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 [ rot drop 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' )
74 [let | kn [ alu length ] |
75 len [ k + kn mod alu nth-unsafe ] "" map-as print
79 : write-repeat-fasta ( n alu desc id -- )
81 [let | k! [ 0 ] alu [ ] |
82 [| len | k len alu make-repeat-fasta k! ] split-lines
86 homo-sapiens make-cumulative
88 [let | homo-sapiens-floats [ ]
89 homo-sapiens-chars [ ]
94 seed [ initial-seed ] |
97 n 2 * ALU "Homo sapiens alu" "ONE" write-repeat-fasta
100 n 3 * homo-sapiens-chars homo-sapiens-floats "IUB ambiguity codes" "TWO" write-random-fasta
101 n 5 * IUB-chars IUB-floats "Homo sapiens frequency" "THREE" write-random-fasta
107 : run-fasta ( -- ) 2500000 reverse-complement-in fasta ;