bwa_perf/fmt_index.cpp

#include <iostream>
#include <stdint.h>
#include <stdlib.h>
#include <vector>
#include <sys/time.h>
#include <string>
#include <stdio.h>
#include <algorithm>
#include <string.h>

#include "util.h"
#include "fmt_index.h"
using namespace std;

const char BASE[4] = {'A', 'C', 'G', 'T'};

// 求反向互补序列
string calc_reverse_seq(string &seq)
{
	string rseq(seq.size(), '0');
	for (int i = 0; i < seq.size(); ++i)
	{
		if (seq[i] == 'A')
			rseq[i] = 'T';
		else if (seq[i] == 'C')
			rseq[i] = 'G';
		else if (seq[i] == 'G')
			rseq[i] = 'C';
		else if (seq[i] == 'T')
			rseq[i] = 'A';
	}
	std::reverse(rseq.begin(), rseq.end());
	return rseq;
}

// 打印32位整型数据中包含的pre-bwt：bwt
void print_base_uint32(uint32_t p)
{
	for (int i = 30; i > 0; i -= 4)
	{
		int b1 = p >> i & 3;
		int b2 = p >> (i - 2) & 3;
		cout << BASE[b1] << BASE[b2] << endl;
	}
}

// 随机生成长度为len的序列
string generate_rand_seq(int len)
{
	string seq(len, 'A');
	for (int i = 0; i < len; ++i)
	{
		seq[i] = BASE[rand() % 4];
	}
	return seq;
}
// 创建bwt矩阵
void create_bwt_mtx(string &seq)
{
	bwtint_t seq_len = seq.size() + 1;
	string sarr[seq_len];
	sarr[0] = seq + '$';
	for (int i = 1; i < seq_len; ++i)
	{
		sarr[i] = sarr[0].substr(i) + sarr[0].substr(0, i);
	}
	std::sort(sarr, sarr + seq_len);

	// print bwt matrix
	// for (int i = 0; i < seq_len; ++i)
	//{
	//	// cout << i << ' ' << sarr[i] << endl;
	//	cout << sarr[i] << endl;
	//}
	// cout << "bwt string" << endl;
	// for (int i = 0; i < seq_len; ++i)
	// {
	// 	cout << sarr[i].back();
	// }
	// cout << endl;
	// cout << "pre bwt string" << endl;
	// for (int i = 0; i < seq_len; ++i)
	// {
	// 	cout << sarr[i][seq_len - 2];
	// }
	// cout << endl;
}

// fmt-index的count table，4对应着bwt碱基的累积量，0,1,2,3分别对应着bwt是A,C,G,T，pre-bwt的累积量
void fmt_gen_cnt_table(FMTIndex *fmt)
{
	int i, j, k;
	for (i = 0; i != 256; ++i) // 遍历单个字节的各种情况
	{
		uint32_t x = 0;
		for (j = 0; j != 4; ++j) // 一个字节有8位，每个碱基是2位，所以一个字节包含4个碱基，从右向左数，第一个和第三个碱基数据bwt，第二个和第四个是对应的pre-bwt
			x |= (((i & 3) == j) + ((i >> 4 & 3) == j)) << (j << 3); // 高位存pre-bwt，挨着存bwt，一一对应
		fmt->cnt_table[4][i] = x; // 保存单个字节中bwt碱基个数，每8位对应一个碱基的个数，从左到右依次是TGCA

		for (k = 0; k < 4;++k) // bwt碱基
		{
			x = 0; // for [A,C,G,T][A,C,G,T]
			for (j = 0; j != 4; ++j) // pre-bwt碱基
				x |= (((i >> 6 & 3) == j && (i >> 4 & 3) == k) + ((i >> 2 & 3) == j && (i & 3) == k)) << (j << 3); 
			fmt->cnt_table[k][i] = x;
		}
	}
}

