Bug fixes in Sequence and make_trapezoid

[FrankZijlstra]

Sequence.waveforms_and_times:

• Fixed missing implementation for recompressing trapezoids from arbitrary gradients. Now proceeds without recompressing and added a TODO to implement the MATLAB code.
• Fixed issue where rounding errors could create non-monotonically increasing times. Also moved the non-monotonicity check outside the loop.

make_trapezoid:

• Fixed bug that would give an 'requested area is too large' error when explicitly specifiying a triangular gradient with duration, rise_time, and area.
• Fixed check for area == 0, while an unspecified area would be None

[btasdelen]

Hi @FrankZijlstra, thanks for the PR!

Would it be possible for you to open a separate PR for make_trapezoid fixes? The reason I ask is, it is a direct merge, but the Sequence fixes require closer inspection, so I did not want to hold it as it is unrelated to the rest of the patch.

Taking the non-monotonicity check out of the loop seems logical, as it is not dependent on the loop variable, but I wonder if the changing variables wave_data and wave_cnt has any unintended side effects. Have you verified that all the checking inside the loop is redundant? I will also have a look to double-check. If it is really redundant, we should also send a PR to the upstream Pulseq, as they also have it inside the loop.

[FrankZijlstra]

It's possible, but the fixes in make_trapezoid are a few edge cases you have to try pretty hard to run into. So I think they can wait.

Maxim Zaitsev already fixed and took the non-monotonicity check out of the loop in the local Pulseq repo here in Freiburg, based on me finding the bug that would give this warning last week (it's a rounding error where a+b+c != a+(b+c)). The reason the check was inside the loop in the first place was because he wanted to find out exactly when it was happening, but as a side effect it was spitting out hundreds of warnings once the check was fixed.

The loop is nothing more than a fancy way to concatenate the shape_pieces together. If it wasn't for the check whether the last element of the previous shape_piece has approximately the same timepoint as the first element of the current shape_piece, it could be implemented as:
wave_data[j] = np.concatenate([sp for sp in shape_pieces[j] if sp is not None], axis=1)
I'm sure the timepoint check could also be done in a more pythonic way. Another thing that could be done here is if we change the shape_pieces to a python list and use shape_pieces[j].append(np.array(...)) instead of shape_pieces[j, block_counter] = np.array(...), then you don't even have to check for None and just concatenate shape_pieces[j] directly.

[btasdelen]

@FrankZijlstra That makes sense, got it.

Yes, I agree that using numpy array for shape_pieces is counter-intuitive. We should do the conversion to a Python list and get rid of None checks in another PR.

I am also questioning the time check entirely. I feel like the only reason that check is needed is due to the usage of floating points and adding eps, which is not a good idea (as discussed here pulseq/pulseq#25). I am planning to convert the time representation to integer in nanoseconds, which should solve many of such problems.

The only other utility of time check is if the user supplied an out-of-order time vector. In that case, check should be done at the make_** functions.

[FrankZijlstra]

Yes, I agree. Because the shape_pieces are constructed in the loop before that, it seems like there are only a few possibilities for the times to go wrong outside of rounding errors, and I don't think they should:

• Regular trapezoid times are constructed using np.cumsum, and checks that the times added are > eps.
• Arbitrary gradients and RF pulses construct the time vector using np.arange and a timestep.
• Making extended trapezoids explicitly checks that the times are ascending.
  That only leaves utility functions like add_gradients to double-check.

That said, I think having the check is still nice to make sure no breaking changes are introduced. But in this case I think it might make more sense to have it as an assert instead of a warning.