src/test/bench/shootout-fasta-redux.rs

   1 // The Computer Language Benchmarks Game
   2 // http://benchmarksgame.alioth.debian.org/
   3 //
   4 // contributed by the Rust Project Developers
   5
   6 // Copyright (c) 2013-2014 The Rust Project Developers
   7 //
   8 // All rights reserved.
   9 //
  10 // Redistribution and use in source and binary forms, with or without
  11 // modification, are permitted provided that the following conditions
  12 // are met:
  13 //
  14 // - Redistributions of source code must retain the above copyright
  15 //   notice, this list of conditions and the following disclaimer.
  16 //
  17 // - Redistributions in binary form must reproduce the above copyright
  18 //   notice, this list of conditions and the following disclaimer in
  19 //   the documentation and/or other materials provided with the
  20 //   distribution.
  21 //
  22 // - Neither the name of "The Computer Language Benchmarks Game" nor
  23 //   the name of "The Computer Language Shootout Benchmarks" nor the
  24 //   names of its contributors may be used to endorse or promote
  25 //   products derived from this software without specific prior
  26 //   written permission.
  27 //
  28 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
  29 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
  30 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
  31 // FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
  32 // COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
  33 // INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
  34 // (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
  35 // SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
  36 // HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
  37 // STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
  38 // ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
  39 // OF THE POSSIBILITY OF SUCH DAMAGE.
  40
  41 use std::cmp::min;
  42 use std::env;
  43 use std::io;
  44 use std::io::prelude::*;
  45
  46 const LINE_LEN: usize = 60;
  47 const LOOKUP_SIZE: usize = 4 * 1024;
  48 const LOOKUP_SCALE: f32 = (LOOKUP_SIZE - 1) as f32;
  49
  50 // Random number generator constants
  51 const IM: u32 = 139968;
  52 const IA: u32 = 3877;
  53 const IC: u32 = 29573;
  54
  55 const ALU: &'static str = "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTG\
  56                             GGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGA\
  57                             GACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAA\
  58                             AATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAAT\
  59                             CCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAAC\
  60                             CCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTG\
  61                             CACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";
  62
  63 const NULL_AMINO_ACID: AminoAcid = AminoAcid { c: ' ' as u8, p: 0.0 };
  64
  65 static IUB: [AminoAcid;15] = [
  66     AminoAcid { c: 'a' as u8, p: 0.27 },
  67     AminoAcid { c: 'c' as u8, p: 0.12 },
  68     AminoAcid { c: 'g' as u8, p: 0.12 },
  69     AminoAcid { c: 't' as u8, p: 0.27 },
  70     AminoAcid { c: 'B' as u8, p: 0.02 },
  71     AminoAcid { c: 'D' as u8, p: 0.02 },
  72     AminoAcid { c: 'H' as u8, p: 0.02 },
  73     AminoAcid { c: 'K' as u8, p: 0.02 },
  74     AminoAcid { c: 'M' as u8, p: 0.02 },
  75     AminoAcid { c: 'N' as u8, p: 0.02 },
  76     AminoAcid { c: 'R' as u8, p: 0.02 },
  77     AminoAcid { c: 'S' as u8, p: 0.02 },
  78     AminoAcid { c: 'V' as u8, p: 0.02 },
  79     AminoAcid { c: 'W' as u8, p: 0.02 },
  80     AminoAcid { c: 'Y' as u8, p: 0.02 },
  81 ];
  82
  83 static HOMO_SAPIENS: [AminoAcid;4] = [
  84     AminoAcid { c: 'a' as u8, p: 0.3029549426680 },
  85     AminoAcid { c: 'c' as u8, p: 0.1979883004921 },
  86     AminoAcid { c: 'g' as u8, p: 0.1975473066391 },
  87     AminoAcid { c: 't' as u8, p: 0.3015094502008 },
  88 ];
  89
  90 // FIXME: Use map().
  91 fn sum_and_scale(a: &'static [AminoAcid]) -> Vec<AminoAcid> {
  92     let mut result = Vec::new();
  93     let mut p = 0f32;
  94     for a_i in a {
  95         let mut a_i = *a_i;
  96         p += a_i.p;
  97         a_i.p = p * LOOKUP_SCALE;
  98         result.push(a_i);
  99     }
 100     let result_len = result.len();
 101     result[result_len - 1].p = LOOKUP_SCALE;
 102     result
 103 }
 104
 105 #[derive(Copy, Clone)]
 106 struct AminoAcid {
 107     c: u8,
 108     p: f32,
 109 }
 110
 111 struct RepeatFasta<'a, W:'a> {
 112     alu: &'static str,
 113     out: &'a mut W
 114 }
 115
 116 impl<'a, W: Write> RepeatFasta<'a, W> {
 117     fn new(alu: &'static str, w: &'a mut W) -> RepeatFasta<'a, W> {
 118         RepeatFasta { alu: alu, out: w }
 119     }
 120
 121     fn make(&mut self, n: usize) -> io::Result<()> {
 122         let alu_len = self.alu.len();
 123         let mut buf = vec![0; alu_len + LINE_LEN];
 124         let alu: &[u8] = self.alu.as_bytes();
 125
 126         for (slot, val) in buf.iter_mut().zip(alu) {
 127             *slot = *val;
 128         }
 129         let buf_len = buf.len();
 130         for (slot, val) in buf[alu_len..buf_len].iter_mut().zip(&alu[..LINE_LEN]) {
 131             *slot = *val;
 132         }
 133
 134         let mut pos = 0;
 135         let mut bytes;
 136         let mut n = n;
 137         while n > 0 {
 138             bytes = min(LINE_LEN, n);
 139             try!(self.out.write_all(&buf[pos..pos + bytes]));
 140             try!(self.out.write_all(&[b'\n']));
 141             pos += bytes;
 142             if pos > alu_len {
 143                 pos -= alu_len;
 144             }
 145             n -= bytes;
 146         }
 147         Ok(())
 148     }
 149 }
 150
 151 fn make_lookup(a: &[AminoAcid]) -> [AminoAcid;LOOKUP_SIZE] {
 152     let mut lookup = [ NULL_AMINO_ACID;LOOKUP_SIZE ];
 153     let mut j = 0;
 154     for (i, slot) in lookup.iter_mut().enumerate() {
 155         while a[j].p < (i as f32) {
 156             j += 1;
 157         }
 158         *slot = a[j];
 159     }
 160     lookup
 161 }
 162
 163 struct RandomFasta<'a, W:'a> {
 164     seed: u32,
 165     lookup: [AminoAcid;LOOKUP_SIZE],
 166     out: &'a mut W,
 167 }
 168
 169 impl<'a, W: Write> RandomFasta<'a, W> {
 170     fn new(w: &'a mut W, a: &[AminoAcid]) -> RandomFasta<'a, W> {
 171         RandomFasta {
 172             seed: 42,
 173             out: w,
 174             lookup: make_lookup(a),
 175         }
 176     }
 177
 178     fn rng(&mut self, max: f32) -> f32 {
 179         self.seed = (self.seed * IA + IC) % IM;
 180         max * (self.seed as f32) / (IM as f32)
 181     }
 182
 183     fn nextc(&mut self) -> u8 {
 184         let r = self.rng(1.0);
 185         for a in &self.lookup[..] {
 186             if a.p >= r {
 187                 return a.c;
 188             }
 189         }
 190         0
 191     }
 192
 193     fn make(&mut self, n: usize) -> io::Result<()> {
 194         let lines = n / LINE_LEN;
 195         let chars_left = n % LINE_LEN;
 196         let mut buf = [0;LINE_LEN + 1];
 197
 198         for _ in 0..lines {
 199             for i in 0..LINE_LEN {
 200                 buf[i] = self.nextc();
 201             }
 202             buf[LINE_LEN] = '\n' as u8;
 203             try!(self.out.write(&buf));
 204         }
 205         for i in 0..chars_left {
 206             buf[i] = self.nextc();
 207         }
 208         self.out.write_all(&buf[..chars_left])
 209     }
 210 }
 211
 212 fn main() {
 213     let mut args = env::args();
 214     let n = if args.len() > 1 {
 215         args.nth(1).unwrap().parse::<usize>().unwrap()
 216     } else {
 217         5
 218     };
 219
 220     let mut out = io::stdout();
 221
 222     out.write_all(b">ONE Homo sapiens alu\n").unwrap();
 223     {
 224         let mut repeat = RepeatFasta::new(ALU, &mut out);
 225         repeat.make(n * 2).unwrap();
 226     }
 227
 228     out.write_all(b">TWO IUB ambiguity codes\n").unwrap();
 229     let iub = sum_and_scale(&IUB);
 230     let mut random = RandomFasta::new(&mut out, &iub);
 231     random.make(n * 3).unwrap();
 232
 233     random.out.write_all(b">THREE Homo sapiens frequency\n").unwrap();
 234     let homo_sapiens = sum_and_scale(&HOMO_SAPIENS);
 235     random.lookup = make_lookup(&homo_sapiens);
 236     random.make(n * 5).unwrap();
 237
 238     random.out.write_all(b"\n").unwrap();
 239 }