void dump_fmt(const char *fn, const FMTIndex *fmt)
{
	FILE *fp;
	fp = xopen(fn, "wb");
	err_fwrite(&fmt->primary, sizeof(bwtint_t), 1, fp);
	err_fwrite(&fmt->sec_primary, sizeof(bwtint_t), 1, fp);
	err_fwrite(&fmt->sec_bcp, sizeof(uint8_t), 1, fp);
	err_fwrite(&fmt->first_base, sizeof(uint8_t), 1, fp);
	err_fwrite(&fmt->last_base, sizeof(uint8_t), 1, fp);
	err_fwrite(fmt->L2 + 1, sizeof(bwtint_t), 4, fp);
	err_fwrite(fmt->bwt, 4, fmt->bwt_size, fp);
	err_fflush(fp);
	err_fclose(fp);
}

FMTIndex *restore_fmt(const char *fn)
{
	FMTIndex *fmt;
	fmt = (FMTIndex *)calloc(1, sizeof(FMTIndex));
	FILE *fp = fopen(fn, "rb");

	fseek(fp, 0, SEEK_END);
	fmt->bwt_size = (ftell(fp) - sizeof(bwtint_t) * 6 - 3) >> 2; // 以32位word为单位计算的size
	fmt->bwt = (uint32_t *)calloc(fmt->bwt_size, 4);
	fseek(fp, 0, SEEK_SET);
	fread(&fmt->primary, sizeof(bwtint_t), 1, fp);
	fread(&fmt->sec_primary, sizeof(bwtint_t), 1, fp);
	fread(&fmt->sec_bcp, sizeof(uint8_t), 1, fp);
	fread(&fmt->first_base, sizeof(uint8_t), 1, fp);
	fread(&fmt->last_base, sizeof(uint8_t), 1, fp);
	fread(fmt->L2 + 1, sizeof(bwtint_t), 4, fp);
	fread_fix(fp, fmt->bwt_size << 2, fmt->bwt);
	fmt->seq_len = fmt->L2[4];
	fclose(fp);
	fmt_gen_cnt_table(fmt); // 字节所能表示的各种碱基组合中，各个碱基的累积数量
	return fmt;
}

