From 1662806e483a1aaa06069c8af000d468ebf1c788 Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Fri, 24 Jul 2026 12:38:06 +0200 Subject: [PATCH 1/6] fix: correct dedx_get_inverse_stp() branch selection for monotone STP curves (#121) find_min() searched a hardcoded x in [0.01, 10] instead of the program's actual tabulated energy range, and side < 0 was the only branch that ever selected the low-energy path, so side == 0 and side == 1 were indistinguishable. In practice dedx_get_inverse_stp() failed with DEDX_ERR_ENERGY_OUT_OF_RANGE for ordinary proton/water queries on both PSTAR and ICRU49, regardless of side. Replace find_min()/find_min_stp_func() with find_stp_peak(), which samples STP on a log-spaced grid over the real [emin, emax] to locate the Bragg-peak energy, then bisects the physically correct branch: the full descending range when there is no interior peak, side == 0 -> ascending [emin, e_peak] (falling back to the descending branch when the STP is below the ascending branch's floor), side == 1 -> descending [e_peak, emax]. Out-of-range STP values now return a clean error instead of looping to a bogus or negative energy. Also add dedx_get_bragg_peak_stp() per the issue's acceptance criteria, and fix dedx_get_inverse_stp() to guard on config->loaded before calling dedx_load_config(), matching dedx_get_csda()'s existing pattern -- without it, calling dedx_get_inverse_stp() on an already-loaded config corrupted the workspace. Co-Authored-By: Claude Sonnet 5 --- include/dedx_tools.h | 29 +++++-- src/dedx_tools.c | 177 ++++++++++++++++++++++++++++----------- tests/test_inverse_stp.c | 176 ++++++++++++++++++++++++++++++++++++++ 3 files changed, 326 insertions(+), 56 deletions(-) create mode 100644 tests/test_inverse_stp.c diff --git a/include/dedx_tools.h b/include/dedx_tools.h index 9f761f5..4825bc3 100644 --- a/include/dedx_tools.h +++ b/include/dedx_tools.h @@ -36,18 +36,37 @@ double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int /** @brief Find the energy corresponding to a given stopping power value. * - * Inverts the stopping power curve. Since stopping power is non-monotonic, - * the @p side parameter selects which branch to use. + * Inverts the stopping power curve. Stopping power rises from a program's + * minimum tabulated energy to a Bragg-peak energy, then falls off; when the + * requested @p stp corresponds to two energies, the @p side parameter + * selects which branch to use. If the curve is monotonically decreasing + * over the whole tabulated range (no interior peak), @p side is ignored and + * the single descending branch is used. * * @param[in] ws Workspace with a loaded configuration. * @param[in] config Loaded configuration. - * @param[in] stp Target stopping power in MeV cm²/g. - * @param[in] side 0 = low-energy branch, 1 = high-energy branch. + * @param[in] stp Target stopping power in MeV cm²/g. Must lie between + * the STP at the program's max energy and the + * Bragg-peak STP (inclusive); otherwise the value is + * unreachable and an error is returned. + * @param[in] side 0 = low-energy (ascending) branch, 1 = high-energy + * (descending) branch. Ignored when the curve has no + * interior Bragg peak. * @param[out] err Error code; 0 on success. - * @return Energy in MeV/nucl (MeV per nucleon). + * @return Energy in MeV/nucl (MeV per nucleon), or -1 on error. */ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, int side, int *err); +/** @brief Find the Bragg-peak stopping power: the maximum stopping power + * over a program's whole tabulated energy range. + * + * @param[in] ws Workspace with a loaded configuration. + * @param[in] config Loaded configuration. + * @param[out] err Error code; 0 on success. + * @return Maximum stopping power in MeV cm²/g, or -1 on error. + */ +double dedx_get_bragg_peak_stp(dedx_workspace *ws, dedx_config *config, int *err); + /** @brief Find the energy corresponding to a given CSDA range. * * @param[in] ws Workspace with a loaded configuration. diff --git a/src/dedx_tools.c b/src/dedx_tools.c index bb74e48..3e58f17 100644 --- a/src/dedx_tools.c +++ b/src/dedx_tools.c @@ -59,42 +59,53 @@ static double adapt_stp(double energy, dedx_tools_settings *set) { return 1.0 / stp; } -static double find_min_stp_func(double x, dedx_tools_settings *set) { - int err = 0; - double stp = dedx_get_stp(set->ws, set->cfg, x, &err); - if (err != 0 || stp == 0.0) - return INFINITY; - return 1.0 / stp; -} +/* Number of log-spaced samples used to locate the Bragg-peak energy before + * bisecting; matches the proven sampling density from the dedx_web reference + * implementation this was ported from (see #121). */ +#define DEDX_INVERSE_STP_SAMPLES 40 + +/* Number of log-spaced samples used by dedx_get_bragg_peak_stp() — denser + * than DEDX_INVERSE_STP_SAMPLES since it reports the peak value itself + * rather than just using it to pick a bisection branch (see #121). */ +#define DEDX_BRAGG_PEAK_SAMPLES 300 -static double find_min(double (*func)(double x, dedx_tools_settings *set), dedx_tools_settings *set, double acc) { - double x[] = {0.01, 10}; - double f[] = {func(x[0], set), func(x[1], set)}; - int i = 0; - double x_temp = 0; - double f_temp = 0; - double h = 0; - while (fabs(x[1] - x[0]) > acc) { - i = 0; - h = x[1] - x[0]; - if (f[1] > f[0]) - i = 1; - // try flip - - x_temp = x[i] + 2 * h * pow(-1, i); - f_temp = func(x_temp, set); - if (f_temp > f[i]) { - // try shrink - x_temp = x[i] + 0.5 * h * pow(-1, i); - f_temp = func(x_temp, set); - - if (f_temp > f[i]) - return -1; +/* Locate the Bragg-peak energy (the energy of maximum stopping power) within + * [emin, emax] by sampling on a log-spaced grid of n_samples points. Also + * reports the STP at the leftmost sample (emin) so callers can tell whether + * a requested STP lies on the ascending branch. Returns 0 on success, -1 if + * no sample succeeded. */ +static int find_stp_peak(dedx_workspace *ws, + dedx_config *config, + double emin, + double emax, + int n_samples, + double *e_peak, + double *max_stp, + double *stp_at_emin) { + double log_emin = log(emin); + double log_emax = log(emax); + double log_step = (log_emax - log_emin) / (n_samples - 1); + int have_sample = 0; + + *max_stp = 0.0; + *e_peak = emin; + *stp_at_emin = 0.0; + + for (int i = 0; i < n_samples; i++) { + double e = exp(log_emin + i * log_step); + int stp_err = 0; + double s = dedx_get_stp(ws, config, (float) e, &stp_err); + if (stp_err != 0) + continue; + have_sample = 1; + if (i == 0) + *stp_at_emin = s; + if (s > *max_stp) { + *max_stp = s; + *e_peak = e; } - f[i] = f_temp; - x[i] = x_temp; } - return (x[0] + x[1]) / 2; + return have_sample ? 0 : -1; } double dedx_get_inverse_csda(dedx_workspace *ws, dedx_config *config, float range, int *err) { @@ -129,41 +140,105 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, *err = DEDX_ERR_ION_A_REQUIRED; return -1; } - double acc = 1e-5; - dedx_tools_settings set; - set.ws = ws; if (*err != 0) return -1; - dedx_load_config(ws, config, err); - set.cfg = config; - + if (config->loaded == 0) + dedx_load_config(ws, config, err); if (*err != 0) return -1; - double max = find_min(find_min_stp_func, &set, acc * 100); + + double acc = 1e-5; + double emin = dedx_get_min_energy(config->program, config->ion); + double emax = dedx_get_max_energy(config->program, config->ion); + + /* Sample the curve to find the Bragg-peak energy, then bisect the + * physically correct monotone branch: + * - no interior peak (monotone descending over the full range): + * bisect [emin, emax] on the single descending branch. + * - interior peak: side == 0 selects the low/ascending branch + * [emin, e_peak]; side == 1 (or any stp below the ascending + * branch's floor at emin) selects the high/descending branch + * [e_peak, emax]. + * See #121 for why the previous find_min()-based approach failed. */ + double e_peak; + double max_stp; + double stp_at_emin; + if (find_stp_peak(ws, config, emin, emax, DEDX_INVERSE_STP_SAMPLES, &e_peak, &max_stp, &stp_at_emin) != 0) { + *err = DEDX_ERR_ENERGY_OUT_OF_RANGE; + return -1; + } + + int stp_err = 0; + double stp_at_emax = dedx_get_stp(ws, config, (float) emax, &stp_err); + if (stp_err != 0 || max_stp == 0.0 || stp > max_stp || stp < stp_at_emax) { + *err = (stp_err != 0) ? stp_err : DEDX_ERR_ENERGY_OUT_OF_RANGE; + return -1; + } + + double log_step = (log(emax) - log(emin)) / (DEDX_INVERSE_STP_SAMPLES - 1); + int has_peak = e_peak > emin * exp(log_step); + double x1; double x2; - double x_temp; - double f_temp; - if (side < 0) { - x1 = dedx_get_min_energy(config->program, config->ion); - x2 = max; + int ascending; + if (!has_peak) { + x1 = emin; + x2 = emax; + ascending = 0; + } else if (side == 0 && stp >= stp_at_emin) { + x1 = emin; + x2 = e_peak; + ascending = 1; } else { - x2 = max; - x1 = dedx_get_max_energy(config->program, config->ion); + x1 = e_peak; + x2 = emax; + ascending = 0; } + + double x_temp; + double f_temp; while (fabs(x1 - x2) > acc) { x_temp = (x1 + x2) / 2; - f_temp = dedx_get_stp(set.ws, set.cfg, x_temp, err); + f_temp = dedx_get_stp(ws, config, (float) x_temp, err); + if (*err != 0) + return -1; - if (f_temp >= stp) { - x2 = x_temp; + if (ascending) { + if (f_temp <= stp) + x1 = x_temp; + else + x2 = x_temp; } else { - x1 = x_temp; + if (f_temp >= stp) + x1 = x_temp; + else + x2 = x_temp; } } return (x1 + x2) / 2; } +double dedx_get_bragg_peak_stp(dedx_workspace *ws, dedx_config *config, int *err) { + if (*err != 0) + return -1; + if (config->loaded == 0) + dedx_load_config(ws, config, err); + if (*err != 0) + return -1; + + double emin = dedx_get_min_energy(config->program, config->ion); + double emax = dedx_get_max_energy(config->program, config->ion); + + double e_peak; + double max_stp; + double stp_at_emin; + if (find_stp_peak(ws, config, emin, emax, DEDX_BRAGG_PEAK_SAMPLES, &e_peak, &max_stp, &stp_at_emin) != 0) { + *err = DEDX_ERR_ENERGY_OUT_OF_RANGE; + return -1; + } + return max_stp; +} + double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int *err) { if (config->ion_a <= 0) { *err = DEDX_ERR_ION_A_REQUIRED; diff --git a/tests/test_inverse_stp.c b/tests/test_inverse_stp.c new file mode 100644 index 0000000..28b3eda --- /dev/null +++ b/tests/test_inverse_stp.c @@ -0,0 +1,176 @@ +#include +#include +#include +#include +#include +#include + +/* Regression tests for dedx_get_inverse_stp() and dedx_get_bragg_peak_stp() + * (issue #121). The previous implementation used find_min() over a hardcoded + * x in [0.01, 10] to locate the Bragg peak, which does not match the real + * tabulated energy range; it also mapped side < 0 (never 0 or 1) to the + * low-energy branch, so side == 0 and side == 1 were indistinguishable. + * Both bugs made dedx_get_inverse_stp() fail with DEDX_ERR_ENERGY_OUT_OF_RANGE + * for ordinary proton/water queries. + */ + +static int failures = 0; + +static void expect_int(const char *label, long got, long expected) { + if (got != expected) { + fprintf(stderr, "FAIL %s: got %ld expected %ld\n", label, got, expected); + failures++; + } +} + +static void expect_near(const char *label, double got, double expected, double rel) { + double scale = fabs(expected) > 1e-12 ? fabs(expected) : 1.0; + if (fabs(got - expected) > rel * scale) { + fprintf(stderr, "FAIL %s: got %.8g expected %.8g\n", label, got, expected); + failures++; + } +} + +static void expect_true(const char *label, int condition) { + if (!condition) { + fprintf(stderr, "FAIL %s\n", label); + failures++; + } +} + +static dedx_config *make_config(int program, int ion, int ion_a, int target) { + dedx_config *cfg = calloc(1, sizeof(dedx_config)); + cfg->program = program; + cfg->ion = ion; + cfg->target = target; + cfg->ion_a = ion_a; + return cfg; +} + +int main(void) { + int err = 0; + + /* --- Below the ascending branch's floor: only the descending branch can + * reach this STP, exactly the scenario that made the old find_min()-based + * code return DEDX_ERR_ENERGY_OUT_OF_RANGE for proton/water. Both sides + * must now fall back to the same (only reachable) descending-branch + * energy, with no error and a positive result. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + expect_int("workspace alloc err", err, DEDX_OK); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + err = 0; + double e0 = dedx_get_inverse_stp(ws, cfg, 100.0f, 0, &err); + expect_int("low-floor stp side=0 err", err, DEDX_OK); + expect_true("low-floor stp side=0 positive", e0 > 0.0); + + err = 0; + double e1 = dedx_get_inverse_stp(ws, cfg, 100.0f, 1, &err); + expect_int("low-floor stp side=1 err", err, DEDX_OK); + expect_near("low-floor stp side=0/1 agree", e0, e1, 1e-3); + + int verify_err = 0; + double stp_at_e0 = dedx_get_stp(ws, cfg, (float) e0, &verify_err); + expect_int("low-floor roundtrip err", verify_err, DEDX_OK); + expect_near("low-floor roundtrip stp", stp_at_e0, 100.0, 1e-2); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- Interior Bragg peak, STP reachable on both branches: side=0 (low + * energy / ascending) and side=1 (high energy / descending) must select + * distinct branches and each round-trip back to the requested STP. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + const float target_stp = 500.0f; + + err = 0; + double e_low = dedx_get_inverse_stp(ws, cfg, target_stp, 0, &err); + expect_int("branch stp side=0 err", err, DEDX_OK); + + err = 0; + double e_high = dedx_get_inverse_stp(ws, cfg, target_stp, 1, &err); + expect_int("branch stp side=1 err", err, DEDX_OK); + + expect_true("side=0 selects the lower energy", e_low < e_high); + + int verify_err = 0; + double stp_low = dedx_get_stp(ws, cfg, (float) e_low, &verify_err); + expect_int("side=0 roundtrip err", verify_err, DEDX_OK); + expect_near("side=0 roundtrip stp", stp_low, target_stp, 1e-2); + + verify_err = 0; + double stp_high = dedx_get_stp(ws, cfg, (float) e_high, &verify_err); + expect_int("side=1 roundtrip err", verify_err, DEDX_OK); + expect_near("side=1 roundtrip stp", stp_high, target_stp, 1e-2); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- Out-of-range STP values (above the Bragg peak, or below the STP at + * max energy) must return an error rather than looping to a bogus or + * negative energy. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + + err = 0; + double e_too_high = dedx_get_inverse_stp(ws, cfg, 2000.0f, 0, &err); + expect_int("stp above peak err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); + expect_near("stp above peak sentinel", e_too_high, -1.0, 1e-9); + + err = 0; + double e_too_low = dedx_get_inverse_stp(ws, cfg, 0.5f, 0, &err); + expect_int("stp below max-energy floor err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); + expect_near("stp below floor sentinel", e_too_low, -1.0, 1e-9); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- ion_a <= 0 is rejected up front, same contract as before. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 0, DEDX_WATER); + err = 0; + double e = dedx_get_inverse_stp(ws, cfg, 100.0f, 0, &err); + expect_int("ion_a required err", err, DEDX_ERR_ION_A_REQUIRED); + expect_near("ion_a required sentinel", e, -1.0, 1e-9); + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- dedx_get_bragg_peak_stp() reports the true maximum of the curve. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + err = 0; + double peak = dedx_get_bragg_peak_stp(ws, cfg, &err); + expect_int("bragg peak err", err, DEDX_OK); + + /* Brute-force reference: sample far more densely and take the max. */ + float emin = dedx_get_min_energy(cfg->program, cfg->ion); + float emax = dedx_get_max_energy(cfg->program, cfg->ion); + double log_emin = log((double) emin); + double log_emax = log((double) emax); + double reference_peak = 0.0; + const int n = 2000; + for (int i = 0; i < n; i++) { + float e = (float) exp(log_emin + (log_emax - log_emin) * i / (n - 1)); + int sample_err = 0; + double s = dedx_get_stp(ws, cfg, e, &sample_err); + if (sample_err == 0 && s > reference_peak) + reference_peak = s; + } + expect_near("bragg peak matches brute-force reference", peak, reference_peak, 5e-2); + expect_true("bragg peak is a real maximum", peak > 0.0); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + return failures; +} From f91977aa233712f0cff7942551fe31c11b0fc477 Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Fri, 24 Jul 2026 12:48:08 +0200 Subject: [PATCH 2/6] fix: use emin directly for the first Bragg-peak sample exp(log(emin)) is not guaranteed to round-trip to exactly emin; a reconstructed value fractionally below the dataset's real lower bound would be rejected by dedx_get_stp(), leaving stp_at_emin stuck at 0 and making the ascending-branch check in dedx_get_inverse_stp() always pass. Addresses Copilot review comment on #148. Co-Authored-By: Claude Sonnet 5 --- src/dedx_tools.c | 8 +++++++- 1 file changed, 7 insertions(+), 1 deletion(-) diff --git a/src/dedx_tools.c b/src/dedx_tools.c index 3e58f17..0c86625 100644 --- a/src/dedx_tools.c +++ b/src/dedx_tools.c @@ -92,7 +92,13 @@ static int find_stp_peak(dedx_workspace *ws, *stp_at_emin = 0.0; for (int i = 0; i < n_samples; i++) { - double e = exp(log_emin + i * log_step); + /* Use emin itself for the first sample: exp(log(emin)) is not + * guaranteed to round-trip to exactly emin, and a reconstructed + * value fractionally below the dataset's real lower bound would be + * rejected by dedx_get_stp(), leaving stp_at_emin stuck at 0 and + * making the ascending-branch check in dedx_get_inverse_stp() + * always pass. */ + double e = (i == 0) ? emin : exp(log_emin + i * log_step); int stp_err = 0; double s = dedx_get_stp(ws, config, (float) e, &stp_err); if (stp_err != 0) From 27b63b27e1868530e61f319d70300e1840bd6973 Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Fri, 24 Jul 2026 16:18:07 +0200 Subject: [PATCH 3/6] rename: dedx_get_bragg_peak_stp() -> dedx_get_max_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. dedx_get_max_stp() names the quantity precisely instead of borrowing dosimetry terminology for a stopping-power quantity. This diverges from the name used in the #121 issue spec and in dedx_web's wasm/dedx_extra.c reference implementation; dedx_extra.c is slated for removal once the migration to libdedx completes, so it will pick up the new name rather than the other way around. Co-Authored-By: Claude Sonnet 5 --- include/dedx_tools.h | 29 +++++++++++++++++------------ src/dedx_tools.c | 32 ++++++++++++++++---------------- tests/test_inverse_stp.c | 35 ++++++++++++++++++----------------- 3 files changed, 51 insertions(+), 45 deletions(-) diff --git a/include/dedx_tools.h b/include/dedx_tools.h index 4825bc3..ee03ffa 100644 --- a/include/dedx_tools.h +++ b/include/dedx_tools.h @@ -37,35 +37,40 @@ double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int /** @brief Find the energy corresponding to a given stopping power value. * * Inverts the stopping power curve. Stopping power rises from a program's - * minimum tabulated energy to a Bragg-peak energy, then falls off; when the - * requested @p stp corresponds to two energies, the @p side parameter - * selects which branch to use. If the curve is monotonically decreasing - * over the whole tabulated range (no interior peak), @p side is ignored and - * the single descending branch is used. + * minimum tabulated energy to the energy of maximum stopping power, then + * falls off; when the requested @p stp corresponds to two energies, the + * @p side parameter selects which branch to use. If the curve is + * monotonically decreasing over the whole tabulated range (no interior + * peak), @p side is ignored and the single descending branch is used. * * @param[in] ws Workspace with a loaded configuration. * @param[in] config Loaded configuration. * @param[in] stp Target stopping power in MeV cm²/g. Must lie between - * the STP at the program's max energy and the - * Bragg-peak STP (inclusive); otherwise the value is - * unreachable and an error is returned. + * the STP at the program's max energy and the maximum + * stopping power (inclusive, see dedx_get_max_stp()); + * otherwise the value is unreachable and an error is + * returned. * @param[in] side 0 = low-energy (ascending) branch, 1 = high-energy * (descending) branch. Ignored when the curve has no - * interior Bragg peak. + * interior peak. * @param[out] err Error code; 0 on success. * @return Energy in MeV/nucl (MeV per nucleon), or -1 on error. */ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, int side, int *err); -/** @brief Find the Bragg-peak stopping power: the maximum stopping power - * over a program's whole tabulated energy range. +/** @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 + * not computed by this library). * * @param[in] ws Workspace with a loaded configuration. * @param[in] config Loaded configuration. * @param[out] err Error code; 0 on success. * @return Maximum stopping power in MeV cm²/g, or -1 on error. */ -double dedx_get_bragg_peak_stp(dedx_workspace *ws, dedx_config *config, int *err); +double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err); /** @brief Find the energy corresponding to a given CSDA range. * diff --git a/src/dedx_tools.c b/src/dedx_tools.c index 0c86625..315f7a2 100644 --- a/src/dedx_tools.c +++ b/src/dedx_tools.c @@ -59,21 +59,21 @@ static double adapt_stp(double energy, dedx_tools_settings *set) { return 1.0 / stp; } -/* Number of log-spaced samples used to locate the Bragg-peak energy before - * bisecting; matches the proven sampling density from the dedx_web reference - * implementation this was ported from (see #121). */ +/* Number of log-spaced samples used to locate the energy of maximum stopping + * power before bisecting; matches the proven sampling density from the + * dedx_web reference implementation this was ported from (see #121). */ #define DEDX_INVERSE_STP_SAMPLES 40 -/* Number of log-spaced samples used by dedx_get_bragg_peak_stp() — denser - * than DEDX_INVERSE_STP_SAMPLES since it reports the peak value itself - * rather than just using it to pick a bisection branch (see #121). */ -#define DEDX_BRAGG_PEAK_SAMPLES 300 +/* Number of log-spaced samples used by dedx_get_max_stp() — denser than + * DEDX_INVERSE_STP_SAMPLES since it reports the peak value itself rather + * than just using it to pick a bisection branch (see #121). */ +#define DEDX_MAX_STP_SAMPLES 300 -/* Locate the Bragg-peak energy (the energy of maximum stopping power) within - * [emin, emax] by sampling on a log-spaced grid of n_samples points. Also - * reports the STP at the leftmost sample (emin) so callers can tell whether - * a requested STP lies on the ascending branch. Returns 0 on success, -1 if - * no sample succeeded. */ +/* Locate the energy of maximum stopping power within [emin, emax] by + * sampling on a log-spaced grid of n_samples points. Also reports the STP at + * the leftmost sample (emin) so callers can tell whether a requested STP + * lies on the ascending branch. Returns 0 on success, -1 if no sample + * succeeded. */ static int find_stp_peak(dedx_workspace *ws, dedx_config *config, double emin, @@ -157,8 +157,8 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, double emin = dedx_get_min_energy(config->program, config->ion); double emax = dedx_get_max_energy(config->program, config->ion); - /* Sample the curve to find the Bragg-peak energy, then bisect the - * physically correct monotone branch: + /* Sample the curve to find the energy of maximum stopping power, then + * bisect the physically correct monotone branch: * - no interior peak (monotone descending over the full range): * bisect [emin, emax] on the single descending branch. * - interior peak: side == 0 selects the low/ascending branch @@ -224,7 +224,7 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, return (x1 + x2) / 2; } -double dedx_get_bragg_peak_stp(dedx_workspace *ws, dedx_config *config, int *err) { +double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err) { if (*err != 0) return -1; if (config->loaded == 0) @@ -238,7 +238,7 @@ double dedx_get_bragg_peak_stp(dedx_workspace *ws, dedx_config *config, int *err double e_peak; double max_stp; double stp_at_emin; - if (find_stp_peak(ws, config, emin, emax, DEDX_BRAGG_PEAK_SAMPLES, &e_peak, &max_stp, &stp_at_emin) != 0) { + if (find_stp_peak(ws, config, emin, emax, DEDX_MAX_STP_SAMPLES, &e_peak, &max_stp, &stp_at_emin) != 0) { *err = DEDX_ERR_ENERGY_OUT_OF_RANGE; return -1; } diff --git a/tests/test_inverse_stp.c b/tests/test_inverse_stp.c index 28b3eda..f799d5e 100644 --- a/tests/test_inverse_stp.c +++ b/tests/test_inverse_stp.c @@ -5,13 +5,13 @@ #include #include -/* Regression tests for dedx_get_inverse_stp() and dedx_get_bragg_peak_stp() +/* Regression tests for dedx_get_inverse_stp() and dedx_get_max_stp() * (issue #121). The previous implementation used find_min() over a hardcoded - * x in [0.01, 10] to locate the Bragg peak, which does not match the real - * tabulated energy range; it also mapped side < 0 (never 0 or 1) to the - * low-energy branch, so side == 0 and side == 1 were indistinguishable. - * Both bugs made dedx_get_inverse_stp() fail with DEDX_ERR_ENERGY_OUT_OF_RANGE - * for ordinary proton/water queries. + * x in [0.01, 10] to locate the energy of maximum stopping power, which does + * not match the real tabulated energy range; it also mapped side < 0 (never + * 0 or 1) to the low-energy branch, so side == 0 and side == 1 were + * indistinguishable. Both bugs made dedx_get_inverse_stp() fail with + * DEDX_ERR_ENERGY_OUT_OF_RANGE for ordinary proton/water queries. */ static int failures = 0; @@ -78,9 +78,10 @@ int main(void) { dedx_free_workspace(ws, &err); } - /* --- Interior Bragg peak, STP reachable on both branches: side=0 (low - * energy / ascending) and side=1 (high energy / descending) must select - * distinct branches and each round-trip back to the requested STP. */ + /* --- Interior peak in the STP curve, STP reachable on both branches: + * side=0 (low energy / ascending) and side=1 (high energy / descending) + * must select distinct branches and each round-trip back to the + * requested STP. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err); dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); @@ -110,9 +111,9 @@ int main(void) { dedx_free_workspace(ws, &err); } - /* --- Out-of-range STP values (above the Bragg peak, or below the STP at - * max energy) must return an error rather than looping to a bogus or - * negative energy. */ + /* --- Out-of-range STP values (above the maximum stopping power, or + * below the STP at max energy) must return an error rather than looping + * to a bogus or negative energy. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err); dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); @@ -143,13 +144,13 @@ int main(void) { dedx_free_workspace(ws, &err); } - /* --- dedx_get_bragg_peak_stp() reports the true maximum of the curve. */ + /* --- dedx_get_max_stp() reports the true maximum of the curve. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err); dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); err = 0; - double peak = dedx_get_bragg_peak_stp(ws, cfg, &err); - expect_int("bragg peak err", err, DEDX_OK); + double peak = dedx_get_max_stp(ws, cfg, &err); + expect_int("max stp err", err, DEDX_OK); /* Brute-force reference: sample far more densely and take the max. */ float emin = dedx_get_min_energy(cfg->program, cfg->ion); @@ -165,8 +166,8 @@ int main(void) { if (sample_err == 0 && s > reference_peak) reference_peak = s; } - expect_near("bragg peak matches brute-force reference", peak, reference_peak, 5e-2); - expect_true("bragg peak is a real maximum", peak > 0.0); + expect_near("max stp matches brute-force reference", peak, reference_peak, 5e-2); + expect_true("max stp is a real maximum", peak > 0.0); dedx_free_config(cfg, &err); dedx_free_workspace(ws, &err); From 31e6e9aa610b3b5113b00393f5003f3c14ce96e6 Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Fri, 24 Jul 2026 18:54:19 +0200 Subject: [PATCH 4/6] redesign: walk exact tabulated knots instead of sampling at a fixed density STP(E) is not simply unimodal. 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) toward the table's max energy. The previous 40/300-sample log-spaced 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 -- verified directly against the data: stp=2.05 side=0 -> err=101 (before), e=1649 (after) stp=2.05 side=1 -> err=101 (before), e=5904 (after) Replace the sampling with a linear scan over the exact tabulated knots (dedx_tools.c can reach dedx_lookup_data.h since it's a private, same-directory header; dedx_spline.c sets coef[i].a = stopping[i], so the knots are the literal (energy, STP) pairs the spline was built from). The scan finds every monotonic run the curve actually has in one pass, no sampling-density guesswork, and dedx_get_inverse_stp() now brackets and bisects each run that contains the requested STP, returning the lowest-energy (side=0) or highest-energy (side=1) reachable solution -- a strict generalization of the old ascending/descending choice that also covers curves with more than one turning point. dedx_get_max_stp() becomes an exact scan of the same knots (previously a 300-sample approximation). Add the complementary dedx_get_min_stp(), which for tables reaching relativistic energies reports the minimum-ionizing point rather than either tabulated endpoint. tests/test_inverse_stp.c: add a regression test for the relativistic-rise dip (the exact scenario above) and for dedx_get_min_stp(). Co-Authored-By: Claude Sonnet 5 --- include/dedx_tools.h | 39 ++++-- src/dedx_tools.c | 259 ++++++++++++++++++++------------------- tests/test_inverse_stp.c | 121 +++++++++++++----- 3 files changed, 248 insertions(+), 171 deletions(-) diff --git a/include/dedx_tools.h b/include/dedx_tools.h index ee03ffa..3d8d822 100644 --- a/include/dedx_tools.h +++ b/include/dedx_tools.h @@ -36,23 +36,22 @@ double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int /** @brief Find the energy corresponding to a given stopping power value. * - * Inverts the stopping power curve. Stopping power rises from a program's - * minimum tabulated energy to the energy of maximum stopping power, then - * falls off; when the requested @p stp corresponds to two energies, the - * @p side parameter selects which branch to use. If the curve is - * monotonically decreasing over the whole tabulated range (no interior - * peak), @p side is ignored and the single descending branch is used. + * Inverts the stopping power curve. Stopping power is not simply unimodal: + * real tables can rise to a maximum (the Bragg peak) at low energy, fall + * through the minimum-ionizing point, and rise again at relativistic + * energies (e.g. proton tables extending to several GeV), so a requested + * @p stp can be reachable at more than one energy. Among all reachable + * energies, @p side selects which one to return: 0 = the lowest, 1 = the + * highest. When only one energy reaches @p stp, @p side has no effect. * * @param[in] ws Workspace with a loaded configuration. * @param[in] config Loaded configuration. * @param[in] stp Target stopping power in MeV cm²/g. Must lie between - * the STP at the program's max energy and the maximum - * stopping power (inclusive, see dedx_get_max_stp()); - * otherwise the value is unreachable and an error is - * returned. - * @param[in] side 0 = low-energy (ascending) branch, 1 = high-energy - * (descending) branch. Ignored when the curve has no - * interior peak. + * dedx_get_min_stp() and dedx_get_max_stp() (inclusive) + * for this configuration; otherwise the value is + * unreachable and an error is returned. + * @param[in] side 0 = lowest-energy solution, 1 = highest-energy + * solution. * @param[out] err Error code; 0 on success. * @return Energy in MeV/nucl (MeV per nucleon), or -1 on error. */ @@ -72,6 +71,20 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, */ double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err); +/** @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 + * curve can rise again past it — see dedx_get_inverse_stp()). + * + * @param[in] ws Workspace with a loaded configuration. + * @param[in] config Loaded configuration. + * @param[out] err Error code; 0 on success. + * @return Minimum stopping power in MeV cm²/g, or -1 on error. + */ +double dedx_get_min_stp(dedx_workspace *ws, dedx_config *config, int *err); + /** @brief Find the energy corresponding to a given CSDA range. * * @param[in] ws Workspace with a loaded configuration. diff --git a/src/dedx_tools.c b/src/dedx_tools.c index 315f7a2..4213aed 100644 --- a/src/dedx_tools.c +++ b/src/dedx_tools.c @@ -5,6 +5,7 @@ #include #include "dedx_data_access.h" +#include "dedx_lookup_data.h" typedef struct { dedx_workspace *ws; @@ -59,59 +60,39 @@ static double adapt_stp(double energy, dedx_tools_settings *set) { return 1.0 / stp; } -/* Number of log-spaced samples used to locate the energy of maximum stopping - * power before bisecting; matches the proven sampling density from the - * dedx_web reference implementation this was ported from (see #121). */ -#define DEDX_INVERSE_STP_SAMPLES 40 - -/* Number of log-spaced samples used by dedx_get_max_stp() — denser than - * DEDX_INVERSE_STP_SAMPLES since it reports the peak value itself rather - * than just using it to pick a bisection branch (see #121). */ -#define DEDX_MAX_STP_SAMPLES 300 - -/* Locate the energy of maximum stopping power within [emin, emax] by - * sampling on a log-spaced grid of n_samples points. Also reports the STP at - * the leftmost sample (emin) so callers can tell whether a requested STP - * lies on the ascending branch. Returns 0 on success, -1 if no sample - * succeeded. */ -static int find_stp_peak(dedx_workspace *ws, - dedx_config *config, - double emin, - double emax, - int n_samples, - double *e_peak, - double *max_stp, - double *stp_at_emin) { - double log_emin = log(emin); - double log_emax = log(emax); - double log_step = (log_emax - log_emin) / (n_samples - 1); - int have_sample = 0; - - *max_stp = 0.0; - *e_peak = emin; - *stp_at_emin = 0.0; +/* Load (if needed) and return the internal dataset backing a config, so + * callers can walk its exact tabulated (energy, STP) knots directly instead + * of re-sampling the curve at an arbitrary density. dedx_get_stp() evaluates + * a spline built from precisely these knots (dedx_spline.c: coef[i].a = + * stopping[i]), so scanning them finds every feature the public curve can + * actually exhibit -- no synthetic sampling can do better or worse. */ +static dedx_internal_lookup_data *load_and_get_dataset(dedx_workspace *ws, dedx_config *config, int *err) { + if (*err != 0) + return NULL; + if (config->loaded == 0) + dedx_load_config(ws, config, err); + if (*err != 0) + return NULL; + int id = config->cfg_id; + if (id < 0 || id >= ws->active_datasets) { + *err = DEDX_ERR_INVALID_DATASET_ID; + return NULL; + } + return ws->loaded_data[id]; +} - for (int i = 0; i < n_samples; i++) { - /* Use emin itself for the first sample: exp(log(emin)) is not - * guaranteed to round-trip to exactly emin, and a reconstructed - * value fractionally below the dataset's real lower bound would be - * rejected by dedx_get_stp(), leaving stp_at_emin stuck at 0 and - * making the ascending-branch check in dedx_get_inverse_stp() - * always pass. */ - double e = (i == 0) ? emin : exp(log_emin + i * log_step); - int stp_err = 0; - double s = dedx_get_stp(ws, config, (float) e, &stp_err); - if (stp_err != 0) - continue; - have_sample = 1; - if (i == 0) - *stp_at_emin = s; - if (s > *max_stp) { - *max_stp = s; - *e_peak = e; - } +static void min_max_stp_over_table(const dedx_internal_lookup_data *data, double *min_stp, double *max_stp) { + double lo = data->base[0].a; + double hi = data->base[0].a; + for (int i = 1; i < data->n; i++) { + double v = data->base[i].a; + if (v < lo) + lo = v; + if (v > hi) + hi = v; } - return have_sample ? 0 : -1; + *min_stp = lo; + *max_stp = hi; } double dedx_get_inverse_csda(dedx_workspace *ws, dedx_config *config, float range, int *err) { @@ -141,74 +122,37 @@ double dedx_get_inverse_csda(dedx_workspace *ws, dedx_config *config, float rang return (min + max) / 2; } -double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, int side, int *err) { - if (config->ion_a <= 0) { - *err = DEDX_ERR_ION_A_REQUIRED; - return -1; - } - if (*err != 0) - return -1; - if (config->loaded == 0) - dedx_load_config(ws, config, err); - if (*err != 0) +/* Bisect the monotonic run of knots [lo_idx, hi_idx] for the energy where + * the spline equals stp, refining via dedx_get_stp() (not linear knot + * interpolation) so the result matches what callers would measure with the + * public curve. Returns 0 and sets *solution on success; returns -1 without + * touching *err if stp is not bracketed by this run. */ +static int bisect_monotonic_run(dedx_workspace *ws, + dedx_config *config, + const dedx_internal_lookup_data *data, + int lo_idx, + int hi_idx, + float stp, + int *err, + double *solution) { + double lo_x = data->base[lo_idx].x; + double hi_x = data->base[hi_idx].x; + double lo_v = data->base[lo_idx].a; + double hi_v = data->base[hi_idx].a; + double range_min = lo_v < hi_v ? lo_v : hi_v; + double range_max = lo_v > hi_v ? lo_v : hi_v; + if (stp < range_min || stp > range_max) return -1; + int ascending = hi_v >= lo_v; + double x1 = lo_x; + double x2 = hi_x; double acc = 1e-5; - double emin = dedx_get_min_energy(config->program, config->ion); - double emax = dedx_get_max_energy(config->program, config->ion); - - /* Sample the curve to find the energy of maximum stopping power, then - * bisect the physically correct monotone branch: - * - no interior peak (monotone descending over the full range): - * bisect [emin, emax] on the single descending branch. - * - interior peak: side == 0 selects the low/ascending branch - * [emin, e_peak]; side == 1 (or any stp below the ascending - * branch's floor at emin) selects the high/descending branch - * [e_peak, emax]. - * See #121 for why the previous find_min()-based approach failed. */ - double e_peak; - double max_stp; - double stp_at_emin; - if (find_stp_peak(ws, config, emin, emax, DEDX_INVERSE_STP_SAMPLES, &e_peak, &max_stp, &stp_at_emin) != 0) { - *err = DEDX_ERR_ENERGY_OUT_OF_RANGE; - return -1; - } - - int stp_err = 0; - double stp_at_emax = dedx_get_stp(ws, config, (float) emax, &stp_err); - if (stp_err != 0 || max_stp == 0.0 || stp > max_stp || stp < stp_at_emax) { - *err = (stp_err != 0) ? stp_err : DEDX_ERR_ENERGY_OUT_OF_RANGE; - return -1; - } - - double log_step = (log(emax) - log(emin)) / (DEDX_INVERSE_STP_SAMPLES - 1); - int has_peak = e_peak > emin * exp(log_step); - - double x1; - double x2; - int ascending; - if (!has_peak) { - x1 = emin; - x2 = emax; - ascending = 0; - } else if (side == 0 && stp >= stp_at_emin) { - x1 = emin; - x2 = e_peak; - ascending = 1; - } else { - x1 = e_peak; - x2 = emax; - ascending = 0; - } - - double x_temp; - double f_temp; while (fabs(x1 - x2) > acc) { - x_temp = (x1 + x2) / 2; - f_temp = dedx_get_stp(ws, config, (float) x_temp, err); + double x_temp = (x1 + x2) / 2; + double f_temp = dedx_get_stp(ws, config, (float) x_temp, err); if (*err != 0) return -1; - if (ascending) { if (f_temp <= stp) x1 = x_temp; @@ -221,30 +165,95 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, x2 = x_temp; } } - return (x1 + x2) / 2; + *solution = (x1 + x2) / 2; + return 0; } -double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err) { - if (*err != 0) +double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, int side, int *err) { + if (config->ion_a <= 0) { + *err = DEDX_ERR_ION_A_REQUIRED; return -1; - if (config->loaded == 0) - dedx_load_config(ws, config, err); - if (*err != 0) + } + dedx_internal_lookup_data *data = load_and_get_dataset(ws, config, err); + if (data == NULL) return -1; - double emin = dedx_get_min_energy(config->program, config->ion); - double emax = dedx_get_max_energy(config->program, config->ion); + /* Stopping power vs. energy is not simply unimodal: real tables can rise + * to the Bragg peak, fall through the minimum-ionizing point, and rise + * again at relativistic energies (e.g. PSTAR/ICRU protons up to 10 GeV) + * -- an arbitrary requested STP can be reachable on more than one of + * these monotonic runs. Walk the exact tabulated knots once, bisecting + * every run that brackets stp, and keep the lowest- and highest-energy + * solutions found; side == 0 returns the low-energy one, side == 1 the + * high-energy one. For the common single-peak case this is exactly the + * old ascending/descending branch choice. See #121. */ + int found = 0; + double x_min_found = 0; + double x_max_found = 0; + int seg_start = 0; + int prev_dir = 0; + + for (int i = 1; i <= data->n; i++) { + int is_turning = 0; + int dir = 0; + if (i < data->n) { + double delta = (double) data->base[i].a - (double) data->base[i - 1].a; + dir = (delta > 0) ? 1 : (delta < 0 ? -1 : 0); + if (dir != 0) { + if (prev_dir == 0) + prev_dir = dir; + else if (dir != prev_dir) + is_turning = 1; + } + } + if (i == data->n || is_turning) { + int seg_end = (i == data->n) ? (data->n - 1) : (i - 1); + if (seg_end > seg_start) { + double solution; + if (bisect_monotonic_run(ws, config, data, seg_start, seg_end, stp, err, &solution) == 0) { + if (!found || solution < x_min_found) + x_min_found = solution; + if (!found || solution > x_max_found) + x_max_found = solution; + found = 1; + } else if (*err != 0) { + return -1; + } + } + if (is_turning) { + seg_start = i - 1; + prev_dir = dir; + } + } + } - double e_peak; - double max_stp; - double stp_at_emin; - if (find_stp_peak(ws, config, emin, emax, DEDX_MAX_STP_SAMPLES, &e_peak, &max_stp, &stp_at_emin) != 0) { + if (!found) { *err = DEDX_ERR_ENERGY_OUT_OF_RANGE; return -1; } + return (side == 0) ? x_min_found : x_max_found; +} + +double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err) { + dedx_internal_lookup_data *data = load_and_get_dataset(ws, config, err); + if (data == NULL) + return -1; + double min_stp; + double max_stp; + min_max_stp_over_table(data, &min_stp, &max_stp); return max_stp; } +double dedx_get_min_stp(dedx_workspace *ws, dedx_config *config, int *err) { + dedx_internal_lookup_data *data = load_and_get_dataset(ws, config, err); + if (data == NULL) + return -1; + double min_stp; + double max_stp; + min_max_stp_over_table(data, &min_stp, &max_stp); + return min_stp; +} + double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int *err) { if (config->ion_a <= 0) { *err = DEDX_ERR_ION_A_REQUIRED; diff --git a/tests/test_inverse_stp.c b/tests/test_inverse_stp.c index f799d5e..11e2f0e 100644 --- a/tests/test_inverse_stp.c +++ b/tests/test_inverse_stp.c @@ -5,13 +5,18 @@ #include #include -/* Regression tests for dedx_get_inverse_stp() and dedx_get_max_stp() - * (issue #121). The previous implementation used find_min() over a hardcoded - * x in [0.01, 10] to locate the energy of maximum stopping power, which does - * not match the real tabulated energy range; it also mapped side < 0 (never - * 0 or 1) to the low-energy branch, so side == 0 and side == 1 were - * indistinguishable. Both bugs made dedx_get_inverse_stp() fail with - * DEDX_ERR_ENERGY_OUT_OF_RANGE for ordinary proton/water queries. +/* Regression tests for dedx_get_inverse_stp(), dedx_get_max_stp(), and + * dedx_get_min_stp() (issue #121). + * + * The stopping-power-vs-energy curve is not simply unimodal: 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) up to + * the table's max energy. A requested STP in that high-energy dip can be + * reachable at two distinct energies, on top of the usual two solutions + * around the Bragg peak. dedx_get_inverse_stp() finds every monotonic run in + * the exact tabulated knots (no arbitrary sampling density) and returns the + * lowest-energy (side=0) or highest-energy (side=1) reachable solution. */ static int failures = 0; @@ -50,11 +55,9 @@ static dedx_config *make_config(int program, int ion, int ion_a, int target) { int main(void) { int err = 0; - /* --- Below the ascending branch's floor: only the descending branch can - * reach this STP, exactly the scenario that made the old find_min()-based - * code return DEDX_ERR_ENERGY_OUT_OF_RANGE for proton/water. Both sides - * must now fall back to the same (only reachable) descending-branch - * energy, with no error and a positive result. */ + /* --- STP below the Bragg peak's ascending floor: only the post-peak + * descending run reaches it, so side=0 and side=1 must agree on the + * same (only reachable) energy, with no error. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err); expect_int("workspace alloc err", err, DEDX_OK); @@ -78,10 +81,10 @@ int main(void) { dedx_free_workspace(ws, &err); } - /* --- Interior peak in the STP curve, STP reachable on both branches: - * side=0 (low energy / ascending) and side=1 (high energy / descending) - * must select distinct branches and each round-trip back to the - * requested STP. */ + /* --- Around the Bragg peak, STP reachable on both the ascending + * (pre-peak) and descending (post-peak) runs: side=0 and side=1 must + * select distinct energies and each round-trip back to the requested + * STP. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err); dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); @@ -111,22 +114,59 @@ int main(void) { dedx_free_workspace(ws, &err); } - /* --- Out-of-range STP values (above the maximum stopping power, or - * below the STP at max energy) must return an error rather than looping - * to a bogus or negative energy. */ + /* --- Relativistic-rise dip around the minimum-ionizing point (~3000-4000 + * MeV for proton/water PSTAR): a target STP in this dip is reachable at + * two widely separated energies, one on each side of the minimum. This + * is exactly the region the old find_min()-based implementation got + * wrong -- it assumed the whole post-peak range was a single descending + * branch and rejected every STP here as out of range. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + const float target_stp = 2.05f; + + err = 0; + double e_low = dedx_get_inverse_stp(ws, cfg, target_stp, 0, &err); + expect_int("dip stp side=0 err", err, DEDX_OK); + + err = 0; + double e_high = dedx_get_inverse_stp(ws, cfg, target_stp, 1, &err); + expect_int("dip stp side=1 err", err, DEDX_OK); + + expect_true("dip side=0 selects the lower energy", e_low < e_high); + expect_true("dip side=0 energy is in the pre-minimum range", e_low > 100.0 && e_low < 3000.0); + expect_true("dip side=1 energy is in the post-minimum range", e_high > 3000.0 && e_high < 10000.0); + + int verify_err = 0; + double stp_low = dedx_get_stp(ws, cfg, (float) e_low, &verify_err); + expect_int("dip side=0 roundtrip err", verify_err, DEDX_OK); + expect_near("dip side=0 roundtrip stp", stp_low, target_stp, 1e-2); + + verify_err = 0; + double stp_high = dedx_get_stp(ws, cfg, (float) e_high, &verify_err); + expect_int("dip side=1 roundtrip err", verify_err, DEDX_OK); + expect_near("dip side=1 roundtrip stp", stp_high, target_stp, 1e-2); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- Out-of-range STP values (above the global maximum, or below the + * global minimum) must return an error rather than looping to a bogus or + * negative energy. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err); dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); err = 0; double e_too_high = dedx_get_inverse_stp(ws, cfg, 2000.0f, 0, &err); - expect_int("stp above peak err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); - expect_near("stp above peak sentinel", e_too_high, -1.0, 1e-9); + expect_int("stp above global max err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); + expect_near("stp above global max sentinel", e_too_high, -1.0, 1e-9); err = 0; - double e_too_low = dedx_get_inverse_stp(ws, cfg, 0.5f, 0, &err); - expect_int("stp below max-energy floor err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); - expect_near("stp below floor sentinel", e_too_low, -1.0, 1e-9); + double e_too_low = dedx_get_inverse_stp(ws, cfg, 1.0f, 0, &err); + expect_int("stp below global min err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); + expect_near("stp below global min sentinel", e_too_low, -1.0, 1e-9); dedx_free_config(cfg, &err); dedx_free_workspace(ws, &err); @@ -144,30 +184,45 @@ int main(void) { dedx_free_workspace(ws, &err); } - /* --- dedx_get_max_stp() reports the true maximum of the curve. */ + /* --- dedx_get_max_stp() / dedx_get_min_stp() report the true extrema of + * the curve: the Bragg peak and the minimum-ionizing point. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err); dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + err = 0; - double peak = dedx_get_max_stp(ws, cfg, &err); + double max_stp = dedx_get_max_stp(ws, cfg, &err); expect_int("max stp err", err, DEDX_OK); - /* Brute-force reference: sample far more densely and take the max. */ + err = 0; + double min_stp = dedx_get_min_stp(ws, cfg, &err); + expect_int("min stp err", err, DEDX_OK); + + expect_true("min stp is less than max stp", min_stp < max_stp); + expect_true("min stp is positive", min_stp > 0.0); + + /* Brute-force reference: sample far more densely than the table + * itself and take the min/max, as an independent cross-check. */ float emin = dedx_get_min_energy(cfg->program, cfg->ion); float emax = dedx_get_max_energy(cfg->program, cfg->ion); double log_emin = log((double) emin); double log_emax = log((double) emax); - double reference_peak = 0.0; - const int n = 2000; + double reference_max = 0.0; + double reference_min = HUGE_VAL; + const int n = 5000; for (int i = 0; i < n; i++) { float e = (float) exp(log_emin + (log_emax - log_emin) * i / (n - 1)); int sample_err = 0; double s = dedx_get_stp(ws, cfg, e, &sample_err); - if (sample_err == 0 && s > reference_peak) - reference_peak = s; + if (sample_err != 0) + continue; + if (s > reference_max) + reference_max = s; + if (s < reference_min) + reference_min = s; } - expect_near("max stp matches brute-force reference", peak, reference_peak, 5e-2); - expect_true("max stp is a real maximum", peak > 0.0); + expect_near("max stp matches brute-force reference", max_stp, reference_max, 5e-2); + expect_near("min stp matches brute-force reference", min_stp, reference_min, 5e-2); dedx_free_config(cfg, &err); dedx_free_workspace(ws, &err); From d2516611efdb5697a9ed7ef3556a89a48053383d Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Fri, 24 Jul 2026 19:11:08 +0200 Subject: [PATCH 5/6] polish: relative bisection tolerance, no ternaries, more comments - load_and_get_dataset() -> get_loaded_dataset(): still static/file-private (only this translation unit needs it, unlike the dedx_internal_* functions declared in headers for cross-file use); comment now says so explicitly. - Replace the fixed acc = 1e-5 (absolute, in the same MeV units as the bisected energy) with a tolerance relative to the current bracket (DEDX_INVERSE_STP_REL_ACC = 1e-6). 1e-5 MeV is 10 eV, not 10 keV -- and as an absolute value it's a poor fit for a domain spanning many decades (proton tables run ~0.001-10000 MeV/nucleon): far too loose relative to the low end, far tighter than dedx_get_stp()'s float precision can even represent at the high end (it truncates the search value to float before evaluating the spline). A relative tolerance keeps achieved precision consistent across the whole range and matches float's own ~1.19e-7 relative precision instead of overshooting it. Documented inline with the reasoning and the eV-scale sanity check. - Remove every ternary in bisect_monotonic_run() and dedx_get_inverse_stp(), replaced with if/else and comments explaining what each branch means physically (bracket direction, run-boundary detection, ascending vs descending step selection). No behavior change beyond the tolerance switch, which tightens precision at low energies and avoids wasted iterations at high energies; verified against the same relativistic-rise and cross-program checks as before. Co-Authored-By: Claude Sonnet 5 --- src/dedx_tools.c | 96 +++++++++++++++++++++++++++++++++++++++++------- 1 file changed, 83 insertions(+), 13 deletions(-) diff --git a/src/dedx_tools.c b/src/dedx_tools.c index 4213aed..c265d0e 100644 --- a/src/dedx_tools.c +++ b/src/dedx_tools.c @@ -65,8 +65,12 @@ static double adapt_stp(double energy, dedx_tools_settings *set) { * of re-sampling the curve at an arbitrary density. dedx_get_stp() evaluates * a spline built from precisely these knots (dedx_spline.c: coef[i].a = * stopping[i]), so scanning them finds every feature the public curve can - * actually exhibit -- no synthetic sampling can do better or worse. */ -static dedx_internal_lookup_data *load_and_get_dataset(dedx_workspace *ws, dedx_config *config, int *err) { + * actually exhibit -- no synthetic sampling can do better or worse. + * + * `static` (file-private, not declared in dedx_tools.h): only this + * translation unit needs it, unlike the dedx_internal_* functions declared + * in headers such as dedx_data_access.h for use across multiple .c files. */ +static dedx_internal_lookup_data *get_loaded_dataset(dedx_workspace *ws, dedx_config *config, int *err) { if (*err != 0) return NULL; if (config->loaded == 0) @@ -122,6 +126,19 @@ double dedx_get_inverse_csda(dedx_workspace *ws, dedx_config *config, float rang return (min + max) / 2; } +/* Bisection convergence, as a fraction of the current bracket's energy scale + * (see the "relative, not absolute" note in bisect_monotonic_run() below for + * why this is relative). 1e-6 is about ten times coarser than a `float`'s + * own precision (~1.19e-7 relative, i.e. 24 bits of mantissa): dedx_get_stp() + * truncates the search value to a float before evaluating the spline, so + * refining much further than the float itself can represent would just + * spend extra iterations re-measuring rounding noise. 1e-6 still comfortably + * resolves the whole tabulated range: even the widest single monotonic run + * seen in practice (proton/water PSTAR's post-peak run, roughly 0.08 to + * 10000 MeV/nucleon) converges in under 30 bisection steps, since halving + * the bracket is a log2(width/tolerance) process regardless of scale. */ +#define DEDX_INVERSE_STP_REL_ACC 1e-6 + /* Bisect the monotonic run of knots [lo_idx, hi_idx] for the energy where * the spline equals stp, refining via dedx_get_stp() (not linear knot * interpolation) so the result matches what callers would measure with the @@ -139,26 +156,60 @@ static int bisect_monotonic_run(dedx_workspace *ws, double hi_x = data->base[hi_idx].x; double lo_v = data->base[lo_idx].a; double hi_v = data->base[hi_idx].a; - double range_min = lo_v < hi_v ? lo_v : hi_v; - double range_max = lo_v > hi_v ? lo_v : hi_v; + + /* This run's STP values span [range_min, range_max] regardless of + * whether the run is rising or falling; stp must land in that span for + * a solution to exist anywhere on this run. */ + double range_min; + double range_max; + if (lo_v < hi_v) { + range_min = lo_v; + range_max = hi_v; + } else { + range_min = hi_v; + range_max = lo_v; + } if (stp < range_min || stp > range_max) return -1; - int ascending = hi_v >= lo_v; + /* Whether STP rises (ascending) or falls (descending) from lo_x to + * hi_x tells us which half to keep at each bisection step below. */ + int ascending; + if (hi_v >= lo_v) + ascending = 1; + else + ascending = 0; + + /* Standard bisection: x1/x2 bracket the root, i.e. dedx_get_stp(x1) and + * dedx_get_stp(x2) sit on opposite sides of the target stp. Each step + * evaluates the midpoint and keeps whichever half still brackets stp, + * halving the bracket width every iteration. + * + * The stopping test compares the bracket width to x2 itself rather than + * to a fixed value: the tabulated energy range spans many decades (e.g. + * proton tables run from ~0.001 to 10000 MeV/nucleon), so a single fixed + * absolute tolerance would be far too loose relative to energies at the + * low end of that range and far tighter than meaningful at the high end + * (see DEDX_INVERSE_STP_REL_ACC above). A tolerance relative to the + * current bracket keeps the achieved *relative* precision the same + * regardless of where in that range the root happens to fall. */ double x1 = lo_x; double x2 = hi_x; - double acc = 1e-5; - while (fabs(x1 - x2) > acc) { + while (fabs(x1 - x2) > DEDX_INVERSE_STP_REL_ACC * fabs(x2)) { double x_temp = (x1 + x2) / 2; double f_temp = dedx_get_stp(ws, config, (float) x_temp, err); if (*err != 0) return -1; if (ascending) { + /* STP too low at the midpoint -> the root is further up. */ if (f_temp <= stp) x1 = x_temp; else x2 = x_temp; } else { + /* STP still at/above target at the midpoint -> the root is + * further up (we're descending, so STP keeps falling as energy + * rises). */ if (f_temp >= stp) x1 = x_temp; else @@ -174,7 +225,7 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, *err = DEDX_ERR_ION_A_REQUIRED; return -1; } - dedx_internal_lookup_data *data = load_and_get_dataset(ws, config, err); + dedx_internal_lookup_data *data = get_loaded_dataset(ws, config, err); if (data == NULL) return -1; @@ -194,11 +245,21 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, int prev_dir = 0; for (int i = 1; i <= data->n; i++) { + /* dir is the sign of the step from knot i-1 to knot i: +1 rising, + * -1 falling, 0 flat (a flat step doesn't end the current run -- + * prev_dir is left untouched so a flat plateau stays part of + * whichever run it interrupts). Only computed for i < data->n, + * since there is no knot i to compare against once i == data->n. */ int is_turning = 0; int dir = 0; if (i < data->n) { double delta = (double) data->base[i].a - (double) data->base[i - 1].a; - dir = (delta > 0) ? 1 : (delta < 0 ? -1 : 0); + if (delta > 0) + dir = 1; + else if (delta < 0) + dir = -1; + else + dir = 0; if (dir != 0) { if (prev_dir == 0) prev_dir = dir; @@ -206,8 +267,15 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, is_turning = 1; } } + /* A run ends either at a direction reversal (is_turning, and the + * turning knot i-1 is shared with the next run) or at the last + * knot (i == data->n, since there's nowhere further to walk). */ if (i == data->n || is_turning) { - int seg_end = (i == data->n) ? (data->n - 1) : (i - 1); + int seg_end; + if (i == data->n) + seg_end = data->n - 1; + else + seg_end = i - 1; if (seg_end > seg_start) { double solution; if (bisect_monotonic_run(ws, config, data, seg_start, seg_end, stp, err, &solution) == 0) { @@ -231,11 +299,13 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, *err = DEDX_ERR_ENERGY_OUT_OF_RANGE; return -1; } - return (side == 0) ? x_min_found : x_max_found; + if (side == 0) + return x_min_found; + return x_max_found; } double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err) { - dedx_internal_lookup_data *data = load_and_get_dataset(ws, config, err); + dedx_internal_lookup_data *data = get_loaded_dataset(ws, config, err); if (data == NULL) return -1; double min_stp; @@ -245,7 +315,7 @@ double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err) { } double dedx_get_min_stp(dedx_workspace *ws, dedx_config *config, int *err) { - dedx_internal_lookup_data *data = load_and_get_dataset(ws, config, err); + dedx_internal_lookup_data *data = get_loaded_dataset(ws, config, err); if (data == NULL) return -1; double min_stp; From 0202a1b887d9bcb9718204f941926579dab31f63 Mon Sep 17 00:00:00 2001 From: Leszek Grzanka Date: Fri, 24 Jul 2026 20:09:49 +0200 Subject: [PATCH 6/6] address review: log-space bisection, get_loaded_dataset robustness, docs - bisect_monotonic_run(): bisect in log-energy space instead of raw energy. The tabulated grid is predominantly log-spaced (constant energy *ratio* per NIST-derived tables), and an absolute tolerance in log-space is exactly a relative tolerance on the energy result (d(ln x) = dx/x) -- more principled than the previous fabs(x1-x2) > REL_ACC*fabs(x2) approximation. Renamed DEDX_INVERSE_STP_REL_ACC -> ...LOG_ACC to match. - get_loaded_dataset(): reset *err = DEDX_OK on entry instead of treating an incoming nonzero value as a precondition (matches dedx_get_stp()); also reload when config->loaded == 1 but config->cfg_id isn't a valid slot in *this* workspace (e.g. the config was loaded into a workspace that was since freed and replaced), instead of just erroring out. Both from Copilot review on #148. - Doc/comment fixes from the same review: dedx_get_inverse_stp()/ dedx_get_max_stp()/dedx_get_min_stp() docstrings now say the config loads automatically rather than implying a precondition; side is documented as exactly what the code checks (side == 0 vs any other value); the "finds every feature" comment no longer overclaims exhaustiveness over the interpolated spline (only over the tabulated points); min/max-stp docstrings say "computed exactly from the tabulated data points" instead of implying a guarantee about the continuous interpolated curve. - tests/test_inverse_stp.c: derive the out-of-range test's STP bounds from dedx_get_max_stp()/dedx_get_min_stp() instead of hardcoded constants; add regression tests for get_loaded_dataset()'s load-failure propagation (through all three affected functions) and its reload-on-stale-workspace fix. - Remove a dead `else dir = 0;` branch (dir is already 0-initialized). Patch coverage: 92.2% (was 91.8%); the remaining uncovered lines in dedx_tools.c are defensive branches unreachable via the public API given dedx_load_config()'s and the bisection's own invariants (verified, not just assumed) -- consistent with similar defensive branches elsewhere in this codebase. Co-Authored-By: Claude Sonnet 5 --- include/dedx_tools.h | 42 ++++++++----- src/dedx_tools.c | 129 ++++++++++++++++++++++++++------------- tests/test_inverse_stp.c | 72 +++++++++++++++++++++- 3 files changed, 181 insertions(+), 62 deletions(-) diff --git a/include/dedx_tools.h b/include/dedx_tools.h index 3d8d822..5706235 100644 --- a/include/dedx_tools.h +++ b/include/dedx_tools.h @@ -41,17 +41,23 @@ double dedx_get_csda(dedx_workspace *ws, dedx_config *config, float energy, int * through the minimum-ionizing point, and rise again at relativistic * energies (e.g. proton tables extending to several GeV), so a requested * @p stp can be reachable at more than one energy. Among all reachable - * energies, @p side selects which one to return: 0 = the lowest, 1 = the - * highest. When only one energy reaches @p stp, @p side has no effect. + * energies, @p side selects which one to return: 0 selects the lowest; + * any other value (the documented convention is 1) selects the highest. + * When only one energy reaches @p stp, @p side has no effect. * - * @param[in] ws Workspace with a loaded configuration. - * @param[in] config Loaded configuration. - * @param[in] stp Target stopping power in MeV cm²/g. Must lie between - * dedx_get_min_stp() and dedx_get_max_stp() (inclusive) - * for this configuration; otherwise the value is - * unreachable and an error is returned. - * @param[in] side 0 = lowest-energy solution, 1 = highest-energy - * solution. + * @param[in] ws Workspace to hold the loaded configuration. + * @param[in] config Configuration to invert against. Loaded automatically + * if not already (config->loaded == 0): on return, + * config->loaded and config->cfg_id reflect that load, + * same as after an explicit dedx_load_config() call. + * @param[in] stp Target stopping power in MeV cm²/g. In practice must + * lie between dedx_get_min_stp() and dedx_get_max_stp() + * (inclusive) for this configuration -- both computed + * from the same tabulated data points this function + * searches -- or the value is unreachable and an error + * is returned. + * @param[in] side 0 = lowest-energy solution, any other value = + * highest-energy solution. * @param[out] err Error code; 0 on success. * @return Energy in MeV/nucl (MeV per nucleon), or -1 on error. */ @@ -62,10 +68,12 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, * * 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 - * not computed by this library). + * not computed by this library). It is computed exactly from the tabulated + * data points backing the configuration, not sampled or estimated. * - * @param[in] ws Workspace with a loaded configuration. - * @param[in] config Loaded configuration. + * @param[in] ws Workspace to hold the loaded configuration. + * @param[in] config Configuration to inspect. Loaded automatically if not + * already (config->loaded == 0); see dedx_get_inverse_stp(). * @param[out] err Error code; 0 on success. * @return Maximum stopping power in MeV cm²/g, or -1 on error. */ @@ -76,10 +84,12 @@ double dedx_get_max_stp(dedx_workspace *ws, dedx_config *config, int *err); * * For tables that reach relativistic energies this is typically the * minimum-ionizing point, not the value at either tabulated endpoint (the - * curve can rise again past it — see dedx_get_inverse_stp()). + * curve can rise again past it — see dedx_get_inverse_stp()). Computed + * exactly from the tabulated data points, like dedx_get_max_stp(). * - * @param[in] ws Workspace with a loaded configuration. - * @param[in] config Loaded configuration. + * @param[in] ws Workspace to hold the loaded configuration. + * @param[in] config Configuration to inspect. Loaded automatically if not + * already (config->loaded == 0); see dedx_get_inverse_stp(). * @param[out] err Error code; 0 on success. * @return Minimum stopping power in MeV cm²/g, or -1 on error. */ diff --git a/src/dedx_tools.c b/src/dedx_tools.c index c265d0e..adb4d89 100644 --- a/src/dedx_tools.c +++ b/src/dedx_tools.c @@ -64,21 +64,44 @@ static double adapt_stp(double energy, dedx_tools_settings *set) { * callers can walk its exact tabulated (energy, STP) knots directly instead * of re-sampling the curve at an arbitrary density. dedx_get_stp() evaluates * a spline built from precisely these knots (dedx_spline.c: coef[i].a = - * stopping[i]), so scanning them finds every feature the public curve can - * actually exhibit -- no synthetic sampling can do better or worse. + * stopping[i]), so scanning them finds every feature present in the + * tabulated values -- though not necessarily every feature of the + * interpolated curve dedx_get_stp() actually evaluates, since a spline can + * in principle overshoot slightly between two knots. That's a property of + * the interpolation, not something this scan (or any sampling density) can + * see; it hasn't been observed in the tables this library ships. * * `static` (file-private, not declared in dedx_tools.h): only this * translation unit needs it, unlike the dedx_internal_* functions declared * in headers such as dedx_data_access.h for use across multiple .c files. */ static dedx_internal_lookup_data *get_loaded_dataset(dedx_workspace *ws, dedx_config *config, int *err) { - if (*err != 0) - return NULL; + /* Public entrypoints reset *err themselves (e.g. dedx_get_stp()) rather + * than treating an incoming nonzero value as a precondition -- do the + * same here so a stale error code left over in the caller's variable + * from an unrelated earlier call can't make a valid request fail. */ + *err = DEDX_OK; + + /* config->loaded == 1 only means *some* workspace holds this config's + * data at config->cfg_id -- not necessarily this ws (the config could + * have been loaded into a workspace that was since freed and replaced, + * or the caller could be reusing one dedx_config across workspaces). + * Reload whenever the config hasn't been loaded at all, or its cfg_id + * isn't a valid slot in *this* workspace; otherwise skip the reload on + * the common already-loaded-here path. */ + int needs_load = 0; if (config->loaded == 0) + needs_load = 1; + else if (config->cfg_id < 0 || config->cfg_id >= ws->active_datasets) + needs_load = 1; + if (needs_load) dedx_load_config(ws, config, err); if (*err != 0) return NULL; + int id = config->cfg_id; if (id < 0 || id >= ws->active_datasets) { + /* Defensive: dedx_load_config() succeeding is expected to always + * leave a valid cfg_id, so this should be unreachable in practice. */ *err = DEDX_ERR_INVALID_DATASET_ID; return NULL; } @@ -126,18 +149,17 @@ double dedx_get_inverse_csda(dedx_workspace *ws, dedx_config *config, float rang return (min + max) / 2; } -/* Bisection convergence, as a fraction of the current bracket's energy scale - * (see the "relative, not absolute" note in bisect_monotonic_run() below for - * why this is relative). 1e-6 is about ten times coarser than a `float`'s - * own precision (~1.19e-7 relative, i.e. 24 bits of mantissa): dedx_get_stp() - * truncates the search value to a float before evaluating the spline, so - * refining much further than the float itself can represent would just - * spend extra iterations re-measuring rounding noise. 1e-6 still comfortably - * resolves the whole tabulated range: even the widest single monotonic run - * seen in practice (proton/water PSTAR's post-peak run, roughly 0.08 to - * 10000 MeV/nucleon) converges in under 30 bisection steps, since halving - * the bracket is a log2(width/tolerance) process regardless of scale. */ -#define DEDX_INVERSE_STP_REL_ACC 1e-6 +/* Bisection convergence, as an absolute tolerance in *log-energy* space (see + * bisect_monotonic_run() below for why bisecting there instead of in the raw + * energy value). For small differences, d(ln x) = dx/x, i.e. a difference in + * log-space is directly a *relative* difference in the original energy -- + * so this one constant fixes the achieved relative precision on the energy + * result regardless of scale. 1e-6 is about ten times coarser than a + * `float`'s own precision (~1.19e-7 relative, i.e. 24 bits of mantissa): + * dedx_get_stp() truncates the search value to a float before evaluating + * the spline, so refining much further than the float itself can represent + * would just spend extra iterations re-measuring rounding noise. */ +#define DEDX_INVERSE_STP_LOG_ACC 1e-6 /* Bisect the monotonic run of knots [lo_idx, hi_idx] for the energy where * the spline equals stp, refining via dedx_get_stp() (not linear knot @@ -180,43 +202,56 @@ static int bisect_monotonic_run(dedx_workspace *ws, else ascending = 0; - /* Standard bisection: x1/x2 bracket the root, i.e. dedx_get_stp(x1) and - * dedx_get_stp(x2) sit on opposite sides of the target stp. Each step - * evaluates the midpoint and keeps whichever half still brackets stp, - * halving the bracket width every iteration. + /* Bisect in log-energy space rather than raw energy. Two reasons: * - * The stopping test compares the bracket width to x2 itself rather than - * to a fixed value: the tabulated energy range spans many decades (e.g. - * proton tables run from ~0.001 to 10000 MeV/nucleon), so a single fixed - * absolute tolerance would be far too loose relative to energies at the - * low end of that range and far tighter than meaningful at the high end - * (see DEDX_INVERSE_STP_REL_ACC above). A tolerance relative to the - * current bracket keeps the achieved *relative* precision the same - * regardless of where in that range the root happens to fall. */ - double x1 = lo_x; - double x2 = hi_x; - while (fabs(x1 - x2) > DEDX_INVERSE_STP_REL_ACC * fabs(x2)) { - double x_temp = (x1 + x2) / 2; + * 1. The tabulated energy grid is itself predominantly log-spaced (e.g. + * NIST PSTAR/ASTAR-derived tables step by a roughly constant energy + * *ratio*, not a constant energy difference, so that a fixed number + * of points can cover ~0.001 to 10000 MeV/nucleon with even relative + * resolution). Bisecting in the same log-scale as the data matches + * the coordinate the curve actually varies smoothly in. + * + * 2. It makes the convergence tolerance exactly a relative tolerance on + * the energy result (see DEDX_INVERSE_STP_LOG_ACC above), rather than + * an approximation of one. A run spanning many decades (e.g. the + * post-peak run for proton/water PSTAR, ~0.08 to 10000 MeV/nucleon) + * then converges in a number of steps set by log2(ln(hi_x/lo_x) / + * DEDX_INVERSE_STP_LOG_ACC) -- a few dozen steps regardless of scale + * or of where in the run the root happens to fall, since bisection + * always halves the bracket every step no matter which half contains + * the root. + * + * dedx_get_stp() evaluates the actual interpolating spline at any energy + * (not just at knots), so this is correct even on the tail of a table + * where the raw knot spacing itself happens to be linear rather than + * log-spaced (observed for a few tables at their highest energies) -- + * the spline doesn't care how its own knots were spaced, and neither + * does bisecting the search variable in log space. */ + double log_lo = log(lo_x); + double log_hi = log(hi_x); + while (fabs(log_lo - log_hi) > DEDX_INVERSE_STP_LOG_ACC) { + double log_mid = (log_lo + log_hi) / 2; + double x_temp = exp(log_mid); double f_temp = dedx_get_stp(ws, config, (float) x_temp, err); if (*err != 0) return -1; if (ascending) { /* STP too low at the midpoint -> the root is further up. */ if (f_temp <= stp) - x1 = x_temp; + log_lo = log_mid; else - x2 = x_temp; + log_hi = log_mid; } else { /* STP still at/above target at the midpoint -> the root is * further up (we're descending, so STP keeps falling as energy * rises). */ if (f_temp >= stp) - x1 = x_temp; + log_lo = log_mid; else - x2 = x_temp; + log_hi = log_mid; } } - *solution = (x1 + x2) / 2; + *solution = exp((log_lo + log_hi) / 2); return 0; } @@ -235,9 +270,9 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, * -- an arbitrary requested STP can be reachable on more than one of * these monotonic runs. Walk the exact tabulated knots once, bisecting * every run that brackets stp, and keep the lowest- and highest-energy - * solutions found; side == 0 returns the low-energy one, side == 1 the - * high-energy one. For the common single-peak case this is exactly the - * old ascending/descending branch choice. See #121. */ + * solutions found across ALL of them; side == 0 returns the lowest, any + * other side value returns the highest. For the common single-peak case + * this is exactly the old ascending/descending branch choice. See #121. */ int found = 0; double x_min_found = 0; double x_max_found = 0; @@ -251,15 +286,13 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, * whichever run it interrupts). Only computed for i < data->n, * since there is no knot i to compare against once i == data->n. */ int is_turning = 0; - int dir = 0; + int dir = 0; /* stays 0 (flat) unless one of the branches below fires */ if (i < data->n) { double delta = (double) data->base[i].a - (double) data->base[i - 1].a; if (delta > 0) dir = 1; else if (delta < 0) dir = -1; - else - dir = 0; if (dir != 0) { if (prev_dir == 0) prev_dir = dir; @@ -279,6 +312,16 @@ double dedx_get_inverse_stp(dedx_workspace *ws, dedx_config *config, float stp, if (seg_end > seg_start) { double solution; if (bisect_monotonic_run(ws, config, data, seg_start, seg_end, stp, err, &solution) == 0) { + /* This run brackets stp, so it contributes one solution. + * Track the lowest- and highest-energy solutions seen so + * far across every run processed so far in this loop -- + * not just this one run -- so that once the whole table + * has been walked, x_min_found/x_max_found hold the + * overall lowest/highest reachable energy regardless of + * how many separate runs bracketed stp along the way. + * `!found` covers the very first solution: there is + * nothing yet to compare it against, so it always + * becomes both the running min and the running max. */ if (!found || solution < x_min_found) x_min_found = solution; if (!found || solution > x_max_found) diff --git a/tests/test_inverse_stp.c b/tests/test_inverse_stp.c index 11e2f0e..013f033 100644 --- a/tests/test_inverse_stp.c +++ b/tests/test_inverse_stp.c @@ -153,18 +153,30 @@ int main(void) { /* --- Out-of-range STP values (above the global maximum, or below the * global minimum) must return an error rather than looping to a bogus or - * negative energy. */ + * negative energy. The bounds are derived from dedx_get_max_stp()/ + * dedx_get_min_stp() themselves rather than hardcoded constants, so this + * test stays valid if the underlying tabulated data is ever updated. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err); dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); err = 0; - double e_too_high = dedx_get_inverse_stp(ws, cfg, 2000.0f, 0, &err); + double max_stp = dedx_get_max_stp(ws, cfg, &err); + expect_int("out-of-range setup: max stp err", err, DEDX_OK); + err = 0; + double min_stp = dedx_get_min_stp(ws, cfg, &err); + expect_int("out-of-range setup: min stp err", err, DEDX_OK); + + float stp_above_max = (float) (max_stp * 2.0); + float stp_below_min = (float) (min_stp * 0.5); + + err = 0; + double e_too_high = dedx_get_inverse_stp(ws, cfg, stp_above_max, 0, &err); expect_int("stp above global max err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); expect_near("stp above global max sentinel", e_too_high, -1.0, 1e-9); err = 0; - double e_too_low = dedx_get_inverse_stp(ws, cfg, 1.0f, 0, &err); + double e_too_low = dedx_get_inverse_stp(ws, cfg, stp_below_min, 0, &err); expect_int("stp below global min err", err, DEDX_ERR_ENERGY_OUT_OF_RANGE); expect_near("stp below global min sentinel", e_too_low, -1.0, 1e-9); @@ -172,6 +184,60 @@ int main(void) { dedx_free_workspace(ws, &err); } + /* --- A config that fails to load (here: a Bethe-type program with a + * custom target and no density set, DEDX_ERR_RHO_REQUIRED -- the same + * failure test_error_codes.c uses) must propagate that error cleanly out + * of all three functions that route through get_loaded_dataset(), not + * crash or silently return a bogus value. */ + { + dedx_workspace *ws = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_BETHE_EXT00, DEDX_PROTON, 1, 0); + + err = 0; + double e = dedx_get_inverse_stp(ws, cfg, 100.0f, 0, &err); + expect_true("unloadable config: inverse_stp err set", err != DEDX_OK); + expect_near("unloadable config: inverse_stp sentinel", e, -1.0, 1e-9); + + err = 0; + double max_stp = dedx_get_max_stp(ws, cfg, &err); + expect_true("unloadable config: max_stp err set", err != DEDX_OK); + expect_near("unloadable config: max_stp sentinel", max_stp, -1.0, 1e-9); + + err = 0; + double min_stp = dedx_get_min_stp(ws, cfg, &err); + expect_true("unloadable config: min_stp err set", err != DEDX_OK); + expect_near("unloadable config: min_stp sentinel", min_stp, -1.0, 1e-9); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws, &err); + } + + /* --- A config already loaded into a workspace that was since freed and + * replaced (config->loaded == 1, but config->cfg_id is not a valid slot + * in the new workspace) must be transparently reloaded into the new + * workspace rather than erroring out -- get_loaded_dataset() detects the + * stale cfg_id and reloads instead of trusting the loaded flag alone. */ + { + dedx_workspace *ws1 = dedx_allocate_workspace(1, &err); + dedx_config *cfg = make_config(DEDX_PSTAR, DEDX_PROTON, 1, DEDX_WATER); + + err = 0; + double max_stp_ws1 = dedx_get_max_stp(ws1, cfg, &err); + expect_int("stale workspace: first load err", err, DEDX_OK); + expect_true("stale workspace: cfg marked loaded", cfg->loaded != 0); + + dedx_free_workspace(ws1, &err); + + dedx_workspace *ws2 = dedx_allocate_workspace(1, &err); + err = 0; + double max_stp_ws2 = dedx_get_max_stp(ws2, cfg, &err); + expect_int("stale workspace: reload into ws2 err", err, DEDX_OK); + expect_near("stale workspace: same result after reload", max_stp_ws2, max_stp_ws1, 1e-9); + + dedx_free_config(cfg, &err); + dedx_free_workspace(ws2, &err); + } + /* --- ion_a <= 0 is rejected up front, same contract as before. */ { dedx_workspace *ws = dedx_allocate_workspace(1, &err);