ArbEPI-python: Python/pypulseq port of ArbEPI

Python port of ArbEPI (MATLAB/Pulseq), using pypulseq as the Pulseq layer. Generates fast, vendor-agnostic 3D-EPI sequences from arbitrary 2D (ky, kz) sampling masks in the phase-encode-partition plane, exports them to GE (.pge), and preprocesses and reconstructs the resulting scanner data (preprocess/, recon/).

Pipeline

Each of the four sequences keeps its own lane from generation through the scanner. The noise scan sets the whitening matrix W, the EPIcal scan calibrates the readout delay and odd/even phase, the ArbEPI time series is gridded and scattered into k-space, and the shared deGRE supplies coil compression, sensitivity maps, B0 and R2*. The sampling schedule never passes through the scanner: it travels in scan_info.mat, which preprocessing needs to put each echo back at its (ky, kz) location.

Stage	Environment	Command	Writes
Sequences	main (uv sync)	uv run python main.py --plot	output/*.seq, output/scan_info.mat, diagnostic plots
GE export	main	uv run python main.py --ge	output/*.pge
Preprocess	.venv-preprocessing	.venv-preprocessing/bin/python -m preprocess.batch_preprocess <datdir> <seq> ...	<datdir>/recon/<seq>_preprocessed.h5
RSS	.venv-recon (uv sync, see recon/README.md)	UV_PROJECT_ENVIRONMENT=.venv-recon uv run --extra recon python -m recon.rss <datdir> <seq>	<datdir>/recon/rss/
Iterative SENSE	.venv-recon	UV_PROJECT_ENVIRONMENT=.venv-recon uv run --extra recon python -m recon.sense <datdir> <seq> --reg mslr --B0	<datdir>/recon/sense_mslr_b0/
Motion correction (+ brain mask)	main (uv sync --extra analyze, see postprocess/README.md)	uv run python -m postprocess --recon <datdir>/recon/sense_mslr_b0/<seq>_recon.h5 --out <datdir>/recon/sense_mslr_b0/<seq>_mc_recon.h5	<seq>_mc_recon.h5 (same layout as the input) + .par + .json; add --brain-mask-from <datdir>/recon/<seq>_preprocessed.h5 to zero everything outside a deGRE-derived brain mask (brain_mask, _brain_mask.nii.gz)
Activation maps	main (uv sync --extra analyze, see analyze/README.md)	uv run python -m analyze --recon <datdir>/recon/sense_mslr_b0/<seq>_recon.h5 --tr <TR> --onsets ... --duration 20 --scales all --out <dir>	<dir>/<seq>_recon_{sum,scaleK}_{t,z,psc}.nii.gz

Before preprocessing, put the raw ScanArchives in <datdir>/scanarchives/ and copy that session's scan_info.mat to <datdir>/seqs/<seq>/ (see preprocess/README.md).

Without a scanner, simulate_fmri/ simulates the session: from scan_info.mat alone it writes the raw readouts of all four scans (a brain phantom with a susceptibility-derived field map, ramp sampling, readout delay and odd/even phase, coil noise, physiological noise and a known BOLD response to a visual-motor block task) into <datdir>/scanarchives/, which preprocess reads like ScanArchives, and a truth file that python -m analyze.score scores a reconstruction against (.venv-simulate; python -m simulate_fmri.session output/scan_info.mat <datdir>). python -m simulate_fmri.study repeats the experiment for ROC and test-retest reliability analyses (simulate_fmri/analysis.py).

The diagram is generated by docs/pipeline/make_pipeline.py. To change it, edit draw() there, run uv run python docs/pipeline/make_pipeline.py, and commit the script and both SVGs together. tests/test_pipeline_diagram.py fails if the SVGs are stale or if a .py/.jl file the diagram names has been moved or deleted.

Getting started

1. Install dependencies, managed with uv:

   uv sync --extra test
   

   Add --extra lint for ruff. Dependencies are listed in pyproject.toml and pinned in uv.lock; preprocess/ and recon/ each use their own venv (see Architecture); recon/'s is synced from the same lockfile with UV_PROJECT_ENVIRONMENT=.venv-recon uv sync --extra recon --extra test.

2. Edit params.py's load_params() to configure the experiment. The "USER CONFIGURATION" section near the top of that function holds the basics:

   Parameter	Meaning	Default
   scanner	Scanner hardware profile: 'GE_MR750' or 'GE_UHP' (see scanners.py)	'GE_MR750'
   res	Voxel resolution [x, y, z], m	[2.4, 2.4, 2.4] mm
   N	Acquisition matrix size [Nx, Ny, Nz] (fov = N * res)	[90, 90, 60]
   TE	Nominal echo time, s	30 ms
   volume_tr	Time to acquire one full 3D volume, s. Set to the fastest the default protocol allows (per-shot min TR x Nshots, rounded up to the ms); re-measure after changing TE, ETL, R, resolution or excitation	0.506 s
   duration	Total scan duration across all frames, s	60 s
   T1	Tissue T1, s (sets the Ernst-angle flip angle)	1.3 s
   n_frames_discard	Frames to discard at the start of the scan (steady-state warm-up). Sets discard_duration = n_frames_discard * volume_tr (s), the warm-up EPIcal plays as dummy (ADC-off) shots; the frames are kept in the data, preprocess/ records the count	0
   ETL	Echo train length (echoes per shot)	54
   custom_mask_path	Path to your own (ky, kz[, t]) sampling mask .mat file, or None to use sampling_method below (see Using custom ky-kz-t sampling masks)	None
   R	Acceleration factor on the (ky, kz) sampling pattern (ignored if custom_mask_path is set)	10
   sampling_method	ky-kz(-t) sampling pattern: 'pd', 'caipi', 'ticaipi', or 'rand' (ignored if custom_mask_path is set)	'pd'
   seed	Sampling-mask RNG seed: an int for a reproducible mask, or None for a fresh one each run (ignored if custom_mask_path is set)	0
   excitation	Fat suppression: 'fatsat' (fat-sat pulse + crusher before the slab-selective sinc excitation) or 'water' (binomial water excitation; see Differences vs. MATLAB original). 'fatsat' needs volume_tr >= 0.586 s at the default protocol	'water'
   epi_trajectory	Echo-train ordering within each shot: 'laminar' or 'radial' (see Trajectory computation below, and demo.ipynb for a side-by-side comparison) — always required, even with a custom mask	'radial'
   Ncoils	Number of receive coil channels (for the noise prescan)	32

   Everything below that section in load_params() is pre-tuned for this sequence's hardware/PNS/timing constraints (see CLAUDE.md) and typically doesn't need to change.

3. Run main.py to generate all four sequences:

   uv run python main.py
   

   Or step by step:

   from params import load_params
   from sample.gen_sampling_masks import resolve_omegas
   from sequences.ArbEPI import generate_arbepi
   from sequences.EPIcal import generate_epical
   from sequences.deGRE import generate_degre
   from sequences.noise import generate_noise
   
   params = load_params()
   # resolve_omegas returns params.custom_omegas when params.custom_mask_path
   # points at your own mask (see "Using custom ky-kz-t sampling masks" below),
   # else falls back to gen_sampling_masks.
   omegas = resolve_omegas(params)
   generate_arbepi(omegas, params)   # writes output/ArbEPI.seq, output/scan_info.mat
   generate_epical(params)           # writes output/EPIcal.seq
   generate_noise(params)            # writes output/noise.seq
   generate_degre(params)            # writes output/deGRE.seq (dual-echo, for coil sensitivity maps + B0 field map)

   To scan the 20-echo, 2.4 mm monopolar multi-echo GRE instead of the dual-echo deGRE (B0, R2* and coil maps in 267 s; sequences/meGRE.py), run uv run python main.py --gre megre or call generate_megre(params) in place of generate_degre(params); it writes output/meGRE.seq and patches scan_info.mat the same way. Preprocessing can use just the first K of its echoes (--gre-n-echoes K) to find how many are needed.

   generate_epical and generate_noise must run after generate_arbepi — they load output/scan_info.mat. generate_degre can run standalone, but run it after generate_arbepi too: it patches its realized echo-time pair (TE_degre) into scan_info.mat when that file exists. All outputs go to params.output_dir (default output/, gitignored).

   demo.ipynb runs these same steps and shows the diagnostic plots (sampling mask, PSF, radial-vs-laminar trajectory, one-TR pulse diagram, PNS) and the GE feasibility check along the way:

   uv run --with jupyter jupyter lab demo.ipynb
   

4. Add --plot to also write diagnostic plots (mask.png, psf.png, trajectory.png, one_tr.png, PNS_one_tr.png) via plot/plot_last_run.py, and/or --ge to also export each sequence to GE .pge (see GE export below):

   uv run python main.py --plot --ge
   

There is no automated end-to-end pytest suite against MATLAB for .seq generation itself (no MATLAB install was available during that initial port — see tests/ for unit tests on algorithm invariants instead, including an independent check that reads k-space back out of the assembled sequence and confirms it matches the sampling schedule). sampling_method='caipi' is deterministic (no RNG) and is the easiest configuration to sanity-check by hand. (The separate GE .pge export path below was validated against real MATLAB output, once a MATLAB install became available — see the GE export section and CLAUDE.md.)

Using custom ky-kz-t sampling masks

If you already have your own (ky, kz) or (ky, kz, t) sampling pattern — designed by hand, by a collaborator's own pipeline, or by a compressed-sensing/non-Cartesian sample-selection method this repo doesn't itself implement — you can feed it straight into the pipeline instead of one of the built-in sampling_methods ('pd', 'caipi', 'ticaipi', 'rand'). This bypasses gen_sampling_masks entirely and replaces the need to specify R, sampling_method, and seed.

1. Save your mask as a plain v5 .mat file (e.g. MATLAB's save(..., '-v5'), or scipy.io.savemat from Python) holding a 0/1 array under one variable, either:

   • (Ny, Nz) — a single static pattern, reused for every frame, or
   • (Ny, Nz, Nframes) — an already time-resolved mask, one pattern per frame/timepoint.

   Ny/Nz must match params.py's N[1]/N[2]. If you provide the 3D form, its Nframes must match what duration/volume_tr/n_frames_discard compute (round((duration + n_frames_discard * volume_tr) / volume_tr)) — adjust one side or the other so they agree.

2. In params.py's load_params(), point custom_mask_path at that file (and custom_mask_key if your variable isn't named 'samp'):

   custom_mask_path = 'my_custom_mask.mat'
   custom_mask_key = 'samp'

   epi_trajectory ('laminar' or 'radial') is still required either way: mask2epi_{laminar,radial} still partitions whatever mask you provide into Nshots EPI trajectories of length ETL, regardless of where the mask came from.

3. Run main.py exactly as usual (uv run python main.py, optionally with --plot/--ge). load_params() loads and validates the mask up front — every frame must sample the same number of (ky, kz) locations, and that count must divide evenly by ETL, since Nshots = samples_per_frame / ETL (see lib/mask2epi.py's Nshots * ETL == n_samples assertion) — and derives Nshots and an effective R (Ny*Nz / samples_per_frame, for display/scan_info.mat bookkeeping only, no longer a design input) from it. An invalid mask raises a ValueError explaining exactly what to fix before any sequence generation starts.

4. Re-verify PNS and the achieved TE — a custom mask changes the numbers this repo's own defaults were tuned against. The default blip_slew/ro_slew_rise/ro_slew_fall in params.py were tuned to sit close to GE's 80% normal-mode PNS limit against this repo's own Poisson-disc masks (see CLAUDE.md's "PNS finding history"); a custom mask can produce very different consecutive-sample ky/kz steps (lib/mask2epi.py's max_blip_steps), which directly changes blip amplitude and the resulting PNS. Run uv run python main.py --ge (or just check_ge_feasibility, see GE export below) and, if PNS comes back over 80%, lower ro_slew_rise/ro_slew_fall/blip_slew (the 2026-09-30 sweep found the default protocol readout-limited, ~1 percentage point under the line: e.g. rise/fall 150/170 measured ~1.5 points lower at an 8 µs longer echo spacing). A custom mask also changes the per-shot minimum TR, so re-check volume_tr too. Separately, lib/calc_te_tr_delays.py only warns, never raises, if your mask's ky/kz steps make the prescribed TE unreachable — it silently falls back to the minimum achievable TE, which is longer than the one you asked for — so check its printed output (or scan_info.mat's schedules[..., 2]) after your first build with a new mask.

Under the hood this is sample/external_mask.py's load_external_mask/resolve_custom_omegas, called from params.py's load_params() to build the full (Ny, Nz, Nframes) array (params.custom_omegas) and derive Nshots/R from it — see that field's comment in params.py for the details.

Trajectory computation (lib/mask2epi.py)

mask2epi_laminar/mask2epi_radial turn a 2D (ky, kz) sampling mask (Ny × Nz booleans, n of them True) into Nshots EPI echo trains of ETL samples each (Nshots * ETL == n), i.e. a partition of the sampled points into Nshots groups plus, within each group, a visiting order — the sequence in which the readout gradient steps from one (ky, kz) sample to the next. This splits into two genuinely different sub-problems: partitioning (which shot does each point belong to) and ordering (the sequence within a shot). Partitioning is sampling-pattern design and differs completely between the two variants (below); ordering is a path-optimization problem shared by both, and is the more interesting half algorithmically.

Why ordering is a bottleneck problem, not an ordinary shortest-path one

lib/make_readout_grads.py builds each shot's phase-encode "blip" gradients (the small pulses that step k-space between readout lines) as a single unit-amplitude waveform, scaled per step by that step's actual (dky, dkz). Critically, the y- and z-blip waveforms themselves — their duration and peak amplitude — are each sized once, from the single largest step seen on that axis across the whole shot. So the quantity that determines whether a shot is achievable on real gradient hardware (and stays under peripheral nerve stimulation limits) is the worst individual step, not the total distance traveled. In graph terms: given the ETL sampled points as vertices of a complete weighted graph (edge weight = a physical step cost, defined below), we want an open Hamiltonian path minimizing the maximum edge weight, not the path with minimum total edge weight. The two objectives can disagree — a path that revisits a large gap several times can still beat, on the max-edge objective, a path that takes one enormous detour to avoid it once — and the max-edge version is the bottleneck traveling salesman path problem (a.k.a. the bottleneck Hamiltonian path / bottleneck wandering salesperson problem): NP-hard, and — because it isn't a sum of independent terms — not amenable to the usual shortest-path or sum-based TSP toolbox. (For metric edge weights, the best possible polynomial-time approximation ratio is 2; nothing better exists unless P = NP.)

The edge weight used throughout is a physically weighted Chebyshev distance:

weighted_step(dy, dz) = max(|dy| * Δky, |dz| * Δkz),   Δk = (1/FOVy, 1/FOVz)

not raw index distance max(|dy|, |dz|). The FOV is generally anisotropic (216 mm in y, 144 mm in z at the default protocol), so one index-step in kz corresponds to a different physical gradient area than one index-step in ky (1.5× as much here, ~5× at the earlier 216 × 40.5 mm protocol); an unweighted metric would let the optimizer "save" index-distance by trading a cheap ky step for an expensive kz step, optimizing the wrong quantity entirely.

Partitioning

• mask2epi_laminar (ported from the MATLAB original): sweeps kz columns outer, ky rows inner in center-out order (_center_out, an fftshift-based interleave: 0, N-1, 1, N-2, ... roughly, so ky = 0 is spread across every shot rather than clustering all of it in shot 1), filling each shot to exactly ETL samples before starting the next. This guarantees ky is non-decreasing along the resulting echo train — a constraint inherited from the original design, not re-derived here.

• mask2epi_radial (this port's own addition): instead of rows, folds every point's angle about (ky, kz) = (0, 0) into [0, π) (θ mod π), so a point and its 180°-antipodal point land in the same bucket, then sorts by that folded angle and cuts the sorted list into Nshots contiguous chunks of exactly ETL. Each shot ends up an "opposite wedge pair" — a spoke through k-space center — rather than a raster row. This deliberately gives up ky non-decreasing, which is safe here because the Nyquist/odd-even ghost-correction reference scan (sequences/EPIcal.py) is keyed to readout-gradient polarity alternation, not to ky adjacency between consecutive echoes.

  Golden-angle progression through the whole sequence. Angles are measured in physical k-space (index offsets scaled by deltak = (1/fov_y, 1/fov_z)), so on a non-square plane such as 90×60 the golden-angle targets, spoke axes and start ends are golden in physical k-space, not in index space. A linear stretch keeps every sample's ray from the center in the same cyclic order, so the wedge membership barely changes (a few of ~900 points, only exact ties on the calibration grid); what changes is the angle values. And the sequence no longer restarts the golden angle every frame: sequences/ArbEPI.py's _compute_schedules first runs lib.mask2epi.radial_golden_chain (cheap, serial: an angle sort per frame), which cuts frame f's wedges starting one golden angle after the wedge acquired last in frame f−1 and continues the golden sequence from there (mask2epi_radial's tau0 and shot_offset options). Wedges are still consecutive equal-count chunks of the angle-sorted samples, taken circularly from the chosen start sample, so every shot has exactly ETL samples. The effect: the step from the last shot of one frame to the first of the next is the golden angle too (measured within 0.1°), and the wedge set is offset from frame to frame by roughly 0.4 of a wedge pitch (π/Nshots) instead of repeating, i.e. wedge angles interleave across frames. It cannot move a wedge by more than one pitch, because wedges are fixed by where the mask's samples fall; free golden-angle spoke directions need a mask generated along spokes (GOCART-style). Frame 0, and a bare mask2epi_radial(mask, ETL, Nshots) call, keep the per-frame behavior. _compute_schedules needs the complete omegas: a slice would give a different chain than the built sequence.

Ordering: a three-pass optimization

Given a shot's point set (fixed by partitioning above), what visiting order minimizes the bottleneck objective? Three passes, each building on the last:

Pass 1 — min-sum construction. Minimizing the bottleneck directly from scratch is hard to search (accepting only strictly-max-improving moves creates huge plateaus — most candidate moves don't touch the current worst edge at all). So pass 1 first builds a tour minimizing ordinary total weighted length instead — a much better-behaved objective, since almost every candidate move changes it by a nonzero amount, giving local search something to descend on everywhere. This total-length tour also turns out to matter for its own sake, not just as a warm start: a path that's long in total tends to have many needlessly-large steps scattered throughout, most of which don't happen to be the worst step and so are invisible to a bottleneck-only search.

• Radial (an unconstrained 2D point-ordering problem): nearest-neighbor construction, refined by ordinary 2-opt (edge-swap local search on the sum objective) — _sum_optimized_order.
• Laminar is more constrained and admits an exact algorithm here rather than a heuristic: within one ky row, the set of consecutive-gap sizes along kz is the same regardless of which end you start from (reversing a sweep just reverses the order the same gaps are visited in), so the only real decision per row is which of its two ends connects to the next row's chosen entry point. That's a chain of binary decisions where each choice only interacts with its immediate neighbor — a textbook 2-state dynamic program (_dp_row_directions, state = "which end does this row exit from," O(number of rows), exact — not an approximation).

Pass 2 — bottleneck refinement. Starting from pass 1's tour (not from scratch), locally reduce the worst remaining step:

• Radial: best-improvement 2-opt (_bottleneck_2opt_order) where a move is accepted only if it strictly lowers the current maximum edge. A single segment-reversal move only ever changes two edges, so re-evaluating a candidate move in O(1) is possible by tracking just the top three largest edges beforehand (at most two can be removed by any one move, so the third-largest is always still present as a bound on "everything else"). This is also where a cheap optimality certificate appears: the bottleneck value of any Hamiltonian path can never beat the bottleneck value of the minimum bottleneck spanning tree over the same points (a Hamiltonian path is itself a spanning tree), and that lower bound is just Prim's algorithm's largest added edge (_mst_bottleneck). Whenever the achieved path matches that bound, it's provably globally optimal — despite the underlying problem being NP-hard in general.
• Laminar: since each row's decision is binary, pass 2 is a greedy hill-climb that flips one row's direction at a time, keeping the flip only if it strictly reduces the current worst inter-row link (_hillclimb_row_directions_max) — provably never worse than pass 1's tour alone, though (unlike the pass-1 DP) not guaranteed globally optimal, since it's a local search over a warm-started starting point rather than an exhaustive search of the 2^(rows) state space.

Radial additionally pins two things before searching, neither of which laminar needs (its row structure already fixes them implicitly): the sample nearest (ky, kz) = (0, 0) is fixed at echo index (ETL - 1) // 2 — matching where calc_te_tr_delays.py places the nominal TE — and each half (before/after center) separately pins its own farthest-from-center point at the outer end, so the echo train visibly sweeps periphery → center → periphery like an actual radial spoke rather than an unconstrained tour that happens to wander near center at either end. Both passes 1 and 2 treat these as fixed endpoints of an open path and only reorder the interior.

Pass 3 — crossing cleanup (radial only). Passes 1-2 optimize the weighted-Chebyshev metric above, whose unit ball is a square rather than a smooth curve — it is not strictly convex. The classical argument that any self-crossing in a tour can always be "uncrossed" to strictly shorten it depends on strict convexity (the crossing quadrilateral's diagonals must be strictly longer than its sides); under a non-strictly-convex metric this can become an equality instead, so a crossing-removing move is a tie under the weighted metric and a strict-improvement-only search will never take it, even though nothing about the underlying hardware objective favors keeping the crossing. _euclidean_uncross_refine fixes this over the whole assembled shot (not each half independently — some crossings straddle the seam at the pinned center sample) in two stages, both capped so no edge may exceed the shot's already-achieved bottleneck value: first, 2-opt reversal and pairwise-exchange moves minimizing weighted Euclidean length instead (strictly convex, so the same uncrossing argument now goes through cleanly); second, for the rare residual crossing where even the fixing move is an exact tie in Euclidean length too (a swap across the pinned center position, where the two new edges are just a relabeling of the two removed ones), directly detect the geometric crossing and apply whichever available fixing move keeps every edge within the cap.

Full derivation, empirical results (bottleneck-only vs. two-pass vs. three-pass path length and crossing-count comparisons), and pointers to the brute-force-verified unit tests live in lib/mask2epi.py's module and per-function docstrings.

Executing on GE scanners

Sequences are exported to GE's .pge format in pure Python, with no MATLAB install or sibling-repo checkouts required:

from ge.ge_export import export_to_ge
export_to_ge('output/ArbEPI.seq', 'output/ArbEPI', params)

export_to_ge runs a feasibility check (hardware limits, PNS, acoustic-resonance — see below) and raises RuntimeError if it fails, then writes the .pge via ge/seq2ceq.py + ge/writeceq.py. Verified end-to-end on a small test sequence and on a full-scale default-params run (main.py's output at the time): the resulting .pge files matched freshly-regenerated real MATLAB output (seq2ceq/pge2.writeceq) byte-for-byte for two of the four default sequences and to within a single float32 ULP on a derived header field for the other two. That validation predates the single-echo GRE → dual-echo deGRE change and has not been re-run since — see CLAUDE.md for the full record.

Hardware limits are keyed off the scanner variable set in params.py (GE_MR750 or GE_UHP, see scanners.py): load_params() builds params.spec (a ScannerSpec) once and derives sys.max_grad/sys.max_slew for .seq generation from the same instance that ge/check.py/ge/writeceq.py read directly, so they can't drift out of sync with each other. main.py --ge also calls ge.ge_export.check_ge_feasibility() on all four sequences — running the hardware/PNS/acoustics check without writing a .pge — before exporting any of them, so an infeasibility surfaces immediately instead of after several full exports. PNS is a physiological safety limit, not a hardware one — PNSwt (a separate Params field, not part of ScannerSpec, since it's phantom-vs-human scan context) defaults to the IEC 60601-2-33:2022-recommended [0.8, 1.0, 0.7]; [0, 0, 0] disables the PNS check entirely and is only appropriate for phantom/non-human scanning. PNS export-blocking uses a deliberately looser threshold than IEC's own "normal operating mode" line: main.py --ge only fails (raises RuntimeError) above the 100% "first controlled mode" threshold (genuinely predicted stimulation), not IEC's 80% normal-mode line — 80–100% is reported as a WARN, not a hard failure, since (unlike the real scanner's own interlock) a Python check has no operational consequence to blocking a .seq file that's still within a standard-defined, permitted mode; see ge/check.py's PNS_NORMAL_MODE_THRESHOLD/FeasibilityReport.ok for the reasoning. Acoustic-resonance is checked but never blocks export — it's a WARN in the report, matching MATLAB's own check_grad_acoustics.m, which only ever calls warning(...), never error(...), when over threshold.

Default sequences now measure under the 80% normal-mode PNS line. PNSwt was [0, 0, 0] for most of this port's lifetime, so PNS was never actually evaluated in any --ge run (weight zero makes the per-channel contribution zero regardless of the real waveform); turning on the real weights (validated against MATLAB's per-instance pipeline to ~0.02 percentage points via ge/pns.py) revealed EPIcal/ArbEPI/GRE (the single-echo predecessor of deGRE) all at or over the limit at full hardware slew (113.2%/114.7%/100.6% on GE_UHP). The resolution is a PNS-optimized asymmetric EPI readout (POPE, Huber et al. 2026): the readout trapezoid's ramp-up slew is throttled (params.ro_slew_rise) while the ramp-down (ro_slew_fall) runs faster, since nerve-integration PNS models peak at the end of each sustained slew event. The defaults in params.py (rise/fall/blip 155/190/200 T/m/s, slew_derate 120) are the fastest found by an empirical sweep (2026-09-30, ~9000 candidates on the current default protocol) that keeps the full ArbEPI build under ~79% peak PNS: they shortened the echo spacing from 680 to 576 µs and the per-shot minimum TR by ~2.9 ms against the earlier 100/120/105 (tuned on the 240 × 240 × 45 protocol in 2026-08). Those numbers depend on the protocol; for the current build, run main.py --ge or see docs/review-findings.md's "Current baseline" table. See CLAUDE.md's PNS section for the sweep record, including why the fall/rise ratio stays mild here (the y/z blips play centered on the readout reversal, so aggressive fall/blip slews RSS-combine into a 3-channel hotspot on this whole-body gradient). uv run python -m plot.compare_readout_pns rebuilds the symmetric-vs-POPE comparison (two full sequences, identical nominal parameters and sampling masks) and writes per-variant and overlay PNS figures to output/compare_pope/. A regression test (test_arbepi_default_params_peak_pns_under_normal_mode_limit) guards the <80% property; PNS is scan-context-dependent, so re-check after changing scanner, mask seed, R/ETL, or resolution.

Copying to the scanner. For internal UM fMRI lab use: ge/coppe.py is a Python port of ../toppe/+toppe/+utils/coppe.m that copies a folder of .pge files (e.g. output/*.pge) to a scanner over SSH, auto-allocating an unused pge2 entry number for each and printing the resulting filename → entry-number mapping to enter on the scanner console.

uv run python ge/coppe.py

See ge/README.md for usage, SSH key setup, and troubleshooting.

Architecture

Only entry points (main.py, demo.ipynb) and tightly-coupled global config sit at the repo root — mirroring ../ArbEPI having params.m/main.m directly at its repo root — everything else lives in a subpackage grouped by role:

params.py                   Params dataclass + load_params() (replaces params.m)
scanners.py                 ScannerSpec hardware profiles (GE_MR750, GE_UHP)
main.py                     Generate all 4 sequences, mirrors ArbEPI/main.m; --plot/--ge flags
demo.ipynb                  main.py step by step, with the diagnostic plots and the GE feasibility check
sample/                     Sampling mask generators (caipi, ticaipi, pd, rand); external_mask.py
                              is a manual escape hatch for a collaborator-supplied mask, not a
                              fifth params.sampling_method -- see CLAUDE.md
lib/
  mask2epi.py                Core algorithm: partitions a 2D mask into EPI trajectories
                              (mask2epi_laminar / mask2epi_radial, selected by params.epi_trajectory)
  trap4ge.py                 GE-raster trapezoid rounding (from ../PulCeq)
  shot_start.py              Per-shot spoiler draws + fat-sat/excitation block shared by ArbEPI and EPIcal
  epi_shot.py                The pulses, gradients and delays of one EPI shot, built once for ArbEPI and EPIcal
  seq_io.py                  Every generator's epilogue (timing check, .seq write) and scan_info.mat access
  (RF/gradient helpers, TE/TR delay calculation)
sequences/                  ArbEPI, EPIcal, deGRE, noise sequence assembly
  meGRE.py                       20-echo monopolar GRE prototype (B0 + R2* + coil maps); `main.py --gre megre`
plot/
  plotting.py                Diagnostic plots (mask/PSF/trajectory/single-TR/PNS)
  plot_last_run.py           Drives plotting.py against the most recent output/ run
  compare_readout_pns.py     Symmetric-vs-POPE readout comparison (PNS figures + table)
ge/                      Pure-Python GE .pge toolchain (Pulseq -> GE), wired into main.py --ge
  ge_export.py                 GE feasibility checking / .pge export entry points
  ceq.py                       Ceq/ParentBlock/Segment data contract
  blocks.py                    Block type/dynamics/comparison (port of PulCeq's compareblocks.m etc.)
  seq2ceq.py                   Pulseq Sequence -> Ceq (port of seq2ceq.m)
  writeceq.py                  Ceq -> binary .pge (port of writeceq.m)
  read_pge.py                  .pge binary reader, used only for validation
  pns.py                       PNS check (port of pge2.pns.m)
  validate_pns.py              Compares pns.py/acoustics.py against MATLAB reference output (not a pytest test)
  acoustics.py                 Acoustic-resonance check (port of check_grad_acoustics.m)
  check.py                     Combines hardware/PNS/acoustics into one FeasibilityReport
  validate_against_matlab.py   Field-by-field / byte-for-byte comparison vs. MATLAB output
  coppe.py                     Copy .pge files to the scanner over SSH, auto-allocating entry numbers
                                (port of ../toppe/+toppe/+utils/coppe.m; UM lab-internal use)
preprocess/                Raw ScanArchives -> one reconstruction-ready file per sequence
                              (.venv-preprocessing), ported from ../epi-preprocessing -- see preprocess/README.md
  preprocess.py                 PreprocessConfig + preprocess(): gridding (cached), coil compression,
                                 maps, output <seq>_preprocessed.h5
  batch_preprocess.py           Command line: several sequences in one go
  coils.py                      Noise whitening; GCC coil compression
  epi_gridding.py               1D NUFFT ramp-sample regridding (sigpy)
  oephase.py                    Odd/even (Nyquist ghost) phase estimation + correction
  smaps.py                      ESPIRiT sensitivity maps (sigpy) + resize/mask/smooth/normalize
  b0map.py                      B0 field map via julia/b0map.jl (MRIFieldmaps.jl + ROMEO.jl)
  r2star.py                     R2* fit over the deGRE echoes (a placeholder with the dual-echo deGRE)
  grid_resize.py                deGRE grid -> EPI grid
  utils.py                      ScanArchive (GERecon), scan_info.mat and NIfTI I/O; QA figures
  julia/                        Self-contained Julia project (Project.toml + Manifest.toml, pinned)
  demo.ipynb                    Every step on one real dataset, and the settings table
recon/                        Image reconstruction (.venv-recon, uv sync), on PyTorch/mirtorch -- see recon/README.md:
                              min_x 0.5||Ax - y||^2 + g(x)
  operators.py               Encoding operators A: SENSE, SENSE_B0 (B0 phase accrual, time-segmented,
                                 L=32 in production -- swept, see CLAUDE.md), SENSE_B0_R2star (+ R2* decay)
  regularizers.py               g(x): MultiScaleLowRank (multi-scale low-rank, patch SVST), WaveletTV (3D wavelet + TV)
  solvers.py                    pogm_restart (PGM/FPGM/POGM), pdhg (Condat-Vu, via mirtorch), cg
  sense.py                      Iterative SENSE driver: --reg {none,mslr,wavelet-tv}, --B0, --R2star
  rss.py                        Root-sum-of-squares, GPU-batched over frames
  utils.py                      I/O, operator norms, tSNR report, one-off L sweep/benchmark and
                                 validation vs. ../mslr-recon (python -m recon.utils ...)
  demo.ipynb                    Every recon type on one real dataset
postprocess/                 Operations on recon/ output before analyze/ (main .venv, `uv sync --extra analyze`) -- see postprocess/README.md:
  motion.py                     Rigid-body motion correction: FSL MCFLIRT if found, else scipy registration; shared resampling
  brainmask.py                  Brain mask: BET on the T1-weighted deGRE, registered to the EPI mean (flirt); applied to the corrected series
  __main__.py                   `python -m postprocess --recon ... --out ... [--brain-mask-from ...]`
analyze/                     Task activation maps from reconstructions (main .venv, `uv sync --extra analyze`) -- see analyze/README.md:
  design.py                     Task blocks x canonical HRF, DCT drift terms, constant
  glm.py                        GLM with SPM-style AR(1) prewhitening; t -> z, FDR
  mask.py                       Brain mask: FSL bet if on PATH, else brainextractor
  io.py                         recon .h5 (summed image, per-scale components) in, NIfTI out
  run.py, __main__.py           analyze_recon and `python -m analyze`
  compare.py                    Several reconstructions or scales: maps, time courses, spectra (`python -m analyze.compare`)
simulate_fmri/               Simulated scan sessions acquired with ArbEPI's schedules (.venv-simulate, uv sync)
                              -- see simulate_fmri/README.md
  session.py                    Raw readouts of all four scans for preprocess/ to read, plus the truth
  forward.py                    Signal equation: ramp-sampled EPI readouts and deGRE k-space from spins on a fine grid
  b0.py                         B0 map from the head's susceptibility (dipole convolution, air cavities, shim)
  coils.py                      Receive array defined in space; coil noise covariance
  task.py                       Block paradigm, canonical HRF, the BOLD response and its regional delays
  physio.py                     Physiological noise (BOLD-like, cardiac, respiratory, drift)
  study.py                      Repetitions of an experiment through preprocess/ and recon/, resumable
  analysis.py                   Activation maps on the T1w image, ROC curves, mixed-binomial test-retest model
  ideal.py                      SNAKE-fMRI's engine on the Cartesian grid, written straight to recon's format
  sampler.py, engine.py,        ideal.py's SNAKE sampler, engine, activation handler and export
    handlers.py, export.py
  protocol.py, phantom.py       scan_info.mat reader; BrainWeb and analytic phantoms
  demo.ipynb                    A visual-motor study: four runs, spectra, activation maps, ROC, reliability
tests/                       Unit tests (pytest)
docs/review-findings.md      Numbered backlog of code-review findings (open and resolved)
docs/TODO.md                 Planned work
docs/pipeline/               Pipeline diagram: make_pipeline.py generates pipeline-{light,dark}.svg

Index convention: internal computation is 0-based throughout (mask2epi's schedule, sampling masks, etc.). output/scan_info.mat is written with schedules' (ky, kz) channels converted to 1-based (matching what MATLAB-side reconstruction code expects) — parts is already a 1-based shot label with 0 = unsampled, so it needs no conversion.

output/scan_info.mat (kxo/kxe, schedules/parts, and a snapshot of the scan scalars preprocess/ needs — consolidating what used to be three separate files, samp_locs.mat/params.mat/kxoe<Nx>.mat) is written via hdf5storage in MATLAB v7.3 format (HDF5-based), matching the original MATLAB code's save(..., '-v7.3'). scipy.io.savemat/loadmat cannot write or read v7.3 at all — use hdf5storage.loadmat (or h5py directly) to read these files from Python, not scipy.io.loadmat.

See ../ArbEPI/README.md for background on the sampling methods and mask2epi_laminar's original (MATLAB) partitioning design — see Trajectory computation above for the ordering optimization and mask2epi_radial, both this port's own addition with no MATLAB counterpart.

preprocess/ is a separate pyproject.toml optional-dependency group (preprocessing) with its own venv (Python 3.10, since GE's Orchestra SDK is ABI-locked to it) -- see preprocess/README.md for setup, usage, outputs and design. The SDK itself is distributed as a GitHub release, accessible by request: GEHC-External/MR-Orchestra-SDK-Python.

recon/ is likewise its own optional-dependency group (recon: PyTorch, mirtorch, ...) in a dedicated venv (.venv-recon), but unlike preprocess/ it is installed from uv.lock (UV_PROJECT_ENVIRONMENT=.venv-recon uv sync --extra recon --extra test); it stays separate from .venv because uv sync is exact and would remove torch whenever the main environment is synced without the recon extra -- see recon/README.md for setup, commands, inputs and outputs.

simulate_fmri/ is a third optional-dependency group (simulate: SNAKE-fMRI, pinned to a GitHub commit) with its own venv, .venv-simulate, installed from uv.lock together with preprocessing and recon so a simulated session can be preprocessed, reconstructed and scored in one environment (UV_PROJECT_ENVIRONMENT=.venv-simulate uv sync --extra simulate --extra preprocessing --extra recon --extra test; no GERecon needed) -- see simulate_fmri/README.md for what it models, the measurements its sizes come from, and how it relates to simulate_fmri/from_real.py.

Documentation

The documentation site (Read the Docs; configured by .readthedocs.yaml) is built from this README, the per-package READMEs and the docstrings: the guides are these README files themselves, and the API reference is generated by Sphinx autosummary, so it follows the code. To build it locally:

uv run --extra docs sphinx-build -b html docs docs/_build/html

then open docs/_build/html/index.html. See CLAUDE.md's "Documentation site" section for how the pieces fit.

Differences vs. MATLAB original

The sequence-generation side (params.py, main.py, sample/, lib/, sequences/, plot/) is the port of the MATLAB original; ge/, preprocess/ and recon/ port other toolchains (see References). A few things from the original MATLAB repo are handled differently — see below.

• GE .pge export: ge/ is a pure-Python port of the PulCeq/pge2 toolchain (seq2ceq, writeceq, and the pge2.pns/check_grad_acoustics feasibility checks), and is the operative path for ge/ge_export.py/main.py --ge — no MATLAB round trip, no sibling ../pulseq/../toppe/../PulCeq/../ArbEPI checkouts required. Validated field-by-field and, for two of the four default sequences at the time (before GRE became today's dual-echo deGRE), byte-for-byte against real MATLAB output, including end to end through main.py --ge itself (see ge/'s module docstrings and CLAUDE.md). See Executing on GE scanners above.
• Fat-sat RF pulse: a min-phase SLR pulse like the MATLAB original's (toppe.utils.rf.makeslr), designed by lib/slr.py, a small port of SigPy's dzrf. It is longer and narrower than the original's (TBW 2 over 6 ms instead of TBW 3 over 4 ms) so it leaves water alone: the original's parameters tip on-resonance water by 9°, and the Gaussian this port used before by 26°, which caused most of the temporal fluctuation in static phantom scans (docs/review-findings.md item 255).
• Water excitation (not in the MATLAB original): params.excitation = 'water' swaps the fat-sat pulse, its crusher and the sinc excitation for a slab-selective binomial 1-3-3-1 water-excitation train (lib/make_water_excitation.py), 8.0 ms shorter per shot on the default protocol. The default since 2026-09-30; 'fatsat' remains available.
• Plotting (plot/plotting.py) covers the sampling mask, k-space trajectory, point-spread-function, and single-TR pulse-diagram plots. The interactive scroll/slider mask viewer from the MATLAB repo is not ported. The single-TR plot (plot_one_tr) isn't a custom pulse-diagram renderer — it's a thin wrapper around pypulseq's own Sequence.plot(). Mask/PSF/trajectory plots take an optional frame_idx to select one frame out of a multi-frame run; see plot/plotting.py's module docstring for why the per-frame trajectory plot draws through exact-sliced ADC samples rather than the fine continuous line pypulseq's calculate_kspace() returns for the whole sequence. plot/plot_last_run.py drives all of this against the most recent output/ run (uv run python main.py --plot, or standalone via uv run python -m plot.plot_last_run).
• Poisson-disc sampling (sample/pd_sample.py) is a local reimplementation, not a dependency on SigPy (whose sigpy.mri.poisson both ../ArbEPI/lib/pd_sample.m and this module's algorithm are based on). SigPy was tried directly and rejected: its own poisson() has an unbounded while slope_min < slope_max binary-search loop with no iteration cap, which hangs forever on small/coarse grids where no achievable density slope lands within tol of the target acceleration (reproduced independently of any code in this repo — confirmed via sigpy.mri.poisson alone). Separately, both ../ArbEPI/lib/pd_sample.m and this module's own first version had a different bug (not reseeding the point-placement RNG identically on every binary-search iteration, unlike real SigPy), which made the search non-convergent and slow rather than truly infinite. A third bug (a missing nx*ny active-list cap that real SigPy has) let the point-placement active list grow unboundedly in the same radius-floor/dense-center regime. pd_sample.py's docstring has the full writeup of these and of later fixes (the seed point is now marked occupied, the exact-count fill stays inside the crop_corner ellipse, and the calibration rectangle is exempt from that crop); the local implementation fixes all of them and adds a bounded outer-search iteration cap (max_search_iters) that SigPy itself lacks. The point-placement core is still JIT-compiled with numba -- a narrow, single-function dependency (unlike depending on the sigpy package wholesale) -- since even with all three fixes it's an inherently sequential loop that can run hundreds of thousands of iterations for worst-case seeds, and pure Python can't get there without JIT (measured ~1-12s/frame in pure Python vs. ~0.02-0.2s/frame JIT-compiled, at production scale).

References

Source repositories this port is derived from or ports code from:

• rextlfung/ArbEPI — the MATLAB/Pulseq original this repo ports.
• rextlfung/epi-preprocessing — the MATLAB original preprocess/ ports.
• HarmonizedMRI/PulCeq — ge/ ports seq2ceq.m/writeceq.m/pge2.pns.m/check_grad_acoustics.m from here.
• HarmonizedMRI/B0shimming — sequences/deGRE.py is adapted from this repo's writeB0.m.
• ../mslr-recon (Julia, MIRT.jl) — the multi-scale low-rank reconstruction recon/ started as a port of, and was validated against on real scanner data (see CLAUDE.md).

Libraries this repo depends on:

• pypulseq — Pulseq sequence assembly (MIT).
• SigPy — ESPIRiT sensitivity maps and NUFFT ramp-sample gridding in preprocess/ (BSD); sample/pd_sample.py is also an independent reimplementation of sigpy.mri.poisson's Poisson-disc algorithm (see that module's docstring for the bugs found in both it and the MATLAB original that motivated the reimplementation).
• Numba — JIT compilation for pd_sample.py's point-placement core.
• PyTorch and MIRTorch — the linear operators, time-segmented B0 approximation (mri_exp_approx) and primal-dual solver behind recon/ (BSD).
• SNAKE-fMRI — in simulate_fmri/: the phantom of both modes, and the BOLD regressor, acquisition engine, MRD file and parallel driver of the ideal mode (MIT; Comby, Vignaud, Ciuciu, arXiv:2404.08282). Its phantom is the BrainWeb digital brain, downloaded on first use by brainweb-dl.
• PyWavelets — wavelet filter taps for recon/regularizers.py's Wavelet3D (MIT).
• BART — considered for coil compression/whitening/parallel imaging and explicitly dropped in favor of plain numpy + SigPy (see CLAUDE.md); noted here since several docstrings describe what this repo does instead of calling it.
• MRIFieldmaps.jl — B0 field map estimation (MIT), invoked as a subprocess from preprocess/b0map.py via the self-contained Julia project in preprocess/julia/; requires a separate julia install (see preprocess/README.md), not otherwise part of this repo's Python dependency set.
• ROMEO.jl — phase unwrapping (MIT), used by preprocess/julia/b0map.jl to build the field map solve's initial guess; same Julia project as MRIFieldmaps.jl above.

Proprietary, not distributed with this repo:

• GE Orchestra SDK / GERecon — required to read raw GE ScanArchive files in preprocess/utils.py; installed separately by the user under GE's own license terms, not covered by this repo's license.

Algorithms and standards:

• ESPIRiT: Uecker M, Lai P, Murphy MJ, et al. "ESPIRiT — an eigenvalue approach to autocalibrating parallel MRI: where SENSE meets GRAPPA." Magn Reson Med. 2014;71(3):990-1001.
• Shinnar-Le Roux (SLR) pulse design: Pauly J, Le Roux P, Nishimura D, Macovski A. "Parameter relations for the Shinnar-Le Roux selective excitation pulse design algorithm." IEEE Trans Med Imaging. 1991;10(1):53-65. (Implemented by lib/slr.py, ported from SigPy's sigpy.mri.rf.dzrf, for the fat-sat pulse in lib/make_fatsat_rf.py; the MATLAB original's toppe.utils.rf.makeslr implements the same design.)
• Golden-angle ordering: Winkelmann S, Schaeffter T, Koehler T, Eggers H, Doessel O. "An optimal radial profile order based on the Golden Ratio for time-resolved MRI." IEEE Trans Med Imaging. 2007;26(1):68-76. (open-access copy; applied in lib/mask2epi.py's _golden_angle_shot_order, see that function's docstring.)
• Regularized B0 field map estimation: Lin CY, Fessler JA. "Efficient Regularized Field Map Estimation in 3D MRI." IEEE Trans Comput Imaging. 2020;7:60-73. (Implemented by MRIFieldmaps.jl, invoked from preprocess/julia/b0map.jl.)
• Phase unwrapping: Dymerska B, Eckstein K, Bachrata B, Siow B, Trattnig S, Shmueli K, Robinson SD. "Phase unwrapping with a rapid opensource minimum spanning tree algorithm (ROMEO)." Magn Reson Med. 2021;85(4):2294-2308. (Implemented by ROMEO.jl, used in b0map.jl to build the field map solve's initial guess -- see the script's header for why the naive two-point guess isn't safe to use directly.)
• Coil compression: Zhang T, Pauly JM, Vasanawala SS, Lustig M. "Coil compression for accelerated imaging with Cartesian sampling." Magn Reson Med. 2013;69(2):571-582. (GCC, preprocess/coils.py.)
• Multi-scale low rank: Ong F, Lustig M. "Beyond low rank + sparse: multiscale low rank matrix decomposition." IEEE J Sel Top Signal Process. 2016;10(4):672-687, arXiv:1507.08751. (recon/regularizers.py's MultiScaleLowRank and its lambda weights.)
• Time-segmented B0 correction: Sutton BP, Noll DC, Fessler JA. "Fast, iterative image reconstruction for MRI in the presence of field inhomogeneities." IEEE Trans Med Imaging. 2003;22(2):178-188. (recon/operators.py's SENSE_B0 sign convention.)
• PNS-optimized EPI readout (POPE): Huber et al., "PNS Optimized Pulses for EPI," bioRxiv 2026, doi:10.64898/2026.07.22.739360. (lib/make_readout_grads.py's asymmetric readout.)
• Peripheral nerve stimulation: IEC 60601-2-33:2022, the PNS prediction model ge/pns.py implements (ported from PulCeq's pge2.pns.m).