1 #ifndef NOISEPROFILE_H_
2 #define NOISEPROFILE_H_
18 memset(c, 0, sizeof(c));
19 memset(p, 0, sizeof(p));
22 NoiseProfile& operator=(const NoiseProfile&);
25 void updateC(const std::string&);
26 void update(const std::string&, double frac);
28 void calcInitParams();
30 double getProb(const std::string&);
31 double getLogP() { return logp; }
33 void collect(const NoiseProfile&);
38 void startSimulation();
39 std::string simulate(simul*, int);
40 void finishSimulation();
43 static const int NCODES = 5;
46 double c[NCODES]; // counts in N0;
49 double *pc; // for simulation
52 NoiseProfile& NoiseProfile::operator=(const NoiseProfile& rv) {
53 if (this == &rv) return *this;
55 memcpy(c, rv.c, sizeof(rv.c));
56 memcpy(p, rv.p, sizeof(rv.p));
60 void NoiseProfile::init() {
61 memset(p, 0, sizeof(p));
64 void NoiseProfile::updateC(const std::string& readseq) {
65 int len = readseq.size();
66 for (int i = 0; i < len; i++) {
67 ++c[get_base_id(readseq[i])];
71 void NoiseProfile::update(const std::string& readseq, double frac) {
72 int len = readseq.size();
73 for (int i = 0; i < len; i++) {
74 p[get_base_id(readseq[i])] += frac;
78 void NoiseProfile::finish() {
83 for (int i = 0; i < NCODES; i++) sum += (p[i] + c[i]);
84 if (sum <= EPSILON) return;
85 for (int i = 0; i < NCODES; i++) {
86 p[i] = (p[i] + c[i]) / sum;
87 if (c[i] > 0.0) { logp += c[i] * log(p[i]); }
91 void NoiseProfile::calcInitParams() {
96 for (int i = 0; i < NCODES; i++) sum += (1.0 + c[i]);
97 for (int i = 0; i < NCODES; i++) {
98 p[i] = (1.0 + c[i]) / sum;
99 if (c[i] > 0.0) { logp += c[i] * log(p[i]); }
103 double NoiseProfile::getProb(const std::string& readseq) {
105 int len = readseq.size();
107 for (int i = 0; i < len; i++) {
108 prob *= p[get_base_id(readseq[i])];
114 void NoiseProfile::collect(const NoiseProfile& o) {
115 for (int i = 0; i < NCODES; i++)
119 void NoiseProfile::read(FILE *fi) {
122 memset(c, 0, sizeof(c));
123 assert(fscanf(fi, "%d", &tmp_ncodes) == 1);
124 assert(tmp_ncodes == NCODES);
125 for (int i = 0; i < NCODES; i++)
126 assert(fscanf(fi, "%lf", &p[i]) == 1);
129 void NoiseProfile::write(FILE *fo) {
130 fprintf(fo, "%d\n", NCODES);
131 for (int i = 0; i < NCODES - 1; i++) {
132 fprintf(fo, "%.10g ", p[i]);
134 fprintf(fo, "%.10g\n", p[NCODES - 1]);
137 void NoiseProfile::startSimulation() {
138 pc = new double[NCODES];
140 for (int i = 0; i < NCODES; i++) {
142 if (i > 0) pc[i] += pc[i - 1];
146 std::string NoiseProfile::simulate(simul* sampler, int len) {
147 std::string readseq = "";
149 for (int i = 0; i < len; i++) {
150 readseq.push_back(getCharacter(sampler->sample(pc, NCODES)));
155 void NoiseProfile::finishSimulation() {
159 #endif /* NOISEPROFILE_H_ */