diff --git a/cmd_cram_demuxlet.cpp b/cmd_cram_demuxlet.cpp index d8ebe75..c4379c7 100644 --- a/cmd_cram_demuxlet.cpp +++ b/cmd_cram_demuxlet.cpp @@ -535,30 +535,6 @@ int32_t main(int32_t argc, char** argv) { if ( ( writePair && wpair == NULL ) || ( wsingle == NULL ) ) error("[E:%s:%d %s] Cannot create %s.single, %s.pair files",__FILE__,__LINE__,__FUNCTION__,outPrefix.c_str(), outPrefix.c_str()); - - // start finding the next-best matching individual - // here we iterate each cell separately. - // pre-calculate nsnp*nv*nv*9, nv*1*9, 1*nv*9, 1*9 - double* gpAB = new double[scl.nsnps * nv * nv * 9]; - double* gpA0 = new double[scl.nsnps * nv * 9]; - double* gp00 = new double[scl.nsnps * 9]; - - int32_t j, k, l, m, n; - for(i=0; i < scl.nsnps; ++i) { - for(j=0; j < nv; ++j) { - for(k=0; k < nv; ++k) { - gps = scl.snps[i].gps; - for(l=0; l < 3; ++l) { - for(m=0; m < 3; ++m) { - gpAB[i*nv*nv*9 + j*nv*9 + k*9 + l*3 + m] = gps[j*3+l] * gps[k*3+m]; - gpA0[i*nv*9 + j*9 + l*3 + m] = gps[j*3+l] * gp0s[i*3+m]; - gp00[i*9 + l*3 + m] = gp0s[i*3+l] * gp0s[i*3+m]; - } - } - } - } - } - // iterate each barcode int32_t n1 = nv; double* llksAB = new double[n1 * nv * nAlpha]; @@ -668,13 +644,14 @@ int32_t main(int32_t argc, char** argv) { double p; int32_t j, k, l, m, n; - for(j=jbeg; j < jend; ++j) { + gps = scl.snps[isnp].gps; + for(j=jbeg; j < jend; ++j) { // pairwise LLK for(k=0; k < nv; ++k) { std::fill(sumPs.begin(), sumPs.end(), 0); for(l=0; l < 3; ++l) { for(m=0; m < 3; ++m) { - p = gpAB[isnp*nv*nv*9 + j*nv*9 + k*9 + l*3 + m]; + p = gps[j*3+l] * gps[k*3+m]; // Compute this on the go for avoiding explosion of the RAM for(n=0; n < nAlpha; ++n) sumPs[n] += (p * pGs[n*9+l*3+m]); } @@ -687,7 +664,7 @@ int32_t main(int32_t argc, char** argv) { std::fill(sumPs.begin(), sumPs.end(), 0); for(l=0; l < 3; ++l) { for(m=0; m < 3; ++m) { - p = gpA0[isnp*nv*9 + j*9 + l*3 + m]; + p = gps[j*3+l] * gp0s[isnp*3+m]; // Compute this on the go for avoiding explosion of the RAM for(n=0; n < nAlpha; ++n) sumPs[n] += (p * pGs[n*9+l*3+m]); } @@ -700,7 +677,7 @@ int32_t main(int32_t argc, char** argv) { std::fill(sumPs.begin(), sumPs.end(), 0); for(l=0; l < 3; ++l) { for(m=0; m < 3; ++m) { - p = gp00[isnp*9 + l*3 + m]; + p = gp0s[isnp*3+l] * gp0s[isnp*3+m]; // Compute this on the go for avoiding explosion of the RAM for(n=0; n < nAlpha; ++n) sumPs[n] += (p * pGs[n*9+l*3+m]); } @@ -709,6 +686,7 @@ int32_t main(int32_t argc, char** argv) { llks00[n] += log(sumPs[n]); } + int32_t j, k, n; // normalize by max likelihood double maxLLK = -1e300; for(j=jbeg; j < jend; ++j) {