Mathews 2004 parameter set (corrected) + paramset plumbing + dangles=2 exterior stacks

[KowalskiBio]

Summary

Stacked behind #11 (native-tm) and #12 (native-dp) — review/merge those first: this branch builds on native-dp, so until #12 merges this PR's diff also shows #11's and #12's commits. After #12 merges, the effective diff here is exactly the 2 commits at the tip.

This PR adds a complete, corrected Mathews 2004 DNA parameter set (strider/thermo/parameters/mathews2004-dna.json) with the plumbing to use it at every API surface, plus a dangles=2 exterior dangling-end option on the MFE DP. All opt-in: with no parameter_set and dangles=0 (both defaults), behavior is bit-identical to v1.2.x.

What's in here

1. strider/thermo/parameters/mathews2004-dna.json + scripts/generate_mathews2004_params.py

The JSON is generated from ViennaRNA's dna_mathews2004.par at development time (docstring in the generator covers the license argument: numeric thermodynamic values are published scientific data; the generator is a dev tool, no runtime dependency). During review-by-comparison against ViennaRNA's own energy evaluation, three convention bugs were caught in my first draft of the generator — all fixed here, each verified against RNA.eval_structure_verbose:

• Stack key orientation: the param-file matrix's column pair is written reversed (from the j-side). Keys must be outer[0] + inner[1] + inner[0] + outer[1]. Verified: GCGC/CGCG = −2.2, CCGG/GGCC = −1.8, AATT/TTAA = −1.0, ATAT = −0.9, TATA = −0.6.
• bulge_size off-by-one: a leading 0.0 "sentinel" made index 0 (the 1-nt bulge entry) free instead of 2.9 — the root cause of DP-absurd single-bulge hairpin structures. Every larger bulge was also shifted by one size.
• interior_size off-by-one + INF sentinels (same alignment class as bulge).
• dH interior_1_1 / interior_2_2 tables: they were computed by the converter but never emitted, so structure_enthalpy silently fell back to native ΔG values inside the ΔH walk.

Verification (built into the generator, verify_against_vienna under python scripts/generate_mathews2004_params.py): 7/7 test structures match ViennaRNA eval_structure to 0.00 kcal/mol at dangles=0, including three bulge-containing hairpins that would have caught the bulge misalignment (the previous straight-stem-only set passed with the bug present).

2. Paramset plumbing

