perf: optimize every custom periodogram implementation across all backends

.github/workflows/ci.yml

@@ -21,6 +21,11 @@ jobs:
         python: ["3.11", "3.12"]
     steps:
       - uses: actions/checkout@v4
+      # PySide6 (the GUI tests, run offscreen) needs the EGL/OpenGL system
+      # libraries, which the ubuntu-latest image no longer preinstalls.
+      - name: Install Qt system libraries (Linux)
+        if: runner.os == 'Linux'
+        run: sudo apt-get update && sudo apt-get install -y libegl1 libopengl0 libgl1
       - uses: astral-sh/setup-uv@v5
       - name: Create environment
         run: uv venv --python ${{ matrix.python }}

CHANGELOG.md

@@ -6,6 +6,51 @@ All notable changes to cuPeriod are documented here. The format is based on
 
 ## [Unreleased]
 
+### Performance
+
+- **Multicore `numba` CPU kernels for PDM, CE, String-Length, MHAOV, and TLS** (BLS
+  already had one). With the `[fast]` extra installed, `backend="cpu"`/`"auto"` now
+  resolve to them; warm-kernel speedups over the vectorized numpy paths on a 3k-point
+  curve: PDM ~300x, CE ~135x, MHAOV ~57x, TLS ~53x, String-Length ~23x, with parity to
+  the numpy results at or below ~1e-11.
+- **MHAOV is rewritten around harmonic trig sums.** Every entry of the normal equations
+  is analytically a trig sum, so the Gram matrix is now assembled from `C_m`/`S_m`
+  sums computed with the Chebyshev recurrence (one `cos`/`sin` evaluation regardless of
+  the harmonic order) instead of materializing the `(F, N, 2H+1)` design tensor:
+  ~2.7-3.2x faster on numpy/cupy/torch alike, an order of magnitude less transient
+  memory (no more multi-GB tensors at 1e5 points), and gemm-free by construction —
+  the Blackwell cuBLAS workaround branch is gone because nothing calls gemm anymore.
+- **The portable GLS direct path evaluates `cos`/`sin` once instead of six times** —
+  the base-grid sums for the weights and the weighted data share one angle matrix and
+  the doubled-frequency sums follow from the double-angle identities: ~1.8x on
+  torch:cpu and torch:cuda; the frequency chunk is auto-capped by the light-curve
+  length so device memory stays bounded regardless of `N`.
+- **PDM bins each point once** into `n_bins*n_covers` fine bins and regroups every
+  cover exactly from that histogram (vectorized, CUDA, and numba paths): ~2.8x on the
+  numpy path, one third of the shared-memory atomics in the CUDA kernel.
+- **Opt-in `precision="float32"` now reaches the CUDA `RawKernel`s** (BLS, PDM, CE,
+  TLS): on consumer GPUs, whose float64 throughput is 1/64 of float32, the
+  FLOP-bound searches speed up ~8.6x (BLS full segmented run 610 -> 71 ms; TLS
+  34 -> 4 ms on an RTX Blackwell card). float64 stays the default; PDM keeps float64
+  accumulators and CE counts are exact integers at any precision. The float32 BLS
+  guards empty box windows with an `ivar` floor scaled to the total inverse variance
+  (float32 cumsum noise would otherwise fabricate boxes).
+- **String-Length gains a real CUDA kernel**: one block per period bitonic-sorts the
+  (phase, index) pairs in shared memory — stable, matching the CPU tie order — with
+  no `(P, N)` intermediates (~75x over numpy; curves longer than the shared-memory
+  capacity fall back to the vectorized path). The torch path fuses its sort+gather
+  via `torch.sort(stable=True)`.
+- **Faster scatter shims everywhere**: numpy binning goes through buffered
+  `bincount` instead of unbuffered `np.add.at`, cupy through `cupyx.scatter_add`,
+  and the batch kernels scatter row-wise so torch reads a stride-0 view instead of
+  materializing `(P, N)` broadcast copies and flat indices.
+- **BLS uploads each light curve to the GPU once per run** (cupy and torch) via a
+  per-curve device cache shared by the period segments — previously every segment
+  re-uploaded the arrays and synced on a device-side `t.min()` — and returns all
+  seven per-period outputs in one stacked device-to-host copy. The GLS NUFFT paths
+  batch the base-grid pair (`w`, `w*y`) into a single `n_trans=2` transform, and the
+  one-shot cufinufft path assembles the power on the GPU (194 -> 290 curves/s).
+
 ### Added
 
 - **Multi-vendor GPU support via PyTorch and the Python array API.** All seven

