1 files changed, 109 insertions, 0 deletions
diff --git a/likelihood.py b/likelihood.py
new file mode 100644
index 0000000..ae712d8
--- /dev/null
+++ b/likelihood.py
@@ -0,0 +1,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))