Skip to content

Conversation

@Stockless
Copy link
Collaborator

This PR addresses issue #555 by adding support for reading Pulseq files in version 1.5.0.

@Stockless Stockless self-assigned this Sep 17, 2025
@Stockless Stockless marked this pull request as draft September 17, 2025 15:40
Copy link
Member

@cncastillo cncastillo left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please use three correct logic so pulseq 1.4 still works. No need to pass the version into the get_X functions, you can just check the input length of the library.

@codecov
Copy link

codecov bot commented Sep 23, 2025

Codecov Report

❌ Patch coverage is 85.98901% with 51 lines in your changes missing coverage. Please review.
✅ Project coverage is 90.89%. Comparing base (523c275) to head (2efa648).

Files with missing lines Patch % Lines
KomaMRIFiles/src/Sequence/Pulseq.jl 91.22% 25 Missing ⚠️
KomaMRIBase/src/datatypes/sequence/RF.jl 70.83% 7 Missing ⚠️
...ase/src/datatypes/sequence/extensions/AdcLabels.jl 14.28% 6 Missing ⚠️
...Base/src/datatypes/sequence/extensions/LabelInc.jl 25.00% 3 Missing ⚠️
...Base/src/datatypes/sequence/extensions/LabelSet.jl 25.00% 3 Missing ⚠️
...IBase/src/datatypes/sequence/extensions/Trigger.jl 25.00% 3 Missing ⚠️
KomaMRIBase/src/datatypes/Sequence.jl 88.23% 2 Missing ⚠️
KomaMRIBase/src/datatypes/sequence/EXT.jl 0.00% 1 Missing ⚠️
KomaMRIPlots/src/ui/DisplayFunctions.jl 90.90% 1 Missing ⚠️
Additional details and impacted files

Impacted file tree graph

@@            Coverage Diff             @@
##           master     #614      +/-   ##
==========================================
- Coverage   91.22%   90.89%   -0.33%     
==========================================
  Files          57       61       +4     
  Lines        3248     3351     +103     
==========================================
+ Hits         2963     3046      +83     
- Misses        285      305      +20     
Flag Coverage Δ
base 90.86% <63.23%> (-0.63%) ⬇️
core 89.73% <ø> (ø)
files 93.90% <91.22%> (-1.08%) ⬇️
komamri 88.13% <ø> (ø)
plots 90.97% <90.90%> (+0.06%) ⬆️

Flags with carried forward coverage won't be shown. Click here to find out more.

Files with missing lines Coverage Δ
KomaMRIBase/src/KomaMRIBase.jl 100.00% <ø> (ø)
KomaMRIBase/src/datatypes/sequence/Grad.jl 94.54% <100.00%> (+5.15%) ⬆️
KomaMRIBase/src/timing/TimeStepCalculation.jl 100.00% <100.00%> (ø)
KomaMRIFiles/src/KomaMRIFiles.jl 100.00% <ø> (ø)
KomaMRIBase/src/datatypes/sequence/EXT.jl 0.00% <0.00%> (-40.00%) ⬇️
KomaMRIPlots/src/ui/DisplayFunctions.jl 91.89% <90.90%> (+0.06%) ⬆️
KomaMRIBase/src/datatypes/Sequence.jl 93.14% <88.23%> (+0.95%) ⬆️
...Base/src/datatypes/sequence/extensions/LabelInc.jl 25.00% <25.00%> (ø)
...Base/src/datatypes/sequence/extensions/LabelSet.jl 25.00% <25.00%> (ø)
...IBase/src/datatypes/sequence/extensions/Trigger.jl 25.00% <25.00%> (ø)
... and 3 more
🚀 New features to boost your workflow:
  • ❄️ Test Analytics: Detect flaky tests, report on failures, and find test suite problems.
  • 📦 JS Bundle Analysis: Save yourself from yourself by tracking and limiting bundle sizes in JS merges.

@cncastillo
Copy link
Member

cncastillo commented Oct 6, 2025

