you can now specify number of devices and morton ordering bits to threadtest script

2 files changed, 33 insertions, 120 deletions

diff --git a/likelihood.py b/likelihood.py
deleted file mode 100644
index ac100b6..0000000
--- a/likelihood.py
+++ /dev/null
@@ -1,96 +0,0 @@
-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)
-propagate = module.get_function('propagate')
-
-lbne = LBNE()
-lbne.build(bits=8)
-texrefs = lbne.load(module)
-
-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(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(lbne.solids))
-    gpu_bincount = np.bincount(lbne.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(lbne.solids)))
-    for i in range(neval):
-        print '\revent: %i' % (i+1),
-        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(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(lbne.solid_index[triangles])
-        bincount[i][:gpu_bincount.size] = gpu_bincount
-    print
-
-    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()
-
-        probability = h.ueval(event_bincount[i])
-
-        if probability.nominal_value == 0.0:
-            probability = ufloat((0.5/h.nentries, 0.5/h.nentries))
-
-        log_likelihood += umath.log(probability)
-
-    return -log_likelihood
diff --git a/threadtest.py b/threadtest.py
index 6d4e228..11bdaa4 100644
--- a/threadtest.py
+++ b/threadtest.py
@@ -11,16 +11,6 @@ from histogram import Histogram
 jobs = Queue()
 output = Queue()
 
-lbne = LBNE()
-lbne.build(bits=8)
-
-cuda.init()
-
-gputhreads = []
-for i in range(cuda.Device.count()):
-    gputhreads.append(GPUThread(i, lbne, jobs, output))
-    gputhreads[-1].start()
-
 def generate_event(pos=(0,0,0), nphotons=1000):
     origins_float3 = np.tile(gpuarray.vec.make_float3(*pos), (nphotons,1))
     directions = uniform_sphere(nphotons)
@@ -67,18 +57,37 @@ def likelihood(event_bincount, pos=(0,0,0), nphotons=1000, neval=1000):
 
     return -log_likelihood
 
-def stop():
-    for gputhread in gputhreads:
-        gputhread.stop()
-
-    for gputhread in gputhreads:
-        gputhread.join()
-
 if __name__ == '__main__':
-    event_bincount = generate_event()
-
-    for z in np.linspace(-1.0, 1.0, 100):
-        log_likelihood = likelihood(event_bincount, (z,0,0))
-        print 'z = %f, likelihood = %s' % (z, log_likelihood)
-
-    stop()
+    import sys
+    import optparse
+    import time
+
+    parser = optparse.OptionParser('%prog')
+    parser.add_option('-b', type='int', dest='nbits', default=8)
+    parser.add_option('-j', type='int', dest='ndevices', default=1)
+    options, args = parser.parse_args()
+
+    lbne = LBNE()
+    lbne.build(bits=options.nbits)
+
+    cuda.init()
+
+    gputhreads = []
+    for i in range(options.ndevices):
+        gputhreads.append(GPUThread(i, lbne, jobs, output))
+        gputhreads[-1].start()
+
+    try:
+        event_bincount = generate_event()
+
+        for z in np.linspace(-1.0, 1.0, 100):
+            t0 = time.time()
+            log_likelihood = likelihood(event_bincount, (z,0,0))
+            elapsed = time.time() - t0
+            print 'z = %f, likelihood = %s, elapsed %f sec' % (z, log_likelihood, elapsed)
+    finally:
+        for gputhread in gputhreads:
+            gputhread.stop()
+
+        for gputhread in gputhreads:
+            gputhread.join()