清理了一下代码，目前比较纯净，无bug，方便各种设置和测试
This commit is contained in:
zzh
parent
45318fc95b
commit
5a055625c3
|
|
@ -1,6 +1,6 @@
|
|||
# specific for bwa_perf
|
||||
*.txt
|
||||
fmtidx
|
||||
bwa_perf
|
||||
*.fmt
|
||||
|
||||
# ---> C
|
||||
|
|
|
|||
4
Makefile
4
Makefile
|
|
@ -7,7 +7,7 @@ SHOW_PERF= -DSHOW_PERF
|
|||
AR= ar
|
||||
DFLAGS= -DHAVE_PTHREAD $(WRAP_MALLOC) $(SHOW_PERF)
|
||||
AOBJS= util.o sa.o fmt_index.o bwt.o
|
||||
PROG= fmtidx
|
||||
PROG= bwa_perf
|
||||
INCLUDES=
|
||||
LIBS= -lm -lz -lpthread
|
||||
SUBDIRS= .
|
||||
|
|
@ -23,7 +23,7 @@ endif
|
|||
|
||||
all:$(PROG)
|
||||
|
||||
fmtidx:$(AOBJS) main.o
|
||||
$(PROG):$(AOBJS) main.o
|
||||
$(CC) $(CFLAGS) $(LDFLAGS) $(AOBJS) main.o -o $@ -L. $(LIBS)
|
||||
|
||||
clean:
|
||||
|
|
|
|||
66
bwt.cpp
66
bwt.cpp
|
|
@ -15,7 +15,6 @@ const static char BASE[4] = {'A', 'C', 'G', 'T'};
|
|||
// 打印32位整型数据中包含的pre-bwt：bwt
|
||||
static void print_base_uint32(uint32_t p)
|
||||
{
|
||||
const char BASE[4] = {'A', 'C', 'G', 'T'};
|
||||
for (int i = 30; i > 0; i -= 2)
|
||||
{
|
||||
int b1 = p >> i & 3;
|
||||
|
|
@ -23,6 +22,7 @@ static void print_base_uint32(uint32_t p)
|
|||
}
|
||||
}
|
||||
|
||||
// 保存bwt数据
|
||||
void dump_bwt(const char *fn, const bwt_t *bwt)
|
||||
{
|
||||
FILE *fp;
|
||||
|
|
@ -90,11 +90,6 @@ void create_interval_occ_bwt(bwt_t *bwt)
|
|||
// 计算occ，生成naive bwt
|
||||
for (i = k = 0; i < bwt->seq_len; ++i)
|
||||
{
|
||||
if (i % 1000000 == 0)
|
||||
{
|
||||
cout << "line: " << i << '\t' << ((uint64_t)c[0] + c[1] + c[2] + c[3]) << '\t';
|
||||
cout << c[0] << ' ' << c[1] << ' ' << c[2] << ' ' << c[3] << endl;
|
||||
}
|
||||
if (i % OCC_INTERVAL == 0)
|
||||
{
|
||||
memcpy(buf + k, c, sizeof(uint32_t) * 4);
|
||||
|
|
@ -105,7 +100,6 @@ void create_interval_occ_bwt(bwt_t *bwt)
|
|||
++c[bwt_B00(bwt, i)];
|
||||
}
|
||||
// the last element
|
||||
// cout << c[0] << '\t' << c[1] << '\t' << c[2] << '\t' << c[3] << endl;
|
||||
memcpy(buf + k, c, sizeof(uint32_t) * 4);
|
||||
xassert(k + 4 == bwt->bwt_size, "inconsistent bwt_size");
|
||||
// update bwt
|
||||
|
|
@ -137,7 +131,6 @@ bwtint_t bwt_occ(const bwt_t *bwt, bwtint_t k, uint8_t c)
|
|||
|
||||
// retrieve Occ at k/OCC_INTERVAL
|
||||
n = ((bwtint_t *)(p = bwt_occ_intv(bwt, k)))[c];
|
||||
// cout << "bwt_occ - 1: " << (int)c << '\t' << k << '\t' << n << endl;
|
||||
p += sizeof(bwtint_t); // jump to the start of the first BWT cell
|
||||
|
||||
// calculate Occ up to the last k/32
|
||||
|
|
@ -145,12 +138,10 @@ bwtint_t bwt_occ(const bwt_t *bwt, bwtint_t k, uint8_t c)
|
|||
for (; p < end; p += 2)
|
||||
n += __occ_aux((uint64_t)p[0] << 32 | p[1], c);
|
||||
|
||||
// cout << "bwt_occ - 2: " << (int)c << '\t' << k << '\t' << n << endl;
|
||||
// calculate Occ
|
||||
n += __occ_aux(((uint64_t)p[0] << 32 | p[1]) & ~((1ull << ((~k & 31) << 1)) - 1), c);
|
||||
if (c == 0)
|
||||
n -= ~k & 31; // corrected for the masked bits
|
||||
// cout << "bwt_occ - 3: " << (int)c << '\t' << k << '\t' << n << endl;
|
||||
return n;
|
||||
}
|
||||
|
||||
|
|
@ -167,37 +158,27 @@ void bwt_occ4(const bwt_t *bwt, bwtint_t k, bwtint_t cnt[4])
|
|||
}
|
||||
k -= (k >= bwt->primary); // because $ is not in bwt
|
||||
p = bwt_occ_intv(bwt, k);
|
||||
//cout << "k: " << k << "; occ: " << p[0] << ' ' << p[1] << ' ' << p[2] << ' ' << p[3] << endl;
|
||||
cnt[0] = p[0];
|
||||
cnt[1] = p[1];
|
||||
cnt[2] = p[2];
|
||||
cnt[3] = p[3];
|
||||
p += 4; // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t))
|
||||
p += 4; // check point之后的bwt字符串首地址
|
||||
end = p + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4)); // this is the end point of the following loop
|
||||
//end = p + 1;
|
||||
// cout << "k - kbase: " << k - bwt_k_base_line << endl;
|
||||
for (; p < end; ++p)
|
||||
{
|
||||
x += __occ_aux4(bwt, *p);
|
||||
//if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550)
|
||||
// print_base_uint32(*p);
|
||||
}
|
||||
tmp = *p & ~((1U << ((~k & 15) << 1)) - 1);
|
||||
//if (k == 3965453116 || k == 3965453672 || k == 2668688087 || k == 2668688550)
|
||||
// print_base_uint32(tmp);
|
||||
x += __occ_aux4(bwt, tmp) - (~k & 15); // 这里多算了A，要减去
|
||||
cnt[0] += x & 0xff;
|
||||
cnt[1] += x >> 8 & 0xff;
|
||||
cnt[2] += x >> 16 & 0xff;
|
||||
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，有提升，但不大
|
||||
void bwt_2occ4(const bwt_t *bwt, bwtint_t k, bwtint_t l, bwtint_t cntk[4], bwtint_t cntl[4])
|
||||
{
|
||||
// double tm_t = realtime();
|
||||
// bwt_occ4(bwt, k, cntk);
|
||||
// bwt_occ4(bwt, l, cntl);
|
||||
// return;
|
||||
|
|
@ -220,7 +201,7 @@ void bwt_2occ4(const bwt_t *bwt, bwtint_t k, bwtint_t l, bwtint_t cntk[4], bwtin
|
|||
cntk[1] = p[1];
|
||||
cntk[2] = p[2];
|
||||
cntk[3] = p[3];
|
||||
p += 4; // sizeof(bwtint_t) = 4*(sizeof(bwtint_t)/sizeof(uint32_t))
|
||||
p += 4;
|
||||
// prepare cntk[]
|
||||
endk = p + ((k >> 4) - ((k & ~OCC_INTV_MASK) >> 4));
|
||||
endl = p + ((l >> 4) - ((l & ~OCC_INTV_MASK) >> 4));
|
||||
|
|
@ -244,8 +225,6 @@ void bwt_2occ4(const bwt_t *bwt, bwtint_t k, bwtint_t l, bwtint_t cntk[4], bwtin
|
|||
cntl[2] += y >> 16 & 0xff;
|
||||
cntl[3] += y >> 24;
|
||||
}
|
||||
// t5 += realtime() - tm_t;
|
||||
// f1 += 1;
|
||||
}
|
||||
|
||||
void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_back)
|
||||
|
|
@ -253,9 +232,7 @@ void bwt_extend(const bwt_t *bwt, const bwtintv_t *ik, bwtintv_t ok[4], int is_b
|
|||
bwtint_t tk[4], tl[4];
|
||||
int i;
|
||||
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl); // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
||||
// cout << "bwt-1-d: " << ik->x[!is_back] - 1 << '\t' << tk[0] << '\t' << tk[1] << '\t' << tk[2] << '\t' << tk[3] << endl;
|
||||
// cout << "bwt-1-d: " << ik->x[!is_back] - 1 + ik->x[2] << '\t' << tl[0] << '\t' << tl[1] << '\t' << tl[2] << '\t' << tl[3] << endl;
|
||||
// 这里是反向扩展
|
||||
// 这里是反向扩展
|
||||
for (i = 0; i != 4; ++i)
|
||||
{
|
||||
ok[i].x[!is_back] = bwt->L2[i] + 1 + tk[i]; // 起始行位置，互补链
|
||||
|
|
@ -276,7 +253,6 @@ void bwt_search(bwt_t *bwt, const string &q)
|
|||
|
||||
bwt_set_intv(bwt, 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 < (int)q.size(); ++i)
|
||||
{
|
||||
|
|
@ -286,7 +262,6 @@ void bwt_search(bwt_t *bwt, const string &q)
|
|||
bwt_extend(bwt, &ik, ok, 0);
|
||||
ik = ok[c];
|
||||
ik.info = i + 1;
|
||||
// cout << "bwt-1: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
|
@ -295,8 +270,6 @@ void bwt_search(bwt_t *bwt, const string &q)
|
|||
void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back, int c1, int c2)
|
||||
{
|
||||
bwtint_t tk[4], tl[4];
|
||||
// oc1, oc2;
|
||||
// int int1 = 0, int2 = 0;
|
||||
int i;
|
||||
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl); // tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
||||
|
||||
|
|
@ -312,17 +285,8 @@ void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back,
|
|||
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];
|
||||
//cout << (ik->x[!is_back] <= bwt->primary && ik->x[!is_back] + ik->x[2] - 1 >= bwt->primary) << endl;
|
||||
*ik = ok[c1];
|
||||
|
||||
//oc1 = oc2 = 0;
|
||||
//for (i = 3; i > c1; --i)
|
||||
//{
|
||||
// int1 += tl[i] - tk[i];
|
||||
// oc1 += tk[i];
|
||||
// oc2 += tl[i];
|
||||
//}
|
||||
//cout << "bwt-2 interval-1: " << int1 << ' ' << oc1 << ' ' << oc2 << endl;
|
||||
*ik = ok[c1];
|
||||
|
||||
bwt_2occ4(bwt, ik->x[!is_back] - 1, ik->x[!is_back] - 1 + ik->x[2], tk, tl);
|
||||
for (i = 0; i != 4; ++i)
|
||||
|
|
@ -335,15 +299,6 @@ void bwt_extend2(const bwt_t *bwt, bwtintv_t *ik, bwtintv_t ok[4], int is_back,
|
|||
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];
|
||||
|
||||
// oc1 = oc2 = 0;
|
||||
// for (i = 3; i > c2; --i)
|
||||
// {
|
||||
// int2 += tl[i] - tk[i];
|
||||
// oc1 += tk[i];
|
||||
// oc2 += tl[i];
|
||||
// }
|
||||
// cout << "bwt-2 interval-2: " << int2 << ' ' << oc1 << ' ' << oc2 << endl;
|
||||
}
|
||||
|
||||
// 每次扩展两步
|
||||
|
|
@ -354,7 +309,6 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q)
|
|||
|
||||
bwt_set_intv(bwt, bval(q[x]), ik);
|
||||
ik.info = x + 1;
|
||||
//cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
|
||||
for (i = x + 1; i + 1 < (int)q.size(); i += 2)
|
||||
{
|
||||
|
|
@ -364,14 +318,10 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q)
|
|||
c1 = cbval(q[i]);
|
||||
c2 = cbval(q[i + 1]);
|
||||
|
||||
// double tm_t = realtime();
|
||||
//cout << BASE[c2] << BASE[c1] << endl;
|
||||
bwt_extend2(bwt, &ik, ok, 0, c1, c2);
|
||||
// t7 += realtime() - tm_t;
|
||||
|
||||
ik = ok[c2];
|
||||
ik.info = i + 1;
|
||||
//cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
@ -381,13 +331,9 @@ bwtintv_t bwt_search2(bwt_t *bwt, const string &q)
|
|||
if (i < (int)q.size() && bval(q[i]) < 4)
|
||||
{ // 最后一次扩展
|
||||
c1 = cbval(q[i]);
|
||||
// double tm_t = realtime();
|
||||
bwt_extend(bwt, &ik, ok, 0);
|
||||
// t10 += realtime() - tm_t;
|
||||
ik = ok[c1];
|
||||
ik.info = i + 1;
|
||||
//cout << "bwt-2: " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
}
|
||||
// cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
return ik;
|
||||
}
|
||||
6
bwt.h
6
bwt.h
|
|
@ -7,7 +7,7 @@
|
|||
using std::string;
|
||||
|
||||
// occ间隔对应的右移位数，5表示间隔32行保存一次
|
||||
#define OCC_INTV_SHIFT 7
|
||||
#define OCC_INTV_SHIFT 6
|
||||
#define OCC_INTERVAL (1LL << OCC_INTV_SHIFT)
|
||||
#define OCC_INTV_MASK (OCC_INTERVAL - 1)
|
||||
|
||||
|
|
@ -23,10 +23,6 @@ using std::string;
|
|||
// 获取碱基c（待查找序列的首个碱基）和对应的互补碱基对应的行，以及间隔
|
||||
#define bwt_set_intv(bwt, c, ik) ((ik).x[0] = (bwt)->L2[(int)(c)] + 1, (ik).x[2] = (bwt)->L2[(int)(c) + 1] - (bwt)->L2[(int)(c)], (ik).x[1] = (bwt)->L2[3 - (c)] + 1, (ik).info = 0)
|
||||
|
||||
// The following two lines are ONLY correct when OCC_INTERVAL==0x80
|
||||
// #define bwt_bwt(b, k) ((b)->bwt[((k) >> 7 << 4) + sizeof(bwtint_t) + (((k) & 0x7f) >> 4)])
|
||||
// #define bwt_occ_intv(b, k) ((b)->bwt + ((k) >> 7 << 4))
|
||||
///* For general OCC_INTERVAL, the following is correct:
|
||||
// k行对应的bwt str碱基，对应的check point bwt str数据起始地址
|
||||
#define bwt_bwt(b, k) ((b)->bwt[(k) / OCC_INTERVAL * (OCC_INTERVAL / 16 + 4) + 4 + (k) % OCC_INTERVAL / 16])
|
||||
// k行（bwt str行（不包含$））对应的check point occ数据起始地址（小于k且是OCC_INTERVAL的整数倍）
|
||||
|
|
|
|||
|
|
@ -315,13 +315,13 @@ FMTIndex *create_fmt_from_bwt(bwt_t *bwt)
|
|||
pre_and_bwt_seq_2 |= (((pre_bwt_16_seq & (3 << lshift_bit)) << (m * 2)) | ((bwt_16_seq & (3 << lshift_bit)) << (m * 2 - 2)));
|
||||
}
|
||||
|
||||
#ifdef FMT_MID_INTERVAL
|
||||
#ifdef FMT_MID_INTERVAL // 计算前边8+8个碱基的mid interval occ
|
||||
s1 = pre_and_bwt_seq;
|
||||
for (m = 0; m < 4; ++m)
|
||||
mc[m] += cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24];
|
||||
#endif
|
||||
|
||||
#if FMT_MID_INTERVAL == 8
|
||||
#if FMT_MID_INTERVAL == 8 // 如果mid interval是8的话，这里要保存一次
|
||||
for (m = 0; m < 4; ++m)
|
||||
buf[k++] = mc[m];
|
||||
#endif
|
||||
|
|
@ -458,13 +458,7 @@ void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back,
|
|||
// tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
||||
fmt_e2_occ(fmt, ik->x[!is_back] - 1, b1, b2, tk);
|
||||
fmt_e2_occ(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl);
|
||||
|
||||
//fmt_e2_occ_2way(fmt, ik->x[!is_back] - 1, b1, b2, tk);
|
||||
//fmt_e2_occ_2way(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl);
|
||||
// 这里是反向扩展
|
||||
//cout << BASE[b1] << BASE[b2] << endl;
|
||||
//cout << "fmt: interval-1: " << tl[0] - tk[0] << ' ' << tl[0] << ' ' << tk[0] << endl;
|
||||
//cout << "fmt: interval-2: " << tl[2] - tk[2] << ' ' << tl[2] << ' ' << tk[2] << endl;
|
||||
ok->x[!is_back] = fmt->L2[b2] + 1 + tk[3];
|
||||
ok->x[2] = tl[3] - tk[3];
|
||||
// 第一次正向扩展
|
||||
|
|
@ -478,12 +472,10 @@ void fmt_extend2(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back,
|
|||
void fmt_extend1(const FMTIndex *fmt, bwtintv_t *ik, bwtintv_t *ok, int is_back, int b1)
|
||||
{
|
||||
bwtint_t tk[4], tl[4];
|
||||
int b2 = 3;
|
||||
int b2 = 3; // 如果只扩展一次，那么第二个碱基设置成T，可以减小一些计算量，如计算大于b2的累积数量
|
||||
// tk表示在k行之前所有各个碱基累积出现次数，tl表示在l行之前的累积
|
||||
fmt_e2_occ(fmt, ik->x[!is_back] - 1, b1, b2, tk);
|
||||
fmt_e2_occ(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl);
|
||||
//fmt_e2_occ_2way(fmt, ik->x[!is_back] - 1, b1, b2, tk);
|
||||
//fmt_e2_occ_2way(fmt, ik->x[!is_back] - 1 + ik->x[2], b1, b2, tl);
|
||||
// 这里是反向扩展
|
||||
ok->x[!is_back] = fmt->L2[b1] + 1 + tk[1];
|
||||
ok->x[2] = tl[1] - tk[1];
|
||||
|
|
@ -501,21 +493,17 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
|
|||
|
||||
fmt_set_intv(fmt, bval(q[x]), ik);
|
||||
ik.info = x + 1;
|
||||
//cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
|
||||
// 每次扩展两个碱基
|
||||
for (i = x + 1; i + 1 < qlen; i += 2)
|
||||
{
|
||||
if (bval(q[i]) < 4 && bval(q[i + 1]) < 4)
|
||||
{
|
||||
c1 = cbval(q[i]);
|
||||
c2 = cbval(q[i + 1]);
|
||||
//cout << BASE[c2] << BASE[c1] << endl;
|
||||
//double tm_t = realtime();
|
||||
fmt_extend2(fmt, &ik, &ok, 0, c1, c2);
|
||||
ik = ok;
|
||||
// t8 += realtime() - tm_t;
|
||||
ik.info = i + 1;
|
||||
//cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
}
|
||||
else
|
||||
{
|
||||
|
|
@ -525,14 +513,10 @@ bwtintv_t fmt_search(FMTIndex *fmt, const string &q)
|
|||
if (i < qlen && bval(q[i]) < 4)
|
||||
{ // 最后一次扩展
|
||||
c1 = cbval(q[i]);
|
||||
// double tm_t = realtime();
|
||||
fmt_extend1(fmt, &ik, &ok, 0, c1);
|
||||
ik = ok;
|
||||
// t9 += realtime() - tm_t;
|
||||
ik.info = i + 1;
|
||||
//cout << "fmt : " << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
}
|
||||
// cout << ik.x[0] << '\t' << ik.x[1] << '\t' << ik.x[2] << endl;
|
||||
return ik;
|
||||
}
|
||||
|
||||
|
|
@ -544,36 +528,14 @@ uint32_t str2bit(string &str)
|
|||
return pac;
|
||||
}
|
||||
|
||||
#define GEN_BWT_IDX 0
|
||||
#define GEN_FMT_IDX 0
|
||||
#define CMP_FMT_BWT_TIME 1
|
||||
#define CMP_FMT_BWT_RESULT 1
|
||||
|
||||
// argv[1] 应该是索引的前缀
|
||||
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;
|
||||
|
||||
// uint32_t cnt_table[4][256];
|
||||
// fmt_gen_cnt_table(cnt_table);
|
||||
// string ss1("ATATATCAATTCTCTT");
|
||||
// string ss2("ATTGCCTCTGAATTAC");
|
||||
// uint32_t bs1 = str2bit(ss1), bs2 = str2bit(ss2);
|
||||
// uint32_t s1 = bs1;
|
||||
// uint32_t x = 0;
|
||||
// for (int m = 0; m < 4; ++m)
|
||||
// {
|
||||
// x = cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24];
|
||||
// cout << ((x) >> 24 & 0xff) << ' ' << ((x) >> 16 & 0xff) << ' ' << ((x) >> 8 & 0xff) << ' ' << ((x) & 0xff) << ' ' << __fmt_mid_sum(x) << endl;
|
||||
// }
|
||||
// x = 0;
|
||||
// s1 = bs2;
|
||||
// for (int m = 0; m < 4; ++m)
|
||||
// {
|
||||
// x = cnt_table[m][s1 & 0xff] + cnt_table[m][s1 >> 8 & 0xff] + cnt_table[m][s1 >> 16 & 0xff] + cnt_table[m][s1 >> 24];
|
||||
// cout << ((x) >> 24 & 0xff) << ' ' << ((x) >> 16 & 0xff) << ' ' << ((x) >> 8 & 0xff) << ' ' << ((x) & 0xff) << ' ' << __fmt_mid_sum(x) << endl;
|
||||
// }
|
||||
|
||||
// return 0;
|
||||
string prefix = argv[1];
|
||||
string bwt_str = prefix + ".bwt.str";
|
||||
string bwt_idx, fmt_idx;
|
||||
|
|
@ -594,7 +556,7 @@ int main_fmtidx(int argc, char **argv)
|
|||
|
||||
// 生成或读取bwt索引文件
|
||||
double time_read_bwt = realtime();
|
||||
#if 0
|
||||
#if GEN_BWT_IDX
|
||||
bwt = restore_bwt(bwt_str.c_str());
|
||||
create_interval_occ_bwt(bwt);
|
||||
dump_bwt(bwt_idx.c_str(), bwt);
|
||||
|
|
@ -606,7 +568,7 @@ int main_fmtidx(int argc, char **argv)
|
|||
|
||||
// 生成或读取fmt索引文件
|
||||
double time_read_fmt = realtime();
|
||||
#if 0
|
||||
#if GEN_FMT_IDX
|
||||
fmt = create_fmt_from_bwt(bwt);
|
||||
dump_fmt(fmt_idx.c_str(), fmt);
|
||||
#else
|
||||
|
|
@ -618,9 +580,9 @@ int main_fmtidx(int argc, char **argv)
|
|||
// 生成随机序列进行测试
|
||||
int seq_size = 10000000;
|
||||
// int seq_size = 1;
|
||||
int seq_len = 19;
|
||||
int seq_len = 11;
|
||||
vector<string> seq_arr(seq_size);
|
||||
seq_arr[0] = "GCGATACTAAGA";
|
||||
// seq_arr[0] = "GCGATACTAAGA";
|
||||
// seq_arr[0] = "GAGAGCTGTTCCCGTGTTTTCCATGGTTT";
|
||||
#if 1
|
||||
srand(time(NULL));
|
||||
|
|
@ -632,7 +594,7 @@ int main_fmtidx(int argc, char **argv)
|
|||
#endif
|
||||
|
||||
// 对比bwt和fmt搜索的时间
|
||||
#if 1
|
||||
#if CMP_FMT_BWT_TIME
|
||||
double time_bwt_search = realtime();
|
||||
for (int i = 0; i < (int)seq_arr.size(); ++i)
|
||||
bwt_search(bwt, seq_arr[i]);
|
||||
|
|
@ -647,9 +609,9 @@ int main_fmtidx(int argc, char **argv)
|
|||
#endif
|
||||
|
||||
// 对比bwt和fmt搜索的结果
|
||||
#if 1
|
||||
#if CMP_FMT_BWT_RESULT
|
||||
double time_cmp = realtime();
|
||||
for (int i = 0; i < (int)seq_arr.size(); ++i)
|
||||
for (int i = 0; i < (int)seq_arr.size() / 100; ++i)
|
||||
{
|
||||
bwtintv_t p1 = bwt_search2(bwt, seq_arr[i]);
|
||||
bwtintv_t p2 = fmt_search(fmt, seq_arr[i]);
|
||||
|
|
|
|||
|
|
@ -5,7 +5,7 @@
|
|||
#define FMT_OCC_INTERVAL (1LL << FMT_OCC_INTV_SHIFT)
|
||||
#define FMT_OCC_INTV_MASK (FMT_OCC_INTERVAL - 1)
|
||||
|
||||
#define FMT_MID_INTV_SHIFT 7
|
||||
#define FMT_MID_INTV_SHIFT 5
|
||||
#define FMT_MID_INTERVAL (1LL << FMT_MID_INTV_SHIFT)
|
||||
#define FMT_MID_INTV_MASK (FMT_MID_INTERVAL - 1)
|
||||
|
||||
|
|
|
|||
11
main.cpp
11
main.cpp
|
|
@ -10,10 +10,19 @@ using namespace std;
|
|||
int main_sa(int argc, char **argv);
|
||||
int main_fmtidx(int argc, char **argv);
|
||||
|
||||
#define TEST_SA 0
|
||||
#define TEST_FMT 1
|
||||
|
||||
// 整个程序的入口
|
||||
int main(int argc, char **argv)
|
||||
{
|
||||
// main_sa(argc, argv);
|
||||
#if TEST_SA
|
||||
main_sa(argc, argv);
|
||||
#endif
|
||||
|
||||
#if TEST_FMT
|
||||
main_fmtidx(argc, argv);
|
||||
#endif
|
||||
|
||||
return 0;
|
||||
}
|
||||
16
sa.cpp
16
sa.cpp
|
|
@ -8,6 +8,7 @@
|
|||
|
||||
using namespace std;
|
||||
|
||||
// 由33bit表示的bwt行对应的所在序列位置信息
|
||||
inline void bwt_set_sa_33(uint8_t *sa_arr, bwtint_t k, bwtint_t val)
|
||||
{
|
||||
const bwtint_t block_idx = (k >> 3) * 33; // 8个数为一组，共享33个字节
|
||||
|
|
@ -66,6 +67,8 @@ static inline uint64_t get_sa_val_40(uint8_t *sa_arr, uint64_t idx) {
|
|||
return sa_val;
|
||||
}
|
||||
|
||||
#define TEST_RAND_READ 1
|
||||
|
||||
int main_sa(int argc, char **argv)
|
||||
{
|
||||
double timeRead40, timeWrite40,
|
||||
|
|
@ -73,20 +76,25 @@ int main_sa(int argc, char **argv)
|
|||
timeRead401, timeWrite401,
|
||||
timeRead64, timeWrite64;
|
||||
double timeStart;
|
||||
int saLen = 1 << 25;
|
||||
bwtint_t diffPos = 0;
|
||||
int saLen = 1 << 25; // 对应于bwt字符串长度
|
||||
bwtint_t diffPos = 0; // 用来显示结果不一致时候的出错位置
|
||||
vector<bwtint_t> valArr(saLen);
|
||||
vector<int> ri(saLen);
|
||||
uint8_t *sa33 = (uint8_t*)calloc(SA_BYTES_40(saLen), 1);
|
||||
uint8_t *sa33 = (uint8_t*)calloc(SA_BYTES_33(saLen), 1);
|
||||
uint8_t *sa40 = (uint8_t*)calloc(SA_BYTES_40(saLen), 1);
|
||||
uint8_t *sa401 = (uint8_t*)calloc(SA_BYTES_40(saLen), 1);
|
||||
bwtint_t *sa64 = (bwtint_t*)calloc(saLen, sizeof(bwtint_t));
|
||||
|
||||
for (int i=0; i< saLen; ++i) {
|
||||
valArr[i] = rand();
|
||||
valArr[i] <<= 1;
|
||||
#if TEST_RAND_READ
|
||||
ri[i] = rand() % saLen;
|
||||
// ri[i] = i;
|
||||
#else
|
||||
ri[i] = i;
|
||||
#endif
|
||||
}
|
||||
|
||||
// 33 test
|
||||
timeStart = realtime();
|
||||
for (int i=0; i<saLen; ++i) bwt_set_sa_33(sa33, i, valArr[i]);
|
||||
|
|
|
|||
Loading…
Reference in New Issue