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
diff --git a/ChangeLog b/ChangeLog
index bd016930..5e0eb3f0 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -14,6 +14,28 @@
## 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)
+
+- 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)
+
+- 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
+
+- --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)
@@ -46,6 +68,15 @@
## 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)
+
+- 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)
@@ -156,6 +187,99 @@ 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;
+ update the seeded-heuristic depiction for approximate pair output (PR #259)
+
+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;
+ 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
+ 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
+ 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 :
+ + 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
+ 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..7a118b25 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 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)
@@ -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.
@@ -1892,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!
@@ -2012,8 +2016,9 @@ 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, or
+`--model=P --mode=H` for an approximate heuristic ensemble. For example:
```sh
IntaRNA -t GGAGGG -q CCCC --energy=B --acc=N --seedBP=2 \
@@ -2026,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, unseeded 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
@@ -2050,20 +2062,98 @@ 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.
+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`, and call
+pass its optional non-owning pointer to `PredictorMfeEns2dSeedExtension` (stack
+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)
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 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 \
+ --model=P --mode=M --out=bpsvg:pairs.svg --out=bpProb:pairs.csv
+```
+
+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
+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)`,
+`[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,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");
+ 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 {
+ 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");
+ }
+ }
+ }
+ 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,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");
+ }
+ }
+ }
+ out<<"\n\n";
+ bool annotated=false;
+ if(result.collectsSeedPairs()) {
+ 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;
+ text(width/2,legendY-14,"Probability ranges (intermolecular base pairs and intramolecular pairing)","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 (markMfe) {
+ const double mfeLegendY=legendY+(annotated?76:48);
+ out<<"\n";
+ 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";
+ // 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,
+ Z_type targetRT,Z_type queryRT,const Interaction * mfe)
+{
+ 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,const Interaction * mfe)
+{
+ emitFile(filename,svgBlock(result,energy,targetRT,queryRT,mfe));
+}
+
+} // namespace IntaRNA
diff --git a/src/IntaRNA/BasePairProbabilityWriter.cpp b/src/IntaRNA/BasePairProbabilityWriter.cpp
index d0aea7af..b8ee89bf 100644
--- a/src/IntaRNA/BasePairProbabilityWriter.cpp
+++ b/src/IntaRNA/BasePairProbabilityWriter.cpp
@@ -3,6 +3,7 @@
#include
#include
#include
+#include
#include
namespace IntaRNA {
@@ -17,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 file(out==&std::cout || out==&std::cerr?nullptr:out);
diff --git a/src/IntaRNA/BasePairProbabilityWriter.h b/src/IntaRNA/BasePairProbabilityWriter.h
index 48578daf..31ec316c 100644
--- a/src/IntaRNA/BasePairProbabilityWriter.h
+++ b/src/IntaRNA/BasePairProbabilityWriter.h
@@ -6,6 +6,8 @@
#include
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.
*/
@@ -30,7 +32,34 @@ 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 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
+ * @param energy original, unshifted energy/accessibility model used in prediction
+ * @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 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);
+ /** 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
+ * @param targetRT positive finite target accessibility energy scale
+ * @param queryRT positive finite query accessibility energy scale
+ * @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);
private:
+ static std::string svgBlock(const BasePairProbabilities & result,const InteractionEnergy & energy,
+ 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);
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/PartitionArithmetic.h b/src/IntaRNA/PartitionArithmetic.h
index 37f28a89..c5d8e107 100644
--- a/src/IntaRNA/PartitionArithmetic.h
+++ b/src/IntaRNA/PartitionArithmetic.h
@@ -4,7 +4,7 @@
#include
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 9b9bb28e..a03506da 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,84 @@ 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 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));
+ }
+ }
+
+ 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..2b77bb67 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,15 @@ 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.
+ * Collection reverses this class's fillHybridZ recurrence; subclasses
+ * changing that recurrence must also supply their own collection.
*/
PredictorMfeEns2d( const InteractionEnergy & energy
, OutputHandler & output
- , PredictionTracker * predTracker );
+ , PredictionTracker * predTracker
+ , BasePairProbabilities * pairProbabilities = nullptr );
virtual ~PredictorMfeEns2d();
@@ -58,6 +64,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 +79,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/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 4f50a61e..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
@@ -185,11 +368,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);
}
}
@@ -245,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
@@ -290,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)
@@ -304,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
@@ -349,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 {
@@ -364,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; ) {
@@ -379,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
@@ -388,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
@@ -406,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
@@ -450,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}});
}
}
@@ -466,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 dc864b8c..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"
@@ -957,7 +958,8 @@ 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 '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 '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/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())
@@ -1245,8 +1247,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 requires --model=P --mode=M with 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.
@@ -1582,7 +1584,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.");
@@ -2233,7 +2235,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()) {
@@ -2249,7 +2251,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)
@@ -2490,19 +2492,146 @@ 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());
+ if (!noSeedRequired) {
+ std::unique_ptr handler(getSeedHandler(energy));
+ 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());
}
-void CommandLineParsing::writeBasePairProbabilities(const BasePairProbabilities & result,const InteractionEnergy & energy) const {
+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();
- 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()) {
+ 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 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");
- pairSeeds.reset(getSeedHandler(energy));
- if(!pairSeeds->guaranteesStackOnlySeeds())
- throw std::invalid_argument("bpProb requires a stack-only seed handler; bulged seeds are unsupported");
+ 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(mode.val=='M' && !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();
@@ -2616,8 +2747,8 @@ 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 'M' : return new PredictorMfeEns2d( 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));
}
} break;
@@ -2661,7 +2792,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));
@@ -2936,9 +3067,18 @@ getTargetRanges( const InteractionEnergy & energy, const size_t sequenceNumber,
////////////////////////////////////////////////////////////////////////////
+Z_type
+CommandLineParsing::getAccessibilityRT(const char mode) const
+{
+ // Computed Nussinov ED uses RT=1; imported ED/Pu and ViennaRNA use the
+ // configured physical temperature, independently of the interaction model.
+ return energy.val=='B' && mode=='C'?Z_type(1):vrnaHandler.getRT();
+}
+
void
CommandLineParsing::
-writeAccessibility( const Accessibility& acc, const std::string & fileOrStream, const bool writeED ) const
+writeAccessibility( const Accessibility& acc, const std::string & fileOrStream, const bool writeED,
+ const char mode ) const
{
if (fileOrStream.empty())
return;
@@ -2959,7 +3099,7 @@ writeAccessibility( const Accessibility& acc, const std::string & fileOrStream,
} else if (writeED) {
acc.writeRNAplfold_ED_text( *out );
} else {
- acc.writeRNAplfold_Pu_text( *out, vrnaHandler.getRT() );
+ acc.writeRNAplfold_Pu_text( *out, getAccessibilityRT(mode) );
}
// Explicit close propagates compression/file errors before releasing ownership.
diff --git a/src/bin/CommandLineParsing.h b/src/bin/CommandLineParsing.h
index fff52480..f0e796ad 100644
--- a/src/bin/CommandLineParsing.h
+++ b/src/bin/CommandLineParsing.h
@@ -282,21 +282,22 @@ 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
*/
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
*/
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
@@ -416,6 +418,7 @@ class CommandLineParsing {
OP_spotProb,
OP_spotProbAll,
OP_bpProb,
+ OP_bpsvg,
OP_UNKNOWN
};
@@ -442,6 +445,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
@@ -1184,9 +1188,16 @@ class CommandLineParsing {
* Writes the accessibility to file or stream if requested by the user
* @param acc the accessibility data assigned
* @param fileOrStream the name of file/stream to write to
- * @param (true) writes ED values, (false) writes Pu values
+ * @param writeED (true) writes ED values, (false) writes Pu values
+ * @param mode accessibility source mode of the strand (C/N/P/E)
*/
- void writeAccessibility( const Accessibility& acc, const std::string & fileOrStream, const bool writeED ) const;
+ void writeAccessibility( const Accessibility& acc, const std::string & fileOrStream, const bool writeED,
+ const char mode ) const;
+ /** Energy scale for converting a strand's accessibility ED to Pu.
+ * @param mode accessibility source mode of the strand (C/N/P/E)
+ * @return RT used by that accessibility source, independent of interaction RT
+ */
+ Z_type getAccessibilityRT(const char mode) const;
/**
* Adds a generic file prefix for input/output files for the given query
@@ -1628,14 +1639,14 @@ writeQueryAccessibility( const Accessibility & acc ) const
writeAccessibility( acc
// get file name prefixed with sequence number if needed
, getFullFilename(outPrefix2streamName.at(OutPrefixCode::OP_qAcc), NULL, &(acc.getSequence()))
- , true );
+ , true, qAcc.val );
}
if (!outPrefix2streamName.at(OutPrefixCode::OP_qPu).empty()) {
VLOG(2) <<"writing unpaired probabilities for query '"<::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/BasePairProbabilities_test.cpp b/tests/BasePairProbabilities_test.cpp
index 500c306c..9cffec2d 100644
--- a/tests/BasePairProbabilities_test.cpp
+++ b/tests/BasePairProbabilities_test.cpp
@@ -48,3 +48,51 @@ 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,1,1));
+ REQUIRE(pending.str().empty());
+ r.finalize();std::ostringstream svg;
+ for (Z_type rt : {Z_type(0), Z_type(-1), std::numeric_limits::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,1,1));
+ r.fail();std::ostringstream failed;
+ 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));
+ 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,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,1,1));
+ REQUIRE(wrong.str().empty());
+}
diff --git a/tests/Makefile.am b/tests/Makefile.am
index 5a919fb1..7ff2a868 100644
--- a/tests/Makefile.am
+++ b/tests/Makefile.am
@@ -30,6 +30,9 @@ runApiTests_SOURCES = \
SeededNativeBenchmark_test.cpp \
StackSeedDomain_test.cpp \
BasePairProbabilities_test.cpp \
+ PredictorUnseededProbability_test.cpp \
+ PredictorHeuristicProbability_test.cpp \
+ PredictorHeuristicSeededProbability_test.cpp \
testEasyLoggingSetup.icc \
AccessibilityConstraint_test.cpp \
AccessibilityFromStream_test.cpp \
@@ -85,4 +88,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/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