Repository: imr-framework/pypulseq

Issue:
make_trapezoid() changes rise_time even if passed as input argument
**Describe the bug**
Creating a gradient using make_trapezoid(channel = 'x',  area = 62, rise_time = 200e-6) results in a gradient event with a different rise_time than provided. 

**To Reproduce**
```
import pypulseq as pp

grad_area = 62
dG = 200e-6

# Create a gradient with a given area and rise time
grad_y1 = pp.make_trapezoid(channel = 'x',  area = grad_area, rise_time = dG)

# Rise time got changed!
print(grad_y1.rise_time)
pe_duration = pp.calc_duration(grad_y1)

# Create a second gradient with the same area and duration triggers error
grad_y2 = pp.make_trapezoid(channel='x',  area = grad_area, duration = pe_duration, rise_time = dG)
```

**Expected behavior**
I would expect that if I provide a specific rise_time that it will not change. 


**Desktop (please complete the following information):**
 - OS: [e.g. Win11]
- `pypulseq` version: [e.g. 1.4.2]

Discussion:
@btasdelen: For the first case, the combination of `area` + `rise_time` is not a valid combination of inputs. Valid combinations as per the documentation are:

```
    The user must supply as a minimum one of the following sets:
    - area
    - amplitude and duration
    - flat_time and flat_area
    - flat_time and amplitude
    - flat_time, area and rise_time
    Additional options may be supplied with the above.
``` 

That being said, wording of the documentation is as if any additional input to the `area` will be a valid combination, so at minimum, wording should be fixed. Also, these input pairs slipped from the input sanity checks, that should be fixed. 

I see two problems:

1. There is a check if `fall_time` is provided, and `rise_time` is not provided, we return an error. But, there is no check for the vice versa. Input logic around this is very convoluted.

2. We have a check if `flat_time` or `duration` is provided, in both cases you shared, it is not. So we default to calculate the shortest grad from area:

``` python
            amplitude2, rise_time, flat_time, fall_time = calculate_shortest_params_for_area(area, max_slew, max_grad, system.grad_raster_time)
``` 

This function overrides the `rise_time`, hence the discrepancy.


For the second case, you provide 20 us for `rise_time`, and `duration`. `fall_time = rise_time` since you did not provide it. That makes the timing impossible. ~~The logic here seems wrong, if we are provided `duration` and `rise_time`, `fall_time=duration-rise_time`~~. Also, we should check if `duration > rise_time` and throw an error, which is your case here.

I will address issues here and create a PR. 

Solution patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index e880e37..6c042c8 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -3,7 +3,7 @@
 
 repos:
   - repo: https://github.com/pre-commit/pre-commit-hooks
-    rev: v4.6.0
+    rev: v5.0.0
     hooks:
       - id: check-added-large-files
         exclude: ^tests/expected_output/pulseq_calcAdcSeg.txt
@@ -16,14 +16,14 @@ repos:
       - id: mixed-line-ending
 
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.6.9
+    rev: v0.8.6
     hooks:
       - id: ruff # linter
         args: [--fix]
       - id: ruff-format # formatter
 
   - repo: https://github.com/crate-ci/typos
-    rev: v1.25.0
+    rev: dictgen-v0.3.1
     hooks:
       - id: typos
         exclude: ^examples/|paper.md
diff --git a/pyproject.toml b/pyproject.toml
index eb48243..341e6f9 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -120,7 +120,10 @@ exclude = ["examples/**"]
 # PyTest section
 [tool.pytest.ini_options]
 testpaths = ["tests"]
