Thursday, August 15, 2013

Simple In-Place Swap

Beginner programmers are often given the task of swapping the values of two variables. First, they will try to assign to each variable the value of the other.

x := y
y := x

Of course, this doesn't behave as expected -- the first assignment overwrites the original value of $x$. The standard solution to this issue is to use a temporary variable to duplicate the value of one of the input variables.

temp := x
x := y
y := temp

For the more mathematically-inclined student, there's an extension to this puzzle: swap the variables' values without storing intermediate values.

The trick is to notice that you can hold the original values of $x$ and $y$ implicitly as simple arithmetic results. For example, let $x := x + y$. We haven't lost the original value of $x$ -- we still have it 'stored' as $x - y$. Once this observation is discovered, some trivial arithmetic on a few small examples is sufficient to find the solution.

x := x + y
y := x - y
x := x - y

To verify that this algorithm actually works, we can let $y'$ be the final value of $y$ and $x'$ be the final value of $x$. By substituting previous values into the above equations, we see that $y' = (x + y) - y = x$ and $x' = (x + y) - ((x + y) - y) = y$, confirming that the values have indeed been swapped.

This is a particularly good task to give to children who are learning to program. It introduces to them the idea that programming isn't just about menial bit-twiddling, and reinforces the importance of analytical problem solving in computer science. It's also very approachable yet still conceptually interesting.

Monday, June 17, 2013

Capacity Scaling for Fast Max Flow

In my previous post, I presented an implementation of the Ford-Fulkerson method (termed 'method' and not 'algorithm' because it encompasses a number of different algorithms) for finding the maximum flow through a network.

The Ford-Fulkerson method involves repeatedly finding so-called augmenting paths in the residual network, then saturating these paths, terminating only when there are no augmenting paths left. The order in which these paths are found is dependent on the specific algorithm employed. The Edmonds-Karp algorithm, for instance, uses a breadth-first search to find the shortest augmenting path. Another algorithm employs a priority queue to find the maximum-flow augmenting path.

There is another approach, however, which accomplishes the same goal but attacks the problem from a different direction: the capacity scaling algorithm. Capacity-scaling for max flow has woefully little literature on it. I found mentions of it in CS lecture slides from several universities, and also in a TopCoder algorithms tutorial. Apparently it's presented in Network Flows: Theory, Algorithms, and Applications, but I'd neven read about it before. In practice, capacity scaling is just as simple to conceptualise, just as easy to implement, and quite a bit faster than traditional Ford-Fulkerson.

Ford-Fulkerson first finds an augmenting path, then saturates that path. Capacity scaling turns this principle on its head: start by assuming that an augmenting path exists with residual flow at least $r$, then find such a path and add $r$ to its flow. (this only works for integral edge weights by the way!)

In particular, we will begin by choosing a value of $r$ such that $r$ is the largest power of 2 that is less than the maximum edge weight in the graph. We will then push $r$ units of flow out of the source node, and hope that all $r$ units reach the sink. If they do, we try push another $r$ units out. Otherwise, set $r := r/2$ and repeat. Iterate while $r \geq 1$.

max_flow = 0
r = 2 ^ max_edge_weight
while r >= 1:
  do:
    push r units from source to sink
    max_flow += r
  while (successful)
  r /= 2
output max_flow

The algorithm should be fairly intuitive. It's clear that it finds the maximum flow, because it terminates only when it cannot push a single unit of flow i.e. there are no augmenting paths left. Furthermore, it pushes flow down an edge at most $\left \lfloor \lg U \right \rfloor$ times, where $U$ is the maximum edge weight in the graph: because each pass considers a different power of 2, one can interpret the algorithm as effectively 'setting each bit' of the flow amount for each edge.

Now to derive the time complexity. First, I want to make a claim that will help us out. I try to avoid proofs, but this one is important so I'll sketch it out.

EDIT: I'm no longer convinced of this proof's correctness. The problem is that by pushing negative flow, you can effectively undo work you did on a previous iteration. I'll leave it here anyway for readers to pick holes in.

Claim: For a graph $G(V,E)$ the algorithm retains a single value of $r$ for at most $|E|$ iterations of the flow-pushing routine.

Proof: Consider the initial case $r=2^{ \left \lfloor \lg U \right \rfloor}, U=\max \; \{ W(e) \mid e \in E \} $, where by design, no edge has residual flow more than twice $r$. Each time we successfully push flow, we decrease the residual flow of at least one edge $e$ by $2^{ \left \lfloor \lg U \right \rfloor}$. Now $W(e)-2^{ \left \lfloor \lg U \right \rfloor} \leq U - 2^{ \left \lfloor \lg U \right \rfloor}$ (since $U$ is chosen to be the greatest of all $W(e)$'s, and furthermore $U - 2^{ \left \lfloor \lg U \right \rfloor} \leq 2^{ \left \lfloor \lg U \right \rfloor - 1}$. So $e$ now has weight $2^{ \left \lfloor \lg U \right \rfloor - 1} < 2^{ \left \lfloor \lg U \right \rfloor}$ so we cannot push $r$ more units of flow down $e$. This process can only occur once for each edge, so the value of r is retained for at most $|E|$ iterations. Following this, we let $r:=r/2$ and classify edges as follows: those that we augmented on the previous iteration (which now have residual flow not more than twice $r$), those that we could not traverse on the previous iteration (because their values were not more than twice $r$), and those which were unreachable on the previous iteration (either because they're not in the connected component or because they can only be reached by first traversing an edge of residual flow not more than twice $r$ -- and hence can be thought of as having the residual flow of this earlier edge). Observe that this scenario is identical to our initial state, in which all reachable edges have $W(e) \leq 2r$, so by induction the algorithm retains a single value of $r$ until the base case $r < 1$.

Since there are $\left \lfloor \lg U \right \rfloor$ values of $r$ to trial, at most $|E|$ repetitions of the push-routine for each value of $r$, and each push-routine conducts an $O(E)$ search (BFS works well), we arrive at a final complexity of $O(E^2 \log U)$. This outperforms Edmonds-Karp's $O(VE^2)$ and the traditional iterated DFS both in theory and in practice.

