lib/maths/common/CMultimodalPrior.cc (1,142 lines of code) (raw):

/* * Copyright Elasticsearch B.V. and/or licensed to Elasticsearch B.V. under one * or more contributor license agreements. Licensed under the Elastic License * 2.0 and the following additional limitation. Functionality enabled by the * files subject to the Elastic License 2.0 may only be used in production when * invoked by an Elasticsearch process with a license key installed that permits * use of machine learning features. You may not use this file except in * compliance with the Elastic License 2.0 and the foregoing additional * limitation. */ #include <maths/common/CMultimodalPrior.h> #include <core/CLogger.h> #include <core/CMemoryDef.h> #include <core/CSmallVector.h> #include <core/CStatePersistInserter.h> #include <core/CStateRestoreTraverser.h> #include <core/Constants.h> #include <core/RestoreMacros.h> #include <maths/common/CBasicStatistics.h> #include <maths/common/CBasicStatisticsPersist.h> #include <maths/common/CChecksum.h> #include <maths/common/CClusterer.h> #include <maths/common/CClustererStateSerialiser.h> #include <maths/common/CKMeansOnline1d.h> #include <maths/common/CMathsFuncs.h> #include <maths/common/CMathsFuncsForMatrixAndVectorTypes.h> #include <maths/common/CNormalMeanPrecConjugate.h> #include <maths/common/CPriorStateSerialiser.h> #include <maths/common/CRestoreParams.h> #include <maths/common/CSampling.h> #include <maths/common/CSetTools.h> #include <maths/common/CSolvers.h> #include <maths/common/CTools.h> #include <maths/common/CToolsDetail.h> #include <maths/common/MathsTypes.h> #include <maths/common/ProbabilityAggregators.h> #include <cmath> #include <limits> #include <numeric> #include <set> namespace ml { namespace maths { namespace common { namespace { using TSizeDoublePr = std::pair<std::size_t, double>; using TSizeDoublePr2Vec = core::CSmallVector<TSizeDoublePr, 2>; using TSizeDoublePr5Vec = core::CSmallVector<TSizeDoublePr, 5>; using TDoubleDoublePr = std::pair<double, double>; using TSizeSet = std::set<std::size_t>; using TDouble1Vec = core::CSmallVector<double, 1>; using maths_t::TDoubleWeightsAry; using maths_t::TDoubleWeightsAry1Vec; using TMeanAccumulator = CBasicStatistics::SSampleMean<double>::TAccumulator; //! \brief Wrapper to call the -log(c.d.f) of a prior object. class CMinusLogJointCdf { public: template<typename T> bool operator()(const T& prior, const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights, double& lowerBound, double& upperBound) const { return prior->minusLogJointCdf(samples, weights, lowerBound, upperBound); } }; //! \brief Wrapper to call the log(1 - c.d.f) of a prior object. class CMinusLogJointCdfComplement { public: template<typename T> bool operator()(const T& prior, const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights, double& lowerBound, double& upperBound) const { return prior->minusLogJointCdfComplement(samples, weights, lowerBound, upperBound); } }; //! \brief Wrapper of CMultimodalPrior::minusLogJointCdf function //! for use with our solver. class CLogCdf { public: using result_type = double; enum EStyle { E_Lower, E_Upper, E_Mean }; public: CLogCdf(EStyle style, const CMultimodalPrior& prior, const TDoubleWeightsAry& weights) : m_Style(style), m_Prior(&prior), m_Weights(1, weights), m_X(1u, 0.0) {} double operator()(double x) const { m_X[0] = x; double lowerBound, upperBound; if (m_Prior->minusLogJointCdf(m_X, m_Weights, lowerBound, upperBound) == false) { throw std::runtime_error("Unable to compute c.d.f. at " + core::CStringUtils::typeToString(x)); } switch (m_Style) { case E_Lower: return -lowerBound; case E_Upper: return -upperBound; case E_Mean: return -(lowerBound + upperBound) / 2.0; } return -(lowerBound + upperBound) / 2.0; } private: EStyle m_Style; const CMultimodalPrior* m_Prior; TDoubleWeightsAry1Vec m_Weights; //! Avoids creating the vector argument to minusLogJointCdf //! more than once. mutable TDouble1Vec m_X; }; const std::size_t MODE_SPLIT_NUMBER_SAMPLES(50u); const std::size_t MODE_MERGE_NUMBER_SAMPLES(25u); const core::TPersistenceTag CLUSTERER_TAG("a", "clusterer"); const core::TPersistenceTag SEED_PRIOR_TAG("b", "seed_prior"); const core::TPersistenceTag MODE_TAG("c", "mode"); const core::TPersistenceTag NUMBER_SAMPLES_TAG("d", "number_samples"); const core::TPersistenceTag DECAY_RATE_TAG("g", "decay_rate"); const std::string EMPTY_STRING; } //////// CMultimodalPrior Implementation //////// CMultimodalPrior::CMultimodalPrior(maths_t::EDataType dataType, const CClusterer1d& clusterer, const CPrior& seedPrior, double decayRate /*= 0.0*/) : CPrior(dataType, decayRate), m_Clusterer(clusterer.clone()), m_SeedPrior(seedPrior.clone()) { // Register the split and merge callbacks. m_Clusterer->splitFunc(CModeSplitCallback(*this)); m_Clusterer->mergeFunc(CModeMergeCallback(*this)); } CMultimodalPrior::CMultimodalPrior(maths_t::EDataType dataType, const TMeanVarAccumulatorVec& moments, double decayRate /*= 0.0*/) : CPrior(dataType, decayRate), m_SeedPrior( CNormalMeanPrecConjugate::nonInformativePrior(dataType, decayRate).clone()) { using TNormalVec = std::vector<CNormalMeanPrecConjugate>; TNormalVec normals; normals.reserve(moments.size()); for (const auto& moments_ : moments) { normals.emplace_back(dataType, moments_, decayRate); } m_Clusterer = std::make_unique<CKMeansOnline1d>(normals); m_Modes.reserve(normals.size()); for (std::size_t i = 0; i < normals.size(); ++i) { m_Modes.emplace_back(i, TPriorPtr(normals.back().clone())); } } CMultimodalPrior::CMultimodalPrior(maths_t::EDataType dataType, double decayRate, TPriorPtrVec& priors) : CPrior(dataType, decayRate) { m_Modes.reserve(priors.size()); for (std::size_t i = 0; i < priors.size(); ++i) { m_Modes.emplace_back(i, std::move(priors[i])); } // TODO - why does this constructor violate the invariant of m_SeedPrior != nullptr? } CMultimodalPrior::CMultimodalPrior(const SDistributionRestoreParams& params, core::CStateRestoreTraverser& traverser) : CPrior(params.s_DataType, params.s_DecayRate) { if (traverser.traverseSubLevel([&](auto& traverser_) { return this->acceptRestoreTraverser(params, traverser_); }) == false) { traverser.setBadState(); } } bool CMultimodalPrior::acceptRestoreTraverser(const SDistributionRestoreParams& params, core::CStateRestoreTraverser& traverser) { do { const std::string& name = traverser.name(); RESTORE_SETUP_TEARDOWN( DECAY_RATE_TAG, double decayRate{this->decayRate()}, m_Clusterer != nullptr && m_SeedPrior != nullptr && core::CStringUtils::stringToType(traverser.value(), decayRate), this->decayRate(decayRate)) RESTORE(CLUSTERER_TAG, traverser.traverseSubLevel( [&, serialiser = CClustererStateSerialiser{} ](auto& traverser_) { return serialiser(params, m_Clusterer, traverser_); })) RESTORE(SEED_PRIOR_TAG, traverser.traverseSubLevel( [&, serialiser = CPriorStateSerialiser{} ](auto& traverser_) { return serialiser(params, m_SeedPrior, traverser_); })) RESTORE_SETUP_TEARDOWN(MODE_TAG, TMode mode, traverser.traverseSubLevel([&](auto& traverser_) { return mode.acceptRestoreTraverser(params, traverser_); }), m_Modes.push_back(std::move(mode))) RESTORE_SETUP_TEARDOWN(NUMBER_SAMPLES_TAG, double numberSamples, core::CStringUtils::stringToType(traverser.value(), numberSamples), this->numberSamples(numberSamples)) } while (traverser.next()); if (m_Clusterer != nullptr) { // Register the split and merge callbacks. m_Clusterer->splitFunc(CModeSplitCallback(*this)); m_Clusterer->mergeFunc(CModeMergeCallback(*this)); } this->checkRestoredInvariants(); return true; } void CMultimodalPrior::checkRestoredInvariants() const { VIOLATES_INVARIANT_NO_EVALUATION(m_SeedPrior, ==, nullptr); } CMultimodalPrior::CMultimodalPrior(const CMultimodalPrior& other) : CPrior(other.dataType(), other.decayRate()), m_Clusterer(other.m_Clusterer != nullptr ? other.m_Clusterer->clone() : nullptr), m_SeedPrior(other.m_SeedPrior != nullptr ? other.m_SeedPrior->clone() : nullptr) { // Register the split and merge callbacks. if (m_Clusterer != nullptr) { m_Clusterer->splitFunc(CModeSplitCallback(*this)); m_Clusterer->mergeFunc(CModeMergeCallback(*this)); } // Clone all the modes up front so we can implement strong exception safety. TModeVec modes; modes.reserve(other.m_Modes.size()); for (const auto& mode : other.m_Modes) { modes.emplace_back(mode.s_Index, TPriorPtr(mode.s_Prior->clone())); } m_Modes.swap(modes); this->addSamples(other.numberSamples()); } CMultimodalPrior& CMultimodalPrior::operator=(const CMultimodalPrior& rhs) { if (this != &rhs) { CMultimodalPrior copy(rhs); this->swap(copy); } return *this; } void CMultimodalPrior::swap(CMultimodalPrior& other) { this->CPrior::swap(other); std::swap(m_Clusterer, other.m_Clusterer); // The call backs for split and merge should point to the // appropriate priors (we don't swap the "this" pointers // after all). So we need to refresh them after swapping. if (m_Clusterer != nullptr) { m_Clusterer->splitFunc(CModeSplitCallback(*this)); m_Clusterer->mergeFunc(CModeMergeCallback(*this)); } if (other.m_Clusterer != nullptr) { other.m_Clusterer->splitFunc(CModeSplitCallback(other)); other.m_Clusterer->mergeFunc(CModeMergeCallback(other)); } std::swap(m_SeedPrior, other.m_SeedPrior); m_Modes.swap(other.m_Modes); } CMultimodalPrior::EPrior CMultimodalPrior::type() const { return E_Multimodal; } CMultimodalPrior* CMultimodalPrior::clone() const { return new CMultimodalPrior(*this); } void CMultimodalPrior::dataType(maths_t::EDataType value) { this->CPrior::dataType(value); m_Clusterer->dataType(value); for (const auto& mode : m_Modes) { mode.s_Prior->dataType(value); } } void CMultimodalPrior::decayRate(double value) { this->CPrior::decayRate(value); m_Clusterer->decayRate(value); for (const auto& mode : m_Modes) { mode.s_Prior->decayRate(value); } m_SeedPrior->decayRate(value); } void CMultimodalPrior::setToNonInformative(double /*offset*/, double decayRate) { m_Clusterer->clear(); m_Modes.clear(); this->decayRate(decayRate); this->numberSamples(0.0); } bool CMultimodalPrior::needsOffset() const { for (const auto& mode : m_Modes) { if (mode.s_Prior->needsOffset()) { return true; } } return false; } double CMultimodalPrior::adjustOffset(const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights) { double result{0.0}; if (this->needsOffset()) { TSizeDoublePr2Vec clusters; for (std::size_t i = 0; i < samples.size(); ++i) { m_Clusterer->cluster(samples[i], clusters); for (const auto& cluster : clusters) { auto j = std::find_if(m_Modes.begin(), m_Modes.end(), CSetTools::CIndexInSet(cluster.first)); if (j != m_Modes.end()) { result += j->s_Prior->adjustOffset({samples[i]}, {weights[i]}); } } } } return result; } double CMultimodalPrior::offset() const { double offset{0.0}; for (const auto& mode : m_Modes) { offset = std::max(offset, mode.s_Prior->offset()); } return offset; } void CMultimodalPrior::addSamples(const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights) { if (samples.empty()) { return; } if (samples.size() != weights.size()) { LOG_ERROR(<< "Mismatch in samples '" << samples << "' and weights '" << weights << "'"); return; } this->adjustOffset(samples, weights); // This uses a clustering methodology (defined by m_Clusterer) // to assign each sample to a cluster. Each cluster has its own // mode in the distribution, which comprises a weight and prior // which are updated based on the assignment. // // The idea in separating clustering from modeling is to enable // different clustering methodologies to be used as appropriate // to the particular data set under consideration and based on // the runtime available. // // It is intended that approximate Bayesian schemes for multimode // distributions, such as the variational treatment of a mixture // of Gaussians, are implemented as separate priors. // // Flow through this function is as follows: // 1) Assign the samples to clusters. // 2) Find or create corresponding modes as necessary. // 3) Update the mode's weight and prior with the samples // assigned to it. // Declared outside the loop to minimize the number of times it // is initialized. TDouble1Vec sample(1); TDoubleWeightsAry1Vec weight(1); TSizeDoublePr2Vec clusters; try { bool hasSeasonalScale{this->isNonInformative() == false && maths_t::hasSeasonalVarianceScale(weights)}; double mean{hasSeasonalScale ? this->marginalLikelihoodMean() : 0.0}; for (std::size_t i = 0; i < samples.size(); ++i) { double x{samples[i]}; if (CMathsFuncs::isFinite(x) == false || CMathsFuncs::isFinite(weights[i]) == false) { LOG_ERROR(<< "Discarding sample = " << x << ", weights = " << weights[i]); continue; } if (hasSeasonalScale) { x = mean + (x - mean) / std::sqrt(maths_t::seasonalVarianceScale(weights[i])); } sample[0] = x; weight[0] = weights[i]; maths_t::setSeasonalVarianceScale(1.0, weight[0]); clusters.clear(); m_Clusterer->add(x, clusters, maths_t::count(weight[0])); double Z{std::accumulate( m_Modes.begin(), m_Modes.end(), maths_t::count(weight[0]), [](double sum, const TMode& mode) { return sum + mode.weight(); })}; double n{0.0}; for (const auto& cluster : clusters) { auto k = std::find_if(m_Modes.begin(), m_Modes.end(), CSetTools::CIndexInSet(cluster.first)); if (k == m_Modes.end()) { LOG_TRACE(<< "Creating mode with index " << cluster.first); m_Modes.emplace_back(cluster.first, m_SeedPrior); k = m_Modes.end() - 1; } maths_t::setCount(cluster.second, weight[0]); if (maths_t::isWinsorised(weight)) { double ww = maths_t::outlierWeight(weight[0]); double f = (k->weight() + cluster.second) / Z; maths_t::setOutlierWeight(std::max(1.0 - (1.0 - ww) / f, ww * f), weight[0]); } k->s_Prior->addSamples(sample, weight); n += maths_t::countForUpdate(weight[0]); } this->addSamples(n); } } catch (const std::exception& e) { LOG_ERROR(<< "Failed to update likelihood: " << e.what()); } } void CMultimodalPrior::propagateForwardsByTime(double time) { if (CMathsFuncs::isFinite(time) == false || time < 0.0) { LOG_ERROR(<< "Bad propagation time " << time); return; } if (this->isNonInformative()) { // Nothing to be done. return; } // We want to hold the probabilities constant. Since the i'th // probability: // p(i) = w(i) / Sum_j{ w(j) } // // where w(i) is its weight we can achieve this by multiplying // all weights by some factor f in the range [0, 1]. m_Clusterer->propagateForwardsByTime(time); for (const auto& mode : m_Modes) { mode.s_Prior->propagateForwardsByTime(time); } // Remove any mode which is non-informative. while (m_Modes.size() > 1) { // Calling remove with the mode's index triggers a callback // which also removes it from s_Modes, see CModeMergeCallback. auto i = std::find_if(m_Modes.begin(), m_Modes.end(), [](const auto& mode) { return mode.s_Prior->isNonInformative(); }); if (i == m_Modes.end() || m_Clusterer->remove(i->s_Index) == false) { break; } } this->numberSamples(this->numberSamples() * std::exp(-this->decayRate() * time)); LOG_TRACE(<< "numberSamples = " << this->numberSamples()); } TDoubleDoublePr CMultimodalPrior::marginalLikelihoodSupport() const { if (m_Modes.size() == 0) { return {-std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}; } if (m_Modes.size() == 1) { return m_Modes[0].s_Prior->marginalLikelihoodSupport(); } TDoubleDoublePr result{-std::numeric_limits<double>::max(), std::numeric_limits<double>::max()}; // We define this is as the union of the mode supports. for (const auto& mode : m_Modes) { TDoubleDoublePr s{mode.s_Prior->marginalLikelihoodSupport()}; result.first = std::min(result.first, s.first); result.second = std::max(result.second, s.second); } return result; } double CMultimodalPrior::marginalLikelihoodMean() const { if (m_Modes.size() == 0) { return 0.0; } if (m_Modes.size() == 1) { return m_Modes[0].s_Prior->marginalLikelihoodMean(); } // By linearity we have that: // Integral{ x * Sum_i{ w(i) * f(x | i) } } // = Sum_i{ w(i) * Integral{ x * f(x | i) } } // = Sum_i{ w(i) * mean(i) } TMeanAccumulator result; for (const auto& mode : m_Modes) { result.add(mode.s_Prior->marginalLikelihoodMean(), mode.weight()); } return CBasicStatistics::mean(result); } double CMultimodalPrior::nearestMarginalLikelihoodMean(double value) const { if (m_Modes.empty()) { return 0.0; } double mean{m_Modes[0].s_Prior->marginalLikelihoodMean()}; double distance{std::fabs(value - mean)}; double result{mean}; for (std::size_t i = 1; i < m_Modes.size(); ++i) { mean = m_Modes[i].s_Prior->marginalLikelihoodMean(); if (std::fabs(value - mean) < distance) { distance = std::fabs(value - mean); result = mean; } } return result; } double CMultimodalPrior::marginalLikelihoodMode(const TDoubleWeightsAry& weights) const { if (m_Modes.size() == 0) { return 0.0; } if (m_Modes.size() == 1) { return m_Modes[0].s_Prior->marginalLikelihoodMode(weights); } using TMaxAccumulator = CBasicStatistics::SMax<double>::TAccumulator; // We'll approximate this as the maximum likelihood mode (mode). double result{0.0}; double seasonalScale{std::sqrt(maths_t::seasonalVarianceScale(weights))}; double countVarianceScale{maths_t::countVarianceScale(weights)}; // Declared outside the loop to minimize number of times they // are created. TDouble1Vec distributionMode(1); TDoubleWeightsAry1Vec weight{maths_t::countVarianceScaleWeight(countVarianceScale)}; TMaxAccumulator maxLikelihood; for (const auto& mode : m_Modes) { double w{mode.weight()}; const auto& prior = mode.s_Prior; distributionMode[0] = prior->marginalLikelihoodMode(weight[0]); double likelihood; if ((prior->jointLogMarginalLikelihood(distributionMode, weight, likelihood) & (maths_t::E_FpFailed | maths_t::E_FpOverflowed)) != 0) { continue; } if (maxLikelihood.add(std::log(w) + likelihood)) { result = distributionMode[0]; } } if (maths_t::hasSeasonalVarianceScale(weights)) { double mean{this->marginalLikelihoodMean()}; result = mean + seasonalScale * (result - mean); } return result; } CMultimodalPrior::TDouble1Vec CMultimodalPrior::marginalLikelihoodModes(const TDoubleWeightsAry& weights) const { TDouble1Vec result(m_Modes.size()); for (std::size_t i = 0; i < m_Modes.size(); ++i) { result[i] = m_Modes[i].s_Prior->marginalLikelihoodMode(weights); } return result; } double CMultimodalPrior::marginalLikelihoodVariance(const TDoubleWeightsAry& weights) const { if (m_Modes.size() == 0) { return std::numeric_limits<double>::max(); } if (m_Modes.size() == 1) { return m_Modes[0].s_Prior->marginalLikelihoodVariance(weights); } // By linearity we have that: // Integral{ (x - m)^2 * Sum_i{ w(i) * f(x | i) } } // = Sum_i{ w(i) * (Integral{ x^2 * f(x | i) } - m^2) } // = Sum_i{ w(i) * ((mi^2 + vi) - m^2) } double varianceScale{maths_t::seasonalVarianceScale(weights) * maths_t::countVarianceScale(weights)}; double mean{this->marginalLikelihoodMean()}; TMeanAccumulator result; for (const auto& mode : m_Modes) { double w{mode.weight()}; double mm{mode.s_Prior->marginalLikelihoodMean()}; double mv{mode.s_Prior->marginalLikelihoodVariance()}; result.add((mm - mean) * (mm + mean) + mv, w); } return std::max(varianceScale * CBasicStatistics::mean(result), 0.0); } TDoubleDoublePr CMultimodalPrior::marginalLikelihoodConfidenceInterval(double percentage, const TDoubleWeightsAry& weights) const { TDoubleDoublePr support{this->marginalLikelihoodSupport()}; if (this->isNonInformative()) { return support; } if (m_Modes.size() == 1) { return m_Modes[0].s_Prior->marginalLikelihoodConfidenceInterval(percentage, weights); } percentage = CTools::truncate(percentage, 0.0, 100.0); if (percentage == 100.0) { return support; } using TMinMaxAccumulator = CBasicStatistics::CMinMax<double>; // The multimodal distribution confidence interval must lie between the // minimum and maximum of the corresponding modes' confidence intervals. TMinMaxAccumulator lower; TMinMaxAccumulator upper; for (const auto& mode : m_Modes) { auto interval = mode.s_Prior->marginalLikelihoodConfidenceInterval(percentage, weights); lower.add(interval.first); upper.add(interval.second); } percentage /= 100.0; double p1{maths_t::count(weights) * std::log((1.0 - percentage) / 2.0)}; double p2{maths_t::count(weights) * std::log((1.0 + percentage) / 2.0)}; double vs{maths_t::seasonalVarianceScale(weights)}; double mean{vs != 1.0 ? this->marginalLikelihoodMean() : 0.0}; CLogCdf fl{CLogCdf::E_Lower, *this, weights}; CLogCdf fu{CLogCdf::E_Upper, *this, weights}; auto computePercentile = [&](const CLogCdf& f, double p, const TMinMaxAccumulator& bounds) { auto fMinusp = [&f, p](double x) { return f(x) - p; }; double a{bounds.min() < mean ? mean - std::max(std::sqrt(vs), 1.0) * (mean - bounds.min()) : bounds.min()}; double b{bounds.max() > mean ? mean + std::max(std::sqrt(vs), 1.0) * (bounds.max() - mean) : bounds.max()}; double fa{fMinusp(a)}; double fb{fMinusp(b)}; // If we enter either of the following loops, since the c.d.f. is // monotonic, we know that the interval is either right or left of // the required root. We can therefore shift the original interval // since a and b are, respectively, upper and lower bounds for the // root. Each iteration also doubles the original interval width. // Note, we can only enter one or other of these loops so provided // we supply a few more than 10 iterations to the solver it'll have // enough iterations to comfortably converge. for (std::size_t i = 0; fa > 0.0 && i < 10; ++i) { std::tie(a, b) = std::make_pair(b - 3.0 * (b - a), a); std::tie(fa, fb) = std::make_pair(fMinusp(a), fa); } for (std::size_t i = 0; fb < 0.0 && i < 10; ++i) { std::tie(a, b) = std::make_pair(b, a + 3.0 * (b - a)); std::tie(fa, fb) = std::make_pair(fb, fMinusp(b)); } std::size_t maxIterations{25}; CEqualWithTolerance<double> equal{CToleranceTypes::E_AbsoluteTolerance, std::min(std::numeric_limits<double>::epsilon() * b, 1e-3 * p / std::max(fa, fb))}; double percentile; CSolvers::solve(a, b, fa, fb, fMinusp, maxIterations, equal, percentile); LOG_TRACE(<< "p1 = " << p << ", x = " << percentile << ", f(x) = " << f(percentile) << " brackets = [" << a << "," << b << "]"); return percentile; }; TDoubleDoublePr result{mean - std::sqrt(vs) * (mean - lower.min()), mean + std::sqrt(vs) * (upper.max() - mean)}; try { result.first = computePercentile(fl, p1, lower); result.second = computePercentile(fu, p2, upper); } catch (const std::exception& e) { LOG_WARN(<< "Unable to compute percentiles: " << e.what() << ", percentiles = [" << p1 << "," << p2 << "] " << ", vs = " << vs); } return result; } maths_t::EFloatingPointErrorStatus CMultimodalPrior::jointLogMarginalLikelihood(const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights, double& result) const { result = 0.0; if (samples.empty()) { LOG_ERROR(<< "Can't compute likelihood for empty sample set"); return maths_t::E_FpFailed; } if (samples.size() != weights.size()) { LOG_ERROR(<< "Mismatch in samples '" << samples << "' and weights '" << weights << "'"); return maths_t::E_FpFailed; } if (this->isNonInformative()) { // The non-informative likelihood is improper and effectively // zero everywhere. We use minus max double because // log(0) = HUGE_VALUE, which causes problems for Windows. // Calling code is notified when the calculation overflows // and should avoid taking the exponential since this will // underflow and pollute the floating point environment. This // may cause issues for some library function implementations // (see fe*exceptflag for more details). result = -std::numeric_limits<double>::max(); return maths_t::E_FpOverflowed; } if (m_Modes.size() == 1) { return m_Modes[0].s_Prior->jointLogMarginalLikelihood(samples, weights, result); } // The likelihood can be computed from the conditional likelihood // that a sample is from each mode. In particular, the likelihood // of a sample x is: // L(x) = Sum_m{ L(x | m) * p(m) } // // where, // L(x | m) is the likelihood the sample is from the m'th mode, // p(m) is the probability a sample is from the m'th mode. // // We compute the combined likelihood by taking the product of the // individual likelihoods. Note, this brushes over the fact that the // joint marginal likelihood that a collection of samples is from // the i'th mode is not just the product of the likelihoods that the // individual samples are from the i'th mode since we're integrating // over a prior. Really, we should compute likelihoods over all // possible assignments of the samples to the modes and use the fact // that: // P(a) = Product_i{ Sum_m{ p(m) * I{a(i) = m} } } // // where, // P(a) is the probability of a given assignment, // p(m) is the probability a sample is from the m'th mode, // I{.} is the indicator function. // // The approximation is increasingly accurate as the prior distribution // on each mode narrows. result = 0.0; // Declared outside the loop to minimize number of times it is created. TDouble1Vec sample(1); TSizeDoublePr5Vec modeLogLikelihoods; modeLogLikelihoods.reserve(m_Modes.size()); double mean{maths_t::hasSeasonalVarianceScale(weights) ? this->marginalLikelihoodMean() : 0.0}; TDoubleWeightsAry1Vec weight{TWeights::UNIT}; try { for (std::size_t i = 0; i < samples.size(); ++i) { double n{maths_t::countForUpdate(weights[i])}; double seasonalScale{std::sqrt(maths_t::seasonalVarianceScale(weights[i]))}; double logSeasonalScale{seasonalScale != 1.0 ? std::log(seasonalScale) : 0.0}; sample[0] = mean + (samples[i] - mean) / seasonalScale; maths_t::setCountVarianceScale(maths_t::countVarianceScale(weights[i]), weight[0]); // We re-normalize so that the maximum log likelihood is one // to avoid underflow. modeLogLikelihoods.clear(); double maxLogLikelihood{-std::numeric_limits<double>::max()}; for (std::size_t j = 0; j < m_Modes.size(); ++j) { double modeLogLikelihood; maths_t::EFloatingPointErrorStatus status{m_Modes[j].s_Prior->jointLogMarginalLikelihood( sample, weight, modeLogLikelihood)}; if ((status & maths_t::E_FpFailed) != 0) { // Logging handled at a lower level. return status; } if ((status & maths_t::E_FpOverflowed) == 0) { modeLogLikelihoods.emplace_back(j, modeLogLikelihood); maxLogLikelihood = std::max(maxLogLikelihood, modeLogLikelihood); } } if (modeLogLikelihoods.empty()) { // Technically, the marginal likelihood is zero here // so the log would be infinite. We use minus max // double because log(0) = HUGE_VALUE, which causes // problems for Windows. Calling code is notified // when the calculation overflows and should avoid // taking the exponential since this will underflow // and pollute the floating point environment. This // may cause issues for some library function // implementations (see fe*exceptflag for more details). result = -std::numeric_limits<double>::max(); return maths_t::E_FpOverflowed; } LOG_TRACE(<< "modeLogLikelihoods = " << modeLogLikelihoods); double sampleLikelihood{0.0}; double Z{0.0}; for (const auto& modeLogLikelihood : modeLogLikelihoods) { double w{m_Modes[modeLogLikelihood.first].weight()}; // Divide through by the largest value to avoid underflow. sampleLikelihood += w * std::exp(modeLogLikelihood.second - maxLogLikelihood); Z += w; } sampleLikelihood /= Z; double sampleLogLikelihood{n * (std::log(sampleLikelihood) + maxLogLikelihood)}; LOG_TRACE(<< "sample = " << sample << ", maxLogLikelihood = " << maxLogLikelihood << ", sampleLogLikelihood = " << sampleLogLikelihood); result += sampleLogLikelihood - n * logSeasonalScale; } } catch (const std::exception& e) { LOG_ERROR(<< "Failed to compute likelihood: " << e.what()); return maths_t::E_FpFailed; } maths_t::EFloatingPointErrorStatus status{CMathsFuncs::fpStatus(result)}; if ((status & maths_t::E_FpFailed) != 0) { LOG_ERROR(<< "Failed to compute likelihood (" << this->debugWeights() << ")"); LOG_ERROR(<< "samples = " << samples); LOG_ERROR(<< "weights = " << weights); } LOG_TRACE(<< "Joint log likelihood = " << result); return status; } void CMultimodalPrior::sampleMarginalLikelihood(std::size_t numberSamples, TDouble1Vec& samples) const { samples.clear(); if (numberSamples == 0 || this->numberSamples() == 0.0) { return; } samples.clear(); if (m_Modes.size() == 1) { m_Modes[0].s_Prior->sampleMarginalLikelihood(numberSamples, samples); return; } // We sample each mode according to its weight. TDoubleVec normalizedWeights; normalizedWeights.reserve(m_Modes.size()); double Z{0.0}; for (const auto& mode : m_Modes) { double weight{mode.weight()}; normalizedWeights.push_back(weight); Z += weight; } for (auto& weight : normalizedWeights) { weight /= Z; } CSampling::TSizeVec sampling; CSampling::weightedSample(numberSamples, normalizedWeights, sampling); LOG_TRACE(<< "normalizedWeights = " << normalizedWeights << ", sampling = " << sampling); if (sampling.size() != m_Modes.size()) { LOG_ERROR(<< "Failed to sample marginal likelihood"); return; } samples.reserve(numberSamples); TDouble1Vec modeSamples; for (std::size_t i = 0; i < m_Modes.size(); ++i) { m_Modes[i].s_Prior->sampleMarginalLikelihood(sampling[i], modeSamples); LOG_TRACE(<< "modeSamples = " << modeSamples); std::copy(modeSamples.begin(), modeSamples.end(), std::back_inserter(samples)); } LOG_TRACE(<< "samples = " << samples); } bool CMultimodalPrior::minusLogJointCdf(const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights, double& lowerBound, double& upperBound) const { return this->minusLogJointCdfImpl(CMinusLogJointCdf{}, samples, weights, lowerBound, upperBound); } bool CMultimodalPrior::minusLogJointCdfComplement(const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights, double& lowerBound, double& upperBound) const { return this->minusLogJointCdfImpl(CMinusLogJointCdfComplement{}, samples, weights, lowerBound, upperBound); } bool CMultimodalPrior::probabilityOfLessLikelySamples(maths_t::EProbabilityCalculation calculation, const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights, double& lowerBound, double& upperBound, maths_t::ETail& tail) const { lowerBound = 0.0; upperBound = 1.0; tail = maths_t::E_UndeterminedTail; if (samples.empty()) { LOG_ERROR(<< "Can't compute distribution for empty sample set"); return false; } if (this->isNonInformative()) { return true; } if (m_Modes.size() == 1) { return m_Modes[0].s_Prior->probabilityOfLessLikelySamples( calculation, samples, weights, lowerBound, upperBound, tail); } // Ideally we'd find the probability of the set of samples whose // total likelihood is less than or equal to that of the specified // samples, i.e. the probability of the set // R = { y | L(y) < L(x) } // // where, // x = {x(1), x(2), ..., x(n)} is the sample vector. // y is understood to be a vector quantity. // // This is not *trivially* related to the probability that the // probabilities of the sets // R(i) = { y | L(y) < L(x(i)) } // // since the joint conditional likelihood must be integrated over // priors for the parameters. However, we'll approximate this as // the joint probability (of a collection of standard normal R.Vs.) // having probabilities {P(R(i))}. This becomes increasingly accurate // as the prior distribution narrows. // // For the two sided calculation, we use the fact that the likelihood // function decreases monotonically away from the interval [a, b] // whose end points are the leftmost and rightmost modes' modes // since all component likelihoods decrease away from this interval. // // To evaluate the probability in the interval [a, b] we relax // the hard constraint that regions where f > f(x) contribute // zero probability. In particular, we note that we can write // the probability as: // P = Integral{ I(f(s) < f(x)) * f(s) }ds // // and that: // I(f(s) < f(x)) = lim_{k->inf}{ exp(-k * (f(s)/f(x) - 1)) // / (1 + exp(-k * (f(s)/f(x) - 1))) } // // We evaluate a smoother integral, i.e. smaller p, initially // to find out which regions contribute the most to P and then // re-evaluate those regions we need with higher resolution // using the fact that the maximum error in the approximation // of I(f(s) < f(x)) is 0.5. switch (calculation) { case maths_t::E_OneSidedBelow: if (this->minusLogJointCdf(samples, weights, upperBound, lowerBound) == false) { return false; } lowerBound = std::exp(-lowerBound); upperBound = std::exp(-upperBound); tail = maths_t::E_LeftTail; break; case maths_t::E_TwoSided: { static const double EPS{1000.0 * std::numeric_limits<double>::epsilon()}; static const std::size_t MAX_ITERATIONS{20}; CJointProbabilityOfLessLikelySamples lowerBoundCalculator; CJointProbabilityOfLessLikelySamples upperBoundCalculator; TDoubleDoublePr support{this->marginalLikelihoodSupport()}; support.first = (1.0 + (support.first > 0.0 ? EPS : -EPS)) * support.first; support.second = (1.0 + (support.first > 0.0 ? EPS : -EPS)) * support.second; bool hasSeasonalScale{maths_t::hasSeasonalVarianceScale(weights)}; double mean{hasSeasonalScale ? this->marginalLikelihoodMean() : 0.0}; double a{std::numeric_limits<double>::max()}; double b{-std::numeric_limits<double>::max()}; double Z{0.0}; for (const auto& mode : m_Modes) { double m{mode.s_Prior->marginalLikelihoodMode()}; a = std::min(a, m); b = std::max(b, m); Z += mode.weight(); } a = CTools::truncate(a, support.first, support.second); b = CTools::truncate(b, support.first, support.second); LOG_TRACE(<< "a = " << a << ", b = " << b << ", Z = " << Z); // Declared outside the loop to minimize the number of times // it is created. TDoubleWeightsAry1Vec weight(1); int tail_{0}; for (std::size_t i = 0; i < samples.size(); ++i) { double x{samples[i]}; weight[0] = weights[i]; if (CMathsFuncs::isNan(x) || CMathsFuncs::isNan(weight[0])) { continue; } if (hasSeasonalScale) { x = mean + (x - mean) / std::sqrt(maths_t::seasonalVarianceScale(weight[0])); maths_t::setSeasonalVarianceScale(1.0, weight[0]); } double fx; maths_t::EFloatingPointErrorStatus status = this->jointLogMarginalLikelihood({x}, weight, fx); if ((status & maths_t::E_FpFailed) != 0) { LOG_ERROR(<< "Unable to compute likelihood for " << x); return false; } if ((status & maths_t::E_FpOverflowed) != 0) { lowerBound = upperBound = 0.0; return true; } LOG_TRACE(<< "x = " << x << ", f(x) = " << fx); CPrior::CLogMarginalLikelihood logLikelihood{*this, weight}; CTools::CMixtureProbabilityOfLessLikelySample calculator{m_Modes.size(), x, fx, a, b}; for (const auto& mode : m_Modes) { double w{mode.weight() / Z}; double centre{mode.s_Prior->marginalLikelihoodMode(weight[0])}; double spread{ std::sqrt(mode.s_Prior->marginalLikelihoodVariance(weight[0]))}; calculator.addMode(w, centre, spread); tail_ = tail_ | (x < centre ? maths_t::E_LeftTail : maths_t::E_RightTail); } double sampleLowerBound{0.0}; double sampleUpperBound{0.0}; double lb; double ub; double xl; CEqualWithTolerance<double> lequal{CToleranceTypes::E_AbsoluteTolerance, EPS * a}; if (calculator.leftTail(logLikelihood, MAX_ITERATIONS, lequal, xl)) { this->minusLogJointCdf({xl}, weight, lb, ub); sampleLowerBound += std::exp(std::min(-lb, -ub)); sampleUpperBound += std::exp(std::max(-lb, -ub)); } else { this->minusLogJointCdf({xl}, weight, lb, ub); sampleUpperBound += std::exp(std::max(-lb, -ub)); } double xr; CEqualWithTolerance<double> requal{CToleranceTypes::E_AbsoluteTolerance, EPS * b}; if (calculator.rightTail(logLikelihood, MAX_ITERATIONS, requal, xr)) { this->minusLogJointCdfComplement({xr}, weight, lb, ub); sampleLowerBound += std::exp(std::min(-lb, -ub)); sampleUpperBound += std::exp(std::max(-lb, -ub)); } else { this->minusLogJointCdfComplement({xr}, weight, lb, ub); sampleUpperBound += std::exp(std::max(-lb, -ub)); } double p{0.0}; if (a < b) { p = calculator.calculate(logLikelihood, sampleLowerBound); } LOG_TRACE(<< "sampleLowerBound = " << sampleLowerBound << ", sampleUpperBound = " << sampleUpperBound << " p = " << p); lowerBoundCalculator.add(CTools::truncate(sampleLowerBound + p, 0.0, 1.0)); upperBoundCalculator.add(CTools::truncate(sampleUpperBound + p, 0.0, 1.0)); } if (lowerBoundCalculator.calculate(lowerBound) == false || upperBoundCalculator.calculate(upperBound) == false) { LOG_ERROR(<< "Couldn't compute probability of less likely samples:" << " " << lowerBoundCalculator << " " << upperBoundCalculator << " (samples = " << samples << ", weights = " << weights << ")"); return false; } tail = static_cast<maths_t::ETail>(tail_); } break; case maths_t::E_OneSidedAbove: if (this->minusLogJointCdfComplement(samples, weights, upperBound, lowerBound) == false) { return false; } lowerBound = std::exp(-lowerBound); upperBound = std::exp(-upperBound); tail = maths_t::E_RightTail; break; } return true; } bool CMultimodalPrior::isNonInformative() const { return m_Modes.empty() || (m_Modes.size() == 1 && m_Modes[0].s_Prior->isNonInformative()); } void CMultimodalPrior::print(const std::string& indent, std::string& result) const { result += "\n" + indent + "multimodal"; if (this->isNonInformative()) { result += " non-informative"; return; } double Z{0.0}; for (const auto& mode : m_Modes) { Z += mode.weight(); } result += ":"; for (const auto& mode : m_Modes) { double weight{mode.weight() / Z}; std::string indent_{indent + " weight " + core::CStringUtils::typeToStringPretty(weight) + " "}; mode.s_Prior->print(indent_, result); } } std::string CMultimodalPrior::printJointDensityFunction() const { return "Not supported"; } std::uint64_t CMultimodalPrior::checksum(std::uint64_t seed) const { seed = this->CPrior::checksum(seed); seed = CChecksum::calculate(seed, m_Clusterer); seed = CChecksum::calculate(seed, m_SeedPrior); return CChecksum::calculate(seed, m_Modes); } void CMultimodalPrior::debugMemoryUsage(const core::CMemoryUsage::TMemoryUsagePtr& mem) const { mem->setName("CMultimodalPrior"); core::memory_debug::dynamicSize("m_Clusterer", m_Clusterer, mem); core::memory_debug::dynamicSize("m_SeedPrior", m_SeedPrior, mem); core::memory_debug::dynamicSize("m_Modes", m_Modes, mem); } std::size_t CMultimodalPrior::memoryUsage() const { std::size_t mem = core::memory::dynamicSize(m_Clusterer); mem += core::memory::dynamicSize(m_SeedPrior); mem += core::memory::dynamicSize(m_Modes); return mem; } std::size_t CMultimodalPrior::staticSize() const { return sizeof(*this); } void CMultimodalPrior::acceptPersistInserter(core::CStatePersistInserter& inserter) const { inserter.insertLevel(CLUSTERER_TAG, [ this, serialiser = CClustererStateSerialiser{} ](auto& inserter_) { serialiser(*m_Clusterer, inserter_); }); inserter.insertLevel(SEED_PRIOR_TAG, [ this, serialiser = CPriorStateSerialiser{} ](auto& inserter_) { serialiser(*m_SeedPrior, inserter_); }); for (const auto& mode : m_Modes) { inserter.insertLevel(MODE_TAG, [&mode](auto& inserter_) { mode.acceptPersistInserter(inserter_); }); } inserter.insertValue(DECAY_RATE_TAG, this->decayRate(), core::CIEEE754::E_SinglePrecision); inserter.insertValue(NUMBER_SAMPLES_TAG, this->numberSamples(), core::CIEEE754::E_SinglePrecision); } std::size_t CMultimodalPrior::numberModes() const { return m_Modes.size(); } bool CMultimodalPrior::checkInvariants(const std::string& tag) const { bool result{true}; if (m_Modes.size() != m_Clusterer->numberClusters()) { LOG_ERROR(<< tag << "# modes = " << m_Modes.size() << ", # clusters = " << m_Clusterer->numberClusters()); result = false; } double numberSamples{this->numberSamples()}; double modeSamples{0.0}; for (const auto& mode : m_Modes) { if (m_Clusterer->hasCluster(mode.s_Index) == false) { LOG_ERROR(<< tag << "Expected cluster for mode = " << mode.s_Index); result = false; } modeSamples += mode.s_Prior->numberSamples(); } CEqualWithTolerance<double> equal{ CToleranceTypes::E_AbsoluteTolerance | CToleranceTypes::E_RelativeTolerance, 1e-3}; if (equal(modeSamples, numberSamples) == false) { LOG_ERROR(<< tag << "Sum mode samples = " << modeSamples << ", total samples = " << numberSamples); result = false; } return result; } bool CMultimodalPrior::participatesInModelSelection() const { return m_Modes.size() > 1; } double CMultimodalPrior::unmarginalizedParameters() const { return std::max(static_cast<double>(m_Modes.size()), 1.0) - 1.0; } template<typename CDF> bool CMultimodalPrior::minusLogJointCdfImpl(CDF minusLogCdf, const TDouble1Vec& samples, const TDoubleWeightsAry1Vec& weights, double& lowerBound, double& upperBound) const { lowerBound = upperBound = 0.0; if (samples.empty()) { LOG_ERROR(<< "Can't compute c.d.f. for empty sample set"); return false; } if (m_Modes.size() == 1) { return minusLogCdf(m_Modes[0].s_Prior, samples, weights, lowerBound, upperBound); } using TMinAccumulator = CBasicStatistics::COrderStatisticsStack<double, 1>; // The c.d.f. of the marginal likelihood is the weighted sum // of the c.d.fs of each mode since: // cdf(x) = Integral{ L(u) }du // = Integral{ Sum_m{ L(u | m) p(m) } }du // = Sum_m{ Integral{ L(u | m) ) p(m) }du } // Declared outside the loop to minimize the number of times // they are created. TDouble1Vec sample(1); TDoubleWeightsAry1Vec weight{TWeights::UNIT}; TDoubleVec modeLowerBounds; TDoubleVec modeUpperBounds; modeLowerBounds.reserve(m_Modes.size()); modeUpperBounds.reserve(m_Modes.size()); bool success = false; try { double mean{maths_t::hasSeasonalVarianceScale(weights) ? this->marginalLikelihoodMean() : 0.0}; for (std::size_t i = 0; i < samples.size(); ++i) { if (!CMathsFuncs::isFinite(samples[i]) || !CMathsFuncs::isFinite(weights[i])) { continue; } success = true; double n{maths_t::count(weights[i])}; double seasonalScale{std::sqrt(maths_t::seasonalVarianceScale(weights[i]))}; if (this->isNonInformative()) { lowerBound -= n * std::log(CTools::IMPROPER_CDF); upperBound -= n * std::log(CTools::IMPROPER_CDF); continue; } sample[0] = mean + (samples[i] - mean) / seasonalScale; maths_t::setCountVarianceScale(maths_t::countVarianceScale(weights[i]), weight[0]); // We re-normalize so that the maximum log c.d.f. is one // to avoid underflow. TMinAccumulator minLowerBound; TMinAccumulator minUpperBound; modeLowerBounds.clear(); modeUpperBounds.clear(); for (const auto& mode : m_Modes) { double modeLowerBound; double modeUpperBound; if (minusLogCdf(mode.s_Prior, sample, weight, modeLowerBound, modeUpperBound) == false) { return false; } minLowerBound.add(modeLowerBound); minUpperBound.add(modeUpperBound); modeLowerBounds.push_back(modeLowerBound); modeUpperBounds.push_back(modeUpperBound); } TMeanAccumulator sampleLowerBound; TMeanAccumulator sampleUpperBound; for (std::size_t j = 0; j < m_Modes.size(); ++j) { LOG_TRACE(<< "Mode -log(c.d.f.) = [" << modeLowerBounds[j] << "," << modeUpperBounds[j] << "]"); double w{m_Modes[j].weight()}; // Divide through by the largest value to avoid underflow. // Remember we are working with minus logs so the largest // value corresponds to the smallest log. sampleLowerBound.add(std::exp(-(modeLowerBounds[j] - minLowerBound[0])), w); sampleUpperBound.add(std::exp(-(modeUpperBounds[j] - minUpperBound[0])), w); } lowerBound += n * std::max(minLowerBound[0] - std::log(CBasicStatistics::mean(sampleLowerBound)), 0.0); upperBound += n * std::max(minUpperBound[0] - std::log(CBasicStatistics::mean(sampleUpperBound)), 0.0); LOG_TRACE(<< "sample = " << sample << ", sample -log(c.d.f.) = [" << sampleLowerBound << "," << sampleUpperBound << "]"); } } catch (const std::exception& e) { LOG_ERROR(<< "Failed to calculate c.d.f.: " << e.what()); return false; } if (!success) { LOG_ERROR(<< "Unable to compute c.d.f. (samples = " << samples << ", weights = " << weights << ")"); } LOG_TRACE(<< "Joint -log(c.d.f.) = [" << lowerBound << "," << upperBound << "]"); return success; } std::string CMultimodalPrior::debugWeights() const { return TMode::debugWeights(m_Modes); } ////////// CMultimodalPrior::CModeSplitCallback Implementation ////////// CMultimodalPrior::CModeSplitCallback::CModeSplitCallback(CMultimodalPrior& prior) : m_Prior(&prior) { } void CMultimodalPrior::CModeSplitCallback::operator()(std::size_t sourceIndex, std::size_t leftSplitIndex, std::size_t rightSplitIndex) const { LOG_TRACE(<< "Splitting mode with index " << sourceIndex); TModeVec& modes = m_Prior->m_Modes; // Remove the split mode. auto mode = std::find_if(modes.begin(), modes.end(), CSetTools::CIndexInSet(sourceIndex)); double numberSamples{mode != modes.end() ? mode->weight() : 0.0}; modes.erase(mode); double pLeft{m_Prior->m_Clusterer->probability(leftSplitIndex)}; double pRight{m_Prior->m_Clusterer->probability(rightSplitIndex)}; double Z = (pLeft + pRight); if (Z > 0.0) { pLeft /= Z; pRight /= Z; } LOG_TRACE(<< "# samples = " << numberSamples << ", pLeft = " << pLeft << ", pRight = " << pRight); // Create the child modes. LOG_TRACE(<< "Creating mode with index " << leftSplitIndex); modes.emplace_back(leftSplitIndex, TPriorPtr(m_Prior->m_SeedPrior->clone())); { TDoubleVec samples; if (m_Prior->m_Clusterer->sample(leftSplitIndex, MODE_SPLIT_NUMBER_SAMPLES, samples) == false) { LOG_ERROR(<< "Couldn't find cluster for " << leftSplitIndex); } LOG_TRACE(<< "samples = " << samples); double wl{pLeft * numberSamples}; double ws{std::min(wl, 4.0)}; double n{static_cast<double>(samples.size())}; LOG_TRACE(<< "# left = " << wl); TDoubleWeightsAry1Vec weights(samples.size(), maths_t::countWeight(ws / n)); modes.back().s_Prior->addSamples(samples, weights); if (wl > ws) { weights.assign(weights.size(), maths_t::countWeight((wl - ws) / n)); modes.back().s_Prior->addSamples(samples, weights); LOG_TRACE(<< modes.back().s_Prior->print()); } } LOG_TRACE(<< "Creating mode with index " << rightSplitIndex); modes.emplace_back(rightSplitIndex, TPriorPtr(m_Prior->m_SeedPrior->clone())); { TDoubleVec samples; if (m_Prior->m_Clusterer->sample(rightSplitIndex, MODE_SPLIT_NUMBER_SAMPLES, samples) == false) { LOG_ERROR(<< "Couldn't find cluster for " << rightSplitIndex); } LOG_TRACE(<< "samples = " << samples); double wr{pRight * numberSamples}; double ws{std::min(wr, 4.0)}; double n{static_cast<double>(samples.size())}; LOG_TRACE(<< "# right = " << wr); TDoubleWeightsAry1Vec weights(samples.size(), maths_t::countWeight(ws / n)); modes.back().s_Prior->addSamples(samples, weights); if (wr > ws) { weights.assign(weights.size(), maths_t::countWeight((wr - ws) / n)); modes.back().s_Prior->addSamples(samples, weights); LOG_TRACE(<< modes.back().s_Prior->print()); } } if (m_Prior->checkInvariants("SPLIT: ") == false) { LOG_ERROR(<< "# samples = " << numberSamples << ", # modes = " << modes.size() << ", pLeft = " << pLeft << ", pRight = " << pRight); } LOG_TRACE(<< "Split mode"); } ////////// CMultimodalPrior::CModeMergeCallback Implementation ////////// CMultimodalPrior::CModeMergeCallback::CModeMergeCallback(CMultimodalPrior& prior) : m_Prior(&prior) { } void CMultimodalPrior::CModeMergeCallback::operator()(std::size_t leftMergeIndex, std::size_t rightMergeIndex, std::size_t targetIndex) const { LOG_TRACE(<< "Merging modes with indices " << leftMergeIndex << " " << rightMergeIndex); TModeVec& modes{m_Prior->m_Modes}; // Create the new mode. TMode newMode(targetIndex, TPriorPtr(m_Prior->m_SeedPrior->clone())); double wl{0.0}; double wr{0.0}; double w{0.0}; std::size_t nl{0}; std::size_t nr{0}; TDouble1Vec samples; auto leftMode = std::find_if(modes.begin(), modes.end(), CSetTools::CIndexInSet(leftMergeIndex)); if (leftMode != modes.end()) { wl = leftMode->s_Prior->numberSamples(); w += wl; TDouble1Vec leftSamples; leftMode->s_Prior->sampleMarginalLikelihood(MODE_MERGE_NUMBER_SAMPLES, leftSamples); nl = leftSamples.size(); samples.insert(samples.end(), leftSamples.begin(), leftSamples.end()); } else { LOG_ERROR(<< "Couldn't find mode for " << leftMergeIndex); } auto rightMode = std::find_if(modes.begin(), modes.end(), CSetTools::CIndexInSet(rightMergeIndex)); if (rightMode != modes.end()) { wr = rightMode->s_Prior->numberSamples(); w += wr; TDouble1Vec rightSamples; rightMode->s_Prior->sampleMarginalLikelihood(MODE_MERGE_NUMBER_SAMPLES, rightSamples); nr = rightSamples.size(); samples.insert(samples.end(), rightSamples.begin(), rightSamples.end()); } else { LOG_ERROR(<< "Couldn't find mode for " << rightMergeIndex); } if (w > 0.0) { double nl_{static_cast<double>(nl)}; double nr_{static_cast<double>(nr)}; double Z{(nl_ * wl + nr_ * wr) / (nl_ + nr_)}; wl /= Z; wr /= Z; } LOG_TRACE(<< "samples = " << samples); LOG_TRACE(<< "w = " << w << ", wl = " << wl << ", wr = " << wr); double ws{std::min(w, 4.0)}; double n{static_cast<double>(samples.size())}; TDoubleWeightsAry1Vec weights; weights.reserve(samples.size()); weights.resize(nl, maths_t::countWeight(wl * ws / n)); weights.resize(nl + nr, maths_t::countWeight(wr * ws / n)); newMode.s_Prior->addSamples(samples, weights); if (w > ws) { weights.clear(); weights.resize(nl, maths_t::countWeight(wl * (w - ws) / n)); weights.resize(nl + nr, maths_t::countWeight(wr * (w - ws) / n)); newMode.s_Prior->addSamples(samples, weights); } // Remove the merged modes. TSizeSet mergedIndices; mergedIndices.insert(leftMergeIndex); mergedIndices.insert(rightMergeIndex); modes.erase(std::remove_if(modes.begin(), modes.end(), CSetTools::CIndexInSet(mergedIndices)), modes.end()); // Add the new mode. LOG_TRACE(<< "Creating mode with index " << targetIndex); modes.push_back(std::move(newMode)); m_Prior->checkInvariants("MERGE: "); LOG_TRACE(<< "Merged modes"); } } } }