From e27fcdf66bfd2f9339fe703d3a9cbde9c861ed61 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Tue, 26 Mar 2024 17:43:26 +0100 Subject: [PATCH 1/7] Fix make_extended_trapezoid area Improve timing/area calculation by ignoring raster_time at first and taking care of it in final calculation. This avoids issues by the discontinuities causing negative flat_times. Co-authored-by: Felix Zimmermann --- pypulseq/make_extended_trapezoid_area.py | 165 +++++++++++------------ 1 file changed, 78 insertions(+), 87 deletions(-) diff --git a/pypulseq/make_extended_trapezoid_area.py b/pypulseq/make_extended_trapezoid_area.py index 7537578..23c41cd 100644 --- a/pypulseq/make_extended_trapezoid_area.py +++ b/pypulseq/make_extended_trapezoid_area.py @@ -13,20 +13,18 @@ def make_extended_trapezoid_area( area: float, channel: str, grad_end: float, grad_start: float, system: Opts ) -> Tuple[SimpleNamespace, np.array, np.array]: - """ - Makes the shortest possible extended trapezoid with a given area which starts and ends (optionally) as non-zero - gradient values. + """Makes the shortest possible extended trapezoid for given area and gradient start and end point. Parameters ---------- + area : float + Area of extended trapezoid. channel : str Orientation of extended trapezoidal gradient event. Must be one of 'x', 'y' or 'z'. grad_start : float Starting non-zero gradient value. grad_end : float Ending non-zero gradient value. - area : float - Area of extended trapezoid. system: Opts System limits. @@ -35,100 +33,93 @@ def make_extended_trapezoid_area( grad : SimpleNamespace Extended trapezoid event. times : numpy.ndarray + Time points of the extended trapezoid. amplitude : numpy.ndarray + Amplitude values of the extended trapezoid. """ - SR = system.max_slew * 0.99 - - Tp = 0 - obj1 = ( - lambda x: ( - area - __testGA(x, 0, SR, system.grad_raster_time, grad_start, grad_end) + slew_rate = system.max_slew * 0.99 + max_grad = system.max_grad * 0.99 + raster_time = system.grad_raster_time + + def _to_raster(time: float) -> float: + return np.ceil(time / raster_time) * raster_time + + def _calc_area_and_times(grad_amp, flat_time, slew_rate, grad_start, grad_end): + time_ramp_up = _calc_ramp_up_time(grad_amp, slew_rate, grad_start) + time_ramp_down = _calc_ramp_down_time(grad_amp, slew_rate, grad_end) + area = _calc_area(grad_amp, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end) + return area + + def _calc_area(grad_amp, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end): + return ( + 0.5 * time_ramp_up * (grad_amp + grad_start) + + grad_amp * flat_time + + 0.5 * (grad_amp + grad_end) * time_ramp_down ) - ** 2 - ) - arr_res = [ - minimize(fun=obj1, x0=-system.max_grad, method="Nelder-Mead"), - minimize(fun=obj1, x0=0, method="Nelder-Mead"), - minimize(fun=obj1, x0=system.max_grad, method="Nelder-Mead"), - ] - arr_res = np.array([(*res.x, res.fun) for res in arr_res]) - Gp, obj1val = arr_res[:, 0], arr_res[:, 1] - i_min = np.argmin(obj1val) - Gp = Gp[i_min] - obj1val = obj1val[i_min] - - if obj1val > 1e-3 or abs(Gp) > system.max_grad: # Search did not converge - Gp = math.copysign(system.max_grad, Gp) - obj2 = ( - lambda x: ( - area - - __testGA(Gp, x, SR, system.grad_raster_time, grad_start, grad_end) - ) - ** 2 - ) - res2 = minimize(fun=obj2, x0=0, method="Nelder-Mead") - T, obj2val = *res2.x, res2.fun - assert obj2val < 1e-2 - Tp = math.ceil(T / system.grad_raster_time) * system.grad_raster_time - - # Fix the ramps - Tru = ( - math.ceil(abs(Gp - grad_start) / SR / system.grad_raster_time) - * system.grad_raster_time - ) - Trd = ( - math.ceil(abs(Gp - grad_end) / SR / system.grad_raster_time) - * system.grad_raster_time - ) - obj3 = lambda x: (area - __testGA1(x, Tru, Tp, Trd, grad_start, grad_end)) ** 2 + def _calc_ramp_up_time(grad_amp, slew_rate, grad_start): + return abs(grad_amp - grad_start) / slew_rate - res = minimize(fun=obj3, x0=Gp, method="Nelder-Mead") - Gp, obj3val = *res.x, res.fun - assert obj3val < 1e-3 # Did the final search converge? + def _calc_ramp_down_time(grad_amp, slew_rate, grad_end): + return abs(grad_amp - grad_end) / slew_rate - assert Tp >= 0 - - if Tp > 0: - times = cumsum(0, Tru, Tp, Trd) - amplitudes = [grad_start, Gp, Gp, grad_end] + # first try to find grad_amp without a flat time. + # for better convergence, we ignore the raster effect here and use a very small time step + # we will take care of raster effect later on + flat_time = 0 + obj1 = lambda x: (area - _calc_area_and_times(x, flat_time, slew_rate, grad_start, grad_end)) ** 2 + # try different starting points for the optimization + optimization_results = [ + minimize(fun=obj1, x0=-max_grad, method="Nelder-Mead"), + minimize(fun=obj1, x0=0, method="Nelder-Mead"), + minimize(fun=obj1, x0=max_grad, method="Nelder-Mead"), + ] + # choose the best result + optimization_results = np.array([(*res.x, res.fun) for res in optimization_results]) + grad_amp_results, obj1val_results = ( + optimization_results[:, 0], + optimization_results[:, 1], + ) + best_index = np.argmin(obj1val_results) + grad_amp = grad_amp_results[best_index] + obj1_value = obj1val_results[best_index] + + if obj1_value**0.5 > 1e-3 or abs(grad_amp) > max_grad: + # Optimization without flat_time did not converge. Try to find flat_time. + grad_amp = math.copysign(max_grad, grad_amp) + obj2 = lambda x: (area - _calc_area_and_times(grad_amp, x, slew_rate, grad_start, grad_end)) ** 2 + res2 = minimize(fun=obj2, x0=system.grad_raster_time, method="Nelder-Mead") + flat_time = res2.x[0] + # assert flat_time is larger than 0 and that the objective function is small + assert flat_time >= 0, "flat time is smaller than 0, this should not happen" + + # ensure times are on raster + flat_time = _to_raster(flat_time) + time_ramp_up = _to_raster(_calc_ramp_up_time(grad_amp, slew_rate, grad_start)) + time_ramp_down = _to_raster(_calc_ramp_down_time(grad_amp, slew_rate, grad_end)) + + # recalculate grad_amp while taking raster effect into account + obj3 = lambda x: (area - _calc_area(x, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end)) ** 2 + res = minimize(fun=obj3, x0=grad_amp, method="Nelder-Mead") + grad_amp = res.x[0] + + if flat_time > 0: + times = cumsum(0, time_ramp_up, flat_time, time_ramp_down) + amplitudes = [grad_start, grad_amp, grad_amp, grad_end] else: - Tru = ( - math.ceil(abs(Gp - grad_start) / SR / system.grad_raster_time) - * system.grad_raster_time - ) - Trd = ( - math.ceil(abs(Gp - grad_end) / SR / system.grad_raster_time) - * system.grad_raster_time - ) - - if Trd > 0: - if Tru > 0: - times = cumsum(0, Tru, Trd) - amplitudes = np.array([grad_start, Gp, grad_end]) + if time_ramp_down > 0: + if time_ramp_up > 0: + times = cumsum(0, time_ramp_up, time_ramp_down) + amplitudes = np.array([grad_start, grad_amp, grad_end]) else: - times = cumsum(0, Trd) + times = cumsum(0, time_ramp_down) amplitudes = np.array([grad_start, grad_end]) else: - times = cumsum(0, Tru) + times = cumsum(0, time_ramp_up) amplitudes = np.array([grad_start, grad_end]) - grad = make_extended_trapezoid( - channel=channel, system=system, times=times, amplitudes=amplitudes - ) - grad.area = __testGA1(Gp, Tru, Tp, Trd, grad_start, grad_end) + grad = make_extended_trapezoid(channel=channel, system=system, times=times, amplitudes=amplitudes) - assert abs(grad.area - area) < 1e-3 + assert abs(grad.area - area) < 1e-3, "Area of the gradient is not equal to the desired area. Optimization failed." return grad, np.array(times), amplitudes - - -def __testGA(Gp, Tp, SR, dT, Gs, Ge): - Tru = math.ceil(abs(Gp - Gs) / SR / dT) * dT - Trd = math.ceil(abs(Gp - Ge) / SR / dT) * dT - ga = __testGA1(Gp, Tru, Tp, Trd, Gs, Ge) - return ga - - -def __testGA1(Gp, Tru, Tp, Trd, Gs, Ge): - return 0.5 * Tru * (Gp + Gs) + Gp * Tp + 0.5 * (Gp + Ge) * Trd From 24a2f3ca06ca0e17cec22f08653ac8bac009ad69 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Fri, 5 Apr 2024 10:52:19 +0200 Subject: [PATCH 2/7] review comments (before general refactoring) change order of grad_end and grad_start arg define default system as None define single _calc_ramp_time function catch different cases for flat_time > 0 --- pypulseq/make_extended_trapezoid_area.py | 58 ++++++++++++++++-------- 1 file changed, 39 insertions(+), 19 deletions(-) diff --git a/pypulseq/make_extended_trapezoid_area.py b/pypulseq/make_extended_trapezoid_area.py index 23c41cd..68eb6ef 100644 --- a/pypulseq/make_extended_trapezoid_area.py +++ b/pypulseq/make_extended_trapezoid_area.py @@ -1,6 +1,6 @@ import math from types import SimpleNamespace -from typing import Tuple +from typing import Tuple, Union import numpy as np from scipy.optimize import minimize @@ -11,9 +11,13 @@ def make_extended_trapezoid_area( - area: float, channel: str, grad_end: float, grad_start: float, system: Opts -) -> Tuple[SimpleNamespace, np.array, np.array]: - """Makes the shortest possible extended trapezoid for given area and gradient start and end point. + area: float, + channel: str, + grad_start: float, + grad_end: float, + system: Union[Opts, None] = None, +) -> Tuple[SimpleNamespace, np.ndarray, np.ndarray]: + """Make (shortest possible) extended trapezoid for given area and gradient start and end point. Parameters ---------- @@ -25,7 +29,7 @@ def make_extended_trapezoid_area( Starting non-zero gradient value. grad_end : float Ending non-zero gradient value. - system: Opts + system: Opts, optional System limits. Returns @@ -37,6 +41,9 @@ def make_extended_trapezoid_area( amplitude : numpy.ndarray Amplitude values of the extended trapezoid. """ + if not system: + system = Opts() + slew_rate = system.max_slew * 0.99 max_grad = system.max_grad * 0.99 raster_time = system.grad_raster_time @@ -45,8 +52,8 @@ def _to_raster(time: float) -> float: return np.ceil(time / raster_time) * raster_time def _calc_area_and_times(grad_amp, flat_time, slew_rate, grad_start, grad_end): - time_ramp_up = _calc_ramp_up_time(grad_amp, slew_rate, grad_start) - time_ramp_down = _calc_ramp_down_time(grad_amp, slew_rate, grad_end) + time_ramp_up = _calc_ramp_time(grad_amp, grad_start, slew_rate) + time_ramp_down = _calc_ramp_time(grad_amp, grad_end, slew_rate) area = _calc_area(grad_amp, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end) return area @@ -57,17 +64,17 @@ def _calc_area(grad_amp, time_ramp_up, flat_time, time_ramp_down, grad_start, gr + 0.5 * (grad_amp + grad_end) * time_ramp_down ) - def _calc_ramp_up_time(grad_amp, slew_rate, grad_start): - return abs(grad_amp - grad_start) / slew_rate - - def _calc_ramp_down_time(grad_amp, slew_rate, grad_end): - return abs(grad_amp - grad_end) / slew_rate + def _calc_ramp_time(grad_1, grad_2, slew_rate): + return abs(grad_1 - grad_2) / slew_rate # first try to find grad_amp without a flat time. # for better convergence, we ignore the raster effect here and use a very small time step # we will take care of raster effect later on flat_time = 0 - obj1 = lambda x: (area - _calc_area_and_times(x, flat_time, slew_rate, grad_start, grad_end)) ** 2 + + def obj1(x): + return (area - _calc_area_and_times(x, flat_time, slew_rate, grad_start, grad_end)) ** 2 + # try different starting points for the optimization optimization_results = [ minimize(fun=obj1, x0=-max_grad, method="Nelder-Mead"), @@ -87,7 +94,10 @@ def _calc_ramp_down_time(grad_amp, slew_rate, grad_end): if obj1_value**0.5 > 1e-3 or abs(grad_amp) > max_grad: # Optimization without flat_time did not converge. Try to find flat_time. grad_amp = math.copysign(max_grad, grad_amp) - obj2 = lambda x: (area - _calc_area_and_times(grad_amp, x, slew_rate, grad_start, grad_end)) ** 2 + + def obj2(x): + return (area - _calc_area_and_times(grad_amp, x, slew_rate, grad_start, grad_end)) ** 2 + res2 = minimize(fun=obj2, x0=system.grad_raster_time, method="Nelder-Mead") flat_time = res2.x[0] # assert flat_time is larger than 0 and that the objective function is small @@ -95,17 +105,27 @@ def _calc_ramp_down_time(grad_amp, slew_rate, grad_end): # ensure times are on raster flat_time = _to_raster(flat_time) - time_ramp_up = _to_raster(_calc_ramp_up_time(grad_amp, slew_rate, grad_start)) - time_ramp_down = _to_raster(_calc_ramp_down_time(grad_amp, slew_rate, grad_end)) + time_ramp_up = _to_raster(_calc_ramp_time(grad_amp, grad_start, slew_rate)) + time_ramp_down = _to_raster(_calc_ramp_time(grad_amp, grad_end, slew_rate)) # recalculate grad_amp while taking raster effect into account - obj3 = lambda x: (area - _calc_area(x, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end)) ** 2 + def obj3(x): + return (area - _calc_area(x, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end)) ** 2 + res = minimize(fun=obj3, x0=grad_amp, method="Nelder-Mead") grad_amp = res.x[0] if flat_time > 0: - times = cumsum(0, time_ramp_up, flat_time, time_ramp_down) - amplitudes = [grad_start, grad_amp, grad_amp, grad_end] + if time_ramp_down > 0: + if time_ramp_up > 0: + times = cumsum(0, time_ramp_up, flat_time, time_ramp_down) + amplitudes = [grad_start, grad_amp, grad_amp, grad_end] + else: + times = cumsum(0, flat_time, time_ramp_down) + amplitudes = [grad_amp, grad_amp, grad_end] + else: + times = cumsum(0, time_ramp_up, flat_time) + amplitudes = [grad_start, grad_amp, grad_end] else: if time_ramp_down > 0: if time_ramp_up > 0: From 394c5fe32ca5712672440e0e377573b244aeebab Mon Sep 17 00:00:00 2001 From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com> Date: Tue, 9 Apr 2024 15:43:25 +0200 Subject: [PATCH 3/7] - Changed `make_extended_trapezoid_area` to use a grid search to search for solutions - Added tests for `make_extended_trapezoid_area` --- pypulseq/make_extended_trapezoid_area.py | 194 ++++++++++-------- .../test_make_extended_trapezoid_area.py | 137 +++++++++++++ 2 files changed, 242 insertions(+), 89 deletions(-) create mode 100644 pypulseq/tests/test_make_extended_trapezoid_area.py diff --git a/pypulseq/make_extended_trapezoid_area.py b/pypulseq/make_extended_trapezoid_area.py index 68eb6ef..421f317 100644 --- a/pypulseq/make_extended_trapezoid_area.py +++ b/pypulseq/make_extended_trapezoid_area.py @@ -1,9 +1,7 @@ -import math from types import SimpleNamespace from typing import Tuple, Union import numpy as np -from scipy.optimize import minimize from pypulseq.make_extended_trapezoid import make_extended_trapezoid from pypulseq.opts import Opts @@ -15,9 +13,9 @@ def make_extended_trapezoid_area( channel: str, grad_start: float, grad_end: float, - system: Union[Opts, None] = None, -) -> Tuple[SimpleNamespace, np.ndarray, np.ndarray]: - """Make (shortest possible) extended trapezoid for given area and gradient start and end point. + 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. Parameters ---------- @@ -43,100 +41,118 @@ def make_extended_trapezoid_area( """ if not system: system = Opts() - - slew_rate = system.max_slew * 0.99 + + max_slew = system.max_slew * 0.99 max_grad = system.max_grad * 0.99 raster_time = system.grad_raster_time def _to_raster(time: float) -> float: return np.ceil(time / raster_time) * raster_time - def _calc_area_and_times(grad_amp, flat_time, slew_rate, grad_start, grad_end): - time_ramp_up = _calc_ramp_time(grad_amp, grad_start, slew_rate) - time_ramp_down = _calc_ramp_time(grad_amp, grad_end, slew_rate) - area = _calc_area(grad_amp, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end) - return area - - def _calc_area(grad_amp, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end): - return ( - 0.5 * time_ramp_up * (grad_amp + grad_start) - + grad_amp * flat_time - + 0.5 * (grad_amp + grad_end) * time_ramp_down - ) - - def _calc_ramp_time(grad_1, grad_2, slew_rate): - return abs(grad_1 - grad_2) / slew_rate - - # first try to find grad_amp without a flat time. - # for better convergence, we ignore the raster effect here and use a very small time step - # we will take care of raster effect later on - flat_time = 0 - - def obj1(x): - return (area - _calc_area_and_times(x, flat_time, slew_rate, grad_start, grad_end)) ** 2 - - # try different starting points for the optimization - optimization_results = [ - minimize(fun=obj1, x0=-max_grad, method="Nelder-Mead"), - minimize(fun=obj1, x0=0, method="Nelder-Mead"), - minimize(fun=obj1, x0=max_grad, method="Nelder-Mead"), - ] - # choose the best result - optimization_results = np.array([(*res.x, res.fun) for res in optimization_results]) - grad_amp_results, obj1val_results = ( - optimization_results[:, 0], - optimization_results[:, 1], - ) - best_index = np.argmin(obj1val_results) - grad_amp = grad_amp_results[best_index] - obj1_value = obj1val_results[best_index] - - if obj1_value**0.5 > 1e-3 or abs(grad_amp) > max_grad: - # Optimization without flat_time did not converge. Try to find flat_time. - grad_amp = math.copysign(max_grad, grad_amp) - - def obj2(x): - return (area - _calc_area_and_times(grad_amp, x, slew_rate, grad_start, grad_end)) ** 2 - - res2 = minimize(fun=obj2, x0=system.grad_raster_time, method="Nelder-Mead") - flat_time = res2.x[0] - # assert flat_time is larger than 0 and that the objective function is small - assert flat_time >= 0, "flat time is smaller than 0, this should not happen" - - # ensure times are on raster - flat_time = _to_raster(flat_time) - time_ramp_up = _to_raster(_calc_ramp_time(grad_amp, grad_start, slew_rate)) - time_ramp_down = _to_raster(_calc_ramp_time(grad_amp, grad_end, slew_rate)) - - # recalculate grad_amp while taking raster effect into account - def obj3(x): - return (area - _calc_area(x, time_ramp_up, flat_time, time_ramp_down, grad_start, grad_end)) ** 2 + def _calc_ramp_time(grad_1, grad_2): + return _to_raster(abs(grad_1 - grad_2) / max_slew) + + # Find a solution for a given duration (specified in integer amounts of raster_time) + def _find_solution(duration): + # Grid search over all possible ramp-up and ramp-down times + # TODO: It should be possible to prune this search space significantly + # based on the reachable gradient amplitudes and the target area. + time_ramp_up = np.arange(1, min( + max(round(_calc_ramp_time(max_grad, grad_start) / raster_time), + round(_calc_ramp_time(-max_grad, grad_start) / raster_time)) + , duration-1) + +1) + time_ramp_down = np.arange(1, min( + max(round(_calc_ramp_time(max_grad, grad_end) / raster_time), + round(_calc_ramp_time(-max_grad, grad_end) / raster_time)) + , duration-1) + +1) + + time_ramp_up, time_ramp_down = np.meshgrid(time_ramp_up, time_ramp_down) + + # Calculate flat time and filter search space for valid timing + flat_time = duration - time_ramp_up - time_ramp_down + mask = flat_time >= 0 + + time_ramp_up = time_ramp_up[mask] + time_ramp_down = time_ramp_down[mask] + flat_time = flat_time[mask] + + # Calculate gradient strength for given timing + grad_amp = -(time_ramp_up*raster_time*grad_start + time_ramp_down*raster_time*grad_end - 2*area)/((time_ramp_up + 2*flat_time + time_ramp_down)*raster_time) + + # Filter solutions that satisfy max_grad and max_slew + slew_rate1 = abs(grad_start - grad_amp) / (time_ramp_up * raster_time) + slew_rate2 = abs(grad_end - grad_amp) / (time_ramp_down * raster_time) + mask = (abs(grad_amp) <= max_grad + 1e-8) & (slew_rate1 <= max_slew + 1e-8) & (slew_rate2 <= max_slew + 1e-8) + + solutions = np.flatnonzero(mask) + + # Check if there are solutions + if solutions.shape[0] == 0: + return None + + # Find solution with lowest slew rate and return it + ind = np.argmin(slew_rate1[mask] + slew_rate2[mask]) + ind = solutions[ind] + return (time_ramp_up[ind], flat_time[ind], time_ramp_down[ind], grad_amp[ind]) + + + # Perform a linear search + # This is necessary because there can exist a dead space where solutions + # do not exist for some durations longer than the optimal duration. The + # binary search below fails to find the optimum in those cases. + # TODO: Check that this range is sufficient, and maybe find out if the dead + # space can be calculated. + min_duration = max(round(_calc_ramp_time(grad_end, grad_start) / raster_time), 2) + max_duration = max(round(_calc_ramp_time(max_grad, grad_start) / raster_time), + round(_calc_ramp_time(-max_grad, grad_start) / raster_time), + round(_calc_ramp_time(max_grad, grad_end) / raster_time), + round(_calc_ramp_time(-max_grad, grad_end) / raster_time), + min_duration) + 1 + + solution = None + for duration in range(min_duration, max_duration): + solution = _find_solution(duration) + if solution: + break + + # Perform a binary search + if not solution: + max_duration = duration + + # First, find the upper limit on duration where a solution exists by + # exponentially expanding the duration. + while not solution: + max_duration *= 2 + solution = _find_solution(max_duration) + + def binary_search(fun, lower_limit, upper_limit): + if lower_limit == upper_limit - 1: + return fun(upper_limit) + + test_value = (upper_limit + lower_limit)//2 + + if fun(test_value): + return binary_search(fun, lower_limit, test_value) + else: + return binary_search(fun, test_value, upper_limit) + + solution = binary_search(_find_solution, max_duration//2, max_duration) - res = minimize(fun=obj3, x0=grad_amp, method="Nelder-Mead") - grad_amp = res.x[0] + # Get timing and gradient amplitude from solution + time_ramp_up = solution[0] * raster_time + flat_time = solution[1] * raster_time + time_ramp_down = solution[2] * raster_time + grad_amp = solution[3] + # Create extended trapezoid if flat_time > 0: - if time_ramp_down > 0: - if time_ramp_up > 0: - times = cumsum(0, time_ramp_up, flat_time, time_ramp_down) - amplitudes = [grad_start, grad_amp, grad_amp, grad_end] - else: - times = cumsum(0, flat_time, time_ramp_down) - amplitudes = [grad_amp, grad_amp, grad_end] - else: - times = cumsum(0, time_ramp_up, flat_time) - amplitudes = [grad_start, grad_amp, grad_end] + times = cumsum(0, time_ramp_up, flat_time, time_ramp_down) + amplitudes = np.array([grad_start, grad_amp, grad_amp, grad_end]) else: - if time_ramp_down > 0: - if time_ramp_up > 0: - times = cumsum(0, time_ramp_up, time_ramp_down) - amplitudes = np.array([grad_start, grad_amp, grad_end]) - else: - times = cumsum(0, time_ramp_down) - amplitudes = np.array([grad_start, grad_end]) - else: - times = cumsum(0, time_ramp_up) - amplitudes = np.array([grad_start, grad_end]) + 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) diff --git a/pypulseq/tests/test_make_extended_trapezoid_area.py b/pypulseq/tests/test_make_extended_trapezoid_area.py new file mode 100644 index 0000000..bb04fbe --- /dev/null +++ b/pypulseq/tests/test_make_extended_trapezoid_area.py @@ -0,0 +1,137 @@ +import pytest +import numpy as np +import random +import math + +from pypulseq import Opts, make_extended_trapezoid_area, make_extended_trapezoid, calc_duration +from pypulseq.utils.cumsum import cumsum + +system = Opts() + +test_zoo = [ + (0, 0, 1), + (0, 0, 10), + (0, 0, 100), + (0, 0, 10000), + (0, 1000, 100), + (-1000, 1000, 100), + (-1000, 0, 100), + (0, 0, -1), + (0, 0, -10), + (0, 0, -100), + (0, 0, -10000), + (0, 1000, -100), + (-1000, 1000, -100), + (-1000, 0, -100), + (0, system.max_grad*0.99, 10000), + (0, system.max_grad*0.99, -10000), + (0, -system.max_grad*0.99, 1000), + (0, -system.max_grad*0.99, -1000), + (system.max_grad*0.99, 0, 100), + (system.max_grad*0.99, 0, -100), + (-system.max_grad*0.99, 0, 1), + (-system.max_grad*0.99, 0, -1), + (0, 100000, 1), + (0, 100000, -1), + (0, -100000, 1), + (0, -100000, -1), + (0, 90000, 0.45), + (0, 90000, -0.45), + (0, -90000, 0.45), + (0, -90000, -0.45), + (0, 10000, 0.5 * (10000)**2 / (system.max_slew*0.99)), + (0, system.max_grad*0.99, 0.5 * (system.max_grad*0.99)**2 / (system.max_slew*0.99)), + (system.max_grad*0.99, system.max_grad*0.99, 1), + (system.max_grad*0.99, system.max_grad*0.99, -1), + ] + + +@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) + + grad_ok = all(abs(g.waveform) <= system.max_grad) + slew_ok = all(abs(np.diff(g.waveform) / np.diff(g.tt)) <= system.max_slew) + + assert pytest.approx(g.area) == area, 'Result area is not correct' + assert grad_ok, 'Maximum gradient strength violated' + assert slew_ok, 'Maximum slew rate violated' + + +random.seed(0) +test_zoo_random = [((random.random()-0.5) * 2 * system.max_grad * 0.99, + (random.random()-0.5) * 2 * system.max_grad * 0.99, + (random.random()-0.5) * 10000) for _ in range(100)] + +@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) + + grad_ok = all(abs(g.waveform) <= system.max_grad) + slew_ok = all(abs(np.diff(g.waveform) / np.diff(g.tt)) <= system.max_slew) + + assert pytest.approx(g.area) == area, 'Result area is not correct' + assert grad_ok, 'Maximum gradient strength violated' + assert slew_ok, 'Maximum slew rate violated' + + +random.seed(0) +test_zoo_random = [((random.random()-0.5) * 3 * system.max_grad, + (random.random()-0.5) * 3 * system.max_grad, + (random.random()-0.5) * 10 * system.max_grad) for _ in range(100)] + +@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 + + def _calc_ramp_time(grad_amp, slew_rate, grad_start): + return _to_raster(abs(grad_amp - grad_start) / slew_rate) + + grad_start = np.clip(grad_start, -system.max_grad, system.max_grad) + grad_end = np.clip(grad_end, -system.max_grad, system.max_grad) + + # Construct extended gradient based on start, intermediate, and end gradient + # strengths, assuming maximum slew rates. + + # If grad_amp > max_grad, convert to max_grad, and keep total area + # approximately equal + if abs(grad_amp) > system.max_grad*0.99: + grad_amp_new = math.copysign(system.max_grad*0.99, grad_amp) + + flat_time = _to_raster(_calc_ramp_time(grad_amp_new, system.max_slew*0.99, grad_amp) * abs(grad_amp_new + grad_amp) / abs(grad_amp_new)) + grad_amp = grad_amp_new + else: + flat_time = 0 + + # Construct extended gradient + if flat_time == 0: + times = cumsum(0, + _calc_ramp_time(grad_amp, system.max_slew*0.99, grad_start), + _calc_ramp_time(grad_amp, system.max_slew*0.99, grad_end), + ) + amplitudes = (grad_start, grad_amp, grad_end) + else: + times = cumsum(0, + _calc_ramp_time(grad_amp, system.max_slew*0.99, grad_start), + _to_raster(flat_time), + _calc_ramp_time(grad_amp, system.max_slew*0.99, grad_end), + ) + amplitudes = (grad_start, grad_amp, grad_amp, grad_end) + + g_true = make_extended_trapezoid(channel='x', amplitudes=amplitudes, times=times) + + # Recreate gradient using make_extended_trapezoid_area + g,_,_ = make_extended_trapezoid_area(channel='x', grad_start=grad_start, grad_end=grad_end, area=g_true.area, 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) + + assert pytest.approx(g.area) == g_true.area, 'Result area is not correct' + assert grad_ok, 'Maximum gradient strength violated' + assert slew_ok, 'Maximum slew rate violated' + + # Check that new gradient duration is smaller or equal to original gradient duration + d1 = calc_duration(g) + d2 = calc_duration(g_true) + + assert pytest.approx(d1) == d2 or d1 < d2 From c68e2ab1ab0be262a85b465b4891641f06e6f4c5 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Wed, 10 Apr 2024 17:58:58 +0200 Subject: [PATCH 4/7] improve readability and documentation --- pypulseq/make_extended_trapezoid_area.py | 147 ++++++++++++++--------- 1 file changed, 88 insertions(+), 59 deletions(-) diff --git a/pypulseq/make_extended_trapezoid_area.py b/pypulseq/make_extended_trapezoid_area.py index 421f317..1688248 100644 --- a/pypulseq/make_extended_trapezoid_area.py +++ b/pypulseq/make_extended_trapezoid_area.py @@ -13,7 +13,7 @@ def make_extended_trapezoid_area( channel: str, grad_start: float, grad_end: float, - system: Union[Opts, None] = None + 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. @@ -41,7 +41,7 @@ def make_extended_trapezoid_area( """ if not system: system = Opts() - + max_slew = system.max_slew * 0.99 max_grad = system.max_grad * 0.99 raster_time = system.grad_raster_time @@ -49,67 +49,95 @@ def make_extended_trapezoid_area( def _to_raster(time: float) -> float: return np.ceil(time / raster_time) * raster_time - def _calc_ramp_time(grad_1, grad_2): + def _calc_ramp_time(grad_1: float, grad_2: float) -> float: return _to_raster(abs(grad_1 - grad_2) / max_slew) - - # Find a solution for a given duration (specified in integer amounts of raster_time) - def _find_solution(duration): - # Grid search over all possible ramp-up and ramp-down times - # TODO: It should be possible to prune this search space significantly - # based on the reachable gradient amplitudes and the target area. - time_ramp_up = np.arange(1, min( - max(round(_calc_ramp_time(max_grad, grad_start) / raster_time), - round(_calc_ramp_time(-max_grad, grad_start) / raster_time)) - , duration-1) - +1) - time_ramp_down = np.arange(1, min( - max(round(_calc_ramp_time(max_grad, grad_end) / raster_time), - round(_calc_ramp_time(-max_grad, grad_end) / raster_time)) - , duration-1) - +1) - + + def _find_solution(duration: int) -> Union[None, Tuple[int, int, int, float]]: + """Find extended trapezoid gradient waveform for given duration. + + The function performs a grid search over all possible ramp-up, ramp-down and flat times + for the given duration and returns the solution with the lowest slew rate. + + Parameters + ---------- + duration + duration of the gradient in integer multiples of raster_time + + Returns + ------- + Tuple of ramp-up time, flat time, ramp-down time, gradient amplitude or None if no solution was found + + Raises + ------ + ValueError if no solution was found that satisfies the constraints and the desired area. + """ + + # Calculate possible ramp-up times for given duration + ramp_time_pos_start = round(_calc_ramp_time(max_grad, grad_start) / raster_time) + ramp_time_neg_start = round(_calc_ramp_time(-max_grad, grad_start) / raster_time) + max_ramp_time_start = max(ramp_time_pos_start, ramp_time_neg_start) + min_time_start = min(max_ramp_time_start, duration - 1) + time_ramp_up = np.arange(1, min_time_start + 1) + + # Calculate possible ramp-down times for given duration + ramp_time_pos_end = round(_calc_ramp_time(max_grad, grad_end) / raster_time) + ramp_time_neg_end = round(_calc_ramp_time(-max_grad, grad_end) / raster_time) + max_ramp_time_end = max(ramp_time_pos_end, ramp_time_neg_end) + min_time_end = min(max_ramp_time_end, duration - 1) + time_ramp_down = np.arange(1, min_time_end + 1) + + # Create meshgrid of possible ramp-up and ramp-down times time_ramp_up, time_ramp_down = np.meshgrid(time_ramp_up, time_ramp_down) - - # Calculate flat time and filter search space for valid timing + + # calculate corresponding flat times flat_time = duration - time_ramp_up - time_ramp_down - mask = flat_time >= 0 - - time_ramp_up = time_ramp_up[mask] - time_ramp_down = time_ramp_down[mask] - flat_time = flat_time[mask] - - # Calculate gradient strength for given timing - grad_amp = -(time_ramp_up*raster_time*grad_start + time_ramp_down*raster_time*grad_end - 2*area)/((time_ramp_up + 2*flat_time + time_ramp_down)*raster_time) - - # Filter solutions that satisfy max_grad and max_slew + + # filter search space for valid timings (flat time >= 0) + valid_indices = flat_time >= 0 + time_ramp_up = time_ramp_up[valid_indices] + time_ramp_down = time_ramp_down[valid_indices] + flat_time = flat_time[valid_indices] + + # Calculate gradient strength for given timing using analytical solution + grad_amp = -(time_ramp_up * raster_time * grad_start + time_ramp_down * raster_time * grad_end - 2 * area) / ( + (time_ramp_up + 2 * flat_time + time_ramp_down) * raster_time + ) + + # Calculate slew rates for given timings slew_rate1 = abs(grad_start - grad_amp) / (time_ramp_up * raster_time) slew_rate2 = abs(grad_end - grad_amp) / (time_ramp_down * raster_time) - mask = (abs(grad_amp) <= max_grad + 1e-8) & (slew_rate1 <= max_slew + 1e-8) & (slew_rate2 <= max_slew + 1e-8) - - solutions = np.flatnonzero(mask) - - # Check if there are solutions + + # Filter solutions that satisfy max_grad and max_slew constraints + valid_indices = ( + (abs(grad_amp) <= max_grad + 1e-8) & (slew_rate1 <= max_slew + 1e-8) & (slew_rate2 <= max_slew + 1e-8) + ) + solutions = np.flatnonzero(valid_indices) + + # Check if any valid solutions were found if solutions.shape[0] == 0: return None - + # Find solution with lowest slew rate and return it - ind = np.argmin(slew_rate1[mask] + slew_rate2[mask]) + ind = np.argmin(slew_rate1[valid_indices] + slew_rate2[valid_indices]) ind = solutions[ind] - return (time_ramp_up[ind], flat_time[ind], time_ramp_down[ind], grad_amp[ind]) - - + return (int(time_ramp_up[ind]), int(flat_time[ind]), int(time_ramp_down[ind]), float(grad_amp[ind])) + # Perform a linear search # This is necessary because there can exist a dead space where solutions # do not exist for some durations longer than the optimal duration. The # binary search below fails to find the optimum in those cases. - # TODO: Check that this range is sufficient, and maybe find out if the dead - # space can be calculated. + # TODO: Check if range is sufficient, try to calculate the dead space. min_duration = max(round(_calc_ramp_time(grad_end, grad_start) / raster_time), 2) - max_duration = max(round(_calc_ramp_time(max_grad, grad_start) / raster_time), - round(_calc_ramp_time(-max_grad, grad_start) / raster_time), - round(_calc_ramp_time(max_grad, grad_end) / raster_time), - round(_calc_ramp_time(-max_grad, grad_end) / raster_time), - min_duration) + 1 + max_duration = ( + max( + round(_calc_ramp_time(max_grad, grad_start) / raster_time), + round(_calc_ramp_time(-max_grad, grad_start) / raster_time), + round(_calc_ramp_time(max_grad, grad_end) / raster_time), + round(_calc_ramp_time(-max_grad, grad_end) / raster_time), + min_duration, + ) + + 1 + ) solution = None for duration in range(min_duration, max_duration): @@ -117,28 +145,28 @@ def _find_solution(duration): if solution: break - # Perform a binary search + # Perform a binary search if no solution was found if not solution: max_duration = duration - + # First, find the upper limit on duration where a solution exists by # exponentially expanding the duration. while not solution: max_duration *= 2 solution = _find_solution(max_duration) - - def binary_search(fun, lower_limit, upper_limit): + + def binary_search(fun, lower_limit, upper_limit): if lower_limit == upper_limit - 1: return fun(upper_limit) - - test_value = (upper_limit + lower_limit)//2 - + + test_value = (upper_limit + lower_limit) // 2 + if fun(test_value): return binary_search(fun, lower_limit, test_value) else: return binary_search(fun, test_value, upper_limit) - - solution = binary_search(_find_solution, max_duration//2, max_duration) + + solution = binary_search(_find_solution, max_duration // 2, max_duration) # Get timing and gradient amplitude from solution time_ramp_up = solution[0] * raster_time @@ -156,6 +184,7 @@ def binary_search(fun, lower_limit, upper_limit): grad = make_extended_trapezoid(channel=channel, system=system, times=times, amplitudes=amplitudes) - assert abs(grad.area - area) < 1e-3, "Area of the gradient is not equal to the desired area. Optimization failed." + 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 From 682864a48a9d0a82f15eae457563acdf1ff47a40 Mon Sep 17 00:00:00 2001 From: Patrick Schuenke Date: Wed, 10 Apr 2024 18:03:48 +0200 Subject: [PATCH 5/7] fix docstrings --- pypulseq/make_extended_trapezoid_area.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/pypulseq/make_extended_trapezoid_area.py b/pypulseq/make_extended_trapezoid_area.py index 1688248..88eb744 100644 --- a/pypulseq/make_extended_trapezoid_area.py +++ b/pypulseq/make_extended_trapezoid_area.py @@ -38,6 +38,10 @@ def make_extended_trapezoid_area( Time points of the extended trapezoid. amplitude : numpy.ndarray Amplitude values of the extended trapezoid. + + Raises + ------ + ValueError if no solution was found that satisfies the constraints and the desired area. """ if not system: system = Opts() @@ -66,10 +70,6 @@ def _find_solution(duration: int) -> Union[None, Tuple[int, int, int, float]]: Returns ------- Tuple of ramp-up time, flat time, ramp-down time, gradient amplitude or None if no solution was found - - Raises - ------ - ValueError if no solution was found that satisfies the constraints and the desired area. """ # Calculate possible ramp-up times for given duration From 3e3628f480de0ea5cfb17b361a434a7f0da1de18 Mon Sep 17 00:00:00 2001 From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com> Date: Thu, 11 Apr 2024 16:44:14 +0200 Subject: [PATCH 6/7] - Changed linear search limits in `make_extended_trapezoid_area` - Minor cleanups --- pypulseq/make_extended_trapezoid_area.py | 38 +++++++++++------------- 1 file changed, 18 insertions(+), 20 deletions(-) diff --git a/pypulseq/make_extended_trapezoid_area.py b/pypulseq/make_extended_trapezoid_area.py index 88eb744..24f31e9 100644 --- a/pypulseq/make_extended_trapezoid_area.py +++ b/pypulseq/make_extended_trapezoid_area.py @@ -76,23 +76,23 @@ def _find_solution(duration: int) -> Union[None, Tuple[int, int, int, float]]: ramp_time_pos_start = round(_calc_ramp_time(max_grad, grad_start) / raster_time) ramp_time_neg_start = round(_calc_ramp_time(-max_grad, grad_start) / raster_time) max_ramp_time_start = max(ramp_time_pos_start, ramp_time_neg_start) - min_time_start = min(max_ramp_time_start, duration - 1) - time_ramp_up = np.arange(1, min_time_start + 1) + max_time_start = min(max_ramp_time_start, duration - 1) + time_ramp_up = np.arange(1, max_time_start + 1) # Calculate possible ramp-down times for given duration ramp_time_pos_end = round(_calc_ramp_time(max_grad, grad_end) / raster_time) ramp_time_neg_end = round(_calc_ramp_time(-max_grad, grad_end) / raster_time) max_ramp_time_end = max(ramp_time_pos_end, ramp_time_neg_end) - min_time_end = min(max_ramp_time_end, duration - 1) - time_ramp_down = np.arange(1, min_time_end + 1) + max_time_end = min(max_ramp_time_end, duration - 1) + time_ramp_down = np.arange(1, max_time_end + 1) # Create meshgrid of possible ramp-up and ramp-down times time_ramp_up, time_ramp_down = np.meshgrid(time_ramp_up, time_ramp_down) - # calculate corresponding flat times + # Calculate corresponding flat times flat_time = duration - time_ramp_up - time_ramp_down - # filter search space for valid timings (flat time >= 0) + # Filter search space for valid timings (flat time >= 0) valid_indices = flat_time >= 0 time_ramp_up = time_ramp_up[valid_indices] time_ramp_down = time_ramp_down[valid_indices] @@ -128,27 +128,25 @@ def _find_solution(duration: int) -> Union[None, Tuple[int, int, int, float]]: # binary search below fails to find the optimum in those cases. # TODO: Check if range is sufficient, try to calculate the dead space. min_duration = max(round(_calc_ramp_time(grad_end, grad_start) / raster_time), 2) - max_duration = ( - max( - round(_calc_ramp_time(max_grad, grad_start) / raster_time), - round(_calc_ramp_time(-max_grad, grad_start) / raster_time), - round(_calc_ramp_time(max_grad, grad_end) / raster_time), - round(_calc_ramp_time(-max_grad, grad_end) / raster_time), - min_duration, - ) - + 1 - ) + # Calculate duration needed to ramp down gradient to zero. + # From this point onwards, solutions can always be found by extending + # the duration and doing a binary search. + max_duration = max( + round(_calc_ramp_time(0, grad_start) / raster_time), + round(_calc_ramp_time(0, grad_end) / raster_time), + min_duration, + ) + + # Linear search solution = None - for duration in range(min_duration, max_duration): + for duration in range(min_duration, max_duration + 1): solution = _find_solution(duration) if solution: break - # Perform a binary search if no solution was found + # Perform a binary search for duration > max_duration if no solution was found if not solution: - max_duration = duration - # First, find the upper limit on duration where a solution exists by # exponentially expanding the duration. while not solution: From 92b5f5876d748e068a6a3b332590feeed04d826e Mon Sep 17 00:00:00 2001 From: Frank Zijlstra <22915457+FrankZijlstra@users.noreply.github.com> Date: Fri, 12 Apr 2024 14:01:20 +0200 Subject: [PATCH 7/7] Modified `_find_solution` in `make_extended_trapezoid_area` to only try reasonable timings (primarily those using maximum slew rate) --- pypulseq/make_extended_trapezoid_area.py | 74 +++++++++++++++++++----- 1 file changed, 58 insertions(+), 16 deletions(-) diff --git a/pypulseq/make_extended_trapezoid_area.py b/pypulseq/make_extended_trapezoid_area.py index 24f31e9..cd42676 100644 --- a/pypulseq/make_extended_trapezoid_area.py +++ b/pypulseq/make_extended_trapezoid_area.py @@ -72,22 +72,64 @@ def _find_solution(duration: int) -> Union[None, Tuple[int, int, int, float]]: Tuple of ramp-up time, flat time, ramp-down time, gradient amplitude or None if no solution was found """ - # Calculate possible ramp-up times for given duration - ramp_time_pos_start = round(_calc_ramp_time(max_grad, grad_start) / raster_time) - ramp_time_neg_start = round(_calc_ramp_time(-max_grad, grad_start) / raster_time) - max_ramp_time_start = max(ramp_time_pos_start, ramp_time_neg_start) - max_time_start = min(max_ramp_time_start, duration - 1) - time_ramp_up = np.arange(1, max_time_start + 1) - - # Calculate possible ramp-down times for given duration - ramp_time_pos_end = round(_calc_ramp_time(max_grad, grad_end) / raster_time) - ramp_time_neg_end = round(_calc_ramp_time(-max_grad, grad_end) / raster_time) - max_ramp_time_end = max(ramp_time_pos_end, ramp_time_neg_end) - max_time_end = min(max_ramp_time_end, duration - 1) - time_ramp_down = np.arange(1, max_time_end + 1) - - # Create meshgrid of possible ramp-up and ramp-down times - time_ramp_up, time_ramp_down = np.meshgrid(time_ramp_up, time_ramp_down) + # Determine timings to check for possible solutions + ramp_up_times = [] + ramp_down_times = [] + + # First, consider solutions that use maximum slew rate: + # Analytically calculate calculate the point where: + # grad_start + ramp_up_time * max_slew == grad_end + ramp_down_time * max_slew + ramp_up_time = (duration * max_slew * raster_time - grad_start + grad_end) / ( + 2*max_slew*raster_time + ) + ramp_up_time = round(ramp_up_time) + + # 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: + 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: + ramp_down_time = duration - ramp_up_time + + # Add possible solution if timing is valid + if (ramp_up_time > 0 + and ramp_down_time > 0 + and ramp_up_time + ramp_down_time <= duration): + ramp_up_times.append(ramp_up_time) + ramp_down_times.append(ramp_down_time) + + # Analytically calculate calculate the point where: + # grad_start - ramp_up_time * max_slew == grad_end - ramp_down_time * max_slew + ramp_up_time = (duration * max_slew * raster_time + grad_start - grad_end) / ( + 2*max_slew*raster_time + ) + ramp_up_time = round(ramp_up_time) + + # 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: + 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: + ramp_down_time = duration - ramp_up_time + + # Add possible solution if timing is valid + if (ramp_up_time > 0 + and ramp_down_time > 0 + and ramp_up_time + ramp_down_time <= duration): + ramp_up_times.append(ramp_up_time) + ramp_down_times.append(ramp_down_time) + + # Second, try any solution with flat_time == 0 + # This appears to be necessary for many cases, but going through all + # timings here is probably too conservative still. + for ramp_up_time in range(1, duration): + ramp_up_times.append(ramp_up_time) + ramp_down_times.append(duration - ramp_up_time) + + time_ramp_up = np.array(ramp_up_times) + time_ramp_down = np.array(ramp_down_times) # Calculate corresponding flat times flat_time = duration - time_ramp_up - time_ramp_down