Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Load pulseq sequences faster and allow comparison #238

Draft
wants to merge 1 commit into
base: master
Choose a base branch
from

Conversation

gabuzi
Copy link
Contributor

@gabuzi gabuzi commented Dec 3, 2023

This addresses #224. Currently a draft, because it hasn't been extensively tested on other pulseq versions.

The below snippet was used for comparison and also contains some further explanations.
An example sequence file is here bSSFP_FA30deg_TE10ms_TR20ms_2D_(69x64)_pulseq.seq.zip. This small sequence loads quite fast, but for e.g. 3D sequences with many excitations, the performance issues and benefits from this diff are greatly increased.

Proper unit tests should of course be required. What I did here was to implement == operator for
sequences (and constituent structs) to compare structurally, rather than by ref, and base the
sequence equality on that. I assert that results of the existing implementation and this one are equal.
There may be an existing julia package to do that equality through a macro (StructEquality.jl),
rather than by hand, but I wanted to keep changes light.

There seemed to be a rather big refactor going on with KomaMRIFiles subpackage that threw my
environment off. As a result, this is based on an earlier version.

using KomaMRI


seqfile = "bSSFP_FA30deg_TE10ms_TR20ms_2D_(69x64)_pulseq.seq"
# force compilation
read_seq(seqfile)
time_read_seq(f) = @time read_seq(f)
println("baseline")
seq = time_read_seq(seqfile)

# profiler says most time is spent in the sequence inner constructor.
# there's really lots of copying going on.
# the read function goes through pulseq blocks and creates a sequence object for the first block
# then, it creates a new seq object for the next block and adds them together with +
# This is implemented by creating yet another sequence object with the concatenation of first and second sequence'
# elements. Then, the sequence constructor itself contains quite some logic that may involve instantiation of
# sequence sub elements.

# One obvious choice is to implement some sort of cat operator for sequences to concatenate
# multiple blocks at once and then have the sequence constructor be called only once on the concatenated args.
# This helped, but then scanf was a new bottleneck (formatter created in the loop, etc.)
# Dictionaries, indexed via id are another source of redirection, which, according to the pulseq
# specification isn't necessary.
# I thus decided to just read the blocks section directly into one big int array
seq_blocks_no_id_sort = read_seq_via_blocks_as_int_array(seqfile)
time_read_seq_via_blocks_as_int_array(f) = @time read_seq_via_blocks_as_int_array(f)
println("optimized without scanf, without blocks id sorting and dicts")
time_read_seq_via_blocks_as_int_array(seqfile)
plot_seq(seq_blocks_no_id_sort)
@assert seq == seq_blocks_no_id_sort "different sequence results!"

# now, a large part of the time (can be up to 40%) is inferring the Nx Ny Nz that we could
# provide directly in the pulseq file...

@beorostica
Copy link
Contributor

The funtions read_blocks, read_events, read_shapes outputs some dictionaries. Maybe it is a good idea not use dictionaries for increasing reading speed, or at least define them as Dict{String, Vector{Float64}}

@cncastillo
Copy link
Member

I am not closing this yet as I would want to compare the current optimized reader with this. @beorostica could you do a performance comparison using gre.seq?

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.

3 participants