1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
|
import sys
import numpy as np
from pycuda import autoinit
from pycuda.compiler import SourceModule
from pycuda import gpuarray
import pycuda.driver as cuda
from uncertainties import ufloat, umath
from detectors import LBNE
from histogram import Histogram
from photon import uniform_sphere
import layout
print 'device %s' % autoinit.device.name()
source = open(layout.source + '/kernel.cu').read()
module = SourceModule(source, options=['-I' + layout.source], no_extern_c=True, cache_dir=False)
geometry = LBNE()
geometry.build(bits=6)
texrefs = geometry.load(module)
propagate = module.get_function('propagate')
nblocks = 64
def generate_event(z=0.0, nphotons=1000):
origins = np.tile(gpuarray.vec.make_float3(0,0,z), (nphotons,1))
origins_gpu = cuda.to_device(origins)
directions = uniform_sphere(nphotons)
directions_float3 = np.empty(nphotons, dtype=gpuarray.vec.float3)
directions_float3['x'] = directions[:,0]
directions_float3['y'] = directions[:,1]
directions_float3['z'] = directions[:,2]
directions_gpu = cuda.to_device(directions_float3)
dest = np.empty(nphotons, dtype=np.int32)
dest_gpu = cuda.to_device(dest)
propagate(np.int32(nphotons), origins_gpu, directions_gpu, np.int32(geometry.first_leaf), dest_gpu, block=(nblocks,1,1), grid=(nphotons//nblocks+1,1))
cuda.memcpy_dtoh(dest, dest_gpu)
triangles = dest[(dest != -1)]
event_bincount = np.zeros(len(geometry.solids))
gpu_bincount = np.bincount(geometry.solid_index[triangles])
event_bincount[:gpu_bincount.size] = gpu_bincount
return event_bincount
def likelihood(event_bincount, z=0.0, nphotons=1000, neval=1000):
origins = np.tile(gpuarray.vec.make_float3(0,0,z), (nphotons,1))
origins_gpu = cuda.to_device(origins)
bincount = np.zeros((neval, len(geometry.solids)))
for i in range(neval):
print '\reval %i' % i,
sys.stdout.flush()
directions = uniform_sphere(nphotons)
directions_float3 = np.empty(nphotons, dtype=gpuarray.vec.float3)
directions_float3['x'] = directions[:,0]
directions_float3['y'] = directions[:,1]
directions_float3['z'] = directions[:,2]
directions_gpu = cuda.to_device(directions_float3)
dest = np.empty(nphotons, dtype=np.int32)
dest_gpu = cuda.to_device(dest)
propagate(np.int32(nphotons), origins_gpu, directions_gpu, np.int32(geometry.first_leaf), dest_gpu, block=(nblocks,1,1), grid=(nphotons//nblocks+1,1))
cuda.memcpy_dtoh(dest, dest_gpu)
triangles = dest[(dest != -1)]
gpu_bincount = np.bincount(geometry.solid_index[triangles])
bincount[i][:gpu_bincount.size] = gpu_bincount
print
histograms = []
for i in range(len(geometry.solids)):
h = Histogram(100, (-0.5, 99.5))
h.fill(bincount[:,i])
h.normalize()
histograms.append(h)
log_likelihood = ufloat((0,0))
for i, h in enumerate(histograms):
probability = ufloat((h.eval(event_bincount[i]), h.eval_err(event_bincount[i])))
if probability.nominal_value == 0.0:
return None
log_likelihood += umath.log(probability)
return -log_likelihood
if __name__ == '__main__':
event_bincount = generate_event()
for z in np.arange(-2.5,2.5,0.1):
print 'z, likelihood = (%f, %s)' % (z, likelihood(event_bincount,z))
|