• dimer_thermo / dimer_thermo_subopt / dimer_tm accept paramset (instance or name) and use it to predict and score the inter-strand helix, symmetric with the hairpin API.
• hairpin_thermo folds with the same paramset it scores with (previously it folded native and then scored custom — two different models).
• param_context and the four structure_* walkers resolve parameter-set names via a memoized loader (load_parameters string input previously crashed with 'str' object has no attribute 'dG').
• hairpin_tm, fraction_folded accept paramset/dangles; fraction_folded returns 0.0 for non-folding sequences instead of raising ValueError (a non-folding beacon's melt is a flat zero curve).

3. dangles option (default 0 → bit-identical)

Exterior dangling-end stacking at W-matrix stem placements in the linear-kernel DP (fold_mfe, subopt_structures, fold_complex, ThermoEngine(dangles=...)). Rule: best negative single dangle per exterior stem end (equivalent to the d1/d2 outer-loop behavior ViennaRNA reports for single-stem structures).

Empirically verified against ViennaRNA:

• 25-case grid of single-stem hairpins with 0–1 nt tails: single-flank cases match exactly, two-flank cases track the better of the two flanks (±0.2–0.9, and VR's complementary-flank zero anomaly is documented in the code).
• End-to-end: 38/42 sequences within 0.4 kcal/mol of VR(dangles=2); remaining misses are structures requiring G·T wobble pairs (unsupported by design).
• The hairpin ΔG/ΔH walk (_sum_elements) applies the identical exterior-tail dangle, so hairpin_thermo stays energy-consistent with the dangle-aware MFE; dangles=0 keeps every existing energy bit-identical (validated: fold↔walk consistency 13/13).

Why

Matching IDT OligoAnalyzer's published numbers for oligo secondary structure (their analysis family encodes Mathews-2004 DNA NN + dangle stacking at helix termini at 25 °C). Two worked examples from the app where this moves strider from "obviously wrong" to "within noise of IDT":

• reference self-dimer ΔG: −3.54 (native) vs IDT −3.61;
• a previously-broken primer hairpin: −0.79/Tm 30.8 (mathews, dangles=2) vs IDT −0.77/35.9 — versus −2.19 with a bogus bulged structure before the table fix.

Test plan

• Repo suite: 657 passed, 4 pre-existing environment failures (TestNativeViennaParity dimer envelope vs local ViennaRNA — verified identical on the pristine v1.2.x baseline); 10 skipped, 1 xfailed (unchanged baseline).
• Generator self-verification vs ViennaRNA: 7/7 exact.
• dangles=0 regression: energies bit-identical across the test battery (42/42 fold↔walk consistent).

[EmilioVenegas]

Holding off on approval until a few things are settled, but the parameter-set work looks solid. The three convention bugs you caught (stack key column orientation, the bulge_size leading zero that made 1-nt bulges free, and the dH interior tables that were computed but never emitted) all affected real results, and the point that the old straight-stem-only verification passed with the bulge bug present is a good case for the expanded set. Folding hairpin_thermo with the same paramset it scores with is a genuine fix. The provenance note follows the precedent already set for rna_turner2004.par.

Six things before this lands:

1. Tests. Neither code commit adds one: a 12k-line parameter set, a new dangles option on four entry points, paramset threading through six public functions, and a behavior change in fraction_folded. The 658 passing tests are the pre-existing suite, and verify_against_vienna needs ViennaRNA at dev time so it cannot run in CI. Enough to hold it in place would be a few pinned mathews2004-dna energies, the "42/42 dangles=0 bit-identical" claim as a regression test, the paramset-by-name path that used to raise 'str' object has no attribute 'dG', and a fraction_folded zero case.

2. dangles does not reach the partition function. It is threaded into _mfe_native and subopt only, so ThermoEngine(dangles=2) computes its MFE and its ensemble under different models, which affects ensemble defect, the differentiable path and concentration solving. The cache key gained d{self.dangles} for every op including pfunc, which implies an effect that is not there. Either thread it through the McCaskill DP, or scope it explicitly to MFE and drop it from the pfunc key. Your call which, the second is much cheaper.

3. The dangles=2 name. Best single negative dangle per exterior stem is not VR d2, which adds both flanks unconditionally, and your own numbers show it: two-flank cases drift 0.2 to 0.9, 38/42 within 0.4 kcal/mol. The docstring says "like VR md.dangles=2" and the ValueError says "2 (VR dangles=2)". Since strider is judged on agreement with IDT and ViennaRNA, I would rename it or note the divergence at the API surface.

4. The dimer path has no dangles. hairpin_thermo(paramset=..., dangles=2) works, dimer_thermo(paramset=..., dangles=2) raises TypeError, and _sum_dimer_elements never sees the flag, though the motivating example in the PR body is a self-dimer.

5. fraction_folded catches ValueError too broadly. The try wraps the whole hairpin_thermo call, so a mistyped salt_model or bad material returns 0.0 instead of raising. Only the "does not fold into a hairpin" case should map to zero.

6. Citation. THIRD_PARTY_NOTICES cites Mathews et al. 1999, J Mol Biol 288:911-940 for a file named dna_mathews2004.par feeding mathews2004-dna.json. ViennaRNA ships both a 1999 and a 2004 DNA set, so worth confirming which this is before the preprint cites it.

Two smaller ones: _resolve_name is lru_cached and returns the same ParameterSet instance to every caller, so worth confirming nothing mutates one in place given the temperature-adjusted paramset machinery. And ext_e's 5' branch has a redundant (ko - 1) not in nicks check that drops a valid dangle when the flanking base starts a strand, where the 3' branch tests only jo + 1.

hairpin_thermo with a paramset and fraction_folded both change existing behavior, so they want changelog entries. #11 is merged, so this will need a rebase once #12 lands.

[KowalskiBio]

All six items, the two smaller ones, and the changelog are in at 0b0145a, now rebased into a clean stack on merged #11 + the #12 tip (main → 7ef9cf0 → 2b81703 → 736419d → 38f2036 → 04314e4 → c23cb1b → 0b0145a). A full rebase onto main alone will follow the #12 merge.

1. Tests — new tests/test_mathews2004.py, 12/12:

• test_mathews2004_stack_markers_match_the_2004_par_file: pinned identity markers (ATAT=-0.9 excluded the 1999 par, GT/CG=+1.2, GATT, AATT, CGCG).
• test_pinned_hairpin_energies_mathews2004: pinned (ΔG37, Tm, structure) triples on 3 pinned hairpins so the 12k-row JSON cannot silently drift in CI (ViennaRNA still only needed at generation time).
• test_paramset_by_name_matches_instance: name-vs-instance equality, the historical 'str' object has no attribute 'dG' path.
• test_dangles_zero_bit_identical_grid: 42-case grid, dangles=0 vs no-flag bit-equal (structure, energy, pairs).
• test_fraction_folded_zero_for_unfoldable_sequence + test_fraction_folded_bad_stem_for_tan_chen_raises (5-bp stem + salt_model="tan_chen" raises; the old blanket catch swallowed exactly this).
• test_dimer_dangles_kwarg_exists_and_shifts_dg (+ subopt/tm/value-raises variants).

2. dangles scoping — took the cheap option: dangles is now explicitly MFE/suboptimal-scoped and dropped from the pfunc cache key (engine docstring spells it out). New test_pfunc_ignores_dangles and test_mfe_cache_key_tracks_dangles_but_pfunc_does_not pin the contract at both behavior and key level.

3. dangles=2 naming — kept the name but placed the divergence note at the API surface: the ValueError and the ThermoEngine doc now say "best single negative dangle per exterior stem; not identical to ViennaRNA dangles=2, which sums both flanks (two-flank cases deviate 0.2-0.9, 38/42 within 0.4)".

4. Dimer path — dimer_thermo/_subopt/dimer_tm accept dangles=0|2 threaded into _sum_dimer_elements (each terminus gathers adjacent negative stacks; convention documented as sum-per-terminus like the dimer DP, not the MFE best-single rule). Default 0 replaces the un-flagged always-count behavior; the changelog calls out the dG/Tm shift for flanked structures. Bonus latent-domain fix along the same flag: the enthalpy walk now always excludes dangle tables, so the ΔG-only values no longer leak into ΔH. 0.8 kcal/mol effect pinned in tests.

5. fraction_folded — the broad catch is gone; the zero case is decided from the MFE structure (no base pairs → 0.0). Bad salt model / degenerate entropy / multiloop MFE propagate. Deterministic raise-path test included.

6. Citation — re-confirmed at value level: the JSON matches dna_mathews2004.par (ATAT=-0.9 vs -0.8 in the 1999 file; GT/CG=+1.2 vs +1.3), excluding the 1999 set. THIRD_PARTY_NOTICES rewritten to record the marker check; citation stays the Mathews 1999 JMB lineage + NNDB Turner & Mathews 2010 (already cited).

Smaller ones:

• _resolve_name: the shared-instance contract is now documented as read-only; grep-verified the tree has zero in-place writes to a ParameterSet's tables, including the temperature-adjusted machinery.
• ext_e 5' branch: the redundant (ko - 1) not in nicks check dropped (3' branch symmetry). test_strand_start_flank_counted_in_complex pins it: on ['TACCT','AGATCTAGAACC'] the fix changes the selected structure to the more stable, strand-start-dangled fold (-1.57 vs -1.28 default).

Changelog — three entries under [Unreleased]: the fraction_folded behavior change, hairpin paramset + dimer dangles flag (incl. the dH-domain fix), dangles/pfunc scoping, the ext_e fix, and a documentation note for the value-level par confirmation.

Verification — full non-slow suite on the rebased tree: 669 passed, same 4 pre-existing TestNativeViennaParity environment failures as pristine main, no regressions. Parity battery (mathews2004 + native + dp + salt): 88/88.

[EmilioVenegas]

Thanks, the review round is thorough and I am happy with all six. Scoping dangles to MFE and dropping it from the pfunc key was the right pick of the two options, and settling the citation at value level (ATAT -0.9, GT/CG +1.2) is better than the argument I asked for. The DEF hairpin-mismatch fix in d5e8096d belongs here too, that is a real 0.5 kcal/mol error on affected hairpins.

#11 and #12 are both merged now, so this wants a rebase onto main.

On the 10k benchmark

The hairpin ΔG panel is the strongest result in it. strider_m sits near +0.2 with a tight body against strider_sl at +0.7 with a tail out to +5, which puts Mathews level with primer3 and Vienna. The Bland-Altman view explains why: the cyan wedge has a hard edge converging on the origin, which is strider_sl returning no structure where IDT finds a real hairpin. That is a structural miss, not noise, and the parameter set fixes it. Good argument for landing this PR's core.

Three things about the analysis before we lean on the rest of it:

Compare on the common subset. With Mathews at 26.5% NA and SantaLucia at 35.9%, and the plots filtered to finite values, the two methods are being scored on different sequence populations, and the harder cases drop out preferentially. Could you rerun restricted to oligos where all four methods return a value, with the NA rate reported separately as its own metric? As it stands, a change that turns hard cases into NAs reads as an accuracy gain.

The duplex Tm slope is probably not ours. Every method including primer3 drifts from about -2 at IDT Tm 50 to +5 at Tm 78, in near-parallel bands. Four independent implementations do not usually share a slope, so I would treat that as an IDT-side convention and chase only the per-method offsets (strider_m +1.3, strider_sl -2.3, primer3 near 0). Worth ruling out condition mismatch first: Na, Mg, dNTP, oligo concentration and the reporting temperature on both sides.

The hairpin Tm outliers are shared. strider_m reaches +250, but primer3 hits +80 and Vienna +120 on the same set. A failure mode common to all four points at the comparison rather than at strider, most likely oligos where a two-state Tm is not defined. I would pull those specific sequences and look at them before keeping a 200 °C ceiling, since the ceiling hides the cases instead of explaining them.

Also worth noting that homodimer ΔG is where strider is weakest, both variants around +2.2 to +2.6 while primer3 and Vienna sit near zero, and Mathews does not fix it. That is unaffected by this PR except that the dimer dangle default change and the ΔH domain fix will both move it, so the panel needs regenerating afterwards either way.

Request: split the hairpin commits out

70ec8e9c and eb18baf2 change hairpin_thermo defaults for every caller with no opt-out: the internal fold moves from 37 °C to 25 °C, multiloop MFEs are decomposed and the most stable stem-loop is scored, and near-degenerate ΔS or Tm above 200 °C now raise. Each is defensible, but together they change ΔG37, ΔH, ΔS and Tm for every hairpin call, and nothing in the suite pins those values tightly enough to notice, which is why CI stayed green.

The reason I want them separate: the Vallone hairpin set is measured Tm, and it is what established that strider is near-unbiased and that the beacon offset is probe-specific. Everything in the 10k benchmark measures agreement with another prediction tool. Tuning defaults toward IDT can move us away from the experimental agreement the validation rests on, and neither the repo suite nor this benchmark would show it.

So: keep the parameter set, the plumbing, dangles and the DEF fix here, and move the hairpin behavior changes to their own PR with the Vallone set run under old and new behavior side by side. If the new defaults hold up there, they land with evidence behind them.

Two smaller notes on that code for when it moves. The candidate loop in the multiloop path uses except Exception: continue, which is the pattern from item 5 reappearing in new code and will mask a real error behind "no valid stem-loop found in multiloop". And test_degenerate_ds_raises asserts MIN_DS_CAL == 0.5 and MAX_TM_CELSIUS == 200.0, which restates the constants rather than testing behavior.

One question rather than a request: is the 25 °C fold chosen because IDT reports hairpin analysis at 25 °C? If so that is a much better justification than the NA rate, and it argues for exposing it as a named parameter with a documented default so the 37 °C behavior stays reachable.

[KowalskiBio]

Split done. The two hairpin commits (70ec8e9 and eb18baf) have been removed from this PR (force-pushed). They are now in #14, stacked on this branch.

What stays in #13: parameter set, plumbing, dangles, DEF fix, and the pair_probs perf flag.

What moved to #14: 25 C fold, multiloop splitting, degenerate-dS/Tm guards, plus the two review notes fixed:

1. except Exception: continue -> except ValueError: continue in the multiloop candidate loop.
2. test_degenerate_ds_raises now patches the energy walkers to force near-zero dS and asserts ValueError is raised, instead of asserting MIN_DS_CAL == 0.5 and MAX_TM_CELSIUS == 200.0.

On the 25 C question: yes, 25 C is because IDT reports hairpin analysis at 25 C. FOLD_CELSIUS is now exposed as a fold_celsius parameter (default 25.0) on hairpin_thermo / hairpin_tm / fraction_folded, so 37 C stays reachable.

On the Vallone set: it is not in the repo. The 10k parity table is prediction-vs-prediction and cannot substitute. #14 notes this explicitly and recommends running the Vallone set under old and new behavior before it lands.