Python调用cuRandSobol生成Sobol

Sobol低偏差序列在期权的蒙特卡洛估值中非常有用。scipy的qmc模块提供了相关的函数，但是速度比普通的伪随机数要慢。cuRand提供了sobol生成器，速度很快，但是cupy目前还没有集成，scipy下面的qmc模块在cupy中没有对应。Python下想使用cuRandSobol，目前还没有特别好用的模块，所以我们可以用c++编写共享库，然后通过Python的Ctypes方式调用。cuRandSobol和scipy的qmc用的Sobol算法都是JoeKuoD6，支持的最大维度为21201。

首先编写一个gen_sobol.cu文件

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <stdexcept>
#include <vector>

#include <cuda_runtime.h>
#include <curand.h>

// CUDA API error checking
#define CUDA_CHECK(err)                                                        \
  do {                                                                         \
    cudaError_t err_ = (err);                                                  \
    if (err_ != cudaSuccess) {                                                 \
      std::printf("CUDA error %d at %s:%d\n", err_, __FILE__, __LINE__);       \
      throw std::runtime_error("CUDA error");                                  \
    }                                                                          \
  } while (0)

// curand API error checking
#define CURAND_CHECK(err)                                                      \
  do {                                                                         \
    curandStatus_t err_ = (err);                                               \
    if (err_ != CURAND_STATUS_SUCCESS) {                                       \
      std::printf("curand error %d at %s:%d\n", err_, __FILE__, __LINE__);     \
      throw std::runtime_error("curand error");                                \
    }                                                                          \
  } while (0)


using data_type = float;
extern "C"{
void  gen_sobol(const int &n, const unsigned long long &offset,
                 const unsigned int &num_dimensions,                 
                 float* h_data) {

  curandOrdering_t order= CURAND_ORDERING_QUASI_DEFAULT;
  curandRngType_t rng= CURAND_RNG_QUASI_SOBOL32;
  cudaStream_t stream=NULL;
  curandGenerator_t gen=NULL;
  CUDA_CHECK(cudaStreamCreateWithFlags(&stream, cudaStreamNonBlocking));

  /* Create quasi-random number generator */
  CURAND_CHECK(curandCreateGeneratorHost(&gen, CURAND_RNG_QUASI_SOBOL32));

  /* Set cuRAND to stream */
  CURAND_CHECK(curandSetStream(gen, stream));

  /* Set offset */
  CURAND_CHECK(curandSetGeneratorOffset(gen, offset));

  /* Set Dimension */
  CURAND_CHECK(curandSetQuasiRandomGeneratorDimensions(gen, num_dimensions));

  /* Set ordering */
  CURAND_CHECK(curandSetGeneratorOrdering(gen, order));

  /* 可以直接生成正态分布，也可以生成均匀分布 */
 // CURAND_CHECK(curandGenerateUniform(gen, h_data, n*num_dimensions));  
  CURAND_CHECK(curandGenerateNormal(gen, h_data, n*num_dimensions,0,1));

  /* Cleanup */
  CURAND_CHECK(curandDestroyGenerator(gen));
}
}

nvcc -Xcompiler -fPIC -shared -o gen_sobol.so  gen_sobol.cu -lcurand

import numpy as np
import ctypes
from ctypes import *
from scipy.stats import norm
# extract cuda_sum function pointer in the shared object cuda_sum.so
def get_gen_sobol():
    dll = ctypes.CDLL('./gen_sobol.so', mode=ctypes.RTLD_GLOBAL)
    func = dll.gen_sobol
    func.argtypes = [POINTER(c_int), POINTER(c_ulonglong), POINTER(c_uint), POINTER(c_float)]
    return func

# create __cuda_sum function with get_cuda_sum()
__cuda_gen_sobol = get_gen_sobol()

# convenient python wrapper for __cuda_sum
# it does all job with types convertation
# from python ones to C++ ones
def cuda_gen_sobol(n, offset, num_dimensions, h_data):
#    a_p = n.ctypes.data_as(POINTER(c_int))
#    b_p = offset.ctypes.data_as(POINTER(c_ulonglong))
#    c_p = num_dimensions.ctypes.data_as(POINTER(c_uint))
    h_data_p=h_data.ctypes.data_as(POINTER(c_float))

    __cuda_gen_sobol(byref(n), byref(offset), byref(num_dimensions), h_data_p)
# testing, sum of two arrays of ones and output head part of resulting array
if __name__ == '__main__':
    n=65535
    offset=c_ulonglong(n+1)
    num_dimensions=520
    size=(num_dimensions,n)
    h_data = np.zeros(size).astype('float32')

    cuda_gen_sobol(c_int(n), offset, c_uint(num_dimensions), h_data)

    print(h_data)
    #和qmc比较
    dimension = num_undl * num_tp
    soboleng = qmc.Sobol(d=dimension, scramble=False)
    soboleng.fast_forward(n+1)        
    qs =soboleng.random(num_path)
    z=norm.ppf(qs)
    print(z)