`waveeqprocessing.Kirchhoff`: opt-in `paired=True` for zero-offset / common-offset data (one trace per source-receiver pair)

The proposal names Kirchhoff and the kernels _travsrcrec_kirch_[r]matvec and _ampsrcrec_kirch_[r]matvec; the draft implementation and tests are in linked PR #840, including pytests/test_kirchhoff.py and pytests/test_lsm.py. Read the draft PR and these tests first, then check how LSM forwards arguments to Kirchhoff. Done means paired source-receiver data work across the described modes and dimensions, existing unpaired behavior stays bit-identical, and the stated checks pass.

Edited 2026-11-06: Proposal points 3-4 and question 2 now describe the paired argument on the existing
kernels that the draft PR uses; the first version copied the kernels.

Motivation

Kirchhoff models every source with every receiver (ns * nr traces). Zero-offset and common-offset data
have one trace per source-receiver pair, for example ground-penetrating radar profiles. We migrate glacier GPR profiles (500–5000 traces per line) with your
operator, and today there are three ways to do zero offset:

• (A) all pairs + Restriction to the diagonal. Models n² traces to keep n: cost ∝ n² · ni,
  memory ∝ n² · nt. A 515-trace × 3963-sample profile ran for more than 45 min and used 27 GB before we
  stopped it.
• (B) VStack of one single-pair Kirchhoff per trace (your suggestion in #755). Correct for
  dynamic=False, but each sub-operator compiles its own numba kernels (0.2–0.4 s per trace on the
  first call), and with dynamic=True the result changes: each sub-operator sets its own
  maxdist = 1e-2 (max dist_srcs + max dist_recs), so the amplitudes depend on how the survey is split
  (5–25 % max difference below; identical once all sub-operators use the same epsilon).
• (C) mode='byot' with one two-way table and a single dummy receiver (nr = 1,
  trav = 2 * t_src). This already works and is fast, but only kinematically: dynamic=True needs
  separate tables, so it raises. You also have to build the tables yourself, which means calling the
  private Kirchhoff._traveltime_table.

Adjoint timings (2-D, constant velocity, analytic tables, engine='numba' serial as pylops' default, image
n x 200, nt = 400, after JIT; (D) is the proposal below). Every route gives the same result as (A),
except where noted:

traces	(A) all pairs + Restriction	(B) VStack (first call incl. JIT)	(C) byot, nr = 1	(D) paired=True
dynamic=False, 50	0.13 s	0.006 s (8.8 s)	0.002 s	0.003 s
dynamic=False, 100	1.05 s	0.016 s (18.3 s)	0.007 s	0.009 s
dynamic=False, 200	6.06 s	0.045 s (36.1 s)	0.022 s	0.029 s
dynamic=True, 50	0.43 s	0.063 s (19.0 s), differs by 5 %	not possible	0.014 s
dynamic=True, 100	3.91 s	0.25 s (38.2 s), differs by 14 %	not possible	0.051 s
dynamic=True, 200	30.2 s	0.97 s (77.0 s), differs by 25 %	not possible	0.203 s

With NUMBA_NUM_THREADS=8, dynamic=True, 200 traces: (A) 6.1 s, (B) 0.087 s (203 s first call),
(D) 0.043 s. Dependencies as in your uv.lock (uv sync --locked, as in build.yaml), Python 3.12.
(D) is the current version of the draft PR; its adjoint is about 40 % slower than with the copied kernels
of the first version (dynamic, 200 traces: 0.20 s against 0.14 s serial, 0.042 s against 0.030 s with 8
threads, timed alternately in one process); the forward is equally fast.

Proposal

Add a new keyword argument paired: bool = False, placed last so that positional calls are unaffected.
With paired=True, source i is paired with receiver i: srcs and recs must have the same size, and
the data have shape (ns, nt) instead of (ns, nr, nt). Everything else stays the same: modes
(analytic, eikonal, byot with a tuple or a single (ni, ns) table), dynamic, wavfilter,
aperture, angleaperture, numpy/numba, 2-D and 3-D, and LSM (which forwards its keyword arguments
to Kirchhoff).

In the current code the change touches only:

1. __init__: nsnr = ns if paired else ns * nr. This value replaces ns * nr in the Convolve1D
   size, in self.nsnr and in the 3-D six/rix zeros. Also dimsd = (ns, nt) when paired, and
   six/rix get one entry per pair (2-D). That is 6 modified lines. With paired=False each one
   computes exactly what it computes now.
2. Input checks, raised only when paired=True: ValueError if ns != nr, or if a single table does not
   have ns columns.
3. _register_multiplications: NotImplementedError for engine='cuda' when paired=True; with
   paired=True, the selected kernels are bound with functools.partial(kernel, paired=True). The CUDA
   helper is not modified. With a single byot table, the existing _trav_kirch_* kernels already loop
   over pairs, so they are used unchanged.
4. The four kernels _travsrcrec_kirch_[r]matvec and _ampsrcrec_kirch_[r]matvec get a last argument
   paired: bool = False (here the forward kinematic kernel; the other three change in the same way):

         trav_recs: NDArray,
+        paired: bool = False,
     ) -> NDArray:
