From d93200c5cccc94c6a5bb544c927c766c52b5d1d1 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Tue, 28 Jul 2020 16:15:03 +0200
Subject: [PATCH 01/17] Simplify pixel log-likelihood.

* simplify the log-likelihood analytically
* discard constants and factors which are not needed in minimization
---
 ctapipe/image/pixel_likelihood.py | 68 ++++++++++++++++++++-----------
 1 file changed, 44 insertions(+), 24 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index 3a1cd76b8dc..2cc6c44e7fd 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -51,42 +51,62 @@ class PixelLikelihoodError(RuntimeError):
     pass
 
 
-def poisson_likelihood_gaussian(image, prediction, spe_width, ped):
-    """
-    Calculate likelihood of prediction given the measured signal, gaussian approx from
-    de Naurois et al 2009
+def poisson_likelihood_gaussian(image, prediction, spe_width, pedestal):
+    """Calculate negative log likelihood for every pixel.
+
+    Gaussian approximation from de Naurois, p. 22 (between (24) and (25)).
+
+    Simplification:
+
+    .. math::
+
+        θ = σ_p^2 + μ · (1 + σ_γ^2)
+
+        → P = \\frac{1}{\\sqrt{2 π θ}} · \\exp\\left(- \\frac{(s - μ)^2}{2 θ}\\right)
+
+        \\ln{P} = \\ln{\\frac{1}{\\sqrt{2 π θ}}} - \\frac{(s - μ)^2}{2 θ}
+
+                = \\ln{1} - \\ln{\\sqrt{2 π θ}} - \\frac{(s - μ)^2}{2 θ}
+
+                = - \\frac{\\ln{2 π θ}}{2} - \\frac{(s - μ)^2}{2 θ}
+
+                = - \\frac{\\ln{2 π} + \\ln{θ}}{2} - \\frac{(s - μ)^2}{2 θ}
+
+        - \\ln{P} = \\frac{\\ln{2 π} + \\ln{θ}}{2} + \\frac{(s - μ)^2}{2 θ}
+
+    and since we can remove constants and factors in the minimization:
+
+    .. math::
+
+        - \\ln{P} = \\ln{θ} + \\frac{(s - μ)^2}{θ}
+
+
+    References
+    ----------
+    - de Naurois, Rolland
+      'A high performance likelihood reconstruction of gamma-rays
+      for Imaging Atmospheric Cherenkov Telescopes' (2009)
 
     Parameters
     ----------
     image: ndarray
-        Pixel amplitudes from image
+        Pixel amplitudes from image (:math:`s`).
     prediction: ndarray
-        Predicted pixel amplitudes from model
+        Predicted pixel amplitudes from model (:math:`μ`).
     spe_width: ndarray
-        width of single p.e. distributio
-    ped: ndarray
-        width of pedestal
+        Width of single p.e. peak (:math:`σ_γ`).
+    pedestal: ndarray
+        Width of pedestal (:math:`σ_p`).
 
     Returns
     -------
-    ndarray: likelihood for each pixel
+    ndarray: Negative log-likelihood for each pixel.
     """
-    image = np.asarray(image)
-    prediction = np.asarray(prediction)
-    spe_width = np.asarray(spe_width)
-    ped = np.asarray(ped)
+    theta = pedestal ** 2 + prediction * (1 + spe_width ** 2)
 
-    sq = 1.0 / np.sqrt(2 * np.pi * (ped ** 2 + prediction * (1 + spe_width ** 2)))
-
-    diff = (image - prediction) ** 2
-    denom = 2 * (ped ** 2 + prediction * (1 + spe_width ** 2))
-    expo = np.asarray(np.exp(-1 * diff / denom))
-
-    # If we are outside of the range of datatype, fix to lower bound
-    min_prob = np.finfo(expo.dtype).tiny
-    expo[expo < min_prob] = min_prob
+    neg_log_l = np.log(theta) + (image - prediction) ** 2 / theta
 
-    return -2 * np.log(sq * expo)
+    return neg_log_l
 
 
 def poisson_likelihood_full(

From 39c2a294e7954b93c6a3d3493b9c4c89344a9da9 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Thu, 6 Aug 2020 14:16:41 +0200
Subject: [PATCH 02/17] Use reference from bibliography.rst

Reference was already there, so use correct link
and remove reference section in this docstring.
---
 ctapipe/image/pixel_likelihood.py | 8 +-------
 1 file changed, 1 insertion(+), 7 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index 2cc6c44e7fd..ea0c7727e3e 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -54,7 +54,7 @@ class PixelLikelihoodError(RuntimeError):
 def poisson_likelihood_gaussian(image, prediction, spe_width, pedestal):
     """Calculate negative log likelihood for every pixel.
 
-    Gaussian approximation from de Naurois, p. 22 (between (24) and (25)).
+    Gaussian approximation from [denaurois2009]_, p. 22 (between (24) and (25)).
 
     Simplification:
 
@@ -81,12 +81,6 @@ def poisson_likelihood_gaussian(image, prediction, spe_width, pedestal):
         - \\ln{P} = \\ln{θ} + \\frac{(s - μ)^2}{θ}
 
 
-    References
-    ----------
-    - de Naurois, Rolland
-      'A high performance likelihood reconstruction of gamma-rays
-      for Imaging Atmospheric Cherenkov Telescopes' (2009)
-
     Parameters
     ----------
     image: ndarray

