13 const int RSPD_DEFAULT_B = 20;
17 RSPD(bool estRSPD, int B = RSPD_DEFAULT_B) {
18 this->estRSPD = estRSPD;
21 pdf = new double[B + 2];
22 cdf = new double[B + 2];
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++) {
38 RSPD& operator=(const RSPD& rv);
43 void update(int fpos, int fullLen, double frac) {
46 if (fpos >= fullLen) return; // if out of range, do not use this hit
49 double a = fpos * 1.0 / fullLen;
52 for (i = ((long long)fpos) * B / fullLen + 1; i < (((long long)fpos + 1) * B - 1) / fullLen + 1; i++) {
54 pdf[i] += (b - a) * fullLen * frac;
57 b = (fpos + 1.0) / fullLen;
58 pdf[i] += (b - a) * fullLen * frac;
63 double evalCDF(int fpos, int fullLen) {
64 int i = ((long long)fpos) * B / fullLen;
65 double val = fpos * 1.0 / fullLen * B;
67 return cdf[i] + (val - i) * pdf[i + 1];
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) ;
77 void collect(const RSPD&);
82 void startSimulation(int, Refs*);
83 int simulate(simul*, int, int);
84 void finishSimulation();
88 int B; // number of bins
95 RSPD& RSPD::operator=(const RSPD& rv) {
96 if (this == &rv) return *this;
100 pdf = new double[rv.B + 2];
101 cdf = new double[rv.B + 2];
104 memcpy(pdf, rv.pdf, sizeof(double) * (B + 2));
105 memcpy(cdf, rv.cdf, sizeof(double) * (B + 2));
112 memset(pdf, 0, sizeof(double) * (B + 2));
113 memset(cdf, 0, sizeof(double) * (B + 2));
116 void RSPD::finish() {
121 for (int i = 1; i <= B; i++) {
125 for (int i = 1; i <= B; i++) {
127 cdf[i] = cdf[i - 1] + pdf[i];
131 void RSPD::collect(const RSPD& o) {
133 for (int i = 1; i <= B; i++) {
138 void RSPD::read(FILE *fi) {
139 //release default space first
144 assert(fscanf(fi, "%d", &val) == 1);
145 estRSPD = (val != 0);
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];
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++) {
166 cdf[i] = i * 1.0 / B;
171 void RSPD::write(FILE *fo) {
172 fprintf(fo, "%d\n", estRSPD);
174 fprintf(fo, "%d\n", B);
175 for (int i = 1; i < B; i++) {
176 fprintf(fo, "%.10g ", pdf[i]);
178 fprintf(fo, "%.10g\n", pdf[B]);
182 void RSPD::startSimulation(int M, Refs* refs) {
183 if (!estRSPD) return;
185 rspdDists = new double*[M + 1];
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);
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);
200 void RSPD::finishSimulation() {
201 if (!estRSPD) return;
202 for (int i = 1; i <= M; i++) delete[] rspdDists[i];