.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_examples/03_algorithms/04_run_depierro.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_algorithms_04_run_depierro.py: DePierro's algorithm to optimize the Poisson logL with quadratic intensity prior ================================================================================ This example demonstrates the use of DePierro's algorithm to minimize the negative Poisson log-likelihood function combined with a quadratic intensity prior: .. math:: f(x) = \sum_{i=1}^m \bar{y}_i - y_i \log(\bar{y}_i (x)) + \frac{\beta}{2} \|x - z \|^2 subject to .. math:: x \geq 0 using the linear forward model .. math:: \bar{y}(x) = A x + s .. GENERATED FROM PYTHON SOURCE LINES 23-29 .. code-block:: Python import matplotlib.pyplot as plt import numpy as np import parallelproj.operators from parallelproj import to_numpy_array .. GENERATED FROM PYTHON SOURCE LINES 30-37 .. code-block:: Python from array_utils import suggest_array_backend_and_device # To use a specific backend and/or device, replace the None arguments, e.g.: # xp, dev = suggest_array_backend_and_device(backend="numpy", dev="cpu") or by setting xp and dev manually xp, dev = suggest_array_backend_and_device(None, None) .. rst-class:: sphx-glr-script-out .. code-block:: none Using array API: array_api_compat.torch, device: cpu .. GENERATED FROM PYTHON SOURCE LINES 38-47 Setup of the forward model :math:`\bar{y}(x) = A x + s` -------------------------------------------------------- We setup a minimal linear forward operator :math:`A` respresented by a 4x4 matrix and an arbritrary contamination vector :math:`s` of length 4. .. note:: The OSEM implementation below works with all linear operators that subclass :class:`.LinearOperator` (e.g. the high-level projectors). .. GENERATED FROM PYTHON SOURCE LINES 47-65 .. code-block:: Python # setup an arbitrary 4x4 matrix mat = xp.asarray( [ [2.5, 1.2, 0.3, 0.1], [0.4, 3.1, 0.7, 0.2], [0.1, 0.3, 4.1, 2.5], [0.2, 0.5, 0.2, 0.9], [0.3, 0.1, 0.7, 0.2], ], dtype=xp.float64, device=dev, ) op_A = parallelproj.operators.MatrixOperator(mat) # setup an arbitrary contamination vector that has shape op_A.out_shape contamination = xp.asarray([0.3, 0.2, 0.1, 0.4, 0.1], dtype=xp.float64, device=dev) .. GENERATED FROM PYTHON SOURCE LINES 66-71 Setup of ground truth and data simulation ----------------------------------------- We setup an arbitrary ground truth :math:`x_{true}` and simulate noise-free and noisy data :math:`y` by adding Poisson noise. .. GENERATED FROM PYTHON SOURCE LINES 71-87 .. code-block:: Python # ground truth x_true = xp.asarray([5.5, 10.7, 8.2, 7.9], dtype=xp.float64, device=dev) # simulated noise-free data noise_free_data = op_A(x_true) + contamination # add Poisson noise np.random.seed(1) y = xp.asarray( np.random.poisson(to_numpy_array(noise_free_data)), device=dev, dtype=xp.float64, ) .. GENERATED FROM PYTHON SOURCE LINES 88-96 .. code-block:: Python x_prior = xp.full(op_A.in_shape, xp.min(x_true), dtype=xp.float64, device=dev) beta = 0.3 num_iter = 50 # initialize x x = xp.ones(op_A.in_shape, dtype=xp.float64, device=dev) .. GENERATED FROM PYTHON SOURCE LINES 97-106 .. code-block:: Python def cost_function(x): exp = op_A(x) + contamination if (xp.min(exp) < 0) or (xp.min(x) < 0): res = xp.finfo(xp.float64).max else: res = (xp.sum(exp - y * xp.log(exp))) + 0.5 * beta * xp.sum((x - x_prior) ** 2) return res .. GENERATED FROM PYTHON SOURCE LINES 107-122 DePierro update to minimize :math:`f(x)` ---------------------------------------- We apply multiple DePierro updates .. math:: b = A^H 1 - \beta z .. math:: t = x A^H \frac{y}{A x + s} .. math:: x^+ = \frac{2 t}{\sqrt{b^2 + 4 \beta t} + b} to calculate the minimizer of :math:`f(x)` iteratively. See :cite:p:`DePierro1995` for more details. .. GENERATED FROM PYTHON SOURCE LINES 122-142 .. code-block:: Python cost = xp.zeros(num_iter, dtype=xp.float64, device=dev) # "b" - modified sensitivity image mod_adjoint_ones = ( op_A.adjoint(xp.ones(op_A.out_shape, dtype=xp.float64, device=dev)) - beta * x_prior ) for i in range(num_iter): # evaluate the forward model exp = op_A(x) + contamination ratio = y / exp t = x * op_A.adjoint(ratio) x = 2 * t / (xp.sqrt(mod_adjoint_ones**2 + 4 * beta * t) + mod_adjoint_ones) cost[i] = cost_function(x) print(f"Solution after {num_iter} DePierro iterations:") print(x) print(f"cost: {cost[-1]:.6e}") .. rst-class:: sphx-glr-script-out .. code-block:: none Solution after 50 DePierro iterations: tensor([6.2868, 7.2615, 6.7209, 6.4391], dtype=torch.float64) cost: -3.321656e+02 .. GENERATED FROM PYTHON SOURCE LINES 143-154 .. code-block:: Python if xp.__name__.endswith("numpy"): from scipy.optimize import fmin_powell x_ref = fmin_powell(cost_function, x, xtol=1e-6, ftol=1e-6) rel_dist = float(xp.sum((x - x_ref) ** 2)) / float(xp.sum(x_ref**2)) print(f"\nReference solution using Powell optimizer:") print(x_ref) print(f"rel. distance to DePierro solution: {rel_dist:.2e}") print(f"cost: {cost_function(x_ref):.6e}") .. GENERATED FROM PYTHON SOURCE LINES 155-163 .. code-block:: Python fig, ax = plt.subplots(1, 1, tight_layout=True) if xp.__name__.endswith("numpy"): ax.axhline(cost_function(x_ref), color="k", linestyle="--") ax.plot(to_numpy_array(cost)) ax.set_xlabel("iteration") ax.set_ylabel("cost function") ax.grid(ls=":") fig.show() .. image-sg:: /auto_examples/03_algorithms/images/sphx_glr_04_run_depierro_001.png :alt: 04 run depierro :srcset: /auto_examples/03_algorithms/images/sphx_glr_04_run_depierro_001.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 0.106 seconds) .. _sphx_glr_download_auto_examples_03_algorithms_04_run_depierro.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: 04_run_depierro.ipynb <04_run_depierro.ipynb>` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: 04_run_depierro.py <04_run_depierro.py>` .. container:: sphx-glr-download sphx-glr-download-zip :download:`Download zipped: 04_run_depierro.zip <04_run_depierro.zip>` .. only:: html .. rst-class:: sphx-glr-signature `Gallery generated by Sphinx-Gallery `_