RSPD.h

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