Code will follow shortly.


Tuesday, May 28, 2013

More like "hacks-imum flow"! Hahaha!

So I have a hunch that the Python maximum flow implementation on Wikipedia is implemented incorrectly. Their find_path method is as follows:

    def find_path(self, source, sink, path):
        if source == sink:
            return path
        for edge in self.get_edges(source):
            residual = edge.capacity - self.flow[edge]
            if residual > 0 and not (edge,residual) in path:
                result = self.find_path( edge.sink, sink, path + [(edge,residual)] ) 
                if result != None:
                    return result
The following condition prevents the depth-first search from getting itself stuck in cycles, as it searches over the residual network:

            if residual > 0 and not (edge,residual) in path:
That looks fine. However, consider the case where your residual network is a dense directed acyclic graph with all edge weights positive. Furthermore, the source node is the 'root' of this DAG (it's an ancestor of every other node), and there is exactly one path from the source to the sink. The number of outgoing paths from the source will increase exponentially in proportion to the number of nodes!

The following graphic illustrates such a case. The source is A and the sink is B. All black edges have residual flow 1, and the red edge has residual flow 0. As the number of horizontal and vertical 'layers' in the graph increases and the number, the number of cycle-free paths from A to X increases exponentially, and none of these paths lead to the sink. So the algorithm could potentially keep trying these faulty paths before ever stumbling upon a path that leads to B.


Of course, I've almost certainly overlooked something. After all, Wikipedia doesn't have any mistakes, right?
As a plea for your forgiveness, I'll post my max-flow code. It doesn't have any of that silly 'object-oriented' stuff. The network example I use is the one found here.
import copy

N = 7
adj = [[0]*N for _ in range(N)]
adj[0][1] = 3
adj[0][2] = 1
adj[1][3] = 3
adj[2][3] = 5
adj[2][4] = 4
adj[3][6] = 2
adj[4][5] = 2
adj[5][6] = 3
source = 0
sink = 6

# the residual network begins as a copy of the original network
residual = copy.deepcopy(adj)

def get_augpath():
    seen = [False] * N

    def dfs(cur):
        if seen[cur]: return
        seen[cur] = True

        for next in xrange(N):
            if residual[cur][next] > 0:
                if next == sink:
                    # cur has an edge directly to the sink
                    return [cur, next]
                else:
                    # suffix is a path from cur to sink, or [] if none exists
                    suffix = dfs(next)
                    if suffix:
                        return [cur] + suffix
        # no paths found
        return []
    return dfs(source)

def max_flow():
    mf  = 0

    path = get_augpath()
    # keep iterating while a path exists
    while path:
        # get the min-cost edge along our augmenting path
        flow_to_push = min(adj[path[i]][path[i+1]] for i in xrange(len(path)-1))
        mf += flow_to_push
        for i in xrange(len(path) - 1):
            residual[path[i]][path[i+1]] -= flow_to_push
            residual[path[i+1]][path[i]] += flow_to_push
        print 'pushed', flow_to_push, 'along', path

        # try find another augmenting path
        path = get_augpath()
    return mf

print max_flow()

Monday, May 20, 2013

Turtle Tower

Another nice question I encountered recently!
You want to build a tower of turtles. Each turtle has a weight and a strength. A turtle tower is stable if each turtle in the tower can support all the turtles above it; that is, if each turtle's strength is at least as large as the total weight on top of it. How many turtles are in the tallest stable tower you can make?
Dynamic programming problems jump out at you like crazed carnivorous turtles, and this is no exception. It looks like the longest increasing subsequence problem, with a hint of knapsack. The first step would be to establish an ordering of turtles such that each turtle in the ordering can be placed after all other turtles, then we'd use some simple dynamic programming to select some maximal subset!

But hang on. What's the ordering? If the turtles are $(s=5, w=5)$ and $(s=10,w=10)$ then the heavier turtle must go below. If the turtles are $(s=2,w=10)$ and $(s=15,w=5)$ then the lighter turtle must go below! If the turtles are $(s=5,w=5)$ and $(s=10,w=2)$ then the stronger turtle must go below; but if the turtles are $(s=10,w=2)$ and $(s=5,w=15)$ then the weaker turtle must go below! Hmm...

It turns out that you have to sort turtles by the sum of their weight and strength. I don't know the intuition behind why this works, but it does; and proving it is an interesting and surprisingly involved exercise.

Denote the strength and weight of turtle $T$ by $s_T$ and $w_T$, respectively.

Lemma: If a stable tower can be built of height $h$, then a stable tower can be built of height $h$ such that for any pair of turtles $A,B$ where $B$ is below $A$, then $s_B+w_B \geq s_A+w_A$.

Proof: The proof will consist of an exchange argument. We will show that for any adjacent pair of turtles in the tower, we can always arrange them so that the lower turtle has a higher strength plus weight than the upper turtle, while still maintaining the tower's stability. If this is possible, then by a 'bubble sort'-type argument, we can repeatedly swap out-of-order adjacent pairs until the tower is sorted by each turtle $t$'s $s_t+w_t$.

Consider a tower consisting only of a pair of turtles $A,B$. Our task is to prove that if $B$ is placed below $A$, then $s_B+w_B \geq s_A+w_A$.

  • Case 1: assume $B$ supports $A$ but $A$ does not support $B$; that is, $s_B\geq w_A$ and $w_B > s_A$. Adding these inequalities, we obtain $s_B+w_B>s_A+w_A$ as required. Easy peasy!
  • Case 2: assume $B$ supports $A$ and vice versa. Regardless of their ordering, the total weight of their tower will obviously be $w_A+w_B$. Furthermore, if $B$ is placed below $A$, the weight supported by the tower is given by $\min(s_A,s_B-w_A)$ and if $A$ is placed below $B$, then the weight supported by the tower is given by $\min(s_B, s_A-w_B)$. We wish to maximise the weight the tower can hold above it, and we can prove that if by placing $B$ below $A$ where $s_B+w_B\geq s_A+w_A$, this is accomplished; that is, we can show that $\min(s_A,s_B-w_A)\geq \min(s_B,s_A-w_B)$. We will work backwards from the statement we wish to prove until we reach a true statement (not the most elegant proof but whatever).
    \[\begin{aligned} \min(s_A,s_B-w_A) &\geq \min(s_B,s_A-w_B)\\ \min(s_A,s_B-w_A) + w_A + w_B &\geq \min(s_B,s_A-w_B) + w_A + w_B\\ \min(s_A+w_A+w_B,s_B+w_B) &\geq \min(s_B+w_B+w_A,s_A+w_A)\quad \text{(by distributivity of min)}\\ \min(s_A+w_A+w_B,s_B+w_B) &\geq s_A+w_A\quad \text{(by our assumption } s_B+w_B\geq s_A+w_A\text{)} \end{aligned}\] Clearly, both $s_A+w_A+w_B$ and $s_B+w_B$ are greater or equal to $s_A+w_A$. From our premise, we've reached a true statement, therefore the premise ("$B$ should be placed below $A$") is true!
These two cases are exhaustive for all situations where $B$ is placed below $A$, so our proof is complete!

The full solution is left as an exercise to the reader :) it can be be solved in $O(n^2)$ time and $O(n)$ space.

