[SciPy-dev] Monte Carlo package

[Robert Kern]

Ed Schofield wrote:
> Robert Kern wrote:

>>I really would prefer that it not use rand() but use numpy.random instead. What
>>do you need to make this happen?
> 
> Hmmm ... do you mean rk_random() defined in
> numpy/random/mtrand/randomkit.c?

Yes.

> I guess I could link to this instead of
> the system default rand(), but this is probably overkill for this
> application.  I just need a single 30-bit random number per instance, to
> use as a seed.  I could use the system time directly, but at the moment
> I use rand() to scramble it up a little first:
> 
>     /* First seed our faster XORshift RNG. */
>     /* We need to do this in a way that supports non-SVID C platforms.
>      * We use rand() twice, constructing a 30-bit integer out of the
>      * two results. */
>     srand( (unsigned int) time( NULL ));
>     seedhigh = rand();
>     seedlow = rand();
> 
>     /* Construct a 30-bit random integer out of these. */
>     jxr = ((seedhigh & 0x3fff) << 16) + (seedlow & 0xffff);

So there's no way for the user to specify a seed? That's not acceptable for
scientific use. Preferably, I should be able to pass in a
numpy.random.RandomState object, and use the Mersenne Twister generator instead
of the XORshift algorithm. If I can do that, then I don't mind also exposing the
XORshift algorithm, too, *if* I can set the seed for that algorithm.

> I fear that depending on the randomkit.h header would just add an
> unnecessary build dependency.  Do you still see a need for it?

Yes. We can add the appropriate headers to numpy and expose the RandomState
extension object to be used in other extension modules. I should do that
anyways, so this is as good a time as any.

So what do you need from me to make this happen?

-- 
Robert Kern
robert.kern at gmail.com

"In the fields of hell where the grass grows high
 Are the graves of dreams allowed to die."
  -- Richard Harter