1 #ifndef NOISEPROFILE_H_
2 #define NOISEPROFILE_H_
17 memset(c, 0, sizeof(c));
18 memset(p, 0, sizeof(p));
21 NoiseProfile& operator=(const NoiseProfile&);
24 void updateC(const std::string&);
25 void update(const std::string&, double frac);
27 void calcInitParams();
29 double getProb(const std::string&);
30 double getLogP() { return logp; }
32 void collect(const NoiseProfile&);
37 void startSimulation();
38 std::string simulate(simul*, int);
39 void finishSimulation();
42 static const int NCODES = 5;
45 double c[NCODES]; // counts in N0;
48 double *pc; // for simulation
51 NoiseProfile& NoiseProfile::operator=(const NoiseProfile& rv) {
52 if (this == &rv) return *this;
54 memcpy(c, rv.c, sizeof(rv.c));
55 memcpy(p, rv.p, sizeof(rv.p));
59 void NoiseProfile::init() {
60 memset(p, 0, sizeof(p));
63 void NoiseProfile::updateC(const std::string& readseq) {
64 int len = readseq.size();
65 for (int i = 0; i < len; i++) {
66 ++c[get_base_id(readseq[i])];
70 void NoiseProfile::update(const std::string& readseq, double frac) {
71 int len = readseq.size();
72 for (int i = 0; i < len; i++) {
73 p[get_base_id(readseq[i])] += frac;
77 void NoiseProfile::finish() {
82 for (int i = 0; i < NCODES; i++) sum += (p[i] + c[i]);
83 if (sum <= 0.0) return;
84 for (int i = 0; i < NCODES; i++) {
85 p[i] = (p[i] + c[i]) / sum;
86 if (c[i] > 0.0) { logp += c[i] * log(p[i]); }
90 void NoiseProfile::calcInitParams() {
95 for (int i = 0; i < NCODES; i++) sum += (1.0 + c[i]);
96 for (int i = 0; i < NCODES; i++) {
97 p[i] = (1.0 + c[i]) / sum;
98 if (c[i] > 0.0) { logp += c[i] * log(p[i]); }
102 double NoiseProfile::getProb(const std::string& readseq) {
104 int len = readseq.size();
106 for (int i = 0; i < len; i++) {
107 prob *= p[get_base_id(readseq[i])];
113 void NoiseProfile::collect(const NoiseProfile& o) {
114 for (int i = 0; i < NCODES; i++)
118 void NoiseProfile::read(FILE *fi) {
121 memset(c, 0, sizeof(c));
122 fscanf(fi, "%d", &tmp_ncodes);
123 assert(tmp_ncodes == NCODES);
124 for (int i = 0; i < NCODES; i++)
125 fscanf(fi, "%lf", &p[i]);
128 void NoiseProfile::write(FILE *fo) {
129 fprintf(fo, "%d\n", NCODES);
130 for (int i = 0; i < NCODES - 1; i++) {
131 fprintf(fo, "%.10g ", p[i]);
133 fprintf(fo, "%.10g\n", p[NCODES - 1]);
136 void NoiseProfile::startSimulation() {
137 pc = new double[NCODES];
139 for (int i = 0; i < NCODES; i++) {
141 if (i > 0) pc[i] += pc[i - 1];
145 std::string NoiseProfile::simulate(simul* sampler, int len) {
146 std::string readseq = "";
148 for (int i = 0; i < len; i++) {
149 readseq.push_back(getCharacter(sampler->sample(pc, NCODES)));
154 void NoiseProfile::finishSimulation() {
158 #endif /* NOISEPROFILE_H_ */