BUGFIX: missing isAccessible() checks added

BackofenLab · Jan 24, 2024 · b465d29 · b465d29
1 parent 26ef8e7
commit b465d29
Show file tree

Hide file tree

Showing 5 changed files with 107 additions and 33 deletions.
diff --git a/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp
@@ -142,13 +142,22 @@ fillHybridE_right( const size_t sj1, const size_t sj2
 			// init current cell (0 if just left (i1,i2) base pair)
 			curMinE = (sj1==j1 && sj2==j2) ? 0 : E_INF;
 
+			// skip if not accessible
 			// check if complementary
-			if( sj1<j1 && sj2<j2 && energy.areComplementary(j1,j2) ) {
+			if ( energy.isAccessible1(j1)
+				&& energy.isAccessible2(j2)
+				&& sj1<j1
+				&& sj2<j2
+				&& energy.areComplementary(j1,j2) )
+			{
 
 				// left-stacking of j if no-LP
 				if (outConstraint.noLP) {
 					// skip if no stacking possible
-					if (!energy.areComplementary(j1-noLpShift,j2-noLpShift)) {
+					if (  !energy.areComplementary(j1-noLpShift,j2-noLpShift)
+						| !energy.isAccessible1(j1-noLpShift)
+						| !energy.isAccessible2(j2-noLpShift))
+					{
 						continue;
 					}
 					// get stacking energy to avoid recomputation in recursion below
@@ -163,19 +172,19 @@ fillHybridE_right( const size_t sj1, const size_t sj2
 				for (k1=j1-noLpShift; k1-- > sj1; ) {
 					// ensure maximal loop length
 					if (j1-noLpShift-k1 > energy.getMaxInternalLoopSize1()+1) break;
-				for (k2=j2-noLpShift; k2-- > sj2; ) {
-					// ensure maximal loop length
-					if (j2-noLpShift-k2 > energy.getMaxInternalLoopSize2()+1) break;
-					// check if (k1,k2) are valid left boundary
-					if ( E_isNotINF( hybridE_right(k1-sj1,k2-sj2) ) ) {
-						curMinE = std::min( curMinE,
-								(hybridE_right(k1-sj1,k2-sj2) // left part
-								+ energy.getE_interLeft(k1,j1-noLpShift,k2,j2-noLpShift) // loop
-								+ iStackE) // right stack if no-LP
-								);
+					for (k2=j2-noLpShift; k2-- > sj2; ) {
+						// ensure maximal loop length
+						if (j2-noLpShift-k2 > energy.getMaxInternalLoopSize2()+1) break;
+						// check if (k1,k2) are valid left boundary
+						if ( E_isNotINF( hybridE_right(k1-sj1,k2-sj2) ) ) {
+							curMinE = std::min( curMinE,
+									(hybridE_right(k1-sj1,k2-sj2) // left part
+									+ energy.getE_interLeft(k1,j1-noLpShift,k2,j2-noLpShift) // loop
+									+ iStackE) // right stack if no-LP
+									);
+						}
 					}
 				}
-				}
 
 				// update mfe if needed
 				if (E_isNotINF(curMinE)) {
@@ -184,9 +193,9 @@ fillHybridE_right( const size_t sj1, const size_t sj2
 					// update mfe for seed+rightExt
 					updateOptima( si1,j1,si2,j2, seedE + curMinE + energy.getE_init(),true,true);
 				}
-			}
-		}
-	}
+			} // if complementary
+		} // for j2
+	} // for j1
 
 }
 
@@ -222,6 +231,7 @@ fillHybridE_left( const size_t si1, const size_t si2 )
 
 	// iterate over all window starts
 	for (size_t l1=0; l1 < hybridE_left.size1(); l1++) {
+
 		for (size_t l2=0; l2 < hybridE_left.size2(); l2++) {
 			i1 = si1-l1;
 			i2 = si2-l2;
@@ -230,13 +240,21 @@ fillHybridE_left( const size_t si1, const size_t si2 )
 			E_type & curMinE = hybridE_left(l1,l2);
 			// init cell
 			curMinE = (i1==si1 && i2==si2) ? energy.getE_init() : E_INF;
+			// skip if not accessible
 			// check if complementary
-			if( i1<si1 && i2<si2 && energy.areComplementary(i1,i2) ) {
-
+			if (energy.isAccessible1(i1)
+				&& energy.isAccessible2(i2)
+				&& i1<si1
+				&& i2<si2
+				&& energy.areComplementary(i1,i2) )
+			{
 				// left-stacking of j if no-LP
 				if (outConstraint.noLP) {
 					// skip if no stacking possible
-					if (!energy.areComplementary(i1+noLpShift,i2+noLpShift)) {
+					if (  !energy.areComplementary(i1+noLpShift,i2+noLpShift)
+						| !energy.isAccessible1(i1+noLpShift)
+						| !energy.isAccessible2(i2+noLpShift))
+					{
 						continue;
 					}
 					// get stacking energy to avoid recomputation in recursion below

diff --git a/src/IntaRNA/PredictorMfe2dSeedExtension.cpp b/src/IntaRNA/PredictorMfe2dSeedExtension.cpp
@@ -144,8 +144,10 @@ fillHybridE_left( const size_t si1, const size_t si2 )
 	const size_t noLpShift = outConstraint.noLP ? 1 : 0;
 	E_type iStackE = E_type(0);
 
-	// iterate over all window starts j1 (seq1) and j2 (seq2)
+	// iterate over all window starts i1 (seq1) and i2 (seq2)
+	// via their distance l1/l2 to the left seed boundary si1/si2
 	for (size_t l1=0; l1 < hybridE_left.size1(); l1++) {
+
 		for (size_t l2=0; l2 < hybridE_left.size2(); l2++) {
 			i1 = si1-l1;
 			i2 = si2-l2;
@@ -154,13 +156,22 @@ fillHybridE_left( const size_t si1, const size_t si2 )
 			E_type & curE = hybridE_left(l1,l2);
 			// init current cell (e_init if just left (i1,i2) base pair; assuming seed is internally stacked on the left end if noLP)
 			curE = (i1==si1 && i2==si2) ? energy.getE_init() : E_INF;
+			// skip if not accessible
 			// check if complementary
-			if( i1<si1 && i2<si2 && energy.areComplementary(i1,i2) ) {
+			if (energy.isAccessible1(i1)
+				&& energy.isAccessible2(i2)
+				&& i1<si1
+				&& i2<si2
+				&& energy.areComplementary(i1,i2) )
+			{
 
 				// right-stacking of i if no-LP
 				if (outConstraint.noLP) {
 					// skip if no stacking possible
-					if (!energy.areComplementary(i1+noLpShift,i2+noLpShift)) {
+					if (  !energy.areComplementary(i1+noLpShift,i2+noLpShift)
+						| !energy.isAccessible1(i1+noLpShift)
+						| !energy.isAccessible2(i2+noLpShift))
+					{
 						continue;
 					}
 					// get stacking energy to avoid recomputation in recursion below
@@ -177,6 +188,7 @@ fillHybridE_left( const size_t si1, const size_t si2 )
 						// ensure maximal loop length
 						if (k2-i2-noLpShift > energy.getMaxInternalLoopSize2()+1) break;
 						// check if (k1,k2) are valid left boundary
+		if ( energy.isAccessible1(k1) && energy.isAccessible2(k2) )
 						if ( E_isNotINF( hybridE_left(si1-k1,si2-k2) ) ) {
 							curE = std::min( curE,
 									(iStackE
@@ -224,13 +236,22 @@ fillHybridE_right( const size_t sj1, const size_t sj2 )
 			// init current cell (0 if just left (i1,i2) base pair)
 			curE = (sj1==j1 && sj2==j2) ? 0 : E_INF;
 
+			// skip if not accessible
 			// check if complementary
-			if( sj1<j1 && sj2<j2 && energy.areComplementary(j1,j2) ) {
+			if (energy.isAccessible1(j1)
+				&& energy.isAccessible2(j2)
+				&& sj1<j1
+				&& sj2<j2
+				&& energy.areComplementary(j1,j2) )
+			{
 
 				// left-stacking of j if no-LP
 				if (outConstraint.noLP) {
 					// skip if no stacking possible
-					if (!energy.areComplementary(j1-noLpShift,j2-noLpShift)) {
+					if (!energy.areComplementary(j1-noLpShift,j2-noLpShift)
+						| !energy.isAccessible1(j1-noLpShift)
+						| !energy.isAccessible2(j2-noLpShift))
+					{
 						continue;
 					}
 					// get stacking energy to avoid recomputation in recursion below
@@ -249,6 +270,7 @@ fillHybridE_right( const size_t sj1, const size_t sj2 )
 					// ensure maximal loop length
 					if (j2-noLpShift-k2 > energy.getMaxInternalLoopSize2()+1) break;
 					// check if (k1,k2) are valid left boundary
+if ( energy.isAccessible1(k1) && energy.isAccessible2(k2) )
 					if ( E_isNotINF( hybridE_right(k1-sj1,k2-sj2) ) ) {
 						curE = std::min( curE,
 								(hybridE_right(k1-sj1,k2-sj2)

diff --git a/src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp b/src/IntaRNA/PredictorMfe2dSeedExtensionRIblast.cpp
@@ -164,7 +164,10 @@ parallelExtension( PredictorMfe2dSeedExtensionRIblast::ExtendedSeed & seed
 	// extend left
 	while (seed.i1 > 1 && seed.i2 > 1) {
 		// todo acc1/acc2 maxLength() termination
-		if (energy.areComplementary(seed.i1-1,seed.i2-1)) {
+		if (energy.isAccessible1(seed.i1-1)
+			&& energy.isAccessible2(seed.i2-1)
+			&& energy.areComplementary(seed.i1-1,seed.i2-1))
+		{
 			E_type newEnergy = seed.energy + energy.getE_interLeft(seed.i1-1,i1min,seed.i2-1,i2min);
 			if (newEnergy < seed.energy) {
 				seed.energy = newEnergy;
@@ -187,7 +190,10 @@ parallelExtension( PredictorMfe2dSeedExtensionRIblast::ExtendedSeed & seed
 	// extend right
 	while (seed.j1 < max_extension1-1 && seed.j2 < max_extension2-1) {
 		// todo acc1/acc2 maxLength() termination
-		if (energy.areComplementary(seed.j1+1,seed.j2+1)) {
+		if (energy.isAccessible1(seed.j1+1)
+			&& energy.isAccessible2(seed.j2+1)
+			&& energy.areComplementary(seed.j1+1,seed.j2+1))
+		{
 			E_type newEnergy = seed.energy + energy.getE_interLeft(j1min,seed.j1+1,j2min,seed.j2+1);
 			if (newEnergy < seed.energy) {
 				seed.energy = newEnergy;
@@ -230,7 +236,12 @@ fillHybridE_left( const size_t j1, const size_t j2 )
 			hybridE_left(i1,i2) = i1==0 && i2==0 ? energy.getE_init() : E_INF;
 
 			// check if complementary
-			if( i1>0 && i2>0 && energy.areComplementary(j1-i1,j2-i2) ) {
+			if( i1>0
+				&& i2>0
+				&& energy.isAccessible1(j1-i1)
+				&& energy.isAccessible2(j2-i2)
+				&& energy.areComplementary(j1-i1,j2-i2) )
+			{
 				curMinE = E_INF;
 
 				// check all combinations of decompositions into (i1,i2)..(k1,k2)-(j1,j2)
@@ -286,7 +297,12 @@ fillHybridE_right( const size_t i1, const size_t i2 )
 			hybridE_right(j1,j2) = j1==0 && j2==0 ? 0 : E_INF;
 
 			// check if complementary
-			if( j1>i1 && j2>i2 && energy.areComplementary(i1+j1,i2+j2) ) {
+			if( j1>i1
+				&& j2>i2
+				&& energy.isAccessible1(i1+j1)
+				&& energy.isAccessible2(i2+j2)
+				&& energy.areComplementary(i1+j1,i2+j2) )
+			{
 				curMinE = E_INF;
 
 				// check all combinations of decompositions into (i1,i2)..(k1,k2)-(j1,j2)

diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp
@@ -210,12 +210,20 @@ fillHybridZ_left( const size_t si1, const size_t si2 )
 			curZ = (i1==si1 && i2==si2) ? energy.getBoltzmannWeight(energy.getE_init()) : 0.0;
 
 			// check if complementary (use global sequence indexing)
-			if( i1<si1 && i2<si2 && energy.areComplementary(i1,i2) ) {
+			if( i1<si1
+				&& i2<si2
+				&& energy.isAccessible1(i1)
+				&& energy.isAccessible2(i2)
+				&& energy.areComplementary(i1,i2) )
+			{
 
 				// right-stacking of i if no-LP
 				if (outConstraint.noLP) {
 					// skip if no stacking possible
-					if (!energy.areComplementary(i1+noLpShift,i2+noLpShift)) {
+					if (!energy.areComplementary(i1+noLpShift,i2+noLpShift)
+						| !energy.isAccessible1(i1+noLpShift)
+						| !energy.isAccessible2(i2+noLpShift) )
+					{
 						continue;
 					}
 					// get stacking energy to avoid recomputation in recursion below
@@ -317,7 +325,7 @@ fillHybridZ_left( const size_t si1, const size_t si2 )
 
 							}
 							curZ -= correctionTerm;
-							// sanity ensurence
+							// sanity insurance
 							if (curZ < 0) {
 								curZ = Z_type(0.0);
 							}
@@ -365,12 +373,20 @@ fillHybridZ_right( const size_t sj1, const size_t sj2 )
 			curZ = sj1==j1 && sj2==j2 ? energy.getBoltzmannWeight(0.0) : 0.0;
 
 			// check if complementary free base pair
-			if( sj1<j1 && sj2<j2 && energy.areComplementary(j1,j2) ) {
+			if( sj1<j1
+				&& sj2<j2
+				&& energy.isAccessible1(j1)
+				&& energy.isAccessible2(j2)
+				&& energy.areComplementary(j1,j2) )
+			{
 
 				// left-stacking of j if no-LP
 				if (outConstraint.noLP) {
 					// skip if no stacking possible
-					if (!energy.areComplementary(j1-noLpShift,j2-noLpShift)) {
+					if (!energy.areComplementary(j1-noLpShift,j2-noLpShift)
+						| !energy.isAccessible1(j1-noLpShift)
+						| !energy.isAccessible2(j2-noLpShift) )
+					{
 						continue;
 					}
 					// get stacking energy to avoid recomputation in recursion below

diff --git a/src/IntaRNA/general.h b/src/IntaRNA/general.h
@@ -11,6 +11,8 @@
 
 #include "IntaRNA/intarna_config.h"
 
+#include <cassert>
+
 //! central flag whether or not debug mode is enabled
 #define INTARNA_IN_DEBUG_MODE ((defined(_DEBUG)) || (!defined (NDEBUG)))