MAINT: interpolate: move fpknot wrap from Cython to manual #14
Conversation
| py_fpknot(PyObject* self, PyObject *args) | ||
| { | ||
| PyObject *py_x=NULL, *py_t=NULL, *py_residuals=NULL; | ||
| int k; // XXX: python int => "i" in ParseTuple? |
There was a problem hiding this comment.
This is fine as long as you never expect reasonable values of k to vall outside the values of a C int. You also need to check the return value of PyArg_ParseTuple for the error case to catch the OverflowError, but you do that.
| return NULL; | ||
| } | ||
|
|
||
| // XXX: Do not need to DECREF py_x etc from PyArg_ParseTuple? |
There was a problem hiding this comment.
Since it's a borrowed reference, you don't if you never create an owned reference.
That said, borrowed references are really hard to reason about, and you're creating objects based on py_x and py_t below, so I'd personally INCREF here and then DECREF at the end when you're done with them.
There was a problem hiding this comment.
Thanks! Point taken about borrowed vs owned references. Could use an advise though: If the whole trampoline function is literally
PyArg_ParseTupleto get python argumentscheck_array- grab pointers from arrays and delegate to numpy-oblivious C code
(this py_fpknot is what I've in mind) keeping all references borrowed looks simpler. How does it fare in the nogil context / under free-threading though, and does switching to owned references change the equation? From reading the free-threading docs (which you wrote?), I've an impression that owned references are generally safer but not always. Don't have a good gut feeling of what is safe and what is not.
|
|
||
|
|
||
| /* | ||
| // XXX: overkill? Mimic cython's `double[::1] x` etc --- replaced by check_array stanza below |
There was a problem hiding this comment.
It doesn't make a ton of difference if you don't want the API to be flexible and cast from other dtypes if you don't hand it a double array or don't hand it a C-contiguous array.
If you already control all that then what you have with check_array is fine.
| PyArrayObject *a_t = (PyArrayObject *)py_t; | ||
| PyArrayObject *a_residuals = (PyArrayObject *)py_residuals; | ||
|
|
||
| // XXX: npy_intp vs Py_ssize_t vs ssize_t (== ptrdiff_t)? |
There was a problem hiding this comment.
In NumPy 2.0, npy_intp and Py_ssize_t are the same on all platforms. Before NumPy 2.0, npy_intp was 32 bit on windows. It's not quite the same as ptrdiff_t because it's a sized int (to make error checking less awful - negative values always indicate an error).
There was a problem hiding this comment.
PY_SSIZE_T_MAX (which is the same as sys.maxsize in python) is also often used as a sigil value for things that are typed as Py_ssize_t.
There was a problem hiding this comment.
Ugh. I full well realize it's asking for more and more, so if you've bandwidth to spare for me, I could use an advice.
For context:
- The goal is to communicate between python and the C++ extension which does actual computations (https://github.com/ev-br/scipy/blob/bspl_no_fused/scipy/interpolate/src/__fitpack.cc);
- The extension does not use anything from python or numpy; Also does not do any memory allocations and relies on the caller to send in everything including temp buffers.
- The extension uses
ssize_tfor all indexing where the classic C++ would usesize_t, to explicitly avoid any signed/unsigned mess. - Under MSVC on Windows
ssize_tis not defined (apparently), so I just#define ssize_t ptrdiff_t - Have to support Windows, and int32 overflows are not impossible
- I control both C++ and python side.
So when I need to return to python an array of integer indices (offsets: for a 2D array A, nonzeros in row j start at A[ j, offset[j]]), my thinking was
- allocate
np.empty(..., dtype=np.intp)on the python side - use
npy_intpin the python/C glue - keep using
ssize_tin C++
Does this sound reasonable or am I making things hard for myself?
There was a problem hiding this comment.
Ah, it turns out that Py_ssize_t is a typedef for ssize_t where available, so I think using Py_ssize_t on the python side and then ssize_t on the C++ side should be sufficient.
There was a problem hiding this comment.
Which leaves a question of the matching numpy dtype/typecode still?
There was a problem hiding this comment.
So the resolution (as discussed offline) is:
- under numpy 2.x,
npy_intpis interchangeable withPy_ssize_t(== ssize_t == ptrdiff_t). - under numpy 1.x and on Windows,
npy_intpis 32-bit, but indexing in numpy itself is problematic if indices overflow 32-bit sized ints.
So in the end, the best one can do here is to detect numpy 1.x and Windows, check that the offsets are too large and bail out if so.
| static_cast<const double *>(PyArray_DATA(a_residuals)) | ||
| ); | ||
|
|
||
| // XXX: need to DECREF a_* variables? Apparently not unless PyArray_ContiguousFromObject |
There was a problem hiding this comment.
ContiguousFromObject returns a new reference to an array object. That array might just be a view onto the original array, but it's a new PyObject and so needs to be DECREF'd. If you don't go through that function you're creating a borrowed reference, which doesn't need to be DECREF'd. That said, borrowed references are hard to reason about (as you saw!) since most APIs return new references.
| Py_ssize_t nc; | ||
| Py_ssize_t startrow=1; // XXX: optional, keeps the value intact if not given? | ||
|
|
||
| const char *kwlist[] = {"a", "offset", "nc", "y", "startrow", NULL}; |
There was a problem hiding this comment.
just an FYI, depending on what you're doing, if this function is called a lot there's a surprising amount of overhead to create python strings from these strings each time the function is called. For internal-only functions implemented in C, I would look at only using arguments positionally. NumPy goes out of its way to avoid using ParseTupleAndKeywords in hot code paths and also has a whole keyword argument name caching system.
| } | ||
|
|
||
| if (!(check_array(py_a, 2, NPY_DOUBLE) && | ||
| check_array(py_offs, 1, NPY_LONG) && // XXX: typecode for ssize_t (==ptrdiff_t)? |
There was a problem hiding this comment.
If you need to support windows and NumPy 1.x, this might be annoying to deal with if you need specifically a Py_ssize_t, since it may or may not be the same depending on architecture. You could always do a conversion, which will be a no-op if it really is the type you want. Also you probably want NPY_DEFAULT_INT instead of NPY_LONG, which is wrong on NumPy 2.x.
| ); | ||
|
|
||
| //Py_DECREF(a_offs); | ||
| // XXX: a & y modified in-place: need to incref? |
There was a problem hiding this comment.
Nope, not with a borrowed reference. You can always temporarily take ownership for the length of the function to make it easier to reason about.
| /* | ||
| XXX: | ||
| 1. how to check for refleaks / refcounting errors? | ||
| 2. goto fail & crosses initialization errors: a better pattern or copy-paste? |
There was a problem hiding this comment.
You either declare all your variables at the top of the function, or you create a new scope with {} and define the temporary variables in the new scope.
|
|
||
| /* | ||
| XXX: | ||
| 1. how to check for refleaks / refcounting errors? |
There was a problem hiding this comment.
Think really hard. There's also https://pypi.org/project/pytest-leaks/ but it's bitrotting pretty hard at this point. You can make a quick and dirty version with a debug python build using sys.gettotalrefcount() though.
ev-br
force-pushed
the
bspl_c_api
branch
2 times, most recently
from
0568a16 to
4148566
Compare
ev-br
force-pushed
the
bspl_no_fused
branch
3 times, most recently
from
fdeb0a3 to
aec8a4f
Compare
Removes an extra `xp.asarray` inside `xp.concat` which was causing errors due to setting an array with a sequence.
- Use `scipy._lib.array_api_compat.numpy` rather than `numpy` as a default when no array specified - Remove use of in-place item assignment in `_cartesian_product` - Remove use of `xp.newaxis` since Torch doesn't support this - Switch to using positional arguments for `xp.linspace` since Torch uses different keywords - Modify slow tests
- Use same style of typing as `_array_api` for `Array` - Use `np_compat` from `_lib.array_api` rather than `_lib.array_api_compat.numpy` - Default to `np.linalg` when `xp.linalg` doesn't exist - Add `xp` argument to test functions rather than calling `array_namespace` every time - Skip testing `array_like` inputs when the backend is not numpy - Change how `FixedRule`, `NestedRule` and `NestedFixedRule` determine the array namespace
- use a `_concat` helper function to pick whether to use `xp.concat` or `xp.concatenate` - use `xp_test = array_namespace(a)` for using `atan` since `atan` isn't supported in np < 2. - compute subregion coordinates and cartesian products without using `np`-specific code
The failing test is testing np.matrix inputs to RGI. This is allowed for backwards compatibility but does not work with xp_assert_* infrastructure. Thus, convert to arrays for assertions. An alternative is to skip the test when running with SCIPY_ARRAY_API env variable (since then we're not bound by backwards compat). But converting to ndarray seems simpler.
The computational complexity is O(1), so there is no reason for this helper to stay in Cython or move to C++.
[skip cirrus] [skip circle]
- Removes an unnecessary `isinstance(a, list)` check since by default `array_namespace(a_list)` is `np_compat` - Deletes the function `_concat` which was used to fall back to `np.concatenate` if `np.concat` wasn't available, and instead make sure `np_compat` is being used
- Determine the compatible version of the array namespace in each rule
* BUG: special.logsumexp: fix type promotion [skip cirrus] [skip circle] * MAINT: adjust dtype determination logic in xp_broadcast_promote * TST: special.logsumexp: add dtype tests * TST: special.logsumexp: fix tests * MAINT: define mixed-type promotion rule if underlying library doesn't * TST: special.logsumexp: update tests * STY: special.logsumexp: use is_array_api_strict where appropriate [skip ci] --------- Co-authored-by: Matt Haberland <[email protected]>
MAINT: interpolate: move an internal utility from cython to python
…ic_filter`, etc.) (scipy#20811) * add axes support to gradient_magnitude and laplace functions update rank filter axes docstring to match more detailed style of other functions * axes should not be cast using xp.asarray in test case * np.empty_like -> xp.empty_like * use concatenation instead of in-place array modification in test suite for JAX compatibility
* ENH: interpolate: add Floater Hormann weights Reviewed at scipy#21497 Co-authored-by: Pauli Virtanen <[email protected]>
…21784) * add tests of basic matmul involving 1D arrays * fix the bug in A@v
…cipy#21792) * BUG: sparse: fix setdiag for matrices with missing diagonal entries We need to first set existing diagonal entries before adding the missing ones, as inserting new values invalidates the computed offsets. * add tests for fix of csr.setdiag --------- Co-authored-by: Dan Schult <[email protected]>
From python, it's find_interval, not py_find_interval
Otherwise, 32-bit build on linux in CI fails on ./tools/check_pyext_symbol_hiding.sh
…trix Reuse _dierckx.data_matrix, which constructs the same object in a "packed" format, with per-row offsets of non-zero elements instead of per-element CSR indices/indptr arrays. Add `extrapolate=False` parameter to data_matrix to account for design_marix usages.
Co-authored-by: Christian Lorentzen <[email protected]>
use vectors not unique_ptr(malloc) for work arrays
Reference issue
What does this implement/fix?
Additional information