-filterwarnings = ["error"]
+filterwarnings = ["error",
+    # Suppress error in debugpy due to mpl deprecation to debug tests.
+    "ignore::matplotlib._api.deprecation.MatplotlibDeprecationWarning:pydev",
+    ]
 markers = [
   "matlab_seq_comp: comparison with matlab generated sequence",
   "sigpy: tests that require sigpy",
diff --git a/src/pypulseq/make_adiabatic_pulse.py b/src/pypulseq/make_adiabatic_pulse.py
index 4d9f5c6..315be5a 100644
--- a/src/pypulseq/make_adiabatic_pulse.py
+++ b/src/pypulseq/make_adiabatic_pulse.py
@@ -221,17 +221,8 @@ def make_adiabatic_pulse(
         rf.delay = rf.dead_time
 
     if return_gz:
-        if max_grad is not None:
-            max_grad_slice_select = max_grad
-        else:
-            # Set to zero, not None for compatibility with existing make_trapezoid
-            max_grad_slice_select = 0
-
-        if max_slew is not None:
-            max_slew_slice_select = max_slew
-        else:
-            # Set to zero, not None for compatibility with existing make_trapezoid
-            max_slew_slice_select = 0
+        max_grad_slice_select = max_grad
+        max_slew_slice_select = max_slew
 
         if pulse_type == 'hypsec':
             bandwidth = mu * beta / np.pi
diff --git a/src/pypulseq/make_arbitrary_grad.py b/src/pypulseq/make_arbitrary_grad.py
index 40bcdd7..1ca008e 100644
--- a/src/pypulseq/make_arbitrary_grad.py
+++ b/src/pypulseq/make_arbitrary_grad.py
@@ -81,10 +81,10 @@ def make_arbitrary_grad(
     if max(abs(waveform)) >= max_grad:
         raise ValueError(f'Gradient amplitude violation {max(abs(waveform)) / max_grad * 100}')
 
-    if not first:
+    if first is None:
         first = (3 * waveform[0] - waveform[1]) * 0.5  # linear extrapolation
 
-    if not last:
+    if last is None:
         last = (3 * waveform[-1] - waveform[-2]) * 0.5  # linear extrapolation
 
     grad = SimpleNamespace()
diff --git a/src/pypulseq/make_trapezoid.py b/src/pypulseq/make_trapezoid.py
index ff1ddaa..7a917bd 100644
--- a/src/pypulseq/make_trapezoid.py
+++ b/src/pypulseq/make_trapezoid.py
@@ -1,57 +1,72 @@
 import math
+import warnings
 from types import SimpleNamespace
-from typing import Union
-
-import numpy as np
+from typing import Literal, Union
 
+from pypulseq import eps
 from pypulseq.opts import Opts
 from pypulseq.utils.tracing import trace, trace_enabled
 
 
-def calculate_shortest_params_for_area(area, max_slew, max_grad, grad_raster_time):
+def calculate_shortest_params_for_area(area: float, max_slew: float, max_grad: float, grad_raster_time: float):
+    """Calculate the shortest possible rise_time, flat_time, and fall_time for a given area."""
+
+    # Calculate initial rise time constrained by max slew rate
     rise_time = math.ceil(math.sqrt(abs(area) / max_slew) / grad_raster_time) * grad_raster_time
-    if rise_time < grad_raster_time:  # Area was almost 0 maybe
-        rise_time = grad_raster_time
-    amplitude = np.divide(area, rise_time)  # To handle nan
-    t_eff = rise_time
+    rise_time = max(rise_time, grad_raster_time)
+
+    # Calculate initial amplitude
+    amplitude = area / rise_time
+    effective_time = rise_time
 
+    # Adjust for max gradient constraint
     if abs(amplitude) > max_grad:
-        t_eff = math.ceil(abs(area) / max_grad / grad_raster_time) * grad_raster_time
-        amplitude = area / t_eff
+        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
+        rise_time = max(rise_time, grad_raster_time)
 
-        if rise_time == 0:
-            rise_time = grad_raster_time
-
-    flat_time = t_eff - rise_time
+    # Calculate flat and fall times
+    flat_time = effective_time - rise_time
     fall_time = rise_time
 
     return amplitude, rise_time, flat_time, fall_time
 
 
+def calculate_shortest_rise_time(amplitude: float, max_slew: float, grad_raster_time: float):
+    """Calculate the shortest possible rise / fall time for a given amplitude and slew rate."""
+
+    return math.ceil(max(abs(amplitude) / max_slew, grad_raster_time) / grad_raster_time) * grad_raster_time
+
+
 def make_trapezoid(
     channel: str,
-    amplitude: float = 0,
+    amplitude: Union[float, None] = None,
     area: Union[float, None] = None,
     delay: float = 0,
-    duration: float = 0,
-    fall_time: float = 0,
+    duration: Union[float, None] = None,
+    fall_time: Union[float, None] = None,
     flat_area: Union[float, None] = None,
-    flat_time: float = -1,
-    max_grad: float = 0,
-    max_slew: float = 0,
-    rise_time: float = 0,
+    flat_time: Union[float, None] = None,
+    max_grad: Union[float, None] = None,
+    max_slew: Union[float, None] = None,
+    rise_time: Union[float, None] = None,
     system: Union[Opts, None] = None,
 ) -> SimpleNamespace:
     """
     Create a trapezoidal gradient event.
 
-    The user must supply as a minimum one of the following sets:
+    The user must supply any of the following sets of parameters:
+    Area based:
     - area
-    - amplitude and duration
-    - flat_time and flat_area
-    - flat_time and amplitude
+    - area and duration
+    - area and duration and rise_time
     - flat_time, area and rise_time
+    Amplitude based:
+    - amplitude and duration
+    - amplitude and flat_time
+    Flat area based:
+    - flat_area and flat_time
     Additional options may be supplied with the above.
 
     See Also
@@ -63,23 +78,23 @@ def make_trapezoid(
     ----------
     channel : str
         Orientation of trapezoidal gradient event. Must be one of `x`, `y` or `z`.
-    amplitude : float, default=0
+    amplitude : float, default=None
         Peak amplitude (Hz/m).
     area : float, default=None
         Area (1/m).
     delay : float, default=0
         Delay in seconds (s).
-    duration : float, default=0
+    duration : float, default=None
         Duration in seconds (s). Duration is defined as rise_time + flat_time + fall_time.
-    fall_time : float, default=0
+    fall_time : float, default=None
         Fall time in seconds (s).
-    flat_area : float, default=0
+    flat_area : float, default=None
         Flat area (1/m).
-    flat_time : float, default=-1
+    flat_time : float, default=None
         Flat duration in seconds (s). Default is -1 to allow for triangular pulses.
-    max_grad : float, default=0
+    max_grad : float, default=None
         Maximum gradient strength (Hz/m).
-    max_slew : float, default=0
+    max_slew : float, default=None
         Maximum slew rate (Hz/m/s).
     rise_time : float, default=0
         Rise time in seconds (s).
@@ -98,6 +113,8 @@ def make_trapezoid(
         If requested area is too large for this gradient
         If `flat_time`, `duration` and `area` are not supplied.
         Amplitude violation
+    NotImplementedError
+        If an input set that might be valid but not implemented is passed.
     """
     if system is None:
         system = Opts.default
@@ -105,50 +122,42 @@ def make_trapezoid(
     if channel not in ['x', 'y', 'z']:
         raise ValueError(f'Invalid channel. Must be one of `x`, `y` or `z`. Passed: {channel}')
 
-    if max_grad <= 0:
+    if max_grad is None:
         max_grad = system.max_grad
 
-    if max_slew <= 0:
+    if max_slew is None:
         max_slew = system.max_slew
 
-    if rise_time <= 0:
-        rise_time = 0.0
-
-    if fall_time > 0:
-        if rise_time == 0:
-            raise ValueError(
-                'Invalid arguments. Must always supply `rise_time` if `fall_time` is specified explicitly.'
-            )
+    # If either of rise_time or fall_time is not provided, set it to the other.
+    rise_time = rise_time or fall_time
+    fall_time = fall_time or rise_time
+
+    # Check which one of area, flat_area and amplitude is provided, and determine the calculation path accordingly.
+    calc_path: Literal['area', 'flat_area', 'amplitude']
+    if area is not None and flat_area is None and amplitude is None:
+        calc_path = 'area'
+    elif area is None and flat_area is not None and amplitude is None:
+        calc_path = 'flat_area'
+    elif area is None and flat_area is None and amplitude is not None:
+        calc_path = 'amplitude'
+    elif area is None and flat_area is not None and amplitude is not None:
+        raise NotImplementedError('Flat Area + Amplitude input pair is not implemented yet.')
+    elif area is not None and flat_area is None and amplitude is not None:
+        raise NotImplementedError('Amplitude + Area input pair is not implemented yet.')
     else:
-        fall_time = 0.0
-
-    if area is None and flat_area is None and amplitude == 0:
         raise ValueError("Must supply either 'area', 'flat_area' or 'amplitude'.")
 
-    if flat_time != -1:
-        if amplitude != 0:
-            amplitude2 = amplitude
-        elif area is not None and rise_time > 0:
-            # We have rise_time, flat_time and area.
-            amplitude2 = area / (rise_time + flat_time)
-        elif flat_area is not None:
-            amplitude2 = flat_area / flat_time
-        else:
-            raise ValueError(
-                'When `flat_time` is provided, either `flat_area`, '
-                'or `amplitude`, or `rise_time` and `area` must be provided as well.'
-            )
-
-        if rise_time == 0:
-            rise_time = abs(amplitude2) / max_slew
-            rise_time = math.ceil(rise_time / system.grad_raster_time) * system.grad_raster_time
-            if rise_time == 0:
-                rise_time = system.grad_raster_time
-        if fall_time == 0:
-            fall_time = rise_time
-    elif duration > 0:
-        if amplitude == 0:
-            if rise_time == 0:
+    # Check if sufficient timing parameters are provided.
+    if flat_time is not None and flat_area is None and amplitude is None and (rise_time is None or area is None):
+        raise ValueError(
+            'When `flat_time` is provided, either `flat_area`, '
+            'or `amplitude`, or `rise_time` and `area` must be provided as well.'
+        )
+
+    if calc_path == 'area':
+        # We have area and duration.
+        if duration is not None and flat_time is None:
+            if rise_time is None:
                 _, rise_time, flat_time, fall_time = calculate_shortest_params_for_area(
                     area, max_slew, max_grad, system.grad_raster_time
                 )
@@ -161,49 +170,73 @@ def make_trapezoid(
                 dc = 1 / abs(2 * max_slew) + 1 / abs(2 * max_slew)
                 amplitude2 = (duration - math.sqrt(duration**2 - 4 * abs(area) * dc)) / (2 * dc)
             else:
-                if fall_time == 0:
+                if duration <= (rise_time + eps):
+                    raise ValueError('The `duration` is too short for the given `rise_time`.')
+
+                if fall_time is None:
                     fall_time = rise_time
+
                 amplitude2 = area / (duration - 0.5 * rise_time - 0.5 * fall_time)
                 possible = duration >= (rise_time + fall_time) and abs(amplitude2) <= max_grad
                 assert possible, (
                     f'Requested area is too large for this gradient. Probably amplitude is violated '
-                    f'{round(abs(amplitude) / max_grad * 100)}'
+                    f'{round(abs(amplitude2) / max_grad * 100)}'
                 )
-        else:
-            amplitude2 = amplitude
+            flat_time = duration - rise_time - fall_time
+            amplitude2 = area / (rise_time / 2 + fall_time / 2 + flat_time)
 
-        if rise_time == 0:
-            rise_time = math.ceil(abs(amplitude2) / max_slew / system.grad_raster_time) * system.grad_raster_time
-            if rise_time == 0:
-                rise_time = system.grad_raster_time
+        elif flat_time is not None:
+            if rise_time is None:
+                raise ValueError('Must supply `rise_time` when `area` and `flat_time` is provided.')
 
-        if fall_time == 0:
-            fall_time = rise_time
-        flat_time = duration - rise_time - fall_time
+            amplitude2 = area / (rise_time + flat_time)
 
-        if amplitude == 0:
-            # Adjust amplitude (after rounding) to match area
-            amplitude2 = area / (rise_time / 2 + fall_time / 2 + flat_time)
-    else:
-        if area is None:
-            raise ValueError('Must supply area or duration.')
         else:
-            # Find the shortest possible duration.
+            if rise_time is not None or fall_time is not None:
+                warnings.warn('Rise time and fall time is ignored when calculating the shortest duration from `area`.')
+
             amplitude2, rise_time, flat_time, fall_time = calculate_shortest_params_for_area(
                 area, max_slew, max_grad, system.grad_raster_time
             )
 
-    assert (
-        abs(amplitude2) <= max_grad
-    ), f'Refined amplitude ({abs(amplitude2):0.0f} Hz/m) is larger than max ({max_grad:0.0f} Hz/m).'
+    elif calc_path == 'flat_area':
+        # Check and raise invalid input sets
+        if duration is not None:
+            raise NotImplementedError('Flat Area + Duration input pair is not implemented yet.')
+        elif flat_time is not None:
+            amplitude2 = flat_area / flat_time
+
+    elif calc_path == 'amplitude':
+        if rise_time is None:
+            rise_time = abs(amplitude) / max_slew
+            rise_time = math.ceil(rise_time / system.grad_raster_time) * system.grad_raster_time
+            if rise_time == 0:
+                rise_time = system.grad_raster_time
+            fall_time = rise_time
+
+        amplitude2 = amplitude
+        if duration is not None and flat_time is None:
+            flat_time = duration - rise_time - fall_time
+        elif flat_time is not None and duration is None:
+            pass
+        else:
+            raise ValueError('Must supply area or duration.')
+
+    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:
+        raise ValueError(f'Refined amplitude ({abs(amplitude2):0.0f} Hz/m) is larger than max ({max_grad:0.0f} Hz/m).')
 
-    assert (
-        abs(amplitude2) / rise_time <= max_slew
-    ), 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) / rise_time > max_slew:
+        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).'
+        )
 
-    assert (
-        abs(amplitude2) / fall_time <= max_slew
-    ), 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).'
+    if abs(amplitude2) / fall_time > max_slew:
+        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).'
+        )
 
     grad = SimpleNamespace()
     grad.type = 'trap'


