Residual per-shot signal drift past ~270 s (same FP collapse as #780, two more sites)

[ArthurAllilaire]

What happened?

#780 fixed the rising-edge FP-epsilon collapse, but the same hazard remains at the time-rebasing sites. _separate_closing_knot! separates the closing knot by a fixed 1e-14 in block-relative coordinates; get_samples and get_variable_times then rebase to absolute time with T0 .+ t. Once eps(T0 + t[end]) > 1e-14 (T0 ≳ 270 s for ms-scale events), the gap collapses, the closing knot equals the block-end knot, and the knots are removed by unique! leading to the same errors as last time.

Fix

Leave the block-relative _separate_closing_knot! untouched, and re-separate the closing knot on the rebased absolute vectors in get_samples / get_variable_times with an idempotent adaptive gap max(MIN_RISE_TIME, eps(t[end])) that widens past 1e-14 once t[end] is large.

Reproducer

24 identical IR-90° shots, TR = 15 s, no gradients — every |signal| should be identical, but shot 18 (sim_time = 270 s) jumps +52 %:

shot 1..17 (≤255 s): |signal| = 0.476267   stable
shot 18    (270 s):  |signal| = 0.722479   +52 %
shot 19..24(≥285 s): |signal| = 0.677159   new wrong plateau

using KomaMRI, Suppressor, Printf

const B1_AMP  = 20e-6     # 20 µT → 180° in 1 ms
const T_RF    = 1e-3
const T_ADC   = 1e-3
const TI      = 3.0
const TR      = 15.0
const N_SHOTS = 24

obj = Phantom(name="spin", x=[0.0], y=[0.0], z=[0.0],
              T1=[1.0], T2=[1.0], T2s=[1.0], ρ=[1.0], Δw=[0.0])

function rf_block(α, T)
    amp = B1_AMP * (α / π)
    Sequence(reshape([Grad(0.0, T), Grad(0.0, T), Grad(0.0, T)], 3, 1),
             reshape([RF(amp, T)], 1, 1), [ADC(0, 0.0)])
end

function adc_block(T)
    Sequence(reshape([Grad(0.0, T), Grad(0.0, T), Grad(0.0, T)], 3, 1),
             reshape([RF(0.0, T)], 1, 1), [ADC(1, T)])
end

function build_seq(N; TI, TR)
    seq = Sequence()
    tr_pad = TR - TI - T_RF - T_ADC
    for _ in 1:N
        seq += rf_block(π,   T_RF)        # 180° inversion
        seq += Delay(TI - T_RF)           # TI
        seq += rf_block(π/2, T_RF)        # 90° excitation
        seq += adc_block(T_ADC)           # 1-sample ADC
        tr_pad > 1e-9 && (seq += Delay(tr_pad))
    end
    seq
end

seq = build_seq(N_SHOTS; TI=TI, TR=TR)
raw = Suppressor.@suppress simulate(obj, seq, Scanner();
                                    sim_params=Dict{String,Any}("Nthreads" => 1))
vals = [abs(p.data[1, 1]) for p in raw.profiles]
ref  = vals[2]

println("  shot │ sim_time [s] │ |signal|   │ dev vs shot 2")
for (k, v) in enumerate(vals)
    dev  = 100 * (v - ref) / ref
    flag = abs(dev) > 1.0 ? "  ← DRIFT" : ""
    @printf("  %4d │ %12g │ %10.6g │ %+7.3f%%%s\n", k, k*TR, v, dev, flag)
end

Environment

OS x86_64-linux-gnu
Julia 1.12.6
KomaMRIPlots 0.11.0
KomaMRIFiles 0.10.4
KomaMRICore 0.11.2
KomaMRIBase 0.11.2