- int local_N;
- int fr, to, len, id;
- vector<double> arr;
- arr.clear();
-
- if (verbose) printf("Begin to sample reads from their posteriors.\n");
- for (int i = 0; i < nThreads; i++) {
- local_N = hitvs[i]->getN();
- for (int j = 0; j < local_N; j++) {
- fr = hitvs[i]->getSAt(j);
- to = hitvs[i]->getSAt(j + 1);
- len = to - fr + 1;
- arr.resize(len);
- arr[0] = ncpvs[i][j];
- for (int k = fr; k < to; k++) arr[k - fr + 1] = arr[k - fr] + hitvs[i]->getHitAt(k).getConPrb();
- id = (arr[len - 1] < EPSILON ? -1 : sample(arr, len)); // if all entries in arr are 0, let id be -1
- for (int k = fr; k < to; k++) hitvs[i]->getHitAt(k).setConPrb(k - fr + 1 == id ? 1.0 : 0.0);
- }
- }
+ READ_INT_TYPE local_N;
+ HIT_INT_TYPE fr, to, len, id;
+ vector<double> arr;
+ uniform01 rg(engine_type(time(NULL)));
+
+ if (verbose) cout<< "Begin to sample reads from their posteriors."<< endl;
+ for (int i = 0; i < nThreads; i++) {
+ local_N = hitvs[i]->getN();
+ for (READ_INT_TYPE j = 0; j < local_N; j++) {
+ fr = hitvs[i]->getSAt(j);
+ to = hitvs[i]->getSAt(j + 1);
+ len = to - fr + 1;
+ arr.assign(len, 0);
+ arr[0] = ncpvs[i][j];
+ for (HIT_INT_TYPE k = fr; k < to; k++) arr[k - fr + 1] = arr[k - fr] + hitvs[i]->getHitAt(k).getConPrb();
+ id = (arr[len - 1] < EPSILON ? -1 : sample(rg, arr, len)); // if all entries in arr are 0, let id be -1
+ for (HIT_INT_TYPE k = fr; k < to; k++) hitvs[i]->getHitAt(k).setConPrb(k - fr + 1 == id ? 1.0 : 0.0);
+ }
+ }
+
+ if (verbose) cout<< "Sampling is finished."<< endl;