Silent rants of my life.........

[TECH] Enumerating set partitions with $DLX$

2011-05-14T11:06:00.000-07:00

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$.

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

2011-05-06T18:52:00.000-07:00

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.

[THEORY] Is nature a giant guessing algorithm ?

2011-05-05T16:46:00.000-07:00

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.

[TECH] Unbuffered message logging of interactive programs.

2011-04-12T12:46:00.000-07:00

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.

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

2011-03-24T18:59:00.000-07:00

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 .

[TECH] External R-Way Merge sorting.

2011-03-23T13:42:00.000-07:00

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.

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

2010-08-27T16:19:00.000-07:00

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.

[TECH] Super-fast tokenizing of sequence word patterns

2010-04-09T01:49:00.000-07:00

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

[PHIL] My best lines in "The zen and art of motor cycle maintenance"

2010-03-11T14:08:00.000-08:00

We often read a lot. On the other hand we also tend to forget most of what we read. However the essence/flavor/curx what ever you might call may still persist. Some time if were asked to recall about something -- say a novel for example. Then some very distinctively marked paragraphs would normally on top of your head. So whenever I think about "The zen and the art of motorcycle maintenance" the following few paragraphs are always on my mind. You can read the full book here

In time...six months; five years, perhaps...a change could easily begin to take place. He would become less and less satisfied with a kind of dumb, day-to-day shopwork. His creative intelligence, stifled by too much theory and too many grades in college, would now become reawakened by the boredom of the shop. Thousands of hours of frustrating mechanical problems would have made him more interested in machine design. He would like to design machinery himself. He’d think he could do a better job. He would try modifying a few engines, meet with success, look for more success, but feel blocked because he didn’t have the theoretical information. He would discover that when before he felt stupid because of his lack of interest in theoretical information, he’d now find a brand of theoretical information which he’d have a lot of respect for, namely, mechanical engineering.

So he would come back to our degreeless and gradeless school, but with a difference. He’d no longer be a grade-motivated person. He’d be a knowledge-motivated person. He would need no external pushing to learn. His push would come from inside. He’d be a free man. He wouldn’t need a lot of discipline to shape him up. In fact, if the instructors assigned him were slacking on the job he would be likely to shape them up by asking rude questions. He’d be there to learn something, would be paying to learn something and they’d better come up with it.

Motivation of this sort, once it catches hold, is a ferocious force, and in the gradeless, degreeless institution where our student would find himself, he wouldn’t stop with rote engineering information. Physics and mathematics were going to come within his sphere of interest because he’d see he needed them. Metallurgy and electrical engineering would come up for attention. And, in the process of intellectual maturing that these abstract studies gave him, he would he likely to branch out into other theoretical areas that weren’t directly related to machines but had become a part of a newer larger goal. This larger goal wouldn’t be the imitation of education in Universities today, glossed over and concealed by grades and degrees that give the appearance of something happening when, in fact, almost nothing is going on. It would be the real thing.

[TECH] Tutte matrix and Sufficiency of Perfect Matching.

2010-03-10T23:07:00.000-08:00

I know that this is a standard results. But a hands on proof is really useful because it makes things more clear. I knew the fact that any permutation could be broken into disjoint cycles. However I did not know how the sign of the permutation changes if we swap the elements with in a disjoint cycle. It turns out that some of the fundamental properties of the signs of the permutation can be used to prove the sufficiency of a perfect matching in graph. A Tutte Matrix is defined as follows.

Given a Tutte matrix we now prove the sufficiency of perfect matching

Determinant of Tutte matrix is as follows . We know that every permutation can be expressed as disjoint cycles. For example permutation can be expressed as two cycles . Also is the number of swaps we need do required to convert to the permutation . We observe the following important properties which help solve the problem.

PROPERTY-1 If is a cycle in then if we replace with reverse cycle then remains the same - because we did not change the number of swaps required but just changed the order of swaps. For instance, let then cycle now we replace it with . This results in .
PROPERTY-2 Let is a cycle of odd length in and let be the permutation by replacing with its reverse cycle . Then the corresponding terms in the determinant for and are and . For example in the previous example so the corresponding term in the determinant for is . On the other hand the reverse cycle has the following term in the determinant.
PROPERTY-3 If all the cycles in the permutation are even and corresponding term in the determinant (i.e. ) is non-zero then its the set is a valid perfect match in the graph. This is because by the definition of Tutte matrix represents an edge in the graph and the permutation is giving a subset of edges since the product is non-zero, on the other than this sub-set of edges is a match.

"IF has a perfect match THEN has at least one non-zero term"

Let be the perfect match in then we construct a permutation where , this permutation on the other hand has even cycles - because we created a cycle for every pair. Now we apply PROPERTY-3 and the term corresponding to in is non-zero.

IF THEN has no perfect matching

Now consider the case when does not have no perfect matching. Then by PROPERTY-3 all the terms corresponding permutations with even cycles must be zero - otherwise our assumption leads to contradiction. Now consider the permutations with at least one odd cycle, let be that permutation and be the odd cycle. Now we replace with its reverse cycle and create the permutation . By PROPERTY-1 . Now we exploit the skew-symmetry of the Tutte matrix. By PROPERTY-2 the terms corresponding to , in the determinant are as follows. , . But by skew-symmetry , on the other hand since the cycle length is odd this gives that the actual terms corresponding to and are and . This means all the terms corresponding to permutations with at least one odd cycle can be grouped and pairs such as and be cancelled. So the only contribution for the comes from permutations will all even cycles - in this case we have none, so .

[LIFE] Keeping the PACE of productivity is very hard.

2010-03-09T15:31:00.000-08:00

This is the first of the posts I'll start writing some thing different from Technology. Spring break started last Friday (03/05/2010). I'm back to the out-of-core assembly project. Well I had this idea of out-of-core de Bruijn graph construction for a long time and in fact we did demonstrate the efficiency in parallel context. However the main focus has been to format the bi-directed de Bruijn graph in the format required by Velvet. I started of this quest some time mid October and finally I had a complete algorithm which generates the 'Graph' file in Velvet format. I was very happy last night when I got that but soon realized that there is a lot of difference in the number of CONTIGS. Went home with full of despair that I'm totally stupid. Today I came back and started it again. Productivity is something I strive for. KARMA is something I really believe in. However it seems to me that keeping up with the PACE of productivity is something very hard, the harder you try to focus the harder it might become. I strongly feel that continuous non-stop productivity is only possible when it becomes a HABIT rather than MOTIVATION. It might start off initially as MOTIVATION, but cultivating it as a HABIT is something on the lines of KARMA YOGA. Finally a true KARMA YOGI seems to be a person for whom continuous productivity is a HABIT rather than anything else.

[TECH] An efficient External sorting API

2010-02-01T19:09:00.000-08:00

Sorting is one of the most fundamental operation on data. Very efficient algorithms exists to perform this operation. Most of the programming languages are shipped with some sorting API with them. Most of them are in-memory sorting algorithms. However I find all these APIs are so bloated which are unfortunately by produces of abusing and overusing the object oriented concepts. On the other hand the bloated code may run fine as long as the input size is bounded by some constant -- which is most of the case for many applications. Any way I guess I'm getting into a little off-topic but if you are interested see the flames between C++ and C hackers on the Linux kernel mailing list here

On the other hand most of these programming APIs lack a solid External sorting algorithms. External sorting is * THE MOST * fundamental operation especially when you want to build algorithm on monster and massive datasets. You can try my new -- well old but wrapped in a new API, External sort please get it here . I'll try to post some examples how to use if when I get more time.

Rounding An Integer To The Next Maximal Mutliple Of A Given Radix Power

2010-01-30T13:12:00.000-08:00

Let be the modulo- representation of an integer , where each is a symbol/digit corresponding to values . Often we are encountered with problems where we need to find smallest integer such that (i.e. divides without any reminder), where . Some of the very common applications include rounding the number of bits required to represent a data structure to the nearest byte (i.e. power of ). Before we see how to address this problem its worth while to understand the following interesting property of modulo representation. Given a modulo- representation of we can get corresponding modulo- representation of by replacing every group of digits (in modulo- representation), by the corresponding digit (i.e. ) in modulo-. For instance if we would like to convert an integer in binary (i.e. ) representation to hexa-decimal representation (). We start from left to right and replace every bits with the corresponding digit in hexa-decimal system. For example if we see we will replace it with digit , by and so on. So the hexa-decimal system is providing us with a one-one function for every bit string, in fact we can use any one-one function here when we move for modulo representation to modulo . However the one-one function has become a standard for module system.

So coming back to our original problem given the modulo- representation , we would like to round to , where is the smallest multiple of such that . To accomplish this task we need to examine the first (from right) digits in the modulo- representation of . In fact the value of in those bits is the reminder we get when we divide by , so if is the value in those bits then . When we can elegantly use the bit-wise operators to accomplish this. So if some one gives an integer and ask to find a smallest which is a multiple of then we use the following C-statement to accomplish this . Where are the standard bit-wise operators in C.

[TECH] An interesting problem on the mailing list [gdb@gnu.org]

2009-12-24T17:53:00.000-08:00

Some one asked the following question of GDB mailing list.

Hi, 

  I was wondering how come the call to "call cos(.2)" returns the wrong
answer in:

#include "math.h"
double    (*fcn)(double);
main() 
{
  double t;
  fcn = cos;
  t=cos(.2);
  t=cos(-.2);
}

(gdb) call cos(.2)
$16 = 104
(gdb) call fcn(.2)
$17 = 0.98006657784124163

Thanks,
-- 
View this message in context:  http://old.nabble.com/Why-%22call-cos%28.2%29%22-returns-garbage--tp26909516p26909516.html 
Sent from the Gnu - gdb - General mailing list archive at Nabble.com.

Here is my analysis.

Hello,

Looks like 'gdb' is missing debug information about the function 'cos' , 
and the debugger implicitly assumes all the functions return an 'int' ,
 if the function is not compiled in debug mode.

I derive my explanation from the following test case.

1. Create a function 'double test_cos(double)' and compile it in optimized mode with -O3 using gcc
 "$gcc -c -O3 test1.c"

2. Complie the file containing the 'main' function in debug mode
"$gcc -g test.c test1.o -lm"

============[gdb a.out]====================

(gdb) list 0
1       #include "math.h"
2       #include "test1.h"
3       double    (*fcn)(double);
4       main()
5       {
6                double t;
7                 fcn = cos;
8                  t=cos(.2);
9                   t=cos(-.2);
10                      test_cos(.2);
(gdb) call fcn(0.2)
$4 = 0.98006657784124163
(gdb) call cos(0.2)
$5 = -1073792992
(gdb) call test_cos(0.2)
$6 = -1073792992
=====================================

============[test1.h]=================
#ifndef _test_1_h
#define _test_1_h
double test_cos(double);
#endif
====================================

============[test1.c]==================
#include "math.h"
#include "test1.h"
double test_cos(double a){
        return cos(a);
}
====================================

============[test.c]==================
#include "math.h"
#include "test1.h"

double    (*fcn)(double);
main()
{
         double t;
          fcn = cos;
           t=cos(.2);
            t=cos(-.2);
                test_cos(.2);
}
===================================

May be the developers of gdb can you more details.

also
=======================================
(gdb) call cos
$6 = {text variable, no debug info} 0xc325d0 {cos}
=======================================

[TECH] Manipulating DNA sequences with bits

2009-12-23T19:33:00.000-08:00

The alphabet of DNA sequences consists of 4 alphabets {A,T,G,C} and they exists in compliments -- A-T, G-C. Often in sequence assembly algorithms we need to squeeze these into bits for saving space since we have too many of them. For example a typical input to a sequence assembler consists of at least 10 million reads (string of fixed length, say 21). In my last post I briefly described the Sequence Assembly (SA) problem. As I pointed one of the basic operations on these reads would be to obtain a reverse compliment of a given read. For example our read is ATTACCAG , then its reverse compliment would be CTGGTAAT -- reverse the original sequence replace every symbol by its compliment (A-T, G-C). The following routines may be useful when we squeeze these symbols into bits using only 2bits per symbol. While using the bit-wise operators I realized that == operator has more precedence than & (bit-wise or). So if you a&b==b then that will always be true, so we need to be more careful use (a&b)==b.

#define MASK_A 3
#define MASK_T 0
#define MASK_G 1
#define MASK_C 2
/*Lets fix the maximum size of the read as 32 base pairs*/
typedef unsigned long long kmer_t;


/*Takes a read of length K and returns its bit represenation*/
kmer_t CreateBKmer(char *kmer, unsigned long K){
    unsigned long i;
    unsigned char bitc;
    /*A binary k-mer*/
    kmer_t bkmer = (kmer_t) 0; 
    for(i=0; i<K; i++){
        bkmer <<= 2; 
        bkmer |= Char2Bit(kmer[i]);
    }
    return bkmer;
}

/*Return the reverse compliment of the kmer_a*/
kmer_t ReverseCompliment(kmer_t a, unsigned int K){
    kmer_t rc=0;
    unsigned int i=0;
    for(i=0; i<K; i++){
        rc <<=2;
        rc |= (3-(a&3));
        a >>=2;
    }
    return rc;
}
/*Print the ASCII equivalent of the underlying bits*/
void PrintKmer(kmer_t a, unsigned int K, FILE *ptr){
    char ntide;
    kmer_t mask_a = MASK_A;
    kmer_t mask_t = MASK_T;
    kmer_t mask_g = MASK_G;
    kmer_t mask_c = MASK_C;

    mask_a <<= 2*(K-1); mask_t <<= 2*(K-1); 
    mask_g <<= 2*(K-1); mask_c <<= 2*(K-1);

    do{
        if((a&mask_a) == mask_a){
            ntide = 'A';
        }else if((a&mask_c) == mask_c){
            ntide = 'C';
        }else if((a&mask_g) == mask_g){
            ntide = 'G';
        }else{ /*This should be the last because T=0*/
            ntide = 'T';
        }
        fprintf(ptr, "%c", ntide);
        mask_a >>=2; mask_t >>=2;
        mask_g >>=2; mask_c >>=2;
    }while(mask_a); /*MASK_A = 3 */
    fprintf(ptr, "\n");
}

inline char Bit2Char(unsigned char a){
    switch(a){
        case MASK_A:
            return 'A';
        case MASK_G: 
            return 'G'; 
        case MASK_C:
            return 'C';
        case MASK_T:
            return 'T';
    }
}

inline unsigned char Char2Bit(char c){
    switch(c){
        case 'A':
            return MASK_A;
        case 'T':
            return MASK_T; 
        case 'G':
            return MASK_G;
        case 'C':
            return MASK_C; 
    }
}

[TECH] Loops in Bi-directed De Bruijn Graphs

2009-12-18T20:53:00.000-08:00

De Bruijn graphs recently found applications in Sequence Assembly(). problem is close to Shortest Common super String() problem, however there are some fundamental differences. Given a set of input strings strings the problem outputs a string such that every is a sub-string in and . On the other hand the input to the problem is a set of strings, every is a randomly chosen sub-string of a bigger string - called the genome. Unfortunately we do not know what is, so given the problem asks to find the such that we can explain every . Adding to the complexity of problem may contain several repeating regions and the direct appliction of to solve problem would result in collapsing the repeating regions in - so we cannot directly reduce the problem to . Also the string is special string coming from the so its double stranded and each may be originating from the forward or reverse compliments of .

A De Bruijn graph on a set of strings is a graph whose vertices correspond to the unique mers (a string of length ) from all the and whose edges correspond to a some mer in some . Also to consider the double stranded-ness we also include a reverse compliment of every mer into the graph . In the following simple example I show how a loop (i.e. both the forward and reverse compliments overlapping each other by ) can originate into . This interesting situation seems to be referred to as hairpin-loop.

[TECH] Simplified Proof For The Application Of Freivalds' Technique to Verify Matrix Multiplication

2009-12-12T20:52:00.000-08:00

We first give a simple and alternative proof for Theorem- in [Randomized Algorithms, Motwani R. and Raghavan P., pp 162--163]. Later in Theorem 2 we show that the assumption on the uniformness is not necessary.

Proof. Let be a matrix and be the column vectors of . Then . This means that multiplying a vector with a matrix is linear combination of the columns, the coefficient is the component of . Since is a boolean and acts as an indicator variable on the selection of column . So if is chosen from a uniform distribution .

Now let and be the column vectors of , similarly let be the column vectors of . Let , clearly since . Then since . Intuitively this means we select our random vector such that for all , such a selection will always ensure even though .

Proof. Continuing with the proof of Theorem- 1 , .

[TECH] Parallel sorting of binary keys

2009-12-04T22:37:00.000-08:00

I have just released the beta of parallel sorting using MPI. Feel free to use it visit the web page http://par-bin-sort.sourceforge.net/usage.html from more details.

[PHIL] My original quote "The more you learn the more stupid you feel"

2009-11-11T13:48:00.000-08:00

I was drinking coffee today and met a friend who was asking me how I was doing my research. I inadvertently replied her "The more I learn the more stupid I feel". Then I came back to my office and was proving a theorem that the cardinality of the optimal edge cover |C| is exactly |V|-|M|. After that I recalled the statement I made, and restated it as follows "IF you are more learned THEN more stupid you feel". From this its clear that learning more or wiseness is a subset of feeling stupid. So if you never feel stupid and always feel that you the genius who landed on this earth then its sufficient that your complexity class is far away from the complexity class of wiseness.

[TECH] Border Length Minimization Solver

2009-09-20T13:05:00.000-07:00

BLM (Border Length Minimization ) is a combinatorial optimization problem which belongs to the same class the QAP (Quadratic Assignment Problem). We now have a new result which indeed proves that BLM problem is NP-hard on any simple connected grid. We are developing a BLM solver which can give good results in practice. I'm now using Launchpad for all the development of this project. The latest commits can be kept track at https://code.launchpad.net/~vamsi-krishnak/blm-solve/trunk .

Algebraic methods for Combinatorial counting.

2009-09-01T21:50:00.000-07:00

I started reading Knuth's Volume-4 series - you can get it here http://www.cs.utsa.edu/~wagner/knuth/. The chapter Introduction to Combinatorial Searching starts off with the analysis of some of the fundamental combinatorial objects such as Langford pairs. I came across an interesting exercise which asks to prove a statement on how to count these Langford pairs. I infact spend some time on this problem and finally the solution seemed extremely elegant. Before that let me quicky define what a Langford pair is.

Langford Pair: Given a multi set of integers - every integer is repeated exactly twice. A Langford pair is a permutation of , with a property that there are exactly elements between and its copy, . (e.g. for , we can see that and its copy has exactly element between them. Also observe the property of and ) Having defined what a Langford pair is we now ask the question how can we count. The exercise I was talking about initially asks to prove the following.

Let denote all the possible Langford pairs which can be formed from set .
Let
Let
Show that

Looking at the problem first, I was a little unclear what is the intuition behind defining such a function (i.e. ). After careful thought and analysis I understood that the function was defined so cleverly that the total number of Langford pairs can be found algebraically. This quickly reminded me about the generating functions , where we use an algebraic expression to count possible partitions of a natural number. Let me first explain the intuition behind that algebraic formulation of function .

Consider a Langford permutation , now consider the positions of each integer occurs for the first time from the left of the permutation. By the definition of the Langford pair if an integer appears at position in the permutation, then none of the integers except can occur at position - which is towards the right. Also note that the integer can only occur at positions . Let , for example, in any Langford permutation of , the integer can only occur at positions . Similarly integer can occur at positions . Now this explains the intuition behind in the definition of the function .
Let be a function - we call it a position map. Which maps every integer to a index/position in the permutation . Let and . We now see that a position map can generate a Langford permutation only if . In other words the position mapping function should map the integers in such a way that, when we introduce the corresponding pairs - at positions , we should not see any duplicates. Now all possible function mappings of is described algebraically as follows.

Now if we observe carefully every term in the expansion of the above algebraic expression represents every possible (see above). More specifically every term in the above expansion is of the form and . Now we are looking for all the possible terms which have , every such term represents a Langford pair.
We now have seen how a clever algebraic expression can help us in counting difficult combinatorial objects. But we are still not done - why? . Because we need to filter out all the terms which don't satisfy . How do we do that ?, this is where the trick of odd , even functions come into play. Now by making we have a chance to wipe-out all terms which contain odd powers for some . This is done by multiplying the expression by by doing this every term in will have at least one with odd power except the Langford term which will be . Thus the summation on knocks out all the terms except the Langford terms as follows.

Vamsi Kundeti 2009-09-02

[TECH] Maintaining a Load Factor (α) in Dynamic Data structures with constant (Θ(1)) amortized cost per operation.

2009-03-10T18:43:00.000-07:00

Accounting Method is an elegant technique for amortized analysis of the data structures. Most of us should have encountered Dynamic data structures (which need to expand and contract to maintain some invariant) at some point of time. A stack implementation is one example of a dynamic data structure. As we perform pop operation on the stack eventually it become full (over-flow), at this stage we might choose to increase (may be by doubling or some other heuristic) the size of the stack. We also want to handle the under-flows i.e. if the number of elements in the stack is too small we might want to relocate the elements into a new stack ( with smaller size ). In general the dynamic data structures have a parameter called the load factor and our aim is to maintain that. One simple way to maintain this is to handle under-flows and over-flows on demand. Unfortunately this will not give us a amortized cost per operation on the stack. However we can still get amortized cost by avoiding this on demand expansion and compaction.

Let denote the number of elements present on the stack after operations (including push and pop). Let denote the allocated size of the stack after operations. We Define a potential energy function as follows.

Clearly . The amortized cost of the push operation which does not require an expansion of the stack is as follows.

The amortized cost of the push operation which requires doubling of the array is as follows.

The amortized cost of the pop

From the above we can see that that amortized cost of all the operations is when we use the doubling heuristic for the stack.

(b)

Let the load factor and let . We use the following heuristic to handle expansion and compaction of the array.

"When (i.e. the array is full) increase the size of the array by times. When (i.e. the array about to under flow) reduce the size of the array by half ()".

We define the following potential function and finally show that the amortized complexity of each of the operations push, pop is .

We now prove that the amortized cost of each operation is bounded by a constant. From the previous problem (since the potential function is same as the previous problem when ) we can observe that push operation when (expansion) has an amortized cost of . We can also observe that whenever expansion does not occur the amortized cost of the push operation is . So we will concentrate on the other possible situations. We use the notation to indicate the load factor of the array after operations on the array.

Amortized cost of push when

Amortized cost of push when ,

Amortized cost of pop when

Amortized cost of pop when and

Amortized cost of pop when with compaction

If the operation is a pop and we find that the load factor then we apply compaction. Compaction reduces the size of the allocated array by . So after compaction . The amortized cost of the pop and compaction is as follows.

We have proved that the amortized cost of both push and pop with all possible values of the load factor to be .

[TECH] Exponential function (e-μt) is the only unique solution satifying R(x+t) = R(x)R(t)(semi-group property).

2009-03-10T17:31:00.000-07:00

Exponential function keeps popping up every where stochastic analysis. However it has a crucial property which is often exploited during the analysis. The property I'm referring to is R(X+Y) = R(Y)R(Y) (e.g. X and Y) could be random variables). Now I prove a strong statement " R(X+Y) = R(X)R(Y) If and Only If R(x) = e^-μt where μ is some constant. The reverse direction (<==) is trivial because we can easily verify the fact by substituting R(x) with e^-μt. In the following we prove that other direction (==>).

One may wonder how about function a^x which also satisfies the semi-group property by inspection. However the μ takes care of that observe that a^x = e^ln(a)x.

[TECH] Simple comparision of real numbers

2008-12-25T03:03:00.000-08:00

We know that the fixed precision representation of real numbers is done by mapping them into integer space I from real space R. However the underlying integer representation of the 32-bit floats and 64-bit doubles is a sign magnitude representation (that is the reason why we see -0.00 during output of some numerical algorithms). People seem to be adding a constant to convert the underlying sign magnitude form to 2's compliment. However following is a more logical conversion (since we know that 2's compliment is 1's compliment +1 , the motivation of 2's compliment is to avoid the -0 as we know). The following comparison routine compares two floating point numbers (add the tolerance in the terms of the number the real numbers to make it approximate comparison).

