Counting partitions

This post is related to part A of Lars' contest. Your job there is to find a (0,1) matrix that maximizes a certain objective function. A (0, 1) matrix is a matrix whose elements are taken from the set {0, 1}. The row sum of a n x m matrix is a n-vector whose i-th element contains the sum of i-th row in the matrix. Similarly you can define the column sum of a matrix. An interesting (standard) problem is to find the necessary and sufficient conditions such that there exist a (0,1) matrix with given row and column sums; and if there is then construct one. But here I'm going to talk about another problem.

I wanted to partition the search space into matrices with certain row/column sums and then proceed with searching these regions in a certain order. The search for each region is implemented. But now I need to order these "regions" by how interesting they are. My first idea was to generate all possible row/column sums, order them by some score and then take the top R regions. In order to get the top R row/column sums I don't need to remember all possibilities and then order them. But I do need to generate them all (at least if I want to keep the definition of a row/column sum score easy to be changed). Fortunately my first instinct was to say: let's count how many possibilities are there to see if it is feasible to generate them all. It turns out they are too many :(

So what I want to present here is a way of counting possible row partitions of a (0,1) matrix. I will also restrict myself to square matrices as in Lars' contest; but the generalization is straightforward. What exactly do we count? The number of solution of the equation: x1+x2+...+xn=s. Here s is the total number of ones in the matrix and n is the number of rows (and columns) in the matrix. Each x is at least 0 and at most n. I will denote the number we search for as p(n, s, m), where m is the maximum value allowable for x, i.e. I have generalized the problem a bit. To paraphrase Dijkstra, sometimes making the problem more general means making it easier to solve.

My first idea was to count the solutions of the equation with sum s, maximum value per term m and number of terms n by counting first solutions with no zero term, then solutions with one zero term, etc. If you write this down you'll see that after very little manipulation (sorry, I hate to write math in HTML -- much worse than LaTeX) the following program works.

   3:let cache = Hashtbl.create 10000000;;
   4:
   5:let rec count n s m = 
   6:  match (s, m) with
   7:    (0, _) -> 1 (* order is relevant! *)
   8:  | (_, 0) -> 0
   9:  | _ -> if n < 0 || s < 0 || m < 0 then 0 else
  10:    try
  11:      Hashtbl.find cache (n, s, m)
  12:    with Not_found ->
  13:      let r = (count (n-1) s m) + (count n (s-n) (m-1)) in
  14:      if r < 0 then failwith "integer overflow" else
  15:      Hashtbl.add cache (n, s, m) r; r;;

Now we can use the above function. For example if we look at matrices that have roughly half elements equal to one and half equal to zero then thair partitions are given by:

  17:let possib n = count n (n * n / 2) n;;

It turns out that possib 18 is about 100 millions, possib 19 is about 400 millions, and possib 20 overflows 231. So there is no chance to generate all possible row sums for matrices of up to 59 rows :( I need a solution specific to the score I'm using for row/column sums...

Topological sort

This entry is still in the spirit of: "let's do some low level processing" inspired by reading SGB. I'll say again that I consider this type of programming "fun" but too risky for real programming. However the main difference is not that big as it might appear at first sight. In a real program I would use std::list to handle adjacency list but otherwise the code would be pretty similar.

The problem I try to solve is a slightly modified topological sort. We are given a graph in an adjacency list representation.

   6:struct Node;
   7:
   8:struct Transition {
   9:    Transition* next;
  10:    Node* dest;
  11:};
  12:
  13:struct Node {
  14:    ~Node() {
  15:        Transition *p, *t = transitions;
  16:        while (p = t->next) {
  17:            delete t;
  18:            t = p;
  19:        }
  20:    }
  21:    
  22:    string name;
  23:    
  24:    Node* rs; // right in stage list
  25:    Node* ls; // left in stage list
  26:    Node* next; // next in graph
  27:    int tag;
  28:    Transition* transitions;
  29:};
  30:
  31:struct Graph {
  32:    Graph() {
  33:        nodes = new Node;
  34:        nodes->ls = nodes->rs = nodes;
  35:    }
  36:    
  37:    ~Graph() {
  38:        Node *p, *t = nodes;
  39:        while (p = t->next) {
  40:            delete t;
  41:            t = p;
  42:        }
  43:    }
  44:    
  45:    Node* nodes; // this is a dummy node
  46:};
  47:
  48:struct Stage {
  49:    Stage* next;
  50:    Node*  n;
  51:};

We want to partition the nodes of the graph in a list of "stages". Each stage is a list of nodes. We want that (1) if the graph contains an edge from s to d then stage[s] comes before stage[d] in the stage list and (2) the stage list is as short as possible. We can use a greedy approach. Why that works is a good exercise. The (untested :( ) code follows

  57:    Stage *r = new Stage, *pr = r; // r is a dummy
  58:    
  59:    // Reset tags and construct initial list
  60:    Node *n, *d;
  61:    for (n = g.nodes->next; n; n = n->next) {
  62:        n->tag = 0;
  63:        (n->ls = g.nodes->ls)->rs = n;
  64:        (n->rs = g.nodes)->ls = n;
  65:    }
  66:    
  67:    // Count inbound edges
  68:    Transition* t;
  69:    for (n = g.nodes->next; n; n = n->next) {
  70:        for (t = n->transitions; t; t = t->next) {
  71:            d = t->dest;
  72:            ++(d->tag);
  73:            d->ls->rs = d->rs;
  74:            d->rs->ls = d->ls;
  75:            d->ls = d->rs = d;
  76:        }
  77:    }
  78:    
  79:    // NOTE: if there are cycles and there is no topological sort this
  80:    // function will return an incomplete result
  81:    while (g.nodes->ls != g.nodes) {
  82:        // Remember the current stage
  83:        pr->next = new Stage;
  84:        pr = pr->next;
  85:        pr->n = g.nodes->ls;
  86:        g.nodes->ls->rs = g.nodes->rs;
  87:        g.nodes->rs->ls = g.nodes->ls;
  88:        g.nodes->rs = g.nodes->ls = g.nodes;
  89:        
  90:        // Compute next stage
  91:        n = pr->n;
  92:        do {
  93:            for (t = n->transitions; t; t = t->next) {
  94:                d = t->dest;
  95:                if (--(d->tag) == 0) {
  96:                    (d->ls = g.nodes->ls)->rs = d;
  97:                    (d->rs = g.nodes)->ls = d;
  98:                }
  99:            }
 100:            n = n->rs;
 101:        } while (n != pr->n);
 102:    }
 103:    
 104:    pr->next = NULL; // this is the last stage
 105:    
 106:    // the result is r

Paul Graham: Lessons from a company

Paul Graham has recently published an essay about start-ups. As usual there are a lots of goodies in there; but one paragraph I want to post here.

I spent a year working for a software company to pay off my college loans. It was the worst year of my adult life, but I learned, without realizing it at the time, a lot of valuable lessons about the software business. In this case they were mostly negative lessons:

1. don't have a lot of meetings
2. don't have chunks of code that multiple people own
3. don't have a sales guy running the company
4. don't make a high-end product
5. don't let your code get too big
6. don't leave finding bugs to QA people
7. don't go too long between releases
8. don't isolate developers from users
9. don't move from Cambridge to Route 128

But negative lessons are just as valuable as positive ones. Perhaps even more valuable: it's hard to repeat a brilliant performance, but it's straightforward to avoid errors.

I agree with most of them wholeheartedly... with a few exceptions. One of them is the code ownership. I feel like I'm in a chains if I can't correct mistakes or simply improve the code written by someone else. I also don't have any problems if someone does the same to my code. At least I think so because in the current project I can rememeber only one occasion when that happen: people are too busy churning out their own new code. The other lesson I'm not sure I agree with is the one about moving from Cambridge to Route 128. But that's because I'm not sure I get it. Does it mean: "stay close (physically) to areas with high densities of smart people"? If so, then I agree.

Size limited cache

I was reading a bit from Knuth's Stanford GraphBase (SGB). He likes to do a lot of pointer handling. In TAOCP he says something along the lines: people avoid handling list themselves because there is a general perception that it is hard to do so; we'll try to dispell that myth. Also he has a paper entitled Dancing Links. Well, Knuth in general likes to do low-level and deep stuff.

Why am I writing about this? Because reading SGB made me a bit nostalgic: I use the Standard Template Library (STL) and very rarely do low level stuff myself. Out of nowhere I remembered that a few days ago I was thinking about how I'd implement a size limited cache with a Least Recently Used (LRU) eviction policy. You need such a beast when you want to cache the result of the last computetions of some function without eating too much memory. The LRU policy is the simplest and probably easiest to implement. The solution I came up at that time (but never implemented) used a map and a doubly linked list to offer O(lgN) operations. So I thought: OK, let me implement this list by hand instead of using STL's list. it should be easy enough and the result would be some code that I sometimes wanted to have handy, but that I never implemented because the effort didn't seem to be worth it. Now I was in the mood of handling pointers :)

Without much comments, here is the code:

   1:/*
   2:    Implementation of a size limited cache. It is based on STL's map.
   3: */
   4:
   5:#include <cassert>
   6:#include <iostream>
   7:#include <map>
   8:
   9:using namespace std;
  10:
  11:template <typename K, typename V>
  12:class Cache
  13:{
  14:private:
  15:    // --- Internal types ---
  16:    struct Node;
  17:    typedef typename map<K, Node*>::iterator It;
  18:    struct Node {
  19:        It it;
  20:        V v;
  21:        Node *L, *R;
  22:    };
  23:    
  24:public:
  25:    // --- Exception for outside ---
  26:    struct NotFound {};
  27:    
  28:public:
  29:    // --- Construction / Destruction ---
  30:    Cache(int sz) : sz(sz) 
  31:    {
  32:        assert (sz > 0);
  33:        q = new Node;
  34:        q->L = q->R = q;
  35:    }
  36:    
  37:    ~Cache()
  38:    {
  39:        Node *d = q, *t;
  40:        do {
  41:            t = d->R;
  42:            delete d;
  43:            d = t;
  44:        } while (d != q);
  45:    }
  46:    
  47:    // --- Put / Get / Rem interface ---
  48:    void Put(const K& k, const V& v)
  49:    {
  50:        Node* t;
  51:        It it = m.find(k);
  52:        if (it == m.end()) {
  53:            if (int(m.size()) == sz) {
  54:                m.erase(q->R->it);
  55:                q->R = q->R->R;
  56:                delete q->R->L;
  57:                q->R->L = q;
  58:            }
  59:            t = new Node;
  60:            t->it = m.insert(m.begin(), make_pair(k, t)); t->v = v;
  61:        } else {
  62:            t = it->second;
  63:            t->L->R = t->R;
  64:            t->R->L = t->L;
  65:        }
  66:        t->L = q->L;
  67:        t->R = q;
  68:        t->L->R = t;
  69:        t->R->L = t;
  70:    }
  71:    
  72:    V Get(const K& k)
  73:    {
  74:        It it = m.find(k);
  75:        if (it == m.end()) throw NotFound();
  76:        Node* t = it->second;
  77:        t->L->R = t->R;
  78:        t->R->L = t->L;
  79:        t->L = q->L;
  80:        t->R = q;
  81:        t->L->R = t;
  82:        t->R->L = t;
  83:        return t->v;
  84:    }
  85:    
  86:    void Rem(const K& k)
  87:    {
  88:        It it = m.find(k);
  89:        if (it == m.end()) throw NotFound();
  90:        Node* t = it->second;
  91:        t->L->R = t->R;
  92:        t->R->L = t->L;
  93:        delete t;
  94:        m.erase(it);
  95:    }
  96:    
  97:private:
  98:    // --- Internal data ---
  99:    map<K, Node*> m;
 100:    int sz;
 101:    Node* q;
 102:};
 103:
 104:int main()
 105:{
 106:    Cache<char, int> cache(3);
 107:    cache.Put('a', 2);
 108:    cache.Put('b', 'b');
 109:    cache.Put('c', 10);
 110:    cache.Put('d', 20); // a should be removed
 111:    try { cout << cache.Get('c') << endl; } catch (...) {}
 112:    try { cout << cache.Get('a') << endl; } catch (...) {}
 113:    try { cout << cache.Get('b') << endl; } catch (...) {}
 114:    try { cout << cache.Get('d') << endl; } catch (...) {}
 115:    cache.Put('e', 100); // c should be removed
 116:    try { cout << cache.Get('a') << endl; } catch (...) {}
 117:    try { cout << cache.Get('b') << endl; } catch (...) {}
 118:    try { cout << cache.Get('c') << endl; } catch (...) {}
 119:    try { cout << cache.Get('d') << endl; } catch (...) {}
 120:    try { cout << cache.Get('e') << endl; } catch (...) {}
 121:}
 122:

Lars' contest

This is a long-term (one month) programming contest. Definitely worth a look. Your job is to find an answer as good as possible to two hard problems.

1. Find zero-one matrices that maximize (determinant / one_count2);
2. Find a linear zero-one system of equations with a single solution such that the closest root to 1 is as close as possible but not quite 1.

The job of finding the zero-one matrix that maximizes the determinant has been a research math problem for over 100 years.