TC SRM 487 D1-500 KiwiJuice

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

class KiwiJuice {
public:
   int theProfit(int, vector <int>, vector <int>);
};

int dp[51][1<<16];
vector<int> data;
vector<int> price;
int _C;

// f(current volume of the special bottle, bitmask of the used bottles)
int f(int i, int mask) {
   // caching mechanism
   if (dp[i][mask] != -1) return dp[i][mask];
   // base case
   if (mask+1 == (1 << (data.size()))-1) {
      return price[i];
   }
   int res = 0;
   for (int j = 1; j < data.size(); j++) {
      if (mask & (1 << j)) continue;   // already used
      // the next bottle to process
      int next = data[j];
      int d1 = 0;
      if (i + next >= _C) {
         // one of the bottles is full, the other is i+next-_C volume
         d1 = price[_C] + f(i+next-_C, mask | (1 << j));
      } else {
         // change to i + next (and used up the empty for bottle[j])
         d1 = price[0] + f(i+next, mask | (1 << j));
      }
      // don't add more to our current bottle
      // we swap with bottle j and "use" up our bottle
      int d2 = f(next, mask | (1 << j)) + price[i];
      
      // take the maximum of the decision branches
      res >?= max(d1,d2);
   }
   return dp[i][mask] = res;
}

int KiwiJuice::theProfit(int C, vector <int> bottles, vector <int> prices) {
   price = prices;
   int res = 0;
   _C = C;
   memset(dp,-1,sizeof(dp));
   data = bottles;
   return res = f(data[0],0);   // arbitrarily assign bottle 0 to be our special bottle
}

TC SRM 487 D1-250 CarrotJumping

Category: Brute Force, Mathematics
URL: http://www.topcoder.com/stat?c=problem_statement&pm=11022

class CarrotJumping {
public:
   int theJump(int);
};

long long MOD = 1000000007;

// compute 2^n (mod MOD)
long long pow(int v) {
   if (v == 0) return 1;
   if (v == 1) return 2;
   long long x = pow(v/2);
   if (v % 2 != 0) return (x * x * 2) % MOD;
   return (x * x) % MOD;
}

int CarrotJumping::theJump(int init) {
   int res = 100001;
   for (int a = 0; a <= 100000; a++) {
      for (int b = 0; b < 3; b++) {
         long long v = pow(a*3 + b*2);
         long long c = (v * init + (v - 1 + MOD) % MOD) % MOD;
         if (c == 0) {
            res <?= (a + b);
         }
      }
   }
   return res >= 100001 ? -1 : res;
}

SPOJ - Assignments (ASSIGN)

URL: http://www.spoj.pl/problems/ASSIGN/
Category: Dynamic Programming, Bit Manipulation

The problem is quite simple: to calculate the number of different assignments of n topics to n students provided that every student gets exactly 1 topic that he likes. Each student likes different topics which is given in the format of a matrix. The essence of the problem lies in optimisation, a brute force algorithm will clearly time out if that isn't already evident. As a note, the maximum number of assignments for a 20x20 (filled with 1's) matrix is a massive 2,432,902,008,176,640,000. Since we are enumerating we can easily eliminate greedy as a design approach. Naturally we turn to dynamic programming to somehow reduce the running time to an acceptable amount.

As with all dynamic programming problems, we first define the state. The most natural and obvious state is:

D(s,u) = amount of ways to assign up to student s with a bitmask of subjects u used.

The space complexity of this DP state is O(N * 2^N). However, a observation allows us to reduce the memory footprint by a factor of N, giving a complexity of O(2^N). How does this work? The unique property of the problem is that the number of students and subjects are the same, so we can deduce that when we are up to student V, there are also exactly V bits set in the bitmask. Hence, we can simply just use the bitmask to determine which student we are up to!

Now let's define the recurrence relation:

F(u) = sum of all i { F(u | (1 << i) } where student count_of_bits(u) likes subject i

The base case is,

F((1 << total_students)-1) = 1

Note we have done something implicit. That is, we do not define the base case for where we have "left over" subjects that weren't assign to a student - this is implicitly handled by our looping.

We then code a memoised solution using the recurrence above:

long long memo[1<<20];

int count_bits(int n) {
   int r = 0;
   while (n > 0) { r += (n&1); n >>= 1; }
   return r;
}

long long func(int topicMask) {
   int idx = count_bits(topicMask);
   if (idx >= size) {
      if (topicMask != (1 << size)-1) return 0;
      return 1;
   }
   if (memo[topicMask] != -1) return memo[topicMask];
   long long res = 0;
   for (int i = 0; i < size; i++) {
      if (adjmat[idx][i] == 0 || topicMask & (1 << i)) {
         continue;
      }
      res += func(topicMask | (1 << i));
   }
   return memo[topicMask] = res;
}

Unfortunately, this times out in SPOJ, i.e. it exceeds over 20 seconds for all test cases. This is where bottom-up dynamic programming shines, given that our algorithm is probably close to the time limit (you can check by running on your own machine and noting the speed - if it's quite fast but times out in SPOJ, you are within a constant time factor of passing). Now turning it into a bottom-up dynamic programming is a bit more tricky than usual. The problem lies with the order on which to evaluate.

If we set the same base case as before with,

F((1 << total_students)-1) = 1

We need to somehow traverse backwards whilst also ensuring that we have computed all DP states that are one bit away from the one we are calculating. Turns out due to binary number representation, traversing from (1 << total_students)-1 to 0 is sufficient. This is a rather easy proof and hence left to the exercise of the reader (just list them out for 4-bit numbers and you should see the connection).

The implementation is detailed below:

int adjmat[21][21];
int size;
long long memo[1<<20];

int count_bits(int n) {
   int r = 0;
   while (n > 0) { r += (n&1); n >>= 1; }
   return r;
}

int main() {
   memset(adjmat,0,sizeof(adjmat));
   ios_base::sync_with_stdio(false);
   int nT;
   cin >> nT;
   while (nT-- && cin >> size) {
      memset(adjmat,0,sizeof(adjmat));
      memset(memo,0,sizeof(memo));
      for (int i = 0; i < size; i++) {
         for (int j = 0; j < size; j++) {
            int val; cin >> val;
            adjmat[i][j] = val;
         }
      }
      memo[(1<<size)-1] = 1;
      for (int j = (1 << size)-1; j >= 0; j--) {
         int idx = count_bits(j);
         for (int k = 0; k < size; k++) {
            if (adjmat[idx][k] == 0 || (j & (1 << k))) continue;
            memo[j] += memo[j | (1 << k)];
         }
      }
      cout << memo[0] << "\n";
   }
   return 0;
}

We can shave off some extra time by using a constant time bit counting algorithm (you can google it):

#include <iostream>
#include <iostream>
#include <vector>
#include <algorithm>
#include <cstring>

using namespace std;

int adjmat[21][21];
int size;
long long memo[1<<20];

unsigned numbits(unsigned i)
{
    const unsigned MASK1  = 0x55555555;
    const unsigned MASK2  = 0x33333333;
    const unsigned MASK4  = 0x0f0f0f0f;
    const unsigned MASK8  = 0x00ff00ff;
    const unsigned MASK16 = 0x0000ffff;

    i = (i&MASK1 ) + (i>>1 &MASK1 );
    i = (i&MASK2 ) + (i>>2 &MASK2 );
    i = (i&MASK4 ) + (i>>4 &MASK4 );
    i = (i&MASK8 ) + (i>>8 &MASK8 );
    i = (i&MASK16) + (i>>16&MASK16);

    return i;
}

int main() {
   memset(adjmat,0,sizeof(adjmat));
   ios_base::sync_with_stdio(false);
   int nT;
   cin >> nT;
   while (nT-- && cin >> size) {
      memset(adjmat,0,sizeof(adjmat));
      memset(memo,0,sizeof(memo));
      for (int i = 0; i < size; i++) {
         for (int j = 0; j < size; j++) {
            int val; cin >> val;
            adjmat[i][j] = val;
         }
      }
      memo[(1<<size)-1] = 1;
      for (int j = (1 << size)-1; j >= 0; j--) {
         int idx = numbits(j);
         for (int k = 0; k < size; k++) {
            if (adjmat[idx][k] == 0 || (j & (1 << k))) continue;
            memo[j] += memo[j | (1 << k)];
         }
      }
      cout << memo[0] << "\n";
   }
   return 0;
}

TC TCO10 Qualification Round 2

Overview:

This was the second qualification round for TCO10 with around 1300 competitors. As 600 competitors advance through to the next round (provided that at least 600 score a positive score), it provided ample opportunity to qualify for Round 1. The problems were rather simple with two greedy problems and a simple DP problem as the 1000 pointer. Having said that, the system test claimed a decent portion of the 500 pointer solutions due to a corner case which wasn't covered in the examples.

Problem 1 (250 Points): JingleRingle
Category: Greedy
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10896

This problem required us to calculate the maximum profit we can make by buying jingles and selling them to various buyers. We are assumed to have infinite wealth so we can essentially buy all the jingles but we also want to maximise our profits.

The first key observation is to note that we can only use each buyer or seller once. The second key observation is that we only buy jingles from sellers if we can find a suitable buyer in which we can make a net profit (after tax). Now the question that remains is which sellers do we buy from and whom do we sell it to. Since each seller only gives us 1 jingle and we can only sell back 1 jingle at a time, this constraint essentially means that we can match a seller to a buyer (i.e. we are buying at a lower price and re-selling it at a higher price for a net profit). A simple greedy strategy would be then to buy at the lowest price and selling it to the highest price for maximum profit.

Does this greedy strategy work? Let's informally prove it. As we need to match a seller and buyer, if we buy from N people then we also have to sell it to N people. To maximise our profit we need to maximise our revenue, so the optimal choice is to choose the N highest buyers. Then to minimise our expenses we need to choose the N lowest sellers of jingles. It turns out it doesn't matter what order we pair the buyers and sellers as the total profit is the same if we choose the N highest buyers and N lowest sellers. Now all we need is to determine the value of N. Our greedy algorithm keeps pairing the highest available buyer to the lowest available seller until we reach a situation where we no longer profit, this is the value of N. If we assume this is not the value of N, then we have to somehow match a lower buyer to a higher seller or vice versa. But our value of N is at least as good as this on all cases hence our deduced value of N is correct.

A straightforward implementation:

class JingleRingle {
public:
   int profit(vector <int>, vector <int>, int);
};

int calcTax(int p, int tax) {
   return (p * tax) / 100;
}

int JingleRingle::profit(vector <int> buyOffers, 
      vector <int> sellOffers, int tax) {
   int res = 0;
   int buyJ = buyOffers.size()-1;
   sort(buyOffers.begin(),buyOffers.end());
   sort(sellOffers.begin(),sellOffers.end());
   for (int i = 0; i < sellOffers.size(); i++) {
      if (buyJ >= 0 && buyOffers[buyJ] - calcTax(buyOffers[buyJ],tax) 
         - sellOffers[i] >= 0) {
         res += buyOffers[buyJ] - calcTax(buyOffers[buyJ],tax) - sellOffers[i];
      }
      buyJ--;
   }
   return res;
}

Problem 2 (500 Points): FuzzyLife
Category: Greedy, Brute Force
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10911

This problem has us calculate the maximum number of live cells on an infinite 2D plane after one iteration of time. The 2D plane is consisted of cells which are considered either dead or alive at a given slice in time. However, there are unknown cells denoted as '?' and our job is to determine whether each of these cells should be considered alive or dead for the next evolution in time to have the maximum number of alive cells.

The rules are:
- Each live cell with less than 2 or more than 3 live neighbours are dead
- Each dead cell with exactly 3 live neighbours becomes alive
- All other cell remain the same

Each cell has exactly 8 neighbours (even the ones on the "border" of the input because of the infinite plane) - 2 horizontal, 2 vertical and 4 diagonal.

The main observation required for this problem is the constraints of the problem allow us to adopt a greedy strategy and consider the neighbouring cells of a '?' as its own problem subinstance. As no cell will have more than 1 neighbour of type '?' this means that every cell is at most affected by 1 unknown cell. All we need to do is try each '?' cell by placing a '1' and a '0' in its place and simulate how many alive cells we are left with on the next iteration. We then simply choose the maximum of these two choices.

One tricky case is that we also need to take into account the '?' cell, if we place a '1' in it's position does it become alive or dead in the next iteration - vice versa if we place a '0'. We need to add this to our alive count before we compare which one is better. The other tricky case was to add a padding of dead cells around the grid as some of these can turn into alive cells (which need to be included in the answer) - this case however, was covered in the example cases. After we have filled in each of the '?' cells we can then do a single iteration pass and work out the number of alive cells.

The implementation using the ideas above:

class FuzzyLife {
public:
   int survivingCells(vector <string>);
};

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

int FuzzyLife::survivingCells(vector <string> grid) {
   int res = 0;
   // pad the grid with extra 0's around the perimeter
   string str(grid[0].size()+2,'0');
   for (int i = 0; i < grid.size(); i++) {
      grid[i] = '0' + grid[i] + '0';
   }
   grid.insert(grid.begin(),str);
   grid.push_back(str);
   for (int i = 0; i < grid.size(); i++) {
      for (int j = 0; j < grid[i].size(); j++) {
         if (grid[i][j] == '?') {
            // compute options ('?' = 1)
            grid[i][j] = '1';
            int c1 = 0;
            for (int m = 0; m < 8; m++) {
               int neighX = i + dx[m];
               int neighY = j + dy[m];
               if (neighX < 0 || neighY < 0 || neighX >= grid.size() 
               || neighY >= grid[0].size()) continue;
               int alive = 0;
               for (int n = 0; n < 8; n++) {
                  int mX = neighX + dx[n];
                  int mY = neighY + dy[n];
                  if (mX < 0 || mY < 0 || mX >= grid.size() 
                  || mY >= grid[0].size()) continue;
                  if (grid[mX][mY] == '1') alive++;
               }
               if (grid[neighX][neighY] == '1' && (alive == 2 || alive == 3)) c1++;
               else if (grid[neighX][neighY] == '0' && (alive == 3)) c1++;
            }
            // consider the square itself under '?' = 1
            int selfAlive = 0;
            for (int n = 0; n < 8; n++) {
               int selfX = i + dx[n];
               int selfY = j + dy[n];
               if (selfX < 0 || selfY < 0 || selfX >= grid.size() 
               || selfY >= grid[0].size()) continue;
               if (grid[selfX][selfY] == '1') selfAlive++; 
            }
            if (selfAlive == 2 || selfAlive == 3) c1++;
            // compute options ('?' = 0)
            grid[i][j] = '0';
            int c2 = 0;
            for (int m = 0; m < 8; m++) {
               int neighX = i + dx[m];
               int neighY = j + dy[m];
               if (neighX < 0 || neighY < 0 || neighX >= grid.size() 
               || neighY >= grid[0].size()) continue;
               int alive = 0;
               for (int n = 0; n < 8; n++) {
                  int mX = neighX + dx[n];
                  int mY = neighY + dy[n];
                  if (mX < 0 || mY < 0 || mX >= grid.size() 
                  || mY >= grid[0].size()) continue;
                  if (grid[mX][mY] == '1') alive++;
               }
               if (grid[neighX][neighY] == '1' && (alive == 2 || alive == 3)) c2++;
               else if (grid[neighX][neighY] == '0' && (alive == 3)) c2++;
            }
            // consider the square itself under '?' = 0
            selfAlive = 0;
            for (int n = 0; n < 8; n++) {
               int selfX = i + dx[n];
               int selfY = j + dy[n];
               if (selfX < 0 || selfY < 0 || selfX >= grid.size() 
               || selfY >= grid[0].size()) continue;
               if (grid[selfX][selfY] == '1') selfAlive++; 
            }
            if (selfAlive == 3) c2++;
            // choose the one with better live cells
            if (c1 >= c2) grid[i][j] = '1';
         }
      }
   }   
   // iterate one period and determine number of alive
   for (int i = 0; i < grid.size(); i++) {
      for (int j = 0; j < grid[i].size(); j++) {
         int alive = 0;
         for (int n = 0; n < 8; n++) {
            int mX = i + dx[n];
            int mY = j + dy[n];
            if (mX < 0 || mY < 0 || mX >= grid.size() 
            || mY >= grid[0].size()) continue;
            if (grid[mX][mY] == '1') alive++;
         }
         if (grid[i][j] == '1' && (alive == 2 || alive == 3)) res++;
         else if (grid[i][j] == '0' && alive == 3) res++;
      }
   }
   return res;
}

Problem 3 (1000 Points): HandlesSpelling
Category: Dynamic Programming, String Manipulation
URL: http://www.topcoder.com/stat?c=problem_statement&pm=10864

This problem required us to play a game and work out the maximum score based on a set of choices we are forced to make. The game is that we have a string of letters in which we can "stamp" a set of pre-determined strings onto it. We need to maximise the score based on A^2 - B where A is the longest contiguous sequence of letters covered by the stampings and B is the number of uncovered letters.

For example, if the string we are given is HELLO and we are given 3 badges namely: E, HE, L. We can stamp them as follows: {HE}{L}{L}O. We stamped HE with the second badge and the 2 L's with the third badge. The value of A is 4 because we have a contiguous sequence of stampings (HELL) and B is 1 because we didn't have a way to stamp the letter O. Hence the score given this set of operations is 4^2 - 1 = 15, which turns out to be the maximum score we can get.

One strategy is to determine the maximum value of A. Why? A is where our score gets maximised, no other measure increases our score except for the value of A so it's only natural that we want to get A as high as possible. We also note that the maximum contiguous sequence (i.e. maximum value A) is in the optimal solution. Why? Consider two configurations:

xxxxx....----n----....xxxxxxxx
xxxxx--m--.......--m--xxxxxxxx

Here the x's denote similar configurations, so we only need to focus on maximising the centre part of m's, n's and dots. We let n be the maximum contiguous length possible in our string, let that be denoted as length x. Let's assign m to be length x-1. Then our goal is to prove that it is always optimal to choose n over the 2 m's.

Then the score for our first configuration is:
F = (x * x) - (x-1) - (x-1) = x^2 - 2x + 2 (length x ^ 2 minus two sections of x-1 uncovered)

The score for our second configuration is:
G = (x-1) * (x-1) - x = x^2 - 3x + 1 (length (x-1) ^ 2 minus one section of x uncovered)

Hence for all positive values of length x, the score for the first is better than the second. In fact, this is a stricter argument than what is required, as if the sections of n and m were disjoint they'll both be covered, hence there must be some conflict between n and m which forces us to choose between one or the other. Our argument will hold for such a scenario.

Our first task is then to find the maximum contiguous sequence of our string - this can be done easily through basic dynamic programming. Then we are tasked with trying all configurations which have the same maximum contiguous sequence to determine the minimum number of uncovered letters. Again, this can be done via a similar dynamic programming algorithm where we divide our main string into a left and right sub-string and count the maximum number of covered cells we can get. Then we simply subtract this by the size of the left and right sub-string respectively to get the number of uncovered letters.

The implementation which uses the strategy above can be seen as:

class HandlesSpelling {
public:
   int spellIt(vector <string>, vector <string>);
};

#define UNI -12345678
#define INF 12345678

vector<string> _badges;
string _S;

// returns true if str is a prefix to S starting at idx
bool isPrefix(int idx, const string& str, string S) {
   if (idx + str.size() > S.size()) return false;
   for (int i = idx; i < idx + str.size(); i++) {
      if (S[i] != str[i-idx]) return false;
   }
   return true;
}

int prefix[1001][51];   // true if badge j is a prefix at index i
int A[1001];   // keeps track of the longest sequence ending with index i
int B[1001];   // left hand side
int C[1001];   // right hand side
int maxSize;

// maximise coverage using ptr as the memoisation table
int func(int idx, int* ptr) {
   int res = 0;
   if (idx >= maxSize) return 0;
   if (ptr[idx] != -1) return ptr[idx];
   for (int i = 0; i < _badges.size(); i++) {
      if (prefix[idx][i] && idx + _badges[i].size() <= maxSize) {
         res >?= func(idx + _badges[i].size(), ptr) + _badges[i].size();
      }
   }
   res >?= func(idx+1, ptr);
   return ptr[idx] = res;
}

int HandlesSpelling::spellIt(vector <string> parts, vector <string> badges) {
   string data;
   for (int i = 0; i < parts.size(); i++) data += parts[i];
   _badges = badges;
   _S = data;
   // pre-compute the prefixes for speed
   for (int i = 0; i < data.size(); i++) {
      for (int j = 0; j < badges.size(); j++) {
         if (isPrefix(i, badges[j], _S)) { prefix[i][j] = 1; }
      }
   }
   // determine the longest contiguous sequence in our string
   for (int i = _S.size()-1; i >= 0; i--) {
      for (int j = 0; j < _badges.size(); j++) {
         if (isPrefix(i, _badges[j], _S)) {
            A[i] >?= A[i+_badges[j].size()] + _badges[j].size();
         }
      }
   }
   int longest = 0;
   for (int i = 0; i < _S.size(); i++) {
      longest >?= A[i];
   }
   // choose the best score out of all combinations which use A[i] == longest
   int res = UNI;
   for (int i = 0; i < _S.size(); i++) {
      if (A[i] == longest) {
         // branch off to left and right
         string left = _S.substr(0, i);
         string right = _S.substr(i+A[i]);
         // reset memo table for each instance
         memset(B,-1,sizeof(B));
         memset(C,-1,sizeof(C));
         int lOffset = 0, rOffset = 0;
         maxSize = i;
         if (left.size() > 0) lOffset = left.size() - func(0,B);
         maxSize = _S.size();
         if (right.size() > 0) rOffset = right.size() - func(i+A[i],C);
         res >?= (A[i] * A[i]) - (lOffset) - (rOffset);
      }
   }
   return res;
}

USACO Training Gateway - "A Game"

Problem Statement:

"Consider the following two-player game played with a sequence of N positive integers (2 <= N <= 100) laid onto a game board. Player 1 starts the game. The players move alternately by selecting a number from either the left or the right end of the sequence. That number is then deleted from the board, and its value is added to the score of the player who selected it. A player wins if his sum is greater than his opponents.

Write a program that implements the optimal strategy. The optimal strategy yields maximum points when playing against the "best possible" opponent. Your program must further implement an optimal strategy for player 2."

Strategy:

This is a simple dynamic programming problem which involves the theme about games (two players) which is quite common. As no number in the sequence can be skipped, it is sufficient enough to compute one player's score and derive the other player's score by subtracting it from the total sum of the sequence on the board. This allows us to condense our DP state by a factor of 2 by keeping track of only one player's progress.

As with all DP problems we first define the DP state. Here an obvious one can be defined as:

D(n,m) = the best score for player 1 in which the board state is a contiguous subsequence between index n and m. Player 1's turn.
E(n,m) = the best score for player 1 in which the board state is a contiguous subsequence between index n and m. Player 2's turn.

As the maximum board size is only 100, an O(n^2) space complexity easily fits within the memory limits. Next we define our recurrence relation:

For the current player we can choose either the leftmost or the rightmost block, i.e. A[n] or A[m]. The next turn in the game is player 2's turn which is denoted by the E(x,y) relation. The obvious goal of player 2 is to minimise the gain of player 1 which will result in a higher score for him/herself.

Using the two recurrences we can combine them into one to get (using direct substitution from E(n,m) to D(n,m)):

Hence, we have constructed our recurrence for player 1. We keep track of the total sum of the sequence of the board during the input/parsing phase. Player 2's score is simply total sum - D(0, size-1). Lastly, we need to handle our base cases. If n = m, then we are left with one choice - choosing A[n]. If n > m, which means we have finished the game as the current board state is non-existent, then we return 0 as neither player can choose any more integers from the sequence. Hence explicitly the base cases are:

Implementation:

Using the idea above we generate the simple memoisation algorithm:

int dp[101][101];
vector<int> vec;

int func(int n, int m) {
   if (n > m) return 0;
   if (n == m) return dp[n][m] = vec[n];
   if (dp[n][m] != -1) return dp[n][m];
   int res = 0;
   res = max(min(func(n+2,m),func(n+1,m-1))+vec[n],
      min(func(n+1,m-1),func(n,m-2))+vec[m]);
   return dp[n][m] = res;
}

int main() {
   ifstream inFile("game1.in");
   ofstream outFile("game1.out");
   int N; inFile >> N;
   int sum = 0;
   for (int c = 0; c < N; c++) {
      int val; inFile >> val;
      sum += val;
      vec.push_back(val);
   }
   memset(dp,-1,sizeof(dp));
   outFile << func(0,N-1) << " " << sum - func(0,N-1) << "\n";
   return 0;
}

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

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