Random Number Generation

class Random123GeneratorBase:
    """
    Base class for Random123 family generators (Philox, Threefry).
    
    Attributes:
    - context (Context): OpenCL context
    - dtype: Output data type
    """

class PhiloxGenerator(Random123GeneratorBase):
    """
    Philox random number generator with excellent statistical properties.
    
    The Philox generator uses a Feistel-like cipher structure and provides
    high-quality random numbers suitable for parallel Monte Carlo simulations.
    """
    
    def __init__(self, context, dtype=np.float32):
        """
        Create Philox random number generator.
        
        Parameters:
        - context (Context): OpenCL context
        - dtype: Output data type (float32, float64, int32, int64)
        """
    
    def fill_uniform(self, queue, ary, luxury=None):
        """
        Fill array with uniform random numbers in [0, 1).
        
        Parameters:
        - queue (CommandQueue): Command queue
        - ary (Array): Target array to fill
        - luxury (int, optional): Quality level (higher = better quality)
        """
    
    def uniform(self, queue, shape, dtype=None, luxury=None):
        """
        Generate uniform random array.
        
        Parameters:
        - queue (CommandQueue): Command queue
        - shape (tuple[int, ...]): Array shape
        - dtype: Data type (uses generator default if None)
        - luxury (int, optional): Quality level
        
        Returns:
        Array: Array of uniform random numbers in [0, 1)
        """

class ThreefryGenerator(Random123GeneratorBase):
    """
    Threefry random number generator based on the Threefish cipher.
    
    The Threefry generator provides high-quality random numbers with
    strong cryptographic properties and excellent parallel performance.
    """
    
    def __init__(self, context, dtype=np.float32):
        """
        Create Threefry random number generator.
        
        Parameters:
        - context (Context): OpenCL context
        - dtype: Output data type
        """
    
    def fill_uniform(self, queue, ary, luxury=None):
        """
        Fill array with uniform random numbers in [0, 1).
        
        Parameters:
        - queue (CommandQueue): Command queue
        - ary (Array): Target array to fill
        - luxury (int, optional): Quality level
        """
    
    def uniform(self, queue, shape, dtype=None, luxury=None):
        """
        Generate uniform random array.
        
        Parameters:
        - queue (CommandQueue): Command queue  
        - shape (tuple[int, ...]): Array shape
        - dtype: Data type
        - luxury (int, optional): Quality level
        
        Returns:
        Array: Array of uniform random numbers in [0, 1)
        """

def rand(queue, shape, dtype=float, luxury=None, generator=None):
    """
    Generate uniform random array in [0, 1).
    
    Parameters:
    - queue (CommandQueue): Command queue
    - shape (int | tuple[int, ...]): Array shape
    - dtype: Data type (float32, float64)
    - luxury (int, optional): Quality level for generator
    - generator (Random123GeneratorBase, optional): Specific generator to use
    
    Returns:
    Array: Array of uniform random numbers
    """

def fill_rand(result, queue=None, luxury=None, generator=None):
    """
    Fill existing array with uniform random numbers.
    
    Parameters:
    - result (Array): Target array to fill
    - queue (CommandQueue, optional): Command queue
    - luxury (int, optional): Quality level
    - generator (Random123GeneratorBase, optional): Specific generator to use
    """

import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.clrandom import PhiloxGenerator, ThreefryGenerator
import pyopencl.clrandom as clrand
import numpy as np

# Setup
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

# Generate uniform random numbers using convenience function
uniform_data = clrand.rand(queue, (10000,), dtype=np.float32)
print(f"Random data shape: {uniform_data.shape}")
print(f"Random samples: {uniform_data.get()[:5]}")
print(f"Range: [{uniform_data.get().min():.3f}, {uniform_data.get().max():.3f}]")

# Fill existing array with random data
existing_array = cl_array.empty(queue, (5000,), dtype=np.float32)
clrand.fill_rand(existing_array, queue)
print(f"Filled array samples: {existing_array.get()[:5]}")

import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.clrandom import PhiloxGenerator, ThreefryGenerator
import numpy as np

# Setup
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

# Create Philox generator
philox_gen = PhiloxGenerator(ctx, dtype=np.float32)

# Generate random array
philox_data = philox_gen.uniform(queue, (10000,))
print(f"Philox generator data: {philox_data.get()[:5]}")

# Create Threefry generator
threefry_gen = ThreefryGenerator(ctx, dtype=np.float32)

# Generate random array
threefry_data = threefry_gen.uniform(queue, (10000,))
print(f"Threefry generator data: {threefry_data.get()[:5]}")

# Fill existing array
target_array = cl_array.empty(queue, (8000,), dtype=np.float32)
philox_gen.fill_uniform(queue, target_array)
print(f"Filled with Philox: {target_array.get()[:5]}")

import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.clrandom import PhiloxGenerator
from pyopencl.elementwise import ElementwiseKernel
import numpy as np

# Setup
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

# Monte Carlo estimation of π using random points in unit square
n_samples = 1000000

# Create generator
generator = PhiloxGenerator(ctx, dtype=np.float32)

# Generate random points
x_coords = generator.uniform(queue, (n_samples,))
y_coords = generator.uniform(queue, (n_samples,))

# Create kernel to check if points are inside unit circle
inside_circle_kernel = ElementwiseKernel(ctx,
    "__global float *x, __global float *y, __global int *inside",
    "inside[i] = (x[i]*x[i] + y[i]*y[i] <= 1.0f) ? 1 : 0",
    "check_inside_circle")

# Count points inside circle
inside_flags = cl_array.empty(queue, (n_samples,), dtype=np.int32)
inside_circle_kernel(x_coords, y_coords, inside_flags)

