summaryrefslogtreecommitdiff
path: root/tests
diff options
context:
space:
mode:
Diffstat (limited to 'tests')
-rw-r--r--tests/__init__.py0
-rw-r--r--tests/data/ray_intersection.npzbin1920080 -> 0 bytes
-rw-r--r--tests/linalg_test.cu128
-rw-r--r--tests/linalg_test.py214
-rw-r--r--tests/matrix_test.cu152
-rw-r--r--tests/matrix_test.py327
-rw-r--r--tests/rotate_test.cu14
-rw-r--r--tests/rotate_test.py67
-rw-r--r--tests/test_fileio.py74
-rw-r--r--tests/test_generator_photon.py21
-rw-r--r--tests/test_generator_vertex.py23
-rw-r--r--tests/test_pdf.py67
-rw-r--r--tests/test_propagation.py58
-rw-r--r--tests/test_ray_intersection.py27
-rw-r--r--tests/test_rayleigh.py57
-rw-r--r--tests/test_sample_cdf.cu16
-rw-r--r--tests/test_sample_cdf.py64
17 files changed, 0 insertions, 1309 deletions
diff --git a/tests/__init__.py b/tests/__init__.py
deleted file mode 100644
index e69de29..0000000
--- a/tests/__init__.py
+++ /dev/null
diff --git a/tests/data/ray_intersection.npz b/tests/data/ray_intersection.npz
deleted file mode 100644
index fba0e40..0000000
--- a/tests/data/ray_intersection.npz
+++ /dev/null
Binary files differ
diff --git a/tests/linalg_test.cu b/tests/linalg_test.cu
deleted file mode 100644
index 4e9c983..0000000
--- a/tests/linalg_test.cu
+++ /dev/null
@@ -1,128 +0,0 @@
-//-*-c-*-
-
-#include "linalg.h"
-
-extern "C"
-{
-
-__global__ void float3add(float3 *a, float3 *b, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = a[idx] + b[idx];
-}
-
-__global__ void float3addequal(float3 *a, float3 *b)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- a[idx] += b[idx];
-}
-
-__global__ void float3sub(float3 *a, float3 *b, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = a[idx] - b[idx];
-}
-
-__global__ void float3subequal(float3 *a, float3 *b)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- a[idx] -= b[idx];
-}
-
-__global__ void float3addfloat(float3 *a, float c, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = a[idx] + c;
-}
-
-__global__ void float3addfloatequal(float3 *a, float c)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- a[idx] += c;
-}
-
-__global__ void floataddfloat3(float3 *a, float c, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = c + a[idx];
-}
-
-__global__ void float3subfloat(float3 *a, float c, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = a[idx] - c;
-}
-
-__global__ void float3subfloatequal(float3 *a, float c)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- a[idx] -= c;
-}
-
-__global__ void floatsubfloat3(float3 *a, float c, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = c - a[idx];
-}
-
-__global__ void float3mulfloat(float3 *a, float c, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = a[idx]*c;
-}
-
-__global__ void float3mulfloatequal(float3 *a, float c)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- a[idx] *= c;
-}
-
-__global__ void floatmulfloat3(float3 *a, float c, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = c*a[idx];
-}
-
-__global__ void float3divfloat(float3 *a, float c, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = a[idx]/c;
-}
-
-__global__ void float3divfloatequal(float3 *a, float c)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- a[idx] /= c;
-}
-
-__global__ void floatdivfloat3(float3 *a, float c, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = c/a[idx];
-}
-
-__global__ void dot(float3 *a, float3 *b, float *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = dot(a[idx],b[idx]);
-}
-
-__global__ void cross(float3 *a, float3 *b, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = cross(a[idx],b[idx]);
-}
-
-__global__ void norm(float3 *a, float *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = norm(a[idx]);
-}
-
-__global__ void minusfloat3(float3 *a, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = -a[idx];
-}
-
-} // extern "c"
diff --git a/tests/linalg_test.py b/tests/linalg_test.py
deleted file mode 100644
index 31688d9..0000000
--- a/tests/linalg_test.py
+++ /dev/null
@@ -1,214 +0,0 @@
-import os
-import numpy as np
-from pycuda import autoinit
-from pycuda.compiler import SourceModule
-import pycuda.driver as cuda
-from pycuda import gpuarray
-
-float3 = gpuarray.vec.float3
-
-print 'device %s' % autoinit.device.name()
-
-current_directory = os.path.split(os.path.realpath(__file__))[0]
-source_directory = current_directory + '/../src'
-
-source = open(current_directory + '/linalg_test.cu').read()
-
-mod = SourceModule(source, options=['-I' + source_directory], no_extern_c=True, cache_dir=False)
-
-float3add = mod.get_function('float3add')
-float3addequal = mod.get_function('float3addequal')
-float3sub = mod.get_function('float3sub')
-float3subequal = mod.get_function('float3subequal')
-float3addfloat = mod.get_function('float3addfloat')
-float3addfloatequal = mod.get_function('float3addfloatequal')
-floataddfloat3 = mod.get_function('floataddfloat3')
-float3subfloat = mod.get_function('float3subfloat')
-float3subfloatequal = mod.get_function('float3subfloatequal')
-floatsubfloat3 = mod.get_function('floatsubfloat3')
-float3mulfloat = mod.get_function('float3mulfloat')
-float3mulfloatequal = mod.get_function('float3mulfloatequal')
-floatmulfloat3 = mod.get_function('floatmulfloat3')
-float3divfloat = mod.get_function('float3divfloat')
-float3divfloatequal = mod.get_function('float3divfloatequal')
-floatdivfloat3 = mod.get_function('floatdivfloat3')
-dot = mod.get_function('dot')
-cross = mod.get_function('cross')
-norm = mod.get_function('norm')
-minusfloat3 = mod.get_function('minusfloat3')
-
-size = {'block': (256,1,1), 'grid': (1,1)}
-
-a = np.empty(size['block'][0], dtype=float3)
-b = np.empty(size['block'][0], dtype=float3)
-c = np.float32(np.random.random_sample())
-
-a['x'] = np.random.random_sample(size=a.size)
-a['y'] = np.random.random_sample(size=a.size)
-a['z'] = np.random.random_sample(size=a.size)
-
-b['x'] = np.random.random_sample(size=b.size)
-b['y'] = np.random.random_sample(size=b.size)
-b['z'] = np.random.random_sample(size=b.size)
-
-def testfloat3add():
- dest = np.empty(a.size, dtype=float3)
- float3add(cuda.In(a), cuda.In(b), cuda.Out(dest), **size)
- if not np.allclose(a['x']+b['x'], dest['x']) or \
- not np.allclose(a['y']+b['y'], dest['y']) or \
- not np.allclose(a['z']+b['z'], dest['z']):
- assert False
-
-def testfloat3sub():
- dest = np.empty(a.size, dtype=float3)
- float3sub(cuda.In(a), cuda.In(b), cuda.Out(dest), **size)
- if not np.allclose(a['x']-b['x'], dest['x']) or \
- not np.allclose(a['y']-b['y'], dest['y']) or \
- not np.allclose(a['z']-b['z'], dest['z']):
- assert False
-
-def testfloat3addequal():
- dest = np.copy(a)
- float3addequal(cuda.InOut(dest), cuda.In(b), **size)
- if not np.allclose(a['x']+b['x'], dest['x']) or \
- not np.allclose(a['y']+b['y'], dest['y']) or \
- not np.allclose(a['z']+b['z'], dest['z']):
- assert False
-
-def testfloat3subequal():
- dest = np.copy(a)
- float3subequal(cuda.InOut(dest), cuda.In(b), **size)
- if not np.allclose(a['x']-b['x'], dest['x']) or \
- not np.allclose(a['y']-b['y'], dest['y']) or \
- not np.allclose(a['z']-b['z'], dest['z']):
- assert False
-
-def testfloat3addfloat():
- dest = np.empty(a.size, dtype=float3)
- float3addfloat(cuda.In(a), c, cuda.Out(dest), **size)
- if not np.allclose(a['x']+c, dest['x']) or \
- not np.allclose(a['y']+c, dest['y']) or \
- not np.allclose(a['z']+c, dest['z']):
- assert False
-
-def testfloat3addfloatequal():
- dest = np.copy(a)
- float3addfloatequal(cuda.InOut(dest), c, **size)
- if not np.allclose(a['x']+c, dest['x']) or \
- not np.allclose(a['y']+c, dest['y']) or \
- not np.allclose(a['z']+c, dest['z']):
- assert False
-
-def testfloataddfloat3():
- dest = np.empty(a.size, dtype=float3)
- floataddfloat3(cuda.In(a), c, cuda.Out(dest), **size)
- if not np.allclose(c+a['x'], dest['x']) or \
- not np.allclose(c+a['y'], dest['y']) or \
- not np.allclose(c+a['z'], dest['z']):
- assert False
-
-def testfloat3subfloat():
- dest = np.empty(a.size, dtype=float3)
- float3subfloat(cuda.In(a), c, cuda.Out(dest), **size)
- if not np.allclose(a['x']-c, dest['x']) or \
- not np.allclose(a['y']-c, dest['y']) or \
- not np.allclose(a['z']-c, dest['z']):
- assert False
-
-def testfloat3subfloatequal():
- dest = np.copy(a)
- float3subfloatequal(cuda.InOut(dest), c, **size)
- if not np.allclose(a['x']-c, dest['x']) or \
- not np.allclose(a['y']-c, dest['y']) or \
- not np.allclose(a['z']-c, dest['z']):
- assert False
-
-def testfloatsubfloat3():
- dest = np.empty(a.size, dtype=float3)
- floatsubfloat3(cuda.In(a), c, cuda.Out(dest), **size)
- if not np.allclose(c-a['x'], dest['x']) or \
- not np.allclose(c-a['y'], dest['y']) or \
- not np.allclose(c-a['z'], dest['z']):
- assert False
-
-def testfloat3mulfloat():
- dest = np.empty(a.size, dtype=float3)
- float3mulfloat(cuda.In(a), c, cuda.Out(dest), **size)
- if not np.allclose(a['x']*c, dest['x']) or \
- not np.allclose(a['y']*c, dest['y']) or \
- not np.allclose(a['z']*c, dest['z']):
- assert False
-
-def testfloat3mulfloatequal():
- dest = np.copy(a)
- float3mulfloatequal(cuda.InOut(dest), c, **size)
- if not np.allclose(a['x']*c, dest['x']) or \
- not np.allclose(a['y']*c, dest['y']) or \
- not np.allclose(a['z']*c, dest['z']):
- assert False
-
-def testfloatmulfloat3():
- dest = np.empty(a.size, dtype=float3)
- floatmulfloat3(cuda.In(a), c, cuda.Out(dest), **size)
- if not np.allclose(c*a['x'], dest['x']) or \
- not np.allclose(c*a['y'], dest['y']) or \
- not np.allclose(c*a['z'], dest['z']):
- assert False
-
-def testfloat3divfloat():
- dest = np.empty(a.size, dtype=float3)
- float3divfloat(cuda.In(a), c, cuda.Out(dest), **size)
- if not np.allclose(a['x']/c, dest['x']) or \
- not np.allclose(a['y']/c, dest['y']) or \
- not np.allclose(a['z']/c, dest['z']):
- assert False
-
-def testfloat3divfloatequal():
- dest = np.copy(a)
- float3divfloatequal(cuda.InOut(dest), c, **size)
- if not np.allclose(a['x']/c, dest['x']) or \
- not np.allclose(a['y']/c, dest['y']) or \
- not np.allclose(a['z']/c, dest['z']):
- assert False
-
-def testfloatdivfloat3():
- dest = np.empty(a.size, dtype=float3)
- floatdivfloat3(cuda.In(a), c, cuda.Out(dest), **size)
- if not np.allclose(c/a['x'], dest['x']) or \
- not np.allclose(c/a['y'], dest['y']) or \
- not np.allclose(c/a['z'], dest['z']):
- assert false
-
-def testdot():
- dest = np.empty(a.size, dtype=np.float32)
- dot(cuda.In(a), cuda.In(b), cuda.Out(dest), **size)
- if not np.allclose(a['x']*b['x'] + a['y']*b['y'] + a['z']*b['z'], dest):
- assert False
-
-def testcross():
- dest = np.empty(a.size, dtype=float3)
- cross(cuda.In(a), cuda.In(b), cuda.Out(dest), **size)
- for u, v, wdest in zip(a,b,dest):
- w = np.cross((u['x'], u['y'], u['z']),(v['x'],v['y'],v['z']))
- if not np.allclose(wdest['x'], w[0]) or \
- not np.allclose(wdest['y'], w[1]) or \
- not np.allclose(wdest['z'], w[2]):
- print w
- print wdest
- assert False
-
-def testnorm():
- dest = np.empty(a.size, dtype=np.float32)
- norm(cuda.In(a), cuda.Out(dest), **size)
-
- for i in range(len(dest)):
- if not np.allclose(np.linalg.norm((a['x'][i],a['y'][i],a['z'][i])), dest[i]):
- assert False
-
-def testminusfloat3():
- dest = np.empty(a.size, dtype=float3)
- minusfloat3(cuda.In(a), cuda.Out(dest), **size)
- if not np.allclose(-a['x'], dest['x']) or \
- not np.allclose(-a['y'], dest['y']) or \
- not np.allclose(-a['z'], dest['z']):
- assert False
diff --git a/tests/matrix_test.cu b/tests/matrix_test.cu
deleted file mode 100644
index d64cb34..0000000
--- a/tests/matrix_test.cu
+++ /dev/null
@@ -1,152 +0,0 @@
-//-*-c-*-
-
-#include "matrix.h"
-
-__device__ Matrix array2matrix(float *a)
-{
- return make_matrix(a[0], a[1], a[2],
- a[3], a[4], a[5],
- a[6], a[7], a[8]);
-}
-
-__device__ void matrix2array(const Matrix &m, float *a)
-{
- a[0] = m.a00;
- a[1] = m.a01;
- a[2] = m.a02;
- a[3] = m.a10;
- a[4] = m.a11;
- a[5] = m.a12;
- a[6] = m.a20;
- a[7] = m.a21;
- a[8] = m.a22;
-}
-
-extern "C"
-{
-
-__global__ void det(float *a, float *dest)
-{
- Matrix m = array2matrix(a);
- dest[0] = det(m);
-}
-
-__global__ void inv(float *a, float *dest)
-{
- Matrix m = array2matrix(a);
- matrix2array(inv(m), dest);
-}
-
-__global__ void minusmatrix(float *a, float *dest)
-{
- matrix2array(-array2matrix(a), dest);
-}
-
-__global__ void matrixadd(float *a, float *b, float *dest)
-{
- matrix2array(array2matrix(a)+array2matrix(b), dest);
-}
-
-__global__ void matrixsub(float *a, float *b, float *dest)
-{
- matrix2array(array2matrix(a)-array2matrix(b), dest);
-}
-
-__global__ void matrixmul(float *a, float *b, float *dest)
-{
- matrix2array(array2matrix(a)*array2matrix(b), dest);
-}
-
-__global__ void multiply(float *a, float3 *x, float3 *dest)
-{
- dest[0] = array2matrix(a)*x[0];
-}
-
-__global__ void matrixaddfloat(float *a, float c, float *dest)
-{
- matrix2array(array2matrix(a)+c, dest);
-}
-
-__global__ void matrixsubfloat(float *a, float c, float *dest)
-{
- matrix2array(array2matrix(a)-c, dest);
-}
-
-__global__ void matrixmulfloat(float *a, float c, float *dest)
-{
- matrix2array(array2matrix(a)*c, dest);
-}
-
-__global__ void matrixdivfloat(float *a, float c, float *dest)
-{
- matrix2array(array2matrix(a)/c, dest);
-}
-
-__global__ void floataddmatrix(float *a, float c, float *dest)
-{
- matrix2array(c+array2matrix(a), dest);
-}
-
-__global__ void floatsubmatrix(float *a, float c, float *dest)
-{
- matrix2array(c-array2matrix(a), dest);
-}
-
-__global__ void floatmulmatrix(float *a, float c, float *dest)
-{
- matrix2array(c*array2matrix(a), dest);
-}
-
-__global__ void floatdivmatrix(float *a, float c, float *dest)
-{
- matrix2array(c/array2matrix(a), dest);
-}
-
-__global__ void matrixaddequals(float *a, float *b)
-{
- Matrix m = array2matrix(a);
- m += array2matrix(b);
- matrix2array(m,a);
-}
-
-__global__ void matrixsubequals(float *a, float *b)
-{
- Matrix m = array2matrix(a);
- m -= array2matrix(b);
- matrix2array(m,a);
-}
-
-__global__ void matrixaddequalsfloat(float *a, float c)
-{
- Matrix m = array2matrix(a);
- m += c;
- matrix2array(m,a);
-}
-
-__global__ void matrixsubequalsfloat(float *a, float c)
-{
- Matrix m = array2matrix(a);
- m -= c;
- matrix2array(m,a);
-}
-
-__global__ void matrixmulequalsfloat(float *a, float c)
-{
- Matrix m = array2matrix(a);
- m *= c;
- matrix2array(m,a);
-}
-
-__global__ void matrixdivequalsfloat(float *a, float c)
-{
- Matrix m = array2matrix(a);
- m /= c;
- matrix2array(m,a);
-}
-
-__global__ void outer(float3 a, float3 b, float* dest)
-{
- matrix2array(outer(a,b), dest);
-}
-
-} // extern "c"
diff --git a/tests/matrix_test.py b/tests/matrix_test.py
deleted file mode 100644
index c843025..0000000
--- a/tests/matrix_test.py
+++ /dev/null
@@ -1,327 +0,0 @@
-import os
-import numpy as np
-from pycuda import autoinit
-from pycuda.compiler import SourceModule
-import pycuda.driver as cuda
-from pycuda import gpuarray
-
-float3 = gpuarray.vec.float3
-
-print 'device %s' % autoinit.device.name()
-
-current_directory = os.path.split(os.path.realpath(__file__))[0]
-source_directory = current_directory + '/../src'
-
-source = open(current_directory + '/matrix_test.cu').read()
-
-mod = SourceModule(source, options=['-I' + source_directory], no_extern_c=True, cache_dir=False)
-
-det = mod.get_function('det')
-inv = mod.get_function('inv')
-matrixadd = mod.get_function('matrixadd')
-matrixsub = mod.get_function('matrixsub')
-matrixmul = mod.get_function('matrixmul')
-multiply = mod.get_function('multiply')
-matrixaddfloat = mod.get_function('matrixaddfloat')
-matrixsubfloat = mod.get_function('matrixsubfloat')
-matrixmulfloat = mod.get_function('matrixmulfloat')
-matrixdivfloat = mod.get_function('matrixdivfloat')
-floataddmatrix = mod.get_function('floataddmatrix')
-floatsubmatrix = mod.get_function('floatsubmatrix')
-floatmulmatrix = mod.get_function('floatmulmatrix')
-floatdivmatrix = mod.get_function('floatdivmatrix')
-matrixaddequals = mod.get_function('matrixaddequals')
-matrixsubequals = mod.get_function('matrixsubequals')
-matrixaddequalsfloat = mod.get_function('matrixaddequalsfloat')
-matrixsubequalsfloat = mod.get_function('matrixsubequalsfloat')
-matrixmulequalsfloat = mod.get_function('matrixmulequalsfloat')
-matrixdivequalsfloat = mod.get_function('matrixdivequalsfloat')
-outer = mod.get_function('outer')
-minusmatrix = mod.get_function('minusmatrix')
-
-size = {'block': (1,1,1), 'grid': (1,1)}
-
-for i in range(1):
- a = np.random.random_sample(size=9).astype(np.float32)
- b = np.random.random_sample(size=9).astype(np.float32)
- dest = np.empty(1, dtype=np.float32)
- c = np.int32(np.random.random_sample())
-
- print 'testing det...',
-
- det(cuda.In(a), cuda.Out(dest), **size)
-
- if not np.allclose(np.float32(np.linalg.det(a.reshape(3,3))), dest[0]):
- print 'fail'
- print np.float32(np.linalg.det(a.reshape(3,3)))
- print dest[0]
- else:
- print 'success'
-
- print 'testing inv...',
-
- dest = np.empty(9, dtype=np.float32)
-
- inv(cuda.In(a), cuda.Out(dest), **size)
-
- if not np.allclose(np.linalg.inv(a.reshape(3,3)).flatten().astype(np.float32), dest):
- print 'fail'
- print np.linalg.inv(a.reshape(3,3)).flatten().astype(np.float32)
- print dest
- else:
- print 'success'
-
- print 'testing matrixadd...',
-
- matrixadd(cuda.In(a), cuda.In(b), cuda.Out(dest), **size)
-
- if not np.allclose(a+b, dest):
- print 'fail'
- print a+b
- print dest
- else:
- print 'success'
-
- print 'testing matrixsub...',
-
- matrixsub(cuda.In(a), cuda.In(b), cuda.Out(dest), **size)
-
- if not np.allclose(a-b, dest):
- print 'fail'
- print a-b
- print dest
- else:
- print 'sucess'
-
- print 'testing matrixmul...',
-
- matrixmul(cuda.In(a), cuda.In(b), cuda.Out(dest), **size)
-
- if not np.allclose(np.dot(a.reshape(3,3),b.reshape(3,3)).flatten(), dest):
- print 'fail'
- print np.dot(a.reshape(3,3),b.reshape(3,3)).flatten()
- print dest
- else:
- print 'success'
-
- print 'testing multiply...',
-
- x_cpu = np.random.random_sample(size=3).astype(np.float32)
- x_gpu = np.array((x_cpu[0], x_cpu[1], x_cpu[2]), dtype=float3)
-
- dest = np.empty(1, dtype=float3)
-
- multiply(cuda.In(a), cuda.In(x_gpu), cuda.Out(dest), **size)
-
- m = a.reshape(3,3)
-
- if not np.allclose(np.dot(x_cpu,m[0]), dest[0]['x']) or \
- not np.allclose(np.dot(x_cpu,m[1]), dest[0]['y']) or \
- not np.allclose(np.dot(x_cpu,m[2]), dest[0]['z']):
- print 'fail'
- print np.dot(x_cpu,m[0])
- print np.dot(x_cpu,m[1])
- print np.dot(x_cpu,m[2])
- print dest[0]['x']
- print dest[0]['y']
- print dest[0]['z']
- else:
- print 'success'
-
- print 'testing matrixaddfloat...',
-
- dest = np.empty(9, dtype=np.float32)
-
- matrixaddfloat(cuda.In(a), c, cuda.Out(dest), **size)
-
- if not np.allclose(a+c, dest):
- print 'fail'
- print a+c
- print dest
- else:
- print 'success'
-
- print 'testing matrixsubfloat...',
-
- matrixsubfloat(cuda.In(a), c, cuda.Out(dest), **size)
-
- if not np.allclose(a-c, dest):
- print 'fail'
- print a-c
- print dest
- else:
- print 'success'
-
- print 'testing matrixmulfloat...',
-
- matrixmulfloat(cuda.In(a), c, cuda.Out(dest), **size)
-
- if not np.allclose(a*c, dest):
- print 'fail'
- print a-c
- print dest
- else:
- print 'success'
-
- print 'testing matrixdivfloat...',
-
- matrixdivfloat(cuda.In(a), c, cuda.Out(dest), **size)
-
- if not np.allclose(a/c, dest):
- print 'fail'
- print a/c
- print dest
- else:
- print 'success'
-
- print 'testing floataddmatrix...',
-
- floataddmatrix(cuda.In(a), c, cuda.Out(dest), **size)
-
- if not np.allclose(c+a, dest):
- print 'fail'
- print c+a
- print dest
- else:
- print 'success'
-
- print 'testing floatsubmatrix...',
-
- floatsubmatrix(cuda.In(a), c, cuda.Out(dest), **size)
-
- if not np.allclose(c-a, dest):
- print 'fail'
- print c-a
- print dest
- else:
- print 'success'
-
- print 'testing floatmulmatrix...',
-
- floatmulmatrix(cuda.In(a), c, cuda.Out(dest), **size)
-
- if not np.allclose(c*a, dest):
- print 'fail'
- print c*a
- print dest
- else:
- print 'success'
-
- print 'testing floatdivmatrix...',
-
- floatdivmatrix(cuda.In(a), c, cuda.Out(dest), **size)
-
- if not np.allclose(c/a, dest):
- print 'fail'
- print c/a
- print dest
- else:
- print 'success'
-
- print 'testing matrixaddequals...',
-
- dest = np.copy(a)
-
- matrixaddequals(cuda.InOut(dest), cuda.In(b), **size)
-
- if not np.allclose(a+b, dest):
- print 'fail'
- print a+b
- print dest
- else:
- print 'success'
-
- print 'testing matrixsubequals...',
-
- dest = np.copy(a)
-
- matrixsubequals(cuda.InOut(dest), cuda.In(b), **size)
-
- if not np.allclose(a-b, dest):
- print 'fail'
- print a-b
- print dest
- else:
- print 'success'
-
- print 'testing matrixaddequalsfloat...',
-
- dest = np.copy(a)
-
- matrixaddequalsfloat(cuda.InOut(dest), c, **size)
-
- if not np.allclose(a+c, dest):
- print 'fail'
- print a+c
- print dest
- else:
- print 'success'
-
- print 'testing matrixsubequalsfloat...',
-
- dest = np.copy(a)
-
- matrixsubequalsfloat(cuda.InOut(dest), c, **size)
-
- if not np.allclose(a-c, dest):
- print 'fail'
- print a-c
- print dest
- else:
- print 'success'
-
- print 'testing matrixmulequalsfloat...',
-
- dest = np.copy(a)
-
- matrixmulequalsfloat(cuda.InOut(dest), c, **size)
-
- if not np.allclose(a*c, dest):
- print 'fail'
- print a*c
- print dest
- else:
- print 'success'
-
- print 'testing matrixdivequalsfloat...',
-
- dest = np.copy(a)
-
- matrixdivequalsfloat(cuda.InOut(dest), c, **size)
-
- if not np.allclose(a/c, dest):
- print 'fail'
- print a*c
- print dest
- else:
- print 'success'
-
- print 'testing outer...',
-
- x1_cpu = np.random.random_sample(size=3).astype(np.float32)
- x2_cpu = np.random.random_sample(size=3).astype(np.float32)
-
- x1_gpu = np.array((x1_cpu[0], x1_cpu[1], x1_cpu[2]), dtype=float3)
- x2_gpu = np.array((x2_cpu[0], x2_cpu[1], x2_cpu[2]), dtype=float3)
-
- outer(x1_gpu, x2_gpu, cuda.Out(dest), **size)
-
- if not np.allclose(np.outer(x1_cpu, x2_cpu).flatten(), dest):
- print 'fail'
- print np.outer(x1_cpu, x2_cpu).flatten()
- print dest
- else:
- print 'success'
-
- print 'testing minus matrix...',
-
- dest = np.copy(a)
-
- minusmatrix(cuda.In(a), cuda.Out(dest), **size)
-
- if not np.allclose(-a, dest):
- print 'fail'
- print -a
- print dest
- else:
- print 'success'
diff --git a/tests/rotate_test.cu b/tests/rotate_test.cu
deleted file mode 100644
index 6cafc12..0000000
--- a/tests/rotate_test.cu
+++ /dev/null
@@ -1,14 +0,0 @@
-//-*-c-*-
-
-#include "rotate.h"
-
-extern "C"
-{
-
-__global__ void rotate(float3 *a, float *phi, float3 *n, float3 *dest)
-{
- int idx = blockIdx.x*blockDim.x + threadIdx.x;
- dest[idx] = rotate(a[idx], phi[idx], n[idx]);
-}
-
-} // extern "c"
diff --git a/tests/rotate_test.py b/tests/rotate_test.py
deleted file mode 100644
index 92eff84..0000000
--- a/tests/rotate_test.py
+++ /dev/null
@@ -1,67 +0,0 @@
-import os
-import numpy as np
-from pycuda import autoinit
-from pycuda.compiler import SourceModule
-import pycuda.driver as cuda
-from pycuda import gpuarray
-float3 = gpuarray.vec.float3
-
-def rotate(x, phi, n):
- x = np.asarray(x)
- n = np.asarray(n)
-
- r = np.cos(phi)*np.identity(3) + (1-np.cos(phi))*np.outer(n,n) + \
- np.sin(phi)*np.array([[0,n[2],-n[1]],[-n[2],0,n[0]],[n[1],-n[0],0]])
-
- return np.inner(x,r)
-
-print 'device %s' % autoinit.device.name()
-
-current_directory = os.path.split(os.path.realpath(__file__))[0]
-source_directory = current_directory + '/../src'
-
-source = open(current_directory + '/rotate_test.cu').read()
-
-mod = SourceModule(source, options=['-I' + source_directory], no_extern_c=True, cache_dir=False)
-
-rotate_gpu = mod.get_function('rotate')
-
-size = {'block': (100,1,1), 'grid': (1,1)}
-
-a = np.empty(size['block'][0], dtype=float3)
-n = np.empty(size['block'][0], dtype=float3)
-phi = np.random.random_sample(size=a.size).astype(np.float32)
-
-a['x'] = np.random.random_sample(size=a.size)
-a['y'] = np.random.random_sample(size=a.size)
-a['z'] = np.random.random_sample(size=a.size)
-
-n['x'] = np.random.random_sample(size=n.size)
-n['y'] = np.random.random_sample(size=n.size)
-n['z'] = np.random.random_sample(size=n.size)
-
-a['x'] = np.ones(a.size)
-a['y'] = np.zeros(a.size)
-a['z'] = np.zeros(a.size)
-
-n['x'] = np.zeros(n.size)
-n['y'] = np.zeros(n.size)
-n['z'] = np.ones(n.size)
-
-phi = np.array([np.pi/2]*a.size).astype(np.float32)
-
-def testrotate():
- dest = np.empty(a.size, dtype=float3)
- rotate_gpu(cuda.In(a), cuda.In(phi), cuda.In(n), cuda.Out(dest), **size)
- for v, theta, w, rdest in zip(a,phi,n,dest):
- r = rotate((v['x'], v['y'], v['z']), theta, (w['x'], w['y'], w['z']))
- if not np.allclose(rdest['x'], r[0]) or \
- not np.allclose(rdest['y'], r[1]) or \
- not np.allclose(rdest['z'], r[2]):
- print v
- print theta
- print w
- print r
- print rdest
- assert False
-
diff --git a/tests/test_fileio.py b/tests/test_fileio.py
deleted file mode 100644
index 3869a9f..0000000
--- a/tests/test_fileio.py
+++ /dev/null
@@ -1,74 +0,0 @@
-import unittest
-from chroma.fileio import root
-from chroma import event
-import numpy as np
-
-class TestFileIO(unittest.TestCase):
- def test_file_write_and_read(self):
- ev = event.Event(1, event.Vertex('e-', pos=(0,0,1), dir=(1,0,0),
- ke=15.0, pol=(0,1,0)))
-
- photons_beg = root.make_photon_with_arrays(1)
- photons_beg.pos[0] = (1,2,3)
- photons_beg.dir[0] = (4,5,6)
- photons_beg.pol[0] = (7,8,9)
- photons_beg.wavelengths[0] = 400.0
- photons_beg.t[0] = 100.0
- photons_beg.last_hit_triangles[0] = 5
- photons_beg.flags[0] = 20
- ev.photons_beg = photons_beg
-
- photons_end = root.make_photon_with_arrays(1)
- photons_end.pos[0] = (1,2,3)
- photons_end.dir[0] = (4,5,6)
- photons_end.pol[0] = (7,8,9)
- photons_end.wavelengths[0] = 400.0
- photons_end.t[0] = 100.0
- photons_end.last_hit_triangles[0] = 5
- photons_end.flags[0] = 20
- ev.photons_end = photons_end
-
- ev.vertices = [ev.primary_vertex]
-
- channels = event.Channels(hit=np.array([True, False]),
- t=np.array([20.0, 1e9], dtype=np.float32),
- q=np.array([2.0, 0.0], dtype=np.float32),
- flags=np.array([8, 32], dtype=np.uint32))
- ev.channels = channels
-
- filename = '/tmp/chroma-filewritertest.root'
- writer = root.RootWriter(filename)
- writer.write_event(ev)
- writer.close()
-
- # Exercise the RootReader methods
- reader = root.RootReader(filename)
- self.assertEquals(len(reader), 1)
-
- self.assertRaises(StopIteration, reader.prev)
-
- reader.next()
-
- self.assertEqual(reader.index(), 0)
- self.assertRaises(StopIteration, reader.next)
-
- reader.jump_to(0)
-
- # Enough screwing around, let's get the one event in the file
- newev = reader.current()
-
- # Now check if everything is correct in the event
- for attribute in ['id']:
- self.assertEqual(getattr(ev, attribute), getattr(newev, attribute), 'compare %s' % attribute)
-
- for attribute in ['pos', 'dir', 'pol', 'ke', 't0']:
- self.assertTrue(np.allclose(getattr(ev.primary_vertex, attribute), getattr(newev.primary_vertex, attribute)), 'compare %s' % attribute)
-
- for i in range(len(ev.vertices)):
- self.assertTrue(np.allclose(getattr(ev.vertices[i], attribute), getattr(newev.vertices[i], attribute)), 'compare %s' % attribute)
-
- for attribute in ['pos', 'dir', 'pol', 'wavelengths', 't', 'last_hit_triangles', 'flags']:
- self.assertTrue(np.allclose(getattr(ev.photons_beg, attribute),
- getattr(newev.photons_beg, attribute)), 'compare %s' % attribute)
- self.assertTrue(np.allclose(getattr(ev.photons_end, attribute),
- getattr(newev.photons_end, attribute)), 'compare %s' % attribute)
diff --git a/tests/test_generator_photon.py b/tests/test_generator_photon.py
deleted file mode 100644
index 13501fe..0000000
--- a/tests/test_generator_photon.py
+++ /dev/null
@@ -1,21 +0,0 @@
-import unittest
-import itertools
-
-from chroma import generator
-from chroma.generator.vertex import constant_particle_gun
-from chroma import optics
-
-class TestG4ParallelGenerator(unittest.TestCase):
- def test_center(self):
- '''Generate Cherenkov light at the center of the world volume'''
- gen = generator.photon.G4ParallelGenerator(1, optics.water_wcsim)
- vertex = itertools.islice(constant_particle_gun('e-', (0,0,0), (1,0,0), 100), 10)
- for event in gen.generate_events(vertex):
- self.assertGreater(len(event.photons_beg.pos), 0)
-
- def test_off_center(self):
- '''Generate Cherenkov light at (1 m, 0 m, 0 m)'''
- gen = generator.photon.G4ParallelGenerator(1, optics.water_wcsim)
- vertex = itertools.islice(constant_particle_gun('e-', (1,0,0), (1,0,0), 100), 10)
- for event in gen.generate_events(vertex):
- self.assertGreater(len(event.photons_beg.pos), 0)
diff --git a/tests/test_generator_vertex.py b/tests/test_generator_vertex.py
deleted file mode 100644
index cec363d..0000000
--- a/tests/test_generator_vertex.py
+++ /dev/null
@@ -1,23 +0,0 @@
-import unittest
-import itertools
-import numpy as np
-import chroma.generator.vertex
-
-class TestParticleGun(unittest.TestCase):
- def test_constant_particle_gun_center(self):
- '''Generate electron vertices at the center of the world volume.'''
- vertex = chroma.generator.vertex.constant_particle_gun('e-', (0,0,0), (1,0,0), 100)
- for ev in itertools.islice(vertex, 100):
- self.assertEquals(ev.primary_vertex.particle_name, 'e-')
- self.assertTrue(np.allclose(ev.primary_vertex.pos, [0,0,0]))
- self.assertTrue(np.allclose(ev.primary_vertex.dir, [1,0,0]))
- self.assertTrue(np.allclose(ev.primary_vertex.ke, 100))
-
- def test_off_center(self):
- '''Generate electron vertices at (1,0,0) in the world volume.'''
- vertex = chroma.generator.vertex.constant_particle_gun('e-', (1,0,0), (1,0,0), 100)
- for ev in itertools.islice(vertex, 100):
- self.assertEquals(ev.primary_vertex.particle_name, 'e-')
- self.assertTrue(np.allclose(ev.primary_vertex.pos, [1,0,0]))
- self.assertTrue(np.allclose(ev.primary_vertex.dir, [1,0,0]))
- self.assertTrue(np.allclose(ev.primary_vertex.ke, 100))
diff --git a/tests/test_pdf.py b/tests/test_pdf.py
deleted file mode 100644
index 571cbd4..0000000
--- a/tests/test_pdf.py
+++ /dev/null
@@ -1,67 +0,0 @@
-import unittest
-import numpy as np
-import itertools
-
-import chroma.detectors
-from chroma.generator.photon import G4ParallelGenerator
-from chroma.generator.vertex import constant_particle_gun
-from chroma.optics import water_wcsim
-from chroma import gpu
-from chroma.sim import Simulation
-
-class TestPDF(unittest.TestCase):
- def setUp(self):
- self.detector = chroma.detectors.microlbne()
- self.detector.build()
- self.vertex_gen = constant_particle_gun('e-', (0,0,0), (1,0,0), 10)
-
- def testGPUPDF(self):
- '''Create a hit count and (q,t) PDF for 10 MeV events in MicroLBNE'''
-
- g4generator = G4ParallelGenerator(1, water_wcsim)
-
- context = gpu.create_cuda_context()
-
- gpu_geometry = gpu.GPUGeometry(self.detector)
-
- nthreads_per_block, max_blocks = 64, 1024
-
- rng_states = gpu.get_rng_states(nthreads_per_block*max_blocks)
-
- gpu_daq = gpu.GPUDaq(gpu_geometry, max(self.detector.pmtids))
- gpu_pdf = gpu.GPUPDF()
- gpu_pdf.setup_pdf(max(self.detector.pmtids), 100, (-0.5, 999.5), 10, (-0.5, 9.5))
-
- gpu_pdf.clear_pdf()
-
- for ev in g4generator.generate_events(itertools.islice(self.vertex_gen, 10)):
- gpu_photons = gpu.GPUPhotons(ev.photons_beg)
- gpu_photons.propagate(gpu_geometry, rng_states, nthreads_per_block,
- max_blocks)
- gpu_channels = gpu_daq.acquire(gpu_photons, rng_states,
- nthreads_per_block, max_blocks)
- gpu_pdf.add_hits_to_pdf(gpu_channels)
-
- hitcount, pdf = gpu_pdf.get_pdfs()
- self.assertTrue( (hitcount > 0).any() )
- self.assertTrue( (pdf > 0).any() )
-
- # Consistency checks
- for i, nhits in enumerate(hitcount):
- self.assertEqual(nhits, pdf[i].sum())
-
- context.pop()
-
- def testSimPDF(self):
- sim = Simulation(self.detector)
-
- vertex_iter = itertools.islice(self.vertex_gen, 10)
-
- hitcount, pdf = sim.create_pdf(vertex_iter, 100, (-0.5, 999.5), 10, (-0.5, 9.5))
-
- self.assertTrue( (hitcount > 0).any() )
- self.assertTrue( (pdf > 0).any() )
-
- # Consistency checks
- for i, nhits in enumerate(hitcount):
- self.assertEqual(nhits, pdf[i].sum())
diff --git a/tests/test_propagation.py b/tests/test_propagation.py
deleted file mode 100644
index 43ecb9b..0000000
--- a/tests/test_propagation.py
+++ /dev/null
@@ -1,58 +0,0 @@
-import unittest
-import numpy as np
-
-from chroma.geometry import Solid, Geometry
-from chroma.make import box
-from chroma.sim import Simulation
-from chroma.optics import vacuum
-from chroma.event import Photons
-
-class TestPropagation(unittest.TestCase):
- def testAbort(self):
- '''Photons that hit a triangle at normal incidence should not abort.
-
- Photons that hit a triangle at exactly normal incidence can sometimes
- produce a dot product that is outside the range allowed by acos().
- Trigger these with axis aligned photons in a box.
- '''
-
- # Setup geometry
- cube = Geometry(vacuum)
- cube.add_solid(Solid(box(100,100,100), vacuum, vacuum))
- cube.pmtids = [0]
- cube.build(use_cache=False)
- sim = Simulation(cube, geant4_processes=0)
-
- # Create initial photons
- nphotons = 10000
- pos = np.tile([0,0,0], (nphotons,1)).astype(np.float32)
- dir = np.tile([0,0,1], (nphotons,1)).astype(np.float32)
- pol = np.zeros_like(pos)
- phi = np.random.uniform(0, 2*np.pi, nphotons).astype(np.float32)
- pol[:,0] = np.cos(phi)
- pol[:,1] = np.sin(phi)
- t = np.zeros(nphotons, dtype=np.float32)
- wavelengths = np.empty(nphotons, np.float32)
- wavelengths.fill(400.0)
-
- photons = Photons(pos=pos, dir=dir, pol=pol, t=t,
- wavelengths=wavelengths)
-
- # First make one step to check for strangeness
- photons_end = sim.simulate([photons], keep_photons_end=True,
- max_steps=1).next().photons_end
-
- self.assertFalse(np.isnan(photons_end.pos).any())
- self.assertFalse(np.isnan(photons_end.dir).any())
- self.assertFalse(np.isnan(photons_end.pol).any())
- self.assertFalse(np.isnan(photons_end.t).any())
- self.assertFalse(np.isnan(photons_end.wavelengths).any())
-
- # Now let it run the usual ten steps
- photons_end = sim.simulate([photons], keep_photons_end=True,
- max_steps=10).next().photons_end
- aborted = (photons_end.flags & (1 << 31)) > 0
- print 'aborted photons: %1.1f' % \
- (float(np.count_nonzero(aborted)) / nphotons)
- self.assertFalse(aborted.any())
-
diff --git a/tests/test_ray_intersection.py b/tests/test_ray_intersection.py
deleted file mode 100644
index 7d0c53c..0000000
--- a/tests/test_ray_intersection.py
+++ /dev/null
@@ -1,27 +0,0 @@
-import unittest
-import chroma
-import numpy as np
-import os
-from pycuda import gpuarray as ga
-
-class TestRayIntersection(unittest.TestCase):
- def setUp(self):
- self.context = chroma.gpu.create_cuda_context()
- self.module = chroma.gpu.get_cu_module('mesh.h')
- self.gpu_funcs = chroma.gpu.GPUFuncs(self.module)
- self.box = chroma.gpu.GPUGeometry(chroma.build(chroma.make.cube()))
-
- pos, dir = chroma.project.from_film()
- self.pos_gpu = ga.to_gpu(chroma.gpu.to_float3(pos))
- self.dir_gpu = ga.to_gpu(chroma.gpu.to_float3(dir))
-
- testdir = os.path.dirname(os.path.abspath(chroma.tests.__file__))
- self.dx_standard = np.load(os.path.join(testdir,
- 'data/ray_intersection.npz'))
- def test_intersection_distance(self):
- dx = ga.zeros(self.pos_gpu.size, dtype=np.float32)
- self.gpu_funcs.distance_to_mesh(np.int32(self.pos_gpu.size), self.pos_gpu, self.dir_gpu, self.box.gpudata, dx, block=(64,1,1), grid=(self.pos_gpu.size//64+1,1))
- self.assertTrue((dx.get() == self.dx_standard).all())
-
- def tearDown(self):
- self.context.pop()
diff --git a/tests/test_rayleigh.py b/tests/test_rayleigh.py
deleted file mode 100644
index 02ccb41..0000000
--- a/tests/test_rayleigh.py
+++ /dev/null
@@ -1,57 +0,0 @@
-import unittest
-import numpy as np
-
-from chroma.geometry import Solid, Geometry
-from chroma.make import box
-from chroma.sim import Simulation
-from chroma.optics import water_wcsim
-from chroma.event import Photons
-import histogram
-from histogram.root import rootify
-import ROOT
-ROOT.gROOT.SetBatch(1)
-
-class TestRayleigh(unittest.TestCase):
- def setUp(self):
- self.cube = Geometry(water_wcsim)
- self.cube.add_solid(Solid(box(100,100,100), water_wcsim, water_wcsim))
- self.cube.pmtids = [0]
- self.cube.build(use_cache=False)
- self.sim = Simulation(self.cube, geant4_processes=0)
-
- nphotons = 100000
- pos = np.tile([0,0,0], (nphotons,1)).astype(np.float32)
- dir = np.tile([0,0,1], (nphotons,1)).astype(np.float32)
- pol = np.zeros_like(pos)
- phi = np.random.uniform(0, 2*np.pi, nphotons).astype(np.float32)
- pol[:,0] = np.cos(phi)
- pol[:,1] = np.sin(phi)
- t = np.zeros(nphotons, dtype=np.float32)
- wavelengths = np.empty(nphotons, np.float32)
- wavelengths.fill(400.0)
-
- self.photons = Photons(pos=pos, dir=dir, pol=pol, t=t, wavelengths=wavelengths)
-
- def testAngularDistributionPolarized(self):
- # Fully polarized photons
- self.photons.pol[:] = [1.0, 0.0, 0.0]
-
- photons_end = self.sim.simulate([self.photons], keep_photons_end=True, max_steps=1).next().photons_end
- aborted = (photons_end.flags & (1 << 31)) > 0
- self.assertFalse(aborted.any())
-
- # Compute the dot product between initial and final dir
- rayleigh_scatters = (photons_end.flags & (1 << 4)) > 0
- cos_scatter = (self.photons.dir[rayleigh_scatters] * photons_end.dir[rayleigh_scatters]).sum(axis=1)
- theta_scatter = np.arccos(cos_scatter)
- h = histogram.Histogram(bins=100, range=(0, np.pi))
- h.fill(theta_scatter)
- h = rootify(h)
-
- # The functional form for polarized light should be
- # (1 + \cos^2 \theta)\sin \theta according to GEANT4 physics
- # reference manual.
- f = ROOT.TF1("pol_func", "[0]*(1+cos(x)**2)*sin(x)", 0, np.pi)
- h.Fit(f)
- self.assertGreater(f.GetProb(), 1e-3)
-
diff --git a/tests/test_sample_cdf.cu b/tests/test_sample_cdf.cu
deleted file mode 100644
index 1401772..0000000
--- a/tests/test_sample_cdf.cu
+++ /dev/null
@@ -1,16 +0,0 @@
-// -*-c++-*-
-#include "random.h"
-
-extern "C" {
-
-__global__ void test_sample_cdf(int offset, int ncdf,
- float *cdf_x, float *cdf_y, float *out)
-{
- int id = blockDim.x * blockIdx.x + threadIdx.x;
- curandState s;
- curand_init(0, id, offset, &s);
-
- out[id] = sample_cdf(&s, ncdf, cdf_x, cdf_y);
-}
-
-}
diff --git a/tests/test_sample_cdf.py b/tests/test_sample_cdf.py
deleted file mode 100644
index 2baa2ac..0000000
--- a/tests/test_sample_cdf.py
+++ /dev/null
@@ -1,64 +0,0 @@
-import pycuda.autoinit
-import pycuda.driver as cuda
-from pycuda.compiler import SourceModule
-from pycuda import gpuarray
-import numpy as np
-import ROOT
-import os
-
-current_directory = os.path.split(os.path.realpath(__file__))[0]
-source_directory = current_directory + '/../src'
-
-source = open(current_directory + '/test_sample_cdf.cu').read()
-
-mod = SourceModule(source, options=['-I' + source_directory], no_extern_c=True, cache_dir=False)
-
-test_sample_cdf = mod.get_function('test_sample_cdf')
-
-def compare_sampling(hist, reps=10):
- nbins = hist.GetNbinsX();
- xaxis = hist.GetXaxis()
- intg = hist.GetIntegral()
- cdf_y = np.empty(nbins+1, dtype=float)
- cdf_x = np.empty_like(cdf_y)
-
- cdf_x[0] = xaxis.GetBinLowEdge(1)
- cdf_y[0] = 0.0
- for i in xrange(1,len(cdf_x)):
- cdf_y[i] = intg[i]
- cdf_x[i] = xaxis.GetBinUpEdge(i)
-
- cdf_x_gpu = gpuarray.to_gpu(cdf_x.astype(np.float32))
- cdf_y_gpu = gpuarray.to_gpu(cdf_y.astype(np.float32))
- block =(128,1,1)
- grid = (128, 1)
- out_gpu = gpuarray.empty(shape=int(block[0]*grid[0]), dtype=np.float32)
-
- out_h = ROOT.TH1D('out_h', '', hist.GetNbinsX(),
- xaxis.GetXmin(),
- xaxis.GetXmax())
- out_h.SetLineColor(ROOT.kGreen)
-
- for i in xrange(reps):
- test_sample_cdf(np.int32(i),
- np.int32(len(cdf_x_gpu)),
- cdf_x_gpu, cdf_y_gpu, out_gpu, block=block, grid=grid)
- out = out_gpu.get()
- for v in out:
- out_h.Fill(v)
-
- prob = out_h.KolmogorovTest(hist)
- return prob, out_h
-
-def test_sampling():
- '''Verify that the CDF-based sampler on the GPU reproduces a binned
- Gaussian distribution'''
- f = ROOT.TF1('f_gaussian', 'gaus(0)', -5, 5)
- f.SetParameters(1.0/np.sqrt(np.pi * 2), 0.0, 1.0)
- gaussian = ROOT.TH1D('gaussian', '', 100, -5, 5)
- gaussian.Add(f)
-
- prob, out_h = compare_sampling(gaussian, reps=50)
-
- assert prob > 0.01
-