nsz Git - libfirm/blob - ir/be/test/langshootout/fasta.c

   1 /* -*- mode: c -*-
   2  * $Id$
   3  * http://shootout.alioth.debian.org/
   4  *
   5  * by Paul Hsieh
   6  */
   7
   8 #include <stdio.h>
   9 #include <stdlib.h>
  10
  11 #define IM 139968
  12 #define IA   3877
  13 #define IC  29573
  14
  15 double gen_random (double max) {
  16     static long last = 42;
  17     return max * (last = (last * IA + IC) % IM) / IM;
  18 }
  19
  20 struct aminoacids {
  21     char c;
  22     double p;
  23 };
  24
  25 /* Weighted selection from alphabet */
  26
  27 void makeCumulative (struct aminoacids * genelist, int count) {
  28     double cp = 0.0;
  29     int i;
  30
  31     for (i=0; i < count; i++) {
  32         cp += genelist[i].p;
  33         genelist[i].p = cp;
  34     }
  35 }
  36
  37 char selectRandom (const struct aminoacids * genelist, int count) {
  38     double r = gen_random (1);
  39     int i, lo, hi;
  40
  41     if (r < genelist[0].p) return genelist[0].c;
  42
  43     lo = 0;
  44     hi = count-1;
  45
  46     while (hi > lo+1) {
  47         i = (hi + lo) / 2;
  48         if (r < genelist[i].p) hi = i; else lo = i;
  49     }
  50     return genelist[hi].c;
  51 }
  52
  53 /* Generate and write FASTA format */
  54
  55 #define LINE_LENGTH (60)
  56
  57 void makeRandomFasta (const char * id, const char * desc, const struct
  58 aminoacids * genelist, int count, int n) {
  59    int todo = n;
  60    int i, m;
  61
  62    printf (">%s %s\n", id, desc);
  63
  64    for (; todo > 0; todo -= LINE_LENGTH) {
  65        char pick[LINE_LENGTH+1];
  66        if (todo < LINE_LENGTH) m = todo; else m = LINE_LENGTH;
  67        for (i=0; i < m; i++) pick[i] = selectRandom (genelist, count);
  68        pick[m] = '\0';
  69        puts (pick);
  70    }
  71 }
  72
  73 void makeRepeatFasta (const char * id, const char * desc, const char *
  74 s, int n) {
  75    char * ss;
  76    int todo = n, k = 0, kn = strlen (s);
  77    int m;
  78
  79    ss = (char *) malloc (kn + 1);
  80    memcpy (ss, s, kn+1);
  81
  82    printf (">%s %s\n", id, desc);
  83
  84    for (; todo > 0; todo -= LINE_LENGTH) {
  85        if (todo < LINE_LENGTH) m = todo; else m = LINE_LENGTH;
  86
  87        while (m >= kn - k) {
  88            printf ("%s", s+k);
  89            m -= kn - k;
  90            k = 0;
  91        }
  92
  93        ss[k + m] = '\0';
  94        puts (ss+k);
  95        ss[k + m] = s[m+k];
  96        k += m;
  97    }
  98
  99    free (ss);
 100 }
 101
 102 /* Main -- define alphabets, make 3 fragments */
 103
 104 struct aminoacids iub[] = {
 105     { 'a', 0.27 },
 106     { 'c', 0.12 },
 107     { 'g', 0.12 },
 108     { 't', 0.27 },
 109
 110     { 'B', 0.02 },
 111     { 'D', 0.02 },
 112     { 'H', 0.02 },
 113     { 'K', 0.02 },
 114     { 'M', 0.02 },
 115     { 'N', 0.02 },
 116     { 'R', 0.02 },
 117     { 'S', 0.02 },
 118     { 'V', 0.02 },
 119     { 'W', 0.02 },
 120     { 'Y', 0.02 }
 121 };
 122
 123 #define IUB_LEN (sizeof (iub) / sizeof (struct aminoacids))
 124
 125 struct aminoacids homosapiens[] = {
 126     { 'a', 0.3029549426680 },
 127     { 'c', 0.1979883004921 },
 128     { 'g', 0.1975473066391 },
 129     { 't', 0.3015094502008 },
 130 };
 131
 132 #define HOMOSAPIENS_LEN (sizeof (homosapiens) / sizeof (struct aminoacids))
 133
 134 char * alu =
 135    "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTGG"
 136    "GAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGAGA"
 137    "CCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAAAAT"
 138    "ACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAATCCCA"
 139    "GCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAACCCGGG"
 140    "AGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTGCACTCC"
 141    "AGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";
 142
 143 int main (int argc, char * argv[]) {
 144     int n = 1000;
 145
 146     if (argc > 1) sscanf (argv[1], "%d", &n);
 147
 148     makeCumulative (iub, IUB_LEN);
 149     makeCumulative (homosapiens, HOMOSAPIENS_LEN);
 150
 151     makeRepeatFasta ("ONE", "Homo sapiens alu", alu, n*2);
 152     makeRandomFasta ("TWO", "IUB ambiguity codes", iub, IUB_LEN, n*3);
 153     makeRandomFasta ("THREE", "Homo sapiens frequency", homosapiens,
 154 HOMOSAPIENS_LEN, n*5);
 155
 156     return 0;
 157 }