# Estimate π
points_inside = cl_array.sum(inside_flags).get()
pi_estimate = 4.0 * points_inside / n_samples

print(f"Samples: {n_samples}")
print(f"Points inside circle: {points_inside}")
print(f"π estimate: {pi_estimate}")
print(f"Error: {abs(pi_estimate - np.pi):.6f}")

import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.clrandom import PhiloxGenerator, ThreefryGenerator
import numpy as np
import matplotlib.pyplot as plt

# Setup
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

# Generate large sample for statistical testing
n_samples = 100000

# Test both generators
philox_gen = PhiloxGenerator(ctx, dtype=np.float32)
threefry_gen = ThreefryGenerator(ctx, dtype=np.float32)

philox_data = philox_gen.uniform(queue, (n_samples,)).get()
threefry_data = threefry_gen.uniform(queue, (n_samples,)).get()

# Statistical tests
print("Statistical Analysis:")
print(f"Philox - Mean: {np.mean(philox_data):.6f} (expected: 0.5)")
print(f"Philox - Std:  {np.std(philox_data):.6f} (expected: {1/np.sqrt(12):.6f})")
print(f"Threefry - Mean: {np.mean(threefry_data):.6f}")
print(f"Threefry - Std:  {np.std(threefry_data):.6f}")

# Histogram comparison
bins = np.linspace(0, 1, 50)
philox_hist, _ = np.histogram(philox_data, bins)
threefry_hist, _ = np.histogram(threefry_data, bins)

print(f"Philox histogram uniformity (chi-square): {np.var(philox_hist):.2f}")
print(f"Threefry histogram uniformity (chi-square): {np.var(threefry_hist):.2f}")

import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.clrandom import PhiloxGenerator
from pyopencl.elementwise import ElementwiseKernel
import pyopencl.clmath as clmath
import numpy as np

# Setup
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)
generator = PhiloxGenerator(ctx, dtype=np.float32)

# Generate uniform random numbers as base
uniform1 = generator.uniform(queue, (10000,))
uniform2 = generator.uniform(queue, (10000,))

# Box-Muller transform for normal distribution
normal_kernel = ElementwiseKernel(ctx,
    "__global float *u1, __global float *u2, __global float *n1, __global float *n2",
    """
    float r = sqrt(-2.0f * log(u1[i]));
    float theta = 2.0f * M_PI * u2[i];
    n1[i] = r * cos(theta);
    n2[i] = r * sin(theta);
    """,
    "box_muller_transform")

normal1 = cl_array.empty_like(uniform1)
normal2 = cl_array.empty_like(uniform2)
normal_kernel(uniform1, uniform2, normal1, normal2)

print(f"Normal distribution samples: {normal1.get()[:5]}")
print(f"Normal mean: {cl_array.sum(normal1).get() / normal1.size:.6f}")

# Exponential distribution: -log(1-u)
exponential = -clmath.log(1.0 - uniform1)
print(f"Exponential samples: {exponential.get()[:5]}")

# Integer random numbers in range [0, max_val)
max_val = 100
uniform_int_base = generator.uniform(queue, (10000,))
int_kernel = ElementwiseKernel(ctx,
    "__global float *u, __global int *result, int max_val",
    "result[i] = (int)(u[i] * max_val)",
    "uniform_int")

random_ints = cl_array.empty(queue, (10000,), dtype=np.int32)
int_kernel(uniform_int_base, random_ints, np.int32(max_val))
print(f"Random integers [0, {max_val}): {random_ints.get()[:10]}")

import pyopencl as cl
import pyopencl.array as cl_array
from pyopencl.clrandom import PhiloxGenerator, ThreefryGenerator
import numpy as np
import time

# Setup
ctx = cl.create_some_context()
queue = cl.CommandQueue(ctx)

# Large array for performance testing
n_samples = 10000000
n_iterations = 5

# Test Philox generator
philox_gen = PhiloxGenerator(ctx, dtype=np.float32)
philox_times = []

for i in range(n_iterations):
    start_time = time.time()
    data = philox_gen.uniform(queue, (n_samples,))
    queue.finish()
    end_time = time.time()
    philox_times.append(end_time - start_time)

# Test Threefry generator
threefry_gen = ThreefryGenerator(ctx, dtype=np.float32)
threefry_times = []

for i in range(n_iterations):
    start_time = time.time()
    data = threefry_gen.uniform(queue, (n_samples,))
    queue.finish()
    end_time = time.time()
    threefry_times.append(end_time - start_time)

# Compare with NumPy
numpy_times = []
for i in range(n_iterations):
    start_time = time.time()
    data = np.random.random(n_samples).astype(np.float32)
    end_time = time.time()
    numpy_times.append(end_time - start_time)

print(f"Performance Comparison ({n_samples:,} samples):")
print(f"Philox (GPU):   {np.mean(philox_times):.4f}s ± {np.std(philox_times):.4f}s")
print(f"Threefry (GPU): {np.mean(threefry_times):.4f}s ± {np.std(threefry_times):.4f}s")
print(f"NumPy (CPU):    {np.mean(numpy_times):.4f}s ± {np.std(numpy_times):.4f}s")

# Throughput
philox_throughput = n_samples / np.mean(philox_times) / 1e6
threefry_throughput = n_samples / np.mean(threefry_times) / 1e6
numpy_throughput = n_samples / np.mean(numpy_times) / 1e6

print(f"Throughput (Million samples/sec):")
print(f"Philox:   {philox_throughput:.1f}")
print(f"Threefry: {threefry_throughput:.1f}")
print(f"NumPy:    {numpy_throughput:.1f}")