[SciPy-Dev] Generate random variates using Cython

[Matti Picus]

On 16/1/20 6:14 am, Christoph Baumgarten wrote:
>
> Hi,
>
> I recently looked into rejection algorithms to generate random 
> variates that rely on repeatedly generating uniform random variates 
> until a termination condition is satisfied.
> A pure Python implementation is reasonably fast but much slower than a 
> C implementation (e.g. Poisson distribution in numpy). I tried to use 
> Cython to adapt my Python code (I just started to look into Cython, so 
> it could be that I overlook something very basic) but I did not 
> succeed and I would like to ask for help.
>
> The following code works that generates rvs on [0, 1] with density 3 
> x**2 as an example (compiling it generates a warning though: Warning 
> Msg: Using deprecated NumPy API, disable it with #define 
> NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION)  but I struggle to get rid 
> of the Python interaction np.random.rand() to speed it up.
>
> cimport numpy as np
> from numpy cimport ndarray, float64_t
> import numpy as np
> def rvs(int N):
>     cdef double u, v
>     cdef int i
>     cdef np.ndarray[np.float64_t, ndim=1] x = np.empty((N, ), 
> dtype=np.float64)
>     for i in range(N):
>         while (1):
>             u = np.random.rand()
>             v = np.random.rand()
>             if u <= v*v:
>                 break
>         x[i] = v
>     return x
>
> Following https://numpy.org/devdocs/reference/random/c-api.html 
> (Cython API for random), I tried
> cimport numpy.random to use random_standard_uniform.
> However, I get an error: 'numpy\random.pxd' not found
>
> I get similar errors when trying to compile the Cython examles here:
> https://docs.scipy.org/doc/numpy-1.17.0/reference/random/extending.html
>
> Cython version: 0.29.14
> Numpy 1.17.4
> python 3.7.1
>
> Any guidance on how to access the uniform distribution in numpy using 
> Cython would be great. Thanks
>
> Christoph
>

The new random c-api was broken in NumPy 1.17. Please try with NumPy 1.18.

Matti