Commit the recovered probability-branch ensemble memory frontier

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

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 \

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.

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

src/IntaRNA/PredictorMfeEns.h

@@ -8,7 +8,10 @@
 #include "IntaRNA/IndexRangeList.h"
 
 #include <list>
+#include <unordered_map>
 #include <utility>
+#include <vector>
+#include <boost/functional/hash.hpp>
 
 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<std::unordered_map<FrontierBoundary,Z_type,FrontierBoundary::Hash>> 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_ */

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);

src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp

@@ -7,6 +7,27 @@
 
 namespace IntaRNA {
 
+
+std::vector<PredictorMfeEns2dSeedExtension::SeedStart>
+PredictorMfeEns2dSeedExtension::getSeedStarts(size_t max1,size_t max2) const
+{
+	std::vector<SeedStart> 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);

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<SeedStart> getSeedStarts(size_t max1,size_t max2) const;
+
 
 	//! matrix type to hold the partition functions for interaction site starts
 	typedef Matrix<Z_type> Z2dMatrix;

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 <cmath>
 #include <iterator>
@@ -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<PredictorMfeEns2dSeedExtension> 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<PredictorMfeEns2dHeuristicSeedExtension> 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<PredictorMfeEns2dHeuristicSeedExtension> streamed(
+						energy,streamedOut,createSeed());
+				InspectableEnsemblePredictor<PredictorMfeEns2dHeuristicSeedExtension> buffered(
+						energy,bufferedOut,createSeed(),new CountingPredictionTracker(callbacks));
+				check(streamed,buffered);
+			} else {
+				InspectableEnsemblePredictor<PredictorMfeEns2dSeedExtension> streamed(
+						energy,streamedOut,createSeed());
+				InspectableEnsemblePredictor<PredictorMfeEns2dSeedExtension> 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");