Test patch:
diff --git a/tests/test_make_trapezoid.py b/tests/test_make_trapezoid.py
index b7a2946..a420c7c 100644
--- a/tests/test_make_trapezoid.py
+++ b/tests/test_make_trapezoid.py
@@ -1,12 +1,7 @@
-"""Tests for the make_trapezoid module
-
-Will Clarke, University of Oxford, 2023
-"""
-
-from types import SimpleNamespace
+"""Tests for the make_trapezoid module"""
 
 import pytest
-from pypulseq import make_trapezoid
+from pypulseq import Opts, eps, make_trapezoid
 
 
 def test_channel_error():
@@ -14,11 +9,11 @@ def test_channel_error():
         make_trapezoid(channel='p')
 
 
-def test_falltime_risetime_error():
-    with pytest.raises(
-        ValueError, match=r'Invalid arguments. Must always supply `rise_time` if `fall_time` is specified explicitly.'
-    ):
-        make_trapezoid(channel='x', fall_time=10)
+# def test_falltime_risetime_error():
+#     with pytest.raises(
+#         ValueError, match=r'Invalid arguments. Must always supply `rise_time` if `fall_time` is specified explicitly.'
+#     ):
+#         make_trapezoid(channel='x', fall_time=10)
 
 
 def test_area_flatarea_amplitude_error():