// 根据interval-bwt创建fmt-index
FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
{
	FILE *fmt_out = fopen("fmt.txt", "w");
	FMTIndex *fmt = (FMTIndex *)calloc(1, sizeof(FMTIndex));
	fmt_gen_cnt_table(fmt);

	bwtint_t i, j, k, m, n, n_occ, cnt[4], cnt2[4];
	uint32_t c[4], c2[16]; /*c用来保存原来的bwt碱基串的累积值，c2用来保存pre-bwt和bwt碱基对的累计值，如AA..TT*/
	uint32_t *buf; /* 计算之后变成fmt结构中bwt部分 */

	fmt->seq_len = bwt->seq_len; // bwt碱基序列的长度，不包含$字符，也就是该长度比bwt matrix长度少1
	for (i = 0; i < 5; ++i)
		fmt->L2[i] = bwt->L2[i]; // 每个碱基的总累积值
	fmt->primary = bwt->primary; // $在末尾的行，在bwt matrix行中的排序位置

	n_occ = (bwt->seq_len + OCC_INTERVAL - 1) / OCC_INTERVAL + 1; // check point 个数
	fmt->bwt_size = (fmt->seq_len * 2  + 15) >> 4; // 要保存最后两列碱基
	fmt->bwt_size += n_occ * 20; // A,C,G,T和AA,AC.....TG,TT共20个
	buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存

	c[0] = c[1] = c[2] = c[3] = 0;
	// 首行的c2，应该是对应的ACGT对应的行，减去1的occ
	for (i = 0; i < 4; ++i)
	{
		bwtint_t before_first_line = fmt->L2[i];
		bwt_occ4(bwt, before_first_line, cnt);
		for (j = i * 4, k = 0; k < 4; ++j, ++k)
			c2[j] = cnt[k];
	}
	// k表示buf存储的偏移量
	for (i = k = 0; i < bwt->seq_len; ++i)
	{
		// 记录occ
		if (i % OCC_INTERVAL == 0)
		{
			memcpy(buf + k, c, sizeof(uint32_t) * 4); // bwt str中各个碱基的occ
			k += 4;
			memcpy(buf + k, c2, sizeof(uint32_t) * 16); // pre-bwt:bwt碱基对的occ
			k += 16;
		}
		// 每个32位整数保存8个倒数第二列碱基（pre-bwt）和8个倒数第一列(bwt)碱基
		if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16
		{
			uint32_t pre_bwt_16_seq = 0; // 16个pre-bwt碱基串
			uint32_t *bwt_addr = bwt_occ_intv(bwt, i) + 8; // bwt字符串i对应的基准行，因为原始的bwt-cp（check point）包含由4个uint64_t(8个uint32_t)组成的occ信息
			int offset = (i % OCC_INTERVAL) / 16;		   // 每OCC_INTERVAL个碱基共享同一个基准地址，每16个碱基共用一个uint32整型，因此需要偏移量来获取当前碱基串的首地址
			uint32_t bwt_16_seq = *(bwt_addr + offset);    // 待处理的当前16个碱基串的首地址
			for (j = 0; j < 16; ++j) // 对于bwt碱基串，一个一个碱基分别处理
			{
				bwtint_t cur_str_line = i + j;  // 当前碱基在bwt str中的行排序
				if (cur_str_line < bwt->seq_len) // 当前碱基行不应超出bwt str总碱基长度（bwt str长度比bwt matrix长度少1，因为bwt str不包含$）
				{
					uint8_t bwt_base = bwt_B0(bwt, cur_str_line); // 对应行的bwt的碱基
					// 先求出该碱基对应在第一列的行（对应的bwt matrix行）
					bwtint_t cur_mtx_line = cur_str_line;
					if (cur_str_line >= bwt->primary) // 因为bwt序列里除去了$符号，所以，超过$所在行之后，对应的seq位置应该加一，才是真正对应bwt matrix的行
						cur_mtx_line += 1;
					bwt_occ4(bwt, cur_mtx_line, cnt);	  // 获取原来bwt-checkpoint中的occ值	  
					for (m=0; m<4; ++m)
						c[m] = (uint32_t)cnt[m]; // 碱基m在cur_bwt_mtx_line(包含)之前的累积值，直接拷贝原bwt中的occ即可

					cnt[bwt_base] -= 1; // 得到cur_bwt_mtx_line(不包含)之前的累积量，即bwt_occ4(bwt, cur_bwt_mtx_line-1, cnt)
					bwtint_t bwt_base_mtx_line = bwt->L2[bwt_base] + 1 + cnt[bwt_base]; // bwt_base对应的bwt matrix行（LF变换）

					bwt_occ4(bwt, bwt_base_mtx_line, cnt2); // 计算bwt_base_mtx_line之前的occ
					for (n = 0; n < 4; ++n)
					{
						int c2_idx = bwt_base << 2 | n; // bwt base放在前边
						c2[c2_idx] = (uint32_t)cnt2[n]; // pre-bwt:bwt 碱基对的累计值
					}
					bwtint_t bwt_base_str_line = bwt_base_mtx_line; // bwt-str中对应的行排序
					if (bwt_base_str_line >= bwt->primary)			// base_line表示在bwt str中的位置，所以超出$为最尾所在行之后，要减掉1
						bwt_base_str_line -= 1;						// bwt碱基序列行（不包含$）
					uint32_t pre_bwt_base = bwt_B0(bwt, bwt_base_str_line); // bwt列碱基对应的前一个碱基pre-bwt
					// 此时，bwt_base对应的bwt matrix首行，是$排在最尾的行，说明bwt_base就是序列的第一个碱基，
					// 此时计算出来的pre_bwt_base就是primary前一行的bwt base，以此来代替$字符，在后续的计算过程中需要考虑
					if (bwt_base_mtx_line == bwt->primary)
					{
						// 计算sec_bcp
						fmt->sec_bcp = pre_bwt_base << 2 | bwt_base; // 因为把$当成A处理了
						fmt->sec_primary = cur_mtx_line;		 // pre-bwt base为$的行排序（bwt-matrix行）
						fmt->first_base = bwt_base; // 原始序列第一个碱基
						fmt->last_base = pre_bwt_base; // 计算后替代$字符的碱基（应该是primary行上边一行对应的bwt base）
					}
					// 暂存 pre-bwt碱基序列
					pre_bwt_16_seq = pre_bwt_16_seq | (pre_bwt_base << (15-j)*2); // 序列靠前的碱基排在uint32_t数据中的高位

					// 输出调试信息
					// cout << "mtx line: " << cur_mtx_line << ' ' << c[0] << ' ' << c[1] << ' ' << c[2] << ' ' << c[3] << ' ';
					// for (m = 0; m < 16; ++m)
					// 	cout << c2[m] << ' ';
					// cout << endl;
				}
				else
					break;
			}
			// 保存bwt和pre_bwt
			uint32_t pre_and_bwt_seq = 0;
			for (m = 0; m < 8; ++m)
			{
				int lshift_bit = 30 - 2 * m;
				pre_and_bwt_seq |= (((pre_bwt_16_seq & (3 << lshift_bit)) >> (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) >> ((m * 2) + 2)));
			}
			buf[k++] = pre_and_bwt_seq;
			if (j > 8)
			{
				pre_and_bwt_seq = 0;
				for (m = 8; m > 0; --m)
				{
					int lshift_bit = 2 * m - 2;
					pre_and_bwt_seq |= (((pre_bwt_16_seq & (3 << lshift_bit)) << (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) << (m * 2 - 2)));
				}
				buf[k++] = pre_and_bwt_seq;
			}
		}
	}
	// the last element
	memcpy(buf + k, c, sizeof(uint32_t) * 4);
	k += 4;
	memcpy(buf + k, c2, sizeof(uint32_t) * 16);
	k += 16;
	// cout << "n occ: " << n_occ << endl;
	// cout << "size: " << k << '\t' << fmt->bwt_size << endl;
	xassert(k == fmt->bwt_size, "inconsistent bwt_size");
	// update fmt
	fmt->bwt = buf; 
	return fmt;
}

