I'm trying to call the SFMT Mersenne Twister implementation (found at http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/) from Python. I'm doing this because I'd like to be able to quickly sample from a discrete pdf with 4 probabilities. I'm writing some simulations and this is by far the slowest part of my code.
I've managed to write some C code which works, and samples the input PDF using a random number in [0,1] created with the SFMT algorithm.
However I don't know how to initialise the SFMT random number generator properly when calling from Python. I only want to initialise it once of course, and then I need to pass the address of the struct used to intialise it (sfmt) into my calls to sfmt_genrand_real1.
So some example code would be:
// Create a struct which stores the Mersenne Twister's state
sfmt_t sfmt;
// Initialise the MT with a random seed from the system time
// NOTE: Only want to do this once
int seed = time(NULL);
sfmt_init_gen_rand(&sfmt, seed);
// Get a random number
double random_number = sfmt_genrand_real1(&sfmt);
The problem is that I don't know how to only seed the SFMT random number generator once when calling this code from Python. If I were just writing everything in C, I'd do the intialisation in the main() function, then pass the &sfmt argument to all subsequent sfmt_genrand_real1() calls.
But because I'm compiling this C code, then calling it from Python, I can't initialise that variable once. For now, I've stuck the initialisation right before the sfmt_genrand_real1() call, because it's the only way I could get the code to even compile and be able to access the sfmt variable in the call to the random number generator. All my attempts to somehow make the sfmt variable "global" backfired.
So my question is: is there a way to initialise the C SFMT random number generator only once, then access the sfmt struct use for that initialisation in all my subsequent calls to c_random_sample from Python?
Many thanks in advance to anyone who can help with this.
Here's my C code in full. (To compile you'd need to have all the SMFT .c and .h files in the same folder, then compile using python setup.py build_ext --inplace)
#include "Python.h"
#include <stdio.h>
#include <time.h>
#include "SFMT.h"
static double*
get_double_array(PyObject* data)
{
int i, size;
double* out;
PyObject* seq;
seq = PySequence_Fast(data, "expected a sequence");
if (!seq)
return NULL;
size = PySequence_Size(seq);
if (size < 0)
return NULL;
out = (double*) PyMem_Malloc(size * sizeof(double));
if (!out) {
Py_DECREF(seq);
PyErr_NoMemory();
return NULL;
}
for (i=0; i<size; i++) {
out[i] = PyFloat_AsDouble(PySequence_Fast_GET_ITEM(seq, i));
}
Py_DECREF(seq);
if (PyErr_Occurred()) {
PyMem_Free(out);
out = NULL;
}
return out;
}
static PyObject*
c_random_sample(PyObject* self, PyObject* args)
{
int i;
double* pdf;
PyObject* pdf_in;
if (!PyArg_ParseTuple(args, "O:c_random_sample", &pdf_in))
return NULL;
pdf = get_double_array(pdf_in);
if (!pdf)
return NULL;
// Create a struct which stores the Mersenne Twister's state
sfmt_t sfmt;
int seed = time(NULL);
// Initialise the MT with a random seed from the system time
sfmt_init_gen_rand(&sfmt, seed);
// NOTE: This simply re-initialises the random number generator
// on every call. We need to only initialise it once...
double r = sfmt_genrand_real1(&sfmt);
for (i=0; i<4; i++) {
r -= pdf[i];
if (r < 0) {
break;
}
}
PyMem_Free(pdf);
return PyInt_FromLong(i);
}
static PyMethodDef functions[] = {
{"c_random_sample", c_random_sample, METH_VARARGS},
{NULL, NULL}
};
PyMODINIT_FUNC initc_random_sample(void)
{
Py_InitModule4(
"c_random_sample", functions, "Trying to implement random number sampling in C", NULL, PYTHON_API_VERSION
);
}
