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
11 CONSTANT: initial-seed 42
12 CONSTANT: line-length 60
14 : random ( seed -- n seed )
15 >float IA * IC + IM mod [ IM /f ] keep ; inline
17 HINTS: random fixnum ;
19 CONSTANT: ALU "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGGGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"
41 CONSTANT: homo-sapiens
43 { CHAR: a 0.3029549426680 }
44 { CHAR: c 0.1979883004921 }
45 { CHAR: g 0.1975473066391 }
46 { CHAR: t 0.3015094502008 }
49 : make-cumulative ( freq -- chars floats )
51 [ values >double-array ] bi unclip [ + ] accumulate swap suffix ;
53 :: select-random ( seed chars floats -- seed elt )
54 floats seed random -rot
55 [ >= ] curry find drop
56 chars nth-unsafe ; inline
58 : make-random-fasta ( seed len chars floats -- seed )
59 [ iota ] 2dip [ [ drop ] 2dip select-random ] 2curry "" map-as print ; inline
61 : write-description ( desc id -- )
62 ">" write write bl print ; inline
64 :: split-lines ( n quot -- )
66 [ [ line-length quot call ] times ] dip
67 quot unless-zero ; inline
69 : write-random-fasta ( seed n chars floats desc id -- seed )
71 [ make-random-fasta ] 2curry split-lines ; inline
73 :: make-repeat-fasta ( k len alu -- k' )
75 len iota [ k + kn mod alu nth-unsafe ] "" map-as print
78 : write-repeat-fasta ( n alu desc id -- )
83 [| len | k len alu make-repeat-fasta k! ] split-lines
87 homo-sapiens make-cumulative
90 :> ( n out IUB-chars IUB-floats homo-sapiens-chars homo-sapiens-floats )
94 n 2 * ALU "Homo sapiens alu" "ONE" write-repeat-fasta
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
105 : run-fasta ( -- ) 2500000 reverse-complement-in fasta ;