PyO3 / rust-numpy

PyO3-based Rust bindings of the NumPy C-API
BSD 2-Clause "Simplified" License
1.07k stars 98 forks source link

Converting a numpy array into a nagebra matrix returns "None" #413

Closed Ruhrpottpatriot closed 4 months ago

Ruhrpottpatriot commented 4 months ago

I have a library that uses nalgebra internally to do some heavy math stuff. I want to offer the users a python frontend, which does some less computationally intensive and less memory critical stuff.

I use PyO3 to create the bindings. In said bindings the function accepts a numpy array, which is transformed into a nalgebra matrix, which is then used in further computations. However, the conversion between numpy and nalgebra via try_as_matrix always returns None. Here is a MWE:

A = np.zeros((2, 2))
filter = test_lib.foo(A)
#[pyfunction]
pub fn foo(py: Python, Pkk: Py<PyArray<f64, Ix2>> ) -> Option<()> {
    let A =  dbg!(Pkk.as_ref(py));
    let pkk: MatrixView<f64, U2, U2> = dbg!(A.try_as_matrix())?;

    Some(())
}
#[pymodule]
pub(crate) fn test_lib(_py: Python, m: &PyModule) -> PyResult<()> {   
    m.add_function(wrap_pyfunction!(foo, m)?)?;
}

As you can see I'm just passing a 2x2 zero matrix to the function. Inside the function, the first dbg statement correctly prints the contents of the matrix. However, the second dbg statement always returns None.

I have activated the "extension-module" for PyO3 and the "nalgebra" feature for this crate.

adamreichold commented 4 months ago

I think the problem here is that zeros will create a C/row-major memory layout but nalgebra uses Fortan/column-major memory layout. So in your example using

np.zeros((2, 2), order='F')

should make it work.

If you do not control the source array, then calling reshape_with_order should allow you to enforce the correct memory layout at the cost of data copy, e.g.

let array = array.reshape_with_order([2, 2], NPY_ORDER::NPY_FORTRANORDER)?;

(Note that try_as_matrix produces views into the NumPy arrays instead of copying, hence the requirement for them to use the memory layout expected by nalgebra.)

adamreichold commented 4 months ago

As a sort of converse, you should also be able to produce views into C order arrays by using the right strides on the nalgebra side by also passing in stride parameters in addition to just specifying the dimensions.

I did not try it, but something like

let pkk: MatrixView<f64, U2, U2, U2, U1> = A.try_as_matrix()?;

might do the trick.

Ruhrpottpatriot commented 4 months ago

I didn't try the reshape_with_order, but the second approach using C order array by manually setting the stride also works.

IMHO, this information should be added to the documentation for the try_as_matrix function. Becuase I can imagine more than one user stumbling about this.

/edit: However, this doesn't work if I pass a view into the function itself, i.e. assume the same rust code but change the python code to:

A = np.zeros((2, 2, 2))
filter = test_lib.foo(A[:, :, 0])
adamreichold commented 4 months ago

However, this doesn't work if I pass a view into the function itself, i.e. assume the same rust code but change the python code to:

Because the strides are different for the view, so either you do know them at compile time and adjust the type accordingly or you opt into nalgebra's support for dynamic dimensions.

adamreichold commented 4 months ago

Or alternatively, you go via reshape_with_order and force the memory layout assumed by nalgebra by copying if necessary.

Ruhrpottpatriot commented 4 months ago

Or alternatively, you go via reshape_with_order and force the memory layout assumed by nalgebra by copying if necessary. That doesn't work for me. Using reshape_with_order in the following code snippet still returns None:

let array = array
    .as_ref(py)
    .reshape_with_order([2, 2], NPY_ORDER::NPY_FORTRANORDER)
    .unwrap();
let array: MatrixView<f64, U2, U2> = array.try_as_matrix()?;

let I = SMatrix::<f64, 2, 2>::identity();
let r = (array + I).to_pyarray(py).to_owned();
Some(r)
adamreichold commented 4 months ago

That doesn't work for me. Using reshape_with_order in the following code snippet still returns None:

I see, NumPy is probably to clever and elides the copy while using column-major but non-standard strides. I guess using

MatrixView<f64, U2, U2, Dyn, Dyn>

as the view type is the better choice then.

Ruhrpottpatriot commented 4 months ago

Using Dyn in the stride dimensions works perfectly.