diff --git a/ChangeLog b/ChangeLog index 07e8fc1c..8c272627 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,13 @@ ## Interface and handling +- deterministic kinetic seed extension with --model=X --mode=K and configurable + local scoring; complete-energy downhill steps, atomic noLP loop/stack moves + and retained greedy traceback; equilibrium probability output is unsupported +- kinetic mode always uses stacked extensions and trusts handler-provided seeds; + reuse candidate pairing/local energies +- IntaRNAsnap personality enables kinetic seed extension with noLP by default + - apply outDeltaE relative to the sequence pair's best interaction when merging regions; preserve local windows for outPerRegion output (PR #253) @@ -35,6 +42,8 @@ ## Technical changes and Optimizations +- install personality links in out-of-tree builds, including IntaRNAsnap + - BUGFIX : normalize single-pair suboptimal boundaries before traceback and boundary-only output validation (PR #253) @@ -76,6 +85,26 @@ - BUGFIX : multi-threading : IntaRNAsTar was not thread-safe due to shared storage of computed ranges +# IntaRNAsnap + +New personality IntaRNAsnap identifies the fastest "folding path" to a (locally) +optimal interaction after seed formation, rather than finding the global optimum +under the assumption of thermodynamic equilibrium. + +To this end, IntaRNAsnap implements a kinetics-motivated greedy search for the +best interaction that can be reached from the given seed(s) by a series of +complete-energy downhill local moves. +Local moves are direct stack extensions or 2-base-pair loop/stack moves, i.e. +the same moves as used in the IntaRNA noLP model. That way, kinetic energy +barriers posed by the formation of loops and bulges are avoided, which mimics +the continuation of the zipping process after "jumping" over a loop or bulge. + +This personality is only available with --model=X and --mode=K, and equilibrium +probability output is unsupported. The local scoring can be configured with +--kineticScore=A|B|C, and the output can include distinct structures, energies, +restricted partition sums, and trackers. + + # IntaRNAeval New personality of IntaRNA to evaluate predefined RNA-RNA interactions. @@ -97,6 +126,38 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261005 Alexander Mitrofanov + * bin/CommandLineParsing, tests/runKineticSeedExtension.sh : + + add IntaRNAsnap personality with model=X, mode=K and outNoLP=true defaults; + support executable-name and explicit personality selection + * configure.ac : + * find personality declarations via srcdir so out-of-tree installations + create the executable links, including IntaRNAsnap + * README.md, doc/kinetic-seed-extension.md, + doc/recursions/IntaRNAsnap.PredictorSeedExtensionKinetic.svg, doc/Makefile.am : + + document personality usage and depict seed initialization, allowed moves, + complete-energy greedy selection, stopping and prefix reporting + + distribute the benchmark's tutorial sequences with its script and results + * doc/benchmark-kix.py, doc/kix-benchmark-20261005.json : + + compare IntaRNAsnap with default IntaRNA with and without GU-end constraints; + record wall time, peak RSS, MFE energy and maximum covered strand length + * implement https://github.com/BackofenLab/IntaRNA/pull/254#issuecomment-5995041057 + +261005 Alexander Mitrofanov + * IntaRNA/PredictorSeedExtensionKinetic : + * always evaluate single and double stacks plus atomic loop/stack extensions; + trust seed-handler structures/energies, including explicit lonely pairs + * cache each end's candidates, shared complementarity checks and local + energies; rebuild only the chosen end and refresh full opposite-end energy + * bin/CommandLineParsing : + * set --outNoLP=true with INFO for kinetic mode K + * tests/PredictorSeedExtensionKinetic_test.cpp, tests/runKineticSeedExtension.sh : + * update independent oracle and CLI expectations for atomic double stacks; + cover trusted seeds, caching and kinetic mode K + * README.md, doc/kinetic-seed-extension.md : + * document revised semantics in response to + https://github.com/BackofenLab/IntaRNA/pull/254 + 261005 Alexander Mitrofanov * IntaRNA/PredictorMfe, tests/PredictorMfeHeuristicCellState_test.cpp, tests/PredictorMfeEnsRegression_test.cpp : @@ -162,6 +223,23 @@ energies, restricted partition sums, and trackers. + define suboptimal predictions by distinct start/end coordinates on both RNAs * distinguish site prediction from evaluation of supplied structures (PR #253) +261003 Alexander Mitrofanov + * IntaRNA/PredictorSeedExtensionKinetic, src/IntaRNA/Makefile.am : + + greedily extend seeds on either side using complete interaction-energy + differences and deterministic thermodynamic or distance-weighted scores + + support noLP macro-steps, GU restrictions, per-strand loop/span constraints, + explicit-seed validation and cached trajectory traceback + + retain valid visited prefixes and select non-overlapping output from the + complete retained candidate set; reject unsupported ensemble statistics + * bin/CommandLineParsing : + + expose --mode=K exclusively for --model=X and --kineticScore=A|B|C + * tests/PredictorSeedExtensionKinetic_test.cpp, tests/runKineticSeedExtension.sh, + tests/Makefile.am : + + validate local move choices, macro-step energetics, structural constraints, + traceback, output filtering and CLI compatibility + * README.md, doc/kinetic-seed-extension.md, doc/Makefile.am : + + document scientific scope and complete-energy move enumeration + 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 diff --git a/README.md b/README.md index 6466b700..ed4cb9e3 100644 --- a/README.md +++ b/README.md @@ -100,6 +100,7 @@ The following topics are covered by this documentation: - [IntaRNAsTar - optimized for sRNA-target prediction](#IntaRNAsTar) - [IntaRNAseed - identifys and reports seed interactions only](#IntaRNAseed) - [IntaRNAens - ensemble-based prediction and partition function computation](#IntaRNAens) + - [IntaRNAsnap - kinetic seed extension](#IntaRNAsnap) - [IntaRNAeval - evaluate predefined interactions](#IntaRNAeval) - [How to constrain predicted interactions](#constraintSetup) - [Interaction restrictions](#interConstr) @@ -729,6 +730,40 @@ minimum free energy interaction. Putative seed interactions (used by the `H` and `M` mode) can be enumerated and studied using the `S` mode. +### Greedy kinetic seed extension + +`--model=X --mode=K` grows each available seed along a deterministic greedy +path. Every step compares feasible extensions on both sides using the complete +change in interaction energy, including accessibility, terminal penalties and +dangling ends. Only strictly negative changes are accepted. Extensions always +add one stacked pair or two stacked pairs, including across a loop; two-pair +moves are evaluated and committed atomically. The CLI sets `--outNoLP=true` +when absent or false and logs an INFO message. Seeds and their energies are +accepted from the seed handler, including explicit seeds with lonely pairs. + +`--kineticScore` selects the local move ranking: + +| Value | Score minimized for a move with gaps `s1`, `s2` | +| --- | --- | +| `A` (default) | Complete energy change | +| `B` | Complete energy change / `(1+s1+s2)` | +| `C` (C1 in the design) | Complete energy change / `(1+2*max(s1,s2))` | + +Equal scores prefer the left side, then fewer unpaired bases, then smaller +`s1`, then the single-pair move. The denominators also apply to two-pair moves. All reportable visited +states, including seeds, participate in the normal energy-ranked output; +traceback preserves the actual chosen path. `--outNoGUend`, separate query and +target loop/span limits, regions, output energy/accessibility filters and +overlap settings remain applicable. + +The [IntaRNAsnap personality](#IntaRNAsnap) selects this mode with noLP enabled +by default. It is a zippering-inspired heuristic without a calibrated time axis +or a guarantee of the global minimum. Equilibrium probability/partition-sum +outputs are rejected, as are other models and `--noSeed`. Scores B and C are +optional distance preferences, not measured kinetic rates. See the +[algorithm and preliminary benchmark](doc/kinetic-seed-extension.md). + + [![up](doc/figures/icon-up.28.png) back to overview](#overview)

