Skip to content

Commit

Permalink
[ML] Refactor code to use std::lgamma instead of boost::math::lgamma (#…
Browse files Browse the repository at this point in the history 
…132)

Migrate remaining codebase to std::lgamma, std::{isfinite,isinf,isnan}andstd::erf`

fixes #128
  • Loading branch information
Hendrik Muhs committed Jun 20, 2018
1 parent 3d394c6 commit e78785c
Show file tree
Hide file tree
Showing 19 changed files with 260 additions and 263 deletions.
8 changes: 4 additions & 4 deletions include/core/CRapidJsonWriterBase.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -20,10 +20,10 @@
#include <rapidjson/writer.h>

#include <boost/iterator/indirect_iterator.hpp>
#include <boost/math/special_functions/fpclassify.hpp>
#include <boost/unordered/unordered_map.hpp>
#include <boost/unordered/unordered_set.hpp>

#include <cmath>
#include <memory>
#include <stack>

Expand Down Expand Up @@ -166,7 +166,7 @@ class CRapidJsonWriterBase

bool Double(double d) {
// rewrite NaN and Infinity to 0
if (!(boost::math::isfinite)(d)) {
if (std::isfinite(d) == false) {
return TRapidJsonWriterBase::Int(0);
}

Expand Down Expand Up @@ -355,7 +355,7 @@ class CRapidJsonWriterBase
//! Adds a double field with the name fieldname to an object.
//! \p fieldName must outlive \p obj or memory corruption will occur.
void addDoubleFieldToObj(const std::string& fieldName, double value, TValue& obj) const {
if (!(boost::math::isfinite)(value)) {
if (std::isfinite(value) == false) {
LOG_ERROR(<< "Adding " << value << " to the \"" << fieldName
<< "\" field of a JSON document");
// Don't return - make a best effort to add the value
Expand Down Expand Up @@ -489,7 +489,7 @@ class CRapidJsonWriterBase
//! Log a message if we're trying to add nan/infinity to a JSON array
template<typename NUMBER>
void checkArrayNumberFinite(NUMBER val, const std::string& fieldName, bool& considerLogging) const {
if (considerLogging && !(boost::math::isfinite)(val)) {
if (considerLogging && (std::isfinite(val) == false)) {
LOG_ERROR(<< "Adding " << val << " to the \"" << fieldName
<< "\" array in a JSON document");
// Don't return - make a best effort to add the value
Expand Down
6 changes: 2 additions & 4 deletions include/maths/CMathsFuncs.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -38,8 +38,7 @@ namespace maths {
//!
class MATHS_EXPORT CMathsFuncs : private core::CNonInstantiatable {
public:
//! Wrapper around boost::math::isnan() which avoids the need to add
//! cryptic brackets everywhere to deal with macros.
//! Check if value is NaN.
static bool isNan(double val);
//! Check if any of the components are NaN.
template<std::size_t N>
Expand All @@ -54,8 +53,7 @@ class MATHS_EXPORT CMathsFuncs : private core::CNonInstantiatable {
//! Check if an element is NaN.
static bool isNan(const core::CSmallVectorBase<double>& val);

//! Wrapper around boost::math::isinf() which avoids the need to add
//! cryptic brackets everywhere to deal with macros.
//! Check if value is infinite.
static bool isInf(double val);
//! Check if any of the components are infinite.
template<std::size_t N>
Expand Down
67 changes: 34 additions & 33 deletions include/maths/CMultivariateNormalConjugate.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,11 +28,11 @@
#include <maths/CNormalMeanPrecConjugate.h>
#include <maths/CRestoreParams.h>
#include <maths/CSampling.h>
#include <maths/CTools.h>
#include <maths/ProbabilityAggregators.h>

#include <boost/make_unique.hpp>
#include <boost/math/special_functions/beta.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/numeric/conversion/bounds.hpp>
#include <boost/optional.hpp>

Expand Down Expand Up @@ -1079,41 +1079,42 @@ class CMultivariateNormalConjugate : public CMultivariatePrior {
LOG_ERROR(<< "Failed to calculate log det " << wishartScaleMatrixPost);
return maths_t::E_FpFailed;
}
double logGammaPost = 0.0;
double logGammaPrior = 0.0;
double logGammaPostMinusPrior = 0.0;

try {
double logGammaPostMinusPrior = 0.0;
for (std::size_t i = 0u; i < N; ++i) {
logGammaPostMinusPrior +=
boost::math::lgamma(
0.5 * (wishartDegreesFreedomPost - static_cast<double>(i))) -
boost::math::lgamma(0.5 * (wishartDegreesFreedomPrior -
static_cast<double>(i)));
for (std::size_t i = 0u; i < N; ++i) {
if (CTools::lgamma(0.5 * (wishartDegreesFreedomPost - static_cast<double>(i)),
logGammaPost, true) &&
CTools::lgamma(0.5 * (wishartDegreesFreedomPrior - static_cast<double>(i)),
logGammaPrior, true)) {

logGammaPostMinusPrior += logGammaPost - logGammaPrior;
} else {
return maths_t::E_FpOverflowed;
}
LOG_TRACE(<< "numberSamples = " << numberSamples);
LOG_TRACE(<< "logGaussianPrecisionPrior = " << logGaussianPrecisionPrior
<< ", logGaussianPrecisionPost = " << logGaussianPrecisionPost);
LOG_TRACE(<< "wishartDegreesFreedomPrior = " << wishartDegreesFreedomPrior
<< ", wishartDegreesFreedomPost = " << wishartDegreesFreedomPost);
LOG_TRACE(<< "wishartScaleMatrixPrior = " << m_WishartScaleMatrix);
LOG_TRACE(<< "wishartScaleMatrixPost = " << wishartScaleMatrixPost);
LOG_TRACE(<< "logDeterminantPrior = " << logDeterminantPrior
<< ", logDeterminantPost = " << logDeterminantPost);
LOG_TRACE(<< "logGammaPostMinusPrior = " << logGammaPostMinusPrior);
LOG_TRACE(<< "logCountVarianceScales = " << logCountVarianceScales);

double d = static_cast<double>(N);
result = 0.5 * (wishartDegreesFreedomPrior * logDeterminantPrior -
wishartDegreesFreedomPost * logDeterminantPost -
d * (logGaussianPrecisionPost - logGaussianPrecisionPrior) +
(wishartDegreesFreedomPost - wishartDegreesFreedomPrior) *
d * core::constants::LOG_TWO +
2.0 * logGammaPostMinusPrior -
numberSamples * d * core::constants::LOG_TWO_PI -
logCountVarianceScales);
} catch (const std::exception& e) {
LOG_ERROR(<< "Failed to calculate marginal likelihood: " << e.what());
return maths_t::E_FpFailed;
}
LOG_TRACE(<< "numberSamples = " << numberSamples);
LOG_TRACE(<< "logGaussianPrecisionPrior = " << logGaussianPrecisionPrior
<< ", logGaussianPrecisionPost = " << logGaussianPrecisionPost);
LOG_TRACE(<< "wishartDegreesFreedomPrior = " << wishartDegreesFreedomPrior
<< ", wishartDegreesFreedomPost = " << wishartDegreesFreedomPost);
LOG_TRACE(<< "wishartScaleMatrixPrior = " << m_WishartScaleMatrix);
LOG_TRACE(<< "wishartScaleMatrixPost = " << wishartScaleMatrixPost);
LOG_TRACE(<< "logDeterminantPrior = " << logDeterminantPrior
<< ", logDeterminantPost = " << logDeterminantPost);
LOG_TRACE(<< "logGammaPostMinusPrior = " << logGammaPostMinusPrior);
LOG_TRACE(<< "logCountVarianceScales = " << logCountVarianceScales);

double d = static_cast<double>(N);
result = 0.5 * (wishartDegreesFreedomPrior * logDeterminantPrior -
wishartDegreesFreedomPost * logDeterminantPost -
d * (logGaussianPrecisionPost - logGaussianPrecisionPrior) +
(wishartDegreesFreedomPost - wishartDegreesFreedomPrior) *
d * core::constants::LOG_TWO +
2.0 * logGammaPostMinusPrior -
numberSamples * d * core::constants::LOG_TWO_PI - logCountVarianceScales);

return static_cast<maths_t::EFloatingPointErrorStatus>(CMathsFuncs::fpStatus(result));
}

Expand Down
2 changes: 1 addition & 1 deletion include/maths/CTools.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -695,7 +695,7 @@ class MATHS_EXPORT CTools : private core::CNonInstantiatable {
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);
static bool lgamma(double value, double& result, bool checkForFinite = true);
};
}
}
Expand Down
33 changes: 26 additions & 7 deletions lib/maths/CCategoricalTools.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,7 +14,6 @@
#include <maths/MathsTypes.h>

#include <boost/math/distributions/binomial.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/numeric/conversion/bounds.hpp>

#include <cmath>
Expand Down Expand Up @@ -408,8 +407,7 @@ double CCategoricalTools::logBinomialCoefficient(std::size_t n, std::size_t m) {
}
double n_ = static_cast<double>(n);
double m_ = static_cast<double>(m);
return boost::math::lgamma(n_ + 1.0) - boost::math::lgamma(m_ + 1.0) -
boost::math::lgamma(n_ - m_ + 1.0);
return std::lgamma(n_ + 1.0) - std::lgamma(m_ + 1.0) - std::lgamma(n_ - m_ + 1.0);
}

double CCategoricalTools::binomialCoefficient(std::size_t n, std::size_t m) {
Expand Down Expand Up @@ -544,8 +542,19 @@ CCategoricalTools::logBinomialProbability(std::size_t n, double p, std::size_t m

double n_ = static_cast<double>(n);
double m_ = static_cast<double>(m);
result = std::min(boost::math::lgamma(n_ + 1.0) - boost::math::lgamma(m_ + 1.0) -
boost::math::lgamma(n_ - m_ + 1.0) +

double logGammaNPlusOne = 0.0;
double logGammaMPlusOne = 0.0;
double logGammaNMinusMPlusOne = 0.0;

if ((CTools::lgamma(n_ + 1.0, logGammaNPlusOne, true) &&
CTools::lgamma(m_ + 1.0, logGammaMPlusOne, true) &&
CTools::lgamma(n_ - m_ + 1.0, logGammaNMinusMPlusOne, true)) == false) {

return maths_t::E_FpOverflowed;
}

result = std::min(logGammaNPlusOne - logGammaMPlusOne - logGammaNMinusMPlusOne +
m_ * std::log(p) + (n_ - m_) * std::log(1.0 - p),
0.0);
return maths_t::E_FpNoErrors;
Expand All @@ -570,7 +579,11 @@ CCategoricalTools::logMultinomialProbability(const TDoubleVec& probabilities,
}

double n_ = static_cast<double>(n);
double logP = boost::math::lgamma(n_ + 1.0);
double logP = 0.0;

if (CTools::lgamma(n_ + 1.0, logP, true) == false) {
return maths_t::E_FpOverflowed;
}

for (std::size_t i = 0u; i < ni.size(); ++i) {
double ni_ = static_cast<double>(ni[i]);
Expand All @@ -584,7 +597,13 @@ CCategoricalTools::logMultinomialProbability(const TDoubleVec& probabilities,
result = boost::numeric::bounds<double>::lowest();
return maths_t::E_FpOverflowed;
}
logP += ni_ * std::log(pi_) - boost::math::lgamma(ni_ + 1.0);

double logGammaNiPlusOne = 0.0;
if (CTools::lgamma(ni_ + 1.0, logGammaNiPlusOne, true) == false) {
return maths_t::E_FpOverflowed;
}

logP += ni_ * std::log(pi_) - logGammaNiPlusOne;
}
}

Expand Down
51 changes: 26 additions & 25 deletions lib/maths/CGammaRateConjugate.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -27,7 +27,6 @@
#include <boost/math/distributions/beta.hpp>
#include <boost/math/distributions/gamma.hpp>
#include <boost/math/special_functions/digamma.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/math/tools/roots.hpp>
#include <boost/numeric/conversion/bounds.hpp>

Expand Down Expand Up @@ -637,35 +636,37 @@ class CLogMarginalLikelihood : core::CNonCopyable {
double logGammaScaledLikelihoodShape = 0.0;
double scaledImpliedShape = 0.0;

try {
for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
double n = maths_t::countForUpdate(m_Weights[i]);
double varianceScale = maths_t::seasonalVarianceScale(m_Weights[i]) *
maths_t::countVarianceScale(m_Weights[i]);
m_NumberSamples += n;
if (varianceScale != 1.0) {
logVarianceScaleSum -= m_LikelihoodShape / varianceScale *
std::log(varianceScale);
logGammaScaledLikelihoodShape +=
n * boost::math::lgamma(m_LikelihoodShape / varianceScale);
scaledImpliedShape += n * m_LikelihoodShape / varianceScale;
} else {
nResidual += n;
}
for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
double n = maths_t::countForUpdate(m_Weights[i]);
double varianceScale = maths_t::seasonalVarianceScale(m_Weights[i]) *
maths_t::countVarianceScale(m_Weights[i]);
m_NumberSamples += n;
if (varianceScale != 1.0) {
logVarianceScaleSum -= m_LikelihoodShape / varianceScale *
std::log(varianceScale);
logGammaScaledLikelihoodShape +=
n * std::lgamma(m_LikelihoodShape / varianceScale);
scaledImpliedShape += n * m_LikelihoodShape / varianceScale;
} else {
nResidual += n;
}
}

m_ImpliedShape = scaledImpliedShape + nResidual * m_LikelihoodShape + m_PriorShape;
m_ImpliedShape = scaledImpliedShape + nResidual * m_LikelihoodShape + m_PriorShape;

LOG_TRACE(<< "numberSamples = " << m_NumberSamples);
LOG_TRACE(<< "numberSamples = " << m_NumberSamples);

m_Constant = m_PriorShape * std::log(m_PriorRate) -
boost::math::lgamma(m_PriorShape) +
logVarianceScaleSum - logGammaScaledLikelihoodShape -
nResidual * boost::math::lgamma(m_LikelihoodShape) +
boost::math::lgamma(m_ImpliedShape);
} catch (const std::exception& e) {
LOG_ERROR(<< "Error calculating marginal likelihood: " << e.what());
m_Constant = m_PriorShape * std::log(m_PriorRate) - std::lgamma(m_PriorShape) +
logVarianceScaleSum - logGammaScaledLikelihoodShape -
nResidual * std::lgamma(m_LikelihoodShape) +
std::lgamma(m_ImpliedShape);

if (std::isnan(m_ImpliedShape) || std::isnan(m_Constant)) {
LOG_ERROR(<< "Error calculating marginal likelihood: floating point nan");
this->addErrorStatus(maths_t::E_FpFailed);
} else if (std::isinf(m_ImpliedShape) || std::isinf(m_Constant)) {
LOG_ERROR(<< "Error calculating marginal likelihood: floating point overflow");
this->addErrorStatus(maths_t::E_FpOverflowed);
}
}

Expand Down
77 changes: 39 additions & 38 deletions lib/maths/CLogNormalMeanPrecConjugate.cc
View file
Original file line number Diff line number Diff line change
Expand Up @@ -30,7 +30,6 @@
#include <boost/math/distributions/lognormal.hpp>
#include <boost/math/distributions/normal.hpp>
#include <boost/math/distributions/students_t.hpp>
#include <boost/math/special_functions/gamma.hpp>
#include <boost/numeric/conversion/bounds.hpp>

#include <algorithm>
Expand Down Expand Up @@ -487,49 +486,51 @@ class CLogMarginalLikelihood : core::CNonCopyable {
private:
//! Compute all the constants in the integrand.
void precompute() {
try {
double logVarianceScaleSum = 0.0;

if (maths_t::hasSeasonalVarianceScale(m_Weights) ||
maths_t::hasCountVarianceScale(m_Weights)) {
m_Scales.reserve(m_Weights.size());
double r = m_Rate / m_Shape;
double s = std::exp(-r);
for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
double varianceScale = maths_t::seasonalVarianceScale(m_Weights[i]) *
maths_t::countVarianceScale(m_Weights[i]);

// Get the scale and shift of the exponentiated Gaussian.
if (varianceScale == 1.0) {
m_Scales.emplace_back(1.0, 0.0);
} else {
double t = r + std::log(s + varianceScale * (1.0 - s));
m_Scales.emplace_back(t / r, 0.5 * (r - t));
logVarianceScaleSum += std::log(t / r);
}
double logVarianceScaleSum = 0.0;

if (maths_t::hasSeasonalVarianceScale(m_Weights) ||
maths_t::hasCountVarianceScale(m_Weights)) {
m_Scales.reserve(m_Weights.size());
double r = m_Rate / m_Shape;
double s = std::exp(-r);
for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
double varianceScale = maths_t::seasonalVarianceScale(m_Weights[i]) *
maths_t::countVarianceScale(m_Weights[i]);

// Get the scale and shift of the exponentiated Gaussian.
if (varianceScale == 1.0) {
m_Scales.emplace_back(1.0, 0.0);
} else {
double t = r + std::log(s + varianceScale * (1.0 - s));
m_Scales.emplace_back(t / r, 0.5 * (r - t));
logVarianceScaleSum += std::log(t / r);
}
}
}

m_NumberSamples = 0.0;
double weightedNumberSamples = 0.0;
m_NumberSamples = 0.0;
double weightedNumberSamples = 0.0;

for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
double n = maths_t::countForUpdate(m_Weights[i]);
m_NumberSamples += n;
weightedNumberSamples +=
n / (m_Scales.empty() ? 1.0 : m_Scales[i].first);
}
for (std::size_t i = 0u; i < m_Weights.size(); ++i) {
double n = maths_t::countForUpdate(m_Weights[i]);
m_NumberSamples += n;
weightedNumberSamples += n / (m_Scales.empty() ? 1.0 : m_Scales[i].first);
}

double impliedShape = m_Shape + 0.5 * m_NumberSamples;
double impliedPrecision = m_Precision + weightedNumberSamples;
double impliedShape = m_Shape + 0.5 * m_NumberSamples;
double impliedPrecision = m_Precision + weightedNumberSamples;

m_Constant = 0.5 * (std::log(m_Precision) - std::log(impliedPrecision)) -
0.5 * m_NumberSamples * LOG_2_PI - 0.5 * logVarianceScaleSum +
boost::math::lgamma(impliedShape) -
boost::math::lgamma(m_Shape) + m_Shape * std::log(m_Rate);
} catch (const std::exception& e) {
LOG_ERROR(<< "Error calculating marginal likelihood: " << e.what());
m_Constant = 0.5 * (std::log(m_Precision) - std::log(impliedPrecision)) -
0.5 * m_NumberSamples * LOG_2_PI -
0.5 * logVarianceScaleSum + std::lgamma(impliedShape) -
std::lgamma(m_Shape) + m_Shape * std::log(m_Rate);

if (std::isnan(m_Constant)) {
LOG_ERROR(<< "Error calculating marginal likelihood, floating point nan");
this->addErrorStatus(maths_t::E_FpFailed);
} else if (std::isinf(m_Constant)) {
LOG_ERROR(<< "Error calculating marginal likelihood, floating point overflow");
this->addErrorStatus(maths_t::E_FpOverflowed);
}
}

Expand Down Expand Up @@ -1219,7 +1220,7 @@ void CLogNormalMeanPrecConjugate::sampleMarginalLikelihood(std::size_t numberSam
double z = (xq - m_GaussianMean - scale * scale) / scale /
boost::math::double_constants::root_two;

double partialExpectation = mean * (1.0 + boost::math::erf(z)) / 2.0;
double partialExpectation = mean * (1.0 + std::erf(z)) / 2.0;

double sample = static_cast<double>(numberSamples) *
(partialExpectation - lastPartialExpectation) -
Expand Down
Loading

0 comments on commit e78785c

Please sign in to comment.