calculate_shortest_params_for_area double rounding causes longer duration

[nathanielatom]

Bug Description
The helper function calculate_shortest_params_for_area raises an AssertionError for a shortest possible trapezoid (well triangle) when make_trapezoid is given area and duration that is the shortest possible.

To Reproduce
Use the EPI example in the documentation, https://github.com/imr-framework/pypulseq/blob/master/examples/scripts/write_epi.py
and only change the FOV to fov = 240e-3 m and max slew rate to max_slew=200 T/m/s. With just these two reasonable changes, the calculate_shortest_params_for_area raises an AssertionError:

Requested area is too large for this gradient. Minimum required duration is 60 us. The passed duration is 50 us.

Expected behavior
The EPI example calculates the shorted possible gradient triangle using the same method as the calculate_shortest_params_for_area, however it only rounds the duration up to system.grad_raster_time once, at the end of the calculation for total duration of the triangle; whereas calculate_shortest_params_for_area rounds up the rise_time on its own, then sets the fall_time equal to the rise_time, effectively rounding up twice, and leading to a longer than shortest possible duration and the error.

Desktop (please complete the following information):

• OS: macOS Sequoia
• OS Version: 15.3.1 (24D70)
• python version: 3.9.10
• pypulseq version: 1.4.2.post1

I'd be happy to help with a PR at some point when I have availability, but not sure the preferred approach.

[btasdelen]

I am able to reproduce the issue. The question is, which is the correct behavior? You said calculate_shortest_params_for_area leads to a longer than the shortest possible duration. But, the shortest duration here is 50 us, leading to 25 us rise and fall times, which are not on the raster. So, I would argue that actually the shortest possible duration is 60 us, and the duration calculated at EPI example is incorrect. What do you think?

[nathanielatom]

Yeah exactly, that's why I wasn't sure about the preferred approach of modifying the EPI example (and any others) or calculate_shortest_params_for_area.

I can see an argument where a 50 us duration triangle is sampled with 5 points: [0, g/2, g, g/2, 0] where g is the gradient amplitude/height of the triangle. But here it's sampled with an odd number of points, so can't be divided into rise and fall times which are multiples of the raster. Sampling a triangle with 6 points misses the peak and gives an actual trapezoid again. So maybe the simplest thing to do is create another function make_triangle that only returns the shortest possible triangle given area or amplitude and slew rate or system.

Let me know what you think!

[btasdelen]

I see what you mean. In this case, I think using make_extended_trapezoid_area might be more suitable here. Using:

gy, gy_times, gy_amps = pp.make_extended_trapezoid_area(channel='y', system=system, area=delta_k, grad_start=0, grad_end=0)

returns:

gy_times
>>> array([0.e+00, 2.e-05, 5.e-05])
gy_amps
>>> array([     0.        , 166666.66666667,      0.        ])

My impression is calculate_shortest_params_for_area can be modified to handle this case. The actual problem is due to the assumption of duration = rise_time + flat_time + fall_time, make_trapezoid can not handle triangles, so that also needs to be handled. You are welcome to have a crack at it.

[btasdelen]

This is how Matlab Pulseq handles it, and we are basically missing this:

https://github.com/pulseq/pulseq/blob/6de38ff11bd1f712bd036a77887def43b4721a59/matlab/%2Bmr/makeTrapezoid.m#L171-L174