Pages

Monday, August 31, 2009

Binary Search Algorithm

We all know that binary search can be used to efficiently find whether or not a specific element is in a sorted array/sequence. The time complexity of binary search is O(log n) which means it's very fast even for extremely large ranges. However in a generalised sense, we can apply binary search to any monotonically increasing sequence or abstract function.

Let a function f(x) return true or false according to a specific criteria. If for all values of f(x) <= f(x+1) then we can use binary search to find the smallest spot in which the value changes from false to true in O(log n) time. In relation to searching for a specific element in a sorted sequence - we can define f(x,s) as follows:

f(x,s) = false if x is smaller than s (the term we want to search for)
f(x,s) = true otherwise

Then if we attempt to map out the function for all values of x we yield a sequence like:

f(1,6) f(2,6) f(3,6) f(4,6) f(5,6) f(6,6) ...
false, false, false, false, false, true, true, true, ...

Algorithmically, for f(x,s) being false we move up the lower bound to the middle (low + high / 2). Conversely, for f(x,s) being true we move down the lower bound to the middle. Hence we can make a generic binary search function like so:

while (hi > lo + 1) {
   long long mid = (lo + hi) / 2;
   if (func(mid)) {
      hi = mid;
   } else {
      lo = mid;
   }
}

output hi as the answer

As a concrete example, we will solve a problem using this generalised version of binary search:

Problem: Mortgage
Source: Topcoder SRM 189 D2 1000
URL: http://www.topcoder.com/stat?c=problem_statement&pm=2427&rd=4765

This problem asks us to find the minimum monthly payment we need to make to fulfill our loan without exceeding the number of terms we are allowed to repay in. The first thing to look for in a binary search problem is whether or not the monotonicity holds. In this case, if we can pay the loan amount using $x under n terms then we can also do the same for any loan amount greater than $x. We begin by defining our binary search boolean function:

Let F(x) = false if we can't pay the loan amount using $x under n terms
Let F(x) = true otherwise

Then we can see the monotonic sequence as:

false, false, false, ..., false, true, ... true

Our task is then to find the x value in which the value of F(x) changes from false to true. We can use a sequential search for this by iteratively going from 1 to x. However given the fact that the worst case scenario for this problem is at least 2 billion - this would surely TLE. Therefore, we turn to binary search in which we can solve the problem with a handful of iterations.

Using our general binary search template given above, we can easily come up with an efficient implementation. We also note the need for 64-bit integers due to overflow problems. Another minor note is that we need to terminate if what we pay is not enough to cover the interest - as in these cases the number of terms required to pay back is infinite (and hence false). The rounding can simply be done by adding 11999 before we divide - this will ensure that we always round up to the next integer as require by the problem.

The implementation is below:
class Mortgage {
public:
  int monthlyPayment(int, int, int);
};

bool calcTime(long long loan, int interest, long long payment, int terms) {
  long long current = loan;
  while (current > 0 && terms > 0) {
    current -= payment;
    if (current <= 0) return true;
    long long d = current + (current * interest + 11999)/12000;
    if (d > current + payment) return false;
    current = d;
    terms--;
  }
  return false;
}

int Mortgage::monthlyPayment(int loan, int interest, int term) {
  long long lo = 1, hi = 2000000000;
  while (hi > lo + 1) {
    long long mid = (lo + hi) / 2;
    if (calcTime(loan,interest,mid,term*12)) {
      hi = mid;
    } else {
      lo = mid;
    }
  }
  return hi;
}

On a side note, the binary search function can also be applied to real (double) numbers. However, in this case due to rounding and machine precision problems it becomes somewhat risky to use straight arithmetic to determine the stopping case. In such cases, it is simpler and less risky to have a maximum number of iteration the loop can run for. Although this is a bit dodgy, given a suitable amount of iterations the answer will be precise as the binary search algorithm converges very fast. A template of such a binary search is given below:

#define MAX_ITER 1000
#define EPS 1e-09

while (fabs(hi-lo) > EPS && iter < MAX_ITER) {
   double mid = (lo + hi) / 2;
   if (func(mid)) {
      hi = mid;
   } else {
      lo = mid;
   }
   iter++;
}
output hi as the answer

Convex Hull Algorithm

The problem:

Given a set of points in a 2D plane determine the smallest encapsulating convex polygon which includes all the points. Conceptually you can think of this as wrapping a rubber band on the set of points and the line segments which make this rubber band is the convex hull. See the figure below:


Applications:

Normally the convex hull is used as a pre-processing step for many computational geometry problems. For example to determine the minimum area convex polygon required to include a given set of points, we first proceed to calculate the convex hull and use the hull as input to a polygon area computation algorithm. We'll see an application of this later.

Algorithm:

There are many algorithms which solve the convex hull problem. Perhaps the simplest is to do a naive brute force. The only way a point belongs to the convex hull is if the point isn't contained in any triangle based off another set of points. Therefore for each point in the set we can iterate all possible 3 points and check, this takes linear time for the set and O(n^3) to check for all possible triangles - hence this approach takes O(n^4).

