# Milking order compatible with maximum number of observations

I'm using some old USA Computing Olympiad (USACO) problems to help teach me. This is the "Milking order" problem:

Farmer John has $N$ cows ($1≤N≤10^5$), numbered $1…N$. He has made $M$ observations about his cows' social structure ($1≤M≤50,000$). Each observation is an ordered list of some of his cows, indicating that these cows should be milked in the same order in which they appear in this list.

Farmer John's observations are prioritized, so his goal is to maximize the value of $X$ for which his milking order meets the conditions outlined in the first $X$ observations. If multiple milking orders satisfy these conditions, Farmer John would like to use the lexicographically smallest one.

### Sample input

4 3
3 1 2 3
2 4 2
3 3 4 1


### Sample Output

1 4 2 3


### Interpretation

Farmer has four cows and would like to satisfy three partial ordering rules, in decreasing order of importance:

1. should milk cow 1 before cow 2 and cow 2 before cow 3 (the first observation)
2. should milk cow 4 before cow 2 (the second observation)
3. should milk cow 3 before cow 4 and cow 4 before cow 1 (the third observation).

The first two observations can be satisfied simultaneously, but Farmer John cannot meet all of these criteria at once, as to do so would require that cow 1 come before cow 3 and cow 3 before cow 1. Therefore, the third rule should be ignored. This means there are two possible orderings: 1 4 2 3 and 4 1 2 3. Output the first, being lexicographically smaller.

This code works, but I don't really know any Python (or any other coding language). But that means I might be picking up some really bad practices just because I'm figuring out whatever I want on my own. Please let me know what you would change in the format of this, etc.

import glob
import copy

def mooSort(currentSort, fullGraphOUT, fullGraphIN, numCows):
""" Sorts the cows given the current list of conditions, described by fullGraphIN and fullGraphOUT. These are dictionaries dictionaries with vertex keys and and nodes either outgoing or incoming vertices"""
noIncoming = [] #keep a list of vertices with no incoming edges
newSort = [] #the topologically sorted vertices
newGraphIN = copy.deepcopy(fullGraphIN)
newGraphOUT = copy.deepcopy(fullGraphOUT)

for k in range(1,numCows+1):
if not newGraphIN[k]:

noIncoming.append(k)
del newGraphIN[k] #use this to keep track of if the graph has edges at the end

while noIncoming: #take each element with no incoming edges and add it to our sort, remove it from the graph
noIncoming.sort() #sorting so we can make it lexicographic as well
m = noIncoming
noIncoming.pop(0)
newSort.append(m)

for k in newGraphOUT[m]:
newGraphIN[k].remove(m)
if not newGraphIN[k]: #if there are no more incoming edges, we put it into noincoming
noIncoming.append(k)
del newGraphIN[k]
if newGraphIN: #there's a cycle in the graph
return False
else: #no cycles
currentSort[:] = newSort #set the current sort to this one
return True

#now do this in a binary search

inputs = glob.glob("*.in")

for input in inputs:
with open(input) as data:
#turn the lines into a list of list ints
inputData=[list(map(int,line.split())) for line in lines]
numCows = inputData
numConditions = inputData
inputData.pop(0)
fullGraphOUT = {} #describe graph by outgoing neighbors
fullGraphIN = {} #describe graph by incoming neighbors
currentSort = [] #keep track of our current sorting
for i in range(1,numCows+1): #initialize the graph to empty
fullGraphOUT[i] = set()
fullGraphIN[i] = set()
for list in inputData:
list.pop(0)
for k in list:
for j in list[list.index(k)+1:]:
if not mooSort(currentSort, fullGraphOUT, fullGraphIN, numCows): #once this gives false, we can no longer sort using all the conditions
print(currentSort)
break


### 1. Review

1. Python's style guide recommends 4 spaces for each level of indentation, and keeping lines to a maximum of 79 characters. It's not compulsory to follow this style, but it makes it easier for other Python programmers to read the code. For example, if you had stuck to 79 characters per line, we wouldn't have to scroll the code horizontally to read it here on Code Review. Also, with shorter lines it would have been easier for you to spot typos like "dictionaries dictionaries".