From b9ca5d239305a18bced7e32c17c0077643016993 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 11:47:33 +0200
Subject: [PATCH 03/17] Add arxiv link to paper.

---
 docs/bibliography.rst | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/docs/bibliography.rst b/docs/bibliography.rst
index f0c44479e71..fab91f9d291 100644
--- a/docs/bibliography.rst
+++ b/docs/bibliography.rst
@@ -8,7 +8,7 @@ References
 .. [denaurois2009] de Naurois, M., Roland, L. "A high performance
 	likelihood reconstruction of gamma-rays for imaging
 	atmospheric Cherenkov telescopes". Astroparticle Physics
-	32.5, (2009)
+	32.5, (2009), https://arxiv.org/abs/0907.2610
 
 .. [temme2016] Temme T.F., "On the hunt for photons: analysis of Crab Nebula
 		data obtainedby the first G-APD Cherenkov telescope" 2016

From febb0613c69f740216b6a64a683c5c236a1eb4d3 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 12:08:38 +0200
Subject: [PATCH 04/17] Improve numeric approximation

---
 ctapipe/image/pixel_likelihood.py | 95 +++++++------------------------
 1 file changed, 22 insertions(+), 73 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index ea0c7727e3e..7a07bb4f8f7 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -30,11 +30,9 @@
 - Additional terms may be useful to add to the likelihood
 """
 
-import math
-
 import numpy as np
 from scipy.integrate import quad
-from scipy.special import factorial
+from scipy.stats import poisson
 
 __all__ = [
     "poisson_likelihood_gaussian",
@@ -103,96 +101,47 @@ def poisson_likelihood_gaussian(image, prediction, spe_width, pedestal):
     return neg_log_l
 
 
-def poisson_likelihood_full(
-    image, prediction, spe_width, ped, width_fac=3, dtype=np.float32
-):
+def neg_log_likelihood_numeric(image, prediction, spe_width, pedestal, confidence=(0.001, 0.999)):
     """
     Calculate likelihood of prediction given the measured signal,
-    full numerical integration from de Naurois et al 2009.
-    The width factor included here defines  the range over
-    which photo electron contributions are summed, and is
-    defined as a multiple of the expected resolution of
-    the highest amplitude pixel. For most applications
-    the defult of 3 is sufficient.
+    full numerical integration from [denaurois2009]_.
 
     Parameters
     ----------
     image: ndarray
-        Pixel amplitudes from image
+        Pixel amplitudes from image (:math:`s`).
     prediction: ndarray
-        Predicted pixel amplitudes from model
+        Predicted pixel amplitudes from model (:math:`μ`).
     spe_width: ndarray
-        width of single p.e. distribution
-    ped: ndarray
-        width of pedestal
-    width_fac: float
-        Factor to determine range of summation on integral
-    dtype: datatype
-        Data type of output array
+        Width of single p.e. peak (:math:`σ_γ`).
+    pedestal: ndarray
+        Width of pedestal (:math:`σ_p`).
+    confidence: tuple(float, float), 0 < x < 1
+        Confidence interval of poisson integration.
 
     Returns
     -------
-    ndarray: likelihood for each pixel
+    float
     """
 
-    image = np.asarray(image, dtype=dtype)
-    prediction = np.asarray(prediction, dtype=dtype)
-    spe_width = np.asarray(spe_width, dtype=dtype)
-    ped = np.asarray(ped, dtype=dtype)
-
-    if image.shape != prediction.shape:
-        raise PixelLikelihoodError(
-            ("Image and prediction arrays" " have different dimensions"),
-            "Image shape: ",
-            image.shape[0],
-            "Prediction shape: ",
-            prediction.shape[0],
-        )
-    max_val = np.max(image)
-    width = ped * ped + max_val * spe_width * spe_width
-    width = np.sqrt(np.abs(width))  # take abs of width if negative
+    epsilon = np.finfo(np.float).eps
 
-    max_sum = max_val + width_fac * width
-    if max_sum < 5:
-        max_sum = 5
+    prediction = prediction + epsilon
 
-    pe_summed = np.arange(max_sum)  # Need to decide how range is determined
-    pe_factorial = factorial(pe_summed)
+    likelihood = epsilon
 
-    first_term = prediction ** pe_summed[:, np.newaxis] * np.exp(-1 * prediction)
-    first_term /= pe_factorial[:, np.newaxis] * np.sqrt(
-        math.pi * 2 * (ped * ped + pe_summed[:, np.newaxis] * spe_width * spe_width)
+    ns = np.arange(
+        *poisson(np.max(prediction)).ppf(confidence),
     )
 
-    # Throw error if we get NaN in likelihood
-    if np.any(np.isnan(first_term)):
-        raise PixelLikelihoodError(
-            "Likelihood returning NaN,"
-            "likely due to extremely high signal"
-            " deviation. Switch to poisson_likelihood_safe"
-            " implementation or"
-            " increase floating point precision"
-            " e.g. dtype=float64"
-        )
-
-    # Should not have any porblems here with NaN that have not bee seens
-    second_term = (image - pe_summed[:, np.newaxis]) * (
-        image - pe_summed[:, np.newaxis]
-    )
-    second_term_denom = 2 * (
-        ped * ped + spe_width * spe_width * pe_summed[:, np.newaxis]
-    )
-
-    second_term = second_term / second_term_denom
-    second_term = np.exp(-1 * second_term)
-
-    # If we are outside of the range of datatype, fix to lower bound
-    min_prob = np.finfo(second_term.dtype).tiny
-    second_term[second_term < min_prob] = min_prob
+    ns = ns[ns >= 0]
 
-    like = first_term * second_term
+    for n in ns:
+        theta = pedestal ** 2 + n * spe_width ** 2
+        _l = prediction ** n * np.exp(-prediction) / theta * np.exp(-(image - n) ** 2 / (2 * theta))
+        likelihood += _l
 
-    return -2 * np.log(np.sum(like, axis=0))
+    return - np.sum(np.log(likelihood))
 
 
 def poisson_likelihood(

From 22049eacbb9ded96dd25d02e7b2bbca6cfda533c Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 12:09:28 +0200
Subject: [PATCH 05/17] Improve safe/smart likelihood calculation.

* Due to likelihoods returning floats not arrays,
calculate numeric and approximation, then sum both
---
 ctapipe/image/pixel_likelihood.py | 81 +++++++++++--------------------
 1 file changed, 29 insertions(+), 52 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index 7a07bb4f8f7..e1788224dcb 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -144,79 +144,56 @@ def neg_log_likelihood_numeric(image, prediction, spe_width, pedestal, confidenc
     return - np.sum(np.log(likelihood))
 
 
-def poisson_likelihood(
+def neg_log_likelihood(
     image,
     prediction,
     spe_width,
-    ped,
-    pedestal_safety=1.5,
-    width_fac=3,
-    dtype=np.float32,
+    pedestal,
+    prediction_safety=20.0,
 ):
     """
     Safe implementation of the poissonian likelihood implementation,
     adaptively switches between the full solution and the gaussian