EDIT: after years, I've finally seen a satisfying explanation of the sort order. Suppose turtle $i$ supports turtle $j$ supports some stack of turtles. The total weight of the stack is no more than $s_i + w_i$. If $s_j+w_j > s_i+w_i$ then $s_j > (s_i + w_i) - w_j$, i.e. turtle $j$ supports the total weight of the stack excluding its own weight. Hence we can place turtle $j$ at the bottom of the stack.

Sunday, May 5, 2013

N Coloured Blobs

N coloured blobs wander around in groups on an open plain. Each group has its own distinct colour, and all coloured blobs in a group share the same colour. Occasionally, two groups will bump into each other: each of the blobs in the smaller group will take on the colour of the larger group, and the two groups will merge. If two groups of equal size bump into each other, the ’mergee’ is chosen at random.

A researcher monitors a coloured blob population. Initially, each blob is on its own (i.e. in a group of size 1) and has its own distinct colour. After a sufficient amount of time has passed, all the blobs have become part of one huge group of the same colour. Let K be the total number of times that a blob changed its colour.

A lower bound on K is N-1, if every blob is merged individually into the huge group, hence only changing colour once: from its initial colour to the colour of the huge group. Determine an upper bound on K that is as tight as possible.

This is a sweet little problem with a not terribly difficult solution, but it involves some interesting complexity analysis and allows us to derive a beautiful new implementation for the disjoint-set data structure.

Let's follow a single blob for its lifecycle and consider the number of times that it changes colour. Assume that at every merge, it merges with blob-set of the same size and takes on the other blob-set's colour.