docs/guide/backends.md

@@ -39,25 +39,25 @@ The selectors:
   - `cupy`
   - `numba` if installed, else `astropy`
 * - PDM
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 * - CE
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 * - String-Length
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 * - MHAOV
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 * - TLS
-  - `numpy`
+  - `numba` *(with `[fast]`)*, `numpy`
   - `cupy`
-  - `numpy`
+  - `numba` if installed, else `numpy`
 ```
 
 † BLS's `numpy` backend is a GPU-parity *reference*, not the product path — it shares one
@@ -172,11 +172,15 @@ How to read this:
   search). They're so fast on the CPU that the GPU's marginal gain on a single curve is
   modest — but the GPU still wins decisively for **large grids and big catalogs**
   ({doc}`batch`).
-- **PDM, CE, String-Length, MHAOV, TLS run on numpy on the CPU**, so the GPU's
-  data-parallelism delivers roughly 50–190× on a single curve. If you're searching many
-  periods or many stars with these, the GPU is a large win. The portable **`torch:cuda`**
-  column tracks alongside — competitive for these methods, and the same code reaches
-  AMD/Intel/Apple GPUs (see the portable-backend note above).
+- The **CPU times in the table are the vectorized numpy paths**. With the `[fast]` extra
+  installed, PDM, CE, String-Length, MHAOV, and TLS now run multicore **numba** kernels
+  on `"cpu"`/`"auto"` — one to two orders of magnitude faster than the numpy column
+  (e.g. PDM ~300×, CE ~135×, MHAOV ~57× on a 3k-point curve) — which puts a warm CPU
+  within reach of the GPU for a *single* curve. The GPU remains the throughput champion
+  for large grids and catalogs.
+- On **consumer NVIDIA cards** (GeForce), whose float64 throughput is 1/64 of float32,
+  the opt-in `precision="float32"` runs the BLS/TLS CUDA kernels ~8-9× faster at
+  detection-grade accuracy; float64 stays the default.
 - cuPeriod's **CPU** path already beats the established reference tools it was checked
   against (GLS ~3× astropy, PDM ~4× PyAstronomy, BLS ~18× astropy's `BoxLeastSquares`).
 

docs/guide/methods.md

@@ -83,7 +83,8 @@ for peak in pg.best_periods(5, alias_diverse=True):     # alias-aware for box se
 ```
 
 With the `[fast]` extra, BLS's CPU backend is a multicore `numba` search ~20× faster than
-astropy. Key settings ({class}`~cuperiod.BLSSettings`): `min_period_days` /
+astropy (PDM, CE, String-Length, MHAOV, and TLS gain `numba` CPU kernels too — see
+{doc}`backends`). Key settings ({class}`~cuperiod.BLSSettings`): `min_period_days` /
 `max_period_days`, `duration_min_frac` / `duration_max_frac`, `n_durations`, `objective`.
 
 ### TLS — transit least squares

docs/index.md

@@ -84,7 +84,7 @@ All seven share one API, one CLI, and the full single/batch machinery. See
 pip install cuperiod            # CPU (numpy, scipy, astropy, finufft)
 pip install "cuperiod[gpu]"     # + CUDA 12 GPU backends (cupy, cufinufft)
 pip install "cuperiod[torch]"   # + portable PyTorch backend (AMD/Intel/Apple GPUs + CPU)
