Discrete Processes: Binomial Tree Model

After reading a bit from Baxter and Rennie's book about discrete processes, I decided to implement the basic binomial tree described in the book. It's a build onto the binomial branch model, in essence it's the same thing but done on a recursive scale. It is mostly used for evaluating the cost of an option via arbitrage. The main downfall to discrete process modelling computationally is that it becomes exceedingly difficult to scale it for problems of large instances. Unlike continuous processes which are evaluated through approximation, discrete processes models each tiny step and for an accurate representation would require a lot of memory and computational time.

The model suggests having a portfolio of stocks and bonds. This simple combination allows us to synthesize the derivative by replicating desired values. The key point is the probability of up and down moves from a given stock at a given time does not influence our ability to generate the desired values. For example, we can say for certainty that if the stock moved up then down and then up again we would reach our desired value f(n). We use an algebraic relation from the next stock up and down prices to determine our claim value. We evaluate the claim value from the back of the tree where the payoff/desired value is guaranteed and move our way back to our initial node. This claim value at the initial node is what we should price the option at.

Now if the market was to move up or down, we need to adjust our proportion of stocks and bonds to compensate. It's essentially like a board game, when the stock moves we make an appropriate adjustment to ensure that the next move (either up or down) is within our described tree regardless of the outcome. At the end of the time period, we are guaranteed to get the claim value for the given sequence of stock movements from t = 0 to t = n, where n is the number of time ticks (i.e. depth of the binomial tree). For more information refer to Baxter and Rennie's Financial Calculus book.

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

using namespace std;

// represent the stock as a value in time
class Stock {
public:
   double price;
   double prUp;
   double prDown;
   Stock() { }
   Stock(double p, double up, double down) : price(p), prUp(up), prDown(down) { }
};

vector<double> fTable;   // represents the claim tree
vector<Stock> graph;   // represents the recombinant tree
int sz;

#define EPS 1e-09

// claims are calculated as:
//      q = (stock_now - stock_down) / (stock_up - stock_down)
//      f = q * f_up + (1 - q) * f_down
// base conditions can be computed based on valuation of the option at t=0 versus the
// payoff given a sequence of situations (when the stock goes up/down for t=0 to t=n).
double computeClaims(int node) {
   int baseNode = (sz-1) >> 1;
   if (fTable[node] >= 0) {
      // memoisation
      return fTable[node];
   }
   if (node >= baseNode) {
      // base case
      return fTable[node] = max(0.0, graph[node].price - graph[1].price);
   }
   // compute
   double q = (graph[node].price - graph[node*2+1].price) / 
      (graph[node*2].price - graph[node*2+1].price);
   fTable[node] = q * computeClaims(node*2) + (1 - q) * computeClaims(node*2+1);
   return fTable[node];
}

// define an recombinant tree and reconstruct the option value
int main() 
{
   int N;   // the depth of the tree
   cin >> N;
   graph = vector<Stock>((1 << N) + 1);
   
   // read in the recombinant tree (stated in order of time)
   // file format:
   // stockprice prUp
   // stockprice2 prUp2
   // etc.

   int totalCnt = (1 << N);
   for (int i = 0; i < totalCnt-1; i++) {
      double stockPrice, prUp, prDown;
      cin >> stockPrice >> prUp;
      prDown = 1 - prUp;
      graph[i+1] = Stock(stockPrice, prUp, prDown);
   }

   // reconstruct the option price by backtracking
   fTable = vector<double>((1 << N) + 1, -1.0);
   sz = (1 << N) + 1;
   double optionPrice = computeClaims(1);
   cout << "Option price is " << optionPrice << "\n";
   
   // run simulations (optional)
   // file format:
   //   sequence of {up, down}'s equal to N-1.
   // note: stockHolding is calculated by:
   //   s = (f_up - f_down) / (s_up - s_down)
   string status;
   string lastStatus = "-";
   int cnt = 0;
   int curTime = 0;
   double stockHolding = 0.0;
   double bondHolding = 0.0;
   int curBranch = 1;
   while (cin >> status && cnt++ <= N-1) {
      // output the status
      cout << "Time: " << curTime << " Last Jump: " << lastStatus << 
         " Stock Price: " << graph[curBranch].price << 
         " Option Value: " << fTable[curBranch] << " Stock Holding: " 
         << stockHolding << " Bond Holding: " << bondHolding << "\n";
      curTime++;
      lastStatus = status;
      stockHolding = (fTable[curBranch*2] - fTable[curBranch*2+1]) / 
         (graph[curBranch*2].price - graph[curBranch*2+1].price);
      bondHolding = fTable[curBranch] - (stockHolding * graph[curBranch].price);
      if (status == "up") {
         curBranch = (curBranch * 2);
      } else {
         curBranch = (curBranch * 2) + 1;
      }
   }
   if (cnt > 0) {
      cout << "Time: " << curTime << " Last Jump: " << lastStatus << 
         " Stock Price: " << graph[curBranch].price << 
         " Option Value: " << fTable[curBranch] << " Stock Holding: " 
         << stockHolding << " Bond Holding: " << bondHolding << "\n";
   }
   return 0;
}

With an example file input based on the book examples (pipe it in):

4
100 0.75
120 0.75
80 0.25
140 0.75
100 0.25
100 0.75
60 0.25
160 0.75
120 0.25
120 0.75
80 0.25
120 0.75
80 0.25
80 0.75
40 0.25
up
up
down

TC SRM 454 D1-500/D2-1000 (NumbersAndMatches)

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10709

We are given a number which is composed of a sequence of digits. Each digit is represented by a combination of matches represented in a certain position to represent the digit. We are given the option of using up to K different match moves and are asked to calculate the number of different integers we can construct without discarding or adding any new matches to our given number.

First, we decompose the problem into several parts. The first involves somehow representing each digit in our program. A simple way to do this is to use a 2D integer array in which index (i,j) is equal to 1 if the digit i has a match in position j. An alternative way is to use a bitmask with the bit representation of the j-th bit set to 1 if there is a match in position j.

Next, we need to somehow compute the transitional costs between two different digits. For example, we need to quantify how many moves we would use up if we convert say a 6 to a 3. Since we are only allowed a limited number of moves its obvious we need a way to keep track of how many we have used up. One possible way is to determine how many places the two numbers differ. Then the number of moves we need to make is equivalent to exactly half of this.

There is another issue we need to deal with, that is the situation where the number of matches used to construct two digits differ. For example, the digit 6 uses 6 matches whereas digit 3 only uses 5 matches. We need to keep track of how much surplus (or deficit) matches we get from transitioning from one digit to another. This is important as we cannot add or remove matches from the final number - i.e. the number of matches used in constructing a unique integer must be the same as the original integer. Now we need to keep track of two things when we consider transforming one digit to another: the number of moves it takes and the difference in matches we have used so far.

We can define the cost function as follows:

To calculate the surplus/deficit it's quite simple, just calculate the number of matches used for each digit and subtract one from the other depending on which direction the transformation is taken. The above cost function ensures that we don't double count the number of moves required from moving a match from a surplus digit to a deficit digit. A way to visualise this is we "precharge" the surplus with an additional move but we are also granted a free move since this surplus is already counted once so the actual number of moves decreases if we need to import more matches from another transformation.

Lastly, we need to somehow compute all the ways we can construct unique integers given that we keep the number of matches used the same (i.e. at the end of our decision algorithm the number of surplus matches is exactly 0) as well as using less than or equal to K moves (which we subtract along the way). This suggests a recursive algorithm from potentially changing the first digit to something else, depending on the transformation we are left with either more/less/equal number of matches and the number of available moves to be less or equal to what we started with. We must also handle the case where we "borrow" matches from the latter part of the integer, so we must also handle cases where our match balance is negative. This can be done by computing an offset to keep the total number positive (to use for our caching purposes).