// 扩展两个个碱基，计算bwt base为b的pre-bwt str中各个碱基的occ
void fmt_e2_occ4(const FMTIndex *fmt, bwtint_t k, int b, uint32_t cnt1[4], uint32_t cnt2[4])
{
	uint32_t x1, x2;
	uint32_t *p, tmp, *end;
	bwtint_t bwt_k_line = k, bwt_k_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT;
	if (k == (bwtint_t)(-1))
	{
		p = fmt->bwt + 4 + b * 4;
		memset(cnt1, 0, 4 * sizeof(uint32_t));
		memcpy(cnt2, p, 4 * sizeof(uint32_t));
		return;
	}
	k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1）
	p = fmt_occ_intv(fmt, k);
	memcpy(cnt1, p, 4 * sizeof(uint32_t));
	memcpy(cnt2, p + 4 + b * 4, 4 * sizeof(uint32_t));
	p += 20;											// 该地址是bwt和pre_bwt字符串数据的首地址
	end = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3)); // this is the end point of the following loop

	for (x1 = 0, x2 = 0; p < end; ++p)
	{
		x1 += __fmt_occ_e2_aux4(fmt, 4, *p);
		x2 += __fmt_occ_e2_aux4(fmt, b, *p);
	}
	tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
	x1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7);
	x2 += __fmt_occ_e2_aux4(fmt, b, tmp);
	if (b == 0)
		x2 -= ~k & 7;
	// 如果跨过了second_primary,那么可能需要减掉一次累积值
	if (b == fmt->first_base && bwt_k_base_line < fmt->sec_primary && bwt_k_line >= fmt->sec_primary)
	{
		x2 -= 1 << (fmt->last_base << 3);
	}

	cnt1[0] += x1 & 0xff;
	cnt1[1] += x1 >> 8 & 0xff;
	cnt1[2] += x1 >> 16 & 0xff;
	cnt1[3] += x1 >> 24;
	cnt2[0] += x2 & 0xff;
	cnt2[1] += x2 >> 8 & 0xff;
	cnt2[2] += x2 >> 16 & 0xff;
	cnt2[3] += x2 >> 24;

	// cout << "fmt-occ: " << k << '\t' << cnt1[0] << '\t' << cnt1[1] << '\t' << cnt1[2] << '\t' << cnt1[3] << endl;
	// cout << "fmt-occ-2: " << k << '\t' << cnt2[0] << '\t' << cnt2[1] << '\t' << cnt2[2] << '\t' << cnt2[3] << endl;
	// cout << "bwt_k_base_line: " << bwt_k_base_line << endl;
	// cout << "bwt_k_line: " << bwt_k_line << endl;
	// cout << "sec_primary: " << fmt->sec_primary << endl;
}

