From 200aa94ca80d8dbc0895c2e77765b77272b07de0 Mon Sep 17 00:00:00 2001 From: toni Date: Mon, 17 Apr 2017 16:50:56 +0200 Subject: [PATCH] fixed some bugs in jensen shannon and kullback leibler --- main.cpp | 2 +- math/Random.h | 4 +- math/divergence/JensenShannon.h | 7 ++ math/divergence/KullbackLeibler.h | 4 +- sensors/radio/VAPGrouper.h | 2 +- tests/math/divergence/TestJensenShannon.cpp | 94 +++++++++++++++++++++ tests/math/divergence/TestJensenShannon.h | 0 7 files changed, 107 insertions(+), 6 deletions(-) create mode 100644 tests/math/divergence/TestJensenShannon.cpp create mode 100644 tests/math/divergence/TestJensenShannon.h diff --git a/main.cpp b/main.cpp index 2808d16..46f1192 100755 --- a/main.cpp +++ b/main.cpp @@ -30,7 +30,7 @@ int main(int argc, char** argv) { //::testing::GTEST_FLAG(filter) = "*Offline.readWrite*"; - ::testing::GTEST_FLAG(filter) = "*Angle*"; + ::testing::GTEST_FLAG(filter) = "*Jensen*"; //::testing::GTEST_FLAG(filter) = "*Barometer*"; diff --git a/math/Random.h b/math/Random.h index 6dd8f8f..ae8d8c5 100644 --- a/math/Random.h +++ b/math/Random.h @@ -18,10 +18,10 @@ class RandomGenerator : public std::minstd_rand { public: /** ctor with default seed */ - RandomGenerator() : std::minstd_rand(RANDOM_SEED) {;} + RandomGenerator() : std::minstd_rand(RANDOM_SEED) {;} /** ctor with custom seed */ - RandomGenerator(result_type seed) : std::minstd_rand(seed) {;} + RandomGenerator(result_type seed) : std::minstd_rand(seed) {;} }; diff --git a/math/divergence/JensenShannon.h b/math/divergence/JensenShannon.h index 5c9bbb3..df64088 100644 --- a/math/divergence/JensenShannon.h +++ b/math/divergence/JensenShannon.h @@ -18,6 +18,13 @@ public: * @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){ + // normalize + P /= P.sum(); + Q /= Q.sum(); + + Assert::isNear((double)P.sum(), 1.0, 0.01,"Normalization failed.. this shouldn't happen"); + Assert::isNear((double)Q.sum(), 1.0, 0.01, "Normalization failed.. this shouldn't happen"); + Eigen::VectorXd M = 0.5 * (P + Q); return (0.5 * KullbackLeibler::getGeneralFromSamples(P, M, mode)) + (0.5 * KullbackLeibler::getGeneralFromSamples(Q, M, mode)); diff --git a/math/divergence/KullbackLeibler.h b/math/divergence/KullbackLeibler.h index 630e5d2..27124e7 100644 --- a/math/divergence/KullbackLeibler.h +++ b/math/divergence/KullbackLeibler.h @@ -56,8 +56,8 @@ namespace Divergence { // calc PQratio if(Q[i] == 0.0){ - Assert::doThrow("Division by zero is not allowed ;)."); - //PQratio = P[i] / 0.00001; + //Assert::doThrow("Division by zero is not allowed ;)."); + PQratio = P[i] / 0.00001; } else { PQratio = P[i] / Q[i]; } diff --git a/sensors/radio/VAPGrouper.h b/sensors/radio/VAPGrouper.h index b8a7545..98211e2 100644 --- a/sensors/radio/VAPGrouper.h +++ b/sensors/radio/VAPGrouper.h @@ -89,7 +89,7 @@ public: } - Log::add(name, "grouped " + std::to_string(original.entries.size()) + " measurements into " + std::to_string(result.entries.size()), true); + //Log::add(name, "grouped " + std::to_string(original.entries.size()) + " measurements into " + std::to_string(result.entries.size()), true); // done diff --git a/tests/math/divergence/TestJensenShannon.cpp b/tests/math/divergence/TestJensenShannon.cpp new file mode 100644 index 0000000..4bf8cf3 --- /dev/null +++ b/tests/math/divergence/TestJensenShannon.cpp @@ -0,0 +1,94 @@ +#ifdef WITH_TESTS + +#include "../../Tests.h" +#include "../../../math/divergence/JensenShannon.h" +#include "../../../math/Distributions.h" + +#include + + +TEST(JensenShannon, generalFromSamples) { + //ge cov + Eigen::VectorXd mu1(2); + mu1[0] = 1.0; + mu1[1] = 1.0; + + Eigen::VectorXd mu2(2); + mu2[0] = 1.0; + mu2[1] = 1.0; + + Eigen::VectorXd mu3(2); + mu3[0] = 1.0; + mu3[1] = 1.0; + + Eigen::VectorXd mu4(2); + mu4[0] = 9.0; + mu4[1] = 9.0; + + Eigen::MatrixXd cov1(2,2); + cov1(0,0) = 1.0; + cov1(0,1) = 0.0; + cov1(1,0) = 0.0; + cov1(1,1) = 1.0; + + Eigen::MatrixXd cov2(2,2); + cov2(0,0) = 1.0; + cov2(0,1) = 0.0; + cov2(1,0) = 0.0; + cov2(1,1) = 1.0; + + Eigen::MatrixXd cov3(2,2); + cov3(0,0) = 1.0; + cov3(0,1) = 0.0; + cov3(1,0) = 0.0; + cov3(1,1) = 1.0; + + Eigen::MatrixXd cov4(2,2); + cov4(0,0) = 3.0; + cov4(0,1) = 0.0; + cov4(1,0) = 0.0; + cov4(1,1) = 13.0; + + Distribution::NormalDistributionN norm1(mu1, cov1); + Distribution::NormalDistributionN norm2(mu2, cov2); + Distribution::NormalDistributionN norm3(mu3, cov3); + Distribution::NormalDistributionN norm4(mu4, cov4); + + int size = 10000; + Eigen::VectorXd samples1(size); + Eigen::VectorXd samples2(size); + Eigen::VectorXd samples3(size); + Eigen::VectorXd samples4(size); + + //random numbers + std::mt19937_64 rng; + // initialize the random number generator with time-dependent seed + uint64_t timeSeed = std::chrono::high_resolution_clock::now().time_since_epoch().count(); + std::seed_seq ss{uint32_t(timeSeed & 0xffffffff), uint32_t(timeSeed>>32)}; + rng.seed(ss); + // initialize a uniform distribution between 0 and 1 + std::uniform_real_distribution unif(-9, 10); + + //generate samples + for(int i = 0; i < size; ++i){ + + double r1 = unif(rng); + double r2 = unif(rng); + Eigen::VectorXd v(2); + v << r1, r2; + + samples1[i] = norm1.getProbability(v); + samples2[i] = norm2.getProbability(v); + samples3[i] = norm3.getProbability(v); + samples4[i] = norm4.getProbability(v); + } + + double kld12 = Divergence::JensenShannon::getGeneralFromSamples(samples1, samples2, Divergence::LOGMODE::NATURALIS); + double kld34 = Divergence::JensenShannon::getGeneralFromSamples(samples3, samples4, Divergence::LOGMODE::NATURALIS); + std::cout << kld34 << " > " << kld12 << std::endl; + + ASSERT_GE(kld34, kld12); + +} + +#endif diff --git a/tests/math/divergence/TestJensenShannon.h b/tests/math/divergence/TestJensenShannon.h new file mode 100644 index 0000000..e69de29