NoiseProfile.h

   1 #ifndef NOISEPROFILE_H_
   2 #define NOISEPROFILE_H_
   3
   4 #include<cmath>
   5 #include<cstdio>
   6 #include<cstring>
   7 #include<string>
   8 #include<cassert>
   9
  10 #include "utils.h"
  11 #include "RefSeq.h"
  12 #include "simul.h"
  13
  14 class NoiseProfile {
  15 public:
  16         NoiseProfile() {
  17                 logp = 0.0;
  18                 memset(c, 0, sizeof(c));
  19                 memset(p, 0, sizeof(p));
  20         }
  21
  22         NoiseProfile& operator=(const NoiseProfile&);
  23
  24         void init();
  25         void updateC(const std::string&);
  26         void update(const std::string&, double frac);
  27         void finish();
  28         void calcInitParams();
  29
  30         double getProb(const std::string&);
  31         double getLogP() { return logp; }
  32
  33         void collect(const NoiseProfile&);
  34
  35         void read(FILE*);
  36         void write(FILE*);
  37
  38         void startSimulation();
  39         std::string simulate(simul*, int);
  40         void finishSimulation();
  41
  42 private:
  43         static const int NCODES = 5;
  44
  45         double logp;
  46         double c[NCODES]; // counts in N0;
  47         double p[NCODES];
  48
  49         double *pc; // for simulation
  50 };
  51
  52 NoiseProfile& NoiseProfile::operator=(const NoiseProfile& rv) {
  53         if (this == &rv) return *this;
  54         logp = rv.logp;
  55         memcpy(c, rv.c, sizeof(rv.c));
  56         memcpy(p, rv.p, sizeof(rv.p));
  57         return *this;
  58 }
  59
  60 void NoiseProfile::init() {
  61         memset(p, 0, sizeof(p));
  62 }
  63
  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])];
  68         }
  69 }
  70
  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;
  75         }
  76 }
  77
  78 void NoiseProfile::finish() {
  79         double sum;
  80
  81         logp = 0.0;
  82         sum = 0.0;
  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]); }
  88         }
  89 }
  90
  91 void NoiseProfile::calcInitParams() {
  92         double sum;
  93
  94         logp = 0.0;
  95         sum = 0.0;
  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]); }
 100         }
 101 }
 102
 103 double NoiseProfile::getProb(const std::string& readseq) {
 104         double prob = 1.0;
 105         int len = readseq.size();
 106
 107         for (int i = 0; i < len; i++) {
 108                 prob *= p[get_base_id(readseq[i])];
 109         }
 110
 111         return prob;
 112 }
 113
 114 void NoiseProfile::collect(const NoiseProfile& o) {
 115         for (int i = 0; i < NCODES; i++)
 116                 p[i] += o.p[i];
 117 }
 118
 119 void NoiseProfile::read(FILE *fi) {
 120         int tmp_ncodes;
 121
 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);
 127 }
 128
 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]);
 133         }
 134         fprintf(fo, "%.10g\n", p[NCODES - 1]);
 135 }
 136
 137 void NoiseProfile::startSimulation() {
 138         pc = new double[NCODES];
 139
 140         for (int i = 0; i < NCODES; i++) {
 141                 pc[i] = p[i];
 142                 if (i > 0) pc[i] += pc[i - 1];
 143         }
 144 }
 145
 146 std::string NoiseProfile::simulate(simul* sampler, int len) {
 147         std::string readseq = "";
 148
 149         for (int i = 0; i < len; i++) {
 150                 readseq.push_back(getCharacter(sampler->sample(pc, NCODES)));
 151         }
 152         return readseq;
 153 }
 154
 155 void NoiseProfile::finishSimulation() {
 156         delete[] pc;
 157 }
 158
 159 #endif /* NOISEPROFILE_H_ */