Rounding error in add_gradients leads to incorrect gradient waveforms

[felixlandmeyer]

Describe the bug
(Note: This is my first issue report in a public repository, I apologise in advance for any mistakes.)
1.
Line 205 in add_gradients returns an error because of dimension inconsistency.

I changed line 205 from
waveforms[ii] = np.concatenate(([t_delay], waveforms[ii]))
to
waveforms[ii] = np.concatenate((t_delay, waveforms[ii]))
to resolve the error (waveforms[ii] should always have shape (N,) - I think-)

When adding gradients, sometimes too many gradient raster times are added as a delay in add_gradients.
The returned gradient waveform has a jump to zero and back when connecting the individual gradients.
With add_gradients():

Gradients in single blocks:

To Reproduce

run this script:

from math import ceil, isclose
from math import pi as PI
import numpy as np
import pypulseq as pp


# Set up system
system = pp.Opts(
    grad_unit='mT/m', max_grad=50, slew_unit='T/m/s',
    max_slew=180, grad_raster_time=10e-6)

Seq = pp.Sequence(system=system)
Seq1 = pp.Sequence(system=system)

amplitudes = np.linspace(0.003, 0.04, 15)  # [T/m]
readout_dur = 4e-3  # [s]
readout_tps = np.arange(0, readout_dur, system.grad_raster_time) 

for ampl in amplitudes:
    # Create readout gradient
    waveform = ampl * np.cos(2*PI/readout_dur * readout_tps)
    g_read = pp.make_arbitrary_grad(
        channel='x', waveform=waveform*system.gamma, system=system)

    # Create rampup gradient
    ramp_dur = g_read.waveform[0]/system.max_slew
    ramp_dur = ceil(
        ramp_dur/system.grad_raster_time
    ) * system.grad_raster_time
    
    g_ramp = pp.make_extended_trapezoid(
        channel='x', amplitudes=[0, g_read.waveform[-1]],
        times=[0, ramp_dur], system=system)
    
    # Create moment nulling gradient for rampup
    g_nullmom = pp.make_trapezoid(
        channel='x', area=-g_ramp.area, system=system)

    single_grads_duration = (
        pp.calc_duration(g_nullmom)
        + pp.calc_duration(g_ramp)
        + pp.calc_duration(g_read)
    )

    # Add native gradients to sequence    
    Seq.add_block(g_nullmom)
    Seq.add_block(g_ramp)
    Seq.add_block(g_read)
    Seq.add_block(pp.make_delay(1e-3))
    
    # Add gradients
    g_ramp.delay = pp.calc_duration(g_nullmom)
    g_read.delay = pp.calc_duration(g_ramp)
    g_new = pp.add_gradients(
        grads=[g_nullmom, g_ramp, g_read],
        max_slew=np.inf,  # Otherwise slew rate conflict 
        system=system
    )

    # Add gradients to sequences
    Seq1.add_block(g_new)
    Seq1.add_block(pp.make_delay(1e-3))

    # Check if gradient timing has changed
    if single_grads_duration == pp.calc_duration(g_new):
        print(f'Duration does not change with added gradients for amplitude {ampl:.4f}')
    else:
        print(f'Duration DOES change with added gradients for amplitude {ampl:.4f}')
        # add_gradients() calc:
        t_delay = np.arange(0, g_read.delay-0.0, step=system.grad_raster_time)
        if t_delay.size != round(g_read.delay/system.grad_raster_time):
            print('\t GRT added for delay in add_gradients():', t_delay.size)
            print('\t Delay to add in units of GRT:', g_read.delay/system.grad_raster_time)

# Show sequence diagramms
dur_seq = Seq.duration()[0]
dur_seq1 = Seq.duration()[0]

Seq.plot(
    show_blocks=False, time_range=(0, dur_seq),
    time_disp='ms', grad_disp='mT/m', plot_now=True)

Seq1.plot(
    show_blocks=False, time_range=(0, dur_seq),
    time_disp='ms', grad_disp='mT/m', plot_now=True)

Expected behavior
Added gradients without a jump at the connection point.

Desktop (please complete the following information):

• OS: macOS
• OS Version: Sonoma 14.5
• pypulseq version: 1.4.2 and 1.4.0

Additional context
The numpy.arange() can add one zero too many to the gradient waveform in some cases due to rounding errors.
In the script above, this happens to some of the readouts (e.g. readout two and five).

I'm still kind of a newbie to pulseq so it might be an issue with my code instead of pypulseq itself.
Other than that, I have fixed the rounding issue in my pypulseq installation and can create a PR if you want.

[btasdelen]

Hi @felixlandmeyer. Thanks for your example snippet and pointers, I was able to reproduce and understand the bug quickly.

1. That looks like an extraneous dimension added to the time vector, removing the brackets seems to be the correct solution.

2. Loosely related to Hacky fix for waveforms_and_times rf append floating point issue. #114 and Making internal time representation integer in units of nanoseconds. #129. I have two solutions to this, and open to the suggestion on which one to pick.

a. Using t_delay = np.arange(0, round(g.delay - common_delay, 6), step=system.grad_raster_time)
This assumes dwell time to be larger than 1us, which is reasonable for most systems, but an assumption nevertheless.

b. Using t_delay = np.arange(0, g.delay - common_delay - system.grad_raster_time/2, step=system.grad_raster_time)
This subtracts half the dwell time to make sure the end point falls between two raster times, and does not rely on the floating point number being exact.

I am more inclined towards b.

[felixlandmeyer]

Hi @btasdelen, thank you for the quick reply!

I used solution c:
t_delay = round((g.delay-common_delay)/system.grad_raster_time)
w_delay = np.zeros(t_delay)
waveforms[ii] = np.concatenate((w_delay, waveforms[ii]))

I was a little confused by the original version because it adds a time array (t_delay) to a gradient amplitude array (waveform).
As I understand the purpose, It should just add zeros (w_delay) to the waveform to create a delay.
Adding t_delay does not add amplitudes of 0 to the waveform but some values above 0.
Since the internal gradient representation is several tousands of Hz/m this does not have a large effect, but is technically not correct as far as I understand it.

[btasdelen]

You are totally right, I missed what the code tries to do as I focused on the single line.

I tracked the reason for the issue to here in Matlab Pulseq:

https://github.com/pulseq/pulseq/blob/f4fe1ec684968c2ce557b0883b32c14d6127d169/matlab/%2Bmr/addGradients.m#L173-L174

        t_delay = 0:opt.system.gradRasterTime:g.delay-common_delay-opt.system.gradRasterTime;
        waveforms{ii} = [t_delay*0 waveforms{ii}];

Python version mimicked that, but there is a missing *0 in the Python version. The way you do it is cleaner.