-    approx depending on the signal. Pedestal safety parameter 
-    determines cross over point between the two solutions,
-    based on the expected p.e. resolution of the image pixels.
-    Therefore the cross over point will change dependent on 
-    the single p.e. resolution and pedestal levels.
+    approx depending on the prediction. Prediction safety parameter
+    determines cross over point between the two solutions.
 
     Parameters
     ----------
     image: ndarray
-        Pixel amplitudes from image
+        Pixel amplitudes from image (:math:`s`).
     prediction: ndarray
-        Predicted pixel amplitudes from model
+        Predicted pixel amplitudes from model (:math:`μ`).
     spe_width: ndarray
-        width of single p.e. distribution
-    ped: ndarray
-        width of pedestal
-    pedestal_safety: float
-        Decision point to choose between poissonian likelihood 
-        and gaussian approximation (p.e. resolution)
-    width_fac: float
-        Factor to determine range of summation on integral
-    dtype: datatype
-        Data type of output array
+        Width of single p.e. peak (:math:`σ_γ`).
+    pedestal: ndarray
+        Width of pedestal (:math:`σ_p`).
+    prediction_safety: float
+        Decision point to choose between poissonian likelihood
+        and gaussian approximation.
 
     Returns
     -------
-    ndarray: pixel likelihoods
+    float
     """
-    # Convert everything to arrays to begin
-    image = np.asarray(image, dtype=dtype)
-    prediction = np.asarray(prediction, dtype=dtype)
-    spe_width = np.asarray(spe_width, dtype=dtype)
-    ped = np.asarray(ped, dtype=dtype)
 
-    # Calculate photoelectron resolution
+    approx_mask = prediction > prediction_safety
 
-    width = ped * ped + image * spe_width * spe_width
-    width = np.asarray(width)
-    width[width < 0] = 0  # Set width to 0 for negative pixel amplitudes
-    width = np.sqrt(width)
+    neg_log_l = 0
+    neg_log_l += neg_log_likelihood_approx(
+        image[approx_mask],
+        prediction[approx_mask],
+        spe_width,
+        pedestal,
+    )
 
-    like = np.zeros(image.shape)
-    # If larger than safety value use gaussian approx
-    poisson_pix = width <= pedestal_safety
-    gaus_pix = width > pedestal_safety
-
-    if np.any(poisson_pix):
-        like[poisson_pix] = poisson_likelihood_full(
-            image[poisson_pix],
-            prediction[poisson_pix],
-            spe_width,
-            ped,
-            width_fac,
-            dtype,
-        )
-    if np.any(gaus_pix):
-        like[gaus_pix] = poisson_likelihood_gaussian(
-            image[gaus_pix], prediction[gaus_pix], spe_width, ped
-        )
+    neg_log_l += neg_log_likelihood_numeric(
+        image[~approx_mask],
+        prediction[~approx_mask],
+        spe_width,
+        pedestal,
+    )
 
-    return like
+    return neg_log_l
 
 
 def mean_poisson_likelihood_gaussian(prediction, spe_width, ped):

From 2d0d4e6709919c7318a52db5ac120e7e730b0281 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 12:12:25 +0200
Subject: [PATCH 06/17] Simplify mean likelihood

---
 ctapipe/image/pixel_likelihood.py | 23 +++++++++--------------
 1 file changed, 9 insertions(+), 14 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index e1788224dcb..da1a9c9b2d9 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -196,9 +196,8 @@ def neg_log_likelihood(
     return neg_log_l
 
 
-def mean_poisson_likelihood_gaussian(prediction, spe_width, ped):
-    """
-    Calculation of the mean  likelihood for a give expectation
+def mean_poisson_likelihood_gaussian(prediction, spe_width, pedestal):
+    """Calculation of the mean likelihood for a give expectation
     value of pixel intensity in the gaussian approximation.
     This is useful in the calculation of the goodness of fit.
 
@@ -207,22 +206,18 @@ def mean_poisson_likelihood_gaussian(prediction, spe_width, ped):
     prediction: ndarray
         Predicted pixel amplitudes from model
     spe_width: ndarray
-        width of single p.e. distribution
-    ped: ndarray
-        width of pedestal
+        Width of single p.e. distribution
+    pedestal: ndarray
+        Width of pedestal
 
     Returns
     -------
-    ndarray: mean likelihood for give pixel expectation
+    float
     """
