RSPD.h

   1 #ifndef RSPD_H_
   2 #define RSPD_H_
   3
   4 #include<cstdio>
   5 #include<cstring>
   6 #include<cassert>
   7
   8 #include "utils.h"
   9 #include "RefSeq.h"
  10 #include "Refs.h"
  11 #include "simul.h"
  12
  13 const int RSPD_DEFAULT_B = 20;
  14
  15 class RSPD {
  16 public:
  17         RSPD(bool estRSPD, int B = RSPD_DEFAULT_B) {
  18                 this->estRSPD = estRSPD;
  19                 this->B = B;
  20
  21                 pdf = new double[B + 2];
  22                 cdf = new double[B + 2];
  23
  24                 //set initial parameters
  25                 memset(pdf, 0, sizeof(double) * (B + 2)); // use B + 2 for evalCDF
  26                 memset(cdf, 0, sizeof(double) * (B + 2));
  27                 for (int i = 1; i <= B; i++) {
  28                         pdf[i] = 1.0 / B;
  29                         cdf[i] = i * 1.0 / B;
  30                 }
  31         }
  32
  33         ~RSPD() {
  34                 delete[] pdf;
  35                 delete[] cdf;
  36         }
  37
  38         RSPD& operator=(const RSPD& rv);
  39
  40         void init();
  41
  42         //fpos starts from 0
  43         void update(int fpos, int fullLen, double frac) {
  44                 assert(estRSPD);
  45
  46                 if (fpos >= fullLen) return;  // if out of range, do not use this hit
  47
  48                 int i;
  49                 double a = fpos * 1.0 / fullLen;
  50                 double b;
  51
  52                 for (i = ((long long)fpos) * B / fullLen + 1; i < (((long long)fpos + 1) * B - 1) / fullLen + 1; i++) {
  53                         b = i * 1.0 / B;
  54                         pdf[i] += (b - a) * fullLen * frac;
  55                         a = b;
  56                 }
  57                 b = (fpos + 1.0) / fullLen;
  58                 pdf[i] += (b - a) * fullLen * frac;
  59         }
  60
  61         void finish();
  62
  63         double evalCDF(int fpos, int fullLen) {
  64                 int i = ((long long)fpos) * B / fullLen;
  65                 double val = fpos * 1.0 / fullLen * B;
  66
  67                 return cdf[i] + (val - i) * pdf[i + 1];
  68         }
  69
  70         double getAdjustedProb(int fpos, int effL, int fullLen) {
  71                 assert(fpos >= 0 && fpos < fullLen && effL <= fullLen);
  72                 if (!estRSPD) return 1.0 / effL;
  73                 double denom = evalCDF(effL, fullLen);
  74                 return (denom >= EPSILON ? (evalCDF(fpos + 1, fullLen) - evalCDF(fpos, fullLen)) / denom : 0.0) ;
  75         }
  76
  77         void collect(const RSPD&);
  78
  79         void read(FILE*);
  80         void write(FILE*);
  81
  82         void startSimulation(int, Refs*);
  83         int simulate(simul*, int, int);
  84         void finishSimulation();
  85
  86 private:
  87         bool estRSPD;
  88         int B; // number of bins
  89         double *pdf, *cdf;
  90
  91         int M;
  92         double **rspdDists;
  93 };
  94
  95 RSPD& RSPD::operator=(const RSPD& rv) {
  96         if (this == &rv) return *this;
  97         if (B != rv.B) {
  98                 delete[] pdf;
  99                 delete[] cdf;
 100                 pdf = new double[rv.B + 2];
 101                 cdf = new double[rv.B + 2];
 102         }
 103         B = rv.B;
 104         memcpy(pdf, rv.pdf, sizeof(double) * (B + 2));
 105         memcpy(cdf, rv.cdf, sizeof(double) * (B + 2));
 106
 107         return *this;
 108 }
 109
 110 void RSPD::init() {
 111         assert(estRSPD);
 112         memset(pdf, 0, sizeof(double) * (B + 2));
 113         memset(cdf, 0, sizeof(double) * (B + 2));
 114 }
 115
 116 void RSPD::finish() {
 117         double sum = 0.0;
 118
 119         assert(estRSPD);
 120
 121         for (int i = 1; i <= B; i++) {
 122                 sum += pdf[i];
 123         }
 124
 125         for (int i = 1; i <= B; i++) {
 126                 pdf[i] /= sum;
 127                 cdf[i] = cdf[i - 1] + pdf[i];
 128         }
 129 }
 130
 131 void RSPD::collect(const RSPD& o) {
 132         assert(estRSPD);
 133         for (int i = 1; i <= B; i++) {
 134                 pdf[i] += o.pdf[i];
 135         }
 136 }
 137
 138 void RSPD::read(FILE *fi) {
 139         //release default space first
 140         delete[] pdf;
 141         delete[] cdf;
 142
 143         int val;
 144         assert(fscanf(fi, "%d", &val) == 1);
 145         estRSPD = (val != 0);
 146
 147         if (estRSPD) {
 148                 assert(fscanf(fi, "%d", &B) == 1);
 149                 pdf = new double[B + 2];
 150                 cdf = new double[B + 2];
 151                 memset(pdf, 0, sizeof(double) * (B + 2));
 152                 memset(cdf, 0, sizeof(double) * (B + 2));
 153                 for (int i = 1; i <= B; i++) {
 154                         assert(fscanf(fi, "%lf", &pdf[i]) == 1);
 155                         cdf[i] = cdf[i - 1] + pdf[i];
 156                 }
 157         }
 158         else {
 159                 B = RSPD_DEFAULT_B;
 160                 pdf = new double[B + 2];
 161                 cdf = new double[B + 2];
 162                 memset(pdf, 0, sizeof(double) * (B + 2));
 163                 memset(cdf, 0, sizeof(double) * (B + 2));
 164                 for (int i = 1; i <= B; i++) {
 165                         pdf[i] = 1.0 / B;
 166                         cdf[i] = i * 1.0 / B;
 167                 }
 168         }
 169 }
 170
 171 void RSPD::write(FILE *fo) {
 172         fprintf(fo, "%d\n", estRSPD);
 173         if (estRSPD) {
 174                 fprintf(fo, "%d\n", B);
 175                 for (int i = 1; i < B; i++) {
 176                         fprintf(fo, "%.10g ", pdf[i]);
 177                 }
 178                 fprintf(fo, "%.10g\n", pdf[B]);
 179         }
 180 }
 181
 182 void RSPD::startSimulation(int M, Refs* refs) {
 183         if (!estRSPD) return;
 184         this->M = M;
 185         rspdDists = new double*[M + 1];
 186         rspdDists[0] = NULL;
 187         for (int i = 1; i <= M; i++) {
 188                 int fullLen = refs->getRef(i).getFullLen();
 189                 rspdDists[i] = new double[fullLen];
 190                 memset(rspdDists[i], 0, sizeof(double) * fullLen);
 191                 for (int j = 0; j < fullLen; j++) rspdDists[i][j] = evalCDF(j + 1, fullLen);
 192         }
 193 }
 194
 195 int RSPD::simulate(simul *sampler, int sid, int effL) {
 196         if (estRSPD) return (rspdDists[sid][effL - 1] > 0.0 ? sampler->sample(rspdDists[sid], effL) : -1);
 197         return int(sampler->random() * effL);
 198 }
 199
 200 void RSPD::finishSimulation() {
 201         if (!estRSPD) return;
 202         for (int i = 1; i <= M; i++) delete[] rspdDists[i];
 203         delete[] rspdDists;
 204 }
 205
 206 #endif /* RSPD_H_ */