-pip install "cuperiod[fast]"    # + numba multicore box search (~20× astropy BLS on CPU)
+pip install "cuperiod[fast]"    # + numba multicore CPU kernels (all methods, 20-300×)
 pip install "cuperiod[gui]"     # + interactive desktop GUI (cuperiod-gui)
 ```
 

docs/installation.md

@@ -19,23 +19,27 @@ CPU, the command line, and batch processing over a process pool.
 | --- | --- | --- |
 | **gpu** | `pip install "cuperiod[gpu]"` | CUDA-12 GPU backends (`cupy-cuda12x`, `cufinufft`, the NVIDIA runtime wheels) |
 | **torch** | `pip install "cuperiod[torch]"` | the portable **PyTorch** backend — runs every method on AMD (ROCm), Intel (XPU), Apple (MPS), and a CPU path |
-| **fast** | `pip install "cuperiod[fast]"` | a multicore `numba` box search — BLS's CPU default, ~20× faster than astropy's compiled `BoxLeastSquares` |
+| **fast** | `pip install "cuperiod[fast]"` | multicore `numba` CPU kernels for BLS, PDM, CE, String-Length, MHAOV, and TLS — the CPU default when installed, one to two orders of magnitude faster than the fallback CPU paths |
 | **gui** | `pip install "cuperiod[gui]"` | the interactive desktop GUI, `cuperiod-gui` (PySide6 + pyqtgraph) — see {doc}`guide/gui` |
 | **pandas** | `pip install "cuperiod[pandas]"` | pandas `DataFrame` ingestion |
 
 Extras combine, e.g. `pip install "cuperiod[gpu,fast]"` or `"cuperiod[gui,torch]"`.
 
 :::{tip}
-The `[fast]` extra is worth installing even without a GPU: it makes BLS an order of
-magnitude faster on the CPU while matching astropy to floating point. When present it
-becomes BLS's default CPU backend automatically; otherwise BLS falls back to astropy.
+The `[fast]` extra is worth installing even without a GPU: it gives BLS, PDM, CE,
+String-Length, MHAOV, and TLS multicore JIT kernels that are one to two orders of
+magnitude faster than the fallback CPU paths while matching them to floating point.
+When present they become the default CPU backends automatically.
 :::
 
 ### What the `[fast]` extra changes
 
-Installing `numba` flips `backend="cpu"` for BLS from astropy's `BoxLeastSquares` to the
-in-house multicore box search. Nothing else about your code changes — the results match
-astropy to round-off; they just arrive ~20× sooner. (Note that `numba` currently caps
+Installing `numba` flips `backend="cpu"` (and the CPU fallback of `"auto"`) from the
+fallback implementations — astropy's `BoxLeastSquares` for BLS, the vectorized numpy
+kernels for PDM, CE, String-Length, MHAOV, and TLS — to in-house multicore JIT kernels.
+Nothing else about your code changes: the results match the fallbacks to round-off
+(BLS still matches astropy), they just arrive much sooner (~20× for BLS over astropy;
+~25–300× for the others over their numpy paths). (Note that `numba` currently caps
 `numpy < 2.5`, so installing it may downgrade numpy slightly.)
 
 ## GPU requirements

pyproject.toml

@@ -162,6 +162,8 @@ module = [
     "cufinufft",
     "cupy",
     "cupy.*",
+    "cupyx",
+    "cupyx.*",
     "numba",
     "numba.*",
     "pandas",

src/cuperiod/core/_arrayapi.py

@@ -93,7 +93,10 @@ def scatter_add(target: Any, index: Any, values: Any) -> None:
 
     ``target`` is 1-D and ``index``/``values`` are 1-D and aligned. ``numpy.add.at`` is
     not part of the array-API standard and PyTorch has no equivalent *function*, so this
-    dispatches per backend: torch uses ``Tensor.index_add_``; numpy/cupy use ``add.at``.
+    dispatches per backend: torch uses ``Tensor.index_add_``; cupy uses
+    ``cupyx.scatter_add`` (its fast documented scatter); numpy uses ``bincount``, which
+    computes the same float64 sums in one buffered pass — ``np.add.at`` is unbuffered
+    and an order of magnitude slower on large inputs.
 
     ``values.dtype`` must equal ``target.dtype``: torch ``index_add_`` rejects a
     mismatch (numpy/cupy would silently cast), so callers build both at the same
@@ -102,12 +105,72 @@ def scatter_add(target: Any, index: Any, values: Any) -> None:
     if is_torch_array(target):
         target.index_add_(0, index.long(), values)
         return
+    if is_cupy_array(target):
+        import cupyx
+
+        cupyx.scatter_add(target, index, values)
+        return
+    binned = np.bincount(
+        np.asarray(index), weights=np.asarray(values), minlength=target.size
+    )
+    target += binned.astype(target.dtype, copy=False)
+
+
+def scatter_add_rows(target: Any, index: Any, values: Any) -> None:
+    """In-place ``target[p, index[p, j]] += values[j]`` with repeats accumulated.
+
+    The batch-kernel binning primitive: ``target`` is 2-D ``(P, W)``, ``index`` is
+    ``(P, N)`` int64, and ``values`` is one shared row ``(N,)`` (the usual case — the
+    same light curve binned at ``P`` trial periods) or a full ``(P, N)``. Compared to
+    flattening and calling :func:`scatter_add`, this avoids materializing the broadcast
+    ``(P, N)`` values *and* the flat ``(P, N)`` int64 index on torch — the scatter
+    reads a stride-0 expanded view — and gives numpy one fused ``bincount`` pass.
+    """
+    if is_torch_array(target):
+        src = values if values.ndim == 2 else values.unsqueeze(0).expand_as(index)
+        target.scatter_add_(1, index, src)
+        return
+    if is_cupy_array(target):
+        import cupy
+        import cupyx
+
+        rows = cupy.arange(target.shape[0])[:, None]
+        cupyx.scatter_add(target, (rows, index), values)
+        return
+    n_rows, n_cols = target.shape
+    flat = index + (np.arange(n_rows, dtype=np.int64) * n_cols)[:, None]
+    weights = np.broadcast_to(values, index.shape)
+    binned = np.bincount(
+        flat.reshape(-1), weights=weights.reshape(-1), minlength=target.size
+    )
+    target += binned.reshape(target.shape).astype(target.dtype, copy=False)
+
+
+def scatter_counts_rows(target: Any, index: Any) -> None:
+    """In-place ``target[p, index[p, j]] += 1`` — a per-row histogram count.
+
+    Same layout contract as :func:`scatter_add_rows` but without a values array:
+    numpy uses weightless ``bincount`` and torch scatters a stride-0 view of a single
+    one, so no ``(P, N)`` ones array is ever built. ``target`` is float (counts are
+    accumulated in the working dtype for the entropy/variance math downstream).
+    """
+    if is_torch_array(target):
+        import torch
+
+        one = torch.ones((1, 1), dtype=target.dtype, device=target.device)
+        target.scatter_add_(1, index, one.expand_as(index))
+        return
     if is_cupy_array(target):
         import cupy
+        import cupyx
 
-        cupy.add.at(target, index, values)
+        rows = cupy.arange(target.shape[0])[:, None]
+        cupyx.scatter_add(target, (rows, index), target.dtype.type(1))
         return
-    np.add.at(target, index, values)
+    n_rows, n_cols = target.shape
+    flat = index + (np.arange(n_rows, dtype=np.int64) * n_cols)[:, None]
+    binned = np.bincount(flat.reshape(-1), minlength=target.size)
+    target += binned.reshape(target.shape).astype(target.dtype, copy=False)
 
 
 def to_host(a: Any) -> np.ndarray:
@@ -201,6 +264,8 @@ def to_device_array(host: np.ndarray, *, device: str, dtype: Any) -> Any:
     "resolve_precision",
     "resolve_torch_device",
     "scatter_add",
+    "scatter_add_rows",
+    "scatter_counts_rows",
     "to_device_array",
     "to_host",
 ]

src/cuperiod/core/backend.py

@@ -193,7 +193,10 @@ def ensure_shared_memory(
             f"the device allows at most {optin // 1024} KB; reduce {hint}, or use "
             "backend='cpu'."
         )
-    if int(dynamic_bytes) > _DEFAULT_SHARED_MEM:
+    # The opt-in is needed whenever the *total* (dynamic + static) exceeds the
+    # default cap — a launch with 48 KB dynamic still fails if the kernel also has
+    # static shared arrays.
+    if needed > _DEFAULT_SHARED_MEM:
         kernel.max_dynamic_shared_size_bytes = int(dynamic_bytes)
 
 

src/cuperiod/core/config.py

@@ -165,7 +165,7 @@ def _check_bounds(self) -> Self:
         default=20, ge=3, description="Skip if fewer finite points."
     )
     backend: Literal[
-        "auto", "cpu", "gpu", "numpy", "astropy", "cupy", "torch"
+        "auto", "cpu", "gpu", "numba", "numpy", "astropy", "cupy", "torch"
     ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=2048, ge=1, description="Trial periods per vectorized batch."
@@ -213,9 +213,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=3, description="Skip if fewer finite points."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=2048, ge=1, description="Trial periods per vectorized batch."
     )
@@ -261,9 +261,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=5, description="Skip if fewer finite points (need > 2H+1)."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=512, ge=1, description="Trial frequencies per vectorized batch."
     )
@@ -308,9 +308,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=3, description="Skip if fewer finite points."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=1024, ge=1, description="Trial periods per vectorized batch."
     )
@@ -353,9 +353,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=3, description="Skip if fewer finite points."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     batch_periods: int = Field(
         default=1024, ge=1, description="Trial periods per vectorized batch."
     )
@@ -425,9 +425,9 @@ def _check_bounds(self) -> Self:
     min_detections: int = Field(
         default=20, ge=3, description="Skip if fewer finite points."
     )
-    backend: Literal["auto", "cpu", "gpu", "numpy", "cupy", "torch"] = Field(
-        default="auto", description="Compute backend."
-    )
+    backend: Literal[
+        "auto", "cpu", "gpu", "numba", "numpy", "cupy", "torch"
+    ] = Field(default="auto", description="Compute backend.")
     period_batch: int = Field(
         default=256, ge=1, description="Trial periods per vectorized batch."
     )