(1) Update get_events to add compat_helper copy: https://github.com/pulseq/pulseq/blob/master/matlab/%2Bmr/%40Sequence/read.m#L520. get_events cannot be hardcoded for the RF case

(2) Test with reader with (take inspiration from https://juliahealth.org/KomaMRI.jl/dev/explanation/5-seq-events/):

  • arbitrary gradients
  • time shaped gradients
  • trapezoidal gradients
  • arbitrary rfs
  • time shaped rfs
  • block rfs

@Stockless can you give use the MATLAB files you were using for the tests?

(3) Modify to consider new RF uses, to calculate k-space correctly (not guess RF use by > 90.0):

  • get_RF_types() → 1st
  • get_RF_center() → 2nd
  • get_Mk() → 3rd

If use is Undefined(), keep using the current method; otherwise, use the provided value.

If center is not defined (i.e., == 0.0), keep the current method; otherwise, use the provided value.

@gsahonero
Copy link
Member

@pvillacorta , @cncastillo mentioned you are interested in working this PR out. I want to get it done too. Let's meet to sort out the efforts?

@pvillacorta
Copy link
Collaborator

@pvillacorta , @cncastillo mentioned you are interested in working this PR out. I want to get it done too. Let's meet to sort out the efforts?

Yes! Let's join forces :)

@cncastillo
Copy link
Member

@pvillacorta @gsahonero do you guys want to plan a meeting to push this?

@pvillacorta
Copy link
Collaborator

@gsahonero and I are meeting this thursday at 4pm (spanish time). I hope we can plan how to finish this PR and the one for writing Pulseq.

@pvillacorta pvillacorta changed the title Read pulseq 1.5.0 Read/Write Pulseq 1.5.0 Nov 13, 2025
@pvillacorta pvillacorta changed the title Read/Write Pulseq 1.5.0 ReadPulseq 1.5.0 Nov 13, 2025
@pvillacorta pvillacorta changed the title ReadPulseq 1.5.0 Read Pulseq 1.5.0 Nov 13, 2025
- Improve Grad and RF constructors for improved readability and error handling.
- Fix documentation typo about `A` and `T` vectors
- Add SHA and MD5 packages for hash generation/verification
- Improve signature generation/verification
- Add missing signature-related functions
- Simplify `read_RF` and `read_ADC` in the same way as `read_Grad`
Copy link
Contributor

@github-actions github-actions bot left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

KomaMRI Benchmarks

Details
Benchmark suite Current: 2efa648 Previous: 11d3180 Ratio
MRI Lab/Bloch/CPU/2 thread(s) 334561715.5 ns 337633469 ns 0.99
MRI Lab/Bloch/CPU/4 thread(s) 278611014 ns 276282606 ns 1.01
MRI Lab/Bloch/CPU/8 thread(s) 279866990.5 ns 209852165 ns 1.33
MRI Lab/Bloch/CPU/1 thread(s) 559995522 ns 555065232 ns 1.01
MRI Lab/Bloch/GPU/CUDA 20997855 ns 21231134 ns 0.99
MRI Lab/Bloch/GPU/oneAPI 73292888 ns 77053690 ns 0.95
MRI Lab/Bloch/GPU/Metal 92637208 ns 95540333 ns 0.97
MRI Lab/Bloch/GPU/AMDGPU 23963376.5 ns 24763685 ns 0.97
Slice Selection 3D/Bloch/CPU/2 thread(s) 1571567038 ns 1592087059 ns 0.99
Slice Selection 3D/Bloch/CPU/4 thread(s) 888550775 ns 889539516 ns 1.00
Slice Selection 3D/Bloch/CPU/8 thread(s) 567568116 ns 565401528.5 ns 1.00
Slice Selection 3D/Bloch/CPU/1 thread(s) 3038525168.5 ns 3029269071 ns 1.00
Slice Selection 3D/Bloch/GPU/CUDA 31434342 ns 32639703 ns 0.96
Slice Selection 3D/Bloch/GPU/oneAPI 118202439 ns 121489006.5 ns 0.97
Slice Selection 3D/Bloch/GPU/Metal 108188250 ns 111152750 ns 0.97
Slice Selection 3D/Bloch/GPU/AMDGPU 31967231 ns 32636098 ns 0.98

This comment was automatically generated by workflow using github-action-benchmark.

- Fix contructors, auxiliary functions and tests
- Restore previous `read_Grads`, `read_RF`, and `read_ADC` implementations
- Improve signature handling
@pvillacorta
Copy link
Collaborator

Pulseq files from KomaMRIFiles/test/test_files/pulseq_read_comparison must ve reviewed, since I have not been able to retrieve their signature. For the moment, I have just added a warning when the retrieved signature and the expected one do not coincide (obviously, the ideal case would be to throw an error, but tests would not pass for pulseq_read_comparison).

@pvillacorta
Copy link
Collaborator

pvillacorta commented Dec 27, 2025

I have modified and created new MATLAB scripts into a local version of this folder in order to make similar (v1.5.1) sequences to those existing in v1.4.
Testing aspects such as labels and extensions is pending.

- Moved `AdcLabels`, `LabelInc`, `LabelSet`, and `Trigger` into separate files for better organization.
- Updated `EXT.jl` to include new extension files and added functions to retrieve extension types from symbols.
- New `get_scanf_format` function to generate format strings based on struct field types.
- Modified `read_events` and `read_extension` functions to accommodate the new structure and improve extensibility.
- Updated `read_seq` to handle new extension libraries and ensure compatibility with existing sequence files.
@pvillacorta
Copy link
Collaborator

pvillacorta commented Dec 31, 2025

The commit above tries to address #575. The main change is that it is no longer necessary to create a library for each extension type inside the reader.
Other important changes are the simplification of function read_events and the deletion of read_labels and read_extension_blocks.
Now, extensions are managed through three main dictionaries within theread_seq function: extensionInstanceLibrary, extensionTypeLibrary, and extensionSpecLibrary.
If we wanted to define a new extension, we would just need to create its corresponding file within the KomaMRIBase/src/datatypes/sequence/extensions directory (obviously, if we wanted that extension to take effect, for instance, in the simulation, more changes would need to be made).
The following v1.5.1 extensions are not implemented yet:

  • Soft Delay
  • Rotation
  • RF Shimming

Maybe leave them for another PR?
A hint about how to do this for the Soft Delay case (not implemented yet):
KomaMRIBase/src/datatypes/sequence/extensions/SoftDelay.jl:

mutable struct SoftDelay <: Extension
    num::Int
    offset::Int
    factor::Int
    hint::String
end

Let me know what you think, I hope we can merge this soon in order to also address #152 during the following days.

@pvillacorta pvillacorta marked this pull request as ready for review January 4, 2026 09:24
@pvillacorta pvillacorta linked an issue Jan 4, 2026 that may be closed by this pull request
Copy link
Member

@cncastillo cncastillo left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I jut wanted to say, that the introduction of trapLibrary and other dictionaries seem to deviate too much from the MALTAB's Pulseq reader. I would like to follow the same logic, even if inefficient so it is easier to update when new versions are released.

Also, the whole handling of fmt strings for scanf (like the function get_scanf_format) seem overly complicated. The comment I left gives a solution, but there might be an easier one like:

struct MyStructure
	x::Float64
    y::Int
end
t = MyStructure(parse.(fieldtypes(MyStructure), split("1.0 1"))...)

The extensive use of Dict's to map inputs to types seems weird to me, and all the generated arrays of type Any are not very performant.

Finally, I think there is too many functions like f(x::Type{T}) instead of f(x::T) which seems too complicated.

I haven't reviewed the whole PR, but I don't think in the current state in can be merged quickly, there are too many "unnecessary" changes. For example, there are some breaking changes on how to read Pulseq <1.4.

My main point: Let's try to stick to whatever Pulseq is doing, let's only deviate for performance or extensibility.

- Remove `get_scanf_format` and `get_scanf_args`
- Reduce the number of type Any arrays
- Put `get_EXT_type_from_symbol` implementations inside their corresponding file
@pvillacorta
Copy link
Collaborator

pvillacorta commented Jan 6, 2026

I jut wanted to say, that the introduction of trapLibrary and other dictionaries seem to deviate too much from the MALTAB's Pulseq reader. I would like to follow the same logic, even if inefficient so it is easier to update when new versions are released.

I agree that this deviates from the exact structure of the MATLAB reader. The main motivation for introducing trapLibrary was to simplify read_events and read_Grad, which were previously hard-coded:

  • RF events relied on character-level checks (e.g. last character determining RF type).
  • Gradients were distinguished via tag characters (tvs g) scattered throughout the code.

By separating trapezoidal gradients from arbitrary gradients into a dedicated library, the dispatch logic in read_events becomes simpler and more explicit. This reduces the amount of branching and character-based parsing while still preserving the underlying Pulseq concepts.

At this point, I see this as a pragmatic compromise: it deviates slightly from MATLAB’s implementation, but significantly reduces code complexity and improves maintainability. My intention was not to redesign the Pulseq logic, but to make the Julia reader easier to reason about and extend while keeping the mapping to Pulseq concepts clear.

Also, the whole handling of fmt strings for scanf (like the function get_scanf_format) seem overly complicated. The comment I left gives a solution, but there might be an easier one like:

struct MyStructure
	x::Float64
    y::Int
end
t = MyStructure(parse.(fieldtypes(MyStructure), split("1.0 1"))...)

You’re absolutely right here, and I’ve already addressed this in the latest commit.

I removed the generic auxiliary functions (get_scale, get_scanf_args) and instead implemented the parsing logic directly in each extension and inside read_events. This removes unnecessary abstraction and makes the parsing much more straightforward and readable.

The extensive use of Dict's to map inputs to types seems weird to me, and all the generated arrays of type Any are not very performant.

This is a fair concern, especially regarding performance.

Here I intentionally tried to stay close to what was already present in the reader before this PR. Dictionaries were already used extensively, and keeping them allowed me to:

  • Minimize the scope of changes.
  • Stay conceptually close to Pulseq’s EventsLibrary (even though Pulseq uses a struct rather than a Dict).

I fully agree that performance could be improved by using more concrete types (or even arrays directly, as it was stated in #224). However, given time constraints and the goal of not overengineering this PR, I preferred to keep consistency with the existing codebase rather than introduce a larger refactor.

Finally, I think there is too many functions like f(x::Type{T}) instead of f(x::T) which seems too complicated.

The only places where I dispatch on Type{T} are get_scanf_format and get_scale. This was a deliberate choice to integrate the extension handling intoread_events without introducing a separate reading pipeline. Dispatching on types allows read_events to treat extensions in the same way as other events (such as gradients or RFs).

Instantiating a dummy object before parsing its arguments felt less practical (and potentially more expensive) than dispatching on the type directly. While I agree this is not ideal in general, in this context it seemed like a reasonable trade-off to keep the overall logic simpler and more uniform.

I haven't reviewed the whole PR, but I don't think in the current state in can be merged quickly, there are too many "unnecessary" changes. For example, there are some breaking changes on how to read Pulseq <1.4.

I understand this concern, but I’m not entirely sure which specific breaking changes for Pulseq < 1.4 you’re referring to. If you can point out a concrete example or file, I would try to adjust or revert that part.

My main point: Let's try to stick to whatever Pulseq is doing, let's only deviate for performance or extensibility.

I fully agree with your main point. we should stick to Pulseq's logic and only deviate when there is a clear benefit in terms of performance or extensibility. My intention was exactly that, not to redesign the reader functions, but to make the Julia reader usable, maintainable and more extensible. Given the time constraints, I think this is a reasonable compromise, and the upcoming Pulseq writer PR will follow the same direction.
Thank you!

Copy link
Member

@cncastillo cncastillo left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, now I fully reviewed the PR, I think it is pretty good. Most of my comments are just to verify that we match what Pulseq does.

@test first(ev_koma.A) first(ev_pulseq.A) || ev_koma.t[2] ev_koma.t[1]
@test last(ev_koma.A) last(ev_pulseq.A)
pth = joinpath(@__DIR__, "test_files/pulseq/pulseq_read_comparison/")
versions = ["v1.4", "v1.5"]
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

No tests for versions lower than 1.4?

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I currently don’t have an easy way to generate test data for Pulseq < 1.4. Supporting it would require setting up older MATLAB versions and regenerating .seq and .mat files. This is doable, but I wasn’t sure whether <1.4 is still actively used (e.g., by JEMRIS). Happy to add this if it’s considered important.

@JuliaHealth JuliaHealth deleted a comment from cncastillo Jan 11, 2026
- Remove `skip_rf`
- Modify `RF` constructor: `_RF_with_center` -> `_RF_with_center_and_use`
- Adapt `get_RF_center` to what Pulseq does
- Change `rf-pulse.seq` (and .mat) to consider a non-zero delay
- "RFs" labels -> "RF_center"
- Rename all dictionary-reading functions (`read_*`->`get_*`)
- Simplify version handling in `get_Grad`, `get_RF` and `get_ADC`
- Dispatch for unknown extensions
- Display RF center in `plot_seq` (new `cents` aux function)
- New `has_hash` variable
- Dispatch depending on signature existance
- Remove `blockDurations` and `delayInd_tmp` dictionaries
- Recover 't' and 'g' gradient types, delete `trapLibrary`
@pvillacorta
Copy link
Collaborator

pvillacorta commented Jan 11, 2026

@cncastillo I have resolved all the comments that I could address by myself.
I have not resolved yet the ones in which you may decide or give additional opinions.
Thank you!

- Fix `fix_first_last_grads!` to mimic Pulseq's MATLAB implementation
- New `simplify_waveforms` function
- Fix RF center representation in `plot_seq`:
   - Remove useless `cents` function
   - Introduce `cents` functionality into `times` function
   - Add previous blocks timing information
- Better distribution of KomaFiles test files
@pvillacorta
Copy link
Collaborator

pvillacorta commented Jan 14, 2026

I made the last commit in order to (we can talk about it more in detail in a meeting if you consider):

  • Simplify "one-value" waveforms before intantiating the actual KomaMRI events (RF or Grad)
    • I used a new function: simplify_waveforms
    • This made me find a discrepancy between our implementation of fix_first_last_grads! and theirs: they update the gradLibrary dictionary with the obtained first and last values, and we were not doing that. I have made the necessary changes to make them do the same. Again, we can better explore this if you need.
  • Fix bugs in plot_seq. RF were not displaying correctly when more than one RF pulses were active in the sequence.

@cncastillo
Copy link
Member

I think we might need to organize a meeting there's too many changes.

I don't understand why the fix_first_last_grads! needs to be applied to the dictionary. Is there a difference in the obtained first and last points? Before they didn't always match Pulseq's, now they always match?

Also I noticed that some of the tests use inside(x)=x[2:end-1] to compare so we don't consider the first and last element of the vector A, but now, first and last are in the Grad data structure. There seems to be a mismatch.

The functions times and ampls are to get samples for simulation and plotting. I don't really understand the addition of the centers, it doesn't seem to be the right place.

Can you explain the bug you encountered with plot_seq?

- finish `fix_first_last_grads!`
- Remove debugging comments
- Recover previous `get_Grad` function
- Simplify tests (not use `inside` function)
- Simplify plot centers in `plot_seq`
- Add warning when a gradient time shape starts at a non-zero value.
- Update info message regarding v1.5 support
@pvillacorta
Copy link
Collaborator

pvillacorta commented Jan 16, 2026

This comment summarizes the points discussed in yesterday’s meeting with @cncastillo and their current status after the latest commit (9c84075). Some items required code changes (included in the commit), while others were already resolved and have now been verified:

  • Simplify center plotting in plot_seq
    Addressed in the latest commit.
  • format_helper
    Verified as discussed. See this comment for confirmation and resolve it if you agree.
  • Test time_shape_id = -1
    The spiral.seq file (v1.5) already contains gradients with time_shape_id = -1.
    See the file that generates this sequence.
  • Emit a warning if the gradient time shape does not start at 0
    Implemented in the latest commit.
  • Review how Pulseq detects the version in get_block
    Pulseq does it as follows:
        % for versions prior to 1.4.0 blockDurations have not been initialized
        obj.blockDurations=zeros(1,length(obj.blockEvents));
        % scan trhough blocks and calculate durations
        for iB = 1:length(obj.blockEvents)
            b=obj.getBlock(iB);
            if delayInd_tmp(iB) > 0
                b.delay.type = 'delay';
                b.delay.delay = tmp_delayLibrary.data(delayInd_tmp(iB)).array;
            end
            obj.blockDurations(iB)=mr.calcDuration(b);
        end
    They do not include the pulseq version as a argument for getBlock, but check if a delay is present (v<1.4) and add it to the tmp_delayLibrary dictionary.
    Then, inside getBlock, they do it as follows:
                block.blockDuration=obj.blockDurations(index);
    %             % now that delays in v1.4 and later Pulseq revisions are not
    %             % stored, we need to see whether we need a delay to explain the
    %             % current block duration
    %             if (mr.calcDuration(block)~=obj.blockDurations(index))
    %                 tmpDelay.type = 'delay';
    %                 tmpDelay.delay = obj.blockDurations(index);
    %                 block.delay = tmpDelay;
    %             end
    I think that this is not intuitive, and introducing pulseq version as an get_block input argument is totally worth it.

Additionally, I have changed the information message that appears when loading a v1.5 sequence:

┌ Info: Pulseq 1.5.1 is supported, but Soft Delay, Rotation, and RF Shimming extensions are not yet included
└ (see https://github.com/JuliaHealth/KomaMRI.jl/issues/714)

The new issue #714 accounts for Pulseq v1.5 extension support.

@pvillacorta
Copy link
Collaborator

Buildkite is failing due to some problems with CUDA:
Screenshot 2026-01-16 at 13 38 02
Not related to this PR

Copy link
Member

@cncastillo cncastillo left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I don't see any test failing!

I think that this is not intuitive, and introducing pulseq version as an get_block input argument is totally worth it.

Ok!

Additionally, I have changed the information message that appears when loading a v1.5 sequence:

I would only give a warning if the user uses any of the non-implemented extensions. As it will become annoying for people to now use them. Or at least, maxlog=1.

A=reduce(vcat, [usrf(block.rf.A); Inf] for block in seq_samples),
t=reduce(vcat, [usrf(block.rf.t); Inf] for block in seq_samples),
ct=center_times,
cA=abs.(KomaMRIBase.get_rfs(seq, center_times)[1])
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Another important parameter to show at the RF center is the phase.

Comment on lines +366 to +377
p[2j + 5] = scatter_fun(;
x=rf.ct * 1e3,
y=rf.cA * 1e6 * frf,
name="RF_center",
hovertemplate="RF center: %{x:.4f} ms<br>Amplitude: %{y:.2f} μT<extra></extra>",
xaxis=xaxis,
yaxis=yaxis,
legendgroup="RF_center",
showlegend=showlegend,
mode="markers",
marker=attr(; color="#FF0000", symbol="x"),
)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

For very long sequences (> 1000 RFs) can you check if this is too slow?, maybe the center display can be turned off by default.

@cncastillo cncastillo changed the title Read Pulseq 1.5.0 Read Pulseq 1.5.1 Jan 17, 2026
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

Make extensions extensible

6 participants