Write an MRI pulse sequence as a tree of overlapping blocks. Compile it to pypulseq.
A pulseq sequence is a flat list of blocks, and one block holds at most one RF, one ADC and one gradient per axis. Hardware needs that. You do not: the moment two things you think of separately have to happen together — a slice rephaser and a phase-encode blip, a diffusion lobe straddling its refocusing pulse — they collide in one block. So you cut waveforms at boundaries by hand and keep the pieces in step yourself.
sc.LogicBlock is a tree instead:
- Anything may overlap anything. Write what you mean, when you mean it.
sc.compilefinds the legal pulseq blocks — splitting, summing, checking them against the amplifier. - Blocks nest, to any depth. A block holds events or other blocks, each timed against its parent. Scan, shot, repetition, readout: one kind of object all the way down.
- So a component is written once. Nothing has to know its neighbours or where a boundary will
fall — which is what makes
sc.Modulepossible: parameters in, oneLogicBlockout.
Here is one repetition of a spoiled 2D GRE, coloured and boxed by the LogicBlock that wrote each
waveform rather than by event type:
Each box is a module. Excitation covers two lanes, because a slice-selective pulse is one idea.
Three boxes overlap on three axes — that is the overlap rule in use, not described. The outer box is
GRE2DTR, a module built from modules. The strip underneath is the only part nobody wrote:
sc.compile put those boundaries there, straight through the boxes above.
what you write the model the compiler the output
───────────────── ───────────────── ────────────────── ─────────────────
pypulseq events ─► sc.LogicBlock ─► block boundaries ─► pypulseq.Sequence
your own modules a tree of events and same-axis sums seq.write('x.seq')
sc.modules.* blocks, with relative amplifier limits
start times
Three sections follow, and they are the whole tool: a tree by hand, a module, a sequence of modules.
Print any of it with .describe(), draw it with sc.plot_block, measure it with sc.moments.
git clone https://github.com/bughht/SeqCraft.git && cd SeqCraft
pip install -e ".[dev,viz]"Python 3.10+, numpy, and the pinned pypulseq fork (installed automatically). Optional extras:
viz (matplotlib), systems (vendor limits), rf (SLR pulses), sim / recon (the simulation
notebooks).
Important
Set the dead times. pypulseq defaults rf_dead_time, rf_ringdown_time and adc_dead_time
to zero, which is wrong on every real scanner: the sequence compiles cleanly, validates
cleanly, and is refused or silently mangled at the console. They belong to your installation, so
no preset can supply them.
Events are made by pypulseq, exactly as you already do it. seqcraft only says when each one plays:
import math
import pypulseq as pp
import seqcraft as sc
opts = pp.Opts(max_grad=40, grad_unit='mT/m', max_slew=150, slew_unit='T/m/s',
rf_dead_time=100e-6, rf_ringdown_time=30e-6, adc_dead_time=10e-6)
dk = 1e3 / 220.0 # k-space step for a 220 mm FOV, 1/m
rf, gz, gz_reph = pp.make_sinc_pulse(flip_angle=math.radians(15), duration=1e-3,
slice_thickness=5e-3, delay=opts.rf_dead_time,
use='excitation', system=opts, return_gz=True)
gx = pp.make_trapezoid('x', flat_area=64 * dk, flat_time=3.2e-3, system=opts)
adc = pp.make_adc(num_samples=64, duration=3.2e-3, delay=gx.rise_time, system=opts)
gx_pre = pp.make_trapezoid('x', area=-gx.area / 2, duration=1e-3, system=opts)
gy = pp.make_trapezoid('y', area=8 * dk, duration=1e-3, system=opts) # phase-encode line 8
t_rf_center = rf.delay + pp.calc_rf_center(rf)[0] # TE is measured from the pulse centre
echo_in_gx = gx.rise_time + gx.flat_time / 2 # and k = 0 sits here inside the readout
tr = sc.LogicBlock('tr') # a tree, tagged for error messages
tr.add(0.0, rf, gz) # RF and slice-select: one instant
tr.add(pp.calc_duration(gz), gz_reph, gy, gx_pre) # z, y, x together — overlap is free
tr.add(t_rf_center + 5e-3 - echo_in_gx, gx, adc) # TE = 5 ms, by construction
seq = sc.compile(tr, opts) # a pypulseq.Sequence — 4 blocks, 7.25 ms
seq.write('tr.seq') # pypulseq's own writeradd is the only method you need, in two shapes. add(t, *items) — above — puts items t seconds
into the block and returns it, so calls chain. add(rows) takes the whole schedule as [time, *items] rows, which is what you want once the times are computed rather than typed:
>>> table = sc.LogicBlock('tr').add([
... [0.0, rf, gz],
... [pp.calc_duration(gz), gz_reph, gy, gx_pre],
... [t_rf_center + 5e-3 - echo_in_gx, gx, adc],
... ])
>>> [n.start for n in table] == [n.start for n in tr] # the same tree, built two ways
TrueRows are never sorted, so both spellings give the same nodes. And tr.duration is measured from
the children rather than declared: a block cannot claim a length it does not play.
Ask the tree what it looks like:
>>> print(tr.describe())
tr 7.25 ms
+0.0 us rf
+0.0 us trap z
+1260.0 us trap z
+1260.0 us trap y
+1260.0 us trap x
+4010.0 us trap x
+4010.0 us adcSeven events, three add calls, no blocks anywhere. The slice rephaser, the phase blip and the
readout prephaser all start at +1260 µs; the RF shares an instant with its slice-select; the ADC
runs inside the readout gradient. In pulseq, working out which of those may share a block — and
where to cut each waveform — is your job. Here you named times, and sc.compile turned seven events
into four legal blocks. sc.plot_block(tr, opts) draws it.
Section 1 covers one phase-encode line. For the other sixty-three — and to reuse the whole thing
inside a bigger sequence — make it a component: subclass sc.Module, design in __init__, assemble
in build, return one LogicBlock.
Don't write those events twice. Three of the pieces already ship as modules — Excitation,
PhaseEncode and CartesianLine — each carrying arithmetic you would otherwise have to get right
twice: the rephaser after a selective pulse, one blip designed once and scaled per line, a prephaser
that exactly cancels the readout's ramp. Composing them is shorter than the raw events were:
class GRETR(sc.Module):
"""One repetition of a 2D gradient echo, composed from the shipped leaf modules."""
def __init__(self, *, opts, fov_mm=220.0, matrix=64, thickness_mm=5.0,
flip_deg=15.0, te_s=5e-3, bandwidth_hz_px=312.5, tag=None):
super().__init__(opts=opts, tag=tag)
self.exc = sc.modules.Excitation(opts=opts, flip_deg=flip_deg,
thickness_mm=thickness_mm, duration_s=1e-3)
# The blip and the readout prephaser play at the same instant, so the shorter is stretched
# to match: every leaf reports its own minimum and accepts an override.
blip = sc.modules.PhaseEncode(opts=opts, fov_mm=fov_mm, matrix=matrix, axis='y')
read = sc.modules.CartesianLine(opts=opts, fov_mm=fov_mm, matrix=matrix,
bandwidth_hz_px=bandwidth_hz_px)
winder_s = max(blip.min_duration_s, read.prephaser_duration_s)
self.pe = sc.modules.PhaseEncode(opts=opts, fov_mm=fov_mm, matrix=matrix, axis='y',
duration_s=winder_s)
self.ro = sc.modules.CartesianLine(opts=opts, fov_mm=fov_mm, matrix=matrix,
bandwidth_hz_px=bandwidth_hz_px,
prephaser_duration_s=winder_s)
# Where the readout has to start for k = 0 to land at TE. Quantised once: a computed start
# time must sit on the gradient raster, and the compiler raises with this exact fix if not.
self._read_start_s = sc.Raster(opts.grad_raster_time).ceil(
self.exc.time_to_center() + te_s - self.ro.time_to_echo())
def time_to_echo(self) -> float:
"""Seconds from the start of this block to k = 0 — two module times, added."""
return self._read_start_s + self.ro.time_to_echo()
def build(self, *, line: int = 0) -> sc.LogicBlock:
return sc.LogicBlock().add([
[0.0, self.exc()],
[self._read_start_s, self.pe(line=line), self.ro()], # blip on y, readout on x
])Two add calls, and no gradient areas, ramp times or dwell arithmetic anywhere: the leaves own that.
Calling the module runs build and tags the block with the class name. What comes back is a
tree, because each leaf contributed a block of its own:
>>> gre_tr = GRETR(opts=opts)
>>> gre_tr(line=40)
LogicBlock(GRETR, 3 nodes, 7.23 ms)
>>> print(gre_tr(line=40).describe())
GRETR 7.23 ms
+0.0 us Excitation 1.80 ms
+0.0 us rf
+0.0 us trap z
+1260.0 us trap z
+3670.0 us PhaseEncode 0.32 ms
+0.0 us trap y
+3670.0 us CartesianLine 3.56 ms
+0.0 us trap x
+320.0 us trap x
+320.0 us adc
>>> sc.moments(gre_tr(line=32), order=0)['y'] # the centre line needs no blip
0.0Three children where section 1 had seven events, and the blip and the prephaser still land on one
instant — +3670 µs, one on y and one on x. Ask for te_s=3e-3 and the pair slides back to
+1670 µs, beside the slice rephaser still playing on z: three axes at once, five pulseq blocks
becoming three, and nothing in the module changed to allow it.
Three conventions make a module reusable, and all three are above:
__init__designs,buildassembles. Waveforms and timings are computed once; sixty-four lines are sixty-four cheap calls.- Calls are pure. A call must not mutate the module or its events —
PhaseEncodederives a scaled copy of its blip rather than rescaling one (pp.scale_gradplussc.events.derive). Mutating in a per-line method still compiles, and makes line 64 differ from line 1 in a way no check can see. - Declare no duration, declare no position. The block measures itself, and a module that must
say when something happens inside it exposes a time instead:
time_to_echo()is the one question the tree cannot answer, because a tree knows when its events play but not which instant among them is the echo.
Check any module in isolation with sc.compile(sc.LogicBlock('probe').add(0.0, gre_tr(line=40)), opts): a component that only works when something else happens to be beside it is not reusable.
sc.modules.GRE2DTR is this same composition with spoiler gradients, a phase-encode rewinder and an
RF-phase argument added — writing one out like this is how that one was written.
A block may hold blocks, which may hold blocks. So the same add builds every level: an inversion
and the train that follows it make a shot, and a handful of shots at different inversion times
make a T1-mapping scan.
gre_tr = GRETR(opts=opts) # from section 2
inv = sc.modules.IRPrep(opts=opts, thickness_mm=None, spoil_voxel_mm=5.0) # ships with seqcraft
tr_s = 12e-3
raster = sc.Raster(opts.grad_raster_time)
lines = sorted(range(64), key=lambda k: (abs(k - 32), k)) # centric: k = 0 acquired first
def shot(ti_s, seg=lines):
"""One inversion and the train recovering into it — itself just a LogicBlock."""
# TI runs from the inversion's effective centre to the acquisition of k = 0. Both ends are
# module questions; the subtraction is the whole layout.
t_train = raster.ceil(inv.time_to_center() + ti_s - gre_tr.time_to_echo())
rows = [[0.0, inv()]] # the inversion, then the train
rows += [[t_train + i * tr_s, gre_tr(line=k)] for i, k in enumerate(seg)]
return sc.LogicBlock('shot').add(rows)
tis = (100e-3, 300e-3, 700e-3, 1500e-3)
scan = sc.LogicBlock('ir_t1').add([[i * 4.0, shot(ti)] for i, ti in enumerate(tis)]) # 4 s apart
seq = sc.compile(scan, opts, name='ir_t1') # 1547 blocks, 14.263 s
seq.write('ir_t1.seq')inv.time_to_center() is 5.101 ms into its own block — a 10 ms hyperbolic secant inverts at its
centre, not at its start, and referencing TI to the block start would be a 5 ms error in the one
quantity the sequence exists to control. Neither that nor gre_tr.time_to_echo() is measurable from
the tree, which is exactly the division of labour: modules know the physics, the tree knows the
times, the compiler knows pulseq.
Five levels, and every one of them is the same kind of object. Here is one shot cut down to two lines so it fits on the page — the real one is the same shape, 64 repetitions wide:
>>> print(sc.LogicBlock('ir_t1').add(0.0, shot(300e-3, lines[:2])).describe())
ir_t1 318.70 ms
+0.0 us shot 318.70 ms
+0.0 us IRPrep 11.34 ms
+0.0 us rf
+10130.0 us spoiler 1.21 ms
+0.0 us trap z
+299470.0 us GRETR 7.23 ms
+0.0 us Excitation 1.80 ms
+0.0 us rf
+0.0 us trap z
+1260.0 us trap z
+3670.0 us PhaseEncode 0.32 ms
+0.0 us trap y
+3670.0 us CartesianLine 3.56 ms
+0.0 us trap x
+320.0 us trap x
+320.0 us adc
+311470.0 us GRETR 7.23 ms
+0.0 us Excitation 1.80 ms
+0.0 us rf
+0.0 us trap z
+1260.0 us trap z
+3670.0 us PhaseEncode 0.32 ms
+0.0 us trap y
+3670.0 us CartesianLine 3.56 ms
+0.0 us trap x
+320.0 us trap x
+320.0 us adctrap x sits inside CartesianLine, inside GRETR, inside shot, inside the scan — and every
offset is read against its own parent, never against the scan. That is what makes a subtree
portable: the shot does not know it is the second one, so moving it moves everything it contains.
>>> len(scan), sum(1 for _ in sc.flatten(scan)) # four children; 1800 leaf events
(4, 1800)
>>> scan.nodes[2].start += 20e-3 # 450 events later, from one number
>>> sorted({path for _, _, path in sc.flatten(scan)})
[('ir_t1', 'shot', 'GRETR', 'CartesianLine'), ('ir_t1', 'shot', 'GRETR', 'Excitation'), ('ir_t1', 'shot', 'GRETR', 'PhaseEncode'), ('ir_t1', 'shot', 'IRPrep'), ('ir_t1', 'shot', 'IRPrep', 'spoiler')]Those paths are provenance, and nobody wrote them: they are the tags on the way down, and every
warning and error message names one. The rest is plain Python — scan.nodes is a list, sc.flatten
walks it, lb.copy() gives you a variant to retime, and one block added at sixty-four times is
shared rather than copied.
Nothing in here coordinates. inv does not know a train follows it, gre_tr does not know what
preceded it, shot does not know it is one of four, and none of them knows where pulseq's block
boundaries will fall.
sc.modules has nine building blocks, each extracted from a working sequence rather than designed:
Excitation, Refocusing, PhaseEncode, CartesianLine, EPI2D, spoiler, IRPrep, GRE2DTR
and GRE2D. The last two are a whole repetition and a whole scan, so the GRE that section 2
composed also comes ready-made:
gre = sc.modules.GRE2D(opts=opts, fov_mm=220.0, matrix=(64, 64), thickness_mm=5.0)
seq = sc.compile(gre(lines=range(64)), opts, name='gre_2d') # 256 blocks, TE 4.62 msGRE2D takes the list of phase-encode lines rather than an acceleration factor: which lines to
acquire is a sequence-programming choice, and it stays yours. Segmenting the train across several
inversions — an MPRAGE — is the same tree one level deeper, and
examples/mprage_2d/ builds it.
Refocusing is the one that shows what a module is for. A refocusing pulse
conjugates k, so between consecutive refocusing centres every axis's gradient area before the echo
has to equal its area after it — measured to the RF's effective centre, not to the middle of the
block. Getting that wrong leaves a residual that alternates sign echo to echo and reads as a
hardware fault, and both pulseq reference implementations avoid it only by setting their transmit
dead time and ringdown to the same number:
refoc = sc.modules.Refocusing(opts=opts, thickness_mm=6.25, crush_voxel_mm=5.0)
assert refoc.area_to_center_per_m == refoc.area_from_center_per_m # 600.000000 both
assert refoc.time_to_center() == refoc().duration / 2 # exactly, to 0 nsA spin echo and a sixteen-echo turbo spin echo are then the same composition with a longer list,
and examples/fse_2d/ writes both.
EPI2D is the newest, and the one where the arithmetic is the module. It is the whole of k-space in
one shot — alternating readout lobes, blips on the zero crossings, one ADC per echo — and an EPI
ghosts because those lobes alternate. Two rules follow. The sampling window is exactly centred
in its lobe, because only then do the two polarities sample one k grid; write the lobe the way a
single line is correctly written and they sample two grids half a gradient raster apart, which is
an N/2 ghost the sequence made itself. And the blip sits on the readout's zero crossing with
sc.barrier() pinned there, because otherwise the compiler splits every readout lobe in the train
and the only report of it is a merge warning:
epi = sc.modules.EPI2D(opts=opts, fov_mm=220.0, matrix=(128, 128), dwell_s=2.5e-6)
assert 2 * epi.guard_s + epi.num_samples * epi.dwell_s == epi.echo_spacing_s # exactly
assert epi.k_read_per_m[epi.echo_sample(0)] == 0.0 # k = 0 is a sample, not a momentParallel imaging then needs nothing added: blip_lines=R is the acceleration and lines is the
table. The calibration band is more build calls too — but with a table of length one, which
is a Cartesian gradient echo on the EPI's own readout lobe rather than another EPI shot.
Segmentation needs exactly one thing: phase_deg, the receiver phase, because a spoiled EPI
needs its receiver locked to a transmitter that is advancing its carrier. Four shots that differ
from one another for any reason put a replica of the object at Ny/4, and it reads as an
under-sampling artefact.
CartesianLine reads its line more than once when asked, which is a multi-echo gradient echo —
echoes and polarity, where 'monopolar' plays every lobe the same sign with a fly-back between
them and 'bipolar' alternates the sign and pays instead. polarity has no default above one
echo, because the two produce files that differ in the dwell, the echo times, the k ordering of
every second echo and the block count, and look identical in a protocol printout:
mono = sc.modules.CartesianLine(opts=opts, fov_mm=220.0, matrix=128, bandwidth_hz_px=500.0,
echoes=8, polarity='monopolar')
assert len(mono.te_s) == 8 # k = 0, once per echo
assert abs(mono.flyback_area_per_m + float(mono.gx.area)) < 1e-9 # minus the WHOLE lobe -- and
assert abs(mono.flyback_area_per_m) > 1.9 * mono.area_to_echo_per_m # not the pre-echo half of itte_s is not echo_spacing_s and that is the trap worth knowing: k = 0 is a sample, so a
bipolar train's echo times alternate about its period by one dwell either way — and a two-point
field map divided by the period instead is 0.75 % low in every voxel, which is the kind of
wrong no image inspection finds. GRE2DTR and GRE2D forward the same three arguments and gain
no arithmetic, so there is no MEGRE2D: a multi-echo GRE is GRE2D(..., echoes=8, polarity='monopolar'), and examples/megre_2d/ is the one example
directory that defines no class of its own.
sc.compile(tree, opts) returns a pypulseq.Sequence and nothing else. If the tree cannot become a
legal sequence it raises; if it had to change a waveform to make it legal it warns. There is
no report object to unpack, so there is nothing to forget to read.
- Block boundaries and overlap. Different axes at one instant become one block. Two gradients on the same axis are summed, with a warning naming both sources.
- Limits measured on the compiled waveform, which is the only place the truth is: two individually legal gradients can sum to an illegal one, and no component can see that alone.
- RF and ADC conflicts, including when only their dead times touch — which pypulseq would otherwise reject 40 000 blocks later.
- One informational warning you should expect. Three axes ramping together exceed the vector-norm slew bound that per-axis limits imply, routinely and legally, so seqcraft reports it rather than raising.
Errors name the offending number, where it happened, the tag path it came from, and what to change:
HardwareLimitError: slew 189% of the 150 T/m/s limit on axis x.
from : probe.a, probe.b
at : 0.000 ms (block 0)
reached: 284.0 T/m/s
fix
lengthen the lobe, or lower the readout bandwidth
or design that part against sc.opts.derate(opts, slew=0.52)
sc.moments, sc.kspace, sc.sample and sc.pns measure a tree directly, before any file is
written.
Notebooks, each one a sequence that works rather than a feature tour:
examples/01_getting_started.ipynb— blocks,Optsandcompile, the overlap rules and the escape hatches. Uses no modules.examples/gre_2d/— a spoiled 2D GRE three ways, then simulated and reconstructed. Six of the nine shipped modules came out of it.examples/mprage_2d/— segmented and inversion-prepared, with the null point checked in simulation.examples/mp2rage_2d/— two trains, theSETlabel that separates them, and the ratio that cancels the receive field.examples/se_2d/— a spin echo, the area balance that makes it one, and the measurement that says the echo is a T2 echo rather than a T2* one.Refocusingcame out of it.examples/fse_2d/— the same composition at sixteen echoes and then seventy-two: 4.3 minutes becomes 18 seconds, and what that costs is measured — including a ghost that every arithmetic check passes through.examples/gre_epi_2d/— the whole of k-space in one shot, and the two rules that stop it ghosting on its own. Both are built the wrong way too, because both failures compile.EPI2Dcame out of it.examples/se_epi_2d/— the same readout after a refocusing pulse, and the measurement most readers expect to come out the other way: a spin echo does not fix EPI distortion.
All seven simulation notebooks share one phantom — examples/phantom.py —
including the off-resonance map the EPI examples distort against, which is the phantom's own.
Documentation: api_reference.md (every public name, executed by CI),
architecture.md (the layering, and what is deliberately absent),
compiler.md (how boundaries are chosen, what every warning means),
writing_a_module.md (the Module contract in full).
Tests: pytest tests --doctest-modules src/seqcraft.
When not to reach for it. For one sequence, once, raw pypulseq is a fine tool. seqcraft pays for itself on a family of sequences, gradients that must overlap without hand-splitting blocks, or files that have to be reproducible six months later — and when it does not fit, what you are holding is already the pypulseq object.
Vendor data stays out of this repository. Siemens .asc descriptors carry proprietary
coefficients, so sc.hardware.load_hardware() reads them only through $SEQCRAFT_ASC_DIR.
sc.hardware.synthetic_hardware() is a vendor-free stand-in for PNS checks — not a real scanner,
and never to be used to clear a sequence for human scanning.
MIT.