A new PDF evaluation method that does not require storage proportional

to [number of bins] * [number of PMTs].  Instead we accumulate information
as the Monte Carlo runs in order to evaluate the PDFs only at the points
required for the likelihood calculation.

This new interface has been propagated all the way up from the GPU class
through the Simulation class to the Likelihood class.  We have preserved
the full binned histogram implementation in case we need it in the future.

1 files changed, 123 insertions, 1 deletions

diff --git a/gpu.py b/gpu.py
index 4d2197b..ab71358 100644
--- a/gpu.py
+++ b/gpu.py
@@ -115,7 +115,8 @@ class GPU(object):
         self.daq_module = SourceModule(chroma.src.daq, options=cuda_options, no_extern_c=True)
         self.daq_funcs = CUDAFuncs(self.daq_module, 
                               ['reset_earliest_time_int', 'run_daq',
-                               'convert_sortable_int_to_float', 'bin_hits'])
+                               'convert_sortable_int_to_float', 'bin_hits',
+                               'accumulate_pdf_eval'])
 
     def print_device_usage(self):
         print 'device usage:'
@@ -422,5 +423,126 @@ class GPU(object):
         '''Returns the 1D hitcount array and the 3D [channel, time, charge] histogram'''
         return self.hitcount_gpu.get(), self.pdf_gpu.get()
 
+    def setup_pdf_eval(self, event_hit, event_time, event_charge,
+                       min_twidth, trange, min_qwidth, qrange, min_bin_content=10,
+                       time_only=True):
+        '''Setup GPU arrays to compute PDF values for the given event.
+        The pdf_eval calculation allows the PDF to be evaluated at a
+        single point for each channel as the Monte Carlo is run.  The
+        effective bin size will be as small as (`min_twidth`,
+        `min_qwidth`) around the point of interest, but will be large
+        enough to ensure that `min_bin_content` Monte Carlo events
+        fall into the bin.
+
+            event_hit: ndarray
+              Hit or not-hit status for each channel in the detector.
+            event_time: ndarray
+              Hit time for each channel in the detector.  If channel 
+              not hit, the time will be ignored.
+            event_charge: ndarray
+              Integrated charge for each channel in the detector.
+              If channel not hit, the charge will be ignored.
+
+            min_twidth: float
+              Minimum bin size in the time dimension
+            trange: (float, float)
+              Range of time dimension in PDF
+            min_qwidth: float
+              Minimum bin size in charge dimension
+            qrange: (float, float)
+              Range of charge dimension in PDF
+            min_bin_content: int
+              The bin will be expanded to include at least this many events
+            time_only: bool
+              If True, only the time observable will be used in the PDF.
+        '''
+        self.event_hit_gpu = gpuarray.to_gpu(event_hit.astype(np.uint32))
+        self.event_time_gpu = gpuarray.to_gpu(event_time.astype(np.float32))
+        self.event_charge_gpu = gpuarray.to_gpu(event_charge.astype(np.float32))
+
+        self.eval_hitcount_gpu = gpuarray.zeros(len(event_hit), dtype=np.uint32)
+        self.eval_bincount_gpu = gpuarray.zeros(len(event_hit), dtype=np.uint32)
+        self.nearest_mc_gpu = gpuarray.empty(shape=len(event_hit) * min_bin_content, 
+                                             dtype=np.float32)
+        self.nearest_mc_gpu.fill(1e9)
+
+        self.min_twidth = min_twidth
+        self.trange = trange
+        self.min_qwidth = min_qwidth
+        self.qrange = qrange
+        self.min_bin_content = min_bin_content
+        self.time_only = time_only
+
+    def clear_pdf_eval(self):
+        '''Reset PDF evaluation counters to start accumulating new Monte Carlo.'''
+        self.eval_hitcount_gpu.fill(0)
+        self.eval_bincount_gpu.fill(0)
+        self.nearest_mc_gpu.fill(1e9)
+
+    def accumulate_pdf_eval(self):
+        '''Add the most recent results of run_daq() to the PDF evaluation.'''
+        
+        self.daq_funcs.accumulate_pdf_eval(np.int32(self.time_only),
+                                           np.int32(len(self.event_hit_gpu)),
+                                           self.event_hit_gpu,
+                                           self.event_time_gpu,
+                                           self.event_charge_gpu,
+                                           self.earliest_time_gpu,
+                                           self.channel_q_gpu,
+                                           self.eval_hitcount_gpu,
+                                           self.eval_bincount_gpu,
+                                           np.float32(self.min_twidth),
+                                           np.float32(self.trange[0]),
+                                           np.float32(self.trange[1]),
+                                           np.float32(self.min_qwidth),
+                                           np.float32(self.qrange[0]),
+                                           np.float32(self.qrange[1]),
+                                           self.nearest_mc_gpu,
+                                           np.int32(self.min_bin_content),
+                                           block=(self.nthreads_per_block,1,1), 
+                                           grid=(len(self.earliest_time_gpu)//self.nthreads_per_block+1,1))
+
+    def get_pdf_eval(self):
+        evhit = self.event_hit_gpu.get().astype(bool)
+        hitcount = self.eval_hitcount_gpu.get()
+        bincount = self.eval_bincount_gpu.get()
+
+        pdf_value = np.zeros(len(hitcount), dtype=float)
+        pdf_frac_uncert = np.zeros_like(pdf_value)
+
+        # PDF value for high stats bins
+        high_stats = bincount >= self.min_bin_content
+        if high_stats.any():
+            if self.time_only:
+                pdf_value[high_stats] = bincount[high_stats].astype(float) / hitcount[high_stats] / self.min_twidth
+            else:
+                assert Exception('Unimplemented 2D (time,charge) mode!')
+
+            pdf_frac_uncert[high_stats] = 1.0/np.sqrt(bincount[high_stats])
+
+        # PDF value for low stats bins
+        low_stats = ~high_stats & (hitcount > 0) & evhit
+
+        nearest_mc = self.nearest_mc_gpu.get().reshape((len(hitcount), self.min_bin_content))
+
+        # Deal with the case where we did not even get min_bin_content events in the PDF
+        # but also clamp the lower range to ensure we don't index by a negative number in 2 lines
+        last_valid_entry = np.maximum(0, (nearest_mc < 1e9).astype(int).sum(axis=1) - 1)
+        distance = nearest_mc[np.arange(len(last_valid_entry)),last_valid_entry]
+
+        if low_stats.any():
+            if self.time_only:
+                pdf_value[low_stats] = (last_valid_entry[low_stats] + 1).astype(float) / hitcount[low_stats] / distance[low_stats] / 2.0
+            else:
+                assert Exception('Unimplemented 2D (time,charge) mode!')
+
+            pdf_frac_uncert[low_stats] = 1.0/np.sqrt(last_valid_entry[low_stats] + 1)
+
+        # PDFs with no stats got zero by default during array creation
+        
+        print 'high_stats:', high_stats.sum(), 'low_stats', low_stats.sum()
+
+        return hitcount, pdf_value, pdf_value * pdf_frac_uncert
+
     def __del__(self):
         self.context.pop()