To efficiently compute the number of ways we need to memoise or apply bottom up DP to our algorithm. Note that the number of unique states is rather small so computational time isn't an issue. An implementation of the above ideas is shown below (using memoisation):

class NumbersAndMatches {
public:
   long long differentNumbers(long long, int);
};

int mat[10][7] = {
  {1,1,1,0,1,1,1}, // 0
  {0,0,1,0,0,1,1}, // 1
  {1,0,1,1,1,0,1}, // 2
  {1,0,1,1,0,1,1}, // 3
  {0,1,1,1,0,1,0}, // 4
  {1,1,0,1,0,1,1}, // 5
  {1,1,0,1,1,1,1}, // 6
  {1,0,1,0,0,1,0}, // 7
  {1,1,1,1,1,1,1}, // 8
  {1,1,1,1,0,1,1}  // 9
};

// defines the total number of matches used for digit i
int tot[10] = {6, 3, 5, 5, 4, 5, 6, 3, 7, 6};   

long long dp[19][155][155];
pair<int,int> costTab[10][10];
string num;

#define SURPLUS_MOD 72

long long func(int idx, int moves, int surplus) {
   // offset due to the possibility of negative surplus
   int surplusMod = surplus + SURPLUS_MOD;
   if (idx >= num.size()) {
      // requires the surplus amount of matches to be 0 otherwise its not valid
      if (moves >= 0 && surplus == 0) return 1;
      return 0;   
   }
   if (dp[idx][moves][surplusMod] != -1) return dp[idx][moves][surplusMod];
   long long res = 0;
   // try each digit
   int curDig = num[idx] - '0';
   for (int dig = 0; dig <= 9; dig++) {
      int c = costTab[curDig][dig].first;
      int d = costTab[curDig][dig].second;
      if (moves - c >= 0) {
         res += func(idx+1, moves - c, surplus + d);
      }
   }
   return dp[idx][moves][surplusMod] = res;
}

long long NumbersAndMatches::differentNumbers(long long N, int K) {
   long long res = 0;
   memset(dp,-1,sizeof(dp));
   // compute the cost function for every digit -> digit transition
   for (int i = 0; i <= 9; i++) {
      for (int j = 0; j <= 9; j++) {
         int c = 0;
         for (int k = 0; k < 7; k++) {
            if (mat[i][k] != mat[j][k]) c++;
         }
         int s = tot[i] - tot[j];
         costTab[i][j] = make_pair((c+s)/2, s);
      }
   }
   // decompose the number into a string
   long long vN = N;
   while (vN > 0) {
      num += ((vN % 10) + '0');
      vN /= 10;
   }
   reverse(num.begin(),num.end());
   // calculate the number of ways
   return res = func(0, K, 0);
}

Matrix Exponentiation

Matrix Exponentiation is a useful technique that can be applied to a wide variety of problems. Most commonly it is used for efficiently solving linear recurrence problems and as such can be used in any problem that can be represented as a linear recurrence. The main problem lies with the fact that it is inefficient when we naively evaluate them. Let's consider an easy example that you should know about: Fibonacci numbers.

The sequence is defined by:

It's easy enough to increasingly evaluate it 1 by 1 each time until we hit our targetted n-th Fibonacci number. However, let's say n was 1 billion - we need at least 1 billion computations of the sequence to derive the answer. If we represent the sequence in terms of a matrix we discover a clever optimisation which can be applied.

The correctness can be verified by a simple matrix multiplication. What's more important to note is that the next two sequence of Fibonacci numbers can be represented like the one above.

If we expand:

This is equivalent to:

We can continue expanding and eventually we see a remarkable observation. We can always reduce it to the power of our matrix multiplied by our initial conditions. This yields the recurrence below:

 

Since exponentiation can be done in logarithmic time we have essentially reduced our Fibonacci computation speed from linear to logarithmic. A huge difference in speed for large numbers! As an exercise, write a logarithmic Fibonacci term generator - as the numbers do get quite large, feel free to modulo it with an appropriate number (like 10^8).

We now expand our knowledge to a slightly more complicated linear recurrence. This is based on the problem "Number Sequences": http://acm.tju.edu.cn/toj/showp2169.html

To solve this recurrence we simply just do the same thing as we did with Fibonacci with a slight modification. Fibonacci is really a special case of the above recurrence where x = 1 and y = 1 and a0 = a1 = 1. Note the prime observation below:

If we matrix multiply - we yield the correct recurrence for both Fn and Fn-1. Therefore we simply need to change our original Fibonacci matrix of [ 1 1, 1 0 ] to [ x y, 1 0] and the initial conditions from being always 1 and 1 (F1 and F0 respectively) to [a1 a0]. We then simply use matrix exponentiation to calculate the correct term, as always we apply modulo arithmetic to keep the number representable with integers.

The implementation below demonstrates this idea:

class Matrix {
public:
   long long a;
   long long b;
   long long c;
   long long d;
   Matrix() { }
   Matrix(long long _a, long long _b, long long _c, long long _d) : 
      a(_a), b(_b), c(_c), d(_d) { }
};

Matrix multiply(const Matrix& lhs, const Matrix& rhs) {
   long long a = lhs.a * rhs.a + lhs.b * rhs.c;
   long long b = lhs.a * rhs.b + lhs.b * rhs.d;
   long long c = lhs.c * rhs.a + lhs.d * rhs.c;
   long long d = lhs.c * rhs.b + lhs.d * rhs.d;
   return Matrix(a % 100, b % 100, c % 100, d % 100);
}

ostream& operator<<(ostream& out, const Matrix& M) {
   out << "a: " << M.a << " b: " << M.b << " c: " << M.c << " d: " << M.d;
   return out;
}

Matrix power(Matrix M, int n) {
   if (n == 1) return M;
   Matrix tmp = power(M, n/2);
   if (n % 2 != 0) {
      Matrix n = multiply(tmp,tmp);
      Matrix m = multiply(n, M);
      return m;
   }
   Matrix v = multiply(tmp,tmp);
   return v;
}

void output(int n) {
   if (n < 10) cout << "0";
   cout << n << "\n";
}

int main() {
   long long x, y, a0, a1, n;
   while (cin >> x >> y) {
      if (x == 0 && y == 0) break;
      cin >> a0 >> a1 >> n;
      if (n == 0) { output(a0%100); continue; }
      else if (n == 1) { output(a1%100); continue; }
      Matrix m = power(Matrix(x,y,1,0), n-1);
      output((m.a * a1 + m.b * a0) % 100);
   }
   return 0;
}

Now for a more advanced example, we use matrix exponentiation to determine the number of cycles with a length smaller than k in a given directed graph.

Problem: TourCounting
Source: TopCoder SRM 306 D1-1050
URL: http://www.topcoder.com/stat?c=problem_statement&pm=6386

An elegant way to calculate the number of paths from a source vertex u to a destination vertex v that takes exactly k steps is to use matrix multiplication. First let's define A to be the adjacency matrix with each entry representing the number of edges from u to v. The (u,v) entry of A^t is precisely the answer we are looking for. The way this works is similar to the way that Floyd-Warshall algorithm works. The matrix multiplication considers the connectivity between an intermediate node k and attempts to link u to k and k to v. Due to the multiplicative nature of matrix multiplication as opposed to an additive (for shortest paths) it determines the number of paths based on the multiplication principle from discrete mathematics.

To illustrate this principle, consider a snapshot of the matrix A^t and it's adjacency matrix which is represented as A^1:

