//------------------------------------------------------------- // Induced Sorting // Reference: Nong, Ge, Sen Zhang, and Wai Hong Chan. "Two efficient algorithms for linear time suffix array construction.” IEEE Transactions on Computers, 60.10 (2011): 1471-1484. // To compile, perform: gcc -std=c99 SAIS.c -o SAIS // To run, SALS . e.g., SAIS 100 //------------------------------------------------------------- #include #include #include #include #include "MT.h" // Use the Mersenne Twister. //#define DEBUG // これを定義すると #ifdef #endif の中の が実行されるので、デバッグの時だけ動かしたいコードを入れるのに便利です. デバッグしないときはコメントではずしましょう. unsigned char mask[]={0x80, 0x40, 0x20, 0x10, 0x08, 0x04, 0x02, 0x01}; // S と L から構成される列を8文字毎に区切り、tyope S を 1, L を 0 でビット表現して 8ビット(1バイト)とし、配列 t に入れる. たとえば LSLSSLSS は 01011011 にコードする(主記憶を節約する本格的なプログラミング技術) // 配列 t から i 番目の文字の type b を取り出すマクロ // ビット演算 or演算は |, and 演算は &, ビットの反転は ~ である #define tget(i) ( (t[(i)/8] & mask[(i)%8]) ? 1 : 0 ) // i 番目の文字の type が b であることを配列 t に入れるマクロ #define tset(i, b) t[(i)/8] = (b) ? (mask[(i)%8] | t[(i)/8]) : ((~mask[(i)%8]) & t[(i)/8]) // 配列 s からのデータの取り出しを cs が sizeof(int) なら整数として s[i] を取り出し、そうでなければ unsigned char として取り出すマクロ #define chr(i) ( cs==sizeof(int) ? ((int*)s)[i] : ((unsigned char *)s)[i] ) // i 番目の文字がleft-most S-type (LMS) ならば true、そうでなければ false を返すマクロ #define isLMS(i) (i>0 && tget(i) && !tget(i-1)) // find the start or end of each bucket void getBuckets(unsigned char *s, int *bkt, int n, int K, int cs, bool end) { int i, sum=0; // clear all buckets. K は $ を除いた文字の数 for(i=0; i<=K; i++) bkt[i]=0; // compute the size of each bucket for(i=0; i=0 && !tget(j)) SA[bkt[chr(j)]++]=j; // j 番目の文字が L-type (0) ならば、 ^ が置かれた位置 bkt[chr(j)] に j を入れて ^ を右へ1つ移動 } } // compute SAs S-type の文字をソート void induceSAs(unsigned char *t, int *SA, unsigned char *s, int *bkt, int n, int K, int cs, bool end) { int i, j; // find ends of buckets 配列 bkt に最後の位置を入れる getBuckets(s, bkt, n, K, cs, end); for(i=n-1; i>=0; i--) { // i は @ の置かれた位置を表現 j=SA[i]-1; if(j>=0 && tget(j)) SA[--bkt[chr(j)]]=j; // j 番目の文字が S-type (1) ならば、 ^ が置かれた位置 bkt[chr(j)] を1つ左に移動して j を入れる } } // find the suffix array SA of s[0..n-1] in {1..K}n require s[n-1]=0 (the sentinel!), n>=2 // use a working space (excluding s and SA) of at most 2.25n+O(1) for a constant alphabet void SA_IS(unsigned char *s, int *SA, int n, int K, int cs) { // LS-type array in bits タイプ S もしくは L をビット表現して記載する配列を t とする unsigned char *t=(unsigned char *)malloc(n/8+1); int i, j; // classify the type of each character the sentinel must be in s1, important!!! // 各文字が S-type か L-type のどちらであるかを計算して配列 t に入れる tset(n-2, 0); tset(n-1, 1); // 最後から2番めの文字は 0 (L-type) で、最後の文字 $ は 1 (S-type) for(i=n-3; i>=0; i--) // 後ろの文字から埋めてゆく tset(i, (chr(i)0 && (isLMS(pos+d) || isLMS(prev+d))) break; // 次の LMS prefix まで到達し、一致したので diff = falseのまま抜ける if(diff) { name++; prev=pos; } // LMS prefix が一致しなければ新しい順番を name に設定 pos=(pos%2==0) ? pos/2 : (pos-1)/2; // LMS は少なくとも1つおきに出現するので場所をつめる SA[n1+pos]=name-1; // 更新した name の1つ前の順番を付与する } for(i=n-1, j=n-1; i>=n1; i--) // -1 を除きながら LMS prefix をコードした文字列を SA の後半に詰める if(SA[i]>=0) SA[j--]=SA[i]; // solve the reduced problem // recurse if names are not yet unique // s1 は LMS prefix をコードした文字列で SA の後半を使うのに対して、その suffix array SA1 は SA の前半に計算することで主記憶を節約する(主記憶を節約する本格的なプログラミング技術) int *SA1=SA, *s1=SA+n-n1; // 同一順序の LMS prefix が存在するならば、LMS prefix の総数 n1 より name が小さくなるので SA_IS を再帰呼び出しして s1 の suffix array SA1 を決定 if(name=0; i--) { j=SA[i]; SA[i]=-1; SA[--bkt[chr(j)]]=j; // LMS suffix の順番を bucket に入れる } // LMS suffix の順番が bucket 内に入ったので L-type をソートして後に S-type をソートする induceSAl(t, SA, s, bkt, n, K, cs, false); induceSAs(t, SA, s, bkt, n, K, cs, true); free(bkt); free(t); } int less(unsigned char *S, int aI, int aJ, int n){ int i, j; for(i=aI, j=aJ; i S[j]) return(0); } if(aI < aJ) return(0); else return(1); } int main(int argc, char *argv[]) { if(argc != 2){ fprintf(stderr, "Usage: SAIS \n"); exit(EXIT_FAILURE); } int n = atoi(argv[1]); if(n <= 0){ fprintf(stderr, "The argument must be a positive integer\n"); exit(EXIT_FAILURE); } unsigned char *S = (unsigned char *)malloc(sizeof(unsigned char) * n); // Use the heap for the input string. if(S == NULL) exit(EXIT_FAILURE); int *SA = (int *)malloc(sizeof(int)*n); if(SA == NULL){ free(S); exit(EXIT_FAILURE); } // We encode $,A,C,G,T by 0,1,2,3,4 so that $=0 < A=1 < C=2 < G=3 < T=4 int K = 4; for(int i=0; i