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