zzh
/
bwa_perf
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

mid interval支持8和16，代码中含有注释掉的调试信息，发现并解决了一个数据溢出的bug

This commit is contained in:
zzh 2024-02-02 03:59:35 +08:00
parent e8ceb3ff58
commit 5f18167703
4 changed files with 87 additions and 27 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
Makefile
View File

@ -1,7 +1,7 @@
CC= g++ CC= g++
NOWARN= -Wno-unused-result -Wno-unused-function NOWARN= -Wno-unused-result -Wno-unused-function
CFLAGS= #-g -Wall $(NOWARN) #-O2 CFLAGS= #-g -Wall $(NOWARN) #-O2
CPPFLAGS= -g -Wall $(NOWARN) -O3 CPPFLAGS= -g -Wall $(NOWARN) -O2
WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS WRAP_MALLOC=-DUSE_MALLOC_WRAPPERS
SHOW_PERF= -DSHOW_PERF SHOW_PERF= -DSHOW_PERF
AR= ar AR= ar

24
bwt.cpp
View File

@ -10,6 +10,8 @@
using namespace std; using namespace std;
const static char BASE[4] = {'A', 'C', 'G', 'T'};
// 打印32位整型数据中包含的pre-bwt：bwt // 打印32位整型数据中包含的pre-bwt：bwt
static void print_base_uint32(uint32_t p) static void print_base_uint32(uint32_t p)
{ {
@ -165,7 +167,7 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
} }
k -= (k >= bwt->primary); // because $ is not in bwt k -= (k >= bwt->primary); // because $ is not in bwt
p = bwt_occ_intv(bwt, k); p = bwt_occ_intv(bwt, k);
// cout << "k: " << k << "; occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl; //cout << "k: " << k << "; occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl;
cnt[0] = p[0]; cnt[0] = p[0];
cnt[1] = p[1]; cnt[1] = p[1];
cnt[2] = p[2]; cnt[2] = p[2];
@ -176,16 +178,20 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
// cout << "k - kbase: " << k - bwt_k_base_line << endl; // cout << "k - kbase: " << k - bwt_k_base_line << endl;
for (; p < end; ++p) for (; p < end; ++p)
{ {
x += __occ_aux4(bwt, *p); x += __occ_aux4(bwt, *p);
//print_base_uint32(*p); //if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550)
// print_base_uint32(*p);
} }
tmp = *p & ~((1U << ((~k & 15) << 1)) - 1); tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
// print_base_uint32(tmp); //if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550)
// print_base_uint32(tmp);
x += __occ_aux4(bwt, tmp) - (~k & 15); // 这里多算了A，要减去 x += __occ_aux4(bwt, tmp) - (~k & 15); // 这里多算了A，要减去
cnt[0] += x & 0xff; cnt[0] += x & 0xff;
cnt[1] += x >> 8 & 0xff; cnt[1] += x >> 8 & 0xff;
cnt[2] += x >> 16 & 0xff; cnt[2] += x >> 16 & 0xff;
cnt[3] += x >> 24; cnt[3] += x >> 24;
//if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550)
// cout << "final-occ: " << cnt[0] << ' ' << cnt[1] << ' ' << cnt[2] << ' ' << cnt[3] << endl;
} }
// an analogy to bwt_occ4() but more efficient, requiring k <= l // an analogy to bwt_occ4() but more efficient, requiring k <= l
@ -210,8 +216,11 @@ void bwt_2occ4(const bwt_t *bwt, bwtint_t k, bwtint_t l, bwtint_t cntk[4], bwtin
k -= (k >= bwt->primary); // because $ is not in bwt k -= (k >= bwt->primary); // because $ is not in bwt
l -= (l >= bwt->primary); l -= (l >= bwt->primary);
p = bwt_occ_intv(bwt, k); p = bwt_occ_intv(bwt, k);
memcpy(cntk, p, 4 * sizeof(bwtint_t)); cntk[0] = p[0];
p += sizeof(bwtint_t); // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t)) cntk[1] = p[1];
cntk[2] = p[2];
cntk[3] = p[3];
p += 4; // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t))
// prepare cntk[] // prepare cntk[]
endk = p + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4)); endk = p + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4));
endl = p + ((l >> 4) - ((l & ~OCC_INTV_MASK) >> 4)); endl = p + ((l >> 4) - ((l & ~OCC_INTV_MASK) >> 4));
@ -356,12 +365,13 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q)
c2 = cbval(q[i + 1]); c2 = cbval(q[i + 1]);
// double tm_t = realtime(); // double tm_t = realtime();
//cout << BASE[c2] << BASE[c1] << endl;
bwt_extend2(bwt, &ik, ok, 0, c1, c2); bwt_extend2(bwt, &ik, ok, 0, c1, c2);
// t7 += realtime() - tm_t; // t7 += realtime() - tm_t;
ik = ok[c2]; ik = ok[c2];
ik.info = i + 1; ik.info = i + 1;
// cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl; //cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
} }
else else
{ {

84
fmt_index.cpp
View File

@ -19,7 +19,7 @@ double t1 = 0, t2 = 0, t3 = 0, t4 = 0, t5 = 0, t6 = 0, t7 = 0, t8 = 0, t9 = 0, t
long f1 = 0, f2 = 0, f3 = 0, f4 = 0, f5 = 0; long f1 = 0, f2 = 0, f3 = 0, f4 = 0, f5 = 0;
const char BASE[4] = {'A', 'C', 'G', 'T'}; const static char BASE[4] = {'A', 'C', 'G', 'T'};
// 求反向互补序列 // 求反向互补序列
string calc_reverse_seq(string &seq) string calc_reverse_seq(string &seq)
@ -207,6 +207,8 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
bwtint_t last_seq_len = bwt->seq_len % FMT_OCC_INTERVAL; bwtint_t last_seq_len = bwt->seq_len % FMT_OCC_INTERVAL;
mn_occ += (last_seq_len + FMT_MID_INTERVAL - 1) / FMT_MID_INTERVAL - 1; mn_occ += (last_seq_len + FMT_MID_INTERVAL - 1) / FMT_MID_INTERVAL - 1;
fmt->bwt_size += mn_occ * 4; fmt->bwt_size += mn_occ * 4;
//cout << "mn_occ: " << mn_occ << endl;
//exit(0);
#endif #endif
buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存 buf = (uint32_t *)calloc(fmt->bwt_size, 4); // 开辟计算fmt用到的缓存
@ -230,7 +232,9 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
k += 4; k += 4;
memcpy(buf + k, c2, sizeof(uint32_t) * 16); // pre-bwt:bwt碱基对的occ memcpy(buf + k, c2, sizeof(uint32_t) * 16); // pre-bwt:bwt碱基对的occ
k += 16; k += 16;
#ifdef FMT_MID_INTERVAL
mc[0] = mc[1] = mc[2] = mc[3] = 0; mc[0] = mc[1] = mc[2] = mc[3] = 0;
#endif
} }
// 每个32位整数保存8个倒数第二列碱基（pre-bwt）和8个倒数第一列(bwt)碱基 // 每个32位整数保存8个倒数第二列碱基（pre-bwt）和8个倒数第一列(bwt)碱基
if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16 if (i % 16 == 0) // 每个32位整数可以包含16个碱基，每次需要处理16个碱基，也就是间隔最小可以设置为16
@ -297,16 +301,36 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
pre_and_bwt_seq |= (((pre_bwt_16_seq & (3 << lshift_bit)) >> (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) >> ((m * 2) + 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; buf[k++] = pre_and_bwt_seq;
if (j > 8) if (j > 8)
{ {
#if FMT_MID_INTERVAL == 8
for (m = 0; m < 4; ++m)
{
uint32_t s1 = pre_and_bwt_seq;
mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24];
buf[k++] = mc[m];
}
#endif
for (m = 8; m > 0; --m) for (m = 8; m > 0; --m)
{ {
int lshift_bit = 2 * m - 2; int lshift_bit = 2 * m - 2;
pre_and_bwt_seq_2 |= (((pre_bwt_16_seq & (3 << lshift_bit)) << (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) << (m * 2 - 2))); pre_and_bwt_seq_2 |= (((pre_bwt_16_seq & (3 << lshift_bit)) << (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) << (m * 2 - 2)));
} }
buf[k++] = pre_and_bwt_seq_2; buf[k++] = pre_and_bwt_seq_2;
#if FMT_MID_INTERVAL == 8
if ((i + 16) % FMT_OCC_INTERVAL != 0 && j == 16)
for (m = 0; m < 4; ++m)
{
uint32_t s1 = pre_and_bwt_seq_2;
mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24];
buf[k++] = mc[m];
}
#endif
} }
#ifdef FMT_MID_INTERVAL #if FMT_MID_INTERVAL == 16
// 加入中间的check point // 加入中间的check point
for (m = 0; m < 4; ++m) for (m = 0; m < 4; ++m)
{ {
@ -329,7 +353,7 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
memcpy(buf + k, c2, sizeof(uint32_t) * 16); memcpy(buf + k, c2, sizeof(uint32_t) * 16);
k += 16; k += 16;
// cout << "n occ: " << n_occ << endl; // cout << "n occ: " << n_occ << endl;
// cout << "size: " << k << '\t' << fmt->bwt_size << endl; cout << "size: " << k << '\t' << fmt->bwt_size << endl;
xassert(k == fmt->bwt_size, "inconsistent bwt_size"); xassert(k == fmt->bwt_size, "inconsistent bwt_size");
// update fmt // update fmt
fmt->bwt = buf; fmt->bwt = buf;
@ -356,13 +380,17 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]
ti = b1 << 2 | b2; ti = b1 << 2 | b2;
k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1） k -= (k >= fmt->primary); // k由bwt矩阵对应的行转换成bwt字符串对应的行（去掉了$，所以大于$的行，都减掉1）
p = fmt_occ_intv(fmt, k); p = fmt_occ_intv(fmt, k);
// cout << "k-base: " << k << "; occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl; //cout << "k-base: " << k << "; occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl;
for (i = b1 + 1; i < 4; ++i) cnt[0] += p[i]; // 大于b1的碱基的occ之和 for (i = b1 + 1; i < 4; ++i) cnt[0] += p[i]; // 大于b1的碱基的occ之和
cnt[1] = p[b1]; // b1的occ cnt[1] = p[b1]; // b1的occ
q = p + 4 + b1 * 4; q = p + 4 + b1 * 4;
for (i = b2 + 1; i < 4 ; ++i) cnt[2] += q[i]; // 大于b2的occ之和 for (i = b2 + 1; i < 4 ; ++i) cnt[2] += q[i]; // 大于b2的occ之和
cnt[3] = q[b2]; // b2的occ cnt[3] = q[b2]; // b2的occ
p += 20; p += 20;
//if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550) {
// cout << "sec base-occ: " << q[0] << ' ' << q[1] << ' ' << q[2] << ' ' << q[3] << endl;
//}
#ifdef FMT_MID_INTERVAL #ifdef FMT_MID_INTERVAL
// 使用mid interval信息 // 使用mid interval信息
int mk = k % FMT_OCC_INTERVAL; int mk = k % FMT_OCC_INTERVAL;
@ -387,7 +415,7 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]
//} //}
if (n_mintv > 0) // 至少超过了第一个mid interval if (n_mintv > 0) // 至少超过了第一个mid interval
{ {
p += n_mintv * 6 - 4; // 对应的mid interval check point的首地址，即A C G T的局部累积量 p += n_mintv * (4 + (FMT_MID_INTERVAL >> 3)) - 4; // 对应的mid interval check point的首地址，即A C G T的局部累积量
q = p + b1; q = p + b1;
for (i = b1 + 1; i < 4; ++i) for (i = b1 + 1; i < 4; ++i)
x += p[i]; // 大于b1的碱基的occ之和 x += p[i]; // 大于b1的碱基的occ之和
@ -406,18 +434,24 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]
x = 0; x = 0;
p += 4; p += 4;
} }
// cout << "mid-occ: " << cnt[0] << ' ' << cnt[1] << ' ' << cnt[2] << ' ' << cnt[3] << endl;
// cout << "k: " << k << ' ' << cnt[1] << endl; // cout << "k: " << k << ' ' << cnt[1] << endl;
#if FMT_MID_INTERVAL == 16
if ((mk & FMT_MID_INTV_MASK) >> 3) if ((mk & FMT_MID_INTV_MASK) >> 3)
{ {
x += __fmt_occ_e2_aux2(fmt, ti, *p); x += __fmt_occ_e2_aux2(fmt, ti, *p);
// print_base_uint32(*p); // print_base_uint32(*p);
++p; ++p;
} }
#endif
tmp = *p & ~((1U << ((~k & 7) << 2)) - 1); tmp = *p & ~((1U << ((~k & 7) << 2)) - 1);
// print_base_uint32(tmp);
x += __fmt_occ_e2_aux2(fmt, ti, tmp); x += __fmt_occ_e2_aux2(fmt, ti, tmp);
//if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550)
//{
// print_base_uint32(tmp);
// cout << "多的: " << (~k & 7) << endl;
// cout << (x >> 24 & 0xff) << ' ' << (x >> 16 & 0xff) << ' ' << (x >> 8 & 0xff) << ' ' << (x & 0xff) << endl;
//}
#else // 该地址是bwt和pre_bwt字符串数据的首地址 #else // 该地址是bwt和pre_bwt字符串数据的首地址
uint32_t *end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); // this is the end point of the following loop uint32_t *end = p + ((k >> 3) - ((k & ~FMT_OCC_INTV_MASK) >> 3)); // this is the end point of the following loop
// p = end - (end - p) / 8; // p = end - (end - p) / 8;
@ -437,19 +471,30 @@ void fmt_e2_occ(const FMTIndex *fmt, bwtint_t k, int b1, int b2, bwtint_t cnt[4]
x -= (~k & 7) << 8; x -= (~k & 7) << 8;
if (b2 == 0) if (b2 == 0)
x -= (~k & 7) << 24; x -= (~k & 7) << 24;
// cout << "here-1" << endl;
} }
// 如果跨过了second_primary,那么可能需要减掉一次累积值 // 如果跨过了second_primary,那么可能需要减掉一次累积值
if (b1 == fmt->first_base && bwt_k_base_line < fmt->sec_primary && bwt_k_line >= fmt->sec_primary) if (b1 == fmt->first_base && bwt_k_base_line < fmt->sec_primary && bwt_k_line >= fmt->sec_primary)
{ {
if (b2 < fmt->last_base) if (b2 < fmt->last_base)
x -= 1 << 16; cnt[2] -= 1 ;
else if (b2 == fmt->last_base) else if (b2 == fmt->last_base)
x -= 1 << 24; cnt[3] -= 1;
// cout << "here-2" << endl;
} }
// if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550)
// {
// cout << BASE[b1] << BASE[b2] << endl;
// cout << "final-x: " << (x >> 24 & 0xff) << ' ' << (x >> 16 & 0xff) << ' ' << (x >> 8 & 0xff) << ' ' << (x & 0xff) << endl;
// }
cnt[0] += x & 0xff; cnt[0] += x & 0xff;
cnt[1] += x >> 8 & 0xff; cnt[1] += x >> 8 & 0xff;
cnt[2] += x >> 16 & 0xff; cnt[2] += x >> 16 & 0xff;
cnt[3] += x >> 24 & 0xff; cnt[3] += x >> 24 & 0xff;
// if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550)
// cout << "final-occ: " << cnt[0] << ' ' << cnt[1] << ' ' << cnt[2] << ' ' << cnt[3] << endl;
} }
// 扩展两个碱基 // 扩展两个碱基
@ -501,7 +546,7 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
fmt_set_intv(fmt, bval(q[x]), ik); fmt_set_intv(fmt, bval(q[x]), ik);
ik.info = x + 1; ik.info = x + 1;
// cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl; //cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
for (i = x + 1; i + 1 < qlen; i += 2) for (i = x + 1; i + 1 < qlen; i += 2)
{ {
@ -509,12 +554,13 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
{ {
c1 = cbval(q[i]); c1 = cbval(q[i]);
c2 = cbval(q[i + 1]); c2 = cbval(q[i + 1]);
//cout << BASE[c2] << BASE[c1] << endl;
//double tm_t = realtime(); //double tm_t = realtime();
fmt_extend2(fmt, &ik, &ok, 0, c1, c2); fmt_extend2(fmt, &ik, &ok, 0, c1, c2);
ik = ok; ik = ok;
// t8 += realtime() - tm_t; // t8 += realtime() - tm_t;
ik.info = i + 1; ik.info = i + 1;
// cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl; //cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
} }
else else
{ {
@ -576,7 +622,7 @@ int main_fmtidx(int argc, char **argv)
// string bwt_str = string(argv[1]) + ".bwt.str"; // string bwt_str = string(argv[1]) + ".bwt.str";
string bwt_idx = string(argv[1]) + ".128.bwt"; string bwt_idx = string(argv[1]) + ".128.bwt";
//string bwt_idx = string(argv[1]) + ".64.bwt"; // string bwt_idx = string(argv[1]) + ".64.bwt";
// string bwt_idx = string(argv[1]) + ".16.bwt"; // string bwt_idx = string(argv[1]) + ".16.bwt";
// string bwt_idx = string(argv[1]) + ".bwt"; // string bwt_idx = string(argv[1]) + ".bwt";
//bwt_t *bwt = restore_bwt(bwt_str.c_str()); //bwt_t *bwt = restore_bwt(bwt_str.c_str());
@ -584,8 +630,9 @@ int main_fmtidx(int argc, char **argv)
//dump_bwt(bwt_idx.c_str(), bwt); //dump_bwt(bwt_idx.c_str(), bwt);
//string fmt_idx = string(argv[1]) + ".fmt"; //string fmt_idx = string(argv[1]) + ".fmt";
//string fmt_idx = string(argv[1]) + ".256.8.fmt";
string fmt_idx = string(argv[1]) + ".256.fmt"; string fmt_idx = string(argv[1]) + ".256.fmt";
// string fmt_idx = string(argv[1]) + ".128.fmt"; //string fmt_idx = string(argv[1]) + ".128.fmt";
// string fmt_idx = string(argv[1]) + ".64.fmt"; // string fmt_idx = string(argv[1]) + ".64.fmt";
// string fmt_idx = string(argv[1]) + ".32.fmt"; // string fmt_idx = string(argv[1]) + ".32.fmt";
@ -595,18 +642,19 @@ int main_fmtidx(int argc, char **argv)
cout << "[time read bwt:] " << t11 << "s" << endl; cout << "[time read bwt:] " << t11 << "s" << endl;
t10 = realtime(); t10 = realtime();
FMTIndex *fmt = restore_fmt(fmt_idx.c_str()); FMTIndex *fmt = restore_fmt(fmt_idx.c_str());
// FMTIndex *fmt = create_fmt_from_bwt(bwt); //FMTIndex *fmt = create_fmt_from_bwt(bwt);
// dump_fmt(fmt_idx.c_str(), fmt); //dump_fmt(fmt_idx.c_str(), fmt);
t10 = realtime() - t10; t10 = realtime() - t10;
cout << "[time gen fmt:] " << t10 << "s" << endl; cout << "[time gen fmt:] " << t10 << "s" << endl;
vector<string> seed_arr(10000000); vector<string> seed_arr(10000000);
seed_arr[0] = "GCGATACTAAGA"; seed_arr[0] = "GCGATACTAAGA";
seed_arr[0] = "GAGAGCTGTTCCCGTGTTTTCCATGGTTT";
srand(time(NULL)); srand(time(NULL));
t1 = realtime(); t1 = realtime();
for (int i = 0; i < (int)seed_arr.size(); ++i) for (int i = 0; i < (int)seed_arr.size(); ++i)
seed_arr[i] = generate_rand_seq(13); seed_arr[i] = generate_rand_seq(7);
t1 = realtime() - t1; t1 = realtime() - t1;
cout << "[time gen seed:] " << t1 << "s" << endl; cout << "[time gen seed:] " << t1 << "s" << endl;
@ -621,7 +669,7 @@ int main_fmtidx(int argc, char **argv)
t3 = realtime() - t3; t3 = realtime() - t3;
cout << "[time fmt search:] " << t3 << "s" << endl; cout << "[time fmt search:] " << t3 << "s" << endl;
return 0; //return 0;
t4 = realtime(); t4 = realtime();
for (int i = 0; i < (int)seed_arr.size(); ++i) for (int i = 0; i < (int)seed_arr.size(); ++i)

4
fmt_index.h
View File

@ -12,9 +12,11 @@
// 获取碱基c（待查找序列的首个碱基）和对应的互补碱基对应的行，以及间隔 // 获取碱基c（待查找序列的首个碱基）和对应的互补碱基对应的行，以及间隔
#define fmt_set_intv(fmt, c, ik) ((ik).x[0] = (fmt)->L2[(int)(c)] + 1, (ik).x[2] = (fmt)->L2[(int)(c) + 1] - (fmt)->L2[(int)(c)], (ik).x[1] = (fmt)->L2[3 - (c)] + 1, (ik).info = 0) #define fmt_set_intv(fmt, c, ik) ((ik).x[0] = (fmt)->L2[(int)(c)] + 1, (ik).x[2] = (fmt)->L2[(int)(c) + 1] - (fmt)->L2[(int)(c)], (ik).x[1] = (fmt)->L2[3 - (c)] + 1, (ik).info = 0)
// k行（bwt str行（不包含$））对应的check point occ数据起始地址（小于k且是OCC_INTERVAL的整数倍） // k行（bwt str行（不包含$））对应的check point occ数据起始地址（小于k且是OCC_INTERVAL的整数倍）
#ifdef FMT_MID_INTERVAL #if FMT_MID_INTERVAL == 16
// 包含mid interval // 包含mid interval
#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 80)) #define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 80))
#elif FMT_MID_INTERVAL == 8
#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 144))
#else #else
#define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 20)) #define fmt_occ_intv(b, k) ((b)->bwt + (k) / FMT_OCC_INTERVAL * (FMT_OCC_INTERVAL / 8 + 20))
#endif #endif