diff --git a/ChangeLog b/ChangeLog index f40c8b25..ae13208a 100644 --- a/ChangeLog +++ b/ChangeLog @@ -71,6 +71,9 @@ - align probability API invalid-input tests with debug-only precondition checks and avoid an unchecked null seed-mask call in release test runs (PR #258) +- finalize completed legacy seed-extension partitions on the probability + branch and release their boundary storage (issue #260, forward port of PR #261) + - use hashed base-pair keys for sparse stack-seed lookup, assert internal seed geometry in debug builds, and document probability computation steps and contributor conventions (PR #258 review) @@ -200,6 +203,13 @@ energies, restricted partition sums, and trackers. * exercise debug-only region/mask/handler/length preconditions only in debug builds + restore the unsupported-handler test under the matching debug guard * keep ordinary lifecycle, arithmetic and output checks in both configurations + * PredictorMfeEns, PredictorMfeEns2dSeedExtension, + PredictorMfeEns2dHeuristicSeedExtension, PredictorMfeEnsRegression_test.cpp + + recover the validation-only right-boundary frontier and buffered comparisons + * preserve seed order, probability tracing, exact kernels and tracker callbacks + * include the containers and hash helper directly in the public header + * doc/probability-memory-frontier.md, doc/Makefile.am + + document the recovered port, original provenance and coverage 261013 Alexander Mitrofanov * StackSeedDomain, SeededPartitionFunction, PredictorMfeEns2dSeedExtension diff --git a/doc/Makefile.am b/doc/Makefile.am index 9d47b6ab..301b4729 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -4,6 +4,7 @@ ################################################################ EXTRA_DIST = \ + probability-memory-frontier.md \ conda.txt \ analysis/base-pair-probabilities-plan.md \ analysis/seeded-kernel-comparison.md \ diff --git a/doc/probability-memory-frontier.md b/doc/probability-memory-frontier.md new file mode 100644 index 00000000..84d24f95 --- /dev/null +++ b/doc/probability-memory-frontier.md @@ -0,0 +1,44 @@ +# Seed-extension memory frontier on the probability branch + +This is the forward port of [PR #261](https://github.com/BackofenLab/IntaRNA/pull/261) +to the base-pair probability branch in PR #258. It groups partial ensemble +partitions by their target right boundary. Once no remaining seed starts +at or before that boundary, the predictor can finalize the site's energy +and release the bucket. The key stores the other three coordinates. + +The seed list retains the seed handler's traversal order and records a suffix +minimum of target starts. This keeps explicit seeds safe even when their +iteration order is not monotonic. Partition contributions and Zall keep +their existing addition order. Tracker-enabled runs retain the deferred map +and callback order. Each prediction resets the frontier, including empty +regions and predictor reuse. + +The change applies to the legacy exact seeded extension path and the +heuristic seeded extension path, including its optional probability trace. +The canonical exact stack-seed kernel already emits complete boundaries; +its recurrence and reverse probability calculation are unchanged. Region +bounds and explicit/singleton seed admission remain those of PR #258. + +## Why the code remained uncommitted + +The recovered implementation lived in `issue-260-validation/feat-fixed` as +an integration-validation port. PR #261 intentionally targeted master and +included this port as `doc/ensemble-memory-probability-port.patch` rather +than committing it to the probability branch. The old validation worktree +also contained SVG changes that were subsequently committed in PR #259; +those changes are already present in this PR's base and are not duplicated. +This PR promotes the six-file memory port to a reviewable commit and adds +the direct public-header includes used by its containers and hash helper. + +The saved October 12 API run passed 101,247 assertions in 97 cases. That is +historical evidence; current validation results are recorded in this PR's +description. Historical measurements and limitations remain in PR #261. +This port makes no new universal memory or runtime improvement claim. + +## Regression coverage + +The recovered tests compare streamed predictions with a buffered subclass +for partition totals and ranked interaction sites. They cover exact and +heuristic paths, stack/bulged/explicit seeds, LP/noLP, cropped regions, +reused predictors and tracker callbacks. The probability CLI suites also +check that CSV/SVG output still agrees with the supported ensembles. diff --git a/src/IntaRNA/PredictorMfeEns.cpp b/src/IntaRNA/PredictorMfeEns.cpp index 3ed283ea..6e077132 100644 --- a/src/IntaRNA/PredictorMfeEns.cpp +++ b/src/IntaRNA/PredictorMfeEns.cpp @@ -34,6 +34,8 @@ initZ() { // reset storage Z_partition.clear(); + Z_frontier.clear(); + Z_frontierNext = 0; } //////////////////////////////////////////////////////////////////////////// @@ -120,12 +122,12 @@ updateZ( const size_t i1, const size_t j1 } // store partial Z (without ED) - Interaction::Boundary key(i1,j1,i2,j2); - auto [keyEntry, inserted] = Z_partition.try_emplace(key, partZ_noED); - if (!inserted) { - // update entry - keyEntry->second += partZ_noED; - } + auto accumulate = [partZ_noED](auto & partitions,const auto & key) { + auto [keyEntry,inserted] = partitions.try_emplace(key,partZ_noED); + if (!inserted) keyEntry->second += partZ_noED; + }; + if (Z_frontier.empty()) accumulate(Z_partition,Interaction::Boundary(i1,j1,i2,j2)); + else accumulate(Z_frontier.at(j1),FrontierBoundary{i1,i2,j2}); } @@ -180,6 +182,7 @@ void PredictorMfeEns:: updateOptimaUsingZ() { + advanceZFrontier(Z_frontier.size()); for (auto it = Z_partition.begin(); it != Z_partition.end(); ++it) { // if partition function is > 0 @@ -189,6 +192,28 @@ updateOptimaUsingZ() } } +void PredictorMfeEns::beginZFrontier(size_t size1) +{ + if (predTracker == nullptr) Z_frontier.resize(size1); +} + +void PredictorMfeEns::advanceZFrontier(size_t minimumSeedStart) +{ + const size_t end = std::min(minimumSeedStart,Z_frontier.size()); + while (Z_frontierNext < end) { + auto & partitions = Z_frontier[Z_frontierNext]; + for (const auto & [boundary,partition] : partitions) { + if (Z_isNotINF(partition) && partition > 0) { + PredictorMfe::updateOptima(boundary.i1,Z_frontierNext,boundary.i2,boundary.j2, + energy.getE(partition),true,false); + } + } + // Free both nodes and bucket storage as soon as this boundary is complete. + decltype(Z_frontier)::value_type{}.swap(partitions); + ++Z_frontierNext; + } +} + //////////////////////////////////////////////////////////////////////////// void diff --git a/src/IntaRNA/PredictorMfeEns.h b/src/IntaRNA/PredictorMfeEns.h index c51fd654..27d352f0 100644 --- a/src/IntaRNA/PredictorMfeEns.h +++ b/src/IntaRNA/PredictorMfeEns.h @@ -8,7 +8,10 @@ #include "IntaRNA/IndexRangeList.h" #include +#include #include +#include +#include namespace IntaRNA { @@ -75,6 +78,21 @@ class PredictorMfeEns : public PredictorMfe { bool updateZisComplete; + //! a bucket supplies j1, so its key needs only the other three coordinates + struct FrontierBoundary { + size_t i1, i2, j2; + bool operator==(const FrontierBoundary &) const = default; + struct Hash { + size_t operator()(const FrontierBoundary & boundary) const; + }; + }; + + //! partial sums grouped by target right boundary; used only without trackers + std::vector> Z_frontier; + //! first right-boundary bucket that has not yet been finalized + size_t Z_frontierNext = 0; + + /** * Initializes the hybridization partition functions. * Will be called by predict(). @@ -121,6 +139,18 @@ class PredictorMfeEns : public PredictorMfe { , const Z_type partFunct , const bool isHybridZ ); + /** Enable partial-partition buckets by target right boundary. Tracker runs + * keep the original map to preserve callback ordering. Call advanceZFrontier + * only when no remaining seed can start below the supplied target index. + * @param size1 number of target positions in the current prediction region + */ + void beginZFrontier(size_t size1); + + /** Consume completed boundary sums without adding their mass to Zall again. + * @param minimumSeedStart smallest target start of any remaining seed + */ + void advanceZFrontier(size_t minimumSeedStart); + /** * Calls for the stored Z_partition information updateOptima() before * calling reportOptima() from its super class. @@ -161,6 +191,17 @@ class PredictorMfeEns : public PredictorMfe { }; +inline +size_t +PredictorMfeEns::FrontierBoundary::Hash::operator()(const FrontierBoundary & boundary) const +{ + size_t key=0; + boost::hash_combine(key,boundary.i1); + boost::hash_combine(key,boundary.i2); + boost::hash_combine(key,boundary.j2); + return key; +} + } // namespace #endif /* INTARNA_PREDICTORMFEENS_H_ */ diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp index 278eb740..527d4935 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp @@ -85,11 +85,10 @@ predict( const IndexRange & r1, const IndexRange & r2 ) // initialize overall partition function for updates initZ(); - size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; - while( seedHandler.updateToNextSeed(si1,si2 - , 0, range_size1-1 - , 0, range_size2-1) ) - { + beginZFrontier(range_size1); + for (const auto & seed : getSeedStarts(range_size1-1,range_size2-1)) { + advanceZFrontier(seed.minimumRemaining); + const size_t si1=seed.i1, si2=seed.i2; const Z_type seedZ = energy.getBoltzmannWeight( seedHandler.getSeedE(si1, si2) ); const size_t sl1 = seedHandler.getSeedLength1(si1, si2); const size_t sl2 = seedHandler.getSeedLength2(si1, si2); diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp index 0b6eb098..274ce01b 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp @@ -7,6 +7,27 @@ namespace IntaRNA { + +std::vector +PredictorMfeEns2dSeedExtension::getSeedStarts(size_t max1,size_t max2) const +{ + std::vector starts; + size_t i1=RnaSequence::lastPos, i2=RnaSequence::lastPos; + while (seedHandler.updateToNextSeed(i1,i2, + 0,max1,0,max2)) { + starts.push_back({i1,i2,i1}); + } + // A suffix minimum is safe for either traversal order and custom handlers. + size_t minimum=max1; + for (auto seed=starts.rbegin(); seed!=starts.rend(); ++seed) { + minimum=std::min(minimum,seed->i1); + seed->minimumRemaining=minimum; + } + return starts; +} + + + // The heuristic objective contains signed seed-overlap corrections. Trace the // implemented arithmetic directly, retaining pair ownership on every term. // This is allocated only for requested probability output; ordinary prediction @@ -268,11 +289,11 @@ predict( const IndexRange & r1, const IndexRange & r2 ) // initialize overall partition function for updates initZ(); - size_t si1 = RnaSequence::lastPos, si2 = RnaSequence::lastPos; - while( seedHandler.updateToNextSeed(si1,si2 - , 0, range_size1+1-seedHandler.getConstraint().getBasePairs() - , 0, range_size2+1-seedHandler.getConstraint().getBasePairs()) ) - { + beginZFrontier(range_size1); + for (const auto & seed : getSeedStarts(range_size1+1-seedHandler.getConstraint().getBasePairs(), + range_size2+1-seedHandler.getConstraint().getBasePairs())) { + advanceZFrontier(seed.minimumRemaining); + const size_t si1=seed.i1, si2=seed.i2; // get Z and boundaries of seed const Z_type seedZ = energy.getBoltzmannWeight( seedHandler.getSeedE(si1, si2) ); const size_t sl1 = seedHandler.getSeedLength1(si1, si2); diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.h b/src/IntaRNA/PredictorMfeEns2dSeedExtension.h index 0b344e4a..7bb46c1e 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.h +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.h @@ -25,6 +25,16 @@ namespace IntaRNA { class PredictorMfeEns2dSeedExtension: public PredictorMfeEns { protected: + //! seed coordinates in handler order and the minimum remaining target start + struct SeedStart { size_t i1, i2, minimumRemaining; }; + + /** Collect seed starts without changing traversal or partition summation order. + * @param max1 largest target seed start + * @param max2 largest reversed-query seed start + * @return starts with a safe lower bound for all remaining target seed starts + */ + std::vector getSeedStarts(size_t max1,size_t max2) const; + //! matrix type to hold the partition functions for interaction site starts typedef Matrix Z2dMatrix; diff --git a/tests/PredictorMfeEnsRegression_test.cpp b/tests/PredictorMfeEnsRegression_test.cpp index e4f55379..f1a20e3a 100644 --- a/tests/PredictorMfeEnsRegression_test.cpp +++ b/tests/PredictorMfeEnsRegression_test.cpp @@ -13,6 +13,8 @@ #include "IntaRNA/ReverseAccessibility.h" #include "IntaRNA/RnaSequence.h" #include "IntaRNA/SeedHandlerNoBulge.h" +#include "IntaRNA/SeedHandlerMfe.h" +#include "IntaRNA/SeedHandlerExplicit.h" #include #include @@ -26,12 +28,14 @@ class InspectableEnsemblePredictor : public PredictorType { InspectableEnsemblePredictor( const InteractionEnergy & energy, OutputHandler & output, - SeedHandler * seedHandler) - : PredictorType(energy, output, NULL, seedHandler) + SeedHandler * seedHandler, PredictionTracker * tracker=nullptr) + : PredictorType(energy, output, tracker, seedHandler) {} size_t getPartitionCount() const { - return this->Z_partition.size(); + size_t count=this->Z_partition.size(); + for (const auto & bucket : this->Z_frontier) count+=bucket.size(); + return count; } }; @@ -231,7 +235,7 @@ TEST_CASE("ensemble predictor regressions", "[PredictorMfeEns]") { InspectableEnsemblePredictor exact( energy, exactOut, new SeedHandlerNoBulge(energy, seedConstraint)); exact.predict(); - // The new exact stack backend streams each completed boundary once. + // The exact stack backend streams each completed boundary once. REQUIRE(exact.getPartitionCount() == 0); REQUIRE(exact.getZall() > 0); exact.predict(IndexRange(0, 0), IndexRange(0, 0)); @@ -242,12 +246,62 @@ TEST_CASE("ensemble predictor regressions", "[PredictorMfeEns]") { InspectableEnsemblePredictor heuristic( energy, heuristicOut, new SeedHandlerNoBulge(energy, seedConstraint)); heuristic.predict(); - REQUIRE(heuristic.getPartitionCount() > 0); + REQUIRE(heuristic.getPartitionCount() == 0); heuristic.predict(IndexRange(0, 0), IndexRange(0, 0)); REQUIRE(heuristic.getPartitionCount() == 0); REQUIRE(heuristic.getZall() == 0.0); } + SECTION("seed frontier preserves partial sums, ranked sites and filtered regions") { + RnaSequence target("target", "GGAGGGGAGGGGAGGGG"); + RnaSequence query("query", "CCCCACCCC"); + AccessibilityDisabled targetAcc(target, 8, nullptr), queryAcc(query, 8, nullptr); + ReverseAccessibility reversed(queryAcc); + InteractionEnergyBasePair energy(targetAcc,reversed); + for (bool noLP : {false,true}) for (bool heuristic : {false,true}) for (int handler : {0,1,2}) { + const size_t bulges=handler==1?1:0; + SeedConstraint seeds(2,bulges,bulges,bulges,E_INF,Accessibility::ED_UPPER_BOUND,E_INF, + IndexRangeList(),IndexRangeList(),handler==2?"4||&2||,10||&6||":"",false,false,noLP); + auto createSeed=[&]() -> SeedHandler * { + if (handler==1) return new SeedHandlerMfe(energy,seeds); + if (handler==2) return new SeedHandlerExplicit(energy,seeds); + return new SeedHandlerNoBulge(energy,seeds); + }; + OutputConstraint constraint(5,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF, + false,noLP,false,true,false); + OutputHandlerInteractionList streamedOut(constraint,5), bufferedOut(constraint,5); + size_t callbacks=0; + auto check=[&](auto & streamed,auto & buffered) { + for (const auto & region : {IndexRange(0,16),IndexRange(3,14)}) { + streamed.predict(region,IndexRange(1,7)); + buffered.predict(region,IndexRange(1,7)); + REQUIRE(streamed.getPartitionCount()==0); + REQUIRE(buffered.getPartitionCount()>0); + REQUIRE(streamed.getZall()==buffered.getZall()); + REQUIRE(streamedOut.getZ()==bufferedOut.getZ()); + REQUIRE(streamedOut.reported()==bufferedOut.reported()); + auto a=streamedOut.begin(), b=bufferedOut.begin(); + for (; a!=streamedOut.end() && b!=bufferedOut.end(); ++a,++b) REQUIRE(**a==**b); + REQUIRE(a==streamedOut.end()); REQUIRE(b==bufferedOut.end()); + } + REQUIRE(callbacks>0); + }; + if (heuristic) { + InspectableEnsemblePredictor streamed( + energy,streamedOut,createSeed()); + InspectableEnsemblePredictor buffered( + energy,bufferedOut,createSeed(),new CountingPredictionTracker(callbacks)); + check(streamed,buffered); + } else { + InspectableEnsemblePredictor streamed( + energy,streamedOut,createSeed()); + InspectableEnsemblePredictor buffered( + energy,bufferedOut,createSeed(),new CountingPredictionTracker(callbacks)); + check(streamed,buffered); + } + } + } + SECTION("exact complete boundaries stream unless tracker ordering is required") { RnaSequence target("target", "GGGG"); RnaSequence query("query", "CCCC");