diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml index 2eef075..89acd49 100644 --- a/.pre-commit-config.yaml +++ b/.pre-commit-config.yaml @@ -56,8 +56,8 @@ repos: hooks: - id: check-pragma-once name: Verify #pragma once in headers - entry: scripts/check-pragma-once.sh - language: script + entry: python scripts/check-pragma-once.py + language: python files: \.(h|hpp|hxx)$ pass_filenames: true diff --git a/CMakeLists.txt b/CMakeLists.txt index ded76de..1d41c3d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -343,6 +343,17 @@ include_directories(include) # Platform-specific compiler flags # NOTE: MSVC flags are applied per-target (see after add_library below) so that # FetchContent dependencies like GTest are not compiled with /W4. +# +# /utf-8 is the exception: it specifies source and execution character sets as +# UTF-8 globally. This suppresses C4566 warnings that occur when source files +# contain non-ASCII Unicode characters (Greek letters, math symbols, etc.) on +# Windows systems with a non-UTF-8 ANSI code page (e.g. CP1252). It is safe +# to apply to all targets including GTest because it only affects how MSVC +# interprets source encoding — it has no effect on ASCII-only code. +if(MSVC) + add_compile_options(/utf-8) +endif() + if(CMAKE_CXX_COMPILER_ID MATCHES "GNU|Clang") set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -Wall -Wextra -Wpedantic -Wpointer-arith") # Note: -Wno-deprecated-declarations and -Wno-defaulted-function-deleted removed @@ -508,6 +519,7 @@ set(LIBHMM_SOURCES src/common/common.cpp src/common/string_tokenizer.cpp src/common/numerical_stability.cpp + src/common/weighted_stats.cpp src/performance/transcendental_kernels.cpp src/distributions/distribution_base.cpp src/distributions/discrete_distribution.cpp diff --git a/include/libhmm/calculators/calculator.h b/include/libhmm/calculators/calculator.h index 5c7b2e2..5149bf6 100755 --- a/include/libhmm/calculators/calculator.h +++ b/include/libhmm/calculators/calculator.h @@ -72,6 +72,24 @@ class Calculator { * @param observations The new observation set */ void setObservations(const ObservationSet &observations) { observations_ = observations; } + + /// Fills logTrans (row-major) and logTransT (column-major/transposed) from the + /// HMM transition matrix. Resizes both matrices to numStates×numStates. + static void precompute_log_transitions(const Hmm &hmm, std::size_t numStates, Matrix &logTrans, + Matrix &logTransT) noexcept { + const Matrix &trans = hmm.getTrans(); + logTrans.resize(numStates, numStates); + logTransT.resize(numStates, numStates); + for (std::size_t i = 0; i < numStates; ++i) { + for (std::size_t j = 0; j < numStates; ++j) { + const double a = trans(i, j); + const double logA = + (a > 0.0) ? std::log(a) : -std::numeric_limits::infinity(); + logTrans(i, j) = logA; + logTransT(j, i) = logA; + } + } + } }; //class Calculator } // namespace libhmm diff --git a/include/libhmm/distributions/beta_distribution.h b/include/libhmm/distributions/beta_distribution.h index b9be62d..218fbea 100644 --- a/include/libhmm/distributions/beta_distribution.h +++ b/include/libhmm/distributions/beta_distribution.h @@ -99,42 +99,10 @@ class BetaDistribution : public DistributionBase { updateCache(); } - BetaDistribution(const BetaDistribution &other) - : DistributionBase{other}, alpha_{other.alpha_}, beta_{other.beta_}, - logBeta_{other.logBeta_}, alphaMinus1_{other.alphaMinus1_}, - betaMinus1_{other.betaMinus1_}, invBeta_{other.invBeta_} {} - - BetaDistribution &operator=(const BetaDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - alpha_ = other.alpha_; - beta_ = other.beta_; - logBeta_ = other.logBeta_; - alphaMinus1_ = other.alphaMinus1_; - betaMinus1_ = other.betaMinus1_; - invBeta_ = other.invBeta_; - } - return *this; - } - - BetaDistribution(BetaDistribution &&other) noexcept - : DistributionBase{std::move(other)}, alpha_{other.alpha_}, beta_{other.beta_}, - logBeta_{other.logBeta_}, alphaMinus1_{other.alphaMinus1_}, - betaMinus1_{other.betaMinus1_}, invBeta_{other.invBeta_} {} - - BetaDistribution &operator=(BetaDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - alpha_ = other.alpha_; - beta_ = other.beta_; - logBeta_ = other.logBeta_; - alphaMinus1_ = other.alphaMinus1_; - betaMinus1_ = other.betaMinus1_; - invBeta_ = other.invBeta_; - } - return *this; - } - + BetaDistribution(const BetaDistribution &other) = default; + BetaDistribution &operator=(const BetaDistribution &other) = default; + BetaDistribution(BetaDistribution &&other) noexcept = default; + BetaDistribution &operator=(BetaDistribution &&other) noexcept = default; ~BetaDistribution() override = default; /** diff --git a/include/libhmm/distributions/binomial_distribution.h b/include/libhmm/distributions/binomial_distribution.h index 8dced3f..77122ad 100644 --- a/include/libhmm/distributions/binomial_distribution.h +++ b/include/libhmm/distributions/binomial_distribution.h @@ -67,40 +67,10 @@ class BinomialDistribution : public DistributionBase { updateCache(); } - BinomialDistribution(const BinomialDistribution &other) - : DistributionBase{other}, n_{other.n_}, p_{other.p_}, - logFactorialCache_{other.logFactorialCache_}, logP_{other.logP_}, - log1MinusP_{other.log1MinusP_} {} - - BinomialDistribution &operator=(const BinomialDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - n_ = other.n_; - p_ = other.p_; - logFactorialCache_ = other.logFactorialCache_; - logP_ = other.logP_; - log1MinusP_ = other.log1MinusP_; - } - return *this; - } - - BinomialDistribution(BinomialDistribution &&other) noexcept - : DistributionBase{std::move(other)}, n_{other.n_}, p_{other.p_}, - logFactorialCache_{std::move(other.logFactorialCache_)}, logP_{other.logP_}, - log1MinusP_{other.log1MinusP_} {} - - BinomialDistribution &operator=(BinomialDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - n_ = other.n_; - p_ = other.p_; - logFactorialCache_ = std::move(other.logFactorialCache_); - logP_ = other.logP_; - log1MinusP_ = other.log1MinusP_; - } - return *this; - } - + BinomialDistribution(const BinomialDistribution &other) = default; + BinomialDistribution &operator=(const BinomialDistribution &other) = default; + BinomialDistribution(BinomialDistribution &&other) noexcept = default; + BinomialDistribution &operator=(BinomialDistribution &&other) noexcept = default; ~BinomialDistribution() override = default; [[nodiscard]] double getProbability(double value) const override; diff --git a/include/libhmm/distributions/chi_squared_distribution.h b/include/libhmm/distributions/chi_squared_distribution.h index 1db1414..78bcac8 100644 --- a/include/libhmm/distributions/chi_squared_distribution.h +++ b/include/libhmm/distributions/chi_squared_distribution.h @@ -78,40 +78,10 @@ class ChiSquaredDistribution : public DistributionBase { updateCache(); } - ChiSquaredDistribution(const ChiSquaredDistribution &other) - : DistributionBase{other}, degrees_of_freedom_{other.degrees_of_freedom_}, - cached_log_gamma_half_k_{other.cached_log_gamma_half_k_}, - cached_log_normalization_{other.cached_log_normalization_}, - cached_half_k_minus_one_{other.cached_half_k_minus_one_} {} - - ChiSquaredDistribution &operator=(const ChiSquaredDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - degrees_of_freedom_ = other.degrees_of_freedom_; - cached_log_gamma_half_k_ = other.cached_log_gamma_half_k_; - cached_log_normalization_ = other.cached_log_normalization_; - cached_half_k_minus_one_ = other.cached_half_k_minus_one_; - } - return *this; - } - - ChiSquaredDistribution(ChiSquaredDistribution &&other) noexcept - : DistributionBase{std::move(other)}, degrees_of_freedom_{other.degrees_of_freedom_}, - cached_log_gamma_half_k_{other.cached_log_gamma_half_k_}, - cached_log_normalization_{other.cached_log_normalization_}, - cached_half_k_minus_one_{other.cached_half_k_minus_one_} {} - - ChiSquaredDistribution &operator=(ChiSquaredDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - degrees_of_freedom_ = other.degrees_of_freedom_; - cached_log_gamma_half_k_ = other.cached_log_gamma_half_k_; - cached_log_normalization_ = other.cached_log_normalization_; - cached_half_k_minus_one_ = other.cached_half_k_minus_one_; - } - return *this; - } - + ChiSquaredDistribution(const ChiSquaredDistribution &other) = default; + ChiSquaredDistribution &operator=(const ChiSquaredDistribution &other) = default; + ChiSquaredDistribution(ChiSquaredDistribution &&other) noexcept = default; + ChiSquaredDistribution &operator=(ChiSquaredDistribution &&other) noexcept = default; ~ChiSquaredDistribution() override = default; /** @@ -193,14 +163,6 @@ class ChiSquaredDistribution : public DistributionBase { */ double getMode() const noexcept { return std::max(0.0, degrees_of_freedom_ - 2.0); } - /** - * Create distribution from string representation - * @param str String representation - * @return ChiSquaredDistribution object - * @throws std::invalid_argument if string format is invalid - */ - static ChiSquaredDistribution fromString(const std::string &str); - /** * Equality comparison operator * @param other Other distribution to compare with diff --git a/include/libhmm/distributions/discrete_distribution.h b/include/libhmm/distributions/discrete_distribution.h index 54f69d1..5c2265f 100755 --- a/include/libhmm/distributions/discrete_distribution.h +++ b/include/libhmm/distributions/discrete_distribution.h @@ -127,52 +127,10 @@ class DiscreteDistribution : public DistributionBase { reset(); } - DiscreteDistribution(const DiscreteDistribution &other) - : DistributionBase{other}, numSymbols_{other.numSymbols_}, pdf_{other.pdf_}, - cachedSum_{other.cachedSum_}, cachedEntropy_{other.cachedEntropy_}, - cachedLogProbs_{other.cachedLogProbs_}, cachedCDF_{other.cachedCDF_}, - cachedMode_{other.cachedMode_}, cachedMaxProb_{other.cachedMaxProb_}, - nonZeroIndices_{other.nonZeroIndices_} {} - - DiscreteDistribution &operator=(const DiscreteDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - numSymbols_ = other.numSymbols_; - pdf_ = other.pdf_; - cachedSum_ = other.cachedSum_; - cachedEntropy_ = other.cachedEntropy_; - cachedLogProbs_ = other.cachedLogProbs_; - cachedCDF_ = other.cachedCDF_; - cachedMode_ = other.cachedMode_; - cachedMaxProb_ = other.cachedMaxProb_; - nonZeroIndices_ = other.nonZeroIndices_; - } - return *this; - } - - DiscreteDistribution(DiscreteDistribution &&other) noexcept - : DistributionBase{std::move(other)}, numSymbols_{other.numSymbols_}, - pdf_{std::move(other.pdf_)}, cachedSum_{other.cachedSum_}, - cachedEntropy_{other.cachedEntropy_}, cachedLogProbs_{std::move(other.cachedLogProbs_)}, - cachedCDF_{std::move(other.cachedCDF_)}, cachedMode_{other.cachedMode_}, - cachedMaxProb_{other.cachedMaxProb_}, nonZeroIndices_{std::move(other.nonZeroIndices_)} {} - - DiscreteDistribution &operator=(DiscreteDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - numSymbols_ = other.numSymbols_; - pdf_ = std::move(other.pdf_); - cachedSum_ = other.cachedSum_; - cachedEntropy_ = other.cachedEntropy_; - cachedLogProbs_ = std::move(other.cachedLogProbs_); - cachedCDF_ = std::move(other.cachedCDF_); - cachedMode_ = other.cachedMode_; - cachedMaxProb_ = other.cachedMaxProb_; - nonZeroIndices_ = std::move(other.nonZeroIndices_); - } - return *this; - } - + DiscreteDistribution(const DiscreteDistribution &other) = default; + DiscreteDistribution &operator=(const DiscreteDistribution &other) = default; + DiscreteDistribution(DiscreteDistribution &&other) noexcept = default; + DiscreteDistribution &operator=(DiscreteDistribution &&other) noexcept = default; ~DiscreteDistribution() override = default; /** diff --git a/include/libhmm/distributions/distribution_io_utils.h b/include/libhmm/distributions/distribution_io_utils.h new file mode 100644 index 0000000..d0ef347 --- /dev/null +++ b/include/libhmm/distributions/distribution_io_utils.h @@ -0,0 +1,30 @@ +#pragma once + +#include +#include +#include +#include +#include + +namespace libhmm::detail { + +/// Splits a "name=value" token at the first '=', strips whitespace from both +/// parts, and returns {name, value}. +/// Throws std::invalid_argument (with context in the message) if no '=' found. +[[nodiscard]] inline std::pair +parse_named_param(const std::string ¶m, const std::string &context) { + const auto eq = param.find('='); + if (eq == std::string::npos) + throw std::invalid_argument("Invalid " + context + " parameter format"); + std::string name = param.substr(0, eq); + std::string value = param.substr(eq + 1); + const auto trim = [](std::string &s) { + s.erase(std::remove_if(s.begin(), s.end(), [](unsigned char c) { return std::isspace(c); }), + s.end()); + }; + trim(name); + trim(value); + return {std::move(name), std::move(value)}; +} + +} // namespace libhmm::detail diff --git a/include/libhmm/distributions/exponential_distribution.h b/include/libhmm/distributions/exponential_distribution.h index 2abbe25..5c68010 100755 --- a/include/libhmm/distributions/exponential_distribution.h +++ b/include/libhmm/distributions/exponential_distribution.h @@ -88,43 +88,10 @@ class ExponentialDistribution : public DistributionBase { updateCache(); } - /** - * Copy constructor - */ - ExponentialDistribution(const ExponentialDistribution &other) - : DistributionBase{other}, lambda_{other.lambda_}, logLambda_{other.logLambda_}, - invLambda_{other.invLambda_}, negLambda_{other.negLambda_}, - invLambdaSquared_{other.invLambdaSquared_} {} - - ExponentialDistribution &operator=(const ExponentialDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - lambda_ = other.lambda_; - logLambda_ = other.logLambda_; - invLambda_ = other.invLambda_; - negLambda_ = other.negLambda_; - invLambdaSquared_ = other.invLambdaSquared_; - } - return *this; - } - - ExponentialDistribution(ExponentialDistribution &&other) noexcept - : DistributionBase{std::move(other)}, lambda_{other.lambda_}, logLambda_{other.logLambda_}, - invLambda_{other.invLambda_}, negLambda_{other.negLambda_}, - invLambdaSquared_{other.invLambdaSquared_} {} - - ExponentialDistribution &operator=(ExponentialDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - lambda_ = other.lambda_; - logLambda_ = other.logLambda_; - invLambda_ = other.invLambda_; - negLambda_ = other.negLambda_; - invLambdaSquared_ = other.invLambdaSquared_; - } - return *this; - } - + ExponentialDistribution(const ExponentialDistribution &other) = default; + ExponentialDistribution &operator=(const ExponentialDistribution &other) = default; + ExponentialDistribution(ExponentialDistribution &&other) noexcept = default; + ExponentialDistribution &operator=(ExponentialDistribution &&other) noexcept = default; ~ExponentialDistribution() override = default; /** diff --git a/include/libhmm/distributions/gamma_distribution.h b/include/libhmm/distributions/gamma_distribution.h index 3474aa1..80d8e95 100755 --- a/include/libhmm/distributions/gamma_distribution.h +++ b/include/libhmm/distributions/gamma_distribution.h @@ -104,38 +104,10 @@ class GammaDistribution : public DistributionBase { updateCache(); } - GammaDistribution(const GammaDistribution &other) - : DistributionBase{other}, k_{other.k_}, theta_{other.theta_}, logGammaK_{other.logGammaK_}, - kLogTheta_{other.kLogTheta_}, kMinus1_{other.kMinus1_} {} - - GammaDistribution &operator=(const GammaDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - k_ = other.k_; - theta_ = other.theta_; - logGammaK_ = other.logGammaK_; - kLogTheta_ = other.kLogTheta_; - kMinus1_ = other.kMinus1_; - } - return *this; - } - - GammaDistribution(GammaDistribution &&other) noexcept - : DistributionBase{std::move(other)}, k_{other.k_}, theta_{other.theta_}, - logGammaK_{other.logGammaK_}, kLogTheta_{other.kLogTheta_}, kMinus1_{other.kMinus1_} {} - - GammaDistribution &operator=(GammaDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - k_ = other.k_; - theta_ = other.theta_; - logGammaK_ = other.logGammaK_; - kLogTheta_ = other.kLogTheta_; - kMinus1_ = other.kMinus1_; - } - return *this; - } - + GammaDistribution(const GammaDistribution &other) = default; + GammaDistribution &operator=(const GammaDistribution &other) = default; + GammaDistribution(GammaDistribution &&other) noexcept = default; + GammaDistribution &operator=(GammaDistribution &&other) noexcept = default; ~GammaDistribution() override = default; /** @@ -174,6 +146,13 @@ class GammaDistribution : public DistributionBase { void fit(std::span data) override; /** Weighted MOM: θ = weighted_var / weighted_mean, k = weighted_mean² / weighted_var. */ void fit(std::span data, std::span weights) override; + +private: + /// Validates mean/var and applies the MOM estimates. Calls reset() on any + /// failure; otherwise assigns theta_, k_ and invalidates the cache. + void apply_fit_params(double mean, double var); + +public: /** Returns false — Gamma is a continuous distribution. */ [[nodiscard]] bool isDiscrete() const noexcept override { return false; } diff --git a/include/libhmm/distributions/gaussian_distribution.h b/include/libhmm/distributions/gaussian_distribution.h index 9aafafb..a2ecac4 100755 --- a/include/libhmm/distributions/gaussian_distribution.h +++ b/include/libhmm/distributions/gaussian_distribution.h @@ -93,51 +93,10 @@ class GaussianDistribution : public DistributionBase { updateCache(); } - // Rule of Five: atomic in base requires explicit copy/move. - // Parameter data is trivially copyable; base handles the atomic. - GaussianDistribution(const GaussianDistribution &other) - : DistributionBase{other}, mean_{other.mean_}, standardDeviation_{other.standardDeviation_}, - normalizationConstant_{other.normalizationConstant_}, - negHalfSigmaSquaredInv_{other.negHalfSigmaSquaredInv_}, - logStandardDeviation_{other.logStandardDeviation_}, sigmaSqrt2_{other.sigmaSqrt2_}, - invStandardDeviation_{other.invStandardDeviation_} {} - - GaussianDistribution &operator=(const GaussianDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - mean_ = other.mean_; - standardDeviation_ = other.standardDeviation_; - normalizationConstant_ = other.normalizationConstant_; - negHalfSigmaSquaredInv_ = other.negHalfSigmaSquaredInv_; - logStandardDeviation_ = other.logStandardDeviation_; - sigmaSqrt2_ = other.sigmaSqrt2_; - invStandardDeviation_ = other.invStandardDeviation_; - } - return *this; - } - - GaussianDistribution(GaussianDistribution &&other) noexcept - : DistributionBase{std::move(other)}, mean_{other.mean_}, - standardDeviation_{other.standardDeviation_}, - normalizationConstant_{other.normalizationConstant_}, - negHalfSigmaSquaredInv_{other.negHalfSigmaSquaredInv_}, - logStandardDeviation_{other.logStandardDeviation_}, sigmaSqrt2_{other.sigmaSqrt2_}, - invStandardDeviation_{other.invStandardDeviation_} {} - - GaussianDistribution &operator=(GaussianDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - mean_ = other.mean_; - standardDeviation_ = other.standardDeviation_; - normalizationConstant_ = other.normalizationConstant_; - negHalfSigmaSquaredInv_ = other.negHalfSigmaSquaredInv_; - logStandardDeviation_ = other.logStandardDeviation_; - sigmaSqrt2_ = other.sigmaSqrt2_; - invStandardDeviation_ = other.invStandardDeviation_; - } - return *this; - } - + GaussianDistribution(const GaussianDistribution &other) = default; + GaussianDistribution &operator=(const GaussianDistribution &other) = default; + GaussianDistribution(GaussianDistribution &&other) noexcept = default; + GaussianDistribution &operator=(GaussianDistribution &&other) noexcept = default; ~GaussianDistribution() override = default; /** diff --git a/include/libhmm/distributions/log_normal_distribution.h b/include/libhmm/distributions/log_normal_distribution.h index c4a125e..631b459 100755 --- a/include/libhmm/distributions/log_normal_distribution.h +++ b/include/libhmm/distributions/log_normal_distribution.h @@ -90,51 +90,10 @@ class LogNormalDistribution : public DistributionBase { updateCache(); } - /** - * Copy constructor - */ - LogNormalDistribution(const LogNormalDistribution &other) - : DistributionBase{other}, mean_{other.mean_}, standardDeviation_{other.standardDeviation_}, - logNormalizationConstant_{other.logNormalizationConstant_}, - negHalfSigmaSquaredInv_{other.negHalfSigmaSquaredInv_} {} - - /** - * Copy assignment operator - */ - LogNormalDistribution &operator=(const LogNormalDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - mean_ = other.mean_; - standardDeviation_ = other.standardDeviation_; - logNormalizationConstant_ = other.logNormalizationConstant_; - negHalfSigmaSquaredInv_ = other.negHalfSigmaSquaredInv_; - } - return *this; - } - - /** - * Move constructor - */ - LogNormalDistribution(LogNormalDistribution &&other) noexcept - : DistributionBase{std::move(other)}, mean_{other.mean_}, - standardDeviation_{other.standardDeviation_}, - logNormalizationConstant_{other.logNormalizationConstant_}, - negHalfSigmaSquaredInv_{other.negHalfSigmaSquaredInv_} {} - - /** - * Move assignment operator - */ - LogNormalDistribution &operator=(LogNormalDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - mean_ = other.mean_; - standardDeviation_ = other.standardDeviation_; - logNormalizationConstant_ = other.logNormalizationConstant_; - negHalfSigmaSquaredInv_ = other.negHalfSigmaSquaredInv_; - } - return *this; - } - + LogNormalDistribution(const LogNormalDistribution &other) = default; + LogNormalDistribution &operator=(const LogNormalDistribution &other) = default; + LogNormalDistribution(LogNormalDistribution &&other) noexcept = default; + LogNormalDistribution &operator=(LogNormalDistribution &&other) noexcept = default; ~LogNormalDistribution() override = default; /** diff --git a/include/libhmm/distributions/negative_binomial_distribution.h b/include/libhmm/distributions/negative_binomial_distribution.h index e793e4f..29147e3 100644 --- a/include/libhmm/distributions/negative_binomial_distribution.h +++ b/include/libhmm/distributions/negative_binomial_distribution.h @@ -114,57 +114,11 @@ class NegativeBinomialDistribution : public DistributionBase { updateCache(); } - /** - * Copy constructor - */ - NegativeBinomialDistribution(const NegativeBinomialDistribution &other) - : DistributionBase{other}, r_{other.r_}, p_{other.p_}, logP_{other.logP_}, - log1MinusP_{other.log1MinusP_}, logGammaR_{other.logGammaR_}, - logFactorialCache_{other.logFactorialCache_} {} - - /** - * Copy assignment operator - */ - NegativeBinomialDistribution &operator=(const NegativeBinomialDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - r_ = other.r_; - p_ = other.p_; - logP_ = other.logP_; - log1MinusP_ = other.log1MinusP_; - logGammaR_ = other.logGammaR_; - logFactorialCache_ = other.logFactorialCache_; - } - return *this; - } - - /** - * Move constructor - */ - NegativeBinomialDistribution(NegativeBinomialDistribution &&other) noexcept - : DistributionBase{std::move(other)}, r_{other.r_}, p_{other.p_}, logP_{other.logP_}, - log1MinusP_{other.log1MinusP_}, logGammaR_{other.logGammaR_}, - logFactorialCache_{std::move(other.logFactorialCache_)} {} - - /** - * Move assignment operator - */ - NegativeBinomialDistribution &operator=(NegativeBinomialDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - r_ = other.r_; - p_ = other.p_; - logP_ = other.logP_; - log1MinusP_ = other.log1MinusP_; - logGammaR_ = other.logGammaR_; - logFactorialCache_ = std::move(other.logFactorialCache_); - } - return *this; - } - - /** - * Destructor - explicitly defaulted to satisfy Rule of Five - */ + NegativeBinomialDistribution(const NegativeBinomialDistribution &other) = default; + NegativeBinomialDistribution &operator=(const NegativeBinomialDistribution &other) = default; + NegativeBinomialDistribution(NegativeBinomialDistribution &&other) noexcept = default; + NegativeBinomialDistribution & + operator=(NegativeBinomialDistribution &&other) noexcept = default; ~NegativeBinomialDistribution() override = default; /** diff --git a/include/libhmm/distributions/pareto_distribution.h b/include/libhmm/distributions/pareto_distribution.h index 2baba0c..75c7dfa 100755 --- a/include/libhmm/distributions/pareto_distribution.h +++ b/include/libhmm/distributions/pareto_distribution.h @@ -117,58 +117,10 @@ class ParetoDistribution : public DistributionBase { updateCache(); } - /** - * Copy constructor - */ - ParetoDistribution(const ParetoDistribution &other) - : DistributionBase{other}, k_{other.k_}, xm_{other.xm_}, logK_{other.logK_}, - kLogXm_{other.kLogXm_}, kPlus1_{other.kPlus1_}, kXmPowK_{other.kXmPowK_}, - negK_{other.negK_}, logXm_{other.logXm_} {} - - /** - * Copy assignment operator - */ - ParetoDistribution &operator=(const ParetoDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - k_ = other.k_; - xm_ = other.xm_; - logK_ = other.logK_; - kLogXm_ = other.kLogXm_; - kPlus1_ = other.kPlus1_; - kXmPowK_ = other.kXmPowK_; - negK_ = other.negK_; - logXm_ = other.logXm_; - } - return *this; - } - - /** - * Move constructor - */ - ParetoDistribution(ParetoDistribution &&other) noexcept - : DistributionBase{std::move(other)}, k_{other.k_}, xm_{other.xm_}, logK_{other.logK_}, - kLogXm_{other.kLogXm_}, kPlus1_{other.kPlus1_}, kXmPowK_{other.kXmPowK_}, - negK_{other.negK_}, logXm_{other.logXm_} {} - - /** - * Move assignment operator - */ - ParetoDistribution &operator=(ParetoDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - k_ = other.k_; - xm_ = other.xm_; - logK_ = other.logK_; - kLogXm_ = other.kLogXm_; - kPlus1_ = other.kPlus1_; - kXmPowK_ = other.kXmPowK_; - negK_ = other.negK_; - logXm_ = other.logXm_; - } - return *this; - } - + ParetoDistribution(const ParetoDistribution &other) = default; + ParetoDistribution &operator=(const ParetoDistribution &other) = default; + ParetoDistribution(ParetoDistribution &&other) noexcept = default; + ParetoDistribution &operator=(ParetoDistribution &&other) noexcept = default; ~ParetoDistribution() override = default; /** diff --git a/include/libhmm/distributions/poisson_distribution.h b/include/libhmm/distributions/poisson_distribution.h index 01436aa..06dd030 100644 --- a/include/libhmm/distributions/poisson_distribution.h +++ b/include/libhmm/distributions/poisson_distribution.h @@ -91,49 +91,10 @@ class PoissonDistribution : public DistributionBase { updateCache(); } - /** - * Copy constructor - */ - PoissonDistribution(const PoissonDistribution &other) - : DistributionBase{other}, lambda_{other.lambda_}, logLambda_{other.logLambda_}, - expNegLambda_{other.expNegLambda_}, sqrtLambda_{other.sqrtLambda_}, - invSqrtLambda_{other.invSqrtLambda_}, sqrtTwoPiLambda_{other.sqrtTwoPiLambda_}, - smallFactorials_{other.smallFactorials_} {} - - PoissonDistribution &operator=(const PoissonDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - lambda_ = other.lambda_; - logLambda_ = other.logLambda_; - expNegLambda_ = other.expNegLambda_; - sqrtLambda_ = other.sqrtLambda_; - invSqrtLambda_ = other.invSqrtLambda_; - sqrtTwoPiLambda_ = other.sqrtTwoPiLambda_; - smallFactorials_ = other.smallFactorials_; - } - return *this; - } - - PoissonDistribution(PoissonDistribution &&other) noexcept - : DistributionBase{std::move(other)}, lambda_{other.lambda_}, logLambda_{other.logLambda_}, - expNegLambda_{other.expNegLambda_}, sqrtLambda_{other.sqrtLambda_}, - invSqrtLambda_{other.invSqrtLambda_}, sqrtTwoPiLambda_{other.sqrtTwoPiLambda_}, - smallFactorials_{std::move(other.smallFactorials_)} {} - - PoissonDistribution &operator=(PoissonDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - lambda_ = other.lambda_; - logLambda_ = other.logLambda_; - expNegLambda_ = other.expNegLambda_; - sqrtLambda_ = other.sqrtLambda_; - invSqrtLambda_ = other.invSqrtLambda_; - sqrtTwoPiLambda_ = other.sqrtTwoPiLambda_; - smallFactorials_ = std::move(other.smallFactorials_); - } - return *this; - } - + PoissonDistribution(const PoissonDistribution &other) = default; + PoissonDistribution &operator=(const PoissonDistribution &other) = default; + PoissonDistribution(PoissonDistribution &&other) noexcept = default; + PoissonDistribution &operator=(PoissonDistribution &&other) noexcept = default; ~PoissonDistribution() override = default; /** diff --git a/include/libhmm/distributions/rayleigh_distribution.h b/include/libhmm/distributions/rayleigh_distribution.h index 8d09308..e33bf49 100644 --- a/include/libhmm/distributions/rayleigh_distribution.h +++ b/include/libhmm/distributions/rayleigh_distribution.h @@ -115,60 +115,10 @@ class RayleighDistribution : public DistributionBase { updateCache(); } - /** - * Copy constructor - */ - RayleighDistribution(const RayleighDistribution &other) - : DistributionBase{other}, sigma_{other.sigma_}, logSigma_{other.logSigma_}, - invSigma_{other.invSigma_}, invSigmaSquared_{other.invSigmaSquared_}, - negHalfInvSigmaSquared_{other.negHalfInvSigmaSquared_}, - sigmaSquared_{other.sigmaSquared_}, mean_{other.mean_}, variance_{other.variance_} {} - - /** - * Copy assignment operator - */ - RayleighDistribution &operator=(const RayleighDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - sigma_ = other.sigma_; - logSigma_ = other.logSigma_; - invSigma_ = other.invSigma_; - invSigmaSquared_ = other.invSigmaSquared_; - negHalfInvSigmaSquared_ = other.negHalfInvSigmaSquared_; - sigmaSquared_ = other.sigmaSquared_; - mean_ = other.mean_; - variance_ = other.variance_; - } - return *this; - } - - /** - * Move constructor - */ - RayleighDistribution(RayleighDistribution &&other) noexcept - : DistributionBase{std::move(other)}, sigma_{other.sigma_}, logSigma_{other.logSigma_}, - invSigma_{other.invSigma_}, invSigmaSquared_{other.invSigmaSquared_}, - negHalfInvSigmaSquared_{other.negHalfInvSigmaSquared_}, - sigmaSquared_{other.sigmaSquared_}, mean_{other.mean_}, variance_{other.variance_} {} - - /** - * Move assignment operator - */ - RayleighDistribution &operator=(RayleighDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - sigma_ = other.sigma_; - logSigma_ = other.logSigma_; - invSigma_ = other.invSigma_; - invSigmaSquared_ = other.invSigmaSquared_; - negHalfInvSigmaSquared_ = other.negHalfInvSigmaSquared_; - sigmaSquared_ = other.sigmaSquared_; - mean_ = other.mean_; - variance_ = other.variance_; - } - return *this; - } - + RayleighDistribution(const RayleighDistribution &other) = default; + RayleighDistribution &operator=(const RayleighDistribution &other) = default; + RayleighDistribution(RayleighDistribution &&other) noexcept = default; + RayleighDistribution &operator=(RayleighDistribution &&other) noexcept = default; ~RayleighDistribution() override = default; [[nodiscard]] double getProbability(double value) const override; diff --git a/include/libhmm/distributions/student_t_distribution.h b/include/libhmm/distributions/student_t_distribution.h index 4fed5cf..0b0fc9d 100644 --- a/include/libhmm/distributions/student_t_distribution.h +++ b/include/libhmm/distributions/student_t_distribution.h @@ -110,10 +110,10 @@ class StudentTDistribution : public DistributionBase { */ StudentTDistribution(double degrees_of_freedom, double location, double scale); - StudentTDistribution(const StudentTDistribution &other); - StudentTDistribution &operator=(const StudentTDistribution &other); - StudentTDistribution(StudentTDistribution &&other) noexcept; - StudentTDistribution &operator=(StudentTDistribution &&other) noexcept; + StudentTDistribution(const StudentTDistribution &other) = default; + StudentTDistribution &operator=(const StudentTDistribution &other) = default; + StudentTDistribution(StudentTDistribution &&other) noexcept = default; + StudentTDistribution &operator=(StudentTDistribution &&other) noexcept = default; ~StudentTDistribution() override = default; /** @@ -217,14 +217,6 @@ class StudentTDistribution : public DistributionBase { */ std::string toString() const override; - /** - * @brief Create distribution from string representation - * @param str String representation - * @return StudentTDistribution object - * @throws std::invalid_argument if string format is invalid - */ - static StudentTDistribution fromString(const std::string &str); - /** * @brief Equality comparison operator * @param other Other distribution to compare with diff --git a/include/libhmm/distributions/uniform_distribution.h b/include/libhmm/distributions/uniform_distribution.h index 84f8bec..c7d066d 100644 --- a/include/libhmm/distributions/uniform_distribution.h +++ b/include/libhmm/distributions/uniform_distribution.h @@ -55,51 +55,10 @@ class UniformDistribution : public DistributionBase { */ UniformDistribution(double a, double b); - UniformDistribution(const UniformDistribution &other) - : DistributionBase{other}, a_{other.a_}, b_{other.b_}, cached_pdf_{other.cached_pdf_}, - cached_log_pdf_{other.cached_log_pdf_}, cached_range_{other.cached_range_}, - cached_inv_range_{other.cached_inv_range_}, cached_mean_{other.cached_mean_}, - cached_variance_{other.cached_variance_}, cached_std_dev_{other.cached_std_dev_} {} - - UniformDistribution &operator=(const UniformDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - a_ = other.a_; - b_ = other.b_; - cached_pdf_ = other.cached_pdf_; - cached_log_pdf_ = other.cached_log_pdf_; - cached_range_ = other.cached_range_; - cached_inv_range_ = other.cached_inv_range_; - cached_mean_ = other.cached_mean_; - cached_variance_ = other.cached_variance_; - cached_std_dev_ = other.cached_std_dev_; - } - return *this; - } - - UniformDistribution(UniformDistribution &&other) noexcept - : DistributionBase{std::move(other)}, a_{other.a_}, b_{other.b_}, - cached_pdf_{other.cached_pdf_}, cached_log_pdf_{other.cached_log_pdf_}, - cached_range_{other.cached_range_}, cached_inv_range_{other.cached_inv_range_}, - cached_mean_{other.cached_mean_}, cached_variance_{other.cached_variance_}, - cached_std_dev_{other.cached_std_dev_} {} - - UniformDistribution &operator=(UniformDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - a_ = other.a_; - b_ = other.b_; - cached_pdf_ = other.cached_pdf_; - cached_log_pdf_ = other.cached_log_pdf_; - cached_range_ = other.cached_range_; - cached_inv_range_ = other.cached_inv_range_; - cached_mean_ = other.cached_mean_; - cached_variance_ = other.cached_variance_; - cached_std_dev_ = other.cached_std_dev_; - } - return *this; - } - + UniformDistribution(const UniformDistribution &other) = default; + UniformDistribution &operator=(const UniformDistribution &other) = default; + UniformDistribution(UniformDistribution &&other) noexcept = default; + UniformDistribution &operator=(UniformDistribution &&other) noexcept = default; ~UniformDistribution() override = default; [[nodiscard]] double getProbability(double val) const override; diff --git a/include/libhmm/distributions/weibull_distribution.h b/include/libhmm/distributions/weibull_distribution.h index 3de493b..6c1a78a 100644 --- a/include/libhmm/distributions/weibull_distribution.h +++ b/include/libhmm/distributions/weibull_distribution.h @@ -110,53 +110,10 @@ class WeibullDistribution : public DistributionBase { /** * Copy constructor */ - WeibullDistribution(const WeibullDistribution &other) - : DistributionBase{other}, k_{other.k_}, lambda_{other.lambda_}, logK_{other.logK_}, - logLambda_{other.logLambda_}, kMinus1_{other.kMinus1_}, invLambda_{other.invLambda_}, - kOverLambda_{other.kOverLambda_} {} - - /** - * Copy assignment operator - */ - WeibullDistribution &operator=(const WeibullDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - k_ = other.k_; - lambda_ = other.lambda_; - logK_ = other.logK_; - logLambda_ = other.logLambda_; - kMinus1_ = other.kMinus1_; - invLambda_ = other.invLambda_; - kOverLambda_ = other.kOverLambda_; - } - return *this; - } - - /** - * Move constructor - */ - WeibullDistribution(WeibullDistribution &&other) noexcept - : DistributionBase{std::move(other)}, k_{other.k_}, lambda_{other.lambda_}, - logK_{other.logK_}, logLambda_{other.logLambda_}, kMinus1_{other.kMinus1_}, - invLambda_{other.invLambda_}, kOverLambda_{other.kOverLambda_} {} - - /** - * Move assignment operator - */ - WeibullDistribution &operator=(WeibullDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - k_ = other.k_; - lambda_ = other.lambda_; - logK_ = other.logK_; - logLambda_ = other.logLambda_; - kMinus1_ = other.kMinus1_; - invLambda_ = other.invLambda_; - kOverLambda_ = other.kOverLambda_; - } - return *this; - } - + WeibullDistribution(const WeibullDistribution &other) = default; + WeibullDistribution &operator=(const WeibullDistribution &other) = default; + WeibullDistribution(WeibullDistribution &&other) noexcept = default; + WeibullDistribution &operator=(WeibullDistribution &&other) noexcept = default; ~WeibullDistribution() override = default; /** @@ -178,6 +135,12 @@ class WeibullDistribution : public DistributionBase { /** Weighted MOM: same approach using weighted mean and variance. */ void fit(std::span data, std::span weights) override; +private: + /// Validates mean/var and applies MOM estimates via weibull_mom_fit. + /// Calls reset() on any failure; otherwise sets k_, lambda_ and invalidates cache. + void apply_fit_params(double mean, double var); + +public: /** Returns false — Weibull is a continuous distribution. */ [[nodiscard]] bool isDiscrete() const noexcept override { return false; } diff --git a/include/libhmm/math/weighted_stats.h b/include/libhmm/math/weighted_stats.h new file mode 100644 index 0000000..2fb8f8f --- /dev/null +++ b/include/libhmm/math/weighted_stats.h @@ -0,0 +1,26 @@ +#pragma once + +#include +#include +#include + +namespace libhmm::detail { + +/// Result of a two-pass weighted statistics computation. +struct WeightedStats { + double mean{0.0}; ///< Weighted mean. + double variance{0.0}; ///< Weighted variance (biased: sum(wi*(xi-mean)^2) / sumW). +}; + +/// Computes weighted mean and variance via a two-pass algorithm. +/// No data filter is applied — all elements contribute. +/// Returns nullopt if the total weight is non-positive or NaN. +[[nodiscard]] std::optional +compute_weighted_stats(std::span data, std::span weights) noexcept; + +/// Computes the weighted mean only (single pass: sumWX / sumW). +/// Returns nullopt if the total weight is non-positive or NaN. +[[nodiscard]] std::optional compute_weighted_mean(std::span data, + std::span weights) noexcept; + +} // namespace libhmm::detail diff --git a/include/libhmm/training/baum_welch_trainer.h b/include/libhmm/training/baum_welch_trainer.h index 72b1814..fa741a3 100755 --- a/include/libhmm/training/baum_welch_trainer.h +++ b/include/libhmm/training/baum_welch_trainer.h @@ -28,6 +28,28 @@ class BaumWelchTrainer : public Trainer { void train() override; private: + /// Accumulates gamma statistics (emission data/weights, pi numerator, trans + /// denominator) for one observation sequence. + void accumulate_gamma(const Matrix &logAlpha, const Matrix &logBeta, const ObservationSet &obs, + double logP, std::size_t N, std::vector> &emisData, + std::vector> &emisWts, std::vector &piNum, + std::vector &transDen) const noexcept; + + /// Accumulates xi statistics (expected transition counts) for one sequence. + /// Dispatches to sparse or dense path based on hasZeroTransitions. + void accumulate_xi(const double *logAlphaData, const double *logBetaData, + const std::vector &logEmitByTime, + const std::vector &logTransT, double logP, std::size_t T, + std::size_t N, bool hasZeroTransitions, + std::vector &transNumT) const noexcept; + + /// M-step: normalize piNum and update the HMM initial distribution. + static void m_step_pi(Hmm &hmm, std::size_t N, const std::vector &piNum); + + /// M-step: normalize transNumT/transDen and update the HMM transition matrix. + static void m_step_transitions(Hmm &hmm, std::size_t N, const std::vector &transNumT, + const std::vector &transDen); + /// Numerically stable log(exp(a) + exp(b)). static double logSumExp(double a, double b) noexcept; }; // class BaumWelchTrainer diff --git a/include/libhmm/training/trainer.h b/include/libhmm/training/trainer.h index f38812b..633b7cc 100644 --- a/include/libhmm/training/trainer.h +++ b/include/libhmm/training/trainer.h @@ -45,6 +45,31 @@ class Trainer { /// Execute one full training pass, updating the HMM in place. virtual void train() = 0; + +protected: + /// Refits emission distributions for all states after an E-step. + /// States with no observations are reset to defaults. + /// If weights is non-empty the weighted fit overload is used (Baum-Welch); + /// otherwise the unweighted overload is used (Viterbi hard assignment). + static void apply_emission_fits(Hmm &hmm, std::size_t numStates, + const std::vector> &data, + const std::vector> &weights = {}) { + for (std::size_t i = 0; i < numStates; ++i) { + const std::size_t M = data[i].size(); + if (M == 0) { + hmm.getDistribution(i).reset(); + continue; + } + if (weights.empty()) { + hmm.getDistribution(i).fit(std::span(data[i].data(), M)); + } else { + hmm.getDistribution(i).fit(std::span(data[i].data(), M), + std::span(weights[i].data(), M)); + } + } + } + +public: /// @return Reference to the HMM under training. [[nodiscard]] Hmm &getHmmRef() const noexcept { return hmm_ref_.get(); } diff --git a/include/libhmm/training/viterbi_trainer.h b/include/libhmm/training/viterbi_trainer.h index 457cdd5..622fcb8 100755 --- a/include/libhmm/training/viterbi_trainer.h +++ b/include/libhmm/training/viterbi_trainer.h @@ -67,6 +67,15 @@ class ViterbiTrainer : public Trainer { bool maxItersReached_{false}; double lastLogProb_{-std::numeric_limits::infinity()}; + /// Runs ViterbiCalculator on one sequence, accumulates pi/trans/emission stats. + /// Returns the log-probability contribution, or -∞ if the sequence is skipped + /// (non-finite log-prob or any exception). + static double accum_sequence(const Hmm &hmm, const ObservationSet &obs, Vector &pi, + Matrix &trans, + std::vector> &emisData) noexcept; + + /// Normalizes pi and every transition row, then commits both to the HMM. + static void normalize_and_commit(Hmm &hmm, std::size_t N, Vector &pi, Matrix &trans); /// Run one Viterbi iteration; returns total log-probability over all sequences. double runIteration(); }; // class ViterbiTrainer diff --git a/scripts/check-pragma-once.py b/scripts/check-pragma-once.py new file mode 100644 index 0000000..33bc42f --- /dev/null +++ b/scripts/check-pragma-once.py @@ -0,0 +1,29 @@ +#!/usr/bin/env python3 +"""Verify that all C/C++ header files use #pragma once. + +libhmm uses #pragma once as the header guard convention. This hook +prevents regressions when new headers are added. +""" + +import sys +from pathlib import Path + +exit_code = 0 + +for path_str in sys.argv[1:]: + path = Path(path_str) + if not path.exists(): + continue # deleted in this commit + + try: + text = path.read_text(encoding="utf-8", errors="replace") + except OSError: + continue + + if not any(line.startswith("#pragma once") for line in text.splitlines()): + print(f"ERROR: {path} is missing '#pragma once'") + print(" libhmm uses #pragma once as the header guard convention.") + print(" Add '#pragma once' as the first line of the file.") + exit_code = 1 + +sys.exit(exit_code) diff --git a/src/calculators/forward_backward_calculator.cpp b/src/calculators/forward_backward_calculator.cpp index b0ab429..6b5dcb0 100755 --- a/src/calculators/forward_backward_calculator.cpp +++ b/src/calculators/forward_backward_calculator.cpp @@ -10,6 +10,20 @@ namespace libhmm { namespace { constexpr double LOG_ZERO = -std::numeric_limits::infinity(); + +/// Initialises log-α at t=0 from log-π and the first row of log-emissions. +void init_log_forward(double *alphaData, const Vector &pi, const double *emitRow0, + std::size_t N) noexcept { + for (std::size_t i = 0; i < N; ++i) + alphaData[i] = (pi(i) > 0.0) ? std::log(pi(i)) + emitRow0[i] : LOG_ZERO; +} + +/// Sets log-β at t=T-1 to 0 (backward terminal condition: log(1) = 0). +void init_log_backward(double *betaData, std::size_t T, std::size_t N) noexcept { + double *finalRow = betaData + (T - 1) * N; + for (std::size_t i = 0; i < N; ++i) + finalRow[i] = 0.0; +} } // namespace FbRecurrenceMode @@ -108,18 +122,7 @@ void ForwardBackwardCalculator::compute() { // --------------------------------------------------------------------------- void ForwardBackwardCalculator::precomputeLogTransitions() { - const Hmm &hmm = getHmmRef(); - const Matrix &trans = hmm.getTrans(); - logTrans_.resize(numStates_, numStates_); - logTransT_.resize(numStates_, numStates_); - for (std::size_t i = 0; i < numStates_; ++i) { - for (std::size_t j = 0; j < numStates_; ++j) { - const double a = trans(i, j); - const double logA = (a > 0.0) ? std::log(a) : LOG_ZERO; - logTrans_(i, j) = logA; - logTransT_(j, i) = logA; - } - } + precompute_log_transitions(getHmmRef(), numStates_, logTrans_, logTransT_); } void ForwardBackwardCalculator::computeLogForward() { @@ -139,12 +142,7 @@ void ForwardBackwardCalculator::computeLogForwardPairwise() { const double *emitByTimeData = logEmitByTime_.data(); double *alphaData = logAlpha_.data(); - // t = 0. - const double *emitRow0 = emitByTimeData; - for (std::size_t i = 0; i < N; ++i) { - const double logPi = (pi(i) > 0.0) ? std::log(pi(i)) : LOG_ZERO; - alphaData[i] = logPi + emitRow0[i]; - } + init_log_forward(alphaData, pi, emitByTimeData, N); // t > 0. for (std::size_t t = 1; t < T; ++t) { @@ -171,12 +169,7 @@ void ForwardBackwardCalculator::computeLogForwardMaxReduce() { const double *emitByTimeData = logEmitByTime_.data(); double *alphaData = logAlpha_.data(); - // t = 0. - const double *emitRow0 = emitByTimeData; - for (std::size_t i = 0; i < N; ++i) { - const double logPi = (pi(i) > 0.0) ? std::log(pi(i)) : LOG_ZERO; - alphaData[i] = logPi + emitRow0[i]; - } + init_log_forward(alphaData, pi, emitByTimeData, N); // t > 0. for (std::size_t t = 1; t < T; ++t) { @@ -217,11 +210,7 @@ void ForwardBackwardCalculator::computeLogBackwardPairwise() { const double *emitByTimeData = logEmitByTime_.data(); double *betaData = logBeta_.data(); - // t = T - 1. - double *finalBetaRow = betaData + (T - 1) * N; - for (std::size_t i = 0; i < N; ++i) { - finalBetaRow[i] = 0.0; - } + init_log_backward(betaData, T, N); // t < T - 1. if (T > 1) { @@ -251,11 +240,7 @@ void ForwardBackwardCalculator::computeLogBackwardMaxReduce() { const double *emitByTimeData = logEmitByTime_.data(); double *betaData = logBeta_.data(); - // t = T - 1. - double *finalBetaRow = betaData + (T - 1) * N; - for (std::size_t i = 0; i < N; ++i) { - finalBetaRow[i] = 0.0; - } + init_log_backward(betaData, T, N); // t < T - 1. if (T > 1) { diff --git a/src/calculators/viterbi_calculator.cpp b/src/calculators/viterbi_calculator.cpp index ae5c18e..a789788 100755 --- a/src/calculators/viterbi_calculator.cpp +++ b/src/calculators/viterbi_calculator.cpp @@ -69,18 +69,7 @@ StateSequence ViterbiCalculator::decode() { // --------------------------------------------------------------------------- void ViterbiCalculator::precomputeLogTransitions() { - const Hmm &hmm = getHmmRef(); - const Matrix &trans = hmm.getTrans(); - logTrans_.resize(numStates_, numStates_); - logTransT_.resize(numStates_, numStates_); - for (std::size_t i = 0; i < numStates_; ++i) { - for (std::size_t j = 0; j < numStates_; ++j) { - const double a = trans(i, j); - const double logA = (a > 0.0) ? std::log(a) : LOG_ZERO; - logTrans_(i, j) = logA; - logTransT_(j, i) = logA; - } - } + precompute_log_transitions(getHmmRef(), numStates_, logTrans_, logTransT_); } void ViterbiCalculator::runViterbi() { diff --git a/src/common/numerical_stability.cpp b/src/common/numerical_stability.cpp index ede2662..021f9e4 100644 --- a/src/common/numerical_stability.cpp +++ b/src/common/numerical_stability.cpp @@ -477,120 +477,128 @@ ErrorRecovery::RecoveryStrategy ErrorRecovery::getRecommendedStrategy(std::size_ return RecoveryStrategy::ADAPTIVE; // Default to adaptive } +// --------------------------------------------------------------------------- +// Internal helpers +// --------------------------------------------------------------------------- + +namespace { + +/// Classifies a single scalar value into the diagnostic report, updating +/// hasNaN, hasInfinity, hasUnderflow, hasOverflow, minValue, and maxValue. +void classify_value(double val, NumericalDiagnostics::DiagnosticReport &report) noexcept { + if (std::isnan(val)) { + report.hasNaN = true; + } else if (std::isinf(val)) { + report.hasInfinity = true; + } else if (std::isfinite(val)) { + report.minValue = std::min(report.minValue, val); + report.maxValue = std::max(report.maxValue, val); + if (val < NumericalConstants::MIN_PROBABILITY && val > 0.0) + report.hasUnderflow = true; + if (val > 1e100) + report.hasOverflow = true; + } +} + +/// Appends the four standard numerical-issue recommendations (NaN, infinity, +/// underflow, overflow) to oss if the corresponding flags are set. +void append_base_recommendations(const NumericalDiagnostics::DiagnosticReport &report, + const std::string &name, std::ostringstream &oss) { + if (report.hasNaN) + oss << "NaN values detected in " << name << ". "; + if (report.hasInfinity) + oss << "Infinite values detected in " << name << ". "; + if (report.hasUnderflow) + oss << "Underflow detected in " << name << ". Consider using scaled arithmetic. "; + if (report.hasOverflow) + oss << "Overflow detected in " << name << ". Consider using log-space arithmetic. "; +} + +/// Returns true if a report contains any numerical issue worth surfacing. +bool report_has_issues(const NumericalDiagnostics::DiagnosticReport &r) noexcept { + return r.hasNaN || r.hasInfinity || r.hasUnderflow || r.hasOverflow || r.conditionNumber > 1e12; +} + +} // namespace + +// --------------------------------------------------------------------------- // NumericalDiagnostics implementation +// --------------------------------------------------------------------------- + +/** + * @brief Analyzes a vector for numerical health issues. + * + * Classifies each element and sets the corresponding flags in the returned + * report. Generates human-readable recommendations when issues are found. + * + * @param vector The vector to examine. + * @param name Human-readable label used in recommendation strings. + * @return DiagnosticReport with health flags and recommendations populated. + */ NumericalDiagnostics::DiagnosticReport NumericalDiagnostics::analyzeVector(const Vector &vector, const std::string &name) { DiagnosticReport report; report.problemSize = vector.size(); - for (std::size_t i = 0; i < vector.size(); ++i) { - double val = vector(i); - - if (std::isnan(val)) { - report.hasNaN = true; - } else if (std::isinf(val)) { - report.hasInfinity = true; - } else if (std::isfinite(val)) { - report.minValue = std::min(report.minValue, val); - report.maxValue = std::max(report.maxValue, val); - - if (val < NumericalConstants::MIN_PROBABILITY && val > 0.0) { - report.hasUnderflow = true; - } - if (val > 1e100) { - report.hasOverflow = true; - } - } - } + for (std::size_t i = 0; i < vector.size(); ++i) + classify_value(vector(i), report); - // Generate recommendations std::ostringstream oss; - if (report.hasNaN) { - oss << "NaN values detected in " << name << ". "; - } - if (report.hasInfinity) { - oss << "Infinite values detected in " << name << ". "; - } - if (report.hasUnderflow) { - oss << "Underflow detected in " << name << ". Consider using scaled arithmetic. "; - } - if (report.hasOverflow) { - oss << "Overflow detected in " << name << ". Consider using log-space arithmetic. "; - } - + append_base_recommendations(report, name, oss); report.recommendations = oss.str(); return report; } +/** + * @brief Analyzes a matrix for numerical health issues. + * + * Classifies each element and also estimates the condition number. Generates + * human-readable recommendations when issues are found. + * + * @param matrix The matrix to examine. + * @param name Human-readable label used in recommendation strings. + * @return DiagnosticReport with health flags, condition number, and + * recommendations populated. + */ NumericalDiagnostics::DiagnosticReport NumericalDiagnostics::analyzeMatrix(const Matrix &matrix, const std::string &name) { DiagnosticReport report; report.problemSize = matrix.size1() * matrix.size2(); - for (std::size_t i = 0; i < matrix.size1(); ++i) { - for (std::size_t j = 0; j < matrix.size2(); ++j) { - double val = matrix(i, j); - - if (std::isnan(val)) { - report.hasNaN = true; - } else if (std::isinf(val)) { - report.hasInfinity = true; - } else if (std::isfinite(val)) { - report.minValue = std::min(report.minValue, val); - report.maxValue = std::max(report.maxValue, val); - - if (val < NumericalConstants::MIN_PROBABILITY && val > 0.0) { - report.hasUnderflow = true; - } - if (val > 1e100) { - report.hasOverflow = true; - } - } - } - } + for (std::size_t i = 0; i < matrix.size1(); ++i) + for (std::size_t j = 0; j < matrix.size2(); ++j) + classify_value(matrix(i, j), report); - // Estimate condition number (simplified) + // Condition number estimate (ratio of max to min diagonal element). report.conditionNumber = estimateConditionNumber(matrix); - // Generate recommendations std::ostringstream oss; - if (report.hasNaN) { - oss << "NaN values detected in " << name << ". "; - } - if (report.hasInfinity) { - oss << "Infinite values detected in " << name << ". "; - } - if (report.hasUnderflow) { - oss << "Underflow detected in " << name << ". Consider using scaled arithmetic. "; - } - if (report.hasOverflow) { - oss << "Overflow detected in " << name << ". Consider using log-space arithmetic. "; - } - if (report.conditionNumber > 1e12) { + append_base_recommendations(report, name, oss); + if (report.conditionNumber > 1e12) oss << "High condition number (" << report.conditionNumber << ") in " << name << ". Matrix may be ill-conditioned. "; - } - report.recommendations = oss.str(); return report; } +/** + * @brief Generates a consolidated health report from a set of diagnostic reports. + * + * Summarises issue flags, recommendations, average condition number, and + * total problem size across all supplied reports. + * + * @param reports Collection of per-component DiagnosticReport values. + * @return Formatted multi-line health report string. + */ std::string NumericalDiagnostics::generateHealthReport(const std::vector &reports) { std::ostringstream oss; oss << std::fixed << std::setprecision(3); oss << "=== Numerical Health Report ===\n\n"; - bool hasIssues = false; - for (const auto &report : reports) { - if (report.hasNaN || report.hasInfinity || report.hasUnderflow || report.hasOverflow || - report.conditionNumber > 1e12) { - hasIssues = true; - break; - } - } + const bool hasIssues = std::any_of(reports.begin(), reports.end(), report_has_issues); if (!hasIssues) { oss << "✅ All systems healthy - no numerical issues detected.\n"; diff --git a/src/common/weighted_stats.cpp b/src/common/weighted_stats.cpp new file mode 100644 index 0000000..16f27ad --- /dev/null +++ b/src/common/weighted_stats.cpp @@ -0,0 +1,45 @@ +#include "libhmm/math/weighted_stats.h" +#include "libhmm/math/constants.h" + +using namespace libhmm::constants; + +namespace libhmm::detail { + +std::optional compute_weighted_stats(std::span data, + std::span weights) noexcept { + double sumW = 0.0; + for (const double w : weights) + sumW += w; + if (sumW < precision::ZERO || std::isnan(sumW)) + return std::nullopt; + + // Pass 1 — weighted mean + double mean = 0.0; + for (std::size_t i = 0; i < data.size(); ++i) + mean += weights[i] * data[i]; + mean /= sumW; + + // Pass 2 — weighted variance (biased) + double var = 0.0; + for (std::size_t i = 0; i < data.size(); ++i) { + const double d = data[i] - mean; + var += weights[i] * d * d; + } + var /= sumW; + + return WeightedStats{mean, var}; +} + +std::optional compute_weighted_mean(std::span data, + std::span weights) noexcept { + double sumW = 0.0, sumWX = 0.0; + for (std::size_t i = 0; i < data.size(); ++i) { + sumW += weights[i]; + sumWX += weights[i] * data[i]; + } + if (sumW < precision::ZERO || std::isnan(sumW)) + return std::nullopt; + return sumWX / sumW; +} + +} // namespace libhmm::detail diff --git a/src/distributions/beta_distribution.cpp b/src/distributions/beta_distribution.cpp index 5a5740d..1f14e0a 100644 --- a/src/distributions/beta_distribution.cpp +++ b/src/distributions/beta_distribution.cpp @@ -1,4 +1,5 @@ #include "libhmm/distributions/beta_distribution.h" +#include "libhmm/math/weighted_stats.h" #include using namespace libhmm::constants; @@ -210,21 +211,13 @@ void BetaDistribution::fit(std::span data) { } void BetaDistribution::fit(std::span data, std::span weights) { - double sumW = 0.0; - for (const double w : weights) - sumW += w; - if (sumW < precision::ZERO || std::isnan(sumW)) { + const auto stats = detail::compute_weighted_stats(data, weights); + if (!stats) { reset(); return; } - double mean = 0.0; - for (std::size_t i = 0; i < data.size(); ++i) - mean += weights[i] * data[i]; - mean /= sumW; - double var = 0.0; - for (std::size_t i = 0; i < data.size(); ++i) - var += weights[i] * (data[i] - mean) * (data[i] - mean); - var /= sumW; + const double mean = stats->mean; + const double var = stats->variance; if (var <= precision::ZERO || mean <= precision::ZERO || mean >= math::ONE) { reset(); return; diff --git a/src/distributions/chi_squared_distribution.cpp b/src/distributions/chi_squared_distribution.cpp index a3baf09..13e0476 100644 --- a/src/distributions/chi_squared_distribution.cpp +++ b/src/distributions/chi_squared_distribution.cpp @@ -1,4 +1,5 @@ #include "libhmm/distributions/chi_squared_distribution.h" +#include "libhmm/math/weighted_stats.h" #include #include @@ -103,16 +104,12 @@ void ChiSquaredDistribution::fit(std::span data) { } void ChiSquaredDistribution::fit(std::span data, std::span weights) { - double sumW = 0.0, sumWX = 0.0; - for (std::size_t i = 0; i < data.size(); ++i) { - sumW += weights[i]; - sumWX += weights[i] * data[i]; - } - if (sumW < precision::ZERO || std::isnan(sumW)) { + const auto mean = detail::compute_weighted_mean(data, weights); + if (!mean) { reset(); return; } - double est = std::max(MIN_DEGREES_OF_FREEDOM, std::min(MAX_DEGREES_OF_FREEDOM, sumWX / sumW)); + double est = std::max(MIN_DEGREES_OF_FREEDOM, std::min(MAX_DEGREES_OF_FREEDOM, *mean)); setDegreesOfFreedom(est); } @@ -129,40 +126,6 @@ std::string ChiSquaredDistribution::toString() const { return oss.str(); } -ChiSquaredDistribution ChiSquaredDistribution::fromString(const std::string &str) { - // Expected format: "ChiSquared(k=value)" or "ChiSquared(df=value)" - std::string::size_type start = str.find('('); - std::string::size_type end = str.find(')', start); - - if (start == std::string::npos || end == std::string::npos) { - throw std::invalid_argument("Invalid ChiSquared distribution string format"); - } - - std::string params = str.substr(start + 1, end - start - 1); - - // Look for parameter patterns - std::string::size_type eq_pos = params.find('='); - if (eq_pos == std::string::npos) { - throw std::invalid_argument("Invalid ChiSquared parameter format"); - } - - std::string param_name = params.substr(0, eq_pos); - std::string param_value = params.substr(eq_pos + 1); - - // Remove whitespace - param_name.erase(std::remove_if(param_name.begin(), param_name.end(), ::isspace), - param_name.end()); - param_value.erase(std::remove_if(param_value.begin(), param_value.end(), ::isspace), - param_value.end()); - - if (param_name == "k" || param_name == "df") { - double df = std::stod(param_value); - return ChiSquaredDistribution(df); - } else { - throw std::invalid_argument("Unknown ChiSquared parameter: " + param_name); - } -} - bool ChiSquaredDistribution::operator==(const ChiSquaredDistribution &other) const { return std::abs(degrees_of_freedom_ - other.degrees_of_freedom_) < PARAMETER_TOLERANCE; } diff --git a/src/distributions/exponential_distribution.cpp b/src/distributions/exponential_distribution.cpp index 4a5b052..ac9d23e 100755 --- a/src/distributions/exponential_distribution.cpp +++ b/src/distributions/exponential_distribution.cpp @@ -1,5 +1,6 @@ #include "libhmm/distributions/exponential_distribution.h" -#include "libhmm/platform/simd_platform.h" // compile-time SIMD macros + intrinsics +#include "libhmm/math/weighted_stats.h" +#include "libhmm/platform/simd_platform.h" #include #include #include @@ -114,21 +115,12 @@ void ExponentialDistribution::fit(std::span data) { void ExponentialDistribution::fit(std::span data, std::span weights) { // Weighted MLE: λ = 1 / weighted_mean - double sumW = 0.0, sumWX = 0.0; - for (std::size_t i = 0; i < data.size(); ++i) { - sumW += weights[i]; - sumWX += weights[i] * data[i]; - } - if (sumW < precision::ZERO || std::isnan(sumW) || sumWX <= 0.0) { - reset(); - return; - } - const double weightedMean = sumWX / sumW; - if (weightedMean <= 0.0 || !std::isfinite(weightedMean)) { + const auto mean = detail::compute_weighted_mean(data, weights); + if (!mean || *mean <= 0.0 || !std::isfinite(*mean)) { reset(); return; } - const double lam = 1.0 / weightedMean; + const double lam = 1.0 / *mean; if (!std::isfinite(lam) || lam <= 0.0) { reset(); return; diff --git a/src/distributions/gamma_distribution.cpp b/src/distributions/gamma_distribution.cpp index 76cc848..3fbd6bd 100755 --- a/src/distributions/gamma_distribution.cpp +++ b/src/distributions/gamma_distribution.cpp @@ -99,6 +99,21 @@ double GammaDistribution::ligamma(double a, double x) noexcept { * * @param values Vector of observed data points */ +void GammaDistribution::apply_fit_params(double mean, double var) { + if (mean <= precision::ZERO || var <= precision::ZERO) { + reset(); + return; + } + const double newTheta = var / mean, newK = (mean * mean) / var; + if (!std::isfinite(newK) || !std::isfinite(newTheta) || newK <= 0.0 || newTheta <= 0.0) { + reset(); + return; + } + theta_ = newTheta; + k_ = newK; + invalidateCache(); +} + void GammaDistribution::fit(std::span data) { if (data.size() < 2) { reset(); @@ -118,19 +133,7 @@ void GammaDistribution::fit(std::span data) { reset(); return; } - const double var = m2 / (static_cast(count) - 1.0); - if (mean <= precision::ZERO || var <= precision::ZERO) { - reset(); - return; - } - const double newTheta = var / mean, newK = (mean * mean) / var; - if (!std::isfinite(newK) || !std::isfinite(newTheta) || newK <= 0.0 || newTheta <= 0.0) { - reset(); - return; - } - theta_ = newTheta; - k_ = newK; - invalidateCache(); + apply_fit_params(mean, m2 / (static_cast(count) - 1.0)); } void GammaDistribution::fit(std::span data, std::span weights) { @@ -154,19 +157,7 @@ void GammaDistribution::fit(std::span data, std::span #include #include @@ -81,16 +82,12 @@ void PoissonDistribution::fit(std::span data) { void PoissonDistribution::fit(std::span data, std::span weights) { // Weighted MLE: λ = weighted mean - double sumW = 0.0, sumWX = 0.0; - for (std::size_t i = 0; i < data.size(); ++i) { - sumW += weights[i]; - sumWX += weights[i] * data[i]; - } - if (sumW < precision::ZERO || std::isnan(sumW)) { + const auto mean = detail::compute_weighted_mean(data, weights); + if (!mean) { reset(); return; } - lambda_ = std::max(sumWX / sumW, precision::ZERO); + lambda_ = std::max(*mean, precision::ZERO); invalidateCache(); } diff --git a/src/distributions/student_t_distribution.cpp b/src/distributions/student_t_distribution.cpp index f7944d2..42afc73 100644 --- a/src/distributions/student_t_distribution.cpp +++ b/src/distributions/student_t_distribution.cpp @@ -24,64 +24,6 @@ StudentTDistribution::StudentTDistribution(double degrees_of_freedom, double loc updateCache(); } -StudentTDistribution::StudentTDistribution(const StudentTDistribution &other) - : DistributionBase{other}, degrees_of_freedom_(other.degrees_of_freedom_), - location_(other.location_), scale_(other.scale_), - cached_log_gamma_half_nu_plus_one_(other.cached_log_gamma_half_nu_plus_one_), - cached_log_gamma_half_nu_(other.cached_log_gamma_half_nu_), - cached_log_normalization_(other.cached_log_normalization_), - cached_normalization_factor_(other.cached_normalization_factor_), - cached_half_nu_plus_one_(other.cached_half_nu_plus_one_), - cached_half_nu_(other.cached_half_nu_), cached_inv_scale_(other.cached_inv_scale_), - cached_log_scale_(other.cached_log_scale_) {} - -StudentTDistribution &StudentTDistribution::operator=(const StudentTDistribution &other) { - if (this != &other) { - DistributionBase::operator=(other); - degrees_of_freedom_ = other.degrees_of_freedom_; - location_ = other.location_; - scale_ = other.scale_; - cached_log_gamma_half_nu_plus_one_ = other.cached_log_gamma_half_nu_plus_one_; - cached_log_gamma_half_nu_ = other.cached_log_gamma_half_nu_; - cached_log_normalization_ = other.cached_log_normalization_; - cached_normalization_factor_ = other.cached_normalization_factor_; - cached_half_nu_plus_one_ = other.cached_half_nu_plus_one_; - cached_half_nu_ = other.cached_half_nu_; - cached_inv_scale_ = other.cached_inv_scale_; - cached_log_scale_ = other.cached_log_scale_; - } - return *this; -} - -StudentTDistribution::StudentTDistribution(StudentTDistribution &&other) noexcept - : DistributionBase{std::move(other)}, degrees_of_freedom_(other.degrees_of_freedom_), - location_(other.location_), scale_(other.scale_), - cached_log_gamma_half_nu_plus_one_(other.cached_log_gamma_half_nu_plus_one_), - cached_log_gamma_half_nu_(other.cached_log_gamma_half_nu_), - cached_log_normalization_(other.cached_log_normalization_), - cached_normalization_factor_(other.cached_normalization_factor_), - cached_half_nu_plus_one_(other.cached_half_nu_plus_one_), - cached_half_nu_(other.cached_half_nu_), cached_inv_scale_(other.cached_inv_scale_), - cached_log_scale_(other.cached_log_scale_) {} - -StudentTDistribution &StudentTDistribution::operator=(StudentTDistribution &&other) noexcept { - if (this != &other) { - DistributionBase::operator=(std::move(other)); - degrees_of_freedom_ = other.degrees_of_freedom_; - location_ = other.location_; - scale_ = other.scale_; - cached_log_gamma_half_nu_plus_one_ = other.cached_log_gamma_half_nu_plus_one_; - cached_log_gamma_half_nu_ = other.cached_log_gamma_half_nu_; - cached_log_normalization_ = other.cached_log_normalization_; - cached_normalization_factor_ = other.cached_normalization_factor_; - cached_half_nu_plus_one_ = other.cached_half_nu_plus_one_; - cached_half_nu_ = other.cached_half_nu_; - cached_inv_scale_ = other.cached_inv_scale_; - cached_log_scale_ = other.cached_log_scale_; - } - return *this; -} - double StudentTDistribution::getProbability(double value) const { if (!isCacheValid()) updateCache(); @@ -295,40 +237,6 @@ std::string StudentTDistribution::toString() const { return oss.str(); } -StudentTDistribution StudentTDistribution::fromString(const std::string &str) { - // Expected format: "StudentT(ν=value)" or "StudentT(nu=value)" or "StudentT(df=value)" - std::string::size_type start = str.find('('); - std::string::size_type end = str.find(')', start); - - if (start == std::string::npos || end == std::string::npos) { - throw std::invalid_argument("Invalid StudentT distribution string format"); - } - - std::string params = str.substr(start + 1, end - start - 1); - - // Look for parameter patterns - std::string::size_type eq_pos = params.find('='); - if (eq_pos == std::string::npos) { - throw std::invalid_argument("Invalid StudentT parameter format"); - } - - std::string param_name = params.substr(0, eq_pos); - std::string param_value = params.substr(eq_pos + 1); - - // Remove whitespace - param_name.erase(std::remove_if(param_name.begin(), param_name.end(), ::isspace), - param_name.end()); - param_value.erase(std::remove_if(param_value.begin(), param_value.end(), ::isspace), - param_value.end()); - - if (param_name == "ν" || param_name == "nu" || param_name == "df") { - double df = std::stod(param_value); - return StudentTDistribution(df); - } else { - throw std::invalid_argument("Unknown StudentT parameter: " + param_name); - } -} - bool StudentTDistribution::operator==(const StudentTDistribution &other) const { return std::abs(degrees_of_freedom_ - other.degrees_of_freedom_) < PARAMETER_TOLERANCE && std::abs(location_ - other.location_) < PARAMETER_TOLERANCE && diff --git a/src/distributions/uniform_distribution.cpp b/src/distributions/uniform_distribution.cpp index e4f034b..52d4b8f 100644 --- a/src/distributions/uniform_distribution.cpp +++ b/src/distributions/uniform_distribution.cpp @@ -1,4 +1,5 @@ #include "libhmm/distributions/uniform_distribution.h" +#include "libhmm/math/weighted_stats.h" #include #include #include @@ -114,33 +115,18 @@ void UniformDistribution::fit(std::span data, std::spanvariance); if (halfRange < thresholds::MIN_DISTRIBUTION_PARAMETER || !std::isfinite(halfRange)) { reset(); return; } - a_ = mean - halfRange; - b_ = mean + halfRange; + a_ = stats->mean - halfRange; + b_ = stats->mean + halfRange; invalidateCache(); } diff --git a/src/distributions/weibull_distribution.cpp b/src/distributions/weibull_distribution.cpp index 3bfd766..3287d1c 100644 --- a/src/distributions/weibull_distribution.cpp +++ b/src/distributions/weibull_distribution.cpp @@ -79,6 +79,22 @@ static void weibull_mom_fit(double mean, double var, double &k_out, double &lamb lambda_out = mean / gamma_term; } +void WeibullDistribution::apply_fit_params(double mean, double var) { + if (var <= precision::ZERO || mean <= precision::ZERO) { + reset(); + return; + } + double k_est, lambda_est; + weibull_mom_fit(mean, var, k_est, lambda_est); + if (lambda_est > precision::ZERO && lambda_est < thresholds::MAX_DISTRIBUTION_PARAMETER) { + k_ = k_est; + lambda_ = lambda_est; + invalidateCache(); + } else { + reset(); + } +} + void WeibullDistribution::fit(std::span data) { if (data.size() < 2) { reset(); @@ -96,20 +112,7 @@ void WeibullDistribution::fit(std::span data) { mean += delta / static_cast(count); m2 += delta * (val - mean); } - const double var = m2 / (n - math::ONE); - if (var <= precision::ZERO || mean <= precision::ZERO) { - reset(); - return; - } - double k_est, lambda_est; - weibull_mom_fit(mean, var, k_est, lambda_est); - if (lambda_est > precision::ZERO && lambda_est < thresholds::MAX_DISTRIBUTION_PARAMETER) { - k_ = k_est; - lambda_ = lambda_est; - invalidateCache(); - } else { - reset(); - } + apply_fit_params(mean, m2 / (n - math::ONE)); } void WeibullDistribution::fit(std::span data, std::span weights) { @@ -128,20 +131,7 @@ void WeibullDistribution::fit(std::span data, std::span= 0.0 && std::isfinite(data[i]) && weights[i] > 0.0) var += weights[i] * (data[i] - mean) * (data[i] - mean); - var /= sumW; - if (var <= precision::ZERO || mean <= precision::ZERO) { - reset(); - return; - } - double k_est, lambda_est; - weibull_mom_fit(mean, var, k_est, lambda_est); - if (lambda_est > precision::ZERO && lambda_est < thresholds::MAX_DISTRIBUTION_PARAMETER) { - k_ = k_est; - lambda_ = lambda_est; - invalidateCache(); - } else { - reset(); - } + apply_fit_params(mean, var / sumW); } void WeibullDistribution::reset() noexcept { diff --git a/src/io/file_io_manager.cpp b/src/io/file_io_manager.cpp index 48ae5b1..dacfd8f 100644 --- a/src/io/file_io_manager.cpp +++ b/src/io/file_io_manager.cpp @@ -2,9 +2,37 @@ #include #include #include +#include #include #include +namespace { + +/// Opens filepath for writing with the given mode, calls writeOp(file), then +/// validates and closes the stream. Throws std::runtime_error on any failure. +void do_write(const std::filesystem::path &filepath, std::ios::openmode mode, + const std::function &writeOp) { + try { + std::ofstream file(filepath, mode); + if (!file.is_open()) { + throw std::runtime_error("Failed to open file for writing: " + filepath.string()); + } + writeOp(file); + if (file.bad()) { + throw std::runtime_error("Error occurred while writing to file: " + filepath.string()); + } + file.close(); + if (file.fail()) { + throw std::runtime_error("Failed to properly close file: " + filepath.string()); + } + } catch (const std::ios_base::failure &e) { + throw std::runtime_error("I/O error while writing file " + filepath.string() + ": " + + e.what()); + } +} + +} // namespace + namespace libhmm { std::string FileIOManager::readTextFile(const std::filesystem::path &filepath) { @@ -34,34 +62,10 @@ std::string FileIOManager::readTextFile(const std::filesystem::path &filepath) { void FileIOManager::writeTextFile(const std::filesystem::path &filepath, const std::string &content, bool append) { - try { - // Ensure directory exists - if (filepath.has_parent_path()) { - ensureDirectoryExists(filepath.parent_path()); - } - - auto mode = append ? (std::ios::out | std::ios::app) : (std::ios::out | std::ios::trunc); - std::ofstream file(filepath, mode); - - if (!file.is_open()) { - throw std::runtime_error("Failed to open file for writing: " + filepath.string()); - } - - file << content; - - if (file.bad()) { - throw std::runtime_error("Error occurred while writing to file: " + filepath.string()); - } - - file.close(); - if (file.fail()) { - throw std::runtime_error("Failed to properly close file: " + filepath.string()); - } - - } catch (const std::ios_base::failure &e) { - throw std::runtime_error("I/O error while writing file " + filepath.string() + ": " + - e.what()); - } + if (filepath.has_parent_path()) + ensureDirectoryExists(filepath.parent_path()); + const auto mode = append ? (std::ios::out | std::ios::app) : (std::ios::out | std::ios::trunc); + do_write(filepath, mode, [&content](std::ofstream &file) { file << content; }); } std::vector FileIOManager::readLines(const std::filesystem::path &filepath) { @@ -94,36 +98,13 @@ std::vector FileIOManager::readLines(const std::filesystem::path &f void FileIOManager::writeLines(const std::filesystem::path &filepath, const std::vector &lines, bool append) { - try { - // Ensure directory exists - if (filepath.has_parent_path()) { - ensureDirectoryExists(filepath.parent_path()); - } - - auto mode = append ? (std::ios::out | std::ios::app) : (std::ios::out | std::ios::trunc); - std::ofstream file(filepath, mode); - - if (!file.is_open()) { - throw std::runtime_error("Failed to open file for writing: " + filepath.string()); - } - - for (const auto &line : lines) { + if (filepath.has_parent_path()) + ensureDirectoryExists(filepath.parent_path()); + const auto mode = append ? (std::ios::out | std::ios::app) : (std::ios::out | std::ios::trunc); + do_write(filepath, mode, [&lines](std::ofstream &file) { + for (const auto &line : lines) file << line << '\n'; - } - - if (file.bad()) { - throw std::runtime_error("Error occurred while writing to file: " + filepath.string()); - } - - file.close(); - if (file.fail()) { - throw std::runtime_error("Failed to properly close file: " + filepath.string()); - } - - } catch (const std::ios_base::failure &e) { - throw std::runtime_error("I/O error while writing file " + filepath.string() + ": " + - e.what()); - } + }); } void FileIOManager::copyFile(const std::filesystem::path &source, diff --git a/src/training/baum_welch_trainer.cpp b/src/training/baum_welch_trainer.cpp index 96b251a..d98b8f1 100755 --- a/src/training/baum_welch_trainer.cpp +++ b/src/training/baum_welch_trainer.cpp @@ -20,6 +20,93 @@ BaumWelchTrainer::BaumWelchTrainer(Hmm &hmm, const ObservationLists &obsLists) BaumWelchTrainer::BaumWelchTrainer(Hmm *hmm, const ObservationLists &obsLists) : Trainer(hmm, obsLists) {} +// --------------------------------------------------------------------------- +// Private helpers +// --------------------------------------------------------------------------- + +void BaumWelchTrainer::accumulate_gamma(const Matrix &logAlpha, const Matrix &logBeta, + const ObservationSet &obs, double logP, std::size_t N, + std::vector> &emisData, + std::vector> &emisWts, + std::vector &piNum, + std::vector &transDen) const noexcept { + const std::size_t T = obs.size(); + const double *logAlphaData = logAlpha.data(); + const double *logBetaData = logBeta.data(); + for (std::size_t t = 0; t < T; ++t) { + const double *alphaRow = logAlphaData + t * N; + const double *betaRow = logBetaData + t * N; + const double obsVal = obs(t); + for (std::size_t i = 0; i < N; ++i) { + const double g = std::exp(alphaRow[i] + betaRow[i] - logP); + emisData[i].push_back(obsVal); + emisWts[i].push_back(g); + if (t == 0) + piNum[i] += g; + if (t < T - 1) + transDen[i] += g; + } + } +} + +void BaumWelchTrainer::accumulate_xi(const double *logAlphaData, const double *logBetaData, + const std::vector &logEmitByTime, + const std::vector &logTransT, double logP, + std::size_t T, std::size_t N, bool hasZeroTransitions, + std::vector &transNumT) const noexcept { + if (hasZeroTransitions) { + for (std::size_t t = 0; t + 1 < T; ++t) { + const double *alphaRow = logAlphaData + t * N; + const double *betaNextRow = logBetaData + (t + 1) * N; + const double *emitNextRow = logEmitByTime.data() + (t + 1) * N; + for (std::size_t j = 0; j < N; ++j) { + const double emitBetaNext = emitNextRow[j] + betaNextRow[j] - logP; + const double *transCol = logTransT.data() + j * N; + double *transNumCol = transNumT.data() + j * N; + for (std::size_t i = 0; i < N; ++i) { + if (transCol[i] == LOG_ZERO) + continue; + transNumCol[i] += std::exp(alphaRow[i] + transCol[i] + emitBetaNext); + } + } + } + } else { + for (std::size_t t = 0; t + 1 < T; ++t) { + const double *alphaRow = logAlphaData + t * N; + const double *betaNextRow = logBetaData + (t + 1) * N; + const double *emitNextRow = logEmitByTime.data() + (t + 1) * N; + for (std::size_t j = 0; j < N; ++j) { + const double emitBetaNext = emitNextRow[j] + betaNextRow[j] - logP; + const double *transCol = logTransT.data() + j * N; + double *transNumCol = transNumT.data() + j * N; + performance::detail::TranscendentalKernels::accumulate_exp_sum2_bias( + transNumCol, alphaRow, transCol, N, emitBetaNext); + } + } + } +} + +void BaumWelchTrainer::m_step_pi(Hmm &hmm, std::size_t N, const std::vector &piNum) { + double piSum = 0.0; + for (std::size_t i = 0; i < N; ++i) + piSum += piNum[i]; + Vector pi(N); + for (std::size_t i = 0; i < N; ++i) + pi(i) = (piSum > 0.0) ? piNum[i] / piSum : 1.0 / static_cast(N); + hmm.setPi(pi); +} + +void BaumWelchTrainer::m_step_transitions(Hmm &hmm, std::size_t N, + const std::vector &transNumT, + const std::vector &transDen) { + Matrix newTrans(N, N); + for (std::size_t i = 0; i < N; ++i) + for (std::size_t j = 0; j < N; ++j) + newTrans(i, j) = (transDen[i] > 0.0) ? transNumT[j * N + i] / transDen[i] + : 1.0 / static_cast(N); + hmm.setTrans(newTrans); +} + // --------------------------------------------------------------------------- // train() — one full EM pass over all sequences // --------------------------------------------------------------------------- @@ -98,59 +185,12 @@ void BaumWelchTrainer::train() { } } - // Accumulate gamma (per timestep per state) and pi/trans denominators - for (std::size_t t = 0; t < T; ++t) { - const double *alphaRow = logAlphaData + t * N; - const double *betaRow = logBetaData + t * N; - const double obsVal = obs(t); - for (std::size_t i = 0; i < N; ++i) { - const double g = std::exp(alphaRow[i] + betaRow[i] - logP); - emisData[i].push_back(obsVal); - emisWts[i].push_back(g); - if (t == 0) - piNum[i] += g; - if (t < T - 1) - transDen[i] += g; - } - } + accumulate_gamma(logAlpha, logBeta, obs, logP, N, emisData, emisWts, piNum, transDen); // Accumulate xi (transition counts). Dense models take a branch-free - // path; sparse models keep the zero-transition skip. - // Sparse path is intentionally scalar: masking non-zero transitions in - // a SIMD loop costs more than it saves for the typically small fraction - // of non-zero entries in a sparse model. - if (hasZeroTransitions) { - for (std::size_t t = 0; t + 1 < T; ++t) { - const double *alphaRow = logAlphaData + t * N; - const double *betaNextRow = logBetaData + (t + 1) * N; - const double *emitNextRow = logEmitByTime.data() + (t + 1) * N; - for (std::size_t j = 0; j < N; ++j) { - const double emitBetaNext = emitNextRow[j] + betaNextRow[j] - logP; - const double *transCol = logTransT.data() + j * N; - double *transNumCol = transNumT.data() + j * N; - for (std::size_t i = 0; i < N; ++i) { - if (transCol[i] == LOG_ZERO) { - continue; - } - const double logXi = alphaRow[i] + transCol[i] + emitBetaNext; - transNumCol[i] += std::exp(logXi); - } - } - } - } else { - for (std::size_t t = 0; t + 1 < T; ++t) { - const double *alphaRow = logAlphaData + t * N; - const double *betaNextRow = logBetaData + (t + 1) * N; - const double *emitNextRow = logEmitByTime.data() + (t + 1) * N; - for (std::size_t j = 0; j < N; ++j) { - const double emitBetaNext = emitNextRow[j] + betaNextRow[j] - logP; - const double *transCol = logTransT.data() + j * N; - double *transNumCol = transNumT.data() + j * N; - performance::detail::TranscendentalKernels::accumulate_exp_sum2_bias( - transNumCol, alphaRow, transCol, N, emitBetaNext); - } - } - } + // SIMD path; sparse models keep the zero-transition skip. + accumulate_xi(logAlphaData, logBetaData, logEmitByTime, logTransT, logP, T, N, + hasZeroTransitions, transNumT); ++validSeqs; } @@ -160,40 +200,9 @@ void BaumWelchTrainer::train() { "(all had zero probability under the current model)"); } - // ---- M-step: pi ---- - { - double piSum = 0.0; - for (std::size_t i = 0; i < N; ++i) - piSum += piNum[i]; - Vector pi(N); - for (std::size_t i = 0; i < N; ++i) { - pi(i) = (piSum > 0.0) ? piNum[i] / piSum : 1.0 / static_cast(N); - } - hmm.setPi(pi); - } - - // ---- M-step: transition matrix ---- - { - Matrix newTrans(N, N); - for (std::size_t i = 0; i < N; ++i) { - for (std::size_t j = 0; j < N; ++j) { - newTrans(i, j) = (transDen[i] > 0.0) ? transNumT[j * N + i] / transDen[i] - : 1.0 / static_cast(N); - } - } - hmm.setTrans(newTrans); - } - - // ---- M-step: emission distributions ---- - for (std::size_t i = 0; i < N; ++i) { - const std::size_t M = emisData[i].size(); - if (M == 0) { - hmm.getDistribution(i).reset(); - continue; - } - hmm.getDistribution(i).fit(std::span(emisData[i].data(), M), - std::span(emisWts[i].data(), M)); - } + m_step_pi(hmm, N, piNum); + m_step_transitions(hmm, N, transNumT, transDen); + apply_emission_fits(hmm, N, emisData, emisWts); } // --------------------------------------------------------------------------- diff --git a/src/training/viterbi_trainer.cpp b/src/training/viterbi_trainer.cpp index d159bb0..678056f 100755 --- a/src/training/viterbi_trainer.cpp +++ b/src/training/viterbi_trainer.cpp @@ -78,6 +78,62 @@ void ViterbiTrainer::train() { maxItersReached_ = true; } +// --------------------------------------------------------------------------- +// Private helpers +// --------------------------------------------------------------------------- + +double ViterbiTrainer::accum_sequence(const Hmm &hmm, const ObservationSet &obs, Vector &pi, + Matrix &trans, + std::vector> &emisData) noexcept { + try { + ViterbiCalculator vc(hmm, obs); + const double lp = vc.getLogProbability(); + if (!std::isfinite(lp)) + return lp; + const StateSequence &seq = vc.getStateSequence(); + const std::size_t T = obs.size(); + pi(static_cast(seq(0))) += 1.0; + for (std::size_t t = 0; t < T; ++t) { + const std::size_t s = static_cast(seq(t)); + emisData[s].push_back(obs(t)); + if (t + 1 < T) + trans(s, static_cast(seq(t + 1))) += 1.0; + } + return lp; + } catch (...) { + return -std::numeric_limits::infinity(); + } +} + +void ViterbiTrainer::normalize_and_commit(Hmm &hmm, std::size_t N, Vector &pi, Matrix &trans) { + // Normalize pi. + double piSum = 0.0; + for (std::size_t i = 0; i < N; ++i) + piSum += pi(i); + if (piSum > 0.0) { + for (std::size_t i = 0; i < N; ++i) + pi(i) /= piSum; + } else { + for (std::size_t i = 0; i < N; ++i) + pi(i) = 1.0 / static_cast(N); + } + hmm.setPi(pi); + // Normalize each transition row. + for (std::size_t i = 0; i < N; ++i) { + double rowSum = 0.0; + for (std::size_t j = 0; j < N; ++j) + rowSum += trans(i, j); + if (rowSum > 0.0) { + for (std::size_t j = 0; j < N; ++j) + trans(i, j) /= rowSum; + } else { + for (std::size_t j = 0; j < N; ++j) + trans(i, j) = 1.0 / static_cast(N); + } + } + hmm.setTrans(trans); +} + // --------------------------------------------------------------------------- // runIteration() — one Viterbi pass // --------------------------------------------------------------------------- @@ -98,75 +154,18 @@ double ViterbiTrainer::runIteration() { for (const auto &obs : obsLists_) { if (obs.size() == 0) continue; - try { - ViterbiCalculator vc(hmm, obs); - const double lp = vc.getLogProbability(); - if (!std::isfinite(lp)) - continue; - - totalLogProb += lp; - const StateSequence &seq = vc.getStateSequence(); - const std::size_t T = obs.size(); - - pi(static_cast(seq(0))) += 1.0; - - for (std::size_t t = 0; t < T; ++t) { - const std::size_t s = static_cast(seq(t)); - emisData[s].push_back(obs(t)); - if (t + 1 < T) { - const std::size_t sNext = static_cast(seq(t + 1)); - trans(s, sNext) += 1.0; - } - } - ++validSeqs; - } catch (...) { + const double lp = accum_sequence(hmm, obs, pi, trans, emisData); + if (!std::isfinite(lp)) continue; - } + totalLogProb += lp; + ++validSeqs; } if (validSeqs == 0) return lastLogProb_; - // Normalise pi - { - double piSum = 0.0; - for (std::size_t i = 0; i < N; ++i) - piSum += pi(i); - if (piSum > 0.0) { - for (std::size_t i = 0; i < N; ++i) - pi(i) /= piSum; - } else { - for (std::size_t i = 0; i < N; ++i) - pi(i) = 1.0 / static_cast(N); - } - hmm.setPi(pi); - } - - // Normalise transition rows - for (std::size_t i = 0; i < N; ++i) { - double rowSum = 0.0; - for (std::size_t j = 0; j < N; ++j) - rowSum += trans(i, j); - if (rowSum > 0.0) { - for (std::size_t j = 0; j < N; ++j) - trans(i, j) /= rowSum; - } else { - for (std::size_t j = 0; j < N; ++j) - trans(i, j) = 1.0 / static_cast(N); - } - } - hmm.setTrans(trans); - - // Refit emission distributions (unweighted — Viterbi hard assignment) - for (std::size_t i = 0; i < N; ++i) { - const std::size_t M = emisData[i].size(); - if (M == 0) { - hmm.getDistribution(i).reset(); - continue; - } - hmm.getDistribution(i).fit(std::span(emisData[i].data(), M)); - } - + normalize_and_commit(hmm, N, pi, trans); + apply_emission_fits(hmm, N, emisData); return totalLogProb; } diff --git a/tests/distributions/test_chi_squared_distribution.cpp b/tests/distributions/test_chi_squared_distribution.cpp index 6180126..1e7a562 100644 --- a/tests/distributions/test_chi_squared_distribution.cpp +++ b/tests/distributions/test_chi_squared_distribution.cpp @@ -234,23 +234,6 @@ void testStringRepresentation() { assert(str.find("3.5") != std::string::npos); assert(str.find("k") != std::string::npos); - // Test fromString functionality - ChiSquaredDistribution chi_dist_from_str = - ChiSquaredDistribution::fromString("ChiSquared(k=7.25)"); - assert(std::abs(chi_dist_from_str.getDegreesOfFreedom() - 7.25) < 1e-10); - - ChiSquaredDistribution chi_dist_from_str2 = - ChiSquaredDistribution::fromString("ChiSquared(df=12.0)"); - assert(std::abs(chi_dist_from_str2.getDegreesOfFreedom() - 12.0) < 1e-10); - - // Test invalid string formats - try { - ChiSquaredDistribution::fromString("InvalidFormat"); - assert(false); // Should not reach here - } catch (const std::invalid_argument &) { - // Expected behavior - } - std::cout << "String representation: " << str << std::endl; std::cout << "✓ String representation tests passed" << std::endl; } diff --git a/tests/distributions/test_student_t_distribution.cpp b/tests/distributions/test_student_t_distribution.cpp index 2f9971b..599c51f 100644 --- a/tests/distributions/test_student_t_distribution.cpp +++ b/tests/distributions/test_student_t_distribution.cpp @@ -218,21 +218,6 @@ void testStringRepresentation() { assert(str.find("3.5") != std::string::npos); assert(str.find("nu") != std::string::npos); - // Test fromString functionality - StudentTDistribution t_dist_from_str = StudentTDistribution::fromString("StudentT(ν=7.25)"); - assert(std::abs(t_dist_from_str.getDegreesOfFreedom() - 7.25) < 1e-10); - - StudentTDistribution t_dist_from_str2 = StudentTDistribution::fromString("StudentT(df=12.0)"); - assert(std::abs(t_dist_from_str2.getDegreesOfFreedom() - 12.0) < 1e-10); - - // Test invalid string formats - try { - StudentTDistribution::fromString("InvalidFormat"); - assert(false); // Should not reach here - } catch (const std::invalid_argument &) { - // Expected behavior - } - std::cout << "String representation: " << str << std::endl; std::cout << "✓ String representation tests passed" << std::endl; } diff --git a/tools/batch_performance.cpp b/tools/batch_performance.cpp index a26a70e..41481d5 100644 --- a/tools/batch_performance.cpp +++ b/tools/batch_performance.cpp @@ -13,7 +13,7 @@ #include "libhmm/hmm.h" #include "libhmm/calculators/forward_backward_calculator.h" #include "libhmm/calculators/viterbi_calculator.h" -#include "libhmm/distributions/gaussian_distribution.h" +#include "libhmm/distributions/distributions.h" #include #include #include diff --git a/tools/bw_hotspot.cpp b/tools/bw_hotspot.cpp index 7109b2a..4c871c1 100644 --- a/tools/bw_hotspot.cpp +++ b/tools/bw_hotspot.cpp @@ -18,7 +18,7 @@ #include "libhmm/calculators/forward_backward_calculator.h" #include "libhmm/hmm.h" #include "libhmm/distributions/discrete_distribution.h" -#include "libhmm/distributions/gaussian_distribution.h" +#include "libhmm/distributions/distributions.h" #include "libhmm/performance/transcendental_kernels.h" #include diff --git a/tools/fb_contour_sweep.cpp b/tools/fb_contour_sweep.cpp index 5b10626..b875fde 100644 --- a/tools/fb_contour_sweep.cpp +++ b/tools/fb_contour_sweep.cpp @@ -1,5 +1,5 @@ #include "libhmm/hmm.h" -#include "libhmm/distributions/gaussian_distribution.h" +#include "libhmm/distributions/distributions.h" #include #include #include diff --git a/tools/fb_crossover_sweep.cpp b/tools/fb_crossover_sweep.cpp index 10d6e14..c1cc299 100644 --- a/tools/fb_crossover_sweep.cpp +++ b/tools/fb_crossover_sweep.cpp @@ -8,7 +8,7 @@ #include "libhmm/performance/fb_recurrence_policy.h" #include "libhmm/calculators/forward_backward_calculator.h" -#include "libhmm/distributions/gaussian_distribution.h" +#include "libhmm/distributions/distributions.h" #include "libhmm/hmm.h" #include "libhmm/platform/simd_platform.h" diff --git a/tools/hotspot_breakdown.cpp b/tools/hotspot_breakdown.cpp index 7e59c40..24827bc 100644 --- a/tools/hotspot_breakdown.cpp +++ b/tools/hotspot_breakdown.cpp @@ -1,5 +1,5 @@ #include "libhmm/hmm.h" -#include "libhmm/distributions/gaussian_distribution.h" +#include "libhmm/distributions/distributions.h" #include "libhmm/math/constants.h" #include #include diff --git a/tools/simd_inspection.cpp b/tools/simd_inspection.cpp index 2d33463..82913cd 100644 --- a/tools/simd_inspection.cpp +++ b/tools/simd_inspection.cpp @@ -12,7 +12,7 @@ #include "libhmm/hmm.h" #include "libhmm/calculators/forward_backward_calculator.h" #include "libhmm/calculators/viterbi_calculator.h" -#include "libhmm/distributions/gaussian_distribution.h" +#include "libhmm/distributions/distributions.h" #include "libhmm/platform/simd_platform.h" #include #include