[FIXED] numpy implementing a custom RNG

Issue

I'm trying to get numpy to use my own implementation of an RNG in for consistency reasons. My understanding, based on the little documentation I could find, from the numpy docs here and here is that I need to provide a custom BitGenerator class that implements the random_raw method and then initialise using np.random.Generator, so I tried this:

import numpy as np

class TestBitGenerator(np.random.BitGenerator):
  def __init__(self):
    super().__init__(0)
    self.counter = 0
  def random_raw(self, size=None):
    self.counter += 1
    if size is None:
      return self.counter
    return np.full(n, self.counter)  

mc_adapter = TestBitGenerator() 
npgen = np.random.Generator(mc_adapter)

print(npgen.random())

which results in a segfault:

$ python bitgen.py 
Segmentation fault (core dumped)

I assume I'm missing something (from TestBitGenerator?) here, can anyone point me in the right direction?

I tried not subclassing np.random.BitGenerator, and got a different error: object has no attribute 'capsule'

I using numpy 1.19.2 and python 3.8.2

Solution

I did some digging here and Generator's constructor looks like this:

    def __init__(self, bit_generator):
        self._bit_generator = bit_generator

        capsule = bit_generator.capsule
        cdef const char *name = "BitGenerator"
        if not PyCapsule_IsValid(capsule, name):
            raise ValueError("Invalid bit generator. The bit generator must "
                             "be instantiated.")
        self._bitgen = (<bitgen_t *> PyCapsule_GetPointer(capsule, name))[0]
        self.lock = bit_generator.lock

also here explains the requirements for the bitgen_t struct:

struct bitgen:
    void *state
    npy_uint64 (*next_uint64)(void *st) nogil
    uint32_t (*next_uint32)(void *st) nogil
    double (*next_double)(void *st) nogil
    npy_uint64 (*next_raw)(void *st) nogil

ctypedef bitgen bitgen_t

so from what I can make out, it looks like any implementation of BitGenerator must be done in cython, any attempt to implement one in pure python (and probably pybind11 too) just won't work(?)

As I'm not familiar with cython and if/how it could coexist with pybind11, for now I'm just going to ensure each of my (parallel) processes explicitly calls np.random.Generator with a numpy RNG seeded deterministically so that it's independent from my own RNG steam.

Answered By - virgesmith