diff --git a/src/IntaRNA/PredictorMfeEns.cpp b/src/IntaRNA/PredictorMfeEns.cpp index ac16413..3bda524 100644 --- a/src/IntaRNA/PredictorMfeEns.cpp +++ b/src/IntaRNA/PredictorMfeEns.cpp @@ -14,6 +14,7 @@ PredictorMfeEns::PredictorMfeEns( , PredictionTracker * predTracker ) : PredictorMfe(energy,output,predTracker) + , updateZisComplete(false) { } @@ -36,16 +37,17 @@ initZ() //////////////////////////////////////////////////////////////////////////// -void +bool PredictorMfeEns:: -updateZ( const size_t i1, const size_t j1 +addPartitionContribution( const size_t i1, const size_t j1 , const size_t i2, const size_t j2 , const Z_type partZ - , const bool isHybridZ ) + , const bool isHybridZ + , Z_type & partZ_noED ) { // check if something to be done if (Z_equal(partZ,0) || Z_isINF(Zall)) - return; + return false; // Apply the same site filters used for MFE candidates before changing // either the global or boundary-specific partition. @@ -53,17 +55,17 @@ updateZ( const size_t i1, const size_t j1 if (outConstraint.noGUend && (energy.isGU(i1,i2) || energy.isGU(j1,j2))) { - return; + return false; } if (outConstraint.maxED < Accessibility::ED_UPPER_BOUND && (energy.getED1(i1,j1) > outConstraint.maxED || energy.getED2(i2,j2) > outConstraint.maxED)) { - return; + return false; } // handle whether or not partZ includes ED values or not - Z_type partZ_withED = 0, partZ_noED = 0; + Z_type partZ_withED = 0; if (isHybridZ) { #if INTARNA_IN_DEBUG_MODE if ( (std::numeric_limits::max() - (partZ*energy.getBoltzmannWeight(energy.getE(i1,j1,i2,j2, E_type(0))))) <= Zall) { @@ -86,6 +88,33 @@ updateZ( const size_t i1, const size_t j1 // increase overall partition function Zall += partZ_withED; + return true; +} + +//////////////////////////////////////////////////////////////////////////// + +void +PredictorMfeEns:: +updateZ( const size_t i1, const size_t j1 + , const size_t i2, const size_t j2 + , const Z_type partZ + , const bool isHybridZ ) +{ + Z_type partZ_noED = 0; + if (!addPartitionContribution(i1, j1, i2, j2, partZ, isHybridZ, partZ_noED)) { + return; + } + + // Exact 2D prediction finalizes each boundary once. In this scoped mode the + // partition can update the optima immediately instead of entering the map. + // Trackers retain the map path to preserve deferred callback ordering. + if (updateZisComplete && predTracker == NULL) { + if (Z_isNotINF(partZ_noED) && partZ_noED > 0) { + PredictorMfe::updateOptima(i1, j1, i2, j2, + energy.getE(partZ_noED), true, false); + } + return; + } // store partial Z (without ED) Interaction::Boundary key(i1,j1,i2,j2); @@ -101,6 +130,37 @@ updateZ( const size_t i1, const size_t j1 //////////////////////////////////////////////////////////////////////////// +void +PredictorMfeEns:: +updateCompleteZ( const size_t i1, const size_t j1 + , const size_t i2, const size_t j2 + , const Z_type partZ + , const bool isHybridZ ) +{ + class CompleteUpdateScope { + public: + explicit CompleteUpdateScope(bool & mode) + : mode(mode), previous(mode) + { + mode = true; + } + + ~CompleteUpdateScope() + { + mode = previous; + } + + private: + bool & mode; + const bool previous; + }; + + CompleteUpdateScope completeUpdate(updateZisComplete); + updateZ(i1, j1, i2, j2, partZ, isHybridZ); +} + +//////////////////////////////////////////////////////////////////////////// + void PredictorMfeEns:: updateOptimaUsingZ() diff --git a/src/IntaRNA/PredictorMfeEns.h b/src/IntaRNA/PredictorMfeEns.h index b9c9803..2037d52 100644 --- a/src/IntaRNA/PredictorMfeEns.h +++ b/src/IntaRNA/PredictorMfeEns.h @@ -61,6 +61,9 @@ class PredictorMfeEns : public PredictorMfe { //! map storing the partition of Zall for all considered interaction sites std::unordered_map Z_partition; + //! whether updateZ() currently receives one complete boundary partition + bool updateZisComplete; + /** * Initializes the hybridization partition functions. @@ -94,6 +97,20 @@ class PredictorMfeEns : public PredictorMfe { , const Z_type partFunct , const bool isHybridZ ); + /** + * Consumes a complete partition for one interaction boundary. Exact 2D + * prediction finalizes each boundary once and therefore does not need to + * retain it until reportOptima(). The virtual updateZ() hook is called so + * subclasses can observe or alter the update; overrides that delegate to the + * base implementation inherit streaming. Tracker-enabled predictions retain + * the partition so their established callback timing and ordering stay intact. + */ + void + updateCompleteZ( const size_t i1, const size_t j1 + , const size_t i2, const size_t j2 + , const Z_type partFunct + , const bool isHybridZ ); + /** * Calls for the stored Z_partition information updateOptima() before * calling reportOptima() from its super class. @@ -113,7 +130,20 @@ class PredictorMfeEns : public PredictorMfe { private: /** - * Calls updateOptima() for each entry of Z_partition. + * Applies output filters, converts the partition representation and adds + * the accepted contribution to Zall. Returns the partition without ED. + */ + bool + addPartitionContribution( const size_t i1, const size_t j1 + , const size_t i2, const size_t j2 + , const Z_type partFunct + , const bool isHybridZ + , Z_type & partZ_noED ); + + /** + * Calls updateOptima() for each buffered entry of Z_partition. This hook is + * still called by reportOptima(), but complete tracker-free partitions were + * already consumed by updateZ() and are therefore not present in the map. */ virtual void diff --git a/src/IntaRNA/PredictorMfeEns2d.cpp b/src/IntaRNA/PredictorMfeEns2d.cpp index 145d355..03820f3 100644 --- a/src/IntaRNA/PredictorMfeEns2d.cpp +++ b/src/IntaRNA/PredictorMfeEns2d.cpp @@ -194,7 +194,7 @@ fillHybridZ( const size_t j1, const size_t j2 // update mfe if needed if (callUpdateZ) { - updateZ(i1, j1, i2, j2, curZ, true); + updateCompleteZ(i1, j1, i2, j2, curZ, true); } } // complementary base pair diff --git a/src/IntaRNA/PredictorMfeEns2d.h b/src/IntaRNA/PredictorMfeEns2d.h index 12c4c74..df845b5 100644 --- a/src/IntaRNA/PredictorMfeEns2d.h +++ b/src/IntaRNA/PredictorMfeEns2d.h @@ -72,13 +72,14 @@ class PredictorMfeEns2d: public PredictorMfeEns { /** * Computes all entries of the hybridE matrix for interactions ending in - * p=j1 and q=j2 and report all valid interactions to updateZ() + * p=j1 and q=j2 and reports every complete boundary via updateCompleteZ(), + * which invokes the virtual updateZ() hook. * * @param j1 end of the interaction within seq 1 * @param j2 end of the interaction within seq 2 * @param i1init smallest value for i1 * @param i2init smallest value for i2 - * @param callUpdateZ whether or not updateZ() is to be called + * @param callUpdateZ whether or not updateCompleteZ() is to be called * */ virtual diff --git a/tests/PredictorMfeEnsRegression_test.cpp b/tests/PredictorMfeEnsRegression_test.cpp index b70326a..e4c5512 100644 --- a/tests/PredictorMfeEnsRegression_test.cpp +++ b/tests/PredictorMfeEnsRegression_test.cpp @@ -14,6 +14,7 @@ #include "IntaRNA/SeedHandlerNoBulge.h" #include +#include using namespace IntaRNA; @@ -32,6 +33,79 @@ class InspectableEnsemblePredictor : public PredictorType { } }; +class InspectableExactEnsemblePredictor : public PredictorMfeEns2d { +public: + InspectableExactEnsemblePredictor( + const InteractionEnergy & energy, + OutputHandler & output, + PredictionTracker * tracker) + : PredictorMfeEns2d(energy, output, tracker) + {} + + size_t getPartitionCount() const { + return Z_partition.size(); + } +}; + +class CountingPredictionTracker : public PredictionTracker { +public: + explicit CountingPredictionTracker(size_t & calls) + : calls(calls) + {} + + void updateOptimumCalled(const size_t, const size_t, + const size_t, const size_t, const E_type) override + { + ++calls; + } + +private: + size_t & calls; +}; + +class InterceptingExactEnsemblePredictor : public PredictorMfeEns2d { +public: + InterceptingExactEnsemblePredictor( + const InteractionEnergy & energy, + OutputHandler & output) + : PredictorMfeEns2d(energy, output, NULL) + , updateCalls(0) + , throwOnUpdate(false) + {} + + size_t getPartitionCount() const { + return Z_partition.size(); + } + + void initializeForUpdateTest() { + initOptima(); + initZ(); + } + + void addCompleteForUpdateTest() { + updateCompleteZ(0, 0, 0, 0, Z_type(1), true); + } + + void addBufferedForUpdateTest() { + PredictorMfeEns::updateZ(0, 0, 0, 0, Z_type(1), true); + } + + size_t updateCalls; + bool throwOnUpdate; + +protected: + void updateZ(const size_t i1, const size_t j1, + const size_t i2, const size_t j2, + const Z_type partZ, const bool isHybridZ) override + { + ++updateCalls; + if (throwOnUpdate) { + throw std::runtime_error("intercepted complete partition"); + } + PredictorMfeEns::updateZ(i1, j1, i2, j2, partZ, isHybridZ); + } +}; + TEST_CASE("ensemble predictor regressions", "[PredictorMfeEns]") { #include "testEasyLoggingSetup.icc" @@ -135,4 +209,88 @@ TEST_CASE("ensemble predictor regressions", "[PredictorMfeEns]") { REQUIRE(heuristic.getPartitionCount() == 0); REQUIRE(heuristic.getZall() == 0.0); } + + SECTION("exact complete boundaries stream unless tracker ordering is required") { + RnaSequence target("target", "GGGG"); + RnaSequence query("query", "CCCC"); + AccessibilityDisabled targetAcc(target, 0, NULL); + AccessibilityDisabled queryAcc(query, 0, NULL); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc, 2, 2); + + for (int overlapValue = OutputConstraint::OVERLAP_NONE; + overlapValue <= OutputConstraint::OVERLAP_BOTH; ++overlapValue) + { + OutputConstraint constraint(5, + static_cast(overlapValue), + E_INF, E_INF, false, false, false, true, false); + + OutputHandlerInteractionList streamedOut(constraint, 5); + InspectableExactEnsemblePredictor streamed(energy, streamedOut, NULL); + streamed.predict(); + REQUIRE(streamed.getPartitionCount() == 0); + + size_t trackerCalls = 0; + OutputHandlerInteractionList bufferedOut(constraint, 5); + InspectableExactEnsemblePredictor buffered(energy, bufferedOut, + new CountingPredictionTracker(trackerCalls)); + buffered.predict(); + REQUIRE(buffered.getPartitionCount() > 0); + REQUIRE(trackerCalls == buffered.getPartitionCount()); + + REQUIRE(streamed.getZall() == buffered.getZall()); + REQUIRE(streamedOut.getZ() == bufferedOut.getZ()); + REQUIRE(streamedOut.reported() == bufferedOut.reported()); + + auto streamedInteraction = streamedOut.begin(); + auto bufferedInteraction = bufferedOut.begin(); + for (; streamedInteraction != streamedOut.end() + && bufferedInteraction != bufferedOut.end(); + ++streamedInteraction, ++bufferedInteraction) + { + REQUIRE(**streamedInteraction == **bufferedInteraction); + } + REQUIRE(streamedInteraction == streamedOut.end()); + REQUIRE(bufferedInteraction == bufferedOut.end()); + } + } + + SECTION("exact streaming retains the virtual update hook") { + RnaSequence target("target", "GGGG"); + RnaSequence query("query", "CCCC"); + AccessibilityDisabled targetAcc(target, 0, NULL); + AccessibilityDisabled queryAcc(query, 0, NULL); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc, 2, 2); + OutputConstraint constraint(5, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, false, false, true, false); + OutputHandlerInteractionList output(constraint, 5); + InterceptingExactEnsemblePredictor predictor(energy, output); + + predictor.predict(); + + REQUIRE(predictor.updateCalls > 0); + REQUIRE(predictor.getPartitionCount() == 0); + } + + SECTION("complete update mode is restored after an override throws") { + RnaSequence target("target", "G"); + RnaSequence query("query", "C"); + AccessibilityDisabled targetAcc(target, 0, NULL); + AccessibilityDisabled queryAcc(query, 0, NULL); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc); + OutputConstraint constraint(1, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, false, false, true, false); + OutputHandlerInteractionList output(constraint, 1); + InterceptingExactEnsemblePredictor predictor(energy, output); + predictor.initializeForUpdateTest(); + predictor.throwOnUpdate = true; + + REQUIRE_THROWS_AS(predictor.addCompleteForUpdateTest(), std::runtime_error); + + predictor.throwOnUpdate = false; + predictor.addBufferedForUpdateTest(); + REQUIRE(predictor.getPartitionCount() == 1); + } }