diff options
| -rw-r--r-- | gputhread.py | 63 | ||||
| -rw-r--r-- | threadtest.py | 84 | ||||
| -rwxr-xr-x | view.py | 1 |
3 files changed, 148 insertions, 0 deletions
diff --git a/gputhread.py b/gputhread.py new file mode 100644 index 0000000..fe05e4f --- /dev/null +++ b/gputhread.py @@ -0,0 +1,63 @@ +import numpy as np +import pycuda.driver as cuda +from pycuda.compiler import SourceModule +import threading +import layout +from Queue import Empty + +class Job(object): + def __init__(self, origins, directions): + self.origins, self.directions = origins, directions + +class GPUThread(threading.Thread): + def __init__(self, device_id, geometry, jobs, output): + threading.Thread.__init__(self) + + self.device_id = device_id + self.geometry = geometry + self.jobs = jobs + self.output = output + self._stop = threading.Event() + + def stop(self): + self._stop.set() + + def stopped(self): + return self._stop.is_set() + + def run(self): + device = cuda.Device(self.device_id) + context = device.make_context() + 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') + texrefs = self.geometry.load(module) + + while not self.stopped(): + try: + job = self.jobs.get(timeout=2) + except Empty: + continue + + origins_gpu, directions_gpu = cuda.to_device(job.origins), \ + cuda.to_device(job.directions) + + dest = np.empty(job.origins.size, dtype=np.int32) + dest_gpu = cuda.to_device(dest) + + propagate(np.int32(job.origins.size), origins_gpu, directions_gpu, np.int32(self.geometry.node_map.size-1), np.int32(self.geometry.first_node), dest_gpu, block=(64,1,1), grid=(job.origins.size//64+1,1), texrefs=texrefs) + cuda.Context.synchronize() + + cuda.memcpy_dtoh(dest, dest_gpu) + + triangles = dest[(dest != -1)] + + bincount = np.zeros(len(self.geometry.solids)) + gpu_bincount = np.bincount(self.geometry.solid_index[triangles]) + bincount[:gpu_bincount.size] = gpu_bincount + + self.output.put(bincount) + self.jobs.task_done() + + context.pop() diff --git a/threadtest.py b/threadtest.py new file mode 100644 index 0000000..6d4e228 --- /dev/null +++ b/threadtest.py @@ -0,0 +1,84 @@ +from gputhread import * +from Queue import Queue +from detectors import LBNE +from photon import uniform_sphere +import numpy as np +from pycuda import gpuarray +import pycuda.driver as cuda +from uncertainties import ufloat, umath +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) + directions_float3 = np.empty(nphotons, dtype=gpuarray.vec.float3) + directions_float3['x'] = directions[:,0] + directions_float3['y'] = directions[:,1] + directions_float3['z'] = directions[:,2] + + jobs.put(Job(origins_float3, directions_float3)) + + jobs.join() + + return output.get() + +def likelihood(event_bincount, pos=(0,0,0), nphotons=1000, neval=1000): + for i in range(neval): + origins_float3 = np.tile(gpuarray.vec.make_float3(*pos), (nphotons,1)) + 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] + + jobs.put(Job(origins_float3, directions_float3)) + + jobs.join() + + bincount = np.zeros((neval, len(lbne.solids))) + for i in range(neval): + bincount[i] = output.get() + + 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 + +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() @@ -82,6 +82,7 @@ def view(geometry, name=''): def render(): """Render the mesh and display to screen.""" cuda_raytrace(np.int32(pixels.size), origins_gpu, directions_gpu, np.int32(geometry.node_map.size-1), np.int32(geometry.first_node), pixels_gpu, texrefs=texrefs, **gpu_kwargs) + cuda.Context.synchronize() cuda.memcpy_dtoh(pixels, pixels_gpu) pygame.surfarray.blit_array(screen, pixels.reshape(size)) |