Bug fixes and input checks for make_trapezoid

[btasdelen]

Tries to address #211.

Instead of throwing an error, I printed a UserWarning for when area and rise_time is supplied. I am not sure if there is a good case for solving an optimization to find the shortest trapezoid for a fixed rise_time and area, let me know if it is of interest.

For the second case, I believe there is a bug with setting fall_time when duration and rise_time is supplied. I also added a sanity check that covers the case in #211.

If this is a good time to overhaul input sanity checks for this function, let me know, we can brainstorm it. I find the current structure too nested and convoluted.

[schuenke]

Hey. Thanks for taking care of it @btasdelen

I think it would be the right moment to think about and implement all possible input arg combinations and add tests for all of them.

I won't find the time before thursday afternoon / friday, but would be happy to contribute.

[FrankZijlstra]

@btasdelen Can you explain the computation for the second case? Shouldn't it consider a flat time as well? By setting fall_time like this it seems that you implicitly set flat_time = 0.

[btasdelen]

@FrankZijlstra Oh shoot, you are right. Do we go back to setting fall_time=rise_time?

[FrankZijlstra]

I think that is the most reasonable assumption.

As for the structure of the function, maybe it would be a little clearer if all parameters are set to None and it explicitly checks every valid combination?
E.g.

if amplitude is not None and duration is not None:
   ...
elif flat_area is not None and flat_time is not None:
  ...
elif area is not None:
   ...
elif ...:

There would be more repeated code though, and each case would still need some nested if's for checking that no other ambiguous parameters are set and variations of the same case (e.g. rise_time/fall_time set).

[FrankZijlstra]

As for all the possible combinations, I wrote a little script using sympy to enumerate all that are non-ambiguous:

amplitude
amplitude + fall_time
amplitude + flat_time
amplitude + flat_time + fall_time
amplitude + rise_time
amplitude + rise_time + fall_time
amplitude + rise_time + flat_time
amplitude + rise_time + flat_time + fall_time
amplitude + duration
amplitude + duration + fall_time
amplitude + duration + flat_time
amplitude + duration + flat_time + fall_time
amplitude + duration + rise_time
amplitude + duration + rise_time + fall_time
amplitude + duration + rise_time + flat_time
flat_area
flat_area + fall_time
flat_area + flat_time
flat_area + flat_time + fall_time
flat_area + rise_time
flat_area + rise_time + fall_time
flat_area + rise_time + flat_time
flat_area + rise_time + flat_time + fall_time
flat_area + duration
flat_area + duration + fall_time
flat_area + duration + flat_time
flat_area + duration + flat_time + fall_time
flat_area + duration + rise_time
flat_area + duration + rise_time + fall_time
flat_area + duration + rise_time + flat_time
flat_area + amplitude
flat_area + amplitude + fall_time
flat_area + amplitude + rise_time
flat_area + amplitude + rise_time + fall_time
flat_area + amplitude + duration
flat_area + amplitude + duration + fall_time
flat_area + amplitude + duration + rise_time
area
area + fall_time
area + flat_time
area + flat_time + fall_time
area + rise_time
area + rise_time + fall_time
area + rise_time + flat_time
area + rise_time + flat_time + fall_time
area + duration
area + duration + fall_time
area + duration + flat_time
area + duration + flat_time + fall_time
area + duration + rise_time
area + duration + rise_time + fall_time
area + duration + rise_time + flat_time
area + amplitude
area + amplitude + fall_time
area + amplitude + flat_time
area + amplitude + flat_time + fall_time
area + amplitude + rise_time
area + amplitude + rise_time + fall_time
area + amplitude + rise_time + flat_time
area + amplitude + duration
area + amplitude + duration + fall_time
area + amplitude + duration + rise_time
area + flat_area
area + flat_area + fall_time
area + flat_area + flat_time
area + flat_area + flat_time + fall_time
area + flat_area + rise_time
area + flat_area + rise_time + fall_time
area + flat_area + rise_time + flat_time
area + flat_area + duration
area + flat_area + duration + fall_time
area + flat_area + duration + rise_time
area + flat_area + amplitude
area + flat_area + amplitude + fall_time
area + flat_area + amplitude + rise_time

Of course, not all of these make sense to implement... E.g. only specifying amplitude would make a max_slew triangular blip up to the specified amplitude, but I can't think of a reason you would ever want that specified that way.

[btasdelen]

@FrankZijlstra I think a good portion of these options at the end just calculates the missing information, and at the end, there are actually only a couple of way to design the gradient. To avoid too much code repetition, what I have in mind is each condition may fill the missing information, and sets a variable that determines which way will be used (i.e. by area, or by amplitude).

That is way too many combinations. I suggest only implement the ways that are actually useful, but also warn the user for any other case. If someone somewhere has an edge case that needs a specific combination of input, they can let us know and we can add it. Higher the complication, worse the bugs.

I will have time to do these in 2-3 days. In the meantime, please let me know if there are ideas.

[FrankZijlstra]

Looking at the possible combinations, I think the main cases are (where the ones marked with (x) is what the current code more or less covers):

amplitude (x)
amplitude + duration (x)
flat_area (x)
flat_area + duration
flat_area + amplitude
area (x)
area + duration (x)
area + amplitude (x)
area + flat_area

I don't think that is too many per se.

Within these cases, the timing would have to be figured out, depending on which timing parameters are specified. The main question then is whether the timing is fully specified, or whether something will need to be inferred. It's possible there is some common code related to the timing that can be put before checking which case we are using. And if you consider duration as part of the timing, there are only 6 cases left (5, if you exclude the exotic area + flat_area case).

Splitting the case explicitly would also be easier to test. Knowing which of amplitude, area, and flat_area are specified, you know which code block you end up in. In the current code, if I specify amplitude, I can end up in all 3 cases of the if block, with the last one seemingly ignoring the amplitude.

[schuenke]

Thanks for rewriting the function @btasdelen. I will definitely have a look at it and go through evertything in detail, but not sure when I find the time. But I guess it's not super urgent anyway, is it?

[btasdelen]

@schuenke I don't think it is urgent, there are a couple of bug fixes regarding the input checks, but nothing critical. Overall calculation and the logic is the same as before.

A short summary of the rework:

• The tests are more comprehensive and instead of just checking the return types, checks the output values.
• All default inputs are None, rather than 0, which might be a valid input value.
• Input sets that are reasonable, but not implemented raise NotImplementedError.
• More strict input checks.
• Logic flow is reworked to be more intuitive to allow easier expansion in the future.

There are small things left to do:

• Need to check code coverage by the tests to make sure we don't forget any edge cases.
• All tests pass, but I never tried this with real sequences.
• If there are any input set that looks useful to have, we can also try implementing them in this PR.

[schuenke]

Currently, if we provide area=0, which often happens in phase encoding loops etc, a total gradient with duration = 2 * grad_raster_time is returned. Should we catch this in calculate_shortest_params_for_area?
For example using something like:

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."""

    if area == 0:  # typical case in phase encoding loops
        # If the area is 0, return all times and amplitude as 0
        return 0.0, 0.0, 0.0, 0.0
    ...

[btasdelen]

Currently, if we provide area=0, which often happens in phase encoding loops etc, a total gradient with duration = 2 * grad_raster_time is returned. Should we catch this in calculate_shortest_params_for_area? For example using something like:

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."""

    if area == 0:  # typical case in phase encoding loops
        # If the area is 0, return all times and amplitude as 0
        return 0.0, 0.0, 0.0, 0.0
    ...

How about we check in make_trapezoid? Similarly, if any of the area, amplitude or duration is 0, that should result in a zero duration gradient. That can be a separate calc_path.

[FrankZijlstra]

Does the interpreter on the scanner accept zero duration gradients? I don't think I ever tried this. And perhaps there are slew rate checks that would fail in our code?

But even if it is accepted, it should only be returned if the user did not specify any other timing. area=0, rise_time=raster_time should still return the 2*raster_time gradients, for example.

[schuenke]

Does the interpreter on the scanner accept zero duration gradients? I don't think I ever tried this. And perhaps there are slew rate checks that would fail in our code?

But even if it is accepted, it should only be returned if the user did not specify any other timing. area=0, rise_time=raster_time should still return the 2*raster_time gradients, for example.

I was wondering the same. I know that rf pulses with amp = 0 were and probably still are a problem. And blocks with only labels and no duration can also cause problems, no? So maybe it's safer to keep it the way it is, but I will still try out the 0 duration gradients on the scanner next week. Hopefully...

[btasdelen]

This PR contains bug fixes, too. I suggest we merge this, and create a new PR for future discussion, if you agree @schuenke?

[schuenke]

lgtm

[FrankZijlstra]

Looks good to me for now, though it could use some more tests. I think combined with the sequence examples tests in #220 we should have a better idea whether anything broke.