What we will look at is an efficient yet simple approach which takes O(n log n) time. The algorithm is called Graham's Scan and it makes use of a fundamental concept in Computational Geometry: cross products. We can use cross products to determine whether a point is clockwise or counter-clockwise relative to the origin (0,0). Thus using this we can transform all our points to the left and bottom-most point by doing vector subtraction.

The basis of Graham's scan is to consider one point at a time using a sorted ordering of the points. The sorted criteria is the polar angle from our fixed point P (the left/bottom-most point in the set) to each of the other points in the set. Then we iterate through each of the points in this set and see if the next point makes a counter-clockwise or clockwise turn. If it's clockwise, it means our last point was non-optimal and we backtrack back until the next point we add is counter-clockwise.

A visualisation is shown below (from Wikipedia):


Note that in Step 1, we have our convex hull as P-A-B. We proceed to try point C, as it results in a counter-clockwise turn (otherwise known as a left turn from the relative position of line segment AB) - we add this to our hull. It now becomes P-A-B-C and we proceed to try point D. From the perspective of line segment BC, point D actually is clockwise (we have to turn right to reach it - you can think of it as an ant traversing through the hull and the direction we need to turn in order to get to the next point). Therefore, we pop point C off our hull and re-compute whether or not point D is counter-clockwise to line segment AB which it turns out to be. Hence we add point D to the hull - after the removal of point C we are left with a hull of P-A-B-D (the image in Step 3). One more step is to link D back to the original point P - implementation wise, we usually append the starting point to act as a sentinel. Hence our hull would be P-A-B-D-P.

To calculate the cross-product in the 2D space:


which results in:


Given two line segments P1P2 and P1P3 we can calculate the coefficients of the above expression:

P1P2 = (x2 - x1)i + (y2 - y1)j => a1 = x2 - x1 and b1 = y2 - y1
P1P3 = (x3 - x1)i + (y3 - y1)j => a2 = x3 - x1 and b2 = y3 - y1

As we are only dealing with the 2D plane, a3 = 0 and b3 = 0. Given the above expression we are left with:

P = ((x2 - x1) * (y3 - y1) - (y2 - y1) * (x3 - x1)) k

Using the sign of the expression P, we can determine if the points are collinear (P evaluates to 0), counter-clockwise (P evaluates to a positive number) or clockwise (P evaluates to a negative number). Note that the actual value doesn't matter we are only interested in the sign - this way we can avoid floating point precision problems which is a big bonus in geometry problems.

If you are curious to why the sign of the expression P dictates the counter-clockwise/clockwise nature of the two vectors, simply look at the definition of the cross product:


As a x b is the thing we just calculated and ab is the magnitude of the vectors a and b respectively. Hence ab will always be positive. Therefore if P is negative, then theta will also be negative and hence will be clockwise. Otherwise if P is positive, then theta will also be positive and hence will be counter-clockwise. Note that the absolute value of theta is between 0 and 180 degrees inclusive.

