From 084dfd08b714faefaea77cb7dc04d2e93dc04b1d Mon Sep 17 00:00:00 2001
From: Stan Seibert <stan@mtrr.org>
Date: Fri, 16 Sep 2011 14:27:46 -0400
Subject: File reorganization to move toward standard python package layout

---
 likelihood.py | 131 ----------------------------------------------------------
 1 file changed, 131 deletions(-)
 delete mode 100644 likelihood.py

(limited to 'likelihood.py')

diff --git a/likelihood.py b/likelihood.py
deleted file mode 100644
index 0e378b1..0000000
--- a/likelihood.py
+++ /dev/null
@@ -1,131 +0,0 @@
-import numpy as np
-from histogram import HistogramDD
-from uncertainties import ufloat, umath, unumpy
-from math import sqrt
-from itertools import islice, izip, compress, repeat
-from chroma.tools import profile_if_possible
-
-class Likelihood(object):
-    "Class to evaluate likelihoods for detector events."
-    def __init__(self, sim, event=None, tbins=100, trange=(-0.5e-9, 999.5e-9), 
-                 qbins=10, qrange=(-0.5, 9.5), time_only=True):
-        """
-        Args:
-            - sim: chroma.sim.Simulation
-                The simulation object used to simulate events and build pdfs.
-            - event: chroma.event.Event, *optional*
-                The detector event being reconstructed. If None, you must call
-                set_event() before eval().
-            - tbins: int, *optional*
-                Number of time bins in PDF
-            - trange: tuple of 2 floats, *optional*
-                Min and max time to include in PDF
-            - qbins: int, *optional*
-                Number of charge bins in PDF
-            - qrange: tuple of 2 floats, *optional*
-                Min and max charge to include in PDF
-            - time_only: bool, *optional*
-                Only use time observable (instead of time and charge) in computing
-                likelihood.  Defaults to True.
-        """
-        self.sim = sim
-        self.tbins = tbins
-        self.trange = trange
-        self.qbins = qbins
-        self.qrange = qrange
-        self.time_only = time_only
-
-        if event is not None:
-            self.set_event(event)
-
-    def set_event(self, event):
-        "Set the detector event being reconstructed."
-        self.event = event
-
-    @profile_if_possible
-    def eval(self, vertex_generator, nevals, nreps=16, ndaq=50):
-        """
-        Return the negative log likelihood that the detector event set in the
-        constructor or by set_event() was the result of a particle generated
-        by `vertex_generator`.  If `nreps` set to > 1, each set of photon
-        vertices will be propagated `nreps` times.
-        """
-        ntotal = nevals * nreps * ndaq
-
-        vertex_generator = islice(vertex_generator, nevals)
-
-        hitcount, pdf_prob, pdf_prob_uncert = \
-            self.sim.eval_pdf(self.event.channels,
-                              vertex_generator,
-                              0.1e-9, self.trange, 
-                              1, self.qrange,
-                              nreps=nreps,
-                              ndaq=ndaq,
-                              time_only=self.time_only)
-        
-        # Normalize probabilities and put a floor to keep the log finite
-        hit_prob = hitcount.astype(np.float32) / ntotal
-        hit_prob = np.maximum(hit_prob, 0.5 / ntotal)
-
-        # Set all zero or nan probabilities to limiting PDF value
-        bad_value = (pdf_prob <= 0.0) | np.isnan(pdf_prob)
-        if self.time_only:
-            pdf_floor = 1.0 / (self.trange[1] - self.trange[0])
-        else:
-            pdf_floor = 1.0 / (self.trange[1] - self.trange[0]) / (self.qrange[1] - self.qrange[0])
-        pdf_prob[bad_value] = pdf_floor
-        pdf_prob_uncert[bad_value] = pdf_floor
-
-        print 'channels with no data:', (bad_value & self.event.channels.hit).astype(int).sum()
-
-        # NLL calculation: note that negation is at the end
-        # Start with the probabilties of hitting (or not) the channels
-        hit_channel_prob = np.log(hit_prob[self.event.channels.hit]).sum() + np.log(1.0-hit_prob[~self.event.channels.hit]).sum()
-        hit_channel_prob_uncert = ( (ntotal * hit_prob * (1.0 - hit_prob)) / hit_prob**2 ).sum()**0.5
-        log_likelihood = ufloat((hit_channel_prob, 0.0))
-
-        # Then include the probability densities of the observed
-        # charges and times.
-        hit_pdf_ufloat = unumpy.uarray((pdf_prob[self.event.channels.hit],
-                                        pdf_prob_uncert[self.event.channels.hit]))
-        log_likelihood += unumpy.log(hit_pdf_ufloat).sum()
-
-        return -log_likelihood
-
-if __name__ == '__main__':
-    from chroma import detectors
-    from chroma.sim import Simulation
-    from chroma.optics import water
-    from chroma.generator import constant_particle_gun
-    from chroma import tools
-    import time
-
-    tools.enable_debug_on_crash()
-
-    detector = detectors.find('lbne')
-    sim = Simulation(detector, seed=0)
-
-    event = sim.simulate(islice(constant_particle_gun('e-',(0,0,0),(1,0,0),100.0), 1)).next()
-
-    print 'nhit = %i' % np.count_nonzero(event.channels.hit)
-
-    likelihood = Likelihood(sim, event)
-
-    x = np.linspace(-10.0, 10.0, 100)
-    l = []
-
-    for pos in izip(x, repeat(0), repeat(0)):
-        t0 = time.time()
-        ev_vertex_iter = constant_particle_gun('e-',pos,(1,0,0),100.0)
-        l.append(likelihood.eval(ev_vertex_iter, 1000))
-        elapsed = time.time() - t0
-
-        print '(%.1f, %.1f, %.1f), %s (%1.1f sec)' % \
-            (pos[0], pos[1], pos[2], tools.ufloat_to_str(l[-1]), elapsed)
-
-    import matplotlib.pyplot as plt
-
-    plt.errorbar(x, [v.nominal_value for v in l], [v.std_dev() for v in l])
-    plt.xlabel('X Position (m)')
-    plt.ylabel('Negative Log Likelihood')
-    plt.show()
-- 
cgit