ExpertSolver.interpolate: index views as Py_ssize_t, not long

CHANGELOG.md

@@ -2,7 +2,11 @@
 
 ## v1.1.1 (in progress)
 
-*No user-visible changes yet.*
+### Fixed
+
+- **`ExpertSolver.interpolate(mode='nearest')` works on Windows.** The model indices it returns in
+  `I_out`, and accepts back as `I`, are now `np.intp`, the type the nearest-neighbour search produces,
+  where they were C `long`. The two differ in width on Windows.
 
 
 ## v1.1.0 (10 September 2026)

TODO_DEFERRED.md

@@ -50,20 +50,6 @@ Potential arXiv tutorial target — rewrite them with the "surrogate model,
 not Taylor series" framing, and mention that the same method appears in the
 literature under several names (MLS, WLSQM, "diffuse approximation").
 
-## Windows `long` width in ExpertSolver.interpolate
-
-*Cluster: portability · Cost: S · Gate: none · Filed: 2026-04-14*
-
-`ExpertSolver.interpolate()` in `wlsqm/fitter/expert.pyx` backs its
-`I_out` return array with a Cython `long[::1]` view and an allocation of
-`dtype=np.int_`. On Linux/macOS 64-bit, C `long` is 64 bits; on Windows
-64-bit (MSVC) it is 32 bits. The function therefore silently produces
-different element widths across platforms, and cannot address more than
-2**31 local models on Windows. Low priority (no one has a WLSQM setup
-with > 2 billion neighborhoods) but the fix is small: switch both the
-cdef view type and the numpy dtype to a fixed width, e.g. `np.int64_t`
-/ `np.int64`.
-
 ## Weighting function support
 
 *Cluster: features · Cost: M · Gate: none · Filed: 2026-04-14*

tests/test_expert.py

@@ -162,3 +162,34 @@ def test_expert_3d_single_case(rng):
     fi = np.zeros((1, wlsqm.number_of_dofs(3, 2)))
     es.solve(fk=fk, fi=fi)
     assert np.allclose(fi[0], fi_expected, atol=1e-10)
+
+
+def test_expert_interpolate_nearest_returns_and_accepts_model_indices(rng):
+    """`interpolate(mode='nearest')` picks the nearest local model for each point, reports which one in
+    `I_out`, and accepts that array back as `I` to skip the search on a re-interpolation.
+
+    The index arrays are what `scipy.spatial.cKDTree.query` returns, `np.intp`. On Windows that is wider
+    than C `long`, so an index view declared `long` refuses the very array the search hands it.
+    """
+    f, _fi_expected = poly2d_order2()
+    npts = 30
+    # Two local models, one either side of the y axis, each fitted to the same exact polynomial.
+    xi_arr = np.array([[-0.5, 0.0], [0.5, 0.0]])
+    xk_arr = np.stack([xi + rng.uniform(-0.3, 0.3, size=(npts, 2)) for xi in xi_arr])
+    fk_arr = np.stack([f(xk) for xk in xk_arr])
+    fi_arr = np.zeros((2, wlsqm.number_of_dofs(2, 2)))
+
+    es = _make_expert_2d(ncases=2, nk_per_case=npts)
+    es.prepare(xi=xi_arr, xk=xk_arr)
+    es.solve(fk=fk_arr, fi=fi_arr)
+    es.prep_interpolate()
+
+    x = np.array([[-0.6, 0.1], [-0.4, -0.1], [0.4, 0.05], [0.6, -0.05]])
+    out, I_out = es.interpolate(x, mode="nearest", diff=wlsqm.i2_F)
+    assert I_out.dtype == np.intp
+    assert list(I_out) == [0, 0, 1, 1]
+    assert np.allclose(out, f(x), atol=1e-10)
+
+    out_again, I_again = es.interpolate(x, mode="nearest", diff=wlsqm.i2_F, I=I_out)
+    assert np.array_equal(out_again, out)
+    assert np.array_equal(I_again, I_out)

wlsqm/fitter/expert.pyx

@@ -726,7 +726,7 @@ I    : If mode='nearest': override which local model to use for each point in x,
 
 Return value: tuple (out, I_out) where
         out   = function value (or derivative value, depending on "diff")
-        I_out = if mode='nearest', index of local model used for each point in x (array of shape (nx,), dtype np.int_, i.e. C long).
+        I_out = if mode='nearest', index of local model used for each point in x (array of shape (nx,), dtype np.intp).
                 This can be passed back in as "I".
 
                 if mode='continuous', this is always None (i.e. currently not supported).
@@ -756,16 +756,14 @@ Return value: tuple (out, I_out) where
         #
         cdef int nx = x.shape[0]
         cdef double[::1] out = np.empty( (nx,), dtype=np.float64 )
-        cdef long[::1] I_out
+        cdef Py_ssize_t[::1] I_out
         if mode == 'nearest':
             if I is not None:
                 I_out = None  # if 'I' was given, don't bother copying it to I_out in expert_interpolate_nearest()
             else:
-                # np.int_ matches the Cython `long[::1]` view declared above,
-                # and — unlike the long-removed `np.long` — works on both
-                # NumPy 1.x and 2.x. (On Windows, C `long` is 32-bit; this is
-                # a pre-existing portability caveat noted in TODO_DEFERRED.md.)
-                I_out = np.empty( (nx,), dtype=np.int_ )  # I_out[j] = index of the local model (in self.cases) used to produce out[j]
+                # np.intp is what cKDTree.query returns, and matches the `Py_ssize_t` views on every
+                # platform. C `long` does not: it is 32-bit on Windows, where np.intp is 64-bit.
+                I_out = np.empty( (nx,), dtype=np.intp )  # I_out[j] = index of the local model (in self.cases) used to produce out[j]
 
             expert_interpolate_nearest( dimension, self.tree, manager, x, out, I_out, I, diff, ntasks )
 
@@ -827,12 +825,12 @@ cdef int expert_solve_one_iterative( infra.Case* case, double* fi, double[::view
     return impl.solve_iterative( case, fk, sens, do_sens, taskid, max_iter, xkManyD, xk1D )
 
 
-cdef void expert_interpolate_nearest( int dimension, xi_tree, infra.CaseManager* manager, x, double[::1] out, long[::1] I_out, long[::1] I_in, int diff, int ntasks ):
+cdef void expert_interpolate_nearest( int dimension, xi_tree, infra.CaseManager* manager, x, double[::1] out, Py_ssize_t[::1] I_out, Py_ssize_t[::1] I_in, int diff, int ntasks ):
     # For each point in x, find the local model whose origin is nearest (the nearest point in xi).
     #
     # This search takes the majority of the runtime of this function.
     #
-    cdef long[::1] I  # scipy.spatial.cKDTree.query() returns an array of long
+    cdef Py_ssize_t[::1] I  # scipy.spatial.cKDTree.query() returns an array of np.intp
     if I_in is None:  # usual case
         _distances, I = xi_tree.query( x, k=1 )  # TODO: _distances could be a useful quality metric
     else:  # use the caller-specified model indices (useful when re-interpolating an updated model for the same points x)