From dc7e794cf2813c9abd67b76535b5008a10f92a9f Mon Sep 17 00:00:00 2001 From: Juha Jeronen Date: Wed, 7 Oct 2026 16:25:12 +0300 Subject: [PATCH 1/2] Test ExpertSolver.interpolate in nearest mode Nothing exercised it: test_interp covers `interpolate_fit`. The model indices it returns and accepts are what cKDTree.query returns, `np.intp`, which on Windows is wider than the C `long` the views are declared as. Expected to fail on the Windows runner until the views are fixed. Co-Authored-By: Claude Opus 5.5 --- tests/test_expert.py | 31 +++++++++++++++++++++++++++++++ 1 file changed, 31 insertions(+) diff --git a/tests/test_expert.py b/tests/test_expert.py index f1193b6..b3e7423 100644 --- a/tests/test_expert.py +++ b/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) From a7c67646068066654a4b2453353163fd9f14f875 Mon Sep 17 00:00:00 2001 From: Juha Jeronen Date: Wed, 7 Oct 2026 16:27:56 +0300 Subject: [PATCH 2/2] ExpertSolver.interpolate: index views as Py_ssize_t, not long MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The model-index views were C `long`, while cKDTree.query returns `np.intp` and `I_out` was allocated as `np.int_`. On Windows `long` is 32-bit and both of those are 64-bit, so `mode='nearest'` raised "Buffer dtype mismatch, expected 'long' but got 'long long'" there — seen on the Windows runner with the previous commit's test. `Py_ssize_t` views and an `np.intp` allocation match on every platform. Co-Authored-By: Claude Opus 5.5 --- CHANGELOG.md | 6 +++++- TODO_DEFERRED.md | 14 -------------- wlsqm/fitter/expert.pyx | 16 +++++++--------- 3 files changed, 12 insertions(+), 24 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index e24e49e..227a553 100644 --- a/CHANGELOG.md +++ b/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) diff --git a/TODO_DEFERRED.md b/TODO_DEFERRED.md index 5bdeeef..1a0df41 100644 --- a/TODO_DEFERRED.md +++ b/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* diff --git a/wlsqm/fitter/expert.pyx b/wlsqm/fitter/expert.pyx index 251f2d2..d2d5d4b 100644 --- a/wlsqm/fitter/expert.pyx +++ b/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)