test/fasta.cpp

/* The Computer Language Benchmarks Game
   http://shootout.alioth.debian.org/
   contributed by Andrew Moon
*/

#include <stdio.h>
#include <stdlib.h>
#include <string.h>

// limit output, so we do not benchmark speed of printing
void puts_limited(char* x) {
  static int left = 550;
  int len = strlen(x);
  if (len <= left) {
    puts(x);
    left -= len;
    return;
  }
  if (left > 0) {
    x[left] = '\0';
    puts(x);
    x[left] = 'z';
    left = 0;
  }
}

struct Random {
  enum { IM = 139968, IA = 3877, IC = 29573 };
  Random() : last(42) {}
  float get(float max = 1.0f) {
    last = (last * IA + IC) % IM;
    return max * last / IM;
  }

protected:
  unsigned int last;
} rng;

struct IUB {
  int c;
  double p;
  unsigned int pi;
};

struct Cumulative {
  enum {
    slots = 512,
  };

  Cumulative(IUB* start) {
    double p = 0;
    for (IUB* iter = start; iter->c; ++iter) {
      p += iter->p;
      iter->p = p < 1.0 ? p : 1.0;
      iter->pi = (unsigned int)(iter->p * slots);
    }

    for (unsigned int i = 0; i <= slots; i++) {
      while (i > start->pi && start->pi != 0) {
        ++start;
      }

      table[i] = start;
    }
  }

  const char operator[](float pct) const {
    IUB* iter = table[(unsigned int)(pct * slots)];
    while (iter->p < pct)
      ++iter;
    return iter->c;
  }

protected:
  IUB* table[slots + 1];
};

static const size_t lineLength = 60;

struct LineBuffer {
  LineBuffer() : lastN(0) {}
  LineBuffer& genrand(Cumulative& table, size_t N) {
    // assert(N <= lineLength);
    for (size_t i = 0; i < N; i++)
      buffer[i] = table[rng.get()];
    buffer[N] = '\n';
    buffer[N + 1] = '\0';
    lastN = N + 1;
    return *this;
  }
  void writeline() { puts_limited(buffer); }

protected:
  char buffer[lineLength + 2];
  size_t lastN;
};

struct RotatingString {
  RotatingString(const char* in) : pos(0) {
    size = strlen(in);
    buffer = new char[size + lineLength];
    memcpy(buffer, in, size);
    memcpy(buffer + size, in, lineLength);
  }
  ~RotatingString() { delete[] buffer; }
  void write(size_t bytes) {
    char* temp = new char[bytes + 2];
    memcpy(temp, buffer + pos, bytes);
    temp[bytes] = '\n';
    temp[bytes] = '\0';
    puts_limited(temp);
    delete temp;
    pos += bytes;
    if (pos > size)
      pos -= size;
  }

protected:
  char* buffer;
  size_t size, pos;
};

template<class Output>
void makeFasta(const char* id, const char* desc, size_t N, Output& output) {
  while (N) {
    const size_t bytes = N < lineLength ? N : lineLength;
    output.writeline(bytes);
    N -= bytes;
  }
}

struct Repeater {
  Repeater(const char* alu) : rot(alu) {}
  void writeline(size_t bytes) { rot.write(bytes); }
  void run(const char* id, const char* desc, size_t N) {
    makeFasta(id, desc, N, *this);
  }

protected:
  RotatingString rot;
};

struct Randomized {
  Randomized(IUB* start) : table(start) {}
  void writeline(size_t bytes) { line.genrand(table, bytes).writeline(); }
  void run(const char* id, const char* desc, size_t N) {
    makeFasta(id, desc, N, *this);
  }

protected:
  Cumulative table;
  LineBuffer line;
};

IUB iub[] = {
  {'a', 0.27, 0},
  {'c', 0.12, 0},
  {'g', 0.12, 0},
  {'t', 0.27, 0},

  {'B', 0.02, 0},
  {'D', 0.02, 0},
  {'H', 0.02, 0},
  {'K', 0.02, 0},
  {'M', 0.02, 0},
  {'N', 0.02, 0},
  {'R', 0.02, 0},
  {'S', 0.02, 0},
  {'V', 0.02, 0},
  {'W', 0.02, 0},
  {'Y', 0.02, 0},
  {0, 0, 0},
};

IUB homosapiens[] = {
  {'a', 0.3029549426680, 0},
  {'c', 0.1979883004921, 0},
  {'g', 0.1975473066391, 0},
  {'t', 0.3015094502008, 0},
  {0, 0, 0},
};

static const char alu[] = "GGCCGGGCGCGGTGGCTCACGCCTGTAATCCCAGCACTTTG"
                          "GGAGGCCGAGGCGGGCGGATCACCTGAGGTCAGGAGTTCGA"
                          "GACCAGCCTGGCCAACATGGTGAAACCCCGTCTCTACTAAA"
                          "AATACAAAAATTAGCCGGGCGTGGTGGCGCGCGCCTGTAAT"
                          "CCCAGCTACTCGGGAGGCTGAGGCAGGAGAATCGCTTGAAC"
                          "CCGGGAGGCGGAGGTTGCAGTGAGCCGAGATCGCGCCACTG"
                          "CACTCCAGCCTGGGCGACAGAGCGAGACTCCGTCTCAAAAA";

int main(int argc, const char* argv[]) {
  const size_t n = (argc > 1) ? atoi(argv[1]) : 512;

  Repeater(alu).run("ONE", "Homo sapiens alu", n * 2);
  Randomized(iub).run("TWO", "IUB ambiguity codes", n * 3);
  Randomized(homosapiens).run("THREE", "Homo sapiens frequency", n * 5);

  return 0;
}