NoiseProfile.h

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