Feat: added Downsample2D operator

.claude/skills/newop/SKILL.md

@@ -0,0 +1,124 @@
+---
+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 <link>").
+---
+
+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__ = ["<Operator>"]`.
+- 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.<Operator>`.
+
+## 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.
+
+## 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_<file>.py -k <Operator> -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_<operator>.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).

.claude/skills/newop/reference/operator_template.py

@@ -0,0 +1,97 @@
+"""Skeleton for a new PyLops operator.
+
+Copy into ``pylops/<subpackage>/<operatorname>.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

.claude/skills/newop/reference/test_template.py

@@ -0,0 +1,57 @@
+"""Skeleton test for a new PyLops operator.
+
+Merge into the ``pytests/test_<subpackage>.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,))

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** 🤖🤖

docs/source/api/index.rst

@@ -103,6 +103,7 @@ Signal processing
     Interp
     InterpCubicSpline
     Bilinear
+    Downsample2D
     FFT
     FFT2D
     FFTND