// 对k行和l行同时计算occ，如果k和l落在同一个interval区间，可以减少一些计算量和访存
void fmt_e2_2occ4(const FMTIndex *fmt, bwtint_t k, bwtint_t l, int b,
			   uint32_t cntk1[4], uint32_t cntk2[4], uint32_t cntl1[4], uint32_t cntl2[4])
{
	// fmt_e2_occ4(fmt, k, b, cntk1, cntk2);
	// fmt_e2_occ4(fmt, l, b, cntl1, cntl2);
	// return;
	bwtint_t _k, _l;
	_k = k - (k >= fmt->primary); // 换算成了seq的行
	_l = l - (l >= fmt->primary);
	if (_l >> OCC_INTV_SHIFT != _k >> OCC_INTV_SHIFT || k == (bwtint_t)(-1) || l == (bwtint_t)(-1))
	{
		fmt_e2_occ4(fmt, k, b, cntk1, cntk2);
		fmt_e2_occ4(fmt, l, b, cntl1, cntl2);
	}
	else
	{
		uint32_t x1, x2, y1, y2;
		uint32_t *p, tmp, *ek, *el;
		bwtint_t bwt_k_line = k, bwt_l_line = l, bwt_base_line = k >> OCC_INTV_SHIFT << OCC_INTV_SHIFT;

		k -= (k >= fmt->primary); // because $ is not in bwt
		l -= (l >= fmt->primary);
		p = fmt_occ_intv(fmt, k);
		memcpy(cntk1, p, 4 * sizeof(uint32_t));
		memcpy(cntk2, p + 4 + b * 4, 4 * sizeof(uint32_t));
		memcpy(cntl1, cntk1, 4 * sizeof(uint32_t));
		memcpy(cntl2, cntk2, 4 * sizeof(uint32_t));
		p += 20;
		// prepare cntk[]
		ek = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3));
		el = p + ((l >> 3) - ((l & ~OCC_INTV_MASK) >> 3));
		for (x1 = 0, x2 = 0; p < ek; ++p)
		{
			x1 += __fmt_occ_e2_aux4(fmt, 4, *p);
			x2 += __fmt_occ_e2_aux4(fmt, b, *p);
		}
		y1 = x1;
		y2 = x2;
		tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
		x1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7);
		x2 += __fmt_occ_e2_aux4(fmt, b, tmp);
		if (b == 0)
			x2 -= ~k & 7;
		for (; p < el; ++p)
		{
			y1 += __fmt_occ_e2_aux4(fmt, 4, *p);
			y2 += __fmt_occ_e2_aux4(fmt, b, *p);
		}
		tmp = *p & ~((1U << ((~l & 7) << 2)) - 1);
		y1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~l & 7);
		y2 += __fmt_occ_e2_aux4(fmt, b, tmp);
		if (b == 0)
			y2 -= ~l & 7;
		// 如果跨过了second_primary,那么可能需要减掉一次累积值
		if (b == fmt->first_base && bwt_base_line < fmt->sec_primary)
		{
			if (bwt_k_line >= fmt->sec_primary)
				x2 -= 1 << (fmt->last_base << 3);
			if (bwt_l_line >= fmt->sec_primary)
				y2 -= 1 << (fmt->last_base << 3);
		}
		cntk1[0] += x1 & 0xff;
		cntk1[1] += x1 >> 8 & 0xff;
		cntk1[2] += x1 >> 16 & 0xff;
		cntk1[3] += x1 >> 24;
		cntk2[0] += x2 & 0xff;
		cntk2[1] += x2 >> 8 & 0xff;
		cntk2[2] += x2 >> 16 & 0xff;
		cntk2[3] += x2 >> 24;

		cntl1[0] += y1 & 0xff;
		cntl1[1] += y1 >> 8 & 0xff;
		cntl1[2] += y1 >> 16 & 0xff;
		cntl1[3] += y1 >> 24;
		cntl2[0] += y2 & 0xff;
		cntl2[1] += y2 >> 8 & 0xff;
		cntl2[2] += y2 >> 16 & 0xff;
		cntl2[3] += y2 >> 24;

		// cout << "fmt-occ: " << k << '\t' << cntk1[0] << '\t' << cntk1[1] << '\t' << cntk1[2] << '\t' << cntk1[3] << endl;
		// cout << "fmt-occ-2: " << k << '\t' << cntk2[0] << '\t' << cntk2[1] << '\t' << cntk2[2] << '\t' << cntk2[3] << endl;
		// cout << "fmt-occ: " << l << '\t' << cntl1[0] << '\t' << cntl1[1] << '\t' << cntl1[2] << '\t' << cntl1[3] << endl;
		// cout << "fmt-occ-2: " << l << '\t' << cntl2[0] << '\t' << cntl2[1] << '\t' << cntl2[2] << '\t' << cntl2[3] << endl;
	}
}

