Repository: imr-framework/pypulseq

Issue:
Max slew check in trapezoidal gradient check fails because of floating point accuracy
When generating gradients at max slew rate, checks on these lines for [risetime](https://github.com/imr-framework/pypulseq/blame/234f4210773016d3ba996d2cb6e4527abcd32e52/src/pypulseq/make_trapezoid.py#L231) and [falltime](https://github.com/imr-framework/pypulseq/blame/234f4210773016d3ba996d2cb6e4527abcd32e52/src/pypulseq/make_trapezoid.py#L236) fail due to floating point error.

I'll propose a PR + test to sort.

Discussion:
@wtclarke: Looks like a similar issue has been handled here https://github.com/imr-framework/pypulseq/blob/master/src/pypulseq/make_extended_trapezoid.py#L133

Solution patch:
diff --git a/src/pypulseq/make_arbitrary_grad.py b/src/pypulseq/make_arbitrary_grad.py
index 1ca008e..445041e 100644
--- a/src/pypulseq/make_arbitrary_grad.py
+++ b/src/pypulseq/make_arbitrary_grad.py
@@ -3,6 +3,7 @@ from typing import Union
 
 import numpy as np
 
+from pypulseq import eps
 from pypulseq.opts import Opts
 from pypulseq.utils.tracing import trace, trace_enabled
 
@@ -76,9 +77,9 @@ def make_arbitrary_grad(
         raise ValueError(f'Invalid channel. Must be one of x, y or z. Passed: {channel}')
 
     slew_rate = np.diff(waveform) / system.grad_raster_time
-    if max(abs(slew_rate)) >= max_slew:
+    if max(abs(slew_rate)) > max_slew + eps:
         raise ValueError(f'Slew rate violation {max(abs(slew_rate)) / max_slew * 100}')
-    if max(abs(waveform)) >= max_grad:
+    if max(abs(waveform)) > max_grad + eps:
         raise ValueError(f'Gradient amplitude violation {max(abs(waveform)) / max_grad * 100}')
 
     if first is None:
diff --git a/src/pypulseq/make_extended_trapezoid_area.py b/src/pypulseq/make_extended_trapezoid_area.py
index c3ae152..6dfa958 100644
--- a/src/pypulseq/make_extended_trapezoid_area.py
+++ b/src/pypulseq/make_extended_trapezoid_area.py
@@ -3,6 +3,7 @@ from typing import Tuple, Union
 
 import numpy as np
 
+from pypulseq import eps
 from pypulseq.make_extended_trapezoid import make_extended_trapezoid
 from pypulseq.opts import Opts
 from pypulseq.utils.cumsum import cumsum
@@ -14,6 +15,7 @@ def make_extended_trapezoid_area(
     channel: str,
     grad_start: float,
     grad_end: float,
+    convert_to_arbitrary: bool = False,
     system: Union[Opts, None] = None,
 ) -> Tuple[SimpleNamespace, np.array, np.array]:
     """Make the shortest possible extended trapezoid for given area and gradient start and end point.
@@ -28,6 +30,8 @@ def make_extended_trapezoid_area(
         Starting non-zero gradient value.
     grad_end : float
         Ending non-zero gradient value.
+    convert_to_arbitrary : bool, default=False
+        Boolean flag to enable converting the extended trapezoid gradient into an arbitrary gradient.
     system: Opts, optional
         System limits.
 
@@ -84,7 +88,7 @@ def make_extended_trapezoid_area(
 
         # Check if gradient amplitude exceeds max_grad, if so, adjust ramp
         # times for a trapezoidal gradient with maximum slew rate.
-        if grad_start + ramp_up_time * max_slew * raster_time > max_grad:
+        if grad_start + ramp_up_time * max_slew * raster_time > max_grad + eps:
             ramp_up_time = round(_calc_ramp_time(grad_start, max_grad) / raster_time)
             ramp_down_time = round(_calc_ramp_time(grad_end, max_grad) / raster_time)
         else:
@@ -102,7 +106,7 @@ def make_extended_trapezoid_area(
 
         # Check if gradient amplitude exceeds -max_grad, if so, adjust ramp
         # times for a trapezoidal gradient with maximum slew rate.
-        if grad_start - ramp_up_time * max_slew * raster_time < -max_grad:
+        if grad_start - ramp_up_time * max_slew * raster_time < -max_grad - eps:
             ramp_up_time = round(_calc_ramp_time(grad_start, -max_grad) / raster_time)
             ramp_down_time = round(_calc_ramp_time(grad_end, -max_grad) / raster_time)
         else:
@@ -214,7 +218,9 @@ def make_extended_trapezoid_area(
         times = cumsum(0, time_ramp_up, time_ramp_down)
         amplitudes = np.array([grad_start, grad_amp, grad_end])
 
-    grad = make_extended_trapezoid(channel=channel, system=system, times=times, amplitudes=amplitudes)
+    grad = make_extended_trapezoid(
+        channel=channel, amplitudes=amplitudes, convert_to_arbitrary=convert_to_arbitrary, system=system, times=times
+    )
 
     # Overwrite trace
     if trace_enabled():
@@ -223,4 +229,4 @@ def make_extended_trapezoid_area(
     if not abs(grad.area - area) < 1e-8:
         raise ValueError(f'Could not find a solution for area={area}.')
 
-    return grad, np.array(times), amplitudes
+    return grad, grad.tt, grad.waveform
diff --git a/src/pypulseq/make_sigpy_pulse.py b/src/pypulseq/make_sigpy_pulse.py
index 382d59c..4819edf 100644
--- a/src/pypulseq/make_sigpy_pulse.py
+++ b/src/pypulseq/make_sigpy_pulse.py
@@ -187,11 +187,6 @@ def sigpy_n_seq(
         if rfp.delay < (gz.rise_time + gz.delay):
             rfp.delay = gz.rise_time + gz.delay
 
-    if rfp.ringdown_time > 0:
-        t_fill = np.arange(1, round(rfp.ringdown_time / 1e-6) + 1) * 1e-6
-        rfp.t = np.concatenate((rfp.t, rfp.t[-1] + t_fill))
-        rfp.signal = np.concatenate((rfp.signal, np.zeros(len(t_fill))))
-
     # Following 2 lines of code are workarounds for numpy returning 3.14... for np.angle(-0.00...)
     negative_zero_indices = np.where(rfp.signal == -0.0)
     rfp.signal[negative_zero_indices] = 0
diff --git a/src/pypulseq/make_trapezoid.py b/src/pypulseq/make_trapezoid.py
index d41a05a..a0050ca 100644
--- a/src/pypulseq/make_trapezoid.py
+++ b/src/pypulseq/make_trapezoid.py
@@ -20,7 +20,7 @@ def calculate_shortest_params_for_area(area: float, max_slew: float, max_grad: f
     effective_time = rise_time
 
     # Adjust for max gradient constraint
-    if abs(amplitude) > max_grad:
+    if abs(amplitude) > max_grad + eps:
         effective_time = math.ceil(abs(area) / max_grad / grad_raster_time) * grad_raster_time
         amplitude = area / effective_time
         rise_time = math.ceil(abs(amplitude) / max_slew / grad_raster_time) * grad_raster_time
@@ -225,15 +225,15 @@ def make_trapezoid(
     if rise_time is None and fall_time is None:
         rise_time = fall_time = calculate_shortest_rise_time(amplitude2, max_slew, system.grad_raster_time)
 
-    if abs(amplitude2) > max_grad:
+    if abs(amplitude2) > max_grad + eps:
         raise ValueError(f'Refined amplitude ({abs(amplitude2):0.0f} Hz/m) is larger than max ({max_grad:0.0f} Hz/m).')
 
-    if abs(amplitude2) / rise_time > max_slew:
+    if abs(amplitude2) / rise_time > max_slew + eps:
         raise ValueError(
             f'Refined slew rate ({abs(amplitude2) / rise_time:0.0f} Hz/m/s) for ramp up is larger than max ({max_slew:0.0f} Hz/m/s).'
         )
 
-    if abs(amplitude2) / fall_time > max_slew:
+    if abs(amplitude2) / fall_time > max_slew + eps:
         raise ValueError(
             f'Refined slew rate ({abs(amplitude2) / fall_time:0.0f} Hz/m/s) for ramp down is larger than max ({max_slew:0.0f} Hz/m/s).'
         )


Test patch:
diff --git a/tests/test_make_extended_trapezoid_area.py b/tests/test_make_extended_trapezoid_area.py
index 7af8279..8f0007c 100644
--- a/tests/test_make_extended_trapezoid_area.py
+++ b/tests/test_make_extended_trapezoid_area.py
@@ -46,7 +46,7 @@ test_zoo = [
 ]
 
 
-@pytest.mark.parametrize('grad_start,grad_end,area', test_zoo)
+@pytest.mark.parametrize('grad_start, grad_end, area', test_zoo)
 def test_make_extended_trapezoid_area(grad_start, grad_end, area):
     g, _, _ = make_extended_trapezoid_area(
         channel='x', grad_start=grad_start, grad_end=grad_end, area=area, system=system
@@ -71,7 +71,7 @@ test_zoo_random = [
 ]
 
 
-@pytest.mark.parametrize('grad_start,grad_end,area', test_zoo_random)
+@pytest.mark.parametrize('grad_start, grad_end, area', test_zoo_random)
 def test_make_extended_trapezoid_area_random_cases(grad_start, grad_end, area):
     g, _, _ = make_extended_trapezoid_area(
         channel='x', grad_start=grad_start, grad_end=grad_end, area=area, system=system
@@ -85,6 +85,34 @@ def test_make_extended_trapezoid_area_random_cases(grad_start, grad_end, area):
     assert slew_ok, 'Maximum slew rate violated'
 
 
+@pytest.mark.parametrize('grad_start, grad_end, area', test_zoo_random)
+def test_make_extended_trapezoid_area_convert_to_arb(grad_start, grad_end, area):
+    g, _, _ = make_extended_trapezoid_area(
+        channel='x', grad_start=grad_start, grad_end=grad_end, area=area, system=system
+    )
+
+    g_arb, _, _ = make_extended_trapezoid_area(
+        channel='x', grad_start=grad_start, grad_end=grad_end, area=area, convert_to_arbitrary=True, system=system
+    )
+
+    grad_ok = all(abs(g.waveform) <= system.max_grad)
+    slew_ok = all(abs(np.diff(g.waveform) / np.diff(g.tt)) <= system.max_slew)
+
+    grad_arb_ok = all(abs(g_arb.waveform) <= system.max_grad)
+    slew_arb_ok = all(abs(np.diff(g_arb.waveform) / np.diff(g_arb.tt)) <= system.max_slew)
+
+    assert pytest.approx(g.area) == g_arb.area, 'Area of extended trapz and arb gradient do not match'
+    assert pytest.approx(g.shape_dur) == g_arb.shape_dur, 'Duration of extended trapz and arb gradient do not match'
+    assert g.tt.shape[0] <= g_arb.tt.shape[0], (
+        'Extended trapezoid should have less or equal number of points than arb gradient'
+    )
+    assert g.waveform.shape[0] <= g_arb.waveform.shape[0], (
+        'Extended trapezoid should have less or equal number of points than arb gradient'
+    )
+    assert grad_ok == grad_arb_ok, 'Gradient strength violation between extended trapz and arb gradient'
+    assert slew_ok == slew_arb_ok, 'Slew rate violation between extended trapz and arb gradient'
+
+
 random.seed(0)
 test_zoo_random = [
     (
@@ -96,7 +124,7 @@ test_zoo_random = [
 ]
 
 
-@pytest.mark.parametrize('grad_start,grad_end,grad_amp', test_zoo_random)
+@pytest.mark.parametrize('grad_start, grad_end, grad_amp', test_zoo_random)
 def test_make_extended_trapezoid_area_recreate(grad_start, grad_end, grad_amp):
     def _to_raster(time: float) -> float:
         return np.ceil(time / system.grad_raster_time) * system.grad_raster_time
diff --git a/tests/test_sigpy.py b/tests/test_sigpy.py
index 0d17d82..71ab8a4 100644
--- a/tests/test_sigpy.py
+++ b/tests/test_sigpy.py
@@ -62,7 +62,7 @@ def test_slr():
     )
 
     # Check that the number of samples in the pulse is correct
-    assert rfp.signal.shape[0] == pytest.approx((duration + system.rf_ringdown_time) / system.rf_raster_time)
+    assert rfp.signal.shape[0] == pytest.approx(duration / system.rf_raster_time)
 
     # Check that the pulse can be added to a PyPulseq Sequence
     seq = pp.Sequence(system=system)
@@ -125,7 +125,7 @@ def test_sms():
     )
 
     # Check that the number of samples in the pulse is correct
-    assert rfp.signal.shape[0] == pytest.approx((duration + system.rf_ringdown_time) / system.rf_raster_time)
+    assert rfp.signal.shape[0] == pytest.approx(duration / system.rf_raster_time)
 
     # Check that the pulse can be added to a PyPulseq Sequence
     seq = pp.Sequence(system=system)