pyopencl

pyopencl copied to clipboard

The fill_normal() method gives inf values sometimes, probably because the random generators sometimes yield zeros especially with 32 bit floats. The logarithm in Box-Muller then causes the infinity.

The following output shows that although zeros are much less frequent than other values such as ones, they occassionally come up in fill_uniform() , and so do infinities in fill_normal(). No difference between Threefry and Philox.

Output on my system, running pyopencl 2017.2.2:

3 zeros found in uniform
570 ones founds in uniform
2 infs found in normal

The code:

from pyopencl import clrandom
import pyopencl as cl
import numpy as np

context = cl.create_some_context(interactive=False)
queue = cl.CommandQueue(context)

#random_gen = clrandom.PhiloxGenerator(context)
random_gen = clrandom.ThreefryGenerator(context)

shape = (1000000,)
arr = cl.array.empty(queue, shape, dtype=np.float32)

n_zeros = 0
n_ones = 0
n_infs = 0
for i in range(20000):
	random_gen.fill_uniform(arr).wait()
	vals = arr.get()
	n_zeros += np.count_nonzero(vals == 0.0)
	n_ones += np.count_nonzero(vals == 1.0)
	
	random_gen.fill_normal(arr).wait()
	vals = arr.get()
	n_infs += np.count_nonzero(np.isinf(vals))
	
print(n_zeros, 'zeros found in uniform')
print(n_ones, 'ones founds in uniform')
print(n_infs, 'infs found in normal')

I notice the same issue, but only when using np.float32. Using np.float64 does not result in inf values. Happens on Intel and AMD devices with pyopencl version 2019.1.2.

The script below reproduces the inf issue immediatly when using seed=393:

import pyopencl.array as cl_array
from pyopencl.clrandom import PhiloxGenerator
import numpy as np

context = cl.create_some_context()
queue = cl.CommandQueue(context)

noise_buff = cl.array.zeros(queue, (10000, 1000), dtype=np.float32)
rng = PhiloxGenerator(context=context, seed=393)
b_no_inf_found = True
while b_no_inf_found:
    rng.fill_normal(noise_buff, queue=queue, mu=1, sigma=0.00001)
    if np.any(np.abs(noise_buff.get()) > 10000):
        b_no_inf_found = False
    queue.finish()
print('Inf value found at positions {}'.format(str(np.where(np.abs(noise_buff.get()) == np.inf)[0])))

A 64 bit float has 2^32 times more possible values than a 32 bit one so encountering a zero is much less likely. With enough iterations, one will probably encounter zeros with 64 bit floats as well.

I am running into the same issue, and as the original poster notes, this is coming from the fact that the fill_uniform method can generate zeros which the Box-Mueller transform does not handle. This could be handled as in numpy (see https://github.com/numpy/numpy/blob/main/numpy/random/src/distributions/distributions.c line 137), where the special case of taking the log of a random variable (r.v) equal to 0 is transformed to taking the log of (1-0).

Additionally, one might think to extend pyopencl capabilities by adding other distributions. Many of the algorithms to generate other distributions (for example, exponential distribution, gamma distribution) expect the r.v. to be distributed inside the half-open interval [0,1), which is not the case for the r.v generated with fill_uniform where the r.v can be 1 in some cases, because (2^32-1)/2^32 equals 1 in float32. One suggestion to remedy this problem is again to borrow from numpy (line 19 in previous link), where (in the case of float32 precision) the uint32 integer is first downgraded to a 24 bit unsigned integer before being divided by 2^24, yielding r.v. < 1.