.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/01-framework/04-custom-operator.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_04-custom-operator.py: ====================== Writing a New Operator ====================== The scope of this notebook is to show how to add a sequence module TorchSim does not ship -- a preparation, or a readout -- without touching a kernel. An operator is a Python function that returns events and says how long it holds the timeline, and ``@`` composes two into one. Two are written here: a T2 preparation, and a readout that takes both samples an unbalanced repetition can carry. .. GENERATED FROM PYTHON SOURCE LINES 16-20 .. colab-link:: :needs_gpu: 0 !pip install torchsim .. GENERATED FROM PYTHON SOURCE LINES 22-25 The operators the new ones are composed from, and the simulator that plays them. .. GENERATED FROM PYTHON SOURCE LINES 26-115 .. code-block:: Python import torch from torchsim import ( Delay, Dephase, Excitation, Readout, Refocusing, SSFPEchoReadout, SSFPFidReadout, Spoil, ) from torchsim.model import Simulator .. GENERATED FROM PYTHON SOURCE LINES 116-129 Composing existing operators ---------------------------- A T2 preparation tips the magnetization into the transverse plane, lets it decay for a chosen time about a Refocusing pulse, tips what is left back along z, and spoils whatever did not come back. All of those are operators already, so ``@`` is the whole of writing it. Nothing new is taught to the kernels -- what is new is the *arrangement*, and that is what an operator is. The Refocusing pulse is asked for uncrushed: a T2 preparation refocuses rather than dephases, and the crusher pair :func:`~torchsim.Refocusing` adds by default would spoil the echo it exists to form. .. GENERATED FROM PYTHON SOURCE LINES 129-153 .. code-block:: Python def t2_preparation(echo_time_s, *, spoil_s=2e-3): """Return a T2 preparation that weights the magnetization by its own decay. Parameters ---------- echo_time_s: How long the magnetization spends in the transverse plane. spoil_s: The spoiler after the tip-up, which removes what did not return. """ half = 0.5 * echo_time_s return ( Excitation(0.5 * torch.pi) @ Delay(half) @ Refocusing(torch.pi, 0.5 * torch.pi, crushed=False) @ Delay(half) @ Excitation(-0.5 * torch.pi) @ Delay(spoil_s) @ Spoil() ) .. GENERATED FROM PYTHON SOURCE LINES 154-159 Using the operator ------------------ A new operator goes into a layout beside the shipped ones. The preparation leaves the weighted magnetization along z, so the train that follows excites it as it would any other longitudinal magnetization. .. GENERATED FROM PYTHON SOURCE LINES 159-194 .. code-block:: Python class T2PreparedFSE(Simulator): """A T2 preparation, then a refocused train. Parameters ---------- TE_prep : float How long the preparation holds the magnetization transverse, in milliseconds. ESP : float The echo spacing, in milliseconds. ETL : int How many echoes are recorded. """ states = 8 def layout(self, *, TE_prep, ESP, ETL): """Return the operators of the whole train, in order.""" half = Delay(0.5 * ESP * 1e-3) parts = [ t2_preparation(TE_prep * 1e-3), self.operators.excitation(0.5 * torch.pi, 0.5 * torch.pi), ] for _ in range(ETL): parts += [ half, self.operators.refocusing(torch.pi, 0.5 * torch.pi), half, self.operators.readout(0.5 * torch.pi), ] return parts .. GENERATED FROM PYTHON SOURCE LINES 195-202 Validation ---------- A preparation is worth only as much as the weighting it imposes, so we check it rather than assert it: sweep the preparation time and hold the first recorded echo against ``exp(-TE / T2)``, which is what a T2 preparation is for. .. GENERATED FROM PYTHON SOURCE LINES 202-236 .. code-block:: Python T2_MS = torch.tensor([40.0, 80.0, 160.0]) prep_times_ms = torch.linspace(0.0, 120.0, 13) prepared = torch.stack( [ T2PreparedFSE(TE_prep=float(prep_ms), ESP=5.0, ETL=1) .simulate(T1=1000.0, T2=T2_MS)[..., 0] .abs() for prep_ms in prep_times_ms ], dim=-1, ) weighting = prepared / prepared[:, :1] expected = torch.exp(-prep_times_ms / T2_MS[:, None]) .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_04-custom-operator_001.png :alt: simulated, against the closed form (solid) :srcset: /generated/autoexamples/01-framework/images/sphx_glr_04-custom-operator_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none worst departure from exp(-TE/T2): 0.0018985345959663391 .. GENERATED FROM PYTHON SOURCE LINES 237-240 The two follow each other to about 0.2%, and the residual is physics rather than error: the tipped-up magnetization recovers a little across the spoiler that follows it, by more for the longer preparations that leave less behind. .. GENERATED FROM PYTHON SOURCE LINES 242-253 Custom readout -------------- The shipped readouts differ only in what they play around the sample. An unbalanced train winds every order on once per repetition, so a sample taken *before* that winding is a free induction decay after the pulse just played, and a sample taken *after* it sits where the next pulse will refocus the previous excitation -- an echo, and far more strongly T2-weighted. TorchSim ships each of those separately. Taking both in one repetition is a double-echo steady state, and writing it is putting the winding between two samples rather than on one side of them. .. GENERATED FROM PYTHON SOURCE LINES 253-268 .. code-block:: Python def dess_readout(phase_rad=0.0, *, duration_s=0.0): """Return the two samples an unbalanced repetition can carry. Parameters ---------- phase_rad : float, optional The receiver phase both samples are taken at. duration_s : float, optional What is left of the repetition after the second sample. """ return Readout(phase_rad) @ Dephase() @ Readout(phase_rad) @ Delay(duration_s) .. GENERATED FROM PYTHON SOURCE LINES 269-274 Whether that is the right arrangement is not a matter of opinion: the first sample has to be what an SSFP-FID train records and the second what an SSFP-Echo train records, since those are the same two samples taken one at a time. So the check is to run all three. .. GENERATED FROM PYTHON SOURCE LINES 275-321 .. code-block:: Python FLIP_DEG, TR_MS, REPETITIONS = 30.0, 20.0, 64 T2_MS = torch.tensor([40.0, 80.0, 160.0]) class DESS(Simulator): """A steady-state train taking both samples each repetition can carry.""" excitation = Excitation readout = dess_readout states = 24 def layout(self, *, flip, TR): """Return the operators of one repetition, in order.""" return [ self.operators.excitation(torch.deg2rad(torch.as_tensor(flip))), self.operators.readout(duration_s=TR * 1e-3), ] # Naming a different readout is the whole of the difference between the three. class SSFPFid(DESS): readout = SSFPFidReadout class SSFPEcho(DESS): readout = SSFPEchoReadout def played(sequence, t1_ms=1000.0): """Return what one train records, over the three T2 values.""" train = sequence(flip=FLIP_DEG, TR=TR_MS, repetitions=REPETITIONS) return train.simulate(T1=t1_ms, T2=T2_MS) both = played(DESS) fid, echo = both[..., 0::2], both[..., 1::2] .. rst-class:: sphx-glr-script-out .. code-block:: none first sample against SSFPFidReadout: 0.0e+00 second sample against SSFPEchoReadout: 0.0e+00 .. GENERATED FROM PYTHON SOURCE LINES 322-330 Both exactly, which is the whole claim: two samples in one repetition, and each is the sample the sequence that takes it alone would have recorded. What it is for is the ratio between them. The echo has spent a further repetition in the transverse plane, so it carries T2 where the free induction decay carries a mixture -- and the ratio of the two is a T2 contrast that needs no separate measurement to normalize. .. GENERATED FROM PYTHON SOURCE LINES 331-366 .. code-block:: Python ratios = {} for t1_ms in (600.0, 1000.0, 2000.0): recorded = played(DESS, t1_ms) ratios[t1_ms] = recorded[..., 1::2][:, -1].abs() / recorded[..., 0::2][:, -1].abs() .. rst-class:: sphx-glr-horizontal * .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_04-custom-operator_002.png :alt: the two samples one repetition takes (echo dashed) :srcset: /generated/autoexamples/01-framework/images/sphx_glr_04-custom-operator_002.png :class: sphx-glr-multi-img * .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_04-custom-operator_003.png :alt: their ratio rises with T2, and moves far less with T1 :srcset: /generated/autoexamples/01-framework/images/sphx_glr_04-custom-operator_003.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-script-out .. code-block:: none at T2 = 40 ms the ratio moves 0.10 over a 3.3x range in T1 .. GENERATED FROM PYTHON SOURCE LINES 367-370 The ratio rises with T2 at every T1, and moves far less with T1 than with T2 -- which is what makes it usable, and why a DESS T2 measurement at a larger flip angle wants T1 known rather than assumed away. .. GENERATED FROM PYTHON SOURCE LINES 373-385 Limits ------ A preparation, a readout, a shaped or per-channel pulse is written from the shipped operators and reaches the kernels unchanged. What an event stream cannot express is *how much* a gradient dephases. It carries one crusher moment for the whole sequence and dephasing is quantized to whole configuration orders, so a bipolar pair, a velocity-encoding moment of its own, or a crusher of twice its neighbour's area have no spelling here. Those need a per-event gradient moment through the packed layout and every kernel, which is a change to the engine rather than to an operator written on top of it. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.464 seconds) .. _sphx_glr_download_generated_autoexamples_01-framework_04-custom-operator.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/firmlab-pisa/torchsim/gh-pages?urlpath=lab/tree/v0.0.8/examples/generated/autoexamples/01-framework/04-custom-operator.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04-custom-operator.ipynb <04-custom-operator.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04-custom-operator.py <04-custom-operator.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04-custom-operator.zip <04-custom-operator.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_