fix: correct dedx_get_inverse_stp() branch selection for monotone STP curves

[grzanka]

Closes #121 · Phase C of the dedx_extra.c migration epic #118.

Bug

dedx_get_inverse_stp() used find_min() to locate the energy of maximum stopping power, searching a hardcoded x ∈ [0.01, 10] instead of the program's actual tabulated energy range. Separately, the documented side contract is 0 = low-energy branch / 1 = high-energy branch, but the implementation only special-cased side < 0 for the low branch — so side == 0 and side == 1 were indistinguishable, both returning the descending-branch result.

In practice this made the function fail outright for ordinary queries. I reproduced it directly before fixing:

ICRU49 ion=1 stp=100 side=0 -> e=-1 err=203
ICRU49 ion=1 stp=100 side=1 -> e=-1 err=203
PSTAR  ion=1 stp=100 side=0 -> e=-1 err=203
PSTAR  ion=1 stp=100 side=1 -> e=-1 err=203

(err=203 = DEDX_ERR_INVALID_DATASET_ID / DEDX_ERR_ENERGY_OUT_OF_RANGE depending on the code path — either way, broken.)

Fix, round 1: sampled-peak + branch-aware bisection

Ported the sampled-peak + branch-aware bisection strategy proposed in #121 (originally proven out in dedx_web's wasm/dedx_extra.c dedx_get_inverse_stp_flat): sample STP at a fixed log-spaced density to locate the Bragg peak, then bisect the ascending or descending branch. This fixed the reported case but assumed a single interior peak.

Fix, round 2: STP(E) is not simply unimodal

