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