/*Less than or greater than*/
char CompareFloats(float a, float b){
    int aint = *(int *)&a;
    int bint = *(int *)&b;
    aint = (aint & (1UL<<31))? ~(aint^(1UL<<31))+1:aint;
    bint = (bint & (1UL<<31))? ~(bint^(1UL<<31))+1:bint;
    aint -= bint;
    printf("%f %s %f \n",a,(aint <0)?"<":">=",b);
    printf("There are %d real numbers \n",(aint<0)?-aint:aint);
    printf("between %f and %f (in 32 bit float representation)\n",a,b);
}

[TECH] Launchpad and Bazaar

2008-12-18T20:22:00.001-08:00

I have now started to move all my development efforts to launchpad and bazaar. I really found launchpad a great place to maintain and release my work. I tried several alternatives before I moved to launchpad, I tried to use code.google but unfortunately its not as rich as launchpad and bazaar. I have released version 0.9 of libEditScript project whose aim to build a high performance and space efficient sequence alignment library. I have charted out the features for the next release, one thing I'm sure people would like to use is a JNI (Java Native Interface) so that people can use this space efficient code in their java code. I have seen that this problem of the need of edit script arises in several occasions, I have seen people using the space inefficient version O(n^2) of edit script computation. This is were we gain significantly in terms of the space saving. Also I have one more idea to make it more parallel in the sense I want to use some concurrent techniques to make this parallel, especially since I have circular queue I can easily use several threads to make this concurrent which I will add as a next step in the next release..

Checkout the libEditScript at https://launchpad.net/libeditscript. Next in my list would be to make a release of the weightedmatching library.

