diff options
Diffstat (limited to 'likelihood.py')
| -rw-r--r-- | likelihood.py | 49 |
1 files changed, 18 insertions, 31 deletions
diff --git a/likelihood.py b/likelihood.py index ae712d8..b31f8fa 100644 --- a/likelihood.py +++ b/likelihood.py @@ -1,32 +1,26 @@ import sys +import time 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') +lbne = LBNE() +lbne.build(bits=9) +texrefs = lbne.load(module) + nblocks = 64 def generate_event(z=0.0, nphotons=1000): @@ -43,14 +37,15 @@ def generate_event(z=0.0, nphotons=1000): 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)) + propagate(np.int32(nphotons), origins_gpu, directions_gpu, np.int32(lbne.node_map.size-1), np.int32(lbne.first_node), dest_gpu, block=(nblocks,1,1), grid=(nphotons//nblocks+1,1), texrefs=texrefs) + cuda.Context.synchronize() 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 = np.zeros(len(lbne.solids)) + gpu_bincount = np.bincount(lbne.solid_index[triangles]) event_bincount[:gpu_bincount.size] = gpu_bincount return event_bincount @@ -59,9 +54,9 @@ 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))) + bincount = np.zeros((neval, len(lbne.solids))) for i in range(neval): - print '\reval %i' % i, + print '\revent: %i' % (i+1), sys.stdout.flush() directions = uniform_sphere(nphotons) @@ -74,36 +69,28 @@ def likelihood(event_bincount, z=0.0, nphotons=1000, neval=1000): 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)) + propagate(np.int32(nphotons), origins_gpu, directions_gpu, np.int32(lbne.node_map.size-1), np.int32(lbne.first_node), dest_gpu, block=(nblocks,1,1), grid=(nphotons//nblocks+1,1)) + cuda.Context.synchronize() cuda.memcpy_dtoh(dest, dest_gpu) triangles = dest[(dest != -1)] - gpu_bincount = np.bincount(geometry.solid_index[triangles]) + gpu_bincount = np.bincount(lbne.solid_index[triangles]) bincount[i][:gpu_bincount.size] = gpu_bincount print - histograms = [] - for i in range(len(geometry.solids)): + log_likelihood = ufloat((0,0)) + for i in range(len(lbne.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]))) + probability = h.ueval(event_bincount[i]) if probability.nominal_value == 0.0: - return None + probability = ufloat((0.5/h.nentries, 0.5/h.nentries)) 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)) |