.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/02-parameter-inference/04-perk.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_04-perk.py: ============================= PERK: kernel ridge regression ============================= The scope of this notebook is to map a brain slice with PERK, to show what the size of the regression buys, and to read the error bar it reports. PERK never builds a dictionary. It is a kernel regression trained on signals drawn from a prior rather than laid on a grid, and at inference it projects a signal onto a fixed set of random Fourier features and reads the answer off a linear combination of them. Its cost per voxel does not depend on how many parameters are unknown, and the training is paid once. .. GENERATED FROM PYTHON SOURCE LINES 17-21 .. colab-link:: :needs_gpu: 0 !pip install torchsim brainweb-dl cmap .. GENERATED FROM PYTHON SOURCE LINES 23-27 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 28-174 .. code-block:: Python import time import numpy as np import torch import torchsim from torchsim import ( Subspace, ) from torchsim.estimators import PERK, DictionaryMatcher from torchsim.simulators import MRFSimulator .. GENERATED FROM PYTHON SOURCE LINES 175-184 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 185-225 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_001.png :alt: BrainWeb subject 0, slice 90 :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Text(0.5, 0.9883529057497298, 'BrainWeb subject 0, slice 90') .. GENERATED FROM PYTHON SOURCE LINES 226-233 Sequence -------- Four hundred repetitions after an inversion, at a fixed repetition time and a flip angle that varies smoothly along the train. A schedule that jumped about would give trajectories differing by noise rather than by physics. .. GENERATED FROM PYTHON SOURCE LINES 234-243 .. code-block:: Python CONTRASTS = 400 TR_MS = 10.0 TI_MS = 20.0 repetition = torch.arange(CONTRASTS, dtype=torch.float32) flip = 10.0 + 50.0 * torch.sin(torch.pi * repetition / CONTRASTS) ** 2 simulator = MRFSimulator(flip=flip, TR=TR_MS, TI=TI_MS, states=20, M0=1.0) .. GENERATED FROM PYTHON SOURCE LINES 244-249 The readouts wind the states on rather than rewinding them, so nothing returns transverse magnetization to the imaginary axis: the trajectory comes back real to within 3e-8, which halves both the dictionary and the arithmetic that searches it. .. GENERATED FROM PYTHON SOURCE LINES 250-266 .. code-block:: Python fingerprints = simulator.simulate( T1=torch.tensor([500.0, 833.0, 2569.0]), T2=torch.tensor([70.0, 83.0, 329.0]) ).real .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_002.png :alt: the schedule, fingerprints :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 267-272 The measurement, with noise at 2% of the peak fingerprint. One number sets it, and the same number is what the estimators are told to expect -- an estimator trained for more noise than the scan has learns to distrust the data and answers with the prior instead. .. GENERATED FROM PYTHON SOURCE LINES 273-316 .. code-block:: Python NOISE_STD = float(0.02 * fingerprints.max()) truth = { "T1": torch.as_tensor(T1_true[mask].copy()), "T2": torch.as_tensor(T2_true[mask].copy()), } density = torch.as_tensor(M0_true[mask].copy()) clean = simulator.simulate(**truth).real * density[:, None] generator = torch.Generator().manual_seed(42) measured = clean + NOISE_STD * torch.randn(clean.shape, generator=generator) .. GENERATED FROM PYTHON SOURCE LINES 317-325 Problem statement ----------------- What is unknown, over what range, and at what noise level. Both relaxation times span more than a decade, so the prior is drawn logarithmically: uniform sampling would spend most of the budget on long T1, where the trajectories are nearly parallel. .. GENERATED FROM PYTHON SOURCE LINES 326-331 .. code-block:: Python T1_RANGE = (200.0, 5000.0) T2_RANGE = (20.0, 600.0) SAMPLES = 20_000 prior = torch.Generator().manual_seed(11) .. GENERATED FROM PYTHON SOURCE LINES 332-339 Subspace basis -------------- Four hundred contrasts do not span four hundred directions. A basis fitted to simulated trajectories says how many they do span; one minus the energy it keeps is the relative squared error of projecting through it and back. .. GENERATED FROM PYTHON SOURCE LINES 340-374 .. code-block:: Python training_signals, _, _ = ( PERK(simulator) .fit( T1=log_uniform(*T1_RANGE, SAMPLES), T2=log_uniform(*T2_RANGE, SAMPLES), noise_std=NOISE_STD, seed=0, ) .training_set(SAMPLES) ) training_signals = training_signals.real RANK = 4 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_003.png :alt: what projecting through the basis loses :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 375-384 Four directions leave less outside the basis than the noise puts in, so rank 4 is used from here on. Training -------- Twenty thousand parameter pairs drawn from the prior, simulated, and given the noise the scan has. Fitting is a linear solve against the features. .. GENERATED FROM PYTHON SOURCE LINES 385-450 .. code-block:: Python FEATURES = 1000 perk = PERK(simulator, n_features=FEATURES, regularization=1e-6, normalize=True).fit( T1=log_uniform(*T1_RANGE, SAMPLES), T2=log_uniform(*T2_RANGE, SAMPLES), noise_std=NOISE_STD, seed=0, rank=RANK, samples=SAMPLES, ) maps = perk(measured) # {"T1": ..., "T2": ...}, one value per voxel .. rst-class:: sphx-glr-script-out .. code-block:: none dictionary: 20000 atoms at rank 4; regression: 1000 features from 20000 training draws .. GENERATED FROM PYTHON SOURCE LINES 451-459 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. The dictionary row is the reference point. .. GENERATED FROM PYTHON SOURCE LINES 460-499 .. rst-class:: sphx-glr-script-out .. code-block:: none method train map model peak T1 T2 ------------------------------------------------------------------------------ match, rank 4 1.0s 1.54s 1.4 MiB -- 0.6% 1.6% PERK, 500 features 1.2s 0.03s 0.0 MiB -- 0.9% 13.2% PERK, 1000 features 1.6s 0.04s 0.1 MiB -- 1.3% 5.7% PERK, 4000 features 8.6s 0.10s 0.4 MiB -- 0.9% 5.5% .. GENERATED FROM PYTHON SOURCE LINES 500-503 Maps ---- .. GENERATED FROM PYTHON SOURCE LINES 504-562 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_004.png :alt: truth, match, PERK :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_004.png :class: sphx-glr-multi-img * .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_005.png :alt: Δ match, Δ PERK :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_005.png :class: sphx-glr-multi-img .. GENERATED FROM PYTHON SOURCE LINES 563-571 Uncertainty ----------- ``uncertainty=True`` returns a second set of maps: how far the answer is expected to sit from the truth. PERK learns that at training, from the residuals of its own fit, so reporting it is a matrix multiply rather than a rerun of the volume. .. GENERATED FROM PYTHON SOURCE LINES 572-574 .. code-block:: Python maps, spread = perk(measured, uncertainty=True) .. GENERATED FROM PYTHON SOURCE LINES 575-579 Read against the Cramer-Rao bound, the lowest standard deviation an unbiased estimate could reach from this train at this noise. The bound belongs to the sequence, so the gap is what the method loses. .. GENERATED FROM PYTHON SOURCE LINES 580-623 .. code-block:: Python _signal, sensitivity = simulator.jacobian("T1 T2".split(), **truth) sensitivity = sensitivity.real * density[:, None, None] floor = torchsim.crlb(sensitivity, noise_variance=NOISE_STD**2, singular="infinite") bound = {"T1": floor[:, 0].sqrt(), "T2": floor[:, 1].sqrt()} .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_006.png :alt: PERK, CRLB, PERK, relative :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_04-perk_006.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none PERK CRLB PERK (median over the brain) T1 8.4 ms 2.4 ms 1.2% T2 10.2 ms 0.5 ms 12.8% .. GENERATED FROM PYTHON SOURCE LINES 624-637 Each row uses its own parameter's colormap. The two absolute panels share a scale; the third is the same spread as a percentage of the relaxation time, which is what says whether ten milliseconds is tight. Both are largest in CSF, whose long T1 this train resolves least. The gap between them is not: T1 sits within a small multiple of the bound, T2 several times above it, and it is T2 whose error moved when the feature count was swept. One is a sequence to redesign, the other a regression to enlarge. The number is not the noise alone. A regression trained on a prior answers with the prior where the data is weak, and is wrong the same way in every realization, so repeating the scan would never show that part. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 27.461 seconds) .. _sphx_glr_download_generated_autoexamples_02-parameter-inference_04-perk.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/04-perk.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04-perk.ipynb <04-perk.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04-perk.py <04-perk.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04-perk.zip <04-perk.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_