import numpy as np
from histogram import HistogramDD
from uncertainties import ufloat, umath, unumpy
from math import sqrt
from itertools import izip, compress
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=1):
        """
        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.
        """

        hitcount, pdf_prob, pdf_prob_uncert = self.sim.eval_pdf(self.event.channels,
                                                                nevals, vertex_generator,
                                                                2.0e-9, self.trange, 
                                                                1, self.qrange,
                                                                nreps=nreps,
                                                                time_only=self.time_only)
        
        # Normalize probabilities and put a floor to keep the log finite
        hit_prob = hitcount.astype(np.float32) / (nreps * nevals)
        hit_prob = np.maximum(hit_prob, 0.5 / (nreps*nevals))

        # 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 = ( (nevals * nreps * 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.tools import enable_debug_on_crash
    import time

    enable_debug_on_crash()

    detector = detectors.find('lbne')
    sim = Simulation(detector, water)

    event = sim.simulate(1, constant_particle_gun('e-',(0,0,0),(1,0,0),100.0)).next()

    print 'nhit = %i' % np.count_nonzero(event.channels.hit)

    likelihood = Likelihood(sim, event)

    for x in np.linspace(-10.0, 10.0, 2*10*5+1):
        start = time.time()
        l = likelihood.eval(constant_particle_gun('e-',(0,0,0),(1,0,0),100.0),100, nreps=200)
        print 'x = %5.1f, %s (%1.1f sec)' % (x, l, time.time() - start)