If we were to calculate the number of paths that start from vertex 0 to itself (i.e. a cycle) that takes exactly t+1 steps. Note that we can access vertex 0 from either vertex 1 or vertex 2. Since at state t, we can reach vertex 1 from vertex 0 with 2 paths and we can reach vertex 2 from vertex 0 with 3 paths. It stands that we can now reach vertex 0 to itself from 2+3=5 paths. Note that this is a direct consequence of the matrix multiplication process - validate this yourself for the first row and column. If the adjacency matrix was changed such that there are two edges from vertex 1 to vertex 0 like the following:

Then we can reach vertex 0 in a total of 2*2 + 1*3 = 7 ways since we have an option of choosing one of the two edges from vertex 1 to vertex 0. Again, this is consequence of the matrix multiplication algorithm. The argument holds for every other cell in the matrix. So what's the advantage of this method? The prime factor is efficiency - we can calculate matrix powers in logarithmic number of matrix multiplications. We accomplish this through a similar algorithm to fast powering.

So now we know how to calculate the number of paths from a source vertex u to a destination vertex v that takes exactly k steps. The actual problem requires us to compute the number of cycles that are present with less than k steps. So for a given instance of the matrix A^t we want to sum the diagonal to yield the number of paths with exactly t steps. To calculate it for less than k steps we would normally have to compute each power of t and sum all the diagonals up, however this is not fast enough as the number of steps can be as much as 1 million. We can reduce this using a similar method to our matrix powering algorithm. First define a function f:

f(n) = represents the number of ways we can get from all vertices to other vertices using less than n steps (matrix form)

The main optimisation is reducing it from linear (summing all matrices less than n) to logarithmic. The huge hint in logarithmic is dividing the input space by a (usually) constant factor. Note that if we divide the space into two f(n/2) parts we are close to our answer but not quite there. The problem is that upper n/2 part is not symmetrical to the lower n/2 part. But wait! Given the simple mathematical property of powers a^(b+n) = a^b * a^n we can simply just multiply one of the f(n/2) parts by A^(n/2) to yield the correct upper half! With this our recurrence becomes:

Like our fast powering algorithm we need to distinguish between odd and even states. The easy way to make an odd number even is to subtract it by 1. We then just calculate one instance of the odd power and make all successively calls even. This maintains the correctness of the recurrence whilst also maintaining the time complexity.

Then it's simply a matter of implementation which is the easy part. Just note we need to use sensible modulo arithmetic to ensure that we don't overflow any of our computations.

class TourCounting {
public:
   int countTours(vector <string>, int, int);
};

class Matrix {
public:
   vector<vector<long long> > data;
   Matrix() { }
   Matrix(int n, int m) {
      data = vector<vector<long long> >(n, vector<long long>(m, 0));
   }
};

long long MOD;

Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
   Matrix res(lhs.data.size(), rhs.data[0].size());
   for (int i = 0; i < lhs.data.size(); i++) {
      for (int j = 0; j < lhs.data[i].size(); j++) {
         for (int k = 0; k < rhs.data[0].size(); k++) {
            res.data[i][k] = (res.data[i][k] + lhs.data[i][j] * rhs.data[j][k])%MOD;
         }
      }
   }
   return res;
}

Matrix operator+(const Matrix& lhs, const Matrix& rhs) {
   Matrix res(lhs.data.size(),lhs.data[0].size());
   for (int i = 0; i < lhs.data.size(); i++) {
      for (int j = 0; j < lhs.data[0].size(); j++) {
         res.data[i][j] = (lhs.data[i][j] + rhs.data[i][j])%MOD;
      }
   }
   return res;
}

Matrix power(const Matrix& lhs, int P) {
   if (P == 1) return lhs;
   Matrix tmp = power(lhs, P/2);
   if (P % 2 != 0) {
      Matrix n = tmp * tmp;
      Matrix m = n * lhs;
      return m;
   }
   Matrix v = tmp * tmp;
   return v;
}

long long compute(const Matrix& m) {
   long long res = 0;
   for (int i = 0; i < m.data.size(); i++) res = (res + m.data[i][i]) % MOD;
   return res % MOD;
}

Matrix dp[1000001];
int visited[1000001];

Matrix f(const Matrix& ref, int n) {
   if (visited[n]) return dp[n];
   visited[n] = 1;
   if (n == 1) return dp[n] = ref;
   if (n % 2 != 0) {
      // odd
      Matrix v = f(ref, n-1) + power(ref, n);
      return dp[n] = v;
   }
   // even 
   Matrix vE = f(ref, n/2) * power(ref, n/2);
   Matrix vEE = vE + f(ref, n/2);
   return dp[n] = vEE;
}

int TourCounting::countTours(vector <string> g, int k, int m) {
   Matrix mat((int)g.size(),(int)g.size());
   for (int i = 0; i < g.size(); i++)
      for (int j = 0; j < g[i].size(); j++)
         mat.data[i][j] = g[i][j] == 'Y' ? 1 : 0;
   MOD = m;
   long long res = 0;
   Matrix rr = f(mat, k-1);
   res = compute(rr);
   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);
}

Google Code Jam Round 2

Links:
Google Code Jam Statistics: http://www.go-hero.net/jam/
Google Code Jam Scoreboard: http://code.google.com/codejam/contest/scoreboard?c=204113

Overview
Round 2 proved to be a bit trickier than Round 1, especially in terms of advancement as only the top 500 go through. The time setting didn't help either for many contestants. Unfortunately I didn't make it through but congratulations to everyone who did! It has been a major blast of fun :)

So what happened? There are several key points I've learnt from participating this round.

1. Don't focus too much on constraints!

For problem A, my first instinct was a greedy based algorithm which was rather close to the correct algorithm. However, I quickly convinced myself (wrongly) that it would fail for cases in which you would possibly make a move which turns the remaining subproblem infeasible. This does not happen though as the strictly increasing ordering of the final output prevents this from happening. This was further misled by looking at the large constraint which had a maximum of 40. Surely a greedy algorithm can't work right? Wrong! :)

