C++ - an improvement over FUZxxl solution
I deserve absolutely no credit for the computation method itself, and if no better approach shows up in the mean time, the bounty should go to FUZxxl by right.
#define _CRT_SECURE_NO_WARNINGS // a Microsoft thing about strcpy security issues
#include <vector>
#include <string>
#include <fstream>
#include <iostream>
#include <cstdlib>
#include <ctime>
#include <cstring>
using namespace std;
#include "divsufsort.h"
// graceful exit of sorts
void panic(const char * msg)
{
cerr << msg;
exit(0);
}
// approximative timing of various steps
struct tTimer {
time_t begin;
tTimer() { begin = time(NULL); }
void print(const char * msg)
{
time_t now = time(NULL);
cerr << msg << " in " << now - begin << "s\n";
begin = now;
}
};
// load input pattern
unsigned char * read_sequence (const char * filename, int& len)
{
ifstream file(filename);
if (!file) panic("could not open file");
string str;
std::string line;
while (getline(file, line)) str += line;
unsigned char * res = new unsigned char[str.length() + 1];
len = str.length()+1;
strcpy((char *)res, str.c_str());
return res;
}
#ifdef FUZXXL_METHOD
/*
* Compute for all k the number of k-mers. kmers will contain at index i the
* number of (i + 1) mers. The count is computed as an array of differences,
* where kmers[i] == kmersum[i] - kmersum[i-1] + 1 and then summed up by the
* caller. This algorithm is a little bit unclear, but when you write subsequent
* suffixes of an array on a piece of paper, it's easy to see how and why it
* works.
*/
static void count(const unsigned char *buf, const int *sa, int *kmers, int n)
{
int i, cl, cp;
/* the first item needs special treatment */
/*
kuroi neko: since SA now includes the null string, kmers[0] is indeed 0 instead of 1
*/
// kmers[0]++;
for (i = 1; i < n; i++) {
/* The longest common prefix of the two suffixes */
cl = n - (sa[i - 1] > sa[i] ? sa[i - 1] : sa[i]);
#ifdef ICAP
cl = (cl > ICAP ? ICAP : cl);
#endif
for (cp = 0; cp < cl; cp++)
if (buf[sa[i - 1] + cp] != buf[sa[i] + cp])
break;
/* add new prefixes to the table */
kmers[cp]++;
}
}
#else // Kasai et al. method
// compute kmer cumulative count using Kasai et al. LCP construction algorithm
void compute_kmer_cumulative_sums(const unsigned char * t, const int * sa, int * kmer, int len)
{
// build inverse suffix array
int * isa = new int[len];
for (int i = 0; i != len; i++) isa[sa[i]] = i;
// enumerate common prefix lengths
memset(kmer, 0, len*sizeof(*kmer));
int lcp = 0;
int limit = len - 1;
for (int i = 0; i != limit; i++)
{
int k = isa[i];
int j = sa[k - 1];
while (t[i + lcp] == t[j + lcp]) lcp++;
// lcp now holds the kth longest commpn prefix length, which is just what we need to compute our kmer counts
kmer[lcp]++;
if (lcp > 0) lcp--;
}
delete[] isa;
}
#endif // FUZXXL_METHOD
int main (int argc, char * argv[])
{
if (argc != 2) panic ("missing data file name");
tTimer timer;
int blen;
unsigned char * sequence;
sequence = read_sequence(argv[1], blen);
timer.print("input read");
vector<int>sa;
sa.assign(blen, 0);
if (divsufsort(sequence, &sa[0], blen) != 0) panic("divsufsort failed");
timer.print("suffix table constructed");
vector<int>kmers;
kmers.assign(blen,0);
#ifdef FUZXXL_METHOD
count(sequence, &sa[0], &kmers[0], blen);
timer.print("FUZxxl count done");
#else
compute_kmer_cumulative_sums(sequence, &sa[0], &kmers[0], blen);
timer.print("Kasai count done");
#endif
/* sum up kmers differences */
for (int i = 1; i < blen; i++) kmers[i] += kmers[i - 1] - 1;
timer.print("sum done");
/* human output */
if (blen>10) blen = 10; // output limited to the first few values to avoid cluttering display or saturating disks
for (int i = 0; i != blen; i++) printf("%d ", kmers[i]);
return 0;
}
I simply used Kasai et al. algorithm to compute LCPs in O(n).
The rest is a mere adaptation of FUZxxl code, using more concise C++ features here and there.
I left original computation code to allow comparisons.
Since the slowest processes are SA construction and LCP count, I removed most of the other optimizations to avoid cluttering the code for neglectible gain.
I extended SA table to include the zero-length prefix. That makes LCP computation easier.
I did not provide a length limitation option, the slowest process being now SA computation that cannot be downsized (or at least I don't see how it could be).
I also removed binary output option and limited display to the first 10 values.
I assume this code is just a proof of concept, so no need to clutter displays or saturate disks.
Building the executable
I had to compile the whole project (including the lite version of divsufsort) for x64 to overcome Win32 2Gb allocation limit.
divsufsort code throws a bunch of warnings due to heavy use of ints instead of size_ts, but that will not be an issue for inputs under 2Gb (which would require 26Gb of RAM anyway :D).
Linux build
compile main.cpp and divsufsort.c using the commands:
g++ -c -O3 -fomit-frame-pointer divsufsort.c
g++ -O3 main.cpp -o kmers divsufsort.o
I assume the regular divsufsort library should work fine on native Linux, as long as you can allocate a bit more than 3Gb.
Performances
The Kasai algorithm requires the inverse SA table, which eats up 4 more bytes per character for a total of 13 (instead of 9 with FUZxxl method).
Memory consumption for reference input is thus above 3Gb.
On the other hand, computation time is dramatically improved, and the whole algorithm is now in O(n):
input read in 5s
suffix table constructed in 37s
FUZxxl count done in 389s
Kasai count done in 27s
14 92 520 2923 15714 71330 265861 890895 2482912 5509765 (etc.)
Further improvements
SA construction is now the slowest process.
Some bits of the divsufsort algorithm are meant to be parallelized with whatever builtin feature of a compiler unknown to me, but if necessary the code should be easy to adapt to more classic multithreading (à la C++11, for instance).
The lib has also a truckload of parameters, including various bucket sizes and the number of distinct symbols in the input string. I only had a cursory look at them, but I suspect compressing the alphabet might be worth a try if your strings are endless ACTG permutations (and you're desperate for performances).
There exist some parallelizable methods for computing LCP from SA too, but since the code should run under one minute on a processor slightly more powerful than my puny [email protected] and the whole algorithm is in O(n), I doubt this would be worth the effort.
J, a naïve solution would simply be `[:~.]` but guess that won't cut it. \$\endgroup\$