[Pyrex] using numpy arrays and random number generators in Pyrex

[Robert Kern]

Mark P. Miller wrote:
> Thanks again Robert:
> 
> Robert Kern wrote:
>> If you have to do many iterations and generate one random number each iterate
>> with different parameters each time (and unknown until you get to that iterate),
>> then pyrexGSL might be useful. MCMC is one such application. Of course, you can
>> also extract the C and Pyrex declarations from numpy/random/mtrand/ for the same
>> purpose. Depending on your needs, that might be a better option: that code is
>> BSD-licensed and smaller and more self-contained than the GSL. Of course, I'm
>> biased; I wrote the non-uniform generators for it.
> 
> Working somehow with mtrand makes some sense to me...but the specific 
> implementation is not yet clear.  In general, my pyrex functions look 
> something like this:
> 
> import numpy as np
> 
> def function1(int var1, int var2, float var3, float mean):
>     cdef float a
>     cdef float result
>     a = np.random.exponential(mean)
>     #do a bunch of additional calculations with a and other passed
>     #function parameters
>     return result
> 
> There are many functions like this that get called at different places 
> throughout the code.  In general, it is not known in advance when a 
> particular function will get called, nor is it always known exactly how 
> many times each function will get called.  For this reason, I have shied 
> away from generating lots of random numbers in advance (that and the 
> fact that I don't completely understand how to efficiently iterate 
> through and use elements from a returned array of random numbers).

This is precisely the situation I describe in the above paragraph where pyrexGSL
or stealing the source code from mtrand would be useful. The recipe for the
latter would be something like this: steal the C source files randomkit.[ch],
distributions.[ch], initarray.[ch]; steal the Pyrex declarations for those
functions from mtrand.pyx; ignore numpy.random.

You will want to organize your module such that all of the functions that call
the random number generators are actually methods on a class that maintains an
rk_state structure.


cdef class MyClass:
    cdef rk_state state

    def function1(self, int var1, int var2, double var3, double mean):
        cdef double a
        cdef double result
        a = rk_exponential(self.state, mean)
        #do a bunch of additional calculations with a and other passed
        #function parameters
        return result

-- 
Robert Kern

"I have come to believe that the whole world is an enigma, a harmless enigma
 that is made terrible by our own mad attempt to interpret it as though it had
 an underlying truth."
  -- Umberto Eco