[SciPy-Dev] Generate random variates using Cython

[Andrew Nelson]

>
> After upgrading to Numpy 1.18, I could indeed compile the code
> https://docs.scipy.org/doc/numpy-1.17.0/reference/random/extending.html
> and it gives a substantial speedup.
> I am still wondering how to integrate such an approach into SciPy. To
> generate uniform samples in the _rvs method of a continuous distribution i
> scipy.stats, I would write:
>
> self._random_state.random_sample
>

`np.random.Generator` doesn't have a `random_sample` method. But it does
have a `uniform` method. So if `_random_state` is either a `RandomState` or
a `Generator` instance you should use the `uniform` method for it to work
across both.

If I want to use Cython, how can I ensure that it relies on the same random
> state?
>

`Generator` and `RandomState` will produce different streams of random
numbers from the same seed. In the PR linked to issue #11402 I'm trying to
ensure that the `rvs` methods will work with both `Generator` and
`RandomState`. The latter is current default, and one would want to use
that to retain back compatibility of random numbers. If you're trying to
use Cython to speedup a particular rvs method you'll need to have two
branches. One branch is given a RandomState instance, the other a Generator
instance.


import numpy as np
> cimport numpy as np
> cimport cython
> from cpython.pycapsule cimport PyCapsule_IsValid, PyCapsule_GetPointer
> from numpy.random cimport bitgen_t
> from numpy.random import PCG64
>
> @cython.boundscheck(False)
> @cython.wraparound(False)
> def rvs(Py_ssize_t n):
>     cdef Py_ssize_t i
>     cdef bitgen_t *rng
>     cdef const char *capsule_name = "BitGenerator"
>     cdef double u, v
>     cdef double[::1] random_values
>     x = PCG64()
>     capsule = x.capsule
>     # Optional check that the capsule if from a BitGenerator
>     if not PyCapsule_IsValid(capsule, capsule_name):
>         raise ValueError("Invalid pointer to anon_func_state")
>     # Cast the pointer
>     rng = <bitgen_t *> PyCapsule_GetPointer(capsule, capsule_name)
>     random_values = np.empty(n, dtype='float64')
>     with x.lock, nogil:
>         for i in range(n):
>             while (1):
>                 u = rng.next_double(rng.state)
>                 v = rng.next_double(rng.state)
>                 if u <= v*v:
>                     break
>             random_values[i] = v
>     randoms = np.asarray(random_values)
>     return randoms
>
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.python.org/pipermail/scipy-dev/attachments/20200123/049dc42a/attachment.html>