learn_bwt_indexing_vanilla.c

## learn_bwt_indexing_vanilla.c
#include <stdio.h>
#include <stdlib.h>
#include <string.h>

#define L 14 // Length of the text.

// Global variables.
char TXT[L]     = "GATGCGAGAGATG$";
int  SA[L]      = {0,1,2,3,4,5,6,7,8,9,10,11,12,13};
char BWT[L]     = {0};
int  C[4]       = {0};
int  OCC[4][L]  = {0};

int NUM[256] = {
  ['A'] = 0, ['C'] = 1, ['G'] = 2, ['T'] = 3,
  ['a'] = 0, ['c'] = 1, ['g'] = 2, ['t'] = 3,
};

int compare_suffixes (const void * a, const void * b) {
   const char * suff_a = &TXT[*(int *)a];
   const char * suff_b = &TXT[*(int *)b];
   return strcmp(suff_a, suff_b);
}

void construct_BWT (void) {
   for (int i = 0 ; i < L ; i++) {
      if (SA[i] > 0) BWT[i] = TXT[SA[i]-1];
      else           BWT[i] = '$';
   }
}

void construct_C_and_OCC (void) {
   for (int i = 0 ; i < L ; i++) {
      if (BWT[i] != '$') C[NUM[BWT[i]]]++;
      for (int j = 0 ; j < 4 ; j++) OCC[j][i] = C[j];
   }
   int tmp = 1;
   for (int j = 0 ; j < 4 ; j++) {
      tmp += C[j];
      C[j] = tmp - C[j];
   }
}

void backward_search (char * query) {
   int bot = 1;
   int top = L-1;
   for (int pos = strlen(query)-1 ; pos > -1 ; pos--) {
      int c = NUM[query[pos]];
      bot = C[c] + OCC[c][bot-1];
      top = C[c] + OCC[c][top]-1;
      if (top < bot) break;
   }
   for (int i = bot ; i <= top ; i++) {
      printf("%s: %d\n", query, SA[i]);
   }
}

int main(void) {
   qsort(SA, L, sizeof(int), compare_suffixes);
   construct_BWT();
   construct_C_and_OCC();
   backward_search("GAGA");
}