After the first merge, its blob-set is size 2. After the second merge, its blob-set is size 4. We extrapolate to see that after the $k$th merge, its blob-set is size $2^k$; so if there were $N$ blobs initially, then a single blob can experience at most $\log_2 N$ merges and thus at most $\log_2 N$ colour changes. To give our final upper-bound, we assume that each blob experiences this number of colour changes, meaning that in the order of $O(N \log N)$ colour changes occur overall. (In fact, we can reduce the upper bound to $\frac{N \log_2 N}{2}$, but being a computer scientist, I'm not concerned about constant factors. seanieg89 from the boredofstudies.org forum proved this upper bound with some fancy maths, and that proof is available here.)

This demonstrates an interesting implementation of the disjoint-set data structure that I hadn't ever considered before. Give each set a label and keep track of the size of each set. When it comes time to merge two sets, overwrite the labels of the smaller set with those of the larger set. This is worst-case time $O(n \log n)$, as opposed to the union-find implementation which is marginally faster: $O(n \cdot a(n))$, where $a(n)$ is the inverse Ackermann function. However, in many of the places where we'd normally use disjoint sets, like Kruskal's algorithm, we need to sort stuff proportional to $n$ anyway, meaning the final complexity of the algorithm becomes $O(n \log n)$ anyway.

So our blob-merging operation, which has virtually no overhead and so is blisteringly fast, can actually replace the standard union-find algorithm in many situations! I've seen plenty of olympiad questions which require fast disjoint sets, but have additional complications which mean that the union-find algorithm doesn't work very nicely. In these situations, it might be preferable to merge some blobs as above.

Monday, April 1, 2013

Pirate Gold

Been super busy, haven't updated in ages! I feel very guilty. Here's a classic interview problem, plus my solution. (my regards to my algorithms lecturer, Aleks Ignjatovic, for the write-up)
There are five pirates who have to split 100 bars of gold. They all line up and proceed as follows:
  1. The first pirate in line gets to propose a way to split up the gold (for example: everyone gets 20 bars)
  2. The pirates, including the one who proposed, vote on whether to accept the proposal. If the proposal is rejected, the prate who made the proposal is killed.
  3. The next pirate in line then makes his proposal, and the 4 pirates vote again. If the vote is tied (2 vs 2) then the proposing pirate is still killed. Only majority can accept a proposal. The process continues until a proposal is accepted or there is only one pirate left. Assume that every pirate :
  4. above all wants to live;
  5. given that he will be alive he wants to get as much gold as possible;
  6. given maximal possible amount of gold, he wants to see any other pirate killed, just for fun;
  7. each pirate knows his exact position in line;
  8. all of the pirates are excellent puzzle solvers.
Like many other problems of the sort, it's helpful to start by considering a very small case -- for example, when $n=2$.

Call the five pirates $A,B,C,D,E$. The leftmost pirate in the line goes first.

2 Pirates: $D,E$.

$D$ proposes. $D$ votes for himself.
$E$ can always vote against $D$, resulting in a tie, meaning $E$ gets all the money.
$D$: dead
$E$: 100

3 Pirates: $C,D,E$.

$C$ proposes. $C$ votes for himself.
$D$ knows that if $C$ dies, then $D$ will also die in the following round (as in the 2 pirates case). So $D$ will avoid death by voting for $C$ regardless of $C$'s proposition.
$C$ wins with $2/3$ votes. 
$C$: 100 
$D$: 0 
$E$: 0

4 Pirates: $B,C,D,E$.

$B$ proposes. $B$ votes for himself.
$C$ can get 100 bars by killing $B$, and gains nothing if $B$ lives, so $C$ will vote against $B$.
$D$ will get nothing if $B$ dies, and as pirates always want to kill pirates for fun, he will kill $B$ unless $B$ gives him some incentive to vote for $B$'s proposition. So if $B$ gives $D$ 1 bar, $D$ will vote for $B$.
Similarly for $E$ -- so $B$ gives $E$ 1 bar in exchange for $E$'s vote. $B$ wins with 3/4 votes.
$B$: 98
$C$: 0
$D$: 1
$E$: 1


5 Pirates: $A,B,C,D,E$.

$A$ proposes. $A$ votes for himself.
$B$ can get 98 bars by killing $A$, and gains nothing if $A$ lives, so $B$ will vote against $A$.
$C$ will get nothing if $A$ dies, and as pirates always want to kill pirates for fun, he will kill $A$ unless $A$ gives him some incentive to vote for $A$'s proposition. So if $A$ gives $C$ 1 bar, $C$ will vote for $A$.
$D$ will get 1 bar if $A$ dies, and as pirates always want to kill pirates for fun, he will kill $A$ unless $A$ gives him some incentive to vote for $A$'s proposition. So if $A$ gives $C$ 2 bars, $C$ will vote for $A$.
$A$ wins with 3/5 votes.
$A$: 97
$B$: 0
$C$: 1
$D$: 2
$E$: 0
Therefore $A$ should propose that $A$ gets 97 gold bars, $C$ gets 1 gold bar and $D$ gets 2 gold bars.

General Case

This is a simple recursive procedure -- the $n$th pirate has to bribe just enough other pirates so that he gets enough votes, while also maximising his own money. This can be done by finding those pirates who would get the least amount of gold bars in the case of $n-1$ pirates, then bribe them with one extra gold bar more than they would otherwise have gotten. Here is a $O(n\log n)$ Python implementation, and its output.
def best_strategy(n):
    if n == 2:
        return [-1, 100]
    else:
        moneys = [100] + best_strategy(n-1)
        copy = sorted(zip(moneys, range(n-1)))
        next = [100] + [0] * (n-1)
        for i in xrange(int((n-1) / 2)):
            their_money, person = copy[i]
            next[person] = their_money + 1
            next[0] -= next[person]
        return next

for i in xrange(2, 15):
    print i, best_strategy(i)

Output

2 [-1, 100]
3 [100, 0, 0]
4 [99, 0, 1, 0]
5 [97, 0, 1, 2, 0]
6 [97, 0, 1, 2, 0, 0]
7 [96, 0, 1, 2, 0, 1, 0]
8 [96, 0, 1, 2, 0, 1, 0, 0]
9 [95, 0, 1, 2, 0, 1, 0, 1, 0]
10 [95, 0, 1, 2, 0, 1, 0, 1, 0, 0]
11 [94, 0, 1, 2, 0, 1, 0, 1, 0, 1, 0]
12 [94, 0, 1, 2, 0, 1, 0, 1, 0, 1, 0, 0]
13 [93, 0, 1, 2, 0, 1, 0, 1, 0, 1, 0, 1, 0]
14 [93, 0, 1, 2, 0, 1, 0, 1, 0, 1, 0, 1, 0, 0]



Monday, September 17, 2012

Finding Circuits in Points on a Plane

The Japanese Olympiad in Informatics Open Contest, held today, was absurdly hard. Of around 140 contestants, two solved Q1. No one solved Q3. 24 contestants, including myself, solved Q2. This post takes a formal approach to solving Q2.
There are N (1 ≤ N ≤ 105) points positioned at integer coordinates on a plane. Starting at any point, find a route that connects all points exactly once via line segments, without any intersections, and returns to the starting point.
Two examples of valid routes through a given set of points. Note that the points will not always be spaced out in a grid.

Breaking Down the Problem

First, let's consider a simpler problem: ignore the requirement of having to return to the start point. All we have to do is connect up all the points with straight lines, without intersections, without "taking the pen off the page". Is it always possible to do this?

Intersections arise when you double back on yourself. If you only moved in one direction -- say, to the right -- then you'd never re-encounter a line you drew, so you'd never create an intersection. If there was a way to connect all points such that your line segments extend only vertically or to the right, then this route would be intersection-less. Clearly, we have to be somewhat clever about how we do this; for example, the case below cannot be completed if we're at position 9.

We can't get to the 3 or 2 without intersecting something.

Concept & Proof

What if we started at the leftmost point, and processed all the points sequentially from left to right? In theory, this would never result in a case such as the one above. Here's a handwavy inductive proof outlining why we can pick all points sequentially:

Lemma: All possible line segments representing connections between points, lie inside the convex hull of those points. This is the inverse of one one of the definitions of the convex hull ("A set of points is defined to be convex if it contains the line segments connecting each pair of its points." -- Wikipedia)
  1. First, we must establish a starting point. If there is a point with a unique minimum x-coordinate, then we can pick that. If there are several points that share a minimum x-coordinate, then we cannot pick a 'middle' point, or else we can't access those above and below without either travelling left later or retracing a vertical line. So we will pick a point on a vertical extremity, the bottom-most left-most point.
  2. Assume we've connected all points to the left and vertically below this point, and none vertically above or to the right. The next sequential point (sorting first by x ascending, then y ascending) is either above us or to the right somewhere. If it is above us, we can always simply draw the next line segment upwards to reach it; because all other line segments only lie inside the convex hull of points vertically below or to the left (by the lemma), there is no risk of intersecting this hull, and thus no risk of intersecting any of the connections inside this hull. By the same logic, if the next point is to the right, we can also connect it. In the new scenario, all points to the left and below are still connected, and none are connected above or to the right, but we've connected one more point.
By this reasoning, all points can be connected from left to right. The first and last points in this ordering will be the first and last points in our route, respectively. However, as the diagram shows, this algorithm often won't allow us to draw a non-intersecting line that returns to the start point. At this point, if you don't have the full solution yet, I encourage you to stop and think about it a bit before moving on.

Connecting the Dots (heh)

To solve the problem, there is one crucial observations to be made:

A circuit can be formed by partitioning the set of nodes into two subsets. One subset represents the points you take while moving right, the other represents those you take while moving left. The 'pivot' point, at which you stop going right and begin going left, can be part of either subset. However, you'd need to make sure that the connections of one subset don't intersect with the connections of another. A simple way to ensure this is to ensure their convex hulls don't intersect (as per the lemma above) -- in other words, partition them into two disjoint convex sets.

This is easily accomplished by simply drawing a line through the plane, and dividing the points into subsets based on what side of the line they're on. The line must be placed so that you can get from the rightmost point in the 'going right' subset to the rightmost point in the 'going left' subset, and also the leftmost point in the 'going left' subset to the leftmost point in the 'going right' subset, without intersecting any other connections.

Failing to fulfil the latter requirement above. No way to get from 7 to 12, or 1 to 6.
This latter requirement is interesting. As the diagram above demonstrates, we need to avoid situations where connections in the 'going right' subset (subset containing the point 1) cannot access the rightmost point in the 'going left' subset (in this case, point 12 by our sort).

The fundamental insight, at least for me, came when I considered what happens when we first use the all-connected-without-circuit algorithm outlined above, then we try to connect the rightmost point to the leftmost point. This final line will usually intersect with many of the previous segments drawn in by the algorithm. However, each of the line segments that intersect with it connect two points, and these points must be on opposing sides of the line. If we were able to reroute these connections so they only connected pairs of points above the line or pairs of points below the line, then no connections would intersect our rightmost-leftmost line.

Now, the inductive proof above tells us that it is possible to connect all the points above the line from left to right, and similarly for the points below the line.

Instead of acting as a connection in itself, the leftmost-rightmost line is acting as a guide. Having connected all the points on top of the line, and also all the points below the line, all we have to do now is connect the top set of points, the bottom set of points, and the leftmost and rightmost points, without creating any more intersections.

The leftmost-rightmost line ensures the convex hulls of the two sets are disjoint -- thus, we don't have to worry about intersections between the top and bottom sets. Consider the potential connection we wish to forge between the set of points above the line, and the rightmost point (that lies on the line). The rightmost element in the set fulfils the condition that all other points in its set are connected, and lie directly below or to the left of it. The rightmost point on the line, being the rightmost point overall, also lies to the right of our rightmost set point. So we can definitely establish a connection between the rightmost point of our above-line set and the rightmost line point without intersecting any segments occuring below the leftmost-rightmost line or any of the segments to the left or directly below our current point in this set; this covers all existing connections. Hence, we can connect our above-line set to the rightmost point on the line without creating any intersections.

An identical argument can be made for the leftmost and rightmost points for the above-line and below-line sets. None of them interfere with each other, due to the disjointedness of their convex hulls. This creates a full circuit, thus solving the general case of the problem.

Here's the 100%-scoring code I wrote for it in the competition:

#include <cstdio>
#include <vector>
#include <utility>
#include <algorithm>
using namespace std;
 
typedef long long ll;
 
typedef pair<int, int> point;
#define x first
#define y second
#define SIZE(z) ((int)z.size())
 
int N;
pair<point, int> points[100005];
 

int pos(point p1, point p2, point p3) {
    double a = p2.y - p1.y;
    double b = -(p2.x - p1.x);
    double c = ((ll)(p2.x - p1.x)*(ll)p1.y - (ll)(p2.y-p1.y) * (ll)p1.x);
    
    if (a*p3.x + b*p3.y + c > 0) return -1;
    if (a*p3.x + b*p3.y + c < 0) return 1;
    return 0; 
}
 

int marked[100005];
int main() {
    scanf("%d", &N);
    for (int i=0 ; i<N ; i++) {
        scanf("%d %d", &points[i].first.x, &points[i].first.y);
        points[i].second = i+1;
    }
 
    sort(points, points+N);
 
    point first_point = points[0].first;
    point last_point = points[N-1].first;
    vector<int> line_points, above_points, below_points;

    bool all_on_line = true;
    for (int i=0 ; i<N ; i++) {

        // Get relative position (above, below, or on line)
        int position = pos(first_point, last_point, points[i].first);
        
        if (position == 0) {
            // Point lies on the line
            above_points.push_back(i);
        } else if (position == 1) {
            // Point lies above the line
            above_points.push_back(i);
            all_on_line = false;
        } else {
            // Point lies below the line
            below_points.push_back(i);
            all_on_line = false;
        }
    }
    if (all_on_line) {
        // All points are collinear. Impossible.
        printf("0\n");
        return 0;
    }
 
    for (int i=0 ; i<SIZE(above_points) ; i++) {
        printf("%d\n", points[above_points[i]].second);
    }
    for (int i=SIZE(below_points)-1 ; i>=0 ; i--) {
        printf("%d\n", points[below_points[i]].second);
    }
    
    return 0;
}

Thursday, May 10, 2012

Quad-tree Range Queries are O(sqrt(S)) Time

I was told my whole informatics career that quad-tree operations were all logarithmically proportional to the size of the tree's range. This is true for insertion and point queries, but it is not true for range queries. In the worst case, a range query will recurse down to 1x1 blocks around the edge of your rectangle -- and there are O(sqrt(S)) of these.

Wednesday, April 18, 2012

Digital Extraction and Rotation

This will be a short one. In the Google Code Jam qualification round, Problem C involved rotating the digits of a number to the left, e.g., 12405 -> 24051 -> 40512 -> 5124.

Python lets us do this with easy string <-> int conversions, but these aren't particularly fast. Instead, we can do it mathematically by extracting digits with % and /. This requires us knowing the number of digits in the number beforehand, but this can be calculated with the log10 function (which is relatively slow, because computers store numbers in binary). In the Code Jam problem the number of digits didn't vary, so I could just calculate the length once then pass it in to my shifting function.

>>> def rot(n, length):
     return n % (10**(length-1))*10 + n / (10**(length - 1))

 
>>> rot(12405, 5)
24051
>>> 
>>> rot(24051, 5)
40512
>>> rot(23829239, 8)
38292392

The % / trick is really great, especially in languages such as C where string <-> int conversions are less straight-forward. It works on two principles: integer division with / by a power of 10 allows us to strip off digits from the right. You can visualise this as moving the decimal point to the left, then deleting everything that follows it. For example, 14231 / 102 gives 142.31, which with int division is 142. Similarly, applying the modulo operator % with a power of 10 allows us to strip digits from the right. 14231 % 102 is the remainder when 14231 is divided by 100, which as we saw before is 31.

Can someone please confirm that mathematical shifting is faster than converting it to a string with str() and manipulating that? And also, that log10 is relatively slow?

Hacky but Effective Gradient Descent

Here's an interesting situation I encountered recently:

Given a L of N (1 ≤ N ≤ 100000) integers, where L0..k is sorted in descending order and Lk..N is sorted in ascending order for some secret k, determine the minimum value in the list. Additionally, this list has the property that:
  • In L0..k, if a < b < c then La-Lb ≥ Lb-Lc.
  • In Lk..N, if a < b < c then Lc-Lb ≥ Lb-La.
You may only look at roughly 1000 elements in the list.

The vague bound on permitted number of lookups is due to the way the original problem was phrased. In the original problem, there was a function with discrete domain and range instead of a list, and the function was slow to evaluate. As long as my program ran in under 3 seconds for some test case, it was marked correct for that test case; knowing the speed of the judging computer, this allowed for roughly 1000 function calls.

To rephrase the problem, the list's values form a sort of discrete parabola (though there are no guarantees on symmetry), and we have to find its minimum.

EDIT: I later found out that ternary search was a correct solution. However, one of my friends used the parabola-like property to do a standard binary search over the difference between successive points: i.e., its derivative, (if such a term is defined for this situation), which is a monotonic function. Kudos to Charlie for that beautiful solution).

