|
0
|
1 #include <stdlib.h>
|
|
|
2 #include <string.h>
|
|
|
3 #include <stdio.h>
|
|
|
4 #include "bwt_lite.h"
|
|
|
5
|
|
|
6 int is_sa(const uint8_t *T, uint32_t *SA, int n);
|
|
|
7 int is_bwt(uint8_t *T, int n);
|
|
|
8
|
|
|
9 bwtl_t *bwtl_seq2bwtl(int len, const uint8_t *seq)
|
|
|
10 {
|
|
|
11 bwtl_t *b;
|
|
|
12 int i;
|
|
|
13 b = (bwtl_t*)calloc(1, sizeof(bwtl_t));
|
|
|
14 b->seq_len = len;
|
|
|
15
|
|
|
16 { // calculate b->bwt
|
|
|
17 uint8_t *s;
|
|
|
18 b->sa = (uint32_t*)calloc(len + 1, 4);
|
|
|
19 is_sa(seq, b->sa, len);
|
|
|
20 s = (uint8_t*)calloc(len + 1, 1);
|
|
|
21 for (i = 0; i <= len; ++i) {
|
|
|
22 if (b->sa[i] == 0) b->primary = i;
|
|
|
23 else s[i] = seq[b->sa[i] - 1];
|
|
|
24 }
|
|
|
25 for (i = b->primary; i < len; ++i) s[i] = s[i + 1];
|
|
|
26 b->bwt_size = (len + 15) / 16;
|
|
|
27 b->bwt = (uint32_t*)calloc(b->bwt_size, 4);
|
|
|
28 for (i = 0; i < len; ++i)
|
|
|
29 b->bwt[i>>4] |= s[i] << ((15 - (i&15)) << 1);
|
|
|
30 free(s);
|
|
|
31 }
|
|
|
32 { // calculate b->occ
|
|
|
33 uint32_t c[4];
|
|
|
34 b->n_occ = (len + 15) / 16 * 4;
|
|
|
35 b->occ = (uint32_t*)calloc(b->n_occ, 4);
|
|
|
36 memset(c, 0, 16);
|
|
|
37 for (i = 0; i < len; ++i) {
|
|
|
38 if (i % 16 == 0)
|
|
|
39 memcpy(b->occ + (i/16) * 4, c, 16);
|
|
|
40 ++c[bwtl_B0(b, i)];
|
|
|
41 }
|
|
|
42 memcpy(b->L2+1, c, 16);
|
|
|
43 for (i = 2; i < 5; ++i) b->L2[i] += b->L2[i-1];
|
|
|
44 }
|
|
|
45 { // generate cnt_table
|
|
|
46 for (i = 0; i != 256; ++i) {
|
|
|
47 u_int32_t j, x = 0;
|
|
|
48 for (j = 0; j != 4; ++j)
|
|
|
49 x |= (((i&3) == j) + ((i>>2&3) == j) + ((i>>4&3) == j) + (i>>6 == j)) << (j<<3);
|
|
|
50 b->cnt_table[i] = x;
|
|
|
51 }
|
|
|
52 }
|
|
|
53 return b;
|
|
|
54 }
|
|
|
55 inline uint32_t bwtl_occ(const bwtl_t *bwt, uint32_t k, uint8_t c)
|
|
|
56 {
|
|
|
57 uint32_t n, b;
|
|
|
58 if (k == bwt->seq_len) return bwt->L2[c+1] - bwt->L2[c];
|
|
|
59 if (k == (uint32_t)(-1)) return 0;
|
|
|
60 if (k >= bwt->primary) --k; // because $ is not in bwt
|
|
|
61 n = bwt->occ[k/16<<2|c];
|
|
|
62 b = bwt->bwt[k/16] & ~((1U<<((15-(k&15))<<1)) - 1);
|
|
|
63 n += (bwt->cnt_table[b&0xff] + bwt->cnt_table[b>>8&0xff]
|
|
|
64 + bwt->cnt_table[b>>16&0xff] + bwt->cnt_table[b>>24]) >> (c<<3) & 0xff;
|
|
|
65 if (c == 0) n -= 15 - (k&15); // corrected for the masked bits
|
|
|
66 return n;
|
|
|
67 }
|
|
|
68 inline void bwtl_occ4(const bwtl_t *bwt, uint32_t k, uint32_t cnt[4])
|
|
|
69 {
|
|
|
70 uint32_t x, b;
|
|
|
71 if (k == (uint32_t)(-1)) {
|
|
|
72 memset(cnt, 0, 16);
|
|
|
73 return;
|
|
|
74 }
|
|
|
75 if (k >= bwt->primary) --k; // because $ is not in bwt
|
|
|
76 memcpy(cnt, bwt->occ + (k>>4<<2), 16);
|
|
|
77 b = bwt->bwt[k>>4] & ~((1U<<((~k&15)<<1)) - 1);
|
|
|
78 x = bwt->cnt_table[b&0xff] + bwt->cnt_table[b>>8&0xff]
|
|
|
79 + bwt->cnt_table[b>>16&0xff] + bwt->cnt_table[b>>24];
|
|
|
80 x -= 15 - (k&15);
|
|
|
81 cnt[0] += x&0xff; cnt[1] += x>>8&0xff; cnt[2] += x>>16&0xff; cnt[3] += x>>24;
|
|
|
82 }
|
|
|
83 inline void bwtl_2occ4(const bwtl_t *bwt, uint32_t k, uint32_t l, uint32_t cntk[4], uint32_t cntl[4])
|
|
|
84 {
|
|
|
85 bwtl_occ4(bwt, k, cntk);
|
|
|
86 bwtl_occ4(bwt, l, cntl);
|
|
|
87 }
|
|
|
88 void bwtl_destroy(bwtl_t *bwt)
|
|
|
89 {
|
|
|
90 if (bwt) {
|
|
|
91 free(bwt->occ); free(bwt->bwt); free(bwt->sa);
|
|
|
92 free(bwt);
|
|
|
93 }
|
|
|
94 }
|
