imr-framework / pypulseq

Pulseq in Python

imr-framework

pypulseq's Issues

Gridding a sequence containing trapezoidal gradients results in 'cut off' corners and a time-grid shift of grad_raster_time/2

Describe the bug
When creating a sequence that contains trapezoidal gradient lobes using make_trapezoid , the result of Sequence.gradient_waveforms() shows an incorrectly gridded waveform. When gridded (by calling the function Sequence.gridded_wave_forms() which in subsequently calls the function `pypulseq.points_to_waveform'), the waveform timing starts at grad_raster_time/2. When defining trapezoids such that the interval border of the piecewise linear function lie exactly on a regular grid with spacing grad_raster_time, this results in 'capped off' corners, a non-zero start value as well as a bend waveform in the last interval.

Not sure why the function pypulseq.points_to_waveform performs the time-grid shift of grad_raster_time/2 but for this case it should not be applied.

So, if this is intentional in some cases this conditional needs to be revised to catch the correct case for the trapezoidal gradients.

To Reproduce

import pypulseq
import matplotlib.pyplot as plt
import numpy as np

system = pypulseq.Opts(max_grad=80, grad_unit='mT/m', max_slew=200,
                       slew_unit='mT/m/ms', grad_raster_time=0.01)
seq = pypulseq.Sequence(system=system)
l1 = pypulseq.make_trapezoid(channel="x", flat_time=5.6, amplitude=1., rise_time=0.9)

seq.add_block(l1)
wf = seq.gradient_waveforms()
t = np.arange(0, seq.duration()[0]+seq.grad_raster_time, seq.grad_raster_time)

f, a = plt.subplots(1, 1)
a.plot(t, wf.T)
a.grid(True)

Expected behavior
My expectation is, that the reference-points of normalized amplitude [0, 1, 1, 0] at times [0, rise_time, rise_time+flat_time, rise_time+flat_time+fall_time] should be met exactly if all time-points match the grad_raster_time as shown in the figure below.

Desktop (please complete the following information):

• OS: [Windows 11 & Ubuntu 20.04]
• pypulseq version: 1.3.1.post1

Additional context
To catch errors like this i propose to add a test similar to the following:

import unittest
import pypulseq
import numpy as np

system = pypulseq.Opts(max_grad=80, grad_unit='mT/m', max_slew=200,
                       slew_unit='mT/m/ms', grad_raster_time=0.01)


class TestGradientWaveform(unittest.TestCase):
    def test_correct_start_trap(self):
        seq = pypulseq.Sequence(system=system)
        l1 = pypulseq.make_trapezoid(channel="x", flat_time=5.6, amplitude=1., rise_time=0.9,
                                     delay=0)
        seq.add_block(l1)
        wf = seq.gradient_waveforms()
        self.assertTrue(np.allclose(wf[:, 0], 0, rtol=1e-10),
                        msg="First gridded sample is not zero")
        self.assertTrue(np.allclose(wf[:, -1], 0, rtol=1e-10),
                        msg="Last gridded sample is not zero")

Return values of make_xxx_pulse functions

Since 7542b39 the make-pulse-functions like make_sinc_pulse.py and make_gauss_pulse.py return a different number of variables as before when using return_gz = False.

For compatibility with previous versions I suggest to return Nones for gz and gzr as it was the case before the release of v1.3.1. However, I think using the _return_gz' flag instead of the try block is a good idea and I suggest to change the return lines only. For example the make_sinc_pulse.py could be modified to:

if return_gz:
    return rf, gz, gzr
else:
    return rf, None, None

This would also solve issue #46

Let me know if I should include this in the PRs for issue #46 and issue #48.

define __version__ in top level __init__.py

Is your feature request related to a problem? Please describe.

It is best practice to define a __version__ for your module at the top level. See PEP 396

Describe the solution you'd like
This can either be done manually or using an automated tool like versioneer.

Describe alternatives you've considered

Additional context

Shape check for arbitrary rf/grad

Is your feature request related to a problem? Please describe.
If a waveform of shape (1,x) is used to construct arbitrary rf or gradient events they cannot be successfully added to a block.

Describe the solution you'd like
Perform shape checking

broken CONTRIBUTING link in README

Describe the bug
The link to CONTRIBUTING.md in the README is broken (.md extension is missing)

Indexing mistake in seq read helper function

Describe the bug
Around line 200 in sequence.py, in the function rf_from_lib_data() , there has been a change between version 1.2.1 and 1.2.0 that seems to fix the error message but reads the sequence incorrectly.

# 1.2.1 (gives error message in cases that lib_data has length less than or equal to 6)
        rf.delay = lib_data[3]
        rf.freq_offset = lib_data[4]
        rf.phase_offset = lib_data[5]
        
        if max(lib_data.shape) < 6:
            lib_data = np.append(lib_data, 0)
        rf.dead_time = lib_data[6]

        if max(lib_data.shape) < 7:
            lib_data = np.append(lib_data, 0)
        rf.ringdown_time = lib_data[7]

        if max(lib_data.shape) < 8:
            lib_data = np.append(lib_data, 0)

        use_cases = {1: 'excitation', 2: 'refocusing', 3: 'inversion'}
        if lib_data[8] in use_cases:
            rf.use = use_cases[lib_data[8]]

# 1.2.0 (current version; reads the data incorrectly)
        rf.delay = lib_data[3]
        rf.freq_offset = lib_data[4]
        rf.phase_offset = lib_data[5] ## <<< note this 

        if max(lib_data.shape) < 6:
            lib_data = np.append(lib_data, 0)
        rf.dead_time = lib_data[5] ## <<< and note this

        if max(lib_data.shape) < 7:
            lib_data = np.append(lib_data, 0)
        rf.ringdown_time = lib_data[6]

        if max(lib_data.shape) < 8:
            lib_data = np.append(lib_data, 0)

        use_cases = {1: 'excitation', 2: 'refocusing', 3: 'inversion'}
        if lib_data[7] in use_cases:
            rf.use = use_cases[lib_data[7]]

For reference, the fields contained in lib_data in the case of an RF pulse are:
lib_data[0] : RF amplitude
lib_data[1] : RF magnitude waveform ID
lib_data[2] : RF phase waveform ID
lib_data[3] : RF delay
lib_data[4] : RF frequency offset
lib_data[5] : RF phase offset

Right now (1.2.0), the code is reading the phase offset as the dead time.

Desktop (please complete the following information):

• OS: Windows 10 Pro
• Version: pypulseq 1.2.0

list dependencies in README and/or docs

1.) a JOSS criteria asks whether dependencies are clearly stated.

I can find them by looking in setup.py, but it would be more apparent to the user if these were listed briefly in the top-level README near where the pip install commands are.

2.) Also, not strictly related to the dependencies, but the API docs are a bit hidden within the README. It may be worth adding a link to the readthedocs page (and optionally some keywords to aid discoverability) at the top of the GitHub page. This is done on GitHub itself (i.e. I think you will have an "Edit" button near where it currently says "No description, website, or topics provided. ")

UnboundLocalError in make_sinc_pulse.py

Since 7542b39 I get the error UnboundLocalError: local variable 'gz' referenced before assignment when creating an rf pulse with make_sinc_pulse.py.

This happens when I set return_gz = False or when I don't set any return_gz value and thus use the False default value.

matplotlib backend TkAgg in sequence problematic in other backends

The hardcoded setting of TkAgg right on top of pypulseq.Sequence.sequence with mpl.use('TkAgg') can lead to problems with another backend like QT5. It raises:
ImportError: Cannot load backend 'TkAgg' which requires the 'tk' interactive framework, as 'qt5' is currently running

In many cases, this can be helped by manually changing the framework to 'tk', if you read up on it. However, we like to read sequence files in a QT5 framework (a Qt5 based application). You can reproduce it by setting the matplotlib framewrk to 'qt5' or other than 'tk'.

Is the TkAgg backend absolutely necessary? If so, would it be possible to separate the backend switch (and maybe the plot function) from the rest of the sequence code?

It would be best if the plot is not dependent on the backend, but it would suffice for us it the functions that do not depend on matplotlib can be imported without the backend switch.

Thanks for your help!

typo in makearbitrary_rf.py

system.rfRasterTime instead of system.rf_raster_time (Line 44)

JOSS Review: Sequence version out of sync with release version.

Current release version is 1.2.2 (as seen here)

...but the generated pulse sequence files report version 1.2.1, because of the hardcoded values here:

pypulseq/pypulseq/Sequence/sequence.py

Lines 20 to 24 in 2268bf6

	class Sequence:
	def __init__(self, system=Opts()):
	self.version_major = 1
	self.version_minor = 2
	self.version_revision = 1

So you might want to (1) fix this release sync issue, and (2) find a better way to update this value. Another point is that since you are always pushing directly to master, people area likely to have the same release version but different versions of the code, if they clone from github.

At a bare minimum, you might want to block devs from pushing directly to master (there's a setting on github to do that), and only develop on branches, where you can then have guidelines to pull request reviewers to check this value whenever reviewing pull requests. But ideally, you would have this updated in an automated way, but I don't know if that's possible.

Error during sequence write when Sequence contains no gradients

Describe the bug
When writing a Sequence which doesn't contain gradients, an error occurs:

~\Anaconda3\envs\mri\lib\site-packages\pypulseq\Sequence\write_seq.py in write(self, file_name)
83
---> 84 if any(arb_grad_mask):
85 output_file.write('# Format of arbitrary gradients:\n')
86 output_file.write('# id amplitude shape_id delay\n')`

TypeError: 'bool' object is not iterable

This happens, because when a Sequence doesn't contain gradients, the EventLibrary self.grad_library is empty and in Line 80 of write_seq.py the variable grad_lib_values will contain a numpy array object as:
array([], dtype=float64)

The empty numpy array can be compared (boolean) element-wise with objects which can be interpreted by numpy as a float data type, for example the integer number 1 or the float 3.3. This will result in an empty boolean array as expected:
array([], dtype=bool)

However, when comparing the grad_lib_values array with a string / character 'g' or 't', numpy can not interpret the string as a float and the '==' operator will be applied as an object comparison, which in turn results in a boolean 'False'.

any(False) will throw the error shown above.

Suggested fix
To overcome this problem, one suggestion is to create the boolean array grad_lib_values inside a try-catch-block using np.equal (which will then throw an exception numpy.core._exceptions.UFuncTypeError: ufunc 'equal' did not contain a loop with signature matching types (dtype('<U32'), dtype('<U32')) -> dtype('bool')) or to use the if any(...) inside a try-catch-block and skip / notify if the exception above occurs.

Desktop

• OS: Windows 10
• pypulseq version: 1.3.1.post1

Default params for make_sinc_pulse

def make_sinc_pulse(flip_angle: float, system: Opts = Opts(), duration: float = 0, freq_offset: float = 0,
phase_offset: float = 0, time_bw_product: float = 4, apodization: float = 0,
center_pos: float = 0.5, max_grad: float = 0, max_slew: float = 0, slice_thickness: float = 0,
delay: float = 0, use: str = None):

duration cannot default to 0

ValueError in make_block_pulse.py when not difining the use parameter

Since 7542b39 the make_block_pulse.py throws the ValueError ValueError: Invalid use parameter. Must be one of 'excitation', 'refocusing' or 'inversion'. Passed: when use is not defined. The reason is that the default for use is an empty string now, but in line 64 you check for a None. This can be fixed by changing line 64 to:

if use != '' and use not in valid_use_pulses:

Non-selective block pulse not playing

Interestingly, the block pulse without any gradient does not play but with the gradient goes through.
System - Siemens 3T Prisma Fit, 32/130, different combinations of ringdown time and rf deadtime

Wrong time unit in function documentation

The time unit in the description of more or less all event block objects (such as pulses, gradients,...) is wrong. The time unit in PyPulseq is seconds [s] instead of milliseconds [ms] as stated in the descriptions. This might be confusing especially for new users.

Bug in block length when using arbitrary RF pulses or gradients

Describe the bug
Sequences won't play when any of the blocks have arbitrary RF or gradients with the same length. Based on trial/error experiments, seems like the gradient waveforms have to have a slightly longer duration compared to the other components of the block (RF pulse signal, adc).

To Reproduce
Create an arbitrary gradient waveform using make_arbitrary_grad
Design ADC to have the same number of samples as the gradient created
seq.add_block(gradient,adc)

Expected behavior
The sequence won't play with the last version of the interpreter.

Screenshots
If applicable, add screenshots to help explain your problem.

Desktop (please complete the following information):

• pypulseq version: 1.2.0.post1

t_factor is missing in the ADC plot

The t_factor for scaling the plot of ADC points is missing in sequence.plot. This results in a wrong scaling of the ADC plot for 'ms' or 'us' representations

feature: set sequence plots to share X axis

Is your feature request related to a problem? Please describe.

The plot method of the Sequence class would be more useful if the
sharex attribute of the subplots was set. This would make it so that if
the user zooms in on one subplot, the others would also be zoomed by the
same amount.

Without this feature, it is hard for the user to inspect the relative
timing of events.

Describe the solution you'd like
If unfamiliar, a figure with subplots sharing an x axis can be created use notation such as:
fig, axes = plt.subplots(3, 1, sharex=True)
axes[0].plot(t, gx)
axes[1].plot(t, gy)
axes[2].plot(t, gz)

Describe alternatives you've considered
the figure add_subplot method can also take a sharex argument if you are adding the subplots one by one instead of creating them upfront as in the example above

Additional context

Integer definitions are missed in `write_seq.py`

Describe the bug
Originally reported by @mavel101 from #45:
The set_definition function accepts an integer as value, but in the write_seq an integer is only accepted as part of a list, tuple or ndarray not as a single integer. Changing this line in the write_seq:

elif isinstance(values[block_counter], float):
                output_file.write(f'{values[block_counter]:0.9g} ')

to

elif isinstance(values[block_counter], (int,float)):
                output_file.write(f'{values[block_counter]:0.9g} ')

would solve this.

To Reproduce

• Share the pulse sequence code, if possible.
• If not, isolate the bug to a particular pulse sequence event by iteratively writing/plotting the sequence one event at a time.
• If the sequence successfully writes/plots with individual pulse sequence events, isolate the bug by writing/plotting the sequence by successively adding events.

Desktop (please complete the following information):

• pypulseq version: 1.3.1

Dwell time on write_epi_rs.py

Describe the bug
The current implementation does not check if the dwell time complies with the hardware requirements - i.e. multiple of 10us

The code is also not easy to read as some of the system requirements are hard-coded instead of being defined with the Opts function.

feature: define a top-level API

Is your feature request related to a problem? Please describe.
It may be worth thinking about the API from a user perspective. If I
import pypulseq, the only things that is defined in the top-level namespace
is a name string. This means that the user has to specifically know where
each function lives to use the library. For example, from write_epi.py the
imports at the top look like:

from pypulseq.Sequence.sequence import Sequence
from pypulseq.make_adc import make_adc
from pypulseq.make_sinc_pulse import make_sinc_pulse
from pypulseq.make_trap_pulse import make_trapezoid

Describe the solution you'd like
It would be easier for the user if these were available at the top level as

from pypulseq import Sequence, make_ads, make_sinc_pulse, make_trapezoid

This would also allow interactive users in Ipython to type pypulseq. and get an autocompleted list of the available functions/classes.

Describe alternatives you've considered

Additional context

Typo in make_arbitrary_rf

Line 72: signal = signal / bp.abs(np.sum(signal * system.rf_raster_time)) * flip_angle / (2 * np.pi)

Change bp to np

float compatibility problems

For our sequence design of long sequences, we're experiencing problems with the np.float16 declaration in the compress_shape function (line 32 of compress_shape.py).

Could you change it to np.float32 or np.float_? That resolves the issue.

Mismatch in the section name for gradients in seq.read and seq.write

Describe the bug
The seq.write function generates files that cannot be read by seq.read because the gradient section name is different between the two ([GRAD] according to the write function and [GRADIENTS] as expected by the read function).

To Reproduce
Run the following code (under Pypulseq 1.2.1):

from pypulseq.Sequence.sequence import Sequence
from pypulseq.make_sinc_pulse import make_sinc_pulse
from pypulseq.make_extended_trapezoid import make_extended_trapezoid
from pypulseq.opts import Opts
from math import pi

system = Opts()
seq = Sequence()
g = make_extended_trapezoid(channel='x',times=[0,0.01],amplitudes=[1,2],system=system)
seq.add_block(g)
seq.write('test.seq')
seq2 = Sequence()
seq2.read('test.seq')

Expected behavior
The sequence should be read without generating an error.

Error message
When run, this error appears:

  File "(personal system path)\pypulseq\Sequence\read_seq.py", line 57, in read
    raise ValueError(f'Unknown section code: {section}')
ValueError: Unknown section code: [GRADIENTS]

Desktop (please complete the following information):

• OS: Microsoft Windows 10 Pro
• Version: pypulseq 1.2.1

Adding delay is breaking seq.report

Describe the bug
seq.report() is unable to detect a Cartesian trajectory for a spin echo implementation

To Reproduce

• Share the pulse sequence code, if possible.
  for i in range(Ny):

  gy_pre = make_trapezoid(channel='y', area=phase_areas[i], duration=t_pe, system=system)
  #Excitation blocks
  seq.add_block(rf_ex, gz_ex)
  seq.add_block(gz_reph, gx_pre, gy_pre)
  #seq.add_block(delay_TE1)

  #Refocusing and phase encoding
  seq.add_block(rf_ref, gz_ref)

  #Delay and readout
  #seq.add_block(delay_TE2)
  seq.add_block(gx,adc)

report = seq.report()
print(report)

Expected behavior
Needs to show that it is 2D Cartesian with spatial resolutions only in x and y

Screenshots Output from report API - wth delay commented
Number of blocks: 512
Number of events:
RF: 256
Gx: 256
Gy: 128
Gz: 384
ADC: 128
Delay: 0
Sequence duration: 1.840640 s
TE: 0.009571 s
TR: 0.014380 s
Flip angle: 90.00 180.00 deg
Unique k-space positions (aka cols, rows, etc.): 128 128 1
Dimensions: 3
Spatial resolution: 1.95 mm
Spatial resolution: 1.95 mm
Spatial resolution: 31.25 mm
Repetitions/slices/contrasts: 1.0
Cartesian encoding trajectory detected
Block timing check passed successfully
Max gradient: 107246 103226 666667 Hz/m == 2.52 2.42 15.66 mT/m
Max slew rate: 5362275148 5333333333 5333333333 Hz/m/s == 125.95 125.27 125.27 mT/m/s
Max absolute gradient: 666667 Hz/m == 15.66 mT/m
Max absolute slew rate: 7.90635e+09 Hz/m/s == 185.70 T/m/s

Desktop (please complete the following information):

• macOS
• 10.14.6

Additional context
If I uncomment the delays, the report cannot evaluate block timings;
Block timing check failed. Error listing follows:

ValueError if triangular gradient is defined with rise_time, flat_time and amplitude

If a triangular gradient is defined with make_trap_pulse and only amplitude, rise_time and flat_time (=0) are provided, a ValueError is raised, stating that an area has to be provided. However, the triangular gradient is fully described with rise_time and amplitude.
This behaviour only occurs in the newest version, where the default value for flat_time was changed from 'None' to '0'.

Seq.write() changes the seq object

Describe the bug
Using seq.write() to output the .seq file causes numerical values in the seq object to change due to unit conversion. The unexpected outcome is that when seq.write() is run twice, the second generated sequence is incorrect.

To Reproduce

1. Go to the write_spgr.ipynb notebook and run all cells ahead of section 6.
2. Run this cell in section 6 twice
3. Compare the 2 seq files output from the 2 runs.

Expected behavior
There should be scaling differences in the timing of gradients, for example.

pulseq 1.3.0 compatibility

I am working on an application where I will probably need to read and write version 1.3 seq files (or have compatibility for both).

Are you working on the implementation of that compatibility? I am happy to contribute in any case, just wanted to check with you first.

test_report method of Sequence should return the report instead of None

Describe the bug
I assume the test_report method of Sequence should return the string corresponding to the report?

To Reproduce
The example file write_gre.py generates a report via report = seq.test_report(). The variable report in this example ends up as None

Expected behavior
returns a string containing the report

Screenshots

Desktop (please complete the following information):

• OS: [e.g. macOS]
• Version [e.g. 1.2.1]

Additional context

Functions to make combining gradients easier

Is your feature request related to a problem? Please describe.
When I am coding oblique gradients, and with multiple types of gradients being played together, there is a need to add gradients together, scale gradients, or change the channels of an already made gradient object.

• A typical example is when readout, phase encoding, and slice refocusing gradients overlap. It's straightforward when orthogonal encoding is used - each type would corresponding to one channel (either 'x', 'y', or 'z'). However, when you want general oblique encoding, each type has components in x, y, and z. The add_block function, however, only accepts one gradient per channel (it overwrites the previous one if more than one is used).
• Scaling: in radial k-space readouts, for example, you'd want to scale the x and y gradients while maintaining their shape. Right now, existing sequences simply create new gradients each iteration. It would be nice to have functions dedicated to scaling gradients and changing their channels as part of PyPulseq.

Describe the solution you'd like
Include these utility functions:

1. Gradient scaling / channel modification. Returns a changed copy of the input gradient. Detects violations of max gradient or max slew rate and raises corresponding errors.
2. Combining any number of gradients with arbitrary channels into 3: gx,gy, and gz. Preserve trapezoidal gradients if timings match. Otherwise, return general gradients.
3. Detect zero gradients during add_block: don't include them to reduce storage waste

Describe alternatives you've considered
The alternative is for each user to perform these actions in their custom code, individually. I believe they are generally useful enough to be part of Pypulseq.

Folder structure

The import statements require a "pulseq" folder underneath the pypulseq folder, which is missing or the import statements should read "pypulseq"

Previous working TSE script generates error

Describe the bug
I tried to run the TSE notebook that previously worked under the version 1.2.0post4; with the current 1.3.1post1, it produces an error at the first make_sinc_pulse

To Reproduce

• Go to write_tse.ipynb: https://github.com/tonggehua/mr-sequence-validation/blob/main/TSE/write_tse.ipynb
• Run the code with PyPulseq 1.3.1post1 or master at installation

Expected behavior
This error should happen at the very first make_sinc_pulse:

    rf_ref, gz, _ = make_sinc_pulse(flip_angle=flip_ref, system=system, duration=t_ref, slice_thickness=thk,
TypeError: cannot unpack non-iterable types.SimpleNamespace object

Desktop (please complete the following information):

• OS: Windows 10
• pypulseq version: 1.3.1post1

Question related to SE and TSE

Hello,
I am now in this field, I have a question and request related to some examples.

1. In the example of https://github.com/imr-framework/pypulseq/blob/master/pypulseq/seq_examples/notebooks/write_t2_se.ipynb, can you please tell me the TE and TR is defined in the beginning as the actual echo and repetition time or TE=delay1+delay2, TR=delay_TR is the actual echo and repetition time?

2. If we want to change the BW then should we need to change the readout_time right? BW=1/readout_time?

3. can you please also change this TSE example https://github.com/imr-framework/pypulseq/blob/master/pypulseq/seq_examples/scripts/write_tse.py with a different TF (Turbo factor ) parameter?

best,
Hafiz

Function to appends the md5_hash signature to the .seq file to make it compatible with pulseq Version 1.4

Is your feature request related to a problem? Please describe.
Pulseq Version 1.4 requires a md5_hash signature at the end of the .seq file

Describe the solution you'd like
A function that appends the signature to make the .seq file compatible with Pulseq V1.4.
Calling the function "add_signature.add_md5_sig(file_name)" after the "seq.write()" function will append the required signature.

Describe alternatives you've considered
An alternative would be to include the function in the seq.write function.
This would be the preferred solution in the future.

There is a pull request for the described add_signature.py function

adding a test suite and continuous integration

Is your feature request related to a problem? Please describe.
In considering the JOSS review criteria, I think the lack of a test suite or any type of continuous integration to verify correct operation across platforms / software versions is the biggest weakness.

There are fortunately a number of demos that the user can run to verify the basic operation of the package, but I did see any automated way to run tests.

Describe the solution you'd like

Ideally some automated unit tests would be created. The pytest package is a popular choice for this.

Once unit tests have been defined, services like Travis-CI, Circle-CI, Appveyor and/or Azure can be used to set up automated cross-platform testing on each pull request, etc. This may be overkill at the moment for a project of this size, but something to keep in mind.

Describe alternatives you've considered

Additional context

Plot Sequence with trigger

Describe the bug
I'm using the make_trigger function to create a trigger object, this way:

trig = make_trigger(channel='physio1', delay=100e-6, duration=Tdelay_trig, system=system)

But then when I tried to plot the sequence with seq.plot() I got this error:

To Reproduce
I used your gre_walkthrough to make sure there was not another cofounding factor, you can see the code here: https://colab.research.google.com/drive/19MruYk0W0AMulph78q_yqywCemuLK0jW?usp=sharing

Desktop (please complete the following information):

• pypulseq version: 1.3.1.post1

Thanks a lot!

broken link to pypulseq spec

There are several broken links in the README.md as rendered by GitHub. For example, the Pulseq spec directs me to:
https://github.com/imr-framework/pypulseq/blob/master/pulseq-spec
which does not exist

However, at the bottom of the raw file it looks like a working link is correctly specified:
http://pulseq.github.io/specification.pdf

I think the issue is that Markdown is expecting a square bracket around the link variable name.
e.g. [Pulseq specification](pulseq-spec) should be [Pulseq specification][pulseq-spec], etc.

seq.read does not support read of previous versions of '.seq' files

Regarding the 'seq.read()' function we realze that we could not read '.seq' files that were created with previous versions of the pypulseq toolbox. This means that we had to generate the '.seq' file again with the latest version of the pypulseq toolbox, in order to read it.

Loosen pinning of numpy & matlplotlib versions in setup.py

Is your feature request related to a problem? Please describe.

Dependencies are pinned too tightly. It would be best to not force users to install a specific micro release of both NumPy and Matplotlib to install pypulseq. As Numpy is at the base of the scientific python ecosystem, we should attempt not to require users to change their version unless strictly necessary.

Similarly, I see that the package is listed as Python 3.6 only in setup.py, but it seems to also work for 3.7, so that should probably be added as well. (I see there are some modern code features like type annotations that may prevent use of earlier versions).

Describe the solution you'd like
It would be better to identify a minimum version (or version range) that are compatible with pypulseq. Matplotlib seems to be only needed for the plotting capabilities and could potentially be converted to an optional runtime requirement.

Describe alternatives you've considered
A clear and concise description of any alternative solutions or features you've considered.

Additional context
Add any other context or screenshots about the feature request here.

RF pulse phase offset not taken into account when plotting the sequence

Changing the phase of the RF pulse using the "phase_offset" argument doesn't seem to change the corresponding sequence plot.
This value should also be considered when plotting the phase of the RF pulse.

Missmatch between compress_shape function in Python and MATLAB

In 2268bf6 you changed the quantization factor in compress_shape.py from 1e-7 (as in MATLAB) to 1e-8. This leads to deviations between seq-files created in MATLAB and python.

For our pulseq-cest project it is important that the seq-files created with MATLAB and python match as good as possible. Therefore, we would prefer to change the factor back to 1e-7.

Is/was there a specific reason why you decided to deviate from the MATLAB implementation in this case and is there any change that you change it back to the old value?

Unfortunately this issue is pretty urgent for us, because we need to create a release for a publication.

Thanks again for your time and help!

Best,
Patrick

Error generating Trapezoid gradient

When trying to generate a trapezoid gradient by providing only its axis, amplitude and duration I get the following error:

UnboundLocalError: local variable 'amplitude_temp' referenced before assignment

Output sequence waveforms from .seq files

Enhancement Description
It would be helpful to be able to output entire waveforms (i.e. all blocks concatenated together) for:

• RF Magnitude
• RF Phase
• Gx
• Gy
• Gz
• ADC points
• Time

Describe the solution you'd like
A Class function for seq objects.
It would be nice to be able to specify (1) Time range and (2) Block index range for custom outputs for examination and debugging.

Describe alternatives you've considered
N/A

Indexing issues happen during seq.read ()

In line 158 - 171 of sequence.py, the code tries to access an index that is out of range in lib_data and causes an error when one calls seq.read() on a file created by the same PyPulseq version (1.2.1).

Those indices should change from n to n - 1.

more pythonic file organization

Is your feature request related to a problem? Please describe.

This is a low priority issue and ultimately a matter of taste, but I thought I would offer it for consideration.

The files in pypulseq seem to be nearly a 1:1 mapping of the corresponding Matlab files from pulseq. It would perhaps be more Pythonic to group related functions into a single file rather than having one file per function.

Describe the solution you'd like

As a concrete example, the RF pulse generating functions:
make_arbitrary_rf
make_gauss_pulse
make_sinc_pulse
make_block_pulse
could be grouped within a single module called rf.py.

Similarly, several gradient related files could be merged into something like grad.py

Describe alternatives you've considered

Additional context

Bug in makearbitrary_rf.py

No option available when nargout=1 (no gz as output)

vds_2d not part of pip install package

vds_2d.py is part of the GitHub repository but not available when pip installing pypulseq.

To Reproduce
After pip install pypulseq:
from pypulseq.utils.vds_2d import vds_2d

Error message:
ModuleNotFoundError: No module named 'pypulseq.utils'

'split_gradient_at' function does not work well

Describe the bug
split_gradient_at function does not split the gradient defined as 'trap' correctly.

Expected behavior
For example, if the gradient maximum value is gy.amplitude=138888 and the rise and fall times are of 3xraster times, when you split in two equal parts the values of, ie., the fall part are array([115740.74074074, 69444.44444444, 23148.14814815]). See that the maximum values is not 138888 anymore, and that it does not end in zero either.

Indexing problem prevents reading of rf pulse in .seq file

Describe the bug
Indexing error appears when trying to load a regular .seq file (FLASH sequence, 2D Cartesian)

To Reproduce

1. Download the .seq file (https://drive.google.com/file/d/1mFMWPrEwwsw086W6KgPRy7Wh-EorLwsJ/view?usp=sharing) and put it in your current folder for running Python
2. pip install the latest PyPulseq, and run the following:

from pypulseq.Sequence.sequence import Sequence
seq = Sequence()
seq.read('FLASH.seq')
seq.get_block(1)

Expected behavior
The following error appears:

Traceback (most recent call last):
  File "<input>", line 1, in <module>
  File "C:\Users\tongg\PycharmProjects\py2jemris\venv\lib\site-packages\pypulseq\Sequence\sequence.py", line 115, in get_block
    return block.get_block(self, block_index)
  File "C:\Users\tongg\PycharmProjects\py2jemris\venv\lib\site-packages\pypulseq\Sequence\block.py", line 177, in get_block
    block.rf = self.rf_from_lib_data(self.rf_library.data[event_ind[1]])
  File "C:\Users\tongg\PycharmProjects\py2jemris\venv\lib\site-packages\pypulseq\Sequence\sequence.py", line 160, in rf_from_lib_data
    rf.dead_time = lib_data[6]
IndexError: index 6 is out of bounds for axis 0 with size 6

Desktop (please complete the following information):

• OS: Windows 10 Pro
• Version: Python 3.6.5, PyPulseq 1.2.1