Bug in: 'def gradient_waveforms' in 'sequence.py'

[Nikbert]

Hello, seems like we have a bug in 'def gradient_waveforms' in 'sequence.py'

I get the following error: " /usr/local/lib/python3.8/dist-packages/pypulseq/Sequence/sequence.py in gradient_waveforms(self) 637 grad_waveforms = np.zeros((3, maxt_i)) 638 --> 639 grad_waveforms[0, gx_i:(gx_i+gx.shape[0])] = gx 640 grad_waveforms[1, gy_i:(gy_i+gy.shape[0])] = gy 641 grad_waveforms[2, gz_i:(gz_i+gz.shape[0])] = gz

ValueError: could not broadcast input array from shape (766316,) into shape (766315,) "

Apparently the error can be resolved by changing line 632 in sequence.py from: maxt_i = int(np.max((tx_e, ty_e, tz_e))/dt) to: maxt_i = int(np.round(np.max((tx_e, ty_e, tz_e))/dt))

[btasdelen]

Hi @Nikbert. np.max((tx_e, ty_e, tz_e))/dt line assumes that time points will always be a multiple of the raster time (dt). It is strange that it was not the case here. I feel like it points to a bug in waveforms_and_times() function.

Would it be possible for you to provide a minimally working sequence that reproduces the same behavior?

[sravan953]

Thanks @bilal-tasdelen, for catching up with this (and multiple other issues!!!). @Nikbert , I too would request you to share code to reproduce the bug.

P.S: bilal, request you to join our Slack.

[Nikbert]

Hi, here is a link to the Jupyter notebook where we encountered the bug.: https://colab.research.google.com/drive/1klQFUabNEv9CmjgclRxTZYV887JSpEmC#scrollTo=ZSg_VUcvgZlo

[schuenke]

Hey @Nikbert, unfortunately the notebook requires access rights. Do you want to enable access to everbody or should we request permission individually?

[Nikbert]

Sorry, I changed the permissions.

[Nikbert]

I deleted most of the junk in the example and attached the python script here. It is not jet minimal, but maybe easier to work with. The example was run with "pypulseq.git@dev" version.

nw_splitgresequence4demo.py.zip

[schuenke]

I can reproduce the error and the suggested solution seems to fix it. I think it's a precision problem and int(x) always rounds to the smaller integer. Using round() or int(np.round()) both should work fine:

edit: In most cases it seems to be no problem because maxt_i is larger than g<x,y,z>.shape anyway.

[btasdelen]

@schuenke @Nikbert The issue does not happen with my fork, but happens with the dev branch. Meaning, I fixed it some time ago and forgot to create a pull request. I will find the difference and create a PR.

[btasdelen]

Found it! The bug was introduced in #87. I will re add np.rint() as a quick fix. But our assumption that time points must be on raster should hold. Long term, we need to find the real problem.

PS: Scratch my comment. I realized it is a floating point precision problem. Rounding is a proper fix.