[BUG] Bounds clamp in OGC_AmsterdamUMC_Bayesian_biexp.ivim_fit() writes to index 0 instead of i

[Devguru-codes]

Is there an existing issue for this?

• I have searched the existing issues

Describe the bug

In OGC_AmsterdamUMC_Bayesian_biexp.ivim_fit(), the loop that pulls the segmented pre-fit inside bounds before it's used as the Bayesian starting point writes to index 0 rather than i:

# src/standardized/OGC_AmsterdamUMC_Bayesian_biexp.py, lines 91-95
fit_results = np.array(fit_results + (1,))          # [D, f, Dp, S0]
for i in range(4):
    if fit_results[i] < bounds[0][i] : fit_results[0] = bounds[0][i]+epsilon
    if fit_results[i] > bounds[1][i] : fit_results[0] = bounds[1][i]-epsilon

So when any parameter is out of bounds, D gets overwritten instead of that parameter. The same clamp in ivim_fit_full_volume() further down (lines 137-145) uses fit_results[i, below], so the two paths in the file don't agree.

What makes this more than cosmetic is that flat_neg_log_prior returns 1e10 outside its ranges - so the point of this clamp, and of epsilon, is to nudge x0 strictly inside the prior's support. With the wrong index the parameter that's actually outside never moves, and the MAP minimisation starts on a flat 1e10 plateau. On the Liver signal with S0 bounds [1.1, 1.3]:

segmented pre-fit   x0 = [0.0015   0.1099 0.0993 1.0     ]   prior = 1e10
current clamp       x0 = [0.005    0.1099 0.0993 1.1     ]   prior = 1e10   <- after scipy clips x0
with [i] instead    x0 = [0.0015   0.1099 0.0993 1.100001]   prior = 0

I should be upfront that this is dormant under default settings. fit_segmented() returns (D, f, Dp) already constrained by curve_fit to these same bounds, so those three can never trip the check - only S0 can, and it isn't fitted, it's appended as the literal 1 on line 92. With the default S0 ∈ [0.7, 1.3] the clamp never fires: 0/14 tissues in generic.json, and 0/280 with added noise at SNR 100/30/15.

But with S0 ∈ [1.1, 1.3] - which is what test_bounds uses - it fires on all 14 tissues and every fitted value shifts. Muscle, for example, goes from f = 0.155, Dp = 0.20 to f = 0.100, Dp = 0.026, and the worst relative error in D across the set is about 64%.

Screenshots [optional]

No response

Steps To Reproduce

import json
import numpy as np
from src.original.fitting.OGC_AmsterdamUMC.LSQ_fitting import fit_segmented, flat_neg_log_prior

with open("tests/IVIMmodels/unit_tests/generic.json") as f:
    data = json.load(f)
bvals = np.array(data.pop("config")["bvalues"])

B = {"S0": [1.1, 1.3], "f": [0, 1.0], "Dp": [0.005, 0.2], "D": [0, 0.005]}   # S0 range excludes 1.0
bounds = ([B["D"][0], B["f"][0], B["Dp"][0], B["S0"][0]],
          [B["D"][1], B["f"][1], B["Dp"][1], B["S0"][1]])
prior = flat_neg_log_prior([B["D"][0], B["D"][1]], [B["f"][0], B["f"][1]],
                           [B["Dp"][0], B["Dp"][1]], [B["S0"][0], B["S0"][1]])

sig = np.asarray(data["Liver"]["data"], float)
sig = sig / sig[0]
pre = np.array(fit_segmented(bvals, sig, bounds=bounds, cutoff=np.array([200]),
                             p0=[0.001, 0.1, 0.01, 1.0]) + (1,))

eps = 1e-6
current, corrected = pre.copy(), pre.copy()
for i in range(4):
    if pre[i] < bounds[0][i]:
        current[0] = bounds[0][i] + eps      # what the code does today
        corrected[i] = bounds[0][i] + eps    # what [i] would do
    if pre[i] > bounds[1][i]:
        current[0] = bounds[1][i] - eps
        corrected[i] = bounds[1][i] - eps

clip = lambda x: np.clip(x, bounds[0], bounds[1])   # scipy.minimize does this internally
print("current  x0:", np.round(clip(current), 6),   "prior =", prior(clip(current)))
print("with [i] x0:", np.round(clip(corrected), 6), "prior =", prior(clip(corrected)))

current  x0: [0.005    0.109872 0.099329 1.1     ] prior = 10000000000.0
with [i] x0: [0.0015   0.109872 0.099329 1.100001] prior = 0

scipy.minimize silently clips an out-of-bounds x0 into bounds, and S0 = 1.1 is still rejected by the prior because the comparison is a strict < - which is what the epsilon is there for.

Expected behavior

The clamp should move the offending parameter inside its own bounds, so x0 ends up strictly within the prior's support:

for i in range(4):
    if fit_results[i] < bounds[0][i] : fit_results[i] = bounds[0][i]+epsilon
    if fit_results[i] > bounds[1][i] : fit_results[i] = bounds[1][i]-epsilon

which is what ivim_fit_full_volume() already does.

Additional context

I went looking at the original code before writing this up. doc/code_contributions_record.csv credits the algorithm to Oliver Gurney-Champion sir / Sebastiano Barbieri sir (DOI 10.1002/mrm.28852) and names fit_bayesian_array as the contributed function. In IVIMNET that function reads:

Dt0, Fp0, Dp0, S00 = paramslsq            # all four come from the LSQ pre-fit
x0 = [Dt0[i], Fp0[i], Dp0[i], S00[i]]     # S0 taken from the fit, not hard-coded

So there's no clamping loop in the original at all - it's specific to this wrapper, and none of this is a problem with the published implementation. IVIMNET also never hard-codes S0; its docstring says arg.fitS0 --> False fixes S0 to 1, True fits S0. This wrapper runs with fitS0=True but starts S0 from a literal 1, because fit_segmented() only returns (D, f, Dp). That may well be the real root cause here: if the pre-fit estimated S0, it would normally be in range and the clamp would rarely fire at all.

Worth noting the current tests can't catch this - test_bounds only checks that results land inside the supplied bounds, and fit_bayesian passes those bounds to scipy.minimize, so that's true regardless of x0.

A few questions before I open a PR, since these interact and I'd rather not guess:

1. Is [0] just a typo for [i], or was something else intended?
2. Should the clamp apply to S0 at all? The full-volume path deliberately skips it (if i == 3: fit_results[i] = np.random.normal(1, 0.2, ...)).
3. Should the segmented pre-fit estimate S0 instead of hard-coding it to 1, the way IVIMNET does?

Happy to submit just the [0] → [i] change, or something broader covering 2 and 3 - whichever you'd prefer.

Tested on main @ f91b53c, Python 3.11.3, scipy 1.17.0.

Are you working on this?

Yes

[Devguru-codes]

Hello, @oliverchampion sir. Thank you.

[oliverchampion]

Hey Devguru,
Thanks for pointing these things out.

Feel free to fix (1).
For (2) there is a difference between the array fitting and the voxel-wise fitting. In the voxel-wise fitting, the prior is defined in the initailisation phase, and the values are only used for initial guess. In this case, initial guesses should fall within the fit bound (flat prior bounds). When a full array is fitted in fit volume, I use a fit to the full image to determine the priors. In this case, the prior is not flat. For S0, I hardcoded a Gaussian prior. As segmented fit always finds S0=1, I need to build my own prior for S0.
And (3), segmented fit, per definition, has S0=1. It is at its heart. Otherwise it cannot estimate f.