@@ -60,10 +55,41 @@ def test_no_area_no_duration_error():
 def test_amplitude_too_large_error():
     errstr = r'Refined amplitude \(\d+ Hz/m\) is larger than max \(\d+ Hz/m\).'
 
-    with pytest.raises(AssertionError, match=errstr):
+    with pytest.raises(ValueError, match=errstr):
         make_trapezoid(channel='x', amplitude=1e10, duration=1)
 
 
+def test_duration_too_short_error():
+    errstr = 'The `duration` is too short for the given `rise_time`.'
+
+    with pytest.raises(ValueError, match=errstr):
+        make_trapezoid(channel='x', area=1, duration=0.1, rise_time=0.1)
+
+
+def test_notimplemented_input_pairs():
+    # flat_area + duration
+    with pytest.raises(NotImplementedError, match=r'Flat Area \+ Duration input pair is not implemented yet.'):
+        make_trapezoid(channel='x', flat_area=1, duration=1)
+    # flat_area + amplitude
+    with pytest.raises(NotImplementedError, match=r'Flat Area \+ Amplitude input pair is not implemented yet.'):
+        make_trapezoid(channel='x', flat_area=1, amplitude=1)
+    # area + amplitude
+    with pytest.raises(NotImplementedError, match=r'Amplitude \+ Area input pair is not implemented yet.'):
+        make_trapezoid(channel='x', area=1, amplitude=1)
+        # compare_trap_out(trap, 1, 2e-5, 0, 2e-5)
+
+
+def round2raster(x, raster=1e-5):
+    return round(x / raster) * raster
+
+
+def compare_trap_out(trap, amplitude, rise_time, flat_time, fall_time):
+    assert abs(trap.amplitude - amplitude) < eps
+    assert abs(trap.rise_time - rise_time) < eps
+    assert abs(trap.flat_time - flat_time) < eps
+    assert abs(trap.fall_time - fall_time) < eps
+
+
 def test_generation_methods():
     """Test minimum input cases
     Cover:
@@ -73,12 +99,41 @@ def test_generation_methods():
         - flat_time and amplitude
         - flat_time, area and rise_time
     """
