TC SRM 459 D2-1000 (ParkAmusement)

Category: Dynamic Programming, Graph Theory
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10723

The problem describes having N landings in the amusement park. The landings can be terminating (i.e. the ride ends at this landing) which are either proper exits or crocodile ponds. The landings can also be non-terminating which means that you will slide down further or until you are stalled with no valid adjacent landings. All of the N landings have different heights and we can only slide down from a taller landing to a shorter one (by the laws of gravity). We start at a random landing and we wish to calculate the probability that we started at a specific landing given that we only went through exactly K different pipes to get safely home (i.e. at a proper exit and not at either a crocodile pond or a dead-end landing).

This problem is best split into two parts. The first is to calculate the probability we reach home from a specific starting point. The second is to calculate the probability that out of all the viable starting points we chose the specific starting landing. We are already given the adjacency matrix in terms of a vector/array of strings which we proceed to convert into a more usable format. We first keep track of exits - i.e. whether or not they are proper exits (home) or crocodile ponds. This is done by specifically checking whether or not landings[i][i] is '0' or not (based on the problem definition). If it's not '0', then the current vertex is an exit and we proceed to determine which type. Otherwise, we add the edges into our adjacency list. This is shown below:

vector<vector<int> > adjlist;   // adjacency list
int table[51];   // the exit table for each vertex, 0 if it's not an exit, 1 if it's home exit, 2 if it's a pond exit

