More xp scipy

CHANGELOG.md

@@ -8,14 +8,22 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 ## [Unreleased]
 
 ### Added
+- MIT license (`LICENSE.txt`, referenced by `pyproject.toml`).
 - `xp.optimize.newton`: Newton's method or the secant method on many
   independent scalar equations at once, with the steps of
   `scipy.optimize.newton` for an array `x0`, on the array's backend (no host
   copies, fixed shapes, one synchronization per iteration).
 - The fake CuPy has `cupyx.scipy`: SciPy behind the CuPy rules (host arrays
   rejected, fake device arrays and proxy objects returned), so `xp.scipy` code
-  runs its CuPy path without a GPU. `cupyx.scipy.interpolate` has only the
-  names CuPy provides. `cupyx.empty_pinned`/`zeros_pinned` return NumPy arrays.
+  runs its CuPy path without a GPU. Every `cupyx.scipy` subpackage has only
+  the names CuPy v14.2.0 provides, so e.g. `special.jv` fails as on a GPU.
+  `cupyx.empty_pinned`/`zeros_pinned` return NumPy arrays.
+- `xp.optimize.fsolve`, `root`, `minimize` and `xp.integrate.quad`, `odeint`
+  (`cupyx.scipy` has neither): SciPy on the host with SciPy's arguments, the
+  callbacks get their arrays on the backend of `x0`/`y0`, and results come
+  back on that backend.
+- `xp.scipy.linalg.solve_circulant` on the CuPy backend, which `cupyx.scipy`
+  lacks: SciPy's FFT solve and arguments, on the device.
 - Solvers guide: splines on the device with `xp.scipy.interpolate`, and
   batched root finding.
 - Compiled Pyccel regression coverage for F-ordered arrays and column-block

LICENSE.txt

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2025 Max Lindqvist
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

docs/source/api.md

@@ -49,7 +49,8 @@ The submodules are named so that they do not hide a NumPy name (`rng`, not
 | [`cunumpy.cuda`](reference/cuda.md) | CUDA only | device selection and memory, `stream`, streams/events, `pin_memory`, debug mode, CUDA headers and source tools (`cuda_include_dir`, `parse_cuda_signature`) |
 | [`cunumpy.rng`](reference/rng.md) | both | `random_streams`, `get_rng`, `philox_*` |
 | [`cunumpy.algorithms`](reference/algorithms.md) | both | `morton_*`, `sort_by_key`, `segment_sum`, `compact_by_mask` |
-| [`cunumpy.optimize`](reference/optimize.md) | both | `newton` (batched Newton/secant), `NewtonResult` |
+| [`cunumpy.optimize`](reference/optimize.md) | both | `newton` (batched Newton/secant), `NewtonResult`; `fsolve`, `root`, `minimize` (SciPy on the host, callbacks on the device) |
+| [`cunumpy.integrate`](reference/integrate.md) | both | `quad`, `odeint` (SciPy on the host, callbacks on the device) |
 | [`cunumpy.mpi`](reference/mpi.md) | both | `mpi_buffer`, `MPIStaging`, CUDA-aware MPI detection, `synchronize_for_mpi`, and the MPI stand-ins of maybempi (`get_mpi`, `local_rank`, ...) |
 | [`cunumpy.profiling`](reference/profiling.md) | both | `timed_region`, `nvtx_range`, `count_transfers`, `assert_no_transfers`, `TransferBudget` |
 | [`cunumpy.memory`](reference/memory.md) | both | `HostStaging`, `HostCopy`, `DeviceMirror` |

docs/source/guides/solvers.md

@@ -44,6 +44,12 @@ def periodic_poisson(rho, length):
 * `cupyx.scipy` has only part of SciPy. A missing name raises
   `AttributeError` saying which backend lacks it; check with
   `xp.scipy.special.available("erfcx")` where a fallback is possible.
+  Commonly missing: `special.jv` (only `j0`, `j1`, `yn`), `special.erfi`,
+  `linalg.solve`/`inv` (use `xp.linalg`), `RectBivariateSpline`, and all of
+  `integrate` and `optimize` (see `xp.optimize` and `xp.integrate` below).
+* cunumpy fills in `xp.scipy.linalg.solve_circulant` on the CuPy backend: an
+  FFT solve on the device, with SciPy's arguments (a circulant preconditioner
+  per iteration stays on the GPU).
 * Keyword arguments can differ between the two: SciPy's `cg` takes `rtol`
   (since SciPy 1.12), CuPy's takes `tol`. Pass tolerances only after checking
   both, or through a small wrapper.
