From a25ffb2ba35a71253a2f1c050920e832c8a162b1 Mon Sep 17 00:00:00 2001
From: toni <toni@localhost.localdomain>
Date: Fri, 24 Mar 2017 17:28:57 +0100
Subject: [PATCH] added jensen shannon, fixe kullback leibler when P=0

---
 math/divergence/JensenShannon.h   | 28 ++++++++++++++++++++++++++++
 math/divergence/KullbackLeibler.h | 11 ++++++++---
 2 files changed, 36 insertions(+), 3 deletions(-)
 create mode 100644 math/divergence/JensenShannon.h

diff --git a/math/divergence/JensenShannon.h b/math/divergence/JensenShannon.h
new file mode 100644
index 0000000..a76686d
--- /dev/null
+++ b/math/divergence/JensenShannon.h
@@ -0,0 +1,28 @@
+#ifndef JENSENSHANNON_H
+#define JENSENSHANNON_H
+
+#include "KullbackLeibler.h"
+
+#include "../../Assertions.h"
+#include <string>
+
+
+namespace Divergence {
+
+template <typename Scalar> class JensenShannon {
+
+    /** Calculate the Jensen Shannon Divergece from a set of sample densities
+     *  Info: https://en.wikipedia.org/wiki/Jensen–Shannon_divergence
+     * @param P is the vector containing the densities of a set of samples
+     * @param Q is a vector containg the densities of the same samples set then P
+    */
+    static inline Scalar getGeneralFromSamples(Eigen::VectorXd P, Eigen::VectorXd Q, LOGMODE mode){
+        Eigen::VectorXd M = 0.5 * (P + Q);
+
+        return (0.5 * KullbackLeibler::getGeneralFromSamples(P, M, mode)) + (0.5 * KullbackLeibler::getGeneralFromSamples(Q, M, mode));
+    }
+};
+
+}
+
+#endif // JENSENSHANNON_H
diff --git a/math/divergence/KullbackLeibler.h b/math/divergence/KullbackLeibler.h
index aab08be..630e5d2 100644
--- a/math/divergence/KullbackLeibler.h
+++ b/math/divergence/KullbackLeibler.h
@@ -45,13 +45,18 @@ namespace Divergence {
                 //sum up the logarithmic difference between P and Q, also called kullback leibler...
                 for(int i = 0; i < P.size(); ++i){
 
+                    // if both prob are near zero we assume a 0 distance since lim->0 = 0
                     if((P[i] == 0.0) && (Q[i] == 0.0)){
-                            dist += 0.0;
-                    } else {
+                        dist += 0.0;
+                    }
+                    //if prob of P is 0 we also assume a 0 distance since lim->0 0 * log(0) = 0
+                    else if ((P[i] == 0.0) && (Q[i] > 0.0)){
+                        dist += 0.0;
+                    } else{
 
                         // calc PQratio
                         if(Q[i] == 0.0){
-                            Assert::doThrow("Division by zero is not allowed ;). TODO: What if the densities are to small?");
+                            Assert::doThrow("Division by zero is not allowed ;).");
                             //PQratio = P[i] / 0.00001;
                         } else {
                             PQratio = P[i] / Q[i];