From: "andrew cooke" <andrew@...>
Date: Sat, 20 Sep 2008 18:54:39 -0400 (CLT)
I have just finished implementing the main core of the algorithm outlined here - http://www.acooke.org/cute/Identifyin0.html - directly in Java and it runs in about 8 seconds! There were two main problems. First, inferring how Postgres did an efficient search and, second, implementing that without using too much memory (my first attempt exhausted the heap so I now have a slight tradeoff, which uses a sort to avoid creating more memory structures and so adds a log term to the big-O). It's easiest to describe both together, by outlining the final solution, but in practice I the development had two distinct steps. So, as in the prototype code, I generate candidate pairs by matching small fragments of the DNA. More exactly: I take 25 fragments, each 8 bits, from each individual and I categorise two individuals as a candidate pair if they have at least 3 fragments in common. So, in psuedocode, I do the following: generate a table of fragments[individual_idx][fragment_idx] generate a table of counts[individual1_idx][individual2_idx] = 0 for each column of fragments in turn: sort the table column containing the fragments; scan the sorted column: for all fragments with the same value: increment the counts associated with the pairs of individuals that share that fragment value; if any count == 3: if the "bit distance" between the pair is < 3000: add the pair for that count to the graph; And I need to repeat this 6 times with different sets of fragments (the number of identified pairs after each set is 8116, 9623, 9935, 9988, 9998, 9999). Instead of sorting each column of fragments the scan could be direct, but you would need to have a separate memory structure to record which individuals were associated with which values (for this amount of data I suspect the log pays for itself in the simplification (reduced constant cost) that the sorted data introduces). Also, Java has no direct support for sorting bytes (the fragments) with keys. I could have wrapped everything in objects, but it was more compact (and probably faster) to bit-pack the DNA fragment and the individual index together in a single integer (obviously the DNA has to occupy the more significant bits for the sorting to give the corrected order). I am going to look for a graph library now to finish this off. 8 seconds is pretty good. When I started out I was looking at many hours; even the optimized Python/SQL code took 30 min... Andrew PS My initial attempt at searching the hashes was to do a depth first search trying to find common fragments for each pair in turn. While this would have fitted well within a constraint programming framework (see my posts here over the last week or two when I was looking at Choco and Gecode) it was, in retrospect, completely stupid - a huge amount of time is spent exhaustively searching irrelevant pairs. The direct scan described above is much more efficient, but it's not yet clear to me how the two approaches are related. Is there some way in which the direct scan with counting is a dual of the search? Or does some kind of optimisation of the search eventually reduce it to the scan? I don't see how either of those pan out, but haven't looked at the CP techniques in any detail yet.
Core Routine
From: "andrew cooke" <andrew@...>
Date: Sat, 20 Sep 2008 19:07:14 -0400 (CLT)
public int search()
{
byte[] counts = new byte[GenomePair.hashSize(population.size())];
int[] scratch = new int[population.size()];
// for each fragment in turn:
for (int column = 0; column < nHashes; ++column) {
// pack into an integer
for (int row = 0; row < population.size(); ++row) {
scratch[row] = pack(hashes[row][column], row);
}
// group individuals with the same hash are together
Arrays.sort(scratch);
// for each group
for (int row = 0; row < population.size();) {
// get the hash for the group
byte hash = unpackHash(scratch[row]);
Set<Integer> allMatching = new HashSet<Integer>();
// note the first individual
allMatching.add(unpackRow(scratch[row]));
// for each additional individual
while (++row < population.size() &&
unpackHash(scratch[row]) == hash) {
int higher = unpackRow(scratch[row]);
// for each pair
for (int lower: allMatching) {
GenomePair pair = new GenomePair(lower, higher);
// if we have sufficient hits, check the distance
if (++counts[pair.hashCode()] == nMatches
&& population.connected(pair, cutoff)) {
graph.add(pair);
}
}
// extend the current set so that we generate all pairs
allMatching.add(higher);
}
}
}
// this should tend to to population.size()-1
return graph.size();
}
Perfect Hash
From: "andrew cooke" <andrew@...>
Date: Sat, 20 Sep 2008 19:10:22 -0400 (CLT)
I should explain that I am abusing GenomePair.hashCode() - the implementation returns a continguous index from 0 over all possible pairs. So all pairs are distinct and there are no gaps. The total number of values is given by hashSize(). At some point I'll change the name. Originally I was using HashMaps of these... Andrew
Same Results
From: "andrew cooke" <andrew@...>
Date: Sun, 21 Sep 2008 20:05:00 -0400 (CLT)
I added the final graph code (using JGraphT, which seems quite capable) and the results are, as expected, identical to the earlier work. I also tried some variations on the numbers of matches and hashes (but not the fragment size, which is hard coded at 8 bits (ie bytes) in this version) - the code is much more stable than the Python/SQL implementation to these changes (I now suspect Postgres was switching algorithms depending on predicted memory usage), and the values chosen aren't particularly critical. I'm considering sending it off to the company that posted the problem, but they only accept submissions that are employment applications, so it seems a bit silly (I'm not looking for a job, and won't move to Boston...). Andrew