## The “BitVectors” Problem

### (From ITA Problem Bank)

Statement of Problem:

The BitVectors are an ancient and immortal race of 10,000, each with a 10,000 bit genome. The race evolved from a single individual by the following process: 9,999 times a BitVector chosen at random from amongst the population was cloned using an error-prone process that considers each bit independently, and flips it with 20% probability. Given a set of BitVectors, determine their reproductive history - i.e. The "Progenitor" BitVector and the BitVector Family tree. Your program's output should be, for each input line, the 0-based line number of that individual's parent, or -1 if it is the progenitor. Balance performance against probability of mistakes as you see fit. The randomly-ordered file bitvectors-genes.data.gz contains a 10,000 bit line for each individual. (However, for testing purposes, I have also included code that can be used to generate other families of BitVectors so that my program's output can be compared to ground truth.

Solution:

My solution strives to both minimize the number of comparisons between BitVectors and to make each comparison as efficient as possible.

The progenitor BitVector is the ancestor of the entire population and has had the most opportunity to clone itself. So, I consider the bit-wise mode of the population to be an idealized version of the progenitor (e.g. the bit-wise mode of [0,0,0], [1,1,0] and [0,1,1] is [0,1,0]) and sort the BitVectors by their Hamming distance from it. The closest BitVector to this platonic ideal is assumed to be the progenitor. Also, I have observed that using only about 500 bits to perform this calculation resulted in a significant increase in speed with only a 0.001% - 0.002% decrease in accuracy.

Bit-flipping during cloning is a Bernoulli process. Each bit flips (or not) independently of the others. The total number of bits flipped after n clonings is closely approximated by a normal distribution whose mean and variance are determined by the Bernoulli parameters. These distributions tend to be fairly tightly grouped about their means, so the number of flipped bits between two BitVectors is a good indication of the number of clonings between them. However, as the number of mutations increases, both the difference between neighboring means decreases and the variance about those means increases. At 10,000 bits, the expected number of mutations for n and n+1 clonings begins to be closer than 2 standard deviations after about 8 clonings. Beyond this distance, the number of clonings between two BitVectors is not determined with extreme precision by their Hamming distance.

When considering the list of sorted BitVectors, it now is evident that with very high probability, the assumed progenitor's immediate children are at the top of the list and need only 1 additional comparison to establish their parentage. This creates the first layer of what can be considered the BitVector's family tree. BitVectors further away from the assumed progenitor, can descend down through the tree by finding their most probable ancestor at each level (i.e. closest BitVector). Then, only children of their respective most probable ancestors need to be checked to see which is either a parent of the descending BitVector or just the most likely ancestor among its siblings. This process continues until either a parent or a non-parent leaf is found for the descending BitVector.

Not all the BitVectors find their parent on this first pass, because the further one gets from the progenitor, the more likely a most probable ancestor is not an actual ancestor, which causes an error when searching for a parent. So, taking the orphan BitVectors, another search begins 4 generations deeper in the BitVector family tree. This is because 4 generations represents about half the distance over which one can reliably convert the distance between two 10,000 bit BitVectors to the number of generations between them. Then, an attempt is made to find most likely ancestors of the orphan BitVectors from these 4th generation BitVectors and to repeat the insertion process. It starts by sorting all the remaining vectors according to their distance from their most probable ancestors and then allowing them to descend into the tree just as before. This process repeats until all the all the BitVectors have been inserted, or the bottom of the tree is reached.

If the bottom of the tree is reached and there are still orphan BitVectors, then a new search begins, starting at the 3rd generation. If another pass through the tree is needed, the next starts at the 2nd generation, with the final (if necessary) pass starting at the 1st generation. If there are still unsorted nodes, a last ditch, brute force BitVector-by-BitVector comparison search is made at each level of the tree. Since this should only involve a very small number of unsorted nodes, the time cost is should not be prohibitive.

Comparisons are made more efficient by taking advantage of the fact that both the BitVectors and the computer's internal logic are binary in nature. An XOR comparison is made between two BitVectors to identify all the flipped bits in parallel. Next, a look up table is used to count the bits, by looking upon each byte as an integer and looking up the number of 1's in the binary expression of that integer. The sum over all the bytes in this XOR-BitVector gives the Hamming distance between the two compared BitVectors. These comparisons are further sped up by using an externally supplied C extension for manipulating bit arrays. Minor modifications to this extension did not result in any appreciable improvement.

Using these techniques, I average 99+% accuracy for problems of the stated size. When running on a 64-bit version of Ubuntu virtual machine on a Mac with a 2.4GHz processor, my code runs in less than 2 minutes. To improve this code, I need to introduce concurrency. I am currently considering both PyCuda and Stackless Python. After the initial sorting of bit vectors, subsequent sortings may be thought of as occuring in separate BitVector family trees. As such, these sortings can occur simultaneously, which should result in significant speed-up.

Using The Code:

Make certain to replace _bitarray.c with the version included in the tar'ed package (I made an minor change to the code as part of an effort to speed up this code and learn about C-extensions to Python). Then install bitarray following the directions given in the package.

ITA's original problem just gave a file with the bitvectors whose genealogy was desired. However, in order to test, my code, I needed to create test cases. So, to do this, run genTestSets.py as follows:

\$ python genTestSets.py [Options]
-h
Help - List/explains options
-b
Name of input file containing bitvectors
-g
Name of output file recording bitvector parents (ground truth)
-t
Number of bitvector populations to be generated for testing
-v
Number of bits per bitvector
-p
Number of bitvectors in the population
-m
Probability of mutation during cloning

The default settings will generate 5 pairs of test files. The files containing the bitvectors will have a default name of bv-raw-genes_XXXX.dat. The XXXX corresponds to the pair of test files to which the bv-raw-genes_XXXX.dat file belongs. Similarly, default name for the file containing the parent of each bitvector is given in a file bv-parents_XXXX.dat. These files will, by default, contain 1000 bitvectors of 1000 bits each.

With a test case in hand, the code used to test my solution is ITAGenes.py. From the command line, just run:
\$ python ITAGenes.py [Options]

-h
Help - List/explains options
-b
Name of input file containing bitvectors (Default: "bv-raw-genes_0000.dat")
-o
Name of output file recording reconstructed bitvector parents (Default: "recon-bitvectors-parents.data.txt")
-p
Number of bits per bitvector used to to estimate progenitor bitvector (Default: 500)
-c
Number of bits per bitvector used to estimate generations between bitvectors (Default: 1000)
-m
Probability of mutation during cloning (Default: 0.2)
-g
Name of output file recording bitvector parents (ground truth) (Default: "" - No comparison made between reconstructed list of bitvector parents and ground truth list of bitvector parents
In addition to the output file of bitvector parents, ITAGenes.py will also output a trace of how many bitvectors it is inserting starting from a given reference generation and a final tally of how many bitvectors were inserted at each level of the bitvector family tree.