11 //from 33 to 126 to encode 0 to 93
15 memset(p_init, 0, sizeof(p_init));
16 memset(p_tran, 0, sizeof(p_tran));
19 QualDist& operator=(const QualDist&);
21 void update(const std::string&);
24 double getProb(const std::string&);
29 void startSimulation();
30 std::string simulate(simul*, int);
31 void finishSimulation();
34 static const int SIZE = 100;
37 double p_tran[SIZE][SIZE]; //p_tran[a][b] = p(b|a)
39 int c2q(char c) { assert(c >= 33 && c <= 126); return c - 33; }
41 double *qc_init, (*qc_trans)[SIZE];
42 char q2c(int qval) { return (char)(qval + 33); }
45 QualDist& QualDist::operator=(const QualDist& rv) {
46 if (this == &rv) return *this;
48 memcpy(p_init, rv.p_init, sizeof(rv.p_init));
49 memcpy(p_tran, rv.p_tran, sizeof(rv.p_tran));
54 void QualDist::update(const std::string& qual) {
55 int len = qual.size();
58 ++p_init[c2q(qual[0])];
60 for (int i = 1; i < len; i++) {
61 ++p_tran[c2q(qual[i - 1])][c2q(qual[i])];
65 void QualDist::finish() {
69 for (int i = 0; i < SIZE; i++) sum += p_init[i];
70 for (int i = 0; i < SIZE; i++) p_init[i] /= sum;
72 for (int i = 0; i < SIZE; i++) {
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;
81 double QualDist::getProb(const std::string& qual) {
82 int len = qual.size();
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])];
94 void QualDist::read(FILE *fi) {
97 fscanf(fi, "%d", &tmp_size);
98 assert(tmp_size == SIZE);
100 for (int i = 0; i < SIZE; i++) { fscanf(fi, "%lf", &p_init[i]); }
101 for (int i = 0; i < SIZE; i++) {
102 for (int j = 0; j < SIZE; j++) { fscanf(fi, "%lf", &p_tran[i][j]); }
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]);
116 void QualDist::startSimulation() {
117 qc_init = new double[SIZE];
118 qc_trans = new double[SIZE][SIZE];
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];
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];
132 std::string QualDist::simulate(simul* sampler, int len) {
134 std::string qual = "";
136 qval = sampler->sample(qc_init, SIZE);
137 qual.push_back(q2c(qval));
138 for (int i = 1; i < len; i++) {
140 qval = sampler->sample(qc_trans[old_qval], SIZE);
141 qual.push_back(q2c(qval));
147 void QualDist::finishSimulation() {
151 #endif /* QUALDIST_H_ */