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 specialized-arrays.double ;
10 CONSTANT: initial-seed 42
11 CONSTANT: line-length 60
13 : random ( seed -- n seed )
14 >float IA * IC + IM mod [ IM /f ] keep ; inline
16 HINTS: random fixnum ;
18 CONSTANT: ALU "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"
40 CONSTANT: homo-sapiens
42 { CHAR: a 0.3029549426680 }
43 { CHAR: c 0.1979883004921 }
44 { CHAR: g 0.1975473066391 }
45 { CHAR: t 0.3015094502008 }
48 : make-cumulative ( freq -- chars floats )
50 [ values >double-array ] bi unclip [ + ] accumulate swap suffix ;
52 :: select-random ( seed chars floats -- seed elt )
53 floats seed random -rot
54 [ >= ] curry find drop
55 chars nth-unsafe ; inline
57 : make-random-fasta ( seed len chars floats -- seed )
58 [ rot drop select-random ] 2curry "" map-as print ; inline
60 : write-description ( desc id -- )
61 ">" write write bl print ; inline
63 :: split-lines ( n quot -- )
65 [ [ line-length quot call ] times ] dip
66 dup zero? [ drop ] quot if ; inline
68 : write-random-fasta ( seed n chars floats desc id -- seed )
70 [ make-random-fasta ] 2curry split-lines ; inline
72 :: make-repeat-fasta ( k len alu -- k' )
73 [let | kn [ alu length ] |
74 len [ k + kn mod alu nth-unsafe ] "" map-as print
78 : write-repeat-fasta ( n alu desc id -- )
80 [let | k! [ 0 ] alu [ ] |
81 [| len | k len alu make-repeat-fasta k! ] split-lines
85 homo-sapiens make-cumulative
87 [let | homo-sapiens-floats [ ]
88 homo-sapiens-chars [ ]
93 seed [ initial-seed ] |
96 n 2 * ALU "Homo sapiens alu" "ONE" write-repeat-fasta
99 n 3 * homo-sapiens-chars homo-sapiens-floats "IUB ambiguity codes" "TWO" write-random-fasta
100 n 5 * IUB-chars IUB-floats "Homo sapiens frequency" "THREE" write-random-fasta
106 : run-fasta ( -- ) 2500000 reverse-complement-in fasta ;