Cython implementation multiply and from_attributes

[oarcelus]

Following the discussion in jonathf/chaospy#433 here the PR for the Cythong implementation of multiply and from_attributes, allowing to bypass inefficient key lookups to numpy structured arrays.

[jonathf]

Okay, now I have a moment to spare.

I am familiar with pybind11 and numpy integration, but I have no experience with Cython specifically. Perhaps you can clarify the approach here.

Is the recommended idea here that we pre-compile binaries with something like cibuildwheel, or de we ship the code as-is, and let the user "Cythonize" their own code? If the latter, is Cython self contained on all major platforms, or do the user need something like gcc or clang to make this work?

Second thing: I don't seem to get the code working on my machine. I did a clean install on my M3 Mac. Seems like it generated darwin shared objects, but the code fails:

>>> import numpoly
>>> x, y = numpoly.variable(2)
>>> x * y
Fatal Python error: Segmentation fault

Thread 0x000000016cc47000 (most recent call first):
  File "/Users/jonathf/.pyenv/versions/3.12.7/lib/python3.12/threading.py", line 355 in wait
  File "/Users/jonathf/.pyenv/versions/3.12.7/lib/python3.12/threading.py", line 655 in wait
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/core/history.py", line 903 in run
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/core/history.py", line 61 in only_when_enabled
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/decorator.py", line 232 in fun
  File "/Users/jonathf/.pyenv/versions/3.12.7/lib/python3.12/threading.py", line 1075 in _bootstrap_inner
  File "/Users/jonathf/.pyenv/versions/3.12.7/lib/python3.12/threading.py", line 1032 in _bootstrap

Current thread 0x00000001ed224f40 (most recent call first):
  File "/Users/jonathf/personal/numpoly/numpoly/array_function/multiply.py", line 112 in multiply
  File "/Users/jonathf/personal/numpoly/numpoly/baseclass.py", line 242 in __array_ufunc__
  File "<ipython-input-4-e32109319f52>", line 1 in <module>
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3577 in run_code
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3517 in run_ast_nodes
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3334 in run_cell_async
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/core/async_helpers.py", line 128 in _pseudo_sync_runner
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3130 in _run_cell
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/core/interactiveshell.py", line 3075 in run_cell
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/terminal/interactiveshell.py", line 910 in interact
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/terminal/interactiveshell.py", line 917 in mainloop
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/terminal/ipapp.py", line 317 in start
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/traitlets/config/application.py", line 1075 in launch_instance
  File "/Users/jonathf/.venv/numpoly/312/lib/python3.12/site-packages/IPython/__init__.py", line 130 in start_ipython
  File "/Users/jonathf/.venv/numpoly/312/bin/ipython", line 8 in <module>

Extension modules: numpy._core._multiarray_umath, numpy.linalg._umath_linalg, numpoly.cfunctions.cmultiply, numpoly.cfunctions.cfrom_attributes, numpoly.cfunctions.cvalues, numpoly.cfunctions.cdispach (total: 6)
fish: Job 1, 'ipython --nosep' terminated by signal SIGSEGV (Address boundary error)

Is this expected? Any pitfalls that I should be aware of that would causes this segfault?

[oarcelus]

Hi, it should work now. I forgot to change some things in multiply.py. I tested the generate expansion function for different hyperbolic truncation values and orders. But maybe a deeper test is require to explore edge cases?

The approach of the cython implemented functions is to bypass the key lookup for numpy structured array. For example in cmultiply.pyx you can see calls to functions cadd_values and cset_values. These functions check for the input array dimensions and type first and branch to a cdef function for the particular data. I implemented functions for input array dimensions 0, 1, and 2, and dtypes np.int64 and np.float64, because these are the cases that I have seen in my tests. But I do not know if the array of coefficients can take more dimensions or not. In this case the error came from the fact that the program branched into a function that request the strides of a ndarray that has ndim==0.

About the possibility of distributing the code. I don't really know. Cython is not completely self contained so I would say that the safest option is to use cibuildwheels.

[jonathf]

That works. I've fixed some numpy v2 bugs. I've tested master branch, and everything works except save/load. That shouldn't be relevant for this work.

Locally running pytest --doctest-modules numpoly/ test/ on this branch reveals that there are a few things that are broken.

I am trying to get CI to work, but seems like there is an issue there. I will take a second look at that soon.

[oarcelus]

I am still getting additional errors apart from save/load errors, when reverting to the non-cython implementations.

Things related to the doctest expected outputs. I guess it is not really important?

____________________________________________ [doctest] numpoly.array_function.any.any ____________________________________________
043             result as dimensions with size one. With this option, the result
044             will broadcast correctly against the input array.
045
046     Return:
047         A new boolean or `ndarray` is returned unless `out` is specified,
048         in which case a reference to `out` is returned.
049
050     Example:
051         >>> q0 = numpoly.variable()
052         >>> numpoly.any(q0)
Expected:
    True
Got:
    np.True_

[jonathf]

On master or on commit b11f7a4?

Yeah, don't worry about cosmetic errors. Get the real problems sorted, and I can iron out the rest.

[oarcelus]

I have made all the changes neccessary. All tests pass correctly, except the load/save ones and the ones related to format changes due to changes in Numpy 2.

Regarding the implementation, it is not so optimal because I am duplicating code for each supported type (bool, uint32, int64, float64, and complex128). I am pretty sure that I can improve upon it. Anyways tests pass, it is still as fast, and my own tests are also correct.

[jonathf]

I've made a beta release. Install it with pip install numpoly==1.3.0b1. I am going to test it out upstream a bit. If it works as expected, will merge this.

For now I deprecated save and load. Pickle still works, and same goes for savetxt and loadtxt.