Вызов быстрой реализации C Mersenne Twister (SFMT) из Python

Question

Я пытаюсь вызвать реализацию SFMT Mersenne Twister (находится по адресу http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/) из Python. Я делаю это, потому что я хотел бы иметь возможность быстро выбрать дискретный PDF с 4 вероятностями. Я пишу некоторые симуляции, и это самая медленная часть моего кода.

Мне удалось написать некоторый код на C, который работает, и сэмплировать входной PDF, используя случайное число в [0,1], созданное с помощью алгоритма SFMT.

Однако я не знаю, как правильно инициализировать генератор случайных чисел SFMT при вызове из Python. Я хочу инициализировать его только один раз, а затем мне нужно передать адрес структуры, используемой для инициализации (sfmt) в мои звонки sfmt_genrand_real1,

Так что пример кода будет:

// 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);

Проблема в том, что я не знаю, как запустить генератор случайных чисел SFMT только один раз при вызове этого кода из Python. Если бы я просто писал все на C, я бы сделал инициализацию в функции main(), а затем передал бы &sfmt аргумент всему последующему sfmt_genrand_real1() звонки.

Но поскольку я компилирую этот код C, а затем вызываю его из Python, я не могу инициализировать эту переменную один раз. Пока я застрял инициализацию прямо перед sfmt_genrand_real1() вызов, потому что это единственный способ, которым я мог заставить код даже скомпилировать и иметь доступ к sfmt переменная в обращении к генератору случайных чисел. Все мои попытки как-то сделать sfmt переменная "глобальная" имеет обратный эффект.

Итак, мой вопрос: есть ли способ инициализировать генератор случайных чисел C SFMT только один раз, а затем получить доступ к sfmt Использование структуры для этой инициализации во всех моих последующих вызовах c_random_sample из питона?

Большое спасибо заранее всем, кто может помочь с этим.

Вот мой C-код в полном объеме. (Для компиляции вам нужно иметь все SMFT .c а также .h файлы в той же папке, а затем скомпилировать с помощью 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
    );
}

Answer 1

Это не совсем то, что вы ищете, но если вы не возражаете против зависимостей (в основном Numpy и Cython), вам может пригодиться библиотека https://github.com/bashtage/ng-numpy-randomstate, которая содержит много быстрых генераторов случайных чисел (в частности, dSFMT).

Если вы можете использовать его, это избавит вас от написания оболочек на языке C.

Кроме того, этот модуль (или, по крайней мере, некоторые функции?) Может быть объединен с Numpy в будущем (см. Открытый выпуск Numpy № 6967).

Answer 2

sfmt идет в глобальном масштабе в C:

static sfmt_t sfmt; /* outside any functions */

Инициализация идет в вашем модуле init функция - вызывается один раз, когда модуль впервые импортируется в Python:

PyMODINIT_FUNC initc_random_sample(void)
{
    /* Py_InitModule4 to initialize the module as before */ 

    int seed = time(NULL);
    sfmt_init_gen_rand(&sfmt, seed);
}