diff --git a/src/IntaRNA/HelixHandlerNoBulgeMax.cpp b/src/IntaRNA/HelixHandlerNoBulgeMax.cpp index 7c98530..7d89a50 100644 --- a/src/IntaRNA/HelixHandlerNoBulgeMax.cpp +++ b/src/IntaRNA/HelixHandlerNoBulgeMax.cpp @@ -58,9 +58,7 @@ fillHelix(const size_t i1min, const size_t i1max, const size_t i2min, const size i2 = i2start+o; // check if valid base pair - if( energy.isAccessible1( i1 ) - && energy.isAccessible2( i2 ) - && energy.areComplementary( i1, i2 )) + if (energy.areComplementary( i1, i2 )) { // start new canonical helix information if (E_isINF(curHelixE)) { @@ -232,9 +230,7 @@ fillHelixSeed(const size_t i1min, const size_t i1max, const size_t i2min, const j1 = seedEnd1+trailingL; j2 = seedEnd2+trailingL; // check if trailing based pairs are possible, otherwise stop computation - if (!(energy.isAccessible1(j1) - && energy.isAccessible2(j2) - && energy.areComplementary(j1,j2))) + if (!energy.areComplementary(j1,j2)) { break; } @@ -258,9 +254,7 @@ fillHelixSeed(const size_t i1min, const size_t i1max, const size_t i2min, const // check if leading based pairs are possible, otherwise stop computation if (leadingBP > 0 - &&!(energy.isAccessible1(i1) - && energy.isAccessible2(i2) - && energy.areComplementary(i1,i2))) + && !energy.areComplementary(i1,i2)) { break; } diff --git a/src/IntaRNA/PredictorMfe2d.cpp b/src/IntaRNA/PredictorMfe2d.cpp index 2c4db30..3edacf8 100644 --- a/src/IntaRNA/PredictorMfe2d.cpp +++ b/src/IntaRNA/PredictorMfe2d.cpp @@ -131,9 +131,7 @@ fillHybridE( const size_t j1, const size_t j2 hybridE_pq(i1,i2) = E_INF; // check if this cell is to be computed (!=E_INF) - if( energy.isAccessible1(i1) - && energy.isAccessible2(i2) - && energy.areComplementary(i1,i2) + if( energy.areComplementary(i1,i2) ) { // w2 = interaction width in seq2 @@ -162,9 +160,7 @@ fillHybridE( const size_t j1, const size_t j2 } else { // no lp allowed // check if right-side stacking of (i1,i2) is possible - if (energy.isAccessible1(i1+noLpShift) - && energy.isAccessible2(i2+noLpShift) - && energy.areComplementary(i1+noLpShift,i2+noLpShift)) + if (energy.areComplementary(i1+noLpShift,i2+noLpShift)) { // get stacking term to avoid recomputation iStackE = energy.getE_interLeft(i1,i1+noLpShift,i2,i2+noLpShift); diff --git a/src/IntaRNA/PredictorMfe2dHelixBlockHeuristic.cpp b/src/IntaRNA/PredictorMfe2dHelixBlockHeuristic.cpp index b8c0d03..e7bafe7 100644 --- a/src/IntaRNA/PredictorMfe2dHelixBlockHeuristic.cpp +++ b/src/IntaRNA/PredictorMfe2dHelixBlockHeuristic.cpp @@ -85,9 +85,7 @@ predict( const IndexRange & r1 for (i2=0; i2 1 && seed.i2 > 1) { // todo acc1/acc2 maxLength() termination - if (energy.isAccessible1(seed.i1-1) - && energy.isAccessible2(seed.i2-1) - && energy.areComplementary(seed.i1-1,seed.i2-1)) + if (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) { @@ -190,9 +188,7 @@ parallelExtension( PredictorMfe2dSeedExtensionRIblast::ExtendedSeed & seed // extend right while (seed.j1 < max_extension1-1 && seed.j2 < max_extension2-1) { // todo acc1/acc2 maxLength() termination - if (energy.isAccessible1(seed.j1+1) - && energy.isAccessible2(seed.j2+1) - && energy.areComplementary(seed.j1+1,seed.j2+1)) + if (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) { @@ -238,8 +234,6 @@ fillHybridE_left( const size_t j1, const size_t j2 ) // check if complementary if( i1>0 && i2>0 - && energy.isAccessible1(j1-i1) - && energy.isAccessible2(j2-i2) && energy.areComplementary(j1-i1,j2-i2) ) { curMinE = E_INF; @@ -299,8 +293,6 @@ fillHybridE_right( const size_t i1, const size_t i2 ) // check if complementary if( j1>i1 && j2>i2 - && energy.isAccessible1(i1+j1) - && energy.isAccessible2(i2+j2) && energy.areComplementary(i1+j1,i2+j2) ) { curMinE = E_INF; diff --git a/src/IntaRNA/PredictorMfeEns2d.cpp b/src/IntaRNA/PredictorMfeEns2d.cpp index 145d355..2bc4b6e 100644 --- a/src/IntaRNA/PredictorMfeEns2d.cpp +++ b/src/IntaRNA/PredictorMfeEns2d.cpp @@ -128,9 +128,7 @@ fillHybridZ( const size_t j1, const size_t j2 hybridZ(i1,i2) = Z_type(0.0); // check if this cell is to be computed (!=E_INF) - if( energy.isAccessible1(i1) - && energy.isAccessible2(i2) - && energy.areComplementary(i1,i2) + if( energy.areComplementary(i1,i2) ) { // w2 = interaction width in seq2 @@ -159,9 +157,7 @@ fillHybridZ( const size_t j1, const size_t j2 } else { // no lp allowed // check if right-side stacking of (i1,i2) is possible - if (energy.isAccessible1(i1+noLpShift) - && energy.isAccessible2(i2+noLpShift) - && energy.areComplementary(i1+noLpShift,i2+noLpShift)) + if (energy.areComplementary(i1+noLpShift,i2+noLpShift)) { // get stacking term to avoid recomputation iStackZ = energy.getBoltzmannWeight(energy.getE_interLeft(i1,i1+noLpShift,i2,i2+noLpShift)); diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp index 1383be4..49cdb37 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp @@ -104,17 +104,13 @@ fillHybridZ() curCellEtotal = E_INF; // check if positions can form interaction - if ( energy.isAccessible1(i1) - && energy.isAccessible2(i2) - && energy.areComplementary(i1,i2) ) + if (energy.areComplementary(i1,i2) ) { // no lp allowed if (noLpShift != 0) { // check if right-side stacking of (i1,i2) is possible if ( i1+noLpShift < energy.size1() && i2+noLpShift < energy.size2() - && energy.isAccessible1(i1+noLpShift) - && energy.isAccessible2(i2+noLpShift) && energy.areComplementary(i1+noLpShift,i2+noLpShift)) { // get stacking term to avoid recomputation diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp index 05a4f85..d5adb5c 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp @@ -213,17 +213,13 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) // check if complementary (use global sequence indexing) if( i1= energy.getED1( i1,i1 ) && seedConstraint.getMaxED() >= energy.getED2( i2,i2 ) diff --git a/tests/SeedHandlerNoBulge_test.cpp b/tests/SeedHandlerNoBulge_test.cpp index 2e2fa14..d993742 100644 --- a/tests/SeedHandlerNoBulge_test.cpp +++ b/tests/SeedHandlerNoBulge_test.cpp @@ -10,6 +10,39 @@ using namespace IntaRNA; +namespace { + +class CountingInteractionEnergyBasePair : public InteractionEnergyBasePair { +public: + using InteractionEnergyBasePair::InteractionEnergyBasePair; + + bool + isAccessible1( const size_t i ) const override + { + accessibilityCalls1++; + return InteractionEnergy::isAccessible1(i); + } + + bool + isAccessible2( const size_t i ) const override + { + accessibilityCalls2++; + return InteractionEnergy::isAccessible2(i); + } + + void + resetAccessibilityCalls() const + { + accessibilityCalls1 = 0; + accessibilityCalls2 = 0; + } + + mutable size_t accessibilityCalls1 = 0; + mutable size_t accessibilityCalls2 = 0; +}; + +} // namespace + TEST_CASE( "SeedHandlerNoBulge", "[SeedHandlerNoBulge]" ) { // setup easylogging++ stuff if not already done @@ -170,3 +203,33 @@ TEST_CASE( "SeedHandlerNoBulge", "[SeedHandlerNoBulge]" ) { } } + +TEST_CASE( "SeedHandler feasibility delegates accessibility once", + "[SeedHandlerNoBulge][AP1Feasibility]" ) +{ + RnaSequence rna1("target", "GA"); + RnaSequence rna2("query", "AC"); + AccessibilityDisabled acc1(rna1, rna1.size(), NULL); + AccessibilityDisabled acc2(rna2, rna2.size(), NULL); + ReverseAccessibility rAcc2(acc2); + CountingInteractionEnergyBasePair energy(acc1, rAcc2); + SeedConstraint constraint(2, 0, 0, 0, E_INF, + Accessibility::ED_UPPER_BOUND, E_INF, IndexRangeList(), + IndexRangeList(), "", false, false, false); + SeedHandlerNoBulge seedHandler(energy, constraint); + + energy.resetAccessibilityCalls(); + REQUIRE(seedHandler.isFeasibleSeedBasePair(0, 0)); + REQUIRE(energy.accessibilityCalls1 == 1); + REQUIRE(energy.accessibilityCalls2 == 1); + + energy.resetAccessibilityCalls(); + REQUIRE_FALSE(seedHandler.isFeasibleSeedBasePair(1, 1)); + REQUIRE(energy.accessibilityCalls1 == 0); + REQUIRE(energy.accessibilityCalls2 == 0); + + energy.resetAccessibilityCalls(); + REQUIRE_FALSE(seedHandler.isFeasibleSeedBasePair(energy.size1(), 0)); + REQUIRE(energy.accessibilityCalls1 == 0); + REQUIRE(energy.accessibilityCalls2 == 0); +}