Now back to the actual Graham Scan algorithm: we need to firstly sort the points in accordance to our fixed P (which will always be in the convex hull as it's the lowest and left-most point). We simply override the comparison operator for the Point class which will automatically sort based on angle without actually calculating. To see this, note that the relationship between the counter-clockwise/clockwise is transitive. As in, if point C is CCW to point B and point B is CCW to point A - this also implies that point C is CCW to point A. Hence this is a sufficient sorting criterion.

We process points in this order and calculate whether or not the point is clockwise or counter-clockwise compared to the last line segment in our convex hull. We stop when we reach our starting point (the sentinel) in which case we end up with the hull.

Implementation of the convex hull is shown below:

#define EPS 1e-09

class Point {
public:
  double x;
  double y;
  Point() { }
  Point(double a, double b) : x(a), y(b) { }
};

Point first_point;

double cross(const Point& a, const Point& b, const Point& c) {
  double res = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
  return res;
}

bool collinear(const Point& a, const Point& b, const Point& c) {
  return fabs(cross(a,b,c)) <= EPS;
}

bool ccw(const Point& a, const Point& b, const Point& c) {
  return cross(a,b,c) > EPS;
}

bool cw(const Point& a, const Point& b, const Point& c) {
  return cross(a,b,c) < -EPS;
}

double distance(const Point& a, const Point& b) {
  return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

// ensure first_point is set before calling it for the sort
bool operator<(const Point& p1, const Point& p2) {
  if (collinear(first_point, p1, p2)) {
    if (distance(first_point, p1) <= distance(first_point, p2)) return true;
    return false;
  }
  if (ccw(first_point,p1,p2)) return true;
  return false;
}

using namespace std;

vector<Point> convexHull(vector<Point> ps) {
  vector<Point> res;
  if (ps.size() <= 3) return res;
  int N = ps.size(), k = 0;
  first_point.x = 1e50;  // sentinel value
  for (int i = 0; i < N; i++) {
    if (ps[i].x < first_point.x || (ps[i].x == first_point.x 
      && ps[i].y < first_point.y)) {
      first_point = ps[i];
      k = i;
    }
  }
  ps.erase(ps.begin()+k);
  res.push_back(first_point);
  sort(ps.begin(),ps.end());
  res.push_back(ps[0]);
  ps.push_back(first_point);
  for (int i = 1; i < N; i++) {
    while (res.size() >= 2 && !ccw(res[res.size()-2],res[res.size()-1],ps[i])) 
      res.pop_back();
    res.push_back(ps[i]);
  }
  return res;
}

Applying the Algorithm (SCUD Busters):

To apply and test the convex hull problem, we take a look at SCUD Busters. This problem can be accessed via: http://acm.pku.edu.cn/JudgeOnline/problem?id=1264

The problem has a set of kingdoms which have a single power station and a set of houses. A set of missiles are then fired which may or may not hit a kingdom, it hits a kingdom if the missile is inside the polygon defined from the set of points belonging to the kingdom. If the missile hits a kingdom - it loses all its power source. We are asked to calculate the total area in which power is lost for all the kingdoms.

First, we need to calculate the convex hull for each kingdom - this is defined by the statement "The wall surrounding a kingdom is the minimal-perimeter wall that completely surrounds all the houses and the power station that comprise a kingdom". Let's say we have calculated the convex hull for each kingdom K, then given a set of missile locations we need to focus on calculating which (if any) kingdom gets hit by the missile. This is determined by iterating through the list of kingdoms and determining if the missile location is inside the convex polygon defined by the convex hull. If it's inside then the kingdom loses power - if it isn't the kingdom is safe. If it's hit we simply need to calculate the area of the convex polygon defined for the kingdom and we need to sum these up to derive the answer. Note that if a kingdom is hit by 2 or more missiles - it doesn't "lose" more power so a checking mechanism is required to make sure we don't double count the areas. A simple visited array suffices.

The implementation is shown below:

#define EPS 1e-09

class Point {
public:
  double x;
  double y;
  Point() { }
  Point(double a, double b) : x(a), y(b) { }
};

Point first_point;

double cross(const Point& a, const Point& b, const Point& c) {
  double res = (b.x - a.x) * (c.y - a.y) - (b.y - a.y) * (c.x - a.x);
  return res;
}

bool collinear(const Point& a, const Point& b, const Point& c) {
  return fabs(cross(a,b,c)) <= EPS;
}

bool ccw(const Point& a, const Point& b, const Point& c) {
  return cross(a,b,c) > EPS;
}

bool cw(const Point& a, const Point& b, const Point& c) {
  return cross(a,b,c) < -EPS;
}

double distance(const Point& a, const Point& b) {
  return sqrt((a.x - b.x) * (a.x - b.x) + (a.y - b.y) * (a.y - b.y));
}

// ensure first_point is set before calling it for the sort
bool operator<(const Point& p1, const Point& p2) {
  if (collinear(first_point, p1, p2)) {
    if (distance(first_point, p1) <= distance(first_point, p2)) return true;
    return false;
  }
  if (ccw(first_point,p1,p2)) return true;
  return false;
}

using namespace std;

vector<Point> convexHull(vector<Point> ps) {
  vector<Point> res;
  if (ps.size() <= 3) return res;
  int N = ps.size(), k = 0;
  first_point.x = 1e50;  // sentinel value
  for (int i = 0; i < N; i++) {
    if (ps[i].x < first_point.x || (ps[i].x == first_point.x 
      && ps[i].y < first_point.y)) {
      first_point = ps[i];
      k = i;
    }
  }
  ps.erase(ps.begin()+k);
  res.push_back(first_point);
  sort(ps.begin(),ps.end());
  res.push_back(ps[0]);
  ps.push_back(first_point);
  for (int i = 1; i < N; i++) {
    while (res.size() >= 2 && !ccw(res[res.size()-2],res[res.size()-1],ps[i])) 
      res.pop_back();
    res.push_back(ps[i]);
  }
  return res;
}

// assumes that the polygon ends with the starting point
double polygonArea(vector<Point> p) {
  double res = 0.0;
  for (int i = 0; i+1 < p.size(); i++) res += (p[i].x*p[i+1].y - p[i+1].x*p[i].y);
  return res / 2.0;
}

bool insideConvex(const vector<Point>& hull, Point p) {
  bool oddNodes = false;
  for (int i = 0; i+1 < hull.size(); i++) {
    if ((hull[i].y < p.y && hull[i+1].y >= p.y) || 
      (hull[i+1].y < p.y && hull[i].y >= p.y)) {
      if (hull[i].x + (p.y - hull[i].y) / (hull[i+1].y - hull[i].y) * 
        (hull[i+1].x - hull[i].x) < p.x) {
        oddNodes = !oddNodes;
      }
    }
  }
  return oddNodes;
}

int destroyed[111];

int main() {
  vector<vector<Point> > points;
  int K;
  // process the kingdoms and add the convex hull to points
  while (cin >> K && K != -1) {
    int x, y;
    vector<Point> temp;
    for (int i = 0; i < K; i++) {
      cin >> x >> y;
      temp.push_back(Point(x,y));
    }
    points.push_back(convexHull(temp));
  }
  // process the missiles
  int mx, my;
  double res = 0.0;
  while (cin >> mx >> my) {
    Point p(mx,my);
    for (int j = 0; j < points.size(); j++) {
      // if its not destroyed already and inside the kingdom - destroy it
      if (!destroyed[j] && insideConvex(points[j], p)) {
        res += polygonArea(points[j]);
        destroyed[j] = 1;
        break;
      }
    }
  }
  printf("%.2f\n",res);
  return 0;
}

Friday, August 28, 2009

TC SRM447 D1 500 (PeopleYouMayKnow)

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

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

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

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

Implementation below:


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

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

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

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

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

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

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

Thursday, August 27, 2009

TC SRM 447 D2 900 (ImportsList)

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

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

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

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

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


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

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

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

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

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

Wednesday, August 26, 2009

Bipartite Matching Algorithm

The problem:

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

Applications:

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

Algorithm:

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

For example, take:

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

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

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

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

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

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

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

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

Matching F and G would force:

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

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

Other Notes:

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

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

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

Applying the Algorithm (SPOJ TAXI):

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

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

The implementation is shown below:


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

using namespace std;

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

int N;
int M;

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

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

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

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

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

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

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

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

Friday, August 21, 2009

TC MPR1 D1 550 (CantorDust)

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

This problem is pretty much exactly the same as the D2 1000. However, there is one additional case we need to handle - patterns which are empty. This is fairly simple to translate from our previous algorithm/implementation. We simply need to check if the pattern supplied is empty (i.e. no X's) and if so make a special case computation.

Now the question is, how do we extend our original algorithm to cover this case? Take note that if there were no X's in the fractal then the number of different occurrences of our blank pattern would be:

(size - row size of the pattern + 1) * (size - column size of the pattern + 1)

However there are X's which act as obstacles from reaching this upper bound. So we simply calculate the number of occurrences where our pattern does not fit into the fractal. This is rather simple, just calculate the number of occurrences in which our pattern does fit and subtract it from the total number of occurrences for the associated row/column.

The total number of occurrences for a row is: size - row size of the pattern + 1
The total number of occurrences for a column is: size - column size of the pattern + 1
The number of occurrences where the pattern doesn't fit (per row): (size - row size of the pattern + 1) - pattern does fit in the row
The number of occurrences where the pattern doesn't fit (per col): (size - column size of the pattern + 1) - pattern does fit in the column

Using the same multiplication principle, the total number of occurrences where the pattern doesn't fit in the 2D space is simply the product of the two. Lastly, we just need to exclude this many occurrences from our upper bound calculation to yield the answer:

Total number of occurrences of an empty pattern: ((size - row size of the pattern + 1) * (size - column size of the pattern + 1)) - ((size - row size of the pattern + 1) - pattern does fit in the row) * ((size - column size of the pattern + 1) - pattern does fit in the column)

The implementation is as follows:


class CantorDust {
public:
int occurrencesNumber(vector <string>, int);
};

string tab[10];
int visited[10];

string build(int num) {
if (num == 0) return "X";
if (visited[num]) return tab[num];
visited[num] = 1;
return tab[num] = build(num-1) + string(pow(3.0,(double)num-1),'.') + build(num-1);
}

int CantorDust::occurrencesNumber(vector <string> pattern, int t) {
int res = 0;
string str = build(t);
string emptyRow = string(pattern[0].size(),'.');
string emptyCol = string(pattern.size(),'.');
string rows = "";
int size = pow(3.0,(double)t);
for (int i = 0; i < pattern.size(); i++) {
if (pattern[i] == emptyRow) continue;
if (rows != "" && rows != pattern[i]) return 0;
rows = pattern[i];
}
string cols = "";
for (int i = 0; i < pattern[0].size(); i++) {
string col;
for (int j = 0; j < pattern.size(); j++) {
col += pattern[j][i];
}
if (col == emptyCol) continue;
if (cols != "" && cols != col) return 0;
cols = col;
}
bool isEmpty = false;
if (cols == "" && rows == "") {
cols = emptyCol;
rows = emptyRow;
isEmpty = true;
}
int numMatchRow = 0, numMatchCol = 0;
string::size_type pos = 0;
while ((pos = str.find(rows,pos)) != string::npos) {
numMatchRow++;
pos++;
}
pos = 0;
while ((pos = str.find(cols,pos)) != string::npos) {
numMatchCol++;
pos++;
}
if (isEmpty) {
int exclude = (size - emptyRow.size() + 1 -
numMatchRow) * (size - emptyCol.size() + 1 - numMatchCol);
return (size - emptyRow.size() + 1) *
(size - emptyCol.size() + 1) - exclude;
}
return res = numMatchRow * numMatchCol;
}

Thursday, August 20, 2009

TC Member Pilot Round 1 Division 2

This is the first member-driven "Single Round Match" for Topcoder (which also isn't rated). Division 2 had an interesting set of problems which proved to be slightly more difficult than the norm.

Division 2 250-pointer (TriFibonacci)
This problem can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10587&rd=13935

Basically this problem asks us to "fill in the blank" of a "Tri-Fibonacci" sequence which is defined as F(a) = F(a-1) + F(a-2) + F(a-3). The small catch here is that there can be sequences where this isn't well-defined (i.e. invalid). All the integers in the sequence must be strictly positive (i.e. zeroes are not allowed).

To solve this problem there are several options:
1. Brute force the blank and check whether or not it leads to a valid Tri-Fibonacci sequence. Due to the constraints this would take around 20 million operations which should run in time. If we can't find a valid sequence then we return -1.

2. Compute the blank based on the sequence given. Roughly speaking, we can divide this into two cases. One is where the number is within the first three terms and hence need to deduce the value by taking the sum of the other 2 and subtracting it from the next Tri-Fibonacci number.

For example, given {2, -1, 3, 7}. We compute 2+3=5, and then we can deduce the blank by noting 5+x=7 (by definition of the recurrence).

The other case is when the missing term is on the fourth term or beyond. Then it's easy to just use the given recurrence and compute the missing number from the last 3 numbers in the sequence before the blank. We then proceed to do a check on whether or not the unique number we obtained from the sequence is a valid Tri-Fibonacci. If it is, we accept and return the value. Otherwise we return -1.

The solution below uses approach 2:


class TriFibonacci {
public:
int complete(vector <int>);
};

bool checkFib(const vector<int>& A) {
for (int i = 3; i < A.size(); i++) {
if (A[i] != A[i-1] + A[i-2] + A[i-3]) return 0;
}
return 1;
}

int TriFibonacci::complete(vector <int> A) {
int pos = 0;
for (int i = 0; i < A.size(); i++) {
if (A[i] < 0) { pos = i; break; }
}
if (pos < 3) {
int s = 0;
for (int i = 0; i < 3; i++) if (A[i] >= 0) s += A[i];
if (A[3] - s <= 0) return -1;
A[pos] = A[3] - s;
return checkFib(A) ? A[pos] : -1;
}
A[pos] = A[pos-1] + A[pos-2] + A[pos-3];
return checkFib(A) ? A[pos] : -1;
}


Division 2 500-pointer (ErdosNumber)
This problem can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10377&rd=13935

This problem asks us to compute the Erdos Number based on a given set of publications. If you are interested in this number there are some facts on Wikipedia: http://en.wikipedia.org/wiki/Erd%C5%91s_number.

Again like most problems there are several ways to solve this. The bulk of the problem is really just parsing the string into a suitable data structure. This can be done using a map which maps a string (author names) and converts it into an integer index (for the adjacency matrix). We simply read each publication at a time and if it's a new author name we insert it into the map. Once we have read all the authors for a given publication we iterate in pairs and assign a weight of 1 to each edge between the authors. Then we simply need to run a shortest path algorithm from Erdos to all the other authors. There are several options available:

1. Most simplest is to use a Floyd-Warshall algorithm. This runs in O(n^3) time and since there are at most 100 vertices, it will easily run in time.
2. Use a breadth-first search from Erdos keeping track of the depth. This works as all the edges have a weight of 1. Reasonably simple to code and should be of choice if you don't know the Floyd-Warshall's All-Pairs shortest path algorithm. It's also more efficient.
3. Alternative methods include using Dijkstra's from Erdos to all vertices (Warshall's is more efficient and easier to code). You can also use Bellman-Ford's algorithm which will get the answer from the Erdos vertex to all other vertices.

A small implementation note is to keep track if a vertex is not reachable from the Erdos vertex.

See the below implementation which uses Floyd-Warshall's algorithm:


class ErdosNumber {
public:
vector <string> calculateNumbers(vector <string>);
};

#define INF (1<<20)
int adjmat[101][101];

vector <string> ErdosNumber::calculateNumbers(vector <string>
publications) {
vector<string> res;
map<string,int> authorMapping;
int authorIdx = 0;
for (int i = 0; i < 101; i++) for (int j = 0; j < 101; j++)
adjmat[i][j] = i == j ? 0 : INF;
for (int i = 0; i < publications.size(); i++) {
istringstream iss(publications[i]);
string author;
vector<string> assocs;
while (iss >> author) {
if (authorMapping.count(author) == 0) {
authorMapping.insert(make_pair(author,authorIdx++));
}
assocs.push_back(author);
}
for (int j = 0; j < assocs.size(); j++) {
for (int k = j+1; k < assocs.size(); k++) {
int first = authorMapping[assocs[j]];
int second = authorMapping[assocs[k]];
adjmat[first][second] = adjmat[second][first] = 1;
}
}
}
int erdosIdx = authorMapping["ERDOS"];
for (int k = 0; k < authorMapping.size(); k++) {
for (int i = 0; i < authorMapping.size(); i++) {
for (int j = 0; j < authorMapping.size(); j++) {
adjmat[i][j] <?= adjmat[i][k] + adjmat[k][j];
}
}
}
for (map<string,int>::iterator it = authorMapping.begin();
it != authorMapping.end(); it++) {
ostringstream oss;
oss << it->first;
int dist = adjmat[erdosIdx][it->second];
if (dist < INF) oss << " " << dist;
res.push_back(oss.str());
}
return res;
}


Division 2 1000-pointer (CantorDustEasy)
This problem can be viewed:
http://www.topcoder.com/stat?c=problem_statement&pm=10588&rd=13935

Finally the Division 2 Hard problem! This is a somewhat tricky problem and takes some hand-drawn diagrams to really see the picture. The problem is as given: given a recursive fractal algorithm determine the number of occurrences of a specified pattern. Given the low constraints of the pattern (50x50 at most) we can try and calculate the fractal up to 3^4 = 81 > 50 and attempt to recursively build on the answer. Why do we only need up to a time of 4? As each black square "sub-section" is divided into the four corners there is an empty space between them made up of 3^(t-1) blocks. This means that no pattern can possible "cross" over and given the constraint there must be at least 1 black square in the pattern - it makes this approach feasible. This is somewhat difficult to implement though so we attempt to make it easier.

On a side note, you can't naively enumerate all the patterns in a given fractal. The size at 9 is 3^18 which is nearly 400 million. We can't store this much let alone compute pattern matching on each cell which would certainly TLE.

So there is an easier approach but requires some effort to discover. If we draw a few diagrams we can see that all rows are either the same or they are empty (made up entirely of white cells). This is due to the symmetric property of the fractal, likewise - the columns are either the same or they are empty.

But how does this help us count the number of occurrences of a pattern? Firstly it helps eliminate infeasible patterns - if the non-empty rows/columns don't match then we can immediately discard it as a feasible pattern. The next observation we need to note is that our pattern now has either rows and columns that are the same or they are empty. Let's take a non-empty row R, since we know that all of the non-empty rows in the actual fractal F are the same. We can just count the number of R's we see in F. How does this help?

To illustrate the concept, take "X.X...X.X" as a sample row from F. Our pattern's row is: ".X.", using a string finding algorithm we find that it occurs twice in F. Hence if our original (non-collapsed) pattern only had 1 row of ".X." then we immediately retrieve the answer we have. Note this is a 1 dimensional problem - we need to convert this to our 2 dimensional problem. So we do the same thing with columns and determine how many occurrences. Let's say the count we have for the rows is "x" and the count we have for columns is "y". Then by the multiplication principle, the total number of occurrences in the 2D space is simply x * y (draw some diagrams to convince yourself - if we have non-empty rows that have x occurrences each then the total number is just how many times these x occurrences occur in a vertical fashion from the column pattern).

Combining this together we yield a simple implementation:


class CantorDustEasy {
public:
int occurrencesNumber(vector <string>, int);
};

string tab[10];
int visited[10];

string build(int num) {
if (num == 0) return "X";
if (visited[num]) return tab[num];
visited[num] = 1;
return tab[num] = build(num-1) + string(pow(3.0,(double)num-1),'.') + build(num-1);
}

int CantorDustEasy::occurrencesNumber(vector <string> pattern, int t) {
int res = 0;
string str = build(t);
string emptyRow = string(pattern[0].size(),'.');
string emptyCol = string(pattern.size(),'.');
string rows = "";
for (int i = 0; i < pattern.size(); i++) {
if (pattern[i] == emptyRow) continue;
if (rows != "" && rows != pattern[i]) return 0;
rows = pattern[i];
}
string cols = "";
for (int i = 0; i < pattern[0].size(); i++) {
string col;
for (int j = 0; j < pattern.size(); j++) {
col += pattern[j][i];
}
if (col == emptyCol) continue;
if (cols != "" && cols != col) return 0;
cols = col;
}
int numMatchRow = 0, numMatchCol = 0;
string::size_type pos = 0;
while ((pos = str.find(rows,pos)) != string::npos) {
numMatchRow++;
pos++;
}
pos = 0;
while ((pos = str.find(cols,pos)) != string::npos) {
numMatchCol++;
pos++;
}
return res = numMatchRow * numMatchCol;
}

Monday, August 17, 2009

TC SRM443 D1 600 (BinaryFlips)

This problem can be accessed via:
http://www.topcoder.com/stat?c=problem_statement&pm=10387

A fairly standard problem with somewhat tricky implementation due to the problem constraints. The problem gives us "A" number of zeroes and "B" and number of ones and using a swap of exactly "K" digits in each turn, determine the minimum number of moves to get it into a state of all ones (and returning -1 if it's impossible to do so). Looking at the constraints of around 100,000 means we can't naively try all combinations. A simple approach to take is to Breadth-First searching the moves and making sure we don't re-visit already processed nodes with lower number of moves (as they will never be optimal).

We need to make several observations to help simplify the problem. As we can choose any K digits to swap over in a given turn, this means we can view the problem as a set of 0's followed by a set of 1's for simplicity. This helps us determine that the state required for the BFS is simply the number of zeroes (or ones) and not a string of disordered 0's and 1's. A further optimisation is to note that the number of zeroes can be derived from the number of ones as the total size of the string does not change. Hence the maximum state space is simply 200,000 so we are fine for memory usage. The only problem left is the transitions between states inside the inner cycle of the BFS. As the K is as large as 100,000 we will run into time limit issues if we aren't efficient in which states we check/traverse to.

This is perhaps the trickest part of the problem - the other parts are fairly standard and simple. Now we ask ourselves, "is there a way to make sure we only consider non-visited nodes?". Using a set data structure which keeps all its elements unique allows us to only traverse the unvisited states. To make another optimisation we need to ensure that we only visit the reachable unvisited states - we calculate the state range in which we can reach and do a subset traversal of the set (instead of going through the whole set and determining whether it's a valid move). To determine the range of states we need to see that if we can traverse to x-number of 1's and we can also traverse to y-number of 1's then we can also traverse to all number of 1's between x and y with the same parity. To see this, arrange a state in a sequence of 0's and 1's and take an arbitrary K (K = 4 in this case):

Move to 9 1's: 0[0000]11111 -> 0[1111]11111
Move to 7 1's: 00[0001]1111 -> 00[1110]1111
Move to 5 1's: 000[0011]111 -> 000[1100]111
Move to 3 1's: 0000[0111]11 -> 0000[1000]11
Move to 1 1's: 00000[1111]1 -> 00000[0000]1

As you can see, we can only move a specific parity (odd/even) between the minimum and maximum range. This can be seen above where we use a sliding "window" mechanism for the inversion. The parity is determined by a combination of the window size (i.e. K) and the parity of the state. To see this connection, we build an informal proof:

Let the K-sized window be denoted by W.
Let P be the parity of the state S, where S is a sequence of 0's followed by 1's.
Then we define x to be the number of 0's in S and y to be the number of 1's in S.

By definition, P = y (mod 2). Furthermore we need to define x' and y' to be the number of 0's and 1's in the window W respectively. We take x'' and y'' to be the number of 0's and 1's in the state S excluding the digits inside window W. Then it follows that:

x'' = x - x' (mod 2)
y'' = y - y' (mod 2)
x' + y' = K

Parity of the inversed window W = x' (mod 2) since it's the number of 0's in the window (as they turn into 1's after the transform which becomes the parity). We now calculate P', the parity of S after being transformed by an arbitrary K-sized window:

P' = x' (mod 2) + y'' (mod 2) (parity of the inversed window plus the region of S - W)
P' = K - y' (mod 2) + y'' (mod 2)
P' = K - y' + y'' (mod 2)
P' = K - y' + 2y' + y'' (mod 2) (add 2y' - this can be done since it doesn't alter the odd/even-ness of the equation)
P' = K + y' + y'' (mod 2)
P' = K + y (mod 2)

Therefore the next state's parity is P' = K + y (mod 2). So we only need to compute this to determine whether the parity of the next state is odd/even. So it suffices to just calculate the minimum and maximum number of 1's and then traverse based on the parity. Using all the concepts above gives us an implementation which runs at around 100ms worst case:


class BinaryFlips {
public:
int minimalMoves(int, int, int);
};

int visited[200011];

int BinaryFlips::minimalMoves(int A, int B, int K) {
if (A == 0) return 0;
if (A + B < K) return -1;
queue<int> q;
int N = A + B;
memset(visited,-1,sizeof(visited));
set<int> evenVisit, oddVisit;
for (int i = 0; i <= N; i++) {
if (i & 1) oddVisit.insert(i);
else evenVisit.insert(i);
}
if (B & 1) oddVisit.erase(oddVisit.find(B));
else evenVisit.erase(evenVisit.erase(B));
q.push(B);
visited[B] = 0;
while (!q.empty()) {
int t = q.front(); q.pop();
if (t == N) return visited[t];
int numZeroes = N - t;
int numOnes = t;
int x = numOnes < K ? K - numOnes : numOnes - K;
int y = numZeroes < K ? numOnes + numZeroes * 2 - K : numOnes + K;
int s = min(x,y);
int st = max(x,y);
set<int>* visitSet = (t+K) & 1 ? &oddVisit : &evenVisit;
if (visitSet->size() == 0) continue;
set<int>::iterator it = visitSet->lower_bound(s);
while (it != visitSet->end() && *it <= st) {
if (visited[*it] < 0) {
q.push(*it);
visited[*it] = 1 + visited[t];
set<int>::iterator itTemp = it; itTemp++;
visitSet->erase(it);
it = itTemp;
} else {
it++;
}
}
}
return -1;
}

Wednesday, August 12, 2009

TC SRM446 D1 500 (AntOnGraph)

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

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

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

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


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


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

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


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


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


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


Then we can begin collapsing:


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

etc.

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


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


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


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

long long tab[51][51];

#define INF (1LL<<50)

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

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

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

Friday, August 7, 2009

TC SRM358 D1 1000 (SharksDinner)

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

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

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

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


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

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

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

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

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

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

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

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

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

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

Tuesday, August 4, 2009

DP: State Compression

Memory constraints is one of the main problems with dynamic programming, usually given an appropriate state the space complexity is fine but sometimes it isn't sufficient enough. Perhaps we are a small factor off from it being under the memory constraints so we can try and optimise our memory usage by compressing the states together and reusing the previous unneeded table entries.

Let's start with an example, the variant of the subset sum problem can be described as follows:

Given a set of strictly positive numbers which contain possible duplicates which sum up to M, for a given K <= M, determine whether or not there is a subset of the numbers such that they sum up to exactly K.


Obviously our DP table will depend on the maximum sum available so we have to be rather reasonable in this dimension. First, let's construct a possible DP state that we may come up with: Let D(i, k) determine whether or not it is possible to construct a sum of exactly k using item 1 to item i. Then we can appropriately define D(i, k) = 1 if it's possible and D(i, k) = 0 if it's impossible. Defining the base cases as D(0, 0) = 1 (possible to construct a sum of 0 by not choosing anything), and D(0, x) = 0 for x > 0 (as it's not possible to construct a positive sum without using any items).

As our base cases and states are defined, as usual we move to define the recurrence relation. This is rather simple:

F(i, k) = max { F(i-1, k - A[i]), F(i-1, k) }

This shows that we are given 2 decisions, we can either use item i into the set and the only way that this is feasible is that if F(i-1, k-A[i]) was also feasible. Take note that we must ensure that k-A[i] >= 0 otherwise we will run into trouble. The other decision is to not use item i into the set (which is perfectly fine since we are allowed any subset of items from 1 to i). So if we can reach the sum of k using item 1 to item i-1 then we can certainly get to the same sum by choosing not to use item i.

The C++ code is listed below:


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

#define M 10000
#define N 1000

int dp[N+1][M+1];

using namespace std;

int main() {
vector<int> vec;
int val;
int num;
int totalSum = 0;
cin >> num;
while (num-- && cin >> val) { vec.push_back(val); totalSum += val; }
// check for size and sum constraint
if (vec.size() > N || totalSum > M) {
cerr << "Exceeds size and/or sum constraint\n";
return 1;
}
for (int i = 0; i < M+1; i++) dp[0][i] = 0;
dp[0][0] = 1;
for (int i = 1; i <= vec.size(); i++) {
for (int j = 0; j <= totalSum; j++) {
dp[i][j] = dp[i-1][j];
if (j - vec[i-1] >= 0) {
dp[i][j] = max(dp[i][j],dp[i-1][j-vec[i-1]]);
}
}
}
// begin queries
while (cin >> val) {
if (val > M) { cerr << "Query too large. Skipping.\n"; continue; }
string outcome = dp[vec.size()][val] ? "Possible" : "Impossible";
cout << outcome << "\n";
}
return 0;
}



However there is an obvious optimisation we can make to the memory usage. This optimisation can be applied to all DP problems - first we observe the recurrence relation and take note that it only depends on the previous item (i-1) and not anything before that. So we can compress the table by just using 2 rows (or columns depending on how you visualise arrays) and swap them around each iteration and process them again. Alternatively, one can use modulus operations to simulate this effect and this is especially fast for recurrences with a "depth" of 2, as we can simply use a bitwise-AND operation to determine whether its odd or even.

As another example, if our recurrence depended on i-1, i-3, and i-4. Then we need to keep track of a depth of at least 5. We can then use modulus to alternate where we calculate from and where we store newly calculated values to.

Such a simple optimisation reduces our memory footprint dramatically and it's very simple to apply, now our DP space complexity only scales with the maximum sum and has no relation with the number of items we are using.

The C++ code is listed below which uses this optimisation:


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

#define M 10000

int dp[2][M+1];

using namespace std;

int main() {
vector<int> vec;
int val;
int num;
int totalSum = 0;
cin >> num;
while (num-- && cin >> val) { vec.push_back(val); totalSum += val; }
// check for sum constraint
if (totalSum > M) {
cerr << "Exceeds sum constraint\n";
return 1;
}
for (int i = 1; i < M+1; i++) dp[0][i] = 0;
dp[0][0] = 1;
for (int i = 1; i <= vec.size(); i++) {
for (int j = 0; j <= totalSum; j++) {
dp[i&1][j] = dp[(i-1)&1][j];
if (j - vec[i-1] >= 0) {
dp[i&1][j] = max(dp[i&1][j],dp[(i-1)&1][j-vec[i-1]]);
}
}
}
// begin queries
while (cin >> val) {
if (val > M) { cerr << "Query too large. Skipping.\n"; continue; }
string outcome = dp[vec.size()&1][val] ? "Possible" : "Impossible";
cout << outcome << "\n";
}
return 0;
}


Sometimes there is an less obvious space optimisation, such as in this problem. We can effectively further reduce the memory footprint by 2 by only using a 1D array. To see this, it is often useful to draw order dependency diagrams which depict the direction in which dependencies are made in the DP table. Obviously this becomes harder as you reach DP tables with dimensions higher than 3. But as our DP table is only 2D it is easy to visualise on paper:



As you can see the column we wish to calculate (i) is dependent on things equal to or directly below it in row in column i-1. So we can traverse from top to bottom on each column and still be assured that all the entries below the one we are calculating are still effectively column i-1. Hence our correctness is assured and we have essentially compressed the two columns in 1 reducing the memory footprint by a factor of 2.

The C++ code is listed below which uses this optimisation:


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

#define M 10000

int dp[M+1];

using namespace std;

int main() {
vector<int> vec;
int val;
int num;
int totalSum = 0;
cin >> num;
while (num-- && cin >> val) { vec.push_back(val); totalSum += val; }
// check for sum constraint
if (totalSum > M) {
cerr << "Exceeds sum constraint\n";
return 1;
}
for (int i = 1; i < M+1; i++) dp[i] = 0;
dp[0] = 1;
for (int i = 1; i <= vec.size(); i++) {
for (int j = totalSum; j >= 0; j--) {
if (j - vec[i-1] >= 0) {
dp[j] = max(dp[j],dp[j-vec[i-1]]);
}
}
}
// begin queries
while (cin >> val) {
if (val > M) { cerr << "Query too large. Skipping.\n"; continue; }
string outcome = dp[val] ? "Possible" : "Impossible";
cout << outcome << "\n";
}
return 0;
}