Get a k-mer count from the fermi2 index

occ1.c

/* To compile this program, you need rld0.{h,c} from ropebwt2 or fermi2:
 *
 *   gcc -g -O2 -Wall -o occ1 rld0.c occ1.c
 */

#include <string.h>
#include <stdint.h>
#include <stdio.h>
#include "rld0.h"

static unsigned char seq_nt6_table[128] = { // lookup table
    0, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 1, 5, 2,  5, 5, 5, 3,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  4, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 1, 5, 2,  5, 5, 5, 3,  5, 5, 5, 5,  5, 5, 5, 5,
    5, 5, 5, 5,  4, 5, 5, 5,  5, 5, 5, 5,  5, 5, 5, 5
};

int main(int argc, char *argv[])
{
	int i, j;
	rld_t *e;
	if (argc < 3) {
		fprintf(stderr, "Usage: occ1 <index.fmd> <string1> [string2 [...]]\n");
		return 1;
	}
	e = rld_restore_mmap(argv[1]);
	for (i = 2; i < argc; ++i) {
		char *s = argv[i];
		uint64_t k = 0, l = e->mcnt[0]; // e->mcnt[0]: the total number of symbols
		for (j = strlen(s) - 1; j >= 0; --j) {
			int c = s[j] >= 0 && s[j] <= 127? seq_nt6_table[(int)s[j]] : 5; // nt6 base encoding
			k = e->cnt[c] + rld_rank11(e, k, c); // using rld_rank21() will be faster
			l = e->cnt[c] + rld_rank11(e, l, c);
			if (k == l) break; // no occurrences
		}
		printf("%lld\t%s\n", (long long)(l - k), s);
	}
	rld_destroy(e);
	return 0;
}