+        # paired=True: only receiver isrc for source isrc (one trace per pair)
+        nrs, ipair = (1, 1) if paired else (nr, 0)
         for isrc in prange(ns):
             travisrc = trav_srcs[:, isrc]
-            for irec in range(nr):
-                travirec = trav_recs[:, irec]
+            for irec in range(nrs):
+                travirec = trav_recs[:, irec + ipair * isrc]
 ...
-                        y[isrc * nr + irec, itravii] += x[ii] * (1 - travdii)
-                        y[isrc * nr + irec, itravii + 1] += x[ii] * travdii
+                        y[isrc * nrs + irec, itravii] += x[ii] * (1 - travdii)
+                        y[isrc * nrs + irec, itravii + 1] += x[ii] * travdii

With paired=False the kernels are called with the same arguments as now, and the indices reduce to the
current ones (nrs = nr, ipair = 0).

In the operator the diff is +83 / −33 lines (20 of the added lines are docstring). The first
version of the draft PR added four copies of the kernels instead (276 lines). Codacy flagged them (8 new
Lizard issues: parameter count, cyclomatic complexity, length), so the PR now uses the argument.

Verification (on dev 06dc276)

We checked 112 configurations: 2-D and 3-D × analytic / eikonal / byot (tuple) / byot
(single table) × dynamic on and off × numpy / numba × (no apertures, no wavelet filter) and
(aperture=2, angleaperture=30, wavfilter=True in 2-D) × offset 0 and 2. The grid has dx != dz.
Each configuration was run serial and with NUMBA_NUM_THREADS=8:

• paired=False gives forward and adjoint bit-identical to unpatched pylops (np.array_equal), with the
  same dimsd.
• paired=True is bit-identical to Restriction(diagonal) @ Kirchhoff(all pairs), both forward and
  adjoint, for zero offset and common offset.
• paired=True passes dottest (rtol 1e-10).
• The three new errors are raised as intended.

Not tested: engine='cuda' (no GPU here). With paired=False the CUDA path receives the same arguments as
before.

The first version of the patch (copied kernels) was also bit-identical to the paired copy of your operator
that we run in production (zero offset, analytic and eikonal, dynamic on and off, both engines). The
current version gives bit-identical forward and adjoint to the first version in all 112 configurations.

Out of scope

• When srcs == recs, the source and receiver tables are computed twice. Sharing them would halve the
  eikonal time and the table memory. This is a possible follow-up and is not part of this change.
• Grid-spacing issues in byot + dynamic gradients and in the offset aperture: reported in #838. The paired option keeps pylops' current conventions, so the results stay identical to
  (A).

Questions

1. Is paired a good name, or would you prefer something like geometry="pairs"?
2. Is the paired argument on the existing kernels fine, or would you prefer separate kernels (they
   would need restructuring to stay within Codacy's Lizard limits)?
3. Should the kinematic-only route (C) also go into the docs?

A draft PR with the patch, tests in pytests/test_kirchhoff.py and pytests/test_lsm.py, a docstring note,
a CHANGELOG entry and a zero-offset section in the LSM tutorial is open: #840. We will adapt it to your
answers before marking it ready.

How AI was used

The analysis of the existing routes, the patch, the checks and the benchmark were done with Claude Code
(Anthropic) in a pair-programming setup: I chose the design (keyword instead of a new class, no CUDA
kernels in a first step), reviewed the diff and the results, and own them.