QualDist.h

   1 #ifndef QUALDIST_H_
   2 #define QUALDIST_H_
   3
   4 #include<cstdio>
   5 #include<cstring>
   6 #include<cassert>
   7 #include<string>
   8
   9 #include "simul.h"
  10
  11 //from 33 to 126 to encode 0 to 93
  12 class QualDist {
  13 public:
  14         QualDist() {
  15                 memset(p_init, 0, sizeof(p_init));
  16                 memset(p_tran, 0, sizeof(p_tran));
  17         }
  18
  19         QualDist& operator=(const QualDist&);
  20
  21         void update(const std::string&);
  22         void finish();
  23
  24         double getProb(const std::string&);
  25
  26         void read(FILE*);
  27         void write(FILE*);
  28
  29         void startSimulation();
  30         std::string simulate(simul*, int);
  31         void finishSimulation();
  32
  33 private:
  34         static const int SIZE = 100;
  35
  36         double p_init[SIZE];
  37         double p_tran[SIZE][SIZE]; //p_tran[a][b] = p(b|a)
  38
  39         int c2q(char c) { assert(c >= 33 && c <= 126); return c - 33; }
  40
  41         double *qc_init, (*qc_trans)[SIZE];
  42         char q2c(int qval) { return (char)(qval + 33); }
  43 };
  44
  45 QualDist& QualDist::operator=(const QualDist& rv) {
  46         if (this == &rv) return *this;
  47
  48         memcpy(p_init, rv.p_init, sizeof(rv.p_init));
  49         memcpy(p_tran, rv.p_tran, sizeof(rv.p_tran));
  50
  51         return *this;
  52 }
  53
  54 void QualDist::update(const std::string& qual) {
  55         int len = qual.size();
  56
  57         assert(len > 0);
  58         ++p_init[c2q(qual[0])];
  59
  60         for (int i = 1; i < len; i++) {
  61                 ++p_tran[c2q(qual[i - 1])][c2q(qual[i])];
  62         }
  63 }
  64
  65 void QualDist::finish() {
  66         double sum;
  67
  68         sum = 0.0;
  69         for (int i = 0; i < SIZE; i++) sum += p_init[i];
  70         for (int i = 0; i < SIZE; i++) p_init[i] /= sum;
  71
  72         for (int i = 0; i < SIZE; i++) {
  73                 sum = 0.0;
  74                 for (int j = 0; j < SIZE; j++) sum += p_tran[i][j];
  75                 if (sum <= 0.0) continue;
  76                 //if (isZero(sum)) continue;
  77                 for (int j = 0; j < SIZE; j++) p_tran[i][j] /= sum;
  78         }
  79 }
  80
  81 double QualDist::getProb(const std::string& qual) {
  82         int len = qual.size();
  83         double prob = 1.0;
  84
  85         assert(len > 0);
  86         prob *= p_init[c2q(qual[0])];
  87         for (int i = 1; i < len; i++) {
  88                 prob *= p_tran[c2q(qual[i - 1])][c2q(qual[i])];
  89         }
  90
  91         return prob;
  92 }
  93
  94 void QualDist::read(FILE *fi) {
  95         int tmp_size;
  96
  97         assert(fscanf(fi, "%d", &tmp_size) == 1);
  98         assert(tmp_size == SIZE);
  99
 100         for (int i = 0; i < SIZE; i++) { assert(fscanf(fi, "%lf", &p_init[i]) == 1); }
 101         for (int i = 0; i < SIZE; i++) {
 102                 for (int j = 0; j < SIZE; j++) { assert(fscanf(fi, "%lf", &p_tran[i][j]) == 1); }
 103         }
 104 }
 105
 106 void QualDist::write(FILE *fo) {
 107         fprintf(fo, "%d\n", SIZE);
 108         for (int i = 0; i < SIZE - 1; i++) { fprintf(fo, "%.10g ", p_init[i]); }
 109         fprintf(fo, "%.10g\n", p_init[SIZE - 1]);
 110         for (int i = 0; i < SIZE; i++) {
 111                 for (int j = 0; j < SIZE -1 ; j++) fprintf(fo, "%.10g ", p_tran[i][j]);
 112                 fprintf(fo, "%.10g\n", p_tran[i][SIZE - 1]);
 113         }
 114 }
 115
 116 void QualDist::startSimulation() {
 117         qc_init = new double[SIZE];
 118         qc_trans = new double[SIZE][SIZE];
 119
 120         for (int i = 0; i < SIZE; i++) {
 121                 qc_init[i] = p_init[i];
 122                 if (i > 0) qc_init[i] += qc_init[i - 1];
 123         }
 124
 125         for (int i = 0; i < SIZE; i++)
 126                 for (int j = 0; j < SIZE; j++) {
 127                         qc_trans[i][j] = p_tran[i][j];
 128                         if (j > 0) qc_trans[i][j] += qc_trans[i][j - 1];
 129                 }
 130 }
 131
 132 std::string QualDist::simulate(simul* sampler, int len) {
 133         int qval, old_qval;
 134         std::string qual = "";
 135
 136         qval = sampler->sample(qc_init, SIZE);
 137         qual.push_back(q2c(qval));
 138         for (int i = 1; i < len; i++) {
 139                 old_qval = qval;
 140                 qval = sampler->sample(qc_trans[old_qval], SIZE);
 141                 qual.push_back(q2c(qval));
 142         }
 143
 144         return qual;
 145 }
 146
 147 void QualDist::finishSimulation() {
 148         delete[] qc_init;
 149         delete[] qc_trans;
 150 }
 151 #endif /* QUALDIST_H_ */