.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/02-parameter-inference/03-nonlinear-least-squares.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_02-parameter-inference_03-nonlinear-least-squares.py: ===================================== T2 mapping by nonlinear least squares ===================================== The scope of this notebook is to map T2 from a multi-echo spin echo by nonlinear least squares, and to compare it against a dictionary match on the same slice. What decides the comparison is a nuisance parameter. A magnitude reconstruction sits on a noise floor, so the decay does not reach zero, and unlike the proton density that offset does not divide out of a normalized match. Ignore it and T2 is biased; put it on the grid and the grid multiplies. A fit pays one more column of the Jacobian instead, and gives up the guarantee that it found the global minimum. .. GENERATED FROM PYTHON SOURCE LINES 19-23 .. colab-link:: :needs_gpu: 0 !pip install torchsim brainweb-dl cmap .. GENERATED FROM PYTHON SOURCE LINES 25-29 The problem is stated over a simulator carrying the sequence and filled in by an estimator. :func:`~torchsim.execution` decides where that work runs, and the timings below are taken inside it. .. GENERATED FROM PYTHON SOURCE LINES 30-173 .. code-block:: Python import time import numpy as np import torch import torchsim from torchsim.estimators import DictionaryMatcher, NonlinearLeastSquares from torchsim.simulators import MultiEchoSimulator .. GENERATED FROM PYTHON SOURCE LINES 174-183 Phantom ------- BrainWeb subject 0, slice 90: an axial slice at 1 mm through the lateral ventricles. BrainWeb publishes fuzzy memberships rather than labels, so each voxel holds a fraction of each tissue, and the relaxation times are weighted by those fractions. A third of the voxels are mixtures, so the truth is a continuum and not four values. .. GENERATED FROM PYTHON SOURCE LINES 184-217 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_03-nonlinear-least-squares_001.png :alt: BrainWeb subject 0, slice 90 :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_03-nonlinear-least-squares_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 0.9878790357573323, 'BrainWeb subject 0, slice 90') .. GENERATED FROM PYTHON SOURCE LINES 218-226 Measurement and noise floor --------------------------- Sixteen echoes out to 200 ms, scaled by the proton density and offset by a constant. A magnitude reconstruction rectifies the noise, so the late echoes measure the floor rather than zero, and a fit that does not model it absorbs it into T2. .. GENERATED FROM PYTHON SOURCE LINES 227-238 .. code-block:: Python ECHOES = 16 TE = torch.linspace(10.0, 200.0, ECHOES) NOISE_STD = 0.005 FLOOR = 0.05 simulator = MultiEchoSimulator(TE=TE) clean = simulator.simulate(T2=truth, M0=density, offset=FLOOR) generator = torch.Generator().manual_seed(42) measured = clean + NOISE_STD * torch.randn(clean.shape, generator=generator) .. GENERATED FROM PYTHON SOURCE LINES 239-242 The short-T2 voxel is on the floor by the fourth echo, so most of its train says nothing about T2; the long-T2 one never reaches it. .. GENERATED FROM PYTHON SOURCE LINES 243-294 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_03-nonlinear-least-squares_002.png :alt: what is measured :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_03-nonlinear-least-squares_002.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 295-306 Nonlinear fit ------------- What is unknown, over what range, and at what noise level. Every voxel steps together in the same pass, carries its own damping, accepts or rejects on its own, and drops out when it converges. The bounds are not clipping. Each is kept by fitting a transformed variable, so no iterate leaves the interval and no bound sits on the answer. That also puts every parameter on one scale, which is what the damping assumes. .. GENERATED FROM PYTHON SOURCE LINES 307-349 .. code-block:: Python BOUNDS = {"T2": (10.0, 500.0), "M0": (0.1, 2.0), "offset": (0.0, 0.2)} START = {"T2": 100.0, "M0": 1.0, "offset": 0.02} fit = NonlinearLeastSquares(simulator, bounds=BOUNDS, initial=START).fit( BOUNDS, noise_std=NOISE_STD, seed=0 ) maps = fit.map(measured) # {"T2": ..., "M0": ..., "offset": ...} .. rst-class:: sphx-glr-script-out .. code-block:: none fitted floor, median 0.0500 against 0.05 .. GENERATED FROM PYTHON SOURCE LINES 350-358 Dictionary match ---------------- The same problem given to a dictionary. Its grid needs no proton density: a match normalizes both sides, so any positive scale is free and one of the three parameters is gone before the grid is built. The offset survives normalization, so modelling it means putting it on the grid. .. GENERATED FROM PYTHON SOURCE LINES 359-363 .. code-block:: Python T2_GRID = torch.logspace(1.0, np.log10(500.0), 400) match = DictionaryMatcher(simulator.bind(M0=1.0, offset=0.0)).fit(T2=T2_GRID, seed=0) .. GENERATED FROM PYTHON SOURCE LINES 364-366 With the floor on the grid, the grid is the product of the two. .. GENERATED FROM PYTHON SOURCE LINES 367-402 .. code-block:: Python offsets = torch.linspace(0.0, 0.15, 40) grid_t2, grid_offset = torch.meshgrid(T2_GRID, offsets, indexing="ij") wide = DictionaryMatcher(simulator.bind(M0=1.0)).fit( T2=grid_t2.reshape(-1), offset=grid_offset.reshape(-1), seed=0 ) .. GENERATED FROM PYTHON SOURCE LINES 403-410 Cost and accuracy ----------------- Best of three passes each, after a warm-up. **model** is what the fitted estimator carries between volumes; **peak** is the high-water mark on the card while the slice was mapped, and a dash on a machine with no card. .. GENERATED FROM PYTHON SOURCE LINES 411-449 .. rst-class:: sphx-glr-script-out .. code-block:: none method atoms train map model peak T2 ---------------------------------------------------------------------------------- match, T2 only 400 0.0s 0.01s 0.1 MiB -- 17.9% match, T2 x 10 offsets 4000 0.0s 0.10s 0.5 MiB -- 1.8% match, T2 x 20 offsets 8000 0.0s 0.20s 1.1 MiB -- 1.8% match, T2 x 40 offsets 16000 0.0s 0.41s 2.1 MiB -- 1.8% match, T2 x 80 offsets 32000 0.0s 0.83s 4.3 MiB -- 1.7% fit, T2 + M0, floor known -- 0.0s 0.47s 0.0 MiB -- 0.6% fit, T2 + M0 + offset -- 0.0s 0.78s 0.0 MiB -- 1.8% .. GENERATED FROM PYTHON SOURCE LINES 450-468 Interpretation -------------- The first row is the trap. A T2-only match is the quickest thing here and the most wrong, because the model it matched was not the model that produced the data. Nothing in the estimator says so: the residual it minimized is small, at the wrong T2. With the offset on the grid the match recovers, and the cost of recovering is the point: ten values of one nuisance is ten times the atoms and ten times the memory, for a parameter that cost the fit one column. The fit is the slower of the two in wall clock and stays that way here: a Levenberg-Marquardt loop is tens of passes where a match is one. What it does not do is grow. Each nuisance multiplies the grid again and adds one column to the Jacobian, so where the two cross is arithmetic; the memory has crossed already. .. GENERATED FROM PYTHON SOURCE LINES 469-496 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_03-nonlinear-least-squares_003.png :alt: time, memory :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_03-nonlinear-least-squares_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none .. GENERATED FROM PYTHON SOURCE LINES 497-500 Maps ---- .. GENERATED FROM PYTHON SOURCE LINES 501-532 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_03-nonlinear-least-squares_004.png :alt: truth, match, T2, match + offset, fit, all three :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_03-nonlinear-least-squares_004.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 533-547 Limits ------ There is no guarantee. The fit started every voxel at 100 ms and landed on the right answer everywhere, which is a property of an exponential: the residual of a single decay has one minimum. A model whose residual has several -- a fingerprinting train, a fat-water fit at a wrong field map -- can be started in the wrong basin and stay there. A match cannot, because it scores every atom. Equality constraints belong in the model, not in the bounds. Two fractions that must sum to one are written with one as the unknown and the other as ``1 - f`` inside the model, so the constraint holds at every iterate. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 10.433 seconds) .. _sphx_glr_download_generated_autoexamples_02-parameter-inference_03-nonlinear-least-squares.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/02-parameter-inference/03-nonlinear-least-squares.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 03-nonlinear-least-squares.ipynb <03-nonlinear-least-squares.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 03-nonlinear-least-squares.py <03-nonlinear-least-squares.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 03-nonlinear-least-squares.zip <03-nonlinear-least-squares.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_