From 20af34770ceb0c76bcaadd31184ee34f02bc43fe Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Fri, 9 Oct 2026 11:59:50 +0200 Subject: [PATCH 1/9] feat: add SVG base-pair probability plots (step 8) --- ChangeLog | 19 ++ README.md | 60 ++++- src/IntaRNA/BasePairProbabilities.cpp | 10 +- src/IntaRNA/BasePairProbabilities.h | 15 +- src/IntaRNA/BasePairProbabilitySvg.cpp | 210 ++++++++++++++++++ src/IntaRNA/BasePairProbabilityWriter.cpp | 8 +- src/IntaRNA/BasePairProbabilityWriter.h | 20 ++ src/IntaRNA/Makefile.am | 1 + .../PredictorMfeEns2dSeedExtension.cpp | 9 +- src/bin/CommandLineParsing.cpp | 21 +- src/bin/CommandLineParsing.h | 4 +- tests/BasePairProbabilities_test.cpp | 41 ++++ tests/Makefile.am | 2 +- tests/runBasePairProbabilities.sh | 1 + tests/testBasePairProbabilitySvg.py | 157 +++++++++++++ 15 files changed, 559 insertions(+), 19 deletions(-) create mode 100644 src/IntaRNA/BasePairProbabilitySvg.cpp create mode 100644 tests/testBasePairProbabilitySvg.py diff --git a/ChangeLog b/ChangeLog index bd016930..d04228f5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,10 @@ ## Interface and handling +- --out=bpsvg:FILE renders actual-pair probabilities with nucleotide/accessibility + frames, a shared CSS palette, signed-coordinate guides, hover values and admitted + seed outlines (PR #258, step 8) + - --out=bpProb:FILE writes actual-pair conditional probabilities for exact model P/mode M stack-seed ensembles, with disjoint-region merging, original coordinates and success-only publication (issue #257, implementation step 6) @@ -156,6 +160,21 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261009 Alexander Mitrofanov + * BasePairProbabilityWriter, BasePairProbabilitySvg, CommandLineParsing : + + add success-only SVG dot plots on files, gzip and synchronized standard + streams; share CSV normalization and use original sequence coordinates + + show sequence labels, 1-nt accessibility, seven CSS probability classes, + legends, signed 10/50-nt guides, missing-zero boundary and hover values + * BasePairProbabilities, PredictorMfeEns2dSeedExtension : + + optionally retain admitted seed-cell masks with completed regional masses; + annotate seed cells without changing probability or ranked-output semantics + * tests, README.md, src/IntaRNA/Makefile.am, tests/Makefile.am : + + document the SVG contract and register the writer and regression coverage + for coordinates, seed admission, accessibility, empty/failed results and I/O + + investigate the referenced dev-frank profile tracker files: their current + versions contain profile output but no SVG implementation (PR #258 step 8) + 261008 Alexander Mitrofanov * tests, README.md, doc/analysis, doc/seeded-probabilities-example.cpp : + finish probability regression and native performance gates, document the diff --git a/README.md b/README.md index 32f77132..fcb80537 100644 --- a/README.md +++ b/README.md @@ -1876,6 +1876,7 @@ targeted file/stream name: - `qSpotProb:`/`tSpotProb:` [query/target's spot probability profile](#profileSpotProb) (CSV format), respectively - `spotProb:` [all spot probabilities](#spotProb) (CSV format) +- `bpsvg:` [SVG base-pair probability dot plot](#bpsvg), with accessibility frames and seed outlines - `bpProb:` [actual base-pair probabilities](#bpProb) for exact stack-seeded ensembles (CSV matrix) - `qMinE:`/`tMinE:` [the query/target's minimal interaction energy profile](profileMinE) (CSV format), respectively - `pMinE:` [minimal interaction energy for all query-target index pairs](pairMinE) (CSV format) @@ -1883,7 +1884,7 @@ targeted file/stream name: - `qPu:`/`tPu:` the [query/target's unpaired probabilities](#accessibility) (RNAplfold format; rounded!!), respectively Note, for *multiple sequences* in FASTA input, the provided file names are suffixed by -- `-t#q#` : for `bpProb`, `*spotProb` and `*minE` output, and +- `-t#q#` : for `bpProb`, `bpsvg`, `*spotProb` and `*minE` output, and - `-s#` : for `*Acc` and `*Pu` output, where `#` denotes the according target/query sequence number within the input where numbering starts with 1. @@ -2064,6 +2065,63 @@ constraint must enable `needZall`. The [standalone API example](doc/seeded-proba and [implementation validation report](doc/analysis/seeded-probabilities-validation.md) show the public API and measured limits. + + +### SVG base-pair probability dot plots + +`--out=bpsvg:FILE` renders the same actual-pair probabilities as `bpProb` in a +standalone SVG. Both outputs can be requested together. The same exact +`--model=P --mode=M` stack-seed restrictions, regional ensemble, normalization +and success-only publication rules apply. + +```sh +IntaRNA -t GGAGGG -q CCCC --energy=B --acc=N --seedBP=2 \ + --model=P --mode=M --out=bpsvg:pairs.svg --out=bpProb:pairs.csv +``` + +The query is the horizontal axis and the target the vertical axis, both in +original 5' to 3' order. Sequence names label the axes; nucleotides appear on all +four sides. Their colored backgrounds show **single-nucleotide unpaired +probabilities**, computed from the prediction accessibility model's opening +energies. These are distinct from the interaction-conditional probabilities in +the matrix. Disabled accessibility gives probability one at unconstrained bases. +Each probability rectangle has a hover title with its index (pair), type and +full-precision value. An empty interaction ensemble uses gray `NA` cells while +its accessibility frame remains defined. + +A shared seven-step, white-to-dark-blue palette encodes both probability types: +`[0,0.01)`, `[0.01,0.1)`, `[0.1,0.25)`, `[0.25,0.5)`, `[0.5,0.75)`, +`[0.75,0.9)`, and `[0.9,1]`. The legend appears below the plot. Edit the `.p0` +through `.p6` fill rules in the SVG's ` +)SVG"; + out<<"\n"; + auto text=[&](double tx,double ty,const std::string & value,const std::string & cls,const char * anchor="middle") { + out<<""<\n"; + }; + text(width/2,30,"Base-pair probabilities","heading"); + text(x+w/2,y-65,q.getId()+" (5' to 3')","axis"); + out<<""<\n"; + auto rect=[&](double rx,double ry,Z_type value,bool missing,const std::string & attributes,const std::string & title) { + out<<""<<xml(title)<<": "; + if (missing) out<<"NA (empty interaction ensemble)"; else out<<value; + out<<"\n"; + }; + out<<"\n"; + for(size_t i=0;i\n\n"; + for(bool target:{false,true}) { + const auto & seq=target?t:q; + const auto & acc=target?at:aq; + for(size_t i=0;i=.75?"nt nt-light":"nt"); + } + } + } + out<<"\n\n"; + for(bool target:{false,true}) { + const auto & seq=target?t:q; + for(size_t i=0;i0 && seq.getInOutIndex(i-1)==-1; + if (idx%10!=0 && idx!=1) continue; + // Negative multiples: before the labelled nt (-11|-10); positive: + // after it (10|11). At the missing zero: -1|+1. + const double offset=cell*(i+(idx>0 && !origin?1:0)); + if (idx%10==0 || origin) out<<"\n"; + const std::string label=idx==1?"+1":std::to_string(idx); + if(target) { + text(x-26,y+cell*(i+.5)+4,label,"index","end"); + text(x+w+26,y+cell*(i+.5)+4,label,"index","start"); + } else { + text(x+cell*(i+.5),y-28,label,"index"); + text(x+cell*(i+.5),y+h+38,label,"index"); + } + } + } + out<<"\n\n"; + bool annotated=false; + if(result.collectsSeedPairs()) { + out<<"\n"; + for(size_t i=0;i\n"; + } + out<<"\n"; + } + const double legendX=(width-7*86)/2, legendY=y+h+85; + text(width/2,legendY-14,"Probability ranges (base pairs and single-nucleotide accessibility)","legend"); + out<<"\n"; + for(size_t k=0;k\n"; + std::ostringstream range; + range.imbue(std::locale::classic()); + range<<'['<\n"; + if (annotated) { + out<<"\n"; + text(legendX+28,legendY+61,"Orange outline: pair in an admitted seed within a searched region","legend","start"); + } + if(empty) text(width/2,legendY+90,"Gray / NA: empty interaction ensemble; accessibility remains defined","legend"); + out<<"\n"; + return out.str(); +} + +void BasePairProbabilityWriter::writeSvg(std::ostream & out,const BasePairProbabilities & result,const InteractionEnergy & energy) +{ + emit(out,svgBlock(result,energy)); +} + +void BasePairProbabilityWriter::writeSvgFile(const std::string & filename,const BasePairProbabilities & result,const InteractionEnergy & energy) +{ + emitFile(filename,svgBlock(result,energy)); +} + +} // namespace IntaRNA diff --git a/src/IntaRNA/BasePairProbabilityWriter.cpp b/src/IntaRNA/BasePairProbabilityWriter.cpp index d0aea7af..6d679fdd 100644 --- a/src/IntaRNA/BasePairProbabilityWriter.cpp +++ b/src/IntaRNA/BasePairProbabilityWriter.cpp @@ -36,7 +36,7 @@ void BasePairProbabilityWriter::emit(std::ostream & out,const std::string & data { try { out< file(out==&std::cout || out==&std::cerr?nullptr:out); diff --git a/src/IntaRNA/BasePairProbabilityWriter.h b/src/IntaRNA/BasePairProbabilityWriter.h index 48578daf..41d7c1a2 100644 --- a/src/IntaRNA/BasePairProbabilityWriter.h +++ b/src/IntaRNA/BasePairProbabilityWriter.h @@ -6,6 +6,7 @@ #include namespace IntaRNA { +class InteractionEnergy; /** Explicit writer for a completed actual-pair matrix, independent of trackers. * Whole blocks are synchronized for shared streams. No destructor emits data. */ @@ -30,7 +31,26 @@ class BasePairProbabilityWriter { */ static void writeFile(const std::string & filename,const BasePairProbabilities & result, const RnaSequence & target,const RnaSequence & query,const std::string & separator=";"); + /** Write a complete SVG dot plot with original target rows/query columns. + * Includes 1-nt accessibility, signed coordinate guides and a CSS color legend. + * Seed outlines are included if the result collected admitted seed pairs. + * @param out destination stream; an I/O failure throws + * @param result explicitly finalized raw result + * @param energy original, unshifted energy/accessibility model used in prediction + */ + static void writeSvg(std::ostream & out,const BasePairProbabilities & result, + const InteractionEnergy & energy); + /** Write a complete SVG to a file or STDOUT/STDERR, including gzip support. + * Validates the whole document before opening the destination. + * @param filename destination (already expanded for multi-FASTA) + * @param result explicitly finalized raw result + * @param energy original, unshifted energy/accessibility model + */ + static void writeSvgFile(const std::string & filename,const BasePairProbabilities & result, + const InteractionEnergy & energy); private: + static std::string svgBlock(const BasePairProbabilities & result,const InteractionEnergy & energy); + static void emitFile(const std::string & name,const std::string & data); static std::string block(const BasePairProbabilities & result,const RnaSequence & target, const RnaSequence & query,const std::string & separator); static void emit(std::ostream & out,const std::string & data); diff --git a/src/IntaRNA/Makefile.am b/src/IntaRNA/Makefile.am index 4fc5db3e..178691d6 100644 --- a/src/IntaRNA/Makefile.am +++ b/src/IntaRNA/Makefile.am @@ -106,6 +106,7 @@ libIntaRNA_a_SOURCES = \ StackSeedDomain.cpp \ BasePairProbabilities.cpp \ BasePairProbabilityWriter.cpp \ + BasePairProbabilitySvg.cpp \ SeededPartitionFunction.cpp \ $(libmylib_a_HEADERS) \ Accessibility.cpp \ diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp index 4f50a61e..fb7c8929 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp @@ -185,11 +185,18 @@ PredictorMfeEns2dSeedExtension::predictStackSeeds(const IndexRange & r1,const In const auto first=energy.getBasePair(0,0), last=energy.getBasePair(n-1,m-1); const IndexRange target(first.first,last.first), query(last.second,first.second); Matrix original(n,m,0); + const bool annotate=pairProbabilities->collectsSeedPairs(); + Matrix seedPairs(annotate?n:0,annotate?m:0,0); for(size_t i=0;iaddRegion(target,query,result.z,original); + pairProbabilities->addRegion(target,query,result.z,original,annotate?&seedPairs:nullptr); } } diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index dc864b8c..77455d1e 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -957,6 +957,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) "\n 'tAcc:' (target) ED accessibility values ('tPu'-like format OR '.agz' file ending for IntaRNA's binary format)." "\n 'tPu:' (target) unpaired probability values (RNAplfold format OR '.agz' file ending for IntaRNA's binary format)." "\n 'pMinE:' (target+query) for each index pair the minimal energy of any interaction covering the pair (CSV format)" + "\n 'bpsvg:' (target+query) SVG actual-pair probability plot with accessibility frame and seed annotations; same requirements as bpProb" "\n 'bpProb:' (target+query) actual-pair conditional probabilities for exact --model=P --mode=M stack-seed ensembles (CSV matrix, original strand order; NA for an empty ensemble)" "\n 'spotProb:' (target+query) tracks for a given set of interaction spots their probability to be covered by an interaction. If no spots are provided, probabilities for all index combinations are computed. Spots are encoded by comma-separated 'idxT&idxQ' pairs (target-query). For each spot a probability is provided in concert with the probability that none of the spots (encoded by '0&0') is covered (CSV format). The spot encoding is followed colon-separated by the output stream/file name, eg. '--out=\"spotProb:3&76,59&2:STDERR\"'. NOTE: value has to be quoted due to '&' symbol!" "\nFor each, provide a file name or STDOUT/STDERR to write to the respective output stream." @@ -1246,7 +1247,7 @@ parse(int argc, char** argv) if (hasBasePairProbabilityOutput()) { if (model.val!='P' || mode.val!='M' || noSeedRequired || !rri.empty() || windowWidth.val!=0) - throw error("bpProb requires --model=P --mode=M with stack-only seeds and no windows"); + throw error("bpProb/bpsvg requires --model=P --mode=M with stack-only seeds and no windows"); outNeedsZall=true; } // K needs a seed even before the usual --noSeed model normalization. @@ -1582,7 +1583,7 @@ parse(int argc, char** argv) for (auto c2=c1; noDuplicate && (++c2)!=outPrefix2streamName.end();) { if ( ! c2->second.empty() && boost::iequals( c1->second, c2->second ) ) { // Pair matrices emit synchronized complete blocks on shared standard streams. - if ((c1->first==OP_bpProb || c2->first==OP_bpProb) + if ((c1->first==OP_bpProb || c2->first==OP_bpProb || c1->first==OP_bpsvg || c2->first==OP_bpsvg) && (boost::iequals(c1->second,"STDOUT") || boost::iequals(c1->second,"STDERR"))) continue; noDuplicate = false; throw error("--out argument shows multiple times '"+toString(c1->second)+"' as target file/stream."); @@ -2490,19 +2491,23 @@ getTemperature() const { //////////////////////////////////////////////////////////////////////////// bool CommandLineParsing::hasBasePairProbabilityOutput() const { - return !outPrefix2streamName.at(OutPrefixCode::OP_bpProb).empty(); + return !outPrefix2streamName.at(OutPrefixCode::OP_bpProb).empty() + || !outPrefix2streamName.at(OutPrefixCode::OP_bpsvg).empty(); } BasePairProbabilities * CommandLineParsing::getBasePairProbabilityResult(const InteractionEnergy & energy) const { if (!hasBasePairProbabilityOutput()) return nullptr; std::unique_ptr handler(getSeedHandler(energy)); if (!handler->guaranteesStackOnlySeeds()) - throw std::invalid_argument("bpProb requires a stack-only seed handler; bulged seeds are unsupported"); - return new BasePairProbabilities(energy.size1(),energy.size2()); + throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); + return new BasePairProbabilities(energy.size1(),energy.size2(),!outPrefix2streamName.at(OP_bpsvg).empty()); } void CommandLineParsing::writeBasePairProbabilities(const BasePairProbabilities & result,const InteractionEnergy & energy) const { const auto & t=energy.getAccessibility1().getSequence(); const auto & q=energy.getAccessibility2().getAccessibilityOrigin().getSequence(); - BasePairProbabilityWriter::writeFile(getFullFilename(outPrefix2streamName.at(OutPrefixCode::OP_bpProb),&t,&q),result,t,q,outSep); + if (!outPrefix2streamName.at(OP_bpProb).empty()) + BasePairProbabilityWriter::writeFile(getFullFilename(outPrefix2streamName.at(OP_bpProb),&t,&q),result,t,q,outSep); + if (!outPrefix2streamName.at(OP_bpsvg).empty()) + BasePairProbabilityWriter::writeSvgFile(getFullFilename(outPrefix2streamName.at(OP_bpsvg),&t,&q),result,energy); } Predictor* @@ -2512,10 +2517,10 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output, BasePair std::unique_ptr pairSeeds; if (pairProbabilities) { if(model.val!='P' || mode.val!='M' || noSeedRequired || !rri.empty() || windowWidth.val!=0) - throw std::invalid_argument("bpProb requires --model=P --mode=M with stack-only seeds and no windows"); + throw std::invalid_argument("bpProb/bpsvg requires --model=P --mode=M with stack-only seeds and no windows"); pairSeeds.reset(getSeedHandler(energy)); if(!pairSeeds->guaranteesStackOnlySeeds()) - throw std::invalid_argument("bpProb requires a stack-only seed handler; bulged seeds are unsupported"); + throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); } // set up hub for prediction tracking (if needed) PredictionTrackerHub * predTracker = new PredictionTrackerHub(); diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index fff52480..f0216b37 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -289,7 +289,7 @@ class CommandLineParsing { Predictor* getPredictor( const InteractionEnergy & energy , OutputHandler & output, BasePairProbabilities * pairProbabilities = nullptr ) const; - /** @return whether actual base-pair matrix output was requested */ + /** @return whether actual base-pair CSV or SVG output was requested */ bool hasBasePairProbabilityOutput() const; /** Allocate a pair result after validating the actual seed capability. * @return owned result or nullptr when pair output is not requested @@ -416,6 +416,7 @@ class CommandLineParsing { OP_spotProb, OP_spotProbAll, OP_bpProb, + OP_bpsvg, OP_UNKNOWN }; @@ -442,6 +443,7 @@ class CommandLineParsing { if (prefLC == "tacc") { return OutPrefixCode::OP_tAcc; } else if (prefLC == "qpu") { return OutPrefixCode::OP_qPu; } else if (prefLC == "tpu") { return OutPrefixCode::OP_tPu; } else + if (prefLC == "bpsvg") { return OutPrefixCode::OP_bpsvg; } else if (prefLC == "bpprob") { return OutPrefixCode::OP_bpProb; } else if (prefLC == "spotprob") { return OutPrefixCode::OP_spotProb; } else // not known diff --git a/tests/BasePairProbabilities_test.cpp b/tests/BasePairProbabilities_test.cpp index 500c306c..669071cb 100644 --- a/tests/BasePairProbabilities_test.cpp +++ b/tests/BasePairProbabilities_test.cpp @@ -48,3 +48,44 @@ TEST_CASE("pair writer validates before emitting and checks stream failure", "[B BasePairProbabilityWriter::write(na,empty,t,q); REQUIRE(na.str()=="bpProb;U_1;C_2\nG_1;NA;NA\nA_2;NA;NA\n"); } + +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +TEST_CASE("seed annotations commit with regional masses and SVG requires success", "[BasePairProbabilities]") { + #include "testEasyLoggingSetup.icc" + RnaSequence t("target<&", "GGGG"),q("query", "CCCCC"); + AccessibilityDisabled at(t,0,nullptr),aq(q,0,nullptr); + ReverseAccessibility ar(aq); InteractionEnergyBasePair energy(at,ar); + BasePairProbabilities r(4,5,true); + Matrix mass(2,2,0);mass(0,1)=2;mass(1,0)=2; + Matrix seeds(2,2,0);seeds(0,1)=1; + r.addRegion({0,1},{1,2},2,mass,&seeds); + r.addRegion({2,3},{3,4},2,mass,&seeds); + REQUIRE(r.collectsSeedPairs());REQUIRE(r.seedPairs()(0,2)==1); + REQUIRE(r.seedPairs()(2,4)==1);REQUIRE(r.seedPairs()(1,1)==0); + std::ostringstream pending; + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(pending,r,energy)); + REQUIRE(pending.str().empty()); + r.finalize();std::ostringstream svg; + BasePairProbabilityWriter::writeSvg(svg,r,energy); + REQUIRE(svg.str().find("target<&")!=std::string::npos); + REQUIRE(svg.str().find("Orange outline:")!=std::string::npos); + std::ostringstream broken;broken.setstate(std::ios::badbit); + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(broken,r,energy)); + r.fail();std::ostringstream failed; + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(failed,r,energy)); + REQUIRE(failed.str().empty()); + BasePairProbabilities missing(4,5,true); + REQUIRE_THROWS(missing.addRegion({0,1},{1,2},2,mass)); + REQUIRE(missing.status()==BasePairProbabilities::Status::failed); + REQUIRE(missing.getZ()==0);REQUIRE(missing.seedPairs()(0,2)==0); + BasePairProbabilities plain(4,5); + REQUIRE_FALSE(plain.collectsSeedPairs());REQUIRE(plain.seedPairs().size1()==0); + plain.finalize();std::ostringstream empty; + BasePairProbabilityWriter::writeSvg(empty,plain,energy); + REQUIRE(empty.str().find("data-probability=\"NA\"")!=std::string::npos); + REQUIRE(empty.str().find("Orange outline:")==std::string::npos); + BasePairProbabilities wrongSize(1,1);wrongSize.finalize();std::ostringstream wrong; + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(wrong,wrongSize,energy)); + REQUIRE(wrong.str().empty()); +} diff --git a/tests/Makefile.am b/tests/Makefile.am index 5a919fb1..4d184d6b 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -85,4 +85,4 @@ runApiTests_CXXFLAGS = -I$(top_builddir)/src \ -EXTRA_DIST = tools/export-seeded-theory.py testBasePairProbabilities.py +EXTRA_DIST = tools/export-seeded-theory.py testBasePairProbabilities.py testBasePairProbabilitySvg.py diff --git a/tests/runBasePairProbabilities.sh b/tests/runBasePairProbabilities.sh index be6d002c..b68221c6 100644 --- a/tests/runBasePairProbabilities.sh +++ b/tests/runBasePairProbabilities.sh @@ -2,3 +2,4 @@ set -euo pipefail script_dir=$(CDPATH= cd -- "$(dirname -- "$0")" && pwd) python3 "$script_dir/testBasePairProbabilities.py" "$INTARNABINPATH/src/bin/IntaRNA" +python3 "$script_dir/testBasePairProbabilitySvg.py" "$INTARNABINPATH/src/bin/IntaRNA" diff --git a/tests/testBasePairProbabilitySvg.py b/tests/testBasePairProbabilitySvg.py new file mode 100644 index 00000000..29047ef3 --- /dev/null +++ b/tests/testBasePairProbabilitySvg.py @@ -0,0 +1,157 @@ +"""SVG semantics, original coordinates, seed provenance and output lifecycle.""" +import csv +import gzip +import math +from pathlib import Path +import re +import subprocess +import sys +import tempfile +import xml.etree.ElementTree as ET + +BIN = sys.argv[1] +COMMON = ['--energy=B', '--acc=N', '--model=P', '--mode=M', '--seedBP=2', + '--threads=1', '--outNoLP=false', '--outNoGUend=false', '--outNumber=0', + '--default-log-file=/dev/null'] +NS = {'s': 'http://www.w3.org/2000/svg'} + + +def run(t='GGAGGG', q='CCCC', options=(), success=True): + keys = {o.split('=')[0] for o in options} + p = subprocess.run([BIN, f'--target={t}', f'--query={q}', + *[o for o in COMMON if o.split('=')[0] not in keys], *options], + text=True, capture_output=True) + assert (p.returncode == 0) == success, (options, p.stdout, p.stderr) + return p + + +def load(path): + return ET.fromstring(gzip.decompress(path.read_bytes()) if str(path).endswith('.gz') else path.read_bytes()) + + +def pairs(root): + return {(int(r.attrib['data-target']), int(r.attrib['data-query'])): r + for r in root.findall('.//s:rect[@data-type="base-pair"]', NS)} + + +def seeds(root): + return {(int(r.attrib['data-target']), int(r.attrib['data-query'])) + for r in root.findall('.//s:rect[@class="seed"][@data-target]', NS)} + + +with tempfile.TemporaryDirectory(prefix='intarna-bpsvg-') as tmp: + d = Path(tmp) + svg, table = d/'pairs.svg', d/'pairs.csv' + + def plot(t='GGAGGG', q='CCCC', opts=()): + run(t, q, [f'--out=bpsvg:{svg}', *opts]) + return load(svg) + + root = plot(opts=[f'--out=bpProb:{table}']) + cells = pairs(root) + rows = list(csv.reader(table.open(), delimiter=';')) + assert len(cells) == 24 + assert '+1' in [e.text for e in root.findall('.//s:text[@class="index"]', NS)] + assert not root.findall('.//s:line[@data-index="1"]', NS) + for i, row in enumerate(rows[1:], 1): + for j, value in enumerate(row[1:], 1): + r = cells[i, j] + assert float(r.attrib['data-probability']) == float(value) + bounds = [0.01, 0.1, 0.25, 0.5, 0.75, 0.9] + expected_color = sum(float(value) >= boundary for boundary in bounds) + assert f'p{expected_color}' in r.attrib['class'].split() + assert 'Base-pair probability' in r.find('s:title', NS).text + assert all(float(cells[3, j].attrib['data-probability']) == 0 for j in range(1, 5)) + frame = root.findall('.//s:rect[@data-type="unpaired"]', NS) + assert len(frame) == 2*(6+4) + assert all(float(r.attrib['data-probability']) == 1 for r in frame) + assert all('Unpaired probability' in r.find('s:title', NS).text for r in frame) + css = root.find('s:style', NS).text + assert all(f'.p{i} {{ fill:' in css for i in range(7)) + assert len(root.findall('.//s:g[@class="probability-legend"]/s:rect', NS)) == 7 + # All admitted computed two-pair stacks, in original query orientation. + expected = {(i+k+1, 4-j-k) for i in (0, 3, 4) for j in range(3) for k in range(2)} + assert seeds(root) == expected, seeds(root) + # Requesting SVG must not affect CSV or ensemble/ranked output. + baseline = run(options=['--outNumber=1', '--outMode=C', '--outCsvCols=Zall,Eall']) + plotted = run(options=['--outNumber=1', '--outMode=C', '--outCsvCols=Zall,Eall', f'--out=bpsvg:{svg}']) + assert baseline.stdout == plotted.stdout + # Explicit singleton/mixed seed masks exclude extension-only cells. + opts = ['--seedTQ=1|&4|,2||&2||'] + mixed = plot('GGGG', 'CCCC', opts) + assert seeds(mixed) == {(1,4), (2,3), (3,2)} + restricted = plot('GGGG', 'CCCC', [*opts, '--tRegion=2-4', '--qRegion=1-3']) + assert seeds(restricted) == {(2,3), (3,2)} + # A seed crossing a searched region boundary is not annotated. + cut = plot('GGGG', 'CCCC', ['--seedTQ=2||&2||', '--tRegion=1-2,3-4']) + assert not seeds(cut) + assert all(r.attrib['data-probability'] == 'NA' for r in pairs(cut).values()) + # Four disjoint rectangles merge both probabilities and masks. + regional = plot('GGGG', 'CCCC', ['--tRegion=1-2,3-4', '--qRegion=1-2,3-4']) + assert seeds(regional) == {(1,2),(2,1),(1,4),(2,3),(3,2),(4,1),(3,4),(4,3)} + # Signed guides follow displayed coordinates; there is no zero nucleotide. + signed = plot('A'*72, 'A'*18, ['--tIdxPos0=-15', '--qIdxPos0=-15']) + signed_cells = pairs(signed) + for strand, length in (('target',72), ('query',18)): + guides = signed.findall(f'.//s:line[@data-strand="{strand}"]', NS) + expected_indices = {-10, 1} | ({10,20,30,40,50} if length == 72 else set()) + assert {int(g.attrib['data-index']) for g in guides} == expected_indices + assert not any(0 in k for k in signed_cells) + for g in guides: + idx = int(g.attrib['data-index']) + r = signed_cells[idx,-15] if strand == 'target' else signed_cells[-15,idx] + axis = 'y' if strand == 'target' else 'x' + expected_pos = float(r.attrib[axis]) + (20 if idx > 1 else 0) + assert float(g.attrib[axis+'1']) == expected_pos + assert ('major' in g.attrib['class']) == (idx % 50 == 0) + assert ('origin' in g.attrib['class']) == (idx == 1) + negative = plot('A'*13, 'A'*13, ['--tIdxPos0=-55', '--qIdxPos0=45']) + assert 'major' in negative.find('.//s:line[@data-strand="target"][@data-index="-50"]', NS).attrib['class'] + # User sequence identifiers must round-trip through XML escaping. + name = '' + named = plot(opts=[f'--tId={name}']) + assert name in named.find('s:title', NS).text + assert name+" (5' to 3')" in [e.text for e in named.findall('.//s:text', NS)] + # Native accessibility values must use original (not reversed) query order. + pu_t, pu_q = d/'target.pu', d/'query.pu' + native = plot('GGGGAAAACCCC', 'GCGCGAAAACGCGC', + ['--energy=V', '--acc=C', f'--out=tPu:{pu_t}', f'--out=qPu:{pu_q}']) + for strand, path in (('target', pu_t), ('query', pu_q)): + expected_pu = [float(line.split()[1]) for line in path.read_text().splitlines() + if line.strip() and not line.lstrip().startswith('#')] + assert any(0 < v < .9 for v in expected_pu) + for r in native.findall(f'.//s:rect[@data-type="unpaired"][@data-strand="{strand}"]', NS): + assert math.isclose(float(r.attrib['data-probability']), expected_pu[int(r.attrib['data-index'])-1], + rel_tol=2e-5, abs_tol=1e-6) + for t,q in (('AAAA','AAAA'), ('G','C')): + empty = plot(t,q) + assert all(r.attrib['data-probability'] == 'NA' for r in pairs(empty).values()) + assert not seeds(empty) + assert 'empty interaction ensemble' in ''.join(empty.itertext()) + compressed = d/'pairs.svg.gz' + run(options=[f'--out=bpsvg:{compressed}']) + assert len(pairs(load(compressed))) == 24 + # Validation/numerical failures must not publish a partial document. + for opts in (['--mode=H'], ['--noSeed'], ['--seedMaxUP=1'], ['--model=X', '--mode=K'], + ['--windowWidth=3'], ['--seedTQ=1|&4|','--energyAdd=900']): + failure = d/'failure.svg' + run(options=[f'--out=bpsvg:{failure}', *opts], success=False) + assert not failure.exists() + run(options=['--out=bpsvg:/dev/full'], success=False) + run(options=[f'--out=bpsvg:{svg}', f'--out=bpProb:{svg}'], success=False) + # Multi-FASTA creates one SVG per pair; shared streams emit indivisible roots. + fasta=d/'targets.fa'; fasta.write_text('>a\nGGGG\n>b\nGAGG\n') + for threads in (1,2): + prefix=d/f'multi{threads}.svg' + run(str(fasta), 'CCCC', [f'--threads={threads}', f'--out=bpsvg:{prefix}']) + files=sorted(d.glob(f'multi{threads}*.svg')) + assert len(files)==2 + assert all(len(pairs(load(file)))==16 for file in files) + for stream in ('STDOUT', 'STDERR'): + shared=run(str(fasta), 'CCCC', [f'--threads={threads}', f'--out=bpsvg:{stream}']) + data=shared.stdout if stream=='STDOUT' else shared.stderr + blocks=re.findall(r'', data, re.S) + assert len(blocks)==2 + assert all(len(pairs(ET.fromstring(block)))==16 for block in blocks) + +print('base-pair probability SVG regressions passed') From fa2f84d6cdbb1e120a98eab9d2dfe3f8154ebad5 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Fri, 9 Oct 2026 12:11:59 +0200 Subject: [PATCH 2/9] docs: depict ensemble and pair-probability recurrences (step 9) --- ChangeLog | 13 ++++ doc/recursions/IntaRNA.PredictorMfeEns2d.svg | 60 +++++++++++++++++ .../IntaRNA.PredictorMfeEns2dHeuristic.svg | 51 ++++++++++++++ ...redictorMfeEns2dHeuristicSeedExtension.svg | 51 ++++++++++++++ ....PredictorMfeEns2dSeedExtension.legacy.svg | 61 +++++++++++++++++ ...IntaRNA.PredictorMfeEns2dSeedExtension.svg | 54 +++++++++++++++ .../IntaRNA.PredictorMfeEnsSeedOnly.svg | 41 ++++++++++++ ...IntaRNA.SeededPartitionFunction.inside.svg | 59 ++++++++++++++++ ...ntaRNA.SeededPartitionFunction.outside.svg | 67 +++++++++++++++++++ 9 files changed, 457 insertions(+) create mode 100644 doc/recursions/IntaRNA.PredictorMfeEns2d.svg create mode 100644 doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg create mode 100644 doc/recursions/IntaRNA.PredictorMfeEns2dHeuristicSeedExtension.svg create mode 100644 doc/recursions/IntaRNA.PredictorMfeEns2dSeedExtension.legacy.svg create mode 100644 doc/recursions/IntaRNA.PredictorMfeEns2dSeedExtension.svg create mode 100644 doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg create mode 100644 doc/recursions/IntaRNA.SeededPartitionFunction.inside.svg create mode 100644 doc/recursions/IntaRNA.SeededPartitionFunction.outside.svg diff --git a/ChangeLog b/ChangeLog index d04228f5..613764ff 100644 --- a/ChangeLog +++ b/ChangeLog @@ -50,6 +50,9 @@ ## Technical changes and Optimizations +- depict ensemble predictor recurrences and the exact base-pair inside/outside + computation in internal SVG documentation (PR #258, step 9) + - validate exact pair probabilities on independent mathematical/native oracles, release/debug suites and installed consumers; document measured runtime and storage, and supply the required logging link dependency (issue #257, step 7) @@ -160,6 +163,16 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261009 Alexander Mitrofanov + * doc/recursions/IntaRNA.PredictorMfeEns*.svg, + doc/recursions/IntaRNA.SeededPartitionFunction.*.svg : + + add eight standalone internal SVG depictions for unseeded exact/heuristic, + seed-only, exact stack-seeded, legacy and heuristic seed-extension paths + + show noLP transitions, seed admission and energy ownership, the production + suffix inside/outside recurrences, actual-pair marking and regional merging + + include implementation anchors and retain editable vector text; validate + XML, glyph coverage and rendered layouts (PR #258 step 9, PR #259) + 261009 Alexander Mitrofanov * BasePairProbabilityWriter, BasePairProbabilitySvg, CommandLineParsing : + add success-only SVG dot plots on files, gzip and synchronized standard diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2d.svg b/doc/recursions/IntaRNA.PredictorMfeEns2d.svg new file mode 100644 index 00000000..674ba774 --- /dev/null +++ b/doc/recursions/IntaRNA.PredictorMfeEns2d.svg @@ -0,0 +1,60 @@ + +Unseeded ensemble · fixed right boundary +PredictorMfeEns2d | --model=P --mode=M --noSeed. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Unseeded ensemble · fixed right boundary +PredictorMfeEns2d | --model=P --mode=M --noSeed + +Pair indices i=(i₁,i₂), k and j increase on both strands; RNA 2 is reversed. +i ≺ k means strict increase in both coordinates. d=BW(E_init), t(i,k)=BW(E_interLeft(i,k)). +Hⱼ(i) = hybrid partition for all admitted chains from i to a fixed j + +Initiation + +Direct pair-to-pair step + +Recursive continuation + + + +i = j +d + + + +i + +j +t(i,j) · d + +Hⱼ(k) + + + +i + +k + +j +t(i,k) · Hⱼ(k) +Lonely pairs allowed +Hⱼ(j)=d; Hⱼ(i)=t(i,j)d + ∑ t(i,k)Hⱼ(k) (i ≺ k ≺ j) +All other or noncomplementary states are zero. Sum only loop-size-admissible steps. +For every valid right pair j, fill i in decreasing order within the span limits. + +noLP: the left boundary must stack to i⁺ = i + (1,1) +Hⱼ(j)=0; Hⱼ(i)=t(i,i⁺) · [ (d if i⁺=j, else Hⱼ(i⁺)) + ∑ t(i⁺,k)Hⱼ(k) ] +The sum requires i⁺ ≺ k ≺ j and k ≠ i⁺+(1,1): the direct stack is already in Hⱼ(i⁺). +Require a complementary i⁺. Unreachable suffixes are zero; the adjacent two-pair case contributes t(i,j)d. +Thus initialization is a stack, and each nonstack continuation starts another stacked suffix. +Zall += B(i,j) · Hⱼ(i); E_site = energy.getE(i,j, energy.getE(Hⱼ(i))) +B(i,j)=BW(energy.getE(i₁,j₁,i₂,j₂,0)) for admitted boundaries; otherwise B=0. +This adds full-interval ED and terminal terms once; updateCompleteZ submits each boundary once. + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +PredictorMfeEns2d.cpp: predict(), fillHybridZ() +PredictorMfeEns.cpp: updateCompleteZ(), addPartitionContribution() + diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg new file mode 100644 index 00000000..fd048c48 --- /dev/null +++ b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg @@ -0,0 +1,51 @@ + +Unseeded ensemble · heuristic continuation +PredictorMfeEns2dHeuristic | --model=P --mode=H --noSeed. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Unseeded ensemble · heuristic continuation +PredictorMfeEns2dHeuristic | --model=P --mode=H --noSeed + +C(i)=(vᵢ,jᵢ) retains one hybrid weight and right boundary for each left pair i. +Use d=BW(E_init), t(i,k)=BW(E_interLeft(i,k)); all pair coordinates increase internally. + +Enumerate candidates at i, visiting i in decreasing order + +vₖ + + + +i + +k + +jₖ +Lonely pairs allowed: + initial candidate (d, i) + continuation (t(i,k) · vₖ, jₖ) for each reachable k ≻ i + within loop limits and the complete span limits. + + +Accumulate every considered candidate +updateZ(i,j,v,true) +Adds B(i,j)v to Zall and v to the boundary map. +B is the admitted complete-boundary Boltzmann factor. + +Retain only the best continuation at i +C(i) ← argmin E_full(i,j,v) +Compare complete energy for each candidate. +Strict improvement replaces vᵢ and jᵢ. + +noLP candidate set +Let i⁺=i+(1,1) be complementary and s=t(i,i⁺). +Initialize (s·d, i⁺); add direct continuation (s·vᵢ⁺, jᵢ⁺). +Other candidates are (s·t(i⁺,k)·vₖ, jₖ), with k ≻ i⁺ and k ≠ i⁺+(1,1). +The adjacent continuation is excluded from the loop sum because it was handled directly. +The boundary map sums contributions, but the DP cell keeps one right continuation: this is a restricted ensemble. + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +PredictorMfeEns2dHeuristic.cpp: fillHybridZ() +PredictorMfeEns.cpp: updateZ(), updateOptimaUsingZ() + diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristicSeedExtension.svg b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristicSeedExtension.svg new file mode 100644 index 00000000..6f7c7fa8 --- /dev/null +++ b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristicSeedExtension.svg @@ -0,0 +1,51 @@ + +Seeded ensemble · one selected right extension +PredictorMfeEns2dHeuristicSeedExtension | --model=P --mode=H. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Seeded ensemble · one selected right extension +PredictorMfeEns2dHeuristicSeedExtension | --model=P --mode=H + +For anchor seed a…b: S=BW(E_seed), d=BW(E_init), L=Lₐ(i), R=R_b(j). +L and R use the inherited legacy extension recurrences, even when the seed handler is stack-only. + +1 Compute R, submit true right extensions, and select j* +updateZ(a,j, d · S · R_b(j), true) for j ≻ b +j* = argmin E_full(a,j, d · S · R_b(j)) over right-end-admissible candidates +updateOptRightZ compares complete energy and the right GU-end restriction; ties retain the first candidate. + + +2 Compute corrected L and submit left-side combinations + +L + + + + +i + +a + + +b + +j* +updateZ(i,b, Lₐ(i) · S, true) +Includes the seed-only case i=a. +updateZ(i,j*, Lₐ(i) · S · R_b(j*), true) +Only for a true left extension and a true selected right extension, +with both complete spans within their limits. + +Candidate family per anchor +All true right-only extensions + all corrected left-only extensions + left extensions paired with just j*. +The second phase runs if noGUend is off, or a valid right optimum exists, or anchor end b is not GU. +Initialize j*=b and E_right_opt=∞. A positive selected R is required before combining both sides. +No full Cartesian product of left and right boundaries is retained, and bpProb / bpsvg are unsupported. +Every update is filtered and accumulated by PredictorMfeEns::updateZ; ranked reporting is a separate step. + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +PredictorMfeEns2dHeuristicSeedExtension.cpp: predict(), fillHybridZ_right(), fillHybridZ_left() +PredictorMfeEns2dHeuristicSeedExtension.h: updateOptRightZ(); inherited legacy fillHybridZ kernels + diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2dSeedExtension.legacy.svg b/doc/recursions/IntaRNA.PredictorMfeEns2dSeedExtension.legacy.svg new file mode 100644 index 00000000..f19dc623 --- /dev/null +++ b/doc/recursions/IntaRNA.PredictorMfeEns2dSeedExtension.legacy.svg @@ -0,0 +1,61 @@ + +Legacy seeded extensions · anchor factorization +Bulged-seed fallback in mode M; also the extension kernels used by the seeded mode H predictor. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Legacy seeded extensions · anchor factorization +Bulged-seed fallback in mode M; also the extension kernels used by the seeded mode H predictor + +For one retained anchor seed a…b, use S=BW(E_seed(a)), d=BW(E_init), t(i,k)=BW(E_interLeft(i,k)). +Lₐ(i) extends left to a (includes initiation); R_b(j) extends right from b (excludes initiation). + +Lₐ(i) + + + + +i + +a + + + + +b + +j + + +R_b(j) + + + + +hybrid contribution = Lₐ(i) · S · R_b(j) + +Lonely pairs allowed +Lₐ(a)=d; Lₐ(i)=∑ t(i,k)Lₐ(k) − Cₐ(i) (i ≺ k ≼ a) +R_b(b)=1; R_b(j)=∑ R_b(k)t(k,j) (b ≼ k ≺ j) +Only nonzero, complementary states and loop-size-admissible terms participate. Sums use componentwise bounds. + +noLP: outward boundary stacks are explicit +Lₐ(i)=t(i,i⁺) · [Lₐ(i⁺) + ∑ t(i⁺,k)Lₐ(k)] − Cₐ(i) (i⁺ ≺ k ≼ a) +R_b(j)=t(j⁻,j) · [1{j⁻=b}R_b(b) + ∑ R_b(k)t(k,j⁻)] (b ≼ k ≺ j⁻) +i⁺=i+(1,1), j⁻=j−(1,1); require the extra stacked pair. Anchor bases still have Lₐ(a)=d and R_b(b)=1. +These are the legacy loop bounds. The exact stack-seed path instead uses the suffix-state recurrence. + +Left-seed correction Cₐ(i) +At a competing seed start i: process it if its end is left of a, or it loop-overlaps the anchor. +If a leftmost loop-overlapping seed h is found, subtract BW(getNonOverlappingEnergy(i,h)) · Lₐ(h). +Otherwise subtract its seed weight times the continuation at its end; noLP adds explicit loop continuations. +Apply the implementation’s overlap/loop guards, then clamp a negative corrected L value to zero. +updateZ(i,j, Lₐ(i) · S · R_b(j), true) for each nonzero span-admissible combination +Complete-site filters and B(i,j) are applied by updateZ. The anchor’s hybrid seed energy is charged once. +This legacy subtraction-based path has no actual-pair outside pass; do not substitute it for the exact seeded union. + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +PredictorMfeEns2dSeedExtension.cpp: predict() legacy branch, fillHybridZ_left(), fillHybridZ_right() +PredictorMfeEns2dSeedExtension.cpp: getNonOverlappingEnergy(); PredictorMfeEns.cpp: updateZ() + diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2dSeedExtension.svg b/doc/recursions/IntaRNA.PredictorMfeEns2dSeedExtension.svg new file mode 100644 index 00000000..b04a6fe1 --- /dev/null +++ b/doc/recursions/IntaRNA.PredictorMfeEns2dSeedExtension.svg @@ -0,0 +1,54 @@ + +Seeded ensemble · exact stack-seed adapter +PredictorMfeEns2dSeedExtension | --model=P --mode=M with stack-only seed capability. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Seeded ensemble · exact stack-seed adapter +PredictorMfeEns2dSeedExtension | --model=P --mode=M with stack-only seed capability + + +SeedHandler: guaranteesStackOnlySeeds() +true → predictStackSeeds(): unique seeded-chain partition +false → legacy anchor extensions; no pair-probability sink + +Capability is geometric +Computed or explicit seeds; mixed lengths and singletons. +Preserve handler admission and per-start selection rules. + +1 Compile active seed occurrences and eligible left boundaries + + + + +p + +a + + +b + +q +max(0, bₖ − Wₖ + 1) ≤ pₖ ≤ aₖ (k=1,2) +Admitted stack a…b must fit both the region and span limits W. +Enumerate the union of these p-boxes once, then require a valid p. +At p, clip the working grid to both inclusive spans. + +2 Give the suffix kernel one consistent energy objective +d(p)=BW(E_init); t(r,q)=BW(E_interLeft(r₁,q₁,r₂,q₂)) +B(p,q)=BW(energy.getE(p₁,q₁,p₂,q₂,0)) if the complete site is admitted +Invalid / inaccessible / noncomplementary vertices and forbidden transitions cannot contribute. +B owns complete-span ED, native dangles, terminal penalties and energyAdd. Seed energies only admit seeds. + +3 Commit a completed region +Z_region = ∑ B(p,q) · T₁,p(q) +Complete boundary callback and outside initialization share exactly the same B. +updateExactCompleteZ preserves the partition objective. +Pair masses and optional seed masks map through energy.getBasePair(). +Only a successful computation reaches BasePairProbabilities::addRegion(). + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +PredictorMfeEns2dSeedExtension.cpp: predict(), predictStackSeeds(), exactBoundaryWeight() +StackSeedDomain.cpp; SeededPartitionFunction.cpp; PredictorMfeEns.cpp: updateExactCompleteZ() + diff --git a/doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg b/doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg new file mode 100644 index 00000000..0a10ebe9 --- /dev/null +++ b/doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg @@ -0,0 +1,41 @@ + +Seed-only ensemble · retained seed enumeration +PredictorMfeEnsSeedOnly | --model=P --mode=S. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Seed-only ensemble · retained seed enumeration +PredictorMfeEnsSeedOnly | --model=P --mode=S + + +1 Fill the seed handler for the active region +Retain the handler-selected seed at each admitted start a; its end is e. +Both endpoints must lie in the region. The handler supplies seed pairs and hybrid energy E_seed(a). +A filled handler with no admitted seeds reports an empty result. + + +2 One complete contribution per retained seed + + + + +a + + + + +e +H_seed = E_seed(a) + E_init +E_full = energy.getE(a₁,e₁,a₂,e₂,H_seed) +Zall += BW(E_full) when needZall is enabled + +3 Report and trace the retained seed +updateOptima(..., isHybridE=true, incrementZall=true) applies the complete-site filters and ranking. +traceBackSeed() supplies its interior pairs; initiation is charged exactly once. +There is no extension recurrence, no sum over alternative seed structures at the same start, and no bpProb output. + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +PredictorMfeEnsSeedOnly.cpp: predict(), traceBack() +PredictorMfe.cpp: updateOptima(); Predictor.h: updateZall() (inherited update path) + diff --git a/doc/recursions/IntaRNA.SeededPartitionFunction.inside.svg b/doc/recursions/IntaRNA.SeededPartitionFunction.inside.svg new file mode 100644 index 00000000..b944efbb --- /dev/null +++ b/doc/recursions/IntaRNA.SeededPartitionFunction.inside.svg @@ -0,0 +1,59 @@ + +Exact seeded inside recurrence · suffix automaton +SeededPartitionFunction::compute() | fixed left boundary p; visit q in increasing pair order. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Exact seeded inside recurrence · suffix automaton +SeededPartitionFunction::compute() | fixed left boundary p; visit q in increasing pair order + +C=max(2, longest admitted seed length). ℓ is the current consecutive stack-suffix length, capped at C. +e(q)=shortest admitted seed ending at q, or 0. A run accepts a seed iff 0 < e(q) ≤ ℓ. + +Uℓ: still seed-free +ℓ=1…C; no seed has occurred. +Stack increments ℓ; a loop resets it. + +S₁: seeded, singleton suffix +A seed has occurred; current run is length 1. +Under noLP, must stack before leaving. + +S₂: seeded, usable suffix +Under noLP: current run has length ≥2. +Otherwise: all seeded suffixes use S₂. + + + +first seed accepted on a stack +stack → S₂ + + +1 Initialize and extend every state +H(p,c₀)=d(p); other H=0. c₀ = accept(p,1) ? (noLP ? S₁ : S₂) : U₁ +Stack r=q−(1,1): H(q,D(c,stack)) += t(r,q) · H(r,c) +For Uℓ: let ℓ′=min(C,ℓ+1); D=Uℓ′ unless the ending seed is accepted, then D=S₂. +For S₁ or S₂: D=S₂. Once seeded, a later seed never creates a second copy of the path. + +2 Aggregate eligible departures, then handle nonstack edges +T₀(q)=∑ H(q,Uℓ) over ℓ≥(noLP ? 2 : 1) +T₁(q)=H(q,S₂) + (noLP ? 0 : H(q,S₁)) +For each nonstack r→q and b∈{0,1}: H(q,D(b,loop)) += t(r,q) · T_b(r) +Nonstack means both coordinates increase, within loop limits, excluding r=q−(1,1). +Reset the suffix to 1; remain seeded if b=1, or accept a singleton seed at q if b=0. + +3 Complete seeded paths at every q +Z += B(p,q) · T₁(q) +Seed-free states never terminate an interaction. +With noLP, singleton runs cannot terminate or cross a nonstack edge. +Such states remain available for a later stack. +Loop destinations: seeded → S₁ under noLP, otherwise S₂; unseeded → U₁ unless a singleton seed accepts. +The Boolean seeded state counts each ordered pair chain once, including paths with several seed witnesses. +Termination-only GU / ED filters belong to B; they do not discard prefixes that can still extend. +d, t and B are the initiation, transition and complete-boundary factors supplied by Weights. + + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +SeededPartitionFunction.cpp: accepts, dest, nonstack, h, total; StackSeedDomain::seedEndingLength() +Code states: Uℓ ↔ c=ℓ−1; S₁ ↔ c=cap; S₂ ↔ c=cap+1. Checked nonnegative arithmetic throughout. + diff --git a/doc/recursions/IntaRNA.SeededPartitionFunction.outside.svg b/doc/recursions/IntaRNA.SeededPartitionFunction.outside.svg new file mode 100644 index 00000000..c98b9271 --- /dev/null +++ b/doc/recursions/IntaRNA.SeededPartitionFunction.outside.svg @@ -0,0 +1,67 @@ + +Actual base-pair probabilities · outside recurrence +Reverse differentiation of the suffix inside computation; same fixed p and same boundary weights B. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Actual base-pair probabilities · outside recurrence +Reverse differentiation of the suffix inside computation; same fixed p and same boundary weights B + +H and T are forward state/aggregate weights. Their adjoints are A=∂Z/∂H and G=∂Z/∂T. +M(q) is the raw mass of complete chains containing actual pair q; it is not site-coverage mass. + +H(r,c) + + + +p + +r + +q + +terminal + +Reverse the incoming edge r → q. +Only actual pair q receives this edge’s mass. + + +1 Initialize the reverse objective +A(q,S₂)=B(p,q); A(q,S₁)=B(p,q) only without noLP; other A and G = 0 +Initialize at reachable positive seeded totals; a zero terminal coefficient contributes nothing. + +2 Reverse q in decreasing order; first scatter aggregate adjoints +A(q,Uℓ) += G₀(q) for the same departure-eligible ℓ as in the forward sum +A(q,S₂) += G₁(q); A(q,S₁) += G₁(q) only without noLP +This scatter must precede reversing the incoming edges at q. + +3a Reverse a stack edge r→q +a = A(q,D(c,stack)) · t(r,q) +A(r,c) += a +M(q) += a · H(r,c) +Repeat for every forward source state c. +The transition marker belongs to destination q. + +3b Reverse a nonstack edge r→q +a = A(q,D(b,loop)) · t(r,q) +Gᵦ(r) += a +M(q) += a · T_b(r) +Repeat for b=0 and b=1 using the same forward totals. +Do not mark the unpaired nucleotides in the loop. + +4 Mark the source pair after the reverse sweep +M(p) += A(p,c₀) · d(p) +Every other pair was marked when its incoming edge was reversed. + +5 Merge successful disjoint search regions before normalization +Z_total = ∑ Z_region; M_total(i,j) = ∑ M_region(i,j); P(i,j)=M_total(i,j)/Z_total +Map regional pairs with energy.getBasePair() into original target/query coordinates, then finalize explicitly. +Z_total=0 → NA. Failed/cancelled computation → no matrix. No unbound-state weight is added. +Seed admission and B are fixed while differentiating. Masses sum across every eligible left boundary p. + + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +SeededPartitionFunction.cpp: dh (A), dt (G), mass and initial-pair contribution +PredictorMfeEns2dSeedExtension::predictStackSeeds(); BasePairProbabilities::addRegion(), finalize(), probabilities() + From 202651e792e01ed188c0253a5251decb193a2b77 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Fri, 9 Oct 2026 12:29:07 +0200 Subject: [PATCH 3/9] fix: address probability output and recurrence audit findings --- ChangeLog | 24 ++++ README.md | 8 +- .../seeded-probabilities-validation.md | 113 ++++++++++++++++++ .../IntaRNA.PredictorMfeEns2dHeuristic.svg | 9 +- .../IntaRNA.PredictorMfeEnsSeedOnly.svg | 3 +- ...IntaRNA.SeededPartitionFunction.inside.svg | 12 +- ...ntaRNA.SeededPartitionFunction.outside.svg | 2 +- src/IntaRNA/BasePairProbabilitySvg.cpp | 25 ++-- src/IntaRNA/BasePairProbabilityWriter.cpp | 3 +- src/IntaRNA/BasePairProbabilityWriter.h | 12 +- src/bin/CommandLineParsing.cpp | 16 ++- src/bin/CommandLineParsing.h | 19 ++- tests/BasePairProbabilities_test.cpp | 19 ++- tests/testBasePairProbabilitySvg.py | 53 ++++++-- 14 files changed, 267 insertions(+), 51 deletions(-) diff --git a/ChangeLog b/ChangeLog index 613764ff..44268e50 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,9 @@ ## Interface and handling +- BUGFIX : use each strand's accessibility energy scale for SVG frames and + unpaired-probability exports, including mixed computed/imported inputs + - --out=bpsvg:FILE renders actual-pair probabilities with nucleotide/accessibility frames, a shared CSS palette, signed-coordinate guides, hover values and admitted seed outlines (PR #258, step 8) @@ -50,6 +53,9 @@ ## Technical changes and Optimizations +- transfer completed probability-output buffers without copying; independently + audit performance, accessibility conversion and recurrence depictions (PR #259) + - depict ensemble predictor recurrences and the exact base-pair inside/outside computation in internal SVG documentation (PR #258, step 9) @@ -163,6 +169,24 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261009 Alexander Mitrofanov + * BasePairProbabilityWriter, BasePairProbabilitySvg, CommandLineParsing : + * pass and validate separate target/query accessibility RT values for SVG + frames; imported probabilities no longer inherit the interaction model RT + * share the source-specific conversion with tPu/qPu, correcting their earlier + physical-temperature conversion of computed base-pair accessibility (RT=1) + * move completed SVG/CSV stream buffers instead of copying large documents + * tests, README.md, doc/analysis/seeded-probabilities-validation.md : + + cover imported Pu/ED, mixed sources, nondefault temperatures and invalid RT; + compare against input probabilities and independent ED conversion + + document independent review, controlled runtime/RSS comparisons and final + release/debug, installation and public-header checks (PR #259) + * doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg, + doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg, + doc/recursions/IntaRNA.SeededPartitionFunction.*.svg : + * show the heuristic initial right-end GU gate and complete-site filters; + place nonstack arrivals before inside aggregation and clarify outside notation + 261009 Alexander Mitrofanov * doc/recursions/IntaRNA.PredictorMfeEns*.svg, doc/recursions/IntaRNA.SeededPartitionFunction.*.svg : diff --git a/README.md b/README.md index fcb80537..d627414e 100644 --- a/README.md +++ b/README.md @@ -1893,6 +1893,9 @@ The column separator within tabular CSV output (defaulting to `;`) can be change using `--outSep`, e.g. to produce tab-separated `.tsv` output. Note further, `qPu:`|`tPu:` will report unpaired probability values based on rounded accessibility (ED) values. +The ED-to-probability conversion uses each strand's accessibility energy scale: +computed base-pair accessibility uses RT=1, while ViennaRNA and imported +accessibility use the configured temperature, independently of the interaction model. Thus, these values will most likely differ from values eg. produced by the program RNAplfold. We therefore strongly recommend to store `qAcc:`|`tAcc:` values when you want to use them as input for subsequent IntaRNA calls! @@ -2116,7 +2119,10 @@ a standalone SVG on standard output, redirect ordinary output, for example produce consecutive SVG documents; use filenames for separate viewable files. Library callers can use `BasePairProbabilityWriter::writeSvg()` or -`writeSvgFile()` with a finalized result and the original energy model. Construct +`writeSvgFile()` with a finalized result, the original energy model, and separate +target/query accessibility RT values. These must match the ED conversion scales +of the accessibility sources, which can differ from the interaction model's RT +(for example, imported Pu with a base-pair interaction model). Construct `BasePairProbabilities(targetLength, queryLength, true)` to retain seed annotations; the default CSV-only result does not allocate that mask. SVG creation buffers the document and uses one rectangle per pair plus the frames, so its size grows diff --git a/doc/analysis/seeded-probabilities-validation.md b/doc/analysis/seeded-probabilities-validation.md index a08848d0..3449a5ba 100644 --- a/doc/analysis/seeded-probabilities-validation.md +++ b/doc/analysis/seeded-probabilities-validation.md @@ -157,3 +157,116 @@ OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 /usr/bin/time -v src/bin/IntaRNA \ Both runs report Eall=-6.52 kcal/mol with the native presentation convention. These measurements quantify added output work; they are not a speedup claim over the legacy seeded ensemble, whose multi-anchor partition is a different result. + + +## Steps 8–9: independent output and documentation audit (2026-10-09) + +Separate agents reviewed code correctness, recurrence depictions, and performance +before the audit fixes were pushed to PR #259. The audit found and corrected: + +- SVG accessibility frames used the interaction model's RT for imported Pu/ED. + The writer now requires separate target/query accessibility RT values. CLI + conversion selects RT=1 for computed base-pair accessibility and the configured + physical temperature for ViennaRNA/imported accessibility. The same selection + fixes the pre-existing `tPu`/`qPu` conversion of computed base-pair ED. +- SVG and CSV stream buffers were copied when returned; moving the completed + strings removes that duplicate allocation without changing document bytes. +- The heuristic depiction omitted its initial right-end GU gate, and the inside + depiction placed aggregation before incoming nonstack contributions. These + were corrected, along with complete-site filtering and outside-index notation. + +The added CLI regressions fail on the pre-audit binary and pass after the fix. +They cover Pu/ED imports, mixed computed/imported strands, 22/37°C, original +query orientation, independent input-Pu round trips (allowing ED quantization), +and computed-base-pair ED conversion at RT=1. API tests reject nonpositive or +nonfinite accessibility scales without writing output. GCC 14.4 release and +debug `make tests -j2` each pass all seven suites, with 91,394 API assertions in +86 cases. Installed public headers and an installed SVG consumer compile; the +consumer produces the expected matrix and seed annotations. Corrected SVG +recurrence documents parse, render, and were visually inspected. Apple Clang +was not available locally. + +### Performance method and existing-mode comparison + +Compared base `48efc541e23571ce1bab88e2754e1db7445cbc89` with steps 8–9 at +`fa2f84d6cdbb1e120a98eab9d2dfe3f8154ebad5`, then compared that saved binary with +the audit fixes. All builds use identical GCC 14.4 C++23 `-O3 +-fno-strict-aliasing -fopenmp`, ViennaRNA 2.7.2, Boost 1.85 and Kokkos mdspan +settings on Linux x86-64 / Ryzen 5 7530U. Processes are pinned to CPU 0 with +`--threads=1`, `OMP_NUM_THREADS=1`, `OPENBLAS_NUM_THREADS=1`, and +`OMP_DYNAMIC=FALSE`. Builds and tests are stopped during measurements. + +There is one warmup per configuration, then seven fresh-process runs in shuffled, +interleaved order for the base comparison. Time is wall time around the GNU time +subprocess; RSS is GNU time `%M`, including executable/libraries, prediction and +serialization. Output goes to the same local filesystem. Initial median seconds: + +| Input | Base ensemble only | Steps 8–9 ensemble only | Base CSV | Steps 8–9 CSV | Steps 8–9 SVG | Steps 8–9 both | +| --- | ---: | ---: | ---: | ---: | ---: | ---: | +| OxyS/fhlA | 0.0975 | 0.0978 | 0.1710 | 0.1703 | 0.1910 | 0.1947 | +| dense40 | 0.3876 | 0.3913 | 0.8958 | 0.9026 | 0.9090 | 0.9060 | +| sparse400 | 0.1191 | 0.1202 | 0.1616 | 0.1633 | 0.7326 | 0.8260 | +| sparse1000 | 1.5816 | 1.5801 | 1.8284 | 1.8372 | 4.0202 | 4.2578 | + +Existing-mode medians differ by −0.4% to +1.1%, within observed timing variation. +Ensemble outputs and CSVs are byte-identical across base/candidate and modes; +all SVG pair values equal their CSV values. Sparse SVG timings have wider +variation than the compute-heavy cases; these are local workload measurements, +not a general throughput guarantee. Full SVG size scales with all target/query +pairs: the sparse1000 document is 220,685,460 bytes despite having few seed cells. + +Reproduce each process with the following command from the corresponding build; +add either/both output options for the CSV/SVG configurations: + +```sh +OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 OMP_DYNAMIC=FALSE \ + /usr/bin/time -f '%e %U %S %M' taskset -c 0 src/bin/IntaRNA \ + --target=doc/handson/fhlA.fasta --query=doc/handson/OxyS.fasta \ + --model=P --mode=M --seedBP=4 --intLenMax=30 \ + --qIntLoopMax=3 --tIntLoopMax=3 --threads=1 \ + --outNoLP=false --outNoGUend=false --outNumber=0 --outMode=E \ + --default-log-file=/dev/null --out=ensemble.txt +# Optional: --out=bpProb:pairs.csv --out=bpsvg:plot.svg +``` + +For dense40, replace the inputs with `G`×40 and `C`×40 and add `--tAcc=N +--qAcc=N`. For sparse length n (400 or 1000), use target `A`×(n/2) + `G`×8 + +`A`×(n/2−8), query with `C`×8 instead of `G`×8, and the same disabled-accessibility +options. All other flags remain identical. + + +### Audit-fix comparison + +Five interleaved fresh-process runs per configuration compare the retained +`fa2f84d` binary with the rebuilt audit fix. All complete SVG, CSV and ensemble +outputs match byte-for-byte on these native-energy inputs. Time is median +[min, max] seconds; RSS is the maximum over the five runs, in KiB. + +| Input/output | Before time | Fixed time | Before RSS | Fixed RSS | +| --- | ---: | ---: | ---: | ---: | +| OxyS SVG | 0.1938 [0.1932, 0.1968] | 0.1934 [0.1918, 0.1938] | 22,132 | 20,720 | +| OxyS both | 0.1965 [0.1945, 0.1976] | 0.1959 [0.1955, 0.1967] | 22,004 | 20,596 | +| sparse400 SVG | 0.5020 [0.5003, 0.7993] | 0.7491 [0.4806, 0.9672] | 87,876 | 84,548 | +| sparse400 both | 0.5553 [0.5441, 0.8527] | 0.7544 [0.5118, 0.8899] | 87,872 | 84,420 | +| sparse1000 SVG | 4.0260 [3.9535, 5.8240] | 3.8527 [3.8373, 5.5779] | 470,084 | 301,128 | +| sparse1000 both | 4.2407 [4.2148, 6.2242] | 4.0592 [4.0216, 4.1580] | 469,960 | 301,260 | + +Transferring the document buffer reduces 1000×1000 SVG peak RSS by **35.9%**. +There is no general elapsed-time speedup claim: filesystem writes varied widely, +and the sparse400 median increased even though its median process CPU time fell +from 0.44 to 0.42 seconds (SVG) and 0.48 to 0.46 seconds (both). Dense40 controls +remain comparable: ensemble-only 0.3848→0.3902 seconds (+1.4%) and CSV +0.8982→0.8869 seconds (−1.3%), with unchanged output bytes. The output is still +buffered as a full document, so very large matrices remain memory intensive. + +One comparison process received an unexplained SIGTERM after the shorter cases; +partial long-case warmups were discarded and the long case completed in a new +process. The table includes only complete five-run groups. + +To separate serializer cost from file-write stalls, a final sparse400 control +sent SVG to `/dev/null` and ensemble output to STDOUT redirected to `/dev/null`. +Seven interleaved runs gave wall medians 0.412419→0.390979 seconds (ranges +0.411222–0.413401 and 0.388902–0.392267), CPU medians 0.40→0.38 seconds, and peak +RSS 87,744→84,544 KiB. The slowdown did not reproduce without filesystem writes; +the fixed serializer was 5.2% faster in this specific control. Use the same +sparse400 command with `--out=STDOUT --out=bpsvg:/dev/null > /dev/null` to repeat. diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg index fd048c48..aceeb06c 100644 --- a/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg +++ b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.svg @@ -10,6 +10,7 @@ C(i)=(vᵢ,jᵢ) retains one hybrid weight and right boundary for each left pair i. Use d=BW(E_init), t(i,k)=BW(E_interLeft(i,k)); all pair coordinates increase internally. +B(i,j)=BW(energy.getE(i₁,j₁,i₂,j₂,0)) is the complete-boundary Boltzmann factor. Enumerate candidates at i, visiting i in decreasing order @@ -23,15 +24,15 @@ jₖ Lonely pairs allowed: - initial candidate (d, i) + initial candidate (d, i), only if !noGUend or i is not GU continuation (t(i,k) · vₖ, jₖ) for each reachable k ≻ i within loop limits and the complete span limits. Accumulate every considered candidate updateZ(i,j,v,true) -Adds B(i,j)v to Zall and v to the boundary map. -B is the admitted complete-boundary Boltzmann factor. +For admitted complete sites: add B(i,j)v to Zall +and v to the boundary map; otherwise discard the update. Retain only the best continuation at i C(i) ← argmin E_full(i,j,v) @@ -40,7 +41,7 @@ noLP candidate set Let i⁺=i+(1,1) be complementary and s=t(i,i⁺). -Initialize (s·d, i⁺); add direct continuation (s·vᵢ⁺, jᵢ⁺). +Initialize (s·d, i⁺) only if !noGUend or i⁺ is not GU; direct continuation is (s·vᵢ⁺, jᵢ⁺). Other candidates are (s·t(i⁺,k)·vₖ, jₖ), with k ≻ i⁺ and k ≠ i⁺+(1,1). The adjacent continuation is excluded from the loop sum because it was handled directly. The boundary map sums contributions, but the DP cell keeps one right continuation: this is a restricted ensemble. diff --git a/doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg b/doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg index 0a10ebe9..0eabab07 100644 --- a/doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg +++ b/doc/recursions/IntaRNA.PredictorMfeEnsSeedOnly.svg @@ -28,7 +28,8 @@ e H_seed = E_seed(a) + E_init E_full = energy.getE(a₁,e₁,a₂,e₂,H_seed) -Zall += BW(E_full) when needZall is enabled +Zall += BW(E_full) +Only when needZall is enabled and complete-site filters pass. 3 Report and trace the retained seed updateOptima(..., isHybridE=true, incrementZall=true) applies the complete-site filters and ranking. diff --git a/doc/recursions/IntaRNA.SeededPartitionFunction.inside.svg b/doc/recursions/IntaRNA.SeededPartitionFunction.inside.svg index b944efbb..e41542a1 100644 --- a/doc/recursions/IntaRNA.SeededPartitionFunction.inside.svg +++ b/doc/recursions/IntaRNA.SeededPartitionFunction.inside.svg @@ -35,12 +35,12 @@ For Uℓ: let ℓ′=min(C,ℓ+1); D=Uℓ′ unless the ending seed is accepted, then D=S₂. For S₁ or S₂: D=S₂. Once seeded, a later seed never creates a second copy of the path. -2 Aggregate eligible departures, then handle nonstack edges -T₀(q)=∑ H(q,Uℓ) over ℓ≥(noLP ? 2 : 1) -T₁(q)=H(q,S₂) + (noLP ? 0 : H(q,S₁)) -For each nonstack r→q and b∈{0,1}: H(q,D(b,loop)) += t(r,q) · T_b(r) -Nonstack means both coordinates increase, within loop limits, excluding r=q−(1,1). -Reset the suffix to 1; remain seeded if b=1, or accept a singleton seed at q if b=0. +2 Add nonstack arrivals, then aggregate eligible departures +For each nonstack r→q and b∈{0,1}: H(q,D(b,loop)) += t(r,q) · T_b(r) +T₀(q)=∑ H(q,Uℓ) over ℓ≥(noLP ? 2 : 1) +T₁(q)=H(q,S₂) + (noLP ? 0 : H(q,S₁)) +T_b(r) is already computed. Nonstack: strict increase within loop limits, excluding r=q−(1,1). +Reset the suffix to 1; remain seeded if b=1, or accept a singleton seed at q if b=0. 3 Complete seeded paths at every q Z += B(p,q) · T₁(q) diff --git a/doc/recursions/IntaRNA.SeededPartitionFunction.outside.svg b/doc/recursions/IntaRNA.SeededPartitionFunction.outside.svg index c98b9271..fe546876 100644 --- a/doc/recursions/IntaRNA.SeededPartitionFunction.outside.svg +++ b/doc/recursions/IntaRNA.SeededPartitionFunction.outside.svg @@ -45,7 +45,7 @@ 3b Reverse a nonstack edge r→q a = A(q,D(b,loop)) · t(r,q) -Gᵦ(r) += a +G_b(r) += a M(q) += a · T_b(r) Repeat for b=0 and b=1 using the same forward totals. Do not mark the unpaired nucleotides in the loop. diff --git a/src/IntaRNA/BasePairProbabilitySvg.cpp b/src/IntaRNA/BasePairProbabilitySvg.cpp index 1b9e5b91..fda2d1cd 100644 --- a/src/IntaRNA/BasePairProbabilitySvg.cpp +++ b/src/IntaRNA/BasePairProbabilitySvg.cpp @@ -7,6 +7,7 @@ #include #include #include +#include namespace IntaRNA { namespace { @@ -42,10 +43,10 @@ size_t colorClass(Z_type probability) return breaks.size()-2; } -Z_type unpaired(const Accessibility & accessibility,size_t i,const InteractionEnergy & energy) +Z_type unpaired(const Accessibility & accessibility,size_t i,Z_type RT) { const E_type ed=accessibility.getED(i,i); - const Z_type p=E_isINF(ed)?Z_type(0):energy.getBoltzmannWeight(ed); + const Z_type p=E_isINF(ed)?Z_type(0):Z_exp(-E_2_Z(ed)/RT); if (!std::isfinite(p) || p<0 || p>1) throw std::runtime_error("bpsvg: invalid single-nucleotide unpaired probability"); return p; @@ -54,8 +55,11 @@ Z_type unpaired(const Accessibility & accessibility,size_t i,const InteractionEn } // namespace std::string -BasePairProbabilityWriter::svgBlock(const BasePairProbabilities & result,const InteractionEnergy & energy) +BasePairProbabilityWriter::svgBlock(const BasePairProbabilities & result,const InteractionEnergy & energy, + Z_type targetRT,Z_type queryRT) { + if (!std::isfinite(targetRT) || targetRT<=0 || !std::isfinite(queryRT) || queryRT<=0) + throw std::invalid_argument("bpsvg: accessibility energy scales must be positive and finite"); const auto & at=energy.getAccessibility1(); const auto & aq=energy.getAccessibility2().getAccessibilityOrigin(); const auto & t=at.getSequence(); @@ -131,7 +135,7 @@ text { font-family: sans-serif; fill: #243247; } const auto & seq=target?t:q; const auto & acc=target?at:aq; for(size_t i=0;i\n"; - return out.str(); + // Transfer the potentially large document instead of copying its buffer. + return std::move(out).str(); } -void BasePairProbabilityWriter::writeSvg(std::ostream & out,const BasePairProbabilities & result,const InteractionEnergy & energy) +void BasePairProbabilityWriter::writeSvg(std::ostream & out,const BasePairProbabilities & result,const InteractionEnergy & energy, + Z_type targetRT,Z_type queryRT) { - emit(out,svgBlock(result,energy)); + emit(out,svgBlock(result,energy,targetRT,queryRT)); } -void BasePairProbabilityWriter::writeSvgFile(const std::string & filename,const BasePairProbabilities & result,const InteractionEnergy & energy) +void BasePairProbabilityWriter::writeSvgFile(const std::string & filename,const BasePairProbabilities & result,const InteractionEnergy & energy, + Z_type targetRT,Z_type queryRT) { - emitFile(filename,svgBlock(result,energy)); + emitFile(filename,svgBlock(result,energy,targetRT,queryRT)); } } // namespace IntaRNA diff --git a/src/IntaRNA/BasePairProbabilityWriter.cpp b/src/IntaRNA/BasePairProbabilityWriter.cpp index 6d679fdd..d38a744b 100644 --- a/src/IntaRNA/BasePairProbabilityWriter.cpp +++ b/src/IntaRNA/BasePairProbabilityWriter.cpp @@ -3,6 +3,7 @@ #include #include #include +#include #include namespace IntaRNA { @@ -25,7 +26,7 @@ std::string BasePairProbabilityWriter::block(const BasePairProbabilities & resul for(size_t j=0;j::infinity(), + std::numeric_limits::quiet_NaN()}) { + std::ostringstream invalid; + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(invalid,r,energy,rt,1)); + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(invalid,r,energy,1,rt)); + REQUIRE(invalid.str().empty()); + } + BasePairProbabilityWriter::writeSvg(svg,r,energy,1,1); REQUIRE(svg.str().find("target<&")!=std::string::npos); REQUIRE(svg.str().find("Orange outline:")!=std::string::npos); std::ostringstream broken;broken.setstate(std::ios::badbit); - REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(broken,r,energy)); + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(broken,r,energy,1,1)); r.fail();std::ostringstream failed; - REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(failed,r,energy)); + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(failed,r,energy,1,1)); REQUIRE(failed.str().empty()); BasePairProbabilities missing(4,5,true); REQUIRE_THROWS(missing.addRegion({0,1},{1,2},2,mass)); @@ -82,10 +89,10 @@ TEST_CASE("seed annotations commit with regional masses and SVG requires success BasePairProbabilities plain(4,5); REQUIRE_FALSE(plain.collectsSeedPairs());REQUIRE(plain.seedPairs().size1()==0); plain.finalize();std::ostringstream empty; - BasePairProbabilityWriter::writeSvg(empty,plain,energy); + BasePairProbabilityWriter::writeSvg(empty,plain,energy,1,1); REQUIRE(empty.str().find("data-probability=\"NA\"")!=std::string::npos); REQUIRE(empty.str().find("Orange outline:")==std::string::npos); BasePairProbabilities wrongSize(1,1);wrongSize.finalize();std::ostringstream wrong; - REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(wrong,wrongSize,energy)); + REQUIRE_THROWS(BasePairProbabilityWriter::writeSvg(wrong,wrongSize,energy,1,1)); REQUIRE(wrong.str().empty()); } diff --git a/tests/testBasePairProbabilitySvg.py b/tests/testBasePairProbabilitySvg.py index 29047ef3..55e5e576 100644 --- a/tests/testBasePairProbabilitySvg.py +++ b/tests/testBasePairProbabilitySvg.py @@ -18,6 +18,8 @@ def run(t='GGAGGG', q='CCCC', options=(), success=True): keys = {o.split('=')[0] for o in options} + if keys & {'--tAcc', '--qAcc'}: + keys.add('--acc') p = subprocess.run([BIN, f'--target={t}', f'--query={q}', *[o for o in COMMON if o.split('=')[0] not in keys], *options], text=True, capture_output=True) @@ -112,17 +114,48 @@ def plot(t='GGAGGG', q='CCCC', opts=()): named = plot(opts=[f'--tId={name}']) assert name in named.find('s:title', NS).text assert name+" (5' to 3')" in [e.text for e in named.findall('.//s:text', NS)] - # Native accessibility values must use original (not reversed) query order. + # Accessibility uses its own per-strand scale and original query order. pu_t, pu_q = d/'target.pu', d/'query.pu' - native = plot('GGGGAAAACCCC', 'GCGCGAAAACGCGC', - ['--energy=V', '--acc=C', f'--out=tPu:{pu_t}', f'--out=qPu:{pu_q}']) - for strand, path in (('target', pu_t), ('query', pu_q)): - expected_pu = [float(line.split()[1]) for line in path.read_text().splitlines() - if line.strip() and not line.lstrip().startswith('#')] - assert any(0 < v < .9 for v in expected_pu) - for r in native.findall(f'.//s:rect[@data-type="unpaired"][@data-strand="{strand}"]', NS): - assert math.isclose(float(r.attrib['data-probability']), expected_pu[int(r.attrib['data-index'])-1], - rel_tol=2e-5, abs_tol=1e-6) + + def one_nt(path): + return [float(line.split()[1]) for line in path.read_text().splitlines() + if line.strip() and not line.lstrip().startswith('#')] + + def check_frame(root, target_pu, query_pu): + for strand, values in (('target', target_pu), ('query', query_pu)): + assert any(0 < v < .9 for v in values) + for r in root.findall(f'.//s:rect[@data-type="unpaired"][@data-strand="{strand}"]', NS): + assert math.isclose(float(r.attrib['data-probability']), values[int(r.attrib['data-index'])-1], + rel_tol=2e-5, abs_tol=1e-6), (strand, r.attrib, values) + + for temperature in (22, 37): + acc_options = [f'--temperature={temperature}', f'--out=tPu:{pu_t}', f'--out=qPu:{pu_q}'] + native = plot('GGGGAAAACCCC', 'GCGCGAAAACGCGC', ['--energy=V', '--acc=C', *acc_options]) + check_frame(native, one_nt(pu_t), one_nt(pu_q)) + input_t, input_q = d/'input-target.pu', d/'input-query.pu' + input_t.write_bytes(pu_t.read_bytes()); input_q.write_bytes(pu_q.read_bytes()) + for t_mode, q_mode in (('P', 'P'), ('P', 'C'), ('C', 'P'), ('E', 'E'), ('C', 'C')): + ed_t, ed_q = d/'target.ed', d/'query.ed' + options = ['--energy=B', f'--tAcc={t_mode}', f'--qAcc={q_mode}', *acc_options] + if t_mode in ('P', 'E'): + options.append(f'--tAccFile={input_t if t_mode == "P" else ed_t}') + if q_mode in ('P', 'E'): + options.append(f'--qAccFile={input_q if q_mode == "P" else ed_q}') + # The first imported-Pu run also supplies ED fixtures for the E case. + if (t_mode, q_mode) in (('P', 'P'), ('C', 'C')): + options += [f'--out=tAcc:{ed_t}', f'--out=qAcc:{ed_q}'] + converted = plot('GGGGAAAACCCC', 'GCGCGAAAACGCGC', options) + check_frame(converted, one_nt(pu_t), one_nt(pu_q)) + for mode, input_pu, output_pu in ((t_mode, input_t, pu_t), (q_mode, input_q, pu_q)): + if mode in ('P', 'E'): + # Independent input reference; ED quantization loses <2% + # at either temperature, but using interaction RT loses far more. + assert all(math.isclose(actual, original, rel_tol=.02, abs_tol=1e-10) + for actual, original in zip(one_nt(output_pu), one_nt(input_pu))) + if t_mode == q_mode == 'C': + # Nussinov accessibility has RT=1, independent of temperature. + check_frame(converted, [math.exp(-ed) for ed in one_nt(ed_t)], + [math.exp(-ed) for ed in one_nt(ed_q)]) for t,q in (('AAAA','AAAA'), ('G','C')): empty = plot(t,q) assert all(r.attrib['data-probability'] == 'NA' for r in pairs(empty).values()) From b492ed695b0ba66ebcc84988347c4bcdd9bf7edc Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Fri, 9 Oct 2026 22:14:58 +0200 Subject: [PATCH 4/9] feat: revise probability SVG layout and highlight best interaction --- ChangeLog | 13 ++++ README.md | 28 +++++---- src/IntaRNA/BasePairProbabilitySvg.cpp | 84 ++++++++++++++++++------- src/IntaRNA/BasePairProbabilityWriter.h | 11 ++-- src/bin/CommandLineParsing.cpp | 9 +-- src/bin/CommandLineParsing.h | 8 ++- src/bin/IntaRNA.cpp | 8 ++- tests/testBasePairProbabilitySvg.py | 65 ++++++++++++++++--- 8 files changed, 169 insertions(+), 57 deletions(-) diff --git a/ChangeLog b/ChangeLog index 44268e50..19fb5983 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,9 @@ ## Interface and handling +- revise SVG axes, sequence labels and rounded tooltips; show intramolecular + pairing (1-Pu) frames and the best interaction region (PR #259 review) + - BUGFIX : use each strand's accessibility energy scale for SVG frames and unpaired-probability exports, including mixed computed/imported inputs @@ -169,6 +172,16 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261009 Alexander Mitrofanov + * BasePairProbabilitySvg, BasePairProbabilityWriter, CommandLineParsing, IntaRNA : + * place target origin at bottom left and query name below the matrix; include + nondefault sequence names in a subtitle and names in rounded hover values + * show 1-Pu intramolecular pairing probabilities on the nucleotide frame + + outline the best reported interaction region in green, including when + ranked output is disabled; preserve ordinary output counts (PR #259 review) + * README.md, tests/testBasePairProbabilitySvg.py : + * document and cover coordinates, frame probabilities, tooltips and MFE bounds + 261009 Alexander Mitrofanov * BasePairProbabilityWriter, BasePairProbabilitySvg, CommandLineParsing : * pass and validate separate target/query accessibility RT values for SVG diff --git a/README.md b/README.md index d627414e..ca099554 100644 --- a/README.md +++ b/README.md @@ -1876,7 +1876,7 @@ targeted file/stream name: - `qSpotProb:`/`tSpotProb:` [query/target's spot probability profile](#profileSpotProb) (CSV format), respectively - `spotProb:` [all spot probabilities](#spotProb) (CSV format) -- `bpsvg:` [SVG base-pair probability dot plot](#bpsvg), with accessibility frames and seed outlines +- `bpsvg:` [SVG base-pair probability dot plot](#bpsvg), with pairing-probability frames and seed/MFE outlines - `bpProb:` [actual base-pair probabilities](#bpProb) for exact stack-seeded ensembles (CSV matrix) - `qMinE:`/`tMinE:` [the query/target's minimal interaction energy profile](profileMinE) (CSV format), respectively - `pMinE:` [minimal interaction energy for all query-target index pairs](pairMinE) (CSV format) @@ -2082,15 +2082,18 @@ IntaRNA -t GGAGGG -q CCCC --energy=B --acc=N --seedBP=2 \ --model=P --mode=M --out=bpsvg:pairs.svg --out=bpProb:pairs.csv ``` -The query is the horizontal axis and the target the vertical axis, both in -original 5' to 3' order. Sequence names label the axes; nucleotides appear on all -four sides. Their colored backgrounds show **single-nucleotide unpaired -probabilities**, computed from the prediction accessibility model's opening -energies. These are distinct from the interaction-conditional probabilities in -the matrix. Disabled accessibility gives probability one at unconstrained bases. -Each probability rectangle has a hover title with its index (pair), type and -full-precision value. An empty interaction ensemble uses gray `NA` cells while -its accessibility frame remains defined. +The query runs left to right and the target bottom to top, both in original +5' to 3' order, starting at the bottom-left corner. Sequence names label the +axes, with the query name below the matrix. A subtitle names both sequences +when either name differs from the defaults `query` and `target`. Nucleotides +appear on all four sides. Their colored backgrounds show **intramolecular +pairing probability**, `1-Pu`, computed from the accessibility model's opening +energies. These are distinct from the interaction-conditional pair probabilities +in the matrix. Disabled accessibility gives zero at unconstrained bases. +Each probability rectangle has a hover title with sequence names, indices, +probability type and three decimal places (scientific notation with two decimal +places below 0.001). Full-precision values remain in `data-probability` attributes. +An empty interaction ensemble uses gray `NA` cells while its frame remains defined. A shared seven-step, white-to-dark-blue palette encodes both probability types: `[0,0.01)`, `[0.01,0.1)`, `[0.1,0.25)`, `[0.25,0.5)`, `[0.5,0.75)`, @@ -2110,7 +2113,10 @@ contained in searched regions, including overlapping, singleton and mixed-length explicit seeds. This annotation uses the predictor's retained seed family and span limits; it is independent of ranked interaction reporting and can include zero-probability cells. Extension-only pairs have no outline. A seed legend is -included whenever an outline is present. +included whenever an outline is present. A green rectangle encloses the best +reported interaction's site, with a matching legend. The best site is retained +for SVG even with `--outNumber=0`, without adding ranked text output. If no site +passes the reporting filters, the rectangle and its legend are omitted. Filenames support `.gz` compression and the usual `-t#q#` multi-FASTA suffixes. `STDOUT` and `STDERR` emit each complete SVG as one synchronized block. To obtain diff --git a/src/IntaRNA/BasePairProbabilitySvg.cpp b/src/IntaRNA/BasePairProbabilitySvg.cpp index fda2d1cd..d350547f 100644 --- a/src/IntaRNA/BasePairProbabilitySvg.cpp +++ b/src/IntaRNA/BasePairProbabilitySvg.cpp @@ -1,5 +1,6 @@ #include "IntaRNA/BasePairProbabilityWriter.h" #include "IntaRNA/InteractionEnergy.h" +#include "IntaRNA/Interaction.h" #include #include #include @@ -56,7 +57,7 @@ Z_type unpaired(const Accessibility & accessibility,size_t i,Z_type RT) std::string BasePairProbabilityWriter::svgBlock(const BasePairProbabilities & result,const InteractionEnergy & energy, - Z_type targetRT,Z_type queryRT) + Z_type targetRT,Z_type queryRT,const Interaction * mfe) { if (!std::isfinite(targetRT) || targetRT<=0 || !std::isfinite(queryRT) || queryRT<=0) throw std::invalid_argument("bpsvg: accessibility energy scales must be positive and finite"); @@ -72,21 +73,35 @@ BasePairProbabilityWriter::svgBlock(const BasePairProbabilities & result,const I const auto p=empty?Matrix():result.probabilities(); const double w=cell*q.size(), h=cell*t.size(); // Leave space for names even for very short sequences; keep a fixed nt pitch. - const double x=std::max(130.0,(std::max(820.0,9.0*q.getId().size()+80)-w)/2); + const double nameWidth=std::max(9.0*q.getId().size()+80,7.0*(t.getId().size()+q.getId().size())+100); + const double x=std::max(130.0,(std::max(820.0,nameWidth)-w)/2); const double y=std::max(120.0,(9.0*t.getId().size()-h)/2+60); - const double width=2*x+w, height=2*y+h+160; + const double width=2*x+w, height=2*y+h+185; + // SVG y increases downwards, while target coordinates increase upwards. + auto targetY=[&](size_t i) { return y+cell*(t.size()-1-i); }; + const bool markMfe=mfe && !mfe->isEmpty(); + size_t tStart=t.size(),tEnd=0,qStart=q.size(),qEnd=0; + if (markMfe) { + for (const auto & bp:mfe->basePairs) { + if (bp.first>=t.size() || bp.second>=q.size()) + throw std::invalid_argument("bpsvg: MFE interaction coordinates exceed sequence dimensions"); + tStart=std::min(tStart,bp.first); tEnd=std::max(tEnd,bp.first); + qStart=std::min(qStart,bp.second); qEnd=std::max(qEnd,bp.second); + } + } std::ostringstream out; out.imbue(std::locale::classic()); out<::max_digits10); out<<"\n" <<"Base-pair probabilities: "<\n" - <<"Target rows and query columns run 5' to 3'. Pair probabilities are conditional on an allowed interaction. " - <<"The frame shows single-nucleotide unpaired probabilities from the accessibility model.\n" + <<"Target rows run 5' to 3' upwards and query columns run 5' to 3' to the right. Pair probabilities are conditional on an allowed interaction. " + <<"The frame shows single-nucleotide intramolecular pairing probabilities (1 - Pu) from the accessibility model.\n" < text { font-family: sans-serif; fill: #243247; } .heading { font-size: 18px; font-weight: bold; } .axis { font-size: 15px; } +.subtitle { font-size: 13px; } .index, .legend { font-size: 11px; } .nt { font-family: monospace; font-size: 13px; pointer-events: none; } .nt-light { fill: #fff; } @@ -104,6 +119,7 @@ text { font-family: sans-serif; fill: #243247; } .origin { stroke: #888; stroke-width: 1.2; stroke-dasharray: 3 2; } .border { fill: none; stroke: #aaa; stroke-width: 0.7; pointer-events: none; } .seed { fill: none; stroke: #c65d00; stroke-width: 1.8; pointer-events: none; } +.mfe { fill: none; stroke: #16823b; stroke-width: 2.5; pointer-events: none; } )SVG"; out<<"\n"; @@ -111,7 +127,9 @@ text { font-family: sans-serif; fill: #243247; } out<<""<\n"; }; text(width/2,30,"Base-pair probabilities","heading"); - text(x+w/2,y-65,q.getId()+" (5' to 3')","axis"); + if (t.getId()!="target" || q.getId()!="query") + text(width/2,53,t.getId()+" / "+q.getId(),"subtitle"); + text(x+w/2,y+h+65,q.getId()+" (5' to 3')","axis"); out<<""<\n"; auto rect=[&](double rx,double ry,Z_type value,bool missing,const std::string & attributes,const std::string & title) { @@ -120,28 +138,33 @@ text { font-family: sans-serif; fill: #243247; } <<" data-probability=\""; if (missing) out<<"NA"; else out<"<<xml(title)<<": "; - if (missing) out<<"NA (empty interaction ensemble)"; else out<<value; + if (missing) out<<"NA (empty interaction ensemble)"; + else { + if (value<0.001) out<<std::scientific<<std::setprecision(2); + else out<<std::fixed<<std::setprecision(3); + out<<value<<std::defaultfloat<<std::setprecision(std::numeric_limits<Z_type>::max_digits10); + } out<<"\n"; }; out<<"\n"; for(size_t i=0;i\n\n"; for(bool target:{false,true}) { const auto & seq=target?t:q; const auto & acc=target?at:aq; for(size_t i=0;i=.75?"nt nt-light":"nt"); } } @@ -158,12 +181,12 @@ text { font-family: sans-serif; fill: #243247; } const double offset=cell*(i+(idx>0 && !origin?1:0)); if (idx%10==0 || origin) out<<"\n"; + <<"\" x1=\""<<(target?x-cell:x+offset)<<"\" y1=\""<<(target?y+h-offset:y-cell) + <<"\" x2=\""<<(target?x+w+cell:x+offset)<<"\" y2=\""<<(target?y+h-offset:y+h+cell)<<"\"/>\n"; const std::string label=idx==1?"+1":std::to_string(idx); if(target) { - text(x-26,y+cell*(i+.5)+4,label,"index","end"); - text(x+w+26,y+cell*(i+.5)+4,label,"index","start"); + text(x-26,targetY(i)+cell/2+4,label,"index","end"); + text(x+w+26,targetY(i)+cell/2+4,label,"index","start"); } else { text(x+cell*(i+.5),y-28,label,"index"); text(x+cell*(i+.5),y+h+38,label,"index"); @@ -177,12 +200,20 @@ text { font-family: sans-serif; fill: #243247; } for(size_t i=0;i\n"; + <<"\" x=\""<\n"; } out<<"\n"; } - const double legendX=(width-7*86)/2, legendY=y+h+85; - text(width/2,legendY-14,"Probability ranges (base pairs and single-nucleotide accessibility)","legend"); + if (markMfe) { + out<<"\n"; + } + const double legendX=(width-7*86)/2, legendY=y+h+110; + text(width/2,legendY-14,"Probability ranges (intermolecular base pairs and intramolecular pairing)","legend"); out<<"\n"; for(size_t k=0;k\n"; @@ -196,6 +227,11 @@ text { font-family: sans-serif; fill: #243247; } out<<"\n"; text(legendX+28,legendY+61,"Orange outline: pair in an admitted seed within a searched region","legend","start"); } + if (markMfe) { + const double mfeLegendY=legendY+(annotated?76:48); + out<<"\n"; + text(legendX+28,mfeLegendY+13,"Green outline: minimum-free-energy interaction region","legend","start"); + } if(empty) text(width/2,legendY+90,"Gray / NA: empty interaction ensemble; accessibility remains defined","legend"); out<<"\n"; // Transfer the potentially large document instead of copying its buffer. @@ -203,15 +239,15 @@ text { font-family: sans-serif; fill: #243247; } } void BasePairProbabilityWriter::writeSvg(std::ostream & out,const BasePairProbabilities & result,const InteractionEnergy & energy, - Z_type targetRT,Z_type queryRT) + Z_type targetRT,Z_type queryRT,const Interaction * mfe) { - emit(out,svgBlock(result,energy,targetRT,queryRT)); + emit(out,svgBlock(result,energy,targetRT,queryRT,mfe)); } void BasePairProbabilityWriter::writeSvgFile(const std::string & filename,const BasePairProbabilities & result,const InteractionEnergy & energy, - Z_type targetRT,Z_type queryRT) + Z_type targetRT,Z_type queryRT,const Interaction * mfe) { - emitFile(filename,svgBlock(result,energy,targetRT,queryRT)); + emitFile(filename,svgBlock(result,energy,targetRT,queryRT,mfe)); } } // namespace IntaRNA diff --git a/src/IntaRNA/BasePairProbabilityWriter.h b/src/IntaRNA/BasePairProbabilityWriter.h index 56acabf3..64773536 100644 --- a/src/IntaRNA/BasePairProbabilityWriter.h +++ b/src/IntaRNA/BasePairProbabilityWriter.h @@ -7,6 +7,7 @@ namespace IntaRNA { class InteractionEnergy; +class Interaction; /** Explicit writer for a completed actual-pair matrix, independent of trackers. * Whole blocks are synchronized for shared streams. No destructor emits data. */ @@ -32,7 +33,7 @@ class BasePairProbabilityWriter { static void writeFile(const std::string & filename,const BasePairProbabilities & result, const RnaSequence & target,const RnaSequence & query,const std::string & separator=";"); /** Write a complete SVG dot plot with original target rows/query columns. - * Includes 1-nt accessibility, signed coordinate guides and a CSS color legend. + * Includes intramolecular pairing, signed coordinate guides and a CSS color legend. * Seed outlines are included if the result collected admitted seed pairs. * @param out destination stream; an I/O failure throws * @param result explicitly finalized raw result @@ -40,9 +41,10 @@ class BasePairProbabilityWriter { * @param targetRT positive finite energy scale used to convert target ED to Pu * @param queryRT positive finite energy scale used to convert query ED to Pu; * accessibility scales can differ from the interaction model's RT + * @param mfe optional best interaction in original coordinates for a region outline */ static void writeSvg(std::ostream & out,const BasePairProbabilities & result, - const InteractionEnergy & energy,Z_type targetRT,Z_type queryRT); + const InteractionEnergy & energy,Z_type targetRT,Z_type queryRT,const Interaction * mfe = nullptr); /** Write a complete SVG to a file or STDOUT/STDERR, including gzip support. * Validates the whole document before opening the destination. * @param filename destination (already expanded for multi-FASTA) @@ -50,12 +52,13 @@ class BasePairProbabilityWriter { * @param energy original, unshifted energy/accessibility model * @param targetRT positive finite target accessibility energy scale * @param queryRT positive finite query accessibility energy scale + * @param mfe optional best interaction in original coordinates for a region outline */ static void writeSvgFile(const std::string & filename,const BasePairProbabilities & result, - const InteractionEnergy & energy,Z_type targetRT,Z_type queryRT); + const InteractionEnergy & energy,Z_type targetRT,Z_type queryRT,const Interaction * mfe = nullptr); private: static std::string svgBlock(const BasePairProbabilities & result,const InteractionEnergy & energy, - Z_type targetRT,Z_type queryRT); + Z_type targetRT,Z_type queryRT,const Interaction * mfe); static void emitFile(const std::string & name,const std::string & data); static std::string block(const BasePairProbabilities & result,const RnaSequence & target, const RnaSequence & query,const std::string & separator); diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 756d8b0f..7f6f1423 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -2234,7 +2234,7 @@ getEnergyHandler( const Accessibility& accTarget, const ReverseAccessibility& ac OutputConstraint CommandLineParsing:: -getOutputConstraint( const InteractionEnergy & energy ) const +getOutputConstraint( const InteractionEnergy & energy, bool forPrediction ) const { checkIfParsed(); if (!rri.empty()) { @@ -2250,7 +2250,7 @@ getOutputConstraint( const InteractionEnergy & energy ) const default : throw std::runtime_error("CommandLineParsing::getOutputConstraint() : unsupported outOverlap value "+toString(outOverlap.val)); } return OutputConstraint( - outNumber.val + (forPrediction && !outPrefix2streamName.at(OP_bpsvg).empty() ? std::max(size_t(1),size_t(outNumber.val)) : outNumber.val) , overlap , Ekcal_2_E(outMaxE.val) , Ekcal_2_E(outDeltaE.val) @@ -2501,14 +2501,15 @@ BasePairProbabilities * CommandLineParsing::getBasePairProbabilityResult(const I throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); return new BasePairProbabilities(energy.size1(),energy.size2(),!outPrefix2streamName.at(OP_bpsvg).empty()); } -void CommandLineParsing::writeBasePairProbabilities(const BasePairProbabilities & result,const InteractionEnergy & energy) const { +void CommandLineParsing::writeBasePairProbabilities(const BasePairProbabilities & result,const InteractionEnergy & energy, + const Interaction * mfe) const { const auto & t=energy.getAccessibility1().getSequence(); const auto & q=energy.getAccessibility2().getAccessibilityOrigin().getSequence(); if (!outPrefix2streamName.at(OP_bpProb).empty()) BasePairProbabilityWriter::writeFile(getFullFilename(outPrefix2streamName.at(OP_bpProb),&t,&q),result,t,q,outSep); if (!outPrefix2streamName.at(OP_bpsvg).empty()) BasePairProbabilityWriter::writeSvgFile(getFullFilename(outPrefix2streamName.at(OP_bpsvg),&t,&q),result,energy, - getAccessibilityRT(tAcc.val),getAccessibilityRT(qAcc.val)); + getAccessibilityRT(tAcc.val),getAccessibilityRT(qAcc.val),mfe); } Predictor* diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h index 14e010f5..f0e796ad 100644 --- a/src/bin/CommandLineParsing.h +++ b/src/bin/CommandLineParsing.h @@ -282,7 +282,7 @@ class CommandLineParsing { * parameters * @param energy the interaction energy handler to be used * @param output the output handler to be used - * @param pairProbabilities optional non-owning raw result sink for exact stack seeds + * @param pairProbabilities optional non-owning raw result sink for supported ensemble predictors * @return the newly allocated Predictor object to be deleted by the calling * function */ @@ -296,7 +296,8 @@ class CommandLineParsing { */ BasePairProbabilities * getBasePairProbabilityResult(const InteractionEnergy & energy) const; /** Publish a successful sequence-pair result using ordinary filename rules. */ - void writeBasePairProbabilities(const BasePairProbabilities & result,const InteractionEnergy & energy) const; + void writeBasePairProbabilities(const BasePairProbabilities & result,const InteractionEnergy & energy, + const Interaction * mfe = nullptr) const; /** * Provides the seed constraint according to the user settings @@ -331,9 +332,10 @@ class CommandLineParsing { /** * The constraints to be applied to the interaction output generation * @param energy the interaction energy to be used for computation + * @param forPrediction retain a best site for SVG even with zero ranked output * @return the output constraints to be applied */ - OutputConstraint getOutputConstraint( const InteractionEnergy & energy ) const; + OutputConstraint getOutputConstraint( const InteractionEnergy & energy, bool forPrediction = false ) const; /** * The stream to write the interaction output to diff --git a/src/bin/IntaRNA.cpp b/src/bin/IntaRNA.cpp index bb971e29..e67b4876 100644 --- a/src/bin/IntaRNA.cpp +++ b/src/bin/IntaRNA.cpp @@ -239,9 +239,9 @@ int main(int argc, char **argv){ // setup collecting output handler to ensure // k-best output per query-target combination // and not per region combination if not requested - OutputHandlerInteractionList bestInteractions( parameters.getOutputConstraint(*energy), + OutputHandlerInteractionList bestInteractions( parameters.getOutputConstraint(*energy,true), (parameters.reportBestPerRegion() ? std::numeric_limits::max() : 1 ) - * parameters.getOutputConstraint(*energy).reportMax ); + * parameters.getOutputConstraint(*energy,true).reportMax ); // One owner spans all disjoint searched regions for this sequence pair. std::unique_ptr pairProbabilities(parameters.getBasePairProbabilityResult(*energy)); @@ -356,12 +356,14 @@ int main(int argc, char **argv){ : (*bestInteractions.begin())->energy + bestInteractions.getOutputConstraint().deltaE; // The collector is sorted by energy, so all later entries are worse. for( const Interaction * inter : bestInteractions) { + if (output->getOutputConstraint().reportMax==0) break; if (inter->energy > maxMergedE) { break; } output->add(*inter); } } - if (pairProbabilities) parameters.writeBasePairProbabilities(*pairProbabilities,*energy); + if (pairProbabilities) parameters.writeBasePairProbabilities(*pairProbabilities,*energy, + bestInteractions.empty()?nullptr:*bestInteractions.begin()); } // successful sequence-pair computation #if INTARNA_MULITHREADING diff --git a/tests/testBasePairProbabilitySvg.py b/tests/testBasePairProbabilitySvg.py index 55e5e576..597afd6a 100644 --- a/tests/testBasePairProbabilitySvg.py +++ b/tests/testBasePairProbabilitySvg.py @@ -55,6 +55,15 @@ def plot(t='GGAGGG', q='CCCC', opts=()): assert len(cells) == 24 assert '+1' in [e.text for e in root.findall('.//s:text[@class="index"]', NS)] assert not root.findall('.//s:line[@data-index="1"]', NS) + assert root.find('.//s:text[@class="subtitle"]', NS) is None + query_axis = root.find('.//s:text[@class="axis"][@x]', NS) + assert query_axis.text == "query (5' to 3')" + assert float(query_axis.attrib['y']) > max(float(r.attrib['y']) for r in cells.values())+20 + # Original target coordinate 1 is at the bottom; query coordinate 1 at left. + for i in range(1, 7): + for j in range(1, 5): + assert float(cells[i, j].attrib['y']) == float(cells[1, 1].attrib['y'])-20*(i-1) + assert float(cells[i, j].attrib['x']) == float(cells[1, 1].attrib['x'])+20*(j-1) for i, row in enumerate(rows[1:], 1): for j, value in enumerate(row[1:], 1): r = cells[i, j] @@ -62,22 +71,57 @@ def plot(t='GGAGGG', q='CCCC', opts=()): bounds = [0.01, 0.1, 0.25, 0.5, 0.75, 0.9] expected_color = sum(float(value) >= boundary for boundary in bounds) assert f'p{expected_color}' in r.attrib['class'].split() - assert 'Base-pair probability' in r.find('s:title', NS).text + title = r.find('s:title', NS).text + assert title.startswith(f'Base-pair probability (target {i}, query {j}): ') + formatted = f'{float(value):.2e}' if float(value) < .001 else f'{float(value):.3f}' + assert title.endswith(': '+formatted), title assert all(float(cells[3, j].attrib['data-probability']) == 0 for j in range(1, 5)) - frame = root.findall('.//s:rect[@data-type="unpaired"]', NS) + frame = root.findall('.//s:rect[@data-type="intramolecular-pairing"]', NS) assert len(frame) == 2*(6+4) - assert all(float(r.attrib['data-probability']) == 1 for r in frame) - assert all('Unpaired probability' in r.find('s:title', NS).text for r in frame) + assert all(float(r.attrib['data-probability']) == 0 for r in frame) + assert all('Intramolecular pairing probability' in r.find('s:title', NS).text for r in frame) css = root.find('s:style', NS).text assert all(f'.p{i} {{ fill:' in css for i in range(7)) assert len(root.findall('.//s:g[@class="probability-legend"]/s:rect', NS)) == 7 # All admitted computed two-pair stacks, in original query orientation. expected = {(i+k+1, 4-j-k) for i in (0, 3, 4) for j in range(3) for k in range(2)} assert seeds(root) == expected, seeds(root) + for r in root.findall('.//s:rect[@class="seed"][@data-target]', NS): + bp = cells[int(r.attrib['data-target']), int(r.attrib['data-query'])] + assert float(r.attrib['x']) == float(bp.attrib['x'])+1 + assert float(r.attrib['y']) == float(bp.attrib['y'])+1 + for r in frame: + idx = int(r.attrib['data-index']) + axis, bp = ('y', cells[idx, 1]) if r.attrib['data-strand'] == 'target' else ('x', cells[1, idx]) + assert float(r.attrib[axis]) == float(bp.attrib[axis]) # Requesting SVG must not affect CSV or ensemble/ranked output. baseline = run(options=['--outNumber=1', '--outMode=C', '--outCsvCols=Zall,Eall']) plotted = run(options=['--outNumber=1', '--outMode=C', '--outCsvCols=Zall,Eall', f'--out=bpsvg:{svg}']) assert baseline.stdout == plotted.stdout + # The green rectangle covers the reported MFE region, including when + # ranked output is disabled. Query coordinates retain their original order. + for target, query, options in ( + ('GGAGGG', 'CCCC', []), + ('AAGGGGAA', 'AACCCCAA', []), + ('GGAGGG', 'CCCC', ['--tIdxPos0=-4', '--qIdxPos0=-2']), + ): + reference = run(target, query, [*options, '--outNumber=1', '--outMode=C', + '--outCsvCols=start1,end1,start2,end2']) + row = next(csv.DictReader(reference.stdout.splitlines(), delimiter=';')) + for count in (0, 1): + mfe_plot = plot(target, query, [*options, f'--outNumber={count}']) + outline = mfe_plot.find('.//s:rect[@class="mfe"][@data-target-start]', NS) + assert outline is not None + bounds = tuple(int(outline.attrib[k]) for k in ( + 'data-target-start', 'data-target-end', 'data-query-start', 'data-query-end')) + assert bounds == tuple(int(row[k]) for k in ('start1', 'end1', 'start2', 'end2')) + ts, te, qs, qe = bounds + matrix = pairs(mfe_plot) + assert float(outline.attrib['x']) == float(matrix[ts, qs].attrib['x']) + assert float(outline.attrib['y']) == float(matrix[te, qs].attrib['y']) + assert float(outline.attrib['width']) == float(matrix[ts, qe].attrib['x'])-float(matrix[ts, qs].attrib['x'])+20 + assert float(outline.attrib['height']) == float(matrix[ts, qs].attrib['y'])-float(matrix[te, qs].attrib['y'])+20 + assert 'Green outline:' in ''.join(mfe_plot.itertext()) # Explicit singleton/mixed seed masks exclude extension-only cells. opts = ['--seedTQ=1|&4|,2||&2||'] mixed = plot('GGGG', 'CCCC', opts) @@ -103,7 +147,7 @@ def plot(t='GGAGGG', q='CCCC', opts=()): idx = int(g.attrib['data-index']) r = signed_cells[idx,-15] if strand == 'target' else signed_cells[-15,idx] axis = 'y' if strand == 'target' else 'x' - expected_pos = float(r.attrib[axis]) + (20 if idx > 1 else 0) + expected_pos = float(r.attrib[axis]) + (0 if idx > 1 else 20) if strand == 'target' else float(r.attrib[axis]) + (20 if idx > 1 else 0) assert float(g.attrib[axis+'1']) == expected_pos assert ('major' in g.attrib['class']) == (idx % 50 == 0) assert ('origin' in g.attrib['class']) == (idx == 1) @@ -113,6 +157,10 @@ def plot(t='GGAGGG', q='CCCC', opts=()): name = '' named = plot(opts=[f'--tId={name}']) assert name in named.find('s:title', NS).text + assert named.find('.//s:text[@class="subtitle"]', NS).text == name+' / query' + assert all(name in r.find('s:title', NS).text for r in pairs(named).values()) + assert all(name in r.find('s:title', NS).text for r in named.findall( + './/s:rect[@data-type="intramolecular-pairing"][@data-strand="target"]', NS)) assert name+" (5' to 3')" in [e.text for e in named.findall('.//s:text', NS)] # Accessibility uses its own per-strand scale and original query order. pu_t, pu_q = d/'target.pu', d/'query.pu' @@ -124,9 +172,9 @@ def one_nt(path): def check_frame(root, target_pu, query_pu): for strand, values in (('target', target_pu), ('query', query_pu)): assert any(0 < v < .9 for v in values) - for r in root.findall(f'.//s:rect[@data-type="unpaired"][@data-strand="{strand}"]', NS): - assert math.isclose(float(r.attrib['data-probability']), values[int(r.attrib['data-index'])-1], - rel_tol=2e-5, abs_tol=1e-6), (strand, r.attrib, values) + for r in root.findall(f'.//s:rect[@data-type="intramolecular-pairing"][@data-strand="{strand}"]', NS): + assert math.isclose(float(r.attrib['data-probability']), 1-values[int(r.attrib['data-index'])-1], + rel_tol=2e-5, abs_tol=5e-6), (strand, r.attrib, values) for temperature in (22, 37): acc_options = [f'--temperature={temperature}', f'--out=tPu:{pu_t}', f'--out=qPu:{pu_q}'] @@ -160,6 +208,7 @@ def check_frame(root, target_pu, query_pu): empty = plot(t,q) assert all(r.attrib['data-probability'] == 'NA' for r in pairs(empty).values()) assert not seeds(empty) + assert not empty.findall('.//s:rect[@class="mfe"]', NS) assert 'empty interaction ensemble' in ''.join(empty.itertext()) compressed = d/'pairs.svg.gz' run(options=[f'--out=bpsvg:{compressed}']) From 322031fc124eec6d01667fbf1f8bedf0ff35f038 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Fri, 9 Oct 2026 22:16:07 +0200 Subject: [PATCH 5/9] feat: add exact unseeded base-pair probabilities (step 10) --- ChangeLog | 12 ++ README.md | 17 +- src/IntaRNA/PredictorMfeEns2d.cpp | 122 ++++++++++- src/IntaRNA/PredictorMfeEns2d.h | 24 ++- src/bin/CommandLineParsing.cpp | 30 +-- tests/Makefile.am | 1 + tests/PredictorUnseededProbability_test.cpp | 217 ++++++++++++++++++++ tests/testBasePairProbabilities.py | 15 +- tests/testBasePairProbabilitySvg.py | 7 +- 9 files changed, 418 insertions(+), 27 deletions(-) create mode 100644 tests/PredictorUnseededProbability_test.cpp diff --git a/ChangeLog b/ChangeLog index 19fb5983..e818a899 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,9 @@ ## Interface and handling +- enable base-pair CSV/SVG probabilities for exact unseeded model P/mode M + predictions (PR #259, step 10) + - revise SVG axes, sequence labels and rounded tooltips; show intramolecular pairing (1-Pu) frames and the best interaction region (PR #259 review) @@ -172,6 +175,15 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261009 Alexander Mitrofanov + * PredictorMfeEns2d, CommandLineParsing : + + compute actual-pair masses by reversing the existing unseeded recurrence; + preserve site filters, noLP and numerical cutoffs, and merge disjoint regions + * enable CSV/SVG with --noSeed; fix default-ended region offset truncation + * tests/PredictorUnseededProbability_test.cpp, tests/Makefile.am, README.md : + + compare against independent chain enumeration, filters, regions and failures; + document the public optional sink (PR #259 step 10) + 261009 Alexander Mitrofanov * BasePairProbabilitySvg, BasePairProbabilityWriter, CommandLineParsing, IntaRNA : * place target origin at bottom left and query name below the matrix; include diff --git a/README.md b/README.md index ca099554..888f143f 100644 --- a/README.md +++ b/README.md @@ -2016,8 +2016,8 @@ or `STDERR` instead of a file name. ## Actual base-pair probabilities `--out=bpProb:FILE` writes the probability that each target/query nucleotide pair -actually pairs, conditional on an allowed seeded interaction. It requires -`--model=P --mode=M` and a stack-only seed handler. For example: +actually pairs, conditional on an allowed interaction. It requires +`--model=P --mode=M` with `--noSeed` or a stack-only seed handler. For example: ```sh IntaRNA -t GGAGGG -q CCCC --energy=B --acc=N --seedBP=2 \ @@ -2031,7 +2031,7 @@ Computed seeds preserve their handler-specific admission thresholds. Explicit stacked seeds may have different lengths, including singletons; competing explicit patterns at the same start still retain only the selected minimum-energy pattern. Other predictors retain their existing behavior. Heuristic, kinetic, -evaluation, seed-only, unseeded and bulged-seed `bpProb` requests are rejected. +evaluation, seed-only and bulged-seed `bpProb` requests are rejected. Rows follow the original target and columns the original query, both 5' to 3'. The header starts with `bpProb`; labels contain each nucleotide and its display @@ -2054,14 +2054,19 @@ bulge can have positive site coverage and zero actual-pair probability. There is no added unbound-state weight. Individual pair events can coexist, so `1 - sum(pair probabilities)` is not a probability of no interaction. -The backend checks nonnegative arithmetic and reports numerical range failures +The unseeded exact predictor is enabled by `--noSeed` with the same CSV/SVG +outputs. Its outside pass follows the existing fixed-right-boundary recurrence, +including its numerical zero cutoffs, without changing ordinary prediction. + +The seeded backend checks nonnegative arithmetic and reports numerical range failures as errors, including underflow, overflow and NaN. Positive subnormals are rejected; there is currently no scaling fallback. Failed or cancelled regional work never produces a successful-looking probability matrix. Ordinary floating-point rounding remains; "exact" describes the enumerated ensemble, not real arithmetic. For library use, keep a `BasePairProbabilities` owner alive across predictions, -pass its optional non-owning pointer to `PredictorMfeEns2dSeedExtension`, and call +pass its optional non-owning pointer to `PredictorMfeEns2dSeedExtension` (stack +seeds) or `PredictorMfeEns2d` (unseeded), and call `finalize()` only after all requested disjoint regions succeed. A writer is called explicitly after finalization; its destructor publishes nothing. The output constraint must enable `needZall`. The [standalone API example](doc/seeded-probabilities-example.cpp) @@ -2074,7 +2079,7 @@ show the public API and measured limits. `--out=bpsvg:FILE` renders the same actual-pair probabilities as `bpProb` in a standalone SVG. Both outputs can be requested together. The same exact -`--model=P --mode=M` stack-seed restrictions, regional ensemble, normalization +`--model=P --mode=M` unseeded/stack-seed restrictions, regional ensemble, normalization and success-only publication rules apply. ```sh diff --git a/src/IntaRNA/PredictorMfeEns2d.cpp b/src/IntaRNA/PredictorMfeEns2d.cpp index 9b9bb28e..490bfe9b 100644 --- a/src/IntaRNA/PredictorMfeEns2d.cpp +++ b/src/IntaRNA/PredictorMfeEns2d.cpp @@ -1,5 +1,6 @@ #include "IntaRNA/PredictorMfeEns2d.h" +#include "IntaRNA/PartitionArithmetic.h" #include @@ -11,10 +12,14 @@ PredictorMfeEns2d:: PredictorMfeEns2d( const InteractionEnergy & energy , OutputHandler & output - , PredictionTracker * predTracker ) + , PredictionTracker * predTracker + , BasePairProbabilities * pairProbabilities ) : PredictorMfeEns(energy,output,predTracker) + , pairProbabilities(pairProbabilities) , hybridZ( 0,0 ) { + if (pairProbabilities && !output.getOutputConstraint().needZall) + throw std::invalid_argument("unseeded pair probabilities require needZall"); } @@ -33,6 +38,13 @@ void PredictorMfeEns2d:: predict( const IndexRange & r1 , const IndexRange & r2 ) +{ + try { predictRegion(r1,r2); } + catch (...) { if (pairProbabilities) pairProbabilities->fail(); throw; } +} + +void +PredictorMfeEns2d::predictRegion(const IndexRange & r1,const IndexRange & r2) { #if INTARNA_MULITHREADING #pragma omp critical(intarna_omp_logOutput) @@ -47,20 +59,34 @@ predict( const IndexRange & r1 throw std::runtime_error("PredictorMfeEns2d::predict("+toString(r1)+","+toString(r2)+") is not sane"); #endif + if (pairProbabilities && (!r1.isAscending() || !r2.isAscending() + || r1.from>=energy.getAccessibility1().getSequence().size() + || r2.from>=energy.getAccessibility2().getSequence().size())) + throw std::invalid_argument("unseeded pair probabilities: invalid region"); + if (pairProbabilities && pairProbabilities->status()!=BasePairProbabilities::Status::pending) + throw std::logic_error("unseeded pair probabilities: accumulator is not pending"); + // set index offset energy.setOffset1(r1.from); energy.setOffset2(r2.from); // resize matrix hybridZ.resize( std::min( energy.size1() - , (r1.to==RnaSequence::lastPos?energy.size1()-1:r1.to)-r1.from+1 ) + , r1.to==RnaSequence::lastPos?energy.size1():r1.to-r1.from+1 ) , std::min( energy.size2() - , (r2.to==RnaSequence::lastPos?energy.size2()-1:r2.to)-r2.from+1 ) ); + , r2.to==RnaSequence::lastPos?energy.size2():r2.to-r2.from+1 ) ); // initialize mfe interaction for updates initOptima(); // initialize overall partition function for updates initZ(); + // No probability buffers or reverse work are needed for ordinary prediction. + Z2dMatrix outside, masses; + Z_type denominator=0; + if (pairProbabilities) { + outside.resize(hybridZ.size1(),hybridZ.size2()); + masses.resize(hybridZ.size1(),hybridZ.size2()); + } // for all right ends j1 for (size_t j1 = hybridZ.size1(); j1-- > 0; ) { @@ -78,12 +104,26 @@ predict( const IndexRange & r1 // fill matrix and store best interaction fillHybridZ( j1, j2, 0, 0, true ); + if (pairProbabilities) + accumulatePairMasses(j1,j2,outside,masses,denominator); } } + if (pairProbabilities && denominator!=PartitionArithmetic::check(Zall)) + throw std::logic_error("unseeded pair probabilities: updateZ override changed the partition objective"); // report mfe interaction reportOptima(); + if (pairProbabilities) { + const size_t n=hybridZ.size1(), m=hybridZ.size2(); + const auto first=energy.getBasePair(0,0), last=energy.getBasePair(n-1,m-1); + // An unseeded prediction has no seed annotations, even for an API sink + // that requested a mask for use by a shared writer. + const bool annotate=pairProbabilities->collectsSeedPairs(); + Matrix seeds(annotate?n:0,annotate?m:0,0); + pairProbabilities->addRegion(IndexRange(first.first,last.first), + IndexRange(last.second,first.second),denominator,masses,annotate?&seeds:nullptr); + } } //////////////////////////////////////////////////////////////////////////// @@ -201,5 +241,81 @@ fillHybridZ( const size_t j1, const size_t j2 //////////////////////////////////////////////////////////////////////////// +void +PredictorMfeEns2d::accumulatePairMasses(size_t j1,size_t j2,Z2dMatrix & outside, + Z2dMatrix & masses,Z_type & denominator) const +{ + using A=PartitionArithmetic; + const bool noLP=output.getOutputConstraint().noLP; + const size_t shift=noLP?1:0; + // Use exactly the rectangle and epsilon gates of fillHybridZ/updateZ: this + // reverses the represented ensemble, including its existing zero cutoff. + const size_t begin1=j1-std::min(j1,energy.getAccessibility1().getMaxLength()+1); + const size_t begin2=j2-std::min(j2,energy.getAccessibility2().getMaxLength()+1); + for(size_t i1=j1+1;i1-- > begin1;) for(size_t i2=j2+1;i2-- > begin2;) { + outside(i1,i2)=0; + const Z_type h=A::check(hybridZ(i1,i2)); + if (!Z_equal(h,0) && isValidOutputSite(i1,j1,i2,j2)) { + const Z_type b=A::check(energy.getBoltzmannWeight(energy.getE(i1,j1,i2,j2,E_type(0)))); + outside(i1,i2)=b; + denominator=A::add(denominator,A::multiply(h,b)); + } + } + + const auto first=energy.getBasePair(0,0); + const auto last=energy.getBasePair(masses.size1()-1,masses.size2()-1); + auto addMass=[&](size_t i1,size_t i2,Z_type value) { + const auto bp=energy.getBasePair(i1,i2); + auto & mass=masses(bp.first-first.first,bp.second-last.second); + mass=A::add(mass,value); + }; + const Z_type initiation=A::check(energy.getBoltzmannWeight(energy.getE_init())); + // Every suffix state contains its left pair once. Reverse topological + // traversal gathers all enclosing left boundaries before visiting a suffix. + for(size_t i1=begin1;i1<=j1;++i1) for(size_t i2=begin2;i2<=j2;++i2) { + const Z_type a=outside(i1,i2), h=hybridZ(i1,i2); + if (a==0 || h==0) continue; + addMass(i1,i2,A::multiply(a,h)); + const size_t w1=j1-i1+1, w2=j2-i2+1; + if (w1<=1 || w2<=1) continue; + + Z_type stack=1; + if (!noLP) { + const Z_type edge=A::multiply(a, + energy.getBoltzmannWeight(energy.getE_interLeft(i1,j1,i2,j2))); + outside(j1,j2)=A::add(outside(j1,j2),edge); + } else { + if (!energy.areComplementary(i1+1,i2+1)) continue; + stack=A::check(energy.getBoltzmannWeight(energy.getE_interLeft(i1,i1+1,i2,i2+1))); + const Z_type edge=A::multiply(a,stack); + if (w1==2 && w2==2) { + // The noLP terminal has no singleton H state: its second pair + // is implicit in stack * initiation and must be counted here. + addMass(j1,j2,A::multiply(edge,initiation)); + } else { + outside(i1+1,i2+1)=A::add(outside(i1+1,i2+1),edge); + } + } + if (w1<=2 || w2<=2) continue; + for(size_t k1=std::min(j1-1,i1+energy.getMaxInternalLoopSize1()+1+shift);k1>i1+shift;--k1) { + for(size_t k2=std::min(j2-1,i2+energy.getMaxInternalLoopSize2()+1+shift);k2>i2+shift;--k2) { + if (noLP && k1-1==i1+shift && k2-1==i2+shift) continue; + if (Z_equal(hybridZ(k1,k2),0)) continue; + const Z_type loop=A::check(energy.getBoltzmannWeight( + energy.getE_interLeft(i1+shift,k1,i2+shift,k2))); + const Z_type edge=A::multiply(a,A::multiply(stack,loop)); + outside(k1,k2)=A::add(outside(k1,k2),edge); + if (noLP) { + // This decomposition consumes (i,i+) before jumping over + // the loop to k; only i and k have explicit suffix states. + addMass(i1+1,i2+1,A::multiply(edge,hybridZ(k1,k2))); + } + } + } + } +} + +//////////////////////////////////////////////////////////////////////////// + } // namespace diff --git a/src/IntaRNA/PredictorMfeEns2d.h b/src/IntaRNA/PredictorMfeEns2d.h index 57b7f5c3..169b9881 100644 --- a/src/IntaRNA/PredictorMfeEns2d.h +++ b/src/IntaRNA/PredictorMfeEns2d.h @@ -3,6 +3,7 @@ #define INTARNA_PREDICTORMFEENS2D_H_ #include "IntaRNA/PredictorMfeEns.h" +#include "IntaRNA/BasePairProbabilities.h" #include "IntaRNA/Interaction.h" #include "IntaRNA/Matrix.h" @@ -34,10 +35,13 @@ class PredictorMfeEns2d: public PredictorMfeEns { * @param predTracker the prediction tracker to be used or NULL if no * tracking is to be done; if non-NULL, the tracker gets deleted * on this->destruction. + * @param pairProbabilities optional non-owning sink; requires needZall and + * must remain pending until all disjoint regions have succeeded */ PredictorMfeEns2d( const InteractionEnergy & energy , OutputHandler & output - , PredictionTracker * predTracker ); + , PredictionTracker * predTracker + , BasePairProbabilities * pairProbabilities = nullptr ); virtual ~PredictorMfeEns2d(); @@ -58,6 +62,9 @@ class PredictorMfeEns2d: public PredictorMfeEns { protected: + //! optional sequence-pair owner, never deleted or finalized by this predictor + BasePairProbabilities * pairProbabilities; + //! access to the interaction energy handler of the super class using PredictorMfeEns::energy; @@ -70,6 +77,21 @@ class PredictorMfeEns2d: public PredictorMfeEns { protected: + /** Compute one region, allowing predict() to fail the sink on any exception. */ + void predictRegion(const IndexRange & r1,const IndexRange & r2); + + /** Reverse the current fixed-right hybrid recurrence. Adjoint seeds obey the + * same site filters and numerical-zero tests as updateZ(). Implicit noLP + * stacking partners are counted in addition to the explicit DP-state pairs. + * @param j1 right boundary in local target coordinates + * @param j2 right boundary in local reversed-query coordinates + * @param outside scratch adjoints, same dimensions as hybridZ + * @param masses accumulated numerators in original regional orientation + * @param denominator accumulated boundary contributions, checked against Zall + */ + void accumulatePairMasses(size_t j1,size_t j2,Z2dMatrix & outside, + Z2dMatrix & masses,Z_type & denominator) const; + /** * Computes all entries of the hybridE matrix for interactions ending in * p=j1 and q=j2 and reports every complete boundary via updateCompleteZ(), diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 7f6f1423..d5a0f30e 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -958,7 +958,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) "\n 'tPu:' (target) unpaired probability values (RNAplfold format OR '.agz' file ending for IntaRNA's binary format)." "\n 'pMinE:' (target+query) for each index pair the minimal energy of any interaction covering the pair (CSV format)" "\n 'bpsvg:' (target+query) SVG actual-pair probability plot with accessibility frame and seed annotations; same requirements as bpProb" - "\n 'bpProb:' (target+query) actual-pair conditional probabilities for exact --model=P --mode=M stack-seed ensembles (CSV matrix, original strand order; NA for an empty ensemble)" + "\n 'bpProb:' (target+query) actual-pair conditional probabilities for --model=P --mode=M with --noSeed or stack-only seeds (CSV matrix, original strand order; NA for an empty ensemble)" "\n 'spotProb:' (target+query) tracks for a given set of interaction spots their probability to be covered by an interaction. If no spots are provided, probabilities for all index combinations are computed. Spots are encoded by comma-separated 'idxT&idxQ' pairs (target-query). For each spot a probability is provided in concert with the probability that none of the spots (encoded by '0&0') is covered (CSV format). The spot encoding is followed colon-separated by the output stream/file name, eg. '--out=\"spotProb:3&76,59&2:STDERR\"'. NOTE: value has to be quoted due to '&' symbol!" "\nFor each, provide a file name or STDOUT/STDERR to write to the respective output stream." ).c_str()) @@ -1246,8 +1246,8 @@ parse(int argc, char** argv) outSep = unescaped_string::getUnescaped( outSep ); if (hasBasePairProbabilityOutput()) { - if (model.val!='P' || mode.val!='M' || noSeedRequired || !rri.empty() || windowWidth.val!=0) - throw error("bpProb/bpsvg requires --model=P --mode=M with stack-only seeds and no windows"); + if (model.val!='P' || mode.val!='M' || !rri.empty() || windowWidth.val!=0) + throw error("bpProb/bpsvg requires --model=P --mode=M with --noSeed or stack-only seeds and no windows"); outNeedsZall=true; } // K needs a seed even before the usual --noSeed model normalization. @@ -2496,10 +2496,12 @@ bool CommandLineParsing::hasBasePairProbabilityOutput() const { } BasePairProbabilities * CommandLineParsing::getBasePairProbabilityResult(const InteractionEnergy & energy) const { if (!hasBasePairProbabilityOutput()) return nullptr; - std::unique_ptr handler(getSeedHandler(energy)); - if (!handler->guaranteesStackOnlySeeds()) - throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); - return new BasePairProbabilities(energy.size1(),energy.size2(),!outPrefix2streamName.at(OP_bpsvg).empty()); + if (!noSeedRequired) { + std::unique_ptr handler(getSeedHandler(energy)); + if (!handler->guaranteesStackOnlySeeds()) + throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); + } + return new BasePairProbabilities(energy.size1(),energy.size2(),!noSeedRequired && !outPrefix2streamName.at(OP_bpsvg).empty()); } void CommandLineParsing::writeBasePairProbabilities(const BasePairProbabilities & result,const InteractionEnergy & energy, const Interaction * mfe) const { @@ -2518,11 +2520,13 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output, BasePair { std::unique_ptr pairSeeds; if (pairProbabilities) { - if(model.val!='P' || mode.val!='M' || noSeedRequired || !rri.empty() || windowWidth.val!=0) - throw std::invalid_argument("bpProb/bpsvg requires --model=P --mode=M with stack-only seeds and no windows"); - pairSeeds.reset(getSeedHandler(energy)); - if(!pairSeeds->guaranteesStackOnlySeeds()) - throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); + if(model.val!='P' || mode.val!='M' || !rri.empty() || windowWidth.val!=0) + throw std::invalid_argument("bpProb/bpsvg requires --model=P --mode=M with --noSeed or stack-only seeds and no windows"); + if (!noSeedRequired) { + pairSeeds.reset(getSeedHandler(energy)); + if(!pairSeeds->guaranteesStackOnlySeeds()) + throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); + } } // set up hub for prediction tracking (if needed) PredictionTrackerHub * predTracker = new PredictionTrackerHub(); @@ -2624,7 +2628,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output, BasePair case 'P' : { switch ( mode.val ) { case 'H' : return new PredictorMfeEns2dHeuristic( energy, output, predTracker ); - case 'M' : return new PredictorMfeEns2d( energy, output, predTracker ); + case 'M' : return new PredictorMfeEns2d( energy, output, predTracker, pairProbabilities ); default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } } break; diff --git a/tests/Makefile.am b/tests/Makefile.am index 4d184d6b..a500da8e 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -30,6 +30,7 @@ runApiTests_SOURCES = \ SeededNativeBenchmark_test.cpp \ StackSeedDomain_test.cpp \ BasePairProbabilities_test.cpp \ + PredictorUnseededProbability_test.cpp \ testEasyLoggingSetup.icc \ AccessibilityConstraint_test.cpp \ AccessibilityFromStream_test.cpp \ diff --git a/tests/PredictorUnseededProbability_test.cpp b/tests/PredictorUnseededProbability_test.cpp new file mode 100644 index 00000000..7124d6e5 --- /dev/null +++ b/tests/PredictorUnseededProbability_test.cpp @@ -0,0 +1,217 @@ +#include "catch.hpp" + +#include "IntaRNA/AccessibilityConstraint.h" +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/BasePairProbabilities.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/OutputHandlerInteractionList.h" +#include "IntaRNA/PredictorMfeEns2d.h" +#include "IntaRNA/ReverseAccessibility.h" + +#include +#include +#include +#include +#include + +using namespace IntaRNA; + +namespace { + +struct UnseededOracle { + Z_type z=0; + Matrix mass; + UnseededOracle(size_t n,size_t m) : mass(n,m,0) {} +}; + +// Enumerate increasing pair chains directly, without the predictor's suffix +// recurrence. Filter lonely pairs by inspecting both neighbours in the chain. +UnseededOracle enumerateUnseeded(const InteractionEnergy & energy, + const OutputConstraint & constraint,const IndexRange & r1,const IndexRange & r2) +{ + UnseededOracle result(energy.size1(),energy.size2()); + const size_t end1=r1.to==RnaSequence::lastPos?energy.size1()-1:r1.to; + const size_t end2=r2.to==RnaSequence::lastPos?energy.size2()-1:r2.to; + std::vector> chain; + std::function visit=[&](E_type hybrid) { + const auto [i1,i2]=chain.front(); + const auto [j1,j2]=chain.back(); + bool accepted=true; + if (constraint.noLP) { + for(size_t k=0;k0 && chain[k-1].first+1==chain[k].first + && chain[k-1].second+1==chain[k].second; + const bool right=k+1constraint.maxED || energy.getED2(i2,j2)>constraint.maxED) accepted=false; + const E_type complete=energy.getE(i1,j1,i2,j2,hybrid); + if (accepted && E_isNotINF(complete)) { + const Z_type weight=energy.getBoltzmannWeight(complete); + result.z+=weight; + for (const auto & p:chain) { + const auto bp=energy.getBasePair(p.first,p.second); + result.mass(bp.first,bp.second)+=weight; + } + } + for(size_t k1=j1+1;k1<=end1;++k1) for(size_t k2=j2+1;k2<=end2;++k2) { + if (!energy.isAccessible1(k1) || !energy.isAccessible2(k2) + || !energy.areComplementary(k1,k2)) continue; + const E_type loop=energy.getE_interLeft(j1,k1,j2,k2); + if (E_isINF(loop)) continue; + chain.emplace_back(k1,k2);visit(hybrid+loop);chain.pop_back(); + } + }; + for(size_t i1=r1.from;i1<=end1;++i1) for(size_t i2=r2.from;i2<=end2;++i2) { + if (!energy.isAccessible1(i1) || !energy.isAccessible2(i2) + || !energy.areComplementary(i1,i2)) continue; + chain.emplace_back(i1,i2);visit(energy.getE_init());chain.pop_back(); + } + return result; +} + +void compareUnseeded(const InteractionEnergy & energy,const OutputConstraint & constraint, + const IndexRange & r1,const IndexRange & r2,bool compareOracle=true) +{ + OutputHandlerInteractionList ordinary(constraint,10), annotated(constraint,10); + PredictorMfeEns2d baseline(energy,ordinary,nullptr); + BasePairProbabilities result(energy.size1(),energy.size2(),true); + PredictorMfeEns2d predictor(energy,annotated,nullptr,&result); + baseline.predict(r1,r2); + predictor.predict(r1,r2); + result.finalize(); + REQUIRE(predictor.getZall()==baseline.getZall()); + REQUIRE(result.getZ()==baseline.getZall()); + REQUIRE(std::distance(ordinary.begin(),ordinary.end())==std::distance(annotated.begin(),annotated.end())); + auto a=ordinary.begin(), b=annotated.begin(); + for(;a!=ordinary.end();++a,++b) { + REQUIRE((*a)->energy==(*b)->energy); + REQUIRE((*a)->basePairs==(*b)->basePairs); + } + for(size_t i=0;igetMaxLength()+1?ED_UPPER_BOUND:E_type((to-from+1)*perNt); + } +private: + E_type perNt; +}; + +class TinyInitiationEnergy : public InteractionEnergyBasePair { +public: + using InteractionEnergyBasePair::InteractionEnergyBasePair; + E_type getE_init() const override { return Ekcal_2_E(60); } +}; + +class ChangedObjectivePredictor : public PredictorMfeEns2d { +public: + using PredictorMfeEns2d::PredictorMfeEns2d; +protected: + void updateZ(size_t i1,size_t j1,size_t i2,size_t j2,Z_type z,bool hybrid) override { + PredictorMfeEns::updateZ(i1,j1,i2,j2,2*z,hybrid); + } +}; + +} + +TEST_CASE("exact unseeded pair probabilities enumerate physical pair chains", "[UnseededProbability][PredictorTinyOracle]") { + #include "testEasyLoggingSetup.icc" + for (const auto & sequence:std::vector>{ + {"GGGGGG","CCCCCC"},{"GUGCA","UGCGC"},{"AAGGGAA","ACCCUAA"}, + {"AAAA","AAAA"},{"G","C"}}) { + for(bool noLP:{false,true}) for(bool noGU:{false,true}) for(size_t loop:{0,2}) { + CAPTURE(sequence.first,sequence.second,noLP,noGU,loop); + RnaSequence target("target",sequence.first), query("query",sequence.second); + AccessibilityDisabled targetAcc(target,0,nullptr), queryAcc(query,0,nullptr); + ReverseAccessibility reverse(queryAcc); + InteractionEnergyBasePair energy(targetAcc,reverse,loop,loop,false,0.7,Ekcal_2_E(-0.8),3,Ekcal_2_E(0.13)); + OutputConstraint constraint(4,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP,noGU,true,true); + compareUnseeded(energy,constraint,{0,RnaSequence::lastPos},{0,RnaSequence::lastPos}); + } + } +} + +TEST_CASE("unseeded pair probabilities preserve filtered offset ensembles", "[UnseededProbability]") { + #include "testEasyLoggingSetup.icc" + RnaSequence target("target","GUGGGCG"),query("query","CUCCCCC"); + for(bool noLP:{false,true}) for(bool noGU:{false,true}) { + ProbabilityWidthAccessibility targetAcc(target,3,Ekcal_2_E(0.10)); + ProbabilityWidthAccessibility queryAcc(query,4,Ekcal_2_E(0.15)); + ReverseAccessibility reverse(queryAcc); + InteractionEnergyBasePair energy(targetAcc,reverse,1,2,false,1,Ekcal_2_E(-1),3,0,true,false); + OutputConstraint constraint(5,OutputConstraint::OVERLAP_NONE,E_INF,E_INF,false,noLP,noGU,true,true,Ekcal_2_E(0.45)); + compareUnseeded(energy,constraint,{1,5},{2,6}); + compareUnseeded(energy,constraint,{2,RnaSequence::lastPos},{1,RnaSequence::lastPos}); + } + AccessibilityConstraint blocked(target,"...b...",0,"","",""); + AccessibilityDisabled targetAcc(target,0,&blocked),queryAcc(query,0,nullptr); + ReverseAccessibility reverse(queryAcc); + InteractionEnergyBasePair energy(targetAcc,reverse,2,2); + OutputConstraint constraint(1,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,false,false,true,true); + compareUnseeded(energy,constraint,{1,5},{1,5}); + + BasePairProbabilities result(target.size(),query.size()); + OutputHandlerInteractionList output(constraint,4); + PredictorMfeEns2d predictor(energy,output,nullptr,&result); + predictor.predict({0,2},{0,2}); + predictor.predict({4,6},{3,6}); + result.finalize(); + const auto first=enumerateUnseeded(energy,constraint,{0,2},{0,2}); + const auto second=enumerateUnseeded(energy,constraint,{4,6},{3,6}); + REQUIRE(result.getZ()==Approx(first.z+second.z).epsilon(2e-12)); + for(size_t i=0;i0 for v in covered[3][1:]);assert all(v==0 for v in with_spot[2]) compressed=d/'pairs.csv.gz';run(options=[f'--out=bpProb:{compressed}']);same(matrix(compressed)[2],values) # Unsupported requests fail explicitly and never create a probability file. - for opts in (['--mode=H'],['--mode=S'],['--model=X','--mode=K'],['--noSeed'],['--windowWidth=3','--windowOverlap=2'],['--seedMaxUP=1'],['--seedTQ=1|.|&2||']): + for opts in (['--mode=H'],['--mode=S'],['--model=X','--mode=K'],['--windowWidth=3','--windowOverlap=2'],['--seedMaxUP=1'],['--seedTQ=1|.|&2||']): fail=d/'unsupported.csv';run(options=[f'--out=bpProb:{fail}',*opts],success=False);assert not fail.exists() fail=d/'range.csv';run('GGGG','CCCC',[f'--out=bpProb:{fail}',*seedopt,'--energyAdd=900'],success=False);assert not fail.exists() run(options=['--out=bpProb:/dev/full'],success=False) diff --git a/tests/testBasePairProbabilitySvg.py b/tests/testBasePairProbabilitySvg.py index 597afd6a..0cbb64b6 100644 --- a/tests/testBasePairProbabilitySvg.py +++ b/tests/testBasePairProbabilitySvg.py @@ -122,6 +122,11 @@ def plot(t='GGAGGG', q='CCCC', opts=()): assert float(outline.attrib['width']) == float(matrix[ts, qe].attrib['x'])-float(matrix[ts, qs].attrib['x'])+20 assert float(outline.attrib['height']) == float(matrix[ts, qs].attrib['y'])-float(matrix[te, qs].attrib['y'])+20 assert 'Green outline:' in ''.join(mfe_plot.itertext()) + unseeded = plot(opts=['--noSeed', f'--out=bpProb:{table}']) + assert not seeds(unseeded) + rows = list(csv.reader(table.open(), delimiter=';')) + assert all(float(pairs(unseeded)[i,j].attrib['data-probability']) == float(rows[i][j]) + for i in range(1,7) for j in range(1,5)) # Explicit singleton/mixed seed masks exclude extension-only cells. opts = ['--seedTQ=1|&4|,2||&2||'] mixed = plot('GGGG', 'CCCC', opts) @@ -214,7 +219,7 @@ def check_frame(root, target_pu, query_pu): run(options=[f'--out=bpsvg:{compressed}']) assert len(pairs(load(compressed))) == 24 # Validation/numerical failures must not publish a partial document. - for opts in (['--mode=H'], ['--noSeed'], ['--seedMaxUP=1'], ['--model=X', '--mode=K'], + for opts in (['--mode=H'], ['--seedMaxUP=1'], ['--model=X', '--mode=K'], ['--windowWidth=3'], ['--seedTQ=1|&4|','--energyAdd=900']): failure = d/'failure.svg' run(options=[f'--out=bpsvg:{failure}', *opts], success=False) From ae0eb78da255c76765e6745d2dc4040766f7b791 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Fri, 9 Oct 2026 22:30:47 +0200 Subject: [PATCH 6/9] feat: add approximate heuristic pair probabilities (step 11) --- ChangeLog | 22 ++ README.md | 32 ++- .../seeded-probabilities-validation.md | 48 ++++ .../IntaRNA.PredictorMfeEns2d.outside.svg | 69 ++++++ ...edictorMfeEns2dHeuristic.probabilities.svg | 80 ++++++ src/IntaRNA/BasePairProbabilities.cpp | 4 + src/IntaRNA/BasePairProbabilities.h | 8 + src/IntaRNA/BasePairProbabilitySvg.cpp | 7 +- src/IntaRNA/BasePairProbabilityWriter.cpp | 3 +- src/IntaRNA/PartitionArithmetic.h | 12 +- src/IntaRNA/PredictorMfeEns2d.cpp | 5 +- src/IntaRNA/PredictorMfeEns2d.h | 4 +- src/IntaRNA/PredictorMfeEns2dHeuristic.cpp | 99 +++++++- src/IntaRNA/PredictorMfeEns2dHeuristic.h | 36 ++- ...redictorMfeEns2dHeuristicSeedExtension.cpp | 41 ++- .../PredictorMfeEns2dHeuristicSeedExtension.h | 4 +- .../PredictorMfeEns2dSeedExtension.cpp | 211 +++++++++++++++- src/IntaRNA/PredictorMfeEns2dSeedExtension.h | 22 ++ src/bin/CommandLineParsing.cpp | 18 +- tests/Makefile.am | 2 + tests/PredictorHeuristicProbability_test.cpp | 234 ++++++++++++++++++ ...edictorHeuristicSeededProbability_test.cpp | 201 +++++++++++++++ tests/PredictorUnseededProbability_test.cpp | 6 + tests/testBasePairProbabilities.py | 12 +- tests/testBasePairProbabilitySvg.py | 11 +- 25 files changed, 1137 insertions(+), 54 deletions(-) create mode 100644 doc/recursions/IntaRNA.PredictorMfeEns2d.outside.svg create mode 100644 doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.probabilities.svg create mode 100644 tests/PredictorHeuristicProbability_test.cpp create mode 100644 tests/PredictorHeuristicSeededProbability_test.cpp diff --git a/ChangeLog b/ChangeLog index e818a899..eaea2fe1 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,12 @@ ## Interface and handling +- enable approximate actual-pair CSV/SVG probabilities for unseeded and seeded + heuristic ensembles, labelled explicitly as approximate (PR #259, step 11) + +- reject nonfinite partition totals and finite-energy boundary underflow before + publishing exact unseeded or heuristic pair probabilities (PR #259 review) + - enable base-pair CSV/SVG probabilities for exact unseeded model P/mode M predictions (PR #259, step 10) @@ -175,6 +181,22 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261009 Alexander Mitrofanov + * PredictorMfeEns2dHeuristic, PredictorMfeEns2dHeuristicSeedExtension, + PredictorMfeEns2dSeedExtension, BasePairProbabilities, probability writers : + + reverse retained candidate chains and signed seed-extension recurrences to + obtain marginals of the existing heuristic ensemble; preserve denominator + and seed-overlap corrections, validate all masses before publication + + expose approximation metadata in CSV/SVG and retain seed provenance + * bound cropped noLP continuations by the active region and correct default + region offsets; reject unsupported singleton heuristic seeds explicitly + * CommandLineParsing, tests, README.md, doc/recursions : + + enable model P/mode H probability output; add independent tiny-chain and + signed-extension checks plus CLI parity/metadata coverage (PR #259 step 11) + * PredictorMfeEns2d, probability validation tests and report : + * reject finite-energy boundary underflow and nonfinite heuristic denominators + before publication; record independent review, full suites and performance + 261009 Alexander Mitrofanov * PredictorMfeEns2d, CommandLineParsing : + compute actual-pair masses by reversing the existing unseeded recurrence; diff --git a/README.md b/README.md index 888f143f..f3089346 100644 --- a/README.md +++ b/README.md @@ -2017,7 +2017,8 @@ or `STDERR` instead of a file name. `--out=bpProb:FILE` writes the probability that each target/query nucleotide pair actually pairs, conditional on an allowed interaction. It requires -`--model=P --mode=M` with `--noSeed` or a stack-only seed handler. For example: +`--model=P --mode=M` with `--noSeed` or a stack-only seed handler, or +`--model=P --mode=H` for an approximate heuristic ensemble. For example: ```sh IntaRNA -t GGAGGG -q CCCC --energy=B --acc=N --seedBP=2 \ @@ -2030,12 +2031,19 @@ is also used without `bpProb`, so requesting the matrix does not change `Zall`. Computed seeds preserve their handler-specific admission thresholds. Explicit stacked seeds may have different lengths, including singletons; competing explicit patterns at the same start still retain only the selected minimum-energy -pattern. Other predictors retain their existing behavior. Heuristic, kinetic, -evaluation, seed-only and bulged-seed `bpProb` requests are rejected. +pattern. With `--mode=H`, probabilities describe the predictor's retained +candidate ensemble, not all possible interactions. Unseeded H retains one +continuation per left pair; seeded H retains its seed-extension candidates and +selected right extension, including the existing seed-overlap corrections. +Seeded H supports computed or explicit seeds with at least two pairs, including +bulged seeds. These heuristic ensembles can retain path multiplicities and +must not be interpreted as the exact unique-chain ensemble of mode M. +Kinetic, evaluation, seed-only and exact bulged-seed probability requests remain +unsupported. Rows follow the original target and columns the original query, both 5' to 3'. -The header starts with `bpProb`; labels contain each nucleotide and its display -index (including `tIdxPos0`/`qIdxPos0` shifts). Entries use sufficient significant +The header starts with `bpProb` (exact) or `bpProbApproximate` (heuristic); labels +contain each nucleotide and its display index (including `tIdxPos0`/`qIdxPos0` shifts). Entries use sufficient significant digits for the partition type, with scientific notation for small probabilities. A successful empty ensemble contains `NA`; pairs absent from a nonempty ensemble have probability zero. The output uses `outSep`, the usual multi-FASTA filename @@ -2061,12 +2069,14 @@ including its numerical zero cutoffs, without changing ordinary prediction. The seeded backend checks nonnegative arithmetic and reports numerical range failures as errors, including underflow, overflow and NaN. Positive subnormals are rejected; there is currently no scaling fallback. Failed or cancelled regional work never -produces a successful-looking probability matrix. Ordinary floating-point -rounding remains; "exact" describes the enumerated ensemble, not real arithmetic. +produces a successful-looking probability matrix. Heuristic reverse passes likewise +validate finite arithmetic and reject invalid masses before publication. Ordinary +floating-point rounding remains; "exact" describes the enumerated ensemble, not +real arithmetic. For library use, keep a `BasePairProbabilities` owner alive across predictions, pass its optional non-owning pointer to `PredictorMfeEns2dSeedExtension` (stack -seeds) or `PredictorMfeEns2d` (unseeded), and call +seeds), `PredictorMfeEns2d` (unseeded), or their heuristic subclasses, and call `finalize()` only after all requested disjoint regions succeed. A writer is called explicitly after finalization; its destructor publishes nothing. The output constraint must enable `needZall`. The [standalone API example](doc/seeded-probabilities-example.cpp) @@ -2078,9 +2088,9 @@ show the public API and measured limits. ### SVG base-pair probability dot plots `--out=bpsvg:FILE` renders the same actual-pair probabilities as `bpProb` in a -standalone SVG. Both outputs can be requested together. The same exact -`--model=P --mode=M` unseeded/stack-seed restrictions, regional ensemble, normalization -and success-only publication rules apply. +standalone SVG. Both outputs can be requested together. The same supported +model/mode combinations, regional ensemble, normalization and success-only +publication rules apply. Heuristic SVGs are explicitly labelled approximate. ```sh IntaRNA -t GGAGGG -q CCCC --energy=B --acc=N --seedBP=2 \ diff --git a/doc/analysis/seeded-probabilities-validation.md b/doc/analysis/seeded-probabilities-validation.md index 3449a5ba..56c83716 100644 --- a/doc/analysis/seeded-probabilities-validation.md +++ b/doc/analysis/seeded-probabilities-validation.md @@ -270,3 +270,51 @@ Seven interleaved runs gave wall medians 0.412419→0.390979 seconds (ranges RSS 87,744→84,544 KiB. The slowdown did not reproduce without filesystem writes; the fixed serializer was 5.2% faster in this specific control. Use the same sparse400 command with `--out=STDOUT --out=bpsvg:/dev/null > /dev/null` to repeat. + +## PR #259 review: SVG revisions and steps 10–11 (2026-10-09) + +Martin's review adds unseeded and heuristic probability outputs to the original +SVG/documentation scope. Exact unseeded mode M now reverses each fixed-right-end +inside table. Unseeded mode H reverses its retained continuation chains; seeded +mode H traces its signed left/right extension arithmetic, including the existing +seed-overlap corrections and selected right extensions. Heuristic output is +explicitly approximate (`bpProbApproximate` and an approximate SVG title). Its +candidate ensemble can retain path multiplicities and is not the exact seeded +unique-chain ensemble. Existing numerical-zero cutoffs remain in the legacy +unseeded/heuristic objectives. + +The SVG now starts both strands at the bottom-left origin, places the query name +below the matrix, uses sequence names and rounded values in tooltips, displays +intramolecular pairing as `1-Pu`, and outlines the best reported interaction region +in green. Machine-readable probability attributes keep full precision. The MFE +outline is available with `--outNumber=0`; ordinary ranked output stays suppressed. +An outline is omitted if no interaction passes the reporting filters. + +Independent agents implemented and reviewed the recurrences, checked actual-pair +ownership, and evaluated performance before push. The new API tests compare exact +and unseeded heuristic results with independently enumerated tiny-chain ensembles. +Seeded heuristic tests change one pair's Boltzmann fugacity and compare the change +in the scalar partition function with that pair's numerator. They freeze the +baseline selected right extension and concrete computed seed family, so the test +checks marginals of the implemented conditional ensemble rather than a different +family selected after perturbation. Cases include overlapping/disjoint and bulged +seeds, LP/noLP, regions, seed annotations, successful empty results, singleton seed +rejection, and failure publication. Separate symbolic checks of 500 randomized +small legacy seeded extension polynomials found no negative or over-total pair +masses in the sampled cases. + +The review found and fixed active-region offset/dimension errors, finite-boundary +Boltzmann underflow masquerading as an empty result, and a nonfinite seeded-H +denominator bypassing a relative-tolerance comparison. These fixes have regression +coverage. CSV/SVG CLI tests cover approximation metadata, numeric parity, seed +masks, MFE geometry, signed coordinates and unchanged ordinary output. + +The final GCC 14.4 release and debug runs each pass all seven `make tests -j2` +suites and 100,947 API assertions in 97 cases. All 74 installed public headers +compile independently. Installed pkg-config consumers for exact/heuristic, +seeded/unseeded predictors compile, link and produce CSV/SVG with the expected +approximation metadata. `make dist` contains the new predictor tests +and implementation. All revised text is UTF-8/LF and `git diff --check` passes. +Both new recurrence SVGs parse, render and passed independent equation review; +representative output SVGs were rendered and visually inspected. Apple Clang was +not available locally. diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2d.outside.svg b/doc/recursions/IntaRNA.PredictorMfeEns2d.outside.svg new file mode 100644 index 00000000..a44c7205 --- /dev/null +++ b/doc/recursions/IntaRNA.PredictorMfeEns2d.outside.svg @@ -0,0 +1,69 @@ + +Unseeded pair probabilities · fixed-right outside pass +PredictorMfeEns2d | --model=P --mode=M --noSeed | reverse the existing hybrid recurrence. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Unseeded pair probabilities · fixed-right outside pass +PredictorMfeEns2d | --model=P --mode=M --noSeed | reverse the existing hybrid recurrence + +A pair i=(i₁,i₂) uses local energy coordinates; RNA 2 is reversed. i ≺ k increases both coordinates. +Fix right pair j. Hⱼ(i) is the forward hybrid weight; Aⱼ(i)=∂Z/∂Hⱼ(i) is its outside adjoint. +M(i) sums complete-interaction weights containing actual pair i, over every right boundary j. +d=BW(E_init); t(i,k)=BW(E_interLeft(i,k)); i⁺=i+(1,1); s=t(i,i⁺) + +1 Initialize the same completed-boundary objective +Aⱼ(i) ← B(i,j); Z += B(i,j)Hⱼ(i) +B(i,j)=BW(energy.getE(i₁,j₁,i₂,j₂,0)) only when !Z_equal(Hⱼ(i),0) and isValidOutputSite() pass. +Otherwise B=0. Complete ED / GU-end filters and native terminal terms match updateZ(). + +2 Reverse i in increasing order; credit its explicit pair +M(i) += Aⱼ(i)Hⱼ(i) +Every represented suffix contains its left pair i exactly once. Process only nonzero Aⱼ(i) and Hⱼ(i). +Then propagate Aⱼ(i) through each forward term below; later suffix states receive all enclosing mass. + +3a Lonely pairs allowed +Hⱼ(j)=d +Aⱼ(k) += Aⱼ(i)t(i,k) +Direct terminal: k=j; its pair is counted at j. +Intermediate k: i ≺ k ≺ j, allowed loop, +and !Z_equal(Hⱼ(k),0) as in the forward sum. + +3b noLP: direct stacking +Aⱼ(i⁺) += Aⱼ(i)s if i⁺ ≠ j +M(j) += Aⱼ(i)s d if i⁺ = j +A noLP terminal has Hⱼ(j)=0. +Its second pair j is implicit in the two-pair +initialization and is credited separately. + +3c noLP: stack, then a nonstack loop + +Hⱼ(k) + + + +i + +i⁺ + +k + +j +a = Aⱼ(i)s t(i⁺,k) +Aⱼ(k) += a +M(i⁺) += a Hⱼ(k) +Require complementary i⁺; i⁺ ≺ k ≺ j; k ≠ i⁺+(1,1); the same loop limits and epsilon gate on Hⱼ(k). +Here i⁺ is an actual pair skipped by the DP-state jump. Unpaired loop positions receive no mass. + +4 Commit successful regions and normalize only once +P(i,j) = ∑ M_region(i,j) / ∑ Z_region +Check the accumulated objective equals the unchanged Zall; map pairs with energy.getBasePair(). +Original target/query orientation is stored in the sink. Z=0 → NA; failed arithmetic → no result. +Conditional on an allowed interaction; no unbound-state weight. No outside buffers when output is disabled. +Exact for the existing forward recurrence, including its numerical-zero gates; no alternative ensemble is substituted. + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +PredictorMfeEns2d.cpp: fillHybridZ(), accumulatePairMasses(), predictRegion() +PredictorMfeEns.cpp: addPartitionContribution(); BasePairProbabilities: addRegion(), finalize() + diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.probabilities.svg b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.probabilities.svg new file mode 100644 index 00000000..3a11d160 --- /dev/null +++ b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristic.probabilities.svg @@ -0,0 +1,80 @@ + +Unseeded heuristic probabilities · candidate flow +PredictorMfeEns2dHeuristic | --model=P --mode=H --noSeed | approximate retained ensemble. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. + + + +INTARNA / RECURRENCE NOTES +Unseeded heuristic probabilities · candidate flow +PredictorMfeEns2dHeuristic | --model=P --mode=H --noSeed | approximate retained ensemble + +Each retained cell C(i) represents one selected pair chain, with hybrid weight vᵢ and right boundary jᵢ. +Pair coordinates increase on both strands internally; RNA 2 is reversed. Set i⁺=i+(1,1). +d=BW(E_init), t(i,k)=BW(E_interLeft(i,k)), and B(i,j)=BW(energy.getE(i₁,j₁,i₂,j₂,0)). + +Approximate means the continuation ensemble is pruned +Visit i in decreasing order. Generate candidates from already selected child chains C(k). +Every admitted candidate contributes to Z; only the candidate with the lowest complete site energy +is retained as C(i) for later extensions. Strict improvement wins; native energy rounding and ties are preserved. + +1 Record which new pairs each candidate owns +Candidate +Hybrid weight v +Owned pairs S +Child k +LP: initial pair +d +{i} +none +LP: loop / stack to k +t(i,k)vₖ +{i} +k +noLP: initial stack +s d, s=t(i,i⁺) +{i, i⁺} +none +noLP: direct stacking +s vᵢ⁺ +{i} +i⁺ +noLP: nonstack loop to k +s t(i⁺,k)vₖ +{i, i⁺} +k +Native pair, loop/span and right-end GU constraints apply; noLP loops exclude the adjacent child k=i⁺+(1,1). + +2 Initialize mass and flow from every admitted candidate +w = B(i,j)v; Z += w; M(p) += w for p ∈ S; F(k) += w if a child exists +M starts at zero; F is accumulated incoming weight for a retained child chain. This includes unselected candidates. +Use exactly updateZ() admission: !Z_equal(v,0), complete-site GU / ED filters, and finite arithmetic. +The current candidate owns S; every remaining pair belongs to its selected child continuation. + +3 Propagate incoming mass through selected chains + +C(k) + + + +i + +k + +jₖ +Sweep i in increasing order. Let Sᵢ and kᵢ describe C(i). +M(p) += F(i) for every p ∈ Sᵢ; +F(kᵢ) += F(i) when C(i) has a child. +Every edge increases both indices, so all parents arrive first. +No division by vᵢ is needed: a retained cell represents one chain, so its full incoming weight follows every pair. + +4 Publish marginals of the retained candidate ensemble +P(i,j) = ∑ M_region(i,j) / ∑ Z_region +Map energy coordinates back to original strand order, merge disjoint regions, then validate and normalize. +The denominator must equal the unchanged heuristic Zall. Mark the result approximate, including empty results. +Missing continuation chains are not recovered. These marginals differ from the exact --mode=M ensemble. +Only actual paired nucleotides receive mass; positions merely enclosed by a retained interaction remain unmarked. + +Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion +PredictorMfeEns2dHeuristic.cpp: fillHybridZ(), recordProbabilityCandidate(), commitProbabilityRegion() +BasePairProbabilities: markApproximate(), addRegion(), finalize(); PredictorMfeEns.cpp: updateZ() + diff --git a/src/IntaRNA/BasePairProbabilities.cpp b/src/IntaRNA/BasePairProbabilities.cpp index 0496e534..4d06f327 100644 --- a/src/IntaRNA/BasePairProbabilities.cpp +++ b/src/IntaRNA/BasePairProbabilities.cpp @@ -4,6 +4,10 @@ #include namespace IntaRNA { +void BasePairProbabilities::markApproximate() { + if (state!=Status::pending) throw std::logic_error("pair probabilities: accumulator is not pending"); + approximate=true; +} BasePairProbabilities::BasePairProbabilities(size_t n,size_t m,bool collectSeedPairs) : mass(n,m,0), collectSeeds(collectSeedPairs), seeds(collectSeedPairs?n:0,collectSeedPairs?m:0,0) {} void BasePairProbabilities::addRegion(const IndexRange & t,const IndexRange & q, diff --git a/src/IntaRNA/BasePairProbabilities.h b/src/IntaRNA/BasePairProbabilities.h index bae84b8d..9e62d048 100644 --- a/src/IntaRNA/BasePairProbabilities.h +++ b/src/IntaRNA/BasePairProbabilities.h @@ -49,6 +49,12 @@ class BasePairProbabilities { const Matrix & rawMasses() const; /** @return whether this sink requests admitted seed-cell annotations */ bool collectsSeedPairs() const; + /** Mark an ensemble restricted by heuristic continuation choices. + * Must be called while pending, before publication. + */ + void markApproximate(); + /** @return whether the represented interaction ensemble is heuristic */ + bool isApproximate() const; /** @return seed-cell mask in original coordinates, empty when not requested */ const Matrix & seedPairs() const; /** @return validated matrix; only available for finalized nonempty results. @@ -57,6 +63,7 @@ class BasePairProbabilities { Matrix probabilities() const; private: Status state=Status::pending; + bool approximate=false; Z_type z=0; Matrix mass; const bool collectSeeds; @@ -67,6 +74,7 @@ inline BasePairProbabilities::Status BasePairProbabilities::status() const { ret inline Z_type BasePairProbabilities::getZ() const { return z; } inline const Matrix & BasePairProbabilities::rawMasses() const { return mass; } inline bool BasePairProbabilities::collectsSeedPairs() const { return collectSeeds; } +inline bool BasePairProbabilities::isApproximate() const { return approximate; } inline const Matrix & BasePairProbabilities::seedPairs() const { return seeds; } inline void BasePairProbabilities::fail() noexcept { state=Status::failed; } } diff --git a/src/IntaRNA/BasePairProbabilitySvg.cpp b/src/IntaRNA/BasePairProbabilitySvg.cpp index d350547f..5b7620ee 100644 --- a/src/IntaRNA/BasePairProbabilitySvg.cpp +++ b/src/IntaRNA/BasePairProbabilitySvg.cpp @@ -94,9 +94,10 @@ BasePairProbabilityWriter::svgBlock(const BasePairProbabilities & result,const I out<::max_digits10); out<<"\n" - <<"Base-pair probabilities: "<\n" + <<""<<(result.isApproximate()?"Approximate base-pair probabilities: ":"Base-pair probabilities: ")<<xml(t.getId())<<" / "<<xml(q.getId())<<"\n" <<"Target rows run 5' to 3' upwards and query columns run 5' to 3' to the right. Pair probabilities are conditional on an allowed interaction. " - <<"The frame shows single-nucleotide intramolecular pairing probabilities (1 - Pu) from the accessibility model.\n" + <<"The frame shows single-nucleotide intramolecular pairing probabilities (1 - Pu) from the accessibility model. " + <<(result.isApproximate()?"Approximate probabilities describe the retained heuristic ensemble.":"Probabilities describe the exact predictor ensemble.")<<"\n" < text { font-family: sans-serif; fill: #243247; } .heading { font-size: 18px; font-weight: bold; } @@ -126,7 +127,7 @@ text { font-family: sans-serif; fill: #243247; } auto text=[&](double tx,double ty,const std::string & value,const std::string & cls,const char * anchor="middle") { out<<""<\n"; }; - text(width/2,30,"Base-pair probabilities","heading"); + text(width/2,30,result.isApproximate()?"Approximate base-pair probabilities":"Base-pair probabilities","heading"); if (t.getId()!="target" || q.getId()!="query") text(width/2,53,t.getId()+" / "+q.getId(),"subtitle"); text(x+w/2,y+h+65,q.getId()+" (5' to 3')","axis"); diff --git a/src/IntaRNA/BasePairProbabilityWriter.cpp b/src/IntaRNA/BasePairProbabilityWriter.cpp index d38a744b..b8ee89bf 100644 --- a/src/IntaRNA/BasePairProbabilityWriter.cpp +++ b/src/IntaRNA/BasePairProbabilityWriter.cpp @@ -18,7 +18,8 @@ std::string BasePairProbabilityWriter::block(const BasePairProbabilities & resul // Validate every normalized value before opening a file or writing a header. const auto p=empty?Matrix():result.probabilities(); std::ostringstream data; - data<::max_digits10)<<"bpProb"; + data<::max_digits10) + <<(result.isApproximate()?"bpProbApproximate":"bpProb"); for(size_t j=0;j namespace IntaRNA { -/** Checked nonnegative arithmetic for exact seeded partitions only. +/** Checked nonnegative arithmetic for pair-probability partitions. * Structural zero is exactly zero. Positive subnormals are rejected because * their relative precision is not guaranteed. Ordinary sum rounding remains. */ @@ -24,7 +24,7 @@ class PartitionArithmetic { inline Z_type PartitionArithmetic::check(const Z_type x) { if (!(x>=0 && x<=std::numeric_limits::max()) || (x>0 && x::min())) - throw std::range_error("exact seeded partition: numerical range exceeded"); + throw std::range_error("partition arithmetic: numerical range exceeded"); return x; } inline Z_type PartitionArithmetic::add(const Z_type a,const Z_type b) { @@ -34,19 +34,19 @@ inline Z_type PartitionArithmetic::multiply(const Z_type a,const Z_type b) { check(a); check(b); if (a==0 || b==0) return 0; const Z_type r=check(a*b); - if (r==0) throw std::range_error("exact seeded partition: product underflow"); + if (r==0) throw std::range_error("partition arithmetic: product underflow"); return r; } inline Z_type PartitionArithmetic::divide(const Z_type a,const Z_type b) { check(a); check(b); - if (b==0) throw std::range_error("exact seeded partition: zero denominator"); + if (b==0) throw std::range_error("partition arithmetic: zero denominator"); const Z_type r=check(a/b); - if (a>0 && r==0) throw std::range_error("exact seeded partition: quotient underflow"); + if (a>0 && r==0) throw std::range_error("partition arithmetic: quotient underflow"); return r; } inline Z_type PartitionArithmetic::exp(const Z_type exponent) { const Z_type r=check(Z_exp(exponent)); - if (r==0) throw std::range_error("exact seeded partition: exponential underflow"); + if (r==0) throw std::range_error("partition arithmetic: exponential underflow"); return r; } } diff --git a/src/IntaRNA/PredictorMfeEns2d.cpp b/src/IntaRNA/PredictorMfeEns2d.cpp index 490bfe9b..a03506da 100644 --- a/src/IntaRNA/PredictorMfeEns2d.cpp +++ b/src/IntaRNA/PredictorMfeEns2d.cpp @@ -256,7 +256,10 @@ PredictorMfeEns2d::accumulatePairMasses(size_t j1,size_t j2,Z2dMatrix & outside, outside(i1,i2)=0; const Z_type h=A::check(hybridZ(i1,i2)); if (!Z_equal(h,0) && isValidOutputSite(i1,j1,i2,j2)) { - const Z_type b=A::check(energy.getBoltzmannWeight(energy.getE(i1,j1,i2,j2,E_type(0)))); + const E_type boundaryEnergy=energy.getE(i1,j1,i2,j2,E_type(0)); + const Z_type b=A::check(energy.getBoltzmannWeight(boundaryEnergy)); + if (b==0 && E_isNotINF(boundaryEnergy)) + throw std::range_error("unseeded pair probabilities: boundary weight underflow"); outside(i1,i2)=b; denominator=A::add(denominator,A::multiply(h,b)); } diff --git a/src/IntaRNA/PredictorMfeEns2d.h b/src/IntaRNA/PredictorMfeEns2d.h index 169b9881..2b77bb67 100644 --- a/src/IntaRNA/PredictorMfeEns2d.h +++ b/src/IntaRNA/PredictorMfeEns2d.h @@ -36,7 +36,9 @@ class PredictorMfeEns2d: public PredictorMfeEns { * tracking is to be done; if non-NULL, the tracker gets deleted * on this->destruction. * @param pairProbabilities optional non-owning sink; requires needZall and - * must remain pending until all disjoint regions have succeeded + * must remain pending until all disjoint regions have succeeded. + * Collection reverses this class's fillHybridZ recurrence; subclasses + * changing that recurrence must also supply their own collection. */ PredictorMfeEns2d( const InteractionEnergy & energy , OutputHandler & output diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp index 806ea98a..93b57cbf 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristic.cpp @@ -1,5 +1,6 @@ #include "IntaRNA/PredictorMfeEns2dHeuristic.h" +#include "IntaRNA/PartitionArithmetic.h" #include @@ -11,8 +12,9 @@ PredictorMfeEns2dHeuristic:: PredictorMfeEns2dHeuristic( const InteractionEnergy & energy , OutputHandler & output - , PredictionTracker * predTracker ) - : PredictorMfeEns2d(energy,output,predTracker) + , PredictionTracker * predTracker + , BasePairProbabilities * pairProbabilities ) + : PredictorMfeEns2d(energy,output,predTracker,pairProbabilities) { } @@ -33,6 +35,14 @@ PredictorMfeEns2dHeuristic:: predict( const IndexRange & r1 , const IndexRange & r2 ) +{ + try { predictRegionHeuristic(r1,r2); } + catch (...) { if (pairProbabilities) pairProbabilities->fail(); throw; } +} + +void +PredictorMfeEns2dHeuristic:: +predictRegionHeuristic(const IndexRange & r1,const IndexRange & r2) { #if INTARNA_MULITHREADING #pragma omp critical(intarna_omp_logOutput) @@ -48,15 +58,28 @@ predict( const IndexRange & r1 #endif + if (!r1.isAscending() || !r2.isAscending() + || r1.from>=energy.getAccessibility1().getSequence().size() + || r2.from>=energy.getAccessibility2().getSequence().size()) + throw std::invalid_argument("heuristic ensemble prediction: invalid region"); + if (pairProbabilities && pairProbabilities->status()!=BasePairProbabilities::Status::pending) + throw std::logic_error("heuristic pair probabilities: accumulator is not pending"); + if (pairProbabilities) pairProbabilities->markApproximate(); // set index offset energy.setOffset1(r1.from); energy.setOffset2(r2.from); // resize matrix hybridZ.resize( std::min( energy.size1() - , (r1.to==RnaSequence::lastPos?energy.size1()-1:r1.to)-r1.from+1 ) + , r1.to==RnaSequence::lastPos?energy.size1():r1.to-r1.from+1 ) , std::min( energy.size2() - , (r2.to==RnaSequence::lastPos?energy.size2()-1:r2.to)-r2.from+1 ) ); + , r2.to==RnaSequence::lastPos?energy.size2():r2.to-r2.from+1 ) ); + if (pairProbabilities) { + probabilityContinuation = Matrix(hybridZ.size1(),hybridZ.size2()); + probabilityFlow = Matrix(hybridZ.size1(),hybridZ.size2(),0); + probabilityMass = Matrix(hybridZ.size1(),hybridZ.size2(),0); + probabilityDenominator = 0; + } // init mfe for later updates initOptima(); @@ -68,6 +91,7 @@ predict( const IndexRange & r1 // trace back and output handler update reportOptima(); + if (pairProbabilities) commitProbabilityRegion(); } @@ -109,8 +133,8 @@ fillHybridZ() // no lp allowed if (noLpShift != 0) { // check if right-side stacking of (i1,i2) is possible - if ( i1+noLpShift < energy.size1() - && i2+noLpShift < energy.size2() + if ( i1+noLpShift < hybridZ.size1() + && i2+noLpShift < hybridZ.size2() && energy.areComplementary(i1+noLpShift,i2+noLpShift)) { // get stacking term to avoid recomputation @@ -130,6 +154,10 @@ fillHybridZ() curCellEtotal = energy.getE(i1,i1+noLpShift, i2,i2+noLpShift ,energy.getE(curCell->val)); // update overall partition function information for initial bps only updateZ( i1,curCell->j1, i2,curCell->j2, curCell->val, true ); + if (pairProbabilities) { + probabilityContinuation(i1,i2)={RnaSequence::lastPos,RnaSequence::lastPos,noLpShift!=0}; + recordProbabilityCandidate(i1,curCell->j1,i2,curCell->j2,curCell->val,probabilityContinuation(i1,i2)); + } } @@ -150,6 +178,8 @@ fillHybridZ() // update overall partition function information for current right extension updateZ( i1,rightExt->j1, i2,rightExt->j2, curZ, true ); + if (pairProbabilities) + recordProbabilityCandidate(i1,rightExt->j1,i2,rightExt->j2,curZ,{i1+noLpShift,i2+noLpShift,false}); // check if this combination yields better energy curEtotal = energy.getE(i1,rightExt->j1, i2,rightExt->j2, energy.getE(curZ)); @@ -164,6 +194,8 @@ fillHybridZ() curCell->val = curZ; // store total energy to avoid recomputation curCellEtotal = curEtotal; + if (pairProbabilities) + probabilityContinuation(i1,i2)={i1+noLpShift,i2+noLpShift,false}; } } } @@ -196,6 +228,8 @@ fillHybridZ() // update overall partition function information for current right extension updateZ( i1,rightExt->j1, i2,rightExt->j2, curZ, true ); + if (pairProbabilities) + recordProbabilityCandidate(i1,rightExt->j1,i2,rightExt->j2,curZ,{i1+noLpShift+w1,i2+noLpShift+w2,noLpShift!=0}); // check if this combination yields better energy curEtotal = energy.getE(i1,rightExt->j1, i2,rightExt->j2, energy.getE(curZ)); @@ -210,6 +244,8 @@ fillHybridZ() curCell->val = curZ; // store total energy to avoid recomputation curCellEtotal = curEtotal; + if (pairProbabilities) + probabilityContinuation(i1,i2)={i1+noLpShift+w1,i2+noLpShift+w2,noLpShift!=0}; } } // w2 @@ -224,4 +260,55 @@ fillHybridZ() //////////////////////////////////////////////////////////////////////////// +void PredictorMfeEns2dHeuristic::recordProbabilityCandidate(size_t i1,size_t j1, + size_t i2,size_t j2,Z_type hybrid,const ProbabilityContinuation & continuation) +{ + // Match the legacy heuristic's numerical-zero and output-site admission. + if (Z_equal(hybrid,0) || !isValidOutputSite(i1,j1,i2,j2)) return; + const E_type boundaryEnergy=energy.getE(i1,j1,i2,j2,E_type(0)); + const Z_type boundary=energy.getBoltzmannWeight(boundaryEnergy); + if (boundary==0 && !E_isINF(boundaryEnergy)) + throw std::range_error("heuristic pair probabilities: boundary weight underflow"); + const Z_type weight=PartitionArithmetic::multiply(hybrid,boundary); + probabilityDenominator=PartitionArithmetic::add(probabilityDenominator,weight); + probabilityMass(i1,i2)=PartitionArithmetic::add(probabilityMass(i1,i2),weight); + if (continuation.extraPair) + probabilityMass(i1+1,i2+1)=PartitionArithmetic::add(probabilityMass(i1+1,i2+1),weight); + if (continuation.next1!=RnaSequence::lastPos) + probabilityFlow(continuation.next1,continuation.next2)=PartitionArithmetic::add( + probabilityFlow(continuation.next1,continuation.next2),weight); +} + +void PredictorMfeEns2dHeuristic::commitProbabilityRegion() +{ + if (probabilityDenominator!=Zall) + throw std::logic_error("heuristic updateZ override did not preserve the pair-probability objective"); + const size_t n=hybridZ.size1(),m=hybridZ.size2(); + // Every retained cell represents exactly one chain. Its incoming candidate + // mass can therefore pass unchanged to the selected child, without division. + for(size_t i=0;i original(n,m,0); + for(size_t i=0;icollectsSeedPairs(); + Matrix seeds(annotate?n:0,annotate?m:0,0); + pairProbabilities->addRegion(target,query,Zall,original,annotate?&seeds:nullptr); +} + +//////////////////////////////////////////////////////////////////////////// + } // namespace diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristic.h b/src/IntaRNA/PredictorMfeEns2dHeuristic.h index 8167c1c5..75e5b935 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristic.h +++ b/src/IntaRNA/PredictorMfeEns2dHeuristic.h @@ -18,6 +18,10 @@ namespace IntaRNA { * all possible interaction ranges. * * This yields a quadratic time and space complexity. + * Optional pair probabilities refer to all admitted candidate chains considered + * by this pruning rule, with the same Zall as ordinary heuristic prediction. + * They approximate the unrestricted interaction ensemble. Collection adds + * quadratic storage and a reverse traversal of the retained continuations. * * @author Martin Raden * @author Frank Gelhausen @@ -40,10 +44,13 @@ class PredictorMfeEns2dHeuristic: public PredictorMfeEns2d { * @param predTracker the prediction tracker to be used or NULL if no * tracking is to be done; if non-NULL, the tracker gets deleted * on this->destruction. + * @param pairProbabilities optional non-owning sink for probabilities of the + * pruned candidate ensemble; requires needZall */ PredictorMfeEns2dHeuristic( const InteractionEnergy & energy , OutputHandler & output - , PredictionTracker * predTracker ); + , PredictionTracker * predTracker + , BasePairProbabilities * pairProbabilities = nullptr ); virtual ~PredictorMfeEns2dHeuristic(); @@ -75,8 +82,35 @@ class PredictorMfeEns2dHeuristic: public PredictorMfeEns2d { //! energy of all interaction hybrids starting in i1,i2 Z2dMatrix hybridZ; + //! The selected chain owns its first pair and, optionally, its next stack pair. + struct ProbabilityContinuation { + size_t next1 = RnaSequence::lastPos, next2 = RnaSequence::lastPos; + bool extraPair = false; + }; + //! Allocated only for probability output; selected continuations form a DAG. + Matrix probabilityContinuation; + //! Complete candidate weight flowing into each retained continuation. + Matrix probabilityFlow; + //! Numerators in local target/reversed-query coordinates until region commit. + Matrix probabilityMass; + //! Independently accumulated objective, required to agree exactly with Zall. + Z_type probabilityDenominator = 0; + protected: + /** Compute one region under predict()'s failure guard. */ + void predictRegionHeuristic(const IndexRange & r1,const IndexRange & r2); + + /** Record precisely the accepted updateZ candidate, crediting newly prepended + * pairs and routing its full weight to the retained child chain. Coordinates + * are local and the query is reversed; lastPos denotes an initial candidate. + */ + void recordProbabilityCandidate(size_t i1,size_t j1,size_t i2,size_t j2, + Z_type hybrid,const ProbabilityContinuation & continuation); + + /** Propagate candidate weights through selected chains and commit this region. */ + void commitProbabilityRegion(); + /** * Computes all entries of the hybridE matrix * and reports all valid interactions via updateOptima() diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp index 4cf8aeaa..278eb740 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.cpp @@ -9,14 +9,16 @@ PredictorMfeEns2dHeuristicSeedExtension( const InteractionEnergy & energy , OutputHandler & output , PredictionTracker * predTracker - , SeedHandler * seedHandlerInstance ) + , SeedHandler * seedHandlerInstance + , BasePairProbabilities * pairProbabilities ) : - PredictorMfeEns2dSeedExtension(energy,output,predTracker,seedHandlerInstance) + PredictorMfeEns2dSeedExtension(energy,output,predTracker,seedHandlerInstance,pairProbabilities) , E_right_opt(E_INF) , j1opt(0) , j2opt(0) { - assert( seedHandler.getConstraint().getBasePairs() > 1 ); + if(seedHandler.getConstraint().getBasePairs()<2 && seedHandler.getConstraint().getExplicitSeeds().empty()) + throw std::invalid_argument("heuristic seed extension requires at least two seed pairs"); } ////////////////////////////////////////////////////////////////////////// @@ -32,6 +34,12 @@ void PredictorMfeEns2dHeuristicSeedExtension:: predict( const IndexRange & r1, const IndexRange & r2 ) { + try { + if(pairProbabilities) pairProbabilities->markApproximate(); + if(!r1.isAscending() || !r2.isAscending() + || r1.from>=energy.getAccessibility1().getSequence().size() + || r2.from>=energy.getAccessibility2().getAccessibilityOrigin().getSequence().size()) + throw std::invalid_argument("heuristic seeded prediction: invalid region"); #if INTARNA_MULITHREADING #pragma omp critical(intarna_omp_logOutput) #endif @@ -52,16 +60,21 @@ predict( const IndexRange & r1, const IndexRange & r2 ) seedHandler.setOffset2(r2.from); const size_t range_size1 = std::min( energy.size1() - , (r1.to==RnaSequence::lastPos?energy.size1()-1:r1.to)-r1.from+1 ); + , r1.to==RnaSequence::lastPos?energy.size1():r1.to-r1.from+1 ); const size_t range_size2 = std::min( energy.size2() - , (r2.to==RnaSequence::lastPos?energy.size2()-1:r2.to)-r2.from+1 ); + , r2.to==RnaSequence::lastPos?energy.size2():r2.to-r2.from+1 ); + + initExtensionProbabilities(range_size1,range_size2); // compute seed interactions for whole range // and check if any seed possible - if (seedHandler.fillSeed( 0, range_size1-1, 0, range_size2-1 ) == 0) { + if ((seedHandler.getConstraint().getExplicitSeeds().empty() + && (range_size1 range_size1 || sj2 > range_size2) + if (sj1 >= range_size1 || sj2 >= range_size2 + || sl1>energy.getAccessibility1().getMaxLength() || sl2>energy.getAccessibility2().getMaxLength()) continue; const size_t maxMatrixLen1 = energy.getAccessibility1().getMaxLength()-sl1+1; const size_t maxMatrixLen2 = energy.getAccessibility2().getMaxLength()-sl2+1; @@ -93,6 +108,8 @@ predict( const IndexRange & r1, const IndexRange & r2 ) j2opt = sj2; E_right_opt = E_INF; + beginExtensionProbabilitySeed(si1,si2); + // ER hybridZ_right.resize( std::min(range_size1-sj1, maxMatrixLen1), std::min(range_size2-sj2, maxMatrixLen2) ); fillHybridZ_right(sj1, sj2, si1, si2); @@ -104,10 +121,13 @@ predict( const IndexRange & r1, const IndexRange & r2 ) fillHybridZ_left(si1, si2); } + finishExtensionProbabilitySeed(); } // si1 / si2 // report mfe interaction + commitExtensionProbabilities(); reportOptima(); + } catch (...) { if(pairProbabilities) pairProbabilities->fail(); throw; } } @@ -141,6 +161,7 @@ fillHybridZ_right( const size_t sj1, const size_t sj2 // update overall partition function information for true right-extensions of the current seed // seed only not covered due to min-val of r1,r2 updateZ(si1, sj1+r1, si2, sj2+r2, seedZ * rightExtZ * initZ, true); + addExtensionProbabilityRoot(si1,sj1+r1,si2,sj2+r2,seedZ*rightExtZ*initZ,false,true); } } @@ -175,6 +196,7 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) // Z( left + seed ); covers seed only updateZ(si1-l1, sj1, si2-l2, sj2, curZ * seedZ, true); + addExtensionProbabilityRoot(si1-l1,sj1,si2-l2,sj2,curZ*seedZ,true,false); // Z( left + seed + rightOpt ) and rightOpt true seed extension if ( l1 > 0 // true left seed extension @@ -186,6 +208,7 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) ) { updateZ(si1-l1, j1opt, si2-l2, j2opt, curZ * seedZ * rightOptZ, true); + addExtensionProbabilityRoot(si1-l1,j1opt,si2-l2,j2opt,curZ*seedZ*rightOptZ,true,true); } } } // i2 diff --git a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.h b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.h index ff8deada..a0ea94c7 100644 --- a/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.h +++ b/src/IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.h @@ -32,13 +32,15 @@ class PredictorMfeEns2dHeuristicSeedExtension: public PredictorMfeEns2dSeedExten * @param predTracker the prediction tracker to be used or NULL if no * tracking is to be done; if non-NULL, the tracker gets deleted * on this->destruction. + * @param pairProbabilities optional non-owning sink for this approximate ensemble * @param seedHandler the seed handler to be used */ PredictorMfeEns2dHeuristicSeedExtension( const InteractionEnergy & energy , OutputHandler & output , PredictionTracker * predTracker - , SeedHandler * seedHandler ); + , SeedHandler * seedHandler + , BasePairProbabilities * pairProbabilities = nullptr ); /** diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp index fb7c8929..7fd217f0 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.cpp @@ -2,9 +2,188 @@ #include "IntaRNA/PredictorMfeEns2dSeedExtension.h" #include "IntaRNA/SeededPartitionFunction.h" #include "IntaRNA/PartitionArithmetic.h" +#include +#include namespace IntaRNA { +// The heuristic objective contains signed seed-overlap corrections. Trace the +// implemented arithmetic directly, retaining pair ownership on every term. +// This is allocated only for requested probability output; ordinary prediction +// follows its original arithmetic and storage path. +struct PredictorMfeEns2dSeedExtension::ExtensionProbabilityTrace { + using Pair=Interaction::BasePair; + struct Edge { size_t child; Z_type weight; std::vector pairs; }; + struct Tape { + size_t columns=0; + std::vector> edges; + std::vector adjoint; + void reset(size_t n,size_t m); + void add(size_t p1,size_t p2,size_t c1,size_t c2,Z_type weight,std::vector pairs); + void reverse(const Matrix & values,ExtensionProbabilityTrace & trace); + }; + Matrix mass,absoluteMass; + Matrix seeds; + Tape left,right; + std::vector anchor; + size_t si1=0,si2=0,sj1=0,sj2=0; + Z_type seedZ=0,z=0; + explicit ExtensionProbabilityTrace(size_t n,size_t m); + static Z_type finite(Z_type x); + static Z_type product(Z_type a,Z_type b); + void addMass(const Pair & bp,Z_type value); +}; + +void +PredictorMfeEns2dSeedExtension::ExtensionProbabilityTrace::Tape::reset(size_t n,size_t m) +{ + columns=m; + edges.assign(n*m,{}); + adjoint.assign(n*m,0); +} + +void +PredictorMfeEns2dSeedExtension::ExtensionProbabilityTrace::Tape::add( + size_t p1,size_t p2,size_t c1,size_t c2,Z_type weight,std::vector pairs) +{ + finite(weight); + edges.at(p1*columns+p2).push_back({c1*columns+c2,weight,std::move(pairs)}); +} + +void +PredictorMfeEns2dSeedExtension::ExtensionProbabilityTrace::Tape::reverse( + const Matrix & values,ExtensionProbabilityTrace & trace) +{ + for(size_t parent=edges.size();parent-- >0;) { + finite(values(parent/columns,parent%columns)); + for(const auto & edge:edges[parent]) { + const Z_type a=product(adjoint[parent],edge.weight); + adjoint[edge.child]=finite(adjoint[edge.child]+a); + const Z_type mass=product(a,values(edge.child/columns,edge.child%columns)); + for(const auto & bp:edge.pairs) trace.addMass(bp,mass); + } + } +} + +PredictorMfeEns2dSeedExtension::ExtensionProbabilityTrace::ExtensionProbabilityTrace(size_t n,size_t m) + : mass(n,m,0),absoluteMass(n,m,0),seeds(n,m,0) +{} + +Z_type +PredictorMfeEns2dSeedExtension::ExtensionProbabilityTrace::finite(Z_type x) +{ + if (!std::isfinite(x)) throw std::range_error("heuristic seeded probabilities: nonfinite arithmetic"); + return x; +} + +Z_type +PredictorMfeEns2dSeedExtension::ExtensionProbabilityTrace::product(Z_type a,Z_type b) +{ + finite(a); finite(b); + const Z_type value=finite(a*b); + if(a!=0 && b!=0 && (value==0 || std::abs(value)::min())) + throw std::range_error("heuristic seeded probabilities: product underflow"); + return value; +} + +void +PredictorMfeEns2dSeedExtension::ExtensionProbabilityTrace::addMass(const Pair & bp,Z_type value) +{ + mass(bp.first,bp.second)=finite(mass(bp.first,bp.second)+value); + absoluteMass(bp.first,bp.second)=finite(absoluteMass(bp.first,bp.second)+std::abs(value)); +} + +void PredictorMfeEns2dSeedExtension::initExtensionProbabilities(size_t n,size_t m) +{ + if(pairProbabilities) extensionProbabilityTrace=std::make_unique(n,m); +} + +std::vector +PredictorMfeEns2dSeedExtension::extensionSeedPrefix(size_t i1,size_t i2,size_t before) const +{ + Interaction seed(energy.getAccessibility1().getSequence(),energy.getAccessibility2().getAccessibilityOrigin().getSequence()); + seed.basePairs.push_back(energy.getBasePair(i1,i2)); + seedHandler.traceBackSeed(seed,i1,i2); + seed.basePairs.push_back(energy.getBasePair(i1+seedHandler.getSeedLength1(i1,i2)-1, + i2+seedHandler.getSeedLength2(i1,i2)-1)); + std::vector pairs; + for(const auto & bp:seed.basePairs) { + const size_t p=energy.getIndex1(bp),q=energy.getIndex2(bp); + if(p::epsilon()*Z_type(n+m+1); + if(std::abs(Zall-trace.z)>tolerance*std::max(std::abs(trace.z),std::abs(Zall))) + throw std::logic_error("heuristic seeded probabilities: outside objective differs from reported partition"); + const auto first=energy.getBasePair(0,0),last=energy.getBasePair(n-1,m-1); + const IndexRange target(first.first,last.first),query(last.second,first.second); + Matrix original(n,m,0); + Matrix annotations(n,m,0); + for(size_t i=0;iaddRegion(target,query,trace.z,original,&annotations); + extensionProbabilityTrace.reset(); +} + ////////////////////////////////////////////////////////////////////////// PredictorMfeEns2dSeedExtension:: @@ -21,8 +200,8 @@ PredictorMfeEns2dSeedExtension( , hybridZ_left( 0,0 ) , hybridZ_right( 0,0 ) { - if (pairProbabilities && (!seedHandler.guaranteesStackOnlySeeds() || !output.getOutputConstraint().needZall)) - throw std::invalid_argument("pair probabilities require a stack-only exact predictor and needZall"); + if (pairProbabilities && !output.getOutputConstraint().needZall) + throw std::invalid_argument("pair probabilities require needZall"); } ////////////////////////////////////////////////////////////////////////// @@ -43,6 +222,10 @@ predict( const IndexRange & r1, const IndexRange & r2 ) catch (...) { if (pairProbabilities) pairProbabilities->fail(); throw; } return; } + if (pairProbabilities) { + pairProbabilities->fail(); + throw std::invalid_argument("exact bulged-seed pair probabilities are unsupported"); + } if (seedHandler.getConstraint().getBasePairs()<2) throw std::invalid_argument("legacy bulged seed extension requires at least two seed pairs"); #if INTARNA_MULITHREADING @@ -252,6 +435,7 @@ void PredictorMfeEns2dSeedExtension:: fillHybridZ_left( const size_t si1, const size_t si2 ) { + if(extensionProbabilityTrace) extensionProbabilityTrace->left.reset(hybridZ_left.size1(),hybridZ_left.size2()); // temporary access const OutputConstraint & outConstraint = output.getOutputConstraint(); #if INTARNA_IN_DEBUG_MODE @@ -297,6 +481,7 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) iStackZ = energy.getBoltzmannWeight(energy.getE_interLeft(i1,i1+noLpShift,i2,i2+noLpShift)); // check just stacked curZ += iStackZ * hybridZ_left(l1-noLpShift,l2-noLpShift); + if(extensionProbabilityTrace) extensionProbabilityTrace->left.add(l1,l2,l1-noLpShift,l2-noLpShift,iStackZ,{{i1,i2}}); } // check all combinations of decompositions into (i1,i2)..(k1,k2)-(j1,j2) @@ -311,6 +496,12 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) curZ += (iStackZ * energy.getBoltzmannWeight(energy.getE_interLeft(i1+noLpShift,k1,i2+noLpShift,k2)) * hybridZ_left(si1-k1,si2-k2)); + if(extensionProbabilityTrace) { + std::vector pairs{{i1,i2}}; + if(noLpShift) pairs.emplace_back(i1+1,i2+1); + extensionProbabilityTrace->left.add(l1,l2,si1-k1,si2-k2, + iStackZ*energy.getBoltzmannWeight(energy.getE_interLeft(i1+noLpShift,k1,i2+noLpShift,k2)),std::move(pairs)); + } } } // k2 } // k1 @@ -356,9 +547,12 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) Z_type correctionTerm = energy.getBoltzmannWeight( nonOverlapE ) * hybridZ_left( si1-si1overlap, si2-si2overlap); curZ -= correctionTerm; + if(extensionProbabilityTrace) extensionProbabilityTrace->left.add(l1,l2,si1-si1overlap,si2-si2overlap, + -energy.getBoltzmannWeight(nonOverlapE),extensionSeedPrefix(i1,i2,si1overlap)); // sanity insurance if (curZ < 0) { curZ = Z_type(0.0); + if(extensionProbabilityTrace) extensionProbabilityTrace->left.edges[l1*extensionProbabilityTrace->left.columns+l2].clear(); } } } else { @@ -371,6 +565,8 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) Z_type correctionTerm = seedZ_rm * hybridZ_left( si1-sj1_rm , si2-sj2_rm ); + if(extensionProbabilityTrace) extensionProbabilityTrace->left.add(l1,l2,si1-sj1_rm,si2-sj2_rm, + -seedZ_rm,extensionSeedPrefix(i1,i2,sj1_rm)); // if noLP : handle explicit loop right of current seed if (outConstraint.noLP) { for (k1=sj1_rm; k1++ < si1; ) { @@ -386,6 +582,8 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) correctionTerm += (seedZ_rm * energy.getBoltzmannWeight(energy.getE_interLeft(sj1_rm,k1,sj2_rm,k2)) * hybridZ_left(si1-k1,si2-k2) ); + if(extensionProbabilityTrace) extensionProbabilityTrace->left.add(l1,l2,si1-k1,si2-k2, + -seedZ_rm*energy.getBoltzmannWeight(energy.getE_interLeft(sj1_rm,k1,sj2_rm,k2)),extensionSeedPrefix(i1,i2,k1)); } } // k2 } // k1 @@ -395,6 +593,7 @@ fillHybridZ_left( const size_t si1, const size_t si2 ) // sanity insurance if (curZ < 0) { curZ = Z_type(0.0); + if(extensionProbabilityTrace) extensionProbabilityTrace->left.edges[l1*extensionProbabilityTrace->left.columns+l2].clear(); } } } // substractThisSeed @@ -413,6 +612,7 @@ void PredictorMfeEns2dSeedExtension:: fillHybridZ_right( const size_t sj1, const size_t sj2 ) { + if(extensionProbabilityTrace) extensionProbabilityTrace->right.reset(hybridZ_right.size1(),hybridZ_right.size2()); // temporary access const OutputConstraint & outConstraint = output.getOutputConstraint(); #if INTARNA_IN_DEBUG_MODE @@ -457,6 +657,7 @@ fillHybridZ_right( const size_t sj1, const size_t sj2 ) // check just stacked seed extension if (j1-noLpShift==sj1 && j2-noLpShift==sj2) { curZ += iStackZ * hybridZ_right(0,0); + if(extensionProbabilityTrace) extensionProbabilityTrace->right.add(j1-sj1,j2-sj2,0,0,iStackZ,{{j1,j2}}); } } @@ -473,6 +674,12 @@ fillHybridZ_right( const size_t sj1, const size_t sj2 ) curZ += ( hybridZ_right(k1-sj1,k2-sj2) * energy.getBoltzmannWeight(energy.getE_interLeft(k1,j1-noLpShift,k2,j2-noLpShift)) * iStackZ ); + if(extensionProbabilityTrace) { + std::vector pairs{{j1,j2}}; + if(noLpShift) pairs.emplace_back(j1-1,j2-1); + extensionProbabilityTrace->right.add(j1-sj1,j2-sj2,k1-sj1,k2-sj2, + energy.getBoltzmannWeight(energy.getE_interLeft(k1,j1-noLpShift,k2,j2-noLpShift))*iStackZ,std::move(pairs)); + } } } // k2 } // k1 diff --git a/src/IntaRNA/PredictorMfeEns2dSeedExtension.h b/src/IntaRNA/PredictorMfeEns2dSeedExtension.h index e046de5b..0b344e4a 100644 --- a/src/IntaRNA/PredictorMfeEns2dSeedExtension.h +++ b/src/IntaRNA/PredictorMfeEns2dSeedExtension.h @@ -6,6 +6,7 @@ #include "IntaRNA/BasePairProbabilities.h" #include "IntaRNA/Matrix.h" #include "IntaRNA/SeedHandlerIdxOffset.h" +#include namespace IntaRNA { @@ -85,6 +86,27 @@ class PredictorMfeEns2dSeedExtension: public PredictorMfeEns { /** Run the selected disjoint stack-seed partition backend. */ void predictStackSeeds(const IndexRange & r1,const IndexRange & r2); + //! Optional signed reverse trace for the legacy anchored extension recurrences. + struct ExtensionProbabilityTrace; + std::unique_ptr extensionProbabilityTrace; + /** Begin a regional trace; no allocation if no probability sink was requested. */ + void initExtensionProbabilities(size_t n,size_t m); + /** Store the current anchor seed's actual pairs, including its endpoints. */ + void beginExtensionProbabilitySeed(size_t i1,size_t i2); + /** Add one accepted left/seed/right objective to the optional reverse trace. + * All coordinates are local, with seed boundaries in increasing DP order. + * @param left whether the left matrix is a factor (otherwise use initiation) + * @param right whether the right matrix is a factor + */ + void addExtensionProbabilityRoot(size_t i1,size_t j1,size_t i2,size_t j2, + Z_type partZ,bool left,bool right); + /** Reverse both extension traces and accumulate actual pair numerators. */ + void finishExtensionProbabilitySeed(); + /** Validate and atomically commit this region in original sequence order. */ + void commitExtensionProbabilities(); + /** Actual seed pairs strictly before the given local target coordinate. */ + std::vector extensionSeedPrefix(size_t i1,size_t i2,size_t before) const; + //! optional sequence-pair owner, never deleted or finalized by the predictor BasePairProbabilities * pairProbabilities; diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index d5a0f30e..623235d3 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -958,7 +958,7 @@ CommandLineParsing::CommandLineParsing( const Personality personality ) "\n 'tPu:' (target) unpaired probability values (RNAplfold format OR '.agz' file ending for IntaRNA's binary format)." "\n 'pMinE:' (target+query) for each index pair the minimal energy of any interaction covering the pair (CSV format)" "\n 'bpsvg:' (target+query) SVG actual-pair probability plot with accessibility frame and seed annotations; same requirements as bpProb" - "\n 'bpProb:' (target+query) actual-pair conditional probabilities for --model=P --mode=M with --noSeed or stack-only seeds (CSV matrix, original strand order; NA for an empty ensemble)" + "\n 'bpProb:' (target+query) actual-pair conditional probabilities for --model=P --mode=M/H (M: --noSeed or stack-only seeds; H: approximate retained ensemble) (CSV matrix, original strand order; NA for an empty ensemble)" "\n 'spotProb:' (target+query) tracks for a given set of interaction spots their probability to be covered by an interaction. If no spots are provided, probabilities for all index combinations are computed. Spots are encoded by comma-separated 'idxT&idxQ' pairs (target-query). For each spot a probability is provided in concert with the probability that none of the spots (encoded by '0&0') is covered (CSV format). The spot encoding is followed colon-separated by the output stream/file name, eg. '--out=\"spotProb:3&76,59&2:STDERR\"'. NOTE: value has to be quoted due to '&' symbol!" "\nFor each, provide a file name or STDOUT/STDERR to write to the respective output stream." ).c_str()) @@ -1246,8 +1246,8 @@ parse(int argc, char** argv) outSep = unescaped_string::getUnescaped( outSep ); if (hasBasePairProbabilityOutput()) { - if (model.val!='P' || mode.val!='M' || !rri.empty() || windowWidth.val!=0) - throw error("bpProb/bpsvg requires --model=P --mode=M with --noSeed or stack-only seeds and no windows"); + if (model.val!='P' || (mode.val!='M' && mode.val!='H') || !rri.empty() || windowWidth.val!=0) + throw error("bpProb/bpsvg requires --model=P --mode=M/H with no windows (M: --noSeed or stack-only seeds)"); outNeedsZall=true; } // K needs a seed even before the usual --noSeed model normalization. @@ -2498,7 +2498,7 @@ BasePairProbabilities * CommandLineParsing::getBasePairProbabilityResult(const I if (!hasBasePairProbabilityOutput()) return nullptr; if (!noSeedRequired) { std::unique_ptr handler(getSeedHandler(energy)); - if (!handler->guaranteesStackOnlySeeds()) + if (mode.val=='M' && !handler->guaranteesStackOnlySeeds()) throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); } return new BasePairProbabilities(energy.size1(),energy.size2(),!noSeedRequired && !outPrefix2streamName.at(OP_bpsvg).empty()); @@ -2520,11 +2520,11 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output, BasePair { std::unique_ptr pairSeeds; if (pairProbabilities) { - if(model.val!='P' || mode.val!='M' || !rri.empty() || windowWidth.val!=0) - throw std::invalid_argument("bpProb/bpsvg requires --model=P --mode=M with --noSeed or stack-only seeds and no windows"); + if(model.val!='P' || (mode.val!='M' && mode.val!='H') || !rri.empty() || windowWidth.val!=0) + throw std::invalid_argument("bpProb/bpsvg requires --model=P --mode=M/H with no windows (M: --noSeed or stack-only seeds)"); if (!noSeedRequired) { pairSeeds.reset(getSeedHandler(energy)); - if(!pairSeeds->guaranteesStackOnlySeeds()) + if(mode.val=='M' && !pairSeeds->guaranteesStackOnlySeeds()) throw std::invalid_argument("bpProb/bpsvg requires a stack-only seed handler; bulged seeds are unsupported"); } } @@ -2627,7 +2627,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output, BasePair // single-site mfe ensemble interactions (contain only interior loops) case 'P' : { switch ( mode.val ) { - case 'H' : return new PredictorMfeEns2dHeuristic( energy, output, predTracker ); + case 'H' : return new PredictorMfeEns2dHeuristic( energy, output, predTracker, pairProbabilities ); case 'M' : return new PredictorMfeEns2d( energy, output, predTracker, pairProbabilities ); default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); } @@ -2672,7 +2672,7 @@ getPredictor( const InteractionEnergy & energy, OutputHandler & output, BasePair // single-site max-prob interactions (contain only interior loops) case 'P' : { switch ( mode.val ) { - case 'H' : return new PredictorMfeEns2dHeuristicSeedExtension( energy, output, predTracker, getSeedHandler( energy ) ); + case 'H' : return new PredictorMfeEns2dHeuristicSeedExtension( energy, output, predTracker, pairSeeds ? pairSeeds.release() : getSeedHandler(energy), pairProbabilities ); case 'M' : return new PredictorMfeEns2dSeedExtension( energy, output, predTracker, pairSeeds ? pairSeeds.release() : getSeedHandler( energy ), pairProbabilities ); case 'S' : return new PredictorMfeEnsSeedOnly( energy, output, predTracker, getSeedHandler(energy) ); default : INTARNA_NOT_IMPLEMENTED("mode "+toString(mode.val)+" not available for model "+toString(model.val)); diff --git a/tests/Makefile.am b/tests/Makefile.am index a500da8e..7ff2a868 100644 --- a/tests/Makefile.am +++ b/tests/Makefile.am @@ -31,6 +31,8 @@ runApiTests_SOURCES = \ StackSeedDomain_test.cpp \ BasePairProbabilities_test.cpp \ PredictorUnseededProbability_test.cpp \ + PredictorHeuristicProbability_test.cpp \ + PredictorHeuristicSeededProbability_test.cpp \ testEasyLoggingSetup.icc \ AccessibilityConstraint_test.cpp \ AccessibilityFromStream_test.cpp \ diff --git a/tests/PredictorHeuristicProbability_test.cpp b/tests/PredictorHeuristicProbability_test.cpp new file mode 100644 index 00000000..96fdf6b8 --- /dev/null +++ b/tests/PredictorHeuristicProbability_test.cpp @@ -0,0 +1,234 @@ +#include "catch.hpp" + +#undef NDEBUG + +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/AccessibilityBasePair.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/InteractionEnergyIdxOffset.h" +#include "IntaRNA/OutputHandlerInteractionList.h" +#include "IntaRNA/PredictorMfeEns2dHeuristic.h" +#include "IntaRNA/ReverseAccessibility.h" + +#include +#include +#include +#include +#include +#include +#include + +using namespace IntaRNA; + +namespace { + +using Pair=std::array; +using Chain=std::vector; + +struct ChainCandidate { + Chain chain; + E_type hybrid; +}; + +struct HeuristicReference { + Z_type z=0; + std::map mass; + size_t candidates=0; +}; + +bool adjacent(Pair p,Pair q) { + return p[0]+1==q[0] && p[1]+1==q[1]; +} + +/** Enumerate complete chains first, then filter by the documented selected + * suffix rule. Unlike the implementation, this stores complete pair lists, + * sums integer chain energies and directly counts pairs in admitted chains. + */ +HeuristicReference enumerateHeuristic(const InteractionEnergy & energy, + const OutputConstraint & constraint,size_t n,size_t m) +{ + std::map> chains; + Chain chain; + std::function enumerate=[&](E_type hybrid) { + const Pair left=chain.front(),right=chain.back(); + if (right[0]-left[0]+1>energy.getAccessibility1().getMaxLength() + || right[1]-left[1]+1>energy.getAccessibility2().getMaxLength()) return; + chains[left].push_back({chain,hybrid}); + for(size_t i=right[0]+1;i selected; + for(size_t i=n;i-- >0;) for(size_t j=m;j-- >0;) { + const Pair left{i,j}; + std::vector,const ChainCandidate *>> admitted; + for(const auto & candidate:chains[left]) { + const auto & c=candidate.chain; + if (constraint.noGUend && energy.isGU(c.back()[0],c.back()[1])) continue; + size_t suffix=1; + std::array order{0,0,0}; + if (constraint.noLP) { + if (c.size()<2 || !adjacent(c[0],c[1])) continue; + if (c.size()==2) suffix=c.size(); + else if (selected[c[1]]==Chain(c.begin()+1,c.end())) order={1,0,0}; + else { + if (adjacent(c[1],c[2])) continue; + suffix=2;order={2,c[2][0],c[2][1]}; + } + } else if (c.size()>1) order={1,c[1][0],c[1][1]}; + if (suffixconstraint.maxED || energy.getED2(j,right[1])>constraint.maxED)) continue; + const Z_type weight=energy.getBoltzmannWeight(full); + result.z+=weight;++result.candidates; + for(const Pair pair:candidate.chain) result.mass[pair]+=weight; + } + } + return result; +} + +class FailingHeuristic : public PredictorMfeEns2dHeuristic { +public: + using PredictorMfeEns2dHeuristic::PredictorMfeEns2dHeuristic; + bool reject=false; +protected: + void updateZ(size_t i1,size_t j1,size_t i2,size_t j2,Z_type z,bool hybrid) override { + if (reject) throw std::runtime_error("test probability failure"); + PredictorMfeEns::updateZ(i1,j1,i2,j2,z,hybrid); + } +}; + +class HeuristicTinyInitiation : public InteractionEnergyBasePair { +public: + using InteractionEnergyBasePair::InteractionEnergyBasePair; + E_type getE_init() const override { return Ekcal_2_E(60); } +}; + +} + +TEST_CASE("heuristic pair probabilities count the retained candidate chains", "[BasePairProbabilities][PredictorHeuristicProbability]") { + #include "testEasyLoggingSetup.icc" + for(const auto & sequences:std::vector>{ + {"GGGGG","CCCCC"},{"GGAGG","CCUCC"},{"GUGGC","GCCUU"}}) { + RnaSequence target("target",sequences.first),query("query",sequences.second); + AccessibilityBasePair at(target,4,nullptr),aq(query,4,nullptr);ReverseAccessibility reverse(aq); + InteractionEnergyBasePair energy(at,reverse,2,1); + for(bool noLP:{false,true}) for(bool noGU:{false,true}) for(unsigned regionKind:{0u,1u,2u}) for(E_type maxED:{Accessibility::ED_UPPER_BOUND,E_type(50)}) { + const IndexRange region=regionKind==0?IndexRange(0,4):regionKind==1?IndexRange(1,3):IndexRange(2,RnaSequence::lastPos); + OutputConstraint constraint(1,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP,noGU,true,true,maxED); + InteractionEnergyIdxOffset shifted(energy,region.from,region.from); + const size_t length=region.to==RnaSequence::lastPos?target.size()-region.from:region.to-region.from+1; + const auto reference=enumerateHeuristic(shifted,constraint,length,length); + OutputHandlerInteractionList enabled(constraint,1),disabled(constraint,1); + BasePairProbabilities result(target.size(),query.size(),true); + PredictorMfeEns2dHeuristic on(energy,enabled,nullptr,&result),off(energy,disabled,nullptr); + on.predict(region,region);off.predict(region,region);result.finalize(); + REQUIRE(result.isApproximate()); + CAPTURE(sequences.first,sequences.second,noLP,noGU,regionKind,maxED); + REQUIRE(on.getZall()==off.getZall()); + REQUIRE(result.getZ()==on.getZall()); + REQUIRE(result.getZ()==Approx(reference.z).epsilon(3e-12)); + Matrix expected(target.size(),query.size(),0); + for(const auto & entry:reference.mass) { + const auto bp=shifted.getBasePair(entry.first[0],entry.first[1]); + expected(bp.first,bp.second)=entry.second; + } + for(size_t i=0;ienergy==(*b)->energy); + REQUIRE((*a)->basePairs==(*b)->basePairs); + } + } + } +} + +TEST_CASE("heuristic pair probabilities handle empty and singleton domains", "[BasePairProbabilities][PredictorHeuristicProbability]") { + #include "testEasyLoggingSetup.icc" + for(const std::string & queryText:{std::string("C"),std::string("A")}) for(bool noLP:{false,true}) { + RnaSequence target("target","G"),query("query",queryText); + AccessibilityDisabled at(target,0,nullptr),aq(query,0,nullptr);ReverseAccessibility reverse(aq); + InteractionEnergyBasePair energy(at,reverse); + OutputConstraint constraint(1,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP,false,true,false); + OutputHandlerInteractionList output(constraint,1); + BasePairProbabilities result(1,1); + PredictorMfeEns2dHeuristic predictor(energy,output,nullptr,&result); + predictor.predict();result.finalize(); + REQUIRE(result.isApproximate()); + const bool empty=noLP || queryText=="A"; + REQUIRE((result.status()==BasePairProbabilities::Status::empty)==empty); + if (!empty) REQUIRE(result.probabilities()(0,0)==1); + } +} + +TEST_CASE("heuristic pair owner merges disjoint regions and rejects failed results", "[BasePairProbabilities][PredictorHeuristicProbability]") { + #include "testEasyLoggingSetup.icc" + RnaSequence target("target","GGGGGG"),query("query","CCCCCC"); + AccessibilityDisabled at(target,0,nullptr),aq(query,0,nullptr);ReverseAccessibility reverse(aq); + InteractionEnergyBasePair energy(at,reverse,2,2); + OutputConstraint constraint(1,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,true,false,true,false); + OutputHandlerInteractionList output(constraint,1); + BasePairProbabilities result(6,6); + FailingHeuristic predictor(energy,output,nullptr,&result); + predictor.predict({0,2},{0,2});const Z_type first=predictor.getZall(); + predictor.predict({3,5},{3,5});result.finalize(); + REQUIRE(result.getZ()==2*first); + REQUIRE(result.probabilities()(0,0)==0); + BasePairProbabilities failed(6,6);OutputHandlerInteractionList failedOutput(constraint,1); + FailingHeuristic failing(energy,failedOutput,nullptr,&failed); + failing.predict({0,2},{0,2});failing.reject=true; + REQUIRE_THROWS(failing.predict({3,5},{3,5})); + REQUIRE(failed.status()==BasePairProbabilities::Status::failed); + REQUIRE_THROWS(failed.finalize()); + BasePairProbabilities invalid(6,6);OutputHandlerInteractionList invalidOutput(constraint,1); + FailingHeuristic badRange(energy,invalidOutput,nullptr,&invalid); + REQUIRE_THROWS(badRange.predict({8,9},{0,2})); + REQUIRE(invalid.status()==BasePairProbabilities::Status::failed); +} + +TEST_CASE("heuristic probability collection preserves numerical-zero admission", "[BasePairProbabilities][PredictorHeuristicProbability]") { + #include "testEasyLoggingSetup.icc" + RnaSequence target("target","GGGGG"),query("query","CCCCC"); + AccessibilityDisabled at(target,0,nullptr),aq(query,0,nullptr);ReverseAccessibility reverse(aq); + HeuristicTinyInitiation tiny(at,reverse,2,2,false,1,Ekcal_2_E(-40)); + InteractionEnergyBasePair vanishing(at,reverse,2,2,false,1,Ekcal_2_E(50)); + for(const InteractionEnergy * energy:{static_cast(&tiny),static_cast(&vanishing)}) for(bool noLP:{false,true}) { + OutputConstraint constraint(1,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP,false,true,false); + OutputHandlerInteractionList baselineOut(constraint,1),probabilityOut(constraint,1); + BasePairProbabilities result(5,5); + PredictorMfeEns2dHeuristic baseline(*energy,baselineOut,nullptr),predictor(*energy,probabilityOut,nullptr,&result); + baseline.predict();predictor.predict();result.finalize(); + REQUIRE(result.getZ()==baseline.getZall()); + REQUIRE(predictor.getZall()==baseline.getZall()); + if (energy==&vanishing) REQUIRE(result.status()==BasePairProbabilities::Status::empty); + } + OutputConstraint constraint(0,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,false,false,true,false); + OutputHandlerInteractionList output(constraint,1); + BasePairProbabilities failed(5,5); + InteractionEnergyBasePair underflow(at,reverse,2,2,false,0.01,Ekcal_2_E(-1),3,Ekcal_2_E(1000)); + PredictorMfeEns2dHeuristic predictor(underflow,output,nullptr,&failed); + REQUIRE_THROWS(predictor.predict()); + REQUIRE(failed.status()==BasePairProbabilities::Status::failed); +} diff --git a/tests/PredictorHeuristicSeededProbability_test.cpp b/tests/PredictorHeuristicSeededProbability_test.cpp new file mode 100644 index 00000000..6a15c116 --- /dev/null +++ b/tests/PredictorHeuristicSeededProbability_test.cpp @@ -0,0 +1,201 @@ +#include "catch.hpp" +#include "IntaRNA/AccessibilityDisabled.h" +#include "IntaRNA/InteractionEnergyBasePair.h" +#include "IntaRNA/OutputHandlerInteractionList.h" +#include "IntaRNA/PredictorMfeEns2dHeuristicSeedExtension.h" +#include "IntaRNA/SeedHandlerExplicit.h" +#include "IntaRNA/SeedHandlerMfe.h" +#include +#include +#include +#include + +using namespace IntaRNA; + +namespace { +// Introduce an independent fugacity x at one actual pair: every loop pays for +// its left pair, and the complete boundary pays for the last pair. Every +// physical chain is multilinear in x, including signed overlap corrections. +// Consequently Z(x)=Z_without+x*M_pair and a finite difference gives M exactly, +// without consulting the reverse tape or its pair-ownership bookkeeping. +class PairFugacityEnergy: public InteractionEnergyBasePair { +public: + using InteractionEnergyBasePair::InteractionEnergyBasePair; + size_t marked1=RnaSequence::lastPos,marked2=RnaSequence::lastPos; + E_type penalty=0; + E_type getE_interLeft(size_t i1,size_t j1,size_t i2,size_t j2) const override { + const E_type e=InteractionEnergyBasePair::getE_interLeft(i1,j1,i2,j2); + return E_isINF(e)?e:e+((i1==marked1 && i2==marked2)?penalty:0); + } + E_type getE(size_t i1,size_t j1,size_t i2,size_t j2,E_type hybrid) const override { + const E_type e=InteractionEnergyBasePair::getE(i1,j1,i2,j2,hybrid); + return E_isINF(e)?e:e+((j1==marked1 && j2==marked2)?penalty:0); + } +}; + +class FixedChoiceSeededHeuristic: public PredictorMfeEns2dHeuristicSeedExtension { +public: + using PredictorMfeEns2dHeuristicSeedExtension::PredictorMfeEns2dHeuristicSeedExtension; + using Key=std::array; + struct Choice { size_t j1,j2; Z_type e; }; + using Choices=std::map; + Choices choices; + const Choices * fixed=nullptr; +protected: + void fillHybridZ_right(size_t sj1,size_t sj2,size_t si1,size_t si2) override { + PredictorMfeEns2dHeuristicSeedExtension::fillHybridZ_right(sj1,sj2,si1,si2); + // Condition on the original heuristic's chosen domain; a fugacity is a + // counting variable and must not change the domain being differentiated. + if(fixed) { + const auto & c=fixed->at({si1,si2}); j1opt=c.j1; j2opt=c.j2; E_right_opt=c.e; + } else choices[{si1,si2}]={j1opt,j2opt,E_right_opt}; + } +}; + +class CorruptSeededHeuristicPartition: public PredictorMfeEns2dHeuristicSeedExtension { +public: + using PredictorMfeEns2dHeuristicSeedExtension::PredictorMfeEns2dHeuristicSeedExtension; + Z_type injected=0; +protected: + void updateZ(size_t i1,size_t j1,size_t i2,size_t j2,Z_type partZ,bool isHybridZ) override; +}; + +void CorruptSeededHeuristicPartition::updateZ(size_t i1,size_t j1,size_t i2,size_t j2, + Z_type partZ,bool isHybridZ) +{ + PredictorMfeEns2dHeuristicSeedExtension::updateZ(i1,j1,i2,j2,partZ,isHybridZ); + Zall=injected; +} + +SeedConstraint explicitSeeds(const std::string & patterns,size_t nominal=2) { + return SeedConstraint(nominal,6,3,3,E_INF,E_INF,E_INF,IndexRangeList(),IndexRangeList(),patterns,true,true,true); +} +OutputConstraint constraints(bool noLP,bool noGU=false) { + return OutputConstraint(0,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF,false,noLP,noGU,true,false); +} +} + +TEST_CASE("seeded heuristic outside agrees with independent pair fugacities", "[BasePairProbabilities][HeuristicSeededProbability]") { + #include "testEasyLoggingSetup.icc" + RnaSequence t("t","GGGGGG"),q("q","CCCCCC"); + AccessibilityDisabled at(t,0,nullptr),aq(q,0,nullptr); ReverseAccessibility ar(aq); + PairFugacityEnergy energy(at,ar,2,2); + // Isolated and competing seeds exercise left/right/both objectives, + // overlapping and nonoverlapping subtraction, and bulged actual pairs. + for(const auto & patterns:{"3||&3||","2||&4||,3||&3||","1||&5||,4||&2||", + "2|.|&3||","1|.||&3|||,3||&3||","2|||&3|||,4||&2||"}) { + const auto sc=explicitSeeds(patterns); + for(bool noLP:{false,true}) { + CAPTURE(patterns,noLP); + const auto oc=constraints(noLP); + energy.penalty=0; + BasePairProbabilities result(t.size(),q.size(),true); + OutputHandlerInteractionList output(oc,0); + FixedChoiceSeededHeuristic predictor(energy,output,nullptr,new SeedHandlerExplicit(energy,sc),&result); + predictor.predict(); result.finalize(); + REQUIRE(result.isApproximate()); + REQUIRE(result.getZ()==Approx(predictor.getZall()).epsilon(2e-12)); + const Z_type baseline=result.getZ(); + for(size_t i=0;i0); + predictor.predict({3,RnaSequence::lastPos},{3,RnaSequence::lastPos});result.finalize(); + REQUIRE(result.getZ()==Approx(2*first).epsilon(2e-12)); + REQUIRE(result.seedPairs()(0,5)==1);REQUIRE(result.seedPairs()(1,4)==1); + REQUIRE(result.seedPairs()(3,2)==1);REQUIRE(result.seedPairs()(4,1)==1); + REQUIRE(result.probabilities()(0,0)==0); + for(const std::string pattern:{"2||&4||","1|||&4|||"}) { + const auto shortSc=explicitSeeds(pattern); + BasePairProbabilities empty(6,6,true);OutputHandlerInteractionList emptyOutput(oc,0); + FixedChoiceSeededHeuristic p(energy,emptyOutput,nullptr,new SeedHandlerExplicit(energy,shortSc),&empty); + p.predict({0,1},{0,1});empty.finalize(); + REQUIRE(empty.status()==BasePairProbabilities::Status::empty); + } + const auto singleton=explicitSeeds("2|&4|",2); + BasePairProbabilities failed(6,6);OutputHandlerInteractionList failedOutput(oc,0); + FixedChoiceSeededHeuristic bad(energy,failedOutput,nullptr,new SeedHandlerExplicit(energy,singleton),&failed); + REQUIRE_THROWS(bad.predict());REQUIRE(failed.status()==BasePairProbabilities::Status::failed); +} + +TEST_CASE("computed bulged seeded heuristic marginals condition on the selected seed family", "[BasePairProbabilities][HeuristicSeededProbability]") { + #include "testEasyLoggingSetup.icc" + RnaSequence t("t","GGAGGCG"),q("q","CCGCACC"); + AccessibilityDisabled at(t,0,nullptr),aq(q,0,nullptr);ReverseAccessibility ar(aq); + PairFugacityEnergy energy(at,ar,2,2); + const SeedConstraint sc(3,2,1,1,E_INF,E_INF,E_INF,IndexRangeList(),IndexRangeList(),"",true,true,true); + // Freeze the concrete MFE seed selected at every start by encoding its + // actual pairs explicitly. Perturbations then cannot select another seed. + SeedHandlerMfe selected(energy,sc);selected.fillSeed(0,t.size()-1,0,q.size()-1); + std::string patterns; + size_t i=RnaSequence::lastPos,j=RnaSequence::lastPos; + while(selected.updateToNextSeed(i,j,0,t.size()-1,0,q.size()-1)) { + const size_t sl1=selected.getSeedLength1(i,j),sl2=selected.getSeedLength2(i,j); + Interaction seed(t,q);seed.basePairs.push_back(energy.getBasePair(i,j)); + selected.traceBackSeed(seed,i,j);seed.basePairs.push_back(energy.getBasePair(i+sl1-1,j+sl2-1)); + const size_t qStart=q.size()-j-sl2; + std::string ts(sl1,'.'),qs(sl2,'.'); + for(const auto & bp:seed.basePairs) { ts[bp.first-i]='|';qs[bp.second-qStart]='|'; } + if(!patterns.empty()) patterns+=','; + patterns+=std::to_string(i+1)+ts+'&'+std::to_string(qStart+1)+qs; + } + REQUIRE_FALSE(patterns.empty()); + const auto fixedSeeds=explicitSeeds(patterns,3); + for(bool noLP:{false,true}) { + const auto oc=constraints(noLP);energy.penalty=0; + BasePairProbabilities result(t.size(),q.size());OutputHandlerInteractionList output(oc,0); + FixedChoiceSeededHeuristic predictor(energy,output,nullptr,new SeedHandlerMfe(energy,sc),&result); + predictor.predict();result.finalize(); + for(size_t a=0;a::infinity(),std::numeric_limits::quiet_NaN()}) { + CAPTURE(injected); + BasePairProbabilities result(t.size(),q.size());OutputHandlerInteractionList output(oc,0); + CorruptSeededHeuristicPartition predictor(energy,output,nullptr,new SeedHandlerExplicit(energy,sc),&result); + predictor.injected=injected; + REQUIRE_THROWS_AS(predictor.predict(),std::range_error); + REQUIRE(result.status()==BasePairProbabilities::Status::failed); + REQUIRE(result.getZ()==0); + REQUIRE_THROWS(result.finalize()); + } +} diff --git a/tests/PredictorUnseededProbability_test.cpp b/tests/PredictorUnseededProbability_test.cpp index 7124d6e5..c45858e4 100644 --- a/tests/PredictorUnseededProbability_test.cpp +++ b/tests/PredictorUnseededProbability_test.cpp @@ -214,4 +214,10 @@ TEST_CASE("unseeded probability collection follows legacy numerical gates and fa PredictorMfeEns2d huge(extreme,output,nullptr,&overflow); REQUIRE_THROWS(huge.predict()); REQUIRE(overflow.status()==BasePairProbabilities::Status::failed); + BasePairProbabilities underflow(5,5); + InteractionEnergyBasePair underflowEnergy(targetAcc,reverse,2,2,false,1, + Ekcal_2_E(-1),3,Ekcal_2_E(1000)); + PredictorMfeEns2d disappearing(underflowEnergy,output,nullptr,&underflow); + REQUIRE_THROWS(disappearing.predict()); + REQUIRE(underflow.status()==BasePairProbabilities::Status::failed); } diff --git a/tests/testBasePairProbabilities.py b/tests/testBasePairProbabilities.py index 09572802..bafed239 100644 --- a/tests/testBasePairProbabilities.py +++ b/tests/testBasePairProbabilities.py @@ -27,7 +27,7 @@ def run(t='GGAGGG', q='CCCC', options=(), success=True): def matrix(path): with (gzip.open(path, 'rt') if str(path).endswith('.gz') else open(path)) as f: rows = list(csv.reader(f, delimiter=';')) - assert rows[0][0] == 'bpProb', rows + assert rows[0][0] in ('bpProb', 'bpProbApproximate'), rows return rows[0][1:], [r[0] for r in rows[1:]], [[None if v=='NA' else float(v) for v in r[1:]] for r in rows[1:]] def oracle(t, q, regions=None, no_lp=False, seeds=None, unseeded=False): @@ -116,8 +116,16 @@ def predict(t='GGAGGG',q='CCCC',extra=()): covered=list(csv.reader(spot.open(),delimiter=';')) assert any(float(v)>0 for v in covered[3][1:]);assert all(v==0 for v in with_spot[2]) compressed=d/'pairs.csv.gz';run(options=[f'--out=bpProb:{compressed}']);same(matrix(compressed)[2],values) + # Heuristic outputs retain their own denominator and declare approximation. + for seed_opts in ([], ['--noSeed'], ['--seedMaxUP=1']): + for no_lp in ('false', 'true'): + opts=['--mode=H', f'--outNoLP={no_lp}', *seed_opts] + predicted, probs=predict(extra=opts) + assert predicted.stdout==run(options=opts).stdout + assert path.read_text().startswith('bpProbApproximate;') + assert all(v is None or 0<=v<=1 for row in probs for v in row) # Unsupported requests fail explicitly and never create a probability file. - for opts in (['--mode=H'],['--mode=S'],['--model=X','--mode=K'],['--windowWidth=3','--windowOverlap=2'],['--seedMaxUP=1'],['--seedTQ=1|.|&2||']): + for opts in (['--mode=S'],['--model=X','--mode=K'],['--windowWidth=3','--windowOverlap=2'],['--seedMaxUP=1'],['--seedTQ=1|.|&2||']): fail=d/'unsupported.csv';run(options=[f'--out=bpProb:{fail}',*opts],success=False);assert not fail.exists() fail=d/'range.csv';run('GGGG','CCCC',[f'--out=bpProb:{fail}',*seedopt,'--energyAdd=900'],success=False);assert not fail.exists() run(options=['--out=bpProb:/dev/full'],success=False) diff --git a/tests/testBasePairProbabilitySvg.py b/tests/testBasePairProbabilitySvg.py index 0cbb64b6..1e7fab9d 100644 --- a/tests/testBasePairProbabilitySvg.py +++ b/tests/testBasePairProbabilitySvg.py @@ -127,6 +127,15 @@ def plot(t='GGAGGG', q='CCCC', opts=()): rows = list(csv.reader(table.open(), delimiter=';')) assert all(float(pairs(unseeded)[i,j].attrib['data-probability']) == float(rows[i][j]) for i in range(1,7) for j in range(1,5)) + for seed_opts in ([], ['--noSeed'], ['--seedMaxUP=1']): + approximate=plot(opts=['--mode=H', f'--out=bpProb:{table}', *seed_opts]) + assert 'Approximate' in approximate.find('s:title', NS).text + rows=list(csv.reader(table.open(), delimiter=';')) + assert rows[0][0]=='bpProbApproximate' + assert all(float(pairs(approximate)[i,j].attrib['data-probability']) == float(rows[i][j]) + for i in range(1,7) for j in range(1,5)) + if '--noSeed' in seed_opts: + assert not seeds(approximate) # Explicit singleton/mixed seed masks exclude extension-only cells. opts = ['--seedTQ=1|&4|,2||&2||'] mixed = plot('GGGG', 'CCCC', opts) @@ -219,7 +228,7 @@ def check_frame(root, target_pu, query_pu): run(options=[f'--out=bpsvg:{compressed}']) assert len(pairs(load(compressed))) == 24 # Validation/numerical failures must not publish a partial document. - for opts in (['--mode=H'], ['--seedMaxUP=1'], ['--model=X', '--mode=K'], + for opts in (['--mode=S'], ['--seedMaxUP=1'], ['--model=X', '--mode=K'], ['--windowWidth=3'], ['--seedTQ=1|&4|','--energyAdd=900']): failure = d/'failure.svg' run(options=[f'--out=bpsvg:{failure}', *opts], success=False) From 5fc10e99ba46af6ff95720753acbb2fe94c57272 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Fri, 9 Oct 2026 22:34:20 +0200 Subject: [PATCH 7/9] docs: record revision audits and probability performance --- ChangeLog | 8 +++ .../seeded-probabilities-validation.md | 65 +++++++++++++++++++ ...redictorMfeEns2dHeuristicSeedExtension.svg | 4 +- 3 files changed, 75 insertions(+), 2 deletions(-) diff --git a/ChangeLog b/ChangeLog index eaea2fe1..e1911f1d 100644 --- a/ChangeLog +++ b/ChangeLog @@ -65,6 +65,9 @@ ## Technical changes and Optimizations +- document independent release/debug, installed-API and performance validation + for the reviewed unseeded and heuristic probability outputs (PR #259) + - transfer completed probability-output buffers without copying; independently audit performance, accessibility conversion and recurrence depictions (PR #259) @@ -181,6 +184,11 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261009 Alexander Mitrofanov + * doc/analysis/seeded-probabilities-validation.md, doc/recursions : + * record final review validation and controlled runtime/memory measurements; + update the seeded-heuristic depiction for approximate pair output (PR #259) + 261009 Alexander Mitrofanov * PredictorMfeEns2dHeuristic, PredictorMfeEns2dHeuristicSeedExtension, PredictorMfeEns2dSeedExtension, BasePairProbabilities, probability writers : diff --git a/doc/analysis/seeded-probabilities-validation.md b/doc/analysis/seeded-probabilities-validation.md index 56c83716..b3c44bd6 100644 --- a/doc/analysis/seeded-probabilities-validation.md +++ b/doc/analysis/seeded-probabilities-validation.md @@ -318,3 +318,68 @@ and implementation. All revised text is UTF-8/LF and `git diff --check` passes. Both new recurrence SVGs parse, render and passed independent equation review; representative output SVGs were rendered and visually inspected. Apple Clang was not available locally. + +### Controlled revision performance + +A separate agent compared the retained `202651e` executable with the final +implementation at `ae0eb78`, using the same GCC 14.4 release configuration, +ViennaRNA 2.7.2, Boost 1.85 and Kokkos mdspan on the machine described above. +Each configuration has one correctness/warmup run and seven deterministically +shuffled, interleaved fresh-process measurements, pinned to CPU 0 and one thread. +No build or test ran during timing. Both ensemble and probability outputs went to +`/dev/null` for timing; separate file-output checks compared complete ensemble +bytes and every CSV/SVG probability cell. + +Wall seconds (median [minimum, maximum], including process startup): + +| Input / predictor | Before, output off | Revised, output off | Revised CSV | Revised SVG | +| --- | ---: | ---: | ---: | ---: | +| OxyS-fhlA / noSeed-M | 0.4116 [0.4077, 0.4262] | 0.4165 [0.4105, 0.4220] | 0.7309 [0.7236, 0.7500] | 0.8490 [0.8415, 0.8668] | +| OxyS-fhlA / noSeed-H | 0.0644 [0.0641, 0.0661] | 0.0646 [0.0637, 0.0668] | 0.0728 [0.0720, 0.0767] | 0.0971 [0.0962, 0.1041] | +| OxyS-fhlA / seeded-H | 0.0517 [0.0510, 0.0532] | 0.0517 [0.0513, 0.0557] | 0.0566 [0.0562, 0.0569] | 0.0772 [0.0771, 0.0784] | +| dense40 / noSeed-M | 0.3409 [0.3396, 0.3412] | 0.3489 [0.3455, 0.3559] | 0.6691 [0.6619, 0.6936] | 0.7231 [0.7143, 0.7377] | +| dense40 / noSeed-H | 0.0195 [0.0194, 0.0207] | 0.0195 [0.0191, 0.0197] | 0.0227 [0.0223, 0.0234] | 0.0264 [0.0261, 0.0271] | +| dense40 / seeded-H | 1.8790 [1.8324, 1.8980] | 1.8767 [1.8584, 1.8861] | 2.7876 [2.7391, 2.8202] | 2.8354 [2.7850, 2.8720] | +| dense32-noLP / noSeed-M | 0.1001 [0.0991, 0.1021] | 0.1014 [0.1000, 0.1042] | 0.1873 [0.1861, 0.1903] | 0.2041 [0.2023, 0.2078] | +| dense32-noLP / noSeed-H | 0.0159 [0.0158, 0.0169] | 0.0157 [0.0155, 0.0159] | 0.0174 [0.0172, 0.0186] | 0.0198 [0.0195, 0.0204] | +| dense32-noLP / seeded-H | 0.5690 [0.5631, 0.5791] | 0.5692 [0.5652, 0.5783] | 0.8398 [0.8290, 0.9217] | 0.8561 [0.8473, 0.8731] | + +Maximum observed process RSS (KiB): + +| Input / predictor | Before, output off | Revised, output off | Revised CSV | Revised SVG | +| --- | ---: | ---: | ---: | ---: | +| OxyS-fhlA / noSeed-M | 17,180 | 17,308 | 17,312 | 20,844 | +| OxyS-fhlA / noSeed-H | 18,208 | 18,204 | 18,844 | 22,328 | +| OxyS-fhlA / seeded-H | 17,184 | 17,180 | 17,184 | 20,848 | +| dense40 / noSeed-M | 16,388 | 16,392 | 16,256 | 16,292 | +| dense40 / noSeed-H | 16,736 | 16,608 | 16,720 | 16,724 | +| dense40 / seeded-H | 34,740 | 34,744 | 35,564 | 35,696 | +| dense32-noLP / noSeed-M | 16,376 | 16,380 | 16,256 | 16,152 | +| dense32-noLP / noSeed-H | 16,564 | 16,564 | 16,608 | 16,604 | +| dense32-noLP / seeded-H | 21,508 | 21,504 | 22,324 | 22,592 | + +Existing output-off medians vary from −1.2% to +2.4% against the pre-revision +binary across these nine cases. All ensemble outputs are byte-identical across +both binaries and all output configurations; every SVG pair value matches CSV. +Probability collection has a measured additional cost: CSV is 1.09–1.92 times +revised output-off wall time here, and SVG is 1.26–2.07 times. These are local +workload measurements, not a universal overhead bound or speedup claim. + +Reproduce one process from each corresponding build (change mode to H for the +unseeded heuristic; remove `--noSeed` for seeded H): + +```sh +OMP_NUM_THREADS=1 OPENBLAS_NUM_THREADS=1 OMP_DYNAMIC=FALSE \ + /usr/bin/time -f '%e %U %S %M' taskset -c 0 src/bin/IntaRNA \ + --target=doc/handson/fhlA.fasta --query=doc/handson/OxyS.fasta \ + --model=P --mode=M --noSeed --seedBP=4 --intLenMax=30 \ + --qIntLoopMax=3 --tIntLoopMax=3 --threads=1 \ + --outNoLP=false --outNoGUend=false --outNumber=0 --outMode=E \ + --default-log-file=/dev/null --out=STDOUT > /dev/null +# Add --out=bpProb:/dev/null or --out=bpsvg:/dev/null. +``` + +For dense40, replace target/query with G×40/C×40 and add `--tAcc=N --qAcc=N`. +For dense32-noLP, use G×32/C×32 with disabled accessibility and +`--outNoLP=true`. Keep all remaining flags identical. Repeat seven times with +interleaved configurations; capture process RSS with GNU time `%M`. diff --git a/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristicSeedExtension.svg b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristicSeedExtension.svg index 6f7c7fa8..2d38e87d 100644 --- a/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristicSeedExtension.svg +++ b/doc/recursions/IntaRNA.PredictorMfeEns2dHeuristicSeedExtension.svg @@ -1,6 +1,6 @@ Seeded ensemble · one selected right extension -PredictorMfeEns2dHeuristicSeedExtension | --model=P --mode=H. Internal recurrence documentation for step 9 of PR 258. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. +PredictorMfeEns2dHeuristicSeedExtension | --model=P --mode=H. Internal recurrence documentation for steps 9 and 11 of PR 258 / PR 259. Equations describe the current implementation; RNA 2 uses reversed energy coordinates. @@ -42,7 +42,7 @@ All true right-only extensions + all corrected left-only extensions + left extensions paired with just j*. The second phase runs if noGUend is off, or a valid right optimum exists, or anchor end b is not GU. Initialize j*=b and E_right_opt=∞. A positive selected R is required before combining both sides. -No full Cartesian product of left and right boundaries is retained, and bpProb / bpsvg are unsupported. +No full left/right Cartesian product; bpProb / bpsvg trace this retained family and are labelled approximate. Every update is filtered and accumulated by PredictorMfeEns::updateZ; ranked reporting is a separate step. Implementation anchors (src/IntaRNA/) · BW(E)=exp(−E/RT) after native energy-unit conversion From c0ef61f1603ffc41c2ae4f8bbcf98b56027880e9 Mon Sep 17 00:00:00 2001 From: AlexanderMitrofanov Date: Mon, 12 Oct 2026 12:41:32 +0200 Subject: [PATCH 8/9] fix: align probability dot plots and trace MFE pair outlines --- ChangeLog | 13 ++ README.md | 15 ++- .../seeded-probabilities-validation.md | 9 +- src/IntaRNA/BasePairProbabilitySvg.cpp | 38 ++++-- src/IntaRNA/BasePairProbabilityWriter.h | 4 +- src/bin/CommandLineParsing.cpp | 124 +++++++++++++++++- tests/testBasePairProbabilitySvg.py | 45 ++++++- 7 files changed, 223 insertions(+), 25 deletions(-) diff --git a/ChangeLog b/ChangeLog index e1911f1d..5e0eb3f0 100644 --- a/ChangeLog +++ b/ChangeLog @@ -14,6 +14,9 @@ ## Interface and handling +- orient probability SVGs from the upper-left corner and highlight every MFE + pair, splitting seed/MFE cell borders into orange and green (PR #259 review) + - enable approximate actual-pair CSV/SVG probabilities for unseeded and seeded heuristic ensembles, labelled explicitly as approximate (PR #259, step 11) @@ -184,6 +187,16 @@ energies, restricted partition sums, and trackers. ################################################################################ ################################################################################ +261012 Alexander Mitrofanov + * BasePairProbabilitySvg, BasePairProbabilityWriter, CommandLineParsing + * restore top-to-bottom target indexing, aligned frames/guides and MFE region + + reconstruct the MFE structure inside the best ensemble site for SVG only + + highlight each MFE pair; shared seed pairs use orange vertical and green + horizontal edges, preserving both annotations (PR #259 review) + * testBasePairProbabilitySvg.py, README.md, seeded-probabilities-validation.md + + verify complete internal MFE pairs, shared borders, signed guides and + output-off parity; document the fixed-boundary traceback and orientation + 261009 Alexander Mitrofanov * doc/analysis/seeded-probabilities-validation.md, doc/recursions : * record final review validation and controlled runtime/memory measurements; diff --git a/README.md b/README.md index f3089346..7a118b25 100644 --- a/README.md +++ b/README.md @@ -2097,8 +2097,8 @@ IntaRNA -t GGAGGG -q CCCC --energy=B --acc=N --seedBP=2 \ --model=P --mode=M --out=bpsvg:pairs.svg --out=bpProb:pairs.csv ``` -The query runs left to right and the target bottom to top, both in original -5' to 3' order, starting at the bottom-left corner. Sequence names label the +The query runs left to right and the target top to bottom, both in original +5' to 3' order, starting at the upper-left corner. Sequence names label the axes, with the query name below the matrix. A subtitle names both sequences when either name differs from the defaults `query` and `target`. Nucleotides appear on all four sides. Their colored backgrounds show **intramolecular @@ -2127,9 +2127,14 @@ Dark orange outlines mark the union of pairs belonging to **admitted seeds** contained in searched regions, including overlapping, singleton and mixed-length explicit seeds. This annotation uses the predictor's retained seed family and span limits; it is independent of ranked interaction reporting and can include -zero-probability cells. Extension-only pairs have no outline. A seed legend is -included whenever an outline is present. A green rectangle encloses the best -reported interaction's site, with a matching legend. The best site is retained +zero-probability cells. A seed legend is included whenever an outline is present. +A green rectangle encloses the best reported ensemble interaction's site, and +green cell borders mark a minimum-energy structure reconstructed within that +site using the seed, loop and noLP constraints. Ensemble predictors ordinarily +report only site boundaries; the SVG requests this additional MFE traceback. +It does not change site ranking or the probability matrix. Where an MFE pair +also belongs to an admitted seed, its left/right edges are orange and its +top/bottom edges are green. Both annotations have legends. The best site is retained for SVG even with `--outNumber=0`, without adding ranked text output. If no site passes the reporting filters, the rectangle and its legend are omitted. diff --git a/doc/analysis/seeded-probabilities-validation.md b/doc/analysis/seeded-probabilities-validation.md index b3c44bd6..924c7048 100644 --- a/doc/analysis/seeded-probabilities-validation.md +++ b/doc/analysis/seeded-probabilities-validation.md @@ -283,10 +283,15 @@ candidate ensemble can retain path multiplicities and is not the exact seeded unique-chain ensemble. Existing numerical-zero cutoffs remain in the legacy unseeded/heuristic objectives. -The SVG now starts both strands at the bottom-left origin, places the query name +The SVG now starts both strands at the upper-left origin, places the query name below the matrix, uses sequence names and rounded values in tooltips, displays intramolecular pairing as `1-Pu`, and outlines the best reported interaction region -in green. Machine-readable probability attributes keep full precision. The MFE +in green. Individual cells mark a minimum-energy structure reconstructed inside +the selected ensemble site. Shared seed/MFE cells have orange vertical edges and +green horizontal edges. This SVG-only fixed-boundary traceback reuses the MFE +recurrences and adds computation proportional to the selected site, without +changing ordinary ensemble output. Machine-readable probability attributes keep +full precision. The MFE outline is available with `--outNumber=0`; ordinary ranked output stays suppressed. An outline is omitted if no interaction passes the reporting filters. diff --git a/src/IntaRNA/BasePairProbabilitySvg.cpp b/src/IntaRNA/BasePairProbabilitySvg.cpp index 5b7620ee..349f2c4f 100644 --- a/src/IntaRNA/BasePairProbabilitySvg.cpp +++ b/src/IntaRNA/BasePairProbabilitySvg.cpp @@ -6,6 +6,7 @@ #include #include #include +#include #include #include #include @@ -77,9 +78,10 @@ BasePairProbabilityWriter::svgBlock(const BasePairProbabilities & result,const I const double x=std::max(130.0,(std::max(820.0,nameWidth)-w)/2); const double y=std::max(120.0,(9.0*t.getId().size()-h)/2+60); const double width=2*x+w, height=2*y+h+185; - // SVG y increases downwards, while target coordinates increase upwards. - auto targetY=[&](size_t i) { return y+cell*(t.size()-1-i); }; + // Both original sequence indices increase from the upper-left corner. + auto targetY=[&](size_t i) { return y+cell*i; }; const bool markMfe=mfe && !mfe->isEmpty(); + std::set mfePairs; size_t tStart=t.size(),tEnd=0,qStart=q.size(),qEnd=0; if (markMfe) { for (const auto & bp:mfe->basePairs) { @@ -87,6 +89,7 @@ BasePairProbabilityWriter::svgBlock(const BasePairProbabilities & result,const I throw std::invalid_argument("bpsvg: MFE interaction coordinates exceed sequence dimensions"); tStart=std::min(tStart,bp.first); tEnd=std::max(tEnd,bp.first); qStart=std::min(qStart,bp.second); qEnd=std::max(qEnd,bp.second); + mfePairs.insert(bp); } } std::ostringstream out; @@ -95,7 +98,7 @@ BasePairProbabilityWriter::svgBlock(const BasePairProbabilities & result,const I out<<"\n" <<""<<(result.isApproximate()?"Approximate base-pair probabilities: ":"Base-pair probabilities: ")<\n" - <<"Target rows run 5' to 3' upwards and query columns run 5' to 3' to the right. Pair probabilities are conditional on an allowed interaction. " + <<"Target rows run 5' to 3' downwards and query columns run 5' to 3' to the right, starting at the upper-left corner. Pair probabilities are conditional on an allowed interaction. " <<"The frame shows single-nucleotide intramolecular pairing probabilities (1 - Pu) from the accessibility model. " <<(result.isApproximate()?"Approximate probabilities describe the retained heuristic ensemble.":"Probabilities describe the exact predictor ensemble.")<<"\n" < @@ -132,7 +135,7 @@ text { font-family: sans-serif; fill: #243247; } text(width/2,53,t.getId()+" / "+q.getId(),"subtitle"); text(x+w/2,y+h+65,q.getId()+" (5' to 3')","axis"); out<<""<\n"; + <<") rotate(90)\" text-anchor=\"middle\">"<\n"; auto rect=[&](double rx,double ry,Z_type value,bool missing,const std::string & attributes,const std::string & title) { out<<"0 && !origin?1:0)); if (idx%10==0 || origin) out<<"\n"; + <<"\" x1=\""<<(target?x-cell:x+offset)<<"\" y1=\""<<(target?y+offset:y-cell) + <<"\" x2=\""<<(target?x+w+cell:x+offset)<<"\" y2=\""<<(target?y+offset:y+h+cell)<<"\"/>\n"; const std::string label=idx==1?"+1":std::to_string(idx); if(target) { text(x-26,targetY(i)+cell/2+4,label,"index","end"); @@ -200,17 +203,35 @@ text { font-family: sans-serif; fill: #243247; } out<<"\n"; for(size_t i=0;i\n" + <<"\n" + <<"\n\n"; + continue; + } out<<"\n"; } out<<"\n"; } if (markMfe) { + out<<"\n"; + for (const auto & bp:mfePairs) { + if (result.collectsSeedPairs() && result.seedPairs()(bp.first,bp.second)) continue; + out<<"\n"; + } + out<<"\n"; out<<"\n"; } const double legendX=(width-7*86)/2, legendY=y+h+110; @@ -231,7 +252,8 @@ text { font-family: sans-serif; fill: #243247; } if (markMfe) { const double mfeLegendY=legendY+(annotated?76:48); out<<"\n"; - text(legendX+28,mfeLegendY+13,"Green outline: minimum-free-energy interaction region","legend","start"); + text(legendX+28,mfeLegendY+13,"Green outline: minimum-free-energy interaction region and base pairs","legend","start"); + if (annotated) text(legendX+28,mfeLegendY+31,"Seed + MFE pair: orange left/right edges, green top/bottom edges","legend","start"); } if(empty) text(width/2,legendY+90,"Gray / NA: empty interaction ensemble; accessibility remains defined","legend"); out<<"\n"; diff --git a/src/IntaRNA/BasePairProbabilityWriter.h b/src/IntaRNA/BasePairProbabilityWriter.h index 64773536..31ec316c 100644 --- a/src/IntaRNA/BasePairProbabilityWriter.h +++ b/src/IntaRNA/BasePairProbabilityWriter.h @@ -41,7 +41,7 @@ class BasePairProbabilityWriter { * @param targetRT positive finite energy scale used to convert target ED to Pu * @param queryRT positive finite energy scale used to convert query ED to Pu; * accessibility scales can differ from the interaction model's RT - * @param mfe optional best interaction in original coordinates for a region outline + * @param mfe optional traced MFE structure in original coordinates for region/pair outlines */ static void writeSvg(std::ostream & out,const BasePairProbabilities & result, const InteractionEnergy & energy,Z_type targetRT,Z_type queryRT,const Interaction * mfe = nullptr); @@ -52,7 +52,7 @@ class BasePairProbabilityWriter { * @param energy original, unshifted energy/accessibility model * @param targetRT positive finite target accessibility energy scale * @param queryRT positive finite query accessibility energy scale - * @param mfe optional best interaction in original coordinates for a region outline + * @param mfe optional traced MFE structure in original coordinates for region/pair outlines */ static void writeSvgFile(const std::string & filename,const BasePairProbabilities & result, const InteractionEnergy & energy,Z_type targetRT,Z_type queryRT,const Interaction * mfe = nullptr); diff --git a/src/bin/CommandLineParsing.cpp b/src/bin/CommandLineParsing.cpp index 623235d3..6650cde4 100644 --- a/src/bin/CommandLineParsing.cpp +++ b/src/bin/CommandLineParsing.cpp @@ -38,6 +38,7 @@ extern "C" { #include "IntaRNA/InteractionEnergyBasePair.h" #include "IntaRNA/InteractionEnergyVrna.h" +#include "IntaRNA/OutputHandlerInteractionList.h" #include "IntaRNA/PredictorMfe2dHeuristic.h" #include "IntaRNA/PredictorMfe2dHelixBlockHeuristic.h" @@ -2503,15 +2504,134 @@ BasePairProbabilities * CommandLineParsing::getBasePairProbabilityResult(const I } return new BasePairProbabilities(energy.size1(),energy.size2(),!noSeedRequired && !outPrefix2streamName.at(OP_bpsvg).empty()); } +namespace { + +// Ensemble output stores boundaries, not a molecular structure. Reuse the +// existing fixed-boundary MFE recurrences solely for the SVG's pair overlay. +class SvgMfeTrace : public PredictorMfe2d { +public: + using PredictorMfe2d::PredictorMfe2d; + void trace(Interaction & interaction); +}; + +void +SvgMfeTrace::trace(Interaction & interaction) +{ + interaction.sort(); + const auto first=interaction.basePairs.front(), last=interaction.basePairs.back(); + energy.setOffset1(first.first); + energy.setOffset2(energy.getAccessibility2().getReversedIndex(first.second)); + hybridE_pq.resize(last.first-first.first+1,first.second-last.second+1); + traceBack(interaction); +} + +class SvgSeededMfeTrace : public PredictorMfe2dSeed { +public: + using PredictorMfe2dSeed::PredictorMfe2dSeed; + bool trace(Interaction & interaction); +}; + +bool +SvgSeededMfeTrace::trace(Interaction & interaction) +{ + interaction.sort(); + const auto first=interaction.basePairs.front(), last=interaction.basePairs.back(); + const size_t offset2=energy.getAccessibility2().getReversedIndex(first.second); + energy.setOffset1(first.first); energy.setOffset2(offset2); + seedHandler.setOffset1(first.first); seedHandler.setOffset2(offset2); + const size_t n=last.first-first.first+1, m=first.second-last.second+1; + hybridE_pq.resize(n,m); hybridE_pq_seed.resize(n,m); + if (seedHandler.fillSeed(0,n-1,0,m-1)==0) return false; + fillHybridE(n-1,m-1,0,0,false); + if (E_isINF(hybridE_pq_seed(0,0))) return false; + interaction.energy=energy.getE(0,n-1,0,m-1,hybridE_pq_seed(0,0)); + traceBack(interaction); + return true; +} + +// Explicit seeds define their own lengths, including singletons. The legacy +// MFE traceback's nominal length check must not reject shorter explicit seeds. +class SvgSeedConstraint : public SeedConstraint { +public: + explicit SvgSeedConstraint(const SeedConstraint & original); +}; + +inline +SvgSeedConstraint::SvgSeedConstraint(const SeedConstraint & original) : SeedConstraint(original) +{ + if (!getExplicitSeeds().empty()) bp=2; +} + +} // namespace + void CommandLineParsing::writeBasePairProbabilities(const BasePairProbabilities & result,const InteractionEnergy & energy, const Interaction * mfe) const { const auto & t=energy.getAccessibility1().getSequence(); const auto & q=energy.getAccessibility2().getAccessibilityOrigin().getSequence(); if (!outPrefix2streamName.at(OP_bpProb).empty()) BasePairProbabilityWriter::writeFile(getFullFilename(outPrefix2streamName.at(OP_bpProb),&t,&q),result,t,q,outSep); - if (!outPrefix2streamName.at(OP_bpsvg).empty()) + if (!outPrefix2streamName.at(OP_bpsvg).empty()) { + std::unique_ptr structure; + if (mfe && !mfe->isEmpty()) { + structure=std::make_unique(*mfe); + const auto original=getOutputConstraint(energy,true); + OutputConstraint constraint(1,OutputConstraint::OVERLAP_BOTH,E_INF,E_INF, + original.bestSeedOnly,original.noLP,original.noGUend,false,true,original.maxED); + OutputHandlerInteractionList unused(constraint,1); + if (noSeedRequired) { + SvgMfeTrace trace(energy,unused,nullptr); + trace.trace(*structure); + } else if (structure->basePairs.size()>1) { + SvgSeedConstraint seedConstraint(getSeedConstraint(energy)); + if (seedConstraint.getExplicitSeeds().empty()) { + SvgSeededMfeTrace trace(energy,unused,nullptr,getSeedHandler(energy)); + if (!trace.trace(*structure)) throw std::runtime_error("bpsvg: no MFE structure in selected seeded site"); + } else { + // Trace each explicit seed separately so mixed lengths and + // alternatives at the same start remain available in this site. + std::vector variants; + const auto first=structure->basePairs.front(), last=structure->basePairs.back(); + const size_t lo2=energy.getAccessibility2().getReversedIndex(first.second); + const size_t hi2=energy.getAccessibility2().getReversedIndex(last.second); + std::istringstream encodings(seedConstraint.getExplicitSeeds()); + for (std::string encoding; std::getline(encodings,encoding,','); ) { + seedConstraint.getExplicitSeeds()=encoding; + SeedHandlerExplicit seed(energy,seedConstraint); + size_t a=RnaSequence::lastPos,b=RnaSequence::lastPos; + if (!seed.updateToNextSeed(a,b,first.first,last.first,lo2,hi2)) continue; + if (seed.getSeedLength1(a,b)>1) { variants.push_back(encoding); continue; } + // A singleton in a multi-pair chain has a neighbor. Expand + // it into each feasible two-pair seed; noLP requires a stack. + for (int direction : {-1,1}) { + const size_t max1=original.noLP?1:energy.getMaxInternalLoopSize1()+1; + const size_t max2=original.noLP?1:energy.getMaxInternalLoopSize2()+1; + for (size_t d1=1; d1<=max1; ++d1) for (size_t d2=1; d2<=max2; ++d2) { + if ((direction<0 && (d1>a-first.first || d2>b-lo2)) + || (direction>0 && (d1>last.first-a || d2>hi2-b))) continue; + const size_t c=direction<0?a-d1:a+d1, d=direction<0?b-d2:b+d2; + if (!energy.areComplementary(c,d)) continue; + variants.push_back(std::to_string(t.getInOutIndex(std::min(a,c)))+"|"+std::string(d1-1,'.')+"|&" + +std::to_string(q.getInOutIndex(energy.getAccessibility2().getReversedIndex(std::max(b,d)))) + +"|"+std::string(d2-1,'.')+"|"); + } + } + } + E_type best=E_INF; + for (const auto & encoding : variants) { + seedConstraint.getExplicitSeeds()=encoding; + SvgSeededMfeTrace trace(energy,unused,nullptr,new SeedHandlerExplicit(energy,seedConstraint)); + Interaction candidate(*mfe); + if (trace.trace(candidate) && candidate.energy max(float(r.attrib['y']) for r in cells.values())+20 - # Original target coordinate 1 is at the bottom; query coordinate 1 at left. + # Both lowest original indices are at the upper-left corner. for i in range(1, 7): for j in range(1, 5): - assert float(cells[i, j].attrib['y']) == float(cells[1, 1].attrib['y'])-20*(i-1) + assert float(cells[i, j].attrib['y']) == float(cells[1, 1].attrib['y'])+20*(i-1) assert float(cells[i, j].attrib['x']) == float(cells[1, 1].attrib['x'])+20*(j-1) for i, row in enumerate(rows[1:], 1): for j, value in enumerate(row[1:], 1): @@ -118,10 +125,36 @@ def plot(t='GGAGGG', q='CCCC', opts=()): ts, te, qs, qe = bounds matrix = pairs(mfe_plot) assert float(outline.attrib['x']) == float(matrix[ts, qs].attrib['x']) - assert float(outline.attrib['y']) == float(matrix[te, qs].attrib['y']) + assert float(outline.attrib['y']) == float(matrix[ts, qs].attrib['y']) assert float(outline.attrib['width']) == float(matrix[ts, qe].attrib['x'])-float(matrix[ts, qs].attrib['x'])+20 - assert float(outline.attrib['height']) == float(matrix[ts, qs].attrib['y'])-float(matrix[te, qs].attrib['y'])+20 + assert float(outline.attrib['height']) == float(matrix[te, qs].attrib['y'])-float(matrix[ts, qs].attrib['y'])+20 assert 'Green outline:' in ''.join(mfe_plot.itertext()) + assert len(mfe_pairs(mfe_plot)) > 2 # includes internal pairs, not just site endpoints + # An unambiguous four-pair stack checks every internal MFE pair against + # independent ordinary MFE output, including the two-color seed overlap. + for seed_options in ([], ['--noSeed'], ['--seedMaxUP=1'], ['--seedTQ=2||&2||']): + for mode in ('M', 'H'): + if '--seedMaxUP=1' in seed_options and mode == 'M': + continue # exact bulged-seed probability output is unsupported + stack = plot('GGGG', 'CCCC', [*seed_options, f'--mode={mode}']) + reference = run('GGGG', 'CCCC', [*seed_options, '--model=X', '--mode=M', + '--outNumber=1', '--outMode=C', '--outCsvCols=bpList']) + row = next(csv.DictReader(reference.stdout.splitlines(), delimiter=';')) + expected = {tuple(map(int, bp)) for bp in re.findall(r'\((\d+),(\d+)\)', row['bpList'])} + assert set(mfe_pairs(stack)) == expected == {(1,4), (2,3), (3,2), (4,1)} + for bp, outline in mfe_pairs(stack).items(): + cell = pairs(stack)[bp] + left, top = float(cell.attrib['x'])+1, float(cell.attrib['y'])+1 + if bp in seeds(stack): + assert outline.tag.endswith('g') + assert outline.find('s:path[@class="seed"]', NS).attrib['d'] == f'M {left:g} {top:g} v 18 M {left+18:g} {top:g} v 18' + assert outline.find('s:path[@class="mfe"]', NS).attrib['d'] == f'M {left:g} {top:g} h 18 M {left:g} {top+18:g} h 18' + else: + assert float(outline.attrib['x']) == left + assert float(outline.attrib['y']) == top + # A sole admitted seed fixes a shorter best site; singleton plotting remains valid. + singleton = plot('GGGG', 'CCCC', ['--seedTQ=4|&1|']) + assert (4,1) in mfe_pairs(singleton) unseeded = plot(opts=['--noSeed', f'--out=bpProb:{table}']) assert not seeds(unseeded) rows = list(csv.reader(table.open(), delimiter=';')) @@ -161,7 +194,7 @@ def plot(t='GGAGGG', q='CCCC', opts=()): idx = int(g.attrib['data-index']) r = signed_cells[idx,-15] if strand == 'target' else signed_cells[-15,idx] axis = 'y' if strand == 'target' else 'x' - expected_pos = float(r.attrib[axis]) + (0 if idx > 1 else 20) if strand == 'target' else float(r.attrib[axis]) + (20 if idx > 1 else 0) + expected_pos = float(r.attrib[axis]) + (20 if idx > 1 else 0) assert float(g.attrib[axis+'1']) == expected_pos assert ('major' in g.attrib['class']) == (idx % 50 == 0) assert ('origin' in g.attrib['class']) == (idx == 1) From 0d61d4b07caca760d64a46ee82afc4e3730ba034 Mon Sep 17 00:00:00 2001 From: Martin Raden Date: Tue, 13 Oct 2026 11:25:31 +0200 Subject: [PATCH 9/9] tmp folder --- .gitignore | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.gitignore b/.gitignore index 7654a03d..ee819cc0 100644 --- a/.gitignore +++ b/.gitignore @@ -55,7 +55,7 @@ doxygen-doc doxygen-doc/* conf3108* test-driver - +tmp/* # temp dir tmp