// 扩展一个碱基，计算bwt str中各个碱基的occ
void fmt_e1_occ4(const FMTIndex *fmt, bwtint_t k, uint32_t cnt[4])
{
	uint32_t x;
	uint32_t *p, tmp, *end;
	if (k == (bwtint_t)(-1))
	{
		memset(cnt, 0, 4 * sizeof(uint32_t));
		return;
	}
	k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1）
	p = fmt_occ_intv(fmt, k);
	memcpy(cnt, p, 4 * sizeof(uint32_t));
	p += 20;											// 该地址是bwt和pre_bwt字符串数据的首地址
	end = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3)); // this is the end point of the following loop

	for (x = 0; p < end; ++p)
	{
		x += __fmt_occ_e2_aux4(fmt, 4, *p);
	}
	tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
	x += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7);
	cnt[0] += x & 0xff;
	cnt[1] += x >> 8 & 0xff;
	cnt[2] += x >> 16 & 0xff;
	cnt[3] += x >> 24;
}

// 对k行和l行同时计算bwt str的occ，如果k和l落在同一个interval区间，可以减少一些计算量和访存
void fmt_e1_2occ4(const FMTIndex *fmt, bwtint_t k, bwtint_t l, uint32_t cntk[4], uint32_t cntl[4])
{
	bwtint_t _k, _l;
	_k = k - (k >= fmt->primary); // 换算成了seq的行
	_l = l - (l >= fmt->primary);
	if (_l >> OCC_INTV_SHIFT != _k >> OCC_INTV_SHIFT || k == (bwtint_t)(-1) || l == (bwtint_t)(-1))
	{
		fmt_e1_occ4(fmt, k, cntk);
		fmt_e1_occ4(fmt, l, cntl);
	}
	else
	{
		uint32_t x1, y1;
		uint32_t *p, tmp, *endk, *endl;

		k -= (k >= fmt->primary); // because $ is not in bwt
		l -= (l >= fmt->primary);
		p = fmt_occ_intv(fmt, k);
		memcpy(cntk, p, 4 * sizeof(uint32_t));
		memcpy(cntl, p, 4 * sizeof(uint32_t));
		p += 20;
		// prepare cntk[]
		endk = p + ((k >> 3) - ((k & ~OCC_INTV_MASK) >> 3));
		endl = p + ((l >> 3) - ((l & ~OCC_INTV_MASK) >> 3));
		for (x1 = 0; p < endk; ++p)
		{
			x1 += __fmt_occ_e2_aux4(fmt, 4, *p);
		}
		y1 = x1;
		tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
		x1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7);
		for (; p < endl; ++p)
		{
			y1 += __fmt_occ_e2_aux4(fmt, 4, *p);
		}
		tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
		y1 += __fmt_occ_e2_aux4(fmt, 4, tmp) - (~k & 7);

		cntk[0] += x1 & 0xff;
		cntk[1] += x1 >> 8 & 0xff;
		cntk[2] += x1 >> 16 & 0xff;
		cntk[3] += x1 >> 24;

		cntl[0] += y1 & 0xff;
		cntl[1] += y1 >> 8 & 0xff;
		cntl[2] += y1 >> 16 & 0xff;
		cntl[3] += y1 >> 24;
	}
}

