extra/benchmark/fasta/fasta.factor

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