Search code examples
pythonnumpyopenclpyopencl

Random number generator with Gaussian distribution for PyOpenCl

Generating Gaussian random numbers using numpy turns out to be the bottleneck in my monte carlo simulation where I make heavy use of PyOpenCl. 

np.random.randn(int(1e9))

Therefore I am looking for a way to generate Gaussian distributed random numbers with PyopenCl too.

I found a 6 years old thread asking a similar question. But I am not sure how to use VexCL library with PyOpenCl: Gaussian distributed random numbers in OpenCL

Any ideas how to implement a good RNG which performs similar as np.random.randn in PyOpenCl?

Solution

  • Seems like pyopencl already includes a random number generator:

    https://github.com/inducer/pyopencl/blob/master/pyopencl/clrandom.py

    Running a simple test shows, that the mean and standard deviation are comparable to numpy's implementation. Also the histogram corresponds to the normal distribution with negligible mean squared error.

    Does anyone know further tests to check the quality of the random number generator?

    Edit: According to https://documen.tician.de/pyopencl/array.html 

    import pytest
from pyopencl.clrandom import RanluxGenerator, PhiloxGenerator
import pyopencl as cl
import pyopencl.array as cl_array
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm


def test_compare_to_numpy_rand():
    context = cl.create_some_context()
    queue = cl.CommandQueue(context)
    rng =  PhiloxGenerator(context=context)
    mu=0
    sigma=3
    size = (int(1e8))
    dtype = np.float32
    numbers_cl = cl_array.zeros(queue, size, dtype=dtype)
    rng.fill_normal(ary=numbers_cl, mu=mu, sigma=sigma)
    numbers_cl_host = numbers_cl.get()
    mu_meas = np.mean(numbers_cl_host)
    sigma_meas = np.std(numbers_cl_host)

    hist_cl, bins_edges = np.histogram(numbers_cl_host, bins=1000, normed=True)
    delta = bins_edges[1] - bins_edges[0]
    bins_center = bins_edges[:-1]+delta/2
    plt.bar(bins_center, hist_cl,width=delta,label='cl')

    numbers_np = mu + sigma * np.random.randn(size)
    hist_np, bins_edges = np.histogram(numbers_np, bins=bins_edges, normed=True)
    plt.bar(bins_center, hist_np, width=delta,alpha=0.8,label='numpy')

    p = norm.pdf(bins_center, mu, sigma)
    plt.plot(bins_center, p, color='red',label='Exact')
    plt.title('Mean squared error: cl={:.5E} np={:.5E}'.format(np.mean(np.abs(p-hist_cl)**2),np.mean(np.abs(p-hist_np)**2)))
    plt.legend()
    plt.show()

    assert pytest.approx(mu_meas, 1e-2) == mu
    assert pytest.approx(sigma_meas, 1e-2) == sigma
    assert pytest.approx(mu_meas, 1e-2) == numbers_np.mean()
    assert pytest.approx(sigma_meas, 1e-2) == numbers_np.std()

    Probability distribution