Assigning numpy data in cython to a view

Question

I am trying to assign the output of linalg inverse function (la.inv) to a view in cython. Unfortunately this does not work. I can always assign the output of la.inv() to a temporary ndarray object and then copy its content to the view.

Is there a better way to do it.

cpdef int testfunc1(np.ndarray[np.float_t, ndim=2] A,
                    double [:,:] B) except -1:

    print("inverse of A:", la.inv(A))
    if np.isnan(A).any():
        return -1
    else:
        B = la.inv(A)
        return 1


cpdef int testfunc2(np.ndarray[np.float_t, ndim=2] A) except -1:
    cdef long p = np.shape(A)[0], status
    cdef B = np.zeros(shape=(p, p), dtype=float)
    cdef double[:,:] BView = B
    print("before inverse. B: ", B)
    status = testfunc1(A, BView)
    print("after inverse. B: ", B)
    if status == -1:
        return -1
    else:
        return 1

The output:

A = np.random.ranf(4).reshape(2, 2)
        status = testfunc2(A)
        if status == -1:
            raise ValueError("nan cell.")
        else:
            print("pass")

('before inverse. B: ', array([[ 0.,  0.],
       [ 0.,  0.]]))
('inverse of A:', array([[ 4.4407987 , -0.10307341],
       [-2.26088593,  1.19604499]]))
('after inverse. B: ', array([[ 0.,  0.],
       [ 0.,  0.]]))

Saullo G. P. Castro · Answer 1 · 2014-05-23T12:15:25.953

You can create a temporary buffer that will receive the value of the la.inv() and then populate the memory view:

import numpy as np
cimport numpy as np
import numpy.linalg as la

cpdef int testfunc1(np.ndarray[np.float_t, ndim=2] A,
                    double [:,:] B) except -1:
    cdef np.ndarray[np.float_t, ndim=2] buff
    cdef int i, j

    print("inverse of A:", la.inv(A))
    if np.isnan(A).any():
        return -1
    else:
        buff = la.inv(A)
        for i in range(buff.shape[0]):
            for j in range(buff.shape[1]):
                B[i, j] = buff[i, j]
        return 1

cpdef int testfunc2(np.ndarray[np.float_t, ndim=2] A) except -1:
    cdef long p = np.shape(A)[0], status
    cdef B = np.zeros(shape=(p, p), dtype=float)
    cdef double[:,:] BView = B
    print("before inverse. B: ", B)
    status = testfunc1(A, BView)
    print("after inverse. B: ", B)
    if status == -1:
        return -1
    else:
        return 1

As pointed out by @MrE, you can use np.copyto() if you use a np.ndarray instead of a MemoryView:

cpdef int testfunc1(np.ndarray[np.float_t, ndim=2] A,
                    np.ndarray[np.float_t, ndim=2] B) except -1:
    cdef int i, j
    print("inverse of A:", la.inv(A))
    if np.isnan(A).any():
        return -1
    else:
        np.copyto(B, la.inv(A))
        return 1

cpdef int testfunc2(np.ndarray[np.float_t, ndim=2] A) except -1:
    cdef long p = np.shape(A)[0], status
    cdef np.ndarray[np.float_t, ndim=2] B, BView
    B = np.zeros(shape=(p, p), dtype=float)
    BView = B
    print("before inverse. B: ", B)
    status = testfunc1(A, BView)
    print("after inverse. B: ", B)
    if status == -1:
        return -1
    else:
        return 1

YXD · Answer 2 · 2014-06-04T17:24:00.607

1

This isn't caused by views or by Cython. B = la.inv(A) creates a new array and gives it the name B in the scope of testfunc1. This does not affect the array with the name B in testfunc2.

Be aware that your code, where the heavy lifting is done by NumPy functions, is unlikely to benefit from Cython.

One way to make this work is to do:

np.copyto(B, la.inv(A))

in testfunc1. @SaulloCastro mentions that this does't work in Cython as B has a memory view type, however you may be able make it work by declaring the argument B as an ndarray (not sure about this). Otherwise without Cython:

>>> import numpy as np
>>> X = np.zeros((5, 5))
>>> B = X[:3, :3]
>>> A = np.ones((3, 3))
>>> np.copyto(B, A)
>>> X
array([[ 1.,  1.,  1.,  0.,  0.],
       [ 1.,  1.,  1.,  0.,  0.],
       [ 1.,  1.,  1.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.],
       [ 0.,  0.,  0.,  0.,  0.]])
>>>

edited Jun 04 '14 at 17:24

answered May 23 '14 at 11:54

YXD

31,741
15
75
115

+1 I liked you approach, but how to make `np.copyto()` work with a memoryvew? Here it returned: `TypeError: argument 1 must be numpy.ndarray, not _stack._memoryviewslice` – Saullo G. P. Castro May 23 '14 at 12:01
Thanks for both the approaches. I had implemented it using a temporary buffer and copying it back to B, but I want to save up on the time required to do the copying. I was however hoping to just redirect B to point to the "contents/data pointer" of what la.inv() returns. I dont mind defining B as a pointer to a double array in that case, but I would want to know how to access just the "data pointer" of la.inv(). Any idea how? – user1971988 May 23 '14 at 12:09
@SaulloCastro updated the answer a little in terms of the description, but I don't know if it'll work with a memoryview. Thanks for the "pointer", will delete the question if it turns out to be useless – YXD May 23 '14 at 12:11

score 1 · Answer 3 · edited May 23 '17 at 12:14

If I create a mememoryview on la.inv(A) I can perform a 1 step, and presumably efficient, memoryview to memoryview copy:

cpdef int testfunc1c(np.ndarray[np.float_t, ndim=2] A,
                    double [:,:] BView) except -1:
    cdef double[:,:] CView
    print("inverse of A:", la.inv(A))
    if np.isnan(A).any():
        return -1
    else:
        CView = la.inv(A)
        BView[...] = CView
        return 1

producing:

In [4]: so23827902.testfunc2(A)
('before inverse. B: ', array([[ 0.,  0.],
       [ 0.,  0.]]))
('inverse of A:', array([[ 1.04082818, -0.14530117],
       [-0.24050511,  1.13292585]]))
('after inverse. B: ', array([[ 1.04082818, -0.14530117],
       [-0.24050511,  1.13292585]]))
Out[4]: 1

I'm guessing that the the memoryview copy will be faster, but the sample arrays are too small for meaningful time tests.

I tested this as part of responding at https://stackoverflow.com/a/30418448/901925

In Python you can reassign the data buffer of an array (albeit at some risk):

B = np.zeros_like(A)
C = la.inv(A)
B.data = C.data

cython raises errors about unsafe pointers during the compile stage with this statement.

Inspired by examples I found for https://stackoverflow.com/a/28855962/901925 using np.PyArray_SimpleNewFromData, I tried using other PyArray... functions to do the same sort of base reassignment:

np.PyArray_SetBaseObject(B, np.PyArray_BASE(la.inv(A)))

Currently I am trying to resolve a AttributeError: 'module' object has no attribute 'PyArray_SetBaseObject' error.

Assigning numpy data in cython to a view

3 Answers3

Linked