Although it is always good to look at constraints (and perhaps I've relied too heavily on it) to dictate what algorithm the problem setter is expecting don't let it "blind" (i.e. let it confirm your incorrect doubts of an algorithm without fully going through it) you though as it did to me.

2. Don't get too focused on a particular methodology and learn to back out of an idea if you can't get anywhere with it!

After solving A-small and D-small, the only way I would advance is to solve both C-small and C-large (or B but fewer people tried to solve it so I also tried C). Instead of wasting time on finding the correct solution to A-large, I proceeded with problem C.

First thought was bipartite matching (and that hunch was correct) as it felt natural (because I've actually solved a very similar problem). I couldn't reduce the graph properly with the correct edge labelling criteria so I completely abandoned this idea (bad move as I did in fact at one point have the correct algorithm but made a false assumption which would of made it incorrect). I eventually deduced it to finding the minimum number of complete induced subgraphs (clique problem). However this is NP-complete, so I really should of backed out of this reduction but I was persistent as this was the only "correct" (although it would time out) solution I had. Sometimes it pays off to let go of an idea if you can't get anywhere with it which is easier said than done!

Despite being slightly disappointed I had a fun time. Unfortunately no Australians (unofficially at this point) have made it through to Round 3. Good luck for all those participating in Round 3 though :)

Official Contest Analysis is available: http://code.google.com/codejam/contest/dashboard?c=204113#s=a

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 SRM424 D1 600 (StrengthOrIntellect)

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10174

This problem asks us to determine the maximum number of missions we can complete given two character attributes: strength and intellect. We meet the pre-requirements to complete a mission if either one of the strength or intellect attributes are at least as high/equal to the mission's strength and intellect requirements respectively.

A greedy approach does not work here. You can informally argue by noting that when we are given two choices to increase strength or intellect, there is no clear decision on which option to choose. The local optimality does not hold for the global optimality, as if we choose one or the other we might miss out on earning a lot of points which will allow us to complete more missions. Another greedy heuristic might be to try both maximising strength and maximising intellect and choose the maximum number of missions from the two possible options. This also does not work because the optimal solution might not be on a extreme point of strength or intellect. Since the decisions we make need to be flexible enough to change between increasing strength and/or then increasing intellect or vice-versa, it's rather clear that no amount of heuristics will pass all the (corner) test cases.

The next obvious approach to consider is brute force. We can keep track of which missions we have accomplished and how much strength and intellect we have achieved so far, then for all possible mission transitions recurse along the mission and adjust the strength and intellect. In effect, this almost works if we memoise the solution. The problem is keeping track of the missions we have completed - the naive approach is using a bitmask which will definitely exceed the memory constraint (2^50 at worst case).

The trickest part of the problem is state compressing our DP algorithm. If we can somehow reduce the memory footprint whilst also keeping the correctness intact, then a viable solution is generated. Let's take out our bitmask and try to find a way to implicitly represent that information. A common reduction in this situation is to somehow find the relationship between our current DP state parameters (in this case, strength and intellect) and the DP state parameter(s) we want to get rid of (in this case, how we can implicitly represent the missions completed so far). In fact we can see that they have a direct relationship! If we are given a particular strength and intellect we can already uniquely determine which missions can be completed, as we take all viable missions as all the points are strictly positive (always something to gain, nothing to lose).

Now our DP state becomes simply:

D(S, I) = the maximum number of missions if we can get to strength S and intellect I.
D(S, I) = 0 if the state is unreachable (not enough points).

We can proceed to pre-compute the maximum number of missions for each state as there are only 1000 * 1000 = 1 million such states. We can also pre-compute whether or not a state is possible or not. A state is impossible if we don't have enough points to cover the increase from strength 1 to S, and intellect 1 to I. Of course, this doesn't fully dictate the feasibility of a state but it's a good start. To compute the number of excess (unused) points we simply sum up all the points gained from completing the missions and subtracting S-1 and I-1 from it (as we expended S-1 strength points and I-1 intellect points to get to the current state). This calculation is done as follows:

pair<int,int> preComp[1011][1011];
...
   for (int i = 0; i <= 1000; i++) {
      for (int j = 0; j <= 1000; j++) {
         int pointsGained = 0, missions = 0;
         for (int k = 0; k < strength.size(); k++) {
            if (strength[k] <= i || intellect[k] <= j) {
               pointsGained += points[k];
               missions++;
            }
         }
         int extraPoints = pointsGained - i - j + 2;
         // extraPoints < 0 --> inaccessible state set to 0 missions
         // the pair is (excess points, number of missions completed)
         preComp[i][j] = make_pair(extraPoints, extraPoints < 0 ? 0 : missions);
      }
   }
...

Next, we need to determine feasibility from our original DP state (1,1) - which is the strength and intellect we start with. We recursively call the next DP state if we have unused points by either deciding to upgrade our strength by 1 or our intellect by 1. Note the number of unique states is only 1 million so this should run in time. Combining this with a caching mechanism (i.e. memoisation) we come up with our final solution:

class StrengthOrIntellect {
public:
   int numberOfMissions(vector <int>, vector <int>, vector <int>);
};

pair<int,int> preComp[1011][1011];
int memo[1011][1011];

int func(int strength, int intellect) {
   // seen it before (caching mechanism)
   if (memo[strength][intellect] != -1) return memo[strength][intellect];
   memo[strength][intellect] = preComp[strength][intellect].second;
   int extra = preComp[strength][intellect].first;
   // if there are unused points, recurse on them otherwise the maximum
   // is simply just the number of missions we have obtained so far
   return memo[strength][intellect] >?= extra <= 0 ? memo[strength]
      [intellect] : max(func(strength+1, intellect), func(strength, intellect+1));
}

int StrengthOrIntellect::numberOfMissions(vector <int> strength, vector <int> 
      intellect, vector <int> points) {
   int res = 0;
   for (int i = 0; i <= 1000; i++) {
      for (int j = 0; j <= 1000; j++) {
         int pointsGained = 0, missions = 0;
         for (int k = 0; k < strength.size(); k++) {
            if (strength[k] <= i || intellect[k] <= j) {
               pointsGained += points[k];
               missions++;
            }
         }
         int extraPoints = pointsGained - i - j + 2;
         // extraPoints < 0 --> inaccessible state set to 0 missions
         preComp[i][j] = make_pair(extraPoints, extraPoints < 0 ? 0 : missions);
      }
   }
   memset(memo,-1,sizeof(memo));
   return res = func(1, 1);
}

Maximum/Minimum Weighted Bipartite Matching using Cycle Cancelling

Problem:

This is an extension to our maximum cardinality bipartite matching problem we introduced earlier. Imagine the same situation, we are given a bipartite graph G = (V,E) in which the vertices can be separated into two disjoint sets such that there are no edges between vertices belonging in the same set. Instead of weightless edges like our previous problem - we are now faced with weighted edges. The problem now is to match the vertices such that the total weight of the matched edges is as small (or as large) as possible whilst still retaining the maximum number of matchings. The figure below depicts this situation:

Applications:

Like the weightless version, we can compute the most optimal pairings between two disjoint sets. For example, let X represent a set of workers and Y to be a set of jobs. There is an edge (u,v) in E of worker u in X to job v in Y of weight W, which gives the productivity of the worker u given job v. If we were tasked to assign the workers to the set of jobs such that each worker is given exactly 1 job and we were also told to maximise the total productivity of the assignments then the answer is simply the maximum weighted matching in the bipartite graph. We can also define the edge weights to mean the converse, for example, the weight for edge (u,v) could mean how difficult worker u would find job v – in which our task could be to minimise the total difficulty whilst ensuring all workers are assigning to exactly 1 job. This does not complicate things as we simply need to reverse our operations to derive the opposite result.

Algorithm:

We need a new approach to tackle this problem. Augmenting the paths isn't enough as it doesn't take into account the weightings on the edge. One idea is to arbitrarily assign the maximum matching without consideration of the weights and then try to successively improve (i.e. maximise) the cost by finding a negative cycle an augmenting the matchings along this cycle. To represent this we calculate a modified graph G’, with edge weights representing the cost difference we will gain by matching a vertex x with another vertex y. In other words, for any u and v, we put a weight W on the edge from u to v, where W is defined as the cost we will gain if we gave the current partner of v to u, i.e. W = cost(u, matching[v]) – cost(u, matching[u]). Now we simply look for a negative cycle in G’ and augment the partners along the cycle to produce a better matching. How does this work? If we can find a negative cycle from a vertex u to itself, it means we can alternate partners to reduce the total cost whilst also keeping the same number of matched partners.

Consider the following bipartite graph:

And the matchings below:

 

The matching on the left has a weight of 15; the optimal minimum matching has a weight of 8 (shown on the right). Algorithmically, let’s assign the vertices on the top-left and bottom-left (A and B respectively) and the vertices on the top-right and bottom-right (C and D respectively). We compute the modified G’ edge weights of (A,B) and (B,A) – the cost of switching the partners around (take note their matched counterparts – i.e. the right side, is calculated implicitly). Proceeding to do this we yield:

W(A,B) = cost(A, matching[B]) – cost(A, matching[A]) = 6 – 5 = 1
W(B,A) = cost(B, matching[A]) – cost(B, matching[B]) = 2 – 10 = -8

Therefore a possible cycle could be from A -> matching[B] -> B -> matching[A] -> A. This yields a cost of W(A,B) + W(B,A) = 1 + (-8) = -7. Since this is a negative cycle we can save a total cost of up to 7 (you can verify this by hand) if we augment along this cycle. From this, A becomes matched to matching[B] (i.e. D) and B becomes matched to matching[A] (i.e. C). The matching on the right hand side diagram highlights this minimal matching. To actually find a negative cycle, an easy way is to use Floyd-Warshall’s algorithm which is illustrated in our implementation below. Additional book-keeping is required to keep track of the parents so we can back-track and augment the negative cycle path.

To summarise, our algorithm is as follows:

- Start with an initial perfect matching M
- While there is a negative cycle C on G’, augment along the negative cycle C to produce a better matching M’
- Return the matching M

Other Notes:

This problem can also be solved using the Hungarian algorithm (a famous combinatorial optimisation algorithm) or using minimum-cost flows. In fact, our cycle cancelling implementation is a specific subset of the minimum cost flow algorithm. The difference between our cycle cancelling algorithm and the one used for the general minimum cost flow problem is that we need to run a maximum flow algorithm to derive an initial network (here we simply assign the vertices to each other as we are guaranteed a perfect matching it’s a complete bipartite graph).

There are many variations on maximum (or minimum) weighted bipartite matching. The first variation is called the assignment problem where we are given an equal number of vertices on each side and the graph itself is complete (i.e. all the vertices link to each other between the two disjoint sets). The second variation we are given an uneven number of vertices on each side but the graph itself is still complete, here we proceed to reduce the problem into the first variation by adding dummy vertices to make the sides even. Any edges which go to these vertices have a weight of 0 and hence do not have an effect on the final answer. Another variation is where there are an uneven number of vertices and the graph isn’t complete – hence a perfect matching may not be possible, this is where the general minimum cost flow algorithm comes in which we will go over later.

Applying the Algorithm:

We briefly demonstrate how the algorithm is used to solve two algorithmic problems.


ACM ICPC Pacific Northwest Regionals 2004 (GoingHome)
URL: http://acm.tju.edu.cn/toj/showp1636.html

The problem can be easily reduced to the assignment problem by constructing a bipartite graph of houses and people. The edge weights are simply the Manhattan distance between house i and person j. We then just run our algorithm over the graph and we obtain our answer!

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

using namespace std;

int cost[111][111];      // cost of assigning house i to person j
int curMatchings[111];   // house i matched to person curMatchings[i]
int delta[111][111];   // the delta/cost graph
int parent[111][111];   // used for backtracking the cycle
int cycle[111];         // path of the cycle
int totalLen;         // length of the cycle

bool augmentingCycle(int N) {
   memset(delta,0,sizeof(delta));
   memset(parent,0,sizeof(parent));
   memset(cycle,-1,sizeof(cycle));
   // derive the delta graph
   for (int i = 0; i < N; i++) {
      for (int j = 0; j < N; j++) {
         delta[i][j] = cost[i][curMatchings[j]] - cost[i][curMatchings[i]];
         parent[i][j] = j;
      }
   }
   for (int k = 0; k < N; k++) {
      for (int i = 0; i < N; i++) {
         for (int j = 0; j < N; j++) {
            if (delta[i][k] + delta[k][j] < delta[i][j]) {
               // minimum cost matching
               delta[i][j] = delta[i][k] + delta[k][j];
               // book-keep the optimal path so far for i->j
               parent[i][j] = parent[i][k];
               // detect for cycle (if vertex i and j are the same vertex)
               if (i == j) {
                  totalLen = 0;
                  // backtrack and construct the cycle path
                  do {
                     cycle[totalLen++] = i;
                     i = parent[i][j];
                  } while (i != j);
                  return true;
               }
            }
         }
      }
   }
   // did not find an augmenting path
   return false;
}

int cycleCancel(int N) {
   int res = 0;
   // pseudo max-flow for assignment problem
   for (int i = 0; i < N; i++) curMatchings[i] = i;
   while (augmentingCycle(N)) {
      // augment the negative cycle
      int r = curMatchings[cycle[0]];
      for (int i = 0; i < totalLen - 1; i++) {
         curMatchings[cycle[i]] = curMatchings[cycle[i+1]];
      }
      curMatchings[cycle[totalLen-1]] = r;
   }
   // compute the final cost
   for (int i = 0; i < N; i++) {
      res += cost[i][curMatchings[i]];
   }
   return res;
}

int main() {
   int N, M;
   while (cin >> N >> M) {
      if (N == 0 && M == 0) break;
      char ch;
      vector<pair<int,int> > houses;
      vector<pair<int,int> > people;
      for (int i = 0; i < N; i++) {
         for (int j = 0; j < M; j++) {
            cin >> ch;
            if (ch == 'H') houses.push_back(make_pair(i,j));
            else if (ch == 'm') people.push_back(make_pair(i,j));
         }
      }
      // compute the weights associated with matching house i to person j
      for (int i = 0; i < houses.size(); i++) {
         for (int j = 0; j < people.size(); j++) {
            cost[i][j] = abs(houses[i].first - people[j].first) + 
               abs(houses[i].second - people[j].second);
         }
      }
      cout << cycleCancel(houses.size()) << "\n";
   }
   return 0;
}

TC SRM387 D2 1000 (MarblesRegroupingHard)
URL: http://www.topcoder.com/stat?c=problem_statement&pm=8538

This is a re-visit of a previous problem we solved using DP. This is a slightly less obvious reduction than the previous problem; we need to compute the total number of colours which exists in the whole dataset. Then we can draw an edge from box i to colour j with a total weight of all the other j-coloured marbles from other boxes minus the ones in box i with colour j (as they don’t need to be moved). Take note that there may be a different number of colours to the number of boxes but we are always guaranteed that the number of colours is less than or equal to the number of boxes, hence a matching is always possible based on the pigeonhole principle. We can add dummy colours to the graph to reduce it to the assignment problem – which is what is done below:

class MarblesRegroupingHard {
public:
   int minMoves(vector <string>);
};

int totColours[15];
int B[51][15];

int cost[55][55];      // cost of assigning house i to person j
int curMatchings[55];   // house i matched to person curMatchings[i]
int delta[55][55];   // the delta/cost graph
int parent[55][55];   // used for backtracking the cycle
int cycle[55];         // path of the cycle
int totalLen;         // length of the cycle

bool augmentingCycle(int N) {
   memset(delta,0,sizeof(delta));
   memset(parent,0,sizeof(parent));
   memset(cycle,-1,sizeof(cycle));
   // derive the delta graph
   for (int i = 0; i < N; i++) {
      for (int j = 0; j < N; j++) {
         delta[i][j] = cost[i][curMatchings[j]] - cost[i][curMatchings[i]];
         parent[i][j] = j;
      }
   }
   for (int k = 0; k < N; k++) {
      for (int i = 0; i < N; i++) {
         for (int j = 0; j < N; j++) {
            if (delta[i][k] + delta[k][j] < delta[i][j]) {
               // minimum cost matching
               delta[i][j] = delta[i][k] + delta[k][j];
               // book-keep the optimal path so far for i->j
               parent[i][j] = parent[i][k];
               // detect for cycle (if vertex i and j are the same vertex)
               if (i == j) {
                  totalLen = 0;
                  // backtrack and construct the cycle path
                  do {
                     cycle[totalLen++] = i;
                     i = parent[i][j];
                  } while (i != j);
                  return true;
               }
            }
         }
      }
   }
   // did not find an augmenting path
   return false;
}

int cycleCancel(int N) {
   int res = 0;
   // pseudo max-flow for assignment problem
   for (int i = 0; i < N; i++) curMatchings[i] = i;
   while (augmentingCycle(N)) {
      // augment the negative cycle
      int r = curMatchings[cycle[0]];
      for (int i = 0; i < totalLen - 1; i++) {
         curMatchings[cycle[i]] = curMatchings[cycle[i+1]];
      }
      curMatchings[cycle[totalLen-1]] = r;
   }
   // compute the final cost
   for (int i = 0; i < N; i++) {
      res += cost[i][curMatchings[i]];
   }
   return res;
}


int MarblesRegroupingHard::minMoves(vector <string> boxes) {
   int res = 0;
   memset(totColours,0,sizeof(totColours));
   memset(B,0,sizeof(B));
   int numColours = 0;
   for (int i = 0; i < boxes.size(); i++) {
      istringstream iss(boxes[i]);
      int colours;
      int k = 0;
      while (iss >> colours) {
         totColours[k] += colours;
         B[i][k++] = colours;
      }
      numColours = k;
   }
   // make the weighted bipartite graph
   // use the cost matrix to fill in
   for (int i = 0; i < boxes.size(); i++) {
      for (int j = 0; j < boxes.size(); j++) {
         if (j >= numColours) cost[i][j] = 0;
         else cost[i][j] = totColours[j] - B[i][j];
      }   
   }
   return res = cycleCancel(boxes.size());
}

Google Code Jam Round 1C

Overview:
This was the last opportunity for qualifiers to advance into Round 2. Congratulations to everyone who has made it to Round 2! Bad luck to those that didn't, there's always next year!

This round was greeted with a set of interesting problems. Problem A was the easy problem in the set, and most competitors were able to solve the small input with some incorrect submissions for the large case. Problem B proved rather interesting and had multiple correct approaches. One of which required a little mathematical manipulation to derive a short and safe implementation, the other had the danger of precision issues but was based on a simpler concept without the use of much mathematics. Problem C was a fairly subtle DP problem in disguise but this did not deter the top competitors to solve all three problems in under 40 minutes.

Links:
Google Code Jam Statistics: http://www.go-hero.net/jam/
Google Code Jam Scoreboard: http://code.google.com/codejam/contest/scoreboard?c=189252
Google Code Jam Unofficial Scoreboard/Online Source Viewer (Thanks to Zibada on TC): http://zibada.ru/gcj/2009/round1a.html, http://zibada.ru/gcj/2009/round1b.html, http://zibada.ru/gcj/2009/round1c.html

Problem A: All Your Base
Category: Greedy
URL: http://code.google.com/codejam/contest/dashboard?c=189252#s=p0

The problem asks us to determine the minimum possible number we can derive from an unknown number representation system. To do this we must assign values to each of the characters present in the symbols. For example let's say we are given "cats" as an input we need to determine what the 'c' means numerically, what the 'a' means numerically etc. Since we want the smallest possible number, it suggests a greedy strategy of assigning what each letter means. In our case, we first see 'c' and therefore wish for this to be as small as possible. However, there is an additional constraint that no number starts with a 0 and since 'c' is the first symbol it cannot be 0. Therefore, the next obvious choice is to give it a value of 1. Next, we process 'a' and assign it a value of 0 since it hasn't been used yet and it's not the first symbol. Following this, 't' gets assigned 2 and 's' gets assigned 3.

Now the next task is to determine the base system the aliens work on, it can be easily shown we want the base as small as possible - yet large enough to allow us to use all the symbols. Therefore, the base system is simply the number of unique characters/letters in the alien number. We are then tasked with converting this to "human" form (i.e. decimal representation) which can be done iteratively like any base conversion algorithm. Implementation wise, we can keep track of our character assignments in a map data structure. When we encounter a letter, we first check whether or not we have seen the letter - if we haven't then we greedily assign it the next available (and smallest) number. If we have, then we simply skip over the character.

#include <iostream>
#include <vector>
#include <string>
#include <sstream>
#include <algorithm>
#include <stack>
#include <queue>
#include <map>

using namespace std;

map<char,int> memo;

int main() {
   int nT; cin >> nT;
   string str;
   int k = 1;
   while (nT-- && cin >> str) {
      memo.clear();
      int next = 0;
      memo[str[0]] = 1;
      for (int i = 1; i < str.size(); i++) {
         if (memo.count(str[i]) == 0) {
            memo[str[i]] = next;
            next++;
            if (next == 1) next++;
         }
      }
      long long val = 0;
      int base = next == 0 ? 2 : next;
      for (int i = 0; i < str.size(); i++) {
         val = val * base + memo[str[i]];
      }
      cout << "Case #" << k << ": " << val << "\n";
      k++;
   }
   return 0;
}

Problem B: Centre of Mass
Category: Maths, Ternary Search
URL: http://code.google.com/codejam/contest/dashboard?c=189252#s=p1

This is a tricky problem to tackle due to precision issues, so it's safer to go for a closed-form solution. Basically, we are given N fireflies each starting off in different positions and each going their own directions with a constant velocity. We need to calculate the distance and time in which the centre of mass of the N fireflies are closest to the origin (0,0,0). We need to make one key observation to heavily simplify this problem. If we look closely, we can collapse all of the N fireflies into one single entity (this is heavily suggested by the notes section in the actual problem statement). If we only have a single point (i.e. the centre of mass) we simply need to find the point in time/distance in which it comes closest to the origin. With some visualisation we can see that the centre of mass movement is divided into two cases:

- There is no movement hence the earliest time and closest distance it is to the origin is at time = 0.
- There is movement (i.e. velocity) then the earliest time and closest distance is unique (as there would be a point in which it gets closer and closer and after the minimum distance it gets larger and larger).

With these two conditions, we have multiple approaches to solve the problem. The first approach is to use a ternary search over the set of possible times, as there is only one minima for the movement case - we are guaranteed that the algorithm will converge to the spot eventually. Ternary search works on functions which have only one extrema that is either strictly increasing then strictly decreasing or the other way around. For more information on this search technique look at: http://en.wikipedia.org/wiki/Ternary_search. The second approach is to differentiate the distance equation - as we know there is only 1 solution we can confidently differentiate and solve for d/dt = 0 to obtain the unique solution. We will now proceed to demonstrate this approach.

We first need to collapse the N fireflies into the centre of mass vector with respect to time. This can be done by a straight summation over the coefficients of the starting positions and velocities. Let S = (Sx, Sy, Sz) be the summation of the starting positions of the N fireflies (this is a straight application of the formula given in the notes section). Let D = (Dx, Dy, Dz) be the summation of the velocities of the N fireflies (similar argument applies to the velocity). Then we yield a distance (with respect to time T) equation of:

D = sqrt((Sx + Dx * T)^2 + (Sy + Dy * T)^2 + (Sz + Dz * T)^2)

We can ignore the sqrt for differentiation purposes as it will get cancelled out later when it's multiplied by 0 (another reason is the function is strictly increasing, and as a result has no effect on the location of the extrema):

D = (Sx + Dx * T)^2 + (Sy + Dy * T)^2 + (Sz + Dz * T)^2

Now let's differentiate (using the Chain Rule):

dD/dT = 2*(Sx + Dx * T)*(Dx) + 2*(Sy + Dy * T)*(Dy) + 2*(Sz + Dz * T)*(Dz)

Let dD/dT = 0,

0 = 2*(Sx + Dx * T)*(Dx) + 2*(Sy + Dy * T)*(Dy) + 2*(Sz + Dz * T)*(Dz)
0 = 2*((Sx + Dx * T)*Dx + (Sy + Dy * T)*Dy + (Sz + Dz * T)*Dz)
0 = (Sx + Dx * T)*Dx + (Sy + Dy * T)*Dy + (Sz + Dz * T)*Dz

Expanding:

0 = Sx*Dx + Dx*Dx*T + Sy*Dy + Dy*Dy*T + Sz*Dz + Dz*Dz*T

Pushing all the T terms into the LHS:

-Dx*Dx*T - Dy*Dy*T - Dz*Dz*T = Sx*Dx + Sy*Dy + Sz*Dz
T(-Dx*Dx - Dy*Dy - Dz*Dz) = Sx*Dx + Sy*Dy + Sz*Dz

Therefore T is equal to:

T = Sx*Dx + Sy*Dy + Sz*Dz / (-Dx*Dx - Dy*Dy - Dz*Dz)

Now we need to be careful with dividing by 0, if (-Dx*Dx - Dy*Dy - Dz*Dz) is 0 this means that there is no movement (i.e. our first condition). We can safely make a separate case for this as the closest time is at t=0. Using these concepts we just need to do it which is pretty simple:

#include <iostream>
#include <string>
#include <algorithm>
#include <cmath>

using namespace std;

int main() {
   int T; cin >> T;
   for (int tt = 0; tt < T; tt++) {
      int N; cin >> N;
      double dX = 0.0, dY = 0.0, dZ = 0.0;
      double sX = 0.0, sY = 0.0, sZ = 0.0;
      // summation of the different components
      for (int i = 0; i < N; i++) {
         double x, y, z, tx, ty, tz;
         cin >> x >> y >> z >> tx >> ty >> tz;
         dX += tx; dY += ty; dZ += tz;
         sX += x; sY += y; sZ += z;
      }
      dX /= N; dY /= N; dZ /= N;
      sX /= N; sY /= N; sZ /= N;
      int before = 0;
      double t;
      // check for division by zero case
      if (dX * dX + dY * dY + dZ * dZ < 1e-09) {
         before = 1;
      } else {
         t = -(sX*dX + sY*dY + sZ*dZ) / (dX*dX + dY*dY + dZ*dZ);
      }
      double dist = 0;
      // time can't be negative, if it is - then the optimal time is on t = 0
      if (t < 0 || before) {
         t = 0;
         dist = sqrt(sX*sX + sY*sY + sZ*sZ);
      } else {
         dist = sqrt((sX+dX*t)*(sX+dX*t)+(sY+dY*t)*(sY+dY*t)+(sZ+dZ*t)*(sZ+dZ*t));
      }
      cout << "Case #" << (tt+1) << ": ";
      cout << dist << " " << t << "\n";
   }
   return 0;
}

Problem C: Bribe the Prisoners
Category: Dynamic Programming
URL: http://code.google.com/codejam/contest/dashboard?c=189252#s=p2

For the final problem, we are given a set of prison cells (from 1 to P) and inside each of these prison cells is one prisoner. Prisoners communicates with its neighbours if the cell is adjacent and there is at least one prisoner inside. As a prisoner is released, the neighbours find out (if any) and this effect gets echoed in both directions until it hits the end of the prison (i.e. left or right wall) or when there is no one in the neighbouring prison (because they were released earlier). We are then tasked to bribe each prisoner who hears about a prisoner getting released 1 gold coin to prevent them from destroying their cell. We need to identify the optimal sequence of prisoners to release as to minimise the number of gold coins spent.

The obvious way is to brute force all possible ordering of prisoners and compute the minimum gold coins spent. This works easily for the small case as there are only 5 prisoners awaiting to be released at the worst case. This obviously does not work for the large case as there can be as many as 100 prisoners to be released. To optimise this we turn to dynamic programming. Seeing there can be as many as 100 prisoners to be released we can't keep a bitmask of which prisoners we have released as it will certainly TLE and/or exceed memory limit (on all current systems). The key observation is to divide the prison into sub-prisons. This can be seen by noting that the ends of the prison cells act as barriers/walls in which the word does not go through. If we release prisoners from cells - the cells we have just freed implicitly acts like the ends of our original prison.

In addition, we need to make a small optimisation as the prisons can be as large as 10,000 in the large case. We need to note that we will never free from a cell which is not in the list of prisoners to be freed. Hence instead of keeping track of the index of the prison cells itself, we keep track of the index of the prisoners to be freed. A simple sort suffices to keep the implementation neat. Hence the space complexity is only O(Q*Q) and since Q is at most only 100 - this is pretty much nothing.

So we define the DP state as:

D(n,m) = the minimum number of gold spent to free prisoners between index n and m of the prisoners to be freed.

Our recurrence relation is then simply defined as:

D(n,m) = min { Cost(i) + D(n,i) + D(i,m) } where n < i < m. (free prisoner i)
D(n,m) = 0 if m - n < 2 (base case)

The base case is for situations where there are no prisoners to be freed between n and m. Cost(i) is defined as the cost to free prisoner index i given that the empty prison to the left is n and the empty prison to the right is m. This can be calculated in constant time as we know which locations index n and m corresponds to. There is a small trick for implementing this - we add explicit left and right barriers to our original prison to make it conform to our recursion. This way we don't need to handle for the original prison special case which would otherwise be unbounded.

Implementing this should be simple with either memoisation (top-down) or a bottom-up approach:

#include <iostream>
#include <vector>
#include <string>
#include <sstream>
#include <algorithm>
#include <stack>
#include <queue>
#include <map>

using namespace std;

int memo[111][111];

vector<int> data;

int func(int leftIdx, int rightIdx) {
   // base case
   if (rightIdx - leftIdx < 2) {
      return 0;
   }
   if (memo[leftIdx][rightIdx] != -1) return memo[leftIdx][rightIdx];
   int res = 1<<30;
   for (int i = leftIdx+1; i < rightIdx; i++) {
      // split on prisoner index i
      int calc = (data[i] - data[leftIdx] - 1) + (data[rightIdx] - data[i] - 1);
      calc = calc + func(leftIdx, i) + func(i, rightIdx);
      res = min(res,calc);
   }
   return memo[leftIdx][rightIdx] = res;
}

int main() {
   int N;
   cin >> N;
   for (int t = 0; t < N; t++) {
      int P, Q; cin >> P >> Q;
      vector<int> pris;
      for (int i = 0; i < Q; i++) {
         int val; cin >> val; pris.push_back(val);
      }
      // add explicit barriers to the ends of the prisons
      pris.push_back(0);
      pris.push_back(P+1);
      // ensure that the prisoner indexes are strictly increasing
      sort(pris.begin(),pris.end());
      data = pris;
      memset(memo,-1,sizeof(memo));
      int best = func(0, data.size()-1);
      cout << "Case #" << t+1 << ": " << best << "\n";
   }
   return 0;
}

TC TCO08 Round 2 D1 500 (CreaturesTraining)

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=8570

This problem asks us to upgrade a set of creatures to maximise the final score after D days. The score is defined as the sum of the power level of the creatures. This is a tricky problem to solve efficiently. The key observation to make is that you never downgrade creatures and you can also skip creatures. Therefore this suggests that we can go from the leftside to rightside in a linear fashion. We can do a rough sketch of the recursive state as:

D(idx, x) = the maximum number of points we can get when we are processing level idx of the creatures and we have x days left in credit to upgrade our creatures.

When we try to formulate our recurrence relation we run into a problem, the state does not fully describe the problem state as it does not keep track of which creatures we have upgraded so far. So we begin to re-formulate our DP state by incorporating the set of monsters we have upgraded so far, i.e. the modifications we have made to the count array. One obvious yet inefficient way is to keep track of the count array as part of the DP state. However, this is too large of a memory footprint and as a result causes our solution to become inefficient. However, it is a good starting point as we can re-use our observation before that we never look back at a previous level. With this in mind, we can simply just keep track of the current extra creatures we have upgraded so far and optimally decide how many of these to keep in each level. Therefore, our DP state becomes:

D(idx, x, carry) = the maximum number of points we can get when we are processing level idx of the creatures and we have x days left in credit to upgrade our creatures. Given that we are also carrying forth carry number of monsters from the previous level.

In each state we can choose between 0 to count[idx] + carry monsters to upgrade. Each of these are valid transitions/upgrades provided we have enough day credit to upgrade them. We need to check whether or not this state space fits within the memory constraints, as there are a maximum of 50 monster levels and 100 days and carries (as we cannot carry more than 100 as this will exceed the upper limit of day credits). So we are faced with a memory footprint of 50 * 100 * 100 = 500,000 64-bit integers which easily fits within the limit.

Lastly, we consider the base/terminating case, this is when we reach the last level of the creatures. We implicitly prune any transitions where we exceed the amount of remaining day credit, so we do not need to consider this in our terminating case.

Implementing the ideas above is rather simple and short:

class CreatureTraining {
public:
   long long maximumPower(vector <int>, vector <int>, int);
};

long long dp[51][111][111];
vector<int> _count;
vector<int> _power;

long long func(int idx, int day, int carry) {
   if (idx >= _count.size()) return dp[idx][day][carry] = 0;
   if (dp[idx][day][carry] != -1) return dp[idx][day][carry];
   long long res = 0;
   int mx = carry + _count[idx];
   for (int i = 0; i <= mx; i++) {
      // move a set of 0 -> mx monsters up to the next level
      // leave mx - i behind
      if (day - i >= 0) {
         res >?= func(idx+1, day - i, i) + (mx-i) * (long long)_power[idx];
      } else {
         break;
      }
   }
   return dp[idx][day][carry] = res;
}

long long CreatureTraining::maximumPower(vector <int> count, vector <int> power, 
      int D) {
   long long res = 0;
   memset(dp,-1,sizeof(dp));
   limitD = D; _count = count; _power = power;
   return res = func(0, D, 0);
}

TC TCO08 Qualification Round 1 D1 600 (PrimeSums)

Category: Dynamic Programming
URL: http://www.topcoder.com/stat?c=problem_statement&pm=8596

This problem continues with our DP-themed run of recent problems. It asks us to compute the number of subsets (including the empty subset) in which the total bag weight is a prime number. The constraints give us an upper bound of 500,000 (50 * 10,000) for the maximum prime. So the first task is to compute all the prime numbers under this number, this is an easy and efficient task if you know the Sieve of Eratosthenes. The next step is determining how many ways we can make a sub-bag of a weight W, we can use this in our last step by adding up all the weights W which are prime to yield our final answer. Before that though, we need to determine the number of sub-bags which totals a weight of W.

We notice that there are overlapping subproblems which screams to us: dynamic programming! To see this, imagine we have a bag of {1, 1, 2, 7}. To compute the number which weights are possible for all subsets: we start off with the first element. We can choose to use this element and place it into our bag, or we can reject this element and not put it into our bag. Either way, we have {1, 2, 7} as a sub-problem however if we don't do any caching or DP we will compute this sub-problem twice. This effect carries on from the first element to the n-th element and the number of sub-problems we are forced to compute grows exponentially. So we can simply memoise our solution, or can we?

Looking at the constraints we need to keep track of our current index inside the bag (so we know which element we are considering) and also the current weight of the bag (which has an upper bound of 500,000). Since the number of elements is 50 are the worst case, using a memoisation table of 50 * 500,000 64-bit integers is too large to fit into the memory limit. There must be another way to optimise our space complexity. As with all DP problems, we only need to keep sub-solutions to as much as our recurrence relation needs them. So let's define our recurrence relation in hopes of a good optimisation:

DP State: F(n,w) = the number of sub-bags we can construct using a set of 1..n bags (1-based index) that sums to exactly a weight of w.
DP Recurrence: F(n,w) = F(n-1,w - b[n]) + F(n-1,w)

The recurrence comes from the fact that we are faced with 2 decisions at a particular state, to take the item and place it into the bag (the first recursive call) and to skip the item and retain our weight (the second recursive call). As you can see from the depth of the recurrence relation is only at most 1 (i.e. we only need to keep track of the n-1'th pre-computed values), then we can reduce our memoisation table to 2 * 500,000 64-bit integers and use modulus arithmetic to implicitly swap the rows around. See the DP: State Compression article for more information on this technique.

As usual, we also need to define a base case. This is pretty simple, we are allowed to have an empty subset so naturally there is 1 way to have a weight of 0 without using any elements in the bag. Hence F(0,0) = 1. However, we encounter a major problem when we try and implement the solution. The fact there are non-unique elements in the bag can cause our dynamic programming algorithm to double-count the same sub-bag. For example, consider the bag {1, 1, 2, 7} if we skip the first bag and decide to use all the other bags we yield {1, 2, 7} with a weight of 10. However, we can also choose to use the first bag and skip the second bag to yield the same sub-bag {1, 2, 7} also with a weight of 10. Our current algorithm does not take this into account and hence will produce incorrect results for certain inputs.

A neat way to tackle this problem is to condense the same weights into one bag and keep track of how many elements we have of that weight. In the example used previously, we can make a mapping of { {1,2}, {2,1}, {7,1} } which has a pair in each element denoting the weight of the element and how many times it is inside our bag. Then we apply the same DP algorithm to this modified mapping and introduce additional decisions to include k = 1 to n of each element, where n is defined as the number of times it is in our bag. So we can choose to use 0, 1 or 2 elements of weight 1 but we are only allowed to choose 0 or 1 elements of either weight 2 and 7. Note that this does not cause us to double count as we only allow a transition of 2 elements of weight 1 once in our algorithm.

Lastly, we just need to add up our DP table for all the prime weights which is pretty self-explanatory.

class PrimeSums {
public:
   long long getCount(vector <int>);
};

int prime[500011];
long long dp[2][500011];

long long PrimeSums::getCount(vector <int> bag) {
   long long res = 0;
   sort(bag.begin(),bag.end());
   // sieve to generate the primes
   memset(prime,1,sizeof(prime));
   prime[0] = prime[1] = 0;
   for (int i = 2; i * i <= 500000; i++) {
      if (prime[i]) {
         for (int j = i * i; j <= 500000; j += i) {
            prime[j] = 0;
         }
      }
   }
   // merge identical bag elements together
   map<int,int> mappings;
   for (int i = 0; i < bag.size(); i++) mappings[bag[i]]++;
   
   int kth = 1;
   
   // base case
   dp[0][0] = 1;
   
   for (map<int,int>::iterator it = mappings.begin(); it != mappings.end(); it++, 
   kth++) {
      for (int i = 0; i <= 500000; i++) {
         dp[kth&1][i] = dp[!(kth&1)][i];
         // loop through possible weight arrangements based on the available set
         for (int j = 1; j <= it->second; j++) {
            if (i-j*it->first < 0) break;
            dp[kth&1][i] += dp[!(kth&1)][i-j*it->first];
         }
      }
   }
   // sum up all the number of sub-bags which have a prime weight
   for (int j = 0; j <= 500000; j++) {
      if (prime[j]) res += dp[!(kth&1)][j];
   }
   return res;
}

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;
}