fasta.go

Documentation: golang.org/x/exp/shootout

     1  //go:build ignore
     2  // +build ignore
     3  
     4  /*
     5  Redistribution and use in source and binary forms, with or without
     6  modification, are permitted provided that the following conditions are met:
     7  
     8      * Redistributions of source code must retain the above copyright
     9      notice, this list of conditions and the following disclaimer.
    10  
    11      * Redistributions in binary form must reproduce the above copyright
    12      notice, this list of conditions and the following disclaimer in the
    13      documentation and/or other materials provided with the distribution.
    14  
    15      * Neither the name of "The Computer Language Benchmarks Game" nor the
    16      name of "The Computer Language Shootout Benchmarks" nor the names of
    17      its contributors may be used to endorse or promote products derived
    18      from this software without specific prior written permission.
    19  
    20  THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
    21  AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
    22  IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
    23  ARE DISCLAIMED.  IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE
    24  LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
    25  CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
    26  SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
    27  INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
    28  CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
    29  ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
    30  POSSIBILITY OF SUCH DAMAGE.
    31  */
    32  
    33  /* The Computer Language Benchmarks Game
    34   * http://shootout.alioth.debian.org/
    35   *
    36   * contributed by The Go Authors.
    37   * Based on C program by by Petr Prokhorenkov.
    38   */
    39  
    40  package main
    41  
    42  import (
    43  	"flag"
    44  	"os"
    45  )
    46  
    47  var out = make(buffer, 0, 32768)
    48  
    49  var n = flag.Int("n", 1000, "length of result")
    50  
    51  const Line = 60
    52  
    53  func Repeat(alu []byte, n int) {
    54  	buf := append(alu, alu...)
    55  	off := 0
    56  	for n > 0 {
    57  		m := n
    58  		if m > Line {
    59  			m = Line
    60  		}
    61  		buf1 := out.NextWrite(m + 1)
    62  		copy(buf1, buf[off:])
    63  		buf1[m] = '\n'
    64  		if off += m; off >= len(alu) {
    65  			off -= len(alu)
    66  		}
    67  		n -= m
    68  	}
    69  }
    70  
    71  const (
    72  	IM = 139968
    73  	IA = 3877
    74  	IC = 29573
    75  
    76  	LookupSize          = 4096
    77  	LookupScale float64 = LookupSize - 1
    78  )
    79  
    80  var rand uint32 = 42
    81  
    82  type Acid struct {
    83  	sym   byte
    84  	prob  float64
    85  	cprob float64
    86  	next  *Acid
    87  }
    88  
    89  func computeLookup(acid []Acid) *[LookupSize]*Acid {
    90  	var lookup [LookupSize]*Acid
    91  	var p float64
    92  	for i := range acid {
    93  		p += acid[i].prob
    94  		acid[i].cprob = p * LookupScale
    95  		if i > 0 {
    96  			acid[i-1].next = &acid[i]
    97  		}
    98  	}
    99  	acid[len(acid)-1].cprob = 1.0 * LookupScale
   100  
   101  	j := 0
   102  	for i := range lookup {
   103  		for acid[j].cprob < float64(i) {
   104  			j++
   105  		}
   106  		lookup[i] = &acid[j]
   107  	}
   108  
   109  	return &lookup
   110  }
   111  
   112  func Random(acid []Acid, n int) {
   113  	lookup := computeLookup(acid)
   114  	for n > 0 {
   115  		m := n
   116  		if m > Line {
   117  			m = Line
   118  		}
   119  		buf := out.NextWrite(m + 1)
   120  		f := LookupScale / IM
   121  		myrand := rand
   122  		for i := 0; i < m; i++ {
   123  			myrand = (myrand*IA + IC) % IM
   124  			r := float64(int(myrand)) * f
   125  			a := lookup[int(r)]
   126  			for a.cprob < r {
   127  				a = a.next
   128  			}
   129  			buf[i] = a.sym
   130  		}
   131  		rand = myrand
   132  		buf[m] = '\n'
   133  		n -= m
   134  	}
   135  }
   136  
   137  func main() {
   138  	defer out.Flush()
   139  
   140  	flag.Parse()
   141  
   142  	iub := []Acid{
   143  		{prob: 0.27, sym: 'a'},
   144  		{prob: 0.12, sym: 'c'},
   145  		{prob: 0.12, sym: 'g'},
   146  		{prob: 0.27, sym: 't'},
   147  		{prob: 0.02, sym: 'B'},
   148  		{prob: 0.02, sym: 'D'},
   149  		{prob: 0.02, sym: 'H'},
   150  		{prob: 0.02, sym: 'K'},
   151  		{prob: 0.02, sym: 'M'},
   152  		{prob: 0.02, sym: 'N'},
   153  		{prob: 0.02, sym: 'R'},
   154  		{prob: 0.02, sym: 'S'},
   155  		{prob: 0.02, sym: 'V'},
   156  		{prob: 0.02, sym: 'W'},
   157  		{prob: 0.02, sym: 'Y'},
   158  	}
   159  
   160  	homosapiens := []Acid{
   161  		{prob: 0.3029549426680, sym: 'a'},
   162  		{prob: 0.1979883004921, sym: 'c'},
   163  		{prob: 0.1975473066391, sym: 'g'},
   164  		{prob: 0.3015094502008, sym: 't'},
   165  	}
   166  
   167  	alu := []byte(
   168  		"GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG" +
   169  			"GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA" +
   170  			"CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT" +
   171  			"ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA" +
   172  			"GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG" +
   173  			"AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC" +
   174  			"AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA")
   175  
   176  	out.WriteString(">ONE Homo sapiens alu\n")
   177  	Repeat(alu, 2**n)
   178  	out.WriteString(">TWO IUB ambiguity codes\n")
   179  	Random(iub, 3**n)
   180  	out.WriteString(">THREE Homo sapiens frequency\n")
   181  	Random(homosapiens, 5**n)
   182  }
   183  
   184  type buffer []byte
   185  
   186  func (b *buffer) Flush() {
   187  	p := *b
   188  	if len(p) > 0 {
   189  		os.Stdout.Write(p)
   190  	}
   191  	*b = p[0:0]
   192  }
   193  
   194  func (b *buffer) WriteString(s string) {
   195  	p := b.NextWrite(len(s))
   196  	copy(p, s)
   197  }
   198  
   199  func (b *buffer) NextWrite(n int) []byte {
   200  	p := *b
   201  	if len(p)+n > cap(p) {
   202  		b.Flush()
   203  		p = *b
   204  	}
   205  	out := p[len(p) : len(p)+n]
   206  	*b = p[:len(p)+n]
   207  	return out
   208  }
   209
View as plain text