From 89ab7051ac3c09beb353ac4fd8baf85024e71a0c Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Sat, 3 Oct 2026 01:59:05 +0200 Subject: [PATCH 1/6] Add deterministic kinetic seed-extension predictor --- ChangeLog | 22 + README.md | 33 ++ doc/Makefile.am | 2 +- doc/kinetic-seed-extension.md | 200 ++++++++ src/IntaRNA/Makefile.am | 2 + src/IntaRNA/PredictorSeedExtensionKinetic.cpp | 435 +++++++++++++++++ src/IntaRNA/PredictorSeedExtensionKinetic.h | 152 ++++++ src/bin/CommandLineParsing.cpp | 35 +- src/bin/CommandLineParsing.h | 2 + tests/Makefile.am | 3 +- tests/PredictorSeedExtensionKinetic_test.cpp | 458 ++++++++++++++++++ tests/runKineticSeedExtension.sh | 119 +++++ 12 files changed, 1458 insertions(+), 5 deletions(-) create mode 100644 doc/kinetic-seed-extension.md create mode 100644 src/IntaRNA/PredictorSeedExtensionKinetic.cpp create mode 100644 src/IntaRNA/PredictorSeedExtensionKinetic.h create mode 100644 tests/PredictorSeedExtensionKinetic_test.cpp create mode 100755 tests/runKineticSeedExtension.sh diff --git a/ChangeLog b/ChangeLog index 62b5f84c..464319e5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,10 @@ ## Interface and handling +- deterministic kinetic seed extension with --model=X --mode=K and configurable + local scoring; complete-energy downhill steps, atomic noLP loop/stack moves + and retained greedy traceback; equilibrium probability output is unsupported + - IntaRNAeval / --rri evaluates predefined RNA-RNA interactions (issue #184) - compressed binary .agz accessibility caches for repeated screens (issue #245) @@ -73,6 +77,24 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261003 Alexander Mitrofanov + * IntaRNA/PredictorSeedExtensionKinetic, src/IntaRNA/Makefile.am : + + greedily extend seeds on either side using complete interaction-energy + differences and deterministic thermodynamic or distance-weighted scores + + support noLP macro-steps, GU restrictions, per-strand loop/span constraints, + explicit-seed validation and cached trajectory traceback + + retain valid visited prefixes and select non-overlapping output from the + complete retained candidate set; reject unsupported ensemble statistics + * bin/CommandLineParsing : + + expose --mode=K exclusively for --model=X and --kineticScore=A|B|C + * tests/PredictorSeedExtensionKinetic_test.cpp, tests/runKineticSeedExtension.sh, + tests/Makefile.am : + + validate local move choices, macro-step energetics, structural constraints, + traceback, output filtering and CLI compatibility + * README.md, doc/kinetic-seed-extension.md, doc/Makefile.am : + + document the council's implementation consensus, scientific scope and + replacement of unsafe loop-only pruning with exact move enumeration + 261001 Alexander Mitrofanov * IntaRNA/PredictorEvalOnly, src/IntaRNA/Makefile.am : + parse colon-separated hybridDB structures with sequence/index/pair validation diff --git a/README.md b/README.md index 746ed498..cc0c72b3 100644 --- a/README.md +++ b/README.md @@ -729,6 +729,39 @@ minimum free energy interaction. Putative seed interactions (used by the `H` and `M` mode) can be enumerated and studied using the `S` mode. +### Greedy kinetic seed extension + +`--model=X --mode=K` grows each available seed along a deterministic greedy +path. Every step compares feasible extensions on both sides using the complete +change in interaction energy, including accessibility, terminal penalties and +dangling ends. Only strictly negative changes are accepted. With `--outNoLP`, +crossing a loop forms its closing pair and the immediately following stack as +one atomic step; a favorable stack can therefore compensate for an unfavorable +loop. Initial seeds must satisfy the selected structural constraints. + +`--kineticScore` selects the local move ranking: + +| Value | Score minimized for a move with gaps `s1`, `s2` | +| --- | --- | +| `A` (default) | Complete energy change | +| `B` | Complete energy change / `(1+s1+s2)` | +| `C` (C1 in the design) | Complete energy change / `(1+2*max(s1,s2))` | + +Equal scores prefer the left side, then fewer unpaired bases, then smaller +`s1`. The denominators also apply to two-pair moves. All reportable visited +states, including seeds, participate in the normal energy-ranked output; +traceback preserves the actual chosen path. `--outNoGUend`, separate query and +target loop/span limits, regions, output energy/accessibility filters and +overlap settings remain applicable. + +This mode is a zippering-inspired heuristic, without a calibrated time axis or +a guarantee of the global minimum. It evaluates all feasible local moves; +loop-only energetic pruning is unsafe for complete loop-plus-stack steps. +Equilibrium probability/partition-sum outputs are rejected, as are other models +and `--noSeed`. Scores B and C are optional distance preferences, not measured +kinetic rates. See the [design and implementation plan](doc/kinetic-seed-extension.md) +for the precise algorithm, scientific limitations and validation cases. + [![up](doc/figures/icon-up.28.png) back to overview](#overview)

diff --git a/doc/Makefile.am b/doc/Makefile.am index 322e9d93..003a02d3 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -5,6 +5,7 @@ EXTRA_DIST = \ conda.txt \ + kinetic-seed-extension.md \ doxygen.cfg \ latex-deps/adjcalc.sty \ latex-deps/adjustbox.sty \ @@ -13,4 +14,3 @@ EXTRA_DIST = \ latex-deps/tocloft.sty \ latex-deps/trimclip.sty \ latex-deps/xtab.sty - diff --git a/doc/kinetic-seed-extension.md b/doc/kinetic-seed-extension.md new file mode 100644 index 00000000..c02acf79 --- /dev/null +++ b/doc/kinetic-seed-extension.md @@ -0,0 +1,200 @@ +# Deterministic kinetic seed extension + +## Council decision and scientific scope + +`PredictorSeedExtensionKinetic`, selected exclusively with `--model=X --mode=K`, +implements a deterministic, seed-conditioned, downhill extension heuristic. +It starts from each seed provided by the selected seed handler, compares moves +at both ends of the current duplex, and commits one move at a time. It does not +sample transition rates, simulate elapsed time, cross uphill barriers, remove +base pairs, or guarantee the global minimum free energy interaction. A +macro-step is a coarse move of this heuristic: favorable combined energy does +not establish a barrier-free physical reaction pathway. + +This scope follows the distinction between gradient walks to local minima and +kinetic dynamics in the ViennaRNA ecosystem. [RNAlocmin](https://www.tbi.univie.ac.at/RNA/BHG/RNAlocmin.html) +uses gradient walks; [Kinfold and treekin](https://www.tbi.univie.ac.at/software/) +use stochastic moves or transition-rate matrices. The name "kinetic" identifies +the proposed extension strategy, not a calibrated kinetic prediction. The +choice of a single seed structure per seed start inherits the selected seed +handler's behavior; this algorithm does not enumerate every seed conformation. + +The second council resolved the first review's open questions as follows. +These decisions supersede contradictory optimization and biological claims in +the supplied `PredSeedExtKinetic.md` proposal. + +## State, energy and accepted moves + +A state stores the complete, ordered chain of intermolecular base pairs, its +inclusive boundaries `(i1,j1,i2,j2)`, its hybridization energy `H`, and its full +interaction energy `E`. Indices use IntaRNA's internal coordinate system: +sequence 2 is reversed and prediction-range offsets are handled by the existing +wrappers. The initial `H` includes the seed's loop energies and `getE_init()` +exactly once. Every accepted step retains its actual pairs for traceback. + +Evaluate every candidate with the active `InteractionEnergy` instance: + +``` +E_current = energy.getE(i1, j1, i2, j2, H_current) +E_next = energy.getE(i1_next, j1_next, i2_next, j2_next, H_next) +delta = E_next - E_current +``` + +This includes both accessibility penalties, both accessibility-weighted dangling +ends, terminal penalties, the configured temperature and parameters, and the +configured additive energy term. Changing one boundary can alter the opposite +end's dangling contribution through its accessibility weight. Consequently, +loop-plus-stack-plus-accessibility differences alone are insufficient. Infinite +states are rejected before subtraction. A constant additive term cancels in the +difference while remaining part of reported energy. + +Only moves with **strictly negative full energy difference** are eligible. +Zero and uphill moves are rejected. Stop when no eligible move remains; the +finite, growing span also guarantees termination. Output energy and accessibility +thresholds filter reported states, rather than introducing additional barriers +in the trajectory. Model-inaccessible states and maximum span limits still +make a move infeasible. + +For each side, enumerate all `0 <= s1 <= m1` and `0 <= s2 <= m2`, where `mk` +is the active energy model's maximum unpaired loop size for strand `k`. +`(s1,s2)=(0,0)` is a stack. Other combinations include bulges and internal loops. +The maximum total skipped length is `m1+m2`, equal to `2*m` only when the +per-strand limits are equal. The proposal's fixed 10 is not a separate limit. + +A normal move adds the complementary loop-closing pair and advances each +boundary by `sk+1`. With `--outNoLP`, a move across any nonzero loop additionally +requires the immediately following outward stack. Evaluate and commit these +two pairs **atomically**, advancing each boundary by `sk+2`. Neither the +isolated closing pair nor its energy is a separately accepted or reported state. +Check both pairs, the intervening loop, and both complete strand spans before +acceptance. Each span must fit that strand's accessibility maximum length and +the requested prediction range; no single ambiguous shared `W` is introduced. + +The complete initial seed must have ordered complementary pairs and finite +per-loop energies. Under `--outNoLP`, every seed pair must already have a direct +stack neighbor. Incompatible explicit seeds are skipped, including seeds with +lonely terminal pairs; this mode does not perform a preliminary seed repair. +Together with atomic macro-steps, this preserves the no-lonely-pair invariant. + +Under `--outNoGUend`, both boundaries of every nonstacking extension loop must +be non-GU. Direct stacking may temporarily expose a GU outer endpoint, as in +existing IntaRNA extension recurrences. A state is reportable only if both outer +endpoints satisfy the flag. The active energy model's internal-loop GU policy +also remains authoritative. A valid earlier state remains available if descent +later stops at an unreportable GU endpoint. + +## Scores and deterministic choice + +`--kineticScore=A|B|C` selects the score; `C` denotes the proposal's C1: + +| Option | Score minimized | Interpretation | +| --- | --- | --- | +| `A` (default) | `delta` | Steepest decrease of the actual modeled interaction energy | +| `B` | `delta/(1+s1+s2)` | Heuristic preference per total skipped length | +| `C` | `delta/(1+2*max(s1,s2))` | Heuristic preference penalizing the longer skipped strand | + +A is the default because it requires no uncalibrated length-to-time assumption. +B and C are retained as explicit alternatives for exploring the supplied +proposal; their denominators are not experimentally calibrated rates or times. +The written `+1` denominator is retained even for a two-pair macro-step: it +counts one candidate move, not the number of pairs added. Scores choose the next +move only. Reported interactions remain ranked by their modeled total energy `E`. + +Compare candidates from **both** sides together. Equal scores are resolved by +left before right, then smaller `s1+s2`, then smaller `s1`. Use sufficiently wide +arithmetic for score comparison, avoiding integer division truncation. Seed +iteration and output comparison follow deterministic existing IntaRNA ordering. + +## Why the proposed pruning is removed + +A loop-only lower bound is not a lower bound for a loop-plus-stack macro-step: +a positive loop may be rescued by the negative following stack. The initial +review reproduced a Turner2004 example at 37 degrees Celsius with a `+0.50` +kcal/mol loop and `-3.30` kcal/mol following stack, totaling `-2.80` kcal/mol. +The full delta also contains terminal and dangling changes absent from the +proposal's filters. Accessibility at the farthest candidate endpoint is not a +certified lower bound for nearer candidates; imported accessibility values need +not obey monotonicity assumptions. + +Therefore the initial implementation exhaustively evaluates all feasible moves +for the current state. It uses neither the six-base-pair-type precomputed Turner +tables nor their `mdspan`/suffix-minimum representation. The suffix-minimum +operation itself is valid, but cannot repair an invalid underlying bound. +Populating an exact full-delta table first would add storage and selection work +without avoiding those energy evaluations. Omitting that cosmetic optimization +is deliberate, rather than replacing one unsafe bound with another. + +There are at most `2*(m1+1)*(m2+1)` candidate shapes per state, with constant-size +local energy updates plus full boundary-energy evaluation for each. No runtime +improvement over other predictors is claimed. Any future pruning must supply a +certified lower bound for the **complete** move under the active energy model +and pass differential tests against this exhaustive implementation. + +## Reporting, traceback and compatibility + +The seed and every atomically committed, reportable prefix are candidates for +normal IntaRNA output. All earlier valid prefixes remain available for suboptimal +and overlap-constrained reporting. Keeping only the last state would lose valid +GU-end prefixes; keeping only the best state for each left boundary would lose +shorter candidates needed after excluding overlap with another report. + +Use existing MFE ordering, output filters, seed annotations and index conversion. +Cache the actual selected pair chain for each retained candidate, including all +macro-step pairs. Traceback must recover that chain rather than run an unrelated +minimum-energy recurrence between its boundaries. Resolve duplicate boundaries +by retaining the lowest total energy, then the lexicographically smallest full +base-pair chain on an energy tie. Reduce duplicates before feeding the ordinary +optimum collector, so no stale energy can select a different cached path. Seed +annotations include only starting seeds that passed validation. Output energies +must agree with an +independent sum over the reported chain and complete boundary terms. + +The cache stores full paths: its memory cost is proportional to the sum of the +retained paths' lengths, in addition to the seed handler's storage. Retaining +prefixes is therefore not a constant-memory walk. It supports reproducible +traceback and shorter alternatives for overlap-constrained output. + +Only `model=X` accepts `mode=K`; seedless operation and `kineticScore` outside K +are rejected. Existing defaults and other models remain unchanged. Ordinary +energy and minimum-energy tracker output are supported. Equilibrium partition +functions and normalized equilibrium probabilities are not defined by this +selected collection of greedy trajectories. Requests needing `Zall`, including +ensemble output and probability trackers, are rejected in CLI validation, with +an API constructor guard for `needZall`. The algorithm does not manufacture an +ensemble by summing repeated prefixes from multiple seeds. + +## Implementation plan and acceptance checklist + +1. Add the public predictor header and implementation, derived from + `PredictorMfe`, with owned seed handler, the existing index-offset wrappers, + validated A/B/C selection and `needZall` rejection. Register both files for + building and installation. +2. Initialize and trace each handler-provided seed; validate its full chain, + seed range, energy, strand spans and active structural constraints. Record + reportable seed states. +3. Enumerate both sides and all feasible loop shapes, form complete normal or + atomic noLP moves, recompute full candidate energy, and select the strictly + downhill winner using the chosen score and the specified tie order. Repeat + until stalling or range exhaustion. +4. Preserve complete committed paths and valid prefixes; integrate ordinary MFE + output and custom candidate lookup for overlap-constrained suboptimals. + Reconstruct seed annotations without changing the retained greedy chain. +5. Wire `--model=X --mode=K` and `--kineticScore=A|B|C` into parsing, help and + factory construction. Reject incompatible model, seedless and ensemble or + probability requests with clear diagnostics. Update README and ChangeLog. +6. Add an independent tiny-sequence reference that enumerates absolute candidate + endpoints and recomputes the entire chain energy. Compare reported energy, + coordinates and traceback against it across scores, constraints and offsets. + Include targeted regressions for a positive-loop/negative-stack rescue, + strict stopping, tie order, nonmonotone accessibility, full boundary-energy + changes, invalid explicit seeds and GU-prefix retention. +7. Run focused API tests, CLI mode/flag compatibility checks, full `make tests`, + debug validation of bounds and ownership, installed standalone-header checks + where supported, and `git diff --check`. Inspect failures; never regenerate + expected outputs merely to hide a change. +8. Review the complete diff and open a pull request documenting this scientific + scope, the deliberate pruning correction, implemented behavior, validation + results and any remaining environmental limitations. + +Test and build results belong in the pull request and development record; this +checklist specifies the required work and does not imply an unrun check passed. diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index 98288166..5f27f746 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -72,6 +72,7 @@ libIntaRNA_a_HEADERS = \ PredictorMfe2dSeed.h \ PredictorMfe2dSeedExtension.h \ PredictorMfe2dSeedExtensionRIblast.h \ + PredictorSeedExtensionKinetic.h \ PredictorMfe2dHeuristic.h \ PredictorMfe2dHeuristicSeed.h \ PredictorMfe2dHelixBlockHeuristic.h \ @@ -134,6 +135,7 @@ libIntaRNA_a_SOURCES = \ PredictorMfe2dSeed.cpp \ PredictorMfe2dSeedExtension.cpp \ PredictorMfe2dSeedExtensionRIblast.cpp \ + PredictorSeedExtensionKinetic.cpp \ PredictorMfe2dHeuristic.cpp \ PredictorMfe2dHeuristicSeed.cpp \ PredictorMfe2dHelixBlockHeuristic.cpp \ diff --git a/src/IntaRNA/PredictorSeedExtensionKinetic.cpp b/src/IntaRNA/PredictorSeedExtensionKinetic.cpp new file mode 100644 index 00000000..4b085b93 --- /dev/null +++ b/src/IntaRNA/PredictorSeedExtensionKinetic.cpp @@ -0,0 +1,435 @@ +#include "IntaRNA/PredictorSeedExtensionKinetic.h" + +#include +#include +#include +#include + +#include + +namespace IntaRNA { +namespace { + +SeedHandler * checkedSeedHandler(SeedHandler * handler) +{ + if (handler == NULL) { + throw std::invalid_argument("PredictorSeedExtensionKinetic requires a seed handler"); + } + return handler; +} + +// Do not add infinity sentinels or overflow the internal integer energy type. +E_type addEnergy(const E_type first, const E_type second) +{ + if (E_isINF(first) || E_isINF(second)) { + return E_INF; + } + const std::int64_t sum = std::int64_t(first) + std::int64_t(second); + return sum >= E_INF || sum < std::numeric_limits::min() + ? E_INF : static_cast(sum); +} + +} // namespace + +////////////////////////////////////////////////////////////////////////// + +PredictorSeedExtensionKinetic::PredictorSeedExtensionKinetic( + const InteractionEnergy & energy, OutputHandler & output, + PredictionTracker * predTracker, SeedHandler * seedHandlerInstance, + const char score) + : PredictorMfe(energy, output, predTracker) + , seedHandler(checkedSeedHandler(seedHandlerInstance)) + , score(score) + , interactions() + , validSeeds() +{ + if (score != 'A' && score != 'B' && score != 'C') { + throw std::invalid_argument("PredictorSeedExtensionKinetic score must be A, B or C"); + } + if (output.getOutputConstraint().needZall) { + throw std::invalid_argument("PredictorSeedExtensionKinetic does not compute an equilibrium partition function"); + } +} + +////////////////////////////////////////////////////////////////////////// + +PredictorSeedExtensionKinetic::~PredictorSeedExtensionKinetic() +{ +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::predict(const IndexRange & r1, const IndexRange & r2) +{ + const size_t size1 = energy.getAccessibility1().getSequence().size(); + const size_t size2 = energy.getAccessibility2().getSequence().size(); + if (!r1.isAscending() || !r2.isAscending() + || r1.from >= size1 || r2.from >= size2) { + throw std::invalid_argument("PredictorSeedExtensionKinetic::predict(): invalid sequence range"); + } + + energy.setOffset1(r1.from); + energy.setOffset2(r2.from); + seedHandler.setOffset1(r1.from); + seedHandler.setOffset2(r2.from); + const size_t last1 = std::min(r1.to, size1 - 1) - r1.from; + const size_t last2 = std::min(r2.to, size2 - 1) - r2.from; + interactions.clear(); + validSeeds.clear(); + initOptima(); + + if (seedHandler.fillSeed(0, last1, 0, last2) != 0) { + size_t i1 = RnaSequence::lastPos, i2 = RnaSequence::lastPos; + while (seedHandler.updateToNextSeed(i1, i2, 0, last1, 0, last2)) { + const size_t length1 = seedHandler.getSeedLength1(i1, i2); + const size_t length2 = seedHandler.getSeedLength2(i1, i2); + if (length1 == 0 || length2 == 0 + || length1 - 1 > last1 - i1 || length2 - 1 > last2 - i2 + || length1 > energy.getAccessibility1().getMaxLength() + || length2 > energy.getAccessibility2().getMaxLength()) { + continue; + } + const size_t j1 = i1 + length1 - 1; + const size_t j2 = i2 + length2 - 1; + if (E_isINF(seedHandler.getSeedE(i1, i2))) { + continue; + } + + Interaction interaction(energy.getAccessibility1().getSequence(), + energy.getAccessibility2().getAccessibilityOrigin().getSequence()); + interaction.basePairs.push_back(energy.getBasePair(i1, i2)); + seedHandler.traceBackSeed(interaction, i1, i2); + if (i1 != j1 || i2 != j2) { + interaction.basePairs.push_back(energy.getBasePair(j1, j2)); + } + interaction.sort(); + E_type hybrid = E_INF; + if (!isValidSeed(interaction, hybrid) + || getBoundary(interaction) != Boundary{i1, j1, i2, j2}) { + continue; + } + interaction.energy = energy.getE(i1, j1, i2, j2, hybrid); + if (E_isINF(interaction.energy)) { + continue; + } + interaction.setSeedRange(interaction.basePairs.front(), + interaction.basePairs.back(), interaction.energy); + validSeeds.emplace(interaction.basePairs.front(), + Interaction::Seed(interaction.basePairs.front(), + interaction.basePairs.back(), interaction.energy)); + extendSeed(interaction, hybrid, last1, last2); + } + } + + // Reduce identical boundaries before updating optima: an earlier, inferior + // path must never be traced back using a later replacement's base pairs. + for (const auto & entry : interactions) { + const Boundary & b = entry.first; + updateOptima(b[0], b[1], b[2], b[3], entry.second.energy, false, false); + } + // The generic reporter assumes a nonempty optimum list. Zero reports can + // still be useful for prediction trackers and must not dereference it. + if (output.getOutputConstraint().reportMax != 0) { + reportOptima(); + } +} + +////////////////////////////////////////////////////////////////////////// + +PredictorSeedExtensionKinetic::Boundary +PredictorSeedExtensionKinetic::getBoundary(const Interaction & interaction) const +{ + return Boundary{energy.getIndex1(interaction.basePairs.front()), + energy.getIndex1(interaction.basePairs.back()), + energy.getIndex2(interaction.basePairs.front()), + energy.getIndex2(interaction.basePairs.back())}; +} + +////////////////////////////////////////////////////////////////////////// + +bool +PredictorSeedExtensionKinetic::isValidSeed(const Interaction & interaction, E_type & hybrid) const +{ + if (interaction.basePairs.empty() || !interaction.isValid()) { + return false; + } + const auto & pairs = interaction.basePairs; + const auto & constraint = output.getOutputConstraint(); + // Reconstruct rather than trust a cached explicit-seed energy: the + // trajectory's energy must correspond to precisely the traced structure. + hybrid = energy.getE_init(); + if (E_isINF(hybrid)) { + return false; + } + for (size_t p = 0; p < pairs.size(); ++p) { + const size_t i1 = energy.getIndex1(pairs[p]); + const size_t i2 = energy.getIndex2(pairs[p]); + if (i1 >= energy.size1() || i2 >= energy.size2() + || !energy.areComplementary(i1, i2)) { + return false; + } + const bool stackedLeft = p > 0 + && pairs[p].first - pairs[p-1].first == 1 + && pairs[p-1].second - pairs[p].second == 1; + const bool stackedRight = p + 1 < pairs.size() + && pairs[p+1].first - pairs[p].first == 1 + && pairs[p].second - pairs[p+1].second == 1; + if (constraint.noLP && !stackedLeft && !stackedRight) { + return false; + } + if (p == 0) { + continue; + } + const size_t previous1 = energy.getIndex1(pairs[p-1]); + const size_t previous2 = energy.getIndex2(pairs[p-1]); + if (!stackedLeft && constraint.noGUend + && (energy.isGU(previous1, previous2) || energy.isGU(i1, i2))) { + return false; + } + hybrid = addEnergy(hybrid, energy.getE_interLeft(previous1, i1, previous2, i2)); + if (E_isINF(hybrid)) { + return false; + } + } + return true; +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::extendSeed(Interaction & interaction, + E_type hybrid, const size_t last1, const size_t last2) +{ + const size_t maxLength1 = energy.getAccessibility1().getMaxLength(); + const size_t maxLength2 = energy.getAccessibility2().getMaxLength(); + while (true) { + retain(interaction); + const Boundary bounds = getBoundary(interaction); + const size_t remaining1 = maxLength1 - (bounds[1] - bounds[0] + 1); + const size_t remaining2 = maxLength2 - (bounds[3] - bounds[2] + 1); + bool found = false; + Candidate best; + for (unsigned int side = 0; side < 2; ++side) { + const bool left = side == 0; + const size_t space1 = std::min(remaining1, left ? bounds[0] : last1 - bounds[1]); + const size_t space2 = std::min(remaining2, left ? bounds[2] : last2 - bounds[3]); + if (space1 == 0 || space2 == 0) { + continue; + } + // A GU boundary can still grow by stacking, but cannot start a + // nonstacking loop when either active constraint forbids it. + const bool stackOnly = (output.getOutputConstraint().noGUend + || !energy.isInternalLoopGUallowed()) + && energy.isGU(bounds[left ? 0 : 1], bounds[left ? 2 : 3]); + const size_t maxGap1 = stackOnly ? 0 + : std::min(energy.getMaxInternalLoopSize1(), space1 - 1); + const size_t maxGap2 = stackOnly ? 0 + : std::min(energy.getMaxInternalLoopSize2(), space2 - 1); + for (size_t s1 = 0; s1 <= maxGap1; ++s1) { + for (size_t s2 = 0; s2 <= maxGap2; ++s2) { + Candidate candidate; + candidate.left = left; + candidate.s1 = s1; + candidate.s2 = s2; + candidate.macro = output.getOutputConstraint().noLP && (s1 != 0 || s2 != 0); + const size_t addedPairs = candidate.macro ? 2 : 1; + if (space1 < addedPairs || space2 < addedPairs + || s1 > space1 - addedPairs || s2 > space2 - addedPairs) { + continue; + } + candidate.bounds = bounds; + if (left) { + candidate.bounds[0] -= s1 + addedPairs; + candidate.bounds[2] -= s2 + addedPairs; + candidate.close1 = bounds[0] - s1 - 1; + candidate.close2 = bounds[2] - s2 - 1; + } else { + candidate.bounds[1] += s1 + addedPairs; + candidate.bounds[3] += s2 + addedPairs; + candidate.close1 = bounds[1] + s1 + 1; + candidate.close2 = bounds[3] + s2 + 1; + } + if (evaluate(candidate, bounds, hybrid, interaction.energy) + && (!found || isBetter(candidate, best))) { + best = candidate; + found = true; + } + } + } + } + if (!found) { + break; + } + const Interaction::BasePair close = energy.getBasePair(best.close1, best.close2); + if (best.left) { + interaction.basePairs.insert(interaction.basePairs.begin(), close); + if (best.macro) { + interaction.basePairs.insert(interaction.basePairs.begin(), + energy.getBasePair(best.bounds[0], best.bounds[2])); + } + } else { + interaction.basePairs.push_back(close); + if (best.macro) { + interaction.basePairs.push_back(energy.getBasePair(best.bounds[1], best.bounds[3])); + } + } + hybrid = best.hybrid; + interaction.energy = best.total; + } +} + +////////////////////////////////////////////////////////////////////////// + +bool +PredictorSeedExtensionKinetic::evaluate(Candidate & candidate, + const Boundary & bounds, const E_type hybrid, const E_type total) const +{ + const size_t old1 = bounds[candidate.left ? 0 : 1]; + const size_t old2 = bounds[candidate.left ? 2 : 3]; + const size_t outer1 = candidate.bounds[candidate.left ? 0 : 1]; + const size_t outer2 = candidate.bounds[candidate.left ? 2 : 3]; + if (!energy.areComplementary(candidate.close1, candidate.close2) + || (candidate.macro && !energy.areComplementary(outer1, outer2))) { + return false; + } + if (output.getOutputConstraint().noGUend && (candidate.s1 != 0 || candidate.s2 != 0) + && (energy.isGU(old1, old2) || energy.isGU(candidate.close1, candidate.close2))) { + return false; + } + const E_type loop = candidate.left + ? energy.getE_interLeft(candidate.close1, old1, candidate.close2, old2) + : energy.getE_interLeft(old1, candidate.close1, old2, candidate.close2); + candidate.hybrid = addEnergy(hybrid, loop); + if (candidate.macro && E_isNotINF(candidate.hybrid)) { + const E_type stack = candidate.left + ? energy.getE_interLeft(outer1, candidate.close1, outer2, candidate.close2) + : energy.getE_interLeft(candidate.close1, outer1, candidate.close2, outer2); + candidate.hybrid = addEnergy(candidate.hybrid, stack); + } + if (E_isINF(candidate.hybrid)) { + return false; + } + const Boundary & b = candidate.bounds; + candidate.total = energy.getE(b[0], b[1], b[2], b[3], candidate.hybrid); + if (E_isINF(candidate.total)) { + return false; + } + candidate.delta = std::int64_t(candidate.total) - std::int64_t(total); + return candidate.delta < 0; +} + +////////////////////////////////////////////////////////////////////////// + +bool +PredictorSeedExtensionKinetic::isBetter(const Candidate & candidate, const Candidate & best) const +{ + // A 32-bit energy difference times a 65-bit gap denominator fits in + // 128 bits, including for public-API loop limits beyond the CLI limits. + using Wide = boost::multiprecision::int128_t; + const auto denominator = [this](const Candidate & c) -> Wide { + if (score == 'B') { + return Wide(1) + Wide(c.s1) + Wide(c.s2); + } + if (score == 'C') { + return Wide(1) + 2 * Wide(std::max(c.s1, c.s2)); + } + return Wide(1); + }; + const Wide lhs = Wide(candidate.delta) * denominator(best); + const Wide rhs = Wide(best.delta) * denominator(candidate); + if (lhs != rhs) { + return lhs < rhs; + } + if (candidate.left != best.left) { + return candidate.left; + } + const Wide size = Wide(candidate.s1) + Wide(candidate.s2); + const Wide bestSize = Wide(best.s1) + Wide(best.s2); + return size != bestSize ? size < bestSize : candidate.s1 < best.s1; +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::retain(const Interaction & interaction) +{ + const Boundary b = getBoundary(interaction); + const auto & constraint = output.getOutputConstraint(); + if (interaction.energy >= E_MAX + || (constraint.noGUend && (energy.isGU(b[0], b[2]) || energy.isGU(b[1], b[3]))) + || energy.getED1(b[0], b[1]) > constraint.maxED + || energy.getED2(b[2], b[3]) > constraint.maxED) { + return; + } + auto existing = interactions.find(b); + if (existing == interactions.end()) { + interactions.emplace(b, interaction); + } else if (interaction.energy < existing->second.energy + || (interaction.energy == existing->second.energy + && interaction.basePairs < existing->second.basePairs)) { + existing->second = interaction; + } +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::traceBack(Interaction & interaction) +{ + if (interaction.basePairs.empty()) { + return; + } + const auto path = interactions.find(getBoundary(interaction)); + if (path == interactions.end() || path->second.energy != interaction.energy) { + throw std::runtime_error("PredictorSeedExtensionKinetic::traceBack(): no matching greedy path"); + } + interaction = path->second; + seedHandler.addSeeds(interaction); + // The generic annotator can recognize explicit seeds that were rejected + // as starting states (e.g. a lonely seed end stacked only by extension). + // Retain only validated starts and use their reconstructed energies. + if (interaction.seed != NULL) { + Interaction::SeedSet validAnnotations; + for (const Interaction::Seed & seed : *interaction.seed) { + const auto valid = validSeeds.find(seed.bp_i); + if (valid != validSeeds.end() && valid->second.bp_j == seed.bp_j) { + validAnnotations.insert(valid->second); + } + } + *interaction.seed = std::move(validAnnotations); + } +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::getNextBest(Interaction & interaction) +{ + const E_type previousEnergy = interaction.energy; + const Interaction * best = NULL; + for (const auto & entry : interactions) { + const Boundary & b = entry.first; + const Interaction & candidate = entry.second; + if (candidate.energy < previousEnergy + || reportedInteractions.first.overlaps(IndexRange(b[0], b[1])) + || reportedInteractions.second.overlaps(IndexRange(b[2], b[3]))) { + continue; + } + if (best == NULL || candidate < *best) { + best = &candidate; + } + } + if (best == NULL) { + interaction.clear(); + interaction.energy = E_INF; + return; + } + interaction = *best; + const Interaction::BasePair right = interaction.basePairs.back(); + interaction.basePairs.resize(interaction.basePairs.size() == 1 ? 1 : 2); + interaction.basePairs.back() = right; + INTARNA_CLEANUP(interaction.seed); +} + +} // namespace IntaRNA diff --git a/src/IntaRNA/PredictorSeedExtensionKinetic.h b/src/IntaRNA/PredictorSeedExtensionKinetic.h new file mode 100644 index 00000000..e9c3f3d7 --- /dev/null +++ b/src/IntaRNA/PredictorSeedExtensionKinetic.h @@ -0,0 +1,152 @@ +#ifndef INTARNA_PREDICTORSEEDEXTENSIONKINETIC_H_ +#define INTARNA_PREDICTORSEEDEXTENSIONKINETIC_H_ + +#include "IntaRNA/PredictorMfe.h" +#include "IntaRNA/SeedHandlerIdxOffset.h" + +#include +#include +#include + +namespace IntaRNA { + +/** + * Deterministic, strictly downhill extension of every feasible seed. + * + * At each step both ends compete using the complete interaction-energy + * difference, including accessibility, terminal penalties and both weighted + * dangling ends. Score A uses this difference directly; B divides by + * 1+s1+s2; C (C1) divides by 1+2*max(s1,s2). These are heuristic move rankings, + * not physical rates or a calibrated folding-time model. Only negative + * energy differences are accepted, including for the normalized scores. + * + * With noLP, every starting seed must already contain no lonely pair and a + * nonstacking extension adds its closing pair and the following stack + * atomically. Ties prefer left, smaller s1+s2, then smaller s1. Every valid + * visited prefix is eligible for normal MFE/suboptimal reporting; traceback + * reproduces the actual greedy path. Equilibrium partition-function output + * is unsupported. + * + * Candidate enumeration uses the active energy model and separate loop/span + * limits for both RNAs. It deliberately has no loop-only energy pruning: + * such bounds omit favorable mandatory stacks and changes to the opposite + * dangling end, and accessibility differences need not be monotone. + */ +class PredictorSeedExtensionKinetic : public PredictorMfe { +public: + + /** + * Constructs a predictor, taking ownership of tracker and seed handler. + * @param energy energy model, which must outlive this predictor + * @param output output handler, which must outlive this predictor + * @param predTracker owned tracker, or NULL + * @param seedHandler owned, non-NULL seed handler + * @param score move ranking: A, B, or C (the C1 formula) + * @throws std::invalid_argument for a NULL seed handler, unknown score or + * output requiring an equilibrium partition function + */ + PredictorSeedExtensionKinetic(const InteractionEnergy & energy, + OutputHandler & output, PredictionTracker * predTracker, + SeedHandler * seedHandler, const char score = 'A'); + + /** Frees the owned seed handler and prediction tracker. */ + virtual ~PredictorSeedExtensionKinetic(); + + /** + * Extends all feasible seeds within inclusive, zero-based sequence ranges. + * Sequence 2 uses the reversed indexing of the energy model. Repeated + * calls reset all trajectories, cached output and index offsets. + * @param r1 permitted range in sequence 1 + * @param r2 permitted range in reversed sequence 2 + * @throws std::invalid_argument for an invalid or empty input range + */ + void predict(const IndexRange & r1 = IndexRange(0, RnaSequence::lastPos), + const IndexRange & r2 = IndexRange(0, RnaSequence::lastPos)) override; + +protected: + /** + * Restores the exact stored greedy path and annotates its contained seeds. + * @param interaction interaction boundaries to expand + */ + void traceBack(Interaction & interaction) override; + + /** + * Finds the best cached prefix disjoint from already reported intervals. + * @param interaction current report, replaced by the next report or E_INF + */ + void getNextBest(Interaction & interaction) override; + +private: + //! Inclusive boundaries (i1,j1,i2,j2), using local energy indices. + using Boundary = std::array; + //! Best actual path for each visited, reportable set of boundaries. + using InteractionCache = std::map; + + /** A feasible move; delta is widened before subtraction. */ + struct Candidate { + Boundary bounds = {}; + size_t close1 = 0; + size_t close2 = 0; + size_t s1 = 0; + size_t s2 = 0; + bool left = true; + bool macro = false; + E_type hybrid = E_INF; + E_type total = E_INF; + std::int64_t delta = 0; + }; + + //! Owned seed handler with offsets matching this->energy. + SeedHandlerIdxOffset seedHandler; + //! The selected deterministic move-ranking formula. + const char score; + //! Paths retained independently of the number of requested reports. + InteractionCache interactions; + //! Valid starting seeds with recomputed energies, keyed by original indices. + std::map validSeeds; + + /** @return local, inclusive boundaries of a nonempty interaction */ + Boundary getBoundary(const Interaction & interaction) const; + + /** + * Checks the complete seed path, including noLP and loop GU constraints. + * @param interaction seed path, in original sequence coordinates + * @param hybrid receives the traced path's loop energies plus initiation + * @return whether every pair and adjacent loop is feasible + */ + bool isValidSeed(const Interaction & interaction, E_type & hybrid) const; + + /** + * Runs a complete greedy trajectory and retains its reportable prefixes. + * @param interaction initial seed, modified to its final state + * @param hybrid initial seed hybridization energy including initiation + * @param last1 last permitted local index in sequence 1 + * @param last2 last permitted local index in reversed sequence 2 + */ + void extendSeed(Interaction & interaction, E_type hybrid, + size_t last1, size_t last2); + + /** + * Checks one geometrically bounded move and evaluates its complete energy. + * @param candidate move indices/side/gaps; energies are filled on success + * @param bounds current interaction boundaries + * @param hybrid current hybridization energy, including initiation + * @param total current full interaction energy + * @return whether the entire move is feasible and strictly downhill + */ + bool evaluate(Candidate & candidate, const Boundary & bounds, + E_type hybrid, E_type total) const; + + /** @return whether a move wins by exact score and deterministic ties */ + bool isBetter(const Candidate & candidate, const Candidate & best) const; + + /** + * Reduces a valid prefix by boundaries, total energy and full-path ties. + * @param interaction actual path to consider for reporting + */ + void retain(const Interaction & interaction); +}; + +} // namespace IntaRNA + +#endif /* INTARNA_PREDICTORSEEDEXTENSIONKINETIC_H_ */ diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 583b61e3..d39c2b0c 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -50,6 +50,7 @@ extern "C" { #include "IntaRNA/PredictorMfe2dSeed.h" #include "IntaRNA/PredictorMfe2dSeedExtension.h" #include "IntaRNA/PredictorMfe2dSeedExtensionRIblast.h" +#include "IntaRNA/PredictorSeedExtensionKinetic.h" #include "IntaRNA/PredictorMfe2dHeuristicSeedExtension.h" #include "IntaRNA/PredictorMfeEnsSeedOnly.h" @@ -183,7 +184,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) temperature("temperature",0,100,37), model("model", "SPBX", 'X'), - mode("mode", "HMSR", 'H'), // R for RIblast heuristic only + mode("mode", "HMSRK", 'H'), // R for RIblast heuristic only + kineticScore("kineticScore", "ABC", 'A'), #if INTARNA_MULITHREADING threads("threads", 0, omp_get_max_threads(), 1), #endif @@ -803,8 +805,16 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , std::string("prediction mode : " "\n 'H' = heuristic (fast and low memory), " "\n 'M' = exact (slow), " - "\n 'S' = seed-only" + "\n 'S' = seed-only, " + "\n 'K' = downhill greedy seed extension (requires --model=X; no time or rate prediction)" ).c_str()) + (kineticScore.name.c_str() + , value(&(kineticScore.val)) + ->default_value(kineticScore.def) + ->notifier(boost::bind(&CommandLineParsing::validate_charArgument,this,kineticScore,_1)) + , "candidate score for --model=X --mode=K: 'A' = complete interaction energy change, " + "'B' = change/(1+s1+s2), 'C' = change/(1+2*max(s1,s2)), where s1/s2 are skipped bases. " + "All modes accept strictly negative energy changes only; B/C are heuristic scores.") (model.name.c_str() , value(&(model.val)) ->default_value(model.def) @@ -1222,6 +1232,24 @@ parse(int argc, char** argv) // parsing escape literals outSep = unescaped_string::getUnescaped( outSep ); + // K needs a seed even before the usual --noSeed model normalization. + if (mode.val == 'K') { + if (model.val != 'X') throw error("--mode=K is available only with --model=X"); + if (noSeedRequired) throw error("--mode=K requires seeds and is incompatible with --noSeed"); + // A selected set of greedy paths is not an equilibrium ensemble. + if (outMode.val == 'E' + || (outMode.val == 'C' && OutputHandlerCsv::needsZall(OutputHandlerCsv::string2list(outCsvCols))) + || !outPrefix2streamName.at(OutPrefixCode::OP_spotProb).empty() + || !outPrefix2streamName.at(OutPrefixCode::OP_spotProbAll).empty() + || !outPrefix2streamName.at(OutPrefixCode::OP_qSpotProb).empty() + || !outPrefix2streamName.at(OutPrefixCode::OP_tSpotProb).empty()) + { + throw error("--mode=K does not support equilibrium ensemble or interaction-probability output"); + } + } else if (vm.count(kineticScore.name) && !vm.at(kineticScore.name).defaulted()) { + throw error("--kineticScore requires --model=X --mode=K"); + } + // open output stream { // open according stream @@ -1669,7 +1697,7 @@ CommandLineParsing::prepareEvaluation( boost::program_options::variables_map & v // Erase before notify(), so even out-of-sequence regions and seed encodings // cannot constrain or invalidate evaluation. Energy/accessibility options stay. const std::set ignored = { - "model", "mode", "noSeed", "intLenMax", "qIntLenMax", "tIntLenMax", + "model", "mode", "kineticScore", "noSeed", "intLenMax", "qIntLenMax", "tIntLenMax", "intLoopMax", "qIntLoopMax", "tIntLoopMax", "qRegion", "tRegion", "qRegionLenMax", "tRegionLenMax", "windowWidth", "windowOverlap", "outNumber", "outOverlap", "outMaxE", "outDeltaE", "outMinPu", @@ -2554,6 +2582,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const case 'H' : return new PredictorMfe2dHeuristicSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); case 'M' : return new PredictorMfe2dSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); case 'R' : return new PredictorMfe2dSeedExtensionRIblast( energy, output, predTracker, getSeedHandler( energy ) ); + case 'K' : return new PredictorSeedExtensionKinetic( energy, output, predTracker, getSeedHandler( energy ), kineticScore.val ); case 'S' : return new PredictorMfeSeedOnly( energy, output, predTracker, getSeedHandler( energy ) ); default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented"); return NULL; } diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index 3731406b..115645b1 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -721,6 +721,8 @@ class CommandLineParsing { CharParameter model; //! the prediction mode (heuristic, space-efficient, exact) CharParameter mode; + //! greedy seed-extension score: A=energy, B=Manhattan, C=asymmetry + CharParameter kineticScore; #if INTARNA_MULITHREADING //! number of threads = number of parallel predictors running NumberParameter threads; diff --git a/tests/Makefile.am b/tests/Makefile.am index 571deef5..565a5b57 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -15,7 +15,7 @@ TEST_EXTENSIONS = $(EXEEXT) .sh SH_LOG_COMPILER = $(SHELL) # the script needed for tests -dist_check_SCRIPTS = runIntaRNA.sh runAccessibilityBinary.sh runIntaRNAeval.sh +dist_check_SCRIPTS = runIntaRNA.sh runAccessibilityBinary.sh runIntaRNAeval.sh runKineticSeedExtension.sh # the program to build check_PROGRAMS = runApiTests @@ -52,6 +52,7 @@ runApiTests_SOURCES = \ PredictorTinyOracle_test.cpp \ PredictorEvalOnly_test.cpp \ PredictorSeedOracle_test.cpp \ + PredictorSeedExtensionKinetic_test.cpp \ Matrix_test.cpp \ NussinovHandler_test.cpp \ RnaSequence_test.cpp \ diff --git a/tests/PredictorSeedExtensionKinetic_test.cpp b/tests/PredictorSeedExtensionKinetic_test.cpp new file mode 100644 index 00000000..cbd47d4d --- /dev/null +++ b/tests/PredictorSeedExtensionKinetic_test.cpp @@ -0,0 +1,458 @@ +#include "catch.hpp" + +#undef NDEBUG + +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/InteractionEnergyVrna.h" +#include "IntaRNA/OutputHandler.h" +#include "IntaRNA/PredictorSeedExtensionKinetic.h" +#include "IntaRNA/SeedHandlerExplicit.h" +#include "IntaRNA/VrnaHandler.h" + +#include +#include +#include +#include +#include +#include +#include + +using namespace IntaRNA; + +namespace { + +using Pair = std::pair; +using Chain = std::vector; +using Bounds = std::array; + +class KineticOutput final : public OutputHandler { +public: + explicit KineticOutput(const OutputConstraint & c); + void add(const Interaction & i) override; + std::vector interactions; +}; + +KineticOutput::KineticOutput(const OutputConstraint & c) : OutputHandler(c) {} + +void KineticOutput::add(const Interaction & i) { + if (!i.basePairs.empty()) interactions.push_back(i); + ++reportedInteractions; +} + +// A nonmonotone table deliberately also models imported accessibility data. +class KineticAccessibility final : public AccessibilityDisabled { +public: + KineticAccessibility(const RnaSequence & s, size_t maxLength = 0); + E_type getED(size_t from, size_t to) const override; + std::map values; +}; + +KineticAccessibility::KineticAccessibility(const RnaSequence & s, size_t maxLength) + : AccessibilityDisabled(s, maxLength, NULL) {} + +E_type KineticAccessibility::getED(size_t from, size_t to) const { + const E_type base = AccessibilityDisabled::getED(from, to); + if (base == ED_UPPER_BOUND) return base; + auto i = values.find({from, to}); + return i == values.end() ? 0 : i->second; +} + +class KineticEnergy final : public InteractionEnergyBasePair { +public: + KineticEnergy(const Accessibility & a, const ReverseAccessibility & b, + size_t m1 = 3, size_t m2 = 3); + E_type getE_interLeft(size_t i, size_t j, size_t k, size_t l) const override; + E_type getE(size_t i, size_t j, size_t k, size_t l, E_type h) const override; + bool customLoops = false; + std::map loops; + std::map boundaryTerms; +}; + +KineticEnergy::KineticEnergy(const Accessibility & a, const ReverseAccessibility & b, + size_t m1, size_t m2) + : InteractionEnergyBasePair(a, b, m1, m2, false, 1., -100, 3, 0, false) {} + +E_type KineticEnergy::getE_interLeft(size_t i, size_t j, size_t k, size_t l) const { + if (!isValidInternalLoop(i, j, k, l)) return E_INF; + if (!customLoops) return InteractionEnergyBasePair::getE_interLeft(i, j, k, l); + auto p = loops.find({i,j,k,l}); + return p == loops.end() ? E_INF : p->second; +} + +E_type KineticEnergy::getE(size_t i, size_t j, size_t k, size_t l, E_type h) const { + const E_type base = InteractionEnergyBasePair::getE(i,j,k,l,h); + auto p = boundaryTerms.find({i,j,k,l}); + return E_isINF(base) ? E_INF : base + (p == boundaryTerms.end() ? 0 : p->second); +} + +struct KineticFixture { + RnaSequence first, second; + KineticAccessibility acc1, acc2; + ReverseAccessibility reversed; + KineticEnergy energy; + KineticFixture(size_t n = 7, size_t m1 = 3, size_t m2 = 3, + size_t span1 = 0, size_t span2 = 0); +}; + +KineticFixture::KineticFixture(size_t n, size_t m1, size_t m2, size_t span1, size_t span2) + : first("target", std::string(n, 'G')), second("query", std::string(n, 'C')), + acc1(first,span1), acc2(second,span2), reversed(acc2), energy(acc1,reversed,m1,m2) {} + +bool stacked(const Pair & a, const Pair & b) { + return a.first+1 == b.first && a.second+1 == b.second; +} + +Bounds bounds(const Chain & chain) { + return {chain.front().first,chain.back().first,chain.front().second,chain.back().second}; +} + +// Independent reference: reconstruct every candidate's complete chain, check +// its topology and recompute all loop and boundary contributions from scratch. +E_type chainEnergy(const InteractionEnergy & energy, const Chain & chain, + const OutputConstraint & out) +{ + const auto b = bounds(chain); + if (b[1]-b[0]+1 > energy.getAccessibility1().getMaxLength() + || b[3]-b[2]+1 > energy.getAccessibility2().getMaxLength()) return E_INF; + E_type h = energy.getE_init(); + for (size_t p = 0; p < chain.size(); ++p) { + if (!energy.areComplementary(chain[p].first,chain[p].second)) return E_INF; + const bool before = p != 0 && stacked(chain[p-1],chain[p]); + const bool after = p+1 < chain.size() && stacked(chain[p],chain[p+1]); + if (out.noLP && !before && !after) return E_INF; + if (p == 0) continue; + if (out.noGUend && !before && (energy.isGU(chain[p-1].first,chain[p-1].second) + || energy.isGU(chain[p].first,chain[p].second))) return E_INF; + const E_type loop = energy.getE_interLeft(chain[p-1].first,chain[p].first, + chain[p-1].second,chain[p].second); + if (E_isINF(loop)) return E_INF; + h += loop; + } + return energy.getE(b[0],b[1],b[2],b[3],h); +} + +std::vector oracle(const InteractionEnergy & energy, Chain chain, + const OutputConstraint & out, char score, + const IndexRange & r1, const IndexRange & r2) +{ + std::vector result; + if (E_isINF(chainEnergy(energy,chain,out))) return result; + while (true) { + const auto b = bounds(chain); + const E_type current = chainEnergy(energy,chain,out); + if ((!out.noGUend || (!energy.isGU(b[0],b[2]) && !energy.isGU(b[1],b[3]))) + && energy.getED1(b[0],b[1]) <= out.maxED + && energy.getED2(b[2],b[3]) <= out.maxED) result.push_back(chain); + bool found = false; + Chain best; + std::tuple bestKey; + // Enumerate absolute endpoints in deliberately different order from + // production's side/loop-size enumeration. + for (size_t x = r1.from; x <= r1.to; ++x) { + for (size_t y = r2.from; y <= r2.to; ++y) { + const bool left = x < b[0] && y < b[2]; + const bool right = x > b[1] && y > b[3]; + if (!left && !right) continue; + const size_t s1 = left ? b[0]-x-1 : x-b[1]-1; + const size_t s2 = left ? b[2]-y-1 : y-b[3]-1; + if (s1 > energy.getMaxInternalLoopSize1() || s2 > energy.getMaxInternalLoopSize2()) continue; + Chain trial = chain; + trial.push_back({x,y}); + if (out.noLP && s1+s2 != 0) { + if (left && (x == r1.from || y == r2.from)) continue; + if (right && (x == r1.to || y == r2.to)) continue; + trial.push_back(left ? Pair{x-1,y-1} : Pair{x+1,y+1}); + } + std::sort(trial.begin(),trial.end()); + const E_type total = chainEnergy(energy,trial,out); + if (E_isINF(total) || total >= current) continue; + const size_t denominator = score == 'A' ? 1 : score == 'B' ? 1+s1+s2 : 1+2*std::max(s1,s2); + const auto key = std::make_tuple(static_cast(total-current)/denominator, + left ? 0 : 1,s1+s2,s1); + if (!found || key < bestKey) { found = true; bestKey = key; best = trial; } + } + } + if (!found) return result; + chain = best; + } +} + +std::string seedEncoding(const Chain & seed, size_t n2) { + const auto b = bounds(seed); + std::string a(b[1]-b[0]+1,'.'), bReverse(b[3]-b[2]+1,'.'); + for (const auto & p : seed) { a[p.first-b[0]] = '|'; bReverse[p.second-b[2]] = '|'; } + std::reverse(bReverse.begin(),bReverse.end()); + return std::to_string(b[0]+1)+a+"&"+std::to_string(n2-b[3])+bReverse; +} + +SeedConstraint seedConstraint(const std::string & explicitSeed) { + return SeedConstraint(2,20,20,20,E_INF,Accessibility::ED_UPPER_BOUND,E_INF, + IndexRangeList(""),IndexRangeList(""),explicitSeed,false,false,false); +} + +std::vector predict(const InteractionEnergy & energy, const Chain & seed, + const OutputConstraint & out, char score = 'A', + const IndexRange & r1 = IndexRange(0,RnaSequence::lastPos), + const IndexRange & r2 = IndexRange(0,RnaSequence::lastPos)) +{ + const auto sc = seedConstraint(seedEncoding(seed,energy.size2())); + KineticOutput output(out); + PredictorSeedExtensionKinetic predictor(energy,output,NULL,new SeedHandlerExplicit(energy,sc),score); + predictor.predict(r1,r2); + return output.interactions; +} + +Chain internalChain(const InteractionEnergy & energy, const Interaction & interaction) { + Chain result; + for (const auto & p : interaction.basePairs) result.push_back({energy.getIndex1(p),energy.getIndex2(p)}); + return result; +} + +void checkOracle(const InteractionEnergy & energy, const Chain & seed, + const OutputConstraint & out, char score, + IndexRange r1 = IndexRange(0,RnaSequence::lastPos), + IndexRange r2 = IndexRange(0,RnaSequence::lastPos)) +{ + r1.to = std::min(r1.to,energy.size1()-1); r2.to = std::min(r2.to,energy.size2()-1); + auto expected = oracle(energy,seed,out,score,r1,r2); + expected.erase(std::remove_if(expected.begin(),expected.end(),[&](const Chain & c) { + return chainEnergy(energy,c,out) >= out.maxE; + }),expected.end()); + std::sort(expected.begin(),expected.end(),[&](const Chain & a,const Chain & b) { + return chainEnergy(energy,a,out) < chainEnergy(energy,b,out); + }); + const auto actual = predict(energy,seed,out,score,r1,r2); + REQUIRE(actual.size() == std::min(expected.size(),out.reportMax)); + for (size_t p = 0; p < actual.size(); ++p) { + REQUIRE(internalChain(energy,actual[p]) == expected[p]); + REQUIRE(actual[p].energy == chainEnergy(energy,expected[p],out)); + } +} + +} // namespace + +TEST_CASE("Kinetic extension agrees with independent whole-chain oracle", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + const Chain seed{{2,2},{3,3}}; + for (char score : {'A','B','C'}) { + for (bool noLP : {false,true}) { + KineticFixture fixture(8,1,3,6,7); + // Nonmonotone values distinguish exact candidate EDs from a + // farthest-end penalty and change full boundary energy deltas. + fixture.acc1.values = {{{0,3},700},{{1,3},30},{{2,5},500},{{2,6},20}}; + fixture.acc2.values = {{{2,5},450},{{1,5},20}}; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP); + checkOracle(fixture.energy,seed,out,score); + checkOracle(fixture.energy,seed,out,score,IndexRange(1,6),IndexRange(0,7)); + } + } +} + +TEST_CASE("Kinetic scoring and strict downhill acceptance", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + const Chain seed{{2,2},{3,3}}; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + SECTION("A B and C select distinct justified moves") { + KineticFixture f(8); + f.energy.customLoops = true; + f.energy.loops = {{{2,3,2,3},-100},{{3,6,3,4},-420},{{3,5,3,5},-390},{{3,4,3,4},-120}}; + const std::map expected{{'A',{6,4}},{'B',{6,4}},{'C',{5,5}}}; + for (char score : {'A','B','C'}) { + checkOracle(f.energy,seed,out,score); + const auto actual = predict(f.energy,seed,out,score); + REQUIRE(internalChain(f.energy,actual.front()).back() == expected.at(score)); + } + // B prefers the immediate stack when its per-distance gain wins. + f.energy.loops[{3,4,3,4}] = -150; + REQUIRE(internalChain(f.energy,predict(f.energy,seed,out,'B').front()).back() == Pair(4,4)); + } + SECTION("zero and positive moves stop even when the local loop is favorable") { + for (E_type delta : {0,1,200}) { + KineticFixture f(5); + f.energy.customLoops = true; + f.energy.loops = {{{2,3,2,3},-100},{{3,4,3,4},-100}}; + f.energy.boundaryTerms[{2,4,2,4}] = 100+delta; + const auto actual = predict(f.energy,seed,out); + REQUIRE(actual.size() == 1); + REQUIRE(internalChain(f.energy,actual.front()) == seed); + } + } + SECTION("ties prefer left then smaller total gap then smaller first gap") { + KineticFixture f(7); + f.energy.customLoops = true; + f.energy.loops = {{{2,3,2,3},-100},{{0,2,1,2},-100},{{1,2,0,2},-100},{{1,2,1,2},-100},{{3,4,3,4},-100}}; + // Limit each span to three: whichever move wins blocks the other side. + KineticFixture shortF(7,3,3,3,3); + shortF.energy.customLoops = true; shortF.energy.loops = f.energy.loops; + REQUIRE(internalChain(shortF.energy,predict(shortF.energy,seed,out).front()).front() == Pair(1,1)); + f.energy.loops.erase({1,2,1,2}); + // Equal total gap: s1=0 (pair 1,0) wins over s1=1 (pair 0,1). + const auto paths = oracle(f.energy,seed,out,'A',IndexRange(0,6),IndexRange(0,6)); + REQUIRE(paths.at(1).front() == Pair(1,0)); + checkOracle(f.energy,seed,out,'A'); + } +} + +TEST_CASE("Kinetic noLP macro-steps and explicit seed validation", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(6,1,1); + f.energy.customLoops = true; + f.energy.loops = {{{0,1,0,1},-100},{{1,3,1,3},50},{{3,4,3,4},-330}}; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,true); + const Chain seed{{0,0},{1,1}}; + checkOracle(f.energy,seed,out,'A'); + const auto actual = predict(f.energy,seed,out); + REQUIRE(actual.size() == 2); + REQUIRE(actual.front().basePairs.size() == 4); + REQUIRE(actual.front().energy == -480); + SECTION("B and C count one macro move in their denominator") { + KineticFixture ranked(7,1,1,5,5); + ranked.energy.customLoops = true; + ranked.energy.loops = {{{2,3,2,3},-100},{{1,2,1,2},-80}, + {{3,5,3,5},50},{{5,6,5,6},-330}}; + const Chain middle{{2,2},{3,3}}; + for (char score : {'B','C'}) { + checkOracle(ranked.energy,middle,out,score); + const auto rankedResult = predict(ranked.energy,middle,out,score); + // -280/3 beats -80. Incorrectly counting the two added pairs + // would instead produce -280/4 and choose the left stack. + REQUIRE(rankedResult.front().basePairs.size() == 4); + REQUIRE(internalChain(ranked.energy,rankedResult.front()).back() == Pair(6,6)); + } + } + SECTION("lonely explicit seeds are skipped") { + REQUIRE(predict(f.energy,Chain{{1,1},{3,3}},out).empty()); + } + SECTION("both additional pairs must fit both strand spans") { + KineticFixture shortF(6,1,1,4,6); + shortF.energy.customLoops = true; shortF.energy.loops = f.energy.loops; + const auto limited = predict(shortF.energy,seed,out); + REQUIRE(limited.size() == 1); + REQUIRE(limited.front().basePairs.size() == 2); + } +} + +TEST_CASE("Kinetic noGU output retains a valid prefix", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + RnaSequence first("target","GGG"), second("query","UCC"); + AccessibilityDisabled a(first,0,NULL), b(second,0,NULL); + ReverseAccessibility reversed(b); + InteractionEnergyBasePair energy(a,reversed,1,1); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,false,true); + const Chain seed{{0,0},{1,1}}; + checkOracle(energy,seed,out,'A'); + const auto actual = predict(energy,seed,out); + REQUIRE(actual.size() == 1); + REQUIRE(internalChain(energy,actual.front()) == seed); + SECTION("a transient GU endpoint can become internal after another stack") { + RnaSequence t("target","GGGG"), q("query","CUCC"); + AccessibilityDisabled at(t,0,NULL), aq(q,0,NULL); + ReverseAccessibility rq(aq); + InteractionEnergyBasePair model(at,rq,1,1); + checkOracle(model,seed,out,'A'); + const auto recovered = predict(model,seed,out); + REQUIRE(recovered.size() == 2); + REQUIRE(recovered.front().basePairs.size() == 4); + REQUIRE(recovered.front().energy == -400); + REQUIRE(model.isGU(2,2)); + } +} + +TEST_CASE("Kinetic full ViennaRNA energy agrees with whole-chain oracle", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + RnaSequence first("target","AGCGACGCA"), second("query","UGCGUCGCU"); + KineticAccessibility a(first), b(second); + a.values = {{{0,4},140},{{1,4},20},{{1,5},350},{{1,6},30},{{2,7},130}}; + b.values = {{{0,4},20},{{1,4},60},{{1,5},20},{{2,6},90}}; + ReverseAccessibility reversed(b); + VrnaHandler vrna(37,"Turner04",false,false); + const Chain seed{{2,2},{3,3}}; + for (bool dangles : {false,true}) { + InteractionEnergyVrna energy(a,reversed,vrna,3,2,false,37,dangles); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + for (char score : {'A','B','C'}) checkOracle(energy,seed,out,score); + } +} + +TEST_CASE("Kinetic refuses undefined ensemble output and invalid scores", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f; + const auto sc = seedConstraint("3||&4||"); + OutputConstraint ensemble(1,OutputConstraint::OVERLAP_BOTH,0,E_INF,false,false,false,true); + KineticOutput ensembleOut(ensemble); + REQUIRE_THROWS_AS(PredictorSeedExtensionKinetic(f.energy,ensembleOut,NULL, + new SeedHandlerExplicit(f.energy,sc)),std::invalid_argument); + OutputConstraint ordinary; + KineticOutput ordinaryOut(ordinary); + REQUIRE_THROWS_AS(PredictorSeedExtensionKinetic(f.energy,ordinaryOut,NULL, + new SeedHandlerExplicit(f.energy,sc),'Z'),std::invalid_argument); +} + +TEST_CASE("Kinetic Turner loop is rescued by its atomic stack", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + // A single query-strand A bulge prevents a direct stack after the seed. + // This is the Turner2004 witness reproduced by the first council: + // +0.50 kcal/mol loop plus -3.30 kcal/mol outward stack. + RnaSequence first("target","CCCC"), second("query","GGAGG"); + AccessibilityDisabled a(first,0,NULL), b(second,0,NULL); + ReverseAccessibility reversed(b); + VrnaHandler vrna(37,"Turner04",false,false); + InteractionEnergyVrna energy(a,reversed,vrna,1,1,false,0,false); + const E_type loop = energy.getE_interLeft(1,2,1,3); + const E_type stack = energy.getE_interLeft(2,3,3,4); + REQUIRE(loop == 50); + REQUIRE(stack == -330); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,true); + const Chain seed{{0,0},{1,1}}; + checkOracle(energy,seed,out,'A'); + const auto actual = predict(energy,seed,out); + REQUIRE(actual.size() == 2); + REQUIRE((internalChain(energy,actual.front()) == Chain{{0,0},{1,1},{2,3},{3,4}})); + REQUIRE(actual.front().energy - actual.back().energy == loop+stack); +} + +TEST_CASE("Kinetic nonoverlap reporting recovers shorter prefixes", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(8,0,0,5,5); + f.energy.customLoops = true; + for (size_t p = 0; p+1 < 8; ++p) f.energy.loops[{p,p+1,p,p+1}] = p < 3 ? -100 : -200; + const auto sc = seedConstraint("1||&7||,4||&4||"); + for (bool needBPs : {false,true}) { + OutputConstraint out(3,OutputConstraint::OVERLAP_NONE,E_INF,E_INF,false,false,false,false,needBPs); + KineticOutput output(out); + PredictorSeedExtensionKinetic predictor(f.energy,output,NULL,new SeedHandlerExplicit(f.energy,sc)); + predictor.predict(); + REQUIRE(output.interactions.size() == 2); + REQUIRE(output.interactions[0].energy == -900); + REQUIRE(output.interactions[1].energy == -300); + REQUIRE((bounds(internalChain(f.energy,output.interactions[0])) == Bounds{{3,7,3,7}})); + REQUIRE((bounds(internalChain(f.energy,output.interactions[1])) == Bounds{{0,2,0,2}})); + if (needBPs) { + REQUIRE(output.interactions[0].basePairs.size() == 5); + REQUIRE(output.interactions[1].basePairs.size() == 3); + } + // Reuse the instance with a different range and ensure old candidates + // and seed metadata cannot leak into the second prediction. + output.interactions.clear(); + predictor.predict(IndexRange(0,2),IndexRange(0,2)); + REQUIRE(output.interactions.size() == 1); + REQUIRE(output.interactions.front().energy == -300); + REQUIRE((bounds(internalChain(f.energy,output.interactions.front())) == Bounds{{0,2,0,2}})); + } +} + +TEST_CASE("Kinetic annotations exclude rejected lonely explicit seeds", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(5,0,0); + const auto sc = seedConstraint("2||&3||,1|&5|"); + OutputConstraint out(1,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,true); + KineticOutput output(out); + PredictorSeedExtensionKinetic predictor(f.energy,output,NULL,new SeedHandlerExplicit(f.energy,sc)); + predictor.predict(); + REQUIRE(output.interactions.size() == 1); + const auto & result = output.interactions.front(); + REQUIRE(result.basePairs.size() == 5); + REQUIRE(result.seed != NULL); + REQUIRE(result.seed->size() == 1); +} diff --git a/tests/runKineticSeedExtension.sh b/tests/runKineticSeedExtension.sh new file mode 100755 index 00000000..1d98dc67 --- /dev/null +++ b/tests/runKineticSeedExtension.sh @@ -0,0 +1,119 @@ +#!/usr/bin/env bash +# Exercise kinetic dispatch, greedy paths, output constraints and energy accounting. +set -euo pipefail +bin="$INTARNABINPATH/src/bin/IntaRNA" +tmp=$(mktemp -d) +trap 'rm -rf "$tmp"' EXIT +common=(--target=GGGGGG --query=CCCCCC --energy=B --acc=N + --threads=1 --default-log-file=/dev/null) +kinetic=(--model=X --mode=K '--seedTQ=3||&3||') +csv=(--outMode=C --outCsvCols=hybridDB,E) + +# Six base pairs have independently known energy -6 in the base-pair model. +printf 'hybridDB;E\n1||||||&1||||||;-6\n' > "$tmp/expected" +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" > "$tmp/default" +cmp "$tmp/expected" "$tmp/default" +for score in A B C; do + for noLP in false true; do + "$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" \ + --kineticScore="$score" --outNoLP="$noLP" --outNoGUend > "$tmp/score" + cmp "$tmp/expected" "$tmp/score" + done +done + +# Ties choose the left extension, including the antiparallel query coordinates. +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" --intLenMax=3 > "$tmp/left" +printf 'hybridDB;E\n2|||&3|||;-3\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/left" + +# Boundary-only output and CLI range/index conversions use the same trajectory. +"$bin" "${common[@]}" "${kinetic[@]}" --outMode=C \ + --outCsvCols=start1,end1,start2,end2,E > "$tmp/boundaries" +printf 'start1;end1;start2;end2;E\n1;6;1;6;-6\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/boundaries" +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" \ + --tRegion=2-5 --qRegion=2-5 > "$tmp/ranges" +printf 'hybridDB;E\n2||||&2||||;-4\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/ranges" +# Signed IntaRNA coordinates skip zero: target position three is labeled 1. +"$bin" "${common[@]}" --model=X --mode=K '--seedTQ=1||&7||' "${csv[@]}" \ + --tIdxPos0=-2 --qIdxPos0=5 > "$tmp/shifted" +printf 'hybridDB;E\n-2||||||&5||||||;-6\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/shifted" + +# Every committed prefix is available for suboptimal output, in energy order. +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" --outNumber=10 > "$tmp/prefixes" +printf 'hybridDB;E\n1||||||&1||||||;-6\n1|||||&2|||||;-5\n1||||&3||||;-4\n2|||&3|||;-3\n3||&3||;-2\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/prefixes" + +# A zero report count still permits energy tracking and never needs traceback. +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" --outNumber=0 \ + --out="tMinE:$tmp/min-energy" > "$tmp/zero" +printf 'hybridDB;E\n' > "$tmp/empty" +cmp "$tmp/empty" "$tmp/zero" +test -s "$tmp/min-energy" +grep -q -- '-6' "$tmp/min-energy" + +# Output constraints apply to explicit seeds too. +"$bin" "${common[@]}" --model=X --mode=K '--seedTQ=3|&3|' \ + "${csv[@]}" --outNoLP > "$tmp/lonely" +cmp "$tmp/empty" "$tmp/lonely" +"$bin" --target=GG --query=UU --energy=B --acc=N --model=X --mode=K \ + '--seedTQ=1||&1||' "${csv[@]}" --outNoGUend --default-log-file=/dev/null > "$tmp/gu" +cmp "$tmp/empty" "$tmp/gu" + +# Config-file selection must behave just like command-line selection. +printf 'model=X\nmode=K\nkineticScore=C\nseedTQ=3||&3||\n' > "$tmp/parameters" +"$bin" "${common[@]}" "${csv[@]}" --parameterFile="$tmp/parameters" > "$tmp/configured" +cmp "$tmp/default" "$tmp/configured" + +# Re-evaluate all reported structures independently through IntaRNAeval. +# This checks the complete energy, including accessibility, terminal and dangle terms. +thermo=(--target=AGCGACGCA --query=UGCGUCGCU --accW=0 --accL=0 --temperature=25 + --energyAdd=1.2 --outMode=C --outMaxE=100 --outNumber=5 + --outCsvCols=hybridDB,E,ED1,ED2,E_init,E_loops,E_dangleL,E_dangleR,E_endL,E_endR,E_add + --threads=1 --default-log-file=/dev/null) +for score in A B C; do + for dangles in false true; do + "$bin" "${thermo[@]}" --model=X --mode=K '--seedTQ=3|||&5|||' \ + --kineticScore="$score" --energyNoDangles="$dangles" --outNoLP --outNoGUend \ + > "$tmp/predicted" + test "$(wc -l < "$tmp/predicted")" -gt 1 + structures=$(awk -F';' 'NR>1 {printf "%s%s", sep, $1; sep=":"}' "$tmp/predicted") + "$bin" "${thermo[@]}" --rri="$structures" --energyNoDangles="$dangles" > "$tmp/evaluated" + cmp "$tmp/predicted" "$tmp/evaluated" + done +done + +expect_error() { + local expected="$1" + shift + local status=0 + "$bin" "$@" > "$tmp/bad.out" 2> "$tmp/bad.err" || status=$? + test "$status" -eq 1 || test "$status" -eq 255 + grep -q -- "$expected" "$tmp/bad.out" "$tmp/bad.err" +} +for model in S P B; do + expect_error 'only with --model=X' "${common[@]}" --model="$model" --mode=K +done +expect_error 'incompatible with --noSeed' "${common[@]}" "${kinetic[@]}" --noSeed +# Explicitly supplying the default A must be rejected outside K as well. +for score in A B C; do + expect_error 'kineticScore requires' "${common[@]}" --mode=H --kineticScore="$score" +done +expect_error 'kineticScore' "${common[@]}" "${kinetic[@]}" --kineticScore=D +expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" --outMode=E +for column in Zall Eall EallTotal P_E; do + expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" \ + --outMode=C --outCsvCols="E,$column" +done +for output in spotProb qSpotProb tSpotProb; do + expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" --out="$output:$tmp/probability" +done +expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" --out="spotProb:1&1:$tmp/probability" + +# Evaluation ignores prediction controls, including invalid kinetic scores. +"$bin" "${common[@]}" "${csv[@]}" '--rri=1||||||&1||||||' \ + --model=S --mode=K --kineticScore=D > "$tmp/evaluation" +cmp "$tmp/default" "$tmp/evaluation" +echo 'Kinetic seed-extension CLI checks passed' From 0f0736aeddb831050f43b3e861b13290fedc5913 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Mon, 5 Oct 2026 14:34:49 +0200 Subject: [PATCH 2/6] Revise kinetic extensions and add experimental pruning mode --- ChangeLog | 22 ++ README.md | 30 +- doc/Makefile.am | 2 +- doc/benchmark-kinetic.py | 77 ++++ doc/kinetic-seed-extension.md | 347 ++++++++---------- src/IntaRNA/InteractionEnergyVrna.h | 16 + src/IntaRNA/Makefile.am | 2 + src/IntaRNA/PredictorSeedExtensionKinetic.cpp | 281 +++++++------- src/IntaRNA/PredictorSeedExtensionKinetic.h | 69 ++-- .../PredictorSeedExtensionKineticPruned.cpp | 109 ++++++ .../PredictorSeedExtensionKineticPruned.h | 65 ++++ src/bin/CommandLineParsing.cpp | 25 +- tests/PredictorSeedExtensionKinetic_test.cpp | 248 +++++++++++-- tests/runKineticSeedExtension.sh | 42 ++- 14 files changed, 898 insertions(+), 437 deletions(-) create mode 100644 doc/benchmark-kinetic.py create mode 100644 src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp create mode 100644 src/IntaRNA/PredictorSeedExtensionKineticPruned.h diff --git a/ChangeLog b/ChangeLog index 464319e5..f89ff44c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -17,6 +17,8 @@ - deterministic kinetic seed extension with --model=X --mode=K and configurable local scoring; complete-energy downhill steps, atomic noLP loop/stack moves and retained greedy traceback; equilibrium probability output is unsupported +- kinetic modes always use stacked extensions and trust handler-provided seeds; + reuse candidate pairing/local energies and add experimental pruning mode L - IntaRNAeval / --rri evaluates predefined RNA-RNA interactions (issue #184) @@ -77,6 +79,26 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261005 Alexander Mitrofanov + * IntaRNA/PredictorSeedExtensionKinetic : + * always evaluate single and double stacks plus atomic loop/stack extensions; + trust seed-handler structures/energies, including explicit lonely pairs + * cache each end's candidates, shared complementarity checks and local + energies; rebuild only the chosen end and refresh full opposite-end energy + * IntaRNA/PredictorSeedExtensionKineticPruned, InteractionEnergyVrna : + + experimental root-pair local-energy suffix bounds from active parameters; + prune before complementarity checks assuming monotone ED, omitting endpoint + changes; preserve exact full-energy evaluation of surviving moves + * bin/CommandLineParsing, src/IntaRNA/Makefile.am : + + expose subclass as --mode=L; set --outNoLP=true with INFO for K/L + * tests/PredictorSeedExtensionKinetic_test.cpp, tests/runKineticSeedExtension.sh : + * update independent oracle and CLI expectations for atomic double stacks; + cover trusted seeds, caching, pruning limitations and both CLI modes + * README.md, doc/kinetic-seed-extension.md, doc/benchmark-kinetic.py, + doc/Makefile.am : + * document revised semantics and reproducible performance evaluation in + response to https://github.com/BackofenLab/IntaRNA/pull/254 + 261003 Alexander Mitrofanov * IntaRNA/PredictorSeedExtensionKinetic, src/IntaRNA/Makefile.am : + greedily extend seeds on either side using complete interaction-energy diff --git a/README.md b/README.md index cc0c72b3..09a1b483 100644 --- a/README.md +++ b/README.md @@ -734,10 +734,11 @@ and studied using the `S` mode. `--model=X --mode=K` grows each available seed along a deterministic greedy path. Every step compares feasible extensions on both sides using the complete change in interaction energy, including accessibility, terminal penalties and -dangling ends. Only strictly negative changes are accepted. With `--outNoLP`, -crossing a loop forms its closing pair and the immediately following stack as -one atomic step; a favorable stack can therefore compensate for an unfavorable -loop. Initial seeds must satisfy the selected structural constraints. +dangling ends. Only strictly negative changes are accepted. Extensions always +add one stacked pair or two stacked pairs, including across a loop; two-pair +moves are evaluated and committed atomically. The CLI sets `--outNoLP=true` +when absent or false and logs an INFO message. Seeds and their energies are +accepted from the seed handler, including explicit seeds with lonely pairs. `--kineticScore` selects the local move ranking: @@ -748,19 +749,24 @@ loop. Initial seeds must satisfy the selected structural constraints. | `C` (C1 in the design) | Complete energy change / `(1+2*max(s1,s2))` | Equal scores prefer the left side, then fewer unpaired bases, then smaller -`s1`. The denominators also apply to two-pair moves. All reportable visited +`s1`, then the single-pair move. The denominators also apply to two-pair moves. All reportable visited states, including seeds, participate in the normal energy-ranked output; traceback preserves the actual chosen path. `--outNoGUend`, separate query and target loop/span limits, regions, output energy/accessibility filters and overlap settings remain applicable. -This mode is a zippering-inspired heuristic, without a calibrated time axis or -a guarantee of the global minimum. It evaluates all feasible local moves; -loop-only energetic pruning is unsafe for complete loop-plus-stack steps. -Equilibrium probability/partition-sum outputs are rejected, as are other models -and `--noSeed`. Scores B and C are optional distance preferences, not measured -kinetic rates. See the [design and implementation plan](doc/kinetic-seed-extension.md) -for the precise algorithm, scientific limitations and validation cases. +`--mode=L` provides an experimental subclass with root-pair-specific local +loop-plus-stack tables and early pruning before complementarity checks. It +uses active energy parameters, assumes monotone accessibility costs, and +ignores terminal/dangling changes in its pruning estimate. It can choose +different paths from K and is available for comparative benchmarking. + +Both modes are zippering-inspired heuristics, without a calibrated time axis +or a guarantee of the global minimum. Equilibrium probability/partition-sum +outputs are rejected, as are other models and `--noSeed`. Scores B and C are +optional distance preferences, not measured kinetic rates. See the +[algorithm and benchmark documentation](doc/kinetic-seed-extension.md) for +precise semantics, pruning limitations and validation cases. [![up](doc/figures/icon-up.28.png) back to overview](#overview) diff --git a/doc/Makefile.am b/doc/Makefile.am index 003a02d3..c3bb4670 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -5,7 +5,7 @@ EXTRA_DIST = \ conda.txt \ - kinetic-seed-extension.md \ + kinetic-seed-extension.md benchmark-kinetic.py \ doxygen.cfg \ latex-deps/adjcalc.sty \ latex-deps/adjustbox.sty \ diff --git a/doc/benchmark-kinetic.py b/doc/benchmark-kinetic.py new file mode 100644 index 00000000..d9efda91 --- /dev/null +++ b/doc/benchmark-kinetic.py @@ -0,0 +1,77 @@ +#!/usr/bin/env python3 +"""Reproducible end-to-end kinetic benchmark; emits timings and output hashes. + +Use an otherwise identical binary with both ends rebuilt after every move as +--uncached to isolate candidate reuse. All runs are single-threaded. Timings +include startup, accessibility, seed generation, table setup and prediction. +""" +import argparse +import hashlib +import json +from pathlib import Path +import random +import statistics +import subprocess +import sys +import tempfile +import time + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument('--cached', required=True, type=Path) + parser.add_argument('--uncached', type=Path) + parser.add_argument('--repeat', type=int, default=5) + args = parser.parse_args() + if args.repeat < 1: + parser.error('--repeat must be positive') + rng = random.Random(254) + target = ''.join(rng.choices('ACGU', k=600)) + query = ''.join(rng.choices('ACGU', k=80)) + cases = [('random-no-ED', target, query, ['--acc=N']), + ('random-folded', target, query, ['--acc=C', '--accW=150', '--accL=100']), + ('stack-rich', 'G'*100, 'C'*30, ['--acc=N'])] + variants = [('cached-K', str(args.cached.resolve()), 'K'), + ('pruned-L', str(args.cached.resolve()), 'L')] + if args.uncached: + variants.insert(0, ('uncached-K', str(args.uncached.resolve()), 'K')) + results = [] + with tempfile.TemporaryDirectory(prefix='intarna-kinetic-bench-') as temp: + for case, t, q, extra in cases: + common = ['--target='+t, '--query='+q, '--model=X', '--seedBP=7', + '--intLenMax=60', '--intLoopMax=10', '--threads=1', + '--outNoLP', '--outNumber=10', '--outMode=C', + '--outCsvCols=hybridDB,E', '--default-log-file=/dev/null'] + extra + measurements = {label: [] for label, _, _ in variants} + outputs = {} + # Warm up every variant, then rotate run order across repetitions. + for iteration in range(args.repeat+1): + order = variants[iteration % len(variants):] + variants[:iteration % len(variants)] + for label, binary, mode in order: + stats = Path(temp)/'time.txt' + command = [binary, '--mode='+mode] + common + start = time.perf_counter() + run = subprocess.run(['/usr/bin/time', '-f', '%M', '-o', str(stats)] + command, + check=True, capture_output=True) + elapsed = time.perf_counter()-start + print(f"{case}/{label} run {iteration}: {elapsed:.3f}s", file=sys.stderr, flush=True) + if label in outputs and outputs[label] != run.stdout: + raise RuntimeError(f'Non-deterministic output: {case}/{label}') + outputs[label] = run.stdout + if iteration: + measurements[label].append((elapsed, int(stats.read_text().strip()))) + if 'uncached-K' in outputs and outputs['uncached-K'] != outputs['cached-K']: + raise RuntimeError(f'Candidate caching changed the result: {case}') + for label, _, _ in variants: + samples = measurements[label] + results.append(dict(case=case, variant=label, + seconds_median=statistics.median(x[0] for x in samples), + rss_KiB_max=max(x[1] for x in samples), + samples_seconds=[x[0] for x in samples], + sha256=hashlib.sha256(outputs[label]).hexdigest(), + equals_K=outputs[label] == outputs['cached-K'])) + print(json.dumps(results, indent=2)) + + +if __name__ == '__main__': + main() diff --git a/doc/kinetic-seed-extension.md b/doc/kinetic-seed-extension.md index c02acf79..fe574978 100644 --- a/doc/kinetic-seed-extension.md +++ b/doc/kinetic-seed-extension.md @@ -1,200 +1,159 @@ # Deterministic kinetic seed extension -## Council decision and scientific scope - -`PredictorSeedExtensionKinetic`, selected exclusively with `--model=X --mode=K`, -implements a deterministic, seed-conditioned, downhill extension heuristic. -It starts from each seed provided by the selected seed handler, compares moves -at both ends of the current duplex, and commits one move at a time. It does not -sample transition rates, simulate elapsed time, cross uphill barriers, remove -base pairs, or guarantee the global minimum free energy interaction. A -macro-step is a coarse move of this heuristic: favorable combined energy does -not establish a barrier-free physical reaction pathway. - -This scope follows the distinction between gradient walks to local minima and -kinetic dynamics in the ViennaRNA ecosystem. [RNAlocmin](https://www.tbi.univie.ac.at/RNA/BHG/RNAlocmin.html) -uses gradient walks; [Kinfold and treekin](https://www.tbi.univie.ac.at/software/) -use stochastic moves or transition-rate matrices. The name "kinetic" identifies -the proposed extension strategy, not a calibrated kinetic prediction. The -choice of a single seed structure per seed start inherits the selected seed -handler's behavior; this algorithm does not enumerate every seed conformation. - -The second council resolved the first review's open questions as follows. -These decisions supersede contradictory optimization and biological claims in -the supplied `PredSeedExtKinetic.md` proposal. - -## State, energy and accepted moves - -A state stores the complete, ordered chain of intermolecular base pairs, its -inclusive boundaries `(i1,j1,i2,j2)`, its hybridization energy `H`, and its full -interaction energy `E`. Indices use IntaRNA's internal coordinate system: -sequence 2 is reversed and prediction-range offsets are handled by the existing -wrappers. The initial `H` includes the seed's loop energies and `getE_init()` -exactly once. Every accepted step retains its actual pairs for traceback. - -Evaluate every candidate with the active `InteractionEnergy` instance: +## Scope and modes + +`--model=X --mode=K` follows a deterministic downhill path from every seed +provided by the selected seed handler. It compares both ends and commits the +best strictly favorable complete move. This is a zippering-inspired heuristic: +it has no calibrated transition rates or time axis, does not cross barriers +between committed states, and does not guarantee a global minimum. A favorable +two-pair move does not establish a barrier-free physical reaction pathway. + +`--mode=L` selects the experimental subclass +`PredictorSeedExtensionKineticPruned`. It applies local-energy and accessibility +pruning **before** checking complementarity and evaluating full energies. +Its additional assumptions can change the selected path; K remains the +reference for measuring those changes. Both modes accept scores A/B/C and +support ordinary energy trackers. They reject seedless operation, other models, +and requests for equilibrium partition functions or probabilities. + +These semantics incorporate the October 5 review of +[PR #254](https://github.com/BackofenLab/IntaRNA/pull/254#issuecomment-5991039020). + +## Seeds, states and allowed extensions + +The predictor trusts the seed handler's structure and hybridization energy. +It does not repeat seed complementarity, noLP or GU-loop checks. Explicit seeds +may contain lonely pairs, including at their ends. Finite energy, prediction +ranges and per-strand span limits still apply. Seed annotations come from the +seed handler without an additional predictor-specific whitelist. + +A state stores its complete ordered base-pair chain, inclusive boundaries +`(i1,j1,i2,j2)`, hybridization energy `H`, and full interaction energy `E`. +Sequence 2 uses reversed energy indices; existing wrappers handle prediction +range offsets and conversion to original coordinates. Initially, +`H = seedHandler.getSeedE(i1,i2) + energy.getE_init()`. + +Extensions **always** use the no-lonely-pair strategy, independent of the API +output constraint. The CLI sets `--outNoLP=true` when absent or false and emits +an INFO message using the normal logging destination. This applies to new +extensions, not to revalidation of handler-provided seeds. Allowed moves are: + +- One pair stacked directly onto the current boundary (`|SEED`). +- Two successive stacked pairs (`||SEED`), evaluated atomically. +- A loop-closing pair followed immediately by its outward stack (`//.SEED`). + +For each strand, `sk` denotes the number of skipped bases between the current +boundary and the closing pair. A single-pair move has `s1=s2=0`. Two-pair moves +include all gap pairs from zero through the separate strand loop limits, +subject to available sequence range and maximum interaction span. They advance +each boundary by `sk+2`. Neither an isolated closing pair nor an intermediate +state of a two-pair move is separately committed or reported. + +Under `--outNoGUend`, a newly formed nonstacking loop must have non-GU closing +pairs. Stacking can temporarily expose a GU outer endpoint; a state is only +reported if its outer endpoints satisfy the flag. The energy model's own +internal-loop GU policy also applies. Earlier reportable prefixes remain +available if a trajectory stops at an unreportable endpoint. + +## Complete energies and deterministic scores + +Every candidate surviving geometric and optional pruning checks is evaluated +with the active energy model: ``` -E_current = energy.getE(i1, j1, i2, j2, H_current) -E_next = energy.getE(i1_next, j1_next, i2_next, j2_next, H_next) -delta = E_next - E_current +H_next = H_current + E_loop + E_stack_if_two_pairs +E_next = energy.getE(i1_next, j1_next, i2_next, j2_next, H_next) +delta = E_next - E_current ``` -This includes both accessibility penalties, both accessibility-weighted dangling -ends, terminal penalties, the configured temperature and parameters, and the -configured additive energy term. Changing one boundary can alter the opposite -end's dangling contribution through its accessibility weight. Consequently, -loop-plus-stack-plus-accessibility differences alone are insufficient. Infinite -states are rejected before subtraction. A constant additive term cancels in the -difference while remaining part of reported energy. - -Only moves with **strictly negative full energy difference** are eligible. -Zero and uphill moves are rejected. Stop when no eligible move remains; the -finite, growing span also guarantees termination. Output energy and accessibility -thresholds filter reported states, rather than introducing additional barriers -in the trajectory. Model-inaccessible states and maximum span limits still -make a move infeasible. - -For each side, enumerate all `0 <= s1 <= m1` and `0 <= s2 <= m2`, where `mk` -is the active energy model's maximum unpaired loop size for strand `k`. -`(s1,s2)=(0,0)` is a stack. Other combinations include bulges and internal loops. -The maximum total skipped length is `m1+m2`, equal to `2*m` only when the -per-strand limits are equal. The proposal's fixed 10 is not a separate limit. - -A normal move adds the complementary loop-closing pair and advances each -boundary by `sk+1`. With `--outNoLP`, a move across any nonzero loop additionally -requires the immediately following outward stack. Evaluate and commit these -two pairs **atomically**, advancing each boundary by `sk+2`. Neither the -isolated closing pair nor its energy is a separately accepted or reported state. -Check both pairs, the intervening loop, and both complete strand spans before -acceptance. Each span must fit that strand's accessibility maximum length and -the requested prediction range; no single ambiguous shared `W` is introduced. - -The complete initial seed must have ordered complementary pairs and finite -per-loop energies. Under `--outNoLP`, every seed pair must already have a direct -stack neighbor. Incompatible explicit seeds are skipped, including seeds with -lonely terminal pairs; this mode does not perform a preliminary seed repair. -Together with atomic macro-steps, this preserves the no-lonely-pair invariant. - -Under `--outNoGUend`, both boundaries of every nonstacking extension loop must -be non-GU. Direct stacking may temporarily expose a GU outer endpoint, as in -existing IntaRNA extension recurrences. A state is reportable only if both outer -endpoints satisfy the flag. The active energy model's internal-loop GU policy -also remains authoritative. A valid earlier state remains available if descent -later stops at an unreportable GU endpoint. - -## Scores and deterministic choice - -`--kineticScore=A|B|C` selects the score; `C` denotes the proposal's C1: - -| Option | Score minimized | Interpretation | -| --- | --- | --- | -| `A` (default) | `delta` | Steepest decrease of the actual modeled interaction energy | -| `B` | `delta/(1+s1+s2)` | Heuristic preference per total skipped length | -| `C` | `delta/(1+2*max(s1,s2))` | Heuristic preference penalizing the longer skipped strand | - -A is the default because it requires no uncalibrated length-to-time assumption. -B and C are retained as explicit alternatives for exploring the supplied -proposal; their denominators are not experimentally calibrated rates or times. -The written `+1` denominator is retained even for a two-pair macro-step: it -counts one candidate move, not the number of pairs added. Scores choose the next -move only. Reported interactions remain ranked by their modeled total energy `E`. - -Compare candidates from **both** sides together. Equal scores are resolved by -left before right, then smaller `s1+s2`, then smaller `s1`. Use sufficiently wide -arithmetic for score comparison, avoiding integer division truncation. Seed -iteration and output comparison follow deterministic existing IntaRNA ordering. - -## Why the proposed pruning is removed - -A loop-only lower bound is not a lower bound for a loop-plus-stack macro-step: -a positive loop may be rescued by the negative following stack. The initial -review reproduced a Turner2004 example at 37 degrees Celsius with a `+0.50` -kcal/mol loop and `-3.30` kcal/mol following stack, totaling `-2.80` kcal/mol. -The full delta also contains terminal and dangling changes absent from the -proposal's filters. Accessibility at the farthest candidate endpoint is not a -certified lower bound for nearer candidates; imported accessibility values need -not obey monotonicity assumptions. - -Therefore the initial implementation exhaustively evaluates all feasible moves -for the current state. It uses neither the six-base-pair-type precomputed Turner -tables nor their `mdspan`/suffix-minimum representation. The suffix-minimum -operation itself is valid, but cannot repair an invalid underlying bound. -Populating an exact full-delta table first would add storage and selection work -without avoiding those energy evaluations. Omitting that cosmetic optimization -is deliberate, rather than replacing one unsafe bound with another. - -There are at most `2*(m1+1)*(m2+1)` candidate shapes per state, with constant-size -local energy updates plus full boundary-energy evaluation for each. No runtime -improvement over other predictors is claimed. Any future pruning must supply a -certified lower bound for the **complete** move under the active energy model -and pass differential tests against this exhaustive implementation. - -## Reporting, traceback and compatibility - -The seed and every atomically committed, reportable prefix are candidates for -normal IntaRNA output. All earlier valid prefixes remain available for suboptimal -and overlap-constrained reporting. Keeping only the last state would lose valid -GU-end prefixes; keeping only the best state for each left boundary would lose -shorter candidates needed after excluding overlap with another report. - -Use existing MFE ordering, output filters, seed annotations and index conversion. -Cache the actual selected pair chain for each retained candidate, including all -macro-step pairs. Traceback must recover that chain rather than run an unrelated -minimum-energy recurrence between its boundaries. Resolve duplicate boundaries -by retaining the lowest total energy, then the lexicographically smallest full -base-pair chain on an energy tie. Reduce duplicates before feeding the ordinary -optimum collector, so no stale energy can select a different cached path. Seed -annotations include only starting seeds that passed validation. Output energies -must agree with an -independent sum over the reported chain and complete boundary terms. - -The cache stores full paths: its memory cost is proportional to the sum of the -retained paths' lengths, in addition to the seed handler's storage. Retaining -prefixes is therefore not a constant-memory walk. It supports reproducible -traceback and shorter alternatives for overlap-constrained output. - -Only `model=X` accepts `mode=K`; seedless operation and `kineticScore` outside K -are rejected. Existing defaults and other models remain unchanged. Ordinary -energy and minimum-energy tracker output are supported. Equilibrium partition -functions and normalized equilibrium probabilities are not defined by this -selected collection of greedy trajectories. Requests needing `Zall`, including -ensemble output and probability trackers, are rejected in CLI validation, with -an API constructor guard for `needZall`. The algorithm does not manufacture an -ensemble by summing repeated prefixes from multiple seeds. - -## Implementation plan and acceptance checklist - -1. Add the public predictor header and implementation, derived from - `PredictorMfe`, with owned seed handler, the existing index-offset wrappers, - validated A/B/C selection and `needZall` rejection. Register both files for - building and installation. -2. Initialize and trace each handler-provided seed; validate its full chain, - seed range, energy, strand spans and active structural constraints. Record - reportable seed states. -3. Enumerate both sides and all feasible loop shapes, form complete normal or - atomic noLP moves, recompute full candidate energy, and select the strictly - downhill winner using the chosen score and the specified tie order. Repeat - until stalling or range exhaustion. -4. Preserve complete committed paths and valid prefixes; integrate ordinary MFE - output and custom candidate lookup for overlap-constrained suboptimals. - Reconstruct seed annotations without changing the retained greedy chain. -5. Wire `--model=X --mode=K` and `--kineticScore=A|B|C` into parsing, help and - factory construction. Reject incompatible model, seedless and ensemble or - probability requests with clear diagnostics. Update README and ChangeLog. -6. Add an independent tiny-sequence reference that enumerates absolute candidate - endpoints and recomputes the entire chain energy. Compare reported energy, - coordinates and traceback against it across scores, constraints and offsets. - Include targeted regressions for a positive-loop/negative-stack rescue, - strict stopping, tie order, nonmonotone accessibility, full boundary-energy - changes, invalid explicit seeds and GU-prefix retention. -7. Run focused API tests, CLI mode/flag compatibility checks, full `make tests`, - debug validation of bounds and ownership, installed standalone-header checks - where supported, and `git diff --check`. Inspect failures; never regenerate - expected outputs merely to hide a change. -8. Review the complete diff and open a pull request documenting this scientific - scope, the deliberate pruning correction, implemented behavior, validation - results and any remaining environmental limitations. - -Test and build results belong in the pull request and development record; this -checklist specifies the required work and does not imply an unrun check passed. +`getE()` includes both accessibility penalties, terminal terms, both weighted +dangling ends, and the configured additive term. Updating one end can change +the opposite end's dangling weight, so complete energies must be refreshed +even when local loop energies are reused. Infinite values are excluded before +subtraction. Only `delta < 0` is accepted. Output energy and accessibility +thresholds filter reports, without imposing additional trajectory barriers. + +| Score | Quantity minimized | +| --- | --- | +| A (default) | `delta` | +| B | `delta / (1+s1+s2)` | +| C | `delta / (1+2*max(s1,s2))` | + +B/C are optional uncalibrated distance preferences. Their denominators count +one move, including two-pair moves. Equal scores prefer left, smaller total +gap, smaller first-strand gap, then a single-pair move. Wide integer cross +products avoid rounding during comparison. Output remains ranked by full `E`. +Two-pair stacking can skip a prefix that the previous implementation visited; +only states actually committed by the revised walk are retained. + +## Candidate storage and reuse + +Each end has a contiguous rectangular table of two-pair candidates plus its +single-stack candidate. A shared position-pair table records unknown, possible +or impossible complementarity. The closing pair of one candidate may be the +outer pair of another; each such check is resolved once per unchanged end. +The first pass resolves pairing and GU feasibility before the energy pass. +The second pass caches local loop-plus-stack energies and updates the best +candidate as it evaluates full energies, without a separate selection pass. + +After committing a move, only that end's tables are rebuilt. The opposite +end's pair checks, feasibility and local energies remain valid. Its total +energy, span eligibility and optional pruning decision are refreshed. An +initially uphill candidate is retained because it may become downhill when +the opposite end changes. Experimentally pruned entries remain unresolved and +can be reconsidered in a later state. Tables use O((m1+2)(m2+2)) space per end. + +## Experimental pruning in L + +For ViennaRNA, precompute local loop-plus-stack minima for each of the six +oriented root base-pair types, both extension sides and each gap pair. Use the +**active temperature-scaled parameter set**, including custom parameter files. +Minimize over the closing and outer pair types and all four adjacent nucleotide +identities. Relaxing consistency between these identities can only make this +local estimate more optimistic. For the base-pair model, two added pairs have +twice its configured base-pair energy. Unknown energy subclasses and API loop +limits above the CLI maximum of 30 fall back to unpruned K enumeration. + +Apply componentwise suffix minima to the gap tables. Each cell then bounds +the local loop-plus-stack term for that gap and every larger gap pair. Before +checking candidate pairs, compare: + +``` +local_suffix_bound + ED1_next + ED2_next - ED1_current - ED2_current >= 0 +``` + +If true, skip that candidate and all componentwise larger gaps. Within the +ordered rectangular traversal this rejects suffixes without more ED, pairing +or loop-energy lookups. A single stacked pair is always evaluated exactly. +Surviving moves still require strictly downhill **complete** energy changes. + +This is deliberately an experiment, not a certified bound on the complete +energy change. It assumes ED is monotone as intervals grow and omits changes in +terminal and dangling contributions. Imported/nonmonotone accessibility and +favorable endpoint changes can invalidate the filter; tests include a concrete +case where L stops while K continues. The local table includes the mandatory +stack, so it avoids the original bare-loop error: a Turner2004 loop of +0.50 +kcal/mol can be rescued by a -3.30 kcal/mol stack. + +Tables require O(12(m1+1)(m2+1)) space and parameter enumeration at construction. +Setup and extra ED lookups may outweigh pruning benefits on small or short-path +inputs. Thus L remains a separate subclass/mode for later real-world evaluation. +See [the reproducible benchmark](benchmark-kinetic.py) and measurements below. + +## Reporting and validation + +Retain each reportable visited prefix and its actual pair chain. For duplicate +boundaries, keep the lowest full energy, then lexicographically smallest chain. +Reduce duplicates before the normal optimum collector. Overlap-constrained +output can select shorter retained prefixes; traceback restores the selected +path directly. Memory for retained paths is proportional to their total length, +not constant per seed. Repeated predictions reset trajectory and ED caches. + +The tests compare K with an independent absolute-endpoint oracle that rebuilds +and reevaluates whole chains. They cover scores/ties, single and double stacks, +positive-loop rescue, strict stopping, separate spans and regions, nonmonotone +ED, GU restrictions, retained prefixes, explicit seeds and annotations, cache +reuse and repeated calls. L has differential and limitation tests. CLI tests +exercise both modes, automatic noLP INFO logging, incompatible requests, and +independent reevaluation of predicted structures through `--rri`. diff --git a/src/IntaRNA/InteractionEnergyVrna.h b/src/IntaRNA/InteractionEnergyVrna.h index 5374f1fb..4dd73c5b 100644 --- a/src/IntaRNA/InteractionEnergyVrna.h +++ b/src/IntaRNA/InteractionEnergyVrna.h @@ -266,6 +266,13 @@ class InteractionEnergyVrna: public InteractionEnergy { Z_type getRT() const; + /** + * Read-only access to the active, temperature-scaled nearest-neighbor + * parameters, e.g. for precomputing local extension-energy estimates. + * @return parameters owned by this energy model, valid for its lifetime + */ + const vrna_param_t & getVrnaParams() const; + protected: @@ -332,6 +339,15 @@ class InteractionEnergyVrna: public InteractionEnergy { //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// +inline +const vrna_param_t & +InteractionEnergyVrna::getVrnaParams() const +{ + return *foldParams; +} + +//////////////////////////////////////////////////////////////////////////// + inline E_type InteractionEnergyVrna:: diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index 5f27f746..f4ec6f61 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -73,6 +73,7 @@ libIntaRNA_a_HEADERS = \ PredictorMfe2dSeedExtension.h \ PredictorMfe2dSeedExtensionRIblast.h \ PredictorSeedExtensionKinetic.h \ + PredictorSeedExtensionKineticPruned.h \ PredictorMfe2dHeuristic.h \ PredictorMfe2dHeuristicSeed.h \ PredictorMfe2dHelixBlockHeuristic.h \ @@ -136,6 +137,7 @@ libIntaRNA_a_SOURCES = \ PredictorMfe2dSeedExtension.cpp \ PredictorMfe2dSeedExtensionRIblast.cpp \ PredictorSeedExtensionKinetic.cpp \ + PredictorSeedExtensionKineticPruned.cpp \ PredictorMfe2dHeuristic.cpp \ PredictorMfe2dHeuristicSeed.cpp \ PredictorMfe2dHelixBlockHeuristic.cpp \ diff --git a/src/IntaRNA/PredictorSeedExtensionKinetic.cpp b/src/IntaRNA/PredictorSeedExtensionKinetic.cpp index 4b085b93..c7990619 100644 --- a/src/IntaRNA/PredictorSeedExtensionKinetic.cpp +++ b/src/IntaRNA/PredictorSeedExtensionKinetic.cpp @@ -41,7 +41,6 @@ PredictorSeedExtensionKinetic::PredictorSeedExtensionKinetic( , seedHandler(checkedSeedHandler(seedHandlerInstance)) , score(score) , interactions() - , validSeeds() { if (score != 'A' && score != 'B' && score != 'C') { throw std::invalid_argument("PredictorSeedExtensionKinetic score must be A, B or C"); @@ -76,7 +75,6 @@ PredictorSeedExtensionKinetic::predict(const IndexRange & r1, const IndexRange & const size_t last1 = std::min(r1.to, size1 - 1) - r1.from; const size_t last2 = std::min(r2.to, size2 - 1) - r2.from; interactions.clear(); - validSeeds.clear(); initOptima(); if (seedHandler.fillSeed(0, last1, 0, last2) != 0) { @@ -104,20 +102,13 @@ PredictorSeedExtensionKinetic::predict(const IndexRange & r1, const IndexRange & interaction.basePairs.push_back(energy.getBasePair(j1, j2)); } interaction.sort(); - E_type hybrid = E_INF; - if (!isValidSeed(interaction, hybrid) - || getBoundary(interaction) != Boundary{i1, j1, i2, j2}) { - continue; - } + const E_type hybrid = addEnergy(seedHandler.getSeedE(i1, i2), energy.getE_init()); interaction.energy = energy.getE(i1, j1, i2, j2, hybrid); if (E_isINF(interaction.energy)) { continue; } interaction.setSeedRange(interaction.basePairs.front(), interaction.basePairs.back(), interaction.energy); - validSeeds.emplace(interaction.basePairs.front(), - Interaction::Seed(interaction.basePairs.front(), - interaction.basePairs.back(), interaction.energy)); extendSeed(interaction, hybrid, last1, last2); } } @@ -149,50 +140,9 @@ PredictorSeedExtensionKinetic::getBoundary(const Interaction & interaction) cons ////////////////////////////////////////////////////////////////////////// bool -PredictorSeedExtensionKinetic::isValidSeed(const Interaction & interaction, E_type & hybrid) const +PredictorSeedExtensionKinetic::prune(const Candidate &, const Boundary &) const { - if (interaction.basePairs.empty() || !interaction.isValid()) { - return false; - } - const auto & pairs = interaction.basePairs; - const auto & constraint = output.getOutputConstraint(); - // Reconstruct rather than trust a cached explicit-seed energy: the - // trajectory's energy must correspond to precisely the traced structure. - hybrid = energy.getE_init(); - if (E_isINF(hybrid)) { - return false; - } - for (size_t p = 0; p < pairs.size(); ++p) { - const size_t i1 = energy.getIndex1(pairs[p]); - const size_t i2 = energy.getIndex2(pairs[p]); - if (i1 >= energy.size1() || i2 >= energy.size2() - || !energy.areComplementary(i1, i2)) { - return false; - } - const bool stackedLeft = p > 0 - && pairs[p].first - pairs[p-1].first == 1 - && pairs[p-1].second - pairs[p].second == 1; - const bool stackedRight = p + 1 < pairs.size() - && pairs[p+1].first - pairs[p].first == 1 - && pairs[p].second - pairs[p+1].second == 1; - if (constraint.noLP && !stackedLeft && !stackedRight) { - return false; - } - if (p == 0) { - continue; - } - const size_t previous1 = energy.getIndex1(pairs[p-1]); - const size_t previous2 = energy.getIndex2(pairs[p-1]); - if (!stackedLeft && constraint.noGUend - && (energy.isGU(previous1, previous2) || energy.isGU(i1, i2))) { - return false; - } - hybrid = addEnergy(hybrid, energy.getE_interLeft(previous1, i1, previous2, i2)); - if (E_isINF(hybrid)) { - return false; - } - } - return true; + return false; } ////////////////////////////////////////////////////////////////////////// @@ -201,66 +151,18 @@ void PredictorSeedExtensionKinetic::extendSeed(Interaction & interaction, E_type hybrid, const size_t last1, const size_t last2) { - const size_t maxLength1 = energy.getAccessibility1().getMaxLength(); - const size_t maxLength2 = energy.getAccessibility2().getMaxLength(); + std::array sides; + Boundary bounds = getBoundary(interaction); + buildCandidates(sides[0], bounds, true, last1, last2); + buildCandidates(sides[1], bounds, false, last1, last2); while (true) { retain(interaction); - const Boundary bounds = getBoundary(interaction); - const size_t remaining1 = maxLength1 - (bounds[1] - bounds[0] + 1); - const size_t remaining2 = maxLength2 - (bounds[3] - bounds[2] + 1); - bool found = false; - Candidate best; - for (unsigned int side = 0; side < 2; ++side) { - const bool left = side == 0; - const size_t space1 = std::min(remaining1, left ? bounds[0] : last1 - bounds[1]); - const size_t space2 = std::min(remaining2, left ? bounds[2] : last2 - bounds[3]); - if (space1 == 0 || space2 == 0) { - continue; - } - // A GU boundary can still grow by stacking, but cannot start a - // nonstacking loop when either active constraint forbids it. - const bool stackOnly = (output.getOutputConstraint().noGUend - || !energy.isInternalLoopGUallowed()) - && energy.isGU(bounds[left ? 0 : 1], bounds[left ? 2 : 3]); - const size_t maxGap1 = stackOnly ? 0 - : std::min(energy.getMaxInternalLoopSize1(), space1 - 1); - const size_t maxGap2 = stackOnly ? 0 - : std::min(energy.getMaxInternalLoopSize2(), space2 - 1); - for (size_t s1 = 0; s1 <= maxGap1; ++s1) { - for (size_t s2 = 0; s2 <= maxGap2; ++s2) { - Candidate candidate; - candidate.left = left; - candidate.s1 = s1; - candidate.s2 = s2; - candidate.macro = output.getOutputConstraint().noLP && (s1 != 0 || s2 != 0); - const size_t addedPairs = candidate.macro ? 2 : 1; - if (space1 < addedPairs || space2 < addedPairs - || s1 > space1 - addedPairs || s2 > space2 - addedPairs) { - continue; - } - candidate.bounds = bounds; - if (left) { - candidate.bounds[0] -= s1 + addedPairs; - candidate.bounds[2] -= s2 + addedPairs; - candidate.close1 = bounds[0] - s1 - 1; - candidate.close2 = bounds[2] - s2 - 1; - } else { - candidate.bounds[1] += s1 + addedPairs; - candidate.bounds[3] += s2 + addedPairs; - candidate.close1 = bounds[1] + s1 + 1; - candidate.close2 = bounds[3] + s2 + 1; - } - if (evaluate(candidate, bounds, hybrid, interaction.energy) - && (!found || isBetter(candidate, best))) { - best = candidate; - found = true; - } - } - } - } - if (!found) { + const Candidate * left = updateCandidates(sides[0], bounds, hybrid, interaction.energy); + const Candidate * right = updateCandidates(sides[1], bounds, hybrid, interaction.energy); + if (left == NULL && right == NULL) { break; } + const Candidate best = left != NULL && (right == NULL || isBetter(*left, *right)) ? *left : *right; const Interaction::BasePair close = energy.getBasePair(best.close1, best.close2); if (best.left) { interaction.basePairs.insert(interaction.basePairs.begin(), close); @@ -276,47 +178,128 @@ PredictorSeedExtensionKinetic::extendSeed(Interaction & interaction, } hybrid = best.hybrid; interaction.energy = best.total; + bounds = best.bounds; + // The opposite end keeps its geometry, pair checks and loop energies. + // Its full energy must still be refreshed (ED and BOTH dangles change). + buildCandidates(sides[best.left ? 0 : 1], bounds, best.left, last1, last2); } } ////////////////////////////////////////////////////////////////////////// -bool -PredictorSeedExtensionKinetic::evaluate(Candidate & candidate, - const Boundary & bounds, const E_type hybrid, const E_type total) const +void +PredictorSeedExtensionKinetic::buildCandidates(SideCandidates & side, + const Boundary & bounds, const bool left, const size_t last1, const size_t last2) const { - const size_t old1 = bounds[candidate.left ? 0 : 1]; - const size_t old2 = bounds[candidate.left ? 2 : 3]; - const size_t outer1 = candidate.bounds[candidate.left ? 0 : 1]; - const size_t outer2 = candidate.bounds[candidate.left ? 2 : 3]; - if (!energy.areComplementary(candidate.close1, candidate.close2) - || (candidate.macro && !energy.areComplementary(outer1, outer2))) { - return false; - } - if (output.getOutputConstraint().noGUend && (candidate.s1 != 0 || candidate.s2 != 0) - && (energy.isGU(old1, old2) || energy.isGU(candidate.close1, candidate.close2))) { - return false; + side.moves.clear(); + const size_t space1 = std::min(energy.getAccessibility1().getMaxLength() - (bounds[1]-bounds[0]+1), + left ? bounds[0] : last1-bounds[1]); + const size_t space2 = std::min(energy.getAccessibility2().getMaxLength() - (bounds[3]-bounds[2]+1), + left ? bounds[2] : last2-bounds[3]); + if (space1 == 0 || space2 == 0) { + return; } - const E_type loop = candidate.left - ? energy.getE_interLeft(candidate.close1, old1, candidate.close2, old2) - : energy.getE_interLeft(old1, candidate.close1, old2, candidate.close2); - candidate.hybrid = addEnergy(hybrid, loop); - if (candidate.macro && E_isNotINF(candidate.hybrid)) { - const E_type stack = candidate.left - ? energy.getE_interLeft(outer1, candidate.close1, outer2, candidate.close2) - : energy.getE_interLeft(candidate.close1, outer1, candidate.close2, outer2); - candidate.hybrid = addEnergy(candidate.hybrid, stack); + const bool stackOnly = (output.getOutputConstraint().noGUend || !energy.isInternalLoopGUallowed()) + && energy.isGU(bounds[left ? 0 : 1], bounds[left ? 2 : 3]); + const size_t maxGap1 = space1 < 2 || stackOnly ? 0 : std::min(energy.getMaxInternalLoopSize1(), space1-2); + const size_t maxGap2 = space2 < 2 || stackOnly ? 0 : std::min(energy.getMaxInternalLoopSize2(), space2-2); + side.columns = maxGap2+2; + side.complementary.assign((maxGap1+2)*side.columns, -1); + const auto append = [&](size_t s1, size_t s2, bool macro) { + Candidate c; + c.left = left; c.s1 = s1; c.s2 = s2; c.macro = macro; + c.bounds = bounds; + const size_t pairs = macro ? 2 : 1; + if (left) { + c.close1 = bounds[0]-s1-1; c.close2 = bounds[2]-s2-1; + c.bounds[0] -= s1+pairs; c.bounds[2] -= s2+pairs; + } else { + c.close1 = bounds[1]+s1+1; c.close2 = bounds[3]+s2+1; + c.bounds[1] += s1+pairs; c.bounds[3] += s2+pairs; + } + side.moves.push_back(c); + }; + append(0, 0, false); + if (space1 >= 2 && space2 >= 2) { + for (size_t s1 = 0; s1 <= maxGap1; ++s1) { + for (size_t s2 = 0; s2 <= maxGap2; ++s2) { + append(s1, s2, true); + } + } } - if (E_isINF(candidate.hybrid)) { - return false; +} + +////////////////////////////////////////////////////////////////////////// + +const PredictorSeedExtensionKinetic::Candidate * +PredictorSeedExtensionKinetic::updateCandidates(SideCandidates & side, + const Boundary & bounds, const E_type hybrid, const E_type total) const +{ + // Phase one: each position pair is tested at most once per unchanged end, + // even when it is the closing pair of one move and outer pair of another. + size_t stopGap2 = std::numeric_limits::max(); + for (Candidate & c : side.moves) { + c.bounds[c.left ? 1 : 0] = bounds[c.left ? 1 : 0]; + c.bounds[c.left ? 3 : 2] = bounds[c.left ? 3 : 2]; + c.active = c.bounds[1]-c.bounds[0]+1 <= energy.getAccessibility1().getMaxLength() + && c.bounds[3]-c.bounds[2]+1 <= energy.getAccessibility2().getMaxLength(); + if (c.active && c.macro) { + if (c.s2 >= stopGap2) { + c.active = false; + } else if (prune(c, bounds)) { + // A suffix bound rejects this rectangle of larger gaps without + // further ED, complementarity or loop-energy lookups. + stopGap2 = c.s2; + c.active = false; + } + } + if (!c.active || c.topologyKnown) { + continue; + } + const auto complementary = [&](size_t s1, size_t s2) { + signed char & cached = side.complementary[s1*side.columns+s2]; + if (cached < 0) { + cached = energy.areComplementary(c.left ? bounds[0]-s1-1 : bounds[1]+s1+1, + c.left ? bounds[2]-s2-1 : bounds[3]+s2+1); + } + return cached != 0; + }; + c.topologyKnown = true; + c.topologyAllowed = complementary(c.s1, c.s2) + && (!c.macro || complementary(c.s1+1, c.s2+1)); + if (c.topologyAllowed && output.getOutputConstraint().noGUend && (c.s1 != 0 || c.s2 != 0)) { + c.topologyAllowed = !energy.isGU(c.close1, c.close2); + } } - const Boundary & b = candidate.bounds; - candidate.total = energy.getE(b[0], b[1], b[2], b[3], candidate.hybrid); - if (E_isINF(candidate.total)) { - return false; + const Candidate * best = NULL; + for (Candidate & c : side.moves) { + if (!c.active || !c.topologyAllowed) { + continue; + } + if (!c.localKnown) { + c.localKnown = true; + c.local = c.left ? energy.getE_interLeft(c.close1, bounds[0], c.close2, bounds[2]) + : energy.getE_interLeft(bounds[1], c.close1, bounds[3], c.close2); + if (c.macro && E_isNotINF(c.local)) { + c.local = addEnergy(c.local, c.left + ? energy.getE_interLeft(c.bounds[0], c.close1, c.bounds[2], c.close2) + : energy.getE_interLeft(c.close1, c.bounds[1], c.close2, c.bounds[3])); + } + } + c.hybrid = addEnergy(hybrid, c.local); + if (E_isINF(c.hybrid)) { + continue; + } + c.total = energy.getE(c.bounds[0], c.bounds[1], c.bounds[2], c.bounds[3], c.hybrid); + if (E_isINF(c.total)) { + continue; + } + c.delta = std::int64_t(c.total)-std::int64_t(total); + if (c.delta < 0 && (best == NULL || isBetter(c, *best))) { + best = &c; + } } - candidate.delta = std::int64_t(candidate.total) - std::int64_t(total); - return candidate.delta < 0; + return best; } ////////////////////////////////////////////////////////////////////////// @@ -346,7 +329,10 @@ PredictorSeedExtensionKinetic::isBetter(const Candidate & candidate, const Candi } const Wide size = Wide(candidate.s1) + Wide(candidate.s2); const Wide bestSize = Wide(best.s1) + Wide(best.s2); - return size != bestSize ? size < bestSize : candidate.s1 < best.s1; + if (size != bestSize) return size < bestSize; + if (candidate.s1 != best.s1) return candidate.s1 < best.s1; + // Identical shape/score: retain the shorter move first. + return !candidate.macro && best.macro; } ////////////////////////////////////////////////////////////////////////// @@ -386,19 +372,6 @@ PredictorSeedExtensionKinetic::traceBack(Interaction & interaction) } interaction = path->second; seedHandler.addSeeds(interaction); - // The generic annotator can recognize explicit seeds that were rejected - // as starting states (e.g. a lonely seed end stacked only by extension). - // Retain only validated starts and use their reconstructed energies. - if (interaction.seed != NULL) { - Interaction::SeedSet validAnnotations; - for (const Interaction::Seed & seed : *interaction.seed) { - const auto valid = validSeeds.find(seed.bp_i); - if (valid != validSeeds.end() && valid->second.bp_j == seed.bp_j) { - validAnnotations.insert(valid->second); - } - } - *interaction.seed = std::move(validAnnotations); - } } ////////////////////////////////////////////////////////////////////////// diff --git a/src/IntaRNA/PredictorSeedExtensionKinetic.h b/src/IntaRNA/PredictorSeedExtensionKinetic.h index e9c3f3d7..4af54a87 100644 --- a/src/IntaRNA/PredictorSeedExtensionKinetic.h +++ b/src/IntaRNA/PredictorSeedExtensionKinetic.h @@ -7,6 +7,7 @@ #include #include #include +#include namespace IntaRNA { @@ -20,9 +21,10 @@ namespace IntaRNA { * not physical rates or a calibrated folding-time model. Only negative * energy differences are accepted, including for the normalized scores. * - * With noLP, every starting seed must already contain no lonely pair and a - * nonstacking extension adds its closing pair and the following stack - * atomically. Ties prefer left, smaller s1+s2, then smaller s1. Every valid + * Extensions always add one stacked pair or two stacked pairs, including + * across a loop. Seeds and their energies are trusted as supplied by the seed + * handler, even when an explicit seed contains lonely pairs. Ties prefer left, + * smaller s1+s2, smaller s1, then the single-pair move. Every valid * visited prefix is eligible for normal MFE/suboptimal reporting; traceback * reproduces the actual greedy path. Equilibrium partition-function output * is unsupported. @@ -76,7 +78,6 @@ class PredictorSeedExtensionKinetic : public PredictorMfe { */ void getNextBest(Interaction & interaction) override; -private: //! Inclusive boundaries (i1,j1,i2,j2), using local energy indices. using Boundary = std::array; //! Best actual path for each visited, reportable set of boundaries. @@ -91,31 +92,43 @@ class PredictorSeedExtensionKinetic : public PredictorMfe { size_t s2 = 0; bool left = true; bool macro = false; + bool topologyKnown = false; + bool topologyAllowed = false; + bool active = false; + bool localKnown = false; + E_type local = E_INF; E_type hybrid = E_INF; E_type total = E_INF; std::int64_t delta = 0; }; + /** + * Optional filter before pairing and loop-energy lookups. The default + * enumerates every move; subclasses may implement heuristic pruning. + * @param candidate geometrically valid extension + * @param bounds current boundaries + * @return whether to skip this two-pair move and all moves with both gaps + * at least as large, in the current state (requires monotone ED) + */ + virtual bool prune(const Candidate & candidate, const Boundary & bounds) const; + +private: + /** Geometry, shared pair checks and local energies for one unchanged end. */ + struct SideCandidates { + std::vector moves; + std::vector complementary; + size_t columns = 0; + }; + //! Owned seed handler with offsets matching this->energy. SeedHandlerIdxOffset seedHandler; //! The selected deterministic move-ranking formula. const char score; //! Paths retained independently of the number of requested reports. InteractionCache interactions; - //! Valid starting seeds with recomputed energies, keyed by original indices. - std::map validSeeds; - /** @return local, inclusive boundaries of a nonempty interaction */ Boundary getBoundary(const Interaction & interaction) const; - /** - * Checks the complete seed path, including noLP and loop GU constraints. - * @param interaction seed path, in original sequence coordinates - * @param hybrid receives the traced path's loop energies plus initiation - * @return whether every pair and adjacent loop is feasible - */ - bool isValidSeed(const Interaction & interaction, E_type & hybrid) const; - /** * Runs a complete greedy trajectory and retains its reportable prefixes. * @param interaction initial seed, modified to its final state @@ -127,14 +140,26 @@ class PredictorSeedExtensionKinetic : public PredictorMfe { size_t last1, size_t last2); /** - * Checks one geometrically bounded move and evaluates its complete energy. - * @param candidate move indices/side/gaps; energies are filled on success - * @param bounds current interaction boundaries - * @param hybrid current hybridization energy, including initiation - * @param total current full interaction energy - * @return whether the entire move is feasible and strictly downhill + * Builds the rectangular gap table plus the single-stack move for one end. + * @param side table to replace + * @param bounds current boundaries + * @param left whether this is the left end + * @param last1 last permitted local index in sequence 1 + * @param last2 last permitted local index in reversed sequence 2 + */ + void buildCandidates(SideCandidates & side, const Boundary & bounds, + bool left, size_t last1, size_t last2) const; + + /** + * Resolves shared pair checks before energies, then refreshes complete + * energies and selects the best downhill candidate during that traversal. + * @param side candidate and pairing cache for one end + * @param bounds current boundaries (the opposite end may have changed) + * @param hybrid current hybridization energy including initiation + * @param total current complete interaction energy + * @return best move within side, or NULL if none is downhill */ - bool evaluate(Candidate & candidate, const Boundary & bounds, + const Candidate * updateCandidates(SideCandidates & side, const Boundary & bounds, E_type hybrid, E_type total) const; /** @return whether a move wins by exact score and deterministic ties */ diff --git a/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp b/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp new file mode 100644 index 00000000..e35f9f53 --- /dev/null +++ b/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp @@ -0,0 +1,109 @@ +#include "IntaRNA/PredictorSeedExtensionKineticPruned.h" + +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/InteractionEnergyVrna.h" + +#include +#include +#include + +namespace IntaRNA { + +PredictorSeedExtensionKineticPruned::PredictorSeedExtensionKineticPruned( + const InteractionEnergy & model, OutputHandler & output, + PredictionTracker * predTracker, SeedHandler * seedHandler, const char score) + : PredictorSeedExtensionKinetic(model, output, predTracker, seedHandler, score) +{ + // Parameter-based bounds cannot describe arbitrary overrides of getE or + // getE_interLeft. Do not silently apply them to derived energy models. + const bool vrna = typeid(model) == typeid(InteractionEnergyVrna); + const bool basePair = typeid(model) == typeid(InteractionEnergyBasePair); + if ((!vrna && !basePair) || model.getMaxInternalLoopSize1() > 30 + || model.getMaxInternalLoopSize2() > 30) { + return; + } + rows = std::min(model.getMaxInternalLoopSize1(), model.size1() > 2 ? model.size1()-2 : 0)+1; + columns = std::min(model.getMaxInternalLoopSize2(), model.size2() > 2 ? model.size2()-2 : 0)+1; + // E_IntLoop takes a mutable pointer in supported ViennaRNA versions. + // Work on a copy and retain no pointer to the model's parameter storage. + vrna_param_t params; + if (vrna) params = static_cast(model).getVrnaParams(); + constexpr int reverse[] = {0,2,1,4,3,6,5}; + for (size_t side = 0; side < 2; ++side) { + for (int root = 1; root <= 6; ++root) { + auto & table = lowerBounds[side*6+root-1]; + table.assign(rows*columns, E_INF); + for (size_t s1 = 0; s1 < rows; ++s1) { + for (size_t s2 = 0; s2 < columns; ++s2) { + E_type & bound = table[s1*columns+s2]; + if (basePair) { + const std::int64_t two = 2*std::int64_t(model.getE_init()); + bound = static_cast(std::clamp(two, std::int64_t(std::numeric_limits::min()), std::int64_t(E_INF))); + continue; + } + for (int close = 1; close <= 6; ++close) { + E_type stack = E_INF; + for (int outer = 1; outer <= 6; ++outer) { + stack = std::min(stack, static_cast(side == 0 + ? params.stack[outer][reverse[close]] : params.stack[close][reverse[outer]])); + } + // Relax sequence consistency of the four mismatch bases. + // This can only lower the local loop/stack minimum. + for (int a = 1; a <= 4; ++a) for (int b = 1; b <= 4; ++b) + for (int c = 1; c <= 4; ++c) for (int d = 1; d <= 4; ++d) { + const E_type loop = E_IntLoop(s1, s2, + side == 0 ? close : root, reverse[side == 0 ? root : close], a,b,c,d, ¶ms); + if (E_isNotINF(loop) && E_isNotINF(stack)) { + const std::int64_t local = std::int64_t(loop)+stack; + bound = std::min(bound, static_cast(std::clamp(local, + std::int64_t(std::numeric_limits::min()), std::int64_t(E_INF)))); + } + } + } + } + } + // Componentwise suffix minima make the bound nondecreasing along + // either gap. With monotone ED this permits rectangle pruning. + for (size_t i = rows; i-- > 0;) for (size_t j = columns; j-- > 0;) { + E_type & value = table[i*columns+j]; + if (i+1 < rows) value = std::min(value, table[(i+1)*columns+j]); + if (j+1 < columns) value = std::min(value, table[i*columns+j+1]); + } + } + } +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKineticPruned::predict(const IndexRange & r1, const IndexRange & r2) +{ + edKnown = false; + PredictorSeedExtensionKinetic::predict(r1, r2); +} + +////////////////////////////////////////////////////////////////////////// + +bool +PredictorSeedExtensionKineticPruned::prune(const Candidate & candidate, const Boundary & bounds) const +{ + if (rows == 0) return false; + // Use global coordinates for the ED memo, including repeated predict() + // calls with different range offsets but identical local boundaries. + const Boundary global{bounds[0]+energy.getOffset1(), bounds[1]+energy.getOffset1(), + bounds[2]+energy.getOffset2(), bounds[3]+energy.getOffset2()}; + if (!edKnown || edBounds != global) { + edKnown = true; edBounds = global; + currentED = std::int64_t(energy.getED1(bounds[0],bounds[1])) + energy.getED2(bounds[2],bounds[3]); + } + const int root = BP_pair[energy.getAccessibility1().getSequence().asCodes().at(global[candidate.left ? 0 : 1])] + [energy.getAccessibility2().getSequence().asCodes().at(global[candidate.left ? 2 : 3])]; + if (root < 1 || root > 6 || candidate.s1 >= rows || candidate.s2 >= columns) return false; + const E_type bound = lowerBounds[(candidate.left ? 0 : 6)+root-1][candidate.s1*columns+candidate.s2]; + const Boundary & b = candidate.bounds; + const E_type ed1 = energy.getED1(b[0],b[1]), ed2 = energy.getED2(b[2],b[3]); + return E_isINF(bound) || ed1 >= Accessibility::ED_UPPER_BOUND || ed2 >= Accessibility::ED_UPPER_BOUND + || std::int64_t(bound)+ed1+ed2-currentED >= 0; +} + +} // namespace IntaRNA diff --git a/src/IntaRNA/PredictorSeedExtensionKineticPruned.h b/src/IntaRNA/PredictorSeedExtensionKineticPruned.h new file mode 100644 index 00000000..342fa302 --- /dev/null +++ b/src/IntaRNA/PredictorSeedExtensionKineticPruned.h @@ -0,0 +1,65 @@ +#ifndef INTARNA_PREDICTORSEEDEXTENSIONKINETICPRUNED_H_ +#define INTARNA_PREDICTORSEEDEXTENSIONKINETICPRUNED_H_ + +#include "IntaRNA/PredictorSeedExtensionKinetic.h" + +namespace IntaRNA { + +/** + * Experimental kinetic extension with root-pair-specific loop/stack bounds. + * + * Before complementarity checks, a suffix-minimum table and the current ED + * increment can reject a rectangle of larger two-pair moves. Local estimates + * use the active ViennaRNA parameters (or the base-pair energy). Pruning + * assumes monotone ED and ignores changes in terminal/dangling contributions: + * it is intentionally heuristic and can change the path compared with K. + * Every surviving move is still checked with its complete energy. Single + * stacks are always evaluated. Unknown energy subclasses and loop limits + * above 30 fall back to exhaustive K enumeration. + */ +class PredictorSeedExtensionKineticPruned : public PredictorSeedExtensionKinetic { +public: + /** + * Constructs the experimental predictor and precomputes local bounds. + * @param energy energy model that must outlive this predictor + * @param output output handler that must outlive this predictor + * @param predTracker owned prediction tracker, or NULL + * @param seedHandler owned, non-NULL seed handler + * @param score move ranking A, B or C + */ + PredictorSeedExtensionKineticPruned(const InteractionEnergy & energy, + OutputHandler & output, PredictionTracker * predTracker, + SeedHandler * seedHandler, char score = 'A'); + + /** + * Resets the ED memo and runs kinetic prediction within the given ranges. + * @param r1 permitted inclusive range in sequence 1 + * @param r2 permitted inclusive range in reversed sequence 2 + */ + void predict(const IndexRange & r1 = IndexRange(0, RnaSequence::lastPos), + const IndexRange & r2 = IndexRange(0, RnaSequence::lastPos)) override; + +protected: + /** + * Tests the local suffix bound plus ED increase, omitting terminal and + * dangling changes. A true result rejects all componentwise larger gaps. + * @param candidate geometrically valid two-pair extension + * @param bounds current interaction boundaries + * @return whether to prune the candidate and its larger-gap rectangle + */ + bool prune(const Candidate & candidate, const Boundary & bounds) const override; + +private: + //! Left/right tables for the six oriented canonical base-pair types. + std::array, 12> lowerBounds; + //! Rectangular gap-table dimensions; zero disables the heuristic. + size_t rows = 0, columns = 0; + //! Current-state ED is shared by all candidates at both ends. + mutable Boundary edBounds = {}; + mutable bool edKnown = false; + mutable std::int64_t currentED = 0; +}; + +} // namespace IntaRNA + +#endif /* INTARNA_PREDICTORSEEDEXTENSIONKINETICPRUNED_H_ */ diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index d39c2b0c..0ecf9206 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -51,6 +51,7 @@ extern "C" { #include "IntaRNA/PredictorMfe2dSeedExtension.h" #include "IntaRNA/PredictorMfe2dSeedExtensionRIblast.h" #include "IntaRNA/PredictorSeedExtensionKinetic.h" +#include "IntaRNA/PredictorSeedExtensionKineticPruned.h" #include "IntaRNA/PredictorMfe2dHeuristicSeedExtension.h" #include "IntaRNA/PredictorMfeEnsSeedOnly.h" @@ -184,7 +185,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) temperature("temperature",0,100,37), model("model", "SPBX", 'X'), - mode("mode", "HMSRK", 'H'), // R for RIblast heuristic only + mode("mode", "HMSRKL", 'H'), // R for RIblast heuristic only kineticScore("kineticScore", "ABC", 'A'), #if INTARNA_MULITHREADING threads("threads", 0, omp_get_max_threads(), 1), @@ -806,13 +807,14 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) "\n 'H' = heuristic (fast and low memory), " "\n 'M' = exact (slow), " "\n 'S' = seed-only, " - "\n 'K' = downhill greedy seed extension (requires --model=X; no time or rate prediction)" + "\n 'K' = downhill greedy seed extension (requires --model=X; always noLP extensions), " + "\n 'L' = experimental K with local-energy/monotone-ED pruning (may change paths; no time or rate prediction)" ).c_str()) (kineticScore.name.c_str() , value(&(kineticScore.val)) ->default_value(kineticScore.def) ->notifier(boost::bind(&CommandLineParsing::validate_charArgument,this,kineticScore,_1)) - , "candidate score for --model=X --mode=K: 'A' = complete interaction energy change, " + , "candidate score for --model=X --mode=K or L: 'A' = complete interaction energy change, " "'B' = change/(1+s1+s2), 'C' = change/(1+2*max(s1,s2)), where s1/s2 are skipped bases. " "All modes accept strictly negative energy changes only; B/C are heuristic scores.") (model.name.c_str() @@ -1232,10 +1234,14 @@ parse(int argc, char** argv) // parsing escape literals outSep = unescaped_string::getUnescaped( outSep ); - // K needs a seed even before the usual --noSeed model normalization. - if (mode.val == 'K') { - if (model.val != 'X') throw error("--mode=K is available only with --model=X"); - if (noSeedRequired) throw error("--mode=K requires seeds and is incompatible with --noSeed"); + // K/L need a seed even before the usual --noSeed model normalization. + if (mode.val == 'K' || mode.val == 'L') { + if (model.val != 'X') throw error("--mode=K/L is available only with --model=X"); + if (noSeedRequired) throw error("--mode=K/L requires seeds and is incompatible with --noSeed"); + if (!outNoLP) { + LOG(INFO) << "--mode=" << mode.val << " uses no-lonely-pair extensions: setting --outNoLP=true (handler-provided seeds are unchanged)"; + outNoLP = true; + } // A selected set of greedy paths is not an equilibrium ensemble. if (outMode.val == 'E' || (outMode.val == 'C' && OutputHandlerCsv::needsZall(OutputHandlerCsv::string2list(outCsvCols))) @@ -1244,10 +1250,10 @@ parse(int argc, char** argv) || !outPrefix2streamName.at(OutPrefixCode::OP_qSpotProb).empty() || !outPrefix2streamName.at(OutPrefixCode::OP_tSpotProb).empty()) { - throw error("--mode=K does not support equilibrium ensemble or interaction-probability output"); + throw error("--mode=K/L does not support equilibrium ensemble or interaction-probability output"); } } else if (vm.count(kineticScore.name) && !vm.at(kineticScore.name).defaulted()) { - throw error("--kineticScore requires --model=X --mode=K"); + throw error("--kineticScore requires --model=X --mode=K or L"); } // open output stream @@ -2583,6 +2589,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const case 'M' : return new PredictorMfe2dSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); case 'R' : return new PredictorMfe2dSeedExtensionRIblast( energy, output, predTracker, getSeedHandler( energy ) ); case 'K' : return new PredictorSeedExtensionKinetic( energy, output, predTracker, getSeedHandler( energy ), kineticScore.val ); + case 'L' : return new PredictorSeedExtensionKineticPruned( energy, output, predTracker, getSeedHandler( energy ), kineticScore.val ); case 'S' : return new PredictorMfeSeedOnly( energy, output, predTracker, getSeedHandler( energy ) ); default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented"); return NULL; } diff --git a/tests/PredictorSeedExtensionKinetic_test.cpp b/tests/PredictorSeedExtensionKinetic_test.cpp index cbd47d4d..6e2ef0b1 100644 --- a/tests/PredictorSeedExtensionKinetic_test.cpp +++ b/tests/PredictorSeedExtensionKinetic_test.cpp @@ -7,6 +7,7 @@ #include "IntaRNA/InteractionEnergyVrna.h" #include "IntaRNA/OutputHandler.h" #include "IntaRNA/PredictorSeedExtensionKinetic.h" +#include "IntaRNA/PredictorSeedExtensionKineticPruned.h" #include "IntaRNA/SeedHandlerExplicit.h" #include "IntaRNA/VrnaHandler.h" @@ -67,6 +68,7 @@ class KineticEnergy final : public InteractionEnergyBasePair { bool customLoops = false; std::map loops; std::map boundaryTerms; + mutable std::map loopCalls; }; KineticEnergy::KineticEnergy(const Accessibility & a, const ReverseAccessibility & b, @@ -74,6 +76,7 @@ KineticEnergy::KineticEnergy(const Accessibility & a, const ReverseAccessibility : InteractionEnergyBasePair(a, b, m1, m2, false, 1., -100, 3, 0, false) {} E_type KineticEnergy::getE_interLeft(size_t i, size_t j, size_t k, size_t l) const { + ++loopCalls[{i,j,k,l}]; if (!isValidInternalLoop(i, j, k, l)) return E_INF; if (!customLoops) return InteractionEnergyBasePair::getE_interLeft(i, j, k, l); auto p = loops.find({i,j,k,l}); @@ -118,12 +121,7 @@ E_type chainEnergy(const InteractionEnergy & energy, const Chain & chain, E_type h = energy.getE_init(); for (size_t p = 0; p < chain.size(); ++p) { if (!energy.areComplementary(chain[p].first,chain[p].second)) return E_INF; - const bool before = p != 0 && stacked(chain[p-1],chain[p]); - const bool after = p+1 < chain.size() && stacked(chain[p],chain[p+1]); - if (out.noLP && !before && !after) return E_INF; if (p == 0) continue; - if (out.noGUend && !before && (energy.isGU(chain[p-1].first,chain[p-1].second) - || energy.isGU(chain[p].first,chain[p].second))) return E_INF; const E_type loop = energy.getE_interLeft(chain[p-1].first,chain[p].first, chain[p-1].second,chain[p].second); if (E_isINF(loop)) return E_INF; @@ -146,7 +144,7 @@ std::vector oracle(const InteractionEnergy & energy, Chain chain, && energy.getED2(b[2],b[3]) <= out.maxED) result.push_back(chain); bool found = false; Chain best; - std::tuple bestKey; + std::tuple bestKey; // Enumerate absolute endpoints in deliberately different order from // production's side/loop-size enumeration. for (size_t x = r1.from; x <= r1.to; ++x) { @@ -157,20 +155,25 @@ std::vector oracle(const InteractionEnergy & energy, Chain chain, const size_t s1 = left ? b[0]-x-1 : x-b[1]-1; const size_t s2 = left ? b[2]-y-1 : y-b[3]-1; if (s1 > energy.getMaxInternalLoopSize1() || s2 > energy.getMaxInternalLoopSize2()) continue; - Chain trial = chain; - trial.push_back({x,y}); - if (out.noLP && s1+s2 != 0) { - if (left && (x == r1.from || y == r2.from)) continue; - if (right && (x == r1.to || y == r2.to)) continue; - trial.push_back(left ? Pair{x-1,y-1} : Pair{x+1,y+1}); + if (s1+s2 != 0 && (out.noGUend || !energy.isInternalLoopGUallowed()) + && (energy.isGU(left ? b[0] : b[1], left ? b[2] : b[3]) || energy.isGU(x,y))) continue; + for (bool macro : {false,true}) { + if (!macro && s1+s2 != 0) continue; + Chain trial = chain; + trial.push_back({x,y}); + if (macro) { + if (left && (x == r1.from || y == r2.from)) continue; + if (right && (x == r1.to || y == r2.to)) continue; + trial.push_back(left ? Pair{x-1,y-1} : Pair{x+1,y+1}); + } + std::sort(trial.begin(),trial.end()); + const E_type total = chainEnergy(energy,trial,out); + if (E_isINF(total) || total >= current) continue; + const size_t denominator = score == 'A' ? 1 : score == 'B' ? 1+s1+s2 : 1+2*std::max(s1,s2); + const auto key = std::make_tuple(static_cast(total-current)/denominator, + left ? 0 : 1,s1+s2,s1,macro); + if (!found || key < bestKey) { found = true; bestKey = key; best = trial; } } - std::sort(trial.begin(),trial.end()); - const E_type total = chainEnergy(energy,trial,out); - if (E_isINF(total) || total >= current) continue; - const size_t denominator = score == 'A' ? 1 : score == 'B' ? 1+s1+s2 : 1+2*std::max(s1,s2); - const auto key = std::make_tuple(static_cast(total-current)/denominator, - left ? 0 : 1,s1+s2,s1); - if (!found || key < bestKey) { found = true; bestKey = key; best = trial; } } } if (!found) return result; @@ -256,8 +259,8 @@ TEST_CASE("Kinetic scoring and strict downhill acceptance", "[PredictorSeedExten SECTION("A B and C select distinct justified moves") { KineticFixture f(8); f.energy.customLoops = true; - f.energy.loops = {{{2,3,2,3},-100},{{3,6,3,4},-420},{{3,5,3,5},-390},{{3,4,3,4},-120}}; - const std::map expected{{'A',{6,4}},{'B',{6,4}},{'C',{5,5}}}; + f.energy.loops = {{{2,3,2,3},-100},{{3,6,3,4},-420},{{3,5,3,5},-390},{{3,4,3,4},-120},{{6,7,4,5},0},{{5,6,5,6},0}}; + const std::map expected{{'A',{7,5}},{'B',{7,5}},{'C',{6,6}}}; for (char score : {'A','B','C'}) { checkOracle(f.energy,seed,out,score); const auto actual = predict(f.energy,seed,out,score); @@ -286,15 +289,17 @@ TEST_CASE("Kinetic scoring and strict downhill acceptance", "[PredictorSeedExten KineticFixture shortF(7,3,3,3,3); shortF.energy.customLoops = true; shortF.energy.loops = f.energy.loops; REQUIRE(internalChain(shortF.energy,predict(shortF.energy,seed,out).front()).front() == Pair(1,1)); - f.energy.loops.erase({1,2,1,2}); - // Equal total gap: s1=0 (pair 1,0) wins over s1=1 (pair 0,1). - const auto paths = oracle(f.energy,seed,out,'A',IndexRange(0,6),IndexRange(0,6)); + f.energy.loops = {{{3,4,3,4},-100},{{2,3,1,3},-50},{{1,2,0,1},-50}, + {{1,3,2,3},-50},{{0,1,1,2},-50}}; + const Chain middle{{3,3},{4,4}}; + // Equal total gap: s1=0 (outer pair 1,0) wins over s1=1 (0,1). + const auto paths = oracle(f.energy,middle,out,'A',IndexRange(0,6),IndexRange(0,6)); REQUIRE(paths.at(1).front() == Pair(1,0)); - checkOracle(f.energy,seed,out,'A'); + checkOracle(f.energy,middle,out,'A'); } } -TEST_CASE("Kinetic noLP macro-steps and explicit seed validation", "[PredictorSeedExtensionKinetic]") { +TEST_CASE("Kinetic always uses stacked extensions and trusts explicit seeds", "[PredictorSeedExtensionKinetic]") { #include "testEasyLoggingSetup.icc" KineticFixture f(6,1,1); f.energy.customLoops = true; @@ -321,8 +326,8 @@ TEST_CASE("Kinetic noLP macro-steps and explicit seed validation", "[PredictorSe REQUIRE(internalChain(ranked.energy,rankedResult.front()).back() == Pair(6,6)); } } - SECTION("lonely explicit seeds are skipped") { - REQUIRE(predict(f.energy,Chain{{1,1},{3,3}},out).empty()); + SECTION("lonely explicit seeds are accepted") { + REQUIRE_FALSE(predict(f.energy,Chain{{1,1},{3,3}},out).empty()); } SECTION("both additional pairs must fit both strand spans") { KineticFixture shortF(6,1,1,4,6); @@ -425,12 +430,12 @@ TEST_CASE("Kinetic nonoverlap reporting recovers shorter prefixes", "[PredictorS predictor.predict(); REQUIRE(output.interactions.size() == 2); REQUIRE(output.interactions[0].energy == -900); - REQUIRE(output.interactions[1].energy == -300); + REQUIRE(output.interactions[1].energy == -200); REQUIRE((bounds(internalChain(f.energy,output.interactions[0])) == Bounds{{3,7,3,7}})); - REQUIRE((bounds(internalChain(f.energy,output.interactions[1])) == Bounds{{0,2,0,2}})); + REQUIRE((bounds(internalChain(f.energy,output.interactions[1])) == Bounds{{0,1,0,1}})); if (needBPs) { REQUIRE(output.interactions[0].basePairs.size() == 5); - REQUIRE(output.interactions[1].basePairs.size() == 3); + REQUIRE(output.interactions[1].basePairs.size() == 2); } // Reuse the instance with a different range and ensure old candidates // and seed metadata cannot leak into the second prediction. @@ -442,7 +447,7 @@ TEST_CASE("Kinetic nonoverlap reporting recovers shorter prefixes", "[PredictorS } } -TEST_CASE("Kinetic annotations exclude rejected lonely explicit seeds", "[PredictorSeedExtensionKinetic]") { +TEST_CASE("Kinetic annotations retain handler-provided lonely explicit seeds", "[PredictorSeedExtensionKinetic]") { #include "testEasyLoggingSetup.icc" KineticFixture f(5,0,0); const auto sc = seedConstraint("2||&3||,1|&5|"); @@ -454,5 +459,182 @@ TEST_CASE("Kinetic annotations exclude rejected lonely explicit seeds", "[Predic const auto & result = output.interactions.front(); REQUIRE(result.basePairs.size() == 5); REQUIRE(result.seed != NULL); - REQUIRE(result.seed->size() == 1); + REQUIRE(result.seed->size() == 2); +} + +TEST_CASE("Kinetic two-stack moves cross a local barrier atomically", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(4,0,0); + f.energy.customLoops = true; + f.energy.loops = {{{0,1,0,1},-100},{{1,2,1,2},50},{{2,3,2,3},-200}}; + const Chain seed{{0,0},{1,1}}; + for (bool noLP : {false,true}) { + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP); + checkOracle(f.energy,seed,out,'A'); + const auto actual = predict(f.energy,seed,out); + REQUIRE(actual.size() == 2); + REQUIRE(actual.front().basePairs.size() == 4); + REQUIRE(actual.front().energy == -350); + } +} + +TEST_CASE("Kinetic exact move ties prefer the single stack", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(4,0,0); + f.energy.customLoops = true; + f.energy.loops = {{{0,1,0,1},-100},{{1,2,1,2},-100},{{2,3,2,3},0}}; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + for (char score : {'A','B','C'}) { + const auto actual = predict(f.energy,Chain{{0,0},{1,1}},out,score); + REQUIRE(actual.size() == 2); + REQUIRE(actual.front().basePairs.size() == 3); + REQUIRE(actual.front().energy == -300); + } +} + +TEST_CASE("Kinetic caches the unchanged end and trusts the seed energy", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(8,0,0); + f.energy.customLoops = true; + for (size_t i = 0; i+1 < 8; ++i) f.energy.loops[{i,i+1,i,i+1}] = i < 2 ? -200 : -100; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + const auto result = predict(f.energy,Chain{{2,2},{3,3}},out); + REQUIRE(result.front().basePairs.size() == 8); + // Seed energy is evaluated by SeedHandlerExplicit just once, never by + // the predictor. The right first stack is shared by single/two-pair moves + // and is not reevaluated after the left end grows. + REQUIRE(f.energy.loopCalls.at(Bounds{{2,3,2,3}}) == 1); + REQUIRE(f.energy.loopCalls.at(Bounds{{3,4,3,4}}) == 2); +} + +TEST_CASE("Kinetic oracle covers varied sequences and endpoint energies", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + unsigned int random = 254; + const auto next = [&]() { random = random*1664525u+1013904223u; return random; }; + for (size_t sample = 0; sample < 32; ++sample) { + std::string t(9,'G'), q(9,'C'); + for (size_t i = 0; i < 9; ++i) { t[i] = "ACGU"[(next() >> 16)%4]; q[i] = "ACGU"[(next() >> 16)%4]; } + t[3] = t[4] = 'G'; q[4] = q[5] = 'C'; + RnaSequence first("t",t), second("q",q); + KineticAccessibility a(first), b(second); + for (size_t i = 0; i < 9; ++i) for (size_t j = i; j < 9; ++j) { + a.values[{i,j}] = (next() >> 16)%300; + b.values[{i,j}] = (next() >> 16)%300; + } + ReverseAccessibility reversed(b); + InteractionEnergyBasePair energy(a,reversed,2,3); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,sample%2,sample%3 == 0); + for (char score : {'A','B','C'}) checkOracle(energy,Chain{{3,3},{4,4}},out,score); + } +} + +TEST_CASE("Pruned kinetic mode preserves complete energies and resets offsets", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + RnaSequence t("t","GGGGGGGG"), q("q","CCCCCCCC"); + KineticAccessibility a(t), b(q); + for (size_t i = 0; i < 8; ++i) for (size_t j = i; j < 8; ++j) { + a.values[{i,j}] = 20*(j-i); b.values[{i,j}] = 25*(j-i); + } + ReverseAccessibility reversed(b); + VrnaHandler vrna(25,"Turner99",false,false); + InteractionEnergyVrna energy(a,reversed,vrna,2,3,false,37,false); + const Chain seed{{3,3},{4,4}}; + const auto sc = seedConstraint(seedEncoding(seed,8)); + for (char score : {'A','B','C'}) { + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + KineticOutput output(out); + PredictorSeedExtensionKineticPruned predictor(energy,output,NULL,new SeedHandlerExplicit(energy,sc),score); + for (const IndexRange & range : {IndexRange(0,7),IndexRange(1,6),IndexRange(2,7),IndexRange(0,7)}) { + output.interactions.clear(); predictor.predict(range,range); + const auto expected = predict(energy,seed,out,score,range,range); + REQUIRE(output.interactions.size() == expected.size()); + for (size_t i = 0; i < expected.size(); ++i) { + REQUIRE(output.interactions[i].basePairs == expected[i].basePairs); + REQUIRE(output.interactions[i].energy == chainEnergy(energy,internalChain(energy,expected[i]),out)); + } + } + } +} + +TEST_CASE("Pruned kinetic mode documents its nonmonotone ED tradeoff", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(6,2,2); + f.acc1.values = {{{0,2},200},{{0,3},300}}; + InteractionEnergyBasePair energy(f.acc1,f.reversed,2,2); + const Chain seed{{0,0},{1,1}}; + const auto sc = seedConstraint(seedEncoding(seed,6)); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + KineticOutput output(out); + PredictorSeedExtensionKineticPruned predictor(energy,output,NULL,new SeedHandlerExplicit(energy,sc)); + predictor.predict(); + REQUIRE(output.interactions.size() == 1); + REQUIRE(output.interactions.front().basePairs.size() == 2); + REQUIRE(predict(energy,seed,out).front().basePairs.size() > 2); + // The subclass has no bound for custom energy overrides: safely use K. + KineticOutput customOutput(out); + PredictorSeedExtensionKineticPruned custom(f.energy,customOutput,NULL,new SeedHandlerExplicit(f.energy,sc)); + custom.predict(); + REQUIRE(customOutput.interactions.front().basePairs == predict(f.energy,seed,out).front().basePairs); +} + +namespace { + +class KineticPruningProbe : public PredictorSeedExtensionKineticPruned { +public: + using PredictorSeedExtensionKineticPruned::PredictorSeedExtensionKineticPruned; + bool skips(const Bounds & before, const Bounds & after, bool left, size_t s1, size_t s2) const; +}; + +bool KineticPruningProbe::skips(const Bounds & before, const Bounds & after, + bool left, size_t s1, size_t s2) const +{ + Candidate c; + c.bounds = after; c.left = left; c.s1 = s1; c.s2 = s2; c.macro = true; + return prune(c,before); +} + +} // namespace + +TEST_CASE("Pruning tables bound local moves at all oriented root types", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + RnaSequence t("t","CGGUAUCGGUAU"); + std::string reversedQuery = "GCUGUAGCUGUA"; + std::reverse(reversedQuery.begin(),reversedQuery.end()); + RnaSequence q("q",reversedQuery); + KineticAccessibility a(t), b(q); + ReverseAccessibility reversed(b); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + KineticOutput output(out); + size_t checked = 0; + for (const std::string & parameters : {"Turner04","Turner99"}) { + for (double temperature : {20.,37.}) { + VrnaHandler vrna(temperature,parameters,false,false); + InteractionEnergyVrna energy(a,reversed,vrna,2,3,false,0,false); + const auto sc = seedConstraint(seedEncoding(Chain{{0,0},{1,1}},12)); + KineticPruningProbe probe(energy,output,NULL,new SeedHandlerExplicit(energy,sc)); + for (bool left : {false,true}) for (size_t root = 0; root < 12; ++root) + for (size_t s1 = 0; s1 <= 2; ++s1) for (size_t s2 = 0; s2 <= 3; ++s2) { + if ((left && root < std::max(s1,s2)+2) + || (!left && root+std::max(s1,s2)+2 >= 12)) continue; + const size_t c1 = left ? root-s1-1 : root+s1+1; + const size_t c2 = left ? root-s2-1 : root+s2+1; + const size_t o1 = left ? c1-1 : c1+1, o2 = left ? c2-1 : c2+1; + const E_type loop = left ? energy.getE_interLeft(c1,root,c2,root) + : energy.getE_interLeft(root,c1,root,c2); + const E_type stack = left ? energy.getE_interLeft(o1,c1,o2,c2) + : energy.getE_interLeft(c1,o1,c2,o2); + if (E_isINF(loop) || E_isINF(stack) || loop+stack >= 0) continue; + const Bounds before{root,root,root,root}; + const Bounds after = left ? Bounds{o1,root,o2,root} : Bounds{root,o1,root,o2}; + a.values.clear(); + // A real local move with delta=-1 must never be rejected by + // a valid local suffix bound (no endpoint terms involved). + a.values[{after[0],after[1]}] = -(loop+stack)-1; + REQUIRE_FALSE(probe.skips(before,after,left,s1,s2)); + ++checked; + } + a.values.clear(); + } + } + REQUIRE(checked > 50); } diff --git a/tests/runKineticSeedExtension.sh b/tests/runKineticSeedExtension.sh index 1d98dc67..f33a3ade 100755 --- a/tests/runKineticSeedExtension.sh +++ b/tests/runKineticSeedExtension.sh @@ -6,7 +6,8 @@ tmp=$(mktemp -d) trap 'rm -rf "$tmp"' EXIT common=(--target=GGGGGG --query=CCCCCC --energy=B --acc=N --threads=1 --default-log-file=/dev/null) -kinetic=(--model=X --mode=K '--seedTQ=3||&3||') +mode="${KINETIC_MODE:-K}" +kinetic=(--model=X --mode="$mode" '--seedTQ=3||&3||') csv=(--outMode=C --outCsvCols=hybridDB,E) # Six base pairs have independently known energy -6 in the base-pair model. @@ -36,14 +37,14 @@ cmp "$tmp/expected" "$tmp/boundaries" printf 'hybridDB;E\n2||||&2||||;-4\n' > "$tmp/expected" cmp "$tmp/expected" "$tmp/ranges" # Signed IntaRNA coordinates skip zero: target position three is labeled 1. -"$bin" "${common[@]}" --model=X --mode=K '--seedTQ=1||&7||' "${csv[@]}" \ +"$bin" "${common[@]}" --model=X --mode="$mode" '--seedTQ=1||&7||' "${csv[@]}" \ --tIdxPos0=-2 --qIdxPos0=5 > "$tmp/shifted" printf 'hybridDB;E\n-2||||||&5||||||;-6\n' > "$tmp/expected" cmp "$tmp/expected" "$tmp/shifted" # Every committed prefix is available for suboptimal output, in energy order. "$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" --outNumber=10 > "$tmp/prefixes" -printf 'hybridDB;E\n1||||||&1||||||;-6\n1|||||&2|||||;-5\n1||||&3||||;-4\n2|||&3|||;-3\n3||&3||;-2\n' > "$tmp/expected" +printf 'hybridDB;E\n1||||||&1||||||;-6\n1||||&3||||;-4\n3||&3||;-2\n' > "$tmp/expected" cmp "$tmp/expected" "$tmp/prefixes" # A zero report count still permits energy tracking and never needs traceback. @@ -54,16 +55,17 @@ cmp "$tmp/empty" "$tmp/zero" test -s "$tmp/min-energy" grep -q -- '-6' "$tmp/min-energy" -# Output constraints apply to explicit seeds too. -"$bin" "${common[@]}" --model=X --mode=K '--seedTQ=3|&3|' \ +# Handler-provided explicit seeds can contain lonely pairs. +"$bin" "${common[@]}" --model=X --mode="$mode" '--seedTQ=3|&3|' \ "${csv[@]}" --outNoLP > "$tmp/lonely" -cmp "$tmp/empty" "$tmp/lonely" -"$bin" --target=GG --query=UU --energy=B --acc=N --model=X --mode=K \ +printf 'hybridDB;E\n1|||||&1|||||;-5\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/lonely" +"$bin" --target=GG --query=UU --energy=B --acc=N --model=X --mode="$mode" \ '--seedTQ=1||&1||' "${csv[@]}" --outNoGUend --default-log-file=/dev/null > "$tmp/gu" cmp "$tmp/empty" "$tmp/gu" # Config-file selection must behave just like command-line selection. -printf 'model=X\nmode=K\nkineticScore=C\nseedTQ=3||&3||\n' > "$tmp/parameters" +printf 'model=X\nmode=%s\nkineticScore=C\nseedTQ=3||&3||\n' "$mode" > "$tmp/parameters" "$bin" "${common[@]}" "${csv[@]}" --parameterFile="$tmp/parameters" > "$tmp/configured" cmp "$tmp/default" "$tmp/configured" @@ -75,7 +77,7 @@ thermo=(--target=AGCGACGCA --query=UGCGUCGCU --accW=0 --accL=0 --temperature=25 --threads=1 --default-log-file=/dev/null) for score in A B C; do for dangles in false true; do - "$bin" "${thermo[@]}" --model=X --mode=K '--seedTQ=3|||&5|||' \ + "$bin" "${thermo[@]}" --model=X --mode="$mode" '--seedTQ=3|||&5|||' \ --kineticScore="$score" --energyNoDangles="$dangles" --outNoLP --outNoGUend \ > "$tmp/predicted" test "$(wc -l < "$tmp/predicted")" -gt 1 @@ -94,7 +96,7 @@ expect_error() { grep -q -- "$expected" "$tmp/bad.out" "$tmp/bad.err" } for model in S P B; do - expect_error 'only with --model=X' "${common[@]}" --model="$model" --mode=K + expect_error 'only with --model=X' "${common[@]}" --model="$model" --mode="$mode" done expect_error 'incompatible with --noSeed' "${common[@]}" "${kinetic[@]}" --noSeed # Explicitly supplying the default A must be rejected outside K as well. @@ -114,6 +116,22 @@ expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" --out="spotPr # Evaluation ignores prediction controls, including invalid kinetic scores. "$bin" "${common[@]}" "${csv[@]}" '--rri=1||||||&1||||||' \ - --model=S --mode=K --kineticScore=D > "$tmp/evaluation" + --model=S --mode="$mode" --kineticScore=D > "$tmp/evaluation" cmp "$tmp/default" "$tmp/evaluation" -echo 'Kinetic seed-extension CLI checks passed' +# A missing or false flag is promoted once, with a visible INFO message. +logging=("${common[@]}") +logging[${#logging[@]}-1]="--default-log-file=$tmp/info.log" +for setting in absent false true; do + : > "$tmp/info.log" + flags=() + if [ "$setting" != absent ]; then flags=(--outNoLP="$setting"); fi + "$bin" "${logging[@]}" "${kinetic[@]}" "${csv[@]}" "${flags[@]}" > "$tmp/info.out" 2> "$tmp/info.err" + if [ "$setting" = true ]; then + ! grep -q 'setting --outNoLP=true' "$tmp/info.log" + else + grep -q 'INFO.*setting --outNoLP=true' "$tmp/info.log" + fi +done +# Run the same API/CLI contract against the experimental subclass. +if [ "$mode" = K ]; then KINETIC_MODE=L bash "$0"; fi +echo "Kinetic seed-extension CLI checks passed ($mode)" From bae5007d3fd78d0b3b90f47348349a967bc91867 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Mon, 5 Oct 2026 14:42:28 +0200 Subject: [PATCH 3/6] Validate pruning bounds and record kinetic benchmark results --- ChangeLog | 6 +- doc/Makefile.am | 4 +- doc/benchmark-kinetic.py | 1 + doc/kinetic-benchmark-20261005.json | 146 ++++++++++++++++++ doc/kinetic-seed-extension.md | 59 ++++++- .../PredictorSeedExtensionKineticPruned.cpp | 7 +- .../PredictorSeedExtensionKineticPruned.h | 4 + tests/PredictorSeedExtensionKinetic_test.cpp | 10 +- 8 files changed, 226 insertions(+), 11 deletions(-) create mode 100644 doc/kinetic-benchmark-20261005.json diff --git a/ChangeLog b/ChangeLog index 031bf170..13960b1c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -113,13 +113,17 @@ energies, restricted partition sums, and trackers. + experimental root-pair local-energy suffix bounds from active parameters; prune before complementarity checks assuming monotone ED, omitting endpoint changes; preserve exact full-energy evaluation of surviving moves + + include all nucleotide mismatch codes in bounds, including N within loops * bin/CommandLineParsing, src/IntaRNA/Makefile.am : + expose subclass as --mode=L; set --outNoLP=true with INFO for K/L * tests/PredictorSeedExtensionKinetic_test.cpp, tests/runKineticSeedExtension.sh : * update independent oracle and CLI expectations for atomic double stacks; cover trusted seeds, caching, pruning limitations and both CLI modes + + check root orientations, temperatures and parameter sets for local bounds * README.md, doc/kinetic-seed-extension.md, doc/benchmark-kinetic.py, - doc/Makefile.am : + doc/kinetic-benchmark-20261005.json, doc/Makefile.am : + * record measured candidate reuse gains and pruning overhead with identical + reported results; keep mode L experimental for future real-world comparison * document revised semantics and reproducible performance evaluation in response to https://github.com/BackofenLab/IntaRNA/pull/254 diff --git a/doc/Makefile.am b/doc/Makefile.am index 4a0c1039..abf7704c 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -7,7 +7,9 @@ EXTRA_DIST = \ analysis/out-overlap.md \ analysis/out-overlap/reproduce.py \ conda.txt \ - kinetic-seed-extension.md benchmark-kinetic.py \ + kinetic-seed-extension.md \ + benchmark-kinetic.py \ + kinetic-benchmark-20261005.json \ doxygen.cfg \ latex-deps/adjcalc.sty \ latex-deps/adjustbox.sty \ diff --git a/doc/benchmark-kinetic.py b/doc/benchmark-kinetic.py index d9efda91..8bcb0283 100644 --- a/doc/benchmark-kinetic.py +++ b/doc/benchmark-kinetic.py @@ -68,6 +68,7 @@ def main(): seconds_median=statistics.median(x[0] for x in samples), rss_KiB_max=max(x[1] for x in samples), samples_seconds=[x[0] for x in samples], + reported_rows=max(0, outputs[label].count(b'\n')-1), sha256=hashlib.sha256(outputs[label]).hexdigest(), equals_K=outputs[label] == outputs['cached-K'])) print(json.dumps(results, indent=2)) diff --git a/doc/kinetic-benchmark-20261005.json b/doc/kinetic-benchmark-20261005.json new file mode 100644 index 00000000..e08f07fb --- /dev/null +++ b/doc/kinetic-benchmark-20261005.json @@ -0,0 +1,146 @@ +[ + { + "case": "random-no-ED", + "variant": "uncached-K", + "seconds_median": 0.01733074593357742, + "rss_KiB_max": 15872, + "samples_seconds": [ + 0.023107814020477235, + 0.016839975956827402, + 0.01652474596630782, + 0.023497022921219468, + 0.01733074593357742 + ], + "reported_rows": 10, + "sha256": "216c3d299c10761bdf0d5ff6633f5cea3803c57e624d24c0932d0df3c86d550b", + "equals_K": true + }, + { + "case": "random-no-ED", + "variant": "cached-K", + "seconds_median": 0.016544974059797823, + "rss_KiB_max": 15872, + "samples_seconds": [ + 0.016544974059797823, + 0.01675323396921158, + 0.016520498087629676, + 0.015950259985402226, + 0.0281936579849571 + ], + "reported_rows": 10, + "sha256": "216c3d299c10761bdf0d5ff6633f5cea3803c57e624d24c0932d0df3c86d550b", + "equals_K": true + }, + { + "case": "random-no-ED", + "variant": "pruned-L", + "seconds_median": 0.06901142804417759, + "rss_KiB_max": 16000, + "samples_seconds": [ + 0.0696510859997943, + 0.06871996296104044, + 0.08893389801960438, + 0.06901142804417759, + 0.06856663594953716 + ], + "reported_rows": 10, + "sha256": "216c3d299c10761bdf0d5ff6633f5cea3803c57e624d24c0932d0df3c86d550b", + "equals_K": true + }, + { + "case": "random-folded", + "variant": "uncached-K", + "seconds_median": 0.26241819199640304, + "rss_KiB_max": 20988, + "samples_seconds": [ + 0.2738717310130596, + 0.2925433471100405, + 0.26241819199640304, + 0.26068354200106114, + 0.26009569107554853 + ], + "reported_rows": 10, + "sha256": "0dfccb866e58f5ce60f2195028b1ed62d252a61142e09f57d0b8175563a6ffae", + "equals_K": true + }, + { + "case": "random-folded", + "variant": "cached-K", + "seconds_median": 0.25387142691761255, + "rss_KiB_max": 20968, + "samples_seconds": [ + 0.25387142691761255, + 0.25619977607857436, + 0.2523735810536891, + 0.2511888009030372, + 0.2838866349775344 + ], + "reported_rows": 10, + "sha256": "0dfccb866e58f5ce60f2195028b1ed62d252a61142e09f57d0b8175563a6ffae", + "equals_K": true + }, + { + "case": "random-folded", + "variant": "pruned-L", + "seconds_median": 0.2807472669519484, + "rss_KiB_max": 21180, + "samples_seconds": [ + 0.2961554420180619, + 0.2807472669519484, + 0.3186804640572518, + 0.27823623293079436, + 0.2778694370063022 + ], + "reported_rows": 10, + "sha256": "0dfccb866e58f5ce60f2195028b1ed62d252a61142e09f57d0b8175563a6ffae", + "equals_K": true + }, + { + "case": "stack-rich", + "variant": "uncached-K", + "seconds_median": 0.6640087210107595, + "rss_KiB_max": 22784, + "samples_seconds": [ + 0.6644330089911819, + 0.6638223639456555, + 0.6640087210107595, + 0.6615705629810691, + 0.680608322029002 + ], + "reported_rows": 10, + "sha256": "988a1a76fcae4f2dd6bf40b47f2f79a456e9a70dea3dc1bf7a6447a46023cb9d", + "equals_K": true + }, + { + "case": "stack-rich", + "variant": "cached-K", + "seconds_median": 0.582312888931483, + "rss_KiB_max": 22784, + "samples_seconds": [ + 0.5762485179584473, + 0.6047428960446268, + 0.5761672409716994, + 0.582312888931483, + 0.6001239109318703 + ], + "reported_rows": 10, + "sha256": "988a1a76fcae4f2dd6bf40b47f2f79a456e9a70dea3dc1bf7a6447a46023cb9d", + "equals_K": true + }, + { + "case": "stack-rich", + "variant": "pruned-L", + "seconds_median": 0.6505173300392926, + "rss_KiB_max": 23144, + "samples_seconds": [ + 0.6505173300392926, + 0.6436482360586524, + 0.6449339290848002, + 0.6763418060727417, + 0.6781582959229127 + ], + "reported_rows": 10, + "sha256": "988a1a76fcae4f2dd6bf40b47f2f79a456e9a70dea3dc1bf7a6447a46023cb9d", + "equals_K": true + } +] diff --git a/doc/kinetic-seed-extension.md b/doc/kinetic-seed-extension.md index fe574978..406663fe 100644 --- a/doc/kinetic-seed-extension.md +++ b/doc/kinetic-seed-extension.md @@ -110,7 +110,7 @@ For ViennaRNA, precompute local loop-plus-stack minima for each of the six oriented root base-pair types, both extension sides and each gap pair. Use the **active temperature-scaled parameter set**, including custom parameter files. Minimize over the closing and outer pair types and all four adjacent nucleotide -identities. Relaxing consistency between these identities can only make this +identities, including unknown nucleotide code 0 within loops. Relaxing consistency between these identities can only make this local estimate more optimistic. For the base-pair model, two added pairs have twice its configured base-pair energy. Unknown energy subclasses and API loop limits above the CLI maximum of 30 fall back to unpruned K enumeration. @@ -157,3 +157,60 @@ ED, GU restrictions, retained prefixes, explicit seeds and annotations, cache reuse and repeated calls. L has differential and limitation tests. CLI tests exercise both modes, automatic noLP INFO logging, incompatible requests, and independent reevaluation of predicted structures through `--rri`. + +## Benchmark record (2026-10-05) + +Linux x86-64, AMD Ryzen 5 7530U, GCC 14.4 release (`-O3`), ViennaRNA 2.7.2, +Boost 1.85, Kokkos mdspan. These are synthetic measurements, not a real-world +screening benchmark. Each cell is the median of five single-threaded process +runs after one warm-up. Run order rotates. Timing includes startup, folding, +seed enumeration, pruning-table setup and prediction. The script records each +sample, peak RSS and output hashes in +[the raw results](kinetic-benchmark-20261005.json). + +The comparison binary uses the same revised move rules, scoring and seed +semantics. Its only change is rebuilding **both** end tables after each move. +It still shares complementarity checks within an update. To reproduce it in +a separate build, replace this line in `extendSeed()`: + +```cpp +buildCandidates(sides[best.left ? 0 : 1], bounds, best.left, last1, last2); +``` + +with: + +```cpp +buildCandidates(sides[0], bounds, true, last1, last2); +buildCandidates(sides[1], bounds, false, last1, last2); +``` + +Build both versions with identical release flags, then run: + +```sh +python3 doc/benchmark-kinetic.py --cached /path/to/revised/IntaRNA \ + --uncached /path/to/rebuild-both/IntaRNA --repeat 5 > measurements.json +``` + +The random cases use deterministic Python seed 254, a 600-nt target and 80-nt +query; the folded case uses `accW=150`, `accL=100`. The stack-rich case uses +100 Gs against 30 Cs without accessibility costs. All use the ViennaRNA energy +model, seven-pair seeds, `intLenMax=60`, `intLoopMax=10`, score A and ten reports. + +| Input | Rebuild both ends K (s) | Cached K (s) | Experimental L (s) | +| --- | ---: | ---: | ---: | +| random-no-ED | 0.0173 | 0.0165 | 0.0690 | +| random-folded | 0.2624 | 0.2539 | 0.2807 | +| stack-rich | 0.6640 | 0.5823 | 0.6505 | + +Every variant produced the same ten reported structures and energies for these +inputs. Candidate reuse reduced the stack-rich median by about 12%; the short +random/folded runs showed only small gains. Peak RSS was about 15.5 MiB for +random/no-ED, 20.5–20.7 MiB for folded input and 22.3–22.6 MiB for stack-rich +input, without a meaningful memory improvement. Runtime gains vary with the +input and host load; this does not establish a general speedup over other +IntaRNA predictors. + +L was slower than cached K on all three samples. Precomputation and extra ED +lookups outweighed any saved candidate work. Together with the endpoint and +monotonicity limitations, this supports keeping L as an explicit experimental +subclass for later real-world benchmarking, rather than enabling it by default. diff --git a/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp b/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp index e35f9f53..b55c3db9 100644 --- a/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp +++ b/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp @@ -48,9 +48,10 @@ PredictorSeedExtensionKineticPruned::PredictorSeedExtensionKineticPruned( ? params.stack[outer][reverse[close]] : params.stack[close][reverse[outer]])); } // Relax sequence consistency of the four mismatch bases. - // This can only lower the local loop/stack minimum. - for (int a = 1; a <= 4; ++a) for (int b = 1; b <= 4; ++b) - for (int c = 1; c <= 4; ++c) for (int d = 1; d <= 4; ++d) { + // Include unknown nucleotide code 0, which can occur inside + // loops. This can only lower the local loop/stack minimum. + for (int a = 0; a <= 4; ++a) for (int b = 0; b <= 4; ++b) + for (int c = 0; c <= 4; ++c) for (int d = 0; d <= 4; ++d) { const E_type loop = E_IntLoop(s1, s2, side == 0 ? close : root, reverse[side == 0 ? root : close], a,b,c,d, ¶ms); if (E_isNotINF(loop) && E_isNotINF(stack)) { diff --git a/src/IntaRNA/PredictorSeedExtensionKineticPruned.h b/src/IntaRNA/PredictorSeedExtensionKineticPruned.h index 342fa302..50d48f33 100644 --- a/src/IntaRNA/PredictorSeedExtensionKineticPruned.h +++ b/src/IntaRNA/PredictorSeedExtensionKineticPruned.h @@ -3,6 +3,10 @@ #include "IntaRNA/PredictorSeedExtensionKinetic.h" +#include +#include +#include + namespace IntaRNA { /** diff --git a/tests/PredictorSeedExtensionKinetic_test.cpp b/tests/PredictorSeedExtensionKinetic_test.cpp index 6e2ef0b1..2a7308ce 100644 --- a/tests/PredictorSeedExtensionKinetic_test.cpp +++ b/tests/PredictorSeedExtensionKinetic_test.cpp @@ -597,8 +597,8 @@ bool KineticPruningProbe::skips(const Bounds & before, const Bounds & after, TEST_CASE("Pruning tables bound local moves at all oriented root types", "[PredictorSeedExtensionKinetic]") { #include "testEasyLoggingSetup.icc" - RnaSequence t("t","CGGUAUCGGUAU"); - std::string reversedQuery = "GCUGUAGCUGUA"; + RnaSequence t("t","CGGUAUNCGGUAU"); + std::string reversedQuery = "GCUGUANGCUGUA"; std::reverse(reversedQuery.begin(),reversedQuery.end()); RnaSequence q("q",reversedQuery); KineticAccessibility a(t), b(q); @@ -610,12 +610,12 @@ TEST_CASE("Pruning tables bound local moves at all oriented root types", "[Predi for (double temperature : {20.,37.}) { VrnaHandler vrna(temperature,parameters,false,false); InteractionEnergyVrna energy(a,reversed,vrna,2,3,false,0,false); - const auto sc = seedConstraint(seedEncoding(Chain{{0,0},{1,1}},12)); + const auto sc = seedConstraint(seedEncoding(Chain{{0,0},{1,1}},t.size())); KineticPruningProbe probe(energy,output,NULL,new SeedHandlerExplicit(energy,sc)); - for (bool left : {false,true}) for (size_t root = 0; root < 12; ++root) + for (bool left : {false,true}) for (size_t root = 0; root < t.size(); ++root) for (size_t s1 = 0; s1 <= 2; ++s1) for (size_t s2 = 0; s2 <= 3; ++s2) { if ((left && root < std::max(s1,s2)+2) - || (!left && root+std::max(s1,s2)+2 >= 12)) continue; + || (!left && root+std::max(s1,s2)+2 >= t.size())) continue; const size_t c1 = left ? root-s1-1 : root+s1+1; const size_t c2 = left ? root-s2-1 : root+s2+1; const size_t o1 = left ? c1-1 : c1+1, o2 = left ? c2-1 : c2+1; From 90d0ae55c63fc75bd3245f39ac570b4a97217728 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Mon, 5 Oct 2026 23:01:18 +0200 Subject: [PATCH 4/6] Add IntaRNAkix personality and remove experimental pruning --- ChangeLog | 48 +- README.md | 45 +- configure.ac | 2 +- doc/Makefile.am | 12 +- doc/benchmark-kinetic.py | 78 -- doc/benchmark-kix.py | 149 ++++ doc/kinetic-benchmark-20261005.json | 146 ---- doc/kinetic-seed-extension.md | 193 +++-- doc/kix-benchmark-20261005.json | 727 ++++++++++++++++++ ...taRNAkix.PredictorSeedExtensionKinetic.svg | 86 +++ src/IntaRNA/InteractionEnergyVrna.h | 16 - src/IntaRNA/Makefile.am | 2 - src/IntaRNA/PredictorSeedExtensionKinetic.cpp | 19 - src/IntaRNA/PredictorSeedExtensionKinetic.h | 17 +- .../PredictorSeedExtensionKineticPruned.cpp | 110 --- .../PredictorSeedExtensionKineticPruned.h | 69 -- src/bin/CommandLineParsing.cpp | 30 +- src/bin/CommandLineParsing.h | 2 + tests/PredictorSeedExtensionKinetic_test.cpp | 112 --- tests/runKineticSeedExtension.sh | 48 +- 20 files changed, 1183 insertions(+), 728 deletions(-) delete mode 100644 doc/benchmark-kinetic.py create mode 100644 doc/benchmark-kix.py delete mode 100644 doc/kinetic-benchmark-20261005.json create mode 100644 doc/kix-benchmark-20261005.json create mode 100644 doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg delete mode 100644 src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp delete mode 100644 src/IntaRNA/PredictorSeedExtensionKineticPruned.h diff --git a/ChangeLog b/ChangeLog index 13960b1c..47227e40 100644 --- a/ChangeLog +++ b/ChangeLog @@ -17,8 +17,9 @@ - deterministic kinetic seed extension with --model=X --mode=K and configurable local scoring; complete-energy downhill steps, atomic noLP loop/stack moves and retained greedy traceback; equilibrium probability output is unsupported -- kinetic modes always use stacked extensions and trust handler-provided seeds; - reuse candidate pairing/local energies and add experimental pruning mode L +- kinetic mode always uses stacked extensions and trusts handler-provided seeds; + reuse candidate pairing/local energies +- IntaRNAkix personality enables kinetic seed extension with noLP by default - apply outDeltaE relative to the sequence pair's best interaction when merging regions; preserve local windows for outPerRegion output (PR #253) @@ -41,6 +42,8 @@ ## Technical changes and Optimizations +- install personality links in out-of-tree builds, including IntaRNAkix + - BUGFIX : normalize single-pair suboptimal boundaries before traceback and boundary-only output validation (PR #253) @@ -103,29 +106,37 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261005 Alexander Mitrofanov + * bin/CommandLineParsing, tests/runKineticSeedExtension.sh : + + add IntaRNAkix personality with model=X, mode=K and outNoLP=true defaults; + support executable-name and explicit personality selection + * configure.ac : + * find personality declarations via srcdir so out-of-tree installations + create the executable links, including IntaRNAkix + * README.md, doc/kinetic-seed-extension.md, + doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg, doc/Makefile.am : + + document personality usage and depict seed initialization, allowed moves, + complete-energy greedy selection, stopping and prefix reporting + + distribute the benchmark's tutorial sequences with its script and results + * doc/benchmark-kix.py, doc/kix-benchmark-20261005.json : + + compare IntaRNAkix with default IntaRNA with and without GU-end constraints; + record wall time, peak RSS, MFE energy and maximum covered strand length + * implement https://github.com/BackofenLab/IntaRNA/pull/254#issuecomment-5995041057 + 261005 Alexander Mitrofanov * IntaRNA/PredictorSeedExtensionKinetic : * always evaluate single and double stacks plus atomic loop/stack extensions; trust seed-handler structures/energies, including explicit lonely pairs * cache each end's candidates, shared complementarity checks and local energies; rebuild only the chosen end and refresh full opposite-end energy - * IntaRNA/PredictorSeedExtensionKineticPruned, InteractionEnergyVrna : - + experimental root-pair local-energy suffix bounds from active parameters; - prune before complementarity checks assuming monotone ED, omitting endpoint - changes; preserve exact full-energy evaluation of surviving moves - + include all nucleotide mismatch codes in bounds, including N within loops - * bin/CommandLineParsing, src/IntaRNA/Makefile.am : - + expose subclass as --mode=L; set --outNoLP=true with INFO for K/L + * bin/CommandLineParsing : + * set --outNoLP=true with INFO for kinetic mode K * tests/PredictorSeedExtensionKinetic_test.cpp, tests/runKineticSeedExtension.sh : * update independent oracle and CLI expectations for atomic double stacks; - cover trusted seeds, caching, pruning limitations and both CLI modes - + check root orientations, temperatures and parameter sets for local bounds - * README.md, doc/kinetic-seed-extension.md, doc/benchmark-kinetic.py, - doc/kinetic-benchmark-20261005.json, doc/Makefile.am : - * record measured candidate reuse gains and pruning overhead with identical - reported results; keep mode L experimental for future real-world comparison - * document revised semantics and reproducible performance evaluation in - response to https://github.com/BackofenLab/IntaRNA/pull/254 + cover trusted seeds, caching and kinetic mode K + * README.md, doc/kinetic-seed-extension.md : + * document revised semantics in response to + https://github.com/BackofenLab/IntaRNA/pull/254 261005 Alexander Mitrofanov * IntaRNA/PredictorMfe, tests/PredictorMfeHeuristicCellState_test.cpp, @@ -207,8 +218,7 @@ energies, restricted partition sums, and trackers. + validate local move choices, macro-step energetics, structural constraints, traceback, output filtering and CLI compatibility * README.md, doc/kinetic-seed-extension.md, doc/Makefile.am : - + document the council's implementation consensus, scientific scope and - replacement of unsafe loop-only pruning with exact move enumeration + + document scientific scope and complete-energy move enumeration 261002 Alexander Mitrofanov * doc/analysis/out-overlap.md, doc/analysis/out-overlap/reproduce.py : diff --git a/README.md b/README.md index f022c19b..075a87e2 100644 --- a/README.md +++ b/README.md @@ -100,6 +100,7 @@ The following topics are covered by this documentation: - [IntaRNAsTar - optimized for sRNA-target prediction](#IntaRNAsTar) - [IntaRNAseed - identifys and reports seed interactions only](#IntaRNAseed) - [IntaRNAens - ensemble-based prediction and partition function computation](#IntaRNAens) + - [IntaRNAkix - kinetic seed extension](#IntaRNAkix) - [IntaRNAeval - evaluate predefined interactions](#IntaRNAeval) - [How to constrain predicted interactions](#constraintSetup) - [Interaction restrictions](#interConstr) @@ -755,18 +756,13 @@ traceback preserves the actual chosen path. `--outNoGUend`, separate query and target loop/span limits, regions, output energy/accessibility filters and overlap settings remain applicable. -`--mode=L` provides an experimental subclass with root-pair-specific local -loop-plus-stack tables and early pruning before complementarity checks. It -uses active energy parameters, assumes monotone accessibility costs, and -ignores terminal/dangling changes in its pruning estimate. It can choose -different paths from K and is available for comparative benchmarking. - -Both modes are zippering-inspired heuristics, without a calibrated time axis +The [IntaRNAkix personality](#IntaRNAkix) selects this mode with noLP enabled +by default. It is a zippering-inspired heuristic without a calibrated time axis or a guarantee of the global minimum. Equilibrium probability/partition-sum outputs are rejected, as are other models and `--noSeed`. Scores B and C are optional distance preferences, not measured kinetic rates. See the -[algorithm and benchmark documentation](doc/kinetic-seed-extension.md) for -precise semantics, pruning limitations and validation cases. +[algorithm and preliminary benchmark](doc/kinetic-seed-extension.md). + [![up](doc/figures/icon-up.28.png) back to overview](#overview) @@ -1042,6 +1038,37 @@ IntaRNA --mode=S ... [![up](doc/figures/icon-up.28.png) back to overview](#overview) +### IntaRNAkix + +**IntaRNAkix** (kinetic seed extension) grows every handler-provided seed by +choosing the most favorable complete energy change at either end. It sets +`--model=X --mode=K --outNoLP=true`; other defaults are those of IntaRNA. +Each move adds one stacked pair, two stacked pairs, or a loop-closing pair +plus its outward stack. Two-pair moves are evaluated and committed together. +Seeds themselves may contain lonely pairs when supplied by the seed handler. + +The following calls are equivalent: + +```sh +IntaRNAkix -t target.fasta -q query.fasta +IntaRNA --personality=IntaRNAkix -t target.fasta -q query.fasta +IntaRNA --model=X --mode=K --outNoLP=true -t target.fasta -q query.fasta +``` + +Only strictly downhill moves are accepted. The default score A chooses the +largest energy decrease; `--kineticScore=B|C` adds distance preferences. +The reported MFE is the best visited, reportable interaction across seeds; +this greedy search has no global-optimum or physical folding-time guarantee. + +![IntaRNAkix initialization, allowed extensions, greedy update and stopping rule](doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg) + +See the [algorithm and preliminary benchmark](doc/kinetic-seed-extension.md) +for time, peak-memory, energy and interaction-length comparisons with default +IntaRNA, with and without GU-end restrictions. + +[![up](doc/figures/icon-up.28.png) back to overview](#overview) + + ### IntaRNAeval **IntaRNAeval** evaluates predefined RNA-RNA interactions with the selected diff --git a/configure.ac b/configure.ac index 3b52ea77..ae2df6d9 100644 --- a/configure.ac +++ b/configure.ac @@ -514,7 +514,7 @@ AS_IF([test "$DEPENDENCYNOTFOUND" = "1"], [ ########################################################################## # get available personalities -AC_SUBST([PERSONALITIES],[`grep 'return.*"IntaRNA..*"' ./src/bin/CommandLineParsing.h | grep -o '".*"' | tr -d '"' | tr "\n" " "`]) +AC_SUBST([PERSONALITIES],[`grep 'return.*"IntaRNA..*"' "$srcdir/src/bin/CommandLineParsing.h" | grep -o '".*"' | tr -d '"' | tr "\n" " "`]) ########################################################################## # Keep strict diagnostics focused on IntaRNA sources. Boost and ViennaRNA diff --git a/doc/Makefile.am b/doc/Makefile.am index abf7704c..b8c4b6bf 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -8,8 +8,16 @@ EXTRA_DIST = \ analysis/out-overlap/reproduce.py \ conda.txt \ kinetic-seed-extension.md \ - benchmark-kinetic.py \ - kinetic-benchmark-20261005.json \ + benchmark-kix.py \ + kix-benchmark-20261005.json \ + handson/README.md \ + handson/fhlA.fasta \ + handson/OxyS.fasta \ + handson/phoB.fasta \ + handson/GcvB.fasta \ + handson/ilvE.fasta \ + handson/GcvB.ST.fasta \ + recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg \ doxygen.cfg \ latex-deps/adjcalc.sty \ latex-deps/adjustbox.sty \ diff --git a/doc/benchmark-kinetic.py b/doc/benchmark-kinetic.py deleted file mode 100644 index 8bcb0283..00000000 --- a/doc/benchmark-kinetic.py +++ /dev/null @@ -1,78 +0,0 @@ -#!/usr/bin/env python3 -"""Reproducible end-to-end kinetic benchmark; emits timings and output hashes. - -Use an otherwise identical binary with both ends rebuilt after every move as ---uncached to isolate candidate reuse. All runs are single-threaded. Timings -include startup, accessibility, seed generation, table setup and prediction. -""" -import argparse -import hashlib -import json -from pathlib import Path -import random -import statistics -import subprocess -import sys -import tempfile -import time - - -def main(): - parser = argparse.ArgumentParser(description=__doc__) - parser.add_argument('--cached', required=True, type=Path) - parser.add_argument('--uncached', type=Path) - parser.add_argument('--repeat', type=int, default=5) - args = parser.parse_args() - if args.repeat < 1: - parser.error('--repeat must be positive') - rng = random.Random(254) - target = ''.join(rng.choices('ACGU', k=600)) - query = ''.join(rng.choices('ACGU', k=80)) - cases = [('random-no-ED', target, query, ['--acc=N']), - ('random-folded', target, query, ['--acc=C', '--accW=150', '--accL=100']), - ('stack-rich', 'G'*100, 'C'*30, ['--acc=N'])] - variants = [('cached-K', str(args.cached.resolve()), 'K'), - ('pruned-L', str(args.cached.resolve()), 'L')] - if args.uncached: - variants.insert(0, ('uncached-K', str(args.uncached.resolve()), 'K')) - results = [] - with tempfile.TemporaryDirectory(prefix='intarna-kinetic-bench-') as temp: - for case, t, q, extra in cases: - common = ['--target='+t, '--query='+q, '--model=X', '--seedBP=7', - '--intLenMax=60', '--intLoopMax=10', '--threads=1', - '--outNoLP', '--outNumber=10', '--outMode=C', - '--outCsvCols=hybridDB,E', '--default-log-file=/dev/null'] + extra - measurements = {label: [] for label, _, _ in variants} - outputs = {} - # Warm up every variant, then rotate run order across repetitions. - for iteration in range(args.repeat+1): - order = variants[iteration % len(variants):] + variants[:iteration % len(variants)] - for label, binary, mode in order: - stats = Path(temp)/'time.txt' - command = [binary, '--mode='+mode] + common - start = time.perf_counter() - run = subprocess.run(['/usr/bin/time', '-f', '%M', '-o', str(stats)] + command, - check=True, capture_output=True) - elapsed = time.perf_counter()-start - print(f"{case}/{label} run {iteration}: {elapsed:.3f}s", file=sys.stderr, flush=True) - if label in outputs and outputs[label] != run.stdout: - raise RuntimeError(f'Non-deterministic output: {case}/{label}') - outputs[label] = run.stdout - if iteration: - measurements[label].append((elapsed, int(stats.read_text().strip()))) - if 'uncached-K' in outputs and outputs['uncached-K'] != outputs['cached-K']: - raise RuntimeError(f'Candidate caching changed the result: {case}') - for label, _, _ in variants: - samples = measurements[label] - results.append(dict(case=case, variant=label, - seconds_median=statistics.median(x[0] for x in samples), - rss_KiB_max=max(x[1] for x in samples), - samples_seconds=[x[0] for x in samples], - reported_rows=max(0, outputs[label].count(b'\n')-1), - sha256=hashlib.sha256(outputs[label]).hexdigest(), - equals_K=outputs[label] == outputs['cached-K'])) - print(json.dumps(results, indent=2)) - - -if __name__ == '__main__': - main() diff --git a/doc/benchmark-kix.py b/doc/benchmark-kix.py new file mode 100644 index 00000000..bbe3faa0 --- /dev/null +++ b/doc/benchmark-kix.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python3 +"""Small, single-thread comparison of IntaRNAkix and default IntaRNA. + +Requires Python 3, GNU time and a release IntaRNA binary. Run from any directory: + python3 doc/benchmark-kix.py /path/to/IntaRNA --output results.json +""" + +import argparse +import csv +import hashlib +import io +import json +import os +from pathlib import Path +import platform +import statistics +import subprocess +import sys +import tempfile +import time + + +PAIRS = [("fhlA", "OxyS"), ("phoB", "GcvB"), ("ilvE", "GcvB.ST")] +COLUMNS = "start1,end1,start2,end2,E,hybridDB" + + +def read_fasta(path): + lines = path.read_text().splitlines() + assert sum(line.startswith(">") for line in lines) == 1, path + sequence = "".join(line.strip() for line in lines if not line.startswith(">")) + return {"file": "handson/" + path.name, "header": lines[0][1:], + "length_nt": len(sequence), "sequence": sequence, + "sha256": hashlib.sha256(path.read_bytes()).hexdigest()} + + +def prediction(stdout): + rows = list(csv.DictReader(io.StringIO(stdout), delimiter=";")) + if not rows: + return None + assert len(rows) == 1, rows + row = rows[0] + result = {key: int(row[key]) for key in ("start1", "end1", "start2", "end2")} + result.update(E_kcal_mol=float(row["E"]), hybridDB=row["hybridDB"]) + # All inputs use the default one-based, ascending coordinates. + result["length_nt"] = max(result["end1"] - result["start1"] + 1, + result["end2"] - result["start2"] + 1) + return result + + +def distribution(samples): + return {"median": statistics.median(samples), "min": min(samples), "max": max(samples)} + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("binary", type=Path) + parser.add_argument("--output", required=True, type=Path) + parser.add_argument("--repetitions", type=int, default=5) + parser.add_argument("--warmups", type=int, default=1) + parser.add_argument("--time", type=Path, default=Path("/usr/bin/time")) + parser.add_argument("--build-description", default="unspecified") + args = parser.parse_args() + if args.repetitions < 1 or args.warmups < 0: + parser.error("repetitions must be positive and warmups nonnegative") + binary = args.binary.resolve() + timer = args.time.resolve() + fixtures = Path(__file__).resolve().parent / "handson" + common = ["--threads=1", "--outMode=C", "--outCsvCols=" + COLUMNS, + "--outNumber=1", "--default-log-file=/dev/null"] + env = dict(os.environ, OMP_NUM_THREADS="1") + report = { + "schema": 1, + "binary_sha256": hashlib.sha256(binary.read_bytes()).hexdigest(), + "version": subprocess.check_output([str(binary), "--version"], text=True).strip(), + "build": args.build_description, + "platform": platform.platform(), "machine": platform.machine(), + "python": platform.python_version(), + "timer": subprocess.check_output([str(timer), "--version"], text=True).splitlines()[0], + "common_arguments": common, + "repetitions": args.repetitions, "warmups_per_configuration": args.warmups, + "timing": "perf_counter wall seconds around GNU time + binary; includes startup, folding, seed search and prediction", + "memory": "GNU time %M: maximum resident set size of each child process, KiB", + "noGU": "--outNoGUend=true; seedNoGU remains false for both programs", + "baseline": "default IntaRNA (model X, mode H, outNoLP false)", + "comparison": "IntaRNAkix defaults (model X, mode K, outNoLP true, kineticScore A)", + "length": "max(end1-start1+1, end2-start2+1), nt; default one-based coordinates", + "deviations": "signed IntaRNAkix minus default IntaRNA, within each GU setting; not a global-optimum error bound", + "cases": [], + } + with tempfile.TemporaryDirectory(prefix="intarna-kix-benchmark-") as directory: + rss_file = Path(directory) / "time.txt" + for target, query in PAIRS: + inputs = ["--target=" + str(fixtures / (target + ".fasta")), + "--query=" + str(fixtures / (query + ".fasta"))] + case = {"name": target + "/" + query, + "target": read_fasta(fixtures / (target + ".fasta")), + "query": read_fasta(fixtures / (query + ".fasta")), "runs": [], "comparisons": []} + configurations = [] + for no_gu in (False, True): + for program in ("IntaRNA", "IntaRNAkix"): + options = [] if program == "IntaRNA" else ["--personality=IntaRNAkix"] + options.append("--outNoGUend=" + str(no_gu).lower()) + configurations.append((program, no_gu, options)) + case["runs"].append({"program": program, "noGU": no_gu, + "arguments": options, "samples": []}) + # Warm each configuration; rotate execution order every round. + for iteration in range(args.warmups + args.repetitions): + for offset in range(len(configurations)): + index = (iteration + offset) % len(configurations) + program, no_gu, options = configurations[index] + run = case["runs"][index] + command = [str(binary), *common, *inputs, *options] + start = time.perf_counter() + process = subprocess.run([str(timer), "-f", "%M", "-o", str(rss_file), + *command], env=env, text=True, capture_output=True, check=True) + elapsed = time.perf_counter() - start + result = prediction(process.stdout) + if "prediction" in run: + assert run["prediction"] == result, (case["name"], options, process.stdout) + run["prediction"] = result + if iteration >= args.warmups: + run["samples"].append({"wall_s": elapsed, "peak_rss_kib": int(rss_file.read_text())}) + for run in case["runs"]: + run["wall_s"] = distribution([sample["wall_s"] for sample in run["samples"]]) + run["peak_rss_kib"] = distribution([sample["peak_rss_kib"] for sample in run["samples"]]) + # Outside timing, reevaluate each reported structure independently. + if run["prediction"] is not None: + evaluated = subprocess.check_output([str(binary), *common, *inputs, + "--rri=" + run["prediction"]["hybridDB"]], + env=env, text=True) + assert prediction(evaluated) == run["prediction"], (case["name"], run) + run["structure_reevaluation_matches"] = True + for index, no_gu in ((0, False), (2, True)): + baseline, kix = case["runs"][index:index + 2] + available = baseline["prediction"] is not None and kix["prediction"] is not None + case["comparisons"].append({ + "noGU": no_gu, + "delta_E_kcal_mol": round(kix["prediction"]["E_kcal_mol"] - baseline["prediction"]["E_kcal_mol"], 2) if available else None, + "delta_length_nt": kix["prediction"]["length_nt"] - baseline["prediction"]["length_nt"] if available else None, + "wall_ratio_kix_over_default": kix["wall_s"]["median"] / baseline["wall_s"]["median"], + "rss_ratio_kix_over_default": kix["peak_rss_kib"]["median"] / baseline["peak_rss_kib"]["median"], + }) + report["cases"].append(case) + print(case["name"] + ": " + json.dumps(case["comparisons"]), file=sys.stderr, flush=True) + args.output.write_text(json.dumps(report, indent=2) + "\n") + + +if __name__ == "__main__": + main() diff --git a/doc/kinetic-benchmark-20261005.json b/doc/kinetic-benchmark-20261005.json deleted file mode 100644 index e08f07fb..00000000 --- a/doc/kinetic-benchmark-20261005.json +++ /dev/null @@ -1,146 +0,0 @@ -[ - { - "case": "random-no-ED", - "variant": "uncached-K", - "seconds_median": 0.01733074593357742, - "rss_KiB_max": 15872, - "samples_seconds": [ - 0.023107814020477235, - 0.016839975956827402, - 0.01652474596630782, - 0.023497022921219468, - 0.01733074593357742 - ], - "reported_rows": 10, - "sha256": "216c3d299c10761bdf0d5ff6633f5cea3803c57e624d24c0932d0df3c86d550b", - "equals_K": true - }, - { - "case": "random-no-ED", - "variant": "cached-K", - "seconds_median": 0.016544974059797823, - "rss_KiB_max": 15872, - "samples_seconds": [ - 0.016544974059797823, - 0.01675323396921158, - 0.016520498087629676, - 0.015950259985402226, - 0.0281936579849571 - ], - "reported_rows": 10, - "sha256": "216c3d299c10761bdf0d5ff6633f5cea3803c57e624d24c0932d0df3c86d550b", - "equals_K": true - }, - { - "case": "random-no-ED", - "variant": "pruned-L", - "seconds_median": 0.06901142804417759, - "rss_KiB_max": 16000, - "samples_seconds": [ - 0.0696510859997943, - 0.06871996296104044, - 0.08893389801960438, - 0.06901142804417759, - 0.06856663594953716 - ], - "reported_rows": 10, - "sha256": "216c3d299c10761bdf0d5ff6633f5cea3803c57e624d24c0932d0df3c86d550b", - "equals_K": true - }, - { - "case": "random-folded", - "variant": "uncached-K", - "seconds_median": 0.26241819199640304, - "rss_KiB_max": 20988, - "samples_seconds": [ - 0.2738717310130596, - 0.2925433471100405, - 0.26241819199640304, - 0.26068354200106114, - 0.26009569107554853 - ], - "reported_rows": 10, - "sha256": "0dfccb866e58f5ce60f2195028b1ed62d252a61142e09f57d0b8175563a6ffae", - "equals_K": true - }, - { - "case": "random-folded", - "variant": "cached-K", - "seconds_median": 0.25387142691761255, - "rss_KiB_max": 20968, - "samples_seconds": [ - 0.25387142691761255, - 0.25619977607857436, - 0.2523735810536891, - 0.2511888009030372, - 0.2838866349775344 - ], - "reported_rows": 10, - "sha256": "0dfccb866e58f5ce60f2195028b1ed62d252a61142e09f57d0b8175563a6ffae", - "equals_K": true - }, - { - "case": "random-folded", - "variant": "pruned-L", - "seconds_median": 0.2807472669519484, - "rss_KiB_max": 21180, - "samples_seconds": [ - 0.2961554420180619, - 0.2807472669519484, - 0.3186804640572518, - 0.27823623293079436, - 0.2778694370063022 - ], - "reported_rows": 10, - "sha256": "0dfccb866e58f5ce60f2195028b1ed62d252a61142e09f57d0b8175563a6ffae", - "equals_K": true - }, - { - "case": "stack-rich", - "variant": "uncached-K", - "seconds_median": 0.6640087210107595, - "rss_KiB_max": 22784, - "samples_seconds": [ - 0.6644330089911819, - 0.6638223639456555, - 0.6640087210107595, - 0.6615705629810691, - 0.680608322029002 - ], - "reported_rows": 10, - "sha256": "988a1a76fcae4f2dd6bf40b47f2f79a456e9a70dea3dc1bf7a6447a46023cb9d", - "equals_K": true - }, - { - "case": "stack-rich", - "variant": "cached-K", - "seconds_median": 0.582312888931483, - "rss_KiB_max": 22784, - "samples_seconds": [ - 0.5762485179584473, - 0.6047428960446268, - 0.5761672409716994, - 0.582312888931483, - 0.6001239109318703 - ], - "reported_rows": 10, - "sha256": "988a1a76fcae4f2dd6bf40b47f2f79a456e9a70dea3dc1bf7a6447a46023cb9d", - "equals_K": true - }, - { - "case": "stack-rich", - "variant": "pruned-L", - "seconds_median": 0.6505173300392926, - "rss_KiB_max": 23144, - "samples_seconds": [ - 0.6505173300392926, - 0.6436482360586524, - 0.6449339290848002, - 0.6763418060727417, - 0.6781582959229127 - ], - "reported_rows": 10, - "sha256": "988a1a76fcae4f2dd6bf40b47f2f79a456e9a70dea3dc1bf7a6447a46023cb9d", - "equals_K": true - } -] diff --git a/doc/kinetic-seed-extension.md b/doc/kinetic-seed-extension.md index 406663fe..c879a7d4 100644 --- a/doc/kinetic-seed-extension.md +++ b/doc/kinetic-seed-extension.md @@ -9,16 +9,14 @@ it has no calibrated transition rates or time axis, does not cross barriers between committed states, and does not guarantee a global minimum. A favorable two-pair move does not establish a barrier-free physical reaction pathway. -`--mode=L` selects the experimental subclass -`PredictorSeedExtensionKineticPruned`. It applies local-energy and accessibility -pruning **before** checking complementarity and evaluating full energies. -Its additional assumptions can change the selected path; K remains the -reference for measuring those changes. Both modes accept scores A/B/C and -support ordinary energy trackers. They reject seedless operation, other models, -and requests for equilibrium partition functions or probabilities. +The `IntaRNAkix` personality (kinetic seed extension) selects +`--model=X --mode=K --outNoLP=true`. It can be invoked through the installed +`IntaRNAkix` executable link or `IntaRNA --personality=IntaRNAkix`. Other defaults +remain those of IntaRNA. Ordinary energy trackers are supported; seedless +operation, other interaction models and equilibrium partition/probability +requests are rejected for mode K. -These semantics incorporate the October 5 review of -[PR #254](https://github.com/BackofenLab/IntaRNA/pull/254#issuecomment-5991039020). +![Kinetic seed extension recursion](recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg) ## Seeds, states and allowed extensions @@ -35,8 +33,9 @@ range offsets and conversion to original coordinates. Initially, `H = seedHandler.getSeedE(i1,i2) + energy.getE_init()`. Extensions **always** use the no-lonely-pair strategy, independent of the API -output constraint. The CLI sets `--outNoLP=true` when absent or false and emits -an INFO message using the normal logging destination. This applies to new +output constraint. Direct `--mode=K` calls promote a missing or false +`--outNoLP` to true with an INFO message using the normal logging destination; +IntaRNAkix already defaults to true. This applies to new extensions, not to revalidation of handler-provided seeds. Allowed moves are: - One pair stacked directly onto the current boundary (`|SEED`). @@ -58,7 +57,7 @@ available if a trajectory stops at an unreportable endpoint. ## Complete energies and deterministic scores -Every candidate surviving geometric and optional pruning checks is evaluated +Every geometrically feasible, complementary candidate is evaluated with the active energy model: ``` @@ -99,47 +98,9 @@ candidate as it evaluates full energies, without a separate selection pass. After committing a move, only that end's tables are rebuilt. The opposite end's pair checks, feasibility and local energies remain valid. Its total -energy, span eligibility and optional pruning decision are refreshed. An +energy and span eligibility are refreshed. An initially uphill candidate is retained because it may become downhill when -the opposite end changes. Experimentally pruned entries remain unresolved and -can be reconsidered in a later state. Tables use O((m1+2)(m2+2)) space per end. - -## Experimental pruning in L - -For ViennaRNA, precompute local loop-plus-stack minima for each of the six -oriented root base-pair types, both extension sides and each gap pair. Use the -**active temperature-scaled parameter set**, including custom parameter files. -Minimize over the closing and outer pair types and all four adjacent nucleotide -identities, including unknown nucleotide code 0 within loops. Relaxing consistency between these identities can only make this -local estimate more optimistic. For the base-pair model, two added pairs have -twice its configured base-pair energy. Unknown energy subclasses and API loop -limits above the CLI maximum of 30 fall back to unpruned K enumeration. - -Apply componentwise suffix minima to the gap tables. Each cell then bounds -the local loop-plus-stack term for that gap and every larger gap pair. Before -checking candidate pairs, compare: - -``` -local_suffix_bound + ED1_next + ED2_next - ED1_current - ED2_current >= 0 -``` - -If true, skip that candidate and all componentwise larger gaps. Within the -ordered rectangular traversal this rejects suffixes without more ED, pairing -or loop-energy lookups. A single stacked pair is always evaluated exactly. -Surviving moves still require strictly downhill **complete** energy changes. - -This is deliberately an experiment, not a certified bound on the complete -energy change. It assumes ED is monotone as intervals grow and omits changes in -terminal and dangling contributions. Imported/nonmonotone accessibility and -favorable endpoint changes can invalidate the filter; tests include a concrete -case where L stops while K continues. The local table includes the mandatory -stack, so it avoids the original bare-loop error: a Turner2004 loop of +0.50 -kcal/mol can be rescued by a -3.30 kcal/mol stack. - -Tables require O(12(m1+1)(m2+1)) space and parameter enumeration at construction. -Setup and extra ED lookups may outweigh pruning benefits on small or short-path -inputs. Thus L remains a separate subclass/mode for later real-world evaluation. -See [the reproducible benchmark](benchmark-kinetic.py) and measurements below. +the opposite end changes. Tables use O((m1+2)(m2+2)) space per end. ## Reporting and validation @@ -148,69 +109,87 @@ boundaries, keep the lowest full energy, then lexicographically smallest chain. Reduce duplicates before the normal optimum collector. Overlap-constrained output can select shorter retained prefixes; traceback restores the selected path directly. Memory for retained paths is proportional to their total length, -not constant per seed. Repeated predictions reset trajectory and ED caches. +not constant per seed. Repeated predictions reset trajectory and candidate caches. The tests compare K with an independent absolute-endpoint oracle that rebuilds and reevaluates whole chains. They cover scores/ties, single and double stacks, positive-loop rescue, strict stopping, separate spans and regions, nonmonotone ED, GU restrictions, retained prefixes, explicit seeds and annotations, cache -reuse and repeated calls. L has differential and limitation tests. CLI tests -exercise both modes, automatic noLP INFO logging, incompatible requests, and +reuse and repeated calls. CLI tests exercise the personality name and option, +explicit parameter overrides, automatic noLP INFO logging, incompatible requests, and independent reevaluation of predicted structures through `--rri`. -## Benchmark record (2026-10-05) - -Linux x86-64, AMD Ryzen 5 7530U, GCC 14.4 release (`-O3`), ViennaRNA 2.7.2, -Boost 1.85, Kokkos mdspan. These are synthetic measurements, not a real-world -screening benchmark. Each cell is the median of five single-threaded process -runs after one warm-up. Run order rotates. Timing includes startup, folding, -seed enumeration, pruning-table setup and prediction. The script records each -sample, peak RSS and output hashes in -[the raw results](kinetic-benchmark-20261005.json). - -The comparison binary uses the same revised move rules, scoring and seed -semantics. Its only change is rebuilding **both** end tables after each move. -It still shares complementarity checks within an update. To reproduce it in -a separate build, replace this line in `extendSeed()`: - -```cpp -buildCandidates(sides[best.left ? 0 : 1], bounds, best.left, last1, last2); -``` - -with: - -```cpp -buildCandidates(sides[0], bounds, true, last1, last2); -buildCandidates(sides[1], bounds, false, last1, last2); -``` - -Build both versions with identical release flags, then run: +## Preliminary benchmark against default IntaRNA + +The small panel uses the repository's tutorial sequences: fhlA/OxyS +(112/108 nt), phoB/GcvB (299/201 nt), and ilvE/GcvB.ST (299/200 nt). +These are the pairs in [hands-on examples 3.2, 3.4 and 3.5](handson/README.md). +No experimental seed, region or accessibility constraints from those examples +are applied here. The raw record includes the sequences and input hashes. + +The comparison uses the actual personality defaults: IntaRNA has model X, +mode H and `outNoLP=false`; IntaRNAkix has model X, mode K, score A and +`outNoLP=true`. Thus energy and length deviations reflect both the search and +the different noLP defaults. Default IntaRNA is itself a heuristic, so the +energy deviation is not a certified error from a global optimum. + +The review's “noGU” setting is interpreted as `--outNoGUend=true`, the same flag +for both programs. This prohibits GU at reported interaction ends and at +nonstacking loop ends; it does **not** forbid all internal GU pairs. +`--seedNoGU` stays at its default false. The other setting explicitly uses +`--outNoGUend=false`. + +Each configuration has one warm-up followed by five measured runs. Execution +order rotates across the four configurations on each pair. All runs use one +thread and compute accessibility from the input sequences. Wall time includes +process startup, accessibility, seed search and prediction; GNU time supplies +the child process's peak resident memory in KiB. The table reports medians; +all samples and min/max values are in the raw JSON. Reported structures and +energies are deterministic across repetitions, and every result is independently +reevaluated via `--rri` outside the timed runs. + +For each reported MFE interaction, covered length is +`L = max(end1-start1+1, end2-start2+1)` in the default one-based coordinates. +Signed deviations are `E_kix - E_default` in kcal/mol and `L_kix - L_default` +in nucleotides, compared within the same GU setting. Missing predictions are +recorded as null, never as zero energy or zero length. + +Measured on 2026-10-05 on Linux x86-64, AMD Ryzen 5 7530U, using GCC 14.4.0 +release (`-O3`, C++23), ViennaRNA 2.7.2, Boost 1.85 and Kokkos mdspan. +No builds or other validation jobs ran alongside these measurements. + +| Pair | noGU | Default time (s) | Kix time (s) | Default peak RSS (KiB) | Kix peak RSS (KiB) | +| --- | --- | ---: | ---: | ---: | ---: | +| fhlA/OxyS | off | 0.0569 | 0.0381 | 16,912 | 16,648 | +| fhlA/OxyS | on | 0.0559 | 0.0381 | 16,660 | 16,796 | +| phoB/GcvB | off | 1.0598 | 0.1491 | 18,448 | 18,456 | +| phoB/GcvB | on | 0.7158 | 0.1446 | 18,456 | 18,452 | +| ilvE/GcvB.ST | off | 1.2935 | 0.1471 | 18,440 | 18,512 | +| ilvE/GcvB.ST | on | 1.0233 | 0.1436 | 18,560 | 18,432 | + +| Pair | noGU | Default E | Kix E | ΔE (kcal/mol) | Default L | Kix L | ΔL (nt) | +| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: | +| fhlA/OxyS | off | -5.59 | -5.57 | +0.02 | 24 | 7 | -17 | +| fhlA/OxyS | on | -5.57 | -5.57 | +0.00 | 7 | 7 | +0 | +| phoB/GcvB | off | -15.70 | -13.19 | +2.51 | 47 | 8 | -39 | +| phoB/GcvB | on | -13.19 | -13.19 | +0.00 | 8 | 8 | +0 | +| ilvE/GcvB.ST | off | -14.24 | -9.84 | +4.40 | 55 | 13 | -42 | +| ilvE/GcvB.ST | on | -10.55 | -9.13 | +1.42 | 40 | 10 | -30 | + +On this small panel, Kix uses 0.11–0.68 times the default runtime. Median peak +RSS differs by less than 2%, within the run-to-run variation. Energy deviations +range from 0 to +4.40 kcal/mol, and length deviations from −42 to 0 nt. +With noGU on, both programs report identical interactions for fhlA/OxyS and +phoB/GcvB. These three selected tutorial pairs are a preliminary performance +and output comparison, not a general speedup or biological-accuracy estimate. + +Reproduce from the repository root with a release binary: ```sh -python3 doc/benchmark-kinetic.py --cached /path/to/revised/IntaRNA \ - --uncached /path/to/rebuild-both/IntaRNA --repeat 5 > measurements.json +python3 doc/benchmark-kix.py /path/to/release/src/bin/IntaRNA \ + --repetitions=5 --warmups=1 --output=kix-benchmark.json ``` -The random cases use deterministic Python seed 254, a 600-nt target and 80-nt -query; the folded case uses `accW=150`, `accL=100`. The stack-rich case uses -100 Gs against 30 Cs without accessibility costs. All use the ViennaRNA energy -model, seven-pair seeds, `intLenMax=60`, `intLoopMax=10`, score A and ten reports. - -| Input | Rebuild both ends K (s) | Cached K (s) | Experimental L (s) | -| --- | ---: | ---: | ---: | -| random-no-ED | 0.0173 | 0.0165 | 0.0690 | -| random-folded | 0.2624 | 0.2539 | 0.2807 | -| stack-rich | 0.6640 | 0.5823 | 0.6505 | - -Every variant produced the same ten reported structures and energies for these -inputs. Candidate reuse reduced the stack-rich median by about 12%; the short -random/folded runs showed only small gains. Peak RSS was about 15.5 MiB for -random/no-ED, 20.5–20.7 MiB for folded input and 22.3–22.6 MiB for stack-rich -input, without a meaningful memory improvement. Runtime gains vary with the -input and host load; this does not establish a general speedup over other -IntaRNA predictors. - -L was slower than cached K on all three samples. Precomputation and extra ED -lookups outweighed any saved candidate work. Together with the endpoint and -monotonicity limitations, this supports keeping L as an explicit experimental -subclass for later real-world benchmarking, rather than enabling it by default. +The [script](benchmark-kix.py) uses Python 3 and GNU time. The +[raw results](kix-benchmark-20261005.json) record all samples, predictions, +signed deviations, flags, sequence data, software versions and binary hash. diff --git a/doc/kix-benchmark-20261005.json b/doc/kix-benchmark-20261005.json new file mode 100644 index 00000000..003a6a24 --- /dev/null +++ b/doc/kix-benchmark-20261005.json @@ -0,0 +1,727 @@ +{ + "schema": 1, + "binary_sha256": "3779a676dd6a96e6edf9fe9a7316c2125466d0c33be51723661eb1513dfcdd64", + "version": "IntaRNA 3.4.1\n using Vienna RNA package 2.7.2 and boost 1.85.0", + "build": "GCC 14.4.0 release (-O3), C++23, ViennaRNA 2.7.2, Boost 1.85, Kokkos mdspan; AMD Ryzen 5 7530U; 2026-10-05", + "platform": "Linux-6.8.0-142-generic-x86_64-with-glibc2.39", + "machine": "x86_64", + "python": "3.12.7", + "timer": "time (GNU Time) UNKNOWN", + "common_arguments": [ + "--threads=1", + "--outMode=C", + "--outCsvCols=start1,end1,start2,end2,E,hybridDB", + "--outNumber=1", + "--default-log-file=/dev/null" + ], + "repetitions": 5, + "warmups_per_configuration": 1, + "timing": "perf_counter wall seconds around GNU time + binary; includes startup, folding, seed search and prediction", + "memory": "GNU time %M: maximum resident set size of each child process, KiB", + "noGU": "--outNoGUend=true; seedNoGU remains false for both programs", + "baseline": "default IntaRNA (model X, mode H, outNoLP false)", + "comparison": "IntaRNAkix defaults (model X, mode K, outNoLP true, kineticScore A)", + "length": "max(end1-start1+1, end2-start2+1), nt; default one-based coordinates", + "deviations": "signed IntaRNAkix minus default IntaRNA, within each GU setting; not a global-optimum error bound", + "cases": [ + { + "name": "fhlA/OxyS", + "target": { + "file": "handson/fhlA.fasta", + "header": "fhlA|NC000913|-53..+60|doi:10.1006/jmbi.2000.3942:Fig7", + "length_nt": 112, + "sequence": "AGUUAGUCAAUGACCUUUUGCACCGCUUUGCGGUGCUUUCCUGGAACAACAAAAUGUCAUAUACACCGAUGAGUGAUCUCGGACAACAAGGGUUGUUCGACAUCACUCGGAC", + "sha256": "667e0df26bf5b8635df5fb5b0eca299cb0484d1c07d0bfd93fff674c2cab9c31" + }, + "query": { + "file": "handson/OxyS.fasta", + "header": "OxyS|NC_000913|56..164|doi:10.1006/jmbi.2000.3942:Fig7", + "length_nt": 108, + "sequence": "GAAACGGAGCGGCACCUCUUUUAACCCUUGAAGUCACUGCCCGUUUCGAGAGUUUCUCAACUCGAAUAACUAAAGCCAACGUGAACUUUUGCGGAUCUCCAGGAUCCG", + "sha256": "bd336bf7f10104837bd7f5e56315dd0f03ba2eb9784d88280588b464adf1a30d" + }, + "runs": [ + { + "program": "IntaRNA", + "noGU": false, + "arguments": [ + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 0.054906384088099, + "peak_rss_kib": 16448 + }, + { + "wall_s": 0.055035896017216146, + "peak_rss_kib": 16912 + }, + { + "wall_s": 0.056886381935328245, + "peak_rss_kib": 16912 + }, + { + "wall_s": 0.057665151078253984, + "peak_rss_kib": 17088 + }, + { + "wall_s": 0.05727210489567369, + "peak_rss_kib": 16788 + } + ], + "prediction": { + "start1": 39, + "end1": 60, + "start2": 81, + "end2": 104, + "E_kcal_mol": -5.59, + "hybridDB": "39|||||||.|..|||||..||||&81||||..|||||..|...|||||||", + "length_nt": 24 + }, + "wall_s": { + "median": 0.056886381935328245, + "min": 0.054906384088099, + "max": 0.057665151078253984 + }, + "peak_rss_kib": { + "median": 16912, + "min": 16448, + "max": 17088 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAkix", + "noGU": false, + "arguments": [ + "--personality=IntaRNAkix", + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 0.038167160004377365, + "peak_rss_kib": 16648 + }, + { + "wall_s": 0.03802963194902986, + "peak_rss_kib": 16576 + }, + { + "wall_s": 0.03889605507720262, + "peak_rss_kib": 16652 + }, + { + "wall_s": 0.03797265503089875, + "peak_rss_kib": 16784 + }, + { + "wall_s": 0.03807136998511851, + "peak_rss_kib": 16620 + } + ], + "prediction": { + "start1": 39, + "end1": 45, + "start2": 98, + "end2": 104, + "E_kcal_mol": -5.57, + "hybridDB": "39|||||||&98|||||||", + "length_nt": 7 + }, + "wall_s": { + "median": 0.03807136998511851, + "min": 0.03797265503089875, + "max": 0.03889605507720262 + }, + "peak_rss_kib": { + "median": 16648, + "min": 16576, + "max": 16784 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNA", + "noGU": true, + "arguments": [ + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.05764365091454238, + "peak_rss_kib": 16660 + }, + { + "wall_s": 0.055861442000605166, + "peak_rss_kib": 16660 + }, + { + "wall_s": 0.0561534590087831, + "peak_rss_kib": 16968 + }, + { + "wall_s": 0.05309296608902514, + "peak_rss_kib": 16640 + }, + { + "wall_s": 0.055408544023521245, + "peak_rss_kib": 16836 + } + ], + "prediction": { + "start1": 39, + "end1": 45, + "start2": 98, + "end2": 104, + "E_kcal_mol": -5.57, + "hybridDB": "39|||||||&98|||||||", + "length_nt": 7 + }, + "wall_s": { + "median": 0.055861442000605166, + "min": 0.05309296608902514, + "max": 0.05764365091454238 + }, + "peak_rss_kib": { + "median": 16660, + "min": 16640, + "max": 16968 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAkix", + "noGU": true, + "arguments": [ + "--personality=IntaRNAkix", + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.038256115978583694, + "peak_rss_kib": 16912 + }, + { + "wall_s": 0.03771536401472986, + "peak_rss_kib": 16772 + }, + { + "wall_s": 0.0362653280608356, + "peak_rss_kib": 16960 + }, + { + "wall_s": 0.03812957904301584, + "peak_rss_kib": 16784 + }, + { + "wall_s": 0.038168090977706015, + "peak_rss_kib": 16796 + } + ], + "prediction": { + "start1": 39, + "end1": 45, + "start2": 98, + "end2": 104, + "E_kcal_mol": -5.57, + "hybridDB": "39|||||||&98|||||||", + "length_nt": 7 + }, + "wall_s": { + "median": 0.03812957904301584, + "min": 0.0362653280608356, + "max": 0.038256115978583694 + }, + "peak_rss_kib": { + "median": 16796, + "min": 16772, + "max": 16960 + }, + "structure_reevaluation_matches": true + } + ], + "comparisons": [ + { + "noGU": false, + "delta_E_kcal_mol": 0.02, + "delta_length_nt": -17, + "wall_ratio_kix_over_default": 0.6692527928459269, + "rss_ratio_kix_over_default": 0.9843897824030274 + }, + { + "noGU": true, + "delta_E_kcal_mol": 0.0, + "delta_length_nt": 0, + "wall_ratio_kix_over_default": 0.6825741992589945, + "rss_ratio_kix_over_default": 1.0081632653061225 + } + ] + }, + { + "name": "phoB/GcvB", + "target": { + "file": "handson/phoB.fasta", + "header": "phoB|NC_000913|b1130|-200..+100|genom-subsequence", + "length_nt": 299, + "sequence": "GAGCTATCACGATGGTTGATGAGCTGAAATAAACCTCGTATCAGTGCCGGATGGCGATGCTGTCCGGCCTGCTTATTAAGATTATCCGCTTTTTATTTTTTCACTTTACCTCCCCTCCCCGCTGGTTTATTTAATGTTTACCCCCATAACCACATAATCGCGTTACACTATTTTAATAATTAAGACAGGGAGAAATAAAAATGCGCGTACTGGTTGTTGAAGACAATGCGTTGTTACGTCACCACCTTAAAGTTCAGATTCAGGATGCTGGTCATCAGGTCGATGACGCAGAAGATGCC", + "sha256": "ce42dc25417d696c2402cc1daaf22ee5045e68c095fcd578b4072395a09ffa98" + }, + "query": { + "file": "handson/GcvB.fasta", + "header": "GcvB|NC_000913", + "length_nt": 201, + "sequence": "ACUUCCUGAGCCGGAACGAAAAGUUUUAUCGGAAUGCGUGUUCUGGUGAACUUUUGGCUUACGGUUGUGAUGUUGUGUUGUUGUGUUUGCAAUUGGUCUGCGAUUCAGACCAUGGUAGCAAAGCUACCUUUUUUCACUUCCUGUACAUUUACCCUGUCUGUCCAUAGUGAUUAAUGUAGCACCGCCUAAUUGCGGUGCUUU", + "sha256": "580180ec3b79266f0684c4c60d44b7b86ef6dbba087f67110d7c25c9ef31bcf7" + }, + "runs": [ + { + "program": "IntaRNA", + "noGU": false, + "arguments": [ + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 1.0826977379620075, + "peak_rss_kib": 18528 + }, + { + "wall_s": 1.0389898010762408, + "peak_rss_kib": 18660 + }, + { + "wall_s": 1.0598379949806258, + "peak_rss_kib": 18448 + }, + { + "wall_s": 1.0448508230037987, + "peak_rss_kib": 18436 + }, + { + "wall_s": 1.0702672540210187, + "peak_rss_kib": 18336 + } + ], + "prediction": { + "start1": 222, + "end1": 268, + "start2": 38, + "end2": 82, + "E_kcal_mol": -15.7, + "hybridDB": "222|||||||||.....|||||||.|||..||||||||...|||||||||&38|||||||||||||||||........|||.||||||||||||||||", + "length_nt": 47 + }, + "wall_s": { + "median": 1.0598379949806258, + "min": 1.0389898010762408, + "max": 1.0826977379620075 + }, + "peak_rss_kib": { + "median": 18448, + "min": 18336, + "max": 18660 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAkix", + "noGU": false, + "arguments": [ + "--personality=IntaRNAkix", + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 0.14911555196158588, + "peak_rss_kib": 18632 + }, + { + "wall_s": 0.1491208989173174, + "peak_rss_kib": 18456 + }, + { + "wall_s": 0.14401900197844952, + "peak_rss_kib": 18652 + }, + { + "wall_s": 0.15017399890348315, + "peak_rss_kib": 18356 + }, + { + "wall_s": 0.147206068970263, + "peak_rss_kib": 18068 + } + ], + "prediction": { + "start1": 183, + "end1": 190, + "start2": 152, + "end2": 159, + "E_kcal_mol": -13.19, + "hybridDB": "183||||||||&152||||||||", + "length_nt": 8 + }, + "wall_s": { + "median": 0.14911555196158588, + "min": 0.14401900197844952, + "max": 0.15017399890348315 + }, + "peak_rss_kib": { + "median": 18456, + "min": 18068, + "max": 18652 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNA", + "noGU": true, + "arguments": [ + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.7203351480420679, + "peak_rss_kib": 18328 + }, + { + "wall_s": 0.6980774619150907, + "peak_rss_kib": 18456 + }, + { + "wall_s": 0.7028307840228081, + "peak_rss_kib": 18396 + }, + { + "wall_s": 0.7158402650384232, + "peak_rss_kib": 18456 + }, + { + "wall_s": 0.7222581719979644, + "peak_rss_kib": 18508 + } + ], + "prediction": { + "start1": 183, + "end1": 190, + "start2": 152, + "end2": 159, + "E_kcal_mol": -13.19, + "hybridDB": "183||||||||&152||||||||", + "length_nt": 8 + }, + "wall_s": { + "median": 0.7158402650384232, + "min": 0.6980774619150907, + "max": 0.7222581719979644 + }, + "peak_rss_kib": { + "median": 18456, + "min": 18328, + "max": 18508 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAkix", + "noGU": true, + "arguments": [ + "--personality=IntaRNAkix", + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.14456413604784757, + "peak_rss_kib": 18192 + }, + { + "wall_s": 0.14236011693719774, + "peak_rss_kib": 18364 + }, + { + "wall_s": 0.1507548289373517, + "peak_rss_kib": 18660 + }, + { + "wall_s": 0.14155757101252675, + "peak_rss_kib": 18452 + }, + { + "wall_s": 0.15410517202690244, + "peak_rss_kib": 18656 + } + ], + "prediction": { + "start1": 183, + "end1": 190, + "start2": 152, + "end2": 159, + "E_kcal_mol": -13.19, + "hybridDB": "183||||||||&152||||||||", + "length_nt": 8 + }, + "wall_s": { + "median": 0.14456413604784757, + "min": 0.14155757101252675, + "max": 0.15410517202690244 + }, + "peak_rss_kib": { + "median": 18452, + "min": 18192, + "max": 18660 + }, + "structure_reevaluation_matches": true + } + ], + "comparisons": [ + { + "noGU": false, + "delta_E_kcal_mol": 2.51, + "delta_length_nt": -39, + "wall_ratio_kix_over_default": 0.1406965523672434, + "rss_ratio_kix_over_default": 1.0004336513443193 + }, + { + "noGU": true, + "delta_E_kcal_mol": 0.0, + "delta_length_nt": 0, + "wall_ratio_kix_over_default": 0.20195027174126337, + "rss_ratio_kix_over_default": 0.9997832683138275 + } + ] + }, + { + "name": "ilvE/GcvB.ST", + "target": { + "file": "handson/ilvE.fasta", + "header": "ilvE|NC_003197|STM3903|-200..+100|genom-subsequence", + "length_nt": 299, + "sequence": "GGTTTTCAGGTGTGCTCCATGAATATGGAAGCCGCGACCGATGCGCAGAATATAAATATTGAATTGACCGTTGCCAGTCCCCGGTCGGTCGACTTACTGTTTAGTCAGTTAAGTAAACTGGTAGATGTTGCGCATGTCGCGATCTGCCAGAGCGCTGCCACATCACAACAAATCCGCGCCTGAGCGCAAAAGGAAGAAAAATGACGACGAAAAAAGCTGATTATATTTGGTTCAATGGCGAGATGGTGCGCTGGGAAGACGCGAAGGTTCACGTAATGTCTCACGCGCTGCACTACGGT", + "sha256": "6fd5cab8545dd5cede4b8d2d23a276f7346a56706b3121de3bcf954d60a2068e" + }, + "query": { + "file": "handson/GcvB.ST.fasta", + "header": "GcvB|NC_003197", + "length_nt": 200, + "sequence": "ACUUCCUGAGCCGGAACGAAAAGUUUUAUCGGAAUGCGUGUUCUGAUGGGCUUUUGGCUUACGGUUGUGAUGUUGUGUUGUUGUGUUUGCAAUUGGUCUGCGAUUCAGACCACGGUAGCGAGACUACCCUUUUUCACUUCCUGUACAUUUACCCUGUCUGUCCAUAGUGAUUAAUGUAGCACCGCCAUAUUGCGGUGCUU", + "sha256": "529466282586ededdea0ed04199197f04978e2be14245feed0655bc6850206bb" + }, + "runs": [ + { + "program": "IntaRNA", + "noGU": false, + "arguments": [ + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 1.3832573930267245, + "peak_rss_kib": 18452 + }, + { + "wall_s": 1.286010636948049, + "peak_rss_kib": 18452 + }, + { + "wall_s": 1.3232459670398384, + "peak_rss_kib": 18064 + }, + { + "wall_s": 1.2741687439847738, + "peak_rss_kib": 18440 + }, + { + "wall_s": 1.2935187839902937, + "peak_rss_kib": 18336 + } + ], + "prediction": { + "start1": 33, + "end1": 81, + "start2": 48, + "end2": 102, + "E_kcal_mol": -14.24, + "hybridDB": "33||||||||||..||||||||..|||||...|||||||||.|||||.|||&48|||...|||||..||||||..|||...|||||..|||.|||||.|||||||.|||", + "length_nt": 55 + }, + "wall_s": { + "median": 1.2935187839902937, + "min": 1.2741687439847738, + "max": 1.3832573930267245 + }, + "peak_rss_kib": { + "median": 18440, + "min": 18064, + "max": 18452 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAkix", + "noGU": false, + "arguments": [ + "--personality=IntaRNAkix", + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 0.14518713497091085, + "peak_rss_kib": 18452 + }, + { + "wall_s": 0.1470703890081495, + "peak_rss_kib": 18512 + }, + { + "wall_s": 0.14972956804558635, + "peak_rss_kib": 18512 + }, + { + "wall_s": 0.1487087750574574, + "peak_rss_kib": 18772 + }, + { + "wall_s": 0.1428277890663594, + "peak_rss_kib": 18452 + } + ], + "prediction": { + "start1": 157, + "end1": 169, + "start2": 64, + "end2": 76, + "E_kcal_mol": -9.84, + "hybridDB": "157||.||||||||||&64||||||||||.||", + "length_nt": 13 + }, + "wall_s": { + "median": 0.1470703890081495, + "min": 0.1428277890663594, + "max": 0.14972956804558635 + }, + "peak_rss_kib": { + "median": 18512, + "min": 18452, + "max": 18772 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNA", + "noGU": true, + "arguments": [ + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 1.0233468420337886, + "peak_rss_kib": 18588 + }, + { + "wall_s": 1.0323557329829782, + "peak_rss_kib": 18448 + }, + { + "wall_s": 1.0218827379867435, + "peak_rss_kib": 18648 + }, + { + "wall_s": 1.0279411339433864, + "peak_rss_kib": 18560 + }, + { + "wall_s": 0.9871100430609658, + "peak_rss_kib": 18520 + } + ], + "prediction": { + "start1": 131, + "end1": 170, + "start2": 39, + "end2": 77, + "E_kcal_mol": -10.55, + "hybridDB": "131||||||||||||||.||||||.....|||.|||||.||||&39||||.||||||||.||||||....|||||||||||.|||", + "length_nt": 40 + }, + "wall_s": { + "median": 1.0233468420337886, + "min": 0.9871100430609658, + "max": 1.0323557329829782 + }, + "peak_rss_kib": { + "median": 18560, + "min": 18448, + "max": 18648 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAkix", + "noGU": true, + "arguments": [ + "--personality=IntaRNAkix", + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.14092493208590895, + "peak_rss_kib": 18624 + }, + { + "wall_s": 0.14769695396535099, + "peak_rss_kib": 18432 + }, + { + "wall_s": 0.143615689012222, + "peak_rss_kib": 18500 + }, + { + "wall_s": 0.14164107700344175, + "peak_rss_kib": 18388 + }, + { + "wall_s": 0.1450244919396937, + "peak_rss_kib": 18300 + } + ], + "prediction": { + "start1": 160, + "end1": 169, + "start2": 64, + "end2": 73, + "E_kcal_mol": -9.13, + "hybridDB": "160||||||||||&64||||||||||", + "length_nt": 10 + }, + "wall_s": { + "median": 0.143615689012222, + "min": 0.14092493208590895, + "max": 0.14769695396535099 + }, + "peak_rss_kib": { + "median": 18432, + "min": 18300, + "max": 18624 + }, + "structure_reevaluation_matches": true + } + ], + "comparisons": [ + { + "noGU": false, + "delta_E_kcal_mol": 4.4, + "delta_length_nt": -42, + "wall_ratio_kix_over_default": 0.11369791519722769, + "rss_ratio_kix_over_default": 1.0039045553145336 + }, + { + "noGU": true, + "delta_E_kcal_mol": 1.42, + "delta_length_nt": -30, + "wall_ratio_kix_over_default": 0.14033921160767127, + "rss_ratio_kix_over_default": 0.993103448275862 + } + ] + } + ] +} diff --git a/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg b/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg new file mode 100644 index 00000000..25dfa059 --- /dev/null +++ b/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg @@ -0,0 +1,86 @@ + + IntaRNAkix: deterministic kinetic seed-extension recurrence + For each handler-provided seed, enumerate single stacks, double stacks and loop-plus-stack moves at both ends. Recompute complete interaction energies, choose the best strictly downhill move, and retain reportable visited prefixes. Stop when no downhill move exists. + + + + + + + + + + + + + + + + + + + + + + + + + + IntaRNAkix + Deterministic kinetic seed extension · model X, mode K · noLP extensions + + 1. Initialize one trajectory per seed + + seed + + S₀ = seed; H₀ = Hseed + Einit; E₀ = E(b₀, H₀) + + Trust the seed handler’s pairs and energy, including explicit lonely pairs. + Sₜ = actual pair chain; bₜ = (i₁, j₁, i₂, j₂); Hₜ = hybridization energy including initiation. + Indices of RNA 2 are reversed; left/right refer to the energy-index coordinate system. + + 2. Enumerate feasible extensions at both ends + + Single stack + Double stack + Loop + outward stack + + Left + Sₜ + Sₜ + Sₜ + Right + + + + + SₜSₜSₜ + + + 1 new pair; s₁ = s₂ = 0 + 2 new pairs; s₁ = s₂ = 0 + 2 new pairs; s₁ + s₂ > 0 + + sₖ = skipped bases on RNA k; 0 ≤ sₖ ≤ mₖ (its loop limit). Two-pair moves are atomic. + Enforce pairing, GU-loop policy, prediction ranges and both span limits. No isolated loop closure. + New boundaries advance by sₖ + 1 or sₖ + 2; the opposite boundaries stay fixed. + + 3. Recompute full energies and choose the best downhill move + + H′ = Hₜ + Eloop + Eoutward stack, if two pairs + ΔE(c) = E(b′, H′) − Eₜ; C₋ = { feasible finite moves c : ΔE(c) < 0 } + c* = arg minc ∈ C₋ ΔE(c) / D(c) + + D = 1 (score A, default), 1 + s₁ + s₂ (B), or 1 + 2 max(s₁, s₂) (C). + Ties: left end, smaller s₁ + s₂, smaller s₁, then single pair. Compare exact ratios. + E includes ED₁, ED₂, terminal terms, both weighted dangles and the additive term. + Refresh full energies at both ends, even when pairing and local energies are cached. + + 4. Commit, retain and repeat + |C₋| > 0: (Sₜ₊₁, Hₜ₊₁, Eₜ₊₁) = (Sₜ extended by c*, H′, E(b′, H′)). + Retain reportable visited states, including the seed; repeat step 2. If C₋ is empty, stop. + Deduplicate boundaries, rank by full E, apply output constraints and restore the actual pair chain. + + Greedy trajectory heuristic: no global-MFE guarantee, transition rates or physical time axis. + + diff --git a/src/IntaRNA/InteractionEnergyVrna.h b/src/IntaRNA/InteractionEnergyVrna.h index 4dd73c5b..5374f1fb 100644 --- a/src/IntaRNA/InteractionEnergyVrna.h +++ b/src/IntaRNA/InteractionEnergyVrna.h @@ -266,13 +266,6 @@ class InteractionEnergyVrna: public InteractionEnergy { Z_type getRT() const; - /** - * Read-only access to the active, temperature-scaled nearest-neighbor - * parameters, e.g. for precomputing local extension-energy estimates. - * @return parameters owned by this energy model, valid for its lifetime - */ - const vrna_param_t & getVrnaParams() const; - protected: @@ -339,15 +332,6 @@ class InteractionEnergyVrna: public InteractionEnergy { //////////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////////// -inline -const vrna_param_t & -InteractionEnergyVrna::getVrnaParams() const -{ - return *foldParams; -} - -//////////////////////////////////////////////////////////////////////////// - inline E_type InteractionEnergyVrna:: diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index f4ec6f61..5f27f746 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -73,7 +73,6 @@ libIntaRNA_a_HEADERS = \ PredictorMfe2dSeedExtension.h \ PredictorMfe2dSeedExtensionRIblast.h \ PredictorSeedExtensionKinetic.h \ - PredictorSeedExtensionKineticPruned.h \ PredictorMfe2dHeuristic.h \ PredictorMfe2dHeuristicSeed.h \ PredictorMfe2dHelixBlockHeuristic.h \ @@ -137,7 +136,6 @@ libIntaRNA_a_SOURCES = \ PredictorMfe2dSeedExtension.cpp \ PredictorMfe2dSeedExtensionRIblast.cpp \ PredictorSeedExtensionKinetic.cpp \ - PredictorSeedExtensionKineticPruned.cpp \ PredictorMfe2dHeuristic.cpp \ PredictorMfe2dHeuristicSeed.cpp \ PredictorMfe2dHelixBlockHeuristic.cpp \ diff --git a/src/IntaRNA/PredictorSeedExtensionKinetic.cpp b/src/IntaRNA/PredictorSeedExtensionKinetic.cpp index c7990619..dba6cd3a 100644 --- a/src/IntaRNA/PredictorSeedExtensionKinetic.cpp +++ b/src/IntaRNA/PredictorSeedExtensionKinetic.cpp @@ -139,14 +139,6 @@ PredictorSeedExtensionKinetic::getBoundary(const Interaction & interaction) cons ////////////////////////////////////////////////////////////////////////// -bool -PredictorSeedExtensionKinetic::prune(const Candidate &, const Boundary &) const -{ - return false; -} - -////////////////////////////////////////////////////////////////////////// - void PredictorSeedExtensionKinetic::extendSeed(Interaction & interaction, E_type hybrid, const size_t last1, const size_t last2) @@ -237,22 +229,11 @@ PredictorSeedExtensionKinetic::updateCandidates(SideCandidates & side, { // Phase one: each position pair is tested at most once per unchanged end, // even when it is the closing pair of one move and outer pair of another. - size_t stopGap2 = std::numeric_limits::max(); for (Candidate & c : side.moves) { c.bounds[c.left ? 1 : 0] = bounds[c.left ? 1 : 0]; c.bounds[c.left ? 3 : 2] = bounds[c.left ? 3 : 2]; c.active = c.bounds[1]-c.bounds[0]+1 <= energy.getAccessibility1().getMaxLength() && c.bounds[3]-c.bounds[2]+1 <= energy.getAccessibility2().getMaxLength(); - if (c.active && c.macro) { - if (c.s2 >= stopGap2) { - c.active = false; - } else if (prune(c, bounds)) { - // A suffix bound rejects this rectangle of larger gaps without - // further ED, complementarity or loop-energy lookups. - stopGap2 = c.s2; - c.active = false; - } - } if (!c.active || c.topologyKnown) { continue; } diff --git a/src/IntaRNA/PredictorSeedExtensionKinetic.h b/src/IntaRNA/PredictorSeedExtensionKinetic.h index 4af54a87..fd996428 100644 --- a/src/IntaRNA/PredictorSeedExtensionKinetic.h +++ b/src/IntaRNA/PredictorSeedExtensionKinetic.h @@ -30,9 +30,8 @@ namespace IntaRNA { * is unsupported. * * Candidate enumeration uses the active energy model and separate loop/span - * limits for both RNAs. It deliberately has no loop-only energy pruning: - * such bounds omit favorable mandatory stacks and changes to the opposite - * dangling end, and accessibility differences need not be monotone. + * limits for both RNAs. Every feasible move is evaluated with its complete + * energy change. */ class PredictorSeedExtensionKinetic : public PredictorMfe { public: @@ -78,6 +77,7 @@ class PredictorSeedExtensionKinetic : public PredictorMfe { */ void getNextBest(Interaction & interaction) override; +private: //! Inclusive boundaries (i1,j1,i2,j2), using local energy indices. using Boundary = std::array; //! Best actual path for each visited, reportable set of boundaries. @@ -102,17 +102,6 @@ class PredictorSeedExtensionKinetic : public PredictorMfe { std::int64_t delta = 0; }; - /** - * Optional filter before pairing and loop-energy lookups. The default - * enumerates every move; subclasses may implement heuristic pruning. - * @param candidate geometrically valid extension - * @param bounds current boundaries - * @return whether to skip this two-pair move and all moves with both gaps - * at least as large, in the current state (requires monotone ED) - */ - virtual bool prune(const Candidate & candidate, const Boundary & bounds) const; - -private: /** Geometry, shared pair checks and local energies for one unchanged end. */ struct SideCandidates { std::vector moves; diff --git a/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp b/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp deleted file mode 100644 index b55c3db9..00000000 --- a/src/IntaRNA/PredictorSeedExtensionKineticPruned.cpp +++ /dev/null @@ -1,110 +0,0 @@ -#include "IntaRNA/PredictorSeedExtensionKineticPruned.h" - -#include "IntaRNA/InteractionEnergyBasePair.h" -#include "IntaRNA/InteractionEnergyVrna.h" - -#include -#include -#include - -namespace IntaRNA { - -PredictorSeedExtensionKineticPruned::PredictorSeedExtensionKineticPruned( - const InteractionEnergy & model, OutputHandler & output, - PredictionTracker * predTracker, SeedHandler * seedHandler, const char score) - : PredictorSeedExtensionKinetic(model, output, predTracker, seedHandler, score) -{ - // Parameter-based bounds cannot describe arbitrary overrides of getE or - // getE_interLeft. Do not silently apply them to derived energy models. - const bool vrna = typeid(model) == typeid(InteractionEnergyVrna); - const bool basePair = typeid(model) == typeid(InteractionEnergyBasePair); - if ((!vrna && !basePair) || model.getMaxInternalLoopSize1() > 30 - || model.getMaxInternalLoopSize2() > 30) { - return; - } - rows = std::min(model.getMaxInternalLoopSize1(), model.size1() > 2 ? model.size1()-2 : 0)+1; - columns = std::min(model.getMaxInternalLoopSize2(), model.size2() > 2 ? model.size2()-2 : 0)+1; - // E_IntLoop takes a mutable pointer in supported ViennaRNA versions. - // Work on a copy and retain no pointer to the model's parameter storage. - vrna_param_t params; - if (vrna) params = static_cast(model).getVrnaParams(); - constexpr int reverse[] = {0,2,1,4,3,6,5}; - for (size_t side = 0; side < 2; ++side) { - for (int root = 1; root <= 6; ++root) { - auto & table = lowerBounds[side*6+root-1]; - table.assign(rows*columns, E_INF); - for (size_t s1 = 0; s1 < rows; ++s1) { - for (size_t s2 = 0; s2 < columns; ++s2) { - E_type & bound = table[s1*columns+s2]; - if (basePair) { - const std::int64_t two = 2*std::int64_t(model.getE_init()); - bound = static_cast(std::clamp(two, std::int64_t(std::numeric_limits::min()), std::int64_t(E_INF))); - continue; - } - for (int close = 1; close <= 6; ++close) { - E_type stack = E_INF; - for (int outer = 1; outer <= 6; ++outer) { - stack = std::min(stack, static_cast(side == 0 - ? params.stack[outer][reverse[close]] : params.stack[close][reverse[outer]])); - } - // Relax sequence consistency of the four mismatch bases. - // Include unknown nucleotide code 0, which can occur inside - // loops. This can only lower the local loop/stack minimum. - for (int a = 0; a <= 4; ++a) for (int b = 0; b <= 4; ++b) - for (int c = 0; c <= 4; ++c) for (int d = 0; d <= 4; ++d) { - const E_type loop = E_IntLoop(s1, s2, - side == 0 ? close : root, reverse[side == 0 ? root : close], a,b,c,d, ¶ms); - if (E_isNotINF(loop) && E_isNotINF(stack)) { - const std::int64_t local = std::int64_t(loop)+stack; - bound = std::min(bound, static_cast(std::clamp(local, - std::int64_t(std::numeric_limits::min()), std::int64_t(E_INF)))); - } - } - } - } - } - // Componentwise suffix minima make the bound nondecreasing along - // either gap. With monotone ED this permits rectangle pruning. - for (size_t i = rows; i-- > 0;) for (size_t j = columns; j-- > 0;) { - E_type & value = table[i*columns+j]; - if (i+1 < rows) value = std::min(value, table[(i+1)*columns+j]); - if (j+1 < columns) value = std::min(value, table[i*columns+j+1]); - } - } - } -} - -////////////////////////////////////////////////////////////////////////// - -void -PredictorSeedExtensionKineticPruned::predict(const IndexRange & r1, const IndexRange & r2) -{ - edKnown = false; - PredictorSeedExtensionKinetic::predict(r1, r2); -} - -////////////////////////////////////////////////////////////////////////// - -bool -PredictorSeedExtensionKineticPruned::prune(const Candidate & candidate, const Boundary & bounds) const -{ - if (rows == 0) return false; - // Use global coordinates for the ED memo, including repeated predict() - // calls with different range offsets but identical local boundaries. - const Boundary global{bounds[0]+energy.getOffset1(), bounds[1]+energy.getOffset1(), - bounds[2]+energy.getOffset2(), bounds[3]+energy.getOffset2()}; - if (!edKnown || edBounds != global) { - edKnown = true; edBounds = global; - currentED = std::int64_t(energy.getED1(bounds[0],bounds[1])) + energy.getED2(bounds[2],bounds[3]); - } - const int root = BP_pair[energy.getAccessibility1().getSequence().asCodes().at(global[candidate.left ? 0 : 1])] - [energy.getAccessibility2().getSequence().asCodes().at(global[candidate.left ? 2 : 3])]; - if (root < 1 || root > 6 || candidate.s1 >= rows || candidate.s2 >= columns) return false; - const E_type bound = lowerBounds[(candidate.left ? 0 : 6)+root-1][candidate.s1*columns+candidate.s2]; - const Boundary & b = candidate.bounds; - const E_type ed1 = energy.getED1(b[0],b[1]), ed2 = energy.getED2(b[2],b[3]); - return E_isINF(bound) || ed1 >= Accessibility::ED_UPPER_BOUND || ed2 >= Accessibility::ED_UPPER_BOUND - || std::int64_t(bound)+ed1+ed2-currentED >= 0; -} - -} // namespace IntaRNA diff --git a/src/IntaRNA/PredictorSeedExtensionKineticPruned.h b/src/IntaRNA/PredictorSeedExtensionKineticPruned.h deleted file mode 100644 index 50d48f33..00000000 --- a/src/IntaRNA/PredictorSeedExtensionKineticPruned.h +++ /dev/null @@ -1,69 +0,0 @@ -#ifndef INTARNA_PREDICTORSEEDEXTENSIONKINETICPRUNED_H_ -#define INTARNA_PREDICTORSEEDEXTENSIONKINETICPRUNED_H_ - -#include "IntaRNA/PredictorSeedExtensionKinetic.h" - -#include -#include -#include - -namespace IntaRNA { - -/** - * Experimental kinetic extension with root-pair-specific loop/stack bounds. - * - * Before complementarity checks, a suffix-minimum table and the current ED - * increment can reject a rectangle of larger two-pair moves. Local estimates - * use the active ViennaRNA parameters (or the base-pair energy). Pruning - * assumes monotone ED and ignores changes in terminal/dangling contributions: - * it is intentionally heuristic and can change the path compared with K. - * Every surviving move is still checked with its complete energy. Single - * stacks are always evaluated. Unknown energy subclasses and loop limits - * above 30 fall back to exhaustive K enumeration. - */ -class PredictorSeedExtensionKineticPruned : public PredictorSeedExtensionKinetic { -public: - /** - * Constructs the experimental predictor and precomputes local bounds. - * @param energy energy model that must outlive this predictor - * @param output output handler that must outlive this predictor - * @param predTracker owned prediction tracker, or NULL - * @param seedHandler owned, non-NULL seed handler - * @param score move ranking A, B or C - */ - PredictorSeedExtensionKineticPruned(const InteractionEnergy & energy, - OutputHandler & output, PredictionTracker * predTracker, - SeedHandler * seedHandler, char score = 'A'); - - /** - * Resets the ED memo and runs kinetic prediction within the given ranges. - * @param r1 permitted inclusive range in sequence 1 - * @param r2 permitted inclusive range in reversed sequence 2 - */ - void predict(const IndexRange & r1 = IndexRange(0, RnaSequence::lastPos), - const IndexRange & r2 = IndexRange(0, RnaSequence::lastPos)) override; - -protected: - /** - * Tests the local suffix bound plus ED increase, omitting terminal and - * dangling changes. A true result rejects all componentwise larger gaps. - * @param candidate geometrically valid two-pair extension - * @param bounds current interaction boundaries - * @return whether to prune the candidate and its larger-gap rectangle - */ - bool prune(const Candidate & candidate, const Boundary & bounds) const override; - -private: - //! Left/right tables for the six oriented canonical base-pair types. - std::array, 12> lowerBounds; - //! Rectangular gap-table dimensions; zero disables the heuristic. - size_t rows = 0, columns = 0; - //! Current-state ED is shared by all candidates at both ends. - mutable Boundary edBounds = {}; - mutable bool edKnown = false; - mutable std::int64_t currentED = 0; -}; - -} // namespace IntaRNA - -#endif /* INTARNA_PREDICTORSEEDEXTENSIONKINETICPRUNED_H_ */ diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index bf2e1533..d8fda9a2 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -51,7 +51,6 @@ extern "C" { #include "IntaRNA/PredictorMfe2dSeedExtension.h" #include "IntaRNA/PredictorMfe2dSeedExtensionRIblast.h" #include "IntaRNA/PredictorSeedExtensionKinetic.h" -#include "IntaRNA/PredictorSeedExtensionKineticPruned.h" #include "IntaRNA/PredictorMfe2dHeuristicSeedExtension.h" #include "IntaRNA/PredictorMfeEnsSeedOnly.h" @@ -185,7 +184,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) temperature("temperature",0,100,37), model("model", "SPBX", 'X'), - mode("mode", "HMSRKL", 'H'), // R for RIblast heuristic only + mode("mode", "HMSRK", 'H'), // R for RIblast heuristic only kineticScore("kineticScore", "ABC", 'A'), #if INTARNA_MULITHREADING threads("threads", 0, omp_get_max_threads(), 1), @@ -333,6 +332,12 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) resetParamDefault<>(accL, 0); resetParamDefault<>(intLenMax, 60); break; + case IntaRNAkix : + // deterministic kinetic seed extension + resetParamDefault<>(model, 'X'); + resetParamDefault<>(mode, 'K'); + resetParamDefault<>(outNoLP, true, "outNoLP"); + break; case IntaRNAseed : // seed-only prediction resetParamDefault<>(mode, 'S'); @@ -807,14 +812,13 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) "\n 'H' = heuristic (fast and low memory), " "\n 'M' = exact (slow), " "\n 'S' = seed-only, " - "\n 'K' = downhill greedy seed extension (requires --model=X; always noLP extensions), " - "\n 'L' = experimental K with local-energy/monotone-ED pruning (may change paths; no time or rate prediction)" + "\n 'K' = downhill greedy seed extension (IntaRNAkix; requires --model=X; always noLP extensions)" ).c_str()) (kineticScore.name.c_str() , value(&(kineticScore.val)) ->default_value(kineticScore.def) ->notifier(boost::bind(&CommandLineParsing::validate_charArgument,this,kineticScore,_1)) - , "candidate score for --model=X --mode=K or L: 'A' = complete interaction energy change, " + , "candidate score for --model=X --mode=K: 'A' = complete interaction energy change, " "'B' = change/(1+s1+s2), 'C' = change/(1+2*max(s1,s2)), where s1/s2 are skipped bases. " "All modes accept strictly negative energy changes only; B/C are heuristic scores.") (model.name.c_str() @@ -1237,10 +1241,10 @@ parse(int argc, char** argv) // parsing escape literals outSep = unescaped_string::getUnescaped( outSep ); - // K/L need a seed even before the usual --noSeed model normalization. - if (mode.val == 'K' || mode.val == 'L') { - if (model.val != 'X') throw error("--mode=K/L is available only with --model=X"); - if (noSeedRequired) throw error("--mode=K/L requires seeds and is incompatible with --noSeed"); + // K needs a seed even before the usual --noSeed model normalization. + if (mode.val == 'K') { + if (model.val != 'X') throw error("--mode=K is available only with --model=X"); + if (noSeedRequired) throw error("--mode=K requires seeds and is incompatible with --noSeed"); if (!outNoLP) { LOG(INFO) << "--mode=" << mode.val << " uses no-lonely-pair extensions: setting --outNoLP=true (handler-provided seeds are unchanged)"; outNoLP = true; @@ -1253,10 +1257,10 @@ parse(int argc, char** argv) || !outPrefix2streamName.at(OutPrefixCode::OP_qSpotProb).empty() || !outPrefix2streamName.at(OutPrefixCode::OP_tSpotProb).empty()) { - throw error("--mode=K/L does not support equilibrium ensemble or interaction-probability output"); + throw error("--mode=K does not support equilibrium ensemble or interaction-probability output"); } } else if (vm.count(kineticScore.name) && !vm.at(kineticScore.name).defaulted()) { - throw error("--kineticScore requires --model=X --mode=K or L"); + throw error("--kineticScore requires --model=X --mode=K"); } // open output stream @@ -2615,7 +2619,6 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const case 'M' : return new PredictorMfe2dSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); case 'R' : return new PredictorMfe2dSeedExtensionRIblast( energy, output, predTracker, getSeedHandler( energy ) ); case 'K' : return new PredictorSeedExtensionKinetic( energy, output, predTracker, getSeedHandler( energy ), kineticScore.val ); - case 'L' : return new PredictorSeedExtensionKineticPruned( energy, output, predTracker, getSeedHandler( energy ), kineticScore.val ); case 'S' : return new PredictorMfeSeedOnly( energy, output, predTracker, getSeedHandler( energy ) ); default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented"); return NULL; } @@ -2965,6 +2968,9 @@ getPersonality( int argc, char ** argv ) } // parse personality + if (value == "IntaRNAkix") { + return Personality::IntaRNAkix; + } if (value == "IntaRNAeval") { return Personality::IntaRNAeval; } diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index 36991034..75205e36 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -60,6 +60,7 @@ class CommandLineParsing { IntaRNA3, // default IntaRNA v3 setup IntaRNAens, // ensemble-based prediction IntaRNAeval, // evaluate predefined interactions + IntaRNAkix, // kinetic seed extension IntaRNAsTar, // sRNA-target prediction (optimized parameter) IntaRNAseed, // seed-only predictions IntaRNAhelix, // helix-block-based predictions @@ -84,6 +85,7 @@ class CommandLineParsing { case IntaRNA3 : return "IntaRNA3"; case IntaRNAens : return "IntaRNAens"; case IntaRNAeval : return "IntaRNAeval"; + case IntaRNAkix : return "IntaRNAkix"; case IntaRNAsTar : return "IntaRNAsTar"; case IntaRNAseed : return "IntaRNAseed"; case IntaRNAhelix : return "IntaRNAhelix"; diff --git a/tests/PredictorSeedExtensionKinetic_test.cpp b/tests/PredictorSeedExtensionKinetic_test.cpp index 2a7308ce..9ed293d9 100644 --- a/tests/PredictorSeedExtensionKinetic_test.cpp +++ b/tests/PredictorSeedExtensionKinetic_test.cpp @@ -7,7 +7,6 @@ #include "IntaRNA/InteractionEnergyVrna.h" #include "IntaRNA/OutputHandler.h" #include "IntaRNA/PredictorSeedExtensionKinetic.h" -#include "IntaRNA/PredictorSeedExtensionKineticPruned.h" #include "IntaRNA/SeedHandlerExplicit.h" #include "IntaRNA/VrnaHandler.h" @@ -527,114 +526,3 @@ TEST_CASE("Kinetic oracle covers varied sequences and endpoint energies", "[Pred for (char score : {'A','B','C'}) checkOracle(energy,Chain{{3,3},{4,4}},out,score); } } - -TEST_CASE("Pruned kinetic mode preserves complete energies and resets offsets", "[PredictorSeedExtensionKinetic]") { - #include "testEasyLoggingSetup.icc" - RnaSequence t("t","GGGGGGGG"), q("q","CCCCCCCC"); - KineticAccessibility a(t), b(q); - for (size_t i = 0; i < 8; ++i) for (size_t j = i; j < 8; ++j) { - a.values[{i,j}] = 20*(j-i); b.values[{i,j}] = 25*(j-i); - } - ReverseAccessibility reversed(b); - VrnaHandler vrna(25,"Turner99",false,false); - InteractionEnergyVrna energy(a,reversed,vrna,2,3,false,37,false); - const Chain seed{{3,3},{4,4}}; - const auto sc = seedConstraint(seedEncoding(seed,8)); - for (char score : {'A','B','C'}) { - OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); - KineticOutput output(out); - PredictorSeedExtensionKineticPruned predictor(energy,output,NULL,new SeedHandlerExplicit(energy,sc),score); - for (const IndexRange & range : {IndexRange(0,7),IndexRange(1,6),IndexRange(2,7),IndexRange(0,7)}) { - output.interactions.clear(); predictor.predict(range,range); - const auto expected = predict(energy,seed,out,score,range,range); - REQUIRE(output.interactions.size() == expected.size()); - for (size_t i = 0; i < expected.size(); ++i) { - REQUIRE(output.interactions[i].basePairs == expected[i].basePairs); - REQUIRE(output.interactions[i].energy == chainEnergy(energy,internalChain(energy,expected[i]),out)); - } - } - } -} - -TEST_CASE("Pruned kinetic mode documents its nonmonotone ED tradeoff", "[PredictorSeedExtensionKinetic]") { - #include "testEasyLoggingSetup.icc" - KineticFixture f(6,2,2); - f.acc1.values = {{{0,2},200},{{0,3},300}}; - InteractionEnergyBasePair energy(f.acc1,f.reversed,2,2); - const Chain seed{{0,0},{1,1}}; - const auto sc = seedConstraint(seedEncoding(seed,6)); - OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); - KineticOutput output(out); - PredictorSeedExtensionKineticPruned predictor(energy,output,NULL,new SeedHandlerExplicit(energy,sc)); - predictor.predict(); - REQUIRE(output.interactions.size() == 1); - REQUIRE(output.interactions.front().basePairs.size() == 2); - REQUIRE(predict(energy,seed,out).front().basePairs.size() > 2); - // The subclass has no bound for custom energy overrides: safely use K. - KineticOutput customOutput(out); - PredictorSeedExtensionKineticPruned custom(f.energy,customOutput,NULL,new SeedHandlerExplicit(f.energy,sc)); - custom.predict(); - REQUIRE(customOutput.interactions.front().basePairs == predict(f.energy,seed,out).front().basePairs); -} - -namespace { - -class KineticPruningProbe : public PredictorSeedExtensionKineticPruned { -public: - using PredictorSeedExtensionKineticPruned::PredictorSeedExtensionKineticPruned; - bool skips(const Bounds & before, const Bounds & after, bool left, size_t s1, size_t s2) const; -}; - -bool KineticPruningProbe::skips(const Bounds & before, const Bounds & after, - bool left, size_t s1, size_t s2) const -{ - Candidate c; - c.bounds = after; c.left = left; c.s1 = s1; c.s2 = s2; c.macro = true; - return prune(c,before); -} - -} // namespace - -TEST_CASE("Pruning tables bound local moves at all oriented root types", "[PredictorSeedExtensionKinetic]") { - #include "testEasyLoggingSetup.icc" - RnaSequence t("t","CGGUAUNCGGUAU"); - std::string reversedQuery = "GCUGUANGCUGUA"; - std::reverse(reversedQuery.begin(),reversedQuery.end()); - RnaSequence q("q",reversedQuery); - KineticAccessibility a(t), b(q); - ReverseAccessibility reversed(b); - OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); - KineticOutput output(out); - size_t checked = 0; - for (const std::string & parameters : {"Turner04","Turner99"}) { - for (double temperature : {20.,37.}) { - VrnaHandler vrna(temperature,parameters,false,false); - InteractionEnergyVrna energy(a,reversed,vrna,2,3,false,0,false); - const auto sc = seedConstraint(seedEncoding(Chain{{0,0},{1,1}},t.size())); - KineticPruningProbe probe(energy,output,NULL,new SeedHandlerExplicit(energy,sc)); - for (bool left : {false,true}) for (size_t root = 0; root < t.size(); ++root) - for (size_t s1 = 0; s1 <= 2; ++s1) for (size_t s2 = 0; s2 <= 3; ++s2) { - if ((left && root < std::max(s1,s2)+2) - || (!left && root+std::max(s1,s2)+2 >= t.size())) continue; - const size_t c1 = left ? root-s1-1 : root+s1+1; - const size_t c2 = left ? root-s2-1 : root+s2+1; - const size_t o1 = left ? c1-1 : c1+1, o2 = left ? c2-1 : c2+1; - const E_type loop = left ? energy.getE_interLeft(c1,root,c2,root) - : energy.getE_interLeft(root,c1,root,c2); - const E_type stack = left ? energy.getE_interLeft(o1,c1,o2,c2) - : energy.getE_interLeft(c1,o1,c2,o2); - if (E_isINF(loop) || E_isINF(stack) || loop+stack >= 0) continue; - const Bounds before{root,root,root,root}; - const Bounds after = left ? Bounds{o1,root,o2,root} : Bounds{root,o1,root,o2}; - a.values.clear(); - // A real local move with delta=-1 must never be rejected by - // a valid local suffix bound (no endpoint terms involved). - a.values[{after[0],after[1]}] = -(loop+stack)-1; - REQUIRE_FALSE(probe.skips(before,after,left,s1,s2)); - ++checked; - } - a.values.clear(); - } - } - REQUIRE(checked > 50); -} diff --git a/tests/runKineticSeedExtension.sh b/tests/runKineticSeedExtension.sh index f33a3ade..62b1aaa2 100755 --- a/tests/runKineticSeedExtension.sh +++ b/tests/runKineticSeedExtension.sh @@ -6,8 +6,7 @@ tmp=$(mktemp -d) trap 'rm -rf "$tmp"' EXIT common=(--target=GGGGGG --query=CCCCCC --energy=B --acc=N --threads=1 --default-log-file=/dev/null) -mode="${KINETIC_MODE:-K}" -kinetic=(--model=X --mode="$mode" '--seedTQ=3||&3||') +kinetic=(--model=X --mode=K '--seedTQ=3||&3||') csv=(--outMode=C --outCsvCols=hybridDB,E) # Six base pairs have independently known energy -6 in the base-pair model. @@ -37,7 +36,7 @@ cmp "$tmp/expected" "$tmp/boundaries" printf 'hybridDB;E\n2||||&2||||;-4\n' > "$tmp/expected" cmp "$tmp/expected" "$tmp/ranges" # Signed IntaRNA coordinates skip zero: target position three is labeled 1. -"$bin" "${common[@]}" --model=X --mode="$mode" '--seedTQ=1||&7||' "${csv[@]}" \ +"$bin" "${common[@]}" --model=X --mode=K '--seedTQ=1||&7||' "${csv[@]}" \ --tIdxPos0=-2 --qIdxPos0=5 > "$tmp/shifted" printf 'hybridDB;E\n-2||||||&5||||||;-6\n' > "$tmp/expected" cmp "$tmp/expected" "$tmp/shifted" @@ -56,16 +55,16 @@ test -s "$tmp/min-energy" grep -q -- '-6' "$tmp/min-energy" # Handler-provided explicit seeds can contain lonely pairs. -"$bin" "${common[@]}" --model=X --mode="$mode" '--seedTQ=3|&3|' \ +"$bin" "${common[@]}" --model=X --mode=K '--seedTQ=3|&3|' \ "${csv[@]}" --outNoLP > "$tmp/lonely" printf 'hybridDB;E\n1|||||&1|||||;-5\n' > "$tmp/expected" cmp "$tmp/expected" "$tmp/lonely" -"$bin" --target=GG --query=UU --energy=B --acc=N --model=X --mode="$mode" \ +"$bin" --target=GG --query=UU --energy=B --acc=N --model=X --mode=K \ '--seedTQ=1||&1||' "${csv[@]}" --outNoGUend --default-log-file=/dev/null > "$tmp/gu" cmp "$tmp/empty" "$tmp/gu" # Config-file selection must behave just like command-line selection. -printf 'model=X\nmode=%s\nkineticScore=C\nseedTQ=3||&3||\n' "$mode" > "$tmp/parameters" +printf 'model=X\nmode=K\nkineticScore=C\nseedTQ=3||&3||\n' > "$tmp/parameters" "$bin" "${common[@]}" "${csv[@]}" --parameterFile="$tmp/parameters" > "$tmp/configured" cmp "$tmp/default" "$tmp/configured" @@ -77,7 +76,7 @@ thermo=(--target=AGCGACGCA --query=UGCGUCGCU --accW=0 --accL=0 --temperature=25 --threads=1 --default-log-file=/dev/null) for score in A B C; do for dangles in false true; do - "$bin" "${thermo[@]}" --model=X --mode="$mode" '--seedTQ=3|||&5|||' \ + "$bin" "${thermo[@]}" --model=X --mode=K '--seedTQ=3|||&5|||' \ --kineticScore="$score" --energyNoDangles="$dangles" --outNoLP --outNoGUend \ > "$tmp/predicted" test "$(wc -l < "$tmp/predicted")" -gt 1 @@ -96,7 +95,7 @@ expect_error() { grep -q -- "$expected" "$tmp/bad.out" "$tmp/bad.err" } for model in S P B; do - expect_error 'only with --model=X' "${common[@]}" --model="$model" --mode="$mode" + expect_error 'only with --model=X' "${common[@]}" --model="$model" --mode=K done expect_error 'incompatible with --noSeed' "${common[@]}" "${kinetic[@]}" --noSeed # Explicitly supplying the default A must be rejected outside K as well. @@ -116,7 +115,7 @@ expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" --out="spotPr # Evaluation ignores prediction controls, including invalid kinetic scores. "$bin" "${common[@]}" "${csv[@]}" '--rri=1||||||&1||||||' \ - --model=S --mode="$mode" --kineticScore=D > "$tmp/evaluation" + --model=S --mode=K --kineticScore=D > "$tmp/evaluation" cmp "$tmp/default" "$tmp/evaluation" # A missing or false flag is promoted once, with a visible INFO message. logging=("${common[@]}") @@ -132,6 +131,31 @@ for setting in absent false true; do grep -q 'INFO.*setting --outNoLP=true' "$tmp/info.log" fi done -# Run the same API/CLI contract against the experimental subclass. -if [ "$mode" = K ]; then KINETIC_MODE=L bash "$0"; fi -echo "Kinetic seed-extension CLI checks passed ($mode)" +# Both personality entry points select mode K and noLP by default. +ln -s "$bin" "$tmp/IntaRNAkix" +for invocation in binary option; do + args=() + executable="$tmp/IntaRNAkix" + if [ "$invocation" = option ]; then + executable="$bin" + args=(--personality=IntaRNAkix) + fi + : > "$tmp/info.log" + "$executable" "${args[@]}" "${logging[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ + --outNumber=10 > "$tmp/kix-prefixes" + cmp "$tmp/prefixes" "$tmp/kix-prefixes" + ! grep -q 'setting --outNoLP=true' "$tmp/info.log" + "$executable" "${args[@]}" "${common[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ + --mode=S > "$tmp/seed-only" + printf 'hybridDB;E\n3||&3||;-2\n' > "$tmp/expected" + cmp "$tmp/expected" "$tmp/seed-only" +done +# Explicitly disabling noLP cannot disable the mode K extension invariant. +"$bin" --personality=IntaRNAkix "${logging[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ + --outNoLP=false > "$tmp/kix-noLP" +cmp "$tmp/default" "$tmp/kix-noLP" +grep -q 'setting --outNoLP=true' "$tmp/info.log" +# Evaluation remains available under the personality and ignores its defaults. +"$tmp/IntaRNAkix" "${common[@]}" "${csv[@]}" '--rri=1||||||&1||||||' > "$tmp/kix-eval" +cmp "$tmp/default" "$tmp/kix-eval" +echo 'Kinetic seed-extension and IntaRNAkix CLI checks passed' From 104433645a359c4ad71c6db04175810a31fd7677 Mon Sep 17 00:00:00 2001 From: Martin Raden Date: Wed, 7 Oct 2026 11:37:56 +0200 Subject: [PATCH 5/6] rename IntaRNAkix to IntaRNAsnap --- ChangeLog | 32 +++++++++++++++---- README.md | 14 ++++---- doc/Makefile.am | 2 +- doc/benchmark-kix.py | 10 +++--- doc/kinetic-seed-extension.md | 10 +++--- doc/kix-benchmark-20261005.json | 28 ++++++++-------- ...taRNAkix.PredictorSeedExtensionKinetic.svg | 4 +-- src/bin/CommandLineParsing.cpp | 8 ++--- src/bin/CommandLineParsing.h | 4 +-- tests/runKineticSeedExtension.sh | 12 +++---- 10 files changed, 72 insertions(+), 52 deletions(-) diff --git a/ChangeLog b/ChangeLog index 47227e40..8c272627 100644 --- a/ChangeLog +++ b/ChangeLog @@ -19,7 +19,7 @@ and retained greedy traceback; equilibrium probability output is unsupported - kinetic mode always uses stacked extensions and trusts handler-provided seeds; reuse candidate pairing/local energies -- IntaRNAkix personality enables kinetic seed extension with noLP by default +- IntaRNAsnap personality enables kinetic seed extension with noLP by default - apply outDeltaE relative to the sequence pair's best interaction when merging regions; preserve local windows for outPerRegion output (PR #253) @@ -42,7 +42,7 @@ ## Technical changes and Optimizations -- install personality links in out-of-tree builds, including IntaRNAkix +- install personality links in out-of-tree builds, including IntaRNAsnap - BUGFIX : normalize single-pair suboptimal boundaries before traceback and boundary-only output validation (PR #253) @@ -85,6 +85,26 @@ - BUGFIX : multi-threading : IntaRNAsTar was not thread-safe due to shared storage of computed ranges +# IntaRNAsnap + +New personality IntaRNAsnap identifies the fastest "folding path" to a (locally) +optimal interaction after seed formation, rather than finding the global optimum +under the assumption of thermodynamic equilibrium. + +To this end, IntaRNAsnap implements a kinetics-motivated greedy search for the +best interaction that can be reached from the given seed(s) by a series of +complete-energy downhill local moves. +Local moves are direct stack extensions or 2-base-pair loop/stack moves, i.e. +the same moves as used in the IntaRNA noLP model. That way, kinetic energy +barriers posed by the formation of loops and bulges are avoided, which mimics +the continuation of the zipping process after "jumping" over a loop or bulge. + +This personality is only available with --model=X and --mode=K, and equilibrium +probability output is unsupported. The local scoring can be configured with +--kineticScore=A|B|C, and the output can include distinct structures, energies, +restricted partition sums, and trackers. + + # IntaRNAeval New personality of IntaRNA to evaluate predefined RNA-RNA interactions. @@ -108,18 +128,18 @@ energies, restricted partition sums, and trackers. 261005 Alexander Mitrofanov * bin/CommandLineParsing, tests/runKineticSeedExtension.sh : - + add IntaRNAkix personality with model=X, mode=K and outNoLP=true defaults; + + add IntaRNAsnap personality with model=X, mode=K and outNoLP=true defaults; support executable-name and explicit personality selection * configure.ac : * find personality declarations via srcdir so out-of-tree installations - create the executable links, including IntaRNAkix + create the executable links, including IntaRNAsnap * README.md, doc/kinetic-seed-extension.md, - doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg, doc/Makefile.am : + doc/recursions/IntaRNAsnap.PredictorSeedExtensionKinetic.svg, doc/Makefile.am : + document personality usage and depict seed initialization, allowed moves, complete-energy greedy selection, stopping and prefix reporting + distribute the benchmark's tutorial sequences with its script and results * doc/benchmark-kix.py, doc/kix-benchmark-20261005.json : - + compare IntaRNAkix with default IntaRNA with and without GU-end constraints; + + compare IntaRNAsnap with default IntaRNA with and without GU-end constraints; record wall time, peak RSS, MFE energy and maximum covered strand length * implement https://github.com/BackofenLab/IntaRNA/pull/254#issuecomment-5995041057 diff --git a/README.md b/README.md index 075a87e2..997daa60 100644 --- a/README.md +++ b/README.md @@ -100,7 +100,7 @@ The following topics are covered by this documentation: - [IntaRNAsTar - optimized for sRNA-target prediction](#IntaRNAsTar) - [IntaRNAseed - identifys and reports seed interactions only](#IntaRNAseed) - [IntaRNAens - ensemble-based prediction and partition function computation](#IntaRNAens) - - [IntaRNAkix - kinetic seed extension](#IntaRNAkix) + - [IntaRNAsnap - kinetic seed extension](#IntaRNAsnap) - [IntaRNAeval - evaluate predefined interactions](#IntaRNAeval) - [How to constrain predicted interactions](#constraintSetup) - [Interaction restrictions](#interConstr) @@ -756,7 +756,7 @@ traceback preserves the actual chosen path. `--outNoGUend`, separate query and target loop/span limits, regions, output energy/accessibility filters and overlap settings remain applicable. -The [IntaRNAkix personality](#IntaRNAkix) selects this mode with noLP enabled +The [IntaRNAsnap personality](#IntaRNAsnap) selects this mode with noLP enabled by default. It is a zippering-inspired heuristic without a calibrated time axis or a guarantee of the global minimum. Equilibrium probability/partition-sum outputs are rejected, as are other models and `--noSeed`. Scores B and C are @@ -1038,9 +1038,9 @@ IntaRNA --mode=S ... [![up](doc/figures/icon-up.28.png) back to overview](#overview) -### IntaRNAkix +### IntaRNAsnap -**IntaRNAkix** (kinetic seed extension) grows every handler-provided seed by +**IntaRNAsnap** (kinetic seed extension) grows every handler-provided seed by choosing the most favorable complete energy change at either end. It sets `--model=X --mode=K --outNoLP=true`; other defaults are those of IntaRNA. Each move adds one stacked pair, two stacked pairs, or a loop-closing pair @@ -1050,8 +1050,8 @@ Seeds themselves may contain lonely pairs when supplied by the seed handler. The following calls are equivalent: ```sh -IntaRNAkix -t target.fasta -q query.fasta -IntaRNA --personality=IntaRNAkix -t target.fasta -q query.fasta +IntaRNAsnap -t target.fasta -q query.fasta +IntaRNA --personality=IntaRNAsnap -t target.fasta -q query.fasta IntaRNA --model=X --mode=K --outNoLP=true -t target.fasta -q query.fasta ``` @@ -1060,7 +1060,7 @@ largest energy decrease; `--kineticScore=B|C` adds distance preferences. The reported MFE is the best visited, reportable interaction across seeds; this greedy search has no global-optimum or physical folding-time guarantee. -![IntaRNAkix initialization, allowed extensions, greedy update and stopping rule](doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg) +![IntaRNAsnap initialization, allowed extensions, greedy update and stopping rule](doc/recursions/IntaRNAsnap.PredictorSeedExtensionKinetic.svg) See the [algorithm and preliminary benchmark](doc/kinetic-seed-extension.md) for time, peak-memory, energy and interaction-length comparisons with default diff --git a/doc/Makefile.am b/doc/Makefile.am index b8c4b6bf..bb58e795 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -17,7 +17,7 @@ EXTRA_DIST = \ handson/GcvB.fasta \ handson/ilvE.fasta \ handson/GcvB.ST.fasta \ - recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg \ + recursions/IntaRNAsnap.PredictorSeedExtensionKinetic.svg \ doxygen.cfg \ latex-deps/adjcalc.sty \ latex-deps/adjustbox.sty \ diff --git a/doc/benchmark-kix.py b/doc/benchmark-kix.py index bbe3faa0..34f1d399 100644 --- a/doc/benchmark-kix.py +++ b/doc/benchmark-kix.py @@ -1,5 +1,5 @@ #!/usr/bin/env python3 -"""Small, single-thread comparison of IntaRNAkix and default IntaRNA. +"""Small, single-thread comparison of IntaRNAsnap and default IntaRNA. Requires Python 3, GNU time and a release IntaRNA binary. Run from any directory: python3 doc/benchmark-kix.py /path/to/IntaRNA --output results.json @@ -82,9 +82,9 @@ def main(): "memory": "GNU time %M: maximum resident set size of each child process, KiB", "noGU": "--outNoGUend=true; seedNoGU remains false for both programs", "baseline": "default IntaRNA (model X, mode H, outNoLP false)", - "comparison": "IntaRNAkix defaults (model X, mode K, outNoLP true, kineticScore A)", + "comparison": "IntaRNAsnap defaults (model X, mode K, outNoLP true, kineticScore A)", "length": "max(end1-start1+1, end2-start2+1), nt; default one-based coordinates", - "deviations": "signed IntaRNAkix minus default IntaRNA, within each GU setting; not a global-optimum error bound", + "deviations": "signed IntaRNAsnap minus default IntaRNA, within each GU setting; not a global-optimum error bound", "cases": [], } with tempfile.TemporaryDirectory(prefix="intarna-kix-benchmark-") as directory: @@ -97,8 +97,8 @@ def main(): "query": read_fasta(fixtures / (query + ".fasta")), "runs": [], "comparisons": []} configurations = [] for no_gu in (False, True): - for program in ("IntaRNA", "IntaRNAkix"): - options = [] if program == "IntaRNA" else ["--personality=IntaRNAkix"] + for program in ("IntaRNA", "IntaRNAsnap"): + options = [] if program == "IntaRNA" else ["--personality=IntaRNAsnap"] options.append("--outNoGUend=" + str(no_gu).lower()) configurations.append((program, no_gu, options)) case["runs"].append({"program": program, "noGU": no_gu, diff --git a/doc/kinetic-seed-extension.md b/doc/kinetic-seed-extension.md index c879a7d4..23cb48d1 100644 --- a/doc/kinetic-seed-extension.md +++ b/doc/kinetic-seed-extension.md @@ -9,14 +9,14 @@ it has no calibrated transition rates or time axis, does not cross barriers between committed states, and does not guarantee a global minimum. A favorable two-pair move does not establish a barrier-free physical reaction pathway. -The `IntaRNAkix` personality (kinetic seed extension) selects +The `IntaRNAsnap` personality (kinetic seed extension) selects `--model=X --mode=K --outNoLP=true`. It can be invoked through the installed -`IntaRNAkix` executable link or `IntaRNA --personality=IntaRNAkix`. Other defaults +`IntaRNAsnap` executable link or `IntaRNA --personality=IntaRNAsnap`. Other defaults remain those of IntaRNA. Ordinary energy trackers are supported; seedless operation, other interaction models and equilibrium partition/probability requests are rejected for mode K. -![Kinetic seed extension recursion](recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg) +![Kinetic seed extension recursion](recursions/IntaRNAsnap.PredictorSeedExtensionKinetic.svg) ## Seeds, states and allowed extensions @@ -35,7 +35,7 @@ range offsets and conversion to original coordinates. Initially, Extensions **always** use the no-lonely-pair strategy, independent of the API output constraint. Direct `--mode=K` calls promote a missing or false `--outNoLP` to true with an INFO message using the normal logging destination; -IntaRNAkix already defaults to true. This applies to new +IntaRNAsnap already defaults to true. This applies to new extensions, not to revalidation of handler-provided seeds. Allowed moves are: - One pair stacked directly onto the current boundary (`|SEED`). @@ -128,7 +128,7 @@ No experimental seed, region or accessibility constraints from those examples are applied here. The raw record includes the sequences and input hashes. The comparison uses the actual personality defaults: IntaRNA has model X, -mode H and `outNoLP=false`; IntaRNAkix has model X, mode K, score A and +mode H and `outNoLP=false`; IntaRNAsnap has model X, mode K, score A and `outNoLP=true`. Thus energy and length deviations reflect both the search and the different noLP defaults. Default IntaRNA is itself a heuristic, so the energy deviation is not a certified error from a global optimum. diff --git a/doc/kix-benchmark-20261005.json b/doc/kix-benchmark-20261005.json index 003a6a24..12f91bd6 100644 --- a/doc/kix-benchmark-20261005.json +++ b/doc/kix-benchmark-20261005.json @@ -20,9 +20,9 @@ "memory": "GNU time %M: maximum resident set size of each child process, KiB", "noGU": "--outNoGUend=true; seedNoGU remains false for both programs", "baseline": "default IntaRNA (model X, mode H, outNoLP false)", - "comparison": "IntaRNAkix defaults (model X, mode K, outNoLP true, kineticScore A)", + "comparison": "IntaRNAsnap defaults (model X, mode K, outNoLP true, kineticScore A)", "length": "max(end1-start1+1, end2-start2+1), nt; default one-based coordinates", - "deviations": "signed IntaRNAkix minus default IntaRNA, within each GU setting; not a global-optimum error bound", + "deviations": "signed IntaRNAsnap minus default IntaRNA, within each GU setting; not a global-optimum error bound", "cases": [ { "name": "fhlA/OxyS", @@ -91,10 +91,10 @@ "structure_reevaluation_matches": true }, { - "program": "IntaRNAkix", + "program": "IntaRNAsnap", "noGU": false, "arguments": [ - "--personality=IntaRNAkix", + "--personality=IntaRNAsnap", "--outNoGUend=false" ], "samples": [ @@ -190,10 +190,10 @@ "structure_reevaluation_matches": true }, { - "program": "IntaRNAkix", + "program": "IntaRNAsnap", "noGU": true, "arguments": [ - "--personality=IntaRNAkix", + "--personality=IntaRNAsnap", "--outNoGUend=true" ], "samples": [ @@ -324,10 +324,10 @@ "structure_reevaluation_matches": true }, { - "program": "IntaRNAkix", + "program": "IntaRNAsnap", "noGU": false, "arguments": [ - "--personality=IntaRNAkix", + "--personality=IntaRNAsnap", "--outNoGUend=false" ], "samples": [ @@ -423,10 +423,10 @@ "structure_reevaluation_matches": true }, { - "program": "IntaRNAkix", + "program": "IntaRNAsnap", "noGU": true, "arguments": [ - "--personality=IntaRNAkix", + "--personality=IntaRNAsnap", "--outNoGUend=true" ], "samples": [ @@ -557,10 +557,10 @@ "structure_reevaluation_matches": true }, { - "program": "IntaRNAkix", + "program": "IntaRNAsnap", "noGU": false, "arguments": [ - "--personality=IntaRNAkix", + "--personality=IntaRNAsnap", "--outNoGUend=false" ], "samples": [ @@ -656,10 +656,10 @@ "structure_reevaluation_matches": true }, { - "program": "IntaRNAkix", + "program": "IntaRNAsnap", "noGU": true, "arguments": [ - "--personality=IntaRNAkix", + "--personality=IntaRNAsnap", "--outNoGUend=true" ], "samples": [ diff --git a/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg b/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg index 25dfa059..a79b875f 100644 --- a/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg +++ b/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg @@ -1,5 +1,5 @@ - IntaRNAkix: deterministic kinetic seed-extension recurrence + IntaRNAsnap: deterministic kinetic seed-extension recurrence For each handler-provided seed, enumerate single stacks, double stacks and loop-plus-stack moves at both ends. Recompute complete interaction energies, choose the best strictly downhill move, and retain reportable visited prefixes. Stop when no downhill move exists. @@ -26,7 +26,7 @@ - IntaRNAkix + IntaRNAsnap Deterministic kinetic seed extension · model X, mode K · noLP extensions 1. Initialize one trajectory per seed diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index d8fda9a2..c06ac8b0 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -332,7 +332,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) resetParamDefault<>(accL, 0); resetParamDefault<>(intLenMax, 60); break; - case IntaRNAkix : + case IntaRNAsnap : // deterministic kinetic seed extension resetParamDefault<>(model, 'X'); resetParamDefault<>(mode, 'K'); @@ -812,7 +812,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) "\n 'H' = heuristic (fast and low memory), " "\n 'M' = exact (slow), " "\n 'S' = seed-only, " - "\n 'K' = downhill greedy seed extension (IntaRNAkix; requires --model=X; always noLP extensions)" + "\n 'K' = downhill greedy seed extension (IntaRNAsnap; requires --model=X; always noLP extensions)" ).c_str()) (kineticScore.name.c_str() , value(&(kineticScore.val)) @@ -2968,8 +2968,8 @@ getPersonality( int argc, char ** argv ) } // parse personality - if (value == "IntaRNAkix") { - return Personality::IntaRNAkix; + if (value == "IntaRNAsnap") { + return Personality::IntaRNAsnap; } if (value == "IntaRNAeval") { return Personality::IntaRNAeval; diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index 75205e36..f6658728 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -60,7 +60,7 @@ class CommandLineParsing { IntaRNA3, // default IntaRNA v3 setup IntaRNAens, // ensemble-based prediction IntaRNAeval, // evaluate predefined interactions - IntaRNAkix, // kinetic seed extension + IntaRNAsnap, // kinetic seed extension IntaRNAsTar, // sRNA-target prediction (optimized parameter) IntaRNAseed, // seed-only predictions IntaRNAhelix, // helix-block-based predictions @@ -85,7 +85,7 @@ class CommandLineParsing { case IntaRNA3 : return "IntaRNA3"; case IntaRNAens : return "IntaRNAens"; case IntaRNAeval : return "IntaRNAeval"; - case IntaRNAkix : return "IntaRNAkix"; + case IntaRNAsnap : return "IntaRNAsnap"; case IntaRNAsTar : return "IntaRNAsTar"; case IntaRNAseed : return "IntaRNAseed"; case IntaRNAhelix : return "IntaRNAhelix"; diff --git a/tests/runKineticSeedExtension.sh b/tests/runKineticSeedExtension.sh index 62b1aaa2..e6054c2d 100755 --- a/tests/runKineticSeedExtension.sh +++ b/tests/runKineticSeedExtension.sh @@ -132,13 +132,13 @@ for setting in absent false true; do fi done # Both personality entry points select mode K and noLP by default. -ln -s "$bin" "$tmp/IntaRNAkix" +ln -s "$bin" "$tmp/IntaRNAsnap" for invocation in binary option; do args=() - executable="$tmp/IntaRNAkix" + executable="$tmp/IntaRNAsnap" if [ "$invocation" = option ]; then executable="$bin" - args=(--personality=IntaRNAkix) + args=(--personality=IntaRNAsnap) fi : > "$tmp/info.log" "$executable" "${args[@]}" "${logging[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ @@ -151,11 +151,11 @@ for invocation in binary option; do cmp "$tmp/expected" "$tmp/seed-only" done # Explicitly disabling noLP cannot disable the mode K extension invariant. -"$bin" --personality=IntaRNAkix "${logging[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ +"$bin" --personality=IntaRNAsnap "${logging[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ --outNoLP=false > "$tmp/kix-noLP" cmp "$tmp/default" "$tmp/kix-noLP" grep -q 'setting --outNoLP=true' "$tmp/info.log" # Evaluation remains available under the personality and ignores its defaults. -"$tmp/IntaRNAkix" "${common[@]}" "${csv[@]}" '--rri=1||||||&1||||||' > "$tmp/kix-eval" +"$tmp/IntaRNAsnap" "${common[@]}" "${csv[@]}" '--rri=1||||||&1||||||' > "$tmp/kix-eval" cmp "$tmp/default" "$tmp/kix-eval" -echo 'Kinetic seed-extension and IntaRNAkix CLI checks passed' +echo 'Kinetic seed-extension and IntaRNAsnap CLI checks passed' From 5f3c4c185b1d9c987d7937669171896e43db9338 Mon Sep 17 00:00:00 2001 From: Martin Raden Date: Wed, 7 Oct 2026 11:44:08 +0200 Subject: [PATCH 6/6] dist + readme cleanup --- README.md | 7 +------ doc/Makefile.am | 13 ------------- 2 files changed, 1 insertion(+), 19 deletions(-) diff --git a/README.md b/README.md index 997daa60..ed4cb9e3 100644 --- a/README.md +++ b/README.md @@ -1058,13 +1058,8 @@ IntaRNA --model=X --mode=K --outNoLP=true -t target.fasta -q query.fasta Only strictly downhill moves are accepted. The default score A chooses the largest energy decrease; `--kineticScore=B|C` adds distance preferences. The reported MFE is the best visited, reportable interaction across seeds; -this greedy search has no global-optimum or physical folding-time guarantee. +this heuristic greedy search has no global-optimum or physical folding-time guarantee. -![IntaRNAsnap initialization, allowed extensions, greedy update and stopping rule](doc/recursions/IntaRNAsnap.PredictorSeedExtensionKinetic.svg) - -See the [algorithm and preliminary benchmark](doc/kinetic-seed-extension.md) -for time, peak-memory, energy and interaction-length comparisons with default -IntaRNA, with and without GU-end restrictions. [![up](doc/figures/icon-up.28.png) back to overview](#overview) diff --git a/doc/Makefile.am b/doc/Makefile.am index bb58e795..c28d05ec 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -4,20 +4,7 @@ ################################################################ EXTRA_DIST = \ - analysis/out-overlap.md \ - analysis/out-overlap/reproduce.py \ conda.txt \ - kinetic-seed-extension.md \ - benchmark-kix.py \ - kix-benchmark-20261005.json \ - handson/README.md \ - handson/fhlA.fasta \ - handson/OxyS.fasta \ - handson/phoB.fasta \ - handson/GcvB.fasta \ - handson/ilvE.fasta \ - handson/GcvB.ST.fasta \ - recursions/IntaRNAsnap.PredictorSeedExtensionKinetic.svg \ doxygen.cfg \ latex-deps/adjcalc.sty \ latex-deps/adjustbox.sty \