"""Conversion of a Pulseq sequence into the scanner's segmented binary IR."""
from __future__ import annotations
from collections.abc import Sequence
from dataclasses import dataclass, field
from enum import IntEnum
from pathlib import Path
from typing import Any
import pypulseqpp as pp
from .._accelerators import require
from ..mrd._sequence import read_chain
from ._checks import SarRatio
from ._source import conversion_payload
def _rasters(system: pp.Opts) -> tuple[float, ...]:
"""Return the RF, gradient, ADC and block rasters, in us."""
return (
system.rf_raster_time * 1e6,
system.grad_raster_time * 1e6,
system.adc_raster_time * 1e6,
system.block_duration_raster * 1e6,
)
[docs]
@dataclass(frozen=True)
class Quantity:
"""The format one quantity is stored in, and what one integer step means.
``step`` is the quantity's SI value of one integer step, so a stored value
``v`` means ``v * step``. It is 0 for a float, which is stored in its SI
unit.
"""
format: Format = Format.FLOAT32
step: float = 0.0
[docs]
@dataclass(frozen=True)
class VendorProfile:
"""What a cache holds its numbers as, for the machine that reads it.
A sequencer that plays integers is given integers, already scaled, so it
converts nothing while it plays; a reader that works in SI is given floats.
Every quantity defaults to a float in its SI unit.
The steps belong to the machine and arrive with its other limits. Nothing
here supplies one.
"""
grad_sample: Quantity = field(default_factory=Quantity)
grad_amplitude: Quantity = field(default_factory=Quantity)
rf_sample: Quantity = field(default_factory=Quantity)
rf_amplitude: Quantity = field(default_factory=Quantity)
rf_phase: Quantity = field(default_factory=Quantity)
rf_frequency: Quantity = field(default_factory=Quantity)
[docs]
def as_pairs(self) -> tuple[float, ...]:
"""Return a format and a step per quantity, as the extension takes them."""
return tuple(
value
for quantity in (
self.grad_sample,
self.grad_amplitude,
self.rf_sample,
self.rf_amplitude,
self.rf_phase,
self.rf_frequency,
)
for value in (float(int(quantity.format)), float(quantity.step))
)
[docs]
@dataclass(frozen=True)
class Grouping:
"""Rule by which the blocks of a repetition are grouped into virtual segments.
The repetition is cut into contiguous runs of blocks, each a virtual
segment the interpreter prepares once and plays as segment instances. A
boundary may fall between two blocks only where every gradient is within
:attr:`boundary_gradient_hz_per_m` of zero at the join. A boundary must
fall where the ``NOROT`` or ``PMC`` label changes, since the interpreter
sets one prescription rotation per segment instance; conversion fails when
such a change falls under a gradient. The ``split_*`` attributes add
boundaries.
The default suits an interpreter that begins a segment by setting its
gradients, so that a boundary falls only where they are at rest.
Attributes
----------
boundary_gradient_hz_per_m
Largest gradient amplitude at a join, in Hz/m, at which a boundary may
fall. An interpreter that can begin a segment under a gradient states
a large value.
split_by_pulses, split_by_readouts
Whether runs that play different RF events, or different ADC events,
are different virtual segments.
split_navigators
Whether a navigator readout is a virtual segment of its own.
split_edge_delays
Whether a pure-delay block at the edge of a virtual segment is a
virtual segment of its own. An interpreter that inserts
transmit-to-receive switching time at every segment boundary needs as
few segments as possible, and sets this False.
"""
boundary_gradient_hz_per_m: float = 100.0
split_by_pulses: bool = True
split_by_readouts: bool = True
split_navigators: bool = True
split_edge_delays: bool = True
[docs]
def as_values(self) -> tuple[float, ...]:
"""Return the rule as the extension takes it, in the order C declares it."""
return (
float(self.boundary_gradient_hz_per_m),
float(self.split_by_pulses),
float(self.split_by_readouts),
float(self.split_navigators),
float(self.split_edge_delays),
)
[docs]
@dataclass(frozen=True)
class WaveBudget:
"""What a playout's waveform memory affords the waves.
A property of the playout, not of the sequence: the C library's
``pulseg_wave_budget``. :func:`convert` lays the waves out for it, and a
playout refuses a cache laid out for another budget.
Attributes
----------
max_samples
Samples each gradient axis holds for waves.
raster_us
The playout's gradient raster, in µs per sample.
load_us_per_sample
Time to sample and load one sample on one axis, in µs; 0 leaves the
loading unchecked.
headroom
Share of the playout's time its loading may take.
slots
Slots per segment position a streamed layout rings through, at least
2: one more than the segment instances the loading may run ahead of
the playout.
Raises
------
ValueError
If the raster or the headroom is not positive, the memory or the load
rate is negative, or there are fewer than two slots.
"""
max_samples: int
raster_us: float
load_us_per_sample: float = 0.0
headroom: float = 0.5
slots: int = 2
def __post_init__(self) -> None:
if self.raster_us <= 0.0 or self.headroom <= 0.0:
raise ValueError("the raster and the headroom must be positive")
if self.max_samples < 0 or self.load_us_per_sample < 0.0:
raise ValueError("the memory and the load rate cannot be negative")
if self.slots < 2:
raise ValueError("a streamed layout rings through at least two slots")
def _held(budget: WaveBudget | None) -> tuple[int, float, float, float, int] | None:
"""Return the budget as the extension takes it."""
if budget is None:
return None
return (
budget.max_samples,
budget.raster_us,
budget.load_us_per_sample,
budget.headroom,
budget.slots,
)
[docs]
def cache_path(seq_path: Path | str, cache_ext: str = ".pseg") -> Path:
"""Return the cache file of a sequence: its last suffix replaced by ``cache_ext``."""
return Path(seq_path).with_suffix(cache_ext)
[docs]
def chain(seq_path: Path | str) -> list[Path]:
"""Return the files of the ``NextSequence`` chain starting at a sequence file, in play order.
Raises
------
FileNotFoundError
If a file of the chain does not exist.
ValueError
If the chain names a file it has already played.
"""
return [path for path, _ in read_chain(seq_path)]
[docs]
def convert(
seq_path: Path | str,
system: pp.Opts,
*,
fov_offset: Sequence[float] | None = None,
vendor: int = 0,
label_column_map: Sequence[int] = (0, 1, 2),
cache_ext: str = ".pseg",
verify_signature: bool = True,
sar_ratios: Sequence[SarRatio] | None = None,
wave_budget: WaveBudget | None = None,
profile: VendorProfile | None = None,
grouping: Grouping | None = None,
designed: list[tuple[Path, pp.Sequence]] | None = None,
) -> Path:
"""Segment a sequence file and write its IR cache beside it.
The ``NextSequence`` chain starting at the file is read as the
subsequences of one scan. An existing cache at the destination is
replaced. RF vendor statistics are left at zero: ``vendor`` only tags the
cache for the reader that loads it, which must be built for that vendor.
The files are expected in the logical frame, and each is moved to
``fov_offset`` with :func:`prescribe` before it is segmented. The
prescription's rotation is not applied here: the scanner plays the cache
through its rotation matrix, composed after each block's own rotation.
ADC phase modulation is not stored in the cache: the reconstruction proxy
applies it to the received samples from the sequence
(:class:`~pulserver.proxy.SequenceTable`).
Parameters
----------
seq_path
Text or binary Pulseq file.
system
Limits and rasters the scan is segmented under.
fov_offset
Translation of the field-of-view centre along the logical readout,
phase and slice axes, in metres. None or zero leaves the files as
designed.
vendor
``PULSEG_VENDOR_*`` code; 0 is vendor-neutral.
label_column_map
Pulseq label state indices filling the three ADC label columns:
0 SLC, 1 PHS, 2 REP, 3 AVG, 4 SEG, 5 SET, 6 ECO, 7 PAR, 8 LIN, 9 ACQ.
cache_ext
Extension of the cache file, dot included.
verify_signature
Refuse a file whose contents do not match the signature it carries. A
file carrying none is read either way.
sar_ratios
One per file of the chain, as :func:`sar_ratios` returns them, written
into each subsequence of the cache; zero when None.
profile
What the cache holds its numbers as; every quantity a float in its SI
unit when left out.
grouping
How the blocks of a repetition are grouped into virtual segments.
Left out, ``Grouping()``: a boundary falls only where every gradient
is within 100 Hz/m of zero, and must fall where the ``NOROT`` or
``PMC`` label changes.
wave_budget
The waveform memory of the playout the cache is for, which the cache
lays the waves out in (:func:`plan_waves`); None holds every wave at
once on the gradient raster of the chain's first file.
designed
The chain as :func:`pulserver.mrd.designed_chain` returns it,
segmented in place of reading and verifying the files, and moved to
``fov_offset`` in place.
Returns
-------
Path
The cache file.
Raises
------
ValueError
If a file of the chain cannot be read, verified or segmented,
``sar_ratios`` does not give one per file, or the waves fit neither
layout of ``wave_budget`` or cannot be loaded in time.
OSError
If no cache was written.
"""
cache_path(Path(seq_path), cache_ext).unlink(missing_ok=True)
payload = _payload(Path(seq_path), system, verify_signature, fov_offset, designed)
return _write_cache(
seq_path,
system,
payload,
vendor=vendor,
label_column_map=label_column_map,
cache_ext=cache_ext,
sar_ratios=sar_ratios,
wave_budget=wave_budget,
profile=profile,
grouping=grouping,
)
def _write_cache(
seq_path: Path | str,
system: pp.Opts,
payload: list[dict[str, Any]],
*,
vendor: int = 0,
label_column_map: Sequence[int] = (0, 1, 2),
cache_ext: str = ".pseg",
sar_ratios: Sequence[SarRatio] | None = None,
wave_budget: WaveBudget | None = None,
profile: VendorProfile | None = None,
grouping: Grouping | None = None,
) -> Path:
"""Segment the libraries :func:`_payload` returned and write the cache, as :func:`convert`.
The libraries are copies, so the sequences they were taken from may
change while this runs; the segmentation releases the GIL.
"""
seq_path = Path(seq_path)
target = cache_path(seq_path, cache_ext)
target.unlink(missing_ok=True)
if sar_ratios is not None:
if len(sar_ratios) != len(payload):
raise ValueError(
f"expected one SAR ratio per file of the chain, {len(payload)}; "
f"got {len(sar_ratios)}"
)
for libraries, ratio in zip(payload, sar_ratios, strict=True):
libraries["reserved"]["vop_sar_ratio"] = float(ratio.local_sar)
libraries["reserved"]["vop_global_sar_ratio"] = float(ratio.global_sar)
require("convert_libraries")(
payload,
str(seq_path),
*_rasters(system),
int(vendor),
list(label_column_map),
cache_ext,
_held(wave_budget),
None if profile is None else profile.as_pairs(),
None if grouping is None else grouping.as_values(),
)
if not target.is_file():
raise OSError(f"no cache was written for {seq_path}")
return target
[docs]
def summary(
seq_path: Path | str,
system: pp.Opts,
*,
cache_ext: str | None = None,
label_column_map: Sequence[int] = (0, 1, 2),
) -> dict[str, Any]:
"""Return the segmentation of a sequence: subsequences, segments and readouts.
Each subsequence lists its unique RF definitions under ``rf``: the flip
angle in degrees at the largest amplitude the definition plays, as
``pypulseqpp.Sequence.rf_flip_angles`` gives it; and the bandwidth at half
the spectral peak, the number of bands, each band's offset from the
carrier and the widest band's bandwidth, all in Hz, as
``pypulseqpp.calc_rf_bandwidth`` measures them; and ``b1sq_integral_s``,
the integral of the squared envelope scaled to unit peak, in s:
``pypulseqpp.calc_rf_power``'s energy over its peak power, both summed
over a dynamic pTx pulse's channels.
``vop_sar_ratio`` and ``vop_global_sar_ratio`` are those the cache was
written with, zero when the chain is read and segmented again.
``readout_labels`` lists, per readout in play order, the values of the
three labels ``label_column_map`` selects, as in force at that readout.
``waves`` lists the subsequence's waves, as :func:`play` indexes
them: each one's point count, and its largest magnitude along x, y and z
over that of the largest gradient event it combines. ``grad_raster_us``
is the raster the file's gradient shapes are sampled on, in µs.
With ``cache_ext``, the cache beside the file is loaded instead of the
chain being read and segmented again; this build loads only vendor-neutral
caches, and only when the size recorded in the cache matches the file.
Raises
------
ValueError
If the file cannot be read or the cache cannot be loaded.
"""
seq_path = Path(seq_path)
if cache_ext is None:
return require("summary_from_libraries")(
_payload(seq_path, system, verify_signature=False),
*_rasters(system),
list(label_column_map),
)
return require("summary_from_cache")(
str(cache_path(seq_path, cache_ext)), seq_path.stat().st_size
)
[docs]
def prescribe(sequence: pp.Sequence, fov_offset: Sequence[float]) -> pp.Sequence:
"""Move a logical-frame sequence to a prescribed field-of-view centre, in place.
An offset along the slice axis becomes an RF frequency, an in-plane offset
an RF and ADC phase, each referenced to the excitation it follows; blocks
labelled ``NOPOS`` are exempt. The gradient area is counted from the
sequence's first block. A block that carries a rotation extension is moved
by the gradients it plays: those it draws, turned by that rotation.
Parameters
----------
sequence
One file of a chain, as designed.
fov_offset
Translation along the logical readout, phase and slice axes, in metres.
Returns
-------
pypulseqpp.Sequence
``sequence`` itself.
Raises
------
ValueError
If ``fov_offset`` is not three values.
"""
shift = tuple(float(v) for v in fov_offset)
if len(shift) != 3:
raise ValueError(f"fov_offset takes three values, got {len(shift)}")
if any(shift):
pp.TransformFOV(translation=shift, through_rotation=True).apply_to_sequence(
sequence, in_place=True
)
return sequence
[docs]
def play(
seq_path: Path | str, cache_ext: str = ".pseg", *, waveforms: bool = False
) -> dict[str, Any]:
"""Walk the cache beside a sequence file as the scanner's playout does.
The cache is loaded by the C library a scanner links, and its execution
stream is walked with that library's cursor: one entry per played block,
in play order, across the subsequences of the chain. This build loads only
vendor-neutral caches, and only when the size recorded in the cache
matches the file.
Returns
-------
dict of str to ndarray
One entry per played block in each:
- ``subsequence``, ``segment``: chain file and cache segment indices;
- ``duration_us``: block duration, in µs;
- ``rf_amp_hz``, ``rf_freq_hz``, ``rf_phase_rad``: RF amplitude
(gamma B1) and frequency and phase offsets, with ppm offsets
resolved at the field strength the cache was converted under; 0
without RF;
- ``rf_use``: the ``PULSEG_RF_USE_*`` code of the RF event, 1 for an
excitation and 2 for a refocusing pulse; a pulse the file leaves
unlabelled takes the use pypulseqpp detects when the chain is read,
and an event a cache carries without one is a refocusing pulse at a
flip angle of 162 to 198 degrees and an excitation otherwise; 0
without RF;
- ``rf_delay_us``: RF delay from the block's start, in µs;
- ``rf_channels``: the transmit channels the RF waveform holds, one
after another over one time base for a dynamic pTx pulse; 0 without
RF;
- ``rf_grad_constant``, ``rf_grad_level``: 1 where every instance of
the block's position plays its RF pulse under one gradient,
steady from the pulse's first sample to its last as
``pypulseqpp.Sequence.rf_gradients`` finds it, and the block
carries no rotation, so the scanner may move the excitation by a
carrier offset; and ``(blocks, 3)``, that gradient along x, y and z
over the amplitude of the event playing it there;
- ``gradient_hz_per_m``: ``(blocks, 3)``, the amplitude along the
logical x, y and z axes: the factor on the shape of the gradient
event, or on the wave, the block plays there, each normalised to a
largest magnitude of one;
- ``wave``: the wave the block plays, as indexed by the
subsequence's ``waves`` in :func:`summary`, or -1 for a block that
plays its gradient events as they are. A block plays one at every
segment position where an instance carries a rotation other than
the identity, or plays a gradient definition or shape other than
the one the position's events are prepared with: its gradient
events combined and turned by its rotation, so that the
prescription's rotation is the only one the scanner applies;
- ``norot``, ``nopos``: the block's NOROT and NOPOS flags;
- ``adc``, ``adc_freq_hz``, ``adc_phase_rad``: whether the block
acquires, and its frequency and phase offsets;
- ``adc_delay_us``, ``adc_dwell_ns``, ``adc_samples``: the ADC delay
from the block's start, dwell time and sample count; 0 without ADC;
- ``trid``: the TRID group in force, 0 when ungrouped.
With ``waveforms``, also:
- ``rf_center_us``: the time of the RF centre the design records, from
the block's start, in µs; NaN without RF;
- ``rf_time_us``, ``rf_waveform_hz``: the samples of every played RF
pulse, concatenated in play order: their times from the block's
start, in µs, and the instance's amplitude times the magnitude
shape and ``exp(i phase)`` of the phase shape its definition
carries, complex, in Hz. The pulse plays each sample turned by the
phase offset and advancing at the frequency offset from its start.
The channels of a pTx pulse follow one another, each over the one
time base;
- ``rf_span``: ``(blocks, 2)``, the start and stop of each block's
RF in those arrays; empty without RF;
- ``gradient_time_us``, ``gradient_waveform_hz_per_m``: the corners of
every played gradient, concatenated in play order: their times from
the block's start, in µs, and the gradient there, the amplitude
times the shape or wave the instance plays. A waveform on the
gradient raster holds its end values over the half raster
intervals before its first sample and after its last. A rotated
wave takes the corners of every gradient event it combines, or the
raster centres across them when one is an arbitrary gradient on
the raster;
- ``gradient_span``: ``(blocks, 3, 2)``, the start and stop of each
block's gradient along x, y and z in those arrays; empty without
one;
- ``adc_phase_modulation_rad``: the phase modulation of every played
readout, one phase per sample in radians, concatenated in play
order; the receiver phase of a sample is the ADC phase offset, plus
its frequency offset times the time since the ADC's start, plus
this. Empty for a cache converted here, which leaves the
modulation to the reconstruction proxy;
- ``adc_modulation_span``: ``(blocks, 2)``, the start and stop of each
block's modulation in that array; empty without one.
Raises
------
ValueError
If the cache cannot be loaded.
RuntimeError
With ``waveforms``, if a block's gradient holds another number of
samples than the segment position it plays at is prepared with, or
its wave more points than that position reserves.
"""
seq_path = Path(seq_path)
return require("play_cache")(
str(cache_path(seq_path, cache_ext)), seq_path.stat().st_size, waveforms
)
def _payload(
seq_path: Path,
system: pp.Opts,
verify_signature: bool,
fov_offset: Sequence[float] | None = None,
designed: list[tuple[Path, pp.Sequence]] | None = None,
) -> list[dict[str, Any]]:
"""Read the chain and return each file's libraries, in play order, prescribed to ``fov_offset``.
A file the reader refuses raises ``ValueError``, whatever the reader
itself raised.
"""
if designed is not None:
chain_read = designed
else:
try:
chain_read = read_chain(seq_path, verify=verify_signature)
except RuntimeError as failure:
raise ValueError(f"cannot read {seq_path}: {failure}") from failure
payload = []
for _, sequence in chain_read:
if fov_offset is not None:
prescribe(sequence, fov_offset)
payload.append(conversion_payload(sequence, system))
return payload