Review (h/t @grzanka) pointed out that "Bragg peak" isn't the right term for the STP(E) maximum (that's a depth-dose concept — see the rename below), and — more importantly — questioned whether a single peak/branch model is even correct. It isn't. PSTAR/ICRU proton tables extend to 10 GeV/nucleon, far enough to show the curve rise to the Bragg peak (~0.08 MeV), fall through a minimum-ionizing point (~3000–4000 MeV), and rise again (relativistic rise / Fermi plateau) toward the table's max energy. I confirmed this directly against the tabulated data (both a dense synthetic scan and the raw knots) before changing anything: exactly two turning points for PSTAR/ICRU/ICRU49 proton/water, one turning point for every other program/ion combination I checked (their tables don't reach high enough energy to show the second one).

The previous fixed-sample-density scan assumed a single interior peak and rejected every STP in that high-energy dip as out of range, even though it's reachable at two energies:

stp=2.05 side=0 -> err=101 (before), e=1649 MeV (after)
stp=2.05 side=1 -> err=101 (before), e=5904 MeV (after)

Redesign: walk the exact tabulated knots instead of sampling

dedx_tools.c can reach dedx_lookup_data.h (a private, same-directory header) and read ws->loaded_data[config->cfg_id]->base[i].x / .a directly — these are the literal (energy, STP) pairs the spline was built from (dedx_spline.c: coef[i].a = stopping[i]). Replaced the fixed-sample-count scan with a single linear pass over these exact knots that finds every monotonic run the curve actually has — no sampling-density guesswork, no risk of straddling a feature between synthetic samples, and no arbitrary sample-count constants (DEDX_INVERSE_STP_SAMPLES / DEDX_MAX_STP_SAMPLES are both gone).

1. load_and_get_dataset() (new) loads the config if needed and returns the backing dedx_internal_lookup_data, bounds-checked like dedx_get_stp() already does.
2. dedx_get_inverse_stp() walks the knots once, brackets and bisects every monotonic run that contains the requested STP (via bisect_monotonic_run(), using dedx_get_stp()'s spline evaluation for the actual root-finding, so precision is unaffected), and returns the lowest-energy reachable solution for side == 0 or the highest-energy one for side == 1. This is a strict generalization of the old ascending/descending branch choice — identical behavior for the common single-peak case, but now also correctly resolves the two solutions on either side of the minimum-ionizing point.
3. dedx_get_max_stp() is now an exact scan of the same knots (previously a 300-sample approximation).
4. Added the complementary dedx_get_min_stp() — for tables reaching relativistic energies this is the minimum-ionizing point, not either tabulated endpoint.

Naming: dedx_get_max_stp() instead of dedx_get_bragg_peak_stp()

The issue spec (and the dedx_web reference this was ported from) named this dedx_get_bragg_peak_stp(). "Bragg peak" properly refers to the maximum of a depth-dose curve, which also depends on range straggling — not the maximum of the stopping-power-vs-energy curve this function returns. Renamed to dedx_get_max_stp() to name the quantity precisely; the header docstring calls out the distinction explicitly. wasm/dedx_extra.c is slated for removal once this migration completes, so it'll pick up the new name rather than the other way around. Commented on #121 with the same rationale.

Tests

tests/test_inverse_stp.c (new, auto-registered by the test_*.c glob):

• Below the Bragg peak's ascending floor: only the post-peak run reaches it, side=0/side=1 agree, round-tripped via dedx_get_stp().
• Around the Bragg peak: side=0/side=1 select distinct energies, both round-trip.
• The relativistic-rise dip: side=0/side=1 select the pre- and post-minimum energies respectively, both round-trip — the exact scenario that used to return DEDX_ERR_ENERGY_OUT_OF_RANGE.
• Out-of-range STP (above the global max, below the global min) returns DEDX_ERR_ENERGY_OUT_OF_RANGE, not a bogus/negative energy.
• ion_a <= 0 still rejected up front (unchanged contract).
• dedx_get_max_stp() / dedx_get_min_stp() checked against a 5000-point brute-force reference scan.

Verification

• ctest: 33/33 pass (32 existing + the new test).
• clang-format-19 --Werror: clean on all changed files.
• gcc -Wall -Wextra -Wpedantic -Wshadow: clean (clang-tidy itself isn't available in this environment; CI will run it).
• Cross-checked dedx_get_inverse_stp()/dedx_get_max_stp()/dedx_get_min_stp() against 7 other program/ion/material combinations (ICRU73/carbon, ASTAR/helium, MSTAR/carbon, Bethe, ICRU49/ICRU protons and helium) — all single-peak curves resolve correctly with distinct low/high solutions.

---

Authored via Claude Code as Phase C of the dedx_extra.c → libdedx migration plan.

[codecov]

Codecov Report

❌ Patch coverage is 96.49123% with 4 lines in your changes missing coverage. Please review.
✅ Project coverage is 77.99%. Comparing base (d5bf0cd) to head (0202a1b).

Files with missing lines	Patch %	Lines
src/dedx_tools.c	96.49%	4 Missing ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main     #148      +/-   ##
==========================================
+ Coverage   74.51%   77.99%   +3.47%     
==========================================
  Files          12       12              
  Lines        1711     1777      +66     
  Branches      317      340      +23     
==========================================
+ Hits         1275     1386     +111     
+ Misses        436      391      -45     

☔ View full report in Codecov by Harness.
📢 Have feedback on the report? Share it here.

🚀 New features to boost your workflow:

• ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.

[Copilot]

Pull request overview

This PR fixes dedx_get_inverse_stp() so it selects the correct inverse-branch over the actual tabulated energy range (including monotone STP curves), and adds a new public helper to query the Bragg-peak (maximum) stopping power for a loaded configuration.

Changes:

• Replace the previous hardcoded-range peak search with a log-sampled Bragg-peak finder and branch-aware bisection in dedx_get_inverse_stp().
• Add dedx_get_bragg_peak_stp() to the public API for retrieving the peak (maximum) stopping power over the tabulated range.
• Add a new regression test suite covering monotone fallback behavior, branch selection (side=0/1), and out-of-range STP handling.

Reviewed changes

Copilot reviewed 3 out of 3 changed files in this pull request and generated 1 comment.

File	Description
src/dedx_tools.c	Implements sampled Bragg-peak detection, fixes inverse STP branch selection/bounds handling, and adds dedx_get_bragg_peak_stp().
include/dedx_tools.h	Updates dedx_get_inverse_stp() docstring semantics and declares the new dedx_get_bragg_peak_stp() API.
tests/test_inverse_stp.c	Adds regression tests for inverse STP branching, monotone fallback, out-of-range behavior, and the new peak-STP helper.

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[Copilot]

Pull request overview

Copilot reviewed 3 out of 3 changed files in this pull request and generated 3 comments.

Comments suppressed due to low confidence (1)

src/dedx_tools.c:231

• Same workspace/config mismatch risk as in dedx_get_inverse_stp(): if config was loaded in a different workspace (or ws was recreated), config->loaded==1 can cause this function to skip dedx_load_config() and then all sampling calls into dedx_get_stp() will fail with DEDX_ERR_INVALID_DATASET_ID (src/dedx.c:490-492). Re-load when cfg_id is not valid for this workspace.

    if (config->loaded == 0)
        dedx_load_config(ws, config, err);

[Copilot]

Pull request overview

Copilot reviewed 3 out of 3 changed files in this pull request and generated 3 comments.

Comments suppressed due to low confidence (1)

include/dedx_tools.h:69

• Same as for dedx_get_inverse_stp(): the docstring states the configuration is already loaded, but dedx_get_max_stp() loads it when config->loaded == 0. Update the docs to reflect that the function may load/mutate config.

 *  @param[in]  ws      Workspace with a loaded configuration.
 *  @param[in]  config  Loaded configuration.

[Copilot]

Pull request overview

Copilot reviewed 3 out of 3 changed files in this pull request and generated 3 comments.

Comments suppressed due to low confidence (3)

src/dedx_tools.c:234

• side is documented as 0/1, but the previous implementation also treated side < 0 as selecting the low-energy branch. The new code maps any nonzero side (including negative) to the high-energy solution, which is a behavioral regression for existing callers that may have relied on side < 0 (and is inconsistent with the common “<=0 / >0” pattern). Consider treating side <= 0 as the low-energy solution and side > 0 as the high-energy solution.

    return (side == 0) ? x_min_found : x_max_found;

include/dedx_tools.h:64

• dedx_get_max_stp() is documented as the maximum of the stopping-power-vs-energy curve, but the implementation computes the maximum over tabulated knots only. With natural cubic spline interpolation, a true maximum can occur between knots. Either clarify that this returns the max knot value, or update the implementation to analyze spline extrema.

/** @brief Find the maximum stopping power over a program's whole tabulated
 *  energy range.
 *
 *  This is the peak of the stopping-power-vs-energy curve, not the Bragg
 *  peak of a depth-dose curve (which also depends on range straggling and is

include/dedx_tools.h:78

• Similarly, dedx_get_min_stp() returns the minimum over tabulated knots, but the doc implies it is the minimum of the interpolated curve. With natural cubic spline interpolation, the true minimum can occur between knots. Clarify the definition (knot min vs spline min) or update the implementation accordingly.

/** @brief Find the minimum stopping power over a program's whole tabulated
 *  energy range.
 *
 *  For tables that reach relativistic energies this is typically the
 *  minimum-ionizing point, not the value at either tabulated endpoint (the

[Copilot]

Pull request overview

Copilot reviewed 3 out of 3 changed files in this pull request and generated 1 comment.

[Copilot]

Pull request overview

Copilot reviewed 3 out of 3 changed files in this pull request and generated no new comments.