extra/benchmark/fasta/fasta.factor

   1 ! Based on http://shootout.alioth.debian.org/gp4/benchmark.php?test=fasta&lang=java&id=2
   2 USING: alien.c-types math kernel io io.files locals multiline
   3 assocs sequences sequences.private benchmark.reverse-complement
   4 hints io.encodings.ascii byte-arrays specialized-arrays ;
   5 SPECIALIZED-ARRAY: double
   6 IN: benchmark.fasta
   7
   8 CONSTANT: IM 139968
   9 CONSTANT: IA 3877
  10 CONSTANT: IC 29573
  11 CONSTANT: initial-seed 42
  12 CONSTANT: line-length 60
  13
  14 : random ( seed -- n seed )
  15     >float IA * IC + IM mod [ IM /f ] keep ; inline
  16
  17 HINTS: random fixnum ;
  18
  19 CONSTANT: ALU "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"
  20
  21 CONSTANT: IUB
  22     {
  23         { CHAR: a 0.27 }
  24         { CHAR: c 0.12 }
  25         { CHAR: g 0.12 }
  26         { CHAR: t 0.27 }
  27
  28         { CHAR: B 0.02 }
  29         { CHAR: D 0.02 }
  30         { CHAR: H 0.02 }
  31         { CHAR: K 0.02 }
  32         { CHAR: M 0.02 }
  33         { CHAR: N 0.02 }
  34         { CHAR: R 0.02 }
  35         { CHAR: S 0.02 }
  36         { CHAR: V 0.02 }
  37         { CHAR: W 0.02 }
  38         { CHAR: Y 0.02 }
  39     }
  40
  41 CONSTANT: homo-sapiens
  42     {
  43         { CHAR: a 0.3029549426680 }
  44         { CHAR: c 0.1979883004921 }
  45         { CHAR: g 0.1975473066391 }
  46         { CHAR: t 0.3015094502008 }
  47     }
  48
  49 : make-cumulative ( freq -- chars floats )
  50     [ keys >byte-array ]
  51     [ values >double-array ] bi unclip [ + ] accumulate swap suffix ;
  52
  53 :: select-random ( seed chars floats -- seed elt )
  54     floats seed random -rot
  55     [ >= ] curry find drop
  56     chars nth-unsafe ; inline
  57
  58 : make-random-fasta ( seed len chars floats -- seed )
  59     [ rot drop select-random ] 2curry "" map-as print ; inline
  60
  61 : write-description ( desc id -- )
  62     ">" write write bl print ; inline
  63
  64 :: split-lines ( n quot -- )
  65     n line-length /mod
  66     [ [ line-length quot call ] times ] dip
  67     quot unless-zero ; inline
  68
  69 : write-random-fasta ( seed n chars floats desc id -- seed )
  70     write-description
  71     [ make-random-fasta ] 2curry split-lines ; inline
  72
  73 :: make-repeat-fasta ( k len alu -- k' )
  74     alu length :> kn
  75     len [ k + kn mod alu nth-unsafe ] "" map-as print
  76     k len + ; inline
  77
  78 : write-repeat-fasta ( n alu desc id -- )
  79     write-description
  80     [let
  81         :> alu
  82         0 :> k!
  83         [| len | k len alu make-repeat-fasta k! ] split-lines
  84     ] ; inline
  85
  86 : fasta ( n out -- )
  87     homo-sapiens make-cumulative
  88     IUB make-cumulative
  89     [let
  90         :> ( n out IUB-chars IUB-floats homo-sapiens-chars homo-sapiens-floats )
  91         initial-seed :> seed
  92
  93         out ascii [
  94             n 2 * ALU "Homo sapiens alu" "ONE" write-repeat-fasta
  95
  96             initial-seed
  97             n 3 * homo-sapiens-chars homo-sapiens-floats
  98             "IUB ambiguity codes" "TWO" write-random-fasta
  99             n 5 * IUB-chars IUB-floats
 100             "Homo sapiens frequency" "THREE" write-random-fasta
 101             drop
 102         ] with-file-writer
 103     ] ;
 104
 105 : run-fasta ( -- ) 2500000 reverse-complement-in fasta ;
 106
 107 MAIN: run-fasta