fasta

# The Computer Language Benchmarks Game
# http://shootout.alioth.debian.org/
#
# by Limin Fu

global  alu  =
'GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG'
'GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA'
'CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT'
'ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA'
'GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG'
'AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC'
'AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA';

typedef  tuple<string,list<float>  >  FreqTable;

global  iub  =  (  'acgtBDHKMNRSVWY',  {  0.27,  0.12,  0.12,  0.27  }  +  {  0.02  :  0.0  :  11  }  );

global  homosapiens  = 
(  'acgt',  {  0.3029549426680,  0.1979883004921,  0.1975473066391,  0.3015094502008  }  );

routine  genRandom(  lim  :  float  )
{
     ia  =  3877;
     ic  =  29573;
     im  =  139968;
     imf  =  139968.0;
     seed  =  42;
     while(  1  ){
          seed  =  (seed  *  ia  +  ic)  %  im;
          yield  (lim  *  seed  /  imf);
     }
     return  0.0;
}
global  gen_random  =  @genRandom(1.0);

routine  makeCumulative(  table  :  list<float>  )
{
     for(  i  =  1  :  table.size()  -1  )  table[i]  +=  table[i-1];
}
routine  selectRandom  (  table  :  FreqTable  )
{
     r  =  gen_random  (1);
     (  syms,  freqs  )  =  table;
     if  (r  <  freqs[0])  return  syms[0];
     lo  =  0;
     hi  =  freqs.size()-1;
     while  (hi  >  lo+1)  {
          i  =  (hi  +  lo)  /  2;
          if  (r  <  freqs[i])  hi  =  i;  else  lo  =  i;
     }
     return  syms[hi];
}
routine  randomFasta  (  table  :  FreqTable,  n  =  0  )
{
     line  =  60;
     pick  =  (string)  {  0  :  line  };
     makeCumulative(  table[1]  );
      for  (  i  =  0  :  line  :  n-1  )  {
          m  =  (  i  +  line  <  n  )  ?  line  :  n  -  i;
          for  (j  =  0  :  m-1  )  pick[j]  =  selectRandom  (  table  );
          if(  m  <  line  )  pick  =  pick[  :  m-1  ];
          stdio.println  (pick);
      }
}
routine  repeatFasta(  src='',  n=0  )
{
     line  =  60;
     r  =  src.size();
     s  =  src  +  src;
     for(  j  =  0  :  line  :  n-1){
          m  =  (  j  +  line  <  n  )  ?  line  :  n  -  j;
          stdio.println(  s[  (j  %  r)  :  (j  %  r)  +  m-1  ]  );
     }
}
routine  main(  n  =  100  )
{
     stdio.println(  '>ONE Homo sapiens alu'  );
     repeatFasta(alu,  n*2  );
     stdio.println(  '>TWO IUB ambiguity codes'  );
     randomFasta(iub,  n*3);
     stdio.println(  '>THREE Homo sapiens frequency'  );
     randomFasta(homosapiens,  n*5);
}