bugfix #190 outNoLP : missing recursion cases

ChangeLog

@@ -16,6 +16,8 @@
   (thanks to Sabine Reisser)
 - BUGFIX: explicit seed encodings within last 7 nucleotides (seedBP) were 
   ignored (thanks to Sebastian Holler)
+- BUGFIX: outNoLP option was not correctly implemented (missing recursion cases)
+  and was thus missing interactions
 
 ################################################################################  
 ################################################################################
@@ -33,6 +35,14 @@
    * seedBP now set to getSeedMinBP() if explicit seeds present
  * bin/IntaRNA :
    * exception information now motivates to send input along with report
+ + docs/recursions : recursion depictions for sanity checks of implementations
+   + PredictorMfe2d (+outNoLP)
+   + PredictorMfe2dSeed (+outNoLP)
+ * IntaRNA/PredictorMfe2dSeed :
+ * IntaRNA/PredictorMfe2dHeuristic :
+ * IntaRNA/PredictorMfe2dHeuristicSeed :
+ * IntaRNA/PredictorMfeEns2dHeuristic :
+   * bugfix outNoLP : missing recursion cases 
 
 201210 Martin Raden
  * IntaRNA/AccessibilityFromStream :

src/IntaRNA/PredictorMfe2d.cpp

@@ -100,6 +100,10 @@ fillHybridE( const size_t j1, const size_t j2
 		throw std::runtime_error("PredictorMfe2d::fillHybridE() : i1init > j1 : "+toString(i1init)+" > "+toString(j1));
 	if (i2init > j2)
 		throw std::runtime_error("PredictorMfe2d::fillHybridE() : i2init > j2 : "+toString(i2init)+" > "+toString(j2));
+	if (!energy.isAccessible2(j1))
+		throw std::runtime_error("PredictorMfe2d::fillHybridE() : !energy.isAccessible2(j1) : "+toString(j1));
+	if (!energy.isAccessible2(j2))
+		throw std::runtime_error("PredictorMfe2d::fillHybridE() : !energy.isAccessible2(j2) : "+toString(j2));
 #endif
 
 	// get minimal start indices heeding max interaction length
@@ -166,6 +170,7 @@ fillHybridE( const size_t j1, const size_t j2
 							iStackE = energy.getE_interLeft(i1,i1+noLpShift,i2,i2+noLpShift);
 
 							// init with stacking only
+							// or stacking with right extension
 							curMinE = iStackE + ((w1==2&&w2==2) ? energy.getE_init() : hybridE_pq(i1+noLpShift, i2+noLpShift) );
 						} else {
 							//
@@ -274,7 +279,7 @@ traceBack( Interaction & interaction  )
 			// get stacking term
 			iStackE = energy.getE_interLeft(i1,i1+noLpShift, i2,i2+noLpShift);
 
-			// check just stacking
+			// check just stacking with right extension
 			if ( E_equal( curE, iStackE + hybridE_pq(i1+noLpShift,i2+noLpShift) )) {
 				// store bp
 				interaction.basePairs.push_back( energy.getBasePair(i1+noLpShift,i2+noLpShift) );

src/IntaRNA/PredictorMfe2dHeuristic.cpp

@@ -124,7 +124,8 @@ fillHybridE()
 
 				// set to interaction initiation with according boundary
 				// if valid right boundary
-				if (!outConstraint.noGUend || !energy.isGU(i1+noLpShift,i2+noLpShift))
+				if (E_isNotINF(iStackE)
+						&& (!outConstraint.noGUend || !energy.isGU(i1+noLpShift,i2+noLpShift)))
 				{
 					*curCell = BestInteractionE(iStackE + energy.getE_init(), i1+noLpShift, i2+noLpShift);
 					// current best total energy value (covers to far E_init only)
@@ -133,6 +134,40 @@ fillHybridE()
 					updateZall( i1,curCell->j1,i2,curCell->j2, curCellEtotal, false );
 				}
 
+				/////////////////////////////////////////
+				// check direct extension to the right of the noLP stacking
+				/////////////////////////////////////////
+				if(outConstraint.noLP)
+				{
+
+					// direct cell access (const)
+					rightExt = &(hybridE(i1+noLpShift,i2+noLpShift));
+					// check if right side can pair
+					// check if interaction length is within boundary
+					if (E_isNotINF(rightExt->val)
+						&& (rightExt->j1 +1 -i1) <= energy.getAccessibility1().getMaxLength()
+						&& (rightExt->j2 +1 -i2) <= energy.getAccessibility2().getMaxLength() )
+					{
+						// compute energy direct extension with stacking
+						curE = iStackE + rightExt->val;
+						// check if this combination yields better energy
+						curEtotal = energy.getE(i1,rightExt->j1,i2,rightExt->j2,curE);
+						// update best extension
+						if ( curEtotal < curCellEtotal )
+						{
+							// update current best for this left boundary
+							// copy right boundary
+							*curCell = *rightExt;
+							// set new energy
+							curCell->val = curE;
+							// store total energy to avoid recomputation
+							curCellEtotal = curEtotal;
+						}
+						// update Zall
+						updateZall( i1,rightExt->j1,i2,rightExt->j2, curEtotal, false );
+					}
+				}
+
 				// iterate over all loop sizes w1 (seq1) and w2 (seq2) (minus 1)
 				for (w1=1; w1-1 <= energy.getMaxInternalLoopSize1() && i1+w1+noLpShift<hybridE.size1(); w1++) {
 				for (w2=1; w2-1 <= energy.getMaxInternalLoopSize2() && i2+w2+noLpShift<hybridE.size2(); w2++) {

src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp

@@ -167,12 +167,15 @@ fillHybridE()
 					E_type seedE = seedHandler.getSeedE(i1,i2);
 					if (sj1 < hybridE.size1() && sj2 < hybridE.size2()) {
 
-						// get energy of seed only explicitly
-						curE = seedE + energy.getE_init();
-						// check if this combination yields better energy
-						curEseedtotal = energy.getE(i1,sj1,i2,sj2,curE);
-						// update Zall for seed only (if otherwise stacking enforced)
-						updateZall( i1,sj1,i2,sj2, curEseedtotal, false );
+						// exclude single-bp seeds in noLP mode
+						if (sj1 >= i1+noLpShift) {
+							// get energy of seed only explicitly
+							curE = seedE + energy.getE_init();
+							// check if this combination yields better energy
+							curEseedtotal = energy.getE(i1,sj1,i2,sj2,curE);
+							// update Zall for seed only (if otherwise stacking enforced)
+							updateZall( i1,sj1,i2,sj2, curEseedtotal, false );
+						}
 
 						// for noLP : check for explicit interior loop after seed
 						// assumption: seed fulfills noLP
@@ -252,6 +255,73 @@ fillHybridE()
 
 				// if !noLP or stacking bp possible
 				if (E_isNotINF(iStackE))  {
+
+					if (outConstraint.noLP) {
+						/////////////////////////////////////////
+						// check direct extension to the right of the noLP stacking
+						/////////////////////////////////////////
+
+						//////////////////////////////////////////////////////////
+						// update hybridE without seed constraint
+						//////////////////////////////////////////////////////////
+
+						// direct cell access (const)
+						rightExt = &(hybridE(i1+noLpShift,i2+noLpShift));
+
+						// check if right side can pair
+						// check if interaction length is within boundary
+						if ( E_isNotINF(rightExt->val)
+							&& (rightExt->j1 +1 -i1) <= energy.getAccessibility1().getMaxLength()
+							&& (rightExt->j2 +1 -i2) <= energy.getAccessibility2().getMaxLength() )
+						{
+							// compute energy for direct stack extension
+							curE = iStackE + rightExt->val;
+							// check if this combination yields better energy
+							curEtotal = energy.getE(i1,rightExt->j1,i2,rightExt->j2,curE);
+							if ( curEtotal < curCellEtotal )
+							{
+								// update current best for this left boundary
+								// copy right boundary
+								*curCell = *rightExt;
+								// set new energy
+								curCell->val = curE;
+								// store total energy to avoid recomputation
+								curCellEtotal = curEtotal;
+							}
+						}
+
+						//////////////////////////////////////////////////////////
+						// update hybridE_seed including seed constraint
+						//////////////////////////////////////////////////////////
+
+						// direct cell access to right side end of loop (seed has to be to the right of it)
+						rightExt = &(hybridE_seed(i1+noLpShift,i2+noLpShift));
+						// check if right side of loop can pair
+						// check if interaction length is within boundary
+						if (E_isNotINF(rightExt->val)
+							&& (rightExt->j1 +1 -i1) <= energy.getAccessibility1().getMaxLength()
+							&& (rightExt->j2 +1 -i2) <= energy.getAccessibility2().getMaxLength() )
+						{
+							// compute energy for direct stack extension
+							curE = iStackE + rightExt->val;
+							// check if this combination yields better energy
+							curEseedtotal = energy.getE(i1,rightExt->j1,i2,rightExt->j2,curE);
+							if ( curEseedtotal < curCellSeedEtotal )
+							{
+								// update current best for this left boundary
+								// copy right boundary
+								*curCellSeed = *rightExt;
+								// set new energy
+								curCellSeed->val = curE;
+								// store overall energy
+								curCellSeedEtotal = curEseedtotal;
+							}
+							// avoid Z update; otherwise double-counting of interactions
+							//   --> that way, underestimation of Z
+						}
+					}
+
+
 					// iterate over all loop sizes w1 (seq1) and w2 (seq2) (minus 1)
 					for (w1=1; w1-1 <= energy.getMaxInternalLoopSize1() && i1+noLpShift+w1<hybridE.size1(); w1++) {
 					for (w2=1; w2-1 <= energy.getMaxInternalLoopSize2() && i2+noLpShift+w2<hybridE.size2(); w2++) {
@@ -409,6 +479,26 @@ traceBack( Interaction & interaction )
 
 		const BestInteractionE * curCell = NULL;
 		bool traceNotFound = true;
+
+		// check for stacking on the left with direct right extension
+		if (outConstraint.noLP) {
+			curCell = &(hybridE_seed(i1+noLpShift,i2+noLpShift));
+			// check if right boundary is equal (part of the heuristic)
+			if ( curCell->j1 == j1 && curCell->j2 == j2 &&
+					// and energy is the source of curE
+					E_equal( curE, (iStackE + curCell->val ) ) )
+			{
+				// stop searching
+				traceNotFound = false;
+				// store bp
+				interaction.basePairs.push_back( energy.getBasePair(i1+noLpShift,i2+noLpShift) );
+				// trace right part of stacking
+				i1+=noLpShift;
+				i2+=noLpShift;
+				curE = curCell->val;
+			}
+		}
+
 		// check all combinations of decompositions into (i1,i2)..(k1,k2)-(j1,j2)
 		for (k1=std::min(j1,i1+energy.getMaxInternalLoopSize1()+1+noLpShift); traceNotFound && k1>i1+noLpShift; k1--) {
 		for (k2=std::min(j2,i2+energy.getMaxInternalLoopSize2()+1+noLpShift); traceNotFound && k2>i2+noLpShift; k2--) {
@@ -454,7 +544,7 @@ traceBack( Interaction & interaction )
 				i2 = k2;
 				break;
 			}
-			// trace remaining base pairs
+			// trace right bulge extension of seed in noLP mode
 			if (outConstraint.noLP) {
 				for (size_t l1=std::min(j1-1,k1+energy.getMaxInternalLoopSize1()+1); traceNotFound && l1>k1; l1--) {
 				for (size_t l2=std::min(j2-1,k2+energy.getMaxInternalLoopSize2()+1); traceNotFound && l2>k2; l2--) {
@@ -477,6 +567,7 @@ traceBack( Interaction & interaction )
 				}} // l1 l2
 			}
 			if (traceNotFound){
+				// as to be direct right extension of seed without bulge
 				curCell = &(hybridE(k1,k2));
 				// sanity check
 				assert( E_equal( curE, (seedE+(curCell->val)) )

src/IntaRNA/PredictorMfe2dSeed.cpp

@@ -192,6 +192,7 @@ fillHybridE( const size_t j1, const size_t j2
 							iStackE = energy.getE_interLeft(i1,i1+noLpShift,i2,i2+noLpShift);
 
 							// init with stacking only
+							// or stacking with right extension
 							curMinE = iStackE + ((w1==2&&w2==2) ? energy.getE_init() : hybridE_pq(i1+noLpShift, i2+noLpShift) );
 						} else {
 							//
@@ -214,7 +215,8 @@ fillHybridE( const size_t j1, const size_t j2
 								curMinEseed = std::min( curMinEseed, seedHandler.getSeedE(i1,i2) + hybridE_pq(k1,k2) );
 							} else
 							// just the seed up to right boundary (explicit noLP handling)
-							if (k1 == j1 && k2 == j2) {
+							// ensure minimal seed length in noLP mode
+							if (k1 == j1 && k2 == j2 && k1>=i1+noLpShift) {
 								curMinEseed = std::min( curMinEseed, seedHandler.getSeedE(i1,i2) + energy.getE_init() );
 							}
 							// handle interior loops after seeds in noLP-mode
@@ -238,6 +240,13 @@ fillHybridE( const size_t j1, const size_t j2
 						}
 					}
 
+					// handle direct left-stacking in noLP-mode
+					if ( outConstraint.noLP) {
+						if ( E_isNotINF( hybridE_pq_seed(i1+noLpShift,i2+noLpShift) ) ) {
+							curMinEseed = std::min( curMinEseed, (iStackE + hybridE_pq_seed(i1+noLpShift,i2+noLpShift) ) );
+						}
+					}
+
 					// check all combinations of decompositions into (i1,i2)..(k1,k2)-(j1,j2)
 					if (w1 > 2 && w2 > 2) {
 						for (k1=std::min(j1-1,i1+energy.getMaxInternalLoopSize1()+1+noLpShift); k1>i1+noLpShift; k1--) {
@@ -433,6 +442,26 @@ traceBack( Interaction & interaction )
 				}
 			}
 
+			// explicit check of direct left-end stacking (in noLP mode)
+			if (outConstraint.noLP) {
+				if ( E_isNotINF( hybridE_pq_seed(i1+noLpShift,i2+noLpShift) ) ) {
+					// check if correct split
+					if (E_equal ( curE,
+							(iStackE + hybridE_pq_seed(i1+noLpShift,i2+noLpShift) )
+							) )
+					{
+						// store stacked base pair
+						interaction.basePairs.push_back( energy.getBasePair(i1+noLpShift,i2+noLpShift) );
+						// update trace back boundary
+						i1+=noLpShift;
+						i2+=noLpShift;
+						curE= hybridE_pq_seed(i1,i2);
+						// start next iteration if trace was found
+						continue;
+					}
+				}
+			}
+
 			// check all interval splits
 			if ( (j1-i1) > 1 && (j2-i2) > 1) {
 				// check all combinations of decompositions into (i1,i2)..(k1,k2)-(j1,j2)

src/IntaRNA/PredictorMfeEns2dHeuristic.cpp

@@ -133,6 +133,46 @@ fillHybridZ()
 					curCellEtotal = energy.getE(i1,i1+noLpShift, i2,i2+noLpShift ,energy.getE(curCell->val));
 					// update overall partition function information for initial bps only
 					updateZ( i1,curCell->j1, i2,curCell->j2, curCell->val, true );
+
+					if(outConstraint.noLP) {
+						/////////////////////////////////////////
+						// check direct extension to the right of the noLP stacking
+						/////////////////////////////////////////
+
+						// direct cell access (const)
+						rightExt = &(hybridZ(i1+noLpShift,i2+noLpShift));
+						// check if right side can pair
+						if (Z_equal(rightExt->val, 0.0)) {
+							continue;
+						}
+						// check if interaction length is within boundary
+						if ( (rightExt->j1 +1 -i1) > energy.getAccessibility1().getMaxLength()
+							|| (rightExt->j2 +1 -i2) > energy.getAccessibility2().getMaxLength() )
+						{
+							continue;
+						}
+
+						// compute Z for direct extension with stacking
+						curZ = iStackZ * rightExt->val;
+
+						// update overall partition function information for current right extension
+						updateZ( i1,rightExt->j1, i2,rightExt->j2, curZ, true );
+
+						// check if this combination yields better energy
+						curEtotal = energy.getE(i1,rightExt->j1, i2,rightExt->j2, energy.getE(curZ));
+
+						// update best right extension for (i1,i2) in curCell
+						if ( curEtotal < curCellEtotal )
+						{
+							// update current best for this left boundary
+							// copy right boundary
+							*curCell = *rightExt;
+							// set new partition function
+							curCell->val = curZ;
+							// store total energy to avoid recomputation
+							curCellEtotal = curEtotal;
+						}
+					}
 				}