[TECH] A non-recursive algorithm to compute the Edit Script in O(min(n1,n2)+log(min(n1,n2)) space

2008-12-04T23:24:00.000-08:00

Recently I was working on an area efficient and synthesizable design to compute the edit script (minimum cost sequence of INSERT, DELETE and CHANGE) operations to transform string S₁ to S₂. After a little bit of thought I have the following idea to build a non-recursive version of the Hirschberg's algorithm, since the algorithm is non-recursive we can build an efficient digital circuit with this idea. We use a simple circular queue and apply DFS (Depth First Search) and we can prove that the capacity of this queue at any stage of the algorithm is θ(log(min(n₁,n₂)). The proof is simple we choose the geometrically decreasing string which has smaller length from the given strings(min(n₁,n₂)). Since we do a depth first search and the depth of subproblem tree is θ(log(min(n₁,n₂)), so we will have atmost θ(log(min(n₁,n₂)) subproblems in the circular queue at any stage of the algorithm.

The core non-recursive algorithm starts off with the initial problem and finds the minimum alignment split ('q') between (S₁[1.....n₁/2] and S₂ creates two sub problems and appends in the circular queue (currently it uses BFS). Download the linear space implementation from here Following is the outline of the core non recursive implementation.

     91     CQueue bfs_list;
     92     ESSubProblem es_sub;
     93 
     94     InitializeCQ(&bfs_list,sizeof(ESSubProblem));
     95     /*Put the toplevel subproblem into the CQueue*/
     96     es_sub.sub_s1 = s1; es_sub.sub_n1 = n1;
     97     es_sub.sub_s2 = s2; es_sub.sub_n2 = n2;
     98 
     99     AppendCQ(&bfs_list,&es_sub);
    100     while(DequeueCQ(&bfs_list,&es_sub)){
    101         sub_s1 = es_sub.sub_s1; sub_s2 = es_sub.sub_s2;
    102         sub_n1 = es_sub.sub_n1; sub_n2 = es_sub.sub_n2;
    103 
    104         if(sub_n1 <=1){
    105             /*If sub_n1 is 1 or 0 we cannot split it any more*/
    106             EditDistanceLS(sub_s1,sub_n1,sub_s2,sub_n2);
    107             i = sub_n1; j = sub_n2;
    108             /*Figure out the edit operations*/
    109             while(i || j){
    110                 op = GetOp(i,j,(i?sub_s1[i-1]:'\0'),
    111                     (j?sub_s2[j-1]:'\0'));
    112                 switch(op){
    113                     case 'I':
    114                             EditScriptS2[(&sub_s2[j-1]) - s2] = 'I';
    115                             j--;
    116                             break;
    117                     case 'D':
    118                             EditScriptS1[(&sub_s1[i-1]) - s1] = 'D';
    119                             i--;
    120                             break;
    121                     case 'C':
    122                             EditScriptS1[(&sub_s1[i-1]) - s1] =
    123                                 (sub_s1[i-1] == sub_s2[j-1])?':':'C';
    124                                 i--; j--;
    125                 }
    126             }
    127         }else {
    128             /*Add the two subproblems*/
    129             q = FindMinSplit(sub_s1,sub_n1,sub_s2,sub_n2);
    130             /*SUB PROBLEM 1*/
    131             es_sub.sub_n1 = CELING(sub_n1,2);
    132             es_sub.sub_n2 = q;
    133             AppendCQ(&bfs_list,&es_sub);
    134 
    135             /*SUB PROBLEM 2*/
    136             es_sub.sub_s1 = &(sub_s1[CELING(sub_n1,2)]);
    137             es_sub.sub_s2 = &(sub_s2[q]);
    138             es_sub.sub_n1 = sub_n1/2;
    139             es_sub.sub_n2 = sub_n2-q;
    140             AppendCQ(&bfs_list,&es_sub);
    141         }
    142     }

[TECH] A Graph with θ(log(n)) approximation ratio for greedy vertex cover

2008-10-24T11:26:00.000-07:00

We know that the greedy algorithm (which picks the vertex with next highest degree) to solve the vertex cover problem is not optimal. However we have 2-Approximation algorithm for the vertex cover problem (based on a integer linear programming formulation). Actually it was not very straight forward to for me prove that the greedy algorithm has an approximation ratio greater than 2. However after some thinking the following construction indeed proves that the greedy algorithms approximation ratio is > 2. In fact we can easily see (from the figure) that the greedy algorithm will first pick the node with highest degree i.e n (the magenta). Then it deletes all the edges adjacent on the magenta node and continues and clearly the greedy algorithm will have all the nodes in the blue bounding box. However its clear that if we choose the nodes in the oval region we get a minimum vertex cover of size n. The size of greedy vertex cover is n/2 + n/3 + n/4 .... n/n ≈ θ(nlog(n)) (for large n) and hence the approximation ratio for this graph is θ(log(n)).

[TECH] On The Uniqueness Of a Perfect Match In Undirected Acyclic graphs.

2008-10-16T14:48:00.000-07:00

I was studying about several results related to Tiling Theory. Tiling Theory adds formalism to the question " Given a simple rectangular region (R) can we find a valid tiling by using tiles from a set (T) ". If our tile set (T) is a singleton set with a unit square then we can tile any region. However what if the the tile set contains only a 2x1 rectangle (domino)? what if the tile set contains a set of polymino's (tetris game)?. Although the question may seem very simple there is a rich combinatorial work based on this question. In fact the study about the tiling problem originated from Wang's Conjecture about periodic tiling. Actually the algorithmic research about tiling theory was dormant for quite some time, even the proof for 2x1 dominoes is fairly recent see this . I was reading about the proof and was really impressed by its elegance. Often people consider simple things useless or don't have enough respect for simple facts, however I always find that good results are based on several very simple facts. In fact I would like to highlight simple ideas in every algorithm I study about.

In this proof to prove that the problem is NP-complete for 2x1 domino's a simple fact that " There is exactly only one perfect matching (if at all one exists) in a Tree " is used. This can be proved intuitively as follows.

Let M be the perfect match in the Tree T, if M is not a unique perfect match then there should be some other perfect match lets call it M^'. If we observe the leaf nodes in T which have only one edge adjacent on them so both M and M^' should agree on the matched edges which have leaf nodes of T as the vertices . We can now delete all the leaf nodes and their adjacent edges and continue the argument until we don't have any leaf nodes. Thus from this argument its clear that both M and M^' are the same and hence there is only a unique perfect match if at all if one exists in a tree T.

[TECH] Converting Market Matrix Sparse Representation to Row Compressed Sparse Representation

2008-10-15T18:26:00.000-07:00

I wanted to use some good benchmarks to illustrate the performance gains of the Fast Dual Update Algorithm for pre-ordering sparse matrices. University of Florida Sparse Library has a rich collection of sparse matrices from various domains. Unfortunately all the matrices are in Market Matrix format but not in Row Compressed format which I currently use.

The following is a handy script which converts the Market Matrix format to Row Compressed Format mm2rc.pl. The script uses a simple logic based on external sorting.

[TECH] Computing Eigen Values and Eigen Vectors using LAPACK

2008-09-24T15:20:00.000-07:00

I have used the LAPACK to compute the eigen values, till recently I needed the eigen vectors both left and right. I thought this could be really handy if you need to compute eigen values and eigen vectors. Download both library and this file from here

/*Computing both the eigen values and eigen vectors using LAPACK.
 *vamsi.krishnak (at) gmail (dot) com
 **/
#include "f2c.h"
#include "clapack.h"
#include<stdio.h>
#include<math.h>
#include<malloc.h>
#include<assert.h>

int main(){
    integer dim=0;
    unsigned int i;
    doublereal *A = NULL;
    real *eig_values = NULL;
    real *work;
    char JOBVL='V'; /*Left eigen vectors*/
    char JOBVR='V'; /*No right eigen vectors*/
    integer iwork[1];
    integer lwork,liwork=1,info,lda;
    doublereal *WR=NULL;
    doublereal *WI=NULL;
    doublereal *VL=NULL;
    doublereal *VR=NULL;
    integer LDVR=1;
    integer LDVL=1;
    integer LDA;
    doublereal *WORK=NULL;
    integer LWORK;


    printf("Please Enter the dimension of the array:\n");
    scanf("%d",&dim); LDA=dim; lda=dim; LDVL=dim;LDVR=dim;
    printf("Please enter the Matrix A\n");
    eig_values = (real *) malloc(sizeof(real)*dim);
    lwork = 10*dim; 
    work = (real *) malloc(sizeof(real)*lwork);
    A = (doublereal *) malloc(sizeof(doublereal)*(dim*dim));
    printf("Please Enter the elements of the matrix\n");
    i = 0;
    
    /*Double real and imaginary parts*/
    WR = (doublereal *)malloc(sizeof(doublereal)*dim);
    WI = (doublereal *)malloc(sizeof(doublereal)*dim);
    VL = (doublereal *)malloc(sizeof(doublereal)*(dim*dim));
    VR = (doublereal *)malloc(sizeof(doublereal)*(dim*dim));
    LWORK = 6*dim;
    WORK = (doublereal *)malloc(sizeof(doublereal)*LWORK);  
    assert(WI && WR && VL && VR);
    /**************/
    do{
        scanf("%lf",&A[i]);
        i++;
    }while(i<(dim*dim));
    /*Compute the Eigen values*/
    //ssyev_("N","U",&dim,A,&dim,eig_values,work,&lwork/*,iwork,&liworki*/,&info);
    dgeev_(&JOBVL,&JOBVR,&dim,A,&LDA,WR,WI,VL,&LDVL,VR,&LDVR,WORK,&LWORK,&info);
    if(!info){
        for(i=0;i<dim;i++){
            printf("%lf+i%lf\n",WR[i],WI[i]);
        }
        printf("===Right Eigen Vectors===\n");
        int k; int conj;
        for(i=0;i<dim;i++){
            /*Print i^th right eigen vector*/
            if(fabs(WI[i] > 0)){
                /*eigen vector i and eigen vector i+1*/
                for(conj=0;conj<2;conj++){
                    printf("[ ");
                    for(k=0;k<dim;k++){
                        if(!conj){
                            printf("(%lf)+i(%lf) ",VL[i*dim+k],VL[(i+1)*dim+k]);

                        }else{
                            printf("(%lf)-i(%lf) ",VL[i*dim+k],VL[(i+1)*dim+k]);
                        }
                    }
                    printf("  ]\n");
                }
                i++;
            }else{
                /*real eigen vector*/
                printf("[  ");
                for(k=0;k<dim;k++){
                    printf("%lf ",VL[i*dim+k]);
                }
                printf("] \n");
            }
            printf("\n");
        }
    }
}

[TECH] Cauchy's Inequality [continued...]

2008-09-18T22:04:00.000-07:00

Just adding the solution to one more problem continuing from last post.

[TECH] Cauchy's inequality a great tool for approximation algorithms

2008-09-17T19:42:00.000-07:00

Its been really long since I made a post, all these day's I'm doing a lot of theory and trying to come up with approximation algorithms to the Border Length Minimization Problem. After a long time I did prove the problem is NP-hard by reducing the T.S.P problem to the B.L.M.P problem. I'm now trying to find some approximation algorithms for the problem. I thought if I could strengthen my skills on inequalities it would be really great and I found a great problem book for inequalities by Michael Steele which has a good collections of problems on inequalities. In this post I provide my solutions to some of the problems which seemed interesting.

[TECH] The shopping problem is as hard as TSP.

2008-07-16T18:51:00.000-07:00

The Shopping Problem : Let I={i₁,i₂...i_k} be the items on your shopping list. Let S = {s₁,s₂,....,s_m} be the list of stores each store s_i sells only a subset of the items on your shopping list, and also the cost the items sold by the shops may vary i.e item i_j may be sold with different prices in different shops, so let C[i,j] be the cost of item i in store j. Since the shops are located in different places let D[p,q] be the cost of reaching store q from store p. So our aim is to start at home (H) and buy all the items we want and return to home with minimum possible cost. As I mentioned in the previous post this problem is was mentioned in Google Code Jam practice problems, however I did not see many submissions to this problem.

Firstly the problem is NP-hard because if items sold by each of the stores are disjoint and they are singleton sets then we should visit all the stores in the right order to minimize the total distance traveled, this is nothing by TSP. The following is the dynamic programming formulation for this problem. Let S(V,j) be a optimal solution to buy all the items in set V ending at store j. Then the following is the DP formulation.
S(V,j) = min { S(V-{k},p) + D[p,j] + C[k,j] } ; k \in V, 1<=p<=m Final Answer = min { S(I,j) + D[j,Home] } Initialization of the DP is simple you may want to look at the code below.

Further complexity to the problem can be added by saying that the items are classified as PerishableP and Non-PerishableN. So when ever you buy a perishable item you need to go back home to put in the refrigerator. The above DP formulation can be extended to this easily as follows. Let S(V,j,P) be the optimal sub problem of shopping all the items in the set V ending at the store j and having a Perishable item in the shopping cart, on the same lines S(V,j,N) optimal sub problem of shopping all the items in the set V ending at store j with NO perishable item in the shopping cart. Then the following formulation will give the optimal solution.

S(V,j,P) = min {S(V-{k},p,P) + D[p,Home] + D[Home,j] + C[k,j] } or
              min {S(V-{k},p,N) + D[p,j] + C[k,j] } 
If k is a perishable item

We can write similar formulations for Non-Perishable items also, get the code here implementing the shopping problem.

  1 /*Vamsi Kundeti (vamsi(dot)krishnak(at)gmail*/
  2 #include<stdio.h>
  3 #include<stdlib.h>
  4 #include<string.h>
  5 #include<assert.h>
  6 #include<math.h>
  7 #include "List.h"
  8 unsigned int **GRAPH;
  9 //#define VERBOSE_MODE 1
 10 unsigned int item_bits;
 11 char **item_names;
 12 unsigned int item_subset;
 13 /*has a unsigned int indicating
 14 *if this shop sells the item or not
 15 */
 16 unsigned int *shop_sells;
 17 /*dist[i][j] indicates the distance 
 18  *between i and j*/
 19 double **dist;
 20 /*cost[i][j] indicates the cost of
 21 *item i in shop j
 22 */
 23 double **cost;
 24 /*The number of the items*/
 25 unsigned int icount;
 26 /*The number of test cases*/
 27 unsigned int tcount;
 28 unsigned int caseno;
 29 /*The number of shops**/
 30 unsigned int scount;
 31 typedef struct _pair_ {
 32     int x;
 33     int y;
 34 }Coor;
 35 Coor *shop_coor;
 36 void EnumSubSetsBFS(unsigned int n);
 37 unsigned int LookUpItem(const char *key){
 38     unsigned int i=0;
 39     for(i=0;i<icount;i++){
 40         if(!strcmp(key,item_names[i])){
 41             return i+1;
 42         }
 43     }
 44     /*Ignore this item*/
 45     return 0;
 46 }
 47 unsigned int gas_price;
 48 unsigned int perish;
 49 inline char CheckPerish(unsigned int i){
 50     return (perish & (1<<(i)))?1:0;
 51 }
 52 void ReadInput(void){
 53     unsigned int i,j;
 54     unsigned int len;
 55     unsigned int icost;
 56     char item_buffer[64];
 57     char delim; 
 58     scanf("%u",&tcount);
 59     while(tcount--){
 60         perish=0;
 61         scanf("%u %u %u",&icount,&scount,&gas_price);
 62         assert(icount && scount);
 63         /*Read the items to be bought*/
 64         item_names = (char **) malloc(sizeof(char *)*icount);
 65         for(i=0;i<icount;i++){
 66             scanf("%s",&item_buffer);
 67             len = strlen(item_buffer);
 68             len = (len > 64)?64:len;
 69             item_names[i] = (char *) malloc(sizeof(char)*(len+1));
 70             assert(memcpy(item_names[i],item_buffer,(len+1))==item_names[i]);
 71             if(len && (item_names[i])[len-1]=='!'){
 72                 perish = perish | (1 << (i));
 73                 (item_names[i])[len-1]='\0';
 74             }
 75         }
 76         shop_coor = (Coor *) malloc(sizeof(Coor)*(scount+1));
 77         cost = (double **) malloc(sizeof(double *)*(scount+1));
 78         shop_sells = (unsigned int *)malloc(sizeof(unsigned int)*(scount+1));
 79         shop_coor[0].x = 0.0;
 80         shop_coor[0].y = 0.0;
 81         for(i=1;i<scount+1;i++){
 82             cost[i] = (double *) malloc(sizeof(double)*icount);
 83             scanf("%d %d",&(shop_coor[i].x),&(shop_coor[i].y));
 84             shop_sells[i] = 0;
 85             for(j=0;j<icount;j++){
 86                 cost[i][j] = 0;
 87             }
 88             while(1){
 89                 scanf(" %[^:]:%u",item_buffer,&icost);
 90                 if(j=LookUpItem(item_buffer)){
 91                     cost[i][j-1] = icost;
 92                     /*Also set the sells bit_set*/
 93                     shop_sells[i] = shop_sells[i] | (1<<(j-1));
 94                 }
 95                 scanf("%c",&delim);
 96                 if(delim == '\n'){
 97                     break;
 98                 }
 99             }
100         }
101         /*Now Compute the distance matrix*/
102         dist = (double **)malloc(sizeof(double *)*(scount+1));
103         for(i=0;i<scount+1;i++){
104             dist[i] = (double *)malloc(sizeof(double)*(scount+1));
105             for(j=0;j<scount+1;j++){
106                 dist[i][j] = sqrtf((shop_coor[i].x - shop_coor[j].x)*
107                     (shop_coor[i].x - shop_coor[j].x) +(shop_coor[i].y - shop_coor[j].y)*
108                         (shop_coor[i].y - shop_coor[j].y));
109             }
110         }
111 #if 0
112         printf("Items..\n");
113         for(i=0;i<icount;i++){
114             printf("%s ",item_names[i]);
115         }
116         printf("\n");
117         printf("Cost Matrix\n");
118         for(i=0;i<scount+1;i++){
119             printf("(%d,%d)",shop_coor[i].x,shop_coor[i].y);
120             for(j=0;j<icount;j++){
121                 if(cost[i][j]){
122                     printf(" %s[%c]:%lf",item_names[j],((CheckPerish(j))?'P':'U'),cost[i][j]);
123                 }
124             }
125             printf("\n");
126         }
127         printf("\n");
128 #endif
129         /*Make the call to the dynamic programming routine*/
130         EnumSubSetsBFS(icount);
131 
132         /*Freeup the stuff*/
133         for(i=0;i<icount;i++){
134             free(item_names[i]); 
135         }
136         free(item_names);
137         free(shop_coor);
138         free(shop_sells);
139         for(i=1;i<scount+1;i++){
140             free(cost[i]);
141         }
142         free(cost);
143         for(i=0;i<scount+1;i++){
144             free(dist[i]);
145         }
146         free(dist);
147     }
148 }
149 
150 void CreateGraph(unsigned int *n){
151     unsigned int i,j;
152     printf("Enter the Graph Size\n");
153     scanf("%u",n);
154     printf("Enter the weighted graph\n");
155     GRAPH = (unsigned int **) malloc(sizeof(unsigned int *)*(*n));
156     for(i=0;i<*n;i++){
157         GRAPH[i] = (unsigned int *) malloc(sizeof(unsigned int)*(*n));
158         for(j=0;j<*n;j++){
159             scanf("%u",&(GRAPH[i][j]));
160         }
161     }
162 }
163 
164 
165 unsigned int subset_count=0;
166 int main(){
167     /*Create the Graph*/
168     unsigned int n;
169     caseno=0;
170     ReadInput();
171     //EnumSubSetsBFS(n);
172 
173 }
174 /*Returns the new tail of the list*/
175 inline List* ListAppend(List **tail,void *_data){
176     List *new_node = (List *)malloc(sizeof(new_node)*1);
177     assert(new_node);
178     if(!(*tail)){
179         (*tail) = new_node;
180         (*tail)->_data = _data;
181         (*tail)->next = NULL;
182         return *tail;
183     }
184     new_node->_data = _data;
185     new_node->next = NULL;
186     (*tail)->next = new_node;
187     return new_node;
188 }
189 typedef struct _subset_{
190     /*if i^th bit is set then subset has i^th
191      *element*/
192     unsigned int ebits;
193     unsigned int max; /*Maximum element in the subset*/
194     unsigned int size; /*Size of subset*/
195     /*Pointer to the sub-problem S({},i,0) S({},i,1)*/
196     void *sub_problem;
197 }SubSet;
198 
199 /*Enumerate the Subset using BFS, n is 
200  *the set size [1,2...n] */
201 void EnumSubSetsBFS(unsigned int n){
202     unsigned int check_sub_set,k,i,j;
203     unsigned int current_level = 0;
204     List *bfs_list = NULL;
205     List *tail = NULL; List *free_node,*ptr;
206     SubSet *sub_set,*new_sub_set;
207     SubSet *one_sets = NULL;
208 
209     one_sets = (SubSet *)malloc(sizeof(SubSet)*icount); 
210     assert(one_sets);
211     /*Add the 1-subsets i=items*/
212     for(i=0;i<icount;i++){
213         one_sets[i].ebits = 0; one_sets[i].ebits |= (1<<i);
214         one_sets[i].max = i; one_sets[i].size = 1;
215         /****SUBPROBLEM SPECIFIC*****/
216         one_sets[i].sub_problem = (double **)malloc(sizeof(double *)*2);
217         /* Perishable S({i},X,0) */
218         ((double **)one_sets[i].sub_problem)[0] = (double *) malloc(sizeof(double)*(scount+1));
219         /* Unperishable S({i},X,1) */
220         ((double **)one_sets[i].sub_problem)[1] = (double *) malloc(sizeof(double)*(scount+1));
221         double **sub_problem_i = (double **)one_sets[i].sub_problem;
222 
223         /*j is the shop*/
224         for(j=1;j<scount+1;j++){
225             /*if item i can be bought at shop j
226              *depending on perishable or not initialize sub_problem[0] 
227              *sub_problem[1]
228              */
229              sub_problem_i[0][j] = HUGE_VAL; sub_problem_i[1][j] = HUGE_VAL;
230              if(shop_sells[j]&(1<<i)){
231                  if(CheckPerish(i)){
232                      sub_problem_i[0][j] = dist[0][j]*(gas_price) + cost[j][i];
233                  }else{
234                      sub_problem_i[1][j] = dist[0][j]*(gas_price) + cost[j][i];
235                  }
236              }
237 #ifdef VERBOSE_MODE
238             printf("S({%u},%u,0) = %lf\n",i+1,j+1,sub_problem_i[0][j]);
239             printf("S({%u},%u,1) = %lf\n",i+1,j+1,sub_problem_i[1][j]);
240 #endif
241         }
242         /*****************************/
243         if(!tail){
244             ListAppend(&bfs_list,(void *)&(one_sets[i]));
245             tail = bfs_list;
246             continue;
247         }
248         tail = ListAppend(&tail,(void *)&(one_sets[i]));
249     }
250 
251 
252 
253     /*Now do the BFS*/
254     while(bfs_list){
255         sub_set = (SubSet *)bfs_list->_data; 
256         free_node = bfs_list;
257         bfs_list=bfs_list->next;
258         subset_count++;
259         assert(sub_set->max < icount);
260         for(i=((sub_set->max)+1);i<icount;i++){
261             new_sub_set = (SubSet *) malloc(sizeof(SubSet)*1);
262             new_sub_set->size = (sub_set->size)+1;
263             new_sub_set->max = i;
264             new_sub_set->ebits = (sub_set->ebits | (1<<i));
265             /*This subset enumeration is critical to Subset based D.P*/
266 
267             /*****<begin>**SUBPROBLEM SPECIFIC******/
268             /* S(X,j) = min_k {S(X-k,k) + d[k][j]}*/
269             ptr = free_node;
270             double current_min = HUGE_VAL;
271             double **new_sub_problem;
272             double **sub_problem; double local_min;
273             double local_min_perish,local_min_unperish;
274             unsigned int k1,perish,shop1;
275             List *ptr_prev;
276             
277             new_sub_set->sub_problem = (double **) malloc(sizeof(double *)*2);
278             ((double **)new_sub_set->sub_problem)[0] = (double *) malloc(sizeof(double)*(scount+1));
279             ((double **)new_sub_set->sub_problem)[1] = (double *) malloc(sizeof(double)*(scount+1));
280             new_sub_problem = (double **)new_sub_set->sub_problem;
281 
282             for(perish=0;perish<=1;perish++){
283                 for(j=1;j<scount+1;j++){
284                     /*Objective here is to fill up S({new_sub_set},j,0) 
285                      *and S({new_sub_set},j,1)*/
286                     ptr = free_node;
287                     current_min = HUGE_VAL;
288                     for(k=1;k<=new_sub_set->size;k++){
289                         do{
290                             check_sub_set = ((SubSet *)ptr->_data)->ebits; 
291                             check_sub_set ^= new_sub_set->ebits; 
292                             k1 = (int)log2(check_sub_set);
293                             ptr_prev = ptr;
294                             ptr=ptr->next;
295                         }while(check_sub_set & (check_sub_set-1));
296 
297                         /*j is the shop and k1 is the item which you want
298                          *to buy at this shop, depending on
299                          *perishablity create the value of current min.
300                         */
301                         if(!(shop_sells[j] & check_sub_set)){
302                              continue; /*The shop don't sell this*/
303                         }
304                         sub_problem = (double **)
305                             ((SubSet *)ptr_prev->_data)->sub_problem;
306                         /*item k1 is sold by shop j*/
307                         local_min = HUGE_VAL;
308                         if((perish && !CheckPerish(k1)) 
309                                 || (!perish && CheckPerish(k1))){
310                             for(shop1=1;shop1<scount+1;shop1++){
311                                 /*S{{new_sub_set}-{k1},shop1} + 
312                                 d(shop1,h) + d(h,j) + c[k1][j]*/
313                                 if(j != shop1 || perish){
314                                     local_min_perish = 
315                                         sub_problem[0][shop1] + 
316                                         (dist[shop1][0] + dist[0][j])*
317                                             (gas_price) + cost[j][k1];
318                                 }else if(!perish){
319                                     local_min_perish = sub_problem[0][shop1] + 
320                                         cost[j][k1];
321                                 }
322                                 local_min_unperish = sub_problem[1][shop1] + 
323                                     (dist[shop1][j])*(gas_price) + cost[j][k1];
324                                 local_min_perish = 
325                                  (local_min_perish < local_min_unperish)? 
326                                     local_min_perish:local_min_unperish;
327 
328                                 if(local_min_perish < local_min){
329                                     local_min = local_min_perish;
330                                 }
331                             }
332                         }
333                         if(local_min < current_min){
334                             current_min = local_min; 
335                         }
336                     }
337                     new_sub_problem[perish][j] = current_min;
338                 }
339             }
340 
341 #ifdef VERBOSE_MODE
342             printf("Printing S(V,i) size=%u \n",sub_set->size);
343             /*Print SubSet*/
344             unsigned int i1;
345             for(i1=1;i1<scount+1;i1++){
346                 printf("S({");
347                 unsigned int j1;
348                 /*Get the elements back*/
349                 check_sub_set = new_sub_set->ebits;
350                 for(j1=0;j1<icount;j1++){
351                     if(check_sub_set & (1<<j1)){
352                         printf(" %u ",j1+1);
353                     }
354                 }
355                 printf("},%u,0)=%lf\n",i1+1,((double **)(new_sub_set->sub_problem))[0][i1]);
356                 printf("S({ABOVE},%u,1)=%lf\n",i1+1,((double **)(new_sub_set->sub_problem))[1][i1]);
357             }
358 #endif
359 
360             /**********<END>****************/
361             tail = ListAppend(&tail,(void *)new_sub_set);
362         }
363 
364         if(sub_set->size == icount){
365             double final_answer=HUGE_VAL;
366             double local_min; double **sub_problem;
367             sub_problem = (double **) sub_set->sub_problem;
368             for(j=1;j<scount+1;j++){
369                 local_min = (sub_problem[0][j]<sub_problem[1][j])?
370                     sub_problem[0][j]:sub_problem[1][j];
371                 local_min += (dist[j][0])*(gas_price);
372                 if(local_min < final_answer){
373                     final_answer = local_min;
374                 }
375             }
376             printf("Case #%u: %0.7lf\n",++caseno,final_answer);
377         }
378 
379         if(!(sub_set >= one_sets && sub_set <= &(one_sets[icount-1]))){
380             free(((double **)sub_set->sub_problem)[0]);
381             free(((double **)sub_set->sub_problem)[1]);
382             free(sub_set->sub_problem);
383             free(sub_set); 
384         }
385         free(free_node);
386     }
387     for(i=0;i<icount;i++){
388         free(((double **)one_sets[i].sub_problem)[0]);
389         free(((double **)one_sets[i].sub_problem)[1]);
390         free(one_sets[i].sub_problem);
391     }
392     free(one_sets);
393 //    printf("SubsetCount=%u\n",subset_count);
394 }
395

[TECH] A O(log(n)) time and O(1) space algorithm to find out Nth Fibonacci number.

2008-07-11T14:57:00.000-07:00

We often come across Fibonacci numbers in several contexts, if you see Knuth's "Concrete Mathematics" in the chapter on "Generating functions" one interesting example illustrates the possible brick patterns to build a wall whose width is 'n' using bricks of dimension 2x1. Also see the UVA problem here

We know a simple dynamic programming based algorithm can lead us to a O(n) time algorithm with O(1). But we can do much better than that in fact we can do it in O(log(n)) time and O(1) space. We can rewrite the dynamic programming as matrix formulation as follows.

F_n

F_n-1

1	1
1	0

F_n-1

F_n-2

This is really great we can get the closed form for F_n by repeated substitution. In fact the following is the closed form for F_n is as follows.

F_n

F_n-1

1	1
1	0

^n-1

So the problem of finding the n^th Fibonacci number is reduces to evaluating X^n-1 where X is a 2x2 matrix as shown above. Now the question is how we can evaluate X^n-1 as efficiently as possible. In fact to compute any a^b we need exactly |_{_}log(n)_{_}| multiplications by making use of the bit representation of b = (2^i₁ + 2^i₂ ...) where i_j are the positions of bits which have 1 in the binary representation of b. I used this trick to get a O(log(n)) time and O(1) time algorithm to compute the n^th Fibonacci number. Please see the code below get this here

  1 /*Find the n^th fibonacci number in O(log(n)) 
  2  *time, using the following idea of matrix formulation 
  3  * 
  4  * [   F_n  ]    [ 1      1  ] [F_{n-1}]
  5  * [        ] = [         ] [       ]
  6  * [ F_{n-1}]    [ 1      0  ] [F_{n-2}]
  7  *
  8  * vamsik(at)engr(dot)uconn July 11, 2008
  9  */
 10 #include<stdio.h>
 11 #include<stdlib.h>
 12 /*Fib returns the n^th Fibonacci number*/
 13 unsigned int Fib(unsigned int n){
 14     unsigned int deg=1;
 15     unsigned int fib[2][2] = {{1,1},{1,0}};
 16     unsigned int fibAns[2][2];
 17     unsigned char initAns=0;
 18     unsigned int mul[4];
 19     if(n<=2){ return (n<2)?n:2;}
 20     while(deg < n-1){
 21          mul[0] = fib[0][0]*fib[0][0] + fib[0][1]*fib[1][0];
 22          mul[1] = fib[0][0]*fib[0][1] + fib[0][1]*fib[1][1];
 23          mul[2] = fib[1][0]*fib[0][0] + fib[1][1]*fib[1][0];
 24          mul[3] = fib[1][0]*fib[0][1] + fib[1][1]*fib[1][1];
 25          fib[0][0] = mul[0]; fib[0][1] = mul[1]; 
 26          fib[1][0] = mul[2]; fib[1][1] = mul[3];
 27          deg <<=1;
 28          if(deg & n-1){ /*i^th bit set*/
 29              if(!initAns){
 30                  fibAns[0][0] = fib[0][0]; fibAns[0][1] = fib[0][1]; 
 31                  fibAns[1][0] = fib[1][0]; fibAns[1][1] = fib[1][1];
 32                 initAns=1;
 33              }else{
 34                  mul[0] = fibAns[0][0]*fib[0][0] + fibAns[0][1]*fib[1][0];
 35                  mul[1] = fibAns[0][0]*fib[0][1] + fibAns[0][1]*fib[1][1];
 36                  mul[2] = fibAns[1][0]*fib[0][0] + fibAns[1][1]*fib[1][0];
 37                  mul[3] = fibAns[1][0]*fib[0][1] + fibAns[1][1]*fib[1][1];
 38                  fibAns[0][0] = mul[0]; fibAns[0][1] = mul[1]; 
 39                  fibAns[1][0] = mul[2]; fibAns[1][1] = mul[3];
 40              }
 41          }
 42     }
 43     if(n-1 & 1){
 44         return fibAns[0][0] + fibAns[0][1];
 45     }
 46     return fibAns[0][0];
 47 }
 48 int main(int argc,char **argv){
 49     unsigned int n;
 50     while(1){
 51         scanf("%u",&n);
 52         if(!n){
 53             exit(0);
 54         }
 55         printf("%u\n",n,Fib(n));
 56     }
 57 }
 58

[TECH] Illustration of how to use the BFS based subset enumeration algorithmic framework to solve TSP.

2008-07-08T09:45:00.000-07:00

In the previous post I posted on the idea of using BFS to enumerate the subsets incrementally based on the size (see the previous post). The power of this core framework is so immense let me show how we can implement the TSP based on the framework. The extra code is between the lines 82-90 (TSP Initialization) and 113-175. The same core can be easily adapted to solve several sub-set based Dynamic programming optimizations like the job-scheduling, vehicle routing problem (VRP). If you have not looked into one of the Google Code Jam 2008 (practice) problems the problem#4 which is a shopping problem can be easily solved by this core framework. Check this code TSP.C (use -std=c99 to compile).

  1 /*******************************
  2  * Subset Enumeration Algorithmic 
  3  * Framework and its illustration in 
  4  * finding TSP. 
  5  ******************************/
  6 #include<stdio.h>
  7 #include<stdlib.h>
  8 #include<string.h>
  9 #include<assert.h>
 10 #include<math.h>
 11 #include "List.h"
 12 unsigned int **GRAPH;
 13 unsigned int subset_count=0;
 14 //#define VERBOSE_MODE 1
 15 void EnumSubSetsBFS(unsigned int n);
 16 
 17 void CreateGraph(unsigned int *n){
 18     unsigned int i,j;
 19     printf("Enter the Graph Size\n");
 20     scanf("%u",n);
 21     if(*n>32){
 22         fprintf(stderr,"TODO: The Bitset currently supports only 32 bits\n");
 23         exit(1);
 24     }
 25     printf("Enter the weighted graph\n");
 26     GRAPH = (unsigned int **) malloc(sizeof(unsigned int *)*(*n));
 27     for(i=0;i<*n;i++){
 28         GRAPH[i] = (unsigned int *) malloc(sizeof(unsigned int)*(*n));
 29         for(j=0;j<*n;j++){
 30             scanf("%u",&(GRAPH[i][j]));
 31         }
 32     }
 33 }
 34 
 35 int main(){
 36     /*Create the Graph*/
 37     unsigned int n;
 38     CreateGraph(&n);
 39     EnumSubSetsBFS(n);
 40 
 41 }
 42 /*Returns the new tail of the list*/
 43 inline List* ListAppend(List **tail,void *_data){
 44     List *new_node = (List *)malloc(sizeof(new_node)*1);
 45     assert(new_node);
 46     if(!(*tail)){
 47         (*tail) = new_node;
 48         (*tail)->_data = _data;
 49         (*tail)->next = NULL;
 50         return *tail;
 51     }
 52     new_node->_data = _data;
 53     new_node->next = NULL;
 54     (*tail)->next = new_node;
 55     return new_node;
 56 }
 57 typedef struct _subset_{
 58     /*if i^th bit is set then subset has i^th
 59      *element*/
 60     unsigned int ebits;
 61     unsigned int max; /*Maximum element in the subset*/
 62     unsigned int size; /*Size of subset*/
 63     /*Pointer to the sub-problem */
 64     void *sub_problem;
 65 }SubSet;
 66 /*Enumerate the Subset using BFS, n is 
 67  *the set size [1,2...n] */
 68 void EnumSubSetsBFS(unsigned int n){
 69     unsigned int check_sub_set,k,i,j;
 70     unsigned int current_level = 0;
 71     List *bfs_list = NULL;
 72     List *tail = NULL; List *free_node,*ptr;
 73     SubSet *sub_set,*new_sub_set;
 74     SubSet *one_sets = NULL;
 75 
 76     one_sets = (SubSet *)malloc(sizeof(SubSet)*n); 
 77     assert(one_sets);
 78     /*Add the 1-subsets*/
 79     for(i=1;i<n;i++){
 80         one_sets[i].ebits = 0; one_sets[i].ebits |= (1<<i);
 81         one_sets[i].max = i; one_sets[i].size = 1;
 82         /****SUBPROBLEM SPECIFIC*****/
 83         one_sets[i].sub_problem = (float *)malloc(sizeof(float)*n);
 84         for(j=1;j<n;j++){
 85             if(i==j){ continue; }
 86             /* TSP({i},j) */
 87             ((float *)one_sets[i].sub_problem)[j] = GRAPH[0][i]+GRAPH[i][j];
 88             printf("S({%u},%u) = %f\n",i+1,j+1,((float *)one_sets[i].sub_problem)[j]);
 89         }
 90         /*****************************/
 91         if(!tail){
 92             ListAppend(&bfs_list,(void *)&(one_sets[i]));
 93             tail = bfs_list;
 94             continue;
 95         }
 96         tail = ListAppend(&tail,(void *)&(one_sets[i]));
 97     }
 98     /*Now do the BFS*/
 99     while(bfs_list){
100         sub_set = (SubSet *)bfs_list->_data; 
101         free_node = bfs_list;
102         bfs_list=bfs_list->next;
103         subset_count++;
104         assert(sub_set->max < n);
105         for(i=((sub_set->max)+1);i<n;i++){
106             new_sub_set = (SubSet *) malloc(sizeof(SubSet)*1);
107             new_sub_set->size = (sub_set->size)+1;
108             new_sub_set->max = i;
109             new_sub_set->ebits = (sub_set->ebits | (1<<i));
110             
111             
112             
113             /*****<begin>**SUBPROBLEM SPECIFIC******/
114             /* S(X,j) = min_k {S(X-k,k) + d[k][j]}*/
115             ptr = free_node;
116             float current_min = HUGE_VAL;
117             float *new_sub_problem;
118             float *sub_problem; float local_min;
119             unsigned int k1,jstart,jend;
120             List *ptr_prev;
121             
122             new_sub_set->sub_problem = (float *) malloc(sizeof(float)*n);
123             new_sub_problem = (float *)new_sub_set->sub_problem;
124 
125             if(new_sub_set->size == n-1){
126                 jstart = 0; jend = 1;
127             }else{
128                 jstart = 1; jend = n;
129             }
130 
131             for(j=jstart;j<jend;j++){
132                 ptr = free_node;
133                 current_min = HUGE_VAL;
134                 if(new_sub_set->ebits & (1<<j)){
135                     /*Avoid computation of redundant subproblems like
136                      *S({....,x_i,x_j,x_k,....},x_i) */
137                     new_sub_problem[j]=HUGE_VAL;
138                     continue;
139                 }
140 
141                 for(k=1;k<=new_sub_set->size;k++){
142                     do{
143                         check_sub_set = ((SubSet *)ptr->_data)->ebits; 
144                         check_sub_set ^= new_sub_set->ebits; 
145                         k1 = (int)log2(check_sub_set);
146                         ptr_prev = ptr;
147                         ptr=ptr->next;
148                     }while(check_sub_set & (check_sub_set-1));
149                     sub_problem = (float *)((SubSet *)ptr_prev->_data)->sub_problem;
150                     local_min = sub_problem[k1] + GRAPH[k1][j];
151 
152                     if(local_min < current_min){
153                         current_min = local_min; 
154                         new_sub_problem[j] = current_min;
155                     }
156 
157                 }
158             }
159             /*Print Solution to the Subproblem */
160             unsigned int i1,istart,iend;
161             unsigned int j1;
162             
163             for(i1=jstart;i1<jend;i1++){
164                 check_sub_set = new_sub_set->ebits;
165                 if(check_sub_set & (1<<i1)){ continue; }
166                 printf("S(%u,{ ",i1+1);
167                 /*Get the elements back*/
168                 for(j1=0;j1<n;j1++){
169                     if(check_sub_set & (1<<j1)){
170                         printf(" %u ",j1+1);
171                     }
172                 }
173                 printf("})=%f\n",((float *)(new_sub_set->sub_problem))[i1]);
174             }
175             /**********<END>***SUBPROBLEM SPECIFIC*************/
176 
177 
178             tail = ListAppend(&tail,(void *)new_sub_set);
179         }
180         if(!(sub_set >= one_sets && sub_set <= &(one_sets[n-1]))){
181             free(sub_set->sub_problem);
182             free(sub_set); 
183         }
184         free(free_node);
185     }
186     for(i=0;i<n;i++){
187         free(one_sets[i].sub_problem);
188     }
189     free(one_sets);
190     printf("SubsetCount=%u\n",subset_count);
191 }
192

Sample Input.
---------test5.txt----------
5
0 3 1 5 4 
1 0 5 4 3 
5 4 0 2 1 
3 1 3 0 3 
5 2 4 1 0
---------------------
[vamsik@abadon ~]$ ./a.out < test5.txt 
Enter the Graph Size
Enter the weighted graph
S({2},3) = 8.000000
S({2},4) = 7.000000
S({2},5) = 6.000000
S({3},2) = 5.000000
S({3},4) = 3.000000
S({3},5) = 2.000000
S({4},2) = 6.000000
S({4},3) = 8.000000
S({4},5) = 8.000000
S({5},2) = 6.000000
S({5},3) = 8.000000
S({5},4) = 5.000000
S(4,{  2  3 })=9.000000
S(5,{  2  3 })=8.000000
S(3,{  2  4 })=10.000000
S(5,{  2  4 })=9.000000
S(3,{  2  5 })=10.000000
S(4,{  2  5 })=7.000000
S(2,{  3  4 })=4.000000
S(5,{  3  4 })=6.000000
S(2,{  3  5 })=4.000000
S(4,{  3  5 })=3.000000
S(2,{  4  5 })=6.000000
S(3,{  4  5 })=8.000000
S(5,{  2  3  4 })=7.000000
S(4,{  2  3  5 })=8.000000
S(3,{  2  4  5 })=10.000000
S(2,{  3  4  5 })=4.000000
S(1,{  2  3  4  5 })=5.000000
SubsetCount=15

[TECH] A note on enumerating subsets using BFS.

2008-07-04T19:09:00.000-07:00

I have several things which I want to make a note on my blog, but unfortunately I'm not getting enough time to write the blog posts. Any way.

Recently I came across several Dynamic Programming formulations which involve subsets in the Subproblem for example take the TSP problem, the Subproblem T(S-{j},j) is often used in the D.P formulation which indicates the optimal path starting from node '1' and covering all the nodes in the subset 'S' and ending in node 'j'. Here 'S' is a subset of V (set of all the nodes/vertices's of the given graph.). The final answer for the TSP problem is found by finding a 'k' such that T(V-{k},k) + d(k,1) is minimum.

If we observe we need a small algorithm framework which moves the solution from all the subsets of size '1' to '|V|'. In simple terms we need to enumerate the subsets in the order of their sizes. I found that a BFS can really help in enumerating the Subsets. The following code should be very handy in providing an enumeration framework and will be greatly useful to solve optimization problems which involve subsets in their D.P formulation.

#include<stdio.h> 
#include<stdlib.h> 
#include<string.h> 
#include<assert.h> 
#include<math.h> 
#include "List.h" 
#define VERBOSE_MODE 1 
 
/*Returns the new tail of the list*/ 
inline List* ListAppend(List **tail,void *_data){ 
    List *new_node = (List *)malloc(sizeof(new_node)*1); 
    assert(new_node); 
    if(!(*tail)){ 
        (*tail) = new_node; 
        (*tail)->_data = _data; 
        (*tail)->next = NULL; 
        return *tail; 
    } 
    new_node->_data = _data; 
    new_node->next = NULL; 
    (*tail)->next = new_node; 
    return new_node; 
} 
typedef struct _subset_{ 
    /*if i^th bit is set then subset has i^th 
     *element*/ 
    unsigned int ebits; 
    unsigned int max; /*Maximum element in the subset*/ 
    unsigned int size; /*Size of subset*/ 
}SubSet; 
/*Enumerate the Subset using BFS, n is  
 *the set size [1,2...n] */ 
void EnumSubSetsBFS(unsigned int n){ 
    unsigned int i; 
    List *bfs_list = NULL; 
    List *tail = NULL; List *free_node; 
    SubSet *sub_set,*new_sub_set; 
    SubSet *one_sets = (SubSet *)malloc(sizeof(SubSet)*n); 
    unsigned int current_level = 0; 
    /*Add the 1-subsets*/ 
    for(i=0;i<n;i++){ 
        one_sets[i].ebits = 0; one_sets[i].ebits |= (1<<i); 
        one_sets[i].max = i; one_sets[i].size = 1; 
        if(!tail){ 
            ListAppend(&bfs_list,(void *)&(one_sets[i])); 
            tail = bfs_list; 
            continue; 
        } 
        tail = ListAppend(&tail,(void *)&(one_sets[i])); 
    } 
    /*Now do the BFS*/ 
    while(bfs_list){ 
        sub_set = (SubSet *)bfs_list->_data;  
        free_node = bfs_list; 
        bfs_list=bfs_list->next; 
 
#ifdef VERBOSE_MODE  
        if(sub_set->size != current_level){ 
            current_level = sub_set->size; 
            printf("Printing SubSets of Size [%u] \n",current_level); 
            printf("===============================\n"); 
        } 
        /*print this set*/ 
        printf("[ "); 
        for(i=0;i<n;i++){ 
            if(sub_set->ebits & (1<<i)){ 
                printf("%u ",i); 
            } 
        } 
        printf("]\n"); 
#endif 
 
        assert(sub_set->max < n); 
        for(i=((sub_set->max)+1);i<n;i++){ 
            new_sub_set = (SubSet *) malloc(sizeof(SubSet)*1); 
            new_sub_set->size = (sub_set->size)+1; 
            new_sub_set->max = i; 
            new_sub_set->ebits = (sub_set->ebits | (1<<i)); 
            tail = ListAppend(&tail,(void *)new_sub_set); 
        } 
        if(!(sub_set >= one_sets && sub_set <= &(one_sets[n-1]))){ 
            free(sub_set);  
        } 
        free(free_node); 
    } 
    free(one_sets); 
} 
int main(int argc,char **argv){ 
    EnumSubSetsBFS(5); 
    return 0; 
}

Thanks to this http://puzzleware.net/codehtmler/default.aspx for converting my code.

Cheers! Vamsi.

[LIFE] It feels like I want to be a Grad student all my life.....Why does industry lack Passion?

2008-06-14T15:29:00.000-07:00

When do you enjoy most being a grad student? from my experience I found the following important facts.

Its passion which drives you, not Money.
When you start seeing the beauty behind what you study.
When you have an advisor who is simply great!

I feel that industry industry lacks passion, the passion and drive which a graduate student has is simply missing in the Industy. I don't know the reason why this happens is money in the industry making people Complacent? . I wanted to write my own views in this post.

The first year I'm into my graduate studies. I always felt that that I should get of UCONN as soon as possible just by getting a M.S degree. But as things went on I really started liking to be a graduate student. I started setting up my own goals rather than some one giving me list of things, and as I started get deep into the field of Combinatorics and Algorithms I started to feel that I need more mathematics for my analysis and now I decided to apply for a M.S degree in Mathematics. My mindset had definitely changed I started to feel that "THERE IS NO SHORTCUT TO SUCCESS". My desire for doing quality work has been getting stronger and stronger and now I'm in a situation that I cannot imagine myself without a PhD and some real and original contribution Algorithmic theory. I guess the only time in your life when sky seems a limit is only during graduate student life. If you really want you can be extremely productive and much professional than any guy working in the industry.

I think I have written soo much of code as a graduate student rather than in the industry, because in the industry when you have contributed to your project you feel good, your boss comes and praises you and then you start taking breaks and your productivity falls, you want to take breaks on weekends and I feel that people become very COMPLACENT once they get into industry. I find that passion and drive which a graduate has is totally lacking in industry.

I feel that one should not change if the environment around him keeps changing. I have so much of advice from my friends to apply for a big school for PhD rather than UCONN. I really don't understand why these guys are so afraid of life, I think these people have no confidence on them selfs but they feel secured only with TAG of a big school they think thats going to save them all the life. Lets take an example consider a PhD from the best university in the world call him X and consider a guy who is not at all a PhD and might have some undergraduate degree but has studied all the Volumes of Donald.E.Knuth and worked out all the exercises himself if far far far wiser than a a guy with a PhD from a top school with a mediocre dissertation.

[TECH] The C v C++ Flame which I have been looking for this time from Linus

2008-05-28T08:36:00.000-07:00

I was doing a Google search with keywords "C++ sucks" and found a flame on the git mailing list. Linus was totally angry by some guy who suggests that "git should have been written in c++ rather c". See the flame here

I have been looking for this kind of stuff from long time when I advocate that why I always stick to writing code in things like C,Perl and shell scripts. I feel that I have total control on the program if I write code in C in which I can deal with the virtual address's but unfortunately as the abstraction increases you cannot have the control on the actual code. I guess its not only me but *You cannot convince any hacker to use c++*

[TECH] A simple UNIX based algorithm to run a set of tasks in parallel.

2008-04-29T16:29:00.000-07:00

Given a set of tasks T={t₁,t₂....,t_n} such that these tasks can be run independently in parallel, however we have a limited set of resources and our constraint is K which tells us the maximum number of process's which can run in parallel. Our problem now is to realize this in UNIX using standard system calls like fork and wait.

I came across this recently and I came across the same thing several times before. Some examples of the contexts which would need this are as follows.

Running test cases (regressions) in parallel, suppose we have several thousands of test cases for a product and we wish to run on a machine which can run K tests at any given instance.
Making builds in parallel
Parallel crawling/web indexing (which was my motivation to write this)

Download this code

$#ARGV == 1 or die("USAGE: perl parallel_runner.pl {file_with_commands} {no.of.process's}");
$command_file = $ARGV[0];
open(CMD_FILE,$command_file) or die("$!");
$max_process = $ARGV[1];
($max_process < 256) or die("Too many process's....use less number");
%child_hash;
$process_count = 0;
while(<CMD_FILE>){
 $cmd = $_;
 if($cmd =~ /^\#/){
  next;
 }
 if($process_count < $max_process){
  $pid = fork();
  if($pid < 0 ){
   die("$! : RESOURCE ERRORS");
  }elsif($pid == 0){
   system("$cmd");
   exit(0);
  }else{
   $child_hash{$pid} = 1;
   $process_count++;
  }
 }

 if($process_count==$max_process){
  $pid = wait();
  if($child_hash{$pid} == 1){
   $child_hash{$pid} = 0;
  }elsif($pid == -1){
   print "No childs to wait?? Confused...quitting\n";
   last;
  } 
  $process_count--;
 }
}
print "Done creating process's...now waiting for childs...\n";
$i = 0;
while($i < $max_process){
 $pid = wait();
 if($pid == -1){
  last;
 }
 $i++;
}

[TECH] Computing all the derivatives of a 'n' degree polynomial at point 'x=a' in O(nlog(n)) time.

2008-04-23T17:13:00.000-07:00

I came across this problem recently, basically you are given a polynomial f(x) = a_nxⁿ+a_n-1x^n-1+....a₀ and you need to evaluate all its derivatives at a given point say x=a. I guess one of the suggested algorithms for this problem runs in O(nlog²(n)). But I found an idea to solve this in just O(nlog(n)). The following is the brief description of my algorithm

The idea is that we can post this as a simple multiplication of two n degree polynomials, the coefficients resulting from the multiplication will be the evaluated derivatives at a given point x=k. And as we might know that we can use FFT to multiply two n-degree polynomials in O(nlog(n))

STEP1: Compute the values of k²,k³....,kⁿ in O(n) time. Also compute the values of n!,(n-1)!,....2! incrementally. The runtime of this step is O(n)
STEP2: Create polynomial p(x) = (n!a_n)x^n-1+((n-1)!a_n-1)x^n-2+.....a₁.
STEP 3: Create polynomial g(x) = (k^n-1/(n-1)!)x+(k^n-2/(n-2)!)x²+(k^n-3/(n-3)!)x³+.... +kx^n-1 .
STEP 4: Multiply p(x) and g(x) using FFT in O(nlog(n) time.
STEP 5: Clearly coefficient of x^n-i+1 gives the value of i^th derivative.

[TECH] What does computing the edit distance between a string and its reverse tell us?

2008-04-19T15:12:00.000-07:00

I came across a problem recently which asks for the following "Given a string find out the minimum number of characters to be inserted to make it a palindrome, and also give the string". There might be several ways to solve this problem I have a simple algorithm to solve this problem in O(n²). The observation is the following reccurence.

Let LCS^r[x,y] = longest common sub sequence between S_1...x and S^r_1...y , where S^r is the reverse of the string S.
Let X=x₁x₂....x_k be the be the palindrome with minimum # of inserted characters to make it a palindrome, then X has to be formed from S as follows.
1. Find a index k in S such that the cost of transforming S[1..k] to S[k+1...n] is minimum and the palindrome would be X = S'[1...k]S'^r[k+1...n] or X = S'[1..k-1]S[k]S'^r[k+1..n]
2. The index k such that LCS^r[k][len-k] is maximum.
3. Get the code here

/*Builds the LCS between the forward and backwards and finds
a index 'p' which will give minimum # of operations to transform
forward string to backward string.*/
void FindLCSReverse(char *str,unsigned int len){
    int i,j;
    int current_min,max_lcs;
    int imax,jmax,palindex,palindex_r;
    unsigned int operations;
    char path=0;
    for(i=0;i<len;i++){
        LCS[i][0] = 0;
    }
    for(j=0;j<len;j++){
        LCS[0][j] = 0;
    }
    for(i=1;i<len;i++){
        for(j=1;j<len;j++){
            /*use len-j to access other string*/
            if(buf[i] == buf[len-j]){
                LCS[i][j] = LCS[i-1][j-1]+1;
            }else{
                LCS[i][j] = (LCS[i-1][j]>LCS[i][j-1])?LCS[i-1][j]:LCS[i][j-1];
            }
        }
    }
    /*Compute the p which gives minimum inserts to transform 
     *forward string to backward
     */
    max_lcs=0;
    imax=len-1;jmax=0;
    for(i=len-1;i>=1;i--){
        if(LCS[i][len-i-1] > max_lcs){
            max_lcs = LCS[i][len-i-1];
            imax = i;
            jmax = len-i-1;
        }
        if(i>0 && LCS[i-1][len-i-1] >= max_lcs){
            max_lcs = LCS[i-1][len-i-1];
            imax = i-1;
            jmax = len-i-1;
        }
    }
    /*Now find the actual string*/
    i=imax;
    j=jmax; 

    palindex_r=0;
    while(i!=0 || j!=0){
        if(i>0 && j>0){
            if(buf[i] == buf[len-j]){
                palindrome_rev[palindex_r++] = buf[i];
                i--; j--;
                continue;
            }
            current_min = (LCS[i-1][j] > LCS[i][j-1])?LCS[i-1][j]:LCS[i][j-1];
            if(current_min == LCS[i-1][j]){
                if(current_min == LCS[i][j-1] && buf[i] < buf[len-j]){
                    palindrome_rev[palindex_r++] = buf[len-j];
                    j--;
                }else{
                    palindrome_rev[palindex_r++] = buf[i];
                    i--;
                }
            }else if(current_min == LCS[i][j-1]){
                if(current_min == LCS[i-1][j] && buf[len-j] < buf[i]){
                    palindrome_rev[palindex_r++] = buf[i];
                    i--;
                }else{
                    palindrome_rev[palindex_r++] = buf[len-j];
                    j--;
                }
            }
        }else if(j==0){
            palindrome_rev[palindex_r++] = buf[i];
            i--;
        }else if(i==0){
            palindrome_rev[palindex_r++] = buf[len-j];
            j--;
        }
        /*path=1 (left) path=2 (down) path=3 (diag)*/
    }

    for(i=0;i<palindex_r;i++){
        palindrome[palindex_r-1-i] = palindrome_rev[i];
    }
    palindrome_rev[palindex_r] = '\0';
    palindrome[palindex_r] = '\0';

    /*printf("imax = %u jmax = %u len=%u \n",imax,jmax,len);*/
    if(imax+jmax==len-1){
        printf("%s%s\n",palindrome,palindrome_rev);
    }else{
        printf("%s%c%s\n",palindrome,buf[imax+1],palindrome_rev);
    }
}

[TECH] Reverse Engineering and Creating Crawler BOTS.

2008-04-15T23:02:00.000-07:00

I have my share of pleasure reverse engineering the underlying details of SCOPUS. SCOPUS as you might know is the most popular scholarly database used for citation searching. I was trying to solve this problem "Given a research paper X produce a set of research papers in the same connected component of the CITATION GRAPH", Let me define a CITATION GRAPH (CITE(V,E)) , V = {set of all research papers} and E={(i,j)} set of all directed edges from research paper 'i' to 'j' such that paper 'j' refers paper 'i' in its references. This edge information comes from SCOPUS however SCOPUS gives only one level (depth 1) in the connected component of all related papers, my goal is to get all the related papers (related in the sense fall in the same connected component of the CITE graph).

The reverse engineering the underlying comes handy when we want to automate the process of searching all these from the browser ourself.

I'm too tired to explain the details of the program which I had written using perl+LWP to create a CRAWLER BOT which gets all the related papers but if you need similar stuff sure the code can help click here

Unfortunately I don't get enough time to write blogs but in past few weeks I had some very interesting technical stuff I want to write.

[TECH] A simple algorithm to find the coefficient's of characteristic polynomial.

2008-03-12T22:31:00.000-07:00

Its often the case that the we need coefficients of the characteristic polynomial (xA = λx) rather than just the Eigen values and Eigen vectors of the matrix A, The following simple algorithm which can determine the coefficients of the polynomial from its roots (Eigen values) would be extremely handy, its a simple dynamic programming algorithm and also illustrates how quickly we can code dynamic programming algorithms in Matlab because the matrices resize automatically.



%
% eigen_poly_coeff: input 'A' should be a square matrix
% The program finds the eigen values of A, and constructs
% the caracteristic polynomial 
%
% output is the order of coefficients in the increasing 
% degree order.
%
function retval = eigen_poly_coeff (A)
v = eig(A)';
B = size(v);
coeff_matrix = eye(1,(B(1,2)+1));
coeff_matrix(1,1) = v(1,1)*-1;
coeff_matrix(1,2) = 1;

for i=3:1:(B(1,2)+1)
   coeff_matrix(1,i) = 1;
   for j=(i-1):-1:1
    coeff_matrix(1,j) = (coeff_matrix(1,j)*
                v(1,(i-1))*-1);
    if(j-1 >=1)
      coeff_matrix(1,j) = 
             coeff_matrix(1,j)+coeff_matrix(1,j-1);
    end
   end
end
 retval = coeff_matrix;
endfunction

[TECH] Algorithmic details of UNIX Sort command.

2008-03-09T23:58:00.000-07:00

I happened to look at the algorithmic details of UNIX Sort, a LINUX version of the classic UNIX sort is a part of GNU coreutils-6.9.90. This is classic example of the standard External R-Way merge , to sort a data of size N bytes with a main memory size of M so it creates N/M runs and merges R at a time, the number of passes through the data is log(N/M)/log(R) passes.In fact the lower bound(runtime) for external sorting is Ω((N/M)log(N/M)/log(R)). All the external memory sorting algorithms provided in the literature are optimal so the fight here is minimizing the constant before the number of passes.

UNIX sort treats keys are lines (strings), the algorithm followed by unix sort is in fact the R-Way merge. Let the input file size be IN_SIZE.
1. Choosing Run Size:
--------------------------------
The sizes of the initial runs are chosen from the total physical memory (TOTAL_PHY) and available memory (AVAIL_PHY). RUN_SIZE = (MAX(TOTAL_PHY/8,AVAIL_PHY))/2
maximum of 1/8th of TOTAL_PHY and AVAIL_PHY and divided by 2. See function "default_sort_size (void)" in the code.
2. Creating Runs:
-------------------------
Unix sort creates a temporary file for every run. So it creates IN_SIZE/RUN_SIZE (celing) temporary files. Internally it uses merge sort to sort internally it uses an optimization mentioned in Knuth volume 3 (2nd edition), problem 5.2.4-23.
3. Merging:
----------------
The number of runs merged at any time is hard coded in the program see macro NMERGE , NMERGE is defined to be 16 so it merges exactly 16 runs at any time.

[TECH] Sorting partially sorted sequences.

2008-03-05T16:12:00.000-08:00

Partially Sorted Sequence: A sequence k₁,k₂,...k_n of n keys such that any key k_j differs from its sorted position by atmost d as an example for d=3 we have 3,5,6,1,2,4 I came across this problem recently, turns out that we can solve this problem in several ways in O(nlog(d)) time. I found an interesting and simple way to solve this problem, a simple observation reveals that the smallest element in the entire sequence of n keys exists in the first d keys.


/*Build a heap(min) H with first d elements
 *in O(d) time
 */
 for(i=d;i<n;i++){
    DeleteMin(H); /*output this key*/
    InsertHeap(H,k_i);
 }
/*Note the Heap will have atmost d keys in it
at any time so the Insert and Delete can be done in 
O(log(d)) worst case.
*/

(n-d)log(d) = O(nlog(d))

[TECH] Converting an Las Vegas algorithm from an Expected Time Bound to High Probability bound.

2008-03-04T14:20:00.000-08:00

Some times its not always possible to derive the High Probability bounds (1-n^-α) for Las Vegas algorithms, I found an interesting way to convert any Las Vegas algorithm with expected bounds on the resources to the High Probability bounds. Assume that T_n be the expected runtime(or any resource) for a Las Vegas algorithm. The plain vanilla Las Vegas algorithm have the following structure.


while(1){
  /*Select a random sample*/
  if(answer found) quit;
}

Basically any analysis of the above randomized algorithm should tell how long should we run the loop. So let T_n gives the expected time on how long this loop runs. What we can do is as follows

Let the algorithm run for exactly 2T_n steps, if the algorithm quits before that we are fine, if the algorithm does not quit after T_n times then stop and restart.


counter = 0;
while(1){
 /*Pick a random sample*/
 if(counter++ == 2T_n){
    counter=0; /*restarting*/
 }
 if(answer found) quit; 
}

We can now show that the above algorithm will now give high probability bounds. Let X be the random variable which indicates the runtime of the algorithm. Then using Markov inequality P[X ≥ aT_n] ≤ 1/a

P[X ≥ 2T_n] ≤ 1/2
Lets assume that the algorithm runs for k ≥ 2T_n time steps, then we should have done a reset exactly k/(2T_n) times.
The probability of doing a reset equals P[X ≥ 2T_n] = P[reset].
Lets assume all the events are independent then the probability of resetting k/(2T_n) P[reset k/(2T_n) times] = P[algorithm running k time steps]
P[reset k/(2T_n)] ≤ (1/2)*(1/2)*(1/2)...k/(2T_n) times
P[reset k/(2T_n)] ≤ (1/2)^k/(2T_n)
we want the above probability to be very small (n^-α), to make that (1/2)^k/(2T_n) = n^-α). So the value of k for which this happens is 2T_nαlog(n)
So with high probability after 2T_nαlog(n) time steps if we quit the loop then we give the correct answer.
So we just have taken an expected bound T_n and converted into ∼O(log(n)T_n) time algorithm with high probability rather than expected bounds.

[TECH] Finding max and min in exactly 3n/2-2 Element Comparisions.

2008-02-29T23:17:00.000-08:00

We know that the lower bound on the minimum number of comparisons to find max and min in an array of n. In fact the recurrence T(n) = 2T(n/2)+2 solves exactly to 3n/2-2 the following is an interesting way to find max and min in exactly 3n/2-2 comparisons non-recursively(assume 'n' is even). Also note by comparison we mean the element comparisons as they are significant.


current_min = a[0]; current_max=a[1];
/* 1 comparison here */
if(current_min < current_max){
  current_min = a[1]; current_max = a[0];
}
/*(n-1)-2+1 times*/
for(i=2;i<n-1;i+=2){
  /*3 Element Comparisons here */
  if(a[i] < a[i+1]){
      if(a[i] < current_min) current_min = a[i];
      if(a[i+1] > current_max) current_max = a[i+1];     
  }else{
      if(a[i+1] < current_min) current_min = a[i+1];
      if(a[i] > current_max) current_max = a[i]; 
  }
}
/*return current_max current_min*/

Analysis: Total Comparisions = 1 + (((n-1)-2+1)/2)*3 = 3*n/2-2.

[TECH] Standard I/O Buffers won't get flushed when program crashes...

2008-02-12T15:08:00.000-08:00

Today I had a interesting problem, assume that you have a program printing out some details on to the screen and suddenly the program crashes and you want to capture what ever it has printed what would you do? you would try to redirect the output to a file and look at that file and what if the contents of the redirected file is empty? how do you explain this? try the following program and try to redirect what ever it prints using '>' or a '|' ("./a.out > out" or "./a.out | more")


char *crash_me = "crash this";
int main(){ 
    int i;  

    for(i=0;i<10;i++){
        printf("Before crash line %d\n",i+1);
        if(i==9){
            crash_me[6] = 'C';
        }
    }
}
/*TRY THE FOLLOWING*/
#1$./a.out
#2$./a.out > out
#3$./a.out | more

We can see that the 'out' file and 'more' don't show any thing which was in fact printed if we just run the program normally, how can we explain this? what exactly is happened? actually this is what has happened in the first case the terminal output is line buffered when ever terminal sees '\n' it prints that but in the next two cases the output is written to a file which is not line buffered (but gets written when the buffer is full) and since the program is crashed before the buffer gets full nothing is printed in case of #2 and #3. I guess one should definitely have a signal handler for SIGSEGV and other signals which end the program and do a explicit flush as below.


/*The buffers should be flushed before
 *program crashes.
 **/
void FlushBuffers(int sig){
    fprintf(stderr,"Segmentation Fault\n");
    fflush(stdout);fflush(stderr);
    exit(1);

}
char *crash_me = "crash this";
int main(){ 
    int i;  
    assert(signal(SIGSEGV,FlushBuffers)!=SIG_ERR);
    for(i=0;i<10;i++){
        printf("Before crash line %d\n",i+1);
        if(i==9){
            crash_me[6] = 'C';
        }
    }
}

I'm thinking of making a list of this kind of problems on UNIX

[TECH] SafeRead and SafeWrite

2008-02-10T18:50:00.001-08:00

'read' and 'write' are most frequently used system calls, its really a blunder to have something like the following in the code.


/*Read and Write blunders*/
if(read(fd,buffer,len) < 0 ){
  perror("Read ERROR:");
}

if(write(fd,buffer,len) < 0){
  perror("Write ERROR:");
}

The above code involving 'read' and 'write' syscalls may seem perfectly fine but unfortunately thats not true, NEVER IGNORE THE RETURN VALUE of a syscall, its very common to functions like 'read' and 'write' to return values less than 'len' in several situations like program interrupted by a signal or timeout or if the 'len' is very large say in 'Mb' make sure you have the following 'SafeRead' and 'SafeWrite' in your coding libraries.


ssize_t SafeWrite(int fd,void *buf,size_t wlen){
 ssize_t len; char *bbuf = (char *)buf;
 size_t writeln=0;
 while(writeln < wlen){
   len = write(fd,&(bbuf[writeln]),wlen-writeln);
   if(len<=0) return len;
   writeln += len;
 }
 return writeln;
}

ssize_t SafeRead(int fd,void *buf,size_t rlen){
 ssize_t len;char *bbuf = (char *)buf;
 size_t readln=0;
 while(readln < rlen){
   len = read(fd,&(bbuf[readln]),rlen-readln);
   if(len<=0) return len;
   readln += len;
 }
 return readln;
}

Its quite often that the while loop only executes once, but its always possible that due to some reason the 'read' and 'write' syscalls may return less than the number of bytes you read or write.

[TECH] A Contradictory gcc message.

2008-01-30T22:06:00.000-08:00

I have been experimenting on some of my PDM (parallel disk model) sorting code and wanted to create huge files with billions of keys I wanted to create a file with 1 billion integers (1024*1024*1024*4), the program stopped after some time saying that it "File size limit exceeded" , however in my shell (csh) the 'limit' command showed unlimited.



[vamsik@abadon PDMSorting]$ ./pdm_sort 
Setting RAND_SEED to 1201763856 
Filesize limit exceeded
[vamsik@abadon PDMSorting]$ ls -al key_file.txt 
---------- 1 vamsik fuse 2147483647 Jan 31 02:18 key_file.txt
[vamsik@abadon PDMSorting]$

I saw that the file was created without any permissions, so I tweaked my 'umask' but nothing changed (I was doing a open with (O_WRONLY|O_CREAT)). You might be wondering what exactly I'm trying to say in this post, in fact the story just started, rather than setting my 'umask' to 'umask 22', I have set it to 'umask 755' and as usual I was doing a build with 'make' this is what happened.


[vamsik@abadon PDMSorting]$ make
gcc  -O2  -I../myutil/            -c ExternalSort.c
ExternalSort.c:1: fatal error: can't open /tmp/cct9kuSr.s for writing: Permission denied
compilation terminated.
The bug is not reproducible, so it is likely a hardware or OS problem.
make: *** [ExternalSort.o] Error 1
[vamsik@abadon PDMSorting]$

Looks strange right? , see the funny thing it says "The bug is not reproducible, so it is likely a hardware or OS problem." , its really a stupid error message how can it say it cannot be reproduces? , just set 'umask 755' and do a 'gcc' on any file its going to say the same thing, in fact the message is a utter contradiction because we can reproduce this by changing the 'umask'. I'm going to report this to the 'gcc' maintainers or submit a patch to the maintainers.

[TECH] Ideas for Optimizing Design pattern implementations with Stack Collapsing

2008-01-14T19:25:00.000-08:00

One major drawback I guess with all these object oriented systems I guess is performance, since people write the code so that its extensible it always ends up creating deep stack sizes, recently I was looking into (Jmol) this structure visualization tool supports several input format descriptions for the structures (pdb,cif,molxyz....). This is how they hide format details from the display code
JmolAdapter (abstract class)
==> SmarterJmolAdapter (this contains a interface called AtomSetCollection ). So depending on the input file we have several AtomSetCollection readers like PdbAtomSetCollection etc...

Although I came to a conclusion that to implement the structural alignment algorithm into Jmol, I need to implement this AtomSetCollection which structurally aligns the protein structures, but I guess there are several draw backs of the implementation of object oriented java implementation, the cost of creating an extensible design does not come for free it comes at the cost of performance, I found that a simple command execution could end up creating a stack size of 32 this is what concerns me is there a way we can collapse the stack when we know that the intermediate functions on the stack are just delegating the call to the actual instance, I guess this STACK COLLAPSING technique can always be applied when ever the Adapter design pattern is used. For example see the below stack of size 32, its not doing any thing great its just opening a new script file, because of the way the code is written (use of design patterns) the stack size seem to increase a lot and ULTIMATELY ALL THESE PATTERNS ARE JUST DELEGATING A FUNCTION CALL TO TO FUNCTION ON TOP OF THE STACK, I'm think of ways to improve the delegation since at each level of the stack you are doing a look up and it increases linearly with the size of the stack.

  [1] org.jmol.viewer.ScriptManager.addScript (ScriptManager.java:69)
  [2] org.jmol.viewer.ScriptManager.addScript (ScriptManager.java:52)
  [3] org.jmol.viewer.Viewer.evalStringQuiet (Viewer.java:3,152)
  [4] org.jmol.viewer.Viewer.evalString (Viewer.java:3,119)
  [5] org.jmol.viewer.Viewer.openFile (Viewer.java:1,379)
  [6] org.openscience.jmol.app.Jmol$OpenAction.actionPerformed (Jmol.java:1,435)
  [7] javax.swing.AbstractButton.fireActionPerformed (AbstractButton.java:1,849)
  [8] javax.swing.AbstractButton$Handler.actionPerformed (AbstractButton.java:2,169)
  [9] javax.swing.DefaultButtonModel.fireActionPerformed (DefaultButtonModel.java:420)
  [10] javax.swing.DefaultButtonModel.setPressed (DefaultButtonModel.java:258)
  [11] javax.swing.AbstractButton.doClick (AbstractButton.java:302)
  [12] javax.swing.plaf.basic.BasicMenuItemUI.doClick (BasicMenuItemUI.java:1,051)
  [13] javax.swing.plaf.basic.BasicMenuItemUI$Handler.mouseReleased (BasicMenuItemUI.java:1,092)
  [14] java.awt.Component.processMouseEvent (Component.java:5,517)
  [15] javax.swing.JComponent.processMouseEvent (JComponent.java:3,135)
  [16] java.awt.Component.processEvent (Component.java:5,282)
  [17] java.awt.Container.processEvent (Container.java:1,966)
  [18] java.awt.Component.dispatchEventImpl (Component.java:3,984)
  [19] java.awt.Container.dispatchEventImpl (Container.java:2,024)
  [20] java.awt.Component.dispatchEvent (Component.java:3,819)
  [21] java.awt.LightweightDispatcher.retargetMouseEvent (Container.java:4,212)
  [22] java.awt.LightweightDispatcher.processMouseEvent (Container.java:3,892)
  [23] java.awt.LightweightDispatcher.dispatchEvent (Container.java:3,822)
  [24] java.awt.Container.dispatchEventImpl (Container.java:2,010)
  [25] java.awt.Window.dispatchEventImpl (Window.java:1,791)
  [26] java.awt.Component.dispatchEvent (Component.java:3,819)
  [27] java.awt.EventQueue.dispatchEvent (EventQueue.java:463)
  [28] java.awt.EventDispatchThread.pumpOneEventForHierarchy (EventDispatchThread.java:242)
  [29] java.awt.EventDispatchThread.pumpEventsForHierarchy (EventDispatchThread.java:163)
  [30] java.awt.EventDispatchThread.pumpEvents (EventDispatchThread.java:157)
  [31] java.awt.EventDispatchThread.pumpEvents (EventDispatchThread.java:149)
  [32] java.awt.EventDispatchThread.run (EventDispatchThread.java:110)

[TECH] Association rules of by Data Mining (TM Algorithm) on Cancer Data.

2008-01-05T20:39:00.000-08:00

I have datamined the Cancer data at http://breastscreening.cancer.gov/rfdataset/ using the TM(Transaction Mapping) and FP-Growth algorithm, this is what I have done to mine the association rules.

1. Randomly partition the data into two parts, I partitioned the data into part1 of size = 148458 records, part2 of size = 153897.
2. Used the part1 (148458 records) and found association rules of support >=0.4 and confidence >=0.4 , I got 72 rules from this.
3. For each of the rule (in step 2) I found the support and confidence of each of the rules in part2, it looks like the support and confidence is close to the support and confidence in training data (part1).


The 72 rules of step1 [
http://www.engr.uconn.edu/~vkk06001/CancerDataMining/rules.txt ]

Support and Confidence of each of this rules in part2
[http://www.engr.uconn.edu/~vkk06001/CancerDataMining/training_result.txt ]

I have made the rules human readable removing all the encoding please
see the rules
[
http://www.engr.uconn.edu/~vkk06001/CancerDataMining/human_readable.txt ]

These are in the following format
==============RULE:1=================
SUP:0.402 ,CONF:0.412,TRAIN_SUP:0.404,TRAIN_CONF:0.414
{
Diagnosis of invasive breast cancer within one year of the index
screening mammogram = no,
}
IMPLIES ===>
{
Diagnosis of invasive or ductal carcinoma in situ breast cancer within
one year of the index screening mammogram = no,
menopaus = postmenopausal or age>=55,
hispanic = no,
}
==============RULE:2=================

SUP indicates support of this rule in part2 , CONF indicates confidence of this
rule in part2, TRAIN_SUP indicates the support of this rule in part1 and
TRAIN_CONF indicates the confidence of this rule in part1.

These rules may not make any sense for me but it might make sense for a cancer doctor. There are several useful perl programs for people who want to do some datamining please feel free to use them http://www.engr.uconn.edu/~vkk06001/CancerDataMining , let me know if you have any questions.

[TECH] An observation to merge upper convex hulls in O(log(n)) time rather than O(log^2(n)).

2007-12-31T19:31:00.000-08:00

The well known divide and conquer algorithm for finding convex hull for a given set of points, needs to merge smaller hulls H_1 and H_2 for doing this merge the algorithm uses lemma 3.2 (in Book Fundamentals of Computer Algorithms). This lemma essentially finds the tangent line (u,v) in O(log^2(n)) time.

We can in fact merge the hulls in more efficient manner, the observation is based on the fact that any point with maximum y-coordinate should definitely be on the hull. Let p1 and p2 be the points in hulls H_1 and H_2 with maximum y-coordinate, then we can safely state that either p1 or p2 will be on the combined hull of H_1 and H_2, we can determine which of these two points lie on the merged hull by comparing the y-coordinates of p1,p2. Among these two the point with maximum y-coordinate will definitely be on the merged hull so one end of the tangent line (u,v) is fixed, the other end of the tangent line can be probed in O(log(n)) time, this basically saves the extra O(log(n)) time which is spent previously in lemma 3.2 for finding the other end of tangent in H_2

[TECH] Hanging simulation..

2007-12-25T15:04:00.000-08:00

I came across a interesting problem in which the entire simulation was hanging just because of an extra non-sensitive always block , in the verilog code the statement "always dummy_reg = in;" was causing both VCS and NC-VERILOG hang. Theoretically speaking irrespective of what I do in my design the simulation should stop after 500 steps because I have a "#500 $finish" in the initial block of the module "TestHangMux", although I got around this problem by removing the "dummy_reg" totally in the multiplexer, but I still don't understand why the simulation was hanging irrespective of the "#500 $finish" statement, is it because we have multiple non-sensitive "always" statement? the LRM(Language Reference Manual) says that all the "always" blocks execute in parallel just like parallel processors in that that the statement "#500 $finish" should have executed in parallel and stop the simulation but why it hangs?



===================================
module HangMux(in,sel,out);
input [1:0]in;
input sel; output out;
reg [1:0]dummy_reg;
reg out;

always @(sel) begin
    case(sel)
        0: out = dummy_reg[0];
        1: out = dummy_reg[1];
    endcase 
end

always dummy_reg = in;


endmodule

module TestHangMux(out_net);
output out_net;
wire out_net;
reg [1:0]in_reg;
reg sel;

initial begin
 in_reg = 2'b01;
 sel = 0;
#500 $finish;
end

always begin
#10 sel = ~sel;
end

HangMux hang_me (in_reg,sel,out_net);
endmodule
==============================

[TECH] Finding frequent itemsets faster than FP-Growth algorithm

2007-12-21T20:24:00.000-08:00

I wanted to implement my prime-number mapping intersection algorithm into the TM(Transaction Mapping) Algorithm, I got the FP-Tree implementation from PERL FP Tree and was trying to test this on a data of around 300,000 records, it works till 10,000 records with support and confidence of 0.8 and 0.9 and fails if the number of records > 20,000. I'm not sure if its a perl bug but the following code which is the cause of the problem seems difficult to understand



===================================================================
@association_rules = map $_->[0], (sort {$b->[1] <=> $a->[1]} 
      (map [$_, $_->confidence], @association_rules));      
===================================================================

I tried to bless the reference $_ but it did'nt work. The following is the error it gives.


[vamsik@abadon Tree-FP-0.04]$ !perl
perl test_vamsi.pl
Setting support
Mining rules...
Can't call method "confidence" on an undefined value at /usr/lib/perl5/site_perl/5.8.8/Tree/FP.pm line 397,  line 20001.

I need to fix this and send it to the FPTree maintainer tomorrow. I guess its really a issue to implement the FPTree in perl because of the speed and time especially when the algorithm need to work on huge data just in my case. A good idea is to re-write FPTree in C and extend TM based on that and compare the results.

[TECH] why is ((unsigned long)a_ptr+(unsigned long)b_int) not equals (a_ptr+(unsigned long)b_int)

2007-11-28T17:24:00.000-08:00

I was looking at some code today on SGI super computer, and wanted to make sure that there are no 32-64 bit problems in the code since SGI was super computer. I found the following thing to be interesting


<---------------------->
((unsigned long)a_ptr+(unsigned long)b_int)!=
          (a_ptr+(unsigned long)b_int)
<---------------------->

That is if I cast both the operands of binary + to (unsigned long) the value I get from that addition is different from if I cast just one operand of binary + to (unsigned long), whats the compiler doing here? since if both the operands to binary + are unsigned long and the value(no.of bits) could be larger than sizeof(unsigned long) is it the reason why I get different values here?, I couldn't find this in FAQ.

[TECH] Verilog incompatabilities between 'ncverilog' and 'vcs'

2007-11-23T19:27:00.000-08:00

Thanksgiving vacation may be fun for some but its really not fun if your IT department is taking a break, our Synopsys License server has been down for last few days all my emails are vain no one responds. I used 'vcs' to compile all my verilog code all these days, now I'm helpless.....how can we solve this problem.
I want to get some synthesis metrics for the Sequence Alignment design, although Synopsys License server is down looks like the Cadence License server is up :). Thats good news but whats equivalent to 'vcs' in Cadence....'NCVERILOG' I found 'ncverilog' and try to compile all my design modules now the in-compatibility, NC-Verilog don't like to see any thing in the module description until all the port list is filled up,but vcs scans the code of the entire declaration list before checking if the port list is filled up (I guess thats the way it should be the algorithm should be very general). I like 'vcs' and not impressed by 'NC-Verilog' but unfortunately I have to live with it for few more days, by porting my verilog code...sounds crazy I have ported code on machines different processors but the way these compilers are build is really crazy....

[TECH] Iso-spectral and Non-Isomorphic Graphs (PINGS)

2007-11-21T04:30:00.000-08:00

Having the same spectrum for the adjacency matrix is a Necessary but not sufficient condition for Graph Isomorphism, our recent algorithmic result states a conjecture that every graph is characterized by its Family of Spectra, rather than Spectrum itself.

I have been doing some work on efficient algorithms for generating PINGS (Pair of Isospectral and Non-Isomorphic Graphs), I wrote a program which can search for PINGS in a very efficient manner, I found some interesting results there are no PINGS for n=2,3,4 for n=5 we have one PING (see the picture) it happens that its the only PING for n=5 it has a symmetric spectrum of {-2,0,0,2}. Also I found a very interesting result that there cannot be more than two INGS(Iso-Spectral and non-isomorphic graphs) which share the same spectrum so if INGS exists they exist as PINGS, I was curious if they exist something like XINGS X=P or T or ... but it happens its just P (only a pair). There are 5 PINGS for n=6 for n=7 there are 55 PINGS see the list at http://trinity.engr.uconn.edu/~vamsik/n.7.

[TECH] Design of my first chip (A sequence Aligner O(n) space implementation)

2007-11-10T23:25:00.000-08:00

Designing hardware is totally different from writing software, writing verilog code is totally different from writing software 'C' programs, often a software engineer tries to solve the problems using loops,arrays etc.. unfortunately implementing the same algorithm in hardware is totally different. I had my experience in designing a chip (verilog code) for finding edit distance between two strings this has a well known O(n^2) dynamic programming based algorithm. Now the trick is how do you create a hardware which realizes the two for loops in the O(n^2) algorithm, after quite a bit of thinking I came up with a solution which involves just SHIFTERS,MULTIPLEXERS and ADDERS. I'm not aware of any such hardware implementation of edit distance algorithm till now...I'm tempted to share my design diagram on the blog but I cannot do it until it gets published some where....

My design flow is as follows (verilog)->vcs; (libs+verilog)->Design Compiler; and I will use cadence virtuso (ICFB) for my placement and routing and also extraction; Will use HSPICE/Nanosim and spectre for postlayout simulation.

[TECH] A simple O(n^2) time algorithm to construct the Ultra-metric trees.

2007-10-27T17:14:00.001-07:00

Given a matrix 'D' which is symmetric an Ultra-metric tree 'T' is a binary tree with internal nodes corresponding to D(i,j) and leaves corresponding to rows in the matrix 'D', and a path from the root to the leaf must be strictly decreasing (Ultra-metric tree). Also D(i,j) is the lowest common ancestor for 'i' and 'j'.
The problem is given 'D' we need to figure out corresponding 'T'.

The well know fact about the Ultra-metric trees is that if we take any 3 leaves i,j,k the maximum is not unique (e.x D(i,j)=5, D(j,k)=6,D(i,k)=5) we have D(i,j) and D(i,k) having same value.
How do we test if D is Ultra-metric?

A trivial solution is O(n^3).
The Errata shows construction of Ultra-metric trees in O(n^2) and asks a question if we can check for the Ultra-metric property in O(n^2) ? which my next item answers affirmatively.

My observation is we can partition the leaves of the ultra-metric tree with out any need of building a complete weighted graph as in Gusfields solution, one fact is that every row in D is going to have a maximum (global) we don't need to spend O(n^2) (comparing all the elements in the matrix D) to find one. The algorithm below we perform such a check if an ultra-metric tree exists and also builds one.


L = {1,2,....n};
UltraMetricPartition(L){
 i = select_a_random_element_in(L) ;
 max = find_max(D(i,*)); /*Takes linear time*/
 /*max = D(i,q) && q should be in L*/
 Part_q = {};
 Part_i = {};
 foreach(j in L){
    if(j!=i){
        if(D(j,i)<=max && D(j,q)==max){
             Part_i += {j};
        }else if(D(j,q)<=max && D(j,i)==max){
             Part_q += {j};
        }else{
           printf("D is not ultra-metric");
           return 0;
        }
    }
 }
 UltraMetricPartition(Part_i);
 UltraMetricPartition(Part_q);

}

Running time would be T(n) = T(n-k) + T(k) + O(n), clearly its quadratic O(n^2).

[TECH] Implementing Reliability(FEARLESS) on top of FUSE

2007-10-16T17:50:00.000-07:00

FUSE is a user space file system implementation framework, the file system implementation can run as a user process and interact with the VFS of the kernel through a named pipe /dev/fuse, FUSE comes by default with all kernels > 2.6.8

"Any data which persists on just one disk is very unreliable" , I was talking with Rohit couple of days back and he compared disk as an "Electric Bulb", and u know about the filament in the bulb it can go off any time. So what do you think of your laptop how many disk's does it have? JUST ONE.....thats extremely unreliable until and unless your system is not backed up every day. So how can we get the reliability achieved by RAID (Redundant Array of Inexpensive Disks) ? the answer is persist what is called an "ACTIVE DATA" on some thing like FLASH memory and we can integrate this with a backup mechanism to get good reliability for personal storage devices like laptops. FEARLESS is idea which came up with this reliability for personal storage devices, you can read paper by Dr.Chandy.

Our IDEA is now to realize FEARLESS using FUSE and rsync (for backup), so FEARLESS would be just another file system similar to SSHFS but very reliable

Cheers!
Vamsi.

[BDAY] unsigned char age = age++;

2007-10-08T14:13:00.000-07:00

Today happened to be my birthday, we just need an unsigned char to represent age, can you tell whats wrong with the code in the title (or below?).

static unsigned char age;

void Birthday(){
  age = age++;
}

Looks like the Standard 'C' says the value in the age is undefined ? because the assignment (=) is not a sequence point where the side effects (++,--) are settled. That can make the age of the person remain the same or increase it by one its ambiguous.

Cheers!
Vamsi. google index this

[TECH] O(n^2/(log(n)^2)) time algorithm for finding edit distance. [How does it feel when you know have reinvented the wheel :(]

2007-10-05T14:36:00.000-07:00

Yesterday night during the workout suddenly I thought about what will happen if we encode the t-vector {-1,0,1}in the Four Russian algorithm as an integer.....the idea was a BOOM! I really felt like EUREKA! just like Archimedes thought when he discovered the law of flotation, I felt that suddenly I speedup the edit-distance dynamic programming algorithm which is O(n^2) by a factor of O(log^2(n)), it was really great I spend working on the details all special cases and written up every thing was great till then. I discussed the solution with Dr.Wu and found that it was actually solution to some exercise problem in Gusfield's book which says that in a RAM based computation model the copy of vector 't' can be done in O(log(t)) time.

How do you feel when you know some thing which you are exited just slips from you ? its very painful and its like killing you, but thats the way the world is its full of SMART people and any path which you are taking might have been taken before........

[TECH] Generic Randomized Quicksort algorithm.

2007-09-20T23:44:00.000-07:00

Adding to my own collection of generic data structures, today I added a generic Randomized Quicksort algorithm, I need a good sorting routine in the SAPE algorithm for sorting distances. I want to makesure that it can sort any kind array (array of char's,int's,float's,double's and also structures).


/*A generic Randomized Quicksort:
 *The key idea of being generic comes
 *from the fact that you take a memory 
 *block and size of each element and 
 *sort this memory chunk, so we always
 *deal with address's.
 *
 *
 *Sep 20,2007 vamsik@engr.uconn.edu
 **/
#include
#include
#include
#include
#include
#include "RandQsort.h"

/*INPUT: swap_routine,compare_routine
 * swap_routine: The algorithm gives the 
 * (void *)'s swap_routine of two things
 * which need to be swapped.
 *
 * compare_routine: The algorithm gives the
 * (void *)'s to the compare_routine its 
 */
static char QSORT_RAND_SEED=0;
static void (*swap_routine)(void *,void *) = NULL;
static unsigned char (*compare_routine_l)(void *,void *) = NULL;
static unsigned char (*compare_routine_g)(void *,void *) = NULL;
static time_t t_rqsort;
unsigned int GetPivot(void){
 if(!QSORT_RAND_SEED){
  assert(((time_t)-1)!=time(&t_rqsort));
  srand(t_rqsort);
  fprintf(stdout,"Setting RAND_SEED to %d \n",t_rqsort);
  QSORT_RAND_SEED=1;
 }
 return ((unsigned int)rand());
}
/*Avoid loosing the pivot during swaps*/
void UpdatePivot(void **pivot,void **i,void **j){
 if(*i==*pivot){
  *pivot = *j;
 }else if(*j == *pivot){
  *pivot = *i;
 }
}
/*Take a memory address location and size 
of the data items[Inplace partition]*/
static void PartitionWithPivot(void* start,void* end,
unsigned int dsize){
 void *pivot=NULL;
 void* i=start;
 void* j=end;
 unsigned int p_index = GetPivot();
 if(start==end){
  return;
 }
 p_index %= (((unsigned int)(end-start))/dsize);
 pivot=(void *)((unsigned long)start+
  (unsigned long)((p_index)*dsize));
 while(i=start)){
   j-=dsize;
  }
  /*Swap i,j*/
  if(istart){
  swap_routine(pivot,(void *)(i-dsize));
  PartitionWithPivot(start,(void *)(i-(2*dsize)),
  dsize);
 }
 if((i+dsize)<=end){
  PartitionWithPivot(i,end,dsize);
 }
}
/*Call PartitionWithPivot recursively*/
void RandQsort(void *start,void *end,unsigned int dsize,
void_void_fptr sroutine,bool_void_fptr croutine_l){
 swap_routine = sroutine;
 compare_routine_l = croutine_l;
 assert(swap_routine && compare_routine_l);
 PartitionWithPivot(start,end,dsize);
}
#ifdef UNIT_TEST_RQSORT
/*TEST1:Testing for chars*/
unsigned char CharCompare_L(void *a,void *b){
 return ((*((char *)a))<=(*((char *)b)))?1:0;
}
unsigned char CharCompare_G(void *a,void *b){
 return ((*((char *)a))>=(*((char *)b)))?1:0;
}
void CharSwap(void *a,void *b){
 char temp=*((char *)a);
 *((char *)a) = *((char *)b);
 *((char *)b) = temp;
}
void TestCharSort(){
 char *data = malloc(sizeof(char)*10);
 char test_buffer[64];
 unsigned int i,j,k;
 time_t test_time;
 for(i=0;i<10;i++){
  data[9-i] = '0'+i;
 }
 /*Randomize the input data*/
 time(&test_time);
 srand(test_time);
 for(i=0;i<100;i++){
  j=(rand())%10;
  k=(rand())%10;
  test_buffer[63]=data[j];
  data[j]=data[k];
  data[k]=test_buffer[63];
 }
 /*Print the buffer*/
 if(strncpy(test_buffer,data,10)!=test_buffer){
  fprintf(stderr,"SYSTEM_ISSUE:Copy into buffer failed\n");
 }else{
  test_buffer[10]='\0';
  fprintf(stdout,"%s\n",test_buffer);
 }
 RandQsort((void *)data,(void *)(data+9),1,
 CharSwap,CharCompare_L);
 for(i=0;i<10;i++){
  if(data[i]-('0'+i)){
   printf("Test Failed for Chars\n");
   printf("Expecting %d but found %d\n",i,
    ((unsigned int)data[i]-'0'));
   exit(1);
  }
 }
 test_buffer[0]='\0';
 if(strncpy(test_buffer,data,10)!=test_buffer){
  fprintf(stderr,"FAILED_TEST: UNABLE TO COPY BUFFER\n");
  fprintf(stderr,"MAY_BE_SOME_MEMORY_CORRUPTION_DURING_SORT\n");
  exit(1);
 }else{
  fprintf(stdout,"%s\n",data);
 }
 printf("Test Passed For Chars....\n");
}

unsigned char FloatCompare(void *a,void *b){
 int a_int = (int)((*(float *)a)*10000);
 int b_int = (int)((*(float *)b)*10000);
 return (a_int<=b_int)?1:0;
}
void FloatSwap(void *a,void *b){
 float temp;
 temp = *((float *)a);
 *((float *)a) = *((float *)b);
 *((float *)b) = temp;
}

#include
void TestFloatSort(){
 float *data = malloc(sizeof(float)*10);
 float temp;
 unsigned int i,j,k;
 time_t test_time;
 for(i=0;i<10;i++){
  data[9-i] = (float)sin(i);
 }
 /*Randomize the input data*/
 time(&test_time);
 srand(test_time);
 for(i=0;i<100;i++){
  j=(rand())%10;
  k=(rand())%10;
  temp=data[j];
  data[j]=data[k];
  data[k]=temp;
 }
 /*Print the buffer*/
 for(i=0;i<10;i++){
  printf("%f ",data[i]);
 }
 printf("\n");
 printf("The # of elements %d\n",(&data[9]-&data[0])/4);
 RandQsort((void *)&data[0],(void *)&(data[9]),4,
 FloatSwap,FloatCompare);
 for(i=1;i<10;i++){
  if(!FloatCompare((void *)&data[i-1],(void *)&data[i])){
   printf("Test FAILED comparing %f , %f\n",
   data[i-1],data[i]);
   exit(1);
  }
 }
 for(i=0;i<10;i++){
  printf("%f ",data[i]);
 }
 printf("\n");
 printf("Test Passed For Floats....\n");
}

/*Test for floats*/
int main(){
 TestCharSort();
 TestFloatSort();
 exit(0);
}
#endif

[TECH] Perl template toolkit and Bugzilla

2007-09-08T11:58:00.000-07:00

I have now done with the customization of Bugzilla for our internal needs, to keep track of things. The core of the idea of templates was simple

Seperate the presentation from the code

The good thing about using Perl + TTK is that the split (presentation code and actual code) is very wide compared to any other framework which I'm aware at least from my past experience. I have used the following frameworks previously. Its a very useful add on as a skill learn more about http://www.template-toolkit.org

Have few things to fix today, Add code to SAPE_1 with the new Center of Gravity algorithm.

Cheers!
Vamsi

[TECH] Bugzilla Internals

2007-09-03T09:38:00.000-07:00

Got back to work today !

Was looking at the Bugzilla's perl code to extend it for some of my internal work, I found that it was a expert level perl code (Object oriented) they are using several things which I never used, It was a good learning experience.

I want extract few concepts from this code

How do all these websystem's avoid hard coding of HTML inside the code which generate's HTML dynamically (I know its kind of servlet and JSP/ASP situation), I found that these guys have some mechanism based on templates which avoids hard coding of the HTML inside the dynamic HTML generation code

I guess this is the key idea/concept behind all the code which generates HTML dynamically.

I was reading Sriram's "Advanced Perl Programming" I think its a great book the way he build's up the concepts, I was really impressed by a quote by him in the book
"It is indicative of inflexible procedural design if you find yourself using conditional statements to distinguish between object types"

I want to go out to get some lunch today :)
Cheers!
Vamsi

[NON-TECH] The old and cold feeling comes back........

2007-09-02T04:11:00.000-07:00

I visited India from Aug10-Aug30 all those 20 days were great most of the time I was in Hyderabad, life was great in India except that there was a lot of pollution especially in Hyderabad and also overcrowded places especially because of the booming IT industry.

Some of the things I noticed have changed in India from the time I left

The money spending power of people have increased a lot, I heard that the average salary for a undergrad is now around 5.0 lack/annum ($12.5k ), this was a huge increase compared to average of 2.4 lack/annum ($6k) which is more than double.
The cost of living has gone up from the last year, I remember that I used to eat in a hotel called Kakatiya in Ameerpet a full meal for Rs 26.00/- the same costs Rs 32.00/- now. A bottle of mineral water went up to 13.00/- from 10.00/-
There have been a lot of malls/multiplexes mushrooming every where in Hyderabad. I saw several existing houses being demolished and new multiplexes and malls are taking their place.
Looks to me that almost every thing is available in India now. I was shocked to see PSP in Music World, I guess people are now getting into games market and if its true India will be a great market for selling computer games.
In India cell phone has become a revolution I find that almost every one in the country has a mobile phone, thanks to Reliance which had a vision of making cell phones affordable to every one.

All these things are fine but my old feeling comes back again, I remember all those dilemma which I was having when I just landed in U.S went to riverside and came to UCONN, all those feeling are fresh in my mind. Well that COLD FEELING of missing India comes back again I don't know why but its back.....I feel sorry for myself I don't know what I'm up to in my life.....But I'm transforming my self to be better than what I used to be, and I believe that what ever is happening is just transformation for something better........as usual Connecticut is very lonely very quiet place, I used to go for jogging in India every day around 7-8 a.m in the morning.....

Well its life and lets face it without fear!
Cheer! Vamsi.

[NON-TECH] From Dennis Ritchie's Bio....

2007-08-02T14:55:00.000-07:00

Some times there are far too many things for people to decide on what to do? how to excel? what exactly is our strength ? . The problem here is that people tend to be *good* in several things, but being *good* is not enough if a person who is *good* compares to some of the smartest people in the same field. All this thoughts haunt every one I found the following from Dennis Ritchie's biography here
"My undergraduate experience convinced me that I was not smart enough to be a physicist, and that computers were quite neat. My graduate school experience convinced me that I was not smart enough to be an expert in the theory of algorithms and also that I liked procedural languages better than functional ones."
see that line "I was not smart enough to be an expert in the theory of algorithms....." I guess any person at his position would never say that (even though he is not smart), I really felt this guy is TRULY a GREAT GUY, he had to be really modest and humble to say such things......
Cheers!
Vamsi.

[TECH] My Quick VIM tip of the day

2007-07-03T09:04:00.000-07:00

To add something to every line in a file with vim
:s/^/something/gc

[TECH] FORTRAN awful stories

2007-07-03T00:35:00.000-07:00

Well first time in my life had a chance to look at a messy FORTRAN code. I used some of the CRAY style POINTER(ptr,pointee) stuff and it looked to me its really weird.
There are several crazy things which got me irritated while writing the code

Can you beleive that when you write FORTRAN code you must make sure that every line is no more than 80 chars, some one told me that FORTRAN had this because of the punch cards whose size limit come from there.
Really strange that FORTRAN has no concept of global variables, and its a real pain especially when you have SUBROUTINES you will have to make a COMMON block and copy all the definitions into the SUBROUTINE, real big pain, thanks to INCLUDE which saves us.
I found they there seems to be a huge variations in the standards of FORTRAN, there are several things 'g77' did'nt support, I used 'ifort' intel fortran compiler.
Really had tough time figuring out the format specifiers equivalent to 'printf' for 'WRITE(6,*)'
Last but not least FORTRAN dont have a ';' at the end of statement and its really pain when you are a 'C' programmers.
Last++ IF conditions seems to require THEN , I figured out after writing several lines of code :(, really an Experience of writing programs in legacy languages and a real test for a real programmer.

I forgot I wanted to keep track of this piece of code, to close all the files in a program with calling close(fd) in each of fd's

fp = tmpfile();
fp->_chain = stderr;
fpclose(fp);
fp = NULL;

Take care guys!
Cheers!
Vamsi

[TECH] Memory corruption and format specifiers "%s%c" is different from "%s%1s"

2007-06-28T13:19:00.000-07:00

I have been fixing several issues last few weeks, and the following issue in someone's code is a clear test for understanding differences between strings and char's.

On the first look "%s%c" and "%s%1s" seem very similar, but unfortunately NO! and they can create some nasty runtime bugs corrupting your variables, suppose the code existing in someone's code like this.

void BuggyScanner(){
    char buf[MAX_SIZE];
    int a;
    char b;
    scanf("%d",&a);
    scanf("%3s%1s",buf,&b);
    printf("a=%d buf=%s b=%c\n",a,buf,b);
}

The format specifier is really a blunder as scanf "%1s" is going to write beyond one byte (the extra '\0' which gets padded for strings)at the address of 'b' , since 'buf','a','b' are on the stack writing one byte beyond the address of 'b' can do really nasty stuff.

Just as in this corrupted the variables.
Potentially corrupt the return address of the function, creating a great security bug.

Be Careful guys! Vamsi.

[BOOKS] List of books I want to buy when I'm in india

2007-06-19T00:16:00.000-07:00

1. The Practice of Programming (Paperback) by Brian W. Kernighan (Author), Rob Pike (Author) 2. Programming Pearls (2nd Edition) (Paperback) by Jon Bentley (Author). 3. Algorithms on Strings, Trees and Sequences: Computer Science and Computational Biology (Hardcover)

I'll keep adding to this list...
Life has been through a lot of code in Matlab/Octave, Perl and C. Cheers! Vamsi.

[Perl] Don't nest your regular expressions when you use /g switch

2007-06-09T23:55:00.000-07:00

I have been busy all these days with quite a few things, well the motive of this blog is to document some of the issues I face in my day to day work.

while($line =~ /input .*? name="(.*?)" .*? value="(.*?)"/g){
  $param_name = $1;
  $param_value = $2;
  if($param_name =~ /blah/){
    # Do some stuff
  }else{
    #Do some stuff
  }
}

The above code can go miserably wrong and hard to debug especially when the outer while loop has hundreds of lines of code in it, as you can see when we use '/g' switch in the matching all the matches within a string, make sure that THERE ARE NO NESTED REGULAR EXPRESSIONS WHEN USING /g switch

I want to become expert in one area

2007-04-04T14:16:00.000-07:00

Some times I think what I'm I worth of? what can sell myself?
In fact I did sell myself when I created perl2exe based on my systems knowledge.
few things I think I should concentrate on
1. Parallel Algorithms (they make me to say ahaa!)
2. Systems (Drepper called me Moron! :()
3. EDA (Feels good for me).

[NON-TECH] Puppy

2007-03-21T15:27:00.001-07:00

Pic of puppy after a fashion show in hyderabad!
IMG_0495
Originally uploaded by vamsi.

[Tech] Bitonic Sequences and Bitonic Merge

2007-03-08T22:59:00.000-08:00

Given a Bitonic Sequence of length 2n (a motonically decreasing (inceasing) and montonically increasing (decreasing).
If we do a Bitonic merge we still of 2 bitonic sequences of length 'n' each.

void BitonicMerge(int *array,size_t size){
 int i;
 assert(!(size%2));
 for(i=0;i lessthan size/2;i++){
   if(a[i+(size/2)] < a[i]){
     swap(a[i+(size/2)],a[i]);
   }
}

What we have is now two Bitonic sequences. Doing a Bitonic sort on a parallel machines with 'p' processors and 'n' elements n >> p. Next post I'll post my code to do the Bitonic merge of the sequences.

[FUN] Top 21 things an Indian does after returning from U.S

2007-03-05T18:30:00.000-08:00

(--- Good One!.... Just for fun---)

21. Tries to use credit cards in a road side hotel.

20. Drinks and carries mineral water and always speaks of being health conscious.

19. Sprays deo so that he doesn't need to take bath.

18. Sneezes and says 'Excuse me'.

17. Says "Hey" instead of "Hi".
says "Yogurt" instead of "Curds".
Says "Cab" instead of "Taxi".
Says "Candy" instead of "Chocolate".
Says "Cookie" instead of "Biscuit".
Says "Free Way" instead of "Highway".
Says "got to go" instead of "Have to go".
Says "Oh" instead of "Zero", (for 704, he will say Seven Oh Four Instead of Seven Zero Four)

16. Doesn't forget to crib about the air pollution. Keeps cribbing every time he steps out.

15. Says all the distances in Miles (Not in Kilo Meters), and counts in Millions. (Not in Lakhs)

14. Tries to figure all the prices in Dollars as far as possible (but deep inside multiplies by 43).

13. Tries to see the % of fat on the cover of a milk pocket.

12. When he needs to say Z (zed), he never says Z (Zed), instead repeats "Zee" several times, and if the other person is unable to get it, then says X, Y Zee(but never says Zed)

11. Writes the date in MM/DD/YYYY. On watching traditional DD/MM/YYYY, says "Oh! British Style!!!!"

10. Makes fun of Indian Standard Time and the Indian Road Conditions.

9. Even after 2 months, complaints about "Jet Lag".

8. Avoids eating spicy food.

7. Tries to drink "Diet Coke", instead of Normal Coke. Eats Pizza instead of Dosa.

6. Tries to complain about any thing in India as if he is experiencing it for the first time. Asks questions etc. about India as though its his first visit to India.

5. Pronounces "schedule" as "skejule", and "module" as "mojule".

4. Looks suspiciously towards any Hotel/Dhaba food.Few more important ones:

3. From the luggage bag, does not remove the stickers of the Airways by which he traveled back to India, even after 4 months of arrival.

2. Takes the cabin luggage bag to short visits in India and tries to roll the bag on Indian Roads.The Ultimate one

1. Tries to begin any conversation with "In US ...." or "When I was in US..."

[Tech] How to determine the bus bandwidth.

2007-02-27T17:48:00.000-08:00

Hi recently I had a need to find out the bus bandwidth. 
The bus had 64 lines and assume that the speed of the bus is xMHz.
This implies the bandwidth = (no.of lines) * (speed of each line).
Bus Bandwidth = ((64/8)*x*10^6)/(2^20) MBytes/sec.

A simple use of this metric is suppose you have 
SMP(Shared memory processors) with a common bus
and the bandwidth of the bus is 1 GBytes/sec
and you have several processors which have an 
(memory) Instruction bandwidth of 250 MBytes/sec
you cannot connect more than 4 processors to 
saturate the bus.

[Tech] Pathetic non unix standard matlab editor

2007-02-22T00:00:00.000-08:00

Well I'm really pissed of this damn editor (matlab editor) , some how I cannot execute my program .m file from the command line so I had to use this matlab's editor to run my program, after using vi for long this is really pathetic, its more worse than windows I just did ctrl+c and ctrl+v it pased some junk instead of the code which I intended to paste.

[SPICE2LAYOUT] 40% Complete

2007-02-19T21:51:00.000-08:00

Got back to the SPICE2LAYOUT project updated some Makefiles and added all the skeletal code today. By the mid of the week I should get this up and running, It also needs to have the channel router in that.
Need to get back home and study some stuff for a upcoming Friday jury.

Algorithm for Permuting in place

2007-02-18T15:29:00.000-08:00

Problem: Given a set of numbers S = 1...n and a permutation PI(S), how can you rearrange elements in S in O(n) time with no extra space.

I don't like explaining things just code it., here's code below

#include
#include
#include
int a[10];
int b[10];

void ResetArray(){
    unsigned int i;
    for(i=0;i<10;i++){
        a[i] = i; 
    }
}
void swap_elements(unsigned int i,unsigned int j){
    int temp;
    temp = a[i];
    a[i] = a[j];
    a[j] = temp;
}
/*b is the swap list*/
void InplacePermute(int *perm){
    unsigned int i;
    for(i=0;i<10;i++){
        if(perm[i] < i){
            perm[i] = perm[perm[i]];
        }
        if(a[i] < i){
            perm[a[i]] = perm[i];
        }
        swap_elements(i,perm[i]);
    }
}

void PrintPermutation(void){
    unsigned int i;
    for(i=0;i<10;i++){
        printf("%d ",a[i]);
    }
}

int main(){
    int runs=10;
    unsigned int i;
    do{
        ResetArray();
        printf("Please enter the permutation of size 10, 10 unique digits from 0-9\n");
        for(i=0;i<10;i++){
         if(scanf("%d",&b[i])<=0){
            break;
         }
         assert(b[i] >=0 && b[i] <=10);
        }
        if(i<10){
            break;
        }
        InplacePermute(b);
        printf("The permutation is::");
        PrintPermutation();
    }while(1);
}

The following are testcases

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

Well this week has been a little productive, got the space for http://phaedrus.sourceforge.net a library of randomized algorithms
Have Fun.... Vamsi.

Passing multidimensional arrays to functions

2007-02-15T00:31:00.000-08:00

Got enlightened by this today. And realized the following.

int a[10][20];

int main(){
 print_array(a); /*results in a crash.*/

}

void print_array(int **a,int i,int j){
 printf("%d",a[i][j]); /*How stupid this can be?*/
}

I love PHP

2007-02-12T12:25:00.000-08:00

Have been using PHP with apache to build some stuff really like programming in this I'am in love with PHP.

[Tech] One more reason why MACRO's are Not Preferred over inlines.

2007-01-30T21:12:00.000-08:00

...macro.h....
/***Dangerouse Macro1***/
#define print_int_value_twice(c) do{\
 printf("first time %c second time %c\n",c,c);\
}while(0);

/***Stupid Preprocessor***/
#define my_error(c) do{\
 fprintf(stderr,"Error::");
 fprintf(stderr,c);
 fprintf(stderr,"\n");
 exit(1);
}while(0);

.....macro user ....
j=0;
print_int_value_twice(j++);
/*j incremented twice by the preprocessor*/


/*In this case the preprocessor is stupid 
 not treat the string spanned across several lines
 as a incomplete macro argument and fails compilation.*/
my_error(" My error spans three lines
            line1
            line2");

Avoid MACROS start using inlines , if you define macros the users should use them with caution

[Personal] I was not good and never have been.....

2007-01-29T19:45:00.000-08:00

Yes! today I realize my potential, I was never good in fact never have been....I accept my failure. I'm a looser, I cannot win and never have won and in fact I'm obliged to life for letting me live and eat what ever I want and wear what ever I could....in the midst of bitter shame....I gave up...I don't want to live as a looser...I will soon end my life and this blog may be the last. I was never honest to myself I have things which I don't deserve. And I get emotional for small things...what else a looser can say except bending down head in shame.

Silent rants of my life.........

[TECH] Enumerating set partitions with $DLX$

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

[THEORY] Is nature a giant guessing algorithm ?

[TECH] Unbuffered message logging of interactive programs.

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

[TECH] External R-Way Merge sorting.

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

[TECH] Super-fast tokenizing of sequence word patterns

[PHIL] My best lines in "The zen and art of motor cycle maintenance"

[TECH] Tutte matrix and Sufficiency of Perfect Matching.

[LIFE] Keeping the PACE of productivity is very hard.

[TECH] An efficient External sorting API

Rounding An Integer To The Next Maximal Mutliple Of A Given Radix Power

[TECH] An interesting problem on the mailing list [gdb@gnu.org]

[TECH] Manipulating DNA sequences with bits

[TECH] Loops in Bi-directed De Bruijn Graphs

[TECH] Simplified Proof For The Application Of Freivalds' Technique to Verify Matrix Multiplication

[TECH] Parallel sorting of binary keys

[PHIL] My original quote "The more you learn the more stupid you feel"

[TECH] Border Length Minimization Solver

Algebraic methods for Combinatorial counting.

[TECH] Maintaining a Load Factor (α) in Dynamic Data structures with constant (Θ(1)) amortized cost per operation.

[TECH] Exponential function (e-μt) is the only unique solution satifying R(x+t) = R(x)R(t)(semi-group property).

[TECH] Simple comparision of real numbers

[TECH] Launchpad and Bazaar

[TECH] A non-recursive algorithm to compute the Edit Script in O(min(n1,n2)+log(min(n1,n2)) space

[TECH] A Graph with θ(log(n)) approximation ratio for greedy vertex cover

[TECH] On The Uniqueness Of a Perfect Match In Undirected Acyclic graphs.

[TECH] Converting Market Matrix Sparse Representation to Row Compressed Sparse Representation

[TECH] Computing Eigen Values and Eigen Vectors using LAPACK

[TECH] Cauchy's Inequality [continued...]

[TECH] Cauchy's inequality a great tool for approximation algorithms

[TECH] The shopping problem is as hard as TSP.

[TECH] A O(log(n)) time and O(1) space algorithm to find out Nth Fibonacci number.

[TECH] Illustration of how to use the BFS based subset enumeration algorithmic framework to solve TSP.

[TECH] A note on enumerating subsets using BFS.

[LIFE] It feels like I want to be a Grad student all my life.....Why does industry lack Passion?

[TECH] The C v C++ Flame which I have been looking for this time from *Linus*

[TECH] A simple UNIX based algorithm to run a set of tasks in parallel.

[TECH] Computing all the derivatives of a 'n' degree polynomial at point 'x=a' in O(nlog(n)) time.

[TECH] What does computing the edit distance between a string and its reverse tell us?

[TECH] Reverse Engineering and Creating Crawler BOTS.

[TECH] A simple algorithm to find the coefficient's of characteristic polynomial.

[TECH] Algorithmic details of UNIX Sort command.

[TECH] Sorting partially sorted sequences.

[TECH] Converting an Las Vegas algorithm from an Expected Time Bound to High Probability bound.

[TECH] Finding max and min in exactly 3n/2-2 Element Comparisions.

[TECH] Standard I/O Buffers won't get flushed when program crashes...

[TECH] SafeRead and SafeWrite

[TECH] A Contradictory gcc message.

[TECH] Ideas for Optimizing Design pattern implementations with Stack Collapsing

[TECH] Association rules of by Data Mining (TM Algorithm) on Cancer Data.

[TECH] An observation to merge upper convex hulls in O(log(n)) time rather than O(log^2(n)).

[TECH] Hanging simulation..

[TECH] Finding frequent itemsets faster than FP-Growth algorithm

[TECH] why is ((unsigned long)a_ptr+(unsigned long)b_int) not equals (a_ptr+(unsigned long)b_int)

[TECH] Verilog incompatabilities between 'ncverilog' and 'vcs'

[TECH] Iso-spectral and Non-Isomorphic Graphs (PINGS)

[TECH] Design of my first chip (A sequence Aligner O(n) space implementation)

[TECH] A simple O(n^2) time algorithm to construct the Ultra-metric trees.

[TECH] Implementing Reliability(FEARLESS) on top of FUSE

[BDAY] unsigned char age = age++;

[TECH] O(n^2/(log(n)^2)) time algorithm for finding edit distance. [How does it feel when you know have reinvented the wheel :(]

[TECH] Generic Randomized Quicksort algorithm.

[TECH] Perl template toolkit and Bugzilla

[TECH] Bugzilla Internals

[NON-TECH] The old and cold feeling comes back........

[NON-TECH] From Dennis Ritchie's Bio....

[TECH] My Quick VIM tip of the day

[TECH] FORTRAN awful stories

[TECH] Memory corruption and format specifiers "%s%c" is different from "%s%1s"

[BOOKS] List of books I want to buy when I'm in india

[Perl] Don't nest your regular expressions when you use /g switch

I want to become expert in one area

[NON-TECH] Puppy

[Tech] Bitonic Sequences and Bitonic Merge

[FUN] Top 21 things an Indian does after returning from U.S

[Tech] How to determine the bus bandwidth.

[Tech] Pathetic non unix standard matlab editor

[TECH] The C v C++ Flame which I have been looking for this time from Linus