bench/yarv/bm_so_fasta.rb

   1 # The Computer Language Shootout\r
   2 # http://shootout.alioth.debian.org/\r
   3 # Contributed by Sokolov Yura\r
   4 \r
   5 $last = 42.0\r
   6 def gen_random (max,im=139968,ia=3877,ic=29573)\r
   7     (max * ($last = ($last * ia + ic) % im)) / im\r
   8 end\r
   9 \r
  10 alu =\r
  11    "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG"+\r
  12    "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA"+\r
  13    "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT"+\r
  14    "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA"+\r
  15    "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG"+\r
  16    "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC"+\r
  17    "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA"\r
  18 \r
  19 iub = [\r
  20     ["a", 0.27],\r
  21     ["c", 0.12],\r
  22     ["g", 0.12],\r
  23     ["t", 0.27],\r
  24 \r
  25     ["B", 0.02],\r
  26     ["D", 0.02],\r
  27     ["H", 0.02],\r
  28     ["K", 0.02],\r
  29     ["M", 0.02],\r
  30     ["N", 0.02],\r
  31     ["R", 0.02],\r
  32     ["S", 0.02],\r
  33     ["V", 0.02],\r
  34     ["W", 0.02],\r
  35     ["Y", 0.02],\r
  36 ]\r
  37 homosapiens = [\r
  38     ["a", 0.3029549426680],\r
  39     ["c", 0.1979883004921],\r
  40     ["g", 0.1975473066391],\r
  41     ["t", 0.3015094502008],\r
  42 ]\r
  43 \r
  44 def make_repeat_fasta(id, desc, src, n)\r
  45     puts ">#{id} #{desc}"\r
  46     v = nil\r
  47     width = 60\r
  48     l = src.length\r
  49     s = src * ((n / l) + 1)\r
  50     s.slice!(n, l)\r
  51     puts(s.scan(/.{1,#{width}}/).join("\n"))\r
  52 end\r
  53 \r
  54 def make_random_fasta(id, desc, table, n)\r
  55     puts ">#{id} #{desc}"\r
  56     rand, v = nil,nil\r
  57     width = 60\r
  58     chunk = 1 * width\r
  59     prob = 0.0\r
  60     table.each{|v| v[1]= (prob += v[1])}\r
  61     for i in 1..(n/width)\r
  62         puts((1..width).collect{\r
  63             rand = gen_random(1.0)\r
  64             table.find{|v| v[1]>rand}[0]\r
  65         }.join)\r
  66     end\r
  67     if n%width != 0\r
  68         puts((1..(n%width)).collect{\r
  69             rand = gen_random(1.0)\r
  70             table.find{|v| v[1]>rand}[0]\r
  71         }.join)\r
  72     end\r
  73 end\r
  74 \r
  75 \r
  76 n = (ARGV[0] or 250_000).to_i\r
  77 \r
  78 make_repeat_fasta('ONE', 'Homo sapiens alu', alu, n*2)\r
  79 make_random_fasta('TWO', 'IUB ambiguity codes', iub, n*3)\r
  80 make_random_fasta('THREE', 'Homo sapiens frequency', homosapiens, n*5)\r
  81 \r