...
double ParkAmusement::getProbability(vector <string> landings, int startLanding, 
      int K) {
   adjlist = vector<vector<int> >(landings.size(), vector<int>());
   for (int i = 0; i < landings.size(); i++) {
      if (landings[i][i] != '0') {
         // exit
         if (landings[i][i] == 'P') table[i] = 2;   // pond
         else table[i] = 1; // exit
         continue;
      }
      for (int j = 0; j < landings[i].size(); j++) {
         if (landings[i][j] == '1') adjlist[i].push_back(j);
      }
   }
...

Once we have parsed the input into a better format, we need to find out the probability of reaching home from a fixed starting vertex/landing. We can accomplish this through a DFS-like recursive function by keeping track of how many steps we have left to reach home. This involves keeping track of which vertex/landing we are currently at and how many hops left we must use. The base case is when the number of steps left is zero, if we are at an exit and it's the proper one (i.e. an 'E' in the matrix) then we have reached home - otherwise we have failed in our mission.

We utilise basic probability calculations to compute our chances of succeeding. If we arrive at a landing that has M out-going edges, then the probability of choosing a particular edge is 1/M. Therefore, the probability we will survive is simply the sum of 1/M * f(edge i, K-1) for all out-going edges i. The base cases are obvious, probability of surviving is 1.0 if we reach an 'E' exit otherwise it's 0.0. The function below illustrates this concept:

double calcPr(int startNode, int K) {
   if (K == 0) {
      // check if it's the exit
      if (table[startNode] != 1) return 0.0;
      return 1.0;
   }
   double res = 0.0;
   // note: if there are no more valid out-going edges it defaults to a 
   //    0.0 probability
   int poss = adjlist[startNode].size();
   for (int i = 0; i < adjlist[startNode].size(); i++) {
      res += (1.0 / poss) * calcPr(adjlist[startNode][i], K-1);
   }
   return res;
}

However, we can optimise this by memoising it. Otherwise we may run into trouble with time limits for the larger cases. This is easily done by setting up a cache table and marking it with a sentinel value to determine whether the state has been computed or not. Since no probability can be negative, we can exploit this fact.

double dp[51][51];

double calcPr(int startNode, int K) {
   if (dp[startNode][K] >= 0.0) return dp[startNode][K];
   if (K == 0) {
      if (table[startNode] != 1) return dp[startNode][K] = 0;
      return dp[startNode][K] = 1.0;
   }
   double res = 0.0;
   int poss = adjlist[startNode].size();
   for (int i = 0; i < adjlist[startNode].size(); i++) {
      res += (1.0 / poss) * calcPr(adjlist[startNode][i], K-1);
   }
   return dp[startNode][K] = res;
}

Now, we have finished the first step of the problem. The second step involves calculating out of the many possible landings we have chosen and given that we did arrive home safely, what is the probability we chose startLanding? The easiest way to think of this is to think geometrically. If we let a circle C be the total sum of the probabilities that we arrived safely from all possible landings, we can scale this to 1.0 as we know for certainty that we arrived safely (given in the problem). Now we simply need to calculate how much probability was contributed by starting from startLanding. This is simply calcPr(startLanding, K) / totalSum. This can be visualised as:

double tot = 0.0;
for (int i = 0; i < 51; i++) for (int j = 0; j < 51; j++) dp[i][j] = -1.0;
for (int i = 0; i < landings.size(); i++) {
   tot += calcPr(i, K);
}
return calcPr(startLanding, K) / tot;

Combining all of the above steps, we yield the final code:

class ParkAmusement {
public:
   double getProbability(vector <string>, int, int);
};

vector<vector<int> > adjlist;
int table[51];
double dp[51][51];

double calcPr(int startNode, int K) {
   if (dp[startNode][K] >= 0.0) return dp[startNode][K];
   if (K == 0) {
      if (table[startNode] != 1) return dp[startNode][K] = 0;
      return dp[startNode][K] = 1.0;
   }
   double res = 0.0;
   int poss = adjlist[startNode].size();
   for (int i = 0; i < adjlist[startNode].size(); i++) {
      res += (1.0 / poss) * calcPr(adjlist[startNode][i], K-1);
   }
   return dp[startNode][K] = res;
}

double ParkAmusement::getProbability(vector <string> landings, int startLanding, 
      int K) {
   adjlist = vector<vector<int> >(landings.size(), vector<int>());
   for (int i = 0; i < landings.size(); i++) {
      if (landings[i][i] != '0') {
         // exit
         if (landings[i][i] == 'P') table[i] = 2;   // pond
         else table[i] = 1; // exit
         continue;
      }
      for (int j = 0; j < landings[i].size(); j++) {
         if (landings[i][j] == '1') adjlist[i].push_back(j);
      }
   }
   double tot = 0.0;
   for (int i = 0; i < 51; i++) for (int j = 0; j < 51; j++) dp[i][j] = -1.0;
   for (int i = 0; i < landings.size(); i++) {
      tot += calcPr(i, K);
   }
   return calcPr(startLanding, K) / tot;
}

Applications of Floyd-Warshall's Algorithm

We will expand on the last post on Floyd-Warshall's algorithm by detailing two simple applications. The first is using the algorithm to compute the transitive closure of a graph, the second is determining whether or not the graph has a negative cycle.

Transitive closure is simply a reachability problem (in terms of graph theory) between all pairs of vertices. So if we compute the transitive closure of a graph we can determine whether or not there is a path from vertex x to vertex y in one or more hops. Unlike the shortest path problem we aren't concerned on how long it takes to get there, only whether or not if we can eventually get there. Obviously you can simply not modify our original algorithm and assume if the distance between vertex x to vertex y is at least "infinity" then there is no way to get from vertex x to vertex y, otherwise there is a way to get there. This perfectly solves the reachability problem.

There is cleaner approach of computing the transitive closure using a slight modification of Warshall's algorithm. Instead of letting no edge between x to y be denoted by infinity, we can let it be 0 (or false). Any weighted edge is simply reduced to 1 (or true). Now instead of adding distances we use the binary-AND operation. If graph[x][k] is true and graph[k][y] is also true, then there is a way to connect x to y. Any other combination will result in failing to connect vertex x to vertex y using an intermediate node k. So graph[x][k] & graph[k][y] correctly fits our desired behaviour. As we just need 1 such intermediate node k to connect vertex x to vertex y, we combine it with the binary-OR operation as we don't care which intermediate vertex it requires to get from x to y so long as there is a path that allows us to get there.

The modified algorithm is the following:

// warshall's algorithm for transitive closure
for (int k = 1; k <= N; k++) {
   for (int i = 1; i <= N; i++) {
      for (int j = 1; j <= N; j++) {
         // only need one connection for it to be reachable, hence the 
         //    OR operation.
         // need both adjmat[i][k] and adjmat[k][j] to be connected for i to go 
         //    to j using intermediate node k.
         adjmat[i][j] |= (adjmat[i][k] & adjmat[k][j]);
      }
   }
}
// post-condition if adjmat[i][j] is 1 then there is a path from vertex i to j

We can make an obvious small optimisation by only running the inner loop if adjmat[i][k] is true, otherwise all the inner computations would fail anyway.

// warshall's algorithm for transitive closure - attempt 2
for (int k = 1; k <= N; k++) {
   for (int i = 1; i <= N; i++) {
      // optimise here by reducing the number of inner loops
      if (adjmat[i][k]) {
         for (int j = 1; j <= N; j++) {
            adjmat[i][j] |= (adjmat[i][k] & adjmat[k][j]);
         }
      }
   }
}
// post-condition if adjmat[i][j] is 1 then there is a path from vertex i to j

The second application is using it to detect negative cycles in a graph. How does Floyd-Warshall relate to this? Since the algorithm calculates the shortest path between all pairs of vertices we can use the "shortest path" for vertex i to itself to determine if there is a cycle. If we can reach from vertex i to itself with a cost less than 0 then we can keep looping through to successively generate shorter paths - hence the term negative cycle. Therefore it is sufficient to just check whether or not adjmat[i][i] < 0 for all vertices. This is used in many other algorithms such as the cycle-cancelling algorithm for minimum cost flows.

// warshall's algorithm for detecting negative cycles
bool cycle = false;
for (int k = 1; k <= N; k++) {
   for (int i = 1; i <= N; i++) {
      for (int j = 1; j <= N; j++) {
         adjmat[i][j] = min(adjmat[i][j], adjmat[i][k] + adjmat[k][j]);
      }
      if (adjmat[i][i] < 0) cycle = true;
   }
}

Floyd-Warshall All-Pairs Shortest Path Algorithm

There are many notable algorithms to calculate the shortest path between vertices in a graph. Most are based on single source to a set of destination vertices. Examples of such famous algorithms include Dijkstra's, Bellman-Ford and even breadth first search for weightless graphs. However, sometimes we wish to calculate the shortest paths between all pairs of vertices. An easy way to calculate this is to loop through each possible pair of vertices and run a Dijkstra through it. However, there exists a more efficient way to calculate this in n^3 time complexity.

The basis of this algorithm lies with dynamic programming. As with all dynamic programming algorithms we need to first define the DP state. The key factor in this algorithm is to consider an intermediate node k of which to connect two vertices i and j. We check if the path from i to k and then from k to j improves our current shortest path length from vertex i to j. If it does, we update - otherwise we don't. Using this definition our DP state becomes:

D(i,j,k) = shortest path length from vertex i to j considering a set of intermediate nodes from vertex 1 to k

We don't need to explicitly label our vertices as this is done implicitly through our algorithm. Now we need to define our recurrence relation. For a given instance of D(i,j,k) we can calculate the shortest path by considering all intermediate vertices from 1 to k. The proposed distance is easily calculated as D(i,k,k-1) + D(k,j,k-1) by definition of an intermediate node. Formally,

F(i,j,k) = min { F(i,j,k-1), F(i,k,k-1) + F(k,j,k-1) }

We have to determine whether or not the new intermediate node k actually helps achieve a shortest path by comparing it to the previous computation without considering intermediate node k, i.e. F(i,j,k-1).

As with all recurrence relations, a base case is required. If we don't consider any intermediate nodes at all, then the shortest path between two vertices is simply the direct edge between vertex i and j. In other words, the weight of the (i,j) entry in the adjacency matrix. So now we have fully defined our recurrence relation with the proper base case:

F(i,j,k) = min { F(i,j,k-1), F(i,k,k-1) + F(k,j,k-1) }
F(i,j,0) = w(i,j)

To calculate this in a bottom-up manner is actually quite simple. First we note the order dependency of the recurrence relation. It can be seen that a particular instance D(i,j,k) can depend on D(i,0,k-1) to D(N,j,k-1) where N is the number of vertices. So it's obvious we need to compute the previous intermediate k values first (i.e. 1, 2, .., k-1 before calculating k). The order of i and j does not matter in our case because we would have computed all k-1 (i,j) pairs before. Using this information we yield an extremely elegant solution:

// warshall's algorithm
for (int k = 1; k <= N; k++) {
   for (int i = 1; i <= N; i++) {
      for (int j = 1; j <= N; j++) {
         adjmat[i][j] = min(adjmat[i][j], adjmat[i][k] + adjmat[k][j]);
      }
   }
}

Let's apply this algorithm to a simple graph theory problem.

TopCoder SRM 301 D1-450 (EscapingJail)
URL: http://www.topcoder.com/tc?module=ProblemDetail&rd=9822&pm=6222

The problem states there are a set of prisoners chained up and we need to determine the maximum distance that a pair of prisoners can be from each other. It's obvious that the maximum length that two prisoners can be apart is defined by the shortest chain separating a connected path between the two prisoners. Sound familiar? This is the exact definition of a shortest path problem! We simply need to find the shortest distance between every pair of prisoner and calculate the maximum of all such pairs.

How do we handle cases in which the prisoners are not chained together (i.e. there is no edge between them)? We can set a sentinel value to denote that the prisoners aren't chained, if we set this to 0 it would be wrong because then it would interfere with our minimum distance calculation. A simple option is to set this to a high/infinity value and if we know the distance from prisoner i to prisoner j exceeds or equals this value then we know they are unbounded. For this problem, we need to return a value of -1 for unbounded distance.

The only remaining part of the problem is translating the given string array into the adjacency matrix. This is a simple exercise of converting a character into the suitable value - if it's a space then we set the distance to be infinity since they aren't chained together. All that leaves is to use the Warshall algorithm and compute the maximum of the pairs (and remembering to return -1 if the distance between at least two prisoners is unbounded).

The implementation is as follows:

class EscapingJail {
public:
   int getMaxDistance(vector <string>);
};

#define INF (1 << 27)

int translate(char ch) {
   if (isdigit(ch)) return ch - '0';
   if (islower(ch)) return ch - 'a' + 10;
   if (isupper(ch)) return ch - 'A' + 36;
   return INF;
}

int adjmat[51][51];

int EscapingJail::getMaxDistance(vector <string> chain) {
   int res = 0;
   for (int i = 0; i < chain.size(); i++) {
      for (int j = 0; j < chain[i].size(); j++) {
         adjmat[i][j] = translate(chain[i][j]);
      }
   }
   for (int k = 0; k < chain.size(); k++) {
      for (int i = 0; i < chain.size(); i++) {
         for (int j = 0; j < chain.size(); j++) {
            adjmat[i][j] = min(adjmat[i][j], adjmat[i][k] + adjmat[k][j]);
         }
      }
   }
   for (int i = 0; i < chain.size(); i++) {
      for (int j = 0; j < chain.size(); j++) {
         if (i == j) continue;
         if (adjmat[i][j] >= INF) return -1;
         res = max(res, adjmat[i][j]);
      }
   }
   return res;
}

TC SRM453.5 D1 450 (PlanarGraphShop)

Category: Dynamic Programming, Graph Theory
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10412

This problem asks us to determine the minimum number of simple planar graphs we need to exact use of N gold coins. Each planar graph has an associated cost function based on the number of vertices and edges the graph has, more specifically: cost = (number of vertices)^3 + (number of edges)^2. We should immediately see a knapsack-like algorithm which strongly suggests using dynamic programming. The main problem lies with determining what number of vertices and edges constitutes a valid simple planar graph.

Luckily there is a simple theorem which we can use as an upper bound calculation to determine whether a pairing of (#V,#E) is a valid pair. For graphs that has at least 3 vertices: if e > 3v - 6 then this will imply nonplanarity. Hence, we can deduce that all e <= 3v - 6 has a simple planar graph representation. For graphs that has less than 3 vertices we can simply hand calculate the combinations. Once we know what our valid pairings are we can construct a set of feasible cost moves based on our cost function by iterating through all possible pairings of # vertices and # edges. We then apply a simple DP algorithm to calculate the minimum number of planar graphs needed.

An implementation is shown below:

class PlanarGraphShop {
public:
   int bestCount(int);
};

int dp[50001];
vector<int> moves;

#define INF 1 << 20

int func(int N) {
   if (N == 0) return 0;
   if (N < 0) return INF;
   if (dp[N] != -1) return dp[N];
   int res = INF;
   for (int i = 0; i < moves.size() && moves[i] <= N; i++) 
      res <?= 1 + func(N - moves[i]);
   return dp[N] = res;
}

int PlanarGraphShop::bestCount(int N) {
   set<int> costs;
   costs.insert(1); costs.insert(8); costs.insert(9);
   for (int v = 1; v < 37; v++) 
      for (int e = 0; e <= 3 * v - 6; e++)
         costs.insert(v*v*v + e*e);
   moves = vector<int>(costs.begin(),costs.end());
   memset(dp,-1,sizeof(dp));
   return func(N);
}

Eulerian Circuits and Paths

Problem:
Given an undirected graph G = (V,E) we need to find a path which uses every edge exactly once - we are allowed to visit vertices multiple times to achieve our goal. If the start and ending vertex must be the same, then this is known as an Eulerian Circuit. Otherwise, this is known as an Eulerian path.

Algorithm:
There is a simple existence theorem which we can use to determine whether or not a graph has an Eulerian circuit and/or path. These are given without proofs:

- A graph has an Eulerian circuit if and only if it is both connected (i.e. from a given vertex we can reach all other vertices through some arbitrary path) and every node has an even degree (which is the number of edges a particular vertex has).

- A graph has an Eulerian path if and only if it is both connected and every node except for two has an even degree.

It follows that an Eulerian circuit is a special case of an Eulerian path in which the start and end vertices are the same. We are allowed to have two non-even degree vertices for an Eulerian path as these denote the start and end vertices.

The classic algorithm to solve this problem is called Fleury's Algorithm. Fleury's Algorithm is pretty simple and intuitive:

- Pick any starting vertex.
- Traverse in a DFS-like manner on edges that haven't been traversed yet. The special case is not to traverse a bridge unless that is the only edge remaining (this is because if we traverse through a bridge we can't come back to the other side by definition as we are only allowed to traverse an edge exactly once).
- When the edge is traversed, erase/delete the edge from the graph.
- Recurse on the new vertex and keep repeating until all edges have been traversed.

The problem with Fleury's when translating it to code is detecting bridges. We present another algorithm which is just as simple but resolves this issue implicitly through the use of backtracking. This algorithm is from the USACO training pages which can also be accessed here: http://acm.fzu.edu.cn/reference/Eulerian%20Tour.htm

Basically, we recurse like in Fleury's Algorithm without consideration of bridges. We also construct the Eulerian path/circuit in reverse order by appending book-keeping information post-processing order. For more details about this algorithm - see the above link. You may ask, how does this avoid having to keep track of bridges?

Firstly, let's make some observations. Take a look at the graph below:

Imagine we are at vertex C. We have several possible choices to traverse - either to vertex E (using edge CE), vertex F (using CF) and vertex D (using CD). Note that if we delete the edge we came from, we are essentially creating a subproblem in which we lose one edge and start "again" from the new vertex we just traversed to in the previous step. For example, if we traverse to vertex E, we are left with only two vertices E-B as a subproblem. The interesting thing to note is that the subproblems carry the same properties as our original problem. Namely, all vertices are either even or they are all even except for two.

Now we need to convince ourselves that the order of traversals of our algorithm does not matter as it produces the correct path for all cases. The first case is when our problem (or subproblem) has no bridges, in which case the order does not matter as we will not exclude any vertices/edges being visited later. The second case is when our problem (or subproblem) has bridges. Our graph above has such a bridge (edge CE). If we were to traverse through this edge before visiting the edges CD, CF and DF then we won't have a valid path as we can't go back on edge CE to visit these edges - thus by definition it won't be an Eulerian path.

However our algorithm adds vertices in a backwards manner. If we visit edge CE first, our stack for the subproblem would be: B-E. Then we visit either edge CF or CD (from backtracking to our vertex C call). The second subproblem's stack are both C-F-D (trace it yourself for edges CF and CD). Since we visited edge CE first, our problem's stack becomes B-E-C-F-D. Since we have now deleted edges CD, CE and CF. Vertex C has no neighbours and hence gets added on as well, yielding: B-E-C-F-D-C. Reversing this path we get the Eulerian path of C-D-F-C-E-B.

If we don't visit edge CE first, then we don't get more than 1 new subproblem (i.e. more than one component) as we don't "violate" the bridge (i.e. we visit all the edges on the left side before traversing it). Hence in this case, our algorithm produces the correct answer also. Intutitively this means that if we generate a disconnected component (by using a bridge too early) we actually put these towards the end of the path first and hence as a result their visiting order gets "reversed" implicitly. Therefore, we still generate a valid path and most importantly still visit all the edges, therefore we have an Eulerian path.

So what happens when there are multiple bridges? Luckily our constraint of at most two vertices having odd degrees saves us a lot of trouble. Draw it out yourself by constructing potential counterexamples!

Our particular graph only has an Eulerian path (not a circuit) because there are two edges with odd degrees (namely vertex C and B). However it is trivial to use the same algorithm for a graph with only even degree vertices - just start at any vertex! For graphs such as the one above, we must start our algorithm on one of the odd degree vertices.

Application (USACO Training Pages - Riding the Fences):

Farmer John owns a large number of fences, which he must periodically check for integrity. Farmer John keeps track of his fences by maintaining a list of their intersection points, along with the fences which end at each point. Each fence has two end points, each at an intersection point, although the intersection point may be the end point of only a single fence. Of course, more than two fences might share an endpoint.

Given the fence layout, calculate if there is a way for Farmer John to ride his horse to all of his fences without riding along a fence more than once. Farmer John can start and end anywhere, but cannot cut across his fields (the only way he can travel between intersection points is along a fence). If there is a way, find one way.

This is a simple application of the Eulerian algorithm outlined above:

using namespace std;

int adjmat[511][511];
int visited[511];
int circuit[1111];
int numEdges[511];
int circuitpos;

void euler_circuit(int pos) {
   bool change = true;
   while (change) {
      change = false;
      for (int i = 0; i < 511; i++) {
         if (adjmat[pos][i] > 0) {
            change = true;
            adjmat[pos][i]--;
            adjmat[i][pos]--;
            euler_circuit(i);
         }
      }
   }
   circuit[circuitpos++] = pos;
}

int main() {
   ifstream inFile("fence.in");
   ofstream outFile("fence.out");
   int F; inFile >> F;
   int x, y;
   int startNode = 1 << 20;
   for (int i = 0; i < F; i++) {
      inFile >> x >> y;
      adjmat[x][y]++;
      adjmat[y][x]++;
      numEdges[x]++; numEdges[y]++;
      startNode = min(startNode, min(x, y));
   }
   for (int i = 505; i >= 0; i--) if (numEdges[i] % 2 != 0) { startNode = i; }
   euler_circuit(startNode);
   for (int i = circuitpos-1; i >= 0; i--) outFile << circuit[i] << "\n";
   return 0;
}

TC TCO06 Round 1 D1 350 (SeparateConnections)

Category: Dynamic Programming, Matching in a General Graph
URL: http://www.topcoder.com/stat?c=problem_statement&pm=6095

This problem just asks us to find the maximum number of matches in a general graph (should be a straightforward reduction given the problem constraints). The most efficient way to solve this problem is to use a maximum matching algorithm (like Edmond's). However, this is not straightforward to code and certainly not required for a 350 pointer problem. Given the constraints, one easy way to implement a maximum matching algorithm is to use dynamic programming. We will be using a bitmask to efficiently keep track of the nodes that have already been matched so far.

For those new to bitmasks, here is a very short tutorial on it. As you may know, all decimal integer numbers (base 10) have a binary representation form (base 2). As a result, each bit inside a binary representation can represent up to 1-bit of information (duh). Besides representing numbers we can use each of these bits (0 or 1) as useful book-keeping information - in our case (and in many cases) we use them to mark whether or not we have visited or used a particular resource. So for example, the bit representation 01100 could mean we have used resource 2 and 3 (0-based index). Note that we go from left to right due to the fact that the left-most bit is also the most significant bit (MSB). This bit representation (01100) is simply known as 2^3 + 2^2 = 14 in decimal. Expanding on this, let's say in our problem we have used node 1, node 4 and node 5. This is represented as 2^5 + 2^4 + 2^1 = 32 + 16 + 2 = 50. There are several important bit operations we can use to manipulate and query the bits:

For starters, to get a '1' in the n-th position from the right (remember 0-based index) - you use the shift left operator '<<'. For example, 1 << 5 yields a bit representation equivalent to 0010 0000 (32).

To set a bit, we use the bitwise-OR operator denoted as '|'. The bitwise-OR operator works by applying the OR operation to each of the corresponding bits between two numbers. The OR operation is defined as:

0 OR 0 = 0
0 OR 1 = 1
1 OR 0 = 1
1 OR 1 = 1

For example to set the 5th bit (again 0-index) from the right, we first use the shift left operator to yield the correct number (i.e. 1 << 5). Let's say the bitmask before the operation is defined as 0100 0010 (68). Then applying 1 << 5 which has a representation of 0010 0000 (32) yields:

0100 0010 (68) OR 0010 0000 (32) Equals 0110 0010 (100)

To test whether a bit is set within a bitmask, we make use of the bitwise-AND operator denoted as '&'. The AND operation is defined as:

0 AND 0 = 0
0 AND 1 = 0
1 AND 0 = 0
1 AND 1 = 1

We use a similar method as we did for setting a bit in a particular position, let's say we wanted to test if the 5th bit was set. We use the shift left operator to yield the correct number (i.e. 1 << 5). Then we simply apply the bitwise-AND operator to the bitmask we wanted to test. If our mask was 0110 0010 then it would return 0010 0000. However if our mask was 0100 0010, then it would return 0000 0000 as it has no bits both turned on at the same position. Thus, we can use this to test whether or not a bitmask contains a bit in a particular position, if it the result is positive after the AND operation then it is inside the bitmask, if it's zero then it isn't inside the bitmask.

There are other bitmasking techniques which won't be discussed here as they aren't used for this problem.

Back to the problem, given that there are only 18 nodes at maximum. We can easily create a bitmask that keeps track of which nodes we have currently used up (there are 2^18 = 262144 different combinations). So for the actual DP state, we can iteratively introduce more nodes to the set and choose any adjacent nodes which haven't been used so far. Therefore, the state is D(n,m) where n is the index of which node we are processing and m is the bitmask of which nodes we have already used up. Next, we need to define the recurrence relation as well as any base cases. The recurrence relation can be found by selecting from a set of valid decisions at a particular node.

The two available options are to:
- Not used the node in the maximum matching (in hopes that it increases the available nodes later on)
- Use the node along with another node which is adjacent to it (increases the number of matchings by 1 and decreases the available set)

We define the recurrence as: D(n,m) = max { 1 + D(n+1, m | (1 << k), D(n+1, m) } where k is the set of all unused nodes which are adjacent to n. Take note, we do not consider matching node n, if it's already been matched by the time we reach index n. Then applying a simple memoization we yield a working solution:

class SeparateConnections {
public:
   int howMany(vector <string>);
};

int adjmat[19][19];
int memo[19][1<<19];
int sz;

int func(int idx, int mask) {
   if (idx >= sz) return 0;
   if (memo[idx][mask] != -1) return memo[idx][mask];
   int res = 0;
   if (!(mask & (1 << idx))) {
      for (int i = idx+1; i < sz; i++) {
         if (mask & (1 << i)) continue;
         if (!adjmat[idx][i]) continue;
         res >?= 1 + func(idx+1, mask | (1 << i));
      }
   }
   res >?= func(idx+1, mask);
   return memo[idx][mask] = res;
}

int SeparateConnections::howMany(vector <string> mat) {
   int res = 0;
   for (int i = 0; i < mat.size(); i++) {
      for (int j = 0; j < mat.size(); j++) {
         if (mat[i][j] == 'Y') adjmat[i][j] = adjmat[j][i] = 1;
      }
   }
   sz = mat.size();
   memset(memo,-1,sizeof(memo));
   return res = func(0,0) * 2;
}

However, one can also use a general matching algorithm like Edmond's Blossoming/Shrinking Algorithm to solve the problem. See below for such an implementation (based on Pape and Conradt variation of Edmond's which is incorrect for a certain set of inputs, but it still passed system tests for this problem). No further explanation would be provided for this algorithm as I'll leave this for a later blog entry in which we build up our matching algorithms.

#define MAX_VERTICES 21

int adjmat[MAX_VERTICES+1][MAX_VERTICES+1];
int table[MAX_VERTICES+1][MAX_VERTICES+1];
int mate[MAX_VERTICES+1];
int level[MAX_VERTICES+1];
int grandparent[MAX_VERTICES+1];

bool checkAncestor(int currentNode, int checkNode) {
       int v = grandparent[currentNode];
       if (v < 0) return false;   // not an ancestor
       if (v == checkNode) return true; // ancestor
       return checkAncestor(v, checkNode);
}

int augment(int curNode) {
       // BFS an augmenting path
       memset(level,1,sizeof(level));
       memset(grandparent,-1,sizeof(grandparent));
       queue<int> q;
       q.push(curNode);
       level[curNode] = 0;
       while (!q.empty()) {
               int x = q.front(); q.pop();
               for (int y = 1; y <= MAX_VERTICES; y++) {
                  // not part of an edge in G or seen the odd-vertex in T
                  if (!adjmat[x][y] || !level[y]) continue;   
                  if (mate[y] == 0) {
                     // found an augmenting path
                     int next = 0;
                     while (true) {
                        mate[y] = x;
                        next = mate[x];
                        mate[x] = y;
                        if (grandparent[x] < 0) break;
                        x = grandparent[x];
                        mate[next] = x;
                        y = next;
                     }
                     return 1;
                 } else if (mate[y] != 0 && !checkAncestor(x,y) && 
                         mate[y] != x) {
                     // push the node into the alternating tree T
                     q.push(mate[y]);
                     grandparent[mate[y]] = x;
                     level[y] = 0;
                 }
               }
       }
       return 0;
}

int maxMatching() {
   int res = 0;
   memset(mate,0,sizeof(mate));    // all currently unmatched
   memset(grandparent,-1,sizeof(grandparent));
   // start at the vertex i if it's single
   while (true) {
      int v = 0;
      for (int i = 1; i < MAX_VERTICES; i++) {
         if (!mate[i]) v += augment(i);
      }
      if (v == 0) break;
      res += v;
   }
   return res;
}

int SeparateConnections::howMany(vector <string> mat) {
   for (int i = 0; i < mat.size(); i++) {
      for (int j = 0; j < mat.size(); j++) {
         if (mat[i][j] == 'Y') adjmat[i+1][j+1] = adjmat[j+1][i+1] = 1;
      }
   }
   return maxMatching() * 2;
}

TC SRM447 D1 500 (PeopleYouMayKnow)

The problem statement can be accessed via:
http://www.topcoder.com/stat?c=problem_statement&pm=10580

The problem itself is pretty simple: it asks us to given a friends graph determine the number of friends to remove such that there are no paths from friend A to friend B which has a length less than 3.

The key observation is to remove all friends which have a shortest path between A and B which exceeds 3 as these will never conflict with our objective. To determine intermediate distances we can just use Floyd-Warshall's algorithm to calculate all-pairs shortest paths between friends. Then we can discard any nodes which violate dist[person1][i] + dist[i][person2] <= 3. Note that the problem fixes the constraint to not allow any paths with a length of less than 3 - this allows us to use a simple bipartite matching to determine the minimum cut of the network. For the general problem of having a length less than n, we can simply use a network flow algorithm (which we use for the implementation here).

We can divide the friends into two (one used as an in-node, the other used as an out-node) and run the maximum flow algorithm over the network. This concept is very similar to the one used in TC SRM360 D1 500 (Prince of Persia).

Implementation below:

class PeopleYouMayKnow {
public:
  int maximalScore(vector <string>, int, int);
};

int augment(int cur, int flow, int ret);

int adjmat[201][201];
int used[201][201];
bool visited[201];
int t;

int maxflow() {
  int flow = 0, thisFlow = 0;
  memset(visited,false,sizeof(visited));
  while ((thisFlow = augment(0, INT_MAX-1, 0)) != 0) {
    memset(visited,false,sizeof(visited));
    flow += thisFlow;
  }
  return flow;
}

int augment(int cur, int flow, int ret) {
  if (cur == t-1) return flow;
  visited[cur] = true;
  for (int i = 0; i < t && ret == 0; i++) {
    if (adjmat[cur][i] > used[cur][i] && !visited[i] && 
      (ret = augment(i, min(flow,adjmat[cur][i]-used[cur][i]), 0)) > 0)
      used[cur][i] += ret;
    else if (used[i][cur] > 0 && !visited[i] && 
      (ret = augment(i,min(flow, used[i][cur]), 0)) > 0)
      used[i][cur] -= ret;
  }
  return ret;
}

int distTable[51][51];
int validNodes[51];
#define INF (1<<20)

int PeopleYouMayKnow::maximalScore(vector <string> friends, int A, int B) {
  int res = 0;
  int sz = friends.size();
  for (int i = 0; i < sz; i++)
    for (int j = 0; j < sz; j++)
      if (i == j) distTable[i][j] = 0;
      else distTable[i][j] = friends[i][j] == 'Y' ? 1 : INF;
  for (int k = 0; k < sz; k++)
    for (int i = 0; i < sz; i++)
      for (int j = 0; j < sz; j++)
        distTable[i][j] <?= distTable[i][k] + distTable[k][j];
  // build the network - allowing infinite capacities between 
  // out -> in nodes and restricting capacities between nodes with 1 capacity
  for (int i = 0; i < sz; i++)
    for (int j = 0; j < sz; j++)
      if (friends[i][j] == 'Y')
        adjmat[i+1+sz][j+1] = INF;
  for (int i = 0; i < sz; i++)
    if (distTable[person1][i] + distTable[i][person2] > 3)
      for (int j = 0; j < sz; j++) adjmat[j+1][i+1] = adjmat[i+1][j+1] = 0;
    else
      validNodes[i] = 1;
  for (int i = 0; i < sz; i++)
    if (validNodes[i]) adjmat[i+1][i+1+sz] = 1; else adjmat[i+1][i+1+sz] = 0;
  adjmat[0][person1+1] = INF;
  adjmat[person1+1][person1+1+sz] = INF;
  adjmat[person2+1][person2+1+sz] = INF;
  adjmat[person2+1+sz][2*sz+1] = INF;
  t = 2*sz+2;
  return res = maxflow();
}

TC SRM 447 D2 900 (ImportsList)

The problem statement can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10579

This problem asks us to calculate the minimum number of scripts to be used within an import list. For a given script A, it needs access to a set of scripts S1, S2, ..., Sn. The relationships are transitive, so if script A has access to script B and script B has access to script C then script A also has access to script C. We must choose a combination such that a given script has access to all its required scripts whilst also choosing the minimum number of scripts to be imported.

The first observation to make is that the graph is directed, so access does not flow both ways. In fact, given the current constraints there are no cycles in the graph hence we have a directed acyclic graph (DAG). This gives us a huge hint to topologically sort (or linearise) the graph before doing further operations. A topological sort can be achieved with a simple DFS and recording the finished nodes in reverse order. Alright let's assume we now have a linearised DAG - how do we choose the minimum amount? Recall that for an DAG we can order the vertices such that all edges go from left to right. Therefore for a vertex V, we can greedily consider vertices from left to right and keeping track of which vertices can be accessed if we access vertex W.

This works for several reasons. Let's consider two vertices: W and X. Let's say that vertex X is to the right of vertex W. If vertex X is reachable from vertex W then vertex W can access at least as much as vertex X since we can simply access vertex X from W using an arbitrary set of traversals. Let's say that vertex X is not reachable from vertex W, then there could be some scripts which we can access between vertex W and X which we need to fulfill the requirements. If there isn't, then we can simply not use vertex W and we will consider vertex X eventually. Therefore, the optimality follows through the linearised order since we consider the largest set and choosing so does not cause conflicts later on.

A very simple way to represent what is accessible from a given vertex is to use a bitmask. Then the implementation is as follows:

class ImportsList {
public:
  vector <int> importsCount(vector <string>);
};

#define N 52
int visited[N];
int order[N];
int cycleDetect;
int adjmat[N][N];
int dfsIdx;

void dfs(int node) {
  if (visited[node]) { if (visited[node] == 1) cycleDetect = 1; return; }
  visited[node] = 1;
  for (int i = 0; i < N; i++) {
    if (adjmat[node][i]) {
      dfs(i);
    }
  } 
  visited[node] = 2;
  order[dfsIdx] = node;
  dfsIdx--;
}

long long mask[N];
long long procMask[N];

vector <int> ImportsList::importsCount(vector <string> requires) {
  vector<int> res(requires.size(),0);
  for (int i = 0; i < requires.size(); i++) {
    long long m = 0;
    for (int j = 0; j < requires.size(); j++) {
      if (requires[i][j] == 'Y') { adjmat[i][j] = 1; m += (1LL << j); }
    }
    mask[i] = m;
  }
  dfsIdx = requires.size()-1;
  for (int i = 0; i < requires.size(); i++) {
    if (!visited[i]) {
      dfs(i);
    }
  }
  // greedy based on topological sorting
  for (int i = requires.size()-1; i >= 0; i--) {
    long long m = mask[order[i]];
    long long p = m;
    int num = 0;
    for (int j = i+1; j < requires.size() && m != 0; j++) {
      if (adjmat[order[i]][order[j]] && (m & (1LL << order[j]))) {
        num++;
        p |= procMask[order[j]];
        m = m & ~(procMask[orderj]]);
      }
    }
    procMask[order[i]] = p;
    res[order[i]] = num;
  }
  return res;
}

Bipartite Matching Algorithm

The problem:

Given a graph G=(V,E) where G is a bipartite graph. Then G can be separated into two disjoint vertex sets A and B where there exists no edges (u,v) where both u and v belong to the same vertex set. We define a matching M to be a set of edges from A to B where each vertex V belongs to at most one such edge. Then the maximum matching refers to the highest possible cardinality of M following the rules we defined above.

Applications:

We can define X to represent a set of workers and Y to be a set of jobs. There is an edge (u,v) in E iff worker u in X is sufficiently skilled enough to work in job v in Y. If we were tasked to assign the maximum number of workers to exactly 1 job then the answer is simply the maximum matching of the graph G.

Algorithm:

Define an alternating path P as a path which alternates between the two sets of vertices A and B. More specifically the alternating P follows the form of: (u0, v0), (v0, u1), (u1, v1), (v1, u2) etc. where each uX is in A and vY is in B. An augmenting path is an alternating path with each (uX, vX) edge being matched (let's define this to be colour blue) and each (vX, u(X+1)) edge being unmatched (define this to be colour red). Recall that one vertex can only be matched to one edge. One special note about an augmenting path is that it ends with an unmatched (uX, vX) edge and we can simply alternate the matchings along the path to obtain exactly 1 more matching in the graph.

For example, take:

A --- B ___ C --- D ___ E --- F ___ G --- H

Where --- denotes an unmatched edge and ___ denotes a matched edge. If H isn't in the matching M then we can alternate:

A ___ B --- C ___ D --- E ___ F --- G ___ H

to obtain 1 more matching (3 vs 4 matchings).

Then our algorithm is to keep finding augmenting paths and improving the number of matchings iteratively. When there are no more augmenting paths we are done.

Theorem: If there exists no augmenting paths then the number of matchings is optimal.

(Informal) Proof: Assume there are no augmenting paths in M. To increase the number of matchings we need to match a new edge (that is currently coloured red). As there are no augmenting paths there exists no red edges from a vertex x in A to a vertex y in B where vertex y is currently unmatched (otherwise there would be an augmenting path). Now if we introduce a new matching into M, we will revert the changes done by an augmenting path:

... C ___ D --- E ___ F --- G ___ H

Matching F and G would force:

... C ___ D --- E --- F ___ G --- H

The subset is the same as our previous alternating path which wasn't augmented. Therefore, introducing a new matching with no augmenting path present does not increase the total number of matchings. Therefore, the absence of augmenting paths means the number of matchings is optimal.

Other Notes:

To find the actual augmenting path, a simple DFS will suffice. A simple way to implement this is to keep track of the colours of each edge (inside the adjacency matrix where red colour = 1, blue colour = 2) as well as keeping track of which side we are on to know which colour edges we need to traverse. We only traverse red edges when we are on the left side and we only traverse blue edges when we are on the right side. Note that if there are no blue edges then we have found an augmenting path (as the right side vertex is unmatched).

Another note is that the number of matchings never goes down so once a vertex is matched in iteration i, then it is also matched in iteration i+1 - although the actual vertex it is matched to can change. Furthermore, we can derive which vertices are matched together simply by looking at the matching array. For vertex i (on the left side) it is matched to matching[i] on the right side, -1 if it isn't matched. For more details, see the implementation below (in the class used for SPOJ TAXI).

The maximum flow algorithm can be reduced to the bipartite matching problem, assign a super source to the set of vertices in A and assign a super sink to the set of vertices in B. Let each edge between A and B have a capacity of 1, and each edge between the source and A as well as B and the sink to have a capacity of infinity. Then the maximum matching is simply the maximum flow from the source to the sink.

Applying the Algorithm (SPOJ TAXI):

To apply the matching algorithm, let's take an easy example. This problem statement can be accessed via:
http://spoj.pl/problems/TAXI

The problem asks us to assign a taxi to exactly 1 passenger. There's a time limit where each taxi must reach the passenger by, this criterion is the basis on whether there is an edge between taxi i and passenger j. We asked to find the maximum number of matchings we can make between the taxi to the passengers. All we need to do is construct the bipartite graph and run the matching algorithm to print out the results.

The implementation is shown below:

#include <iostream>
#include <string>
#include <algorithm>
#include <sstream>
#include <vector>

using namespace std;

class BipartiteMatching {
private:
 vector<vector<int> > adjmat;
 vector<int> matched;
 vector<int> visited;

 int N;
 int M;

 bool augment(int node, int side) {
   visited[side*N+node] = 1;
   if (!side) {
     // left side
     for (int i = 0; i < M; i++) {
       if (adjmat[node][i] == 1 && !visited[i+N]) {
         // check for right-side unmatched node
         if (matched[i+N] == -1) {
           adjmat[node][i] = 2;
           matched[i+N] = matched[node] = i;
           return true;
         }
         // otherwise dfs more
         bool f = augment(i, 1-side);
         if (f) {
           // augmenting path found
           adjmat[node][i] = 2;
           matched[i+N] = matched[node] = i;
           return true;
         }
       }
     }
     return false;
   } else {
     // right side
     for (int i = 0; i < N; i++) {
       if (adjmat[i][node] == 2 && !visited[i]) {
         bool f = augment(i, 1-side);
         if (f) {
           // augmenting path found
           adjmat[i][node] = 1;
           return true;
         }
       }
     }
     return false;
   }
   return false;
 }

public:
 BipartiteMatching(int _N, int _M) : N(_N), M(_M) {
   // initialise the data structures
   adjmat = vector<vector<int> >(N, vector<int>(M, 0));
   visited = vector<int>(N+M,0);
   matched = vector<int>(N+M,0);
 }

 void Reset() {
   for (int i = 0; i < adjmat.size(); i++) {
     fill(adjmat[i].begin(),adjmat[i].end(),0);
   }
 }

 void AddEdge(int i, int j) {
   adjmat[i][j] = 1;
 }

 int maxMatching() {
   int res = 0;
   bool valid = true;
   // all vertices are unmatched at the start
   for (int i = 0; i < N + M; i++) matched[i] = -1;
   // while there is at least 1 valid augmenting path keep
   // iterating
   while (valid) {
     valid = false;
     for (int i = 0; i < N; i++) {
       if (matched[i] == -1) {
         fill(visited.begin(),visited.end(),0);
         bool f = augment(i,0);
         if (f) valid = true;
       }
     }
   }
   // count the number of matches
   for (int i = 0; i < N; i++) {
     if (matched[i] != -1) { res++; }
   }
   return res;
 }
};

class point {
public:
 int x;
 int y;
 point() { }
 point(int _x, int _y) : x(_x), y(_y) { }
};

bool operator<(const point& lhs, const point& rhs) {
 if (lhs.x != rhs.x) return lhs.x < rhs.x;
 if (lhs.y != rhs.y) return lhs.y < rhs.y;
 return false;
}

int main(int argc, char** argv) {
 ios_base::sync_with_stdio(false);
 int numTest = 0;
 cin >> numTest;
 int p, t, s, c;
 BipartiteMatching* bpm = new BipartiteMatching(406,206);
 while (numTest-- && cin >> p >> t >> s >> c) {
   vector<point> people;
   vector<point> taxi;
   // people
   int x, y;
   while (p-- && cin >> x >> y) {
     people.push_back(point(x,y));
   }
   while (t-- && cin >> x >> y) {
     taxi.push_back(point(x,y));
   }
   // make graph
   bpm->Reset();
   for (int i = 0; i < people.size(); i++) {
     for (int j = 0; j < taxi.size(); j++) {
       long long dist = abs(people[i].x-taxi[j].x)+
         abs(people[i].y-taxi[j].y);
       if (dist*200 <= c*s) {
         // can reach
         bpm->AddEdge(i,j);
       }
     }
   }
   // calculate maximum bipartite matching
   cout << bpm->maxMatching() << "\n";
 }
 return 0;
}

TC SRM446 D1 500 (AntOnGraph)

The problem statement can be accessed via:
http://www.topcoder.com/stat?c=problem_statement&pm=10516

The problem depicts an ant which travels around a weighted directed graph. The objective of this problem is to determine the maximum score in which the ant can have given that the score is defined by the sum of all the edge weights it travels. The ant is given at most timeLimit seconds to travel (it can use less than this) and each second is defined by a mandatory set of stepsPerSecond adjacent vertex traversals within the graph. An additional constraint is that the ant must end at vertex 1 at the end (the ant starts at vertex 0), it is also possible that the ant can't reach the vertex at the right interval in time or that the vertex is inaccessible.

Given the constraints it's rather obvious we can't simulate it by time (as it'll easily exceed the time limit). There are several clues in this problem which leads to an elegant solution - firstly, we must move exactly stepsPerSecond traversals each second. What does this mean? It means that we can compute the cost of moving exactly this many traversals per second. Let's assume timeLimit = 1 second, then the answer is simply the greatest value (which could be negative) from vertex 0 to vertex 1 using exactly stepsPerSecond vertex traversals.

How do we compute this? We can think of it in a dynamic programming manner. Let's define M(x,i,j) to be the maximum value we can achieve using exactly x steps from vertex i to vertex j. If this is unreachable then we define M(x,i,j) to be negative infinity. So now we define the relationship between step x and step x-1:

M(x,i,j) = max { M(x-1,i,k) + M(1,k,j) }       
     iff (i,k) and (k,j) are edges (not negative infinity)

We can iteratively compute this from 0 to stepsPerSecond traversals and the maximum score we can get is simply M(stepsPerSecond, 0, 1). Now we need to somehow incorporate this into the case where timeLimit > 1. We can now try and relate timeLimit and timeLimit-1 seconds. It can be seen that the optimal/maximal moves does not change inside each individual timeLimit (the stepsPerSecond matrix). To see this, let's say we are on vertex x and we just finished using the timeLimit-1 interval. We now wish to go from vertex x to vertex y on the next interval. Then the optimal score is simply M(stepsPerSecond, x, y) + currentScore. We simply reused the matrix we computed since it forms a "basic indivisible unit" due to the constraint that each second must have exactly stepsPerSecond vertex traversals.

We now attempt to define a recurrence between seconds (and not specific vertex traversals) using the matrix M as a basic unit:

N(t,i,j) = max { M(t-1,i,k) + M(1,k,j) }
     iff (i,k) and (k,j) are edges

There is a problem though, since the maximum timeLimit is 1 billion - this recurrence is too slow. To optimise it we need to make one observation. With basic expansion:

N(t,i,j) = N(1,i,k) + N(1,k,l) + N(1,l,m) + ... + N(1,o,j)   (t additions)

Then we can begin collapsing:

N(t,i,j) = N(2,i,l) + ... + N(2,n,j)   (t/2 additions)

etc.

As you can see we can recursively define N(t,i,j) in terms of additions of smaller values. Note that we can define it much like how you do exponentation squaring:

N(t,i,j) = N(t/2,i,k) + N(t/2,k,j)            iff (t is even)
         = N(t/2,i,k) + N(t/2,k,m) + N(1,m,j) otherwise

With this we can compute the maximal score for exactly "t" seconds. We need to make a minor modification to allow us to compute the correct score in case we don't need to use all of the allotted "t" seconds. Basically we need to add a dummy edge of 0 score going from the last vertex to itself (self loop) so it's a possible decision to terminate early. Given the recurrence and the self-loop modification on the finishing vertex we can simply implement:

class AntOnGraph {
public:
 string maximumBonus(vector <string>, vector <string>, 
   vector <string>, int, int);
};

long long tab[51][51];

#define INF (1LL<<50)

vector<vector<long long> > mat(const vector<vector<long long> 
    >& lhs, const vector<vector<long long> >& rhs) {
 vector<vector<long long> > res = vector<vector<long long> 
     >(lhs.size(), vector<long long>(lhs.size(), -INF));
 for (int i = 0; i < lhs.size(); i++) {
  for (int j = 0; j < lhs[0].size(); j++) {
   for (int k = 0; k < lhs.size(); k++) {
    if (lhs[i][k] > -INF && rhs[k][j] > -INF)
     res[i][j] >?= (lhs[i][k] + rhs[k][j]);
   }
  }
 }
 return res;
}

vector<vector<long long> > pow(const vector<vector<long long> 
    >& vec, long long n) {
 vector<vector<long long> > res;
 if (n == 1) {
  return vec;
 }
 vector<vector<long long> > base = pow(vec, n/2);
 res = mat(base, base);
 if (n % 2 != 0) {
  return mat(res, vec);
 }
 return res;
}

string AntOnGraph::maximumBonus(vector <string> p0, vector <string> p1, 
    vector <string> p2, int stepsPerSecond, int timeLimit) {
 string res;
 vector<vector<long long> > vec = vector<vector<long long> >(p0.size()
   ,vector<long long>(p0.size(),0));
 for (int i = 0; i < p0.size(); i++) {
  for (int j = 0; j < p0[i].size(); j++) {
   vec[i][j] = ((p0[i][j] - '0') * 100 + (p1[i][j] - '0') * 10 + (p2[i][j] 
       - '0')) - 500;
   if (vec[i][j] == -500) vec[i][j] = -INF;
  }
 }
 vector<vector<long long> > M = pow(vec, stepsPerSecond);
 if (M[1][1] < 0) M[1][1] = 0;
 vector<vector<long long> > M2 = pow(M, timeLimit);
 ostringstream oss; oss << M2[0][1];
 res = oss.str();
 return M2[0][1] == -INF ? "IMPOSSIBLE" : res;
}

TC SRM358 D1 1000 (SharksDinner)

The problem statement can be accessed via:
http://www.topcoder.com/stat?c=problem_statement&pm=7834&rd=10768

The problem states that there are sharks that can eat another provided that all their statistics (skill, speed and intelligence) are at least as good as another. We are asked to find the minimum number of sharks that will remain at the end, in other words, this can be viewed as finding the maximum number of sharks being eaten. This problem easily lends itself to a matching problem, if it weren't for the fact that a shark can eat at most 2 other sharks (instead of just one other shark) then this problem is simply a maximum cardinality bipartite matching problem. However, since the constraint is rather low (allowing a shark to eat at most 2 sharks) we can divide a shark into multiple nodes and solve it as a bipartite matching problem.

However we'll solve it using the maximum flow algorithm (which can also be used to solve the bipartite matching) by making a slight modification. For each shark we give it an "in" node and an "out" node and the capacity between these two nodes is exactly 2 (corresponding to the fact that a shark can only eat at most 2). Using this method allows us to scale the maximum number a shark can eat without increasing the number of nodes (beyond the in and out node expansion). One small note we need to be careful of is sharks which can eat each other. To handle this, a simple modification towards matching edges (i.e. shark A can eat shark B) is made - if shark A and shark B has the same statistics then shark A can eat shark B iff it has a lower index. So you resolve the problem of sharks eating each other cyclically.

Now it's just a matter of constructing the network and running the maximum flow algorithm over it. The final answer is simply the number of sharks in the data set minus the ones that got eaten from our matching.

class SharksDinner {
public:
 int minSurvivors(vector <int>, vector <int>, vector <int>);
};

class Shark {
public:
 int size;
 int speed;
 int intel;
 Shark(int a, int b, int c) : size(a), speed(b), intel(c) { }
};

bool operator<(const Shark& lhs, const Shark& rhs) {
 if (lhs.size != rhs.size) return lhs.size < rhs.size;
 if (lhs.speed != rhs.speed) return lhs.speed < rhs.speed;
 return lhs.intel < rhs.intel;
}

int augment(int cur, int flow, int ret);

class node {
public:
 int num;
 int cap;
 node() { }
 node(int n, int c) : num(n), cap(c) { }
};

bool operator<(const node& lhs, const node& rhs) {
 if (lhs.num != rhs.num) return lhs.num < rhs.num;
 if (lhs.cap != rhs.cap) return lhs.cap < rhs.cap;
 return false;
}

#define INF 9999999
int sz; // the size of the graph (vertices)
int cap[202][202]; // capacity matrix
int used[202][202]; // network flow array
bool visited[202];
vector<vector<node> > adjlist;

int maxflow() {
 int flow = 0;
 int thisFlow = 0;
 memset(visited,false,sizeof(visited));
 while ((thisFlow = augment(0, INT_MAX-1, 0)) != 0) {
  memset(visited,false,sizeof(visited));
  flow += thisFlow;
 }
 return flow;
}

int augment(int cur, int flow, int ret) {
 if (cur == sz-1) return flow;
 visited[cur] = true;
 for (int k = 0; k < adjlist[cur].size() && ret == 0; k++) {
  int i = adjlist[cur][k].num;
  int c = adjlist[cur][k].cap;
  if (c > used[cur][i] && !visited[i] && (ret = augment(i, min(flow,
   c-used[cur][i]), 0)) > 0) used[cur][i] += ret;
  else if (used[i][cur] > 0 && !visited[i] && (ret = augment(i, 
   min(flow, used[i][cur]), 0)) > 0)
   used[i][cur] -= ret;
 }
 return ret;
}

int SharksDinner::minSurvivors(vector <int> size, vector <int> speed,
      vector <int> intelligence) {
 int res = 0;
 vector<Shark> sharkData;
 for (int i = 0; i < size.size(); i++) {
  sharkData.push_back(Shark(size[i],speed[i],intelligence[i]));
 }
 adjlist = vector<vector<node> >(sharkData.size() * 3 + 2, vector<node>());
 // source
 for (int i = 0; i < sharkData.size(); i++) {
  adjlist[0].push_back(node(i+1,INF));
 }
 // sub-nodes
 for (int i = 0; i < sharkData.size(); i++) {
  adjlist[i+1].push_back(node(i+1+sharkData.size(),2));
 }
 // end-nodes
 for (int i = 0; i < sharkData.size(); i++) {
  for (int j = 0; j < sharkData.size(); j++) {
   if (i == j) continue;
   if (sharkData[i].size < sharkData[j].size || 
      sharkData[i].speed < sharkData[j].speed || 
      sharkData[i].intel < sharkData[j].intel) continue;
   if (sharkData[i].size == sharkData[j].size && 
      sharkData[i].speed == sharkData[j].speed && 
      sharkData[i].intel == sharkData[j].intel && i > j) continue;
   adjlist[i+1+sharkData.size()].push_back(node(j+1+(sharkData.size()*2),INF));
   adjlist[j+1+(sharkData.size()*2)].push_back(node(i+1+sharkData.size(),0));
  }
 }
 // sink
 for (int i = 0; i < sharkData.size(); i++) {
  adjlist[i+1+(sharkData.size()*2)].push_back(node(3*sharkData.size()+1,1));
 }
 sz = 3 * sharkData.size() + 2;
 // run flow
 int flow = maxflow();
 return res = sharkData.size() - flow;
}

TC SRM360 D1 500 (PrinceOfPersia)

The problem statement can be accessed via:
http://www.topcoder.com/stat?c=problem_statement&pm=7876&rd=10772

This problem lends itself to a graph problem as can be observed by the fact that, the input format is in a grid format and there are obstacles (forbidden cells) and its related to path accessibility between the prince and princess. One way to solve this problem is to use the Maximum Flow algorithm. Take note, although they are required to "meet" in any cell - it is sufficient enough for the prince (or princess) to go to princess (or prince) as it is not limited by time.

How does it relate to the Maximum Flow problem? We are basically asked to find the minimum cut and using the max-flow min-cut theorem we can calculate this by running a maximum flow algorithm through the graph. To see that what we really want is the minimum cut, let's revise the definition of a cut (in terms of graph theory of course): A cut is a partition of vertices of a graph into two disjoint subsets. So if we imagine the princess in one subset and the prince in another subset, the minimum cut between these two disjoint subsets is the number of accessible ways in which one can go to another.

So now we move on to applying the actual algorithm to the problem. If we let each empty cell be a node (or vertex) on the network and assign 1-capacity edges to adjacent empty nodes and apply flow from one "P" (source) to another "P" (sink) do we obtain the correct answer? Well, unfortunately no. The problem lies in which one "bottleneck" node (on the minimum cut set) can have an outgoing flow of more than 1, so it may over-calculate the answer. To see this, consider the example:

# . P
. . .
P . #

The central node in the snippet can have an incoming flow of 2 if both its adjacent nodes gives it a flow of 1. The problem with this is that the bottom-left P (assume it to be the sink) will have a maximum flow of 2 (and hence a minimum cut of 2). So our problem now is to ensure that each empty node will only give out a maximum of 1 out-going flow.

One wrong approach is to ensure that there is only an incoming flow of maximum 1 (i.e. 1 edge going into the node). The problem with this is that it totally breaks down the network layout. The solution is to further break each empty cell into 2 nodes: let's call them A and B. These nodes are used as "filters" or more precisely input and output filters. The edge from A->B will have a capacity of 1 which limits the maximum flow across the node to be 1 whereas the incoming and outgoing arcs of A and B respectively are infinite. This allows the network to retain the correct layout (and hence the correct answer) as well as fixing the original problem. We just need to ensure that all input arcs to the node use the sub-node A, and all output arcs from the node use the sub-node B.

Implementation-wise it's pretty standard if you follow the ideas discussed earlier. A sample implementation from the practice room is given below:

// the maximum number of vertices
#define NN 256

// adjacency matrix (fill this up)
int cap[NN][NN];

// flow network
int fnet[NN][NN];

// BFS
int q[NN], qf, qb, prev[NN];

int fordFulkerson( int n, int s, int t )
{
    memset( fnet, 0, sizeof( fnet ) );
    int flow = 0;
 // find an augmenting path of at least 1
    while( true )
    {
        memset( prev, -1, sizeof( prev ) );
        qf = qb = 0;
        prev[q[qb++] = s] = -2;
        while( qb > qf && prev[t] == -1 )
            for( int u = q[qf++], v = 0; v < n; v++ )
                if( prev[v] == -1 && fnet[u][v] - fnet[v][u] < cap[u][v] )
                    prev[q[qb++] = v] = u;
        if( prev[t] == -1 ) break;
        // get the bottleneck capacity
        int bot = 0x7FFFFFFF;
        for( int v = t, u = prev[v]; u >= 0; v = u, u = prev[v] )
            bot <?= cap[u][v] - fnet[u][v] + fnet[v][u];
        // update the flow network
        for( int v = t, u = prev[v]; u >= 0; v = u, u = prev[v] )
            fnet[u][v] += bot;
        flow += bot;
    }
    return flow;
}

int dx[] = {-1, 1, 0, 0};
int dy[] = {0, 0, -1, 1};

#define INF 999999

int PrinceOfPersia::minObstacles(vector <string> maze) {
 int res = 0;
 bool gotP = false;
 int source = 0, sink = 0;
 for (int i = 0; i < maze.size(); i++) {
  for (int j = 0; j < maze[i].size(); j++) {
      // assign the source and sink
   if (maze[i][j] == 'P' && !gotP) { 
     source = (i * maze[0].size() + j)*2; gotP = true; 
   }
   else if (maze[i][j] == 'P' && gotP) { 
     sink = (i * maze[0].size() + j)*2; 
   }
   // build connections
   if (maze[i][j] == '.') { 
     cap[(i * maze[0].size() + j)*2][(i * maze[0].size() + j)*2+1] = 1; 
   }
   else if (maze[i][j] == 'P') { 
     cap[(i * maze[0].size() + j)*2][(i * maze[0].size() + j)*2+1] = INF; 
   }
   // iterate through adjacent cells
   for (int k = 0; k < 4; k++) {
    int mi = i + dx[k];
    int mj = j + dy[k];
    if (mi < 0 || mj < 0 || mi >= maze.size() || mj >= maze[0].size()) 
      continue;
    // impossible case
    if (maze[mi][mj] == 'P' && maze[i][j] == 'P') return -1;
    if (maze[mi][mj] != '#') { 
     cap[(i * maze[0].size() + j)*2+1][(mi * maze[0].size() + mj)*2] = INF;
    }
   }
  }
 }
 int flow = fordFulkerson(maze[0].size()*maze.size()*2, source, sink);
 return res = flow;
}