Read the sequence description from ISMRMRD waveforms; description's RF raster by default; 20 mm non-selective rule

docs/developer-guide/internals/description.md

@@ -154,20 +154,31 @@ echo flag is set.
 
 Pulserver's reconstruction proxy reads the sequence file a series was played
 from and sends its description after the MRD XML header and before the first
-acquisition, as one `TEXT` message: a JSON object under the key
-`pulserver_sequence_description`, holding one entry per file of the sequence
-chain. {func}`~blochsim.read_mrd_description` reads the stream up to that
-message, returns the XML header and the descriptions in play order, and leaves
-the stream at the first acquisition. A message other than a filename, config,
-header or text before the description is refused.
-
-Each entry holds `subsequence_index`, `tr_duration_us`, `rf_raster_time_s`,
-the events as three base64 arrays -- `type` (int32), `timestamp_us`
-(float64) and `params` (float32, seven per event, in the wire order above) --
-and `rf_definitions`, each with its bandwidth, bands, and Pulseq-compressed
-magnitude, phase and time shapes. The stream carries no shim definitions, no
-per-sample slice-select gradient, no crusher dephasing or voxel size, and no
-event actions.
+acquisition, as ISMRMRD waveforms, which every MRD writer keeps with the series:
+per file of the sequence chain a header (`waveform_id` 1024) and an event list
+(1025), and one waveform per RF definition (1026). Each carries one vector of
+32-bit words led by its length, over as many channels as keep a channel within
+65 535 samples; floats are bit-cast and float64 values take two words, low then
+high.
+
+- Header: format version, subsequence index, number of subsequences, number of
+  events, number of RF definitions, `tr_duration_us`, `rf_raster_time_s`.
+- Event list: subsequence index, number of events, then ten words an event:
+  type, `timestamp_us` (float64), and the seven float32 parameters in the wire
+  order above.
+- RF definition: subsequence index, id, bandwidth, number of bands, eight band
+  offsets, band bandwidth, a word left 0, then for the magnitude, phase and time
+  shapes whether each is given and its uncompressed and stored sample counts,
+  then the samples of those given.
+
+{func}`~blochsim.read_mrd_description` reads the stream up to the last of
+them, passing the physiological waveforms over, returns the XML header and the
+descriptions in play order, and leaves the stream at the first acquisition;
+{func}`~blochsim.description_from_waveforms` reads them off the waveforms of a
+stored series. A message other than a filename, config, header, text or
+waveform before the description is refused. The description
+carries no shim definitions, no per-sample slice-select gradient, no crusher
+dephasing or voxel size, and no event actions.
 
 ## Gradients
 

docs/explanations/description.md

@@ -56,7 +56,7 @@ wire format are in {doc}`../developer-guide/internals/description`.
 - **Decodes the description pulserver sends on an MRD stream, one per file of the sequence chain, in play order.**
 
   :::{dropdown} Show code and tests
-  Code: {func}`~blochsim.read_mrd_description`. Tests: *a written stream reads back into the rows it was built from*; *every file of a chain is read in play order*; *a stream whose data arrives first is refused* (`test_mrd_reader.py`).
+  Code: {func}`~blochsim.read_mrd_description`, {func}`~blochsim.description_from_waveforms`. Tests: *a written stream reads back into the rows it was built from*; *every file of a chain is read in play order*; *a stream whose data arrives first is refused*; *the physiological waveforms ahead of the description are passed over*; *a kept series reads its description from its waveforms*; *a pulse is integrated on the raster its description states* (`test_mrd_reader.py`).
   :::
 - **Reproduces the spoiled gradient-echo steady state from a description read either way.**
 

docs/explanations/epg.md

@@ -32,10 +32,10 @@ kernels hold. The derivations and the operator-by-operator detail are in
   :::{dropdown} Show code and tests
   Code: `transmit_field`, {class}`~blochsim.ShimDefinition`. Tests: *channels in antiphase cancel*; *the field is the complex sum written out* (`test_transmit_array.py`).
   :::
-- **Integrates a shaped pulse into the rotation it performs across the slice, or per voxel when its channel weights vary while it plays.**
+- **Integrates a shaped pulse into the rotation it performs across the slice, or per voxel when its channel weights vary while it plays; a pulse selecting more than `NONSELECTIVE_THICKNESS_M` (20 mm) is taken at its centre rotation, a slab the encoding resolves.**
 
   :::{dropdown} Show code and tests