// 扩展一个碱基
void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int b1)
{
	uint32_t tk[4], tl[4];
	int i;

	fmt_e1_occ4(fmt, ik->x[!is_back] - 1, tk);
	fmt_e1_occ4(fmt, ik->x[!is_back] - 1 + ik->x[2], tl);
	fmt_e1_2occ4(fmt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl);
	for (i = 0; i != 4; ++i)
	{
		ok[i].x[!is_back] = fmt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
		ok[i].x[2] = tl[i] - tk[i];				 // 间隔
	}
	ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= fmt->primary && ik->x[!is_back] + ik->x[2] - 1 >= fmt->primary);
	for (i = 2; i >= b1; --i)
		ok[i].x[is_back] = ok[i + 1].x[is_back] + ok[i + 1].x[2];
	*ik = ok[b1];
}

// 扩展两个碱基
void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int b1, int b2)
{
	uint32_t tk1[4], tl1[4], tk2[4], tl2[4];
	int i;

	// fmt_e2_occ4(fmt, ik->x[!is_back] - 1, b1, tk1, tk2);
	// fmt_e2_occ4(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, tl1, tl2);
	// tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
	fmt_e2_2occ4(fmt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], b1, tk1, tk2, tl1, tl2);

	// cout << "k: " << tk1[0] << '\t' << tk1[1] << '\t' << tk1[2] << '\t' << tk1[3] << endl;
	// cout << "l: " << tl1[0] << '\t' << tl1[1] << '\t' << tl1[2] << '\t' << tl1[3] << endl;
	// cout << "k: " << tk2[0] << '\t' << tk2[1] << '\t' << tk2[2] << '\t' << tk2[3] << endl;
	// cout << "l: " << tl2[0] << '\t' << tl2[1] << '\t' << tl2[2] << '\t' << tl2[3] << endl;

	// 这里是反向扩展
	for (i = 0; i != 4; ++i)
	{
		ok[i].x[!is_back] = fmt->L2[i] + 1 + tk2[i]; // 起始行位置，互补链
		ok[i].x[2] = tl2[i] - tk2[i];  				// 间隔
	}
	// 因为计算的是互补碱基，所以3对应着0,2对应1，下边是正向扩展

	ok[3].x[is_back] = ik->x[is_back] + (ik->x[!is_back] <= fmt->primary && ik->x[!is_back] + ik->x[2] - 1 >= fmt->primary);
	ok[2].x[is_back] = ok[3].x[is_back] + tl1[3] - tk1[3];
	ok[1].x[is_back] = ok[2].x[is_back] + tl1[2] - tk1[2];
	ok[0].x[is_back] = ok[1].x[is_back] + tl1[1] - tk1[1];

	cout << "fmt-d: " << BASE[b1] << '\t' << ok[b1].x[is_back] << '\t' << ok[b1].x[2] << endl;

	ok[3].x[is_back] = ok[b1].x[is_back] + (ok[b1].x[!is_back] <= fmt->primary && ok[b1].x[!is_back] + ok[b1].x[2] - 1 >= fmt->primary);
	ok[2].x[is_back] = ok[3].x[is_back] + ok[3].x[2];
	ok[1].x[is_back] = ok[2].x[is_back] + ok[2].x[2];
	ok[0].x[is_back] = ok[1].x[is_back] + ok[1].x[2];

	*ik = ok[b2];
}
// 利用fmt搜索seed，完整搜索，只需要单向搜索
void fmt_search(FMTIndex *fmt, const string &q)
{
	bwtintv_t ik, ok[4];
	int i, j, c1, c2, x = 0;

	fmt_set_intv(fmt, bval(q[x]), ik);
	ik.info = x + 1;
	cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;

	for (i = x + 1; i + 1 < q.size(); i += 2)
	{
		if (bval(q[i]) < 4 && bval(q[i + 1]) < 4)
		{

			c1 = cbval(q[i]);
			c2 = cbval(q[i + 1]);

			fmt_extend2(fmt, &ik, ok, 0, c1, c2);
			ik.info = i + 1;
			cout << "fmt  : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
		}
		else
		{
			break;
		}
	}
	if (i < q.size() && bval(q[i]) < 4)
	{ // 最后一次扩展
		c1 = cbval(q[i]);
		fmt_extend1(fmt, &ik, ok, 0, c1);
		ik.info = i + 1;
		cout << "fmt  : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
	}
}