-    prediction = np.asarray(prediction)
-    spe_width = np.asarray(spe_width)
-    ped = np.asarray(ped)
-
-    mean_like = 1 + np.log(2 * math.pi)
-    mean_like += np.log(ped * ped + prediction * (1 + spe_width * spe_width))
+    theta = pedestal ** 2 + prediction * (1 + spe_width ** 2)
+    mean_log_likelihood = 1 + np.log(2 * np.pi) + np.log(theta)
 
-    return mean_like
+    return np.sum(mean_log_likelihood)
 
 
 def _integral_poisson_likelihood_full(s, prediction, spe_width, ped):

From b25ee8c07c881d34a7aece7f3ae1d3538776cb72 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 12:13:00 +0200
Subject: [PATCH 07/17] Use link to bibliography

---
 ctapipe/image/pixel_likelihood.py | 6 ++----
 1 file changed, 2 insertions(+), 4 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index da1a9c9b2d9..67a853ecd72 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -1,9 +1,7 @@
 """
 Class for calculation of likelihood of a pixel expectation, given the pixel amplitude,
-the level of noise in the pixel and the photoelectron resolution. This calculation is
-taken from:
-de Naurois & Rolland, Astroparticle Physics, Volume 32, Issue 5, p. 231-252 (2009)
-https://arxiv.org/abs/0907.2610
+the level of noise in the pixel and the photoelectron resolution.
+This calculation is taken from [denaurois2009]_.
 
 The likelihood is essentially a poissonian convolved with a gaussian, at low signal
 a full possonian approach must be adopted, which requires the sum of contibutions

From 1f1d2d8a3e6c68b5ca850db2d91eb9ef53d44fc1 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 12:13:19 +0200
Subject: [PATCH 08/17] Fix doc

* Removed the factorial
---
 ctapipe/image/pixel_likelihood.py | 5 ++---
 1 file changed, 2 insertions(+), 3 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index 67a853ecd72..d4f8567aa90 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -5,9 +5,8 @@
 
 The likelihood is essentially a poissonian convolved with a gaussian, at low signal
 a full possonian approach must be adopted, which requires the sum of contibutions
-over a number of potential contributing photoelectrons (which is slow and can fail
-at high signals due to the factorial which mst be calculated). At high signal this
-simplifies to a gaussian approximation.
+over a number of potential contributing photoelectrons (which is slow).
+At high signal this simplifies to a gaussian approximation.
 
 The full and gaussian approximations are implemented, in addition to a general purpose
 implementation, which tries to intellegently switch 

From d096fc88d015ed08a327118e65f8cfb193a5ad50 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 12:14:23 +0200
Subject: [PATCH 09/17] Return sum of pixel likelihoods

---
 ctapipe/image/pixel_likelihood.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index d4f8567aa90..04d93b220de 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -89,13 +89,13 @@ def poisson_likelihood_gaussian(image, prediction, spe_width, pedestal):
 
     Returns
     -------
-    ndarray: Negative log-likelihood for each pixel.
+    float
     """
     theta = pedestal ** 2 + prediction * (1 + spe_width ** 2)
 
     neg_log_l = np.log(theta) + (image - prediction) ** 2 / theta
 
-    return neg_log_l
+    return np.sum(neg_log_l)
 
 
 def neg_log_likelihood_numeric(image, prediction, spe_width, pedestal, confidence=(0.001, 0.999)):

From c8e0fe813e80913010745ff52e3fd8966698e132 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 12:14:49 +0200
Subject: [PATCH 10/17] Rename functions

---
 ctapipe/image/pixel_likelihood.py | 20 ++++++++++----------
 ctapipe/reco/ImPACT.py            |  4 ++--
 2 files changed, 12 insertions(+), 12 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index 04d93b220de..d5e34167863 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -12,13 +12,13 @@
 implementation, which tries to intellegently switch 
 between the two. Speed tests are below:
 
-poisson_likelihood_gaussian(image, prediction, spe, ped)
+neg_log_likelihood_approx(image, prediction, spe, ped)
 29.8 µs per loop
 
-poisson_likelihood_full(image, prediction, spe, ped)
+neg_log_likelihood_numeric(image, prediction, spe, ped)
 93.4 µs per loop
 
-poisson_likelihood(image, prediction, spe, ped)
+neg_log_likelihood(image, prediction, spe, ped)
 59.9 µs per loop
 
 TODO:
