.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/03-regularization/02-operators-and-solvers.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. .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_examples_03-regularization_02-operators-and-solvers.py: ===================== Operators and solvers ===================== This lesson rebuilds the reconstruction of the previous lessons from its parts -- the encoding operator, the regularization term and the iterative algorithm -- instead of calling a BART application, and shows that the result is identical. :func:`bartorch.apps.pics` builds three objects and runs BART's iteration with them: the SENSE encoding operator :math:`A = PFS`, the regularization terms, and the algorithm. Building them separately with :mod:`bartorch.linop` and :mod:`bartorch.optim` is what a reconstruction BART has no application for requires: an encoding with an additional factor, such as an off-resonance or motion-induced phase, a solver called from an outer loop, an operator defined in Python, or a gradient with respect to the data for training a network. The example builds the Cartesian SENSE encoding of :doc:`../01-basics/02-from-kspace-to-image`, checks it against the definition of the adjoint, solves with it, and compares the result with the application. The phantom, the coil sensitivities and the sampling are that example's; the cell that builds them is hidden on this page and present in the script this page can be downloaded as. **Learning objectives** - Build :class:`bartorch.linop.CartesianSense` and read the plan it was lowered into. - Check an operator's adjoint with the dot-product test. - Prepare data as ``pics`` does and reproduce :func:`bartorch.apps.pics` bit for bit with :class:`bartorch.optim.FISTA`. - Compose operators, include one defined in Python, and differentiate through an application. It follows :doc:`01-regularized-reconstruction`. The next section, :doc:`../04-non-cartesian/01-trajectories-and-transforms`, uses these operators off the Cartesian grid. .. GENERATED FROM PYTHON SOURCE LINES 42-253 .. code-block:: Python import csv from pathlib import Path import brainweb_dl import numpy as np import torch from brainweb_dl import get_mri import bartorch import bartorch.tools as bt from bartorch import apps, linop, optim, priors SIZE = 192 COILS = 8 ACCELERATION = 3 CALIBRATION = 24 .. GENERATED FROM PYTHON SOURCE LINES 254-268 The encoding operator --------------------- :func:`bartorch.linop.CartesianSense` is :math:`A = P F S` as one operator. It takes the sensitivities, the shape of the image it maps from, and the sampling pattern; the shape of the k-space it maps to follows from those. :func:`bartorch.tools.pattern` reads the pattern off the measured data, as for a prospectively undersampled acquisition. ``modulated=True`` selects BART's uncentred sample convention, which its applications iterate in; the default is the centred convention that :func:`bartorch.fft` produces. The two differ by a modulation of the samples and give the same image, so the choice matters only when the operator is applied to data already in one of them, as it is below. .. GENERATED FROM PYTHON SOURCE LINES 269-277 .. code-block:: Python pattern = bt.pattern(kspace) A = linop.CartesianSense(maps.squeeze(1), (SIZE, SIZE), pattern.squeeze(), modulated=True) print(f"{A.ishape} -> {A.oshape}") print(A.plan) print(f"fused: {A.plan.fused}") .. rst-class:: sphx-glr-script-out .. code-block:: none (192, 192) -> (8, 192, 192) Plan(transform=fft, image=sensitivities, normal=transform, coil_batch=1, streamed=coils, executor=slab) fused: True .. GENERATED FROM PYTHON SOURCE LINES 278-290 ``A.plan`` reports the form the operator was lowered into: which transform, what multiplies the image and the samples, and how the normal operator :math:`A^H A` is applied. It is a property of the built operator rather than a prediction, and ``plan.fused`` is false where the composition could not be expressed as one encoding and fell back to a chain, which computes the same numbers more slowly. Applying the adjoint is not the same as applying the transpose, and a reconstruction built on the wrong one converges to the wrong image. The definition :math:`\langle Ax, y\rangle = \langle x, A^H y\rangle` holds for any pair of vectors, and holds for random vectors as readily as for real data, so it is a usable check on an operator. .. GENERATED FROM PYTHON SOURCE LINES 291-300 .. code-block:: Python generator = torch.Generator().manual_seed(0) probe = torch.randn(A.ishape, dtype=torch.complex64, generator=generator) samples = torch.randn(A.oshape, dtype=torch.complex64, generator=generator) forward = (A(probe).conj() * samples).sum() adjoint = (probe.conj() * A.H(samples)).sum() print(f"relative difference {abs(forward - adjoint) / abs(forward):.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none relative difference 1.35e-07 .. GENERATED FROM PYTHON SOURCE LINES 301-310 Solving ------- A solver is called as ``solver(y, A)``. What it is given is not the array the scanner wrote but what ``pics`` iterates on: the sampling pattern applied, the modulation into the uncentred convention, and the data divided by the scaling :func:`bartorch.optim.data_scaling` estimates from the adjoint reconstruction, which is the step that makes a regularization weight transferable from one dataset to the next. .. GENERATED FROM PYTHON SOURCE LINES 311-319 .. code-block:: Python measured = bartorch.fftmod(kspace * pattern, axes=(-1, -2, -3), inverse=True) scale = optim.data_scaling(measured) data = (measured / scale).squeeze(1) term = priors.Wavelet(axes=(-1, -2), weight=0.002) assembled = optim.FISTA(term, maxiter=100)(data, A) .. GENERATED FROM PYTHON SOURCE LINES 320-323 With the same preprocessing the assembled solve and the application are not merely close: they are the same iteration over the same operator, and return the same bits. .. GENERATED FROM PYTHON SOURCE LINES 324-330 .. code-block:: Python tool = apps.pics(kspace, maps, regularizers=term, solver="fista", maxiter=100) print(f"identical to pics: {torch.equal(assembled.squeeze(), tool.squeeze())}") print(f"NRMSE, adjoint {bt.nrmse(image.abs(), A.H(data).abs(), scaled=True):.3f}") print(f"NRMSE, FISTA {bt.nrmse(image.abs(), assembled.abs(), scaled=True):.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none identical to pics: True NRMSE, adjoint 0.119 NRMSE, FISTA 0.032 .. GENERATED FROM PYTHON SOURCE LINES 331-351 .. image-sg:: /auto_examples/03-regularization/images/sphx_glr_02-operators-and-solvers_001.png :alt: reference, adjoint, $A^H y$, FISTA, wavelet, |error|, FISTA :srcset: /auto_examples/03-regularization/images/sphx_glr_02-operators-and-solvers_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 352-369 The adjoint of the encoding is not its inverse: :math:`A^H y` is the sensitivity-weighted coil combination of the zero-filled k-space, and carries the aliasing of the undersampling and the shading of :math:`\sum_c |S_c|^2`, which the solve removes, as the NRMSE printed above shows. The error of the solution, at a tenth of the image peak, is concentrated at the tissue boundaries. Operator algebra ---------------- ``@`` composes, ``+`` adds, ``A.H`` is the adjoint and ``A.gram()`` the normal operator :math:`A^H A`. A composition builds a single BART operator rather than a Python chain, so a solver iterating on it does not return to Python between applications. :func:`bartorch.optim.maxeigen` runs the power iteration on an operator, which is how a gradient step size is chosen: the Lipschitz constant of the least-squares gradient is the largest eigenvalue of :math:`A^H A`. .. GENERATED FROM PYTHON SOURCE LINES 370-373 .. code-block:: Python print(f"largest eigenvalue of A^H A: {optim.maxeigen(A.gram()):.3f}") .. rst-class:: sphx-glr-script-out .. code-block:: none largest eigenvalue of A^H A: 0.999 .. GENERATED FROM PYTHON SOURCE LINES 374-378 An operator defined in Python is composed with BART's through :meth:`~bartorch.linop.LinearOperator.from_callbacks`, which BART applies as a callback. Here it is a spatially varying phase, as an off-resonance or an eddy-current phase would be, placed between the image and the encoding. .. GENERATED FROM PYTHON SOURCE LINES 379-387 .. code-block:: Python field = torch.exp(1j * 0.4 * torch.pi * grid_x).to(torch.complex64) phase = linop.LinearOperator.from_callbacks( (SIZE, SIZE), (SIZE, SIZE), lambda u: field * u, lambda u: field.conj() * u ) composed = A @ phase print(f"{composed.ishape} -> {composed.oshape}, fused: {composed.plan.fused}") .. rst-class:: sphx-glr-script-out .. code-block:: none (192, 192) -> (8, 192, 192), fused: True .. GENERATED FROM PYTHON SOURCE LINES 388-398 Differentiation --------------- Applying an operator to a tensor that requires a gradient records the application for autograd. The gradient torch propagates back through :math:`y = Ax` is :math:`A^H g` rather than :math:`A^T g`, the conjugate Wirtinger convention torch uses for complex tensors. For a real :math:`A`, :math:`A^H = A^T`, so only a complex check distinguishes the two; :doc:`../../explanation/differentiation` describes the backward passes of the solvers. .. GENERATED FROM PYTHON SOURCE LINES 399-408 .. code-block:: Python variable = data.new_zeros(A.ishape).requires_grad_(True) residual = A(variable) - data (residual.abs() ** 2).sum().backward() expected = 2 * A.H(-data) difference = float((variable.grad - expected).abs().max() / expected.abs().max()) print(f"relative difference from 2 A^H (Ax - y): {difference:.2e}") .. rst-class:: sphx-glr-script-out .. code-block:: none relative difference from 2 A^H (Ax - y): 2.04e-07 .. GENERATED FROM PYTHON SOURCE LINES 409-413 The regularization terms are the subject of :mod:`bartorch.priors`, and the iterations of :mod:`bartorch.optim`; :doc:`../../explanation/inverse-problems` states which algorithm applies to which problem. .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 2.284 seconds) .. _sphx_glr_download_auto_examples_03-regularization_02-operators-and-solvers.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 02-operators-and-solvers.ipynb <02-operators-and-solvers.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 02-operators-and-solvers.py <02-operators-and-solvers.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 02-operators-and-solvers.zip <02-operators-and-solvers.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_