@@ -54,8 +60,8 @@ def periodic_poisson(rho, length):
   (`erf`, `erfc`, Bessel functions `i0`, `i1`, `k0`, ...) are in
   `xp.scipy.special` on both backends.
 * On the fake CuPy (tests without a GPU), `xp.scipy` runs SciPy with the CuPy
-  rules: host arrays are rejected, results are fake device arrays, and
-  `interpolate` has only the names CuPy has.
+  rules: host arrays are rejected, results are fake device arrays, and each
+  subpackage has only the names CuPy has.
 
 ## Splines and interpolation
 
@@ -87,7 +93,7 @@ transfers.
 
 SciPy's optimizers run on the host, and `cupyx.scipy` has no `optimize`.
 Sequential scalar problems (one `fsolve`, a `quad` integral, an ODE) belong on
-the host at setup (see `xp.setup_on_host`). Many independent scalar equations
+the host at setup (see the next section and `xp.setup_on_host`). Many independent scalar equations
 (a flux surface on every ray, an inverse mapping at every marker) are one
 array problem: `xp.optimize.newton` runs Newton's method (with `fprime`) or
 the secant method on all of them at once, with the steps of
@@ -104,6 +110,30 @@ r = xp.optimize.newton(residual, xp.full(levels.shape, 0.3))
 Each iteration synchronizes once to test convergence.
 `full_output=True` returns `root`, `converged` and `zero_der` arrays.
 
+## SciPy solvers with device callbacks
+
+For a single problem, `xp.optimize.fsolve`, `root` and `minimize` and
+`xp.integrate.quad` and `odeint` run SciPy on the host with SciPy's
+arguments, so they work on both backends. With `x0` (`y0` for `odeint`) on
+the device, the callbacks get `x` on the device and the result comes back on
+the device, so a residual written with `xp` that uses device data needs no
+host copy of that data:
+
+```python
+B = equilibrium.b_field_on_grid()  # device array
+
+def residual(x):
+    return interpolate(B, x) - b_target
+
+x = xp.optimize.fsolve(residual, xp.asarray([0.5, 0.0]))  # device array
+```
+
+Every evaluation copies `x` to the device and the result to the host, so use
+them at setup, not in a time loop. `quad` calls its integrand with a Python
+float; a 0-d device array returned is copied to the host. Bounds and
+constraint matrices of `minimize` are host data; constraint functions get
+device arrays like `fun`.
+
 ## Fused elementwise updates: `xp.kernels.fuse`
 
 On the GPU, `(gamma - 1) * (E - 0.5 * rho * u**2)` runs as five kernels, each

docs/source/reference/index.md

@@ -15,6 +15,7 @@ cuda
 rng
 algorithms
 optimize
+integrate
 mpi
 profiling
 memory

docs/source/reference/integrate.md

@@ -0,0 +1,8 @@
+# Integration: `cunumpy.integrate`
+
+SciPy's `quad` and `odeint` with callbacks on the device.
+
+```{eval-rst}
+.. automodule:: cunumpy.integrate
+   :members:
+```

docs/source/reference/optimize.md

@@ -1,6 +1,7 @@
-# Root finding: `cunumpy.optimize`
+# Optimization: `cunumpy.optimize`
 