@@ -32,9 +32,9 @@
 from scipy.stats import poisson
 
 __all__ = [
-    "poisson_likelihood_gaussian",
-    "poisson_likelihood_full",
-    "poisson_likelihood",
+    "neg_log_likelihood_approx",
+    "neg_log_likelihood_numeric",
+    "neg_log_likelihood",
     "mean_poisson_likelihood_gaussian",
     "mean_poisson_likelihood_full",
     "PixelLikelihoodError",
@@ -46,10 +46,10 @@ class PixelLikelihoodError(RuntimeError):
     pass
 
 
-def poisson_likelihood_gaussian(image, prediction, spe_width, pedestal):
-    """Calculate negative log likelihood for every pixel.
+def neg_log_likelihood_approx(image, prediction, spe_width, pedestal):
+    """Calculate negative log likelihood for telescope.
 
-    Gaussian approximation from [denaurois2009]_, p. 22 (between (24) and (25)).
+    Gaussian approximation from [denaurois2009]_, p. 22 (equation between (24) and (25)).
 
     Simplification:
 
@@ -222,7 +222,7 @@ def _integral_poisson_likelihood_full(s, prediction, spe_width, ped):
     Wrapper function around likelihood calculation, used in numerical
     integration.
     """
-    like = poisson_likelihood(s, prediction, spe_width, ped)
+    like = neg_log_likelihood(s, prediction, spe_width, ped)
     return like * np.exp(-0.5 * like)
 
 
diff --git a/ctapipe/reco/ImPACT.py b/ctapipe/reco/ImPACT.py
index e49ba6aa9ac..384e0c4a65b 100644
--- a/ctapipe/reco/ImPACT.py
+++ b/ctapipe/reco/ImPACT.py
@@ -18,7 +18,7 @@
     GroundFrame,
     project_to_ground,
 )
-from ctapipe.image import poisson_likelihood_gaussian, mean_poisson_likelihood_gaussian
+from ctapipe.image import neg_log_likelihood, mean_poisson_likelihood_gaussian
 from ctapipe.instrument import get_atmosphere_profile_functions
 from ctapipe.containers import (
     ReconstructedShowerContainer,
@@ -484,7 +484,7 @@ def get_likelihood(
         prediction *= self.template_scale
 
         # Get likelihood that the prediction matched the camera image
-        like = poisson_likelihood_gaussian(self.image, prediction, self.spe, self.ped)
+        like = neg_log_likelihood(self.image, prediction, self.spe, self.ped)
         like[np.isnan(like)] = 1e9
         like *= np.invert(ma.getmask(self.image))
         like = ma.MaskedArray(like, mask=ma.getmask(self.image))

From 61f81fa6c19239aa2f7f5e439b873b215638c231 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <nbiederbeck@tutanota.com>
Date: Mon, 17 Aug 2020 12:39:43 +0200
Subject: [PATCH 11/17] Run black

---
 ctapipe/image/pixel_likelihood.py | 42 +++++++++++++++----------------
 1 file changed, 21 insertions(+), 21 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index d5e34167863..4975abf44c3 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -98,7 +98,9 @@ def neg_log_likelihood_approx(image, prediction, spe_width, pedestal):
     return np.sum(neg_log_l)
 
 
-def neg_log_likelihood_numeric(image, prediction, spe_width, pedestal, confidence=(0.001, 0.999)):
+def neg_log_likelihood_numeric(
+    image, prediction, spe_width, pedestal, confidence=(0.001, 0.999)
+):
     """
     Calculate likelihood of prediction given the measured signal,
     full numerical integration from [denaurois2009]_.
@@ -127,26 +129,25 @@ def neg_log_likelihood_numeric(image, prediction, spe_width, pedestal, confidenc
 
     likelihood = epsilon
 
-    ns = np.arange(
-        *poisson(np.max(prediction)).ppf(confidence),
-    )
+    ns = np.arange(*poisson(np.max(prediction)).ppf(confidence),)
 
     ns = ns[ns >= 0]
 
     for n in ns:
         theta = pedestal ** 2 + n * spe_width ** 2
-        _l = prediction ** n * np.exp(-prediction) / theta * np.exp(-(image - n) ** 2 / (2 * theta))
+        _l = (
+            prediction ** n
+            * np.exp(-prediction)
+            / theta
+            * np.exp(-((image - n) ** 2) / (2 * theta))
+        )
         likelihood += _l
 
-    return - np.sum(np.log(likelihood))
+    return -np.sum(np.log(likelihood))
 
 
 def neg_log_likelihood(
-    image,
-    prediction,
-    spe_width,
-    pedestal,
-    prediction_safety=20.0,
+    image, prediction, spe_width, pedestal, prediction_safety=20.0,
 ):
     """
     Safe implementation of the poissonian likelihood implementation,
@@ -177,17 +178,11 @@ def neg_log_likelihood(
 
     neg_log_l = 0
     neg_log_l += neg_log_likelihood_approx(
-        image[approx_mask],
-        prediction[approx_mask],
-        spe_width,
-        pedestal,
+        image[approx_mask], prediction[approx_mask], spe_width, pedestal,
     )
 
     neg_log_l += neg_log_likelihood_numeric(
-        image[~approx_mask],
-        prediction[~approx_mask],
-        spe_width,
-        pedestal,
+        image[~approx_mask], prediction[~approx_mask], spe_width, pedestal,
     )
 
     return neg_log_l
@@ -260,7 +255,10 @@ def mean_poisson_likelihood_full(prediction, spe_width, ped):
     width = np.sqrt(width)
 
     for p in range(len(prediction)):
-        int_range = (prediction[p] - 10 * width[p], prediction[p] + 10 * width[p])
+        int_range = (
+            prediction[p] - 10 * width[p],
+            prediction[p] + 10 * width[p],
+        )
         mean_like[p] = quad(
             _integral_poisson_likelihood_full,
             int_range[0],
@@ -298,7 +296,9 @@ def chi_squared(image, prediction, ped, error_factor=2.9):
     if image.shape is not prediction.shape:
         PixelLikelihoodError(
             "Image and prediction arrays have different dimensions Image "
-            "shape: {} Prediction shape: {}".format(image.shape, prediction.shape)
+            "shape: {} Prediction shape: {}".format(
+                image.shape, prediction.shape
+            )
         )
 
     chi_square = (image - prediction) * (image - prediction)

From ccc19a11af6a47208f95e3d03b741b11f24c0107 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <noah.biederbeck@tu-dortmund.de>
Date: Tue, 15 Sep 2020 08:27:32 +0200
Subject: [PATCH 12/17] Black fixes

---
 ctapipe/image/pixel_likelihood.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index 4975abf44c3..495a9d1101a 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -9,7 +9,7 @@
 At high signal this simplifies to a gaussian approximation.
 
 The full and gaussian approximations are implemented, in addition to a general purpose
-implementation, which tries to intellegently switch 
+implementation, which tries to intellegently switch
 between the two. Speed tests are below:
 
 neg_log_likelihood_approx(image, prediction, spe, ped)
@@ -129,7 +129,7 @@ def neg_log_likelihood_numeric(
 
     likelihood = epsilon
 
-    ns = np.arange(*poisson(np.max(prediction)).ppf(confidence),)
+    ns = np.arange(*poisson(np.max(prediction)).ppf(confidence))
 
     ns = ns[ns >= 0]
 

From 7b81f6492e88fb4f0e782c4db942b365fc0e28ea Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <noah.biederbeck@tu-dortmund.de>
Date: Tue, 15 Sep 2020 08:32:13 +0200
Subject: [PATCH 13/17] Only calculate likelihood when mask is not empty

---
 ctapipe/image/pixel_likelihood.py | 14 ++++++++------
 1 file changed, 8 insertions(+), 6 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index 495a9d1101a..13bcb5d3e37 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -177,13 +177,15 @@ def neg_log_likelihood(
     approx_mask = prediction > prediction_safety
 
     neg_log_l = 0
-    neg_log_l += neg_log_likelihood_approx(
-        image[approx_mask], prediction[approx_mask], spe_width, pedestal,
-    )
+    if np.sum(approx_mask) > 0:
+        neg_log_l += neg_log_likelihood_approx(
+            image[approx_mask], prediction[approx_mask], spe_width, pedestal,
+        )
 
-    neg_log_l += neg_log_likelihood_numeric(
-        image[~approx_mask], prediction[~approx_mask], spe_width, pedestal,
-    )
+    if np.sum(~approx_mask) > 0:
+        neg_log_l += neg_log_likelihood_numeric(
+            image[~approx_mask], prediction[~approx_mask], spe_width, pedestal,
+        )
 
     return neg_log_l
 

From af77901d75dece6417d9d33abf20942ed0bb540e Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <noah.biederbeck@tu-dortmund.de>
Date: Tue, 15 Sep 2020 08:34:08 +0200
Subject: [PATCH 14/17] Fix tests.

* change names
* use numpy arrays
* use sum, because we calculate telescope-likelihood not pixel-likelihood
---
 ctapipe/image/tests/test_pixel_likelihood.py | 31 ++++++++++----------
 1 file changed, 15 insertions(+), 16 deletions(-)

diff --git a/ctapipe/image/tests/test_pixel_likelihood.py b/ctapipe/image/tests/test_pixel_likelihood.py
index 3ae6aa401ee..4d86f279a97 100644
--- a/ctapipe/image/tests/test_pixel_likelihood.py
+++ b/ctapipe/image/tests/test_pixel_likelihood.py
@@ -1,5 +1,5 @@
 import numpy as np
-from ctapipe.image import poisson_likelihood_full, poisson_likelihood_gaussian
+from ctapipe.image import neg_log_likelihood, neg_log_likelihood_approx
 
 
 def test_full_likelihood():
@@ -11,30 +11,29 @@ def test_full_likelihood():
     spe = 0.5  # Single photo-electron width
     pedestal = 1  # width of the pedestal distribution
 
-    image_small = [0, 1, 2]
-    expectation_small = [1, 1, 1]
+    image_small = np.array([0, 1, 2])
+    expectation_small = np.array([1, 1, 1])
 
-    full_like_small = poisson_likelihood_full(
-        image_small, expectation_small, spe, pedestal
+    full_like_small = neg_log_likelihood(image_small, expectation_small, spe, pedestal)
+    exp_diff = full_like_small - np.sum(
+        np.asarray([2.75630505, 2.62168656, 3.39248449])
     )
-    exp_diff = full_like_small - np.asarray([2.75630505, 2.62168656, 3.39248449])
-    exp_diff = np.sum(np.abs(exp_diff))
+
     # Check against known values
-    print(exp_diff)
     assert exp_diff / np.sum(full_like_small) < 1e-4
 
-    image_large = [40, 50, 60]
-    expectation_large = [50, 50, 50]
-    full_like_large = poisson_likelihood_full(
-        image_large, expectation_large, spe, pedestal
-    )
+    image_large = np.array([40, 50, 60])
+    expectation_large = np.array([50, 50, 50])
+
+    full_like_large = neg_log_likelihood(image_large, expectation_large, spe, pedestal)
     # Check against known values
-    exp_diff = full_like_large - np.asarray([7.45489137, 5.99305388, 7.66226007])
-    exp_diff = np.sum(np.abs(exp_diff))
+    exp_diff = full_like_large - np.sum(
+        np.asarray([7.45489137, 5.99305388, 7.66226007])
+    )
 
     assert exp_diff / np.sum(full_like_large) < 1e-4
 
-    gaus_like_large = poisson_likelihood_gaussian(
+    gaus_like_large = neg_log_likelihood_approx(
         image_large, expectation_large, spe, pedestal
     )
 

From c2d842ff2d21904badc328b9f79d3a311b0223b0 Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <noah.biederbeck@tu-dortmund.de>
Date: Tue, 15 Sep 2020 11:38:12 +0200
Subject: [PATCH 15/17] Use np.any and not np.all

* also black fixes that were somehow skipped earlier
---
 ctapipe/image/pixel_likelihood.py | 21 +++++++--------------
 1 file changed, 7 insertions(+), 14 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index 13bcb5d3e37..e5d61e28761 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -146,9 +146,7 @@ def neg_log_likelihood_numeric(
     return -np.sum(np.log(likelihood))
 
 
-def neg_log_likelihood(
-    image, prediction, spe_width, pedestal, prediction_safety=20.0,
-):
+def neg_log_likelihood(image, prediction, spe_width, pedestal, prediction_safety=20.0):
     """
     Safe implementation of the poissonian likelihood implementation,
     adaptively switches between the full solution and the gaussian
@@ -177,14 +175,14 @@ def neg_log_likelihood(
     approx_mask = prediction > prediction_safety
 
     neg_log_l = 0
-    if np.sum(approx_mask) > 0:
+    if np.any(approx_mask):
         neg_log_l += neg_log_likelihood_approx(
-            image[approx_mask], prediction[approx_mask], spe_width, pedestal,
+            image[approx_mask], prediction[approx_mask], spe_width, pedestal
         )
 
-    if np.sum(~approx_mask) > 0:
+    if not np.all(approx_mask):
         neg_log_l += neg_log_likelihood_numeric(
-            image[~approx_mask], prediction[~approx_mask], spe_width, pedestal,
+            image[~approx_mask], prediction[~approx_mask], spe_width, pedestal
         )
 
     return neg_log_l
@@ -257,10 +255,7 @@ def mean_poisson_likelihood_full(prediction, spe_width, ped):
     width = np.sqrt(width)
 
     for p in range(len(prediction)):
-        int_range = (
-            prediction[p] - 10 * width[p],
-            prediction[p] + 10 * width[p],
-        )
+        int_range = (prediction[p] - 10 * width[p], prediction[p] + 10 * width[p])
         mean_like[p] = quad(
             _integral_poisson_likelihood_full,
             int_range[0],
@@ -298,9 +293,7 @@ def chi_squared(image, prediction, ped, error_factor=2.9):
     if image.shape is not prediction.shape:
         PixelLikelihoodError(
             "Image and prediction arrays have different dimensions Image "
-            "shape: {} Prediction shape: {}".format(
-                image.shape, prediction.shape
-            )
+            "shape: {} Prediction shape: {}".format(image.shape, prediction.shape)
         )
 
     chi_square = (image - prediction) * (image - prediction)

From 9244b5f85d9f3c11d99d25982872580dc408fbfb Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <noah.biederbeck@tu-dortmund.de>
Date: Wed, 16 Sep 2020 15:02:46 +0200
Subject: [PATCH 16/17] Add tests

---
 ctapipe/image/tests/test_pixel_likelihood.py | 42 +++++++++++++++++++-
 1 file changed, 41 insertions(+), 1 deletion(-)

diff --git a/ctapipe/image/tests/test_pixel_likelihood.py b/ctapipe/image/tests/test_pixel_likelihood.py
index 4d86f279a97..9a2b5bff49e 100644
--- a/ctapipe/image/tests/test_pixel_likelihood.py
+++ b/ctapipe/image/tests/test_pixel_likelihood.py
@@ -1,5 +1,45 @@
 import numpy as np
-from ctapipe.image import neg_log_likelihood, neg_log_likelihood_approx
+from ctapipe.image import (
+    neg_log_likelihood,
+    neg_log_likelihood_approx,
+    mean_poisson_likelihood_gaussian,
+    chi_squared,
+    mean_poisson_likelihood_full,
+)
+
+
+def test_chi_squared():
+    image = np.array([20, 20, 20])
+    prediction = np.array([20, 20, 20])
+    bad_prediction = np.array([1, 1, 1])
+
+    ped = 1
+
+    chi = chi_squared(image, prediction, ped)
+    bad_chi = chi_squared(image, bad_prediction, ped)
+
+    assert chi < bad_chi
+
+
+def test_mean_poisson_likelihoood_gaussian():
+    prediction = np.array([1, 1, 1], dtype="float")
+    spe = 0.5
+
+    small_mean_likelihood = mean_poisson_likelihood_gaussian(prediction, spe, 0)
+    large_mean_likelihood = mean_poisson_likelihood_gaussian(prediction, spe, 1)
+
+    assert small_mean_likelihood < large_mean_likelihood
+
+
+def test_mean_poisson_likelihood_full():
+    prediction = np.array([30.0, 30.0])
+
+    spe = np.array([0.5])
+
+    small_mean_likelihood = mean_poisson_likelihood_full(prediction, spe, [0])
+    large_mean_likelihood = mean_poisson_likelihood_full(prediction, spe, [1])
+
+    assert small_mean_likelihood < large_mean_likelihood
 
 
 def test_full_likelihood():

From 49bc6a4539354bee3d7f734f71f92bd91cdb94af Mon Sep 17 00:00:00 2001
From: Noah Biederbeck <noah.biederbeck@tu-dortmund.de>
Date: Wed, 16 Sep 2020 15:04:03 +0200
Subject: [PATCH 17/17] Change wording, simpler loop

---
 ctapipe/image/pixel_likelihood.py | 82 +++++++++++++++----------------
 1 file changed, 39 insertions(+), 43 deletions(-)

diff --git a/ctapipe/image/pixel_likelihood.py b/ctapipe/image/pixel_likelihood.py
index e5d61e28761..5ca857ea1e3 100644
--- a/ctapipe/image/pixel_likelihood.py
+++ b/ctapipe/image/pixel_likelihood.py
@@ -212,12 +212,14 @@ def mean_poisson_likelihood_gaussian(prediction, spe_width, pedestal):
     return np.sum(mean_log_likelihood)
 
 
-def _integral_poisson_likelihood_full(s, prediction, spe_width, ped):
+def _integral_poisson_likelihood_full(image, prediction, spe_width, ped):
     """
     Wrapper function around likelihood calculation, used in numerical
     integration.
     """
-    like = neg_log_likelihood(s, prediction, spe_width, ped)
+    image = np.asarray(image)
+    prediction = np.asarray(prediction)
+    like = neg_log_likelihood(image, prediction, spe_width, ped)
     return like * np.exp(-0.5 * like)
 
 
@@ -234,70 +236,64 @@ def mean_poisson_likelihood_full(prediction, spe_width, ped):
     prediction: ndarray
         Predicted pixel amplitudes from model
     spe_width: ndarray
-        width of single p.e. distribution
-    ped: ndarray
-        width of pedestal
+        Width of single p.e. distribution
+    pedestal: ndarray
+        Width of pedestal
 
     Returns
     -------
-    ndarray: mean likelihood for give pixel expectation
+    float
     """
-    prediction = np.asarray(prediction)
-    spe_width = np.asarray(spe_width)
-
-    if len(spe_width.shape) == 0:
-        spe_width = np.ones(prediction.shape) * spe_width
-    ped = np.asarray(ped)
-    if len(ped.shape) == 0:
-        ped = np.ones(prediction.shape) * ped
-    mean_like = np.zeros(prediction.shape)
-    width = ped * ped + prediction * spe_width * spe_width
+
+    if len(spe_width) == 1:
+        spe_width = np.full_like(prediction, spe_width)
+
+    if len(ped) == 1:
+        ped = np.full_like(prediction, ped)
+
+    mean_like = 0
+
+    width = ped ** 2 + prediction * spe_width ** 2
     width = np.sqrt(width)
 
-    for p in range(len(prediction)):
-        int_range = (prediction[p] - 10 * width[p], prediction[p] + 10 * width[p])
-        mean_like[p] = quad(
+    for pred, w, spe, p in zip(prediction, width, spe_width, ped):
+        lower_integration_bound = pred - 10 * w
+        upper_integration_bound = pred + 10 * w
+
+        integral, *_ = quad(
             _integral_poisson_likelihood_full,
-            int_range[0],
-            int_range[1],
-            args=(prediction[p], spe_width[p], ped[p]),
+            lower_integration_bound,
+            upper_integration_bound,
+            args=(pred, spe, p),
             epsrel=0.05,
-        )[0]
+        )
+
+        mean_like += integral
+
     return mean_like
 
 
-def chi_squared(image, prediction, ped, error_factor=2.9):
+def chi_squared(image, prediction, pedestal, error_factor=2.9):
     """
     Simple chi-squared statistic from Le Bohec et al 2008
 
     Parameters
     ----------
     image: ndarray
-        Pixel amplitudes from image
+        Pixel amplitudes from image (:math:`s`).
     prediction: ndarray
-        Predicted pixel amplitudes from model
-    ped: ndarray
-        width of pedestal
+        Predicted pixel amplitudes from model (:math:`μ`).
+    pedestal: ndarray
+        Width of pedestal (:math:`σ_p`).
     error_factor: float
-        ad hoc error factor
+        Ad-hoc error factor
 
     Returns
     -------
-    ndarray: likelihood for each pixel
+    float
     """
 
-    image = np.asarray(image)
-    prediction = np.asarray(prediction)
-    ped = np.asarray(ped)
-
-    if image.shape is not prediction.shape:
-        PixelLikelihoodError(
-            "Image and prediction arrays have different dimensions Image "
-            "shape: {} Prediction shape: {}".format(image.shape, prediction.shape)
-        )
-
-    chi_square = (image - prediction) * (image - prediction)
-    chi_square /= ped + 0.5 * (image - prediction)
+    chi_square = (image - prediction) ** 2 / (pedestal + 0.5 * (image - prediction))
     chi_square *= 1.0 / error_factor
 
-    return chi_square
+    return np.sum(chi_square)