.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/01-framework/03-custom-signal-model.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code or to run this example in your browser via Binder. .. rst-class:: sphx-glr-example-title .. _sphx_glr_generated_autoexamples_01-framework_03-custom-signal-model.py: ======================= 3. Custom signal model ======================= Lessons 1 and 2 used shipped simulators. This lesson implements a new one, saturation recovery, by subclassing :class:`~blochsim.Simulator`. The resulting class supports the same features as a shipped simulator: batched tissue inputs, derivatives with respect to tissue and sequence parameters, and the additional physics of lesson 2. **Learning objectives** - Define the handlers, which map each kind of event to an operator. - Define the layout, the ordered events of one repetition. - Validate the model against its closed-form solution. - Inspect the generated sequence description. - Differentiate with respect to tissue and sequence parameters. Previous: :doc:`02-advanced-physics`. Next: :doc:`04-description-based-simulation`, where the events are read from a Pulseq file or a scanner instead of the layout. .. GENERATED FROM PYTHON SOURCE LINES 26-31 .. colab-link:: :needs_gpu: 0 !pip install blochsim matplotlib !wget --quiet --no-clobber https://raw.githubusercontent.com/pulserver/blochsim/main/docs/figure_style.py .. GENERATED FROM PYTHON SOURCE LINES 31-43 .. GENERATED FROM PYTHON SOURCE LINES 44-71 Sequence -------- Each block of the sequence is a 90-degree pulse followed by a spoiler, which nulls the magnetization; a saturation time TS of free recovery; and a low-flip-angle spoiled readout. Over a set of saturation times the signal is .. math:: S(T_S) = M_0 \sin\alpha \, \left(1 - e^{-T_S/T_1}\right). which serves as the reference for the simulation. Handlers -------- A sequence consists of RF, ADC and wait events. The *handlers* define the operator that plays each kind of event. They are class attributes of the simulator, one per kind: ``excitation``, ``refocusing``, ``inversion``, ``saturation``, ``readout`` and ``delay`` (free precession with no RF and no ADC). Unset handlers keep their defaults. Saturation recovery uses four. The saturation handler is a 90-degree pulse composed with a spoiler through ``@``. The readout is :func:`~blochsim.SPGRReadout`, a sample followed by ideal transverse spoiling: .. GENERATED FROM PYTHON SOURCE LINES 71-90 .. code-block:: Python import torch from blochsim import Delay, Excitation, SPGRReadout, Spoil from blochsim import Simulator def Saturate(flip_rad=torch.pi / 2, phase_rad=0.0): """A 90-degree pulse, then a spoiler.""" return Excitation(flip_rad, phase_rad) @ Spoil() class SaturationRecovery(Simulator): saturation = Saturate # plays a saturation pulse excitation = Excitation # plays an excitation pulse delay = Delay # plays a wait readout = SPGRReadout # plays a readout, then spoils states = 1 .. GENERATED FROM PYTHON SOURCE LINES 91-112 Since every block starts and ends with spoiling, no transverse coherence survives between blocks, and a single configuration state (``states = 1``) is exact. The handlers also determine how an external description is simulated. A Pulseq file or a scanner stream specifies RF pulses, ADC windows and timing, but blochsim does not take gradients from it. When such a description is simulated (lesson 4), the handlers define the dephasing between events: spoiled, unbalanced or balanced. Layout ------ The *layout* returns the events of one repetition in playing order, built from the handlers through ``self.operators``. It is the default sequence description, generated by the constructor from the sequence parameters. Timestamps are derived from the duration of each operator. The sequence parameters are the keyword arguments of ``layout``: here the saturation times in milliseconds and the readout flip angle in degrees. Operators take SI units (seconds, radians). .. GENERATED FROM PYTHON SOURCE LINES 112-124 .. code-block:: Python class SaturationRecovery(SaturationRecovery): def layout(self, *, TS, flip): angle = torch.deg2rad(torch.as_tensor(flip)) parts = [] for wait in torch.as_tensor(TS) * 1e-3: # ms to s parts += [self.operators.saturation(), self.operators.delay(wait)] parts += [self.operators.excitation(angle), self.operators.readout()] return parts .. GENERATED FROM PYTHON SOURCE LINES 125-129 Simulation ---------- The class is used like a shipped simulator: .. GENERATED FROM PYTHON SOURCE LINES 129-151 .. code-block:: Python TS = torch.logspace(1.3, 3.7, 40) # 20 ms to 5 s T1 = torch.tensor([830.0, 1330.0, 4000.0]) recovery = SaturationRecovery(TS=TS, flip=10.0) signal = recovery.simulate(T1=T1, T2=torch.tensor([80.0, 110.0, 2000.0])) closed_form = torch.sin(torch.deg2rad(torch.tensor(10.0))) * ( 1 - torch.exp(-TS / T1[:, None]) ) print(f"largest difference: {(signal.abs() - closed_form).abs().max():.1e}") .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_03-custom-signal-model_001.png :alt: 03 custom signal model :srcset: /generated/autoexamples/01-framework/images/sphx_glr_03-custom-signal-model_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none largest difference: 3.0e-08 .. GENERATED FROM PYTHON SOURCE LINES 152-161 The simulation matches the closed form to single precision. Since blochsim evaluates the events individually, the agreement validates the layout. Sequence description -------------------- ``describe`` returns the events your layout produced, with their timestamps, as a :class:`~blochsim.SequenceDescription`; ``plot`` displays them. This is the first thing to inspect when a result is unexpected: .. GENERATED FROM PYTHON SOURCE LINES 161-171 .. code-block:: Python description = recovery.describe(TS=[50.0, 100.0, 200.0], flip=10.0) .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_03-custom-signal-model_002.png :alt: 03 custom signal model :srcset: /generated/autoexamples/01-framework/images/sphx_glr_03-custom-signal-model_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 172-181 Three blocks, each with a saturation pulse, a wait of increasing length, a 10-degree excitation and a readout. Lesson 4 reads the same kind of object from a Pulseq file or a scanner. Derivatives ----------- ``jacobian`` works on the new class without additional code. The derivative with respect to T1 is what a T1 fit requires: .. GENERATED FROM PYTHON SOURCE LINES 181-183 .. code-block:: Python signal, dT1 = recovery.jacobian("T1", T1=T1, T2=80.0) .. GENERATED FROM PYTHON SOURCE LINES 184-187 Sequence parameters are PyTorch tensors, so derivatives with respect to them are obtained by reverse-mode autodiff (``backward``). Here, the gradient of the summed signal with respect to each saturation time: .. GENERATED FROM PYTHON SOURCE LINES 187-201 .. code-block:: Python TS_designed = TS.clone().requires_grad_(True) total = recovery.simulate(TS=TS_designed, T1=830.0, T2=80.0).abs().sum() total.backward() .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_03-custom-signal-model_003.png :alt: 03 custom signal model :srcset: /generated/autoexamples/01-framework/images/sphx_glr_03-custom-signal-model_003.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 202-211 This gradient is the basis of sequence optimization; the sequence design Applications use it on complete protocols. Additional physics ------------------ The physics of lesson 2 is controlled by tissue properties and therefore applies to the new class unchanged. With B1 = 0.8 the saturation pulse is 72 degrees, and saturation is incomplete: .. GENERATED FROM PYTHON SOURCE LINES 211-223 .. code-block:: Python low_b1 = recovery.simulate(T1=830.0, T2=80.0, B1=0.8) ideal = recovery.simulate(T1=830.0, T2=80.0) .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_03-custom-signal-model_004.png :alt: 03 custom signal model :srcset: /generated/autoexamples/01-framework/images/sphx_glr_03-custom-signal-model_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 224-241 At short saturation times the residual longitudinal magnetization raises the signal; at long ones the 8-degree readout lowers it. A T1 fit that ignores B1 is biased by both effects. Specification ------------- What this lesson built, stated the way you would ask an agent for it: .. code-block:: text Write a blochsim Simulator for saturation recovery. Handlers: a saturation that is a 90-degree pulse followed by a spoiler, an ideal excitation, a delay, and a spoiled readout; keep one configuration state. The layout takes the saturation times TS (ms) and the readout flip angle (deg), and plays saturate, wait TS, excite, read for each TS. Check it against M0 sin(flip) (1 - exp(-TS/T1)) in white matter, grey matter and CSF. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.617 seconds) .. _sphx_glr_download_generated_autoexamples_01-framework_03-custom-signal-model.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: binder-badge .. image:: images/binder_badge_logo.svg :target: https://mybinder.org/v2/gh/pulserver/blochsim/gh-pages?urlpath=lab/tree/latest/notebooks/generated/autoexamples/01-framework/03-custom-signal-model.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 03-custom-signal-model.ipynb <03-custom-signal-model.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 03-custom-signal-model.py <03-custom-signal-model.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 03-custom-signal-model.zip <03-custom-signal-model.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_