-Batched Newton and secant iterations on either backend.
+Batched Newton and secant iterations on the device, and SciPy's `fsolve`,
+`root` and `minimize` with callbacks on the device.
 
 ```{eval-rst}
 .. automodule:: cunumpy.optimize

src/cunumpy/LLM_GUIDE.md

@@ -100,6 +100,8 @@ https://max-models.github.io/cunumpy/ and in `docs/source/` of the repository.
 | CUDA struct from a Pyccel argument class | `xp.arguments.CudaStruct.from_signature(Cls.__init__, "Name")`, `xp.arguments.write_cuda_header(...)` |
 | SciPy (sparse, sparse.linalg, fft, special, ndimage, ...) on either backend | `xp.scipy.<subpackage>.<name>` (SciPy or `cupyx.scipy`); `xp.scipy.special.available(name)` |
 | splines on the device | `xp.scipy.interpolate.UnivariateSpline`, `make_interp_spline`, `NdBSpline` (no `RectBivariateSpline`: `make_interp_spline` per axis + `NdBSpline`) |
+| circulant solve (FFT preconditioner) on the device | `xp.scipy.linalg.solve_circulant(c, b)` (cunumpy fill-in on CuPy, SciPy's arguments) |
+| one root/minimization/integral/ODE at setup with an `xp` callback using device data | `xp.optimize.fsolve`/`root`/`minimize`, `xp.integrate.quad`/`odeint` (SciPy on the host, callback args and results on `x0`'s backend) |
 | many independent scalar root problems (one per ray, per marker) | `xp.optimize.newton(func, x0_array, fprime=None)` (SciPy's array Newton/secant, on the device) |
 | chain of elementwise operations as one GPU kernel | `@xp.kernels.fuse` (`cupy.fuse` for CuPy arrays, plain call otherwise) |
 | PETSc solve on device arrays without copies | `xp.petsc.petsc_vec(array)` (CUDA/HIP petsc4py for CuPy arrays); `xp.synchronize()` around PETSc calls |

src/cunumpy/__init__.py

@@ -19,6 +19,7 @@
     algorithms,
     arguments,
     cuda,
+    integrate,
     kernels,
     memory,
     mpi,
@@ -110,6 +111,7 @@ def require_version(minimum: str) -> None:
     "get_array_module",
     "get_backend",
     "host_call",
+    "integrate",
     "is_cpu",
     "is_gpu",
     "kernels",

src/cunumpy/__init__.pyi

@@ -11,6 +11,7 @@ from numpy import *
 from cunumpy import algorithms as algorithms
 from cunumpy import arguments as arguments
 from cunumpy import cuda as cuda
+from cunumpy import integrate as integrate
 from cunumpy import kernels as kernels
 from cunumpy import memory as memory
 from cunumpy import mpi as mpi

src/cunumpy/_fake_cupy.py

@@ -16,7 +16,7 @@
   :class:`~cunumpy.arguments.CudaStruct` packing and the argument checks of
   :class:`~cunumpy.kernels.CudaKernel` work;
 * ``cupyx.scipy`` forwards to SciPy under the same rules, with the names of
-  ``cupyx.scipy.interpolate`` only (see :mod:`cunumpy._fake_cupyx`);
+  ``cupyx.scipy`` (see :mod:`cunumpy._fake_cupyx`);
 * CUDA kernels cannot run: ``RawKernel`` and friends raise
   ``NotImplementedError`` when called, and :func:`cunumpy.kernel_testing.requires_cupy`
   skips tests while the fake is active. Inside

src/cunumpy/_fake_cupyx.py

@@ -9,10 +9,10 @@
   in CuPy;
 * array results are fake device arrays, and objects (a fitted spline, a
   factorization) are proxies whose methods follow the same rules;
-* ``cupyx.scipy.interpolate`` has only the names CuPy provides
-  (:data:`INTERPOLATE`), so code using a SciPy name that CuPy lacks (e.g.
-  ``RectBivariateSpline``) fails here as it would on a GPU. The other
-  subpackages forward every SciPy name.
+* each subpackage has only the names CuPy provides (:data:`NAMES`), so code
+  using a SciPy name that CuPy lacks (e.g. ``RectBivariateSpline``,
+  ``special.jv``, ``linalg.solve_circulant``) fails here as it would on a
+  GPU. A CuPy name that the installed SciPy lacks is missing too.
 
 ``cupyx.empty_pinned`` and ``cupyx.zeros_pinned`` return NumPy arrays, like
 CuPy's (which are backed by pinned memory).
@@ -43,41 +43,139 @@
     "stats",
 )
 
-#: The names of ``cupyx.scipy.interpolate`` (CuPy reference, stable).
-INTERPOLATE = frozenset(
-    {
-        "Akima1DInterpolator",
-        "BPoly",
-        "BSpline",
-        "BarycentricInterpolator",
-        "CloughTocher2DInterpolator",
-        "CubicHermiteSpline",
-        "CubicSpline",
-        "InterpolatedUnivariateSpline",
-        "KroghInterpolator",
-        "LSQUnivariateSpline",
-        "LinearNDInterpolator",
-        "NdBSpline",
-        "NdPPoly",
-        "NearestNDInterpolator",
-        "PPoly",
-        "PchipInterpolator",
-        "RBFInterpolator",
-        "RegularGridInterpolator",
-        "UnivariateSpline",
-        "barycentric_interpolate",
-        "interp1d",
-        "interpn",
-        "krogh_interpolate",
-        "make_interp_spline",
-        "make_lsq_spline",
-        "pchip_interpolate",
-        "splantider",
-        "splder",
-    },
-)
 
-_ALLOWED = {"interpolate": INTERPOLATE}
+def _names(text: str) -> frozenset[str]:
+    return frozenset(text.split())
+
+
+#: The public names of each ``cupyx.scipy`` subpackage (CuPy v14.2.0).
+NAMES = {
+    "fft": _names(
+        """
+        dct dctn dst dstn fft fft2 fftfreq fftn fftshift fht get_fft_plan hfft
+        hfft2 hfftn idct idctn idst idstn ifft ifft2 ifftn ifftshift ifht ihfft
+        ihfft2 ihfftn irfft irfft2 irfftn next_fast_len rfft rfft2 rfftfreq
+        rfftn
+        """,
+    ),
+    "fftpack": _names(
+        """
+        fft fft2 fftn get_fft_plan ifft ifft2 ifftn irfft rfft
+        """,
+    ),
+    "interpolate": _names(
+        """
+        Akima1DInterpolator BPoly BSpline BarycentricInterpolator
+        CloughTocher2DInterpolator CubicHermiteSpline CubicSpline
+        InterpolatedUnivariateSpline KroghInterpolator LSQUnivariateSpline
+        LinearNDInterpolator NdBSpline NdPPoly NearestNDInterpolator PPoly
+        PchipInterpolator RBFInterpolator RegularGridInterpolator
+        UnivariateSpline barycentric_interpolate interp1d interpn
+        krogh_interpolate make_interp_spline make_lsq_spline pchip
+        pchip_interpolate splantider splder
+        """,
+    ),
+    "linalg": _names(
+        """
+        bandwidth block_diag circulant companion convolution_matrix dft expm
+        fiedler fiedler_companion hadamard hankel helmert hilbert khatri_rao
+        kron leslie lu lu_factor lu_solve solve_triangular toeplitz
+        """,
+    ),
+    "ndimage": _names(
+        """
+        affine_transform binary_closing binary_dilation binary_erosion
+        binary_fill_holes binary_hit_or_miss binary_opening binary_propagation
+        black_tophat center_of_mass convolve convolve1d correlate correlate1d
+        distance_transform_edt extrema find_objects fourier_ellipsoid
+        fourier_gaussian fourier_shift fourier_uniform gaussian_filter
+        gaussian_filter1d gaussian_gradient_magnitude gaussian_laplace
+        generate_binary_structure generic_filter generic_filter1d
+        generic_gradient_magnitude generic_laplace grey_closing grey_dilation
+        grey_erosion grey_opening histogram iterate_structure label
+        labeled_comprehension laplace map_coordinates maximum maximum_filter
+        maximum_filter1d maximum_position mean median median_filter minimum
+        minimum_filter minimum_filter1d minimum_position morphological_gradient
+        morphological_laplace percentile_filter prewitt rank_filter rotate shift
+        sobel spline_filter spline_filter1d standard_deviation sum sum_labels
+        uniform_filter uniform_filter1d value_indices variance white_tophat zoom
+        """,
+    ),
+    "signal": _names(
+        """
+        BadCoefficients CZT StateSpace TransferFunction ZerosPolesGain ZoomFFT
+        abcd_normalize argrelextrema argrelmax argrelmin band_stop_obj bilinear
+        bilinear_zpk bode buttap butter buttord cheb1ap cheb1ord cheb2ap
+        cheb2ord cheby1 cheby2 check_COLA check_NOLA chirp choose_conv_method
+        coherence cont2discrete convolve convolve2d correlate correlate2d
+        correlation_lags csd cspline1d cspline1d_eval cspline2d cwt czt
+        czt_points dbode decimate deconvolve detrend dfreqresp dimpulse dlsim
+        dlti dstep ellip ellipap ellipord fftconvolve filtfilt find_peaks
+        findfreqs firls firwin firwin2 freqresp freqs freqs_zpk freqz freqz_sos
+        freqz_zpk gammatone gauss_spline gausspulse get_window group_delay
+        hilbert hilbert2 iircomb iirdesign iirfilter iirnotch iirpeak impulse
+        invres invresz istft kaiser_atten kaiser_beta kaiserord lfilter
+        lfilter_zi lfiltic lombscargle lp2bp lp2bp_zpk lp2bs lp2bs_zpk lp2hp
+        lp2hp_zpk lp2lp lp2lp_zpk lsim lti max_len_seq medfilt medfilt2d
+        minimum_phase morlet morlet2 normalize oaconvolve order_filter
+        peak_prominences peak_widths periodogram place_poles qmf qspline1d
+        qspline1d_eval qspline2d resample resample_poly residue residuez ricker
+        savgol_coeffs savgol_filter sawtooth sepfir2d sos2tf sos2zpk sosfilt
+        sosfilt_zi sosfiltfilt sosfreqz spectrogram spline_filter square ss2tf
+        ss2zpk step stft sweep_poly symiirorder1 symiirorder2 tf2sos tf2ss
+        tf2zpk unique_roots unit_impulse upfirdn vectorstrength welch wiener
+        zoom_fft zpk2sos zpk2ss zpk2tf
+        """,
+    ),
+    "sparse": _names(
+        """
+        SparseEfficiencyWarning SparseWarning block_array block_diag bmat
+        coo_array coo_matrix csc_array csc_matrix csr_array csr_matrix dia_array
+        dia_matrix diags diags_array eye eye_array find get_index_dtype hstack
+        identity issparse isspmatrix isspmatrix_coo isspmatrix_csc
+        isspmatrix_csr isspmatrix_dia kron kronsum matrix_transpose permute_dims
+        rand random random_array safely_cast_index_arrays sparray spdiags
+        spmatrix swapaxes tril triu vstack
+        """,
+    ),
+    "sparse.csgraph": _names(
+        """
+        connected_components
+        """,
+    ),
+    "sparse.linalg": _names(
+        """
+        LinearOperator SuperLU aslinearoperator bicgstab cg cgs eigsh factorized
+        gmres lobpcg lsmr lsqr matrix_power minres norm spilu splu spsolve
+        spsolve_triangular svds
+        """,
+    ),
+    "spatial": _names(
+        """
+        Delaunay KDTree distance_matrix
+        """,
+    ),
+    "special": _names(
+        """
+        bdtr bdtrc bdtri beta betainc betaincinv betaln binom boxcox boxcox1p
+        btdtr btdtri cbrt chdtr chdtrc chdtri cosdg cosm1 cotdg digamma
+        ellipeinc ellipj ellipk ellipkinc ellipkm1 entr erf erfc erfcinv erfcx
+        erfinv exp1 exp10 exp2 expi expit expm1 expn exprel fdtr fdtrc fdtri
+        gamma gammainc gammaincc gammainccinv gammaincinv gammaln gammasgn gdtr
+        gdtrc huber i0 i0e i1 i1e inv_boxcox inv_boxcox1p j0 j1 k0 k0e k1 k1e
+        kl_div lambertw log1p log_expit log_ndtr log_softmax loggamma logit
+        logsumexp lpmv multigammaln nbdtr nbdtrc nbdtri ndtr ndtri pdtr pdtrc
+        pdtri poch polygamma pseudo_huber psi radian rel_entr rgamma round
+        shichi sici sinc sindg softmax sph_harm sph_harm_y spherical_yn tandg
+        wright_bessel xlog1py xlogy y0 y1 yn zeta zetac
+        """,
+    ),
+    "stats": _names(
+        """
+        boxcox_llf entropy trim_mean zmap zscore
+        """,
+    ),
+}
 
 
 class _Proxy:
@@ -168,10 +266,10 @@ def call(*args: Any, **kwargs: Any) -> Any:
 
 def _subpackage(path: str, cupy: Any) -> types.ModuleType:
     module = types.ModuleType(f"cupyx.scipy.{path}")
-    allowed = _ALLOWED.get(path)
+    allowed = NAMES[path]
 
     def __getattr__(name: str) -> Any:
-        if name.startswith("__") or (allowed is not None and name not in allowed):
+        if name.startswith("__") or name not in allowed:
             raise AttributeError(
                 f"module 'cupyx.scipy.{path}' has no attribute {name!r}"
             )