After trying to code a couple of inevitably buggy ternary searches, I tried a different approach, based on the notion of gradient descent. The idea here is that we have a function whose minimum we want to find, so we pick some starting point on the slope of the function and sort of 'slide' down the slopes (though really, this 'slide' is more of a 'hop' to another point on the function's slope), stopping once we reach the minimum. This is also related to the notion of simulated annealing . When implementing gradient descent with a continuous function, there are two ways of getting the slope of our current point: either differentiate the function at the point (thus giving an exact result) or just approximate it by picking two close points.

I used the latter approach for this problem. We can get the 'tangent' of some point on the list by comparing two adjacent points. For example, if our current position in the list holds the value 1000 and the subsequent position holds the value 500, we should jump forward by a large amount; but if the subsequent position held the value 990, we're probably quite close to the minimum, so we should jump forward by a small amount.

It turns out that 1000 list lookups is a hell of a lot, so even dodgy heuristics like these work most of the time -- it's a number of lookups vs. speed of convergence trade-off, and since our permitted number of lookups is very generous, we can be confident we'll converge to the correct solution.

Pretty much, all I had to do was implement something like this. Note that this specific code is largely untested.

# this list is small, but the solution works on very large lists as well
L = [123912, 12000, 11000, 9000, 5000, 4600, 4300, 4250, 4230, 4228, 4231, 4235, 4250, 5400, 6900, 10000, 234432, 23423432, 8645632423436]

min_point, min_val = -1, float('inf')

# start at list item 0, but really we could start anywhere
cur = 0

# initialise the jump size to some value that's slightly smaller than the size of L
jump = int((len(L)-1 - cur) / 1.5)

for i in xrange(500):
    # this could theoretically go over the length of the list, but that won't happen
    # due to the pattern of jumps. i don't think you can ever normally get to the
    # rightmost element, but this should really be proved
    adj = cur + 1

    # perform two lookups: one for the current point, one for the adjacent point
    val_cur = L[cur]
    val_adj = L[adj]

    # have we found a new best?
    if val_cur < min_val:
        min_point, min_val = cur, val_cur

    if val_adj < val_cur:
        # point to the right is smaller than point to the left...
        # meaning we are on the left of the min, meaning jump right
        cur += jump
    else:
     # otherwise, jump left
        cur -= jump

    # perform a vague reduction of the jump size that will allow it to reach 1 within
    # 500 steps (this should be dependent on list size, but this is all a huge hack).
    # just make sure it doesn't get to 0 prematurely, we always want to allow adjacent jumps!
    jump = max(0, int(jump * 0.8))

# and, magic!
print min_point, min_val

Pretty cool, huh? This is a very crude version of gradient descent, but hey -- it works, and it's surprisingly effective. It found the precise minimum value in all ~30 cases, including ones that were specifically designed to break heuristic programs like this. Sometimes you've just gotta hack something up.

EDIT 2: Here's an implementation of Charlie's algorithm:

# this list is small, but the solution works on very large lists as well
L = [123912, 12000, 11000, 9000, 5000, 4600, 4300, 4250, 4230, 4228, 4227, 4231, 4235, 4250, 5400, 6900, 10000, 234432, 23423432, 8645632423436]
 
min_point, min_val = -1, float('inf')
 
# start at list item 0, but really we could start anywhere
start, end = 0, len(L) - 1
while start <= end:
    mid = (start + end) / 2
    print start, end, mid
    mid_val, next_val = L[mid], L[mid+1]

    if mid_val < min_val:
        min_point, min_val = mid, mid_val

    if next_val > mid_val:
        # we're on the pos grad. slope
        end = mid - 1
    else:
        start = mid + 1

print min_point, min_val

Friday, February 3, 2012

We Need To Talk About Binary Search.

Although the basic idea of binary search is comparatively straightforward, the details can be surprisingly tricky.
-- Donald Knuth

Why is binary search so damn hard to get right? Why is it that 90% of programmers are unable to code up a binary search on the spot, even though it's easily the most intuitive of the standard algorithms?

  • Firstly, binary search has a lot of potential for off-by-one errors. Do you do inclusive bounds or exclusive bounds? What's your break condition: lo=hi+1, lo=hi, or lo=hi-1? Is the midpoint (lo+hi)/2 or (lo+hi)/2 - 1 or (lo+hi)/2 + 1? And what about the comparison, < or ≤? Certain combinations of these work, but it's easy to pick one that doesn't.
  • Secondly, there are actually two variants of binary search: a lower-bound search and an upper-bound search. Bugs are often caused by a careless programmer accidentally applying a lower-bound search when an upper-bound search was required, or vice versa.
  • Finally, binary search is very easy to underestimate and very hard to debug. You'll get it working on one case, but when you increase the array size by 1 it'll stop working; you'll then fix it for this case, but now it won't work in the original case!

I want to generalise and nail down the binary search, with the goal of introducing a shift in the way the you perceive it. By the end of this post you should be able to code any variant of binary search without hesitation and with complete confidence. But first, back to the start: here is the binary search you were probably taught...

Input: sorted list of elements, query term
Output: the index of the first appearance of the query in the list, or an ERROR value otherwise

I propose an alternative definition: a binary search takes as input a (monotonic) function f(x) and a boolean predicate function p(v), and searches over the finite domain of the function for arguments where the predicate is true for the function's value -- i.e., values x such that p(f(x)) is true.

Based on this definition, here are the definitions of the variants:

Upper-bound: find the maximum argument x such that p(f(x)) is true
Lower-bound: find the minimum argument x such that p(f(x)) is true

The traditional binary search described above is a special case of the general lower-bound search, where f(x) = array[x], the domain of f(x) is the set of integers {0, 1, ..., N-1} (N being the length of array) and p(v) = v ≥ query.

In other words, you're searching for the minimum argument x such that array[x] ≥ query. Hey, this is just what we had before!

Let's face it: this description of binary search isn't very helpful. For example, why use this predicate thing when all we need is a simple array[mid] < query in our binary search?

The advantage of this somewhat convoluted definition comes when either the query is not in the array or it's in the array many times. Say you're searching for the first instance of the number 6 in the following array, using the traditional method:

[1, 1, 2, 4, 5, 5, 5, 6, 6, 6, 6, 8, 10, 10, 11]

It should output 7. Let's try this out...

lo, hi = 0, len(arr) - 1
while lo < hi:
    mid = (lo + hi) / 2
    if arr[mid] >= 6:  hi = mid - 1
    else:              lo = mid + 1
print mid  # should print 7

This looks about right. If mid is greater or equal to 6, it'll keep searching to the left, trying to find smaller values of 6. Otherwise, it'll search to the right of mid. Will mid always be the correct index by the end? I run it... apparently not, it outputs 5. Oh, I know! It's because I used ≥ instead of >! Right now it'll keep searching left after it encounters the first 6. OK, change that to a > and run it again. Now I'm getting 9... oh, it must be because I'm accessing arr[mid] at the end instead of arr[lo]. On the final iteration, arr[mid] would have become too large, but arr[lo] will be just right -- it should be at exactly 7, which is what we want. Hit F5; wtf, 10? Undo that ≥/> change from before but keep the other changes to see if it makes a difference -- nope, now 6? And I haven't even considered the issue of inclusive or exclusive bounds...

No joke, I just messed around with the +'s and -'s and bounds and ≥s and >s and ≤s <s and los and mids and his for about 10 minutes and I can't find a single combination which gets me an answer of 7. This is worse than any other bug because you'll undoubtedly end up in a loop of case-bashing: getting it to work for this case, then finding it fails for another, then fixing it for that case and finding it now fails the original case. This style of debugging never works. Instead, turn off your monitor, grab a pen and paper and plan this out.

OK, are we doing upper-bound or lower-bound search? We're finding the minimum index with the value 6; lower-bound then. What's the predicate? It's a lower-bound search, so we want it to return true while we're larger than our desired index and return false when it's smaller. Easy! p(v) = (v >= 6).

So let's do this right. I now write my own p(v) function, even though the logic is ridiculously simple, and I translate the predicate-based binary search definition into code. Most importantly, I introduce a new variable, the best_so_far variable.

def p(v): return v >= 6

lo, hi = 0, len(arr) - 1
best_so_far = None

while lo <= hi:
    x = (lo + hi) / 2
    if p(arr[x]):
        # we found a potential minimum x, but we should still check to see if any smaller ones work
        best_so_far  = x
        hi = x - 1
    else:
        # the predicate is false, so we need to go right to find true values
        lo = x + 1

print best_so_far 

And it works first time. *whistles*

But the problem definition changes! Your boss tells you that now, you must search for the last occurrence of 6!

But hey, that's cool. Re-evaluate the problem: it's now an upper-bound search. Our predicate must return true for values smaller or equal to 6, but start returning false after we get to 7, so the maximum x that p(f(x)) = true is the index of the last 6.

def p(v): return v <= 6

lo, hi = 0, len(arr) - 1
best_so_far = None

while lo <= hi:
    x = (lo + hi) / 2
    if p(arr[x]):
        # we found a potential maximum x, but we should still check to see if any larger ones work
        best_so_far = x
        lo = x + 1
    else:
        # the predicate is false, so we need to go left to find true values
        hi = x - 1

print best_so_far 

Again, it works first go. Notice that I only changed two things: the predicate function (reversed the sign) and the direction we head when we find a predicate=true (i.e. I swapped the lines lo = x + 1 and hi = x - 1). I did not have to mess with any +'s or -'s. No off-by-ones were introduced during the making of this function.

Notice also that I use lo = x + 1 and not lo = x. Similarly, hi = x - 1 and not hi = x. This is a foolproof way to avoid the nasty infinite loop binary search bug caused by integer division -- it ensures that you're never considering a value of x more than once, so you're always narrowing down your search space by at least 1 each time, hence ensuring termination. The use of max/min_so_far gives us complete control over how we're approaching the solution, meaning that we don't need to mess around trying to work out whether it's lo, hi or mid that contain the return value at the conclusion of the algorithm. I personally find that inclusive bounds work best with this form of binary search, but your mileage may vary. If you use exclusive bounds, I make no guarantees on this strategy's correctness.

Yet again, the requirements change. The list is now in descending order, and you need to find the index of the first item less than 5.

[11, 10, 10, 8, 6, 6, 6, 6, 5, 5, 5, 4, 2, 1, 1]

As usual, you only need to make two decisions. Upper- or lower-bound? It's clearly lower-bound: you're finding the first item that satisfies the predicate. What's the predicate? p(v) = (v < 5). Expected output is 11.

def p(v): return v < 5

lo, hi = 0, len(arr) - 1
best_so_far  = None

while lo <= hi:
    x = (lo + hi) / 2
    if p(arr[x]):
        # we found a potential minimum x, but we should still check to see if any smaller ones work
        best_so_far  = x
        hi = x - 1
    else:
        # the predicate is false, so we need to go right to find true values
        lo = x + 1

print best_so_far 

It prints 11.

Do I expect every programmer to write out a trivial p(v) function for every binary search they write? Of course not. It might help you think about the problem, but it's not required. If you take one thing away from this post, let it be this: in any binary search you ever write, whether it be over a list of strings, or a multi-dimensional space, or over the domain of a function that uses the inclusion-exclusion principle on O(1) cumulative sums of rectangular regions of a ternary predicate mapped over the integer values in a 2D grid (this has happened before), you just need to worry about whether it is upper- or lower-bound and what your predicate is.

BAM. No more bugs in a binary search, ever. You can thank me later.

Tuesday, January 3, 2012

Shortest Superstring Problem

If I give you a list of words, can you find the shortest string that contains all the words?

I can!

It turns out that this problem is equivalent to TSP, which means there are no polynomial-time algorithms to solve it. However, we can tackle it the same way as we tackle TSP: do the dynamic programming solution, which although slow, is the optimal correct algorithm. There's a trick to it though, which caught me out the first time: if one word is a complete subset of another and does not appear at the beginning or the end (as in the case ['abc', 'b']), the edge traversal principle doesn't work. In cases like these, the smaller word must first be removed from the list. You can see this below with ['germanic', 'german', 'germ']. EDIT: And if two words are the same, they'll both be substrings of each other, but we don't want to remove both of them. Thanks to Werner Lemberg for picking this up in the comments.

For once, some neat code. What's going on???

all_words = ['ginger', 'german', 'minutes', 'testing', 'tingling', 'minor', 'testicle', 'manage', 'guilt', 'germanic', 'normal', 'malt', 'german', 'germ']
all_N = len(all_words)

words = []

# remove strings which are substrings of others
for i in xrange(all_N):
    for j in xrange(all_N):
        if i != j and all_words[i] != all_words[j] and all_words[i] in all_words[j]:
            break
    else:
        words.append(all_words[i])

N = len(words)
        
 
# determine the numerical overlap between two strings a,b
def overlap(a, b):
    best = 0
    for i in xrange(1, min(len(a), len(b))+1):
        if b.startswith(a[-i:]):
            best = i
    return best
 
 
cost = [[None] * N for _ in xrange(N)]
# Precompute edge costs
# for every pair of words with indices u,v
for u in xrange(N):
    for v in xrange(N):
        # work out the best compressed concatenation you can make with u then v
        cost[u][v] = len(words[u]) + len(words[v]) - overlap(words[u], words[v])
 
                 
cache = {}
backtrace = {}
                 
# top-down DP
def solve(used, last):
    if (used,last) not in cache:
     
        bestCost = 0
        bestOption = None
 
        # the base case is when used == 0. using no words, the optimal solution is of length 0
        if used != 0:
         
            bestCost = float('inf')
         
            # for each word we can use
            for i in xrange(N):
                 
                # if the word is as yet unused
                if (1 << i) & used:
                 
                    # calc the cost of using it
                    newCost = cost[i][last] + solve(used & ~(1<<i), i)
                     
                    # if we've reached a new best solution, update stuff
                    if newCost < bestCost:
                        bestCost = newCost
                        bestOption = i
         
        # cache stuff
        cache[(used, last)] = bestCost
        backtrace[(used, last)] = bestOption
     
    return cache[(used, last)]
 
 
# run it for all possible starting cases
bestCost = float('inf')
bestOption = None
for i in xrange(N):
    cur = solve(((1<<N)-1) & ~(1<<i), i)
    if cur < bestCost:
        bestCost = cur
        bestOption = i
 
 
# reconstruct the words that were used
soln = []
used = (1<<N) - 1
last = bestOption
while last is not None:
    soln.append(words[last])
    used &= ~(1<<last)
    last = backtrace[(used, last)]
soln.reverse()
 
 
# now compress the words of the solution into the final string
cur = soln[0]
for i in xrange(1, N):
    cur += soln[i][overlap(cur, soln[i]):]
print cur

The output is minutestinglingingermanicminormaltmanageguiltesticle. Fascinating stuff, I know.