Tiny change to last-split probabilities
authorMartin C. Frith
Wed Mar 11 16:58:27 2020 +0900 (2 months ago)
changeset 1051eb537293aa74
parent 1050 cf3b0bf3dac8
child 1052 ba6739cb6e7b
Tiny change to last-split probabilities
src/split/cbrc_split_aligner.cc
test/last-split-test.out
     1.1 --- a/src/split/cbrc_split_aligner.cc	Wed Mar 11 16:35:18 2020 +0900
     1.2 +++ b/src/split/cbrc_split_aligner.cc	Wed Mar 11 16:58:27 2020 +0900
     1.3 @@ -504,9 +504,9 @@
     1.4  	updateInplayAlnIndicesF(sortedAlnPos, oldNumInplay, newNumInplay, j);
     1.5  	unsigned oldInplayPos = 0;
     1.6  	cell(rescales, j) = rescale;
     1.7 -	zF /= rescale;
     1.8 +	zF *= rescale;
     1.9  	double pSum = 0.0;
    1.10 -	double rNew = 1.0;
    1.11 +	double rNew = 0.0;
    1.12  	for (unsigned x = 0; x < newNumInplay; ++x) {
    1.13  	    unsigned i = newInplayAlnIndices[x];
    1.14  	    size_t ij = matrixRowOrigins[i] + j;
    1.15 @@ -517,22 +517,22 @@
    1.16  	    p *= spliceEndProb(ij);
    1.17  	    p += Fmat[ij] * Dexp[ij];
    1.18  	    if (restartProb <= 0 && alns[i].qstart == j) p += begprob;
    1.19 -	    p = p * Aexp[ij] / rescale;
    1.20 +	    p = p * Aexp[ij] * rescale;
    1.21  
    1.22  	    Fmat[ij + 1] = p;
    1.23  	    if (alns[i].qend == j+1) zF += p;
    1.24  	    pSum += p * spliceBegProb(ij + 1);
    1.25  	    rNew += p;
    1.26          }
    1.27 -        begprob /= rescale;
    1.28 -	sumProb = pSum + sumProb / rescale;
    1.29 +        begprob *= rescale;
    1.30 +	sumProb = pSum + sumProb * rescale;
    1.31  	probFromJump = pSum * jumpProb + sumProb * restartProb;
    1.32 -	rescale = rNew;
    1.33 +	rescale = 1 / (rNew + 1);
    1.34      }
    1.35  
    1.36      if (restartProb > 0) zF = sumProb;
    1.37 -    //zF /= cell(rescales, maxEnd);
    1.38 -    cell(rescales, maxEnd) = zF;  // this causes scaled zF to equal 1
    1.39 +    //zF *= cell(rescales, maxEnd);
    1.40 +    cell(rescales, maxEnd) = 1 / zF;  // this causes scaled zF to equal 1
    1.41  }
    1.42  
    1.43  void SplitAligner::backward() {
    1.44 @@ -555,7 +555,7 @@
    1.45  	updateInplayAlnIndicesB(sortedAlnPos, oldNumInplay, newNumInplay, j);
    1.46  	unsigned oldInplayPos = 0;
    1.47  	double rescale = cell(rescales, j);
    1.48 -	//zB /= rescale;
    1.49 +	//zB *= rescale;
    1.50  	double pSum = 0.0;
    1.51  	for (unsigned x = 0; x < newNumInplay; ++x) {
    1.52  	    unsigned i = newInplayAlnIndices[x];
    1.53 @@ -567,19 +567,19 @@
    1.54  	    p *= spliceBegProb(ij);
    1.55  	    p += Bmat[ij] * Dexp[ij];
    1.56  	    if (restartProb <= 0 && alns[i].qend == j) p += endprob;
    1.57 -	    p = p * Aexp[ij - 1] / rescale;
    1.58 +	    p = p * Aexp[ij - 1] * rescale;
    1.59  
    1.60  	    Bmat[ij - 1] = p;
    1.61  	    //if (alns[i].qstart == j-1) zB += p;
    1.62  	    pSum += p * spliceEndProb(ij - 1);
    1.63          }
    1.64 -        endprob /= rescale;
    1.65 -	sumProb = pSum * restartProb + sumProb / rescale;
    1.66 +        endprob *= rescale;
    1.67 +	sumProb = pSum * restartProb + sumProb * rescale;
    1.68  	probFromJump = pSum * jumpProb + sumProb;
    1.69      }
    1.70  
    1.71      //if (restartProb > 0) zB = sumProb;
    1.72 -    //zB /= cell(rescales, minBeg);
    1.73 +    //zB *= cell(rescales, minBeg);
    1.74  }
    1.75  
    1.76  std::vector<double>
    1.77 @@ -591,7 +591,7 @@
    1.78      for (unsigned pos = alnBeg; pos < alnEnd; ++pos) {
    1.79  	size_t ij = matrixRowOrigins[i] + j;
    1.80          if (alns[i].qalign[pos] == '-') {
    1.81 -            double value = Fmat[ij] * Bmat[ij] * Dexp[ij] / cell(rescales, j);
    1.82 +            double value = Fmat[ij] * Bmat[ij] * Dexp[ij] * cell(rescales, j);
    1.83              output.push_back(value);
    1.84          } else {
    1.85              double value = Fmat[ij + 1] * Bmat[ij] / Aexp[ij];
    1.86 @@ -992,7 +992,7 @@
    1.87    assert(rescales.size() == rescalesRev.size());
    1.88    double logOdds = 0;
    1.89    for (unsigned j = 0; j < rescales.size(); ++j) {
    1.90 -    logOdds += std::log(rescales[j] / rescalesRev[j]);
    1.91 +    logOdds += std::log(rescalesRev[j] / rescales[j]);
    1.92    }
    1.93    return logOdds;
    1.94  }
     2.1 --- a/test/last-split-test.out	Wed Mar 11 16:35:18 2020 +0900
     2.2 +++ b/test/last-split-test.out	Wed Mar 11 16:58:27 2020 +0900
     2.3 @@ -46848,7 +46848,7 @@
     2.4  a score=952 mismap=1e-10
     2.5  s chr19 1390869 183 + 58617616 GAGTTCTCTGTGGCCCATGACCTTCGGCCTGGCCTGCTGCGCCGTGGAGATGATGCACATGGCAGCACCCCGCTACGACATGGACCGCTTTGGCGTGGTCTTCCGCGCCAGCccgcgccagtccgACGTCATGATCGTGGCCGGCACACTCA---CCAACAAGATGGCCCCAGCGCTTCGCAAGGT
     2.6  s 102       297 181 -      699 GAGTTCTCTGTGGCCCATGACCTTCGGCCTGGCCTGCTGCGCCGTGGAGATGATGCGCATGGCAGCACCCCGCTACGACATGGACCGC--CGGC--GGTCTTCCGCGCCAGCccgcgccagtccgACGTCATGATCGTGGCCGGCACACTCACCGCCAACAAGATGGCCCCAGCGC-TCGCAAGGT
     2.7 -p                              #)/5:@FLRX^djou{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{uoiiikkkkkqw}~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ysmmmsy~~~~~~~|vpjd^YSMGA;55542.(&$#
     2.8 +p                              #)/5:@FLRX^djou{~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~{uoiiikkkkkqw|~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ysmmmsy~~~~~~~|vpjd^YSMGA;55542.(&$#
     2.9  
    2.10  a score=489 mismap=1e-10
    2.11  s chr19 1388829 108 + 58617616 CAGCACCCAGCCTGCCCTGCCAAAGGCCAG-AGCCGTGGCTCCCAAACCCAGCAGCCGGGGCGAGTATG-TGGTGGCCAAGCTGGATGACCTCGTCAACTGGGCCCGCCG