-    assert isinstance(make_trapezoid(channel='x', area=1), SimpleNamespace)
 
-    assert isinstance(make_trapezoid(channel='x', amplitude=1, duration=1), SimpleNamespace)
-
-    assert isinstance(make_trapezoid(channel='x', flat_time=1, flat_area=1), SimpleNamespace)
+    opts = Opts()
+    # Amplitude specified
+    # amplitude + duration
+    trap = make_trapezoid(channel='x', amplitude=1, duration=1)
+    compare_trap_out(trap, 1, 1e-5, 1 - 2e-5, 1e-5)
+
+    # flat_time + amplitude
+    trap = make_trapezoid(channel='x', flat_time=1, amplitude=1)
+    compare_trap_out(trap, 1, 1e-5, 1, 1e-5)
+
+    # Flat area specified
+    # flat_area + flat_time
+    trap = make_trapezoid(channel='x', flat_time=1, flat_area=1)
+    compare_trap_out(trap, 1, 1e-5, 1, 1e-5)
+
+    # Area specified
+    # area
+    # triangle case
+    trap = make_trapezoid(channel='x', area=1)
+    compare_trap_out(trap, 50000, 2e-5, 0, 2e-5)
+    # trap case
+    trap = make_trapezoid(channel='x', area=opts.max_grad * 2)
+    time_to_max = round2raster(opts.max_grad / opts.max_slew)
+    compare_trap_out(
+        trap, opts.max_grad, time_to_max, (opts.max_grad * 2 - time_to_max * opts.max_grad) / opts.max_grad, time_to_max
+    )
 
-    assert isinstance(make_trapezoid(channel='x', flat_time=1, amplitude=1), SimpleNamespace)
+    # area + duration
+    trap = make_trapezoid(channel='x', area=1, duration=1)
+    compare_trap_out(trap, 1.00002, 2e-5, 1 - 4e-5, 2e-5)
 
-    assert isinstance(make_trapezoid(channel='x', flat_time=0.5, area=1, rise_time=0.1), SimpleNamespace)
+    # area + duration + rise_time
+    trap = make_trapezoid(channel='x', area=1, duration=1, rise_time=0.01)
+    compare_trap_out(trap, 1 / 0.99, 0.01, 0.98, 0.01)
+    # flat_time + area + rise_time
+    trap = make_trapezoid(channel='x', flat_time=0.5, area=1, rise_time=0.1)
+    compare_trap_out(trap, 1 / 0.6, 0.1, 0.5, 0.1)