From 0b20908db9ec1a3d714dc1cfc0f18adbe7cd2bb1 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:09:50 +0200 Subject: [PATCH 01/27] test: document phase 1 audit and expose baseline defects Record the 3.4.1 architecture, numeric contracts, bottlenecks, test gaps, and independently reproduced correctness findings for discussion #232. Add focused regressions before fixes: Interaction self-assignment and null-seed equality, zero-capacity result storage, invalid Nussinov paired intervals, base-pair ES semantics, noLP/terminal-GU partition counting, and target/query accessibility option isolation. These regressions intentionally fail on the audited baseline so each correction can follow independently. --- doc/refactor/1-current-state.md | 150 ++++++++++++++++++ tests/InteractionEnergyBasePair_test.cpp | 6 +- tests/Interaction_test.cpp | 35 ++++ tests/Makefile.am | 4 +- tests/NussinovHandler_test.cpp | 11 ++ tests/OutputHandlerInteractionList_test.cpp | 13 ++ tests/PredictorMfeEnsRegression_test.cpp | 58 +++++++ .../data/energyB-target-acc-options.parameter | 17 ++ .../energyB-target-acc-options.testresult | 8 + 9 files changed, 297 insertions(+), 5 deletions(-) create mode 100644 doc/refactor/1-current-state.md create mode 100644 tests/PredictorMfeEnsRegression_test.cpp create mode 100644 tests/data/energyB-target-acc-options.parameter create mode 100644 tests/data/energyB-target-acc-options.testresult diff --git a/doc/refactor/1-current-state.md b/doc/refactor/1-current-state.md new file mode 100644 index 00000000..a63dd8df --- /dev/null +++ b/doc/refactor/1-current-state.md @@ -0,0 +1,150 @@ +# Phase 1: current state + +This document records the phase-1 audit requested in discussion #232. It is +about scientific behavior and runtime cost, not a style rewrite. The audited +baseline is IntaRNA 3.4.1, commit `3b14bc0`. + +## Reproducible baseline + +The source was built as an optimized OpenMP build with GCC 13.3, Boost 1.85 and +ViennaRNA 2.7.2. Before any source change: + +- the API binary passed 3,859 assertions in 30 test cases; +- all 17 command-line golden cases passed; and +- the golden cases covered only the table-independent base-pair energy model + with disabled accessibility. + +Thus the old suite is a useful output-stability gate, but it is not broad enough +to establish correctness of the ensemble predictors, ViennaRNA paths, +parallel/window reduction, or several public boundary cases. + +## Architecture and data flow + +`src/bin/IntaRNA.cpp` owns one invocation. `CommandLineParsing` parses the +selected personality and works as a factory for sequence, accessibility, +energy, seed/helix, predictor, tracker and output objects. Query accessibility +is cached across targets. For each target/query pair the program computes or +loads accessibility, constructs an interaction-energy facade, decomposes +allowed ranges into optional overlapping windows, runs one predictor per +window, and merges interactions in an `OutputHandlerInteractionList`. + +The library keeps the biological pieces modular: + +| Family | Responsibility | Scientific invariant | +| --- | --- | --- | +| `RnaSequence`, `IndexRange*`, `Interaction*` | sequence and inclusive coordinate value types | external/internal and reversed-query coordinates remain bijective and in range | +| `Accessibility*` | opening energy ED for a sequence interval | `ED=-RT log(Pu)`; unavailable intervals return the documented upper bound | +| `InteractionEnergy*` | initiation, loops, ends, dangles, ED and Boltzmann conversion | total energy uses integer centi-kcal/mol without arithmetic on infinity sentinels | +| `SeedHandler*`, `HelixHandler*` | admissible seeds/helices and traceback | preprocessing bounds and traceback describe the same structure | +| `PredictorMfe*` | exact/heuristic MFE dynamic programs | every reported boundary has a valid recurrence and traceback | +| `PredictorMfeEns*` | interaction partition functions and ensemble representatives | every valid interaction contributes exactly once with multiplicative independent weights | +| `OutputHandler*`, `PredictionTracker*` | filtering, reduction and rendering | filters apply consistently to both reported sites and their requested ensemble | + +The second sequence is exposed in biological orientation and wrapped by +`ReverseAccessibility` for the antiparallel dynamic programs. This convention +is performance-sensitive: offsets can be cached, but changing the convention +would invalidate almost every recurrence and traceback. + +### Numeric model + +MFE values use `E_type=int` in centi-kcal/mol. `E_INF` and `E_MAX` are finite +sentinels chosen to leave arithmetic headroom. Partition values use `double` +unless the optional quadmath build is enabled. A hybrid path with base pairs +`(i1,i2),...,(j1,j2)` has the schematic energy + +``` +E = E_init + sum(E_interLeft) + ED1(i1,j1) + ED2(i2,j2) + + E_dangleLeft + E_dangleRight + E_endLeft + E_endRight + E_add +``` + +and weight `exp(-E/RT)`. Exact energy ties intentionally use integer equality; +loop order and tie-breaking must remain deterministic in optimized code. + +## Runtime structure and bottlenecks + +The default personality uses ViennaRNA accessibility, a bulge-free seed and +`PredictorMfe2dHeuristicSeedExtension`. Important costs are: + +1. ViennaRNA accessibility preprocessing and storage of banded ED tables. +2. Two-dimensional predictor matrices, repeatedly resized and initialized for + windows, seeds, and right boundaries. +3. Nested internal-loop scans in which virtual offset, complementarity, ED and + energy calls are repeated for nearby cells. +4. Exact ensemble prediction stores a hash entry keyed by all four site + boundaries. That can grow as `O(n1^2*n2^2)` even though the recurrence matrix + itself is two-dimensional. +5. Output/top-k and partition updates enter global OpenMP critical sections; + parallel work is selected at one outer dimension rather than represented as + independent pair/window tasks. +6. Automatic accessibility range decomposition is performed from mutable + parser-owned vectors and can be repeated for sequence pairs. + +The first optimization work should therefore remove unnecessary state and +calls from existing recurrences before changing their mathematics. + +## Confirmed defects and regressions + +The following findings were traced in this source tree. Tests are added before +fixes, as required by phase 1. + +| Finding | Evidence / consequence | Phase-1 regression | +| --- | --- | --- | +| `Interaction::operator=` is not self-assignment safe | it clears `basePairs` before reading the same object | preserve pairs, energy and seed metadata across `x=x` | +| interaction equality dereferences asymmetric null seed pointers | `seed == i.seed || *seed == *i.seed` dereferences when only one side is null | seeded and unseeded interactions compare unequal without crashing | +| zero-capacity interaction storage dereferences an empty reverse iterator | `maxToStore=0` reaches `*storage.rbegin()` | count a report but retain no interaction | +| `NussinovHandler::getQb` returns one for an out-of-range paired interval | an invalid pair obtains the multiplicative identity instead of zero weight | `getQb(i,n)==0` | +| base-pair `getES*` includes the empty monomer structure | it stores `-RT log(Q)`, while the API specifies structures containing at least one pair, i.e. `-RT log(Q-1)` | a four-base sequence with one possible pair has `ES=-1`, not `-log(1+e)` | +| noLP heuristic ensemble counts the direct stack twice | the explicit direct continuation and the `w1=w2=1` loop iteration denote the same path | on `GGGG/CCCC`, heuristic `Zall` must not exceed exact `Zall` | +| ensemble `updateZ` bypasses `noGUend` and `maxED` | its direct `Zall`/boundary-map update does not use the filters in `PredictorMfe::updateOptima` | a single terminal GU contributes zero when terminal GU is forbidden | +| target base-pair accessibility uses query limits | target factory reads `qIntLenMax/qAccW` instead of target parameters | asymmetric CLI target-accessibility table retains the target maximum length | + +Additional high-confidence findings are not treated as tiny local fixes because +they change aggregation or recurrence ownership and need the benchmark gates of +later phases: + +- exact ensemble four-boundary storage violates the advertised practical + two-dimensional memory bound; +- seed-extension ensemble code has early-return state reuse, additive where + independent partition factors must be multiplicative, noLP duplication and + unsigned boundary risks; +- windowed/global partition and tracker aggregation can count overlapping + domains more than once; +- query/target range factories mutate parser-owned vectors from `const` + methods, repeat decomposition and can race across pair tasks; +- partition accumulation can overflow in release builds; debug-only warnings + neither prevent nor repair the result; and +- CLI object ownership uses raw factories and `const_cast` cleanup, making + exceptional and parallel paths difficult to reason about. + +## Missing coverage + +The original suite has no direct test for exact MFE, exact seeded MFE, any +ensemble predictor, seed-extension predictors, zero requested output, or +partition/filter agreement. It also lacks: + +- a brute-force small-instance partition oracle; +- reuse tests for stateful predictors; +- asymmetric target/query option tests; +- multiple-region and overlapping-window aggregation tests; +- deterministic threaded-output tests; +- optimized/release invalid-input tests; +- sanitizer coverage; and +- a GCC 14 plus macOS Clang portability build. + +Phase 1 adds focused deterministic cases for the confirmed local defects. +Phase 2 adds compiler/build coverage. The benchmark and differential corpus +introduced before phase 3 supplies executable-level output parity for every +performance change. + +## Change record + +- Created the required `refactoring` branch from `master` at `3b14bc0`. +- Added this architecture, correctness and performance audit. +- Established the clean 30-case / 3,859-assertion API and 17-case CLI baseline. +- Added table-independent regressions for ownership, boundary, partition and + target/query isolation defects before their fixes. + +The discussion once calls the phase-1 document `refactor-changelog.md`; no such +file exists and the same phase otherwise consistently requires +`doc/refactor/1-current-state.md`. This file is therefore both the current-state +analysis and the reproducibility/change log. diff --git a/tests/InteractionEnergyBasePair_test.cpp b/tests/InteractionEnergyBasePair_test.cpp index 27360332..82e810ae 100644 --- a/tests/InteractionEnergyBasePair_test.cpp +++ b/tests/InteractionEnergyBasePair_test.cpp @@ -72,8 +72,10 @@ TEST_CASE( "InteractionEnergyBasePair", "[InteractionEnergyBasePair]" ) { } SECTION("ES computation") { - REQUIRE( E_equal(energy.getES1(0, 3), Ekcal_2_E(-1.313186)) ); - REQUIRE( E_equal(energy.getES2(0, 3), Ekcal_2_E(-1.313186)) ); + // ES covers structures containing at least one intramolecular base pair. + // ACGU has exactly one admissible pair of weight exp(1), hence ES=-1. + REQUIRE( E_equal(energy.getES1(0, 3), Ekcal_2_E(-1.0)) ); + REQUIRE( E_equal(energy.getES2(0, 3), Ekcal_2_E(-1.0)) ); REQUIRE( E_isINF(energy.getES1(0, 2)) ); REQUIRE( E_isINF(energy.getES1(1, 2)) ); REQUIRE( E_isINF(energy.getES2(0, 2)) ); diff --git a/tests/Interaction_test.cpp b/tests/Interaction_test.cpp index a14ead65..c02cd412 100644 --- a/tests/Interaction_test.cpp +++ b/tests/Interaction_test.cpp @@ -90,4 +90,39 @@ TEST_CASE( "Interaction", "[Interaction]" ) { } + SECTION("self assignment preserves owned state") { + + Interaction inter(r,r); + inter.basePairs.push_back( Interaction::BasePair(0,7) ); + inter.basePairs.push_back( Interaction::BasePair(1,6) ); + inter.energy = Ekcal_2_E(-2.0); + inter.seed = new Interaction::SeedSet(); + inter.seed->insert( Interaction::Seed( + Interaction::BasePair(0,7), Interaction::BasePair(1,6), inter.energy) ); + + inter = inter; + + REQUIRE( inter.basePairs.size() == 2 ); + REQUIRE( inter.basePairs.front() == Interaction::BasePair(0,7) ); + REQUIRE( inter.basePairs.back() == Interaction::BasePair(1,6) ); + REQUIRE( inter.energy == Ekcal_2_E(-2.0) ); + REQUIRE( inter.seed != NULL ); + REQUIRE( inter.seed->size() == 1 ); + } + + SECTION("seeded and unseeded interactions compare safely") { + + Interaction unseeded(r,r); + unseeded.basePairs.push_back( Interaction::BasePair(0,7) ); + unseeded.energy = Ekcal_2_E(-1.0); + + Interaction seeded(unseeded); + seeded.seed = new Interaction::SeedSet(); + seeded.seed->insert( Interaction::Seed( + Interaction::BasePair(0,7), Interaction::BasePair(0,7), seeded.energy) ); + + REQUIRE_FALSE( unseeded == seeded ); + REQUIRE_FALSE( seeded == unseeded ); + } + } diff --git a/tests/Makefile.am b/tests/Makefile.am index ff862502..07d0bcab 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -46,6 +46,7 @@ runApiTests_SOURCES = \ PredictionTrackerSpotProb_test.cpp \ PredictorMfe2dHelixBlockHeuristic_test.cpp \ PredictorMfe2dHelixBlockHeuristicSeed_test.cpp \ + PredictorMfeEnsRegression_test.cpp \ NussinovHandler_test.cpp \ RnaSequence_test.cpp \ OutputStreamHandlerSortedCsv_test.cpp \ @@ -65,6 +66,3 @@ LIBS= -L$(top_builddir)/src/IntaRNA -lIntaRNA \ runApiTests_CXXFLAGS = -I$(top_builddir)/src \ @AM_CXXFLAGS@ @CXXFLAGS@ \ -DELPP_NO_LOG_TO_FILE - - - \ No newline at end of file diff --git a/tests/NussinovHandler_test.cpp b/tests/NussinovHandler_test.cpp index cb345b58..af70f566 100644 --- a/tests/NussinovHandler_test.cpp +++ b/tests/NussinovHandler_test.cpp @@ -68,4 +68,15 @@ TEST_CASE("NussinovHandler", "[NussinovHandler]") { REQUIRE(nuss == 3u); } + + SECTION("Out-of-range paired intervals have zero weight") { + + RnaSequence rna("test", "ACGU"); + NussinovHandler::Z2dMatrix Q(rna.size(), rna.size()); + NussinovHandler::Z2dMatrix Qb(rna.size(), rna.size()); + Q.clear(); + Qb.clear(); + + REQUIRE(NussinovHandler::getQb(0, rna.size(), rna, 1.0, 0, Q, Qb) == 0.0); + } } diff --git a/tests/OutputHandlerInteractionList_test.cpp b/tests/OutputHandlerInteractionList_test.cpp index 52da59b8..0ee9674c 100644 --- a/tests/OutputHandlerInteractionList_test.cpp +++ b/tests/OutputHandlerInteractionList_test.cpp @@ -88,4 +88,17 @@ TEST_CASE( "OutputHandlerInteractionList", "[OutputHandlerInteractionList]" ) { } + SECTION("zero storage capacity") { + + Interaction i(r,r); + i.basePairs.push_back( Interaction::BasePair(0,7) ); + i.energy = Ekcal_2_E(-1.0); + + OutputHandlerInteractionList out(oc,0); + out.add(i); + + REQUIRE( out.reported() == 1 ); + REQUIRE( out.empty() ); + } + } diff --git a/tests/PredictorMfeEnsRegression_test.cpp b/tests/PredictorMfeEnsRegression_test.cpp new file mode 100644 index 00000000..600b4498 --- /dev/null +++ b/tests/PredictorMfeEnsRegression_test.cpp @@ -0,0 +1,58 @@ +#include "catch.hpp" + +#undef NDEBUG + +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/OutputHandlerInteractionList.h" +#include "IntaRNA/PredictorMfeEns2d.h" +#include "IntaRNA/PredictorMfeEns2dHeuristic.h" +#include "IntaRNA/ReverseAccessibility.h" +#include "IntaRNA/RnaSequence.h" + +using namespace IntaRNA; + +TEST_CASE("ensemble predictor regressions", "[PredictorMfeEns]") { + + #include "testEasyLoggingSetup.icc" + + SECTION("a pruning heuristic never exceeds the exact noLP partition") { + 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); + + OutputConstraint constraint(1, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, true, false, true, false); + OutputHandlerInteractionList exactOut(constraint, 1); + OutputHandlerInteractionList heuristicOut(constraint, 1); + PredictorMfeEns2d exact(energy, exactOut, NULL); + PredictorMfeEns2dHeuristic heuristic(energy, heuristicOut, NULL); + + exact.predict(); + heuristic.predict(); + + REQUIRE(exact.getZall() > 0.0); + REQUIRE(heuristic.getZall() <= exact.getZall() * (1.0 + 1e-12)); + } + + SECTION("terminal GU filtering also applies to the partition") { + RnaSequence target("target", "G"); + RnaSequence query("query", "U"); + 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, true, true, false); + OutputHandlerInteractionList out(constraint, 1); + PredictorMfeEns2d predictor(energy, out, NULL); + + predictor.predict(); + + REQUIRE(predictor.getZall() == 0.0); + } +} diff --git a/tests/data/energyB-target-acc-options.parameter b/tests/data/energyB-target-acc-options.parameter new file mode 100644 index 00000000..ebd84ebe --- /dev/null +++ b/tests/data/energyB-target-acc-options.parameter @@ -0,0 +1,17 @@ +mode=M +noSeed=true +model=S +energy=B +tAcc=C +qAcc=N +tAccW=4 +tAccL=4 +tIntLenMax=4 +qIntLenMax=3 +target=ACGU +query=GGGG +out=STDERR +out=tAcc:STDOUT +outMode=C +outNumber=1 +threads=1 diff --git a/tests/data/energyB-target-acc-options.testresult b/tests/data/energyB-target-acc-options.testresult new file mode 100644 index 00000000..b17376e3 --- /dev/null +++ b/tests/data/energyB-target-acc-options.testresult @@ -0,0 +1,8 @@ +id1;start1;end1;id2;start2;end2;subseqDP;hybridDP;E +target;2;4;query;1;2;CGU&GG;(.(&));-2 +#ensemble delta energy to unpair a region ED + #i$ l=1 2 3 4 +1 0.000000e+00 NA NA NA +2 0.000000e+00 0.000000e+00 NA NA +3 0.000000e+00 0.000000e+00 0.000000e+00 NA +4 0.000000e+00 0.000000e+00 0.000000e+00 1.313186e+00 From 11ac391e6113f1a7dce6ea0e981e9b35852d3cd2 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:10:39 +0200 Subject: [PATCH 02/27] fix: preserve interactions on self-assignment Return before clearing owned pairing and seed state when the source and destination are the same Interaction. This makes the assignment contract safe without changing normal copy behavior. --- src/IntaRNA/Interaction.cpp | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/IntaRNA/Interaction.cpp b/src/IntaRNA/Interaction.cpp index f043cd48..bce4dcad 100644 --- a/src/IntaRNA/Interaction.cpp +++ b/src/IntaRNA/Interaction.cpp @@ -151,6 +151,9 @@ Interaction & Interaction:: operator= ( const Interaction & toCopy ) { + if (this == &toCopy) { + return *this; + } #if INTARNA_IN_DEBUG_MODE if (!toCopy.isValid()) throw std::runtime_error("Interaction::=("+toString(toCopy)+") not valid!"); From f498659890c3df3273354fe4f4acfa114376fb11 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:11:15 +0200 Subject: [PATCH 03/27] fix: compare optional interaction seeds safely Treat identical seed pointers as equal and compare seed sets only when both interactions own one. Seeded and unseeded interactions now compare unequal instead of dereferencing null. --- src/IntaRNA/Interaction.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/IntaRNA/Interaction.cpp b/src/IntaRNA/Interaction.cpp index bce4dcad..87f55695 100644 --- a/src/IntaRNA/Interaction.cpp +++ b/src/IntaRNA/Interaction.cpp @@ -230,7 +230,8 @@ operator == ( const Interaction &i ) const && s2 == i.s2 && E_equal( energy, i.energy ) && basePairs == i.basePairs - && (seed == i.seed || *seed == *(i.seed)) + && (seed == i.seed + || (seed != NULL && i.seed != NULL && *seed == *(i.seed))) ; } From 8db8b0c6d383af82269bbefa6c85fc05176df3b2 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:11:40 +0200 Subject: [PATCH 04/27] fix: handle zero-capacity interaction storage Continue counting reported candidates when maxToStore is zero, but skip list comparison and insertion. This avoids dereferencing an empty reverse iterator for --outNumber=0 style consumers. --- src/IntaRNA/OutputHandlerInteractionList.cpp | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/src/IntaRNA/OutputHandlerInteractionList.cpp b/src/IntaRNA/OutputHandlerInteractionList.cpp index 0fed214c..034d5652 100644 --- a/src/IntaRNA/OutputHandlerInteractionList.cpp +++ b/src/IntaRNA/OutputHandlerInteractionList.cpp @@ -52,7 +52,9 @@ add( const Interaction & interaction ) { // count interaction reportedInteractions++; - if (storage.size() < maxToStore || lessThan_StorageContainer( &interaction, *(storage.rbegin()) )) { + if (maxToStore > 0 + && (storage.size() < maxToStore + || lessThan_StorageContainer( &interaction, *(storage.rbegin()) ))) { // find where to insert this interaction StorageContainer::iterator insertPos = std::lower_bound( storage.begin(), storage.end(), &interaction, lessThan_StorageContainer ); // check if interaction is NOT already part of the list From 9e2a39edf279c882bb1b21abba2ccef5f1946c0b Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:12:05 +0200 Subject: [PATCH 05/27] fix: give invalid Nussinov pairs zero weight An out-of-range paired state is impossible and must contribute the additive zero of the partition semiring. Returning one incorrectly admitted an invalid pair as a multiplicative identity. --- src/IntaRNA/NussinovHandler.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IntaRNA/NussinovHandler.cpp b/src/IntaRNA/NussinovHandler.cpp index 37727add..15e4a190 100644 --- a/src/IntaRNA/NussinovHandler.cpp +++ b/src/IntaRNA/NussinovHandler.cpp @@ -94,7 +94,7 @@ NussinovHandler::getQb(const size_t i, const size_t j, const RnaSequence &seq, const Z_type bpWeight, const size_t minLoopLength, NussinovHandler::Z2dMatrix &Q, NussinovHandler::Z2dMatrix &Qb) { if (j >= seq.size()) { - return 1.0; + return 0.0; } if (i + minLoopLength >= j) { return 0.0; From 0f901cc9b992fd5fcd00f101d68039cec2fbcc45 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:12:38 +0200 Subject: [PATCH 06/27] fix: exclude the empty structure from base-pair ES getES1/getES2 promise the ensemble of intramolecular structures containing at least one pair. Subtract the unit-weight empty structure before converting Q to energy, while leaving getEall1/getEall2 based on the complete monomer partition. --- src/IntaRNA/InteractionEnergyBasePair.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/IntaRNA/InteractionEnergyBasePair.cpp b/src/IntaRNA/InteractionEnergyBasePair.cpp index 0ba248fe..3b38b887 100644 --- a/src/IntaRNA/InteractionEnergyBasePair.cpp +++ b/src/IntaRNA/InteractionEnergyBasePair.cpp @@ -26,7 +26,10 @@ void InteractionEnergyBasePair::computeES(const RnaSequence &seq, if (Z_equal(q_val, 1.0)) { logQ(i, j) = E_INF; } else { - logQ(i, j) = getE(q_val); + // getES* covers only structures containing at least one base pair. + // The full monomer partition Q also contains the empty structure with + // unit weight, which has to be removed here. + logQ(i, j) = getE(q_val - Z_type(1.0)); } } } From 302bd38d1ccd643bccfc7259eae84bfc3706a8b9 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:13:20 +0200 Subject: [PATCH 07/27] fix: count noLP heuristic ensemble paths once The direct continuation after the mandatory left stack already covers the adjacent (1,1) extension. Skip that one loop transition so the heuristic remains a subset of the exact interaction ensemble instead of assigning duplicate weight. --- src/IntaRNA/PredictorMfeEns2dHeuristic.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp index 8d33466b..df9cd01d 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp @@ -179,6 +179,12 @@ fillHybridZ() // iterate over all loop sizes w1 (seq1) and w2 (seq2) (minus 1) for (w1=1; w1-1 <= energy.getMaxInternalLoopSize1() && i1+w1+noLpShift Date: Sat, 15 Aug 2026 03:13:52 +0200 Subject: [PATCH 08/27] fix: apply site filters before ensemble accumulation Reject terminal-GU and excessive-ED sites in updateZ before they alter Zall or the boundary partition. Reported candidates and their requested interaction ensemble now use the same scientific domain. --- src/IntaRNA/PredictorMfeEns.cpp | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/IntaRNA/PredictorMfeEns.cpp b/src/IntaRNA/PredictorMfeEns.cpp index 20913d71..ac164136 100644 --- a/src/IntaRNA/PredictorMfeEns.cpp +++ b/src/IntaRNA/PredictorMfeEns.cpp @@ -46,6 +46,22 @@ updateZ( const size_t i1, const size_t j1 // check if something to be done if (Z_equal(partZ,0) || Z_isINF(Zall)) return; + + // Apply the same site filters used for MFE candidates before changing + // either the global or boundary-specific partition. + const OutputConstraint & outConstraint = output.getOutputConstraint(); + if (outConstraint.noGUend + && (energy.isGU(i1,i2) || energy.isGU(j1,j2))) + { + return; + } + if (outConstraint.maxED < Accessibility::ED_UPPER_BOUND + && (energy.getED1(i1,j1) > outConstraint.maxED + || energy.getED2(i2,j2) > outConstraint.maxED)) + { + return; + } + // handle whether or not partZ includes ED values or not Z_type partZ_withED = 0, partZ_noED = 0; if (isHybridZ) { From eae6555fbeb1bfb0152e06a690a72e32e70f2b1f Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:14:31 +0200 Subject: [PATCH 09/27] fix: isolate target base-pair accessibility limits Construct target AccessibilityBasePair with tIntLenMax and tAccW. Query settings must not truncate the target ED matrix; asymmetric target/query CLI options now retain their independent biological domains. --- src/bin/CommandLineParsing.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 7238f973..455a0f3f 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -2015,8 +2015,8 @@ getTargetAccessibility( const size_t sequenceNumber ) const case 'B' : // base-pair based accessibility return new AccessibilityBasePair( seq - , std::min( qIntLenMax.val == 0 ? seq.size() : qIntLenMax.val - , qAccW.val == 0 ? seq.size() : qAccW.val ) + , std::min( tIntLenMax.val == 0 ? seq.size() : tIntLenMax.val + , tAccW.val == 0 ? seq.size() : tAccW.val ) , &accConstraint ); @@ -2778,4 +2778,3 @@ getPersonality( int argc, char ** argv ) //////////////////////////////////////////////////////////////////////////// - From a1292058c4baeb4ca3963a08b75fdb0cef043d23 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:18:08 +0200 Subject: [PATCH 10/27] test: correct base-pair accessibility oracle --- tests/data/energyB-target-acc-options.testresult | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/data/energyB-target-acc-options.testresult b/tests/data/energyB-target-acc-options.testresult index b17376e3..bb1a566a 100644 --- a/tests/data/energyB-target-acc-options.testresult +++ b/tests/data/energyB-target-acc-options.testresult @@ -5,4 +5,4 @@ target;2;4;query;1;2;CGU&GG;(.(&));-2 1 0.000000e+00 NA NA NA 2 0.000000e+00 0.000000e+00 NA NA 3 0.000000e+00 0.000000e+00 0.000000e+00 NA -4 0.000000e+00 0.000000e+00 0.000000e+00 1.313186e+00 +4 0.000000e+00 0.000000e+00 0.000000e+00 0.000000e+00 From 18bdac74a50d09605506d6da37ea71b5347d27d8 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:21:05 +0200 Subject: [PATCH 11/27] fix: preserve no-seed ranges with outMinPu --- src/bin/CommandLineParsing.cpp | 5 ++--- tests/data/energyB-accN-noSeed-outMinPu.parameter | 13 +++++++++++++ tests/data/energyB-accN-noSeed-outMinPu.testresult | 2 ++ 3 files changed, 17 insertions(+), 3 deletions(-) create mode 100644 tests/data/energyB-accN-noSeed-outMinPu.parameter create mode 100644 tests/data/energyB-accN-noSeed-outMinPu.testresult diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 455a0f3f..949d3931 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -2634,7 +2634,7 @@ getQueryRanges( const InteractionEnergy & energy, const size_t sequenceNumber, c if (outMinPu.val > Z_type(0) && !Z_equal(outMinPu.val, Z_type(0))) { // decompose ranges based in minimal unpaired probability value per position // since all ranges covering a position will have a lower unpaired probability - acc.decomposeByMaxED( qRegion[sequenceNumber], energy.getE( outMinPu.val ), (noSeedRequired ? RnaSequence::lastPos : seedBP.val ) ); + acc.decomposeByMaxED( qRegion[sequenceNumber], energy.getE( outMinPu.val ), (noSeedRequired ? 1 : seedBP.val ) ); } return qRegion.at(sequenceNumber); @@ -2668,7 +2668,7 @@ getTargetRanges( const InteractionEnergy & energy, const size_t sequenceNumber, if (outMinPu.val > Z_type(0) && !Z_equal(outMinPu.val, Z_type(0))) { // decompose ranges based in minimal unpaired probability value per position // since all ranges covering a position will have a lower unpaired probability - acc.decomposeByMaxED( tRegion[sequenceNumber], energy.getE( outMinPu.val ), (noSeedRequired ? RnaSequence::lastPos : seedBP.val ) ); + acc.decomposeByMaxED( tRegion[sequenceNumber], energy.getE( outMinPu.val ), (noSeedRequired ? 1 : seedBP.val ) ); } return tRegion.at(sequenceNumber); @@ -2777,4 +2777,3 @@ getPersonality( int argc, char ** argv ) //////////////////////////////////////////////////////////////////////////// - diff --git a/tests/data/energyB-accN-noSeed-outMinPu.parameter b/tests/data/energyB-accN-noSeed-outMinPu.parameter new file mode 100644 index 00000000..55be72eb --- /dev/null +++ b/tests/data/energyB-accN-noSeed-outMinPu.parameter @@ -0,0 +1,13 @@ +target=GGGG +query=CCCC +energy=B +tAcc=N +qAcc=N +model=S +mode=M +noSeed=true +outMinPu=0.5 +outMode=C +outCsvCols=id1,start1,end1,id2,start2,end2,E +outNumber=1 +threads=1 diff --git a/tests/data/energyB-accN-noSeed-outMinPu.testresult b/tests/data/energyB-accN-noSeed-outMinPu.testresult new file mode 100644 index 00000000..486d9d08 --- /dev/null +++ b/tests/data/energyB-accN-noSeed-outMinPu.testresult @@ -0,0 +1,2 @@ +id1;start1;end1;id2;start2;end2;E +target;1;4;query;1;4;-4 From eefbda58b81198d4a5592f4ab1ece7b70a8081ba Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:23:11 +0200 Subject: [PATCH 12/27] fix: reject double-counted window partition output --- src/bin/CommandLineParsing.cpp | 7 ++++++- tests/data/energyB-accN-window-zall.parameter | 15 +++++++++++++++ tests/data/energyB-accN-window-zall.testresult | 1 + 3 files changed, 22 insertions(+), 1 deletion(-) create mode 100644 tests/data/energyB-accN-window-zall.parameter create mode 100644 tests/data/energyB-accN-window-zall.testresult diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 949d3931..b554bd86 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -1549,6 +1549,12 @@ parse(int argc, char** argv) if (outNumber.val > 1 && outOverlap.val != 'B') { throw error("window-based computation: non-overlapping subopt output (-n > 1) only supported for --outOverlap=B"); } + const bool windowNeedsZall = outMode.val == 'E' + || (outMode.val == 'C' + && OutputHandlerCsv::needsZall(OutputHandlerCsv::string2list(outCsvCols))); + if (windowNeedsZall) { + throw error("window-based computation cannot provide Zall/Eall output: overlapping windows count interactions more than once"); + } } @@ -2776,4 +2782,3 @@ getPersonality( int argc, char ** argv ) //////////////////////////////////////////////////////////////////////////// - diff --git a/tests/data/energyB-accN-window-zall.parameter b/tests/data/energyB-accN-window-zall.parameter new file mode 100644 index 00000000..1ec3131a --- /dev/null +++ b/tests/data/energyB-accN-window-zall.parameter @@ -0,0 +1,15 @@ +target=GGGGGGGGGGGG +query=CCCCCCCCCCCC +energy=B +tAcc=N +qAcc=N +model=S +mode=M +noSeed=true +tIntLenMax=4 +qIntLenMax=4 +windowWidth=10 +windowOverlap=4 +outMode=C +outCsvCols=Eall +threads=1 diff --git a/tests/data/energyB-accN-window-zall.testresult b/tests/data/energyB-accN-window-zall.testresult new file mode 100644 index 00000000..c069ea4b --- /dev/null +++ b/tests/data/energyB-accN-window-zall.testresult @@ -0,0 +1 @@ +# ERROR : window-based computation cannot provide Zall/Eall output: overlapping windows count interactions more than once From 77e01d378ddeb3d6ec2ac0c929c45f6d14956de9 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:23:59 +0200 Subject: [PATCH 13/27] fix: retain distinct equal-energy seeds --- src/IntaRNA/Interaction.h | 20 +++++++++++++++----- tests/Interaction_test.cpp | 24 ++++++++++++++++++++++++ 2 files changed, 39 insertions(+), 5 deletions(-) diff --git a/src/IntaRNA/Interaction.h b/src/IntaRNA/Interaction.h index ea44204d..6ef98063 100644 --- a/src/IntaRNA/Interaction.h +++ b/src/IntaRNA/Interaction.h @@ -175,16 +175,26 @@ class Interaction { E_type energy; /** - * order definition: first by increasing energy using increasing seq1 - * index as tie breaker. + * order definition: first by increasing energy, then lexicographically + * by all seed boundaries as tie breakers. * @param s the seed to compare to * @return true if this seed is considered smaller than s */ const bool operator < ( const Seed &s ) const { - return ( energy < s.energy - || (E_equal(energy,s.energy) && (bp_i.first < s.bp_i.first)) - ); + if (energy != s.energy) { + return energy < s.energy; + } + if (bp_i.first != s.bp_i.first) { + return bp_i.first < s.bp_i.first; + } + if (bp_i.second != s.bp_i.second) { + return bp_i.second < s.bp_i.second; + } + if (bp_j.first != s.bp_j.first) { + return bp_j.first < s.bp_j.first; + } + return bp_j.second < s.bp_j.second; } /** diff --git a/tests/Interaction_test.cpp b/tests/Interaction_test.cpp index c02cd412..977da5cd 100644 --- a/tests/Interaction_test.cpp +++ b/tests/Interaction_test.cpp @@ -125,4 +125,28 @@ TEST_CASE( "Interaction", "[Interaction]" ) { REQUIRE_FALSE( seeded == unseeded ); } + SECTION("seed ordering preserves distinct equal-energy ranges") { + + const E_type seedEnergy = -100; + const Interaction::Seed seed( + Interaction::BasePair(1,6), Interaction::BasePair(3,4), seedEnergy ); + const Interaction::Seed differentSeq2Right( + Interaction::BasePair(1,7), Interaction::BasePair(3,4), seedEnergy ); + const Interaction::Seed differentSeq1Right( + Interaction::BasePair(1,6), Interaction::BasePair(4,4), seedEnergy ); + const Interaction::Seed differentSeq2Left( + Interaction::BasePair(1,6), Interaction::BasePair(3,3), seedEnergy ); + Interaction::SeedSet seeds; + + REQUIRE( seeds.insert(seed).second ); + REQUIRE( seeds.insert(differentSeq2Right).second ); + REQUIRE( seeds.insert(differentSeq1Right).second ); + REQUIRE( seeds.insert(differentSeq2Left).second ); + REQUIRE( seeds.size() == 4 ); + + // Only an exact duplicate is equivalent in the ordering. + REQUIRE_FALSE( seeds.insert(seed).second ); + REQUIRE( seeds.size() == 4 ); + } + } From b7a4ba159ca376315f2adb7e37e1de9877345bd5 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:26:30 +0200 Subject: [PATCH 14/27] fix: exclude inaccessible split positions from ranges --- src/IntaRNA/Accessibility.cpp | 4 ++-- tests/AccessibilityConstraint_test.cpp | 21 ++++++++++++++++++++- 2 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/IntaRNA/Accessibility.cpp b/src/IntaRNA/Accessibility.cpp index 9c15f59f..c514549e 100644 --- a/src/IntaRNA/Accessibility.cpp +++ b/src/IntaRNA/Accessibility.cpp @@ -213,8 +213,8 @@ decomposeByMaxED( IndexRangeList & ranges, const E_type maxED, const size_t minR for (size_t i= range->from; i <= range->to; i++) { if (E_isINF(getED(i,i)) || (getED(i,i) > maxED && !E_equal(getED(i,i),maxED))) { // check if end of range found and to be stored - if (lastStart < i && minRangeLength <= (i +1 - lastStart)) { - out.push_back(IndexRange(lastStart,i)); + if (lastStart < i && minRangeLength <= (i - lastStart)) { + out.push_back(IndexRange(lastStart,i - 1)); } lastStart = range->to +1; } else { diff --git a/tests/AccessibilityConstraint_test.cpp b/tests/AccessibilityConstraint_test.cpp index 776909b6..53271bee 100644 --- a/tests/AccessibilityConstraint_test.cpp +++ b/tests/AccessibilityConstraint_test.cpp @@ -5,6 +5,7 @@ #undef NDEBUG #include "IntaRNA/AccessibilityConstraint.h" +#include "IntaRNA/AccessibilityDisabled.h" using namespace IntaRNA; @@ -132,6 +133,24 @@ TEST_CASE( "AccessibilityConstraint", "[AccessibilityConstraint]" ) { REQUIRE( vrnaStyle == "..xx..xx|.xx"); } + SECTION("decomposeByMaxED excludes forbidden positions") { + RnaSequence rna("test", "AAAAAAAA"); + AccessibilityConstraint constraint(rna, "...b...b", 0, "", "", ""); + AccessibilityDisabled accessibility(rna, 0, &constraint); + IndexRangeList ranges; + ranges.push_back(IndexRange(0, rna.size() - 1)); -} + accessibility.decomposeByMaxED(ranges, 0, 3); + + REQUIRE(ranges.size() == 2); + REQUIRE(ranges.get(0) == IndexRange(0, 2)); + REQUIRE(ranges.get(1) == IndexRange(4, 6)); + + ranges.clear(); + ranges.push_back(IndexRange(0, rna.size() - 1)); + accessibility.decomposeByMaxED(ranges, 0, 4); + REQUIRE(ranges.empty()); + } + +} From c132a613b9691a1c84f84e268476201e5e779d7e Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:28:20 +0200 Subject: [PATCH 15/27] fix: honor base-pair span in Vienna ES --- src/IntaRNA/InteractionEnergyVrna.cpp | 2 +- tests/AccessibilityVrna_test.cpp | 25 +++++++++++++++++++++++++ 2 files changed, 26 insertions(+), 1 deletion(-) diff --git a/src/IntaRNA/InteractionEnergyVrna.cpp b/src/IntaRNA/InteractionEnergyVrna.cpp index 348d90e3..33931bcc 100644 --- a/src/IntaRNA/InteractionEnergyVrna.cpp +++ b/src/IntaRNA/InteractionEnergyVrna.cpp @@ -112,7 +112,7 @@ computeES( const Accessibility & acc, InteractionEnergyVrna::EsMatrix & esToFill curModel.max_bp_span = -1; } // TODO check if VRNA_OPTION_WINDOW reasonable to speedup - vrna_fold_compound_t * foldData = vrna_fold_compound( sequence, &foldModel, VRNA_OPTION_PF); + vrna_fold_compound_t * foldData = vrna_fold_compound( sequence, &curModel, VRNA_OPTION_PF); // Adding hard constraints from pseudo dot-bracket unsigned int constraint_options = VRNA_CONSTRAINT_DB_DEFAULT; diff --git a/tests/AccessibilityVrna_test.cpp b/tests/AccessibilityVrna_test.cpp index f37f44b5..31b93231 100644 --- a/tests/AccessibilityVrna_test.cpp +++ b/tests/AccessibilityVrna_test.cpp @@ -4,7 +4,10 @@ #undef NDEBUG #include +#include "IntaRNA/AccessibilityDisabled.h" #include "IntaRNA/AccessibilityVrna.h" +#include "IntaRNA/InteractionEnergyVrna.h" +#include "IntaRNA/ReverseAccessibility.h" using namespace IntaRNA; @@ -28,4 +31,26 @@ TEST_CASE("AccessibilityVrna", "[AccessibilityVrna]") { REQUIRE( E_equal( acc.getED(0, 1), 0 ) ); } + + + SECTION("InteractionEnergyVrna respects accessibility base-pair span") { + + const std::string seq = "GGGGAAAACCCC"; + RnaSequence rna("test", seq); + AccessibilityConstraint unrestrictedConstraint(rna.size(), 0, "", "", ""); + AccessibilityConstraint shortSpanConstraint(rna.size(), 3, "", "", ""); + AccessibilityDisabled unrestrictedAcc(rna, rna.size(), &unrestrictedConstraint); + AccessibilityDisabled shortSpanAcc(rna, rna.size(), &shortSpanConstraint); + ReverseAccessibility unrestrictedAccReversed(unrestrictedAcc); + ReverseAccessibility shortSpanAccReversed(shortSpanAcc); + VrnaHandler vrnaHandler(37, "Turner04", false, false); + + InteractionEnergyVrna unrestrictedEnergy( + unrestrictedAcc, unrestrictedAccReversed, vrnaHandler, 16, 16, true); + InteractionEnergyVrna shortSpanEnergy( + shortSpanAcc, shortSpanAccReversed, vrnaHandler, 16, 16, true); + + REQUIRE( E_isNotINF(unrestrictedEnergy.getES1(0, rna.size()-1)) ); + REQUIRE( E_isINF(shortSpanEnergy.getES1(0, rna.size()-1)) ); + } } From 595029ff388b87d7d25709dcd592583f313411da Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:29:34 +0200 Subject: [PATCH 16/27] fix: release ViennaRNA ensemble resources --- src/IntaRNA/InteractionEnergyVrna.cpp | 44 ++++++++++++++++++++------- tests/AccessibilityVrna_test.cpp | 5 +++ 2 files changed, 38 insertions(+), 11 deletions(-) diff --git a/src/IntaRNA/InteractionEnergyVrna.cpp b/src/IntaRNA/InteractionEnergyVrna.cpp index 33931bcc..16a02459 100644 --- a/src/IntaRNA/InteractionEnergyVrna.cpp +++ b/src/IntaRNA/InteractionEnergyVrna.cpp @@ -3,6 +3,8 @@ #include "IntaRNA/AccessibilityVrna.h" #include +#include +#include #include // ES computation @@ -16,6 +18,26 @@ extern "C" { } +namespace { + +struct VrnaAllocatedDeleter { + void operator()( char * data ) const { + free(data); + } +}; + +struct VrnaFoldCompoundDeleter { + void operator()( vrna_fold_compound_t * foldCompound ) const { + vrna_fold_compound_free(foldCompound); + } +}; + +typedef std::unique_ptr VrnaAllocatedPtr; +typedef std::unique_ptr VrnaFoldCompoundPtr; + +} + + namespace IntaRNA { @@ -94,8 +116,10 @@ computeES( const Accessibility & acc, InteractionEnergyVrna::EsMatrix & esToFill const Z_type RT = getRT(); // VRNA compatible data structures - char * sequence = (char *) vrna_alloc(sizeof(char) * (seqLength + 1)); - char * structureConstraint = (char *) vrna_alloc(sizeof(char) * (seqLength + 1)); + VrnaAllocatedPtr sequenceOwner( (char *) vrna_alloc(sizeof(char) * (seqLength + 1)) ); + VrnaAllocatedPtr structureConstraintOwner( (char *) vrna_alloc(sizeof(char) * (seqLength + 1)) ); + char * const sequence = sequenceOwner.get(); + char * const structureConstraint = structureConstraintOwner.get(); for (int i=0; i Date: Sat, 15 Aug 2026 03:41:35 +0200 Subject: [PATCH 17/27] test: expose seed ensemble algebra and reuse defects --- tests/PredictorMfeEnsRegression_test.cpp | 80 ++++++++++++++++++++++++ 1 file changed, 80 insertions(+) diff --git a/tests/PredictorMfeEnsRegression_test.cpp b/tests/PredictorMfeEnsRegression_test.cpp index 600b4498..b70326af 100644 --- a/tests/PredictorMfeEnsRegression_test.cpp +++ b/tests/PredictorMfeEnsRegression_test.cpp @@ -7,11 +7,31 @@ #include "IntaRNA/OutputHandlerInteractionList.h" #include "IntaRNA/PredictorMfeEns2d.h" #include "IntaRNA/PredictorMfeEns2dHeuristic.h" +#include "IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.h" +#include "IntaRNA/PredictorMfeEns2dSeedExtension.h" #include "IntaRNA/ReverseAccessibility.h" #include "IntaRNA/RnaSequence.h" +#include "IntaRNA/SeedHandlerNoBulge.h" + +#include using namespace IntaRNA; +template +class InspectableEnsemblePredictor : public PredictorType { +public: + InspectableEnsemblePredictor( + const InteractionEnergy & energy, + OutputHandler & output, + SeedHandler * seedHandler) + : PredictorType(energy, output, NULL, seedHandler) + {} + + size_t getPartitionCount() const { + return this->Z_partition.size(); + } +}; + TEST_CASE("ensemble predictor regressions", "[PredictorMfeEns]") { #include "testEasyLoggingSetup.icc" @@ -55,4 +75,64 @@ TEST_CASE("ensemble predictor regressions", "[PredictorMfeEns]") { REQUIRE(predictor.getZall() == 0.0); } + + SECTION("noLP seed extension multiplies stacked partition factors") { + RnaSequence target("target", "GGG"); + RnaSequence query("query", "CCC"); + AccessibilityDisabled targetAcc(target, 0, NULL); + AccessibilityDisabled queryAcc(query, 0, NULL); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc); + IndexRangeList targetSeedRange; + targetSeedRange.push_back(IndexRange(1, 2)); + IndexRangeList querySeedRange; + querySeedRange.push_back(IndexRange(1, 2)); + SeedConstraint seedConstraint(2, 0, 0, 0, + E_INF, Accessibility::ED_UPPER_BOUND, E_INF, + targetSeedRange, querySeedRange, "", false, false, true); + OutputConstraint constraint(1, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, true, false, true, false); + OutputHandlerInteractionList out(constraint, 1); + PredictorMfeEns2dSeedExtension predictor( + energy, out, NULL, new SeedHandlerNoBulge(energy, seedConstraint)); + + predictor.predict(); + + // The only seed starts at internal (1,1), contributing exp(2). + // Its sole noLP extension stacks (0,0) to the left, contributing exp(3). + const Z_type expectedZ = std::exp(2.0) + std::exp(3.0); + REQUIRE(predictor.getZall() == Approx(expectedZ).epsilon(1e-12)); + } + + SECTION("seed-extension predictor reuse clears an empty range partition") { + RnaSequence target("target", "GGG"); + RnaSequence query("query", "CCC"); + AccessibilityDisabled targetAcc(target, 0, NULL); + AccessibilityDisabled queryAcc(query, 0, NULL); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc); + SeedConstraint seedConstraint(2, 0, 0, 0, + E_INF, Accessibility::ED_UPPER_BOUND, E_INF, + IndexRangeList(), IndexRangeList(), "", false, false, true); + OutputConstraint constraint(1, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, true, false, true, false); + + OutputHandlerInteractionList exactOut(constraint, 1); + InspectableEnsemblePredictor exact( + energy, exactOut, new SeedHandlerNoBulge(energy, seedConstraint)); + exact.predict(); + REQUIRE(exact.getPartitionCount() > 0); + exact.predict(IndexRange(0, 0), IndexRange(0, 0)); + REQUIRE(exact.getPartitionCount() == 0); + REQUIRE(exact.getZall() == 0.0); + + OutputHandlerInteractionList heuristicOut(constraint, 1); + InspectableEnsemblePredictor heuristic( + energy, heuristicOut, new SeedHandlerNoBulge(energy, seedConstraint)); + heuristic.predict(); + REQUIRE(heuristic.getPartitionCount() > 0); + heuristic.predict(IndexRange(0, 0), IndexRange(0, 0)); + REQUIRE(heuristic.getPartitionCount() == 0); + REQUIRE(heuristic.getZall() == 0.0); + } } From 1ed5819bffdab0403063054737578455bf41de94 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:43:33 +0200 Subject: [PATCH 18/27] fix: multiply stacked seed partition factors --- src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp index 250339aa..e0175d97 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp @@ -229,7 +229,7 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) // get stacking energy to avoid recomputation in recursion below iStackZ = energy.getBoltzmannWeight(energy.getE_interLeft(i1,i1+noLpShift,i2,i2+noLpShift)); // check just stacked - curZ += iStackZ + hybridZ_left(l1-noLpShift,l2-noLpShift); + curZ += iStackZ * hybridZ_left(l1-noLpShift,l2-noLpShift); } // check all combinations of decompositions into (i1,i2)..(k1,k2)-(j1,j2) From a975c4d42893055cb866c0f88ca0c69213fe90df Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:44:53 +0200 Subject: [PATCH 19/27] fix: clear seed ensemble state on empty ranges --- src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp | 1 + src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp | 1 + 2 files changed, 2 insertions(+) diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp index 9fd05b9e..4cf8aeaa 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp @@ -61,6 +61,7 @@ predict( const IndexRange & r1, const IndexRange & r2 ) if (seedHandler.fillSeed( 0, range_size1-1, 0, range_size2-1 ) == 0) { // trigger empty interaction reporting initOptima(); + initZ(); reportOptima(); // stop computation return; diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp index e0175d97..05a4f852 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp @@ -62,6 +62,7 @@ predict( const IndexRange & r1, const IndexRange & r2 ) if (seedHandler.fillSeed( 0, range_size1-1, 0, range_size2-1 ) == 0) { // trigger empty interaction reporting initOptima(); + initZ(); reportOptima(); // stop computation return; From 70920c6d76cab1f1fb39761ada0352b7a31c874f Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 03:46:28 +0200 Subject: [PATCH 20/27] docs: finalize phase 1 audit and corrections --- doc/refactor/1-current-state.md | 122 +++++++++++++++++++++++--------- 1 file changed, 87 insertions(+), 35 deletions(-) diff --git a/doc/refactor/1-current-state.md b/doc/refactor/1-current-state.md index a63dd8df..358afad1 100644 --- a/doc/refactor/1-current-state.md +++ b/doc/refactor/1-current-state.md @@ -2,7 +2,8 @@ This document records the phase-1 audit requested in discussion #232. It is about scientific behavior and runtime cost, not a style rewrite. The audited -baseline is IntaRNA 3.4.1, commit `3b14bc0`. +baseline is IntaRNA 3.4.1, commit `3b14bc0`. The correction and regression +record below is current through `a975c4d`. ## Reproducible baseline @@ -18,6 +19,17 @@ Thus the old suite is a useful output-stability gate, but it is not broad enough to establish correctness of the ensemble predictors, ViennaRNA paths, parallel/window reduction, or several public boundary cases. +After the phase-1 regressions and local corrections through `a975c4d`: + +- the API binary passes 3,895 assertions in 31 test cases; and +- all 20 command-line golden cases pass. + +The added CLI cases cover asymmetric target/query accessibility options, +seed-free `outMinPu` range filtering, and rejection of partition output from +overlapping windows. The API additions cover the local ownership, boundary, +partition-filter, seed-extension, range-splitting and ViennaRNA cases listed +below. + ## Architecture and data flow `src/bin/IntaRNA.cpp` owns one invocation. `CommandLineParsing` parses the @@ -84,19 +96,30 @@ calls from existing recurrences before changing their mathematics. ## Confirmed defects and regressions -The following findings were traced in this source tree. Tests are added before -fixes, as required by phase 1. - -| Finding | Evidence / consequence | Phase-1 regression | -| --- | --- | --- | -| `Interaction::operator=` is not self-assignment safe | it clears `basePairs` before reading the same object | preserve pairs, energy and seed metadata across `x=x` | -| interaction equality dereferences asymmetric null seed pointers | `seed == i.seed || *seed == *i.seed` dereferences when only one side is null | seeded and unseeded interactions compare unequal without crashing | -| zero-capacity interaction storage dereferences an empty reverse iterator | `maxToStore=0` reaches `*storage.rbegin()` | count a report but retain no interaction | -| `NussinovHandler::getQb` returns one for an out-of-range paired interval | an invalid pair obtains the multiplicative identity instead of zero weight | `getQb(i,n)==0` | -| base-pair `getES*` includes the empty monomer structure | it stores `-RT log(Q)`, while the API specifies structures containing at least one pair, i.e. `-RT log(Q-1)` | a four-base sequence with one possible pair has `ES=-1`, not `-log(1+e)` | -| noLP heuristic ensemble counts the direct stack twice | the explicit direct continuation and the `w1=w2=1` loop iteration denote the same path | on `GGGG/CCCC`, heuristic `Zall` must not exceed exact `Zall` | -| ensemble `updateZ` bypasses `noGUend` and `maxED` | its direct `Zall`/boundary-map update does not use the filters in `PredictorMfe::updateOptima` | a single terminal GU contributes zero when terminal GU is forbidden | -| target base-pair accessibility uses query limits | target factory reads `qIntLenMax/qAccW` instead of target parameters | asymmetric CLI target-accessibility table retains the target maximum length | +The first eight findings below were exposed by `0b20908` before their +individual source fixes. Later local findings carry focused API or CLI evidence +in the same change sequence. The corrected accessibility oracle is called out +separately because it was a test-data error, not a production defect. + +| Finding | Evidence / consequence | Phase-1 regression evidence | Correction | +| --- | --- | --- | --- | +| `Interaction::operator=` was not self-assignment safe | it cleared `basePairs` before reading the same object | `x=x` preserves pairs, energy and owned seed metadata | `11ac391` returns immediately for self-assignment | +| interaction equality dereferenced asymmetric optional seeds | `seed == i.seed || *seed == *i.seed` dereferenced when only one pointer was null | seeded and unseeded interactions compare unequal in both operand orders without crashing | `f498659` dereferences only when both seeds are non-null | +| zero-capacity interaction storage dereferenced an empty reverse iterator | `maxToStore=0` reached `*storage.rbegin()` | adding an interaction increments the report count but stores nothing | `8db8b0c` gates all storage access on positive capacity | +| `NussinovHandler::getQb` returned one for an out-of-range paired interval | an invalid pair obtained the multiplicative identity instead of zero weight | `getQb(0,n)==0` on a four-base sequence | `9e2a39e` returns zero for `j>=n` | +| base-pair `getES*` included the empty monomer structure | it stored `-RT log(Q)`, while the API specifies structures containing at least one pair, i.e. `-RT log(Q-1)` | `ACGU` with one admissible pair has `ES=-1`; intervals without a pair have infinite ES | `0f901cc` removes the empty structure's unit weight | +| the noLP heuristic ensemble counted the direct stack twice | the explicit direct continuation and the `w1=w2=1` loop iteration denoted the same paths | on `GGGG/CCCC`, heuristic `Zall` does not exceed exact `Zall` | `302bd38` skips that duplicate loop term when noLP is active | +| ensemble `updateZ` bypassed site filters | its direct global and boundary-partition updates did not apply `noGUend` or `maxED` | a single terminal GU contributes zero when terminal GU is forbidden | `f9612c3` applies both filters before either partition update | +| target base-pair accessibility used query limits | the target factory read `qIntLenMax/qAccW` instead of target parameters | an asymmetric CLI case with target limit 4 and query limit 3 retains the target length-4 accessibility column | `eae6555` uses `tIntLenMax/tAccW` | +| the initial asymmetric-accessibility oracle assigned a nonzero ED to the full `ACGU` target | at the default minimum loop length this four-base target has no admissible intramolecular pair, so its unpaired probability is one | the golden target-accessibility table expects zero ED for every interval, including length 4 | `a129205` corrects the test oracle only | +| seed-free `outMinPu` filtering removed every range | `RnaSequence::lastPos` was used as the minimum resulting length when no seed was required | a `GGGG/CCCC` CLI case with `noSeed=true,outMinPu=0.5` retains and reports the full interaction | `18bdac7` uses minimum length one for seed-free prediction | +| overlapping windows could report a double-counted global partition | the same interaction can occur in more than one window, so `Zall/Eall` cannot be summed safely | a CLI case requesting CSV `Eall` in window mode expects a deterministic error | `eefbda5` rejects ensemble output and CSV columns needing `Zall` in window mode | +| seed ordering collapsed distinct equal-energy seed ranges | the set comparator used only energy and the first sequence-1 boundary | four seeds differing in any remaining boundary all survive; only an exact duplicate is rejected | `77e01d3` uses all four boundaries as lexicographic tie breakers | +| accessibility range splitting included the forbidden split position | the closed prefix ended at `i` and its length counted `i`, although singleton ED at `i` exceeded the threshold | blocked positions 3 and 7 split `[0,7]` into `[0,2]` and `[4,6]`; minimum length four retains neither | `b7a4ba1` measures `i-lastStart` and closes at `i-1` | +| ViennaRNA ES ignored the accessibility base-pair-span model | `computeES` configured `curModel.max_bp_span` but passed the unmodified model to `vrna_fold_compound` | unrestricted `GGGGAAAACCCC` has finite ES while span 3 makes full-range ES infinite | `c132a61` passes the configured model | +| ViennaRNA ensemble temporaries lacked complete scoped ownership | `computeIntraEall` leaked its allocated sequence and fold compound on normal return; `computeES` cleanup was not exception-safe | the API suite exercises span-sensitive ES plus both `getEall1/getEall2` paths; direct leak detection still belongs to sanitizer coverage | `595029f` gives Vienna allocations and fold compounds scoped deleters | +| exact noLP seed extension added independent partition factors | the stack weight and remaining left subensemble were added, violating the sum-product recurrence | with one allowed two-pair seed in a `3x3` complementary grid, the exact partition is `exp(2)+exp(3)` | `1ed5819` multiplies the stack and subensemble weights | +| empty seed-extension ranges retained prior boundary partitions | early return reset scalar optima but did not clear `Z_partition`, so predictor reuse replayed stale sites | exact and heuristic predictors run on a seeded range and then a singleton range; the second partition and boundary store are zero | `a975c4d` calls `initZ()` on both no-seed paths | Additional high-confidence findings are not treated as tiny local fixes because they change aggregation or recurrence ownership and need the benchmark gates of @@ -104,45 +127,74 @@ later phases: - exact ensemble four-boundary storage violates the advertised practical two-dimensional memory bound; -- seed-extension ensemble code has early-return state reuse, additive where - independent partition factors must be multiplicative, noLP duplication and - unsigned boundary risks; -- windowed/global partition and tracker aggregation can count overlapping - domains more than once; +- the two proven seed-extension algebra/state defects are corrected, but the + family still needs a broader recurrence and traceback oracle for overlapping + bulged seeds, noLP corrections and unsigned boundary arithmetic; +- window-mode `Zall/Eall` is now rejected, but aggregation across general + multiple or overlapping input regions and tracker domains still needs an + ownership/counting audit; - query/target range factories mutate parser-owned vectors from `const` methods, repeat decomposition and can race across pair tasks; - partition accumulation can overflow in release builds; debug-only warnings neither prevent nor repair the result; and -- CLI object ownership uses raw factories and `const_cast` cleanup, making +- local ViennaRNA ensemble allocations now have scoped ownership, but CLI + factories still return raw objects and use `const_cast` cleanup, making exceptional and parallel paths difficult to reason about. ## Missing coverage -The original suite has no direct test for exact MFE, exact seeded MFE, any -ensemble predictor, seed-extension predictors, zero requested output, or -partition/filter agreement. It also lacks: - -- a brute-force small-instance partition oracle; -- reuse tests for stateful predictors; -- asymmetric target/query option tests; -- multiple-region and overlapping-window aggregation tests; +At the audited baseline the suite had no direct test for exact MFE, exact +seeded MFE, any ensemble predictor, seed-extension predictors, zero requested +output, or partition/filter agreement. Phase 1 now has focused evidence for +zero-capacity output, one exact-versus-heuristic noLP ensemble comparison, +terminal-GU partition filtering, seed-extension sum-product/reuse behavior, +asymmetric target/query options, seed-free range filtering, window-partition +rejection, accessibility splitting, and the span-sensitive ViennaRNA ES path. + +Coverage still lacks: + +- a brute-force small-instance partition oracle and direct `maxED` partition + filter regression; +- direct exact-MFE and exact-seeded-MFE regressions, plus wider bulged-seed and + heuristic seed-extension oracles; +- reuse tests for other stateful predictor families; +- general multiple-region and overlapping-domain aggregation tests; - deterministic threaded-output tests; - optimized/release invalid-input tests; -- sanitizer coverage; and +- sanitizer/leak coverage for the scoped ViennaRNA ownership change; and - a GCC 14 plus macOS Clang portability build. -Phase 1 adds focused deterministic cases for the confirmed local defects. -Phase 2 adds compiler/build coverage. The benchmark and differential corpus -introduced before phase 3 supplies executable-level output parity for every -performance change. +The local seed-extension sum-product and lifecycle audit is complete; broader +seed/traceback enumeration remains a phase-2 correctness gate. Phase 2 also +adds compiler/build coverage. The benchmark and differential corpus introduced +before phase 3 supplies executable-level output parity for every performance +change. ## Change record - Created the required `refactoring` branch from `master` at `3b14bc0`. - Added this architecture, correctness and performance audit. - Established the clean 30-case / 3,859-assertion API and 17-case CLI baseline. -- Added table-independent regressions for ownership, boundary, partition and - target/query isolation defects before their fixes. +- `0b20908` added table-independent regressions for ownership, boundary, + partition and target/query isolation defects before their fixes. +- `11ac391` and `f498659` made interaction assignment and optional-seed + comparison safe. +- `8db8b0c`, `9e2a39e` and `0f901cc` corrected zero-capacity storage, + invalid paired-interval weight and base-pair ES semantics. +- `302bd38` and `f9612c3` removed the noLP duplicate contribution and made + ensemble accumulation honor site filters. +- `eae6555` isolated target accessibility limits; `a129205` then corrected + the scientific golden oracle used by that regression. +- `18bdac7` retained seed-free ranges under `outMinPu`, and `eefbda5` + rejected undefined global partition output from overlapping windows. +- `77e01d3` retained distinct equal-energy seeds, and `b7a4ba1` excluded + inaccessible split positions from returned ranges. +- `c132a61` honored the configured ViennaRNA base-pair span, and `595029f` + introduced scoped ownership for the ViennaRNA ensemble resources. +- `de135fd` exposed seed-extension algebra and reuse defects before `1ed5819` + restored multiplication and `a975c4d` cleared stale partitions. +- Established the current 31-case / 3,895-assertion API and 20-case CLI gate at + `a975c4d`. The discussion once calls the phase-1 document `refactor-changelog.md`; no such file exists and the same phase otherwise consistently requires From 08d375ed4be9825bc64724116b167a3dc58e4b95 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 05:47:14 +0200 Subject: [PATCH 21/27] test: expose stale heuristic cell state --- tests/Makefile.am | 1 + tests/PredictorMfeHeuristicCellState_test.cpp | 189 ++++++++++++++++++ 2 files changed, 190 insertions(+) create mode 100644 tests/PredictorMfeHeuristicCellState_test.cpp diff --git a/tests/Makefile.am b/tests/Makefile.am index 07d0bcab..70e5c4b8 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -46,6 +46,7 @@ runApiTests_SOURCES = \ PredictionTrackerSpotProb_test.cpp \ PredictorMfe2dHelixBlockHeuristic_test.cpp \ PredictorMfe2dHelixBlockHeuristicSeed_test.cpp \ + PredictorMfeHeuristicCellState_test.cpp \ PredictorMfeEnsRegression_test.cpp \ NussinovHandler_test.cpp \ RnaSequence_test.cpp \ diff --git a/tests/PredictorMfeHeuristicCellState_test.cpp b/tests/PredictorMfeHeuristicCellState_test.cpp new file mode 100644 index 00000000..2e66ccfa --- /dev/null +++ b/tests/PredictorMfeHeuristicCellState_test.cpp @@ -0,0 +1,189 @@ +#include "catch.hpp" + +#undef NDEBUG + +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/OutputHandlerInteractionList.h" +#include "IntaRNA/PredictorMfe2dHeuristic.h" +#include "IntaRNA/PredictorMfe2dHeuristicSeed.h" +#include "IntaRNA/PredictorMfeEns2dHeuristic.h" +#include "IntaRNA/ReverseAccessibility.h" +#include "IntaRNA/RnaSequence.h" +#include "IntaRNA/SeedConstraint.h" +#include "IntaRNA/SeedHandlerNoBulge.h" + +#include + +using namespace IntaRNA; + +namespace { + +class InspectableMfeEns2dHeuristic : public PredictorMfeEns2dHeuristic { +public: + InspectableMfeEns2dHeuristic(const InteractionEnergy & energy, + OutputHandler & output) + : PredictorMfeEns2dHeuristic(energy, output, NULL) + {} + + Z_type getBoundaryZ(const size_t i1, const size_t j1, + const size_t i2, const size_t j2) const + { + const auto entry = Z_partition.find(Interaction::Boundary(i1, j1, i2, j2)); + return entry == Z_partition.end() ? Z_type(0) : entry->second; + } +}; + +void requireThreePairStack(const OutputHandlerInteractionList & output) { + REQUIRE_FALSE(output.empty()); + const Interaction & interaction = **output.begin(); + REQUIRE(interaction.energy == Ekcal_2_E(-3.0)); + REQUIRE(interaction.basePairs.front() == Interaction::BasePair(0, 3)); + REQUIRE(interaction.basePairs.back() == Interaction::BasePair(2, 1)); +} + +} // namespace + +TEST_CASE("heuristic cells reset their incumbent energy", "[PredictorMfeHeuristicCellState]") { + + #include "testEasyLoggingSetup.icc" + + // In reversed query coordinates, the diagonal interaction is GC-GU-GC. + // The GU pair is internal and therefore valid when terminal GU pairs are + // filtered. With stacking-only loops, its energy is exactly -3 kcal/mol. + RnaSequence target("target", "GGG"); + RnaSequence query("query", "CCUC"); + AccessibilityDisabled targetAcc(target, 0, NULL); + AccessibilityDisabled queryAcc(query, 0, NULL); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc, 0, 0); + OutputConstraint constraint(1, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, false, true, true, true); + + SECTION("unseeded MFE retains an internal GU continuation") { + OutputHandlerInteractionList output(constraint, 1); + PredictorMfe2dHeuristic predictor(energy, output, NULL); + + predictor.predict(); + + requireThreePairStack(output); + REQUIRE((**output.begin()).basePairs.size() == 3); + REQUIRE((**output.begin()).basePairs.at(1) + == Interaction::BasePair(1, 2)); + } + + SECTION("ensemble heuristic retains the internal GU boundary partition") { + OutputHandlerInteractionList output(constraint, 1); + PredictorMfeEns2dHeuristic predictor(energy, output, NULL); + + predictor.predict(); + + requireThreePairStack(output); + // Ensemble traceback intentionally reports boundaries only. + REQUIRE((**output.begin()).basePairs.size() == 2); + } + + SECTION("seeded MFE can extend a seed through an internal GU cell") { + // Only the leading GC-GC seed is admissible. Its right extension is + // GC-GU-GC in the unseeded matrix, so the seed cannot bypass the + // poisoned cell via a later seed start. + RnaSequence seedTarget("seedTarget", "GGGG"); + RnaSequence seedQuery("seedQuery", "CCUCC"); + AccessibilityDisabled seedTargetAcc(seedTarget, 0, NULL); + AccessibilityDisabled seedQueryAcc(seedQuery, 0, NULL); + ReverseAccessibility reverseSeedQueryAcc(seedQueryAcc); + InteractionEnergyBasePair seedEnergy( + seedTargetAcc, reverseSeedQueryAcc, 0, 0); + IndexRangeList targetSeedRange; + targetSeedRange.push_back(IndexRange(0, 1)); + IndexRangeList querySeedRange; + querySeedRange.push_back(IndexRange(0, 1)); + SeedConstraint seedConstraint(2, 0, 0, 0, + E_INF, Accessibility::ED_UPPER_BOUND, E_INF, + targetSeedRange, querySeedRange, "", false, false, true); + OutputConstraint seedOutputConstraint(1, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, true, true, true, true); + OutputHandlerInteractionList output(seedOutputConstraint, 1); + PredictorMfe2dHeuristicSeed predictor(seedEnergy, output, NULL, + new SeedHandlerNoBulge(seedEnergy, seedConstraint)); + + predictor.predict(); + + REQUIRE_FALSE(output.empty()); + const Interaction & interaction = **output.begin(); + REQUIRE(interaction.energy == Ekcal_2_E(-4.0)); + REQUIRE(interaction.basePairs.size() == 4); + REQUIRE(interaction.basePairs.front() == Interaction::BasePair(0, 4)); + REQUIRE((**output.begin()).basePairs.at(1) + == Interaction::BasePair(1, 3)); + REQUIRE((**output.begin()).basePairs.at(2) + == Interaction::BasePair(2, 2)); + REQUIRE(interaction.basePairs.back() == Interaction::BasePair(3, 1)); + } +} + +TEST_CASE("ensemble noLP heuristic keeps valid non-direct extensions", + "[PredictorMfeHeuristicCellState][PredictorMfeEns2dHeuristic]") { + + #include "testEasyLoggingSetup.icc" + + SECTION("a GU mandatory stack can be internal to non-GU boundaries") { + // Reversed query CUCC gives the unique best GC-GU-GC path at + // internal boundary (0,2,0,2). The mandatory second pair is GU, + // but the actual interaction right end is the following GC pair. + RnaSequence target("target", "GGG"); + RnaSequence query("query", "CCUC"); + AccessibilityDisabled targetAcc(target, 0, NULL); + AccessibilityDisabled queryAcc(query, 0, NULL); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc, 0, 0); + OutputConstraint constraint(1, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, true, true, true, true); + OutputHandlerInteractionList output(constraint, 1); + InspectableMfeEns2dHeuristic predictor(energy, output); + + predictor.predict(); + + const Z_type pathZ = std::exp(3.0); + const Z_type expectedZall = Z_type(2) * std::exp(2.0) + pathZ; + REQUIRE(predictor.getBoundaryZ(0, 2, 0, 2) + == Approx(pathZ).epsilon(1e-12)); + REQUIRE(predictor.getZall() == Approx(expectedZall).epsilon(1e-12)); + REQUIRE_FALSE(output.empty()); + const Interaction & interaction = **output.begin(); + REQUIRE(interaction.energy == Ekcal_2_E(-3.0)); + REQUIRE(interaction.basePairs.size() == 2); + REQUIRE(interaction.basePairs.front() == Interaction::BasePair(0, 3)); + REQUIRE(interaction.basePairs.back() == Interaction::BasePair(2, 1)); + } + + SECTION("an absent direct extension does not suppress a later bulge") { + // The first GC-GC block cannot continue directly because A-C at + // internal (2,2) is impossible. It can still cross the target A + // bulge to the second GC-GC block via loop offset (2,1). + RnaSequence target("target", "GGAGG"); + RnaSequence query("query", "CCCC"); + AccessibilityDisabled targetAcc(target, 0, NULL); + AccessibilityDisabled queryAcc(query, 0, NULL); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc, 1, 1); + OutputConstraint constraint(1, OutputConstraint::OVERLAP_BOTH, + E_INF, E_INF, false, true, true, true, true); + OutputHandlerInteractionList output(constraint, 1); + InspectableMfeEns2dHeuristic predictor(energy, output); + + predictor.predict(); + + const Z_type pathZ = std::exp(4.0); + const Z_type expectedZall = Z_type(6) * std::exp(2.0) + pathZ; + REQUIRE(predictor.getBoundaryZ(0, 4, 0, 3) + == Approx(pathZ).epsilon(1e-12)); + REQUIRE(predictor.getZall() == Approx(expectedZall).epsilon(1e-12)); + REQUIRE_FALSE(output.empty()); + const Interaction & interaction = **output.begin(); + REQUIRE(interaction.energy == Ekcal_2_E(-4.0)); + REQUIRE(interaction.basePairs.size() == 2); + REQUIRE(interaction.basePairs.front() == Interaction::BasePair(0, 3)); + REQUIRE(interaction.basePairs.back() == Interaction::BasePair(4, 0)); + } +} From 96c14b30d32e310b3da47c7d35ee9963fdb2c85e Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 05:48:18 +0200 Subject: [PATCH 22/27] fix: reset heuristic cell state --- src/IntaRNA/PredictorMfe2dHeuristic.cpp | 1 + src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp | 1 + src/IntaRNA/PredictorMfeEns2dHeuristic.cpp | 29 +++++++++------------ 3 files changed, 15 insertions(+), 16 deletions(-) diff --git a/src/IntaRNA/PredictorMfe2dHeuristic.cpp b/src/IntaRNA/PredictorMfe2dHeuristic.cpp index d0b12253..8a9a99a7 100644 --- a/src/IntaRNA/PredictorMfe2dHeuristic.cpp +++ b/src/IntaRNA/PredictorMfe2dHeuristic.cpp @@ -99,6 +99,7 @@ fillHybridE() // init as invalid boundary *curCell = BestInteractionE(E_INF, RnaSequence::lastPos, RnaSequence::lastPos); + curCellEtotal = E_INF; // check if positions can form interaction if ( energy.isAccessible1(i1) diff --git a/src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp b/src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp index 78b80262..2a32ad76 100644 --- a/src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp +++ b/src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp @@ -116,6 +116,7 @@ fillHybridE() // init as invalid boundary *curCell = BestInteractionE(E_INF, RnaSequence::lastPos, RnaSequence::lastPos); *curCellSeed = BestInteractionE(E_INF, RnaSequence::lastPos, RnaSequence::lastPos); + curCellEtotal = E_INF; // check if positions can form interaction if ( energy.isAccessible1(i1) diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp index df9cd01d..1383be43 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp @@ -101,6 +101,7 @@ fillHybridZ() // init as invalid boundary *curCell = BestInteractionZ(0.0, RnaSequence::lastPos, RnaSequence::lastPos); + curCellEtotal = E_INF; // check if positions can form interaction if ( energy.isAccessible1(i1) @@ -134,24 +135,20 @@ fillHybridZ() // update overall partition function information for initial bps only updateZ( i1,curCell->j1, i2,curCell->j2, curCell->val, true ); - if(outConstraint.noLP) { - ///////////////////////////////////////// - // check direct extension to the right of the noLP stacking - ///////////////////////////////////////// + } - // direct cell access (const) - rightExt = &(hybridZ(i1+noLpShift,i2+noLpShift)); - // check if right side can pair - if (Z_equal(rightExt->val, 0.0)) { - continue; - } - // check if interaction length is within boundary - if ( (rightExt->j1 +1 -i1) > energy.getAccessibility1().getMaxLength() - || (rightExt->j2 +1 -i2) > energy.getAccessibility2().getMaxLength() ) - { - continue; - } + if(outConstraint.noLP) { + ///////////////////////////////////////// + // check direct extension to the right of the noLP stacking + ///////////////////////////////////////// + // direct cell access (const) + rightExt = &(hybridZ(i1+noLpShift,i2+noLpShift)); + // check if right side can pair and interaction length is within boundary + if (!Z_equal(rightExt->val, 0.0) + && (rightExt->j1 +1 -i1) <= energy.getAccessibility1().getMaxLength() + && (rightExt->j2 +1 -i2) <= energy.getAccessibility2().getMaxLength() ) + { // compute Z for direct extension with stacking curZ = iStackZ * rightExt->val; From ade93ddb9c56bd40a77ee78638b69c49e7d668b4 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 05:51:36 +0200 Subject: [PATCH 23/27] docs: record heuristic cell corrections --- doc/refactor/1-current-state.md | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/doc/refactor/1-current-state.md b/doc/refactor/1-current-state.md index 358afad1..20316eea 100644 --- a/doc/refactor/1-current-state.md +++ b/doc/refactor/1-current-state.md @@ -3,7 +3,7 @@ This document records the phase-1 audit requested in discussion #232. It is about scientific behavior and runtime cost, not a style rewrite. The audited baseline is IntaRNA 3.4.1, commit `3b14bc0`. The correction and regression -record below is current through `a975c4d`. +record below is current through `96c14b3`. ## Reproducible baseline @@ -24,6 +24,10 @@ After the phase-1 regressions and local corrections through `a975c4d`: - the API binary passes 3,895 assertions in 31 test cases; and - all 20 command-line golden cases pass. +The late heuristic-cell regressions and correction through `96c14b3` raise the +API gate to 3,927 assertions in 33 test cases. The 20-case CLI result above is +the latest recorded CLI gate. + The added CLI cases cover asymmetric target/query accessibility options, seed-free `outMinPu` range filtering, and rejection of partition output from overlapping windows. The API additions cover the local ownership, boundary, @@ -120,6 +124,7 @@ separately because it was a test-data error, not a production defect. | ViennaRNA ensemble temporaries lacked complete scoped ownership | `computeIntraEall` leaked its allocated sequence and fold compound on normal return; `computeES` cleanup was not exception-safe | the API suite exercises span-sensitive ES plus both `getEall1/getEall2` paths; direct leak detection still belongs to sanitizer coverage | `595029f` gives Vienna allocations and fold compounds scoped deleters | | exact noLP seed extension added independent partition factors | the stack weight and remaining left subensemble were added, violating the sum-product recurrence | with one allowed two-pair seed in a `3x3` complementary grid, the exact partition is `exp(2)+exp(3)` | `1ed5819` multiplies the stack and subensemble weights | | empty seed-extension ranges retained prior boundary partitions | early return reset scalar optima but did not clear `Z_partition`, so predictor reuse replayed stale sites | exact and heuristic predictors run on a seeded range and then a singleton range; the second partition and boundary store are zero | `a975c4d` calls `initZ()` on both no-seed paths | +| heuristic cell incumbents leaked, and ensemble noLP direct control was over-scoped | `curCellEtotal` was not reset in `PredictorMfe2dHeuristic`, `PredictorMfeEns2dHeuristic` or `PredictorMfe2dHeuristicSeed`, so a stale strict-tie incumbent could suppress a valid cell; ensemble direct extension was also nested under singleton `noGUend` acceptance, while its empty/too-long guards continued the whole cell and skipped later loop extensions | `08d375e` records five pre-fix failures: unseeded and ensemble energies `-2` instead of `-3`, seeded energy `-2` instead of `-4`, and missing direct/bulged boundaries of `exp(3)` and `exp(4)` | `96c14b3` resets each incumbent per cell, independently guards the ensemble direct extension, and leaves later loop enumeration reachable | Additional high-confidence findings are not treated as tiny local fixes because they change aggregation or recurrence ownership and need the benchmark gates of @@ -149,12 +154,16 @@ output, or partition/filter agreement. Phase 1 now has focused evidence for zero-capacity output, one exact-versus-heuristic noLP ensemble comparison, terminal-GU partition filtering, seed-extension sum-product/reuse behavior, asymmetric target/query options, seed-free range filtering, window-partition -rejection, accessibility splitting, and the span-sensitive ViennaRNA ES path. +rejection, accessibility splitting, the span-sensitive ViennaRNA ES path, and +per-cell state plus direct/bulged noLP continuation in the three affected +heuristic families. Coverage still lacks: - a brute-force small-instance partition oracle and direct `maxED` partition filter regression; +- systematic small-grid differential coverage for heuristic pruning and ties + beyond the five focused counterexamples; - direct exact-MFE and exact-seeded-MFE regressions, plus wider bulged-seed and heuristic seed-extension oracles; - reuse tests for other stateful predictor families; @@ -193,8 +202,12 @@ change. introduced scoped ownership for the ViennaRNA ensemble resources. - `de135fd` exposed seed-extension algebra and reuse defects before `1ed5819` restored multiplication and `a975c4d` cleared stale partitions. -- Established the current 31-case / 3,895-assertion API and 20-case CLI gate at - `a975c4d`. +- Established the intermediate 31-case / 3,895-assertion API and 20-case CLI + gate at `a975c4d`. +- `08d375e` exposed stale cell incumbents in all three affected heuristic + families and the ensemble direct/loop control-flow defects in five focused + sections; `96c14b3` reset the incumbents and corrected the ensemble guard. +- Established the current 33-case / 3,927-assertion API gate at `96c14b3`. The discussion once calls the phase-1 document `refactor-changelog.md`; no such file exists and the same phase otherwise consistently requires From 0c9cbc1357b8ef1807b59fd14dede001cae7dabf Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 06:04:37 +0200 Subject: [PATCH 24/27] test: expose output hub reporting defects --- tests/Makefile.am | 1 + tests/OutputHandlerHub_test.cpp | 51 +++++++++++++++++++++++++++++++++ 2 files changed, 52 insertions(+) create mode 100644 tests/OutputHandlerHub_test.cpp diff --git a/tests/Makefile.am b/tests/Makefile.am index 70e5c4b8..30e84189 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -51,6 +51,7 @@ runApiTests_SOURCES = \ NussinovHandler_test.cpp \ RnaSequence_test.cpp \ OutputStreamHandlerSortedCsv_test.cpp \ + OutputHandlerHub_test.cpp \ OutputHandlerInteractionList_test.cpp \ SeedHandlerExplicit_test.cpp \ SeedHandlerNoBulge_test.cpp \ diff --git a/tests/OutputHandlerHub_test.cpp b/tests/OutputHandlerHub_test.cpp new file mode 100644 index 00000000..508eadfb --- /dev/null +++ b/tests/OutputHandlerHub_test.cpp @@ -0,0 +1,51 @@ +#include "catch.hpp" + +#undef NDEBUG + +#include "IntaRNA/OutputHandlerHub.h" +#include "IntaRNA/RnaSequence.h" + +using namespace IntaRNA; + +namespace { + +class CountingOutputHandler : public OutputHandler { +public: + explicit CountingOutputHandler(const OutputConstraint & constraint) + : OutputHandler(constraint) + {} + + void add(const Interaction &) override + { + ++reportedInteractions; + } +}; + +} // namespace + +TEST_CASE("OutputHandlerHub reports the largest child count", "[OutputHandlerHub]") { + + #include "testEasyLoggingSetup.icc" + + OutputConstraint constraint; + CountingOutputHandler first(constraint); + CountingOutputHandler second(constraint); + OutputHandlerHub hub(constraint, false); + + REQUIRE(hub.reported() == 0); + hub.addOutputHandler(&first); + hub.addOutputHandler(&second); + + RnaSequence sequence("sequence", "GG"); + Interaction interaction(sequence, sequence); + hub.add(interaction); + + REQUIRE(first.reported() == 1); + REQUIRE(second.reported() == 1); + REQUIRE(hub.reported() == 1); + + second.add(interaction); + REQUIRE(first.reported() == 1); + REQUIRE(second.reported() == 2); + REQUIRE(hub.reported() == 2); +} From ce0599b4057679419b9ff5d617e775d4e40ea666 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 06:05:25 +0200 Subject: [PATCH 25/27] fix: align output hub forwarding signature --- src/IntaRNA/OutputHandlerHub.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/IntaRNA/OutputHandlerHub.h b/src/IntaRNA/OutputHandlerHub.h index 1e38dd61..08d13eef 100644 --- a/src/IntaRNA/OutputHandlerHub.h +++ b/src/IntaRNA/OutputHandlerHub.h @@ -204,11 +204,11 @@ addOutputHandler( OutputHandler * handler ) inline void OutputHandlerHub:: -add( const Interaction & inter, const OutputConstraint & outConstraint ) +add( const Interaction & inter ) { // forward to all in list for (auto it=outList.begin(); it!=outList.end(); it++) { - (*it)->add(inter,outConstraint); + (*it)->add(inter); } } From 87f45f118ed5d25c6a2be774ab27d66541205048 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 06:06:47 +0200 Subject: [PATCH 26/27] fix: report maximum output hub count --- src/IntaRNA/OutputHandlerHub.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/IntaRNA/OutputHandlerHub.h b/src/IntaRNA/OutputHandlerHub.h index 08d13eef..41270add 100644 --- a/src/IntaRNA/OutputHandlerHub.h +++ b/src/IntaRNA/OutputHandlerHub.h @@ -222,7 +222,7 @@ reported() const size_t maxReported = 0; // get maximal number of reports among all handlers for (auto it=outList.begin(); it!=outList.end(); it++) { - maxReported = std::min( maxReported, (*it)->reported() ); + maxReported = std::max( maxReported, (*it)->reported() ); } // return maximum return maxReported; From ba5f9eb41229b5aa47339b6ee954858ee3d6b93b Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 15 Aug 2026 06:07:55 +0200 Subject: [PATCH 27/27] docs: record output hub corrections --- doc/refactor/1-current-state.md | 16 +++++++++++----- 1 file changed, 11 insertions(+), 5 deletions(-) diff --git a/doc/refactor/1-current-state.md b/doc/refactor/1-current-state.md index 20316eea..f5aebe4c 100644 --- a/doc/refactor/1-current-state.md +++ b/doc/refactor/1-current-state.md @@ -3,7 +3,7 @@ This document records the phase-1 audit requested in discussion #232. It is about scientific behavior and runtime cost, not a style rewrite. The audited baseline is IntaRNA 3.4.1, commit `3b14bc0`. The correction and regression -record below is current through `96c14b3`. +record below is current through `87f45f1`. ## Reproducible baseline @@ -25,8 +25,9 @@ After the phase-1 regressions and local corrections through `a975c4d`: - all 20 command-line golden cases pass. The late heuristic-cell regressions and correction through `96c14b3` raise the -API gate to 3,927 assertions in 33 test cases. The 20-case CLI result above is -the latest recorded CLI gate. +API gate to 3,927 assertions in 33 test cases. The subsequent output-hub +regression and corrections through `87f45f1` raise it to 3,934 assertions in +34 test cases. The complete 20-case CLI suite still passes. The added CLI cases cover asymmetric target/query accessibility options, seed-free `outMinPu` range filtering, and rejection of partition output from @@ -125,6 +126,7 @@ separately because it was a test-data error, not a production defect. | exact noLP seed extension added independent partition factors | the stack weight and remaining left subensemble were added, violating the sum-product recurrence | with one allowed two-pair seed in a `3x3` complementary grid, the exact partition is `exp(2)+exp(3)` | `1ed5819` multiplies the stack and subensemble weights | | empty seed-extension ranges retained prior boundary partitions | early return reset scalar optima but did not clear `Z_partition`, so predictor reuse replayed stale sites | exact and heuristic predictors run on a seeded range and then a singleton range; the second partition and boundary store are zero | `a975c4d` calls `initZ()` on both no-seed paths | | heuristic cell incumbents leaked, and ensemble noLP direct control was over-scoped | `curCellEtotal` was not reset in `PredictorMfe2dHeuristic`, `PredictorMfeEns2dHeuristic` or `PredictorMfe2dHeuristicSeed`, so a stale strict-tie incumbent could suppress a valid cell; ensemble direct extension was also nested under singleton `noGUend` acceptance, while its empty/too-long guards continued the whole cell and skipped later loop extensions | `08d375e` records five pre-fix failures: unseeded and ensemble energies `-2` instead of `-3`, seeded energy `-2` instead of `-4`, and missing direct/bulged boundaries of `exp(3)` and `exp(4)` | `96c14b3` resets each incumbent per cell, independently guards the ensemble direct extension, and leaves later loop enumeration reachable | +| `OutputHandlerHub` could neither instantiate nor report child counts correctly | its inline `add` definition retained an obsolete two-argument signature, while `reported()` initialized zero and repeatedly took the minimum, so a nonempty hub always reported zero | `0c9cbc1` first fails to compile against the public declaration; after the forwarding repair, a two-child hub reports zero instead of the required maximum two | `ce0599b` aligns and forwards the one-argument `add`; `87f45f1` aggregates child counts with `max` | Additional high-confidence findings are not treated as tiny local fixes because they change aggregation or recurrence ownership and need the benchmark gates of @@ -156,7 +158,7 @@ terminal-GU partition filtering, seed-extension sum-product/reuse behavior, asymmetric target/query options, seed-free range filtering, window-partition rejection, accessibility splitting, the span-sensitive ViennaRNA ES path, and per-cell state plus direct/bulged noLP continuation in the three affected -heuristic families. +heuristic families, and output-hub forwarding/count aggregation. Coverage still lacks: @@ -207,7 +209,11 @@ change. - `08d375e` exposed stale cell incumbents in all three affected heuristic families and the ensemble direct/loop control-flow defects in five focused sections; `96c14b3` reset the incumbents and corrected the ensemble guard. -- Established the current 33-case / 3,927-assertion API gate at `96c14b3`. +- `0c9cbc1` exposed the unusable output-hub forwarder and zero report count; + `ce0599b` aligned its public forwarding call and `87f45f1` returns the + documented maximum child count. +- Established the current 34-case / 3,934-assertion API and 20-case CLI gate at + `87f45f1`. The discussion once calls the phase-1 document `refactor-changelog.md`; no such file exists and the same phase otherwise consistently requires