Fix make_extended_trapezoid area

[schuenke]

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.

The issue accured for example when trying to use make_extended_trapezoid_area to calculate the rewinder gradients for spiral MRI sequences.

[FrankZijlstra]

Looks good in principle, works better than the original function in all cases I've tested.

In the code as-is, there are edge-cases in the form of area == grad_end * raster_time * N (assuming grad_start == 0). That is, where the requested area can be achieved with a triangle gradient between start and end.

In some cases they give an assertion error on (very small) negative flat_time:
pp.make_extended_trapezoid_area(1, 'x', 100000, 0, system)

In other cases they result in a non-time-optimal solution, due to the rounding to raster times:
pp.make_extended_trapezoid_area(0.45, 'x', 90000, 0, system)
Where the optimal solution is given by:
pp.make_extended_trapezoid('x', times=[0, 1e-5], amplitudes=[0,90000])

Another set of edge cases is in the form of area == 0.5 * grad_end**2 / (system.max_slew * 0.99) (assuming grad_start == 0), where the optimization fails. This is the case where the requested area is equal to the area of ramping up to grad_end with 0.99 times the maximum slew rate (as used in the function).
pp.make_extended_trapezoid_area(0.5 * 100000**2 / (system.max_slew*0.99), 'x', 100000, 0, system)

An even more rare set of edge cases exist where the taking ceil of the ramp up/down and flat times is not optimal. I think the worst case here is that the solution ends up being one or two raster times off the optimum, which is not a big disaster. But the code could in principle consider all 8 solutions (rounding up or down for each of the times).

Besides these edge-cases, an issue I see is that with these changes that remove the discontinuity from the optimization problems, all optimization problems now have closed-form solutions, making all the iterative optimization unnecessary, solutions below:

The first set of optimization:

    r = slew_rate
    s = grad_start
    z = grad_end
    v = area
       
    if (s**2 + 2*r*v + z**2) < 0:
        x1 = np.sqrt(-(s**2 + 2*r*v + z**2)) / np.sqrt(2)
    else:
        x1 = np.sqrt(s**2 + 2*r*v + z**2) / np.sqrt(2)
    
    if (s**2 - 2*r*v + z**2) < 0:
        x2 = np.sqrt(-(s**2 - 2*r*v + z**2)) / np.sqrt(2)
    else:
        x2 = np.sqrt(s**2 - 2*r*v + z**2) / np.sqrt(2)
    
    grad_amp_results = [x1, -x1, x2, -x2]
    obj1val_results = [obj1(x) for x in grad_amp_results]
    print(grad_amp_results, obj1val_results)

Based on https://www.wolframalpha.com/input?i=%28s-x%29%2Fr*+%28s%2Bx%29%2F2+%2B+%28z-x%29%2Fr*+%28x%2Bz%29%2F2+%3D+v+solve+for+x and https://www.wolframalpha.com/input?i=%28x-s%29%2Fr*+%28s%2Bx%29%2F2+%2B+%28x-z%29%2Fr*+%28x%2Bz%29%2F2+%3D+v+solve+for+x
I think here it should be possible to reduce the solutions to two, or maybe even one, solution, based on whether the requested area is larger or smaller than the area of a triangle gradient between grad_start and grad_end. E.g. if the area is larger, then x must be larger than max(grad_start,grad_end).

The second optimization:

        x = grad_amp
        flat_time = -(-2*r*v + s**2 - 2*x**2 + z**2)/(2*r*x)
        print(flat_time)

Based on https://www.wolframalpha.com/input?i=%28s-x%29%2Fr*+%28s%2Bx%29%2F2+%2B+f*x+%2B+%28z-x%29%2Fr*+%28x%2Bz%29%2F2+%3D+v+solve+for+f

And the final optimization:

    a = time_ramp_up
    f = flat_time
    b = time_ramp_down
    grad_amp = -(a*s + b*z - 2*v)/(a+b+2*f)

Based on https://www.wolframalpha.com/input?i=solve+a*%28s%2Bx%29%2F2+%2B+f*x+%2B+b*%28x%2Bz%29%2F2+%3D+v+for+x

These can be easily replaced into the code, obviously with some renaming of the variables I used to make wolfram alpha happy :). These may also solve some of the edge-cases, or make them easier to detect.

[schuenke]

Looks good in principle, works better than the original function in all cases I've tested.

In the code as-is, there are edge-cases in the form of area == grad_end * raster_time * N (assuming grad_start == 0). That is, where the requested area can be achieved with a triangle gradient between start and end.

In some cases they give an assertion error on (very small) negative flat_time: pp.make_extended_trapezoid_area(1, 'x', 100000, 0, system)

In other cases they result in a non-time-optimal solution, due to the rounding to raster times: pp.make_extended_trapezoid_area(0.45, 'x', 90000, 0, system) Where the optimal solution is given by: pp.make_extended_trapezoid('x', times=[0, 1e-5], amplitudes=[0,90000])

Another set of edge cases is in the form of area == 0.5 * grad_end**2 / (system.max_slew * 0.99) (assuming grad_start == 0), where the optimization fails. This is the case where the requested area is equal to the area of ramping up to grad_end with 0.99 times the maximum slew rate (as used in the function). pp.make_extended_trapezoid_area(0.5 * 100000**2 / (system.max_slew*0.99), 'x', 100000, 0, system)

Yeah, the simple ramp-up or ramp-down cases were and are not handled separately atm. We only differentiate between with or without flat time. Probably this is something we should do directly in the beginning. It should fix all cases mentioned above, no ?!

An even more rare set of edge cases exist where the taking ceil of the ramp up/down and flat times is not optimal. I think the worst case here is that the solution ends up being one or two raster times off the optimum, which is not a big disaster. But the code could in principle consider all 8 solutions (rounding up or down for each of the times).

Rounding down the ramp-times should not be possible, should it? Maybe a recalculation of the flat_time due to larger ramping areas makes sense, but I don't think this will ever lead to deviations of more than one raster_time from the optimum.

Besides these edge-cases, an issue I see is that with these changes that remove the discontinuity from the optimization problems, all optimization problems now have closed-form solutions, making all the iterative optimization unnecessary, solutions below:

Yeah, we were thinking about the closed-form solution as well, but catching all edge cases, especially when enforcing the raster, is also not super simple and clean. Thus, we decided to provide a simple fix for the current version and maybe implement an alternative version later on. Would be nice to consider nulling the 1st moment as well as for example in https://onlinelibrary.wiley.com/doi/10.1002/mrm.20489

The first set of optimization:

    r = slew_rate
    s = grad_start
    z = grad_end
    v = area
       
    if (s**2 + 2*r*v + z**2) < 0:
        x1 = np.sqrt(-(s**2 + 2*r*v + z**2)) / np.sqrt(2)
    else:
        x1 = np.sqrt(s**2 + 2*r*v + z**2) / np.sqrt(2)
    
    if (s**2 - 2*r*v + z**2) < 0:
        x2 = np.sqrt(-(s**2 - 2*r*v + z**2)) / np.sqrt(2)
    else:
        x2 = np.sqrt(s**2 - 2*r*v + z**2) / np.sqrt(2)
    
    grad_amp_results = [x1, -x1, x2, -x2]
    obj1val_results = [obj1(x) for x in grad_amp_results]
    print(grad_amp_results, obj1val_results)

Based on https://www.wolframalpha.com/input?i=%28s-x%29%2Fr*+%28s%2Bx%29%2F2+%2B+%28z-x%29%2Fr*+%28x%2Bz%29%2F2+%3D+v+solve+for+x and https://www.wolframalpha.com/input?i=%28x-s%29%2Fr*+%28s%2Bx%29%2F2+%2B+%28x-z%29%2Fr*+%28x%2Bz%29%2F2+%3D+v+solve+for+x I think here it should be possible to reduce the solutions to two, or maybe even one, solution, based on whether the requested area is larger or smaller than the area of a triangle gradient between grad_start and grad_end. E.g. if the area is larger, then x must be larger than max(grad_start,grad_end).

The second optimization:

        x = grad_amp
        flat_time = -(-2*r*v + s**2 - 2*x**2 + z**2)/(2*r*x)
        print(flat_time)

Based on https://www.wolframalpha.com/input?i=%28s-x%29%2Fr*+%28s%2Bx%29%2F2+%2B+f*x+%2B+%28z-x%29%2Fr*+%28x%2Bz%29%2F2+%3D+v+solve+for+f

And the final optimization:

    a = time_ramp_up
    f = flat_time
    b = time_ramp_down
    grad_amp = -(a*s + b*z - 2*v)/(a+b+2*f)

Based on https://www.wolframalpha.com/input?i=solve+a*%28s%2Bx%29%2F2+%2B+f*x+%2B+b*%28x%2Bz%29%2F2+%3D+v+for+x

These can be easily replaced into the code, obviously with some renaming of the variables I used to make wolfram alpha happy :). These may also solve some of the edge-cases, or make them easier to detect.

For the moment, I would prefer to stick with the old style using iterative optimization and just change the other points you mentioned. I will adjust the code and push the changes soon.

[fzimmermann89]

An even more rare set of edge cases exist where the taking ceil of the ramp up/down and flat times is not optimal. I think the worst case here is that the solution ends up being one or two raster times off the optimum, which is not a big disaster. But the code could in principle consider all 8 solutions (rounding up or down for each of the times).

The calculated time point is the earliest possible time to achieve the area while satisfying the slope constraint.
Rounding down would always slightly exceed the maximum slope afaik.
Thus, we have to always round up the end times, right?

[FrankZijlstra]

Yeah, the simple ramp-up or ramp-down cases were and are not handled separately atm. We only differentiate between with or without flat time. Probably this is something we should do directly in the beginning. It should fix all cases mentioned above, no ?!
All, except the final case, where the triangle is achievable with unrounded ramp times, but not with rounded ramp times. In the first closed form solution, this is the case where one of the discriminants is exactly zero. That is, any solution grad_start <= x <= grad_end would yield an objective of exactly 0.

Rounding down the ramp-times should not be possible, should it? Maybe a recalculation of the flat_time due to larger ramping areas makes sense, but I don't think this will ever lead to deviations of more than one raster_time from the optimum.
It's not a big issue, I agree that it will be on the order of one or two raster times off the optimum, but it is possible. If you round up any of the times, the gradient area increases, and grad_amp will need to come down as a result. This lower amplitude can mean a ramp time that was just above the raster time could come down to just below.

Yeah, we were thinking about the closed-form solution as well, but catching all edge cases, especially when enforcing the raster, is also not super simple and clean. Thus, we decided to provide a simple fix for the current version and maybe implement an alternative version later on. Would be nice to consider nulling the 1st moment as well as for example in https://onlinelibrary.wiley.com/doi/10.1002/mrm.20489

I don't think the solutions I provided introduce any new edge cases. In fact, I think it makes it easier to detect the edge cases (i.e. avoid divisions by 0).

For the moment, I would prefer to stick with the old style using iterative optimization and just change the other points you mentioned. I will adjust the code and push the changes soon.

All that needs to happen on the analytical solutions is a little cleanup as I just hacked them together :). But they do work fine!

[FrankZijlstra]

Little update on the analytical solutions:

First optimization, simplified and added a check for the case where the gradient is triangular (still needs to be preceded with a check that the solution is not just a triangular gradient with a ramp time that fits raster_time):

    if area == (s+z)/2 * abs(s-z)/r:
        # Solution is triangular gradient, set grad_amp to violate max_grad, so flat_time optimization follows
        if area > (s+z)/2 * _to_raster(abs(s-z)/r):
            grad_amp = -max_grad*2
        else:
            grad_amp = max_grad*2
    elif area > (s+z)/2 * abs(s-z)/r:
        grad_amp = np.sqrt(s**2 + 2*r*v + z**2) / np.sqrt(2)
    else:
        grad_amp = -np.sqrt(s**2 - 2*r*v + z**2) / np.sqrt(2)

Flat time optimization needed a second case:

        x = grad_amp
        if x < 0:
            flat_time = -(-2*r*v + s**2 - 2*x**2 + z**2)/(2*r*x)
        else:
            flat_time = (2*r*v + s**2 - 2*x**2 + z**2)/(2*r*x)

In addition, I wrote a little pytest script for the function, with some of the things I tried, I suggest you add it to pypulseq/tests, and add test cases for things I have not covered, or if new things come up:
test_make_extended_trapezoid_area.zip (see comment below)

[FrankZijlstra]

Bigger update... Because I was curious about the impact of rounding up or down, I added a test function based on recreating an existing extended gradient. It finds many failure cases, some 1e-5 off on the timing, but some a few that are a multiple of the duration off. I'm not sure what exactly these cases represent, but that the come up in a test with randomly generated values is a bit worrying.

test_make_extended_trapezoid_area.zip

[btasdelen]

@schuenke Just as a heads-up, I don't know how it will perform with the recent changes to make_extended_trapezoid_area, but when I tried using this approach for Spiral rewinders, accumulated errors in k-space were large enough to see some artifacts (for bSSFP), which I believe comes from rounding up to the raster time. My guess is there might be another solution that has less error when round-up, but we won't know as the optimization assumes continuous case.

Another issue was the rewinders were significantly longer than they need to be. At the end of the day, I started using GrOpt, and I am quite happy with the result.

[fzimmermann89]

The goal of the rewrite was also to reduce this error..

Once the different time points are determined, the last step should find the optimal gradient amplitude given the (rounded to raster time) time points to achieve exactly the requested amplitude.

It should also always find the shortest possible gradient up to (worst case) 3 gradient raster steps.

[btasdelen]

@fzimmermann89 Correct me if I'm wrong, but it is the shortest trapezoid, not the shortest gradient, that nulls the M0.

[fzimmermann89]

@fzimmermann89 Correct me if I'm wrong, but it is the shortest trapezoid, not the shortest gradient, that nulls the M0.

Maybe I am wrong, but I'd say it should be the shortest gradient that has the requested area.
Ignoring raster_time for now, the shortest possible should rise with the maximum allowed slew rate. If there was a short gradient that has somewhere a slope smaller than the max_slew rate, you can always make it shorter by rising faster. So the fastest gradient to null the 0-th moment will be the trapezoid. See also, for example, https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7606655/figure/F2/.

The only difference can be of the order of few (~3) gradient_raster steps ..

Or where is my error? How would a shorter gradient look like that has the same area as a max-slew trapezoid?

[btasdelen]

@fzimmermann89 You might be right, and it might be a problem with my implementation back then that I was not able to find the shortest gradient. I can not think of any case that there is a shorter gradient.

[schuenke]

@schuenke Just as a heads-up, I don't know how it will perform with the recent changes to make_extended_trapezoid_area, but when I tried using this approach for Spiral rewinders, accumulated errors in k-space were large enough to see some artifacts (for bSSFP), which I believe comes from rounding up to the raster time. My guess is there might be another solution that has less error when round-up, but we won't know as the optimization assumes continuous case.

Another issue was the rewinders were significantly longer than they need to be. At the end of the day, I started using GrOpt, and I am quite happy with the result.

Thanks for the feedback. I'll keep it in mind, but I would definitely prefer an elegant and efficient solution in pypulseq. At least for the "simple" 0-th moment nulling. Pretty optimistic we can get this done 💪

[schuenke]

    if area == (s+z)/2 * abs(s-z)/r:
        # Solution is triangular gradient, set grad_amp to violate max_grad, so flat_time optimization follows
        if area > (s+z)/2 * _to_raster(abs(s-z)/r):
            grad_amp = -max_grad*2
        else:
            grad_amp = max_grad*2
    elif area > (s+z)/2 * abs(s-z)/r:
        grad_amp = np.sqrt(s**2 + 2*r*v + z**2) / np.sqrt(2)
    else:
        grad_amp = -np.sqrt(s**2 - 2*r*v + z**2) / np.sqrt(2)

Hmm, not sure about this.
Lets define ramp_time = abs(s-z)/r as the ramp time, which is NOT on raster

_to_raster(ramp_time) will always be larger than ramp_time, because we use np.ceil()
This means area > (s+z)/2 * _to_raster(abs(s-z)/r) should never be reachable, no ?!

IMO we only want to use a single triangular gradient (just ONE ramp) if area == (s+z)/2 * _to_raster(abs(s-z)/r), which is a very rare case.

Just to make clear we all have the same things in mind... we expect the following cases:

1. a grad with a single ramp from grad_start to grad_end, which could potentially be of length 1 * raster_time
2. a grad with 2 ramps, but no flat time. This means one ramp from grad_start to grad_amp and another from grad_amp to grad_end, with arbitrary signs of grad_start/grad_amp/grad_end. Note: this might be the same as case 1 (for length >= 2 * raster_time), if the slopes of both ramps are identical.
3. a grad with 2 slopes (ramp_up and ramp_down) and a flat time
4. a grad with ramp_up and a flat_time, but no ramp_down
5. a grad with a flat_time and ramp_down, but no ramp_up

[FrankZijlstra]

Hmm, not sure about this. Lets define ramp_time = abs(s-z)/r as the ramp time, which is NOT on raster

_to_raster(ramp_time) will always be larger than ramp_time, because we use np.ceil() This means area > (s+z)/2 * _to_raster(abs(s-z)/r) should never be reachable, no ?!

IMO we only want to use a single triangular gradient (just ONE ramp) if area == (s+z)/2 * _to_raster(abs(s-z)/r), which is a very rare case.

s+z can be negative. In plain text it might be easier to understand:
If the requested area is equal to the area of a triangle between start and end without rounding, then (assuming we already checked that it wasn't the one special case that you mention) we need to check if the desired area is larger or smaller than the area of a triangle between start and end with rounding, which will inform whether grad_amp should be +grad_max or -grad_max. And then we solve for flat_time.

But I just realize there is still a problem with that case, because that does not necessarily mean flat_time > 0. It will get fixed later by the third optimization, but the ramp times will be too long.

Just to be clear, this is not a problem with just the analytical solution, the same thing happens when using optimization.

Just to make clear we all have the same things in mind... we expect the following cases:

1. a grad with a single ramp from grad_start to grad_end, which could potentially be of length 1 * raster_time

2. a grad with 2 ramps, but no flat time. This means one ramp from grad_start to grad_amp and another from grad_amp to grad_end, with arbitrary signs of grad_start/grad_amp/grad_end. Note: this might be the same as case 1 (for length >= 2 * raster_time), if the slopes of both ramps are identical.

3. a grad with 2 slopes (ramp_up and ramp_down) and a flat time

4. a grad with ramp_up and a flat_time, but no ramp_down

5. a grad with a flat_time and ramp_down, but no ramp_up

On case 1, I assume you mean (integer) N * raster_time. And if you check for case 1, then a situation with equal slopes in case 2 will never happen. And it shouldn't, because the optimization problem is ill-defined in that case, every solution grad_start <= grad_amp <= grad_end will be identical. And with that case excluded, the optimization problem is much narrower with only two possibilities, because we already know whether grad_amp should be larger or smaller than max(grad_start, grad_end).

And in case 4 and 5, I noticed issues with start or end set to max_grad*0.99: the initial grad_amp solution can be max_grad*0.99, which would mean there is no ramp up or ramp down. But the flat_time is not rounded yet, which would change the area, and change grad_amp, and therefore require both ramp up and ramp down.

To be clear, I'm not expecting every edge case to be solved per se, but at the very least we can mark where there are issues. Right now I'm most worried about my second set of tests (recreating existing gradients), giving more errors than with the original code.

[schuenke]

Hmm, not sure about this. Lets define ramp_time = abs(s-z)/r as the ramp time, which is NOT on raster
_to_raster(ramp_time) will always be larger than ramp_time, because we use np.ceil() This means area > (s+z)/2 * _to_raster(abs(s-z)/r) should never be reachable, no ?!
IMO we only want to use a single triangular gradient (just ONE ramp) if area == (s+z)/2 * _to_raster(abs(s-z)/r), which is a very rare case.

s+z can be negative. In plain text it might be easier to understand: If the requested area is equal to the area of a triangle between start and end without rounding, then (assuming we already checked that it wasn't the one special case that you mention) we need to check if the desired area is larger or smaller than the area of a triangle between start and end with rounding, which will inform whether grad_amp should be +grad_max or -grad_max. And then we solve for flat_time.

But I just realize there is still a problem with that case, because that does not necessarily mean flat_time > 0. It will get fixed later by the third optimization, but the ramp times will be too long.

Just to be clear, this is not a problem with just the analytical solution, the same thing happens when using optimization.

Just to make clear we all have the same things in mind... we expect the following cases:

1. a grad with a single ramp from grad_start to grad_end, which could potentially be of length 1 * raster_time

2. a grad with 2 ramps, but no flat time. This means one ramp from grad_start to grad_amp and another from grad_amp to grad_end, with arbitrary signs of grad_start/grad_amp/grad_end. Note: this might be the same as case 1 (for length >= 2 * raster_time), if the slopes of both ramps are identical.

3. a grad with 2 slopes (ramp_up and ramp_down) and a flat time

4. a grad with ramp_up and a flat_time, but no ramp_down

5. a grad with a flat_time and ramp_down, but no ramp_up

On case 1, I assume you mean (integer) N * raster_time. And if you check for case 1, then a situation with equal slopes in case 2 will never happen. And it shouldn't, because the optimization problem is ill-defined in that case, every solution grad_start <= grad_amp <= grad_end will be identical. And with that case excluded, the optimization problem is much narrower with only two possibilities, because we already know whether grad_amp should be larger or smaller than max(grad_start, grad_end).

And in case 4 and 5, I noticed issues with start or end set to max_grad*0.99: the initial grad_amp solution can be max_grad*0.99, which would mean there is no ramp up or ramp down. But the flat_time is not rounded yet, which would change the area, and change grad_amp, and therefore require both ramp up and ramp down.

To be clear, I'm not expecting every edge case to be solved per se, but at the very least we can mark where there are issues. Right now I'm most worried about my second set of tests (recreating existing gradients), giving more errors than with the original code.

Ah, sorry I somehow thought there was an abs() around the s+z... makes sense.

Just to clarify: I really like the idea of having a closed-form solution and I am motivated to implement it. I just want to think of potential edge-cases before the implementation.

Regarding the errors in the test: you also test for existing gradient without max_slew, don't you? There might be cases where using max_slew is not optimal. E.g. if area is slightly larger than the area of the simple single ramp, but using a slightly lower slew rate might do the job?! Or Is there always a solution with max_slew that should result in the same total time?

I will think about this again and try to define some example cases 😄

[FrankZijlstra]

Just to clarify: I really like the idea of having a closed-form solution and I am motivated to implement it. I just want to think of potential edge-cases before the implementation.

Regarding the errors in the test: you also test for existing gradient without max_slew, don't you? There might be cases where using max_slew is not optimal. E.g. if area is slightly larger than the area of the simple single ramp, but using a slightly lower slew rate might do the job?! Or Is there always a solution with max_slew that should result in the same total time?

I will think about this again and try to define some example cases 😄

I have also found myself strangely intrigued by this seemingly simple geometric problem ;). I am currently looking into the cases where the total duration is multiple times off the optimum.
E.g.
g,_,_ = pp.make_extended_trapezoid_area(channel='x', grad_start=-1680000, grad_end=-1000000, area=-100, system=system)
Which has a duration of 69e-5, while I know the optimum should be around 10e-5...

And no, I have not explicitly considered that not using max_slew could be more optimal. If this is the case, then not rounding the ramp times in the initial optimizations may create additional issues.

[FrankZijlstra]

So, it appears the issue at hand is that the first optimization problem can have at least 3 solutions which achieve a given area:

(y = area, x = grad_amp)

These will have near-zero, but not identical objective values, so it's essentially random which one gets picked. But they do not have the same duration. In a bad case, the longest duration gets picked, which is far off the optimum.

In addition, the solver does not always find these solutions. I think the tolerances on the optimization should be lower (at the cost of time), and the 0 initialization is always a stationary point, so that needs to be split into small positive and small negative initializations.

This is actually a good reason to use the closed-form solutions. But I think we need to consider all 4 possible solutions ([x1, -x1, x2, -x2] like in my first comment), and use the one with the shortest total time.

[FrankZijlstra]

In the closed-form solution of the third problem, this fixes any slew rate violations:

    a = time_ramp_up
    f = flat_time
    b = time_ramp_down
    grad_amp_new = -(a*s + b*z - 2*v)/(a+b+2*f)
    
    changed = True
    
    while changed:
        time_ramp_up_new = _to_raster(_calc_ramp_up_time(grad_amp_new, slew_rate, grad_start))
        time_ramp_down_new = _to_raster(_calc_ramp_down_time(grad_amp_new, slew_rate, grad_end))
        if time_ramp_up_new > time_ramp_up or time_ramp_down_new > time_ramp_down:
            a = time_ramp_up_new
            f = flat_time
            b = time_ramp_down_new
            time_ramp_up = time_ramp_up_new
            time_ramp_down = time_ramp_down_new
            
            grad_amp_new = -(a*s + b*z - 2*v)/(a+b+2*f)
        else:
            changed = False
    
    grad_amp = grad_amp_new

Even when only rounding up, grad_amp can increase in some cases. So you need to keep recalculating until the ramp times are no longer increasing.

[FrankZijlstra]

@schuenke I pushed the fixed version. I went back to using a grid search to solve for the timing. It's not as efficient as it can be, but for normal system parameters it is fine. For low slew rate (relative to maximum gradient strength), it can be slow to find a solution (but it will).

[schuenke]

@schuenke I pushed the fixed version. I went back to using a grid search to solve for the timing. It's not as efficient as it can be, but for normal system parameters it is fine. For low slew rate (relative to maximum gradient strength), it can be slow to find a solution (but it will).

Nice work 👍 While reviewing the code, I added some type hints, documentation and introduced some temp vars to improve the readability. Feel free to revert changes you don't like.

[FrankZijlstra]

Looks good. I tightened the bound for the linear search a bit by only considering a ramp down/up to gradient strength 0. After this, the desired area can always be reached by extending the duration (so binary search works there). Also did some minor cleaning.

I had a look at the dead space issue. Without considering the raster time it is possible to calculate the range where the area is unreachable with some ugly quadratic formula. But the bounds it find are not entirely accurate once the raster time is considered. I do think, that if this can be solved accurately, it gives only one or two possible solutions for the entire problem (i.e. find the points where the minimum or maximum area for a given duration is equal to the desired area). But for now I'll leave it as is, we can revisit it if the performance of the function turns out to be too poor for some situations.

[FrankZijlstra]

Made one final change that substantially narrows the timings that are checked by _find_solution.