[ML] Fix #51 - fix lgamma calculation for x-means clustering (#126) (#131)

replacing boost::lgamma with std::lgamma and guarding the output of lgamma for finite values.

fixes #51

Author: Hendrik Muhs

docs/CHANGELOG.asciidoc

@@ -40,6 +40,7 @@ by unnecessary reference counting ({pull}108[#108])
 Age seasonal components in proportion to the fraction of values with which they're updated ({pull}88[#88])
 Persist and restore was missing some of the trend model state ({pull}#99[#99])
 Stop zero variance data generating a log error in the forecast confidence interval calculation ({pull}#107[#107])
+Fix corner case failing to calculate lgamma values and the correspoinding log errors ({pull}#126[#126])
 
 === Regressions
 

include/maths/CTools.h

@@ -693,6 +693,9 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
     //! A custom implementation of \f$\log(1 - x)\f$ which handles the
     //! cancellation error for small x.
     static double logOneMinusX(double x);
+
+    //! A wrapper around lgamma which handles corner cases if requested
+    static bool lgamma(double value, double& result, bool checkForFinite = false);
 };
 }
 }

lib/maths/CTools.cc

@@ -2005,5 +2005,10 @@ double CTools::logOneMinusX(double x) {
 
     return result;
 }
+
+bool CTools::lgamma(double value, double& result, bool checkForFinite) {
+    result = std::lgamma(value);
+    return checkForFinite ? std::isfinite(result) : true;
+}
 }
 }

lib/maths/CXMeansOnline1d.cc

@@ -30,7 +30,6 @@
 #include <boost/math/constants/constants.hpp>
 #include <boost/math/distributions/normal.hpp>
 #include <boost/math/special_functions/digamma.hpp>
-#include <boost/math/special_functions/gamma.hpp>
 
 #include <algorithm>
 #include <cmath>
