Saturday, June 11, 2016

What kind of strings have Kolmogorov complexity equal to their length ?

I have been doing some reading about Information Theory and Kolmogorov complexity lately. I just thought about the following problem. Consider a language of all C programs (strings) which print themselves, formally Let $$L_{self} = \{ w\,\, | w\,\, \mbox{is a C program which prints itself} \}$$.
Following is an example of a C program which prints itself (on a side-note the trick behind the construction is to use the ASCII values for quotes and newlines). Clearly the Kolmogorov complexity of each w is less than or equal to the length of w -- given that we have constructed a program of length same as the string w.
I'm wondering if this is indeed optimal ? -- that is given a w, if its possible to produce a program to generate w and size strictly less than |w|.

char qt=34;
char nl=10;
char *str="char qt=34;%cchar nl=10;%cchar *str=%c%s%c;%cint main() {%c printf(str,nl,nl,qt,str,qt,nl,nl,nl,nl);%c}%c";
int main() {
 printf(str,nl,nl,qt,str,qt,nl,nl,nl,nl);
}

BTW, compile the program with
gcc -w 
to suppress warnings about header inclusion.

Saturday, March 07, 2015

Some Special Cases Of the Sequence Assembly Problem.

Sequence Assembly (SA) is a well studied problem in Bioinformatics. To quickly summarize the problem: let $\Sigma$ be some alphabet and $S =\{ s_1 , s_2 \ldots s_n\}$ be a set of strings from $\Sigma$. The goal of the SA problem is re-construct a string $s_a$ such that every $s_i \in S$ is a sub-string in $s_a$. Although its very tempting to reduce the SA problem to the Shortest Super String Problem (which happens to be $\mathcal{NP}$-hard), the repeats in genome makes the SA problem complex and ambiguous.

