Inaccurate RF pulse duration causing 25% amplitude drift on sequences (due to rise-edge collapse at t > 120s)

[ArthurAllilaire]

What happened?

Summary

Identical RF pulses played at different absolute times produce
different signal once cumulative simulated time crosses ~128 s. The
drift is silent — no warning, no NaN. Reproduced on master
@ d0d5c5fdd.

Root cause: the fixed constant MIN_RISE_TIME = 1e-14 in
KomaMRIBase/src/timing/TimeStepCalculation.jl is used to place a
rise-edge marker at t1 + MIN_RISE_TIME. Past the point where
eps(t)/2 > 1e-14 (i.e. t ≳ 128 s for Float64), that addition
rounds back onto t1 itself. sort!(unique!(t)) removes the
bit-equal duplicate, linear_interpolation returns the leading 0 at
t1, and the integrator loses exactly one Δt_rf of B1 area at
the leading edge of every affected pulse (5 % of pulse area with the
default Δt_rf = 5e-5, T_pulse = 1 ms).

For a single 90°-ish pulse, that 5 % area loss → 7.7° flip angle drop
→ sin(145.61°) / sin(153.27°) = 1.256 → +25.6 % |Mxy| jump. For
the 16-shot IR sequence below the compound cos(α_180)·sin(α_90)
ratio lands at ~1.20 — the observed 22 % shift.

Minimal reproducer

using KomaMRI
Trf = 1e-3
rf = Sequence(
    reshape([Grad(0, Trf), Grad(0, Trf), Grad(0, Trf)], 3, 1),
    reshape([RF(10e-6, Trf)], 1, 1), [ADC(0, 0.0)])
adc = Sequence(
    reshape([Grad(0, Trf), Grad(0, Trf), Grad(0, Trf)], 3, 1),
    reshape([RF(0.0, Trf)], 1, 1), [ADC(1, Trf)])

seq = Sequence()
for _ in 1:4
    seq += Delay(50.0); seq += rf; seq += adc
end

obj = Phantom(name="v", x=[0.], y=[0.], z=[0.],
              T1=[1.], T2=[1.], T2s=[1.], ρ=[1.], Δw=[0.])
raw = simulate(obj, seq, Scanner();
               sim_params=Dict{String,Any}("Nthreads"=>1))

for (k, p) in enumerate(raw.profiles)
    @info "shot $k" t_pulse_s = k*50.002  abs_signal = abs(p.data[1,1])
end

Output (master)

shot 1: t≈50 s   |Mxy| = 0.4497
shot 2: t≈100 s  |Mxy| = 0.4497
shot 3: t≈150 s  |Mxy| = 0.5646   ← +25.6 % jump
shot 4: t≈200 s  |Mxy| = 0.5646

The threshold is time-driven, not shot-count-driven: it sits at the
power-of-2 boundary above which eps(t) > 2·MIN_RISE_TIME (~128 s).

Mechanism in one walk

KomaMRIBase/src/timing/TimeStepCalculation.jl:106:

taux = points_from_key_times(sort([t1, t1 + ϵ, tc, t2 - ϵ, t2]); dt=Δt_rf)

The pair (t1, t1 + ϵ) encodes the rise-edge step. At small t both
entries are distinct Float64s and linear_interpolation evaluates to
B1 = 0 at t1, B1 = full at t1 + ϵ. Past ~128 s,
t1 + ϵ == t1 bit-exactly, sort!(unique!(t)) collapses the pair,
and the integrator integrates B1 = 0 over a full Δt_rf interval
instead of the intended ULP-sized one.

This is not the same as the cumsum/sum 1-ULP path disagreement at
block boundaries; that one is benign (Interpolations.deduplicate_knots!
handles it). The load-bearing failure is the bit-equal collapse of
the rise-edge marker.

Affected scenarios

Any multi-shot acquisition where total simulated time > ~128 s. There is currently no way to detect this from the simulator output. I am using KomaMRI for RL agent training, so I was testing longer sequences (250s).

Why existing tests didn't catch it

Every existing KomaMRICore test runs sequences of milliseconds total
duration (no Delay past tens of ms anywhere in the test suite). The
bug only manifests past 128 s of cumulative simulated time, so the
suite was structurally blind to time-scaling regressions in
get_variable_times.

Fix and PR

PR #780 implements the fix in get_variable_times by:

1. Wrapping the RF key-time list with
   Interpolations.deduplicate_knots!(_; move_knots=true) so any
   bit-equal duplicates introduced by t ± ϵ rounding are perturbed
   by 1 ULP via nextfloat. This reuses the same library function
   already used in 5 places elsewhere on event timelines
   (Sequence.jl:770/793, KeyValuesCalculation.jl:29/30/66).

2. Replacing the sequence-bookend t[1] - ϵ / t[end] + ϵ with
   max(MIN_RISE_TIME, 2·eps(t)) so the pad is distinct at any t
   while preserving a MIN_RISE_TIME floor for tiny t.

The fix is a no-op at normal time scales (markers stay at
t1 + ε / t2 - ε inside the pulse window) and only activates in
the regime where the original design was silently broken.

Regression tests added

Three layers in the PR — all fail on master, all pass after the fix:

• KomaMRIBase (unit): get_variable_times produces sub-Δt_rf
  marker pairs on both sides of every RF pulse, at any absolute t,
  with strictly increasing t.
• KomaMRIBase (semantic): the discrete ∫ B1·dt over the RF
  window equals A · T_pulse at any absolute t (rtol = 1e-6).
• KomaMRICore (integration): four identical 50 s-spaced shots
  produce bit-identical |Mxy|. Crosses the 128 s threshold between
  shots 2 and 3.

Versions reproduced on

• master @ d0d5c5fdd (2026-05-08, "Fix CLI world age dispatch Fix CLI world age dispatch #776")
• KomaMRI v0.10.3 / KomaMRICore v0.11.0 / KomaMRIBase v0.11.0 (Julia 1.11/1.12)
• KomaMRI v0.9.2 / KomaMRICore v0.9.6 (earlier check, same behaviour)
• CPU, single- and multi-thread
• Same drift bit-identical with precision = "f64" and precision = "f32" as time is always Float64.#
• Dependent on Δt_rf (loss scales linearly: 5e-5 → 25.6 %, 1e-5 → 5.3 %, 2e-6 → 1.1 %)

Environment

OS x86_64-linux-gnu
Julia 1.12.1
KomaMRIPlots 0.11.0
KomaMRIFiles 0.10.3
KomaMRICore 0.11.0
KomaMRIBase 0.11.0

[cncastillo]

I haven't experienced this error myself. Thanks for reporting!!

I will try to reproduce and check if the submitted solution seems correct.