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 ;
10 : initial-seed 42 ; inline
11 : line-length 60 ; inline
15 : random ( seed -- n seed )
16 >float IA * IC + IM mod [ IM /f ] keep ; inline
18 HINTS: random fixnum ;
20 : ALU "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA" ; inline
44 { CHAR: a 0.3029549426680 }
45 { CHAR: c 0.1979883004921 }
46 { CHAR: g 0.1975473066391 }
47 { CHAR: t 0.3015094502008 }
50 : make-cumulative ( freq -- chars floats )
52 swap values >float-array unclip [ + ] accumulate swap suffix ;
54 :: select-random ( seed chars floats -- seed elt )
55 floats seed random -rot
56 [ >= ] curry find drop
57 chars nth-unsafe ; inline
59 : make-random-fasta ( seed len chars floats -- seed )
60 [ rot drop select-random ] 2curry B{ } map-as print ; inline
62 : write-description ( desc id -- )
63 ">" write write bl print ; inline
65 :: split-lines ( n quot -- )
67 [ [ line-length quot call ] times ] dip
68 dup zero? [ drop ] quot if ; inline
70 : write-random-fasta ( seed n chars floats desc id -- seed )
72 [ make-random-fasta ] 2curry split-lines ; inline
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
80 : write-repeat-fasta ( n alu desc id -- )
82 [let | k! [ 0 ] alu [ ] |
83 [| len | k len alu make-repeat-fasta k! ] split-lines
87 homo-sapiens make-cumulative
89 [let | homo-sapiens-floats [ ]
90 homo-sapiens-chars [ ]
95 seed [ initial-seed ] |
98 n 2 * ALU "Homo sapiens alu" "ONE" write-repeat-fasta
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
108 : run-fasta ( -- ) 2500000 reverse-complement-in fasta ;