@@ -1003,6 +1038,32 @@ IntaRNA --mode=S ... [![up](doc/figures/icon-up.28.png) back to overview](#overview) +### IntaRNAsnap + +**IntaRNAsnap** (kinetic seed extension) grows every handler-provided seed by +choosing the most favorable complete energy change at either end. It sets +`--model=X --mode=K --outNoLP=true`; other defaults are those of IntaRNA. +Each move adds one stacked pair, two stacked pairs, or a loop-closing pair +plus its outward stack. Two-pair moves are evaluated and committed together. +Seeds themselves may contain lonely pairs when supplied by the seed handler. + +The following calls are equivalent: + +```sh +IntaRNAsnap -t target.fasta -q query.fasta +IntaRNA --personality=IntaRNAsnap -t target.fasta -q query.fasta +IntaRNA --model=X --mode=K --outNoLP=true -t target.fasta -q query.fasta +``` + +Only strictly downhill moves are accepted. The default score A chooses the +largest energy decrease; `--kineticScore=B|C` adds distance preferences. +The reported MFE is the best visited, reportable interaction across seeds; +this heuristic greedy search has no global-optimum or physical folding-time guarantee. + + +[![up](doc/figures/icon-up.28.png) back to overview](#overview) + + ### IntaRNAeval **IntaRNAeval** evaluates predefined RNA-RNA interactions with the selected diff --git a/configure.ac b/configure.ac index 3b52ea77..ae2df6d9 100644 --- a/configure.ac +++ b/configure.ac @@ -514,7 +514,7 @@ AS_IF([test "$DEPENDENCYNOTFOUND" = "1"], [ ########################################################################## # get available personalities -AC_SUBST([PERSONALITIES],[`grep 'return.*"IntaRNA..*"' ./src/bin/CommandLineParsing.h | grep -o '".*"' | tr -d '"' | tr "\n" " "`]) +AC_SUBST([PERSONALITIES],[`grep 'return.*"IntaRNA..*"' "$srcdir/src/bin/CommandLineParsing.h" | grep -o '".*"' | tr -d '"' | tr "\n" " "`]) ########################################################################## # Keep strict diagnostics focused on IntaRNA sources. Boost and ViennaRNA diff --git a/doc/Makefile.am b/doc/Makefile.am index ac9bab1a..c28d05ec 100644 --- a/doc/Makefile.am +++ b/doc/Makefile.am @@ -4,8 +4,6 @@ ################################################################ EXTRA_DIST = \ - analysis/out-overlap.md \ - analysis/out-overlap/reproduce.py \ conda.txt \ doxygen.cfg \ latex-deps/adjcalc.sty \ @@ -15,4 +13,3 @@ EXTRA_DIST = \ latex-deps/tocloft.sty \ latex-deps/trimclip.sty \ latex-deps/xtab.sty - diff --git a/doc/benchmark-kix.py b/doc/benchmark-kix.py new file mode 100644 index 00000000..34f1d399 --- /dev/null +++ b/doc/benchmark-kix.py @@ -0,0 +1,149 @@ +#!/usr/bin/env python3 +"""Small, single-thread comparison of IntaRNAsnap and default IntaRNA. + +Requires Python 3, GNU time and a release IntaRNA binary. Run from any directory: + python3 doc/benchmark-kix.py /path/to/IntaRNA --output results.json +""" + +import argparse +import csv +import hashlib +import io +import json +import os +from pathlib import Path +import platform +import statistics +import subprocess +import sys +import tempfile +import time + + +PAIRS = [("fhlA", "OxyS"), ("phoB", "GcvB"), ("ilvE", "GcvB.ST")] +COLUMNS = "start1,end1,start2,end2,E,hybridDB" + + +def read_fasta(path): + lines = path.read_text().splitlines() + assert sum(line.startswith(">") for line in lines) == 1, path + sequence = "".join(line.strip() for line in lines if not line.startswith(">")) + return {"file": "handson/" + path.name, "header": lines[0][1:], + "length_nt": len(sequence), "sequence": sequence, + "sha256": hashlib.sha256(path.read_bytes()).hexdigest()} + + +def prediction(stdout): + rows = list(csv.DictReader(io.StringIO(stdout), delimiter=";")) + if not rows: + return None + assert len(rows) == 1, rows + row = rows[0] + result = {key: int(row[key]) for key in ("start1", "end1", "start2", "end2")} + result.update(E_kcal_mol=float(row["E"]), hybridDB=row["hybridDB"]) + # All inputs use the default one-based, ascending coordinates. + result["length_nt"] = max(result["end1"] - result["start1"] + 1, + result["end2"] - result["start2"] + 1) + return result + + +def distribution(samples): + return {"median": statistics.median(samples), "min": min(samples), "max": max(samples)} + + +def main(): + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("binary", type=Path) + parser.add_argument("--output", required=True, type=Path) + parser.add_argument("--repetitions", type=int, default=5) + parser.add_argument("--warmups", type=int, default=1) + parser.add_argument("--time", type=Path, default=Path("/usr/bin/time")) + parser.add_argument("--build-description", default="unspecified") + args = parser.parse_args() + if args.repetitions < 1 or args.warmups < 0: + parser.error("repetitions must be positive and warmups nonnegative") + binary = args.binary.resolve() + timer = args.time.resolve() + fixtures = Path(__file__).resolve().parent / "handson" + common = ["--threads=1", "--outMode=C", "--outCsvCols=" + COLUMNS, + "--outNumber=1", "--default-log-file=/dev/null"] + env = dict(os.environ, OMP_NUM_THREADS="1") + report = { + "schema": 1, + "binary_sha256": hashlib.sha256(binary.read_bytes()).hexdigest(), + "version": subprocess.check_output([str(binary), "--version"], text=True).strip(), + "build": args.build_description, + "platform": platform.platform(), "machine": platform.machine(), + "python": platform.python_version(), + "timer": subprocess.check_output([str(timer), "--version"], text=True).splitlines()[0], + "common_arguments": common, + "repetitions": args.repetitions, "warmups_per_configuration": args.warmups, + "timing": "perf_counter wall seconds around GNU time + binary; includes startup, folding, seed search and prediction", + "memory": "GNU time %M: maximum resident set size of each child process, KiB", + "noGU": "--outNoGUend=true; seedNoGU remains false for both programs", + "baseline": "default IntaRNA (model X, mode H, outNoLP false)", + "comparison": "IntaRNAsnap defaults (model X, mode K, outNoLP true, kineticScore A)", + "length": "max(end1-start1+1, end2-start2+1), nt; default one-based coordinates", + "deviations": "signed IntaRNAsnap minus default IntaRNA, within each GU setting; not a global-optimum error bound", + "cases": [], + } + with tempfile.TemporaryDirectory(prefix="intarna-kix-benchmark-") as directory: + rss_file = Path(directory) / "time.txt" + for target, query in PAIRS: + inputs = ["--target=" + str(fixtures / (target + ".fasta")), + "--query=" + str(fixtures / (query + ".fasta"))] + case = {"name": target + "/" + query, + "target": read_fasta(fixtures / (target + ".fasta")), + "query": read_fasta(fixtures / (query + ".fasta")), "runs": [], "comparisons": []} + configurations = [] + for no_gu in (False, True): + for program in ("IntaRNA", "IntaRNAsnap"): + options = [] if program == "IntaRNA" else ["--personality=IntaRNAsnap"] + options.append("--outNoGUend=" + str(no_gu).lower()) + configurations.append((program, no_gu, options)) + case["runs"].append({"program": program, "noGU": no_gu, + "arguments": options, "samples": []}) + # Warm each configuration; rotate execution order every round. + for iteration in range(args.warmups + args.repetitions): + for offset in range(len(configurations)): + index = (iteration + offset) % len(configurations) + program, no_gu, options = configurations[index] + run = case["runs"][index] + command = [str(binary), *common, *inputs, *options] + start = time.perf_counter() + process = subprocess.run([str(timer), "-f", "%M", "-o", str(rss_file), + *command], env=env, text=True, capture_output=True, check=True) + elapsed = time.perf_counter() - start + result = prediction(process.stdout) + if "prediction" in run: + assert run["prediction"] == result, (case["name"], options, process.stdout) + run["prediction"] = result + if iteration >= args.warmups: + run["samples"].append({"wall_s": elapsed, "peak_rss_kib": int(rss_file.read_text())}) + for run in case["runs"]: + run["wall_s"] = distribution([sample["wall_s"] for sample in run["samples"]]) + run["peak_rss_kib"] = distribution([sample["peak_rss_kib"] for sample in run["samples"]]) + # Outside timing, reevaluate each reported structure independently. + if run["prediction"] is not None: + evaluated = subprocess.check_output([str(binary), *common, *inputs, + "--rri=" + run["prediction"]["hybridDB"]], + env=env, text=True) + assert prediction(evaluated) == run["prediction"], (case["name"], run) + run["structure_reevaluation_matches"] = True + for index, no_gu in ((0, False), (2, True)): + baseline, kix = case["runs"][index:index + 2] + available = baseline["prediction"] is not None and kix["prediction"] is not None + case["comparisons"].append({ + "noGU": no_gu, + "delta_E_kcal_mol": round(kix["prediction"]["E_kcal_mol"] - baseline["prediction"]["E_kcal_mol"], 2) if available else None, + "delta_length_nt": kix["prediction"]["length_nt"] - baseline["prediction"]["length_nt"] if available else None, + "wall_ratio_kix_over_default": kix["wall_s"]["median"] / baseline["wall_s"]["median"], + "rss_ratio_kix_over_default": kix["peak_rss_kib"]["median"] / baseline["peak_rss_kib"]["median"], + }) + report["cases"].append(case) + print(case["name"] + ": " + json.dumps(case["comparisons"]), file=sys.stderr, flush=True) + args.output.write_text(json.dumps(report, indent=2) + "\n") + + +if __name__ == "__main__": + main() diff --git a/doc/kinetic-seed-extension.md b/doc/kinetic-seed-extension.md new file mode 100644 index 00000000..23cb48d1 --- /dev/null +++ b/doc/kinetic-seed-extension.md @@ -0,0 +1,195 @@ +# Deterministic kinetic seed extension + +## Scope and modes + +`--model=X --mode=K` follows a deterministic downhill path from every seed +provided by the selected seed handler. It compares both ends and commits the +best strictly favorable complete move. This is a zippering-inspired heuristic: +it has no calibrated transition rates or time axis, does not cross barriers +between committed states, and does not guarantee a global minimum. A favorable +two-pair move does not establish a barrier-free physical reaction pathway. + +The `IntaRNAsnap` personality (kinetic seed extension) selects +`--model=X --mode=K --outNoLP=true`. It can be invoked through the installed +`IntaRNAsnap` executable link or `IntaRNA --personality=IntaRNAsnap`. Other defaults +remain those of IntaRNA. Ordinary energy trackers are supported; seedless +operation, other interaction models and equilibrium partition/probability +requests are rejected for mode K. + +![Kinetic seed extension recursion](recursions/IntaRNAsnap.PredictorSeedExtensionKinetic.svg) + +## Seeds, states and allowed extensions + +The predictor trusts the seed handler's structure and hybridization energy. +It does not repeat seed complementarity, noLP or GU-loop checks. Explicit seeds +may contain lonely pairs, including at their ends. Finite energy, prediction +ranges and per-strand span limits still apply. Seed annotations come from the +seed handler without an additional predictor-specific whitelist. + +A state stores its complete ordered base-pair chain, inclusive boundaries +`(i1,j1,i2,j2)`, hybridization energy `H`, and full interaction energy `E`. +Sequence 2 uses reversed energy indices; existing wrappers handle prediction +range offsets and conversion to original coordinates. Initially, +`H = seedHandler.getSeedE(i1,i2) + energy.getE_init()`. + +Extensions **always** use the no-lonely-pair strategy, independent of the API +output constraint. Direct `--mode=K` calls promote a missing or false +`--outNoLP` to true with an INFO message using the normal logging destination; +IntaRNAsnap already defaults to true. This applies to new +extensions, not to revalidation of handler-provided seeds. Allowed moves are: + +- One pair stacked directly onto the current boundary (`|SEED`). +- Two successive stacked pairs (`||SEED`), evaluated atomically. +- A loop-closing pair followed immediately by its outward stack (`//.SEED`). + +For each strand, `sk` denotes the number of skipped bases between the current +boundary and the closing pair. A single-pair move has `s1=s2=0`. Two-pair moves +include all gap pairs from zero through the separate strand loop limits, +subject to available sequence range and maximum interaction span. They advance +each boundary by `sk+2`. Neither an isolated closing pair nor an intermediate +state of a two-pair move is separately committed or reported. + +Under `--outNoGUend`, a newly formed nonstacking loop must have non-GU closing +pairs. Stacking can temporarily expose a GU outer endpoint; a state is only +reported if its outer endpoints satisfy the flag. The energy model's own +internal-loop GU policy also applies. Earlier reportable prefixes remain +available if a trajectory stops at an unreportable endpoint. + +## Complete energies and deterministic scores + +Every geometrically feasible, complementary candidate is evaluated +with the active energy model: + +``` +H_next = H_current + E_loop + E_stack_if_two_pairs +E_next = energy.getE(i1_next, j1_next, i2_next, j2_next, H_next) +delta = E_next - E_current +``` + +`getE()` includes both accessibility penalties, terminal terms, both weighted +dangling ends, and the configured additive term. Updating one end can change +the opposite end's dangling weight, so complete energies must be refreshed +even when local loop energies are reused. Infinite values are excluded before +subtraction. Only `delta < 0` is accepted. Output energy and accessibility +thresholds filter reports, without imposing additional trajectory barriers. + +| Score | Quantity minimized | +| --- | --- | +| A (default) | `delta` | +| B | `delta / (1+s1+s2)` | +| C | `delta / (1+2*max(s1,s2))` | + +B/C are optional uncalibrated distance preferences. Their denominators count +one move, including two-pair moves. Equal scores prefer left, smaller total +gap, smaller first-strand gap, then a single-pair move. Wide integer cross +products avoid rounding during comparison. Output remains ranked by full `E`. +Two-pair stacking can skip a prefix that the previous implementation visited; +only states actually committed by the revised walk are retained. + +## Candidate storage and reuse + +Each end has a contiguous rectangular table of two-pair candidates plus its +single-stack candidate. A shared position-pair table records unknown, possible +or impossible complementarity. The closing pair of one candidate may be the +outer pair of another; each such check is resolved once per unchanged end. +The first pass resolves pairing and GU feasibility before the energy pass. +The second pass caches local loop-plus-stack energies and updates the best +candidate as it evaluates full energies, without a separate selection pass. + +After committing a move, only that end's tables are rebuilt. The opposite +end's pair checks, feasibility and local energies remain valid. Its total +energy and span eligibility are refreshed. An +initially uphill candidate is retained because it may become downhill when +the opposite end changes. Tables use O((m1+2)(m2+2)) space per end. + +## Reporting and validation + +Retain each reportable visited prefix and its actual pair chain. For duplicate +boundaries, keep the lowest full energy, then lexicographically smallest chain. +Reduce duplicates before the normal optimum collector. Overlap-constrained +output can select shorter retained prefixes; traceback restores the selected +path directly. Memory for retained paths is proportional to their total length, +not constant per seed. Repeated predictions reset trajectory and candidate caches. + +The tests compare K with an independent absolute-endpoint oracle that rebuilds +and reevaluates whole chains. They cover scores/ties, single and double stacks, +positive-loop rescue, strict stopping, separate spans and regions, nonmonotone +ED, GU restrictions, retained prefixes, explicit seeds and annotations, cache +reuse and repeated calls. CLI tests exercise the personality name and option, +explicit parameter overrides, automatic noLP INFO logging, incompatible requests, and +independent reevaluation of predicted structures through `--rri`. + +## Preliminary benchmark against default IntaRNA + +The small panel uses the repository's tutorial sequences: fhlA/OxyS +(112/108 nt), phoB/GcvB (299/201 nt), and ilvE/GcvB.ST (299/200 nt). +These are the pairs in [hands-on examples 3.2, 3.4 and 3.5](handson/README.md). +No experimental seed, region or accessibility constraints from those examples +are applied here. The raw record includes the sequences and input hashes. + +The comparison uses the actual personality defaults: IntaRNA has model X, +mode H and `outNoLP=false`; IntaRNAsnap has model X, mode K, score A and +`outNoLP=true`. Thus energy and length deviations reflect both the search and +the different noLP defaults. Default IntaRNA is itself a heuristic, so the +energy deviation is not a certified error from a global optimum. + +The review's “noGU” setting is interpreted as `--outNoGUend=true`, the same flag +for both programs. This prohibits GU at reported interaction ends and at +nonstacking loop ends; it does **not** forbid all internal GU pairs. +`--seedNoGU` stays at its default false. The other setting explicitly uses +`--outNoGUend=false`. + +Each configuration has one warm-up followed by five measured runs. Execution +order rotates across the four configurations on each pair. All runs use one +thread and compute accessibility from the input sequences. Wall time includes +process startup, accessibility, seed search and prediction; GNU time supplies +the child process's peak resident memory in KiB. The table reports medians; +all samples and min/max values are in the raw JSON. Reported structures and +energies are deterministic across repetitions, and every result is independently +reevaluated via `--rri` outside the timed runs. + +For each reported MFE interaction, covered length is +`L = max(end1-start1+1, end2-start2+1)` in the default one-based coordinates. +Signed deviations are `E_kix - E_default` in kcal/mol and `L_kix - L_default` +in nucleotides, compared within the same GU setting. Missing predictions are +recorded as null, never as zero energy or zero length. + +Measured on 2026-10-05 on Linux x86-64, AMD Ryzen 5 7530U, using GCC 14.4.0 +release (`-O3`, C++23), ViennaRNA 2.7.2, Boost 1.85 and Kokkos mdspan. +No builds or other validation jobs ran alongside these measurements. + +| Pair | noGU | Default time (s) | Kix time (s) | Default peak RSS (KiB) | Kix peak RSS (KiB) | +| --- | --- | ---: | ---: | ---: | ---: | +| fhlA/OxyS | off | 0.0569 | 0.0381 | 16,912 | 16,648 | +| fhlA/OxyS | on | 0.0559 | 0.0381 | 16,660 | 16,796 | +| phoB/GcvB | off | 1.0598 | 0.1491 | 18,448 | 18,456 | +| phoB/GcvB | on | 0.7158 | 0.1446 | 18,456 | 18,452 | +| ilvE/GcvB.ST | off | 1.2935 | 0.1471 | 18,440 | 18,512 | +| ilvE/GcvB.ST | on | 1.0233 | 0.1436 | 18,560 | 18,432 | + +| Pair | noGU | Default E | Kix E | ΔE (kcal/mol) | Default L | Kix L | ΔL (nt) | +| --- | --- | ---: | ---: | ---: | ---: | ---: | ---: | +| fhlA/OxyS | off | -5.59 | -5.57 | +0.02 | 24 | 7 | -17 | +| fhlA/OxyS | on | -5.57 | -5.57 | +0.00 | 7 | 7 | +0 | +| phoB/GcvB | off | -15.70 | -13.19 | +2.51 | 47 | 8 | -39 | +| phoB/GcvB | on | -13.19 | -13.19 | +0.00 | 8 | 8 | +0 | +| ilvE/GcvB.ST | off | -14.24 | -9.84 | +4.40 | 55 | 13 | -42 | +| ilvE/GcvB.ST | on | -10.55 | -9.13 | +1.42 | 40 | 10 | -30 | + +On this small panel, Kix uses 0.11–0.68 times the default runtime. Median peak +RSS differs by less than 2%, within the run-to-run variation. Energy deviations +range from 0 to +4.40 kcal/mol, and length deviations from −42 to 0 nt. +With noGU on, both programs report identical interactions for fhlA/OxyS and +phoB/GcvB. These three selected tutorial pairs are a preliminary performance +and output comparison, not a general speedup or biological-accuracy estimate. + +Reproduce from the repository root with a release binary: + +```sh +python3 doc/benchmark-kix.py /path/to/release/src/bin/IntaRNA \ + --repetitions=5 --warmups=1 --output=kix-benchmark.json +``` + +The [script](benchmark-kix.py) uses Python 3 and GNU time. The +[raw results](kix-benchmark-20261005.json) record all samples, predictions, +signed deviations, flags, sequence data, software versions and binary hash. diff --git a/doc/kix-benchmark-20261005.json b/doc/kix-benchmark-20261005.json new file mode 100644 index 00000000..12f91bd6 --- /dev/null +++ b/doc/kix-benchmark-20261005.json @@ -0,0 +1,727 @@ +{ + "schema": 1, + "binary_sha256": "3779a676dd6a96e6edf9fe9a7316c2125466d0c33be51723661eb1513dfcdd64", + "version": "IntaRNA 3.4.1\n using Vienna RNA package 2.7.2 and boost 1.85.0", + "build": "GCC 14.4.0 release (-O3), C++23, ViennaRNA 2.7.2, Boost 1.85, Kokkos mdspan; AMD Ryzen 5 7530U; 2026-10-05", + "platform": "Linux-6.8.0-142-generic-x86_64-with-glibc2.39", + "machine": "x86_64", + "python": "3.12.7", + "timer": "time (GNU Time) UNKNOWN", + "common_arguments": [ + "--threads=1", + "--outMode=C", + "--outCsvCols=start1,end1,start2,end2,E,hybridDB", + "--outNumber=1", + "--default-log-file=/dev/null" + ], + "repetitions": 5, + "warmups_per_configuration": 1, + "timing": "perf_counter wall seconds around GNU time + binary; includes startup, folding, seed search and prediction", + "memory": "GNU time %M: maximum resident set size of each child process, KiB", + "noGU": "--outNoGUend=true; seedNoGU remains false for both programs", + "baseline": "default IntaRNA (model X, mode H, outNoLP false)", + "comparison": "IntaRNAsnap defaults (model X, mode K, outNoLP true, kineticScore A)", + "length": "max(end1-start1+1, end2-start2+1), nt; default one-based coordinates", + "deviations": "signed IntaRNAsnap minus default IntaRNA, within each GU setting; not a global-optimum error bound", + "cases": [ + { + "name": "fhlA/OxyS", + "target": { + "file": "handson/fhlA.fasta", + "header": "fhlA|NC000913|-53..+60|doi:10.1006/jmbi.2000.3942:Fig7", + "length_nt": 112, + "sequence": "AGUUAGUCAAUGACCUUUUGCACCGCUUUGCGGUGCUUUCCUGGAACAACAAAAUGUCAUAUACACCGAUGAGUGAUCUCGGACAACAAGGGUUGUUCGACAUCACUCGGAC", + "sha256": "667e0df26bf5b8635df5fb5b0eca299cb0484d1c07d0bfd93fff674c2cab9c31" + }, + "query": { + "file": "handson/OxyS.fasta", + "header": "OxyS|NC_000913|56..164|doi:10.1006/jmbi.2000.3942:Fig7", + "length_nt": 108, + "sequence": "GAAACGGAGCGGCACCUCUUUUAACCCUUGAAGUCACUGCCCGUUUCGAGAGUUUCUCAACUCGAAUAACUAAAGCCAACGUGAACUUUUGCGGAUCUCCAGGAUCCG", + "sha256": "bd336bf7f10104837bd7f5e56315dd0f03ba2eb9784d88280588b464adf1a30d" + }, + "runs": [ + { + "program": "IntaRNA", + "noGU": false, + "arguments": [ + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 0.054906384088099, + "peak_rss_kib": 16448 + }, + { + "wall_s": 0.055035896017216146, + "peak_rss_kib": 16912 + }, + { + "wall_s": 0.056886381935328245, + "peak_rss_kib": 16912 + }, + { + "wall_s": 0.057665151078253984, + "peak_rss_kib": 17088 + }, + { + "wall_s": 0.05727210489567369, + "peak_rss_kib": 16788 + } + ], + "prediction": { + "start1": 39, + "end1": 60, + "start2": 81, + "end2": 104, + "E_kcal_mol": -5.59, + "hybridDB": "39|||||||.|..|||||..||||&81||||..|||||..|...|||||||", + "length_nt": 24 + }, + "wall_s": { + "median": 0.056886381935328245, + "min": 0.054906384088099, + "max": 0.057665151078253984 + }, + "peak_rss_kib": { + "median": 16912, + "min": 16448, + "max": 17088 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAsnap", + "noGU": false, + "arguments": [ + "--personality=IntaRNAsnap", + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 0.038167160004377365, + "peak_rss_kib": 16648 + }, + { + "wall_s": 0.03802963194902986, + "peak_rss_kib": 16576 + }, + { + "wall_s": 0.03889605507720262, + "peak_rss_kib": 16652 + }, + { + "wall_s": 0.03797265503089875, + "peak_rss_kib": 16784 + }, + { + "wall_s": 0.03807136998511851, + "peak_rss_kib": 16620 + } + ], + "prediction": { + "start1": 39, + "end1": 45, + "start2": 98, + "end2": 104, + "E_kcal_mol": -5.57, + "hybridDB": "39|||||||&98|||||||", + "length_nt": 7 + }, + "wall_s": { + "median": 0.03807136998511851, + "min": 0.03797265503089875, + "max": 0.03889605507720262 + }, + "peak_rss_kib": { + "median": 16648, + "min": 16576, + "max": 16784 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNA", + "noGU": true, + "arguments": [ + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.05764365091454238, + "peak_rss_kib": 16660 + }, + { + "wall_s": 0.055861442000605166, + "peak_rss_kib": 16660 + }, + { + "wall_s": 0.0561534590087831, + "peak_rss_kib": 16968 + }, + { + "wall_s": 0.05309296608902514, + "peak_rss_kib": 16640 + }, + { + "wall_s": 0.055408544023521245, + "peak_rss_kib": 16836 + } + ], + "prediction": { + "start1": 39, + "end1": 45, + "start2": 98, + "end2": 104, + "E_kcal_mol": -5.57, + "hybridDB": "39|||||||&98|||||||", + "length_nt": 7 + }, + "wall_s": { + "median": 0.055861442000605166, + "min": 0.05309296608902514, + "max": 0.05764365091454238 + }, + "peak_rss_kib": { + "median": 16660, + "min": 16640, + "max": 16968 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAsnap", + "noGU": true, + "arguments": [ + "--personality=IntaRNAsnap", + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.038256115978583694, + "peak_rss_kib": 16912 + }, + { + "wall_s": 0.03771536401472986, + "peak_rss_kib": 16772 + }, + { + "wall_s": 0.0362653280608356, + "peak_rss_kib": 16960 + }, + { + "wall_s": 0.03812957904301584, + "peak_rss_kib": 16784 + }, + { + "wall_s": 0.038168090977706015, + "peak_rss_kib": 16796 + } + ], + "prediction": { + "start1": 39, + "end1": 45, + "start2": 98, + "end2": 104, + "E_kcal_mol": -5.57, + "hybridDB": "39|||||||&98|||||||", + "length_nt": 7 + }, + "wall_s": { + "median": 0.03812957904301584, + "min": 0.0362653280608356, + "max": 0.038256115978583694 + }, + "peak_rss_kib": { + "median": 16796, + "min": 16772, + "max": 16960 + }, + "structure_reevaluation_matches": true + } + ], + "comparisons": [ + { + "noGU": false, + "delta_E_kcal_mol": 0.02, + "delta_length_nt": -17, + "wall_ratio_kix_over_default": 0.6692527928459269, + "rss_ratio_kix_over_default": 0.9843897824030274 + }, + { + "noGU": true, + "delta_E_kcal_mol": 0.0, + "delta_length_nt": 0, + "wall_ratio_kix_over_default": 0.6825741992589945, + "rss_ratio_kix_over_default": 1.0081632653061225 + } + ] + }, + { + "name": "phoB/GcvB", + "target": { + "file": "handson/phoB.fasta", + "header": "phoB|NC_000913|b1130|-200..+100|genom-subsequence", + "length_nt": 299, + "sequence": "GAGCTATCACGATGGTTGATGAGCTGAAATAAACCTCGTATCAGTGCCGGATGGCGATGCTGTCCGGCCTGCTTATTAAGATTATCCGCTTTTTATTTTTTCACTTTACCTCCCCTCCCCGCTGGTTTATTTAATGTTTACCCCCATAACCACATAATCGCGTTACACTATTTTAATAATTAAGACAGGGAGAAATAAAAATGCGCGTACTGGTTGTTGAAGACAATGCGTTGTTACGTCACCACCTTAAAGTTCAGATTCAGGATGCTGGTCATCAGGTCGATGACGCAGAAGATGCC", + "sha256": "ce42dc25417d696c2402cc1daaf22ee5045e68c095fcd578b4072395a09ffa98" + }, + "query": { + "file": "handson/GcvB.fasta", + "header": "GcvB|NC_000913", + "length_nt": 201, + "sequence": "ACUUCCUGAGCCGGAACGAAAAGUUUUAUCGGAAUGCGUGUUCUGGUGAACUUUUGGCUUACGGUUGUGAUGUUGUGUUGUUGUGUUUGCAAUUGGUCUGCGAUUCAGACCAUGGUAGCAAAGCUACCUUUUUUCACUUCCUGUACAUUUACCCUGUCUGUCCAUAGUGAUUAAUGUAGCACCGCCUAAUUGCGGUGCUUU", + "sha256": "580180ec3b79266f0684c4c60d44b7b86ef6dbba087f67110d7c25c9ef31bcf7" + }, + "runs": [ + { + "program": "IntaRNA", + "noGU": false, + "arguments": [ + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 1.0826977379620075, + "peak_rss_kib": 18528 + }, + { + "wall_s": 1.0389898010762408, + "peak_rss_kib": 18660 + }, + { + "wall_s": 1.0598379949806258, + "peak_rss_kib": 18448 + }, + { + "wall_s": 1.0448508230037987, + "peak_rss_kib": 18436 + }, + { + "wall_s": 1.0702672540210187, + "peak_rss_kib": 18336 + } + ], + "prediction": { + "start1": 222, + "end1": 268, + "start2": 38, + "end2": 82, + "E_kcal_mol": -15.7, + "hybridDB": "222|||||||||.....|||||||.|||..||||||||...|||||||||&38|||||||||||||||||........|||.||||||||||||||||", + "length_nt": 47 + }, + "wall_s": { + "median": 1.0598379949806258, + "min": 1.0389898010762408, + "max": 1.0826977379620075 + }, + "peak_rss_kib": { + "median": 18448, + "min": 18336, + "max": 18660 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAsnap", + "noGU": false, + "arguments": [ + "--personality=IntaRNAsnap", + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 0.14911555196158588, + "peak_rss_kib": 18632 + }, + { + "wall_s": 0.1491208989173174, + "peak_rss_kib": 18456 + }, + { + "wall_s": 0.14401900197844952, + "peak_rss_kib": 18652 + }, + { + "wall_s": 0.15017399890348315, + "peak_rss_kib": 18356 + }, + { + "wall_s": 0.147206068970263, + "peak_rss_kib": 18068 + } + ], + "prediction": { + "start1": 183, + "end1": 190, + "start2": 152, + "end2": 159, + "E_kcal_mol": -13.19, + "hybridDB": "183||||||||&152||||||||", + "length_nt": 8 + }, + "wall_s": { + "median": 0.14911555196158588, + "min": 0.14401900197844952, + "max": 0.15017399890348315 + }, + "peak_rss_kib": { + "median": 18456, + "min": 18068, + "max": 18652 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNA", + "noGU": true, + "arguments": [ + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.7203351480420679, + "peak_rss_kib": 18328 + }, + { + "wall_s": 0.6980774619150907, + "peak_rss_kib": 18456 + }, + { + "wall_s": 0.7028307840228081, + "peak_rss_kib": 18396 + }, + { + "wall_s": 0.7158402650384232, + "peak_rss_kib": 18456 + }, + { + "wall_s": 0.7222581719979644, + "peak_rss_kib": 18508 + } + ], + "prediction": { + "start1": 183, + "end1": 190, + "start2": 152, + "end2": 159, + "E_kcal_mol": -13.19, + "hybridDB": "183||||||||&152||||||||", + "length_nt": 8 + }, + "wall_s": { + "median": 0.7158402650384232, + "min": 0.6980774619150907, + "max": 0.7222581719979644 + }, + "peak_rss_kib": { + "median": 18456, + "min": 18328, + "max": 18508 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAsnap", + "noGU": true, + "arguments": [ + "--personality=IntaRNAsnap", + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.14456413604784757, + "peak_rss_kib": 18192 + }, + { + "wall_s": 0.14236011693719774, + "peak_rss_kib": 18364 + }, + { + "wall_s": 0.1507548289373517, + "peak_rss_kib": 18660 + }, + { + "wall_s": 0.14155757101252675, + "peak_rss_kib": 18452 + }, + { + "wall_s": 0.15410517202690244, + "peak_rss_kib": 18656 + } + ], + "prediction": { + "start1": 183, + "end1": 190, + "start2": 152, + "end2": 159, + "E_kcal_mol": -13.19, + "hybridDB": "183||||||||&152||||||||", + "length_nt": 8 + }, + "wall_s": { + "median": 0.14456413604784757, + "min": 0.14155757101252675, + "max": 0.15410517202690244 + }, + "peak_rss_kib": { + "median": 18452, + "min": 18192, + "max": 18660 + }, + "structure_reevaluation_matches": true + } + ], + "comparisons": [ + { + "noGU": false, + "delta_E_kcal_mol": 2.51, + "delta_length_nt": -39, + "wall_ratio_kix_over_default": 0.1406965523672434, + "rss_ratio_kix_over_default": 1.0004336513443193 + }, + { + "noGU": true, + "delta_E_kcal_mol": 0.0, + "delta_length_nt": 0, + "wall_ratio_kix_over_default": 0.20195027174126337, + "rss_ratio_kix_over_default": 0.9997832683138275 + } + ] + }, + { + "name": "ilvE/GcvB.ST", + "target": { + "file": "handson/ilvE.fasta", + "header": "ilvE|NC_003197|STM3903|-200..+100|genom-subsequence", + "length_nt": 299, + "sequence": "GGTTTTCAGGTGTGCTCCATGAATATGGAAGCCGCGACCGATGCGCAGAATATAAATATTGAATTGACCGTTGCCAGTCCCCGGTCGGTCGACTTACTGTTTAGTCAGTTAAGTAAACTGGTAGATGTTGCGCATGTCGCGATCTGCCAGAGCGCTGCCACATCACAACAAATCCGCGCCTGAGCGCAAAAGGAAGAAAAATGACGACGAAAAAAGCTGATTATATTTGGTTCAATGGCGAGATGGTGCGCTGGGAAGACGCGAAGGTTCACGTAATGTCTCACGCGCTGCACTACGGT", + "sha256": "6fd5cab8545dd5cede4b8d2d23a276f7346a56706b3121de3bcf954d60a2068e" + }, + "query": { + "file": "handson/GcvB.ST.fasta", + "header": "GcvB|NC_003197", + "length_nt": 200, + "sequence": "ACUUCCUGAGCCGGAACGAAAAGUUUUAUCGGAAUGCGUGUUCUGAUGGGCUUUUGGCUUACGGUUGUGAUGUUGUGUUGUUGUGUUUGCAAUUGGUCUGCGAUUCAGACCACGGUAGCGAGACUACCCUUUUUCACUUCCUGUACAUUUACCCUGUCUGUCCAUAGUGAUUAAUGUAGCACCGCCAUAUUGCGGUGCUU", + "sha256": "529466282586ededdea0ed04199197f04978e2be14245feed0655bc6850206bb" + }, + "runs": [ + { + "program": "IntaRNA", + "noGU": false, + "arguments": [ + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 1.3832573930267245, + "peak_rss_kib": 18452 + }, + { + "wall_s": 1.286010636948049, + "peak_rss_kib": 18452 + }, + { + "wall_s": 1.3232459670398384, + "peak_rss_kib": 18064 + }, + { + "wall_s": 1.2741687439847738, + "peak_rss_kib": 18440 + }, + { + "wall_s": 1.2935187839902937, + "peak_rss_kib": 18336 + } + ], + "prediction": { + "start1": 33, + "end1": 81, + "start2": 48, + "end2": 102, + "E_kcal_mol": -14.24, + "hybridDB": "33||||||||||..||||||||..|||||...|||||||||.|||||.|||&48|||...|||||..||||||..|||...|||||..|||.|||||.|||||||.|||", + "length_nt": 55 + }, + "wall_s": { + "median": 1.2935187839902937, + "min": 1.2741687439847738, + "max": 1.3832573930267245 + }, + "peak_rss_kib": { + "median": 18440, + "min": 18064, + "max": 18452 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAsnap", + "noGU": false, + "arguments": [ + "--personality=IntaRNAsnap", + "--outNoGUend=false" + ], + "samples": [ + { + "wall_s": 0.14518713497091085, + "peak_rss_kib": 18452 + }, + { + "wall_s": 0.1470703890081495, + "peak_rss_kib": 18512 + }, + { + "wall_s": 0.14972956804558635, + "peak_rss_kib": 18512 + }, + { + "wall_s": 0.1487087750574574, + "peak_rss_kib": 18772 + }, + { + "wall_s": 0.1428277890663594, + "peak_rss_kib": 18452 + } + ], + "prediction": { + "start1": 157, + "end1": 169, + "start2": 64, + "end2": 76, + "E_kcal_mol": -9.84, + "hybridDB": "157||.||||||||||&64||||||||||.||", + "length_nt": 13 + }, + "wall_s": { + "median": 0.1470703890081495, + "min": 0.1428277890663594, + "max": 0.14972956804558635 + }, + "peak_rss_kib": { + "median": 18512, + "min": 18452, + "max": 18772 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNA", + "noGU": true, + "arguments": [ + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 1.0233468420337886, + "peak_rss_kib": 18588 + }, + { + "wall_s": 1.0323557329829782, + "peak_rss_kib": 18448 + }, + { + "wall_s": 1.0218827379867435, + "peak_rss_kib": 18648 + }, + { + "wall_s": 1.0279411339433864, + "peak_rss_kib": 18560 + }, + { + "wall_s": 0.9871100430609658, + "peak_rss_kib": 18520 + } + ], + "prediction": { + "start1": 131, + "end1": 170, + "start2": 39, + "end2": 77, + "E_kcal_mol": -10.55, + "hybridDB": "131||||||||||||||.||||||.....|||.|||||.||||&39||||.||||||||.||||||....|||||||||||.|||", + "length_nt": 40 + }, + "wall_s": { + "median": 1.0233468420337886, + "min": 0.9871100430609658, + "max": 1.0323557329829782 + }, + "peak_rss_kib": { + "median": 18560, + "min": 18448, + "max": 18648 + }, + "structure_reevaluation_matches": true + }, + { + "program": "IntaRNAsnap", + "noGU": true, + "arguments": [ + "--personality=IntaRNAsnap", + "--outNoGUend=true" + ], + "samples": [ + { + "wall_s": 0.14092493208590895, + "peak_rss_kib": 18624 + }, + { + "wall_s": 0.14769695396535099, + "peak_rss_kib": 18432 + }, + { + "wall_s": 0.143615689012222, + "peak_rss_kib": 18500 + }, + { + "wall_s": 0.14164107700344175, + "peak_rss_kib": 18388 + }, + { + "wall_s": 0.1450244919396937, + "peak_rss_kib": 18300 + } + ], + "prediction": { + "start1": 160, + "end1": 169, + "start2": 64, + "end2": 73, + "E_kcal_mol": -9.13, + "hybridDB": "160||||||||||&64||||||||||", + "length_nt": 10 + }, + "wall_s": { + "median": 0.143615689012222, + "min": 0.14092493208590895, + "max": 0.14769695396535099 + }, + "peak_rss_kib": { + "median": 18432, + "min": 18300, + "max": 18624 + }, + "structure_reevaluation_matches": true + } + ], + "comparisons": [ + { + "noGU": false, + "delta_E_kcal_mol": 4.4, + "delta_length_nt": -42, + "wall_ratio_kix_over_default": 0.11369791519722769, + "rss_ratio_kix_over_default": 1.0039045553145336 + }, + { + "noGU": true, + "delta_E_kcal_mol": 1.42, + "delta_length_nt": -30, + "wall_ratio_kix_over_default": 0.14033921160767127, + "rss_ratio_kix_over_default": 0.993103448275862 + } + ] + } + ] +} diff --git a/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg b/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg new file mode 100644 index 00000000..a79b875f --- /dev/null +++ b/doc/recursions/IntaRNAkix.PredictorSeedExtensionKinetic.svg @@ -0,0 +1,86 @@ + + IntaRNAsnap: deterministic kinetic seed-extension recurrence + For each handler-provided seed, enumerate single stacks, double stacks and loop-plus-stack moves at both ends. Recompute complete interaction energies, choose the best strictly downhill move, and retain reportable visited prefixes. Stop when no downhill move exists. + + + + + + + + + + + + + + + + + + + + + + + + + + IntaRNAsnap + Deterministic kinetic seed extension · model X, mode K · noLP extensions + + 1. Initialize one trajectory per seed + + seed + + S₀ = seed; H₀ = Hseed + Einit; E₀ = E(b₀, H₀) + + Trust the seed handler’s pairs and energy, including explicit lonely pairs. + Sₜ = actual pair chain; bₜ = (i₁, j₁, i₂, j₂); Hₜ = hybridization energy including initiation. + Indices of RNA 2 are reversed; left/right refer to the energy-index coordinate system. + + 2. Enumerate feasible extensions at both ends + + Single stack + Double stack + Loop + outward stack + + Left + Sₜ + Sₜ + Sₜ + Right + + + + + SₜSₜSₜ + + + 1 new pair; s₁ = s₂ = 0 + 2 new pairs; s₁ = s₂ = 0 + 2 new pairs; s₁ + s₂ > 0 + + sₖ = skipped bases on RNA k; 0 ≤ sₖ ≤ mₖ (its loop limit). Two-pair moves are atomic. + Enforce pairing, GU-loop policy, prediction ranges and both span limits. No isolated loop closure. + New boundaries advance by sₖ + 1 or sₖ + 2; the opposite boundaries stay fixed. + + 3. Recompute full energies and choose the best downhill move + + H′ = Hₜ + Eloop + Eoutward stack, if two pairs + ΔE(c) = E(b′, H′) − Eₜ; C₋ = { feasible finite moves c : ΔE(c) < 0 } + c* = arg minc ∈ C₋ ΔE(c) / D(c) + + D = 1 (score A, default), 1 + s₁ + s₂ (B), or 1 + 2 max(s₁, s₂) (C). + Ties: left end, smaller s₁ + s₂, smaller s₁, then single pair. Compare exact ratios. + E includes ED₁, ED₂, terminal terms, both weighted dangles and the additive term. + Refresh full energies at both ends, even when pairing and local energies are cached. + + 4. Commit, retain and repeat + |C₋| > 0: (Sₜ₊₁, Hₜ₊₁, Eₜ₊₁) = (Sₜ extended by c*, H′, E(b′, H′)). + Retain reportable visited states, including the seed; repeat step 2. If C₋ is empty, stop. + Deduplicate boundaries, rank by full E, apply output constraints and restore the actual pair chain. + + Greedy trajectory heuristic: no global-MFE guarantee, transition rates or physical time axis. + + diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index 98288166..5f27f746 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -72,6 +72,7 @@ libIntaRNA_a_HEADERS = \ PredictorMfe2dSeed.h \ PredictorMfe2dSeedExtension.h \ PredictorMfe2dSeedExtensionRIblast.h \ + PredictorSeedExtensionKinetic.h \ PredictorMfe2dHeuristic.h \ PredictorMfe2dHeuristicSeed.h \ PredictorMfe2dHelixBlockHeuristic.h \ @@ -134,6 +135,7 @@ libIntaRNA_a_SOURCES = \ PredictorMfe2dSeed.cpp \ PredictorMfe2dSeedExtension.cpp \ PredictorMfe2dSeedExtensionRIblast.cpp \ + PredictorSeedExtensionKinetic.cpp \ PredictorMfe2dHeuristic.cpp \ PredictorMfe2dHeuristicSeed.cpp \ PredictorMfe2dHelixBlockHeuristic.cpp \ diff --git a/src/IntaRNA/PredictorSeedExtensionKinetic.cpp b/src/IntaRNA/PredictorSeedExtensionKinetic.cpp new file mode 100644 index 00000000..dba6cd3a --- /dev/null +++ b/src/IntaRNA/PredictorSeedExtensionKinetic.cpp @@ -0,0 +1,389 @@ +#include "IntaRNA/PredictorSeedExtensionKinetic.h" + +#include +#include +#include +#include + +#include + +namespace IntaRNA { +namespace { + +SeedHandler * checkedSeedHandler(SeedHandler * handler) +{ + if (handler == NULL) { + throw std::invalid_argument("PredictorSeedExtensionKinetic requires a seed handler"); + } + return handler; +} + +// Do not add infinity sentinels or overflow the internal integer energy type. +E_type addEnergy(const E_type first, const E_type second) +{ + if (E_isINF(first) || E_isINF(second)) { + return E_INF; + } + const std::int64_t sum = std::int64_t(first) + std::int64_t(second); + return sum >= E_INF || sum < std::numeric_limits::min() + ? E_INF : static_cast(sum); +} + +} // namespace + +////////////////////////////////////////////////////////////////////////// + +PredictorSeedExtensionKinetic::PredictorSeedExtensionKinetic( + const InteractionEnergy & energy, OutputHandler & output, + PredictionTracker * predTracker, SeedHandler * seedHandlerInstance, + const char score) + : PredictorMfe(energy, output, predTracker) + , seedHandler(checkedSeedHandler(seedHandlerInstance)) + , score(score) + , interactions() +{ + if (score != 'A' && score != 'B' && score != 'C') { + throw std::invalid_argument("PredictorSeedExtensionKinetic score must be A, B or C"); + } + if (output.getOutputConstraint().needZall) { + throw std::invalid_argument("PredictorSeedExtensionKinetic does not compute an equilibrium partition function"); + } +} + +////////////////////////////////////////////////////////////////////////// + +PredictorSeedExtensionKinetic::~PredictorSeedExtensionKinetic() +{ +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::predict(const IndexRange & r1, const IndexRange & r2) +{ + const size_t size1 = energy.getAccessibility1().getSequence().size(); + const size_t size2 = energy.getAccessibility2().getSequence().size(); + if (!r1.isAscending() || !r2.isAscending() + || r1.from >= size1 || r2.from >= size2) { + throw std::invalid_argument("PredictorSeedExtensionKinetic::predict(): invalid sequence range"); + } + + energy.setOffset1(r1.from); + energy.setOffset2(r2.from); + seedHandler.setOffset1(r1.from); + seedHandler.setOffset2(r2.from); + const size_t last1 = std::min(r1.to, size1 - 1) - r1.from; + const size_t last2 = std::min(r2.to, size2 - 1) - r2.from; + interactions.clear(); + initOptima(); + + if (seedHandler.fillSeed(0, last1, 0, last2) != 0) { + size_t i1 = RnaSequence::lastPos, i2 = RnaSequence::lastPos; + while (seedHandler.updateToNextSeed(i1, i2, 0, last1, 0, last2)) { + const size_t length1 = seedHandler.getSeedLength1(i1, i2); + const size_t length2 = seedHandler.getSeedLength2(i1, i2); + if (length1 == 0 || length2 == 0 + || length1 - 1 > last1 - i1 || length2 - 1 > last2 - i2 + || length1 > energy.getAccessibility1().getMaxLength() + || length2 > energy.getAccessibility2().getMaxLength()) { + continue; + } + const size_t j1 = i1 + length1 - 1; + const size_t j2 = i2 + length2 - 1; + if (E_isINF(seedHandler.getSeedE(i1, i2))) { + continue; + } + + Interaction interaction(energy.getAccessibility1().getSequence(), + energy.getAccessibility2().getAccessibilityOrigin().getSequence()); + interaction.basePairs.push_back(energy.getBasePair(i1, i2)); + seedHandler.traceBackSeed(interaction, i1, i2); + if (i1 != j1 || i2 != j2) { + interaction.basePairs.push_back(energy.getBasePair(j1, j2)); + } + interaction.sort(); + const E_type hybrid = addEnergy(seedHandler.getSeedE(i1, i2), energy.getE_init()); + interaction.energy = energy.getE(i1, j1, i2, j2, hybrid); + if (E_isINF(interaction.energy)) { + continue; + } + interaction.setSeedRange(interaction.basePairs.front(), + interaction.basePairs.back(), interaction.energy); + extendSeed(interaction, hybrid, last1, last2); + } + } + + // Reduce identical boundaries before updating optima: an earlier, inferior + // path must never be traced back using a later replacement's base pairs. + for (const auto & entry : interactions) { + const Boundary & b = entry.first; + updateOptima(b[0], b[1], b[2], b[3], entry.second.energy, false, false); + } + // The generic reporter assumes a nonempty optimum list. Zero reports can + // still be useful for prediction trackers and must not dereference it. + if (output.getOutputConstraint().reportMax != 0) { + reportOptima(); + } +} + +////////////////////////////////////////////////////////////////////////// + +PredictorSeedExtensionKinetic::Boundary +PredictorSeedExtensionKinetic::getBoundary(const Interaction & interaction) const +{ + return Boundary{energy.getIndex1(interaction.basePairs.front()), + energy.getIndex1(interaction.basePairs.back()), + energy.getIndex2(interaction.basePairs.front()), + energy.getIndex2(interaction.basePairs.back())}; +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::extendSeed(Interaction & interaction, + E_type hybrid, const size_t last1, const size_t last2) +{ + std::array sides; + Boundary bounds = getBoundary(interaction); + buildCandidates(sides[0], bounds, true, last1, last2); + buildCandidates(sides[1], bounds, false, last1, last2); + while (true) { + retain(interaction); + const Candidate * left = updateCandidates(sides[0], bounds, hybrid, interaction.energy); + const Candidate * right = updateCandidates(sides[1], bounds, hybrid, interaction.energy); + if (left == NULL && right == NULL) { + break; + } + const Candidate best = left != NULL && (right == NULL || isBetter(*left, *right)) ? *left : *right; + const Interaction::BasePair close = energy.getBasePair(best.close1, best.close2); + if (best.left) { + interaction.basePairs.insert(interaction.basePairs.begin(), close); + if (best.macro) { + interaction.basePairs.insert(interaction.basePairs.begin(), + energy.getBasePair(best.bounds[0], best.bounds[2])); + } + } else { + interaction.basePairs.push_back(close); + if (best.macro) { + interaction.basePairs.push_back(energy.getBasePair(best.bounds[1], best.bounds[3])); + } + } + hybrid = best.hybrid; + interaction.energy = best.total; + bounds = best.bounds; + // The opposite end keeps its geometry, pair checks and loop energies. + // Its full energy must still be refreshed (ED and BOTH dangles change). + buildCandidates(sides[best.left ? 0 : 1], bounds, best.left, last1, last2); + } +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::buildCandidates(SideCandidates & side, + const Boundary & bounds, const bool left, const size_t last1, const size_t last2) const +{ + side.moves.clear(); + const size_t space1 = std::min(energy.getAccessibility1().getMaxLength() - (bounds[1]-bounds[0]+1), + left ? bounds[0] : last1-bounds[1]); + const size_t space2 = std::min(energy.getAccessibility2().getMaxLength() - (bounds[3]-bounds[2]+1), + left ? bounds[2] : last2-bounds[3]); + if (space1 == 0 || space2 == 0) { + return; + } + const bool stackOnly = (output.getOutputConstraint().noGUend || !energy.isInternalLoopGUallowed()) + && energy.isGU(bounds[left ? 0 : 1], bounds[left ? 2 : 3]); + const size_t maxGap1 = space1 < 2 || stackOnly ? 0 : std::min(energy.getMaxInternalLoopSize1(), space1-2); + const size_t maxGap2 = space2 < 2 || stackOnly ? 0 : std::min(energy.getMaxInternalLoopSize2(), space2-2); + side.columns = maxGap2+2; + side.complementary.assign((maxGap1+2)*side.columns, -1); + const auto append = [&](size_t s1, size_t s2, bool macro) { + Candidate c; + c.left = left; c.s1 = s1; c.s2 = s2; c.macro = macro; + c.bounds = bounds; + const size_t pairs = macro ? 2 : 1; + if (left) { + c.close1 = bounds[0]-s1-1; c.close2 = bounds[2]-s2-1; + c.bounds[0] -= s1+pairs; c.bounds[2] -= s2+pairs; + } else { + c.close1 = bounds[1]+s1+1; c.close2 = bounds[3]+s2+1; + c.bounds[1] += s1+pairs; c.bounds[3] += s2+pairs; + } + side.moves.push_back(c); + }; + append(0, 0, false); + if (space1 >= 2 && space2 >= 2) { + for (size_t s1 = 0; s1 <= maxGap1; ++s1) { + for (size_t s2 = 0; s2 <= maxGap2; ++s2) { + append(s1, s2, true); + } + } + } +} + +////////////////////////////////////////////////////////////////////////// + +const PredictorSeedExtensionKinetic::Candidate * +PredictorSeedExtensionKinetic::updateCandidates(SideCandidates & side, + const Boundary & bounds, const E_type hybrid, const E_type total) const +{ + // Phase one: each position pair is tested at most once per unchanged end, + // even when it is the closing pair of one move and outer pair of another. + for (Candidate & c : side.moves) { + c.bounds[c.left ? 1 : 0] = bounds[c.left ? 1 : 0]; + c.bounds[c.left ? 3 : 2] = bounds[c.left ? 3 : 2]; + c.active = c.bounds[1]-c.bounds[0]+1 <= energy.getAccessibility1().getMaxLength() + && c.bounds[3]-c.bounds[2]+1 <= energy.getAccessibility2().getMaxLength(); + if (!c.active || c.topologyKnown) { + continue; + } + const auto complementary = [&](size_t s1, size_t s2) { + signed char & cached = side.complementary[s1*side.columns+s2]; + if (cached < 0) { + cached = energy.areComplementary(c.left ? bounds[0]-s1-1 : bounds[1]+s1+1, + c.left ? bounds[2]-s2-1 : bounds[3]+s2+1); + } + return cached != 0; + }; + c.topologyKnown = true; + c.topologyAllowed = complementary(c.s1, c.s2) + && (!c.macro || complementary(c.s1+1, c.s2+1)); + if (c.topologyAllowed && output.getOutputConstraint().noGUend && (c.s1 != 0 || c.s2 != 0)) { + c.topologyAllowed = !energy.isGU(c.close1, c.close2); + } + } + const Candidate * best = NULL; + for (Candidate & c : side.moves) { + if (!c.active || !c.topologyAllowed) { + continue; + } + if (!c.localKnown) { + c.localKnown = true; + c.local = c.left ? energy.getE_interLeft(c.close1, bounds[0], c.close2, bounds[2]) + : energy.getE_interLeft(bounds[1], c.close1, bounds[3], c.close2); + if (c.macro && E_isNotINF(c.local)) { + c.local = addEnergy(c.local, c.left + ? energy.getE_interLeft(c.bounds[0], c.close1, c.bounds[2], c.close2) + : energy.getE_interLeft(c.close1, c.bounds[1], c.close2, c.bounds[3])); + } + } + c.hybrid = addEnergy(hybrid, c.local); + if (E_isINF(c.hybrid)) { + continue; + } + c.total = energy.getE(c.bounds[0], c.bounds[1], c.bounds[2], c.bounds[3], c.hybrid); + if (E_isINF(c.total)) { + continue; + } + c.delta = std::int64_t(c.total)-std::int64_t(total); + if (c.delta < 0 && (best == NULL || isBetter(c, *best))) { + best = &c; + } + } + return best; +} + +////////////////////////////////////////////////////////////////////////// + +bool +PredictorSeedExtensionKinetic::isBetter(const Candidate & candidate, const Candidate & best) const +{ + // A 32-bit energy difference times a 65-bit gap denominator fits in + // 128 bits, including for public-API loop limits beyond the CLI limits. + using Wide = boost::multiprecision::int128_t; + const auto denominator = [this](const Candidate & c) -> Wide { + if (score == 'B') { + return Wide(1) + Wide(c.s1) + Wide(c.s2); + } + if (score == 'C') { + return Wide(1) + 2 * Wide(std::max(c.s1, c.s2)); + } + return Wide(1); + }; + const Wide lhs = Wide(candidate.delta) * denominator(best); + const Wide rhs = Wide(best.delta) * denominator(candidate); + if (lhs != rhs) { + return lhs < rhs; + } + if (candidate.left != best.left) { + return candidate.left; + } + const Wide size = Wide(candidate.s1) + Wide(candidate.s2); + const Wide bestSize = Wide(best.s1) + Wide(best.s2); + if (size != bestSize) return size < bestSize; + if (candidate.s1 != best.s1) return candidate.s1 < best.s1; + // Identical shape/score: retain the shorter move first. + return !candidate.macro && best.macro; +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::retain(const Interaction & interaction) +{ + const Boundary b = getBoundary(interaction); + const auto & constraint = output.getOutputConstraint(); + if (interaction.energy >= E_MAX + || (constraint.noGUend && (energy.isGU(b[0], b[2]) || energy.isGU(b[1], b[3]))) + || energy.getED1(b[0], b[1]) > constraint.maxED + || energy.getED2(b[2], b[3]) > constraint.maxED) { + return; + } + auto existing = interactions.find(b); + if (existing == interactions.end()) { + interactions.emplace(b, interaction); + } else if (interaction.energy < existing->second.energy + || (interaction.energy == existing->second.energy + && interaction.basePairs < existing->second.basePairs)) { + existing->second = interaction; + } +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::traceBack(Interaction & interaction) +{ + if (interaction.basePairs.empty()) { + return; + } + const auto path = interactions.find(getBoundary(interaction)); + if (path == interactions.end() || path->second.energy != interaction.energy) { + throw std::runtime_error("PredictorSeedExtensionKinetic::traceBack(): no matching greedy path"); + } + interaction = path->second; + seedHandler.addSeeds(interaction); +} + +////////////////////////////////////////////////////////////////////////// + +void +PredictorSeedExtensionKinetic::getNextBest(Interaction & interaction) +{ + const E_type previousEnergy = interaction.energy; + const Interaction * best = NULL; + for (const auto & entry : interactions) { + const Boundary & b = entry.first; + const Interaction & candidate = entry.second; + if (candidate.energy < previousEnergy + || reportedInteractions.first.overlaps(IndexRange(b[0], b[1])) + || reportedInteractions.second.overlaps(IndexRange(b[2], b[3]))) { + continue; + } + if (best == NULL || candidate < *best) { + best = &candidate; + } + } + if (best == NULL) { + interaction.clear(); + interaction.energy = E_INF; + return; + } + interaction = *best; + const Interaction::BasePair right = interaction.basePairs.back(); + interaction.basePairs.resize(interaction.basePairs.size() == 1 ? 1 : 2); + interaction.basePairs.back() = right; + INTARNA_CLEANUP(interaction.seed); +} + +} // namespace IntaRNA diff --git a/src/IntaRNA/PredictorSeedExtensionKinetic.h b/src/IntaRNA/PredictorSeedExtensionKinetic.h new file mode 100644 index 00000000..fd996428 --- /dev/null +++ b/src/IntaRNA/PredictorSeedExtensionKinetic.h @@ -0,0 +1,166 @@ +#ifndef INTARNA_PREDICTORSEEDEXTENSIONKINETIC_H_ +#define INTARNA_PREDICTORSEEDEXTENSIONKINETIC_H_ + +#include "IntaRNA/PredictorMfe.h" +#include "IntaRNA/SeedHandlerIdxOffset.h" + +#include +#include +#include +#include + +namespace IntaRNA { + +/** + * Deterministic, strictly downhill extension of every feasible seed. + * + * At each step both ends compete using the complete interaction-energy + * difference, including accessibility, terminal penalties and both weighted + * dangling ends. Score A uses this difference directly; B divides by + * 1+s1+s2; C (C1) divides by 1+2*max(s1,s2). These are heuristic move rankings, + * not physical rates or a calibrated folding-time model. Only negative + * energy differences are accepted, including for the normalized scores. + * + * Extensions always add one stacked pair or two stacked pairs, including + * across a loop. Seeds and their energies are trusted as supplied by the seed + * handler, even when an explicit seed contains lonely pairs. Ties prefer left, + * smaller s1+s2, smaller s1, then the single-pair move. Every valid + * visited prefix is eligible for normal MFE/suboptimal reporting; traceback + * reproduces the actual greedy path. Equilibrium partition-function output + * is unsupported. + * + * Candidate enumeration uses the active energy model and separate loop/span + * limits for both RNAs. Every feasible move is evaluated with its complete + * energy change. + */ +class PredictorSeedExtensionKinetic : public PredictorMfe { +public: + + /** + * Constructs a predictor, taking ownership of tracker and seed handler. + * @param energy energy model, which must outlive this predictor + * @param output output handler, which must outlive this predictor + * @param predTracker owned tracker, or NULL + * @param seedHandler owned, non-NULL seed handler + * @param score move ranking: A, B, or C (the C1 formula) + * @throws std::invalid_argument for a NULL seed handler, unknown score or + * output requiring an equilibrium partition function + */ + PredictorSeedExtensionKinetic(const InteractionEnergy & energy, + OutputHandler & output, PredictionTracker * predTracker, + SeedHandler * seedHandler, const char score = 'A'); + + /** Frees the owned seed handler and prediction tracker. */ + virtual ~PredictorSeedExtensionKinetic(); + + /** + * Extends all feasible seeds within inclusive, zero-based sequence ranges. + * Sequence 2 uses the reversed indexing of the energy model. Repeated + * calls reset all trajectories, cached output and index offsets. + * @param r1 permitted range in sequence 1 + * @param r2 permitted range in reversed sequence 2 + * @throws std::invalid_argument for an invalid or empty input range + */ + void predict(const IndexRange & r1 = IndexRange(0, RnaSequence::lastPos), + const IndexRange & r2 = IndexRange(0, RnaSequence::lastPos)) override; + +protected: + /** + * Restores the exact stored greedy path and annotates its contained seeds. + * @param interaction interaction boundaries to expand + */ + void traceBack(Interaction & interaction) override; + + /** + * Finds the best cached prefix disjoint from already reported intervals. + * @param interaction current report, replaced by the next report or E_INF + */ + void getNextBest(Interaction & interaction) override; + +private: + //! Inclusive boundaries (i1,j1,i2,j2), using local energy indices. + using Boundary = std::array; + //! Best actual path for each visited, reportable set of boundaries. + using InteractionCache = std::map; + + /** A feasible move; delta is widened before subtraction. */ + struct Candidate { + Boundary bounds = {}; + size_t close1 = 0; + size_t close2 = 0; + size_t s1 = 0; + size_t s2 = 0; + bool left = true; + bool macro = false; + bool topologyKnown = false; + bool topologyAllowed = false; + bool active = false; + bool localKnown = false; + E_type local = E_INF; + E_type hybrid = E_INF; + E_type total = E_INF; + std::int64_t delta = 0; + }; + + /** Geometry, shared pair checks and local energies for one unchanged end. */ + struct SideCandidates { + std::vector moves; + std::vector complementary; + size_t columns = 0; + }; + + //! Owned seed handler with offsets matching this->energy. + SeedHandlerIdxOffset seedHandler; + //! The selected deterministic move-ranking formula. + const char score; + //! Paths retained independently of the number of requested reports. + InteractionCache interactions; + /** @return local, inclusive boundaries of a nonempty interaction */ + Boundary getBoundary(const Interaction & interaction) const; + + /** + * Runs a complete greedy trajectory and retains its reportable prefixes. + * @param interaction initial seed, modified to its final state + * @param hybrid initial seed hybridization energy including initiation + * @param last1 last permitted local index in sequence 1 + * @param last2 last permitted local index in reversed sequence 2 + */ + void extendSeed(Interaction & interaction, E_type hybrid, + size_t last1, size_t last2); + + /** + * Builds the rectangular gap table plus the single-stack move for one end. + * @param side table to replace + * @param bounds current boundaries + * @param left whether this is the left end + * @param last1 last permitted local index in sequence 1 + * @param last2 last permitted local index in reversed sequence 2 + */ + void buildCandidates(SideCandidates & side, const Boundary & bounds, + bool left, size_t last1, size_t last2) const; + + /** + * Resolves shared pair checks before energies, then refreshes complete + * energies and selects the best downhill candidate during that traversal. + * @param side candidate and pairing cache for one end + * @param bounds current boundaries (the opposite end may have changed) + * @param hybrid current hybridization energy including initiation + * @param total current complete interaction energy + * @return best move within side, or NULL if none is downhill + */ + const Candidate * updateCandidates(SideCandidates & side, const Boundary & bounds, + E_type hybrid, E_type total) const; + + /** @return whether a move wins by exact score and deterministic ties */ + bool isBetter(const Candidate & candidate, const Candidate & best) const; + + /** + * Reduces a valid prefix by boundaries, total energy and full-path ties. + * @param interaction actual path to consider for reporting + */ + void retain(const Interaction & interaction); +}; + +} // namespace IntaRNA + +#endif /* INTARNA_PREDICTORSEEDEXTENSIONKINETIC_H_ */ diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 4478868a..c06ac8b0 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -50,6 +50,7 @@ extern "C" { #include "IntaRNA/PredictorMfe2dSeed.h" #include "IntaRNA/PredictorMfe2dSeedExtension.h" #include "IntaRNA/PredictorMfe2dSeedExtensionRIblast.h" +#include "IntaRNA/PredictorSeedExtensionKinetic.h" #include "IntaRNA/PredictorMfe2dHeuristicSeedExtension.h" #include "IntaRNA/PredictorMfeEnsSeedOnly.h" @@ -183,7 +184,8 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) temperature("temperature",0,100,37), model("model", "SPBX", 'X'), - mode("mode", "HMSR", 'H'), // R for RIblast heuristic only + mode("mode", "HMSRK", 'H'), // R for RIblast heuristic only + kineticScore("kineticScore", "ABC", 'A'), #if INTARNA_MULITHREADING threads("threads", 0, omp_get_max_threads(), 1), #endif @@ -330,6 +332,12 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) resetParamDefault<>(accL, 0); resetParamDefault<>(intLenMax, 60); break; + case IntaRNAsnap : + // deterministic kinetic seed extension + resetParamDefault<>(model, 'X'); + resetParamDefault<>(mode, 'K'); + resetParamDefault<>(outNoLP, true, "outNoLP"); + break; case IntaRNAseed : // seed-only prediction resetParamDefault<>(mode, 'S'); @@ -803,8 +811,16 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) , std::string("prediction mode : " "\n 'H' = heuristic (fast and low memory), " "\n 'M' = exact (slow), " - "\n 'S' = seed-only" + "\n 'S' = seed-only, " + "\n 'K' = downhill greedy seed extension (IntaRNAsnap; requires --model=X; always noLP extensions)" ).c_str()) + (kineticScore.name.c_str() + , value(&(kineticScore.val)) + ->default_value(kineticScore.def) + ->notifier(boost::bind(&CommandLineParsing::validate_charArgument,this,kineticScore,_1)) + , "candidate score for --model=X --mode=K: 'A' = complete interaction energy change, " + "'B' = change/(1+s1+s2), 'C' = change/(1+2*max(s1,s2)), where s1/s2 are skipped bases. " + "All modes accept strictly negative energy changes only; B/C are heuristic scores.") (model.name.c_str() , value(&(model.val)) ->default_value(model.def) @@ -1225,6 +1241,28 @@ parse(int argc, char** argv) // parsing escape literals outSep = unescaped_string::getUnescaped( outSep ); + // K needs a seed even before the usual --noSeed model normalization. + if (mode.val == 'K') { + if (model.val != 'X') throw error("--mode=K is available only with --model=X"); + if (noSeedRequired) throw error("--mode=K requires seeds and is incompatible with --noSeed"); + if (!outNoLP) { + LOG(INFO) << "--mode=" << mode.val << " uses no-lonely-pair extensions: setting --outNoLP=true (handler-provided seeds are unchanged)"; + outNoLP = true; + } + // A selected set of greedy paths is not an equilibrium ensemble. + if (outMode.val == 'E' + || (outMode.val == 'C' && OutputHandlerCsv::needsZall(OutputHandlerCsv::string2list(outCsvCols))) + || !outPrefix2streamName.at(OutPrefixCode::OP_spotProb).empty() + || !outPrefix2streamName.at(OutPrefixCode::OP_spotProbAll).empty() + || !outPrefix2streamName.at(OutPrefixCode::OP_qSpotProb).empty() + || !outPrefix2streamName.at(OutPrefixCode::OP_tSpotProb).empty()) + { + throw error("--mode=K does not support equilibrium ensemble or interaction-probability output"); + } + } else if (vm.count(kineticScore.name) && !vm.at(kineticScore.name).defaulted()) { + throw error("--kineticScore requires --model=X --mode=K"); + } + // open output stream { // open according stream @@ -1675,7 +1713,7 @@ CommandLineParsing::prepareEvaluation( boost::program_options::variables_map & v // Erase before notify(), so even out-of-sequence regions and seed encodings // cannot constrain or invalidate evaluation. Energy/accessibility options stay. const std::set ignored = { - "model", "mode", "noSeed", "intLenMax", "qIntLenMax", "tIntLenMax", + "model", "mode", "kineticScore", "noSeed", "intLenMax", "qIntLenMax", "tIntLenMax", "intLoopMax", "qIntLoopMax", "tIntLoopMax", "qRegion", "tRegion", "qRegionLenMax", "tRegionLenMax", "windowWidth", "windowOverlap", "outNumber", "outOverlap", "outMaxE", "outDeltaE", "outMinPu", @@ -2580,6 +2618,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output ) const case 'H' : return new PredictorMfe2dHeuristicSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); case 'M' : return new PredictorMfe2dSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); case 'R' : return new PredictorMfe2dSeedExtensionRIblast( energy, output, predTracker, getSeedHandler( energy ) ); + case 'K' : return new PredictorSeedExtensionKinetic( energy, output, predTracker, getSeedHandler( energy ), kineticScore.val ); case 'S' : return new PredictorMfeSeedOnly( energy, output, predTracker, getSeedHandler( energy ) ); default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not implemented"); return NULL; } @@ -2929,6 +2968,9 @@ getPersonality( int argc, char ** argv ) } // parse personality + if (value == "IntaRNAsnap") { + return Personality::IntaRNAsnap; + } if (value == "IntaRNAeval") { return Personality::IntaRNAeval; } diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index 0deec81e..f6658728 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -60,6 +60,7 @@ class CommandLineParsing { IntaRNA3, // default IntaRNA v3 setup IntaRNAens, // ensemble-based prediction IntaRNAeval, // evaluate predefined interactions + IntaRNAsnap, // kinetic seed extension IntaRNAsTar, // sRNA-target prediction (optimized parameter) IntaRNAseed, // seed-only predictions IntaRNAhelix, // helix-block-based predictions @@ -84,6 +85,7 @@ class CommandLineParsing { case IntaRNA3 : return "IntaRNA3"; case IntaRNAens : return "IntaRNAens"; case IntaRNAeval : return "IntaRNAeval"; + case IntaRNAsnap : return "IntaRNAsnap"; case IntaRNAsTar : return "IntaRNAsTar"; case IntaRNAseed : return "IntaRNAseed"; case IntaRNAhelix : return "IntaRNAhelix"; @@ -721,6 +723,8 @@ class CommandLineParsing { CharParameter model; //! the prediction mode (heuristic, space-efficient, exact) CharParameter mode; + //! greedy seed-extension score: A=energy, B=Manhattan, C=asymmetry + CharParameter kineticScore; #if INTARNA_MULITHREADING //! number of threads = number of parallel predictors running NumberParameter threads; diff --git a/tests/Makefile.am b/tests/Makefile.am index 6f580ca6..ea9662b2 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 runOutputOverlap.sh +dist_check_SCRIPTS = runIntaRNA.sh runAccessibilityBinary.sh runIntaRNAeval.sh runKineticSeedExtension.sh runOutputOverlap.sh # the program to build check_PROGRAMS = runApiTests @@ -52,6 +52,7 @@ runApiTests_SOURCES = \ PredictorTinyOracle_test.cpp \ PredictorEvalOnly_test.cpp \ PredictorSeedOracle_test.cpp \ + PredictorSeedExtensionKinetic_test.cpp \ Matrix_test.cpp \ NussinovHandler_test.cpp \ RnaSequence_test.cpp \ diff --git a/tests/PredictorSeedExtensionKinetic_test.cpp b/tests/PredictorSeedExtensionKinetic_test.cpp new file mode 100644 index 00000000..9ed293d9 --- /dev/null +++ b/tests/PredictorSeedExtensionKinetic_test.cpp @@ -0,0 +1,528 @@ +#include "catch.hpp" + +#undef NDEBUG + +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/InteractionEnergyVrna.h" +#include "IntaRNA/OutputHandler.h" +#include "IntaRNA/PredictorSeedExtensionKinetic.h" +#include "IntaRNA/SeedHandlerExplicit.h" +#include "IntaRNA/VrnaHandler.h" + +#include +#include +#include +#include +#include +#include +#include + +using namespace IntaRNA; + +namespace { + +using Pair = std::pair; +using Chain = std::vector; +using Bounds = std::array; + +class KineticOutput final : public OutputHandler { +public: + explicit KineticOutput(const OutputConstraint & c); + void add(const Interaction & i) override; + std::vector interactions; +}; + +KineticOutput::KineticOutput(const OutputConstraint & c) : OutputHandler(c) {} + +void KineticOutput::add(const Interaction & i) { + if (!i.basePairs.empty()) interactions.push_back(i); + ++reportedInteractions; +} + +// A nonmonotone table deliberately also models imported accessibility data. +class KineticAccessibility final : public AccessibilityDisabled { +public: + KineticAccessibility(const RnaSequence & s, size_t maxLength = 0); + E_type getED(size_t from, size_t to) const override; + std::map values; +}; + +KineticAccessibility::KineticAccessibility(const RnaSequence & s, size_t maxLength) + : AccessibilityDisabled(s, maxLength, NULL) {} + +E_type KineticAccessibility::getED(size_t from, size_t to) const { + const E_type base = AccessibilityDisabled::getED(from, to); + if (base == ED_UPPER_BOUND) return base; + auto i = values.find({from, to}); + return i == values.end() ? 0 : i->second; +} + +class KineticEnergy final : public InteractionEnergyBasePair { +public: + KineticEnergy(const Accessibility & a, const ReverseAccessibility & b, + size_t m1 = 3, size_t m2 = 3); + E_type getE_interLeft(size_t i, size_t j, size_t k, size_t l) const override; + E_type getE(size_t i, size_t j, size_t k, size_t l, E_type h) const override; + bool customLoops = false; + std::map loops; + std::map boundaryTerms; + mutable std::map loopCalls; +}; + +KineticEnergy::KineticEnergy(const Accessibility & a, const ReverseAccessibility & b, + size_t m1, size_t m2) + : InteractionEnergyBasePair(a, b, m1, m2, false, 1., -100, 3, 0, false) {} + +E_type KineticEnergy::getE_interLeft(size_t i, size_t j, size_t k, size_t l) const { + ++loopCalls[{i,j,k,l}]; + if (!isValidInternalLoop(i, j, k, l)) return E_INF; + if (!customLoops) return InteractionEnergyBasePair::getE_interLeft(i, j, k, l); + auto p = loops.find({i,j,k,l}); + return p == loops.end() ? E_INF : p->second; +} + +E_type KineticEnergy::getE(size_t i, size_t j, size_t k, size_t l, E_type h) const { + const E_type base = InteractionEnergyBasePair::getE(i,j,k,l,h); + auto p = boundaryTerms.find({i,j,k,l}); + return E_isINF(base) ? E_INF : base + (p == boundaryTerms.end() ? 0 : p->second); +} + +struct KineticFixture { + RnaSequence first, second; + KineticAccessibility acc1, acc2; + ReverseAccessibility reversed; + KineticEnergy energy; + KineticFixture(size_t n = 7, size_t m1 = 3, size_t m2 = 3, + size_t span1 = 0, size_t span2 = 0); +}; + +KineticFixture::KineticFixture(size_t n, size_t m1, size_t m2, size_t span1, size_t span2) + : first("target", std::string(n, 'G')), second("query", std::string(n, 'C')), + acc1(first,span1), acc2(second,span2), reversed(acc2), energy(acc1,reversed,m1,m2) {} + +bool stacked(const Pair & a, const Pair & b) { + return a.first+1 == b.first && a.second+1 == b.second; +} + +Bounds bounds(const Chain & chain) { + return {chain.front().first,chain.back().first,chain.front().second,chain.back().second}; +} + +// Independent reference: reconstruct every candidate's complete chain, check +// its topology and recompute all loop and boundary contributions from scratch. +E_type chainEnergy(const InteractionEnergy & energy, const Chain & chain, + const OutputConstraint & out) +{ + const auto b = bounds(chain); + if (b[1]-b[0]+1 > energy.getAccessibility1().getMaxLength() + || b[3]-b[2]+1 > energy.getAccessibility2().getMaxLength()) return E_INF; + E_type h = energy.getE_init(); + for (size_t p = 0; p < chain.size(); ++p) { + if (!energy.areComplementary(chain[p].first,chain[p].second)) return E_INF; + if (p == 0) continue; + const E_type loop = energy.getE_interLeft(chain[p-1].first,chain[p].first, + chain[p-1].second,chain[p].second); + if (E_isINF(loop)) return E_INF; + h += loop; + } + return energy.getE(b[0],b[1],b[2],b[3],h); +} + +std::vector oracle(const InteractionEnergy & energy, Chain chain, + const OutputConstraint & out, char score, + const IndexRange & r1, const IndexRange & r2) +{ + std::vector result; + if (E_isINF(chainEnergy(energy,chain,out))) return result; + while (true) { + const auto b = bounds(chain); + const E_type current = chainEnergy(energy,chain,out); + if ((!out.noGUend || (!energy.isGU(b[0],b[2]) && !energy.isGU(b[1],b[3]))) + && energy.getED1(b[0],b[1]) <= out.maxED + && energy.getED2(b[2],b[3]) <= out.maxED) result.push_back(chain); + bool found = false; + Chain best; + std::tuple bestKey; + // Enumerate absolute endpoints in deliberately different order from + // production's side/loop-size enumeration. + for (size_t x = r1.from; x <= r1.to; ++x) { + for (size_t y = r2.from; y <= r2.to; ++y) { + const bool left = x < b[0] && y < b[2]; + const bool right = x > b[1] && y > b[3]; + if (!left && !right) continue; + const size_t s1 = left ? b[0]-x-1 : x-b[1]-1; + const size_t s2 = left ? b[2]-y-1 : y-b[3]-1; + if (s1 > energy.getMaxInternalLoopSize1() || s2 > energy.getMaxInternalLoopSize2()) continue; + if (s1+s2 != 0 && (out.noGUend || !energy.isInternalLoopGUallowed()) + && (energy.isGU(left ? b[0] : b[1], left ? b[2] : b[3]) || energy.isGU(x,y))) continue; + for (bool macro : {false,true}) { + if (!macro && s1+s2 != 0) continue; + Chain trial = chain; + trial.push_back({x,y}); + if (macro) { + if (left && (x == r1.from || y == r2.from)) continue; + if (right && (x == r1.to || y == r2.to)) continue; + trial.push_back(left ? Pair{x-1,y-1} : Pair{x+1,y+1}); + } + std::sort(trial.begin(),trial.end()); + const E_type total = chainEnergy(energy,trial,out); + if (E_isINF(total) || total >= current) continue; + const size_t denominator = score == 'A' ? 1 : score == 'B' ? 1+s1+s2 : 1+2*std::max(s1,s2); + const auto key = std::make_tuple(static_cast(total-current)/denominator, + left ? 0 : 1,s1+s2,s1,macro); + if (!found || key < bestKey) { found = true; bestKey = key; best = trial; } + } + } + } + if (!found) return result; + chain = best; + } +} + +std::string seedEncoding(const Chain & seed, size_t n2) { + const auto b = bounds(seed); + std::string a(b[1]-b[0]+1,'.'), bReverse(b[3]-b[2]+1,'.'); + for (const auto & p : seed) { a[p.first-b[0]] = '|'; bReverse[p.second-b[2]] = '|'; } + std::reverse(bReverse.begin(),bReverse.end()); + return std::to_string(b[0]+1)+a+"&"+std::to_string(n2-b[3])+bReverse; +} + +SeedConstraint seedConstraint(const std::string & explicitSeed) { + return SeedConstraint(2,20,20,20,E_INF,Accessibility::ED_UPPER_BOUND,E_INF, + IndexRangeList(""),IndexRangeList(""),explicitSeed,false,false,false); +} + +std::vector predict(const InteractionEnergy & energy, const Chain & seed, + const OutputConstraint & out, char score = 'A', + const IndexRange & r1 = IndexRange(0,RnaSequence::lastPos), + const IndexRange & r2 = IndexRange(0,RnaSequence::lastPos)) +{ + const auto sc = seedConstraint(seedEncoding(seed,energy.size2())); + KineticOutput output(out); + PredictorSeedExtensionKinetic predictor(energy,output,NULL,new SeedHandlerExplicit(energy,sc),score); + predictor.predict(r1,r2); + return output.interactions; +} + +Chain internalChain(const InteractionEnergy & energy, const Interaction & interaction) { + Chain result; + for (const auto & p : interaction.basePairs) result.push_back({energy.getIndex1(p),energy.getIndex2(p)}); + return result; +} + +void checkOracle(const InteractionEnergy & energy, const Chain & seed, + const OutputConstraint & out, char score, + IndexRange r1 = IndexRange(0,RnaSequence::lastPos), + IndexRange r2 = IndexRange(0,RnaSequence::lastPos)) +{ + r1.to = std::min(r1.to,energy.size1()-1); r2.to = std::min(r2.to,energy.size2()-1); + auto expected = oracle(energy,seed,out,score,r1,r2); + expected.erase(std::remove_if(expected.begin(),expected.end(),[&](const Chain & c) { + return chainEnergy(energy,c,out) >= out.maxE; + }),expected.end()); + std::sort(expected.begin(),expected.end(),[&](const Chain & a,const Chain & b) { + return chainEnergy(energy,a,out) < chainEnergy(energy,b,out); + }); + const auto actual = predict(energy,seed,out,score,r1,r2); + REQUIRE(actual.size() == std::min(expected.size(),out.reportMax)); + for (size_t p = 0; p < actual.size(); ++p) { + REQUIRE(internalChain(energy,actual[p]) == expected[p]); + REQUIRE(actual[p].energy == chainEnergy(energy,expected[p],out)); + } +} + +} // namespace + +TEST_CASE("Kinetic extension agrees with independent whole-chain oracle", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + const Chain seed{{2,2},{3,3}}; + for (char score : {'A','B','C'}) { + for (bool noLP : {false,true}) { + KineticFixture fixture(8,1,3,6,7); + // Nonmonotone values distinguish exact candidate EDs from a + // farthest-end penalty and change full boundary energy deltas. + fixture.acc1.values = {{{0,3},700},{{1,3},30},{{2,5},500},{{2,6},20}}; + fixture.acc2.values = {{{2,5},450},{{1,5},20}}; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP); + checkOracle(fixture.energy,seed,out,score); + checkOracle(fixture.energy,seed,out,score,IndexRange(1,6),IndexRange(0,7)); + } + } +} + +TEST_CASE("Kinetic scoring and strict downhill acceptance", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + const Chain seed{{2,2},{3,3}}; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + SECTION("A B and C select distinct justified moves") { + KineticFixture f(8); + f.energy.customLoops = true; + f.energy.loops = {{{2,3,2,3},-100},{{3,6,3,4},-420},{{3,5,3,5},-390},{{3,4,3,4},-120},{{6,7,4,5},0},{{5,6,5,6},0}}; + const std::map expected{{'A',{7,5}},{'B',{7,5}},{'C',{6,6}}}; + for (char score : {'A','B','C'}) { + checkOracle(f.energy,seed,out,score); + const auto actual = predict(f.energy,seed,out,score); + REQUIRE(internalChain(f.energy,actual.front()).back() == expected.at(score)); + } + // B prefers the immediate stack when its per-distance gain wins. + f.energy.loops[{3,4,3,4}] = -150; + REQUIRE(internalChain(f.energy,predict(f.energy,seed,out,'B').front()).back() == Pair(4,4)); + } + SECTION("zero and positive moves stop even when the local loop is favorable") { + for (E_type delta : {0,1,200}) { + KineticFixture f(5); + f.energy.customLoops = true; + f.energy.loops = {{{2,3,2,3},-100},{{3,4,3,4},-100}}; + f.energy.boundaryTerms[{2,4,2,4}] = 100+delta; + const auto actual = predict(f.energy,seed,out); + REQUIRE(actual.size() == 1); + REQUIRE(internalChain(f.energy,actual.front()) == seed); + } + } + SECTION("ties prefer left then smaller total gap then smaller first gap") { + KineticFixture f(7); + f.energy.customLoops = true; + f.energy.loops = {{{2,3,2,3},-100},{{0,2,1,2},-100},{{1,2,0,2},-100},{{1,2,1,2},-100},{{3,4,3,4},-100}}; + // Limit each span to three: whichever move wins blocks the other side. + KineticFixture shortF(7,3,3,3,3); + shortF.energy.customLoops = true; shortF.energy.loops = f.energy.loops; + REQUIRE(internalChain(shortF.energy,predict(shortF.energy,seed,out).front()).front() == Pair(1,1)); + f.energy.loops = {{{3,4,3,4},-100},{{2,3,1,3},-50},{{1,2,0,1},-50}, + {{1,3,2,3},-50},{{0,1,1,2},-50}}; + const Chain middle{{3,3},{4,4}}; + // Equal total gap: s1=0 (outer pair 1,0) wins over s1=1 (0,1). + const auto paths = oracle(f.energy,middle,out,'A',IndexRange(0,6),IndexRange(0,6)); + REQUIRE(paths.at(1).front() == Pair(1,0)); + checkOracle(f.energy,middle,out,'A'); + } +} + +TEST_CASE("Kinetic always uses stacked extensions and trusts explicit seeds", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(6,1,1); + f.energy.customLoops = true; + f.energy.loops = {{{0,1,0,1},-100},{{1,3,1,3},50},{{3,4,3,4},-330}}; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,true); + const Chain seed{{0,0},{1,1}}; + checkOracle(f.energy,seed,out,'A'); + const auto actual = predict(f.energy,seed,out); + REQUIRE(actual.size() == 2); + REQUIRE(actual.front().basePairs.size() == 4); + REQUIRE(actual.front().energy == -480); + SECTION("B and C count one macro move in their denominator") { + KineticFixture ranked(7,1,1,5,5); + ranked.energy.customLoops = true; + ranked.energy.loops = {{{2,3,2,3},-100},{{1,2,1,2},-80}, + {{3,5,3,5},50},{{5,6,5,6},-330}}; + const Chain middle{{2,2},{3,3}}; + for (char score : {'B','C'}) { + checkOracle(ranked.energy,middle,out,score); + const auto rankedResult = predict(ranked.energy,middle,out,score); + // -280/3 beats -80. Incorrectly counting the two added pairs + // would instead produce -280/4 and choose the left stack. + REQUIRE(rankedResult.front().basePairs.size() == 4); + REQUIRE(internalChain(ranked.energy,rankedResult.front()).back() == Pair(6,6)); + } + } + SECTION("lonely explicit seeds are accepted") { + REQUIRE_FALSE(predict(f.energy,Chain{{1,1},{3,3}},out).empty()); + } + SECTION("both additional pairs must fit both strand spans") { + KineticFixture shortF(6,1,1,4,6); + shortF.energy.customLoops = true; shortF.energy.loops = f.energy.loops; + const auto limited = predict(shortF.energy,seed,out); + REQUIRE(limited.size() == 1); + REQUIRE(limited.front().basePairs.size() == 2); + } +} + +TEST_CASE("Kinetic noGU output retains a valid prefix", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + RnaSequence first("target","GGG"), second("query","UCC"); + AccessibilityDisabled a(first,0,NULL), b(second,0,NULL); + ReverseAccessibility reversed(b); + InteractionEnergyBasePair energy(a,reversed,1,1); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,false,true); + const Chain seed{{0,0},{1,1}}; + checkOracle(energy,seed,out,'A'); + const auto actual = predict(energy,seed,out); + REQUIRE(actual.size() == 1); + REQUIRE(internalChain(energy,actual.front()) == seed); + SECTION("a transient GU endpoint can become internal after another stack") { + RnaSequence t("target","GGGG"), q("query","CUCC"); + AccessibilityDisabled at(t,0,NULL), aq(q,0,NULL); + ReverseAccessibility rq(aq); + InteractionEnergyBasePair model(at,rq,1,1); + checkOracle(model,seed,out,'A'); + const auto recovered = predict(model,seed,out); + REQUIRE(recovered.size() == 2); + REQUIRE(recovered.front().basePairs.size() == 4); + REQUIRE(recovered.front().energy == -400); + REQUIRE(model.isGU(2,2)); + } +} + +TEST_CASE("Kinetic full ViennaRNA energy agrees with whole-chain oracle", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + RnaSequence first("target","AGCGACGCA"), second("query","UGCGUCGCU"); + KineticAccessibility a(first), b(second); + a.values = {{{0,4},140},{{1,4},20},{{1,5},350},{{1,6},30},{{2,7},130}}; + b.values = {{{0,4},20},{{1,4},60},{{1,5},20},{{2,6},90}}; + ReverseAccessibility reversed(b); + VrnaHandler vrna(37,"Turner04",false,false); + const Chain seed{{2,2},{3,3}}; + for (bool dangles : {false,true}) { + InteractionEnergyVrna energy(a,reversed,vrna,3,2,false,37,dangles); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + for (char score : {'A','B','C'}) checkOracle(energy,seed,out,score); + } +} + +TEST_CASE("Kinetic refuses undefined ensemble output and invalid scores", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f; + const auto sc = seedConstraint("3||&4||"); + OutputConstraint ensemble(1,OutputConstraint::OVERLAP_BOTH,0,E_INF,false,false,false,true); + KineticOutput ensembleOut(ensemble); + REQUIRE_THROWS_AS(PredictorSeedExtensionKinetic(f.energy,ensembleOut,NULL, + new SeedHandlerExplicit(f.energy,sc)),std::invalid_argument); + OutputConstraint ordinary; + KineticOutput ordinaryOut(ordinary); + REQUIRE_THROWS_AS(PredictorSeedExtensionKinetic(f.energy,ordinaryOut,NULL, + new SeedHandlerExplicit(f.energy,sc),'Z'),std::invalid_argument); +} + +TEST_CASE("Kinetic Turner loop is rescued by its atomic stack", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + // A single query-strand A bulge prevents a direct stack after the seed. + // This is the Turner2004 witness reproduced by the first council: + // +0.50 kcal/mol loop plus -3.30 kcal/mol outward stack. + RnaSequence first("target","CCCC"), second("query","GGAGG"); + AccessibilityDisabled a(first,0,NULL), b(second,0,NULL); + ReverseAccessibility reversed(b); + VrnaHandler vrna(37,"Turner04",false,false); + InteractionEnergyVrna energy(a,reversed,vrna,1,1,false,0,false); + const E_type loop = energy.getE_interLeft(1,2,1,3); + const E_type stack = energy.getE_interLeft(2,3,3,4); + REQUIRE(loop == 50); + REQUIRE(stack == -330); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,true); + const Chain seed{{0,0},{1,1}}; + checkOracle(energy,seed,out,'A'); + const auto actual = predict(energy,seed,out); + REQUIRE(actual.size() == 2); + REQUIRE((internalChain(energy,actual.front()) == Chain{{0,0},{1,1},{2,3},{3,4}})); + REQUIRE(actual.front().energy - actual.back().energy == loop+stack); +} + +TEST_CASE("Kinetic nonoverlap reporting recovers shorter prefixes", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(8,0,0,5,5); + f.energy.customLoops = true; + for (size_t p = 0; p+1 < 8; ++p) f.energy.loops[{p,p+1,p,p+1}] = p < 3 ? -100 : -200; + const auto sc = seedConstraint("1||&7||,4||&4||"); + for (bool needBPs : {false,true}) { + OutputConstraint out(3,OutputConstraint::OVERLAP_NONE,E_INF,E_INF,false,false,false,false,needBPs); + KineticOutput output(out); + PredictorSeedExtensionKinetic predictor(f.energy,output,NULL,new SeedHandlerExplicit(f.energy,sc)); + predictor.predict(); + REQUIRE(output.interactions.size() == 2); + REQUIRE(output.interactions[0].energy == -900); + REQUIRE(output.interactions[1].energy == -200); + REQUIRE((bounds(internalChain(f.energy,output.interactions[0])) == Bounds{{3,7,3,7}})); + REQUIRE((bounds(internalChain(f.energy,output.interactions[1])) == Bounds{{0,1,0,1}})); + if (needBPs) { + REQUIRE(output.interactions[0].basePairs.size() == 5); + REQUIRE(output.interactions[1].basePairs.size() == 2); + } + // Reuse the instance with a different range and ensure old candidates + // and seed metadata cannot leak into the second prediction. + output.interactions.clear(); + predictor.predict(IndexRange(0,2),IndexRange(0,2)); + REQUIRE(output.interactions.size() == 1); + REQUIRE(output.interactions.front().energy == -300); + REQUIRE((bounds(internalChain(f.energy,output.interactions.front())) == Bounds{{0,2,0,2}})); + } +} + +TEST_CASE("Kinetic annotations retain handler-provided lonely explicit seeds", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(5,0,0); + const auto sc = seedConstraint("2||&3||,1|&5|"); + OutputConstraint out(1,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,true); + KineticOutput output(out); + PredictorSeedExtensionKinetic predictor(f.energy,output,NULL,new SeedHandlerExplicit(f.energy,sc)); + predictor.predict(); + REQUIRE(output.interactions.size() == 1); + const auto & result = output.interactions.front(); + REQUIRE(result.basePairs.size() == 5); + REQUIRE(result.seed != NULL); + REQUIRE(result.seed->size() == 2); +} + +TEST_CASE("Kinetic two-stack moves cross a local barrier atomically", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(4,0,0); + f.energy.customLoops = true; + f.energy.loops = {{{0,1,0,1},-100},{{1,2,1,2},50},{{2,3,2,3},-200}}; + const Chain seed{{0,0},{1,1}}; + for (bool noLP : {false,true}) { + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP); + checkOracle(f.energy,seed,out,'A'); + const auto actual = predict(f.energy,seed,out); + REQUIRE(actual.size() == 2); + REQUIRE(actual.front().basePairs.size() == 4); + REQUIRE(actual.front().energy == -350); + } +} + +TEST_CASE("Kinetic exact move ties prefer the single stack", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(4,0,0); + f.energy.customLoops = true; + f.energy.loops = {{{0,1,0,1},-100},{{1,2,1,2},-100},{{2,3,2,3},0}}; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + for (char score : {'A','B','C'}) { + const auto actual = predict(f.energy,Chain{{0,0},{1,1}},out,score); + REQUIRE(actual.size() == 2); + REQUIRE(actual.front().basePairs.size() == 3); + REQUIRE(actual.front().energy == -300); + } +} + +TEST_CASE("Kinetic caches the unchanged end and trusts the seed energy", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + KineticFixture f(8,0,0); + f.energy.customLoops = true; + for (size_t i = 0; i+1 < 8; ++i) f.energy.loops[{i,i+1,i,i+1}] = i < 2 ? -200 : -100; + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF); + const auto result = predict(f.energy,Chain{{2,2},{3,3}},out); + REQUIRE(result.front().basePairs.size() == 8); + // Seed energy is evaluated by SeedHandlerExplicit just once, never by + // the predictor. The right first stack is shared by single/two-pair moves + // and is not reevaluated after the left end grows. + REQUIRE(f.energy.loopCalls.at(Bounds{{2,3,2,3}}) == 1); + REQUIRE(f.energy.loopCalls.at(Bounds{{3,4,3,4}}) == 2); +} + +TEST_CASE("Kinetic oracle covers varied sequences and endpoint energies", "[PredictorSeedExtensionKinetic]") { + #include "testEasyLoggingSetup.icc" + unsigned int random = 254; + const auto next = [&]() { random = random*1664525u+1013904223u; return random; }; + for (size_t sample = 0; sample < 32; ++sample) { + std::string t(9,'G'), q(9,'C'); + for (size_t i = 0; i < 9; ++i) { t[i] = "ACGU"[(next() >> 16)%4]; q[i] = "ACGU"[(next() >> 16)%4]; } + t[3] = t[4] = 'G'; q[4] = q[5] = 'C'; + RnaSequence first("t",t), second("q",q); + KineticAccessibility a(first), b(second); + for (size_t i = 0; i < 9; ++i) for (size_t j = i; j < 9; ++j) { + a.values[{i,j}] = (next() >> 16)%300; + b.values[{i,j}] = (next() >> 16)%300; + } + ReverseAccessibility reversed(b); + InteractionEnergyBasePair energy(a,reversed,2,3); + OutputConstraint out(100,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,sample%2,sample%3 == 0); + for (char score : {'A','B','C'}) checkOracle(energy,Chain{{3,3},{4,4}},out,score); + } +} diff --git a/tests/runKineticSeedExtension.sh b/tests/runKineticSeedExtension.sh new file mode 100755 index 00000000..e6054c2d --- /dev/null +++ b/tests/runKineticSeedExtension.sh @@ -0,0 +1,161 @@ +#!/usr/bin/env bash +# Exercise kinetic dispatch, greedy paths, output constraints and energy accounting. +set -euo pipefail +bin="$INTARNABINPATH/src/bin/IntaRNA" +tmp=$(mktemp -d) +trap 'rm -rf "$tmp"' EXIT +common=(--target=GGGGGG --query=CCCCCC --energy=B --acc=N + --threads=1 --default-log-file=/dev/null) +kinetic=(--model=X --mode=K '--seedTQ=3||&3||') +csv=(--outMode=C --outCsvCols=hybridDB,E) + +# Six base pairs have independently known energy -6 in the base-pair model. +printf 'hybridDB;E\n1||||||&1||||||;-6\n' > "$tmp/expected" +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" > "$tmp/default" +cmp "$tmp/expected" "$tmp/default" +for score in A B C; do + for noLP in false true; do + "$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" \ + --kineticScore="$score" --outNoLP="$noLP" --outNoGUend > "$tmp/score" + cmp "$tmp/expected" "$tmp/score" + done +done + +# Ties choose the left extension, including the antiparallel query coordinates. +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" --intLenMax=3 > "$tmp/left" +printf 'hybridDB;E\n2|||&3|||;-3\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/left" + +# Boundary-only output and CLI range/index conversions use the same trajectory. +"$bin" "${common[@]}" "${kinetic[@]}" --outMode=C \ + --outCsvCols=start1,end1,start2,end2,E > "$tmp/boundaries" +printf 'start1;end1;start2;end2;E\n1;6;1;6;-6\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/boundaries" +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" \ + --tRegion=2-5 --qRegion=2-5 > "$tmp/ranges" +printf 'hybridDB;E\n2||||&2||||;-4\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/ranges" +# Signed IntaRNA coordinates skip zero: target position three is labeled 1. +"$bin" "${common[@]}" --model=X --mode=K '--seedTQ=1||&7||' "${csv[@]}" \ + --tIdxPos0=-2 --qIdxPos0=5 > "$tmp/shifted" +printf 'hybridDB;E\n-2||||||&5||||||;-6\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/shifted" + +# Every committed prefix is available for suboptimal output, in energy order. +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" --outNumber=10 > "$tmp/prefixes" +printf 'hybridDB;E\n1||||||&1||||||;-6\n1||||&3||||;-4\n3||&3||;-2\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/prefixes" + +# A zero report count still permits energy tracking and never needs traceback. +"$bin" "${common[@]}" "${kinetic[@]}" "${csv[@]}" --outNumber=0 \ + --out="tMinE:$tmp/min-energy" > "$tmp/zero" +printf 'hybridDB;E\n' > "$tmp/empty" +cmp "$tmp/empty" "$tmp/zero" +test -s "$tmp/min-energy" +grep -q -- '-6' "$tmp/min-energy" + +# Handler-provided explicit seeds can contain lonely pairs. +"$bin" "${common[@]}" --model=X --mode=K '--seedTQ=3|&3|' \ + "${csv[@]}" --outNoLP > "$tmp/lonely" +printf 'hybridDB;E\n1|||||&1|||||;-5\n' > "$tmp/expected" +cmp "$tmp/expected" "$tmp/lonely" +"$bin" --target=GG --query=UU --energy=B --acc=N --model=X --mode=K \ + '--seedTQ=1||&1||' "${csv[@]}" --outNoGUend --default-log-file=/dev/null > "$tmp/gu" +cmp "$tmp/empty" "$tmp/gu" + +# Config-file selection must behave just like command-line selection. +printf 'model=X\nmode=K\nkineticScore=C\nseedTQ=3||&3||\n' > "$tmp/parameters" +"$bin" "${common[@]}" "${csv[@]}" --parameterFile="$tmp/parameters" > "$tmp/configured" +cmp "$tmp/default" "$tmp/configured" + +# Re-evaluate all reported structures independently through IntaRNAeval. +# This checks the complete energy, including accessibility, terminal and dangle terms. +thermo=(--target=AGCGACGCA --query=UGCGUCGCU --accW=0 --accL=0 --temperature=25 + --energyAdd=1.2 --outMode=C --outMaxE=100 --outNumber=5 + --outCsvCols=hybridDB,E,ED1,ED2,E_init,E_loops,E_dangleL,E_dangleR,E_endL,E_endR,E_add + --threads=1 --default-log-file=/dev/null) +for score in A B C; do + for dangles in false true; do + "$bin" "${thermo[@]}" --model=X --mode=K '--seedTQ=3|||&5|||' \ + --kineticScore="$score" --energyNoDangles="$dangles" --outNoLP --outNoGUend \ + > "$tmp/predicted" + test "$(wc -l < "$tmp/predicted")" -gt 1 + structures=$(awk -F';' 'NR>1 {printf "%s%s", sep, $1; sep=":"}' "$tmp/predicted") + "$bin" "${thermo[@]}" --rri="$structures" --energyNoDangles="$dangles" > "$tmp/evaluated" + cmp "$tmp/predicted" "$tmp/evaluated" + done +done + +expect_error() { + local expected="$1" + shift + local status=0 + "$bin" "$@" > "$tmp/bad.out" 2> "$tmp/bad.err" || status=$? + test "$status" -eq 1 || test "$status" -eq 255 + grep -q -- "$expected" "$tmp/bad.out" "$tmp/bad.err" +} +for model in S P B; do + expect_error 'only with --model=X' "${common[@]}" --model="$model" --mode=K +done +expect_error 'incompatible with --noSeed' "${common[@]}" "${kinetic[@]}" --noSeed +# Explicitly supplying the default A must be rejected outside K as well. +for score in A B C; do + expect_error 'kineticScore requires' "${common[@]}" --mode=H --kineticScore="$score" +done +expect_error 'kineticScore' "${common[@]}" "${kinetic[@]}" --kineticScore=D +expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" --outMode=E +for column in Zall Eall EallTotal P_E; do + expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" \ + --outMode=C --outCsvCols="E,$column" +done +for output in spotProb qSpotProb tSpotProb; do + expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" --out="$output:$tmp/probability" +done +expect_error 'equilibrium ensemble' "${common[@]}" "${kinetic[@]}" --out="spotProb:1&1:$tmp/probability" + +# Evaluation ignores prediction controls, including invalid kinetic scores. +"$bin" "${common[@]}" "${csv[@]}" '--rri=1||||||&1||||||' \ + --model=S --mode=K --kineticScore=D > "$tmp/evaluation" +cmp "$tmp/default" "$tmp/evaluation" +# A missing or false flag is promoted once, with a visible INFO message. +logging=("${common[@]}") +logging[${#logging[@]}-1]="--default-log-file=$tmp/info.log" +for setting in absent false true; do + : > "$tmp/info.log" + flags=() + if [ "$setting" != absent ]; then flags=(--outNoLP="$setting"); fi + "$bin" "${logging[@]}" "${kinetic[@]}" "${csv[@]}" "${flags[@]}" > "$tmp/info.out" 2> "$tmp/info.err" + if [ "$setting" = true ]; then + ! grep -q 'setting --outNoLP=true' "$tmp/info.log" + else + grep -q 'INFO.*setting --outNoLP=true' "$tmp/info.log" + fi +done +# Both personality entry points select mode K and noLP by default. +ln -s "$bin" "$tmp/IntaRNAsnap" +for invocation in binary option; do + args=() + executable="$tmp/IntaRNAsnap" + if [ "$invocation" = option ]; then + executable="$bin" + args=(--personality=IntaRNAsnap) + fi + : > "$tmp/info.log" + "$executable" "${args[@]}" "${logging[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ + --outNumber=10 > "$tmp/kix-prefixes" + cmp "$tmp/prefixes" "$tmp/kix-prefixes" + ! grep -q 'setting --outNoLP=true' "$tmp/info.log" + "$executable" "${args[@]}" "${common[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ + --mode=S > "$tmp/seed-only" + printf 'hybridDB;E\n3||&3||;-2\n' > "$tmp/expected" + cmp "$tmp/expected" "$tmp/seed-only" +done +# Explicitly disabling noLP cannot disable the mode K extension invariant. +"$bin" --personality=IntaRNAsnap "${logging[@]}" '--seedTQ=3||&3||' "${csv[@]}" \ + --outNoLP=false > "$tmp/kix-noLP" +cmp "$tmp/default" "$tmp/kix-noLP" +grep -q 'setting --outNoLP=true' "$tmp/info.log" +# Evaluation remains available under the personality and ignores its defaults. +"$tmp/IntaRNAsnap" "${common[@]}" "${csv[@]}" '--rri=1||||||&1||||||' > "$tmp/kix-eval" +cmp "$tmp/default" "$tmp/kix-eval" +echo 'Kinetic seed-extension and IntaRNAsnap CLI checks passed'