Any way one of my friends showed me this problem yesterday, which is a much simpler special case of the SA problem. Let $s_{\pi}$ be a some permutation (unknown) of $\Sigma$. Now given a set $S' = \{s'_1,s'_2\ldots s'_m\}$ of sub-sequences of $s_{\pi}$, we want to re-construct this un-known $s_{\pi}$. Just like the SA problem for a given $S'$ there could be several answers for $s_{\pi}$ -- in which case we want to find out the lexicographically smallest permutation. If $\Sigma_i |s'_i| = N$ the problem could be solved efficiently by reducing it to a topological sorting algorithm on a DAG. Following is a quick solution to this problem I wrote while watching the RSA vs PAK world cup cricket match -- I'm little disappointed that RSA lost despite a ferocious knock by AB. The construction of the DAG from $S'$ takes $O(Nlog(N)$ operations, however the cost of topological sorting to report the answer is still $O(N)$ operations.

// 03/06/2015: vamsik (at) engr.uconn.edu//
#include <cmath>
#include <cstdio>
#include <map>
#include <list>
#include <unordered_map>
#include <iostream>
#include <algorithm>
#include <cassert>
using namespace std;


struct CDirectedGraph{
    map<unsigned int, map<unsigned int, bool> > m_aEdges;
    unordered_map<unsigned int, unsigned int> m_aInDegree;
    
    
    inline CDirectedGraph() : m_aEdges(), m_aInDegree() {}
    
    inline void add_edge(unsigned int i, unsigned int j){
        if((m_aEdges[i]).find(j) == (m_aEdges[i]).end()){
            if(m_aInDegree.find(j) == m_aInDegree.end()){
                m_aInDegree[j] = 0;
            }
            m_aInDegree[j]++;
            (m_aEdges[i])[j] = true;
        }
    }
    
    inline void top_sort(){
        map<unsigned int, bool> snodes;
        unsigned int v;
        for(auto vItr=m_aEdges.begin(); vItr!=m_aEdges.end(); ++vItr){
            if(!m_aInDegree[vItr->first]){
                snodes[vItr->first] = true;
            }
        }
        
        while(!snodes.empty()){
            auto vItr = snodes.begin();
            v = vItr->first;
            auto eEnd = m_aEdges[v].end();
            auto eBeg = m_aEdges[v].begin();
            
            for(auto eItr=eBeg;  eItr!=eEnd; ++eItr){
                m_aInDegree[eItr->first]--;
                if(!m_aInDegree[eItr->first]){
                    snodes[eItr->first] = true;
                }
            }
            printf("%u ", v);
            snodes.erase(vItr);
        }
    }
};

int main() {
    unsigned int N;
    unsigned int K, pnode, cnode;
    scanf("%u", &N);
    CDirectedGraph diGraph;
    
    for(unsigned int i=0; i<N; i++){
        scanf("%u", &K);
        scanf("%u", &pnode);
        for(unsigned int j=0; j<K-1; j++){
            scanf("%u", &cnode);
            diGraph.add_edge(pnode, cnode);
            pnode = cnode;
        }
    }
    diGraph.top_sort();
    return 0;
}

Sunday, November 16, 2014

Pseudo Polynomial Dynamic Programming

Sub Set Sum and Integer Knapsack are few examples of combinatorial optimization problems with pseudo polynomial running time algorithms. Often we the implementation of these pseudo polynomial dynamic programming algorithms is done using a huge array to hold intermediate sub problems (e.g. $OPT[K][n]$ would indicate a presence of a subset in the first $n$ integer elements ($x_1,x_2\ldots x_n$) which sum to $K$). If $N$ is the bound on the subset value you were looking for then the space complexity would $\Theta(Nn)$ -- which is clearly exponential on the input size which is $O(n\log(N))$ bits. To build robust algorithms -- reduce the dependence on the value of $N$ -- for these pseudo polynomial algorithms, it would be efficient if we could exploit the underlying DAG (Directed Acyclic Graph) nature of any dynamic programming formulation. But that might be a little bit more code to traverse the sub-problems in a topological order. However we can make our life a little easy by using a hash table to keep track of the sub-problems, which might add to the worst case asymptotic runtime by a factor of $O(log(nN))$. On the other hand the average case of sub-problem lookup would be a constant and would run fairly efficiently in practice.

The following code is a quick illustration of using a hash tables to keep track (and also lookup) of sub-problems in pseudo polynomial dynamic programming algorithms. You can read the problem statement here . Notice the obvious generalization of the problem (scheduling with constant number of resources) from the way I have formulated the dynamic program.


// 11/15/2014: vamsik@engr.uconn.edu//
#include 
#include 
#include 
#include 
#include 
#include 
using namespace std;
typedef std::pair KeyType;
struct KeyTypeHasher{
    inline size_t operator()(const KeyType& a) const {
        return (a.first)*1729 + a.second;
    }
};
typedef std::unordered_map SubProbType;
typedef typename SubProbType::iterator SubProbItr;

bool IsFeasible(const std::vector& A, unsigned int G){
    SubProbType sub_problems_even, sub_problems_odd;
    SubProbItr itr;
    
    //init//
    if(A[0] <= G){
        sub_problems_even[KeyType(A[0], 0)] = true;
        sub_problems_even[KeyType(0, A[0])] = true;
    }
    
    for(size_t i=1; i < A.size(); i++){
        SubProbType& curr_problems = 
            (i%2) ? sub_problems_even : sub_problems_odd;
        SubProbType& next_problems =
            (i%2) ? sub_problems_odd : sub_problems_even;
        
        if(curr_problems.empty()){
            return false;
        }
        next_problems.clear();
        //create new sub-problems in the next level//
        for(itr = curr_problems.begin(); 
            itr != curr_problems.end(); ++itr){
            const KeyType &a = itr->first;
            
            if( (a.first + A[i]) <= G){
                next_problems[ KeyType((a.first+A[i]), a.second)] = true;
            }
            
            if( (a.second + A[i]) <= G){
                next_problems[ KeyType(a.first, (a.second + A[i]))] = true;
            }
        }
        curr_problems.clear();
    }
    
    return ((A.size())%2) ? !(sub_problems_even.empty()) : !(sub_problems_odd.empty());
}

int main() {
    /* Enter your code here. Read input from STDIN. Print output to STDOUT */
    size_t T,N;
    unsigned int G, a;
    scanf("%lu", &T);
    for(size_t i=0; i < T; i++){
       std::vector A;
       scanf("%lu %u",&N, &G);
       for(size_t i=0; i < N; i++){
           scanf("%u", &a);
           A.push_back(a);
       }
       printf("%s\n", IsFeasible(A, G) ? "YES" : "NO"); 
    }
    return 0;
}

Monday, August 25, 2014

Planarity Puzzle

I'm obsessed with the untangling of the planar graphs. I'm now able to untangle a 30-vertex planar graph in ~240sec see below (yellow one is the final solution):


Sunday, August 10, 2014

Sunday, January 05, 2014

Useful properties of the function $$f(x) = \frac{x}{x-k}$$

The function f(x) = x/(x-k) has some very useful properties -- especially when we restrict x to Integers. Let me start with a simple question about Maxima and Minima of f(x). Applying elementary Calculus may not be useful since f(x) ↦ +∞ as x ↦ k+ (goes to -∞ from the other side and is discontinuous). A quick plot on WolframAlpha reveals the vertical asymptotes [http://wolfr.am/JSXDdV].



Now lets move our attention to the case where x is restricted to Integers. First lets consider the function f(x) when x takes on integers strictly greater than k (i.e. x > k). Notice that f(x) is decreasing the maximum in this case occurs when x = k+1. This gives us the following useful fact:
$$
f(n) = \frac{n}{n-k} \leq k+1 \,\,\,\,\,\,\, n \in \{ x\, | \, x > k \vee x \in \cal{I} \}
$$

The above inequality is quite handy in your combinatorial analysis tool book. Let me give you a simple application of this inequality by proving that nlog(n) = O(log(n!)) (also recall that log(n!) = O(nlog(n)) which means nlog(n) = Θ(log(n!))).

$$
\begin{array}{lcl}
&& \frac{n}{n-k} \leq k+1 \,\,\,\,\,\, k \geq 1\\
&& \\
\rightarrow && n \leq (n-k)\times (k+1) \\
\rightarrow && n < (n-k) \times (k+k) = 2 (n-k)\times k
&& \\
&&\\
\rightarrow && n < 2(n-1)\times 1\\
\rightarrow && n < 2(n-2)\times 2\\
\rightarrow && n < 2(n-3)\times 3\\
&& \ldots \\
\rightarrow && n < 2(n-(n-1))\times (n-1) = 2(1\times (n-1)\\
&&\\
&&\mbox{multiplying all the above n-1 inequalities you get the following} \\
&&\\
\rightarrow&& n^{n-1} < \displaystyle\Pi_{k=1}^{n-1} 2(n-k)\times k = 2^n ((n-1)!)^2 \\
\rightarrow&& n^n < 2^n (n!)^2 \\
\rightarrow&& \log(n^n) < \log(2^n(n!)^2) = 2\log(n!) + n \\
\rightarrow&& n\log(n) < 2\log(n!) + n < 2\log(n!) + \log(n!) = 3\log(n!) \\
\rightarrow&& n\log(n) = O(\log(n!)) \,\, \Box\\
\end{array}
$$

Saturday, November 30, 2013

Getting around with integer overflows during Combinatorial counting

Many closed form expressions characterizing the number of combinatorial objects involve manipulation of raising and falling factorials (e.g. $$\displaystyle\Pi_{i=1}^{k} (n-i+1)$$). These quantities accumulate very fast we need to use extreme caution during computation. Lets take the example in which we want to compute $${n\choose k} = \frac{n\times (n-1)\ldots (n-k+1)}{k\times k-1\times \ldots 1}$$ which is the ratio of two falling factorials which is always an integer. A naive algorithm which computes numerator and denominator will easily overflow even though the actual value is well within the limits of an integer. If $$n_{max}$$ is the largest value represented by the underlying data-type, the naive algorithm would overflow when: $$n\geq \sqrt[k]{n_{max}}$$. One way to fix this issue is to maintain the ratio implicitly as two integers and knocking off the GCD when ever possible. This reduces the overflow during computation see below for a simple example. Thanks to my friend Chittu for discussing this over lunch.


  public static long computeNchooseK(int n, int k){
    assert n >=1 && k<=n : "Invalid n and k to compute";
    long num=1, denom=1, factor;
    do{
      factor = gcd(n,k);
      num *= (n--/factor); denom *= (k--/factor);
      factor = gcd(num, denom);
      num /= factor; denom /= factor;
    }while(k>=1);
    assert denom == 1 : "Something went wrong in computing NchooseK";
    return num;
  }

Saturday, September 21, 2013

$\sum_{k=1}^{n} H_k = (n+1)(H_{n+1} - 1)$

$$
\frac{1}{1-z} = 1 + z + z^2 + z^3 \dots \\
\int \frac{1}{1-z}dz = \frac{z}{1} + \frac{z^2}{2} + \frac{z^3}{3} \dots = \sum_{k=1}^{\infty} \frac{z}{k} z^k = -\ln(1-z)\\
\\
H_n = \sum_{k=1}^{n} \frac{1}{k} \hspace{0.5in} \mbox{By definition, also notice they are partial sums}\\
\rightarrow \mbox{O.G.F of $H_n$ is}\hspace{0.02in} \frac{1}{1-z}\times (-\ln(1-z)) = \frac{1}{1-z}\ln(\frac{1}{1-z}) \hspace{0.5in}\mbox{from our previous result on partial sums}\\
\rightarrow \frac{1}{1-z}\ln(\frac{1}{1-z}) = H_1\times z + H_2\times z^2 + H_3\times z^3 \ldots \hspace{0.5in} \mbox{(1)}\\
\rightarrow \frac{d}{dz}(\frac{1}{1-z}\ln(\frac{1}{1-z})) = \frac{1}{(1-z)^2}\ln(\frac{1}{1-z}) + \frac{1}{(1-z)^2} = H_1 + 2H_2\times z + 3H_3\times z^2 \ldots \hspace{0.5in}\mbox{(2)}\\
\mbox{Notice that}\hspace{0.02in} \frac{1}{(1-z)^2} = 1 + 2z + 3z^2 \dots \hspace{0.5in} \mbox{(3)}\\
\rightarrow\mbox{(2) - (3)} = \frac{1}{(1-z)^2}\ln(\frac{1}{(1-z)}) = (H_1-1) + (2H_2-2)\times z + (3H_3 -3) \times z^2 \ldots \\
\rightarrow \frac{1}{(1-z)^2}\ln(\frac{1}{(1-z)}) \hspace{0.02in} \mbox{is the O.G.F of sequence} \hspace{0.02in} (n+1)H_{n+1} - (n+1) = (n+1)(H_{n+1}-1) \hspace{0.5in} \mbox{(4)}\\
\\
\mbox{Applying the result of the partial sums to (1) the O.G.F of sequence}\hspace{0.02in} \sum_{k=1}^{n} H_n \hspace{0.02in} \mbox{is} \hspace{0.02in} \frac{1}{(1-z)}\times {\frac{1}{(1-z)}\ln(\frac{1}{1-z})} = \frac{1}{(1-z)^2}\ln(\frac{1}{(1-z)}) \hspace{0.5in} \mbox{(5)}\\

\mbox{From (4) and (5) the generating functions of both the sequences are same, hence the corresponding terms must be the same}





$$

Simple proof for the O.G.F of partial sums

Let $A(z)$ be the O.G.F of the sequence $S = \{a_0, a_1, a_2 \ldots \}$, then we can derive the O.G.F of the sequence $b_n = \sum_{i=0}^{n} a_i$ (i.e. $B = \{a_0, a_0+a_1, a_0+a_1+a_2 \ldots \}$) as follows.

Consider the sequence $S^1 = \{0 , a_0, a_1, a_2, a_3 \ldots \}$ then O.G.F corresponding to $S^1$ is $0\times z^0 + a_0 \times z + a_1 \times z^2 + a_2 \times z^3 \ldots = z\times A(z)$. Similarly the O.G.F corresponding to the sequence $S^2 = \{0, 0, a_0, a_1, a_3 \ldots \}$ is $z^2 \times A(z)$. Now notice that if we add corresponding terms in the sequences $S^1 , S^2, S^3 \ldots$ then this
results in the sequence $B$. Hence the generating function corresponding to $B$ is $A(z) + z\times A(z) + z^2\times A(z) + z^3\times A(z) \ldots = A(z)\times(1 + z + z^2 + z^3 \ldots) = \frac{A(z)}{1-z}$ $\Box$.

Sunday, June 09, 2013

An iterative method to compute the integer linear combination of Greatest Common Divisor

Euclid's iterative method to compute the GCD is very well known. You might have also heard that $GCD(a,b)$ can be expressed as integer linear combination of $a,b$ (i.e. $GCD(a,b) = \alpha\times a + \beta \times b \,\, \alpha,\beta \in \mathcal{I}$). The following is an iterative method -- I wrote on a flight between AUS -> DFW just for fun -- which computes this integer linear combination. Runs in $O(\log(\max\{a,b\}))$ time.


template
 static inline Unit gcd_integer_linear_combination(const Unit& a, const Unit& b, Unit& alpha, Unit& beta ){
  Unit divd = a, divs = b, remd=0, quot;
  Unit a_divd[2] = {1,0}; //linear combination of the divident//
  Unit a_divs[2] = {0,1}, temp[2] = {0,0}; //linear combination of the divisor//
  bool b_flip = false;
  if(divd < divs){
   divd = b; divs = a;
   b_flip = true;
  }
  //Invariant's: divd = a_divd[0]*a + a_divd[1]*b //
  //      divs = a_divs[0]*a + a_divs[1]*b //
  while(divs){
   remd = divd%divs; quot = (divd-remd)/divs;
   divd = divs; divs = remd;
   temp[0] = a_divs[0]; temp[1] = a_divs[1];
   a_divs[0] = a_divd[0] - (quot * a_divs[0]); 
   a_divs[1] = a_divd[1] - (quot * a_divs[1]);
   a_divd[0] = temp[0]; a_divd[1] = temp[1];
  }
  if(!b_flip){
   alpha = a_divd[0]; beta = a_divd[1];
  }else{
   alpha = a_divd[1]; beta = a_divd[0];
  }
  return divd;
 }

The following are some examples $GCD(a,b) = (\alpha)*a + (\beta)*b$.
PASSED: 2 = (1613)*42 + (-8)*8468 
PASSED: 1 = (-1684)*6335 + (1641)*6501 
PASSED: 5 = (344)*9170 + (-551)*5725 
PASSED: 1 = (-1949)*1479 + (308)*9359 
PASSED: 1 = (967)*6963 + (-1508)*4465 
PASSED: 2 = (858)*5706 + (-601)*8146 
PASSED: 6 = (181)*3282 + (-87)*6828 
PASSED: 2 = (121)*9962 + (-2450)*492 
PASSED: 1 = (-596)*2996 + (919)*1943 
PASSED: 1 = (-2348)*4828 + (2085)*5437 
PASSED: 1 = (283)*2392 + (-147)*4605 
PASSED: 1 = (-61)*3903 + (1546)*154 
PASSED: 1 = (122)*293 + (-15)*2383 
PASSED: 1 = (-1528)*7422 + (1301)*8717 
PASSED: 1 = (615)*9719 + (-604)*9896 
PASSED: 1 = (-304)*5448 + (959)*1727 
PASSED: 1 = (-139)*4772 + (431)*1539 
PASSED: 1 = (493)*1870 + (-93)*9913 
PASSED: 4 = (628)*5668 + (-565)*6300 
PASSED: 1 = (3461)*7036 + (-2461)*9895 
PASSED: 4 = (-60)*8704 + (137)*3812 
PASSED: 1 = (77)*1323 + (-305)*334 
PASSED: 3 = (662)*7674 + (-1089)*4665 
PASSED: 2 = (3)*5142 + (-2)*7712 
PASSED: 1 = (1840)*8254 + (-2211)*6869 
PASSED: 1 = (2027)*5548 + (-1471)*7645 
PASSED: 1 = (929)*2663 + (-897)*2758 
PASSED: 2 = (-301)*38 + (4)*2860 
PASSED: 2 = (-823)*8724 + (737)*9742 
PASSED: 1 = (-3)*7530 + (29)*779 
PASSED: 1 = (-1499)*2317 + (1144)*3036 
PASSED: 1 = (662)*2191 + (-787)*1843 
PASSED: 1 = (10)*289 + (-27)*107 
PASSED: 1 = (4289)*9041 + (-4336)*8943

Friday, October 05, 2012

Friday fun with a dynamic program.

Its been almost an year I wrote my personal blog. I missed so many fun things I used to do as student. Well I'm a pursuit of relentless improvement. Today I just wanted to write a dynamic program --  which I have missed for a while -- so I randomly pick one of the hardest DIV-1 problems and got the dynamic program working. This is a very interesting dynamic program -- a two level one. You can read the problem -- because I cannot' post here-- from this link http://community.topcoder.com/stat?c=problem_statement&pm=1996&rd=4710 . Following is my code which passed all the 88 system tests.

/*
 * Solution to  problem: http://community.topcoder.com/stat?c=problem_statement&pm=1996&rd=4710 (SRM-178, DIV-1, hard)
 *
 * Vamsi K. Kundeti 10/05/2012
 *
 **/
#include
#include
#include
#include
#include

using namespace std;
#define MAX_ROWS 50
#define MAX_STROKES 3000
class MinPaint {
   private:
      /*Optimal number of mismatches to color a row 'i' till 'j' with a
       *stoke ending black using 'k' strokes : 'm_aB[j%2][i][k]'
       */
      size_t m_aB[2][MAX_ROWS][MAX_STROKES];
      size_t m_aW[2][MAX_ROWS][MAX_STROKES];
      /*Optimal number of mismatches in the image upto row 'i' using 'k' 
       *strokes*/
      int m_aOPT[2][3001];
      //optimal value is stored in m_aB[0][i][j]
      vector m_aImg;

   public:
      MinPaint(): m_aB(), m_aW(), m_aImg() {}

      void initRow(size_t idx, size_t maxStrokes){
         for(size_t j=0; j<=maxStrokes; j++){
            m_aB[0][idx][j] = !j ? 1 : m_aImg[idx][0] == 'W';
            m_aW[0][idx][j] = !j ? 1 : m_aImg[idx][0] == 'B';
         }
      }
      size_t findMin(size_t a, size_t b){
         return a < b ? a : b;
      }
      void dynamicProgramPerRow(size_t idx, size_t maxStrokes){
         initRow(idx, maxStrokes);
         for(size_t i=1; i < m_aImg[idx].size(); i++){
            for(size_t j=0; j <= maxStrokes; j++){
               if(m_aImg[idx][i] == 'W'){
                  m_aW[i%2][idx][j] = !j ? i+1 : 
                     findMin(m_aW[(i-1)%2][idx][j], m_aB[(i-1)%2][idx][j-1]);
                  m_aB[i%2][idx][j] = !j ? i+1 : 
                     findMin(m_aB[(i-1)%2][idx][j]+1, m_aW[(i-1)%2][idx][j-1]+1);
               }else{
                  m_aB[i%2][idx][j] = !j ? i+1 : 
                     findMin(m_aB[(i-1)%2][idx][j], m_aW[(i-1)%2][idx][j-1]);
                  m_aW[i%2][idx][j] = !j ? i+1 : 
                     findMin(m_aW[(i-1)%2][idx][j]+1, m_aB[(i-1)%2][idx][j-1]+1);
               }
            }
         }
         //choose the min between m_aB and m_aW
         size_t nidx = m_aImg[idx].size()-1;
         for(size_t i=0; i<=maxStrokes; i++){
            m_aB[nidx%2][idx][i] = findMin(m_aB[nidx%2][idx][i], 
               m_aW[nidx%2][idx][i]);
         }
      }

      int leastBad(const vector &Img, const int &maxStrokes){
         m_aImg = Img;
          for(size_t i=0; i< Img.size(); i++){
            dynamicProgramPerRow(i, (size_t)maxStrokes);
         }

         //now find the optimal across rows//
         for(int i=0; i<=maxStrokes; i++){
            m_aOPT[0][i] = m_aB[(Img[0].size()-1)%2][0][i];
         }

         for(size_t j=1; j<Img.size(); j++){
            /*choose number of strokes from [1, maxStrokes-j] to paint row 'j'*/
            size_t rsize_idx = m_aImg[j].size()-1;
            for(int usedStrokes=0;  usedStrokes<=maxStrokes; usedStrokes++){
               int min = numeric_limits::max();
               int current_min;
               int min_k=0;
               for(int k=0;  k<=usedStrokes; k++){
                  //use 'k' strokes to color current row 'j' and 'usedStrokes-j' for the
                  //array below 'j'.
                  current_min = 
                     m_aOPT[(j-1)%2][usedStrokes-k] + (m_aB[(rsize_idx)%2][j][k]);
                  if(current_min < min){
                     min = current_min;
                     min_k = k;
                  }

               }
               m_aOPT[j%2][usedStrokes] = min;
            }
         }
         return m_aOPT[(Img.size()-1)%2][maxStrokes];
      }

      size_t OptPerRow(size_t idx, size_t maxStrokes){
         return m_aB[(m_aImg[idx].size()-1)%2][idx][maxStrokes-1];
      }
};
/*INPUT FILE FORMAT:
 * 
 *
 * X Y
 * string1
 * string2
 * ...
 * ...
 * ...
 * stringX
 *
 *X is the number of rows
 *Y is the maximum number of strokes
 **/
int main(int argc, char **argv){
   char buf[4096];
   vector input;
   size_t imgSize, maxStrokes;
   MinPaint X;

   while(scanf("%lu%lu", &imgSize, &maxStrokes)==2){
      input.clear();
      while(imgSize-- && scanf("%s", &buf) == 1){
         input.push_back(string(buf));
      }
      int opt = X.leastBad(input, maxStrokes);
      printf("%d\n", opt);
   }
}


Have fun -- because you can afford it only on a Friday.

Vamsi.

Saturday, May 14, 2011

[TECH] Enumerating set partitions with $DLX$

I have been a great fan of Dr. Knuth's paper about Dancing Links (Algorithm X) . In that paper Don suggests how the combinatorial search using backtracking can be improved in practice with the original idea from Hitotumatu and Noshita. One strange thing about Dancing Links is that its not Don's original idea, however Don made it popular by writing a paper about it. Don used the $DLX$ idea to solve the exact cover (matrix cover) problem and he gave several examples (e.g. Covering an rectilinear area with Pentominos) where he certain combinatorial problems to exact cover.

I have been waiting to find a problem where I get to apply $DLX$ without reducing it to exact cover problem. I'm really happy now that simultaneous hamming neighborhood problem (given a set $S=\{s_1,s_2\ldots s_n\},\, s_i\in \Sigma^{l}$ of strings compute the set $N_d(S) = \{ t | hamming(s_i,t)\leq d ,\forall s_i \in S\}$) which I have been working lately can be solved by applying $DLX$. As I mentioned in my earlier post the algorithm to compute $N_d(S)$ needs solve a set partition problem -- which further needs an algorithm to enumerate all the subsets of a given size $r$. In my previous post we have seen how to compute $n \choose r$ efficiently with backtracking. I'm now going use $DLX$ to solve the same problem. Notice that $DLX$ is an overkill to this simple problem, however it really makes sense if we look it from the perspective of the partitioning problem -- which is what we want to solve. $DLX$ just plays the role of UNDO operation see below.

/*Implementation of $n \choose r$ using $DLX$*/
enum_count_t EnumSubSetsCountDLX(size_t n, size_t r){
 enum_count_t n_choose_r=0;
 DLXState *dlx_sentinal = GetEnumSubSetsDLXState(n); /*Initialize the DLXState*/
 DLXState *sidx_dat[r+1], **state_idx=sidx_dat; /*item picked in the state*/
 /*total # of items picked in a given state*/
 size_t scidx_dat[r+1], *state_count_idx=scidx_dat; 
 size_t level;

 state_count_idx--; state_idx--; /*1-indexing*/
 level=1; state_count_idx[level]=0; state_idx[level] = dlx_sentinal;
 while(level){
  if(level == r+1){
   n_choose_r++;
   level--; 
   if(level) Dance(state_idx[level]);
  }else if(state_count_idx[level]+r-level > n-1){ /*state_idx[level]*/
   level--; 
   if(level) Dance(state_idx[level]);
  }else{
   state_idx[level] = state_idx[level]->right; 
   Break(state_idx[level]);
   state_count_idx[level]++; 
   /*track forward*/
   state_count_idx[level+1] = state_count_idx[level]; 
   state_idx[level+1] = state_idx[level];
   level++;
  }
 }
 CleanDLXState(dlx_sentinal); /*cleanup*/
 return n_choose_r;
}

/*Break and Dance operations*/
void Break(DLXState *state){
 state->left->right = state->right;
 state->right->left = state->left;
}
void Dance(DLXState *state){
 state->left->right = state;
 state->right->left = state;
}

/*Initializes a new DLX state required to enumerate $n \choose r$*/
EnumSubSetState* GetNewEnumSStateWithDLX(DLXState *dlx_sentinal, 
 size_t n, size_t r){

 EnumSubSetState *estate = malloc(sizeof(EnumSubSetState));
 assert(estate);
 estate->dlx_sentinal = dlx_sentinal;

 estate->state_count_idx = malloc(sizeof(size_t)*(r+1));
 estate->state_idx = malloc(sizeof(DLXState *)*(r+1));
 assert(estate->state_count_idx && estate->state_idx);

 (estate->state_count_idx)--; (estate->state_idx)--;

 estate->level=1; estate->state_count_idx[1]=0; 
 estate->state_idx[1] = estate->dlx_sentinal;
 estate->n = n; estate->r = r;
 return estate;

}



Next I'll post the most interesting part which is to enumerate the partitions using $DLX$.

Friday, May 06, 2011

[TECH] Enumerating subsets with backtracking and some interesting observations on the backtracking enumeration tree.

During my current research on enumerating all the common Hamming neighborhood of a set of strings. I needed an algorithm which can enumerate seven disjoint subsets from a given set. The fundamental building block in that is an algorithm which enumerates all the subsets of size $r$ from a given set of size $n$. I'm very much aware of the standard recursive algorithm which uses the fact ${n \choose r} = {n-1 \choose r} + {n-1 \choose r-1}$. However I wanted a much efficient version because I need to use this several times and cannot afford the constants in stack based implementations. So I wrote the following to backtracking based algorithm to enumerate the subsets. Once I did that I started observing some very interesting patterns in the number of leaves for each internal backtrack node at level $r-1$.

/*Back tracking*/ 
enum_count_t EnumSubSetsPrint(char *set, size_t n, size_t r){ 
    size_t state_idx_data[r+1], i, level; 
    size_t *state_idx = state_idx_data; 
    enum_count_t n_choose_r=0;  

    set--; 
    state_idx--; 
    for(i=1; i<=r+1; i++){ /*initialize*/ 
        state_idx[i] = 0; 
    } 

    level = 1; 
    while(level){ 
        if(level == r+1){
            /*back track*/ 
            for(i=1; i<=r; i++){ 
                printf("%c ", set[state_idx[i]]); 
            } 
            printf("\n"); 
            n_choose_r++; 
            level--; 
        }else if(state_idx[level]+r-level > n-1){ 
            level--; 
        }else{  
            /*track forward*/ 
            state_idx[level]++; /*pick a item*/ 
            state_idx[level+1] = state_idx[level]; 
            level++; 
        } 
    } 
    return n_choose_r; 
} 

By looking at these patterns (at the end of this paragraph) I came up with the following interesting identities for $n \choose 3$ and $n \choose 4$.
IDENTITY-1: ${n\choose 3} = \sum_{i=1}^{n-2} (n-i)(i-1)$
IDENTITY-2: ${n\choose 4} = \sum_{j=1}^{n-3}\sum_{i=1}^{n-2-j} (i)(j) $
I'm actually printing the number of leaves under an internal node at level $r-1$ in the backtracking tree. I'm observing the following pattern.

=================
8 3 0

6 5 4 3 2 1 
5 4 3 2 1 
4 3 2 1 
3 2 1 
2 1 
1 
(8 \choose 3) = 56 

Please enter n and r < n and print_option
 (e.g. 15 14 1 or 15 14 0)
9 3 0

7 6 5 4 3 2 1 
6 5 4 3 2 1 
5 4 3 2 1 
4 3 2 1 
3 2 1 
2 1 
1 
(9 \choose 3) = 84 

Please enter n and r < n and print_option
 (e.g. 15 14 1 or 15 14 0)
10 3 0

8 7 6 5 4 3 2 1 
7 6 5 4 3 2 1 
6 5 4 3 2 1 
5 4 3 2 1 
4 3 2 1 
3 2 1 
2 1 
1 
(10 \choose 3) = 120 
Please enter n and r < n and print_option
 (e.g. 15 14 1 or 15 14 0)
11 3 0

9 8 7 6 5 4 3 2 1 
8 7 6 5 4 3 2 1 
7 6 5 4 3 2 1 
6 5 4 3 2 1 
5 4 3 2 1 
4 3 2 1 
3 2 1 
2 1 
1 
===============================

The following are for $n \choose 4$ for various values of $n$.
===============================
8 4 0

5 4 3 2 1 
4 3 2 1 
3 2 1 
2 1 
1 
4 3 2 1 
3 2 1 
2 1 
1 
3 2 1 
2 1 
1 
2 1 
1 
1 
(8 \choose 4) = 70 

Please enter n and r < n and print_option
 (e.g. 15 14 1 or 15 14 0)
9 4 0

6 5 4 3 2 1 
5 4 3 2 1 
4 3 2 1 
3 2 1 
2 1 
1 
5 4 3 2 1 
4 3 2 1 
3 2 1 
2 1 
1 
4 3 2 1 
3 2 1 
2 1 
1 
3 2 1 
2 1 
1 
2 1 
1 
1 

===============================

Play around with this program Enumerate.c , Enumerate.h.

Finally I enabled $LaTeX$ setting for the blog.

Thursday, May 05, 2011

[THEORY] Is nature a giant guessing algorithm ?

I was reading Dr. Lipton's blog on guessing and was really impressed on how he relates the power of guessing with the question whether P=NP or P!=NP.

Looks like nature on the other had seems to be solving several combinatorial problems almost instantly (e.g. folding of the proteins to minimize its energy), does nature guess the solution to the problem ?. One way to guess a solution to a problem is to look inside the building blocks which formulated the problem. Dr. Lipton gives an example were we were supposed to guess the digits which can open a lock -- the lock we normally use in the GYM. It may be hard to guess the digits which can open a lock without any help. However if we can get an X-ray of the levers we can guess the solution to the problem more easily. So does this mean nature looks at the problems in a totally different perspective ? , is there something very intuitive to the nature which the humans fail to recognize ?.

Why did all the great mathematicians failed to find a polynomial solution 3-SAT ? -- or prove it does not exist. Is there something missing in the way we are looking at a digital circuit ? , can we design guessing algorithms to solve the NP-Complete problems ?. We don't know the answers for all these. But I'm sure that nature is a great guessing algorithm and it may be impossible for us to guess its guessing algorithm.

Tuesday, April 12, 2011

[TECH] Unbuffered message logging of interactive programs.

We generally use one of the following methods to log the output messages of programs which print a lot of messages.

  •  $ ./program >& message.log 
  •  $ ./program |& tee message.log 
Since stdout is line buffered we expect the same to hold when we re-direct stdout or any other stream. However this is not true if you keep monitoring the program status with tail -f message.log you will not notice any output from the program this is because the redirection is not line buffered. And in certain unfortunate cases when the running program aborts the redirection may not even flush the buffers and you find that message.log is empty and there is NO way to figure out what has happened.

One obvious way to fix this is to put fflush(stdout) after every printf statement you have, this is OK for small programs but for very large programs its really tedious. I came up with the following technique to solve this problem. This solution is generic.

[main.c]

int main(int argc, char **argv){

 ForceLogMessages("message.log");
 ....
 ... 
}

void ForceLogMessages(const char *mesg_file){
    int ret_setvbuf;
    /*close the stdout and re-open a text file*/
    fclose(stdout);
    stdout = fopen(name, "w");
    //make stdout unbuffered
    ret_setvbuf = setvbuf(stdout, (char *) NULL, _IOLBF, 0);   
    if(ret_setvbuf){
        fprintf(stderr, "UNABLE TO MAKE stdout UNBUFFERED %s", strerror(errno));
    }   
    assert(stdout);
}

Now if we run the program and we can monitor the status with tail -f message.log since we have made stdout unbuffered we can see a message in message.log every time a printf is executed in the original program and you can know what exactly is happening.

Thursday, March 24, 2011

[TECH] juggernaut-asm: An out-of-core sequence assembler

Development updates on the juggernaut-asm project.I'm now branching off the trunk to implement the feature which avoids using the UNIX sort. To checkout the main branch as follows. This is the version of the code I used to obtain the results on 4million reads which the in-memory algorithm of Velvet could not handle.

cvs -d:pserver:anonymous@juggernaut-asm.cvs.sourceforge.net:/cvsroot/juggernaut-asm co .

The code to replace the unix sort being developed in the branch name 'replace-unix-sort'.

cvs -d:ext:vamsi99@juggernaut-asm.cvs.sourceforge.net:/cvsroot/juggernaut-asm co -r replace-unix-sort .

Wednesday, March 23, 2011

[TECH] External R-Way Merge sorting.

http://lib-ex-sort.sourceforge.net/ is an external sorting program the following are the limitations I wish to remove before I leave in May. Most of them are engineering changes but are very useful for several algorithms based on external sorting.

  • Currently the size of each key is a constant. We need to remove this by adding a key_header for each key.
  • Use MMAP for integer sorting to avoid copy between user and kernel space.
  • Support for sorting when the data spans multiple files.

Friday, August 27, 2010

[TECH] Combinatorial counting under symmetries with Polya's Theorem

Recently I have started looking at the problem of enumerating Graphs. During my journey I was thrilled at the beauty of Polya's theorem. After applying the theorem on several examples I felt that I should write an article which explains the context and need for Polya's theorem.

This is the first blog I'm writing after taking a break over the summer. I had few things which had happened to me recently. I'm going to talk about them some other time as they are more philosophical.

Friday, April 09, 2010

[TECH] Super-fast tokenizing of sequence word patterns

C++ code colored by C++2HTML
#!/bin/perl 
#
# AHO-CORASICK-TOKENIZER on words, will also work on any delemiter.
# INPUT(1): Set of patterns (sequence of words) 'P'
# INPUT(2): Text 'T' 
# OUTPUT: report all the occuring patterns -- we don't even 
# want to know their locations. 
#
# NOTE: We don't want to worry much about the time required to 
# construct the failure function. All we need is to improve the
# scan speed.
#
# PATTERN FILE FORMAT:
# <annotation> $$$ pattern
#
# vamsik@engr.uconn.edu 04/08/2010
#

($#ARGV == 1) or die("USAGE: perl <name> pattern_file text_file");

open(PATTERN_FILE, $ARGV[0]) or die($!);

open(TEXT_FILE, $ARGV[1]) or die($!);
%AC_TREE;
%PATTERN_HASH;

$ANNOT_DELIM = ";;a;;";
$FAILED_DELIM = ";;f;;"; 
$HEIGHT_DELIM = ";;h;;"; #to move within the text

sub BuildKwordTree{

# root of the Aho-Corasick tree #
	my @params = @_;
	my $ac_root = $params[0];

	my $PFILE = $params[1];
	my ($line, $pword_list, $pattern, $curr_node, $pword);

	
	while(<$PFILE>){
		$line = $_; chomp($line);

		if($line =~ /([^\$]+)\$\$\$ (.*)/){
			$annotation = $1;

			$pattern = $2;
# split the pattern into words delimited by any #
# non-newline character. #
			@pword_list = split(/\s+/, $pattern);

			$curr_node = $ac_root; 
			foreach $pword (@pword_list){
				${$curr_node}{$pword} = 
					{} unless exists ${$curr_node}{$pword};

				$curr_node = ${$curr_node}{$pword};
			}
# put-the annotation into the last-node of the key-word tree 
			${$curr_node}{$ANNOT_DELIM} = $annotation;
		}
	}
}

#
# Build the failure function level-by level
#
#
sub BuildFailureFunction{
	my @params = @_;
	my $ac_root = $params[0];

	my ($k, @bfs_list, $keyword); 
#
# Incremental algorithm; For level-1 its the $ac_root itself  
#
	${$ac_root}{$FAILED_DELIM} = $ac_root;

	${$ac_root}{$HEIGHT_DELIM} = 0;
	foreach $k (keys %$ac_root){

		if(($k eq $FAILED_DELIM) or ($k eq $ANNOT_DELIM)
			or ($k eq $HEIGHT_DELIM)){

			next;
		}
		${${$ac_root}{$k}}{$FAILED_DELIM} =  $ac_root;

		${${$ac_root}{$k}}{$HEIGHT_DELIM} = 1;
		push(@bfs_list, ${$ac_root}{$k});
	}

	
	while($#bfs_list >= 0){
		$nref = shift(@bfs_list); #node

		foreach $keyword (keys %$nref){

			if(($keyword eq $FAILED_DELIM) or 
					($keyword eq $ANNOT_DELIM) or 
					($keyword eq $HEIGHT_DELIM)){

				next;
			}

			$nnref = ${$nref}{$keyword};

			$fpref = ${$nref}{$FAILED_DELIM}; 
			while(!(exists ${$fpref}{$keyword})

				and $fpref != $ac_root){
				$fpref = ${$fpref}{$FAILED_DELIM};
			}

			if(exists ${$fpref}{$keyword}){
				${$nnref}{$FAILED_DELIM} = ${$fpref}{$keyword};
			}else{

				${$nnref}{$FAILED_DELIM} = $ac_root; 
			}
			${$nnref}{$HEIGHT_DELIM} = ${$nref}{$HEIGHT_DELIM}+1;

			push(@bfs_list, $nnref);
		}
	}
}
#
# Reads word delimited text from the stream or a file
#
sub MatchWordDelimText{
	my @params = @_;

	my $ac_root = $params[0];
	my $stream = $params[1];

	my ($curr_state, $annot, $sword, $c, $height);

	my $line;

	while(<$stream>){
		$line = $_; chomp($line);

		$curr_state = $ac_root;
		@swords = split(/\s+/, $_);

		$c = 0; 
		while($c <= $#swords){

			$sword = $swords[$c];
# Change the state of the automaton by reading the input 
#			print "I:$sword C:$c 1:$curr_state ";
			if(exists ${$curr_state}{$sword}){

				$curr_state = ${$curr_state}{$sword};
				$c++;
			}else{

				if($curr_state == $ac_root){
					$c++;
				}
				$curr_state = ${$curr_state}{$FAILED_DELIM};
			}

#			print "2:$curr_state \n";
			if(exists ${$curr_state}{$ANNOT_DELIM}){
				$annot = ${$curr_state}{$ANNOT_DELIM};

				print "$annot \n";
			}
		}
	}
}

sub StressTestMatchWordDelimText{
	my @params = @_;

	my $ac_root = $params[0];
	my $stream = $params[1];

	my ($curr_state, $annot, $sword, $c, $height);

	my $line;
	my ($stress_me, $ss);
	$line = <$stream>;

	print "ABSTRACT: $line\n";
	chomp($line);
	$stress_me = 1000000;

	for($s=0; $s< $stress_me; $s++){
		if($s%100000 == 0){

			print "iteration $s \n";
		}
		$curr_state = $ac_root;
		@swords = split(/\s+/, $_);

		$c = 0; 
		while($c <= $#swords){

			$sword = $swords[$c];
# Change the state of the automaton by reading the input 
#			print "I:$sword C:$c 1:$curr_state ";
			if(exists ${$curr_state}{$sword}){

				$curr_state = ${$curr_state}{$sword};
				$c++;
			}else{

				if($curr_state == $ac_root){
					$c++;
				}
				$curr_state = ${$curr_state}{$FAILED_DELIM};
			}

#			print "2:$curr_state \n";
			if(exists ${$curr_state}{$ANNOT_DELIM}){
				$annot = ${$curr_state}{$ANNOT_DELIM};

#				print "$annot \n";
			}
		}
	}
}

#
# Just to see how the tree looks -- for viewing pleasure :)
#
sub PrintKwordTree{
	my @params = @_;

	my $ac_root = $params[0]; 
	my $level = $params[1];

	my ($k, $aref, $annot, $j, $nref, $fref, $height);

	foreach $k (keys %$ac_root){
		if(($k eq $ANNOT_DELIM) 
			or ($k eq $FAILED_DELIM) or ($k eq $HEIGHT_DELIM)){

			next;
		}
		for($j=0; $j<$level; $j++){

			print " ";
		}
		$fref = ${${$ac_root}{$k}}{$FAILED_DELIM};

		$nref = ${$ac_root}{$k};
		$height = ${$nref}{$HEIGHT_DELIM};

		print "-$k-> N=$nref F=$fref h=$height\n";
		PrintKwordTree(${$ac_root}{$k}, $level+1);
	}

	

	if(exists ${$ac_root}{$ANNOT_DELIM}){
		$annot = ${$ac_root}{$ANNOT_DELIM};

		for($j=0; $j<$level; $j++){
			print " ";
		}

		print "; $annot ;\n\n";
		return;
	}
}
#
# A simple unit-test
#
sub main{
	BuildKwordTree(\%AC_TREE, PATTERN_FILE);

#	PrintKwordTree(\%AC_TREE, 0);
	BuildFailureFunction(\%AC_TREE);
#	PrintKwordTree(\%AC_TREE, 0);
# Now build the failure function.
	print "please enter the words\n";
#	MatchWordDelimText(\%AC_TREE, STDIN);
	StressTestMatchWordDelimText(\%AC_TREE, TEXT_FILE);

	print "\n....THANK YOU......\n";
}
main();