从 Python 调用快速 C Mersenne Twister 实现 (SFMT)答案

【问题标题】：Calling fast C Mersenne Twister implementation (SFMT) from Python从 Python 调用快速 C Mersenne Twister 实现 (SFMT)
【发布时间】：2017-07-16 11:21:09
【问题描述】：

我正在尝试从 Python 调用 SFMT Mersenne Twister 实现（位于 http://www.math.sci.hiroshima-u.ac.jp/~m-mat/MT/SFMT/）。我这样做是因为我希望能够从具有 4 个概率的离散 pdf 中快速采样。我正在编写一些模拟，这是迄今为止我的代码中最慢的部分。

我设法编写了一些有效的 C 代码，并使用使用 SFMT 算法创建的 [0,1] 中的随机数对输入 PDF 进行采样。

但是我不知道从 Python 调用时如何正确初始化 SFMT 随机数生成器。我当然只想初始化一次，然后我需要将用于初始化它的结构的地址 (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);

问题是，当从 Python 调用此代码时，我不知道如何只为 SFMT 随机数生成器播种一次。如果我只是用 C 编写所有内容，我会在 main() 函数中进行初始化，然后将 &sfmt 参数传递给所有后续的 sfmt_genrand_real1() 调用。

但是因为我正在编译这个 C 代码，然后从 Python 调用它，所以我无法初始化该变量一次。目前，我已经在sfmt_genrand_real1() 调用之前进行了初始化，因为这是我可以让代码甚至编译并能够在调用随机数生成器时访问sfmt 变量的唯一方法。我所有试图以某种方式使sfmt 变量“全局”的尝试都适得其反。

所以我的问题是：有没有办法只初始化 C SFMT 随机数生成器一次，然后在我随后从 Python 对 c_random_sample 的所有调用中访问用于该初始化的 sfmt 结构？

非常感谢任何可以提供帮助的人。

这是我的完整 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
    );
}

【问题讨论】：

这种 PRNG 可能仍然比较新的 PRNG 慢得多（PCG、xoroshiro 等）。这些甚至在 python 中可用（通过外部库）。但更重要的是：也许你应该首先展示你的代码来推理这种性能。除了 PRNG 之外，我可以想象出一些原因来说明这种缓慢的根源。
感谢您的回复萨沙。我正在编写一个模拟器，它需要从依赖于系统当前状态的分布中抽取随机数。所以我不能轻易地提前缓存样本。本质上，Python 对于我所追求的东西来说太慢了，这就是我转向 C 的原因。不过，我会看看你提到的那些 PRNG，谢谢。
我刚刚对 SFMT 与 PCG 进行了快速测试。 SFMT 创建 10 亿个 uint32_t 随机数的速度似乎快了两倍。
虽然xoroshiro.di.unimi.it 的速度测试似乎表明它可以很快。也许我做错了什么。感谢您向我指出这些，萨沙。通过在 C 中定义一个 Python 类并在我的类的构造函数中初始化随机数生成器，我设法克服了我的初始化问题。
我明白了。我预计离散采样的使用会有些错误（例如，如果使用 O(1) 采样，则每次都要为设置付费）。但我现在看到，您的案例是高度动态的，变化很大。请记住，Python 在 core + numpy 中的 PRNG 实现也是 C 代码。但是，如果做得正确，相信你的基准肯定没有错。

标签： random python-c-api mersenne-twister

【解决方案1】：

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

【讨论】：

【解决方案2】：

这不是您正在寻找的，但如果您不介意依赖项（主要是 Numpy 和 Cython），您可能会发现包含许多快速随机数生成器（特别是 dSFMT）的 ng-numpy-randomstate 库很有用.

如果您可以使用它，它将为您节省编写 C 包装器的工作。

此外，这个模块（或至少是某些功能？）将来可能会合并到 Numpy 中（参见 Numpy 开放问题 #6967）。

【讨论】：