2. Although the problem description has to do with cows and milking, it should be clear that this is a very thin disguise and you're really just finding a topological order on a directed graph. So there is no need for the docstring to refer to "cows", and it's clearer to keep things abstract by naming the function something like topological_order.

3. It would be simpler if the function returned the topological order, instead of updating a list passed as a parameter. I guess that you implemented it this way because the function already returns True if the graph is acyclic and False otherwise. However, in Python it's easy to return multiple results from a function, by returning a tuple. So the code could become:

if newGraphIN: # There's a cycle in the graph.
return False, None
else: # Graph is acyclic.
return True, newSort


and then the caller could use tuple assignment to give names to the two results:

acyclic, order = topological_order(...)
if acyclic:
# Do something with order.
else:
# Handle cyclic case.

4. The code deletes vertices from the graph when there are no more incoming edges:

del newGraphIN[k]


The reason it does this is in order to "keep track of if the graph has edges at the end". That is, in order to be able to write:

if newGraphIN: #there's a cycle in the graph


But I think it would be simpler to leave the vertices in the graph as you go along, and to test at the end for any remaining edges:

if any(newGraphIN.values()):


This saves a couple of lines of code, and also enables a further simplification, described below.

5. Now that the code is no longer modifying newGraphIN during the construction of noIncoming, the argument numCows is not needed, because it is possible to deduce the vertices from the graph. So where the code has:

for k in range(1,numCows+1):
if not newGraphIN[k]:
noIncoming.append(k)


for v, incoming in newGraphIN.items():
if not incoming:
noIncoming.append(v)


or more simply, using a list comprehension:

noIncoming = [v for v, incoming in newGraphIN.items() if not incoming]

6. At each step the code needs to extract the smallest vertex in the noIncoming list (in order to get the lexicographically smallest topological sort). In order to do so, it sorts the whole noIncoming list. But this is inefficient: if there are $n$ items in the list it takes $O(n\log n)$ to sort them, and another $O(n)$ to pop the smallest item.

This can be improved to $O(\log n)$ per item if noIncoming is represented as a heap. A heap is a data structure that is specialized for efficient repeated extraction of the minimum item, which is what is needed here. Python has a heapq module implementing this data structure and the associated algorithms. In this case we start by turning noIncoming into a heap by calling heapq.heapify:

heapify(noIncoming)


then to extract the smallest vertex from the heap, we call heapq.heappop:

m = heappop(noIncoming)


and to add a new vertex to the heap, we call heapq.heappush:

heappush(noIncoming, k)

7. Since the code never modifies newGraphOUT, there is no need to make a copy.

8. It is not necessary for the function to take both fullGraphOUT (the directed graph) and fullGraphIN (the same graph with all edges reversed), since the latter can easily be computed from the former, taking no longer than it would take to make the copy.

9. I think the names could be improved: I would write graph instead of fullGraphOUT, reversed_graph instead of newGraphIN, order instead of newSort, v and w for vertices instead of m and k, and sources instead of noIncoming (because sources and sinks are the terms from graph theory for vertices with no incoming and outgoing edges respectively).

### 2. Revised code

from heapq import heapify, heappop, heappush

def topological_order(graph):
"""Return the lexicographically smallest topological order of the
vertices in graph (if possible).

graph must be represented as a mapping from vertex to an iterable
of its neighbours.

Returns:
acyclic -- True, if graph is acyclic, otherwise False.
order -- The topological order, if graph is acyclic, otherwise None.

"""
# Copy of graph with all edges reversed.
reversed_graph = {v:set() for v in graph}
for v, outgoing in graph.items():
for w in outgoing:

# Topological order.
order = []

# Heap of sources (vertices with no incoming edges).
sources = [v for v, incoming in reversed_graph.items() if not incoming]
heapify(sources)

while sources:
v = heappop(sources)
order.append(v)
for w in graph[v]:
incoming = reversed_graph[w]
incoming.remove(v)
if not incoming:
heappush(sources, w)

if any(reversed_graph.values()):
# Some edges remain, hence there's a cycle in the graph.
return False, None
else:
return True, order


### 3. Analysis

The milking order problem is a disguise for a problem known as incremental topological sorting or online cycle detection. In this problem you have a directed graph and you are given edges one at a time (or in batches as in the USACO problem) to add to the graph, and after adding any edge (or batch of edges) you need to be able to determine if the graph contains a cycle, or to output a topological order on the vertices of the graph.

When the code in the post is given new edges, it updates the directed graph with those edges, and then calls mooSort (or topological_order in the revised code), which runs a modification of Kahn's topological sorting algorithm on the graph. Kahn's algorithm normally takes time proportional to the number of vertices $N$ plus the number of edges $E$, that is, $O(N + E)$, but because we need the lexicographically smallest topological order, we incur a further penalty for finding the smallest vertex at each step, making it $O(N^2\log N + E)$ (in the original code) or $O(N\log N + E)$ (in the revised code).

The topological order gets re-computed after every batch of edges has been added to the graph, so topological_order gets called $O(M)$ times, making the overall runtime $O(MN\log N + ME)$ (in the revised code). The number of edges is $O(\min(N^2, MN))$, as there can be at most one edge between every pair of vertices, and at most $N$ edges are added in each batch, so the overall runtime is $O(MN(\log N + \min(N, M)))$.

We would like to do better than this, but that's not a simple matter, because the incremental topological sorting problem is still an area of active research, and the best known algorithms are complicated. Take a look at this paper for a survey of the state of the art as of 2008:

• Thank you Gareth! I learned a lot from this. I really liked implementing the heap there and learning about it. There were also several tricks and standard of Python I learned from this, my code is definitely better here and hopefully the next puzzle I work out will have look better than this one from the start! – silverware May 31 '18 at 4:03
def mooSort(currentSort, fullGraphOUT, fullGraphIN, numCows):
""" Sorts the cows given the current list of conditions, described by fullGraphIN and fullGraphOUT. These are dictionaries dictionaries with vertex keys and and nodes either outgoing or incoming vertices"""


You can split a docstring over multiple lines to make it easier to read.

It's not immediately obvious what is meant by "nodes".

  noIncoming = [] #keep a list of vertices with no incoming edges
newSort = [] #the topologically sorted vertices
newGraphIN = copy.deepcopy(fullGraphIN)
newGraphOUT = copy.deepcopy(fullGraphOUT)


Python style is to use lower case words with underscores for variables.

I don't see the benefit of creating new names for the defensive copies. IMO it's safer to reuse the name of the parameter, because that way you can't accidentally modify the original dictionary.

The names xIN and xOUT sound like in- and out-parameters. The comment explains their meaning, but I think there are more transparent names, such as predecessors and successors.

  for k in range(1,numCows+1):
if not newGraphIN[k]:

noIncoming.append(k)
del newGraphIN[k] #use this to keep track of if the graph has edges at the end


It's unnecessary to delete from the graph to track whether it has edges: you can easily test that with len(newSort) == numCows.

In fact, I can't see that newGraphIN is necessary at all. The only thing you care about it how many predecessors a vertex still has.

  while noIncoming: #take each element with no incoming edges and add it to our sort, remove it from the graph
noIncoming.sort() #sorting so we can make it lexicographic as well


I won't repeat Gareth's point about using a heap.

   currentSort[:] = newSort #set the current sort to this one


This was unexpected. The docstring said nothing about currentSort, and I wasn't expecting it to be a form of output. IMO it would be sufficient to return newSort if it's complete or None otherwise.

#now do this in a binary search


The code which follows looks like a linear search, not a binary one. Is this comment a statement of intention which was forgotten when the code gave the right answer?

       lines = data.readlines()
#turn the lines into a list of list ints
inputData=[list(map(int,line.split())) for line in lines]
numCows = inputData
numConditions = inputData
inputData.pop(0)
fullGraphOUT = {} #describe graph by outgoing neighbors
fullGraphIN = {} #describe graph by incoming neighbors
currentSort = [] #keep track of our current sorting
for i in range(1,numCows+1): #initialize the graph to empty
fullGraphOUT[i] = set()
fullGraphIN[i] = set()
for list in inputData:
list.pop(0)
for k in list:
for j in list[list.index(k)+1:]:
if not mooSort(currentSort, fullGraphOUT, fullGraphIN, numCows): #once this gives false, we can no longer sort using all the conditions
print(currentSort)
break


This looks like a bit much for a single loop body. Why not split the input processing into a function?

         for k in list:
for j in list[list.index(k)+1:]: