imr-framework / pypulseq

Pulseq in Python
https://pypulseq.readthedocs.io
GNU Affero General Public License v3.0
112 stars 58 forks source link

Hacky fix for waveforms_and_times rf append floating point issue. #114

Closed btasdelen closed 11 months ago

btasdelen commented 1 year ago

When we have RF pulses with non-zero start/end times, waveforms_and_times() function tries to add trailing/leading zeros as necessary to the beginning and end of the RF waveform. To avoid having duplicate times, leading zero is added at time t[0]-eps, trailing zero is added to the t[-1]+eps, where eps is supposedly "smallest floating increment". These are all inline with upstream MATLAB Pulseq implementation.

Trouble begins with the eps. eps is not actually the smallest increment. Floating point have variable increments depending on the number of significant digits. Quoting numpy documentation:

eps: float

The difference between 1.0 and the next smallest representable float larger than 1.0. For example, for 64-bit binary floats in the IEEE-754 standard, eps = 2**-52, approximately 2.22e-16.

It is only true if we are close to 1.0. Time vector is an unbounded vector, so our comparison with eps becomes incorrect as we add eps to larger and larger floats. When we apply np.unique() to these values, they are identified as duplicates, thus we raise:

raise Warning(
                           "Not all elements of the generated time vector are unique."
                        )

I could not find a solution to fix this without changing the way we look at things here, so I just added 50*eps arbitrarily to silence the error. I am open to suggestions on how to solve this properly. I am also not sure why we need to add trailing or leading zeros, maybe the answer lies removing that.


The second issue this PR addresses is, using np.unique is extremely slow, and prohibits the use of the function for large sequences. It seems like it is used as a check to see if we have repeating time points. I would say it is an overkill. We do not expect to have an arbitrarily ordered time vector. If there is a repeat, it is most likely happening at consecutive time points. Thus, I replaced np.unique with a call to np.diff and a comparison with eps. This will make sure that time points are monotonically increasing. This way, I observed around 2.7 times improvement in speed. np.diff is still slow, in the near future, a better approach is needed.

Sorry for the long post, let me know what you think.