-  Code: `exact_slice_profile`, `transition_table`, `dynamic_pair`. Tests: *the slice shows up when the pulse selects one* (`test_exact_profile.py`); *the pair is what exponentiating each sample gives* (`test_dynamic_transmit.py`).
+  Code: `exact_slice_profile`, `transition_table`, `dynamic_pair`. Tests: *the slice shows up when the pulse selects one*; *a slab thicker than the threshold is profiled as its centre* (`test_exact_profile.py`); *the pair is what exponentiating each sample gives* (`test_dynamic_transmit.py`).
   :::
 - **Plays a waveform sample by sample, with relaxation, precession and exchange between samples.**
 

src/blochsim/sequence/__init__.py

@@ -5,6 +5,7 @@
 __all__ = [
     "SequenceDescription",
     "ShimDefinition",
+    "description_from_waveforms",
     "read_mrd_description",
     "FSEReadout",
     "Operator",
@@ -53,7 +54,7 @@
     ideal_rf_definition,  # noqa: F401
     rf_definition,
 )
-from ._mrd import read_mrd_description
+from ._mrd import description_from_waveforms, read_mrd_description
 from ._operators import (
     Delay,
     Dephase,

src/blochsim/sequence/_accelerators.py

@@ -5,6 +5,7 @@
 __all__: list[str] = []
 
 import functools
+import math
 import os
 import weakref
 from collections import OrderedDict
@@ -288,6 +289,40 @@ def _pack_for(
     )
 
 
+#: Thickness in m above which a pulse is taken as non-selective: a slab this
+#: thick is resolved by the encoding, so each position of it sees the pulse's
+#: centre rotation, not a profile across it.
+NONSELECTIVE_THICKNESS_M = 20e-3
+
+
+def _selected_thickness_m(description: SequenceDescription) -> dict[int, float]:
+    """The thinnest slice each RF definition is played over, in m; 0 where unknown.
+
+    A definition's ``bandwidth_hz`` over the event's slice-select gradient. A
+    description that states no gradient or no bandwidth states no thickness,
+    and its pulse keeps the mode its waveform asks for.
+    """
+    thickness: dict[int, float] = {}
+    for event in description.events:
+        if event.type is not EventType.RF:
+            continue
+        key = event.rf_definition_id
+        bandwidth = float(description.rf_definitions[key].bandwidth_hz)
+        gradient = abs(float(event.slice_select_gradient_hz_per_m))
+        here = bandwidth / gradient if bandwidth > 0.0 and gradient > 0.0 else 0.0
+        thickness[key] = min(thickness.get(key, math.inf), here)
+    return thickness
+
+
+def _selective(description: SequenceDescription, key: int) -> bool:
+    """Whether a definition selects a slice thin enough to profile across, or states none."""
+    cached = getattr(description, "_thickness_cache", None)
+    if cached is None:
+        cached = _selected_thickness_m(description)
+        object.__setattr__(description, "_thickness_cache", cached)
+    return cached.get(key, math.inf) <= NONSELECTIVE_THICKNESS_M
+
+
 def _wants_a_table(description: SequenceDescription, slice_profile: Any) -> bool:
     """Whether the pulses are to be integrated into a rotation over the slice.
 
@@ -311,6 +346,7 @@ def _wants_a_table(description: SequenceDescription, slice_profile: Any) -> bool
         }
         cached = any(
             description.rf_definitions[key].rf_mode() is RfMode.PROFILED
+            and _selective(description, key)
             for key in driven
         )
         object.__setattr__(description, "_profiled_cache", cached)
@@ -327,7 +363,9 @@ def _pulse_shapes(
     different shapes get one each, and the event carries which it reads.
 
     Ordered by first use, so the row a definition takes does not depend on how
-    the description happened to number it.
+    the description happened to number it. A pulse selecting more than
+    :data:`NONSELECTIVE_THICKNESS_M` is tabulated with no bandwidth, so every
+    position reads its centre rotation.
 
     Raises
     ------
@@ -339,7 +377,10 @@ def _pulse_shapes(
             continue
         key = event.rf_definition_id
         if key not in used:
-            used[key] = description.rf_definitions[key]
+            definition = description.rf_definitions[key]
+            if not _selective(description, key):
+                definition = replace(definition, bandwidth_hz=0.0)
+            used[key] = definition
     if not used:
         raise ValueError(
             "a slice profile is the rotation a pulse performs across the "