@@ -285,6 +284,7 @@ void BICGain(maths_t::EDataType dataType,
     double mr = mean(dataType, mvr);
     double vr = std::max(variance(dataType, mvr), vmin);
 
+    bool haveGamma = distributions.haveGamma();
     try {
         // Mixture of log-normals (method of moments)
         double sl = std::log(1.0 + vl / CTools::pow2(ml + logNormalOffset));
@@ -299,7 +299,21 @@ void BICGain(maths_t::EDataType dataType,
 
         double log2piv = std::log(boost::math::double_constants::two_pi * v);
         double log2pis = std::log(boost::math::double_constants::two_pi * s);
-        double loggn = boost::math::lgamma(a) - a * std::log(b);
+
+        double loggn = 0.0;
+        double loggnl = 0.0;
+        double loggnr = 0.0;
+
+        if (haveGamma && maths::CTools::lgamma(a, loggn, true) &&
+            maths::CTools::lgamma(al, loggnl, true) &&
+            maths::CTools::lgamma(ar, loggnr, true)) {
+            loggn -= a * std::log(b);
+            loggnl -= al * std::log(bl) + std::log(wl);
+            loggnr -= ar * std::log(br) + std::log(wr);
+        } else {
+            haveGamma = false;
+        }
+
         double log2pivl =
             std::log(boost::math::double_constants::two_pi * vl / CTools::pow2(wl));
         double log2pivr =
@@ -308,8 +322,6 @@ void BICGain(maths_t::EDataType dataType,
             std::log(boost::math::double_constants::two_pi * sl / CTools::pow2(wl));
         double log2pisr =
             std::log(boost::math::double_constants::two_pi * sr / CTools::pow2(wr));
-        double loggnl = boost::math::lgamma(al) - al * std::log(bl) - std::log(wl);
-        double loggnr = boost::math::lgamma(ar) - ar * std::log(br) - std::log(wr);
 
         for (std::size_t i = start; i < split; ++i) {
             double ni = CBasicStatistics::count(categories[i]);
@@ -372,15 +384,15 @@ void BICGain(maths_t::EDataType dataType,
     double ll1 =
         min(distributions.haveNormal() ? ll1n : boost::numeric::bounds<double>::highest(),
             distributions.haveLogNormal() ? ll1l : boost::numeric::bounds<double>::highest(),
-            distributions.haveGamma() ? ll1g : boost::numeric::bounds<double>::highest()) +
+            haveGamma ? ll1g : boost::numeric::bounds<double>::highest()) +
         distributions.parameters() * logn;
     double ll2 =
         min(distributions.haveNormal() ? ll2nl : boost::numeric::bounds<double>::highest(),
             distributions.haveLogNormal() ? ll2ll : boost::numeric::bounds<double>::highest(),
-            distributions.haveGamma() ? ll2gl : boost::numeric::bounds<double>::highest()) +
+            haveGamma ? ll2gl : boost::numeric::bounds<double>::highest()) +
         min(distributions.haveNormal() ? ll2nr : boost::numeric::bounds<double>::highest(),
             distributions.haveLogNormal() ? ll2lr : boost::numeric::bounds<double>::highest(),
-            distributions.haveGamma() ? ll2gr : boost::numeric::bounds<double>::highest()) +
+            haveGamma ? ll2gr : boost::numeric::bounds<double>::highest()) +
         (2.0 * distributions.parameters() + 1.0) * logn;
 
     LOG_TRACE(<< "BIC(1) = " << ll1 << ", BIC(2) = " << ll2);

lib/maths/unittest/CToolsTest.cc

@@ -26,6 +26,8 @@
 #include <boost/optional.hpp>
 #include <boost/range.hpp>
 
+#include <array>
+
 using namespace ml;
 using namespace maths;
 using namespace test;
@@ -1168,6 +1170,54 @@ void CToolsTest::testMiscellaneous() {
     }
 }
 
+void CToolsTest::testLgamma() {
+    std::array<double, 8> testData = {3.5,
+                                      0.125,
+                                      -0.125,
+                                      0.000244140625,
+                                      1.3552527156068805e-20,
+                                      4.95547e+25,
+                                      5.01753e+25,
+                                      8.64197e+25};
+
+    std::array<double, 8> expectedData = {
+        1.2009736023470742248160218814507129957702389154682,
+        2.0194183575537963453202905211670995899482809521344,
+        2.1653002489051702517540619481440174064962195287626,
+        8.3176252939431805089043336920440196990966796875000,
+        45.7477139169563926657247066032141447067260742187500,
+        2.882355039447984e+27,
+        2.919076782442754e+27,
+        5.074673490557339e+27};
+
+    for (std::size_t i = 0u; i < testData.size(); ++i) {
+        double actual;
+        double expected = expectedData[i];
+        CPPUNIT_ASSERT(maths::CTools::lgamma(testData[i], actual, true));
+        CPPUNIT_ASSERT_DOUBLES_EQUAL(expected, actual, 1e-5 * expected);
+    }
+
+    double result;
+    CPPUNIT_ASSERT(maths::CTools::lgamma(0, result));
+    CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
+
+    CPPUNIT_ASSERT((maths::CTools::lgamma(0, result, true) == false));
+    CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
+
+    CPPUNIT_ASSERT((maths::CTools::lgamma(-1, result)));
+    CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
+
+    CPPUNIT_ASSERT((maths::CTools::lgamma(-1, result, true) == false));
+    CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
+
+    CPPUNIT_ASSERT((maths::CTools::lgamma(std::numeric_limits<double>::max() - 1, result)));
+    CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
+
+    CPPUNIT_ASSERT((maths::CTools::lgamma(std::numeric_limits<double>::max() - 1,
+                                          result, true) == false));
+    CPPUNIT_ASSERT_EQUAL(result, std::numeric_limits<double>::infinity());
+}
+
 CppUnit::Test* CToolsTest::suite() {
     CppUnit::TestSuite* suiteOfTests = new CppUnit::TestSuite("CToolsTest");
 
@@ -1187,6 +1237,8 @@ CppUnit::Test* CToolsTest::suite() {
         "CToolsTest::testFastLog", &CToolsTest::testFastLog));
     suiteOfTests->addTest(new CppUnit::TestCaller<CToolsTest>(
         "CToolsTest::testMiscellaneous", &CToolsTest::testMiscellaneous));
+    suiteOfTests->addTest(new CppUnit::TestCaller<CToolsTest>(
+        "CToolsTest::testLgamma", &CToolsTest::testLgamma));
 
     return suiteOfTests;
 }

lib/maths/unittest/CToolsTest.h

@@ -18,6 +18,7 @@ class CToolsTest : public CppUnit::TestFixture {
     void testSpread();
     void testFastLog();
     void testMiscellaneous();
+    void testLgamma();
 
     static CppUnit::Test* suite();
 };