9 const int RSPD_DEFAULT_B = 20;
13 RSPD(bool estRSPD, int B = RSPD_DEFAULT_B) {
14 this->estRSPD = estRSPD;
17 pdf = new double[B + 2];
18 cdf = new double[B + 2];
20 //set initial parameters
21 memset(pdf, 0, sizeof(double) * (B + 2)); // use B + 2 for evalCDF
22 memset(cdf, 0, sizeof(double) * (B + 2));
23 for (int i = 1; i <= B; i++) {
34 RSPD& operator=(const RSPD& rv);
39 void update(int fpos, int fullLen, double frac) {
42 if (fpos >= fullLen) return; // if out of range, do not use this hit
45 double a = fpos * 1.0 / fullLen;
48 for (i = ((long long)fpos) * B / fullLen + 1; i < (((long long)fpos + 1) * B - 1) / fullLen + 1; i++) {
50 pdf[i] += (b - a) * fullLen * frac;
53 b = (fpos + 1.0) / fullLen;
54 pdf[i] += (b - a) * fullLen * frac;
59 double evalCDF(int fpos, int fullLen) {
60 int i = ((long long)fpos) * B / fullLen;
61 double val = fpos * 1.0 / fullLen * B;
63 return cdf[i] + (val - i) * pdf[i + 1];
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) ;
73 void collect(const RSPD&);
78 void startSimulation(int, Refs*);
79 int simulate(simul*, int, int);
80 void finishSimulation();
84 int B; // number of bins
91 RSPD& RSPD::operator=(const RSPD& rv) {
92 if (this == &rv) return *this;
96 pdf = new double[rv.B + 2];
97 cdf = new double[rv.B + 2];
100 memcpy(pdf, rv.pdf, sizeof(double) * (B + 2));
101 memcpy(cdf, rv.cdf, sizeof(double) * (B + 2));
108 memset(pdf, 0, sizeof(double) * (B + 2));
109 memset(cdf, 0, sizeof(double) * (B + 2));
112 void RSPD::finish() {
117 for (int i = 1; i <= B; i++) {
121 for (int i = 1; i <= B; i++) {
123 cdf[i] = cdf[i - 1] + pdf[i];
127 void RSPD::collect(const RSPD& o) {
129 for (int i = 1; i <= B; i++) {
134 void RSPD::read(FILE *fi) {
135 //release default space first
140 fscanf(fi, "%d", &val);
141 estRSPD = (val != 0);
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];
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++) {
162 cdf[i] = i * 1.0 / B;
167 void RSPD::write(FILE *fo) {
168 fprintf(fo, "%d\n", estRSPD);
170 fprintf(fo, "%d\n", B);
171 for (int i = 1; i < B; i++) {
172 fprintf(fo, "%.10g ", pdf[i]);
174 fprintf(fo, "%.10g\n", pdf[B]);
178 void RSPD::startSimulation(int M, Refs* refs) {
179 if (!estRSPD) return;
181 rspdDists = new double*[M + 1];
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);
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);
196 void RSPD::finishSimulation() {
197 if (!estRSPD) return;
198 for (int i = 1; i <= M; i++) delete[] rspdDists[i];