From 78733985367182b1f1a24912b94919dbdb1ac2b3 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 30 Aug 2026 14:31:49 +0100 Subject: [PATCH 1/3] feat: added Downsample2D operator --- .claude/skills/newop/SKILL.md | 126 ++++++++++ .../newop/reference/operator_template.py | 97 ++++++++ .../skills/newop/reference/test_template.py | 57 +++++ AIPOLICY.md | 1 + docs/source/api/index.rst | 1 + examples/plot_downsample.py | 140 +++++++++++ pylops/signalprocessing/__init__.py | 3 + pylops/signalprocessing/downsample2d.py | 235 ++++++++++++++++++ pylops/signalprocessing/radon2d.py | 2 +- pylops/signalprocessing/radon3d.py | 2 +- pytests/test_downsample.py | 156 ++++++++++++ pytests/test_radon.py | 39 +++ tutorials/ctscan.py | 167 ------------- 13 files changed, 857 insertions(+), 169 deletions(-) create mode 100644 .claude/skills/newop/SKILL.md create mode 100644 .claude/skills/newop/reference/operator_template.py create mode 100644 .claude/skills/newop/reference/test_template.py create mode 100644 examples/plot_downsample.py create mode 100644 pylops/signalprocessing/downsample2d.py create mode 100644 pytests/test_downsample.py delete mode 100755 tutorials/ctscan.py diff --git a/.claude/skills/newop/SKILL.md b/.claude/skills/newop/SKILL.md new file mode 100644 index 00000000..08f2517a --- /dev/null +++ b/.claude/skills/newop/SKILL.md @@ -0,0 +1,126 @@ +--- +name: newop +description: Create a new PyLops linear operator following docs/source/adding.rst - class file, docstring, tests, docs entry and example. Use when the user asks to add/implement/port a new operator into PyLops, including porting an existing non-PyLops forward/adjoint implementation from a URL or a local file (e.g. "add a Foo operator", "turn this script into a PyLops operator", "port the operator at "). +--- + +Goal: add a new, PyLops-compliant `LinearOperator` to the library, complete with +docstring, registration, tests, docs entry and a gallery example, following +`docs/source/adding.rst` (the authoritative guide - read it if unsure). + +The operator may be written from scratch (from a mathematical description) or +**ported** from an existing non-PyLops implementation supplied as a **web link** +or a **local file**. Ask for the operator name and the source only if neither is +inferable from the invocation. + +## 0. Get the source material + +- **Web link**: fetch it with `WebFetch` (or the browser tools if the page needs + JS). Extract the actual forward/adjoint code, not the prose. +- **Local file**: read it in full. +- **Neither**: work from the user's mathematical description, and state the + assumed definition of the operator before writing code. + +Then write down explicitly, before touching `pylops/`: + +- what the forward map does, and its input/output shapes; +- whether the source's "adjoint" is a true adjoint (\(\mathbf{A}^H\)) or merely an + inverse/transpose/approximation - **this is the most common porting bug**; +- which source parameters become `__init__` arguments, which become derived + members, and which are irrelevant (e.g. plotting, I/O, CLI args); +- whether the operator is real- or complex-linear, and whether it is `explicit`. + +If the source adjoint is not the true adjoint, say so and implement the correct +adjoint - the dot-test in step 4 will fail otherwise. Never relax the dot-test +tolerance to make a wrong adjoint pass. + +## 1. Place the file + +- One class per file; file named after the class but **lowercase** + (`pylops/basicoperators/diagonal.py` holds `Diagonal`). Choose the subpackage by + theme: `basicoperators`, `signalprocessing`, `waveeqprocessing`, `optimization`, + etc. Create a new subpackage only if nothing fits. +- If the operator is just a composition of existing operators, write a **function** + returning the composed operator instead of a class (see `pylops.Laplacian`). +- Start the file with `__all__ = [""]`. +- Register it: add the import/`__all__` entry in the subpackage `__init__.py` + (and its module-level summary table), plus the top-level `pylops/__init__.py` + if the operator is meant to be user-facing as `pylops.`. + +## 2. Write the class + +Use `reference/operator_template.py` as the skeleton. Key rules: + +- Inherit from `pylops.LinearOperator` and initialize via + `super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name)`. + Prefer `dims`/`dimsd` over setting `shape` directly; `shape` is derived. + Set `explicit=True` only when the operator also exposes a dense matrix `A`. +- Decorate `_matvec`/`_rmatvec` with `@reshaped` when the operator is + n-dimensional, so `x` arrives shaped as `dims` (`dimsd` for `_rmatvec`) and the + return value is flattened for you. +- Use the backend helpers rather than raw NumPy so CuPy/JAX work: + `pylops.utils.backend.get_array_module`, `to_cupy_conditional`, and friends. + Do not `import numpy` for array creation inside `_matvec`/`_rmatvec`. +- Type-annotate with `pylops.utils.typing` (`NDArray`, `DTypeLike`, + `InputDimsLike`). +- Keep a `name` argument (default a short string) for `pylops.utils.describe`. +- Write the `numpydoc` docstring with, at minimum: one-line summary, expanded + description, `Parameters`, `Attributes` (when non-obvious), `Raises` (when the + `__init__` validates inputs), and a `Notes` section giving the maths of forward + and adjoint in `.. math::` blocks. Match the level of detail of neighbouring + operators. +- Add `.. versionadded:: ` to the class docstring for a brand-new + operator (check the current version in `pylops/version.py` / `pyproject.toml`). + +## 3. Add tests + +Add to the existing `pytests/test_*.py` matching the subpackage, or create a new +one following the same header (the `TEST_CUPY_PYLOPS` / `backend` guard block). +Follow `reference/test_template.py`: + +- module-level `par*` dicts, parametrized with `@pytest.mark.parametrize("par", [...])` + covering real/complex and, where relevant, square/over-/under-determined; +- an `assert dottest(Op, nr, nc, rtol=..., complexflag=0 if par["imag"] == 0 else 3, backend=backend)` + in every test of a new configuration; +- a forward check against an independently computed expected result + (e.g. `Op.todense() @ x`, or the original source implementation's output); +- an inversion round-trip with `lsqr` / `Op / y` and `assert_array_almost_equal` + when the operator is invertible; +- error-path tests for anything the `__init__` raises. + +## 4. Run + +Always use `uv`: + +```bash +uv run pytest pytests/test_.py -k -q +make lint_uv +``` + +Iterate until the dot-test and all assertions pass cleanly. + +## 5. Document + +- Add the operator name to the right `autosummary` block in + `docs/source/api/index.rst`. +- Add a gallery example `examples/plot_.py` (or a tutorial in + `tutorials/` for a heavier workflow), following the sphinx-gallery format of + `examples/plot_diagonal.py`: `r"""` title/underline/description `"""` header, + then `###...` comment blocks separating narrative from code, and matplotlib + figures showing forward and adjoint (and inversion, if relevant). + +## 6. Final checklist (from `docs/source/adding.rst`) + +Report back confirming each item: + +- [ ] single class (or function) in its own file, in a suitable `pylops` subpackage +- [ ] `__init__`, `_matvec`, `_rmatvec` implemented (plus `todense`/`matrix` if cheap) +- [ ] operator exported from the subpackage and top-level `__init__.py` +- [ ] numpydoc docstring with `Parameters` and a mathematical `Notes` section +- [ ] test added, `dottest` passes, forward/inverse checked +- [ ] listed in `docs/source/api/index.rst` +- [ ] used in at least one `examples/` script or `tutorials/` script +- [ ] `make lint_uv` clean + +When porting, close with a short note on what differed between the source +implementation and the PyLops version (adjoint correction, shape/flattening +conventions, dtype handling, removed I/O). diff --git a/.claude/skills/newop/reference/operator_template.py b/.claude/skills/newop/reference/operator_template.py new file mode 100644 index 00000000..316e901a --- /dev/null +++ b/.claude/skills/newop/reference/operator_template.py @@ -0,0 +1,97 @@ +"""Skeleton for a new PyLops operator. + +Copy into ``pylops//.py``, rename, and fill in. +Delete anything that does not apply. +""" + +__all__ = ["MyOperator"] + +import numpy as np + +from pylops import LinearOperator +from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import get_array_module, to_cupy_conditional +from pylops.utils.decorators import reshaped +from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray + + +class MyOperator(LinearOperator): + r"""One-line summary of the operator. + + Longer description of what the operator applies in forward mode and what + its adjoint does. + + .. versionadded:: X.Y.Z + + Parameters + ---------- + param : :obj:`numpy.ndarray` + Description of the main parameter. + dims : :obj:`list` or :obj:`int`, optional + Number of samples for each dimension of the model. + axis : :obj:`int`, optional + Axis along which the operator is applied. + dtype : :obj:`str`, optional + Type of elements in input array. + name : :obj:`str`, optional + Name of operator (to be used by :func:`pylops.utils.describe.describe`) + + Attributes + ---------- + shape : :obj:`tuple` + Operator shape. + explicit : :obj:`bool` + Operator contains a matrix that can be solved explicitly (``True``) or + not (``False``). + + Raises + ------ + ValueError + If ``param`` has incompatible size with ``dims``. + + Notes + ----- + In forward mode the operator applies + + .. math:: + y_i = \ldots \quad \forall i=1,2,\ldots,N + + and in adjoint mode + + .. math:: + x_i = \ldots \quad \forall i=1,2,\ldots,M + + """ + + def __init__( + self, + param: NDArray, + dims: int | InputDimsLike | None = None, + axis: int = -1, + dtype: DTypeLike = "float64", + name: str = "M", + ) -> None: + self.param = param + self.axis = axis + dims = param.shape if dims is None else _value_or_sized_to_tuple(dims) + # dimsd is the shape of the data (output of the forward) + dimsd = dims + super().__init__(dtype=np.dtype(dtype), dims=dims, dimsd=dimsd, name=name) + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + if type(self.param) is not type(x): + self.param = to_cupy_conditional(x, self.param) + y = ncp.zeros(self.dimsd, dtype=self.dtype) + # ... forward implementation, y = A x + return y + + @reshaped + def _rmatvec(self, y: NDArray) -> NDArray: + ncp = get_array_module(y) + if type(self.param) is not type(y): + self.param = to_cupy_conditional(y, self.param) + x = ncp.zeros(self.dims, dtype=self.dtype) + # ... adjoint implementation, x = A^H y (conjugate for complex params!) + return x diff --git a/.claude/skills/newop/reference/test_template.py b/.claude/skills/newop/reference/test_template.py new file mode 100644 index 00000000..d125f5b1 --- /dev/null +++ b/.claude/skills/newop/reference/test_template.py @@ -0,0 +1,57 @@ +"""Skeleton test for a new PyLops operator. + +Merge into the ``pytests/test_.py`` that matches the operator. +Keep the CuPy guard header identical to the one already in that file. +""" + +import os + +if int(os.environ.get("TEST_CUPY_PYLOPS", 0)): + import cupy as np + from cupy.testing import assert_array_almost_equal + + backend = "cupy" +else: + import numpy as np + from numpy.testing import assert_array_almost_equal + + backend = "numpy" +import pytest + +from pylops.basicoperators import MyOperator # noqa: F401 (adjust import) +from pylops.optimization.basic import lsqr +from pylops.utils import dottest + +par1 = {"ny": 11, "nx": 11, "imag": 0, "dtype": "float64"} # square real +par2 = {"ny": 21, "nx": 11, "imag": 0, "dtype": "float64"} # overdetermined real +par1j = {"ny": 11, "nx": 11, "imag": 1j, "dtype": "complex128"} # square complex +par2j = {"ny": 21, "nx": 11, "imag": 1j, "dtype": "complex128"} # overdet. complex + + +@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)]) +def test_MyOperator(par): + """Dot-test, forward and inversion for MyOperator""" + param = np.arange(par["nx"]) + 1.0 + par["imag"] * (np.arange(par["nx"]) + 1.0) + + Op = MyOperator(param, dtype=par["dtype"]) + assert dottest( + Op, + par["ny"], + par["nx"], + rtol=1e-6 if par["dtype"] in ("float64", "complex128") else 1e-3, + complexflag=0 if par["imag"] == 0 else 3, + backend=backend, + ) + + x = np.ones(par["nx"]) + par["imag"] * np.ones(par["nx"]) + y = Op * x + assert_array_almost_equal(y, Op.todense() @ x, decimal=6) + + xinv = lsqr(Op, y, x0=np.zeros_like(x), niter=300, show=0)[0] + assert_array_almost_equal(x, xinv, decimal=4) + + +def test_MyOperator_raises(): + """Check input validation of MyOperator""" + with pytest.raises(ValueError): + MyOperator(np.ones(5), dims=(4,)) diff --git a/AIPOLICY.md b/AIPOLICY.md index 04e00845..ab8f6420 100644 --- a/AIPOLICY.md +++ b/AIPOLICY.md @@ -36,5 +36,6 @@ More specifically, we currently provide: - ``.pi/prompts/optest.md`` / ``.claude/skills/optest``: a skill to increase the test coverage of an operator; +- ``.claude/skills/newop``: a skill to create a new operator from a mathematical description or a plain implementation of forward and adjoint from file or URL; 🤖🤖 **This Policy was written by humans and polished by AI** 🤖🤖 diff --git a/docs/source/api/index.rst b/docs/source/api/index.rst index 8084d114..246ec9cb 100755 --- a/docs/source/api/index.rst +++ b/docs/source/api/index.rst @@ -103,6 +103,7 @@ Signal processing Interp InterpCubicSpline Bilinear + Downsample2D FFT FFT2D FFTND diff --git a/examples/plot_downsample.py b/examples/plot_downsample.py new file mode 100644 index 00000000..872f8cc1 --- /dev/null +++ b/examples/plot_downsample.py @@ -0,0 +1,140 @@ +r""" +Downsampling +============ +This example shows how to use the +:py:class:`pylops.signalprocessing.Downsample2D` operator to reduce the size +of a 2-dimensional array along both of its directions. + +Downsampling is performed in two steps: an anti-aliasing Gaussian filter is +first applied to the input array, and the smoothed array is subsequently +subsampled by the required decimation factors. Whilst a naive subsampling of +the input array would fold any energy above the Nyquist wavenumber of the +coarse grid back onto the retained wavenumbers (i.e., aliasing), the Gaussian +filter removes such energy prior to decimation. + +As the operator is linear, its adjoint (and, more interestingly, its inverse) +can also be used to move back from the coarse to the fine grid; the latter +represents a very simple form of *super-resolution*. +""" + +import matplotlib.pyplot as plt +import numpy as np +from scipy import datasets + +import pylops + +plt.close("all") +np.random.seed(0) + +############################################################################### +# Let's start by creating a 2-dimensional input vector containing an image +# from the ``scipy.datasets`` family and downsample it by a factor of 4 in +# both directions. +x = datasets.face()[::2, ::2, 0].astype(np.float64) +nz, nx = x.shape + +Dop = pylops.signalprocessing.Downsample2D((nz, nx), factors=4) +y = Dop @ x + +print(Dop) +print(f"Model size: {Dop.dims}, Data size: {Dop.dimsd}") + +fig, axs = plt.subplots(1, 2, figsize=(10, 4)) +axs[0].imshow(x, cmap="gray") +axs[0].set_title(f"Original {Dop.dims}") +axs[0].axis("tight") +axs[1].imshow(y, cmap="gray") +axs[1].set_title(f"Downsampled {Dop.dimsd}") +axs[1].axis("tight") +plt.tight_layout() + +############################################################################### +# The role of the anti-aliasing Gaussian filter becomes evident if we take +# the input and output in the frequency domain. Note how a simple resampling +# of the input array (i.e., picking one every four samples in each direction) +# would lead to aliasing of the high wavenumbers, which is instead not present +# in the downsampled array. +Fop = pylops.signalprocessing.FFT2D((nz, nx), fftshift_after=True) +F1op = pylops.signalprocessing.FFT2D((nz // 4, nx // 4), fftshift_after=True) + +xf = Fop @ x +yf = F1op @ y +yf1 = F1op @ x[::4, ::4] + +fig, axs = plt.subplots(1, 3, figsize=(10, 4)) +axs[0].imshow( + np.abs(xf)[ + nz // 2 - nz // 8 : nz // 2 + nz // 8, nx // 2 - nx // 8 : nx // 2 + nx // 8 + ], + cmap="jet", + vmin=0, + vmax=0.005 * np.abs(xf).max(), +) +axs[0].set_title("Original (centered)") +axs[0].axis("tight") +axs[1].imshow(np.abs(yf), cmap="jet", vmin=0, vmax=0.005 * np.abs(yf).max()) +axs[1].set_title("Downsampled") +axs[1].axis("tight") +axs[2].imshow(np.abs(yf1), cmap="jet", vmin=0, vmax=0.005 * np.abs(yf1).max()) +axs[2].set_title("Resampled (no filter)") +axs[2].axis("tight") +plt.tight_layout() + +############################################################################### +# Similarly, if we take a synthetic image containing a rapidly oscillating +# pattern, where aliasing is easy to spot, we can see the difference between +# our downsampled image with the one obtained by simply picking one every +# four samples in each direction (which is equivalent to using ``sigma=0``). +nz1, nx1 = 201, 201 +iz, ix = np.meshgrid(np.arange(nz1), np.arange(nx1), indexing="ij") +xosc = np.sin(0.1 * np.sqrt((iz - nz1 // 2) ** 2 + (ix - nx1 // 2) ** 2) ** 2 / 10.0) + +Dop = pylops.signalprocessing.Downsample2D((nz1, nx1), factors=4) +Dop_noaa = pylops.signalprocessing.Downsample2D((nz1, nx1), factors=4, sigma=0.0) + +fig, axs = plt.subplots(1, 3, figsize=(12, 4)) +axs[0].imshow(xosc, cmap="gray") +axs[0].set_title("Original") +axs[0].axis("tight") +axs[1].imshow(Dop_noaa @ xosc, cmap="gray") +axs[1].set_title("Subsampled (aliased)") +axs[1].axis("tight") +axs[2].imshow(Dop @ xosc, cmap="gray") +axs[2].set_title("Downsampled (anti-aliased)") +axs[2].axis("tight") +plt.tight_layout() + +############################################################################### +# Finally, we consider the inverse problem: given the downsampled data, can we +# retrieve the original, finely sampled image? As the operator has many more +# columns than rows, this problem is heavily underdetermined and we must +# regularize it. Here we simply ask for a smooth solution by penalizing the +# Laplacian of the model. We compare the estimated model with the adjoint, +# which spreads each coarse sample back over the fine grid. +x = datasets.face()[400:528:2, 400:528:2, 0].astype(np.float64) +nz, nx = x.shape + +Dop = pylops.signalprocessing.Downsample2D((nz, nx), factors=2) +y = Dop @ x + +xadj = Dop.H @ y +D2op = pylops.Laplacian((nz, nx), weights=(1, 1), dtype="float64") +xinv = pylops.optimization.leastsquares.regularized_inversion( + Dop, y.ravel(), [D2op], epsRs=[np.sqrt(0.1)], **dict(iter_lim=200) +)[0] +xinv = xinv.reshape(nz, nx) + +fig, axs = plt.subplots(1, 4, figsize=(14, 4)) +axs[0].imshow(x, cmap="gray", vmin=0, vmax=255) +axs[0].set_title("Original") +axs[0].axis("tight") +axs[1].imshow(y, cmap="gray", vmin=0, vmax=255) +axs[1].set_title("Downsampled") +axs[1].axis("tight") +axs[2].imshow(xadj, cmap="gray") +axs[2].set_title("Adjoint") +axs[2].axis("tight") +axs[3].imshow(xinv, cmap="gray", vmin=0, vmax=255) +axs[3].set_title("Inverse") +axs[3].axis("tight") +plt.tight_layout() diff --git a/pylops/signalprocessing/__init__.py b/pylops/signalprocessing/__init__.py index c9d58262..de21b40a 100755 --- a/pylops/signalprocessing/__init__.py +++ b/pylops/signalprocessing/__init__.py @@ -18,6 +18,7 @@ Interp Interpolation operator. InterpCubicSpline Cubic Spline Interpolation operator. Bilinear Bilinear interpolation operator. + Downsample2D 2D downsampling operator. FFT One dimensional Fast-Fourier Transform. FFT2D Two dimensional Fast-Fourier Transform. FFTND N-dimensional Fast-Fourier Transform. @@ -59,6 +60,7 @@ from .interp import * from .interpspline import * from .bilinear import * +from .downsample2d import * from .radon2d import * from .radon3d import * from .fourierradon2d import * @@ -96,6 +98,7 @@ "Interp", "InterpCubicSpline", "Bilinear", + "Downsample2D", "Radon2D", "Radon3D", "FourierRadon2D", diff --git a/pylops/signalprocessing/downsample2d.py b/pylops/signalprocessing/downsample2d.py new file mode 100644 index 00000000..f811eb1d --- /dev/null +++ b/pylops/signalprocessing/downsample2d.py @@ -0,0 +1,235 @@ +__all__ = ["Downsample2D"] + +from typing import Literal + +import numpy as np + +from pylops import LinearOperator +from pylops.signalprocessing import Convolve2D +from pylops.utils._internal import _value_or_sized_to_tuple +from pylops.utils.backend import get_array_module, get_normalize_axis_index +from pylops.utils.decorators import reshaped +from pylops.utils.typing import DTypeLike, InputDimsLike, NDArray, SamplingLike + + +def _gaussian_kernel1d(sigma: float, truncate: float) -> NDArray: + """Create a normalized, symmetric 1d Gaussian kernel. + + The kernel is truncated at ``truncate`` standard deviations, leading to a + kernel of size :math:`2 \\lfloor \\text{truncate} \\sigma + 0.5 \\rfloor + 1`. + A unitary kernel (i.e., ``[1.]``) is returned when ``sigma=0``. + """ + if sigma == 0.0: + return np.ones(1) + radius = int(truncate * sigma + 0.5) + x = np.arange(-radius, radius + 1) + h = np.exp(-0.5 * (x / sigma) ** 2) + return h / h.sum() + + +class Downsample2D(LinearOperator): + r"""2D downsampling operator. + + Downsample a two (or more) dimensional array along a pair of ``axes`` by + applying an anti-aliasing Gaussian filter followed by subsampling with + a given decimation factor in each of the two directions. + + Parameters + ---------- + dims : :obj:`list` or :obj:`int` + Number of samples for each dimension. + factors : :obj:`int` or :obj:`tuple`, optional + Decimation factors along each of the two ``axes``. If a single value is + provided, the same factor is used in both directions. + sigma : :obj:`float` or :obj:`tuple`, optional + Standard deviations (in number of samples) of the Gaussian filter along + each of the two ``axes``. If a single value is provided, the same + standard deviation is used in both directions. If ``None``, the + standard deviations are set to ``(factor - 1) / 2`` for each direction. + truncate : :obj:`float`, optional + Number of standard deviations at which the Gaussian filter is + truncated. The filter has ``2 * int(truncate * sigma + 0.5) + 1`` + samples along each direction. + axes : :obj:`tuple`, optional + Axes along which downsampling is applied. + method : :obj:`str`, optional + Method used to calculate the Gaussian filtering (``auto``, ``direct`` + or ``fft``) - see :func:`scipy.signal.convolve` for details. + dtype : :obj:`str`, optional + Type of elements in input array. + name : :obj:`str`, optional + Name of operator (to be used by :func:`pylops.utils.describe.describe`) + + Attributes + ---------- + h : :obj:`numpy.ndarray` + 2d Gaussian filter applied prior to subsampling. + Cop : :obj:`pylops.signalprocessing.Convolve2D` + Gaussian filtering operator. + dims : :obj:`tuple` + Shape of the array after the adjoint, but before flattening. + + For example, ``x_reshaped = (Op.H * y.ravel()).reshape(Op.dims)``. + dimsd : :obj:`tuple` + Shape of the array after the forward, but before flattening. + + For example, ``y_reshaped = (Op * x.ravel()).reshape(Op.dimsd)``. + shape : :obj:`tuple` + Operator shape. + explicit : :obj:`bool` + Operator contains a matrix that can be solved explicitly (``True``) or + not (``False``). + + Raises + ------ + ValueError + If ``dims`` has less than 2 dimensions, if ``axes``, ``factors``, or + ``sigma`` do not contain 2 elements, if any element of ``factors`` is + smaller than 1 or larger than half the size of the corresponding axis, + or if any element of ``sigma`` is negative. + + See Also + -------- + pylops.signalprocessing.Convolve2D : 2D convolution operator + pylops.Restriction : Restriction (or sampling) operator + + Notes + ----- + The Downsample2D operator reduces the size of a two-dimensional array + :math:`\mathbf{x}` of size :math:`n_0 \times n_1` by a factor + :math:`f_0` and :math:`f_1` along the first and second direction, + respectively. Direct subsampling of the input array would however lead to + aliasing of any energy above the Nyquist wavenumber of the coarse grid; + for this reason the array is first smoothed by a separable Gaussian kernel + + .. math:: + h[p, q] = g_{\sigma_0}[p]\, g_{\sigma_1}[q], \qquad + g_\sigma[p] = \frac{e^{-p^2 / (2\sigma^2)}} + {\sum_{p'} e^{-p'^2 / (2\sigma^2)}} + + with :math:`|p| \leq r_0`, :math:`|q| \leq r_1`, and + :math:`r_i = \lfloor \tau \sigma_i + 0.5 \rfloor` where :math:`\tau` is the + ``truncate`` parameter. In forward mode, filtering and subsampling are + applied one after the other + + .. math:: + y[i, j] = \sum_{p=-r_0}^{r_0} \sum_{q=-r_1}^{r_1} + h[p, q] \, x[f_0 i - p, f_1 j - q] + \quad \forall i=0,\ldots,\lceil n_0 / f_0 \rceil - 1, + \; j=0,\ldots,\lceil n_1 / f_1 \rceil - 1 + + where the input array is assumed to be zero-padded outside of its + boundaries. Since the adjoint of subsampling is zero-interleaving and the + adjoint of convolution is correlation, in adjoint mode the data is first + spread over the fine grid and then correlated with the same kernel + + .. math:: + x[k, l] = \sum_{p=-r_0}^{r_0} \sum_{q=-r_1}^{r_1} + h[p, q] \, \tilde{y}[k + p, l + q], \qquad + \tilde{y}[k, l] = + \begin{cases} + y[k / f_0, l / f_1] & k \bmod f_0 = 0 \land l \bmod f_1 = 0\\ + 0 & \text{otherwise} + \end{cases} + + Note that, as the Gaussian kernel is real and symmetric, the operator + is effectively the composition of a self-adjoint smoothing operator and a + restriction operator. + + """ + + def __init__( + self, + dims: InputDimsLike, + factors: int | InputDimsLike = 2, + sigma: float | SamplingLike | None = None, + truncate: float = 4.0, + axes: InputDimsLike = (-2, -1), + method: Literal["auto", "direct", "fft"] | None = "fft", + dtype: DTypeLike = "float64", + name: str = "D", + ) -> None: + # check dims + dims = _value_or_sized_to_tuple(dims) + if len(dims) < 2: + msg = "dims must contain at least 2 dimensions" + raise ValueError(msg) + + # check axes + if len(axes) != 2: + msg = "axes must contain 2 elements" + raise ValueError(msg) + axes = tuple(get_normalize_axis_index()(ax, len(dims)) for ax in axes) + + # check factors + factors = _value_or_sized_to_tuple(factors, repeat=2) + if len(factors) != 2: + msg = "factors must contain 2 elements" + raise ValueError(msg) + + for f, ax in zip(factors, axes, strict=True): + if f < 1: + msg = "factors must be greater or equal to 1" + raise ValueError(msg) + if f > dims[ax] // 2: + msg = ( + f"factor={f} is larger than the half of the " + f"number of samples ({dims[ax]}) along axis={ax}" + ) + raise ValueError(msg) + + if sigma is None: + sigma = tuple((f - 1) / 2.0 for f in factors) + else: + sigma = _value_or_sized_to_tuple(sigma, repeat=2) + if len(sigma) != 2: + msg = "sigma must contain 2 elements" + raise ValueError(msg) + if any(s < 0 for s in sigma): + msg = "sigma must be positive" + raise ValueError(msg) + + self.axes = axes + self.factors = factors + self.sigma = sigma + self.truncate = truncate + + # data dimensions after subsampling + dimsd = list(dims) + for f, ax in zip(factors, axes, strict=True): + dimsd[ax] = int(np.ceil(dims[ax] / f)) + + super().__init__( + dtype=np.dtype(dtype), dims=dims, dimsd=tuple(dimsd), name=name + ) + + # separable gaussian filter and associated convolution operator + h0 = _gaussian_kernel1d(sigma[0], truncate) + h1 = _gaussian_kernel1d(sigma[1], truncate) + self.h = np.outer(h0, h1).astype(self.dtype) + self.Cop = Convolve2D( + dims, + h=self.h, + offset=(h0.size // 2, h1.size // 2), + axes=axes, + method=method, + dtype=dtype, + ) + + # slices used to subsample the filtered model + self.slices = tuple( + slice(None, None, factors[axes.index(ax)]) if ax in axes else slice(None) + for ax in range(len(dims)) + ) + + @reshaped + def _matvec(self, x: NDArray) -> NDArray: + y = self.Cop._matvec(x.ravel()).reshape(self.dims) + return y[self.slices] + + @reshaped + def _rmatvec(self, x: NDArray) -> NDArray: + ncp = get_array_module(x) + y = ncp.zeros(self.dims, dtype=self.dtype) + y[self.slices] = x + return self.Cop._rmatvec(y.ravel()).reshape(self.dims) diff --git a/pylops/signalprocessing/radon2d.py b/pylops/signalprocessing/radon2d.py index 77bd0313..97c75f7b 100644 --- a/pylops/signalprocessing/radon2d.py +++ b/pylops/signalprocessing/radon2d.py @@ -144,7 +144,7 @@ def Radon2D( taxis: NDArray, haxis: NDArray, pxaxis: NDArray, - kind: Literal["linear", "parabolic", "hyperbolic"] = "linear", + kind: Literal["linear", "parabolic", "hyperbolic"] | Callable = "linear", centeredh: bool = True, interp: bool = True, onthefly: bool = False, diff --git a/pylops/signalprocessing/radon3d.py b/pylops/signalprocessing/radon3d.py index 96eabd0f..cbc7a935 100644 --- a/pylops/signalprocessing/radon3d.py +++ b/pylops/signalprocessing/radon3d.py @@ -164,7 +164,7 @@ def Radon3D( hxaxis: NDArray, pyaxis: NDArray, pxaxis: NDArray, - kind: Literal["linear", "parabolic", "hyperbolic"] = "linear", + kind: Literal["linear", "parabolic", "hyperbolic"] | Callable = "linear", centeredh: bool = True, interp: bool = True, onthefly: bool = False, diff --git a/pytests/test_downsample.py b/pytests/test_downsample.py new file mode 100644 index 00000000..2bb7e56a --- /dev/null +++ b/pytests/test_downsample.py @@ -0,0 +1,156 @@ +import os + +if int(os.environ.get("TEST_CUPY_PYLOPS", 0)): + import cupy as np + from cupy.testing import assert_array_almost_equal + + backend = "cupy" +else: + import numpy as np + from numpy.testing import assert_array_almost_equal + + backend = "numpy" + +import numpy as npp +import pytest +from scipy.ndimage import gaussian_filter + +from pylops.optimization.basic import lsqr +from pylops.signalprocessing import Downsample2D +from pylops.utils import dottest + +par1 = { + "ny": 21, + "nx": 15, + "factors": 3, + "imag": 0, + "dtype": "float64", +} # same factor, real +par2 = { + "ny": 20, + "nx": 16, + "factors": (2, 4), + "imag": 0, + "dtype": "float64", +} # different factors, real +par3 = { + "ny": 11, + "nx": 13, + "factors": 1, + "imag": 0, + "dtype": "float64", +} # unitary factor, real +par1j = { + "ny": 21, + "nx": 15, + "factors": 3, + "imag": 1j, + "dtype": "complex128", +} # same factor, complex +par2j = { + "ny": 20, + "nx": 16, + "factors": (2, 4), + "imag": 1j, + "dtype": "complex128", +} # different factors, complex + + +@pytest.mark.parametrize( + "kwargs", + [ + {"dims": (10,)}, + {"dims": (10, 10), "axes": (0,)}, + {"dims": (10, 10), "factors": (2, 2, 2)}, + {"dims": (10, 10), "factors": 0}, + {"dims": (10, 10), "factors": 11}, + {"dims": (10, 10), "sigma": (1.0, 1.0, 1.0)}, + {"dims": (10, 10), "sigma": -1.0}, + ], +) +def test_Downsample2D_raises(kwargs): + """Check input validation of Downsample2D""" + with pytest.raises(ValueError): + Downsample2D(**kwargs) + + +def test_Downsample2D_sigma(): + """Check that a null sigma leads to pure subsampling""" + x = np.random.normal(0.0, 1.0, (12, 9)) + Dop = Downsample2D((12, 9), factors=(3, 3), sigma=0.0) + assert Dop.h.shape == (1, 1) + assert_array_almost_equal(Dop @ x, x[::3, ::3], decimal=10) + + +def test_Downsample2D_ndim(): + """Check that Downsample2D can be applied to a subset of axes of a + 3-dimensional array + """ + Dop = Downsample2D((7, 9, 5), factors=2, axes=(0, 1)) + assert Dop.dimsd == (4, 5, 5) + assert dottest(Dop, *Dop.shape, rtol=1e-6, backend=backend) + + +@pytest.mark.parametrize("par", [(par1), (par2), (par3), (par1j), (par2j)]) +def test_Downsample2D(par): + """Dot-test and shapes for Downsample2D""" + Dop = Downsample2D( + (par["ny"], par["nx"]), factors=par["factors"], dtype=par["dtype"] + ) + factors = ( + (par["factors"], par["factors"]) + if isinstance(par["factors"], int) + else par["factors"] + ) + assert Dop.dimsd == ( + int(npp.ceil(par["ny"] / factors[0])), + int(npp.ceil(par["nx"] / factors[1])), + ) + assert dottest( + Dop, + *Dop.shape, + rtol=1e-6, + complexflag=0 if par["imag"] == 0 else 3, + backend=backend, + ) + + +@pytest.mark.parametrize("par", [(par1), (par2), (par1j), (par2j)]) +def test_Downsample2D_scipy(par): + """Compare Downsample2D forward with scipy + gaussian filtering plus subsampling""" + factors = ( + (par["factors"], par["factors"]) + if isinstance(par["factors"], int) + else par["factors"] + ) + sigma = tuple((f - 1) / 2.0 for f in factors) + + shape = (par["ny"], par["nx"]) + x = np.random.normal(0.0, 1.0, shape) + par["imag"] * np.random.normal( + 0.0, 1.0, shape + ) + Dop = Downsample2D( + (par["ny"], par["nx"]), factors=par["factors"], dtype=par["dtype"] + ) + y = Dop @ x + + xnp = np.asnumpy(x) if backend == "cupy" else x + ynp = gaussian_filter(xnp, sigma=sigma, truncate=4.0, mode="constant")[ + :: factors[0], :: factors[1] + ] + assert_array_almost_equal(y, np.asarray(ynp), decimal=10) + + +@pytest.mark.parametrize("par", [(par3)]) +def test_Downsample2D_inverse(par): + """Invert Downsample2D when no decimation is applied (factors=1) as in + this case the operator is a square, invertible smoothing operator + """ + x = np.random.normal(0.0, 1.0, (par["ny"], par["nx"])) + Dop = Downsample2D( + (par["ny"], par["nx"]), factors=par["factors"], sigma=0.6, dtype=par["dtype"] + ) + y = Dop @ x + xinv = lsqr(Dop, y.ravel(), x0=np.zeros(Dop.shape[1]), niter=500, show=0)[0] + assert_array_almost_equal(x.ravel(), xinv, decimal=3) diff --git a/pytests/test_radon.py b/pytests/test_radon.py index 787fa635..8f299c73 100644 --- a/pytests/test_radon.py +++ b/pytests/test_radon.py @@ -126,6 +126,18 @@ def test_unknown_engine(): _ = Radon3D(None, None, None, None, None, engine="foo") +@pytest.mark.skipif( + int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled" +) +def test_Radon2D_unknown_kind(): + """Check error is raised if unknown (and non-callable) kind is passed""" + t = np.arange(11, dtype=np.float64) * 0.005 + h = np.arange(21, dtype=np.float64) + px = np.linspace(0, 2e-2, 21, dtype=np.float64) + with pytest.raises(NotImplementedError, match="Wrong kind of basis function"): + _ = Radon2D(t, h, px, kind="foo") + + @pytest.mark.skipif( int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled" ) @@ -189,6 +201,33 @@ def test_Radon2D(par, dtype): assert_array_almost_equal(x.ravel(), xinv, decimal=1) +@pytest.mark.skipif( + int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled" +) +def test_Radon2D_callable_kind(): + """Dot-test for Radon2D operator when kind is a custom callable""" + dt, dh = 0.005, 1 + t = np.arange(par1["nt"], dtype=np.float64) * dt + h = np.arange(par1["nhx"], dtype=np.float64) * dh + px = np.linspace(0, par1["pxmax"], par1["npx"], dtype=np.float64) + + def _linear(x, t, px): + return t + px * x + + Rop = Radon2D( + t, + h, + px, + centeredh=par1["centeredh"], + interp=par1["interp"], + kind=_linear, + onthefly=False, + engine="numpy", + dtype=np.float64, + ) + assert dottest(Rop, par1["nhx"] * par1["nt"], par1["npx"] * par1["nt"]) + + @pytest.mark.skipif( int(os.environ.get("TEST_CUPY_PYLOPS", 0)) == 1, reason="Not CuPy enabled" ) diff --git a/tutorials/ctscan.py b/tutorials/ctscan.py deleted file mode 100755 index 529017d2..00000000 --- a/tutorials/ctscan.py +++ /dev/null @@ -1,167 +0,0 @@ -r""" -16. CT Scan Imaging -=================== -This tutorial considers a very well-known inverse problem from the field of -medical imaging. - -First, we will be using the :class:`pylops.signalprocessing.Radon2D` operator -to model a *sinogram*, which is a graphic representation of the raw data -obtained from a CT scan. - -Note that whilst we can twick the Radon2D operator to work in a CT-like style, -this has initially been designed with other applications in mind -(i.e., seismic). We will see that if we use :class:`pylops.medical.CT2D` the produced -sinogram will be very similar in the middle (horizontal and near horizontal lines) but -it will greatly differ at both end (vertical and near vertical lines). The latter lines -are in fact not easy to parametrize using the convention chosen in Radon2D. - -The sinogram created by the :class:`pylops.medical.CT2D` operator is further -inverted using both a L2 solver and a TV-regularized solver like Split-Bregman. -""" - -import matplotlib.pyplot as plt - -# sphinx_gallery_thumbnail_number = 2 -import numpy as np -from numba import jit - -import pylops - -plt.close("all") -np.random.seed(10) - -############################################################################### -# Let's start by loading the Shepp-Logan phantom model. We can then construct -# the sinogram by providing a custom-made function to the -# :func:`pylops.signalprocessing.Radon2D` that samples parametric curves of -# such a type: -# -# .. math:: -# t(r,\theta; x) = \tan(90°-\theta)x + \frac{r}{\sin(\theta)} -# -# where :math:`\theta` is the angle between the x-axis (:math:`x`) and -# the perpendicular to the summation line and :math:`r` is the distance -# from the origin of the summation line. - - -@jit(nopython=True) -def radoncurve(x, r, theta): - return ( - (r - ny // 2) / (np.sin(theta) + 1e-15) - + np.tan(np.pi / 2.0 - theta) * x - + ny // 2 - ) - - -x = np.load("../testdata/optimization/shepp_logan_phantom.npy").T -x = x / x.max() -nx, ny = x.shape - -ntheta = 151 -theta = np.linspace(0.0, np.pi, ntheta, endpoint=False) - -RLop = pylops.signalprocessing.Radon2D( - np.arange(ny), - np.arange(nx), - theta, - kind=radoncurve, - centeredh=True, - interp=False, - engine="numba", - dtype="float64", -) - -y = RLop.H * x - -############################################################################### -# We can now first perform the adjoint, which in the medical imaging literature -# is also referred to as back-projection. -# -# This is the first step of a common reconstruction technique, named filtered -# back-projection, which simply applies a correction filter in the -# frequency domain to the adjoint model. -xrec = RLop * y - -fig, axs = plt.subplots(1, 3, figsize=(10, 4)) -axs[0].imshow(x.T, vmin=0, vmax=1, cmap="gray") -axs[0].set_title("Model") -axs[0].axis("tight") -axs[1].imshow(y.T, cmap="gray") -axs[1].set_title("Data") -axs[1].axis("tight") -axs[2].imshow(xrec.T, cmap="gray") -axs[2].set_title("Adjoint model") -axs[2].axis("tight") -fig.tight_layout() - - -############################################################################### -# Let's now repeat the same exercise, this time using the CT2D operator -Cop = pylops.medical.CT2D((ny, nx), 1.0, ny, theta, engine="cpu") - -y = Cop * x.T -xrec = Cop.H * y - -fig, axs = plt.subplots(1, 3, figsize=(10, 4)) -axs[0].imshow(x.T, vmin=0, vmax=1, cmap="gray") -axs[0].set_title("Model") -axs[0].axis("tight") -axs[1].imshow(np.flipud(y.T), cmap="gray") -axs[1].set_title("Data") -axs[1].axis("tight") -axs[2].imshow(xrec, cmap="gray") -axs[2].set_title("Adjoint model") -axs[2].axis("tight") -fig.tight_layout() - -############################################################################### -# Finally we take advantage of our different solvers and try to invert the -# modelling operator both in a least-squares sense and using TV-reg. -Dop = [ - pylops.FirstDerivative( - (ny, nx), axis=0, edge=True, kind="backward", dtype=np.float64 - ), - pylops.FirstDerivative( - (ny, nx), axis=1, edge=True, kind="backward", dtype=np.float64 - ), -] -D2op = pylops.Laplacian(dims=(ny, nx), edge=True, dtype=np.float64) - -# L2 -xinv_sm = pylops.optimization.leastsquares.regularized_inversion( - Cop, y.ravel(), [D2op], epsRs=[1e1], **dict(iter_lim=20) -)[0] -xinv_sm = np.real(xinv_sm.reshape(ny, nx)).T - -# TV -mu = 1.5 -lamda = [1.0, 1.0] -niter = 3 -niterinner = 4 - -xinv = pylops.optimization.sparsity.splitbregman( - Cop, - y.ravel(), - Dop, - niter_outer=niter, - niter_inner=niterinner, - mu=mu, - epsRL1s=lamda, - tol=1e-4, - tau=1.0, - show=False, - **dict(iter_lim=20, damp=1e-2), -)[0] -xinv = np.real(xinv.reshape(ny, nx)).T - -fig, axs = plt.subplots(1, 3, figsize=(10, 4)) -axs[0].imshow(x.T, vmin=0, vmax=1, cmap="gray") -axs[0].set_title("Model") -axs[0].axis("tight") -axs[1].imshow(xinv_sm.T, vmin=0, vmax=1, cmap="gray") -axs[1].set_title("L2 Inversion") -axs[1].axis("tight") -axs[2].imshow(xinv.T, vmin=0, vmax=1, cmap="gray") -axs[2].set_title("TV-Reg Inversion") -axs[2].axis("tight") -fig.tight_layout() From 4196155438205e0281b3f103829a226b97a6df0d Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 30 Aug 2026 14:50:58 +0100 Subject: [PATCH 2/3] minor: restore ctscan tutorial --- tutorials/ctscan.py | 167 ++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 167 insertions(+) create mode 100755 tutorials/ctscan.py diff --git a/tutorials/ctscan.py b/tutorials/ctscan.py new file mode 100755 index 00000000..529017d2 --- /dev/null +++ b/tutorials/ctscan.py @@ -0,0 +1,167 @@ +r""" +16. CT Scan Imaging +=================== +This tutorial considers a very well-known inverse problem from the field of +medical imaging. + +First, we will be using the :class:`pylops.signalprocessing.Radon2D` operator +to model a *sinogram*, which is a graphic representation of the raw data +obtained from a CT scan. + +Note that whilst we can twick the Radon2D operator to work in a CT-like style, +this has initially been designed with other applications in mind +(i.e., seismic). We will see that if we use :class:`pylops.medical.CT2D` the produced +sinogram will be very similar in the middle (horizontal and near horizontal lines) but +it will greatly differ at both end (vertical and near vertical lines). The latter lines +are in fact not easy to parametrize using the convention chosen in Radon2D. + +The sinogram created by the :class:`pylops.medical.CT2D` operator is further +inverted using both a L2 solver and a TV-regularized solver like Split-Bregman. +""" + +import matplotlib.pyplot as plt + +# sphinx_gallery_thumbnail_number = 2 +import numpy as np +from numba import jit + +import pylops + +plt.close("all") +np.random.seed(10) + +############################################################################### +# Let's start by loading the Shepp-Logan phantom model. We can then construct +# the sinogram by providing a custom-made function to the +# :func:`pylops.signalprocessing.Radon2D` that samples parametric curves of +# such a type: +# +# .. math:: +# t(r,\theta; x) = \tan(90°-\theta)x + \frac{r}{\sin(\theta)} +# +# where :math:`\theta` is the angle between the x-axis (:math:`x`) and +# the perpendicular to the summation line and :math:`r` is the distance +# from the origin of the summation line. + + +@jit(nopython=True) +def radoncurve(x, r, theta): + return ( + (r - ny // 2) / (np.sin(theta) + 1e-15) + + np.tan(np.pi / 2.0 - theta) * x + + ny // 2 + ) + + +x = np.load("../testdata/optimization/shepp_logan_phantom.npy").T +x = x / x.max() +nx, ny = x.shape + +ntheta = 151 +theta = np.linspace(0.0, np.pi, ntheta, endpoint=False) + +RLop = pylops.signalprocessing.Radon2D( + np.arange(ny), + np.arange(nx), + theta, + kind=radoncurve, + centeredh=True, + interp=False, + engine="numba", + dtype="float64", +) + +y = RLop.H * x + +############################################################################### +# We can now first perform the adjoint, which in the medical imaging literature +# is also referred to as back-projection. +# +# This is the first step of a common reconstruction technique, named filtered +# back-projection, which simply applies a correction filter in the +# frequency domain to the adjoint model. +xrec = RLop * y + +fig, axs = plt.subplots(1, 3, figsize=(10, 4)) +axs[0].imshow(x.T, vmin=0, vmax=1, cmap="gray") +axs[0].set_title("Model") +axs[0].axis("tight") +axs[1].imshow(y.T, cmap="gray") +axs[1].set_title("Data") +axs[1].axis("tight") +axs[2].imshow(xrec.T, cmap="gray") +axs[2].set_title("Adjoint model") +axs[2].axis("tight") +fig.tight_layout() + + +############################################################################### +# Let's now repeat the same exercise, this time using the CT2D operator +Cop = pylops.medical.CT2D((ny, nx), 1.0, ny, theta, engine="cpu") + +y = Cop * x.T +xrec = Cop.H * y + +fig, axs = plt.subplots(1, 3, figsize=(10, 4)) +axs[0].imshow(x.T, vmin=0, vmax=1, cmap="gray") +axs[0].set_title("Model") +axs[0].axis("tight") +axs[1].imshow(np.flipud(y.T), cmap="gray") +axs[1].set_title("Data") +axs[1].axis("tight") +axs[2].imshow(xrec, cmap="gray") +axs[2].set_title("Adjoint model") +axs[2].axis("tight") +fig.tight_layout() + +############################################################################### +# Finally we take advantage of our different solvers and try to invert the +# modelling operator both in a least-squares sense and using TV-reg. +Dop = [ + pylops.FirstDerivative( + (ny, nx), axis=0, edge=True, kind="backward", dtype=np.float64 + ), + pylops.FirstDerivative( + (ny, nx), axis=1, edge=True, kind="backward", dtype=np.float64 + ), +] +D2op = pylops.Laplacian(dims=(ny, nx), edge=True, dtype=np.float64) + +# L2 +xinv_sm = pylops.optimization.leastsquares.regularized_inversion( + Cop, y.ravel(), [D2op], epsRs=[1e1], **dict(iter_lim=20) +)[0] +xinv_sm = np.real(xinv_sm.reshape(ny, nx)).T + +# TV +mu = 1.5 +lamda = [1.0, 1.0] +niter = 3 +niterinner = 4 + +xinv = pylops.optimization.sparsity.splitbregman( + Cop, + y.ravel(), + Dop, + niter_outer=niter, + niter_inner=niterinner, + mu=mu, + epsRL1s=lamda, + tol=1e-4, + tau=1.0, + show=False, + **dict(iter_lim=20, damp=1e-2), +)[0] +xinv = np.real(xinv.reshape(ny, nx)).T + +fig, axs = plt.subplots(1, 3, figsize=(10, 4)) +axs[0].imshow(x.T, vmin=0, vmax=1, cmap="gray") +axs[0].set_title("Model") +axs[0].axis("tight") +axs[1].imshow(xinv_sm.T, vmin=0, vmax=1, cmap="gray") +axs[1].set_title("L2 Inversion") +axs[1].axis("tight") +axs[2].imshow(xinv.T, vmin=0, vmax=1, cmap="gray") +axs[2].set_title("TV-Reg Inversion") +axs[2].axis("tight") +fig.tight_layout() From 790ce31f39bcc3ff2a064b138c803855f743ee20 Mon Sep 17 00:00:00 2001 From: mrava87 Date: Sun, 30 Aug 2026 14:53:34 +0100 Subject: [PATCH 3/3] minor: removed unwanted step from newop skill --- .claude/skills/newop/SKILL.md | 2 -- 1 file changed, 2 deletions(-) diff --git a/.claude/skills/newop/SKILL.md b/.claude/skills/newop/SKILL.md index 08f2517a..0f1ae2e5 100644 --- a/.claude/skills/newop/SKILL.md +++ b/.claude/skills/newop/SKILL.md @@ -68,8 +68,6 @@ Use `reference/operator_template.py` as the skeleton. Key rules: `__init__` validates inputs), and a `Notes` section giving the maths of forward and adjoint in `.. math::` blocks. Match the level of detail of neighbouring operators. -- Add `.. versionadded:: ` to the class docstring for a brand-new - operator (check the current version in `pylops/version.py` / `pyproject.toml`). ## 3. Add tests