diff --git a/ChangeLog b/ChangeLog index 62b5f84c..07e8fc1c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,14 @@ ## Interface and handling +- apply outDeltaE relative to the sequence pair's best interaction when merging + regions; preserve local windows for outPerRegion output (PR #253) + +- reject merged regions that cannot enforce the requested output overlap mode; + retain independent selection with --outPerRegion (PR #253) + +- clarify that suboptimal predictions have distinct interaction-site boundaries + - IntaRNAeval / --rri evaluates predefined RNA-RNA interactions (issue #184) - compressed binary .agz accessibility caches for repeated screens (issue #245) @@ -27,6 +35,22 @@ ## Technical changes and Optimizations +- BUGFIX : normalize single-pair suboptimal boundaries before traceback and + boundary-only output validation (PR #253) + +- verify restricted heuristic traces against the independent tiny oracle and + cover output overlap checks after outMinPu decomposition (PR #253) + +- BUGFIX : apply terminal GU and accessibility output limits to later heuristic + results, including seeded and helix predictors (PR #253) + +- BUGFIX : heuristic ensemble suboptimals use finalized site energies and stop + safely when no compatible candidate remains (PR #253) + +- analyse output overlap modes and implement the reviewed regional guards, + ensemble/output-filter fixes and global energy window; retain and document + the one-right-extension limitation (issue #212, PR #253) + - serialize stored accessibility rows directly and load retained rows in place; benchmark direct/generic export and raw/gzip binary I/O (PR #250) - add repository-specific AI coding guidance in AGENTS.md (issue #247) @@ -73,6 +97,82 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261005 Alexander Mitrofanov + * IntaRNA/PredictorMfe, tests/PredictorMfeHeuristicCellState_test.cpp, + tests/PredictorMfeEnsRegression_test.cpp : + * normalize coincident suboptimal boundaries to one base pair before + traceback or output, matching the initial optimum representation + + assert valid single-pair candidates in traced and boundary-only output; + fixes failures exposed by the PR #253 debug/sanitizer regressions + +261005 Alexander Mitrofanov + * tests/PredictorTinyOracle_test.cpp, tests/runOutputOverlap.sh : + + verify every reported heuristic trace and energy against independently + enumerated structures filtered by GU/ED limits; check forbidden intervals + + exercise automatic outMinPu region splitting on both RNAs in all overlap + modes with merged and independent output (PR #253) + +261005 Alexander Mitrofanov + * README.md : + * retain and explain the one-right-extension strategy, including mode=M, + incomplete compatible alternatives and the limits of blocking paired bases + * doc/analysis/out-overlap.md, doc/analysis/out-overlap/reproduce.py : + + record the accepted review decisions separately from historical findings + * capture old output and new regional input rejections without treating + historical bugs as regression expectations (PR #253) + +261005 Alexander Mitrofanov + * bin/IntaRNA, bin/CommandLineParsing, README.md : + * filter merged results by the global minimum plus outDeltaE; keep the + independent energy windows for outPerRegion=true (PR #253) + * tests/runOutputOverlap.sh : + + verify both region orders, zero/inclusive delta boundaries, permitted + overlap modes, independent regional output and empty merged results + +261005 Alexander Mitrofanov + * IntaRNA/PredictorMfe, PredictorMfeEns, PredictorMfe2dHeuristic, + PredictorMfe2dHeuristicSeed, PredictorMfe2dHelixBlockHeuristic, + PredictorMfe2dHelixBlockHeuristicSeed : + * share complete-site GU/ED validation between initial and later candidates + * retain the existing matrix scan order and one-extension heuristic (PR #253) + * tests/PredictorMfeHeuristicCellState_test.cpp : + + cover all four matrix selectors, terminal GU and nonzero ED limits, + all overlap modes, region offsets and boundary/traced output + +261005 Alexander Mitrofanov + * bin/CommandLineParsing, README.md : + + validate manual regions and automatic decomposition against outOverlap + + N requires one region per RNA, T one target region, Q one query region; + B and independent outPerRegion output accept multiple regions (PR #253) + * tests/runOutputOverlap.sh, tests/Makefile.am : + + exercise all overlap modes, manual/automatic regions, per-region output, + unchanged single regions and shifted sequence indices + +261005 Alexander Mitrofanov + * IntaRNA/PredictorMfeEns2dHeuristic : + - remove raw-matrix suboptimal selection and its floating-point sentinel + * reuse validated, finalized per-site energies and E_INF exhaustion from + PredictorMfe; retain one best right extension per left boundary (PR #253) + * tests/PredictorMfeEnsRegression_test.cpp : + + compare N/T/Q energies against B with accessibility and energy cutoffs + + cover exhausted and empty candidate sets with and without traceback + +261005 Alexander Mitrofanov + * README.md : + + define suboptimal predictions by distinct start/end coordinates on both RNAs + * distinguish site prediction from evaluation of supplied structures (PR #253) + +261002 Alexander Mitrofanov + * doc/analysis/out-overlap.md, doc/analysis/out-overlap/reproduce.py : + + analyse all four output overlap modes and the documented enumeration limits + + reproduce regional overlap and energy-window violations, heuristic output + filtering gaps, and invalid/incomplete ensemble candidate energies + + provide 118 CLI observations, 64 controls and three independent site oracles + + propose selection semantics and repairs for consultation; no predictor changes + * addresses analysis request in https://github.com/BackofenLab/IntaRNA/issues/212 + * doc/Makefile.am : + + distribute the analysis and standalone reproduction script + 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..6466b700 100644 --- a/README.md +++ b/README.md @@ -1129,7 +1129,19 @@ and should be in the format `from1-end1,from2-end2,..` using integers. Note, if you want to have predictions individually for each region combination (rather than just the best for each query-target combination) you -want to add `--outPerRegion` to the call. +want to add `--outPerRegion` to the call. Overlap restrictions then apply +independently within each region combination; results from different combinations +can overlap on either RNA. + +With the default `--outPerRegion=false`, region combinations are merged. To +avoid reporting forbidden overlaps, `--outOverlap=N` requires a single region +on each RNA, `T` requires a single target region, and `Q` requires a single query +region. `B` permits any number of regions. For example, multiple target regions +can reuse the same query interval, so they are allowed with `Q` or `B`, but +rejected with `N` or `T`. The same checks apply to regions produced automatically +by `--qRegionLenMax`, `--tRegionLenMax`, or `--outMinPu`. Use +`--outPerRegion=true` to select interactions independently for multiple regions +with any overlap mode. If you are dealing with very long sequences it might be useful to use the *automatic identification of accessible regions*, which dramatically reduces @@ -1666,19 +1678,27 @@ interaction energy = -6.39 kcal/mol Besides the identification of the optimal (e.g. minimum-free-energy) RNA-RNA interaction, IntaRNA enables the enumeration of suboptimal interactions. To this end, the argument `-n N` or `--outNumber=N` can be used to generate up to `N` -interactions for each query-target pair (including the optimal one). +interactions for each query-target pair (including the optimal one). Reported +interactions have distinct interaction-site boundaries: at least one of the four +indices `start1`, `end1`, `start2`, or `end2` differs. Different internal base-pair +patterns with the same boundaries are not enumerated as separate predictions. +This restriction applies to prediction; [interaction evaluation](#intarnaeval) +can evaluate distinct supplied structures with identical boundaries. *Note*: suboptimal interaction enumeration is not exhaustive! That is, for each interaction site (defined by the left- and right-most intermolecular base pair) -only the best interaction is reported! In heuristic prediction mode (default -mode of IntaRNA), this is even less exhaustive, since only for each left-most -interaction boundary one interaction is reported! +only one prediction is reported. Heuristic predictors additionally prune +candidate extensions during computation and can therefore consider fewer +interaction sites. Furthermore, it is possible to *restrict (sub)optimal enumeration* using - `--outMaxE` : maximal energy for any interaction reported - `--outDeltaE` : maximal energy difference of suboptimal interactions' energy - to the minimum free energy interaction + to the minimum free energy interaction for the query-target pair, including + when results from multiple regions are merged. With `--outPerRegion=true`, + the minimum and energy window are determined independently for each region + combination - `--outOverlap` : defines if and where overlapping of reported interaction sites is allowed: - 'N' : no overlap neither in target nor query allowed for reported interactions @@ -1686,11 +1706,21 @@ Furthermore, it is possible to *restrict (sub)optimal enumeration* using - 'T' : overlap allowed for interacting subsequences in target only - 'Q' : overlap allowed for interacting subsequences in query only -*Note*: non-overlapping output (i) is heuristic by considering for each left -interaction site only the best right extension for overlap computation and -(ii) increases runtime. To get optimized results of non-overlapping suboptimals, -rerun IntaRNA and mark the optimal (mfe) interaction region as -[blocked](#accConstraints). +Overlap refers to the entire interval between the outermost intermolecular +base pairs, including unpaired positions inside that interval. + +*Note*: non-overlapping output is heuristic, even with exact prediction +(`--mode=M`), and increases runtime. For each left interaction boundary, only +the best right extension is retained for overlap selection. If that extension +overlaps an earlier result, a shorter compatible extension may already have +been discarded. Thus fewer than `outNumber` results does not imply that no +further compatible interaction exists. This strategy keeps storage bounded. + +Rerunning IntaRNA with the optimal interaction region +[blocked](#accConstraints) may reveal additional alternatives. Blocking prevents +base pairing at those positions; an interaction can still span the blocked +region through an internal loop, so this does not guarantee interval-disjoint +results. diff --git a/doc/Makefile.am b/doc/Makefile.am index 322e9d93..ac9bab1a 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -4,6 +4,8 @@ ################################################################ EXTRA_DIST = \ + analysis/out-overlap.md \ + analysis/out-overlap/reproduce.py \ conda.txt \ doxygen.cfg \ latex-deps/adjcalc.sty \ diff --git a/doc/analysis/out-overlap.md b/doc/analysis/out-overlap.md new file mode 100644 index 00000000..fa48fa81 --- /dev/null +++ b/doc/analysis/out-overlap.md @@ -0,0 +1,459 @@ +# Output overlap analysis for issue 212 + +This report examines `--outOverlap=B,N,T,Q` at master revision +`afede020fe7d1c44cf686684dc7eb40927f91583` (2026-10-02). It answers the +[request for analysis and consultation](https://github.com/BackofenLab/IntaRNA/issues/212#issuecomment-5927769865). +The findings, example outputs, and original proposals below describe that +baseline. The implementation following review is summarized next; the historical +examples are not expected output for the revised program. + +## Implementation following Martin's review + +The [review on PR #253](https://github.com/BackofenLab/IntaRNA/pull/253#issuecomment-5991808858) +requested separate implementation steps, which are recorded in separate commits: + +- Define distinct prediction results by their four interaction-site boundaries + in the README. Alternative internal structures with the same boundaries are + not enumerated separately. +- Reject unsupported merged-region output in `CommandLineParsing`, both for + explicit regions and after automatic decomposition (including `outMinPu`). + `N` requires one region per RNA, `T` one target region, and `Q` one query + region. `B` and explicit `outPerRegion=true` permit multiple regions. + This follows the actual overlap semantics: multiple target regions can reuse + the query, so they are unsafe for `T`; the converse applies to `Q`. +- Use finalized per-site ensemble energies for later `P/H --noSeed` results, + through the base predictor's existing best-site-per-left-boundary storage. + This also removes the floating-point exhaustion sentinel in finding 2. +- Share complete-site terminal GU and maximum ED validation between initial + selection and all four specialized MFE/helix matrix selectors. Ensemble + partition accumulation uses the same check. Debug validation also exposed + coincident boundaries for single-pair suboptimals; these are normalized before + traceback or boundary-only output. +- Apply `outDeltaE` relative to the best result for the sequence pair when + merging regional output. `outPerRegion=true` retains independent local + energy windows and overlap selection. +- Keep the one-right-extension strategy in finding 6, and clarify its limits + in the README, including its use with exact prediction mode. + +The new regional checks deliberately replace the proposed global overlap +selector; no exhaustive global selection or recovery of discarded extensions +is introduced. Existing specialized MFE/helix tie order is retained. Ensemble +restricted selection now uses the base predictor's deterministic tie order. + +The capture script accepts either historical successful output or a regional +input rejection for the formerly unsafe region examples. It remains an +observation tool; `tests/runOutputOverlap.sh` and the predictor API tests enforce +the revised behavior. The original controls and independent site oracles remain +assertions in the capture script. + +## Baseline findings + + +There are reproducible correctness bugs beyond the documented enumeration +heuristic. Region merging can violate the requested overlap rule, and some +heuristic predictors bypass output validation or use invalid energies for later +results. Separately, even exact prediction retains only one right extension per +left boundary for non-overlapping output. Consequently, a missing suboptimal +result does not establish that no compatible interaction exists. + +## Findings and proposed priority + +| Finding | Scope | Assessment | Proposed priority | +| --- | --- | --- | --- | +| Independent region results can overlap on a forbidden sequence | `N`, `T`, `Q`, with multiple region combinations | Correctness bug under the default output per sequence pair | High | +| Exhausted seed-free heuristic ensemble enumeration emits an invalid energy and repeats a site | `P/H --noSeed`, `N/T/Q`, requesting more results than available | Undefined numeric conversion | High | +| Later seed-free heuristic ensemble results omit accessibility contributions | `P/H --noSeed`, `N/T/Q` | Incorrect reported energy, ranking and energy filtering | High | +| Later heuristic MFE results can ignore `outNoGUend` | Confirmed for `S/H`, seeded and seed-free, with `T` | Output constraint violation | High | +| Region merging does not apply a global `outDeltaE` threshold | All four overlap modes, multiple regions | Same local versus global scope problem | Medium | +| Valid non-overlapping alternatives can be discarded before selection | `N/T/Q`, including `mode=M` | Documented heuristic with a substantial completeness limitation | Consult on semantics and cost | + +These priorities are recommendations, not implementation decisions. The +ensemble-specific findings also explain a single-region overlap violation; +regional merging is not the only way to obtain one. + +## What the four settings mean + +Overlap concerns inclusive interaction **intervals**, including unpaired bases +between the outermost intermolecular pairs. It does not mean shared base pairs +or only shared paired nucleotides. Coordinates below use the normal one-based +forward target and query indices. + +| Setting | Target overlap allowed | Query overlap allowed | Reason to reject a later interaction | +| --- | --- | --- | --- | +| `B` | Yes | Yes | Neither interval is excluded by this option | +| `T` | Yes | No | Its query interval intersects an already selected query interval | +| `Q` | No | Yes | Its target interval intersects an already selected target interval | +| `N` | No | No | Either interval intersects a selected interval on that sequence | + +The CLI mapping and `PredictorMfe::reportOptima()` implement this interpretation +correctly within their local selection state. The internally reversed query +coordinates are converted through the energy interface. There is no evidence +here of a simple swapped `T`/`Q` flag or reversed-query indexing error. + +The [README](../../README.md#subopts) promises **up to** `outNumber` results, +including the optimum. `B` does not enumerate every structure: exact prediction +normally retains the best structure for each boundary pair, and heuristic +predictors can consider fewer candidates. Allowing an overlap does not require +one. `N` does not search for the best *set* of simultaneously occupied sites; +the current algorithm chooses an optimum first and then compatible alternatives. + +## Reproduction setup + +Build IntaRNA using the [repository instructions](../../AGENTS.md), then run +from the repository root: + +```bash +python3 doc/analysis/out-overlap/reproduce.py ./src/bin/IntaRNA \ + --output /tmp/intarna-overlap-observations.json +``` + +The script captures the arguments, exit status, output, and forbidden interval +overlaps for 118 invocations. It checks 64 non-exhausting controls: all four +overlap settings for 16 CLI predictor configurations (including RIblast mode). +It independently enumerates every antiparallel structure for three tiny +seed-free examples and verifies their exact `B` site energies. Known defects +are recorded as observations, not installed as regression-test expectations. + +All commands below use a Bash helper to keep the examples short: + +```bash +rna() { + ./src/bin/IntaRNA --threads=1 --outMode=C \ + --outCsvCols=start1,end1,start2,end2,E "$@" +} +``` + +Except where stated, examples use `energy=B` (minus one kcal/mol per pair) and +`acc=N` (zero accessibility penalty) to make energies directly inspectable. These +are demonstrations of program behavior, not biological predictions. Explicit +small seeds and seed-free examples are intentional. `--outDeltaE=100` makes +the energy window non-limiting in these examples. + +## 1 Regional merging permits forbidden overlaps + +```bash +rna -t CCAACC -q GG --energy=B --acc=N --seedBP=2 \ + --tRegion=1-2,5-6 --outOverlap=N -n 10 --outDeltaE=100 +``` + +Actual output: + +```text +start1;end1;start2;end2;E +1;2;1;2;-2 +5;6;1;2;-2 +``` + +The query interval `1..2` is reused. With the default `outPerRegion=false`, +`N` should permit at most one of these two interactions. `T` also produces this +violation; `Q` and `B` permit these rows. Conversely: + +```bash +rna -t CC -q GGAAGG --energy=B --acc=N --seedBP=2 \ + --qRegion=1-2,5-6 --outOverlap=Q -n 10 --outDeltaE=100 +``` + +reports target `1..2` with query `1..2` and `5..6`, violating `Q` and, if selected, +`N`. The permitted cases are `T` and `B`. + +Automatic decomposition reaches the same path: + +```bash +rna -t UAUCGGCC -q GG --energy=B --acc=C --seedBP=2 \ + --tRegionLenMax=4 --outOverlap=N -n 10 --outDeltaE=100 +``` + +This reports `(3..4,1..2,-2)` and `(7..8,1..2,-2)`. With `-v`, the prediction +messages show the two target regions. The problem also occurs with ViennaRNA +energies: replace the first example's target/query with `CCCACCC`/`GGG`, use +`--tRegion=1-3,5-7`, and omit `--energy=B`. It reports `-4.2` and `-3` kcal/mol +while reusing query `1..3`. + +**Cause.** [IntaRNA.cpp](../../src/bin/IntaRNA.cpp) creates a fresh predictor +for every region/window combination. Each predictor clears its own +`reportedInteractions` lists. +[OutputHandlerInteractionList::add](../../src/IntaRNA/OutputHandlerInteractionList.cpp) +merges, sorts, deduplicates and truncates their results but never enforces +`reportOverlap` across combinations. Disjoint target regions do not prevent +their interactions from sharing a query region, and vice versa. + +**Proposed change.** For default output per sequence pair, apply overlap +selection in a common sequence-pair context using original sequence +coordinates. Preserve independent selection when users explicitly request +`outPerRegion=true`, and document that scope. The reproduction with +`outPerRegion=true -n 1` intentionally produces one result per region, so its +combined overlaps are not evidence of a violation of that independent scope. + +A final overlap filter alone can ensure validity, but cannot ensure the best +available compatible output or fill the requested count. A region's local +winner may suppress a local alternative and later be rejected against a better +winner from another region. That suppressed alternative can now be admissible. +Selection therefore needs unpruned candidate streams, or a way to request or +recompute alternatives after global exclusions. Applying a fixed top-k cap +before compatibility filtering has the same shortfall. + +Window mode is explicitly different: the CLI rejects `outNumber>1` with +`N/T/Q` when `windowWidth` is set. The reproduction script checks that guard. +Manual and automatic regions currently lack equivalent protection or global +selection. `--rri` evaluation intentionally ignores overlap and other prediction +filters and is outside this report's prediction contract. + +## 2 Exhausted ensemble enumeration converts infinity to an integer + +```bash +rna -t CC -q GG --energy=B --acc=N --model=P --mode=H --noSeed \ + --outOverlap=N -n 2 --outDeltaE=100 +``` + +Actual output on the checked GCC build: + +```text +start1;end1;start2;end2;E +1;2;1;2;-2.14748e+07 +1;2;1;2;-2 +``` + +Only the `-2` row is valid. The corrupted row is sorted ahead of the real +optimum by the collector. `T` and `Q` reproduce the problem; exact `mode=M` +returns only the valid row. `B` follows a different enumeration path and does +not encounter this exhaustion conversion. + +**Cause.** In +[PredictorMfeEns2dHeuristic::getNextBest](../../src/IntaRNA/PredictorMfeEns2dHeuristic.cpp), +`curBestCellE` has type `Z_type` and starts at `Z_INF`. With no eligible cell, +the function assigns floating-point infinity to the integer `curBest.energy`. +It leaves the previous coordinates intact. The invalid converted energy can +pass the next report-loop checks. This is undefined behavior; the exact printed +number and even the failure mode are not portable. + +A focused rebuild of this unmodified translation unit with GCC 14.4 and +`-fsanitize=undefined,float-cast-overflow -fno-sanitize-recover=all` confirmed: + +```text +PredictorMfeEns2dHeuristic.cpp:291:19: runtime error: +inf is outside the range of representable values of type 'int' +``` + +**Proposed change.** Represent candidate energies with `E_type` and its +`E_INF` sentinel, and explicitly return exhaustion before updating coordinates +or emitting another result. Test zero, one and several remaining candidates +for all restricted overlap modes, with both boundary-only and traced output. +Reusing the final validated candidate path described next would remove the +need for this independent raw-matrix enumeration. + +## 3 Later ensemble results use the wrong energy + +This is distinct from exhaustion; `-n 2` below stops while a real second +candidate exists. + +```bash +rna -t AGAGC -q GAUUC --energy=B --acc=C --model=P --mode=H --noSeed \ + --outOverlap=N -n 2 --outDeltaE=100 --outMaxE=0 +``` + +Actual output: + +```text +start1;end1;start2;end2;E +1;2;3;4;-2 +5;5;1;1;-1 +``` + +The second site has `ED1=0`, `ED2=1.31`, and total energy `0.31`. The same +predictor reports `0.31` for that site with `--outOverlap=B -n 100 +--outMaxE=100`; include `ED1,ED2` in `outCsvCols` to inspect the penalties. +There is only one pair in the site, so no alternative interior structure +explains this difference. It should be excluded by `outMaxE=0`. + +**Cause.** The heuristic ensemble `getNextBest()` uses +`energy.getE(curCell->val)`, the hybrid score from its recursion matrix. +The first-result path instead calls `updateOptimaUsingZ()` on finalized site +partition values and `updateOptima(..., isHybridE=true)`, which adds the site +energy contributions. Thus changing the overlap option changes the energy +definition for later rows. This also invalidates their ranking and their +comparison with `outDeltaE`/`outMaxE`. + +**Proposed change.** Enumerate later results from the same finalized, filtered +site energies used for the first result. Adding ED to a raw matrix score alone +would not establish equivalence with finalized per-site partition aggregation. +Test the energy of identical boundaries across overlap modes, including +nonzero ED and energy cutoffs. Seeded `P/H` uses a different predictor class; +this reproducer and the raw-matrix diagnosis concern `P/H --noSeed`. + +## 4 Later heuristic results bypass terminal pair filtering + +```bash +rna -t UUGA -q CAUU --energy=B --acc=N --model=S --mode=H --noSeed \ + --outNoGUend --outNoLP=false --outOverlap=T -n 10 --outDeltaE=100 +``` + +Actual output: + +```text +start1;end1;start2;end2;E +1;3;1;2;-2 +3;4;3;4;-2 +``` + +The second row pairs target G at position 3 with query U at position 4. This +terminal GU pair violates `outNoGUend`. The row is absent with `B` and with +exact `S/M`. Replacing `--noSeed` by `--seedBP=2 --seedMaxUP=2` reproduces the +same violation in seeded `S/H`. The corresponding seeded `X/H` control does +not reproduce it. + +**Cause.** [PredictorMfe::updateOptima](../../src/IntaRNA/PredictorMfe.cpp) +checks terminal GU and maximum ED before recording output candidates. The +specialized `getNextBest()` functions in +[PredictorMfe2dHeuristic](../../src/IntaRNA/PredictorMfe2dHeuristic.cpp) and +[PredictorMfe2dHeuristicSeed](../../src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp) +read recursion cells directly and bypass that validation. The raw state can +legitimately contain intermediate extensions that are unsuitable as complete +reported interactions. + +**Proposed change.** Share complete-candidate validation and the final energy +calculation between initial and subsequent selection. Audit the analogous +helix and ensemble overrides for all output constraints, including maximum ED; +the GU reproducer above establishes the two `S/H` paths, not every possible +constraint failure in every override. Merely rejecting the bad final row may +still hide a valid alternative at its starting position, so candidate retention +must also be considered. + +## 5 Energy windows are local to regions + +```bash +rna -t CCAAC -q GG --energy=B --acc=N --model=S --noSeed \ + --tRegion=1-2,5-5 --outOverlap=B -n 10 --outDeltaE=0 +``` + +Actual output: + +```text +start1;end1;start2;end2;E +1;2;1;2;-2 +5;5;1;1;-1 +5;5;2;2;-1 +``` + +The pair's best energy is `-2`. The `-1` rows are outside a zero-width energy +window around it. They survive because each predictor uses its own local MFE +as the reference and the collector does not apply `deltaE` again. The script +reproduces this in all four overlap settings; restricted modes change how many +of the local rows survive, not the scope of the threshold. + +**Proposed change.** Under default output per sequence pair, evaluate the +energy window relative to that pair's best eligible candidate when combining +regions. Under explicit output per region, retain and document local windows. +This is coupled to the global selector in finding 1; `B` needs the global +energy-window check even though it needs no overlap rejection. + +## 6 A valid right extension can be discarded before selection + +```bash +rna -t CCCAAU -q GGGUU --energy=B --acc=N --seedBP=2 \ + --model=X --mode=M --outOverlap=N -n 100 --outDeltaE=100 +``` + +Only `(1..3,1..3,-3)` is returned. Yet `(4..5,4..5,-2)`, a two-pair AU seed, +is disjoint on both sequences. It appears with `--outOverlap=B`. It is also +returned by rerunning this particular example with +`--tAccConstr=b:1-3 --qAccConstr=b:1-3`. `mode=H` shows the same omission; +changing `N` to `T` does too. + +**Cause.** [PredictorMfe::updateMfe4leftEnd](../../src/IntaRNA/PredictorMfe.cpp) +stores just one best right end for each `(i1,i2)`. At target start 4 and forward +query end 5, it retains the better `(4..6,1..5,-3)` extension. That extension +overlaps the first result on the query. The shorter `(4..5,4..5,-2)` extension +has already been discarded, so `getNextBest()` cannot recover it. Several +heuristic predictors instead scan matrices with the same one-extension limit. + +The following still smaller seed-free examples establish the limitation for +each restricted mode. The script verifies their complete exact `B` site +energies against an independent structure enumerator. Use `model=S`, `mode=M`, +`noSeed`, `outNoLP=false`, `energy=B`, `acc=N`, `outDeltaE=100`, and `-n 1000`. + +| Mode | Target | Query | Only reported site and energy | Omitted compatible site and energy | +| --- | --- | --- | --- | --- | +| `N` | `CCAU` | `CGGU` | `(1..2,2..3,-2)` | `(3..3,4..4,-1)` | +| `T` | `GUA` | `GUA` | `(1..2,1..2,-2)` | `(2..2,3..3,-1)` | +| `Q` | `GCGC` | `AGUG` | `(2..4,2..4,-3)` | `(1..1,3..3,-1)` | + +This behavior is already acknowledged in the README's warning about +non-overlapping output. It is therefore a documented completeness limitation, +not evidence that every missing suboptimal row is a newly introduced bug. +Calling `mode=M` exact can nevertheless mislead users unless the distinction +between exact candidate calculation and heuristic restricted enumeration is +explicit. These examples are small witnesses; no global minimality claim is +made. + +**Proposed choices for consultation.** Keep and clearly document the heuristic, +or support complete greedy selection by retaining more endpoints or recomputing +the best candidate after excluding previously selected intervals. For `N`, +recomputation must cover all combinations of remaining target and query +intervals; for `T` or `Q`, only the forbidden-overlap sequence is split. Avoid +changing the underlying accessibility model when excluding interaction spans. +Retaining all boundaries can require quartic storage, whereas repeated +prediction trades memory for time. Do not promise exact completeness for a +heuristic predictor merely because output selection has improved. + +Blocking paired nucleotides is not a general replacement for interval exclusion. +For example, `-t CCGG -q CCGG --model=S --mode=M --noSeed --energy=B --acc=N +--tAccConstr=b:2-3 --outOverlap=N -n 2 --outDeltaE=100` still produces target +interval `1..4`, with the blocked positions inside a loop. A span-based selector +must reject crossings over an excluded interval, not just forbid pairs at its +positions. + +## Original recommendations before review + +The recommended contract is: choose the best eligible interaction first, then +greedily choose the best eligible interaction compatible with every previously +selected interval, within a documented sequence-pair or region scope. This +preserves the current interpretation of ranked alternatives. Optimizing the +number of interactions or the total energy of a compatible set is a different +objective that can discard the single best interaction; it should not be +introduced implicitly as an overlap bug fix. + +Agree on that scope, the desired completeness/runtime tradeoff, and tie-breaking +before changing enumeration. `Interaction::operator<`, the base map scan, and +specialized reverse matrix scans currently use different tie rules. Equal-energy +choices can change which later intervals remain available. This report does +not classify that unspecified ordering as a separate bug, but a common selector +needs an explicit rule in original sequence coordinates. + +Suggested implementation order: + +1. Repair the invalid ensemble sentinel and use common validated final energies + for all emitted candidates. Add exhaustion, ED and terminal-pair regressions. +2. Enforce sequence-pair overlap and energy-window scope across regions, with + explicit independent behavior for `outPerRegion`. Include alternatives that + become available when a local winner is globally rejected. +3. Decide whether to retain the documented one-extension heuristic or introduce + more complete restricted enumeration, with measured memory/runtime costs. + +Future regressions should check interval validity, energy identity, output +constraints, maximum count and completeness separately. Cover `B/N/T/Q`, +multiple regions on either sequence, automatic decomposition, nonzero region +offsets, shifted output indices, tied energies, seeds, terminal GU constraints, +nonzero accessibility and exhaustion. Existing overlap fixtures cap interaction +lengths at two bases, so they cannot detect a discarded longer/shorter extension +like finding 6. The independent toy enumerator is a useful oracle for those +regressions; the known-bug observations themselves are not correct expected +output. + +## Baseline validation and limits + +The analysis used a fresh release build of the stated master revision with +GCC 14.4, bundled Kokkos mdspan, ViennaRNA 2.7.2 and Boost 1.85.0 on Linux. +The 118-run script passed all 64 control checks and three independent site +oracles. All 20 existing CLI fixtures passed via `tests/runIntaRNA.sh`, run +from `tests/` with `INTARNABINPATH` set to the repository root. The full API +suite was not rerun for this analysis-only change. Documentation distribution +was checked after regenerating the Autotools files. The focused sanitizer +probe instrumented the ensemble heuristic +translation unit, not the entire application. No production source, CLI +behavior, or existing test expectations were modified. + +These controls establish the reported examples and shared code paths, not +exhaustive correctness of every predictor, energy model or constraint +combination. In particular, the successful `B` controls do not prove exhaustive +structure enumeration, and the platform-specific numeric manifestation of +undefined behavior must not be used as a portable regression expectation. diff --git a/doc/analysis/out-overlap/reproduce.py b/doc/analysis/out-overlap/reproduce.py new file mode 100644 index 00000000..47eb693c --- /dev/null +++ b/doc/analysis/out-overlap/reproduce.py @@ -0,0 +1,198 @@ +#!/usr/bin/env python3 +"""Capture issue 212 observations without changing IntaRNA. + +Historical defects and their repairs are recorded, not blessed as expectations. +Regional examples accept either historical output or the revised input rejection. +Only the small independent oracle and the non-exhausting controls are asserted. +Requires Python 3 and an already built IntaRNA executable. +""" + +import argparse +import csv +import io +import itertools +import json +import pathlib +import subprocess + + +COLS = "start1,end1,start2,end2,E,ED1,ED2" +COMMON = ["--threads=1", "--outMode=C", "--outCsvCols=" + COLS] +TOY = ["--energy=B", "--acc=N", "--outDeltaE=100"] +FORBIDDEN = {"N": "12", "T": "2", "Q": "1", "B": ""} + + +def overlaps(a, b, axis): + return max(a["start" + axis], b["start" + axis]) <= min( + a["end" + axis], b["end" + axis]) + + +def conflicts(rows, mode): + return [ + [i + 1, j + 1, axis] + for (i, a), (j, b) in itertools.combinations(enumerate(rows), 2) + for axis in FORBIDDEN[mode] if overlaps(a, b, axis) + ] + + +def site(row): + return tuple(row[key] for key in ("start1", "end1", "start2", "end2")) + + +def base_pair_oracle(target, query): + """Enumerate every antiparallel structure for the four-base toy example. + + No accessibility, seed or lonely-pair restrictions apply. At this length + every interior loop fits the default loop limit. Each pair contributes -1. + Keep the best energy for each pair of inclusive sequence intervals. + """ + best = {} + pairs = [(i, j) for i, t in enumerate(target, 1) + for j, q in enumerate(query, 1) if t + q in + {"AU", "UA", "CG", "GC", "GU", "UG"}] + + def extend(structure): + i, j = structure[-1] + bounds = (structure[0][0], i, j, structure[0][1]) + best[bounds] = min(best.get(bounds, 0), -len(structure)) + for ni, nj in pairs: + if ni > i and nj < j: + extend(structure + [(ni, nj)]) + + for pair in pairs: + extend([pair]) + return best + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("binary", type=pathlib.Path) + parser.add_argument("--output", type=pathlib.Path, required=True) + args = parser.parse_args() + binary = str(args.binary.resolve()) + observations = {} + + def run(name, options, mode="B", expected_status=0): + command = COMMON + options + ["--outOverlap=" + mode] + process = subprocess.run([binary] + command, capture_output=True, + text=True, timeout=30) + if expected_status is not None and (process.returncode == 0) != (expected_status == 0): + raise RuntimeError(f"{name}: unexpected exit {process.returncode}\n" + + process.stdout + process.stderr) + lines = [line for line in process.stdout.splitlines() + if line and not line.startswith("#")] + rows = [] + if process.returncode == 0: + reader = csv.DictReader(io.StringIO("\n".join(lines)), delimiter=";") + if reader.fieldnames != COLS.split(","): + raise RuntimeError(f"{name}: unexpected CSV header {reader.fieldnames}") + for row in reader: + rows.append({key: (int(value) if key.startswith(("start", "end")) + else float(value)) for key, value in row.items()}) + observations[name] = { + "arguments": command, "returncode": process.returncode, + "stdout": process.stdout, "stderr": process.stderr, + "rows": rows, "forbidden_overlaps": conflicts(rows, mode), + } + return rows + + # Every CLI predictor combination, with enough disjoint two-pair sites to + # avoid the separately investigated ensemble exhaustion defect. + configurations = [] + for model, modes, seeded in [ + ("S", "HM", False), ("S", "HMS", True), ("X", "HMRS", True), + ("P", "HM", False), ("P", "HMS", True), + ("B", "H", False), ("B", "H", True), + ]: + for mode in modes: + configurations.append((model, mode, seeded)) + for model, mode, seeded in configurations: + for overlap in "BNTQ": + name = f"control_{model}_{mode}_{'seed' if seeded else 'noSeed'}_{overlap}" + options = TOY + ["-t", "CCAACC", "-q", "GGAAGG", "-n", "2", + "--intLenMax=2", "--model=" + model, "--mode=" + mode, + "--helixMinBP=2", "--helixMaxBP=2", + "--seedBP=2" if seeded else "--noSeed"] + rows = run(name, options, overlap) + assert len(rows) == 2, (name, rows) + assert all(row["E"] == -2 for row in rows), (name, rows) + assert not conflicts(rows, overlap), (name, rows) + + for overlap in "BNTQ": + run("target_regions_" + overlap, TOY + ["-t", "CCAACC", "-q", "GG", + "--seedBP=2", "--tRegion=1-2,5-6", "-n", "10"], overlap, + expected_status=None if overlap in "NT" else 0) + run("query_regions_" + overlap, TOY + ["-t", "CC", "-q", "GGAAGG", + "--seedBP=2", "--qRegion=1-2,5-6", "-n", "10"], overlap, + expected_status=None if overlap in "NQ" else 0) + run("regional_delta_" + overlap, ["--energy=B", "--acc=N", "-t", "CCAAC", + "-q", "GG", "--noSeed", "--model=S", "--tRegion=1-2,5-5", + "--outDeltaE=0", "-n", "10"], overlap, + expected_status=None if overlap in "NT" else 0) + run("window_rejection_" + overlap, TOY + ["-t", "CCAACCAACCAA", "-q", "GGAAGGAAGGAA", + "--seedBP=2", "--intLenMax=2", "--windowWidth=10", "--windowOverlap=2", + "-n", "2"], overlap, expected_status=0 if overlap == "B" else 1) + + run("automatic_regions", ["--energy=B", "--acc=C", "-t", "UAUCGGCC", "-q", "GG", + "--seedBP=2", "--tRegionLenMax=4", "--outDeltaE=100", "-n", "10"], "N", expected_status=None) + run("regions_vienna", ["--acc=N", "-t", "CCCACCC", "-q", "GGG", "--seedBP=2", + "--tRegion=1-3,5-7", "--outDeltaE=100", "-n", "10"], "N", expected_status=None) + run("explicit_per_region", TOY + ["-t", "CCAACC", "-q", "GG", "--seedBP=2", + "--tRegion=1-2,5-6", "--outPerRegion", "-n", "1"], "N") + + seeded_missing = TOY + ["-t", "CCCAAU", "-q", "GGGUU", "--seedBP=2", "-n", "100"] + for mode in "HM": + for overlap in "BNTQ": + run(f"missing_seed_{mode}_{overlap}", seeded_missing + ["--mode=" + mode], overlap) + run("blocked_seed_" + mode, seeded_missing + ["--mode=" + mode, + "--tAccConstr=b:1-3", "--qAccConstr=b:1-3"], "N") + + for overlap, target, query in [("N", "CCAU", "CGGU"), ("T", "GUA", "GUA"), + ("Q", "GCGC", "AGUG")]: + options = TOY + ["-t", target, "-q", query, "--model=S", "--mode=M", + "--noSeed", "--outNoLP=false", "-n", "1000"] + all_rows = run("oracle_B_" + overlap, options) + oracle = base_pair_oracle(target, query) + assert {site(row): row["E"] for row in all_rows} == oracle + assert len(all_rows) == len(oracle) + selected = run("oracle_" + overlap, options, overlap) + observations["oracle_" + overlap]["omitted_compatible_sites"] = [ + row for row in all_rows if all(not overlaps(row, old, axis) + for old in selected for axis in FORBIDDEN[overlap])] + + for model, mode, seed in [("S", "H", ["--noSeed"]), + ("S", "M", ["--noSeed"]), + ("S", "H", ["--seedBP=2", "--seedMaxUP=2"]), + ("X", "H", ["--seedBP=2", "--seedMaxUP=2"])]: + for overlap in "BT": + name = f"gu_filter_{model}_{mode}_{seed[0]}_{overlap}" + rows = run(name, TOY + ["-t", "UUGA", "-q", "CAUU", "--model=" + model, + "--mode=" + mode, "--outNoGUend", "--outNoLP=false", "-n", "10"] + seed, overlap) + observations[name]["gu_ended_rows"] = [i + 1 for i, row in enumerate(rows) + if any("UUGA"[row[t] - 1] + "CAUU"[row[q] - 1] in ("GU", "UG") + for t, q in [("start1", "end2"), ("end1", "start2")])] + + for mode in "HM": + for overlap in "BNTQ": + run(f"ensemble_exhaustion_{mode}_{overlap}", TOY + ["-t", "CC", "-q", "GG", + "--model=P", "--mode=" + mode, "--noSeed", "-n", "2"], overlap) + for overlap in "BN": + run("ensemble_accessibility_" + overlap, ["--energy=B", "--acc=C", "-t", "AGAGC", + "-q", "GAUUC", "--model=P", "--mode=H", "--noSeed", "--outDeltaE=100", + "-n", "100" if overlap == "B" else "2", "--outMaxE=100" if overlap == "B" + else "--outMaxE=0"], overlap) + run("blocking_can_bridge", TOY + ["-t", "CCGG", "-q", "CCGG", "--noSeed", + "--model=S", "--mode=M", "--tAccConstr=b:2-3", "-n", "2"], "N") + + data = {"version": subprocess.check_output([binary, "--version"], text=True).strip(), + "control_configurations": len(configurations), + "control_runs": len(configurations) * 4, + "independent_oracles": 3, "invocations": len(observations), + "observations": observations} + args.output.write_text(json.dumps(data, indent=2) + "\n", encoding="utf-8") + print(f"Captured {len(observations)} runs; {len(configurations) * 4} controls and " + f"3 independent site oracles passed. Observations: {args.output}") + + +if __name__ == "__main__": + main() diff --git a/src/IntaRNA/PredictorMfe.cpp b/src/IntaRNA/PredictorMfe.cpp index 9cf7277f..91792bea 100644 --- a/src/IntaRNA/PredictorMfe.cpp +++ b/src/IntaRNA/PredictorMfe.cpp @@ -73,6 +73,22 @@ initOptima() //////////////////////////////////////////////////////////////////////////// +bool +PredictorMfe:: +isValidOutputSite( const size_t i1, const size_t j1, + const size_t i2, const size_t j2 ) const +{ + const OutputConstraint & constraint = output.getOutputConstraint(); + if (constraint.noGUend && (energy.isGU(i1,i2) || energy.isGU(j1,j2))) { + return false; + } + return constraint.maxED >= Accessibility::ED_UPPER_BOUND + || (energy.getED1(i1,j1) <= constraint.maxED + && energy.getED2(i2,j2) <= constraint.maxED); +} + +//////////////////////////////////////////////////////////////////////////// + void PredictorMfe:: updateOptima( const size_t i1, const size_t j1 @@ -87,17 +103,7 @@ updateOptima( const size_t i1, const size_t j1 return; } - // check GU ends if needed - if (output.getOutputConstraint().noGUend && (energy.isGU(i1,i2) || energy.isGU(j1,j2)) ) { - return; - } - - // check ED penalties - if (output.getOutputConstraint().maxED < Accessibility::ED_UPPER_BOUND - && (energy.getED1(i1,j1) > output.getOutputConstraint().maxED - || energy.getED2(i2,j2) > output.getOutputConstraint().maxED) - ) - { + if (!isValidOutputSite(i1, j1, i2, j2)) { return; } @@ -228,6 +234,11 @@ reportOptima() && (curBest.energy < mfeDeltaE || E_equal(curBest.energy,mfeDeltaE)) && reported < outConstraint.reportMax ) { + // Selectors return two boundaries, which coincide for a single pair. + // Normalize before either traceback or boundary-only output validation. + if (curBest.basePairs.size() == 2 && curBest.basePairs.front() == curBest.basePairs.back()) { + curBest.basePairs.resize(1); + } // report current best if (outConstraint.needBPs) { // fill interaction with according base pairs diff --git a/src/IntaRNA/PredictorMfe.h b/src/IntaRNA/PredictorMfe.h index 239ce7a9..35501d33 100644 --- a/src/IntaRNA/PredictorMfe.h +++ b/src/IntaRNA/PredictorMfe.h @@ -114,6 +114,19 @@ class PredictorMfe : public Predictor { void initOptima(); + /** + * Checks terminal GU and accessibility constraints for a complete site. + * Recursion cells may be useful extensions without being valid output sites. + * @param i1 inclusive start in sequence 1, relative to the energy offset + * @param j1 inclusive end in sequence 1, relative to the energy offset + * @param i2 inclusive start in reversed sequence 2, relative to its offset + * @param j2 inclusive end in reversed sequence 2, relative to its offset + * @return whether both terminal pairs and ED penalties satisfy output limits + */ + bool + isValidOutputSite( const size_t i1, const size_t j1, + const size_t i2, const size_t j2 ) const; + /** * updates the global optimum to be the mfe interaction if needed * diff --git a/src/IntaRNA/PredictorMfe2dHelixBlockHeuristic.cpp b/src/IntaRNA/PredictorMfe2dHelixBlockHeuristic.cpp index e7bafe7d..537719b6 100644 --- a/src/IntaRNA/PredictorMfe2dHelixBlockHeuristic.cpp +++ b/src/IntaRNA/PredictorMfe2dHelixBlockHeuristic.cpp @@ -375,7 +375,8 @@ getNextBest( Interaction & curBest ) // direct cell access curCell = &(hybridE(i1,i2)); // check if left side can pair - if (E_isINF(curCell->val)) + if (E_isINF(curCell->val) + || !isValidOutputSite(i1,curCell->j1,i2,curCell->j2)) { continue; } diff --git a/src/IntaRNA/PredictorMfe2dHelixBlockHeuristicSeed.cpp b/src/IntaRNA/PredictorMfe2dHelixBlockHeuristicSeed.cpp index bfa37c65..aaa6cef8 100644 --- a/src/IntaRNA/PredictorMfe2dHelixBlockHeuristicSeed.cpp +++ b/src/IntaRNA/PredictorMfe2dHelixBlockHeuristicSeed.cpp @@ -503,7 +503,8 @@ getNextBest( Interaction & curBest ) // direct cell access curCell = &(hybridE_seed(i1,i2)); // check if left side can pair - if (E_isINF(curCell->val)) + if (E_isINF(curCell->val) + || !isValidOutputSite(i1,curCell->j1,i2,curCell->j2)) { continue; } diff --git a/src/IntaRNA/PredictorMfe2dHeuristic.cpp b/src/IntaRNA/PredictorMfe2dHeuristic.cpp index c9398be8..06983318 100644 --- a/src/IntaRNA/PredictorMfe2dHeuristic.cpp +++ b/src/IntaRNA/PredictorMfe2dHeuristic.cpp @@ -378,7 +378,8 @@ getNextBest( Interaction & curBest ) // direct cell access curCell = &(hybridE(i1,i2)); // check if left side can pair - if (E_isINF(curCell->val)) + if (E_isINF(curCell->val) + || !isValidOutputSite(i1,curCell->j1,i2,curCell->j2)) { continue; } diff --git a/src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp b/src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp index 822bb9ba..2c7a29f5 100644 --- a/src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp +++ b/src/IntaRNA/PredictorMfe2dHeuristicSeed.cpp @@ -650,7 +650,8 @@ getNextBest( Interaction & curBest ) // direct cell access curCell = &(hybridE_seed(i1,i2)); // check if left side can pair - if (E_isINF(curCell->val)) + if (E_isINF(curCell->val) + || !isValidOutputSite(i1,curCell->j1,i2,curCell->j2)) { continue; } diff --git a/src/IntaRNA/PredictorMfeEns.cpp b/src/IntaRNA/PredictorMfeEns.cpp index e7391ae1..08b2e2d7 100644 --- a/src/IntaRNA/PredictorMfeEns.cpp +++ b/src/IntaRNA/PredictorMfeEns.cpp @@ -51,16 +51,7 @@ addPartitionContribution( const size_t i1, const size_t j1 // 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 false; - } - if (outConstraint.maxED < Accessibility::ED_UPPER_BOUND - && (energy.getED1(i1,j1) > outConstraint.maxED - || energy.getED2(i2,j2) > outConstraint.maxED)) - { + if (!isValidOutputSite(i1, j1, i2, j2)) { return false; } diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp index 49cdb37f..806ea98a 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp @@ -224,89 +224,4 @@ fillHybridZ() //////////////////////////////////////////////////////////////////////////// -void -PredictorMfeEns2dHeuristic:: -getNextBest( Interaction & curBest ) -{ - - // get original - const Z_type curBestE = curBest.energy; - - // identify cell with next best non-overlapping interaction site - // iterate (decreasingly) over all left interaction starts - size_t i1,i2; - BestInteractionZ * curBestCell = NULL; - Z_type curBestCellE = Z_INF; - Interaction::BasePair curBestCellStart; - BestInteractionZ * curCell = NULL; - Z_type curCellE = Z_INF; - IndexRange r1,r2; - for (i1=hybridZ.size1(); i1-- > 0;) { - // ensure interaction site start is not covered - if (reportedInteractions.first.covers(i1)) { - continue; - } - for (i2=hybridZ.size2(); i2-- > 0;) { - // ensure interaction site start is not covered - if (reportedInteractions.second.covers(i2)) { - continue; - } - // direct cell access - curCell = &(hybridZ(i1,i2)); - // check if left side can pair - if (Z_isINF(curCell->val)) - { - continue; - } - // get overall energy of the interaction - curCellE = energy.getE(curCell->val); - // or energy is too low to be considered - // or energy is higher than current best found so far - if (curCellE < curBestE || curCellE >= curBestCellE ) - { - continue; - } - // ensure site is not overlapping - r1.from = i1; - r1.to = curCell->j1; - if ( reportedInteractions.first.overlaps( r1 )) { - continue; - } - r2.from = i2; - r2.to = curCell->j2; - if ( reportedInteractions.second.overlaps( r2 )) { - continue; - } - //// FOUND THE NEXT BETTER SOLUTION - // overwrite current best found so far - curBestCell = curCell; - curBestCellE = curCellE; - curBestCellStart.first = i1; - curBestCellStart.second = i2; - - } // i2 - } // i1 - - // overwrite curBest - curBest.energy = curBestCellE; - curBest.basePairs.resize(2); - if (E_isNotINF(curBestCellE)) { - curBest.basePairs[0] = energy.getBasePair( curBestCellStart.first, curBestCellStart.second ); - curBest.basePairs[1] = energy.getBasePair( curBestCell->j1, curBestCell->j2 ); - } -} - -//////////////////////////////////////////////////////////////////////////// - -void -PredictorMfeEns2dHeuristic:: -updateMfe4leftEnd(const size_t i1, const size_t j1 - , const size_t i2, const size_t j2 - , const Interaction & curInteraction ) -{ - // do nothing since getNextBest() is based on local data structure -} - -//////////////////////////////////////////////////////////////////////////// - } // namespace diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristic.h b/src/IntaRNA/PredictorMfeEns2dHeuristic.h index cec1c905..8167c1c5 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristic.h +++ b/src/IntaRNA/PredictorMfeEns2dHeuristic.h @@ -85,35 +85,8 @@ class PredictorMfeEns2dHeuristic: public PredictorMfeEns2d { void fillHybridZ(); - /** - * Identifies the next best interaction (containing a seed) - * with an energy equal to or higher - * than the given interaction. The new interaction will not overlap any - * index range stored in reportedInteractions. - * - * @param curBest IN/OUT the current best interaction to be replaced with one - * of equal or higher energy not overlapping with any reported - * interaction so far; an interaction with energy E_INF is set, if - * there is no better interaction left - */ - virtual - void - getNextBest( Interaction & curBest ); - - /** - * Overwrites function of super class to surpress the update. - * - * @param i1 interaction start in seq1 - * @param j1 interaction end in seq1 - * @param i2 interaction start in seq2 - * @param i2 interaction end in seq2 - * @param curInteraction the interaction information to be used for update - */ - virtual - void - updateMfe4leftEnd(const size_t i1, const size_t j1 - , const size_t i2, const size_t j2 - , const Interaction & curInteraction ); + // Restricted output uses PredictorMfe's best finalized site per left boundary. + // Intermediate hybridZ cells do not contain complete site ensemble energies. }; diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 583b61e3..4478868a 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -972,7 +972,9 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) "\n 'N' in none of the sequences, " "\n 'T' in the target only, " "\n 'Q' in the query only, " - "\n 'B' in both sequences").c_str()) + "\n 'B' in both sequences. With merged regions, N requires one region" + " per sequence, T one target region, Q one query region;" + " use --outPerRegion for independent region combinations.").c_str()) ; opts_cmdline_short.add(opts_output); opts_output.add_options() @@ -992,7 +994,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , value(&(outDeltaE.val)) ->default_value(outDeltaE.def) ->notifier(boost::bind(&CommandLineParsing::validate_numberArgument,this,outDeltaE,_1)) - , std::string("suboptimal output : only interactions with E <= (minE+deltaE) are reported" + , std::string("suboptimal output : only interactions with E <= (minE+deltaE) are reported;" + " minE is per sequence pair, or per region combination with --outPerRegion" " (arg in range ["+toString(outDeltaE.min)+","+toString(outDeltaE.max)+"])").c_str()) ("outBestSeedOnly" , value(&(outBestSeedOnly)) @@ -1370,6 +1373,9 @@ parse(int argc, char** argv) // parse region string if available parseRegion( "qRegion", qRegionString, query, qRegion ); parseRegion( "tRegion", tRegionString, target, tRegion ); + for (const auto & ranges : qRegion) { validateRegionOverlap(ranges, true); } + for (const auto & ranges : tRegion) { validateRegionOverlap(ranges, false); } + //////////////// ACCESSIBILITY CONSTRAINTS /////////////////// @@ -1911,6 +1917,26 @@ parseRegion( const std::string & argName, const std::string & value, const RnaSe //////////////////////////////////////////////////////////////////////////// +void +CommandLineParsing:: +validateRegionOverlap( const IndexRangeList & ranges, const bool isQuery ) const +{ + if (isEvaluation() || outPerRegion || ranges.size() <= 1 || outOverlap.val == 'B') { + return; + } + // Splitting one sequence allows separate predictions to reuse the other. + // Thus multiple target regions require query overlap, and vice versa. + if (outOverlap.val == 'N' || outOverlap.val == (isQuery ? 'Q' : 'T')) { + throw boost::program_options::error( + std::string("multiple ") + (isQuery ? "query" : "target") + + " regions cannot be merged with --outOverlap=" + outOverlap.val + + "; use --outPerRegion=true, --outOverlap=B, or a single " + + (isQuery ? "query" : "target") + " region"); + } +} + +//////////////////////////////////////////////////////////////////////////// + const CommandLineParsing::RnaSequenceVec & CommandLineParsing:: getQuerySequences() const @@ -2783,6 +2809,7 @@ getQueryRanges( const InteractionEnergy & energy, const size_t sequenceNumber, c } + validateRegionOverlap(qRegion.at(sequenceNumber), true); return qRegion.at(sequenceNumber); } @@ -2827,6 +2854,7 @@ getTargetRanges( const InteractionEnergy & energy, const size_t sequenceNumber, } } + validateRegionOverlap(tRegion.at(sequenceNumber), false); return tRegion.at(sequenceNumber); } diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index 3731406b..0deec81e 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -1129,6 +1129,16 @@ class CommandLineParsing { , const RnaSequenceVec & sequences , IndexRangeListVec & rangeList ); + /** + * Rejects region merging that cannot enforce the requested overlap rule. + * Applies to explicit regions and to the result of automatic decomposition. + * @param ranges non-overlapping prediction regions for one sequence + * @param isQuery whether these are query (true) or target (false) regions + * @throws boost::program_options::error for unsupported merged output + */ + void + validateRegionOverlap( const IndexRangeList & ranges, const bool isQuery ) const; + /** * Checks whether or not any command line argument were parsed. Throws a * std::runtime_error if not. diff --git a/src/bin/IntaRNA.cpp b/src/bin/IntaRNA.cpp index b7bdbd46..555e0d02 100644 --- a/src/bin/IntaRNA.cpp +++ b/src/bin/IntaRNA.cpp @@ -336,8 +336,14 @@ int main(int argc, char **argv){ {// update final output handler // copy partition function information if available output->incrementZ( bestInteractions.getZ() ); - // forward all reported interactions for all regions to final output handler + // Apply the energy window to the sequence pair's best candidate. + // Independent per-region output retains each region's local window. + const E_type maxMergedE = parameters.reportBestPerRegion() || bestInteractions.empty() + ? E_INF + : (*bestInteractions.begin())->energy + bestInteractions.getOutputConstraint().deltaE; + // The collector is sorted by energy, so all later entries are worse. for( const Interaction * inter : bestInteractions) { + if (inter->energy > maxMergedE) { break; } output->add(*inter); } } diff --git a/tests/Makefile.am b/tests/Makefile.am index 571deef5..6f580ca6 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 runOutputOverlap.sh # the program to build check_PROGRAMS = runApiTests diff --git a/tests/PredictorMfeEnsRegression_test.cpp b/tests/PredictorMfeEnsRegression_test.cpp index 76b0d3fb..faf71ec2 100644 --- a/tests/PredictorMfeEnsRegression_test.cpp +++ b/tests/PredictorMfeEnsRegression_test.cpp @@ -3,6 +3,7 @@ #undef NDEBUG #include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/AccessibilityBasePair.h" #include "IntaRNA/InteractionEnergyBasePair.h" #include "IntaRNA/OutputHandlerInteractionList.h" #include "IntaRNA/PredictorMfeEns2d.h" @@ -14,6 +15,7 @@ #include "IntaRNA/SeedHandlerNoBulge.h" #include +#include #include using namespace IntaRNA; @@ -328,3 +330,87 @@ TEST_CASE("ensemble predictor regressions", "[PredictorMfeEns]") { REQUIRE(predictor.getPartitionCount() == 1); } } + +TEST_CASE("heuristic ensemble suboptimals use finalized site energies", "[PredictorMfeEns][Overlap]") { + #include "testEasyLoggingSetup.icc" + + RnaSequence target("target", "AGAGC"); + RnaSequence query("query", "GAUUC"); + AccessibilityBasePair targetAcc(target, 0, nullptr); + AccessibilityBasePair queryAcc(query, 0, nullptr); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc); + + // B records the finalized boundary partitions, including nonzero ED. + OutputConstraint referenceConstraint(100, OutputConstraint::OVERLAP_BOTH, + Ekcal_2_E(100), Ekcal_2_E(100)); + OutputHandlerInteractionList reference(referenceConstraint, 100); + PredictorMfeEns2dHeuristic referencePredictor(energy, reference, nullptr); + referencePredictor.predict(); + REQUIRE_FALSE(reference.empty()); + + for (auto overlap : {OutputConstraint::OVERLAP_NONE, + OutputConstraint::OVERLAP_SEQ1, OutputConstraint::OVERLAP_SEQ2}) { + for (bool trace : {false, true}) { + for (E_type maxE : {E_type(0), Ekcal_2_E(100)}) { + for (E_type deltaE : {E_type(0), Ekcal_2_E(100)}) { + CAPTURE(overlap, trace, maxE, deltaE); + OutputConstraint constraint(100, overlap, maxE, deltaE, + false, false, false, false, trace); + OutputHandlerInteractionList output(constraint, 100); + PredictorMfeEns2dHeuristic predictor(energy, output, nullptr); + predictor.predict(); + REQUIRE_FALSE(output.empty()); + for (const Interaction * interaction : output) { + REQUIRE(interaction->isValid()); + const Interaction * sameSite = nullptr; + for (const Interaction * site : reference) { + if (site->basePairs.front() == interaction->basePairs.front() + && site->basePairs.back() == interaction->basePairs.back()) { + sameSite = site; + break; + } + } + REQUIRE(sameSite != nullptr); + REQUIRE(interaction->energy == sameSite->energy); + REQUIRE(interaction->energy < maxE); + REQUIRE(interaction->energy <= (*reference.begin())->energy + deltaE); + } + if (overlap == OutputConstraint::OVERLAP_NONE) { + REQUIRE(std::distance(output.begin(), output.end()) == (maxE == 0 || deltaE == 0 ? 1 : 2)); + } + } + } + } + } +} + +TEST_CASE("heuristic ensemble exhaustion stops without repeating sites", "[PredictorMfeEns][Overlap]") { + #include "testEasyLoggingSetup.icc" + + for (auto overlap : {OutputConstraint::OVERLAP_NONE, + OutputConstraint::OVERLAP_SEQ1, OutputConstraint::OVERLAP_SEQ2}) { + for (bool trace : {false, true}) { + for (bool canPair : {false, true}) { + CAPTURE(overlap, trace, canPair); + RnaSequence target("target", "CC"); + RnaSequence query("query", canPair ? "GG" : "CC"); + AccessibilityDisabled targetAcc(target, 0, nullptr); + AccessibilityDisabled queryAcc(query, 0, nullptr); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc); + OutputConstraint constraint(10, overlap, 0, Ekcal_2_E(100), + false, false, false, false, trace); + OutputHandlerInteractionList output(constraint, 10); + PredictorMfeEns2dHeuristic predictor(energy, output, nullptr); + predictor.predict(); + REQUIRE(std::distance(output.begin(), output.end()) == (canPair ? 1 : 0)); + if (canPair) { + REQUIRE((*output.begin())->energy == Ekcal_2_E(-2)); + REQUIRE((*output.begin())->basePairs.front() == Interaction::BasePair(0, 1)); + REQUIRE((*output.begin())->basePairs.back() == Interaction::BasePair(1, 0)); + } + } + } + } +} diff --git a/tests/PredictorMfeHeuristicCellState_test.cpp b/tests/PredictorMfeHeuristicCellState_test.cpp index 2e66ccfa..c2c4cac5 100644 --- a/tests/PredictorMfeHeuristicCellState_test.cpp +++ b/tests/PredictorMfeHeuristicCellState_test.cpp @@ -8,12 +8,17 @@ #include "IntaRNA/PredictorMfe2dHeuristic.h" #include "IntaRNA/PredictorMfe2dHeuristicSeed.h" #include "IntaRNA/PredictorMfeEns2dHeuristic.h" +#include "IntaRNA/PredictorMfe2dHelixBlockHeuristic.h" +#include "IntaRNA/PredictorMfe2dHelixBlockHeuristicSeed.h" +#include "IntaRNA/SeedHandlerMfe.h" #include "IntaRNA/ReverseAccessibility.h" #include "IntaRNA/RnaSequence.h" #include "IntaRNA/SeedConstraint.h" #include "IntaRNA/SeedHandlerNoBulge.h" #include +#include +#include using namespace IntaRNA; @@ -187,3 +192,112 @@ TEST_CASE("ensemble noLP heuristic keeps valid non-direct extensions", REQUIRE(interaction.basePairs.back() == Interaction::BasePair(4, 0)); } } + +namespace { + +std::unique_ptr makeOutputFilterPredictor(const InteractionEnergy & energy, + OutputHandler & output, const bool seeded, const bool helix) { + static const SeedConstraint seed(2, 2, 2, 2, E_INF, Accessibility::ED_UPPER_BOUND, E_INF, + IndexRangeList(), IndexRangeList(), "", false, false, true); + static const HelixConstraint helixConstraint(2, 4, 2, Accessibility::ED_UPPER_BOUND, E_INF, false); + if (helix) { + if (seeded) { + return std::make_unique( + energy, output, nullptr, helixConstraint, new SeedHandlerMfe(energy, seed)); + } + return std::make_unique(energy, output, nullptr, helixConstraint); + } + if (seeded) { + return std::make_unique(energy, output, nullptr, + new SeedHandlerMfe(energy, seed)); + } + return std::make_unique(energy, output, nullptr); +} + +// A small accessibility penalty confined to one end of the sequence. The +// second disjoint two-pair site is favorable but exceeds an ED limit of zero. +class EndPenaltyAccessibility : public AccessibilityDisabled { +public: + EndPenaltyAccessibility(const RnaSequence & sequence, const bool atStart) + : AccessibilityDisabled(sequence, 0, nullptr), atStart(atStart) {} + + E_type getED(const size_t from, const size_t to) const override { + const E_type base = AccessibilityDisabled::getED(from, to); + return base + ((atStart ? from < 2 : to >= 5) ? Ekcal_2_E(0.25) : 0); + } +private: + const bool atStart; +}; + +} // namespace + +TEST_CASE("heuristic suboptimals respect terminal GU constraints", "[PredictorMfeHeuristicCellState][Overlap]") { + #include "testEasyLoggingSetup.icc" + + for (bool seeded : {false, true}) { + for (bool helix : {false, true}) { + for (bool trace : {false, true}) { + for (size_t offset : {size_t(0), size_t(1)}) { + RnaSequence target("target", offset ? "NUUGAN" : "UUGA"); + RnaSequence query("query", offset ? "NCAUUN" : "CAUU"); + AccessibilityDisabled targetAcc(target, 0, nullptr); + AccessibilityDisabled queryAcc(query, 0, nullptr); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc); + for (auto overlap : {OutputConstraint::OVERLAP_NONE, OutputConstraint::OVERLAP_SEQ1, + OutputConstraint::OVERLAP_SEQ2, OutputConstraint::OVERLAP_BOTH}) { + CAPTURE(seeded, helix, trace, offset, overlap); + OutputConstraint constraint(10, overlap, 0, Ekcal_2_E(100), + false, false, true, false, trace); + OutputHandlerInteractionList output(constraint, 10); + auto predictor = makeOutputFilterPredictor(energy, output, seeded, helix); + predictor->predict(IndexRange(offset, offset+3), IndexRange(offset, offset+3)); + REQUIRE_FALSE(output.empty()); + for (const Interaction * interaction : output) { + REQUIRE(interaction->isValid()); + for (auto bp : {interaction->basePairs.front(), interaction->basePairs.back()}) { + const char t = target.asString().at(bp.first), q = query.asString().at(bp.second); + REQUIRE_FALSE((t == 'G' && q == 'U')); + REQUIRE_FALSE((t == 'U' && q == 'G')); + } + } + } + } + } + } + } +} + +TEST_CASE("heuristic suboptimals respect complete-site accessibility limits", "[PredictorMfeHeuristicCellState][Overlap]") { + #include "testEasyLoggingSetup.icc" + + RnaSequence target("target", "CCCAACC"); + RnaSequence query("query", "GGAAGGG"); + EndPenaltyAccessibility targetAcc(target, false), queryAcc(query, true); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc, 0, 0); + for (bool seeded : {false, true}) { + for (bool helix : {false, true}) { + for (bool trace : {false, true}) { + for (auto overlap : {OutputConstraint::OVERLAP_NONE, + OutputConstraint::OVERLAP_SEQ1, OutputConstraint::OVERLAP_SEQ2}) { + CAPTURE(seeded, helix, trace, overlap); + OutputConstraint constraint(10, overlap, 0, Ekcal_2_E(100), + false, false, false, false, trace, 0); + OutputHandlerInteractionList output(constraint, 10); + auto predictor = makeOutputFilterPredictor(energy, output, seeded, helix); + predictor->predict(); + REQUIRE(std::distance(output.begin(), output.end()) == 1); + for (const Interaction * interaction : output) { + REQUIRE(interaction->isValid()); + REQUIRE(targetAcc.getED(interaction->basePairs.front().first, + interaction->basePairs.back().first) == 0); + REQUIRE(queryAcc.getED(interaction->basePairs.back().second, + interaction->basePairs.front().second) == 0); + REQUIRE(interaction->energy == Ekcal_2_E(-3)); + } + } + } + } + } +} diff --git a/tests/PredictorTinyOracle_test.cpp b/tests/PredictorTinyOracle_test.cpp index 45ae2ea9..015b5fb4 100644 --- a/tests/PredictorTinyOracle_test.cpp +++ b/tests/PredictorTinyOracle_test.cpp @@ -7,6 +7,7 @@ #include "IntaRNA/InteractionEnergyBasePair.h" #include "IntaRNA/OutputHandlerInteractionList.h" #include "IntaRNA/PredictorMfe2d.h" +#include "IntaRNA/PredictorMfe2dHeuristic.h" #include "IntaRNA/PredictorMfeEns2d.h" #include "IntaRNA/ReverseAccessibility.h" #include "IntaRNA/RnaSequence.h" @@ -375,3 +376,50 @@ TEST_CASE("tiny exhaustive oracle for exact predictors", "[PredictorTinyOracle]" == Approx(oracle.partition).epsilon(1e-12)); } } + +TEST_CASE("restricted heuristic results belong to the filtered tiny oracle", "[PredictorTinyOracle][Overlap]") { + #include "testEasyLoggingSetup.icc" + + RnaSequence target("target", "UUGA"), query("query", "CAUU"); + WidthAccessibility targetAcc(target, 0, Ekcal_2_E(0.1)); + WidthAccessibility queryAcc(query, 0, Ekcal_2_E(0.2)); + ReverseAccessibility reverseQueryAcc(queryAcc); + InteractionEnergyBasePair energy(targetAcc, reverseQueryAcc); + for (bool noGU : {false, true}) { + for (E_type maxED : {Ekcal_2_E(0.5), Accessibility::ED_UPPER_BOUND}) { + for (auto overlap : {OutputConstraint::OVERLAP_NONE, OutputConstraint::OVERLAP_SEQ1, + OutputConstraint::OVERLAP_SEQ2, OutputConstraint::OVERLAP_BOTH}) { + CAPTURE(noGU, maxED, overlap); + OutputConstraint constraint(20, overlap, 0, Ekcal_2_E(100), + false, false, noGU, false, true, maxED); + const auto oracle = enumerateInteractions(energy, constraint, IndexRange(0,3), IndexRange(0,3)); + OutputHandlerInteractionList output(constraint, 20); + PredictorMfe2dHeuristic predictor(energy, output, nullptr); + predictor.predict(); + REQUIRE_FALSE(output.empty()); + for (const Interaction * interaction : output) { + bool found = false; + for (const auto & candidate : oracle.interactions) { + if (candidate.energy == interaction->energy + && toBasePairs(energy, candidate.chain) == interaction->basePairs) { + found = true; + break; + } + } + REQUIRE(found); + for (const Interaction * other : output) { + if (other == interaction) { break; } + if (overlap == OutputConstraint::OVERLAP_NONE || overlap == OutputConstraint::OVERLAP_SEQ2) { + REQUIRE_FALSE((interaction->basePairs.front().first <= other->basePairs.back().first + && other->basePairs.front().first <= interaction->basePairs.back().first)); + } + if (overlap == OutputConstraint::OVERLAP_NONE || overlap == OutputConstraint::OVERLAP_SEQ1) { + REQUIRE_FALSE((interaction->basePairs.back().second <= other->basePairs.front().second + && other->basePairs.back().second <= interaction->basePairs.front().second)); + } + } + } + } + } + } +} diff --git a/tests/runOutputOverlap.sh b/tests/runOutputOverlap.sh new file mode 100644 index 00000000..580dbac1 --- /dev/null +++ b/tests/runOutputOverlap.sh @@ -0,0 +1,131 @@ +#!/usr/bin/env bash +# Regional overlap guards, including automatic decomposition. +set -euo pipefail +bin="$INTARNABINPATH/src/bin/IntaRNA" +tmp=$(mktemp -d) +trap 'rm -rf "$tmp"' EXIT +region_acc=C +common=(--energy=B --seedBP=2 --threads=1 --outMode=C + --outCsvCols=start1,end1,start2,end2,E --outNumber=10 --default-log-file="$tmp/info.log") + +check_regions() { + local expect=$1 + shift + : > "$tmp/info.log" + if "$bin" "${common[@]}" --acc="$region_acc" "$@" > "$tmp/result" 2> "$tmp/error"; then + if test "$expect" != ok; then + echo "Expected regional-overlap rejection: $*" >&2 + exit 1 + fi + else + if test "$expect" != reject; then cat "$tmp/error" "$tmp/info.log" "$tmp/result" >&2; exit 1; fi + grep -q 'regions cannot be merged with --outOverlap=' "$tmp/error" "$tmp/info.log" "$tmp/result" + fi +} + +for overlap in B N T Q; do + for per_region in false true; do + for target_regions in 1 2; do + for query_regions in 1 2; do + expect=ok + if test "$per_region" = false; then + if test "$target_regions" = 2 && [[ "$overlap" = N || "$overlap" = T ]]; then expect=reject; fi + if test "$query_regions" = 2 && [[ "$overlap" = N || "$overlap" = Q ]]; then expect=reject; fi + fi + tr=1-2; qr=1-2 + if test "$target_regions" = 2; then tr+=,5-6; fi + if test "$query_regions" = 2; then qr+=,5-6; fi + check_regions "$expect" -t CCAACC -q GGAAGG --tRegion="$tr" --qRegion="$qr" \ + --outOverlap="$overlap" --outPerRegion="$per_region" + done + done + # Decomposition really produces multiple regions on the selected RNA. + for side in target query; do + expect=ok + if test "$side" = target; then + args=(-t UAUCGGCC -q GG --tRegionLenMax=4) + if test "$per_region" = false && [[ "$overlap" = N || "$overlap" = T ]]; then expect=reject; fi + else + args=(-t GG -q UAUCGGCC --qRegionLenMax=4) + if test "$per_region" = false && [[ "$overlap" = N || "$overlap" = Q ]]; then expect=reject; fi + fi + check_regions "$expect" "${args[@]}" --outOverlap="$overlap" --outPerRegion="$per_region" + done + done + # An automatic-region option alone is fine if no split is needed. + check_regions ok -t CC -q GG --tRegionLenMax=4 --qRegionLenMax=4 --outOverlap="$overlap" +done + +# outMinPu also decomposes regions; blocked positions supply deterministic gaps. +region_acc=N +for overlap in B N T Q; do + for per_region in false true; do + for side in target query; do + expect=ok + if test "$side" = target; then + args=(-t CCAACC -q GG --tAccConstr=b:3-4) + if test "$per_region" = false && [[ "$overlap" = N || "$overlap" = T ]]; then expect=reject; fi + else + args=(-t CC -q GGAAGG --qAccConstr=b:3-4) + if test "$per_region" = false && [[ "$overlap" = N || "$overlap" = Q ]]; then expect=reject; fi + fi + check_regions "$expect" "${args[@]}" --outMinPu=0.5 --outOverlap="$overlap" --outPerRegion="$per_region" + done + done +done + +region_acc=C + +# Shifted indices still obey the same manual-region rule. +check_regions reject -t CCAACC -q GG --tIdxPos0=10 --qIdxPos0=20 --tRegion=10-11,14-15 --outOverlap=T +check_regions ok -t CCAACC -q GG --tIdxPos0=10 --qIdxPos0=20 --tRegion=10-11,14-15 --outOverlap=Q +printf 'start1;end1;start2;end2;E\n10;11;20;21;-2\n14;15;20;21;-2\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/result" +echo 'Regional output overlap checks passed' + +# Global energy windows must be independent of which region supplies the MFE. +# In this base-pair energy model the optimum is exactly -2; a one-pair site is -1. +for overlap in B T Q; do + for best_first in true false; do + if test "$best_first" = true; then seq=CCAAC; regions=1-2,5-5; + else seq=CAACC; regions=1-1,4-5; fi + if test "$overlap" = T; then + seq=${seq//C/G} + args=(-t CC -q "$seq" --qRegion="$regions") + else + args=(-t "$seq" -q GG --tRegion="$regions") + fi + for per_region in false true; do + for delta in 0 0.99 1; do + "$bin" --energy=B --acc=N --noSeed --model=S --mode=M --threads=1 \ + --outMode=C --outCsvCols=E --outNumber=10 --outOverlap="$overlap" \ + --outPerRegion="$per_region" --outDeltaE="$delta" "${args[@]}" \ + --default-log-file="$tmp/info.log" > "$tmp/energies" + awk -v delta="$delta" -v per_region="$per_region" ' + NR == 2 { if ($0 != -2) exit 1 } + NR > 1 { if ($0 != -2 && $0 != -1) exit 1; if ($0 == -1) weak++ } + END { + if (NR < 2) exit 1 + if (per_region == "false" && delta < 1 && weak) exit 1 + if ((per_region == "true" || delta == 1) && !weak) exit 1 + }' "$tmp/energies" + done + done + done +done + +# With independent regions all four modes retain the weaker region's optimum. +for overlap in B N T Q; do + "$bin" --energy=B --acc=N --noSeed --model=S --threads=1 -t CCAAC -q GG \ + --tRegion=1-2,5-5 --outMode=C --outCsvCols=E --outNumber=10 \ + --outOverlap="$overlap" --outPerRegion --outDeltaE=0 \ + --default-log-file="$tmp/info.log" > "$tmp/energies" + grep -q '^-1$' "$tmp/energies" +done + +# An empty merged result has no minimum to dereference. +"$bin" --energy=B --acc=N --noSeed --threads=1 -t AAAAA -q AAAAA --tRegion=1-2,4-5 \ + --outMode=C --outCsvCols=E --outDeltaE=0 --default-log-file="$tmp/info.log" > "$tmp/energies" +printf 'E\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/energies" +echo 'Global and per-region energy window checks passed'