.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/01-framework/05-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_05-custom-operator.py: ================== 5. Custom operator ================== In lesson 3 you wrote a sequence from the operators that ship with blochsim. In this lesson you write two operators of your own: a T2 preparation, and a readout that records two samples in one repetition. An operator is a Python function that returns events, and ``@`` composes operators into one. No kernel changes, and the new operators work in any simulator. **Learning objectives** - Write a preparation by composing the shipped operators with ``@``. - Play it from the layout of a simulator. - Write a readout that takes two samples per repetition. - Check an operator against a closed form or a shipped sequence. Previous: :doc:`04-description-based-simulation`. .. GENERATED FROM PYTHON SOURCE LINES 23-28 .. 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 28-39 .. GENERATED FROM PYTHON SOURCE LINES 40-52 Preparation ----------- A T2 preparation weights the magnetization by its transverse decay before the readout. It tips the magnetization into the transverse plane, lets it decay for the echo time TE_prep around a refocusing pulse, tips what is left back along z, and spoils whatever did not come back. Each of those steps is an operator that ships with blochsim, so ``@`` is all you write. The refocusing pulse is requested without crushers (``crushed=False``): this refocusing must form an echo, and the crusher pair that :func:`~blochsim.Refocusing` adds by default would spoil it. .. GENERATED FROM PYTHON SOURCE LINES 52-71 .. code-block:: Python import torch from blochsim import Delay, Excitation, Refocusing, Spoil def t2_preparation(echo_time_s, *, spoil_s=2e-3): """A T2 preparation of duration ``echo_time_s``, then a spoiler.""" 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 72-80 Preparation in a sequence ------------------------- An operator goes into a layout like any other. The preparation leaves the weighted magnetization along z, and the fast spin echo that follows excites it as it would any other longitudinal magnetization. The handlers are the defaults, so the layout only uses ``self.operators`` for the pulses and the readout of the train: .. GENERATED FROM PYTHON SOURCE LINES 80-102 .. code-block:: Python from blochsim import Simulator class T2PreparedFSE(Simulator): states = 8 def layout(self, *, TE_prep, ESP, ETL): 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 103-108 Contrast weighting ------------------ Sweep TE_prep and compare the first echo with exp(-TE_prep / T2), which is what a T2 preparation is for. Times are in milliseconds: .. GENERATED FROM PYTHON SOURCE LINES 108-142 .. code-block:: Python T2 = torch.tensor([40.0, 80.0, 160.0]) prep_times = torch.linspace(0.0, 120.0, 13) prepared = torch.stack( [ T2PreparedFSE(TE_prep=float(te), ESP=5.0, ETL=1) .simulate(T1=1000.0, T2=T2)[..., 0] .abs() for te in prep_times ], dim=-1, ) weighting = prepared / prepared[:, :1] expected = torch.exp(-prep_times / T2[:, None]) print(f"largest difference from exp(-TE/T2): {(weighting - expected).abs().max():.1e}") .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_05-custom-operator_001.png :alt: 05 custom operator :srcset: /generated/autoexamples/01-framework/images/sphx_glr_05-custom-operator_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none largest difference from exp(-TE/T2): 1.9e-03 .. GENERATED FROM PYTHON SOURCE LINES 143-162 The two agree to about 0.2%. The remainder is physics, not error: the magnetization tipped back up recovers a little during the spoiler that follows, and more of it does so after the longer preparations, which leave less behind. Readout ------- The shipped readouts differ in what they play around the sample. In an unbalanced steady-state train each repetition winds every configuration order on by one. A sample taken *before* that winding is a free induction decay (FID) after the pulse just played. A sample taken *after* it sits where the next pulse will refocus the previous excitation: an echo, and far more strongly T2-weighted. blochsim ships these as :func:`~blochsim.SSFPFidReadout` and :func:`~blochsim.SSFPEchoReadout`. Recording both in one repetition is a double-echo steady state (DESS), and you write it by putting the winding between two samples: .. GENERATED FROM PYTHON SOURCE LINES 162-170 .. code-block:: Python from blochsim import Dephase, Readout def dess_readout(phase_rad=0.0, *, duration_s=0.0): """Two samples, with the winding between them, then a wait.""" return Readout(phase_rad) @ Dephase() @ Readout(phase_rad) @ Delay(duration_s) .. GENERATED FROM PYTHON SOURCE LINES 171-174 Use it as the readout handler of a steady-state sequence. The sequence is one excitation and one readout per repetition, with the repetition time as the readout's wait: .. GENERATED FROM PYTHON SOURCE LINES 174-189 .. code-block:: Python FLIP, TR, REPETITIONS = 30.0, 20.0, 64 class DESS(Simulator): excitation = Excitation readout = dess_readout states = 24 def layout(self, *, flip, TR): return [ self.operators.excitation(torch.deg2rad(torch.as_tensor(flip))), self.operators.readout(duration_s=TR * 1e-3), ] .. GENERATED FROM PYTHON SOURCE LINES 190-194 Check it against the shipped readouts. The first sample has to be what an FID train records and the second what an echo train records, since those are the same two samples taken one at a time. Changing the readout handler is all that separates the three simulators: .. GENERATED FROM PYTHON SOURCE LINES 194-219 .. code-block:: Python from blochsim import SSFPEchoReadout, SSFPFidReadout class SSFPFid(DESS): readout = SSFPFidReadout class SSFPEcho(DESS): readout = SSFPEchoReadout def played(sequence, T1=1000.0): train = sequence(flip=FLIP, TR=TR, repetitions=REPETITIONS) return train.simulate(T1=T1, T2=T2) both = played(DESS) fid, echo = both[..., 0::2], both[..., 1::2] print( f"first sample against SSFPFidReadout: {(fid - played(SSFPFid)).abs().max():.1e}" ) print( f"second sample against SSFPEchoReadout: {(echo - played(SSFPEcho)).abs().max():.1e}" ) .. 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 220-225 Both samples agree exactly with the corresponding single-sample sequences. The ratio of the two is a T2 contrast. The echo has spent one more repetition in the transverse plane, so it carries T2 where the FID carries a mixture of T1 and T2. Compute it at three T1 values: .. GENERATED FROM PYTHON SOURCE LINES 225-249 .. code-block:: Python ratios = {} for t1 in (600.0, 1000.0, 2000.0): recorded = played(DESS, t1) ratios[t1] = recorded[..., 1::2][:, -1].abs() / recorded[..., 0::2][:, -1].abs() .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_05-custom-operator_002.png :alt: two samples, T1 = 1000 ms, ratio :srcset: /generated/autoexamples/01-framework/images/sphx_glr_05-custom-operator_002.png :class: sphx-glr-single-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 250-281 The ratio rises with T2 at every T1 and moves much less with T1 than with T2. That makes it usable as a T2 contrast. It also means that a DESS T2 measurement at a larger flip angle needs T1 to be known. Limits ------ A preparation, a readout, or a shaped or per-channel pulse is built from the shipped operators and reaches the kernels unchanged. An event cannot say *how much* a gradient dephases. A sequence carries one crusher moment, and dephasing comes in whole configuration orders. A bipolar pair, a velocity-encoding moment of its own, or a crusher of twice the area of its neighbour cannot be written as operators. They need a gradient moment per event through the packed layout and every kernel, which is a change to blochsim itself. Specification ------------- What this lesson built, stated the way you would ask an agent for it: .. code-block:: text Write two blochsim operators by composition with @. A T2 preparation: 90 degree pulse, wait TE/2, uncrushed 180 degree refocusing, wait TE/2, -90 degree pulse, short wait, spoiler. Check the first echo of a spin echo train that follows it against exp(-TE/T2). A DESS readout: sample, dephase, sample, wait. Use it as the readout handler of a steady-state Simulator (30 degrees, TR 20 ms) and check its two samples against the shipped SSFPFidReadout and SSFPEchoReadout. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.385 seconds) .. _sphx_glr_download_generated_autoexamples_01-framework_05-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/pulserver/blochsim/gh-pages?urlpath=lab/tree/latest/notebooks/generated/autoexamples/01-framework/05-custom-operator.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 05-custom-operator.ipynb <05-custom-operator.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 05-custom-operator.py <05-custom-operator.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 05-custom-operator.zip <05-custom-operator.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_