From b465d299e7afcfb4abb7cc76b417556b4ffe8a11 Mon Sep 17 00:00:00 2001 From: martin-raden Date: Wed, 24 Jan 2024 15:49:33 +0100 Subject: [PATCH] BUGFIX: missing isAccessible() checks added --- .../PredictorMfe2dHeuristicSeedExtension.cpp | 56 ++++++++++++------- src/IntaRNA/PredictorMfe2dSeedExtension.cpp | 32 +++++++++-- .../PredictorMfe2dSeedExtensionRIblast.cpp | 24 ++++++-- .../PredictorMfeEns2dSeedExtension.cpp | 26 +++++++-- src/IntaRNA/general.h | 2 + 5 files changed, 107 insertions(+), 33 deletions(-) diff --git a/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfe2dHeuristicSeedExtension.cpp index 1d14cee1..2b4c4d28 100644 --- 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 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 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 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 index f150ff4c..50cc1f45 100644 --- 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 index ec047052..4b195653 100644 --- 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 + //! central flag whether or not debug mode is enabled #define INTARNA_IN_DEBUG_MODE ((defined(_DEBUG)) || (!defined (NDEBUG)))