.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/01-framework/01-first-simulation.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_01-first-simulation.py: ======================== 1. Your first simulation ======================== blochsim computes the signal of an MR sequence for given tissue properties, together with its derivatives with respect to those properties and to the sequence parameters. This lesson uses the fast spin echo (FSE) simulator shipped with blochsim: you simulate three tissues, change the refocusing angle, and compute the derivative of the signal with respect to T2. **Learning objectives** - Construct a shipped simulator and evaluate it for several tissues. - Take the derivative of the signal with respect to a tissue property. - Choose how many configuration states to keep, and where the simulation runs. Next: :doc:`02-advanced-physics`, where the same sequence meets B1, off resonance, magnetization transfer, diffusion and flow. .. 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-40 .. GENERATED FROM PYTHON SOURCE LINES 41-51 Sequence -------- The sequence is a CPMG echo train: a 90-degree excitation followed by 48 refocusing pulses, echo spacing (ESP) 5 ms, repetition time 3 s. Each shipped sequence is a *simulator* class, constructed from its sequence parameters. Times are in milliseconds and angles in degrees. The refocusing angles are given per echo, so variable flip angle trains use the same interface. .. GENERATED FROM PYTHON SOURCE LINES 51-58 .. code-block:: Python import torch from blochsim.simulators import FSESimulator fse = FSESimulator(ESP=5.0, TR=3000.0) flip = torch.full((48,), 180.0) .. GENERATED FROM PYTHON SOURCE LINES 59-67 Tissue ------ Tissue properties are passed at call time; here T1 and T2, in milliseconds. Each simulator accepts scalar and tensor-valued inputs. Tensor-valued inputs are simulated in parallel, over multiple CPU threads or on a GPU if available. The three entries below approximate white matter, grey matter and cerebrospinal fluid, with a common T1 of 1 s: .. GENERATED FROM PYTHON SOURCE LINES 67-81 .. code-block:: Python T2 = torch.tensor([80.0, 110.0, 2000.0]) signal = fse.simulate(flip=flip, T1=1000.0, T2=T2) print(signal.shape) .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_01-first-simulation_001.png :alt: 01 first simulation :srcset: /generated/autoexamples/01-framework/images/sphx_glr_01-first-simulation_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none torch.Size([3, 48]) .. GENERATED FROM PYTHON SOURCE LINES 82-92 The output has one row per tissue and one column per echo. With ideal 180-degree refocusing each train decays as exp(-TE/T2). Sequence parameters ------------------- Clinical FSE protocols often use refocusing angles below 180 degrees to reduce SAR, which scales with the square of the flip angle. At 60 degrees more coherence pathways contribute, so more configuration states are retained (``states``, discussed below): .. GENERATED FROM PYTHON SOURCE LINES 92-108 .. code-block:: Python signal60 = fse.simulate(flip=torch.full((48,), 60.0), T1=1000.0, T2=T2, states=48) .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_01-first-simulation_002.png :alt: 180° refocusing, 60° refocusing :srcset: /generated/autoexamples/01-framework/images/sphx_glr_01-first-simulation_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 109-121 With 60-degree refocusing, part of the magnetization is stored along z and returns as stimulated echoes. The first echoes oscillate before the train approaches a pseudo steady state, and the late echoes decay more slowly than exp(-TE/T2). blochsim models all coherence pathways with the extended phase graph (EPG) formalism; see the :doc:`EPG explanation `. Tissue derivatives ------------------ Model-based T2 fitting, as well as Cramér-Rao bound analysis, needs the derivative of the signal with respect to T2. ``jacobian`` returns the signal and its derivatives with respect to the named properties: .. GENERATED FROM PYTHON SOURCE LINES 121-123 .. code-block:: Python signal, dT2 = fse.jacobian("T2", flip=flip, T1=1000.0, T2=T2) .. GENERATED FROM PYTHON SOURCE LINES 124-127 The derivative is computed in forward mode alongside the signal, not by finite differences. As a check, compare it with a forward difference over a 1 ms step: .. GENERATED FROM PYTHON SOURCE LINES 127-144 .. code-block:: Python step = 1.0 finite = (fse.simulate(flip=flip, T1=1000.0, T2=T2 + step) - signal) / step error = (dT2 - finite).abs().max() / dT2.abs().max() print(f"largest difference: {error:.1e}") .. image-sg:: /generated/autoexamples/01-framework/images/sphx_glr_01-first-simulation_003.png :alt: 01 first simulation :srcset: /generated/autoexamples/01-framework/images/sphx_glr_01-first-simulation_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none largest difference: 7.7e-03 .. GENERATED FROM PYTHON SOURCE LINES 145-160 The two agree to within the truncation error of the finite difference. The sensitivity to T2 peaks at TE close to T2. Passing several names, ``jacobian(("T1", "T2"), ...)``, returns one derivative per property. Accuracy and devices -------------------- Two settings control accuracy and cost. Both can be passed to the constructor or to the call. ``states`` is the number of EPG configuration states retained. A train of N refocusing pulses populates at most N states, so 48 is exact here. The FSE simulator defaults to 10, which is sufficient for near-180-degree trains. Truncation is a systematic error, so check it against the exact number when the refocusing angle is low: .. GENERATED FROM PYTHON SOURCE LINES 160-167 .. code-block:: Python flip60 = torch.full((48,), 60.0) reference = fse.simulate(flip=flip60, T1=1000.0, T2=T2, states=48) for states in (10, 24, 32): approximate = fse.simulate(flip=flip60, T1=1000.0, T2=T2, states=states) error = (approximate - reference).abs().max() / reference.abs().max() print(f"{states:2d} states: {error:.0%} off") .. rst-class:: sphx-glr-script-out .. code-block:: none 10 states: 8% off 24 states: 6% off 32 states: 6% off .. GENERATED FROM PYTHON SOURCE LINES 168-175 The 60-degree train requires the full 48 states, which the previous figure used. ``execution=`` selects the device. By default (``"auto"``) a small problem runs where its tensors are, and a large one is cut into chunks of voxels spread over every GPU, each sized to the card's free memory; the signal comes back where the tensors are. ``"cpu"`` or a named device overrides this: .. GENERATED FROM PYTHON SOURCE LINES 175-177 .. code-block:: Python on_cpu = fse.simulate(flip=flip, T1=1000.0, T2=T2, execution="cpu") .. GENERATED FROM PYTHON SOURCE LINES 178-183 The first call resolves the sequence structure; later calls only rebind numerical values. Construct a simulator once and reuse it inside loops. Each shipped sequence also has a functional form, convenient for one-off calls: .. GENERATED FROM PYTHON SOURCE LINES 183-187 .. code-block:: Python from blochsim import fse_sim signal, dT2 = fse_sim(flip=flip, ESP=5.0, TR=3000.0, T1=1000.0, T2=T2, diff="T2") .. GENERATED FROM PYTHON SOURCE LINES 188-201 Specification ------------- What this lesson did, stated the way you would ask an agent for it: .. code-block:: text With blochsim, simulate a fast spin echo of 48 echoes, 5 ms apart, TR 3 s, in white matter, grey matter and CSF (T1 1 s; T2 80, 110 and 2000 ms). Compare 180 and 60 degree refocusing trains. Return the derivative of the signal with respect to T2, checked against a finite difference, and check that the number of configuration states kept is enough for the 60 degree train. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 1.709 seconds) .. _sphx_glr_download_generated_autoexamples_01-framework_01-first-simulation.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/01-first-simulation.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01-first-simulation.ipynb <01-first-simulation.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 01-first-simulation.py <01-first-simulation.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 01-first-simulation.zip <01-first-simulation.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_