.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "generated/autoexamples/02-parameter-inference/01-dictionary-matching.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_01-dictionary-matching.py: =================== Dictionary matching =================== The scope of this notebook is to map a brain slice by exhaustive dictionary matching, and to show the two ways of making that affordable: working in the low-rank basis the train spans, and clustering the dictionary so that most atoms are never scored. A dictionary spans every combination of the parameters, so its size is the product of the grids and each atom is as long as the train. The two savings are independent and multiply; what each costs and what each gets wrong is read off the same slice. .. GENERATED FROM PYTHON SOURCE LINES 18-22 .. colab-link:: :needs_gpu: 0 !pip install torchsim brainweb-dl cmap .. GENERATED FROM PYTHON SOURCE LINES 24-28 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 29-175 .. code-block:: Python import time import numpy as np import torch import torchsim from torchsim import ( Subspace, ) from torchsim.estimators import DictionaryMatcher from torchsim.simulators import MRFSimulator .. GENERATED FROM PYTHON SOURCE LINES 176-185 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 186-226 .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_01-dictionary-matching_001.png :alt: BrainWeb subject 0, slice 90 :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_01-dictionary-matching_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none Downloading tissues: 0%| | 0/10 [00:00] .. GENERATED FROM PYTHON SOURCE LINES 267-271 The measurement, with noise at 2% of the peak fingerprint. The estimators are told the same number: one trained for more noise than the scan has learns to distrust the data and answers with the prior. .. GENERATED FROM PYTHON SOURCE LINES 272-315 .. 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 316-324 Problem statement ----------------- What is unknown, over what range, and at what noise level. Both relaxation times span more than a decade, so the grid is logarithmic: uniform spacing would spend most of it on long T1, where the trajectories are nearly parallel. .. GENERATED FROM PYTHON SOURCE LINES 325-330 .. 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 331-338 Subspace rank ------------- A basis fitted to simulated trajectories says how much of their energy each rank keeps. One minus that fraction is the relative squared error of projecting through the basis and back. .. GENERATED FROM PYTHON SOURCE LINES 339-373 .. code-block:: Python training_signals, _, _ = ( DictionaryMatcher(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_01-dictionary-matching_003.png :alt: what projecting through the basis loses :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_01-dictionary-matching_003.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-script-out .. code-block:: none [] .. GENERATED FROM PYTHON SOURCE LINES 374-385 Four directions out of four hundred leave less outside the basis than the noise puts in, and every contrast dropped is arithmetic the match avoids. Dictionary ---------- The dictionary spans the parameters jointly, so its size is the product of the grids: twenty thousand atoms for two parameters, on a grid fine enough that the spacing is not what limits the answer. A third parameter multiplies it again. .. GENERATED FROM PYTHON SOURCE LINES 386-406 .. code-block:: Python T1_GRID = torch.logspace(np.log10(T1_RANGE[0]), np.log10(T1_RANGE[1]), 200) T2_GRID = torch.logspace(np.log10(T2_RANGE[0]), np.log10(T2_RANGE[1]), 100) grid_t1, grid_t2 = torch.meshgrid(T1_GRID, T2_GRID, indexing="ij") full = DictionaryMatcher(simulator).fit( T1=grid_t1.reshape(-1), T2=grid_t2.reshape(-1), seed=0 ) maps = full.map(measured) # {"T1": ..., "T2": ...}, one value per voxel .. GENERATED FROM PYTHON SOURCE LINES 407-414 Matching in the subspace ------------------------ ``rank`` is the whole change. The dictionary is fitted, projected and stored in four directions instead of four hundred, and the measurement is projected the same way before scoring. .. GENERATED FROM PYTHON SOURCE LINES 415-429 .. code-block:: Python low = DictionaryMatcher(simulator).fit( T1=grid_t1.reshape(-1), T2=grid_t2.reshape(-1), seed=0, rank=RANK ) .. GENERATED FROM PYTHON SOURCE LINES 430-441 Clustered dictionary -------------------- Compressing shortened every inner product; grouping cuts how many are taken. Neighbouring tissues give nearly parallel signals, so the atoms cluster, and a voxel scored against one representative per group rules out most groups before any atom inside them is touched. The clustering is done in the compressed basis, so a group is entered without leaving the space the measurement is already in. .. GENERATED FROM PYTHON SOURCE LINES 442-479 .. code-block:: Python GROUPS = 32 grouped = DictionaryMatcher(simulator, groups=GROUPS).fit( T1=grid_t1.reshape(-1), T2=grid_t2.reshape(-1), seed=0, rank=RANK, ) .. rst-class:: sphx-glr-script-out .. code-block:: none 32 groups of 625 atoms; 1.4 still open per voxel .. GENERATED FROM PYTHON SOURCE LINES 480-488 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, which decides whether a volume fits or has to be streamed, and a dash on a machine with no card. .. GENERATED FROM PYTHON SOURCE LINES 489-544 .. rst-class:: sphx-glr-script-out .. code-block:: none method train map model peak T1 T2 M0 ---------------------------------------------------------------------------------------- match, 400 contrasts 8.0s 7.58s 122.3 MiB -- 0.6% 1.6% 0.4% match, rank 4 1.0s 1.54s 1.4 MiB -- 0.6% 1.6% 0.4% match, rank 4 + groups 1.0s 0.09s 1.4 MiB -- 0.6% 1.6% 0.4% .. GENERATED FROM PYTHON SOURCE LINES 545-548 Maps ---- .. GENERATED FROM PYTHON SOURCE LINES 549-606 .. rst-class:: sphx-glr-horizontal * .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_01-dictionary-matching_004.png :alt: truth, full, rank 4, + groups :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_01-dictionary-matching_004.png :class: sphx-glr-multi-img * .. image-sg:: /generated/autoexamples/02-parameter-inference/images/sphx_glr_01-dictionary-matching_005.png :alt: Δ full, Δ rank 4, Δ + groups :srcset: /generated/autoexamples/02-parameter-inference/images/sphx_glr_01-dictionary-matching_005.png :class: sphx-glr-multi-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 53.892 seconds) .. _sphx_glr_download_generated_autoexamples_02-parameter-inference_01-dictionary-matching.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/01-dictionary-matching.ipynb :alt: Launch binder :width: 150 px .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 01-dictionary-matching.ipynb <01-dictionary-matching.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 01-dictionary-matching.py <01-dictionary-matching.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 01-dictionary-matching.zip <01-dictionary-matching.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_