int main_fmtidx(int argc, char **argv)
{
	// string seq("ACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTA");
	string seq("ACCCT");
	string rseq = calc_reverse_seq(seq);
	seq = seq + rseq;
	//create_bwt_mtx(seq);
	//cout << seq << endl;

	bwt_t *bwt = restore_bwt(argv[1]); // 读取bwt原始字符串（带ACGT总的累积量）
	// create_interval_occ_bwt(bwt); // 根据bwt字符串创建包含interval occ的bwt（128碱基+ACGT累积量）
	cout << "L2: " << bwt->L2[0] << '\t' << bwt->L2[1] << '\t' << bwt->L2[2] << '\t'
		 << bwt->L2[3] << '\t' << bwt->L2[4] << endl;

	string s = "AACCCTAA";

	srand(time(NULL));
	s = generate_rand_seq(10);
	cout << "seq: " << s << endl;
	// s = "TTC";

	bwt_search(bwt, s);
	bwt_search2(bwt, s);

	// for (int i = 0; i < 120; ++i)
	// {
	// 	cout << i << '\t' << bwt_B0(bwt, i) << endl;
	// }
	// TGGGAT
	FMTIndex *fmt = create_fmt_from_bwt(bwt);
	dump_fmt("ref.fmt", fmt);
	// FMTIndex *fmt = restore_fmt("tiny.fmt");

	fmt_search(fmt, s);
	// cout << bwt->bwt_size << endl;
	// cout << bwt->seq_len << endl;
	cout << "sec_: " << fmt->sec_bcp << '\t' << fmt->sec_primary << endl;

	uint8_t b8 = 2 << 4 | 2;
	cout << "AGAG: " << fmt->cnt_table[2][b8] << endl;
	cout << (((b8 >> 6) == 0 && (b8 >> 4 & 3) == 2) + ((b8 >> 2 & 3) == 0 && (b8 & 3) == 2)) << endl;
	return 0;
}