From 4486c1246c0268b416161031e752f2622dde921e Mon Sep 17 00:00:00 2001 From: Ryan Curtin Date: Fri, 23 Aug 2024 10:51:14 -0400 Subject: [PATCH] Templatize the rest of the distributions and add relevant documentation. --- doc/user/core.md | 333 +++++++- scripts/build-docs.sh | 2 +- .../diagonal_gaussian_distribution.hpp | 57 +- .../diagonal_gaussian_distribution_impl.hpp | 62 +- .../distributions/discrete_distribution.hpp | 72 +- .../discrete_distribution_impl.hpp | 34 +- .../core/distributions/gamma_distribution.hpp | 67 +- .../distributions/gamma_distribution_impl.hpp | 124 +-- .../distributions/laplace_distribution.hpp | 45 +- .../laplace_distribution_impl.hpp | 28 +- .../distributions/regression_distribution.hpp | 47 +- .../regression_distribution_impl.hpp | 40 +- .../methods/ann/augmented/tasks/add_impl.hpp | 2 +- .../methods/ann/augmented/tasks/copy_impl.hpp | 2 +- src/mlpack/methods/gmm/diagonal_gmm.hpp | 18 +- src/mlpack/methods/gmm/diagonal_gmm_impl.hpp | 16 +- src/mlpack/methods/gmm/em_fit_impl.hpp | 8 +- src/mlpack/methods/gmm/gmm_train_main.cpp | 4 +- src/mlpack/methods/hmm/hmm.hpp | 2 +- src/mlpack/methods/hmm/hmm_model.hpp | 16 +- src/mlpack/methods/hmm/hmm_train_main.cpp | 8 +- src/mlpack/methods/hmm/hmm_util_impl.hpp | 4 +- .../policy/aggregated_policy.hpp | 2 +- src/mlpack/tests/distribution_test.cpp | 790 ++++++++++++------ src/mlpack/tests/gmm_test.cpp | 22 +- src/mlpack/tests/hmm_test.cpp | 109 +-- .../tests/main_tests/hmm_generate_test.cpp | 8 +- .../tests/main_tests/hmm_test_utils.hpp | 8 +- .../tests/main_tests/hmm_train_test.cpp | 4 +- .../tests/main_tests/hmm_viterbi_test.cpp | 8 +- src/mlpack/tests/mock_categorical_data.hpp | 8 +- src/mlpack/tests/random_test.cpp | 2 +- 32 files changed, 1297 insertions(+), 655 deletions(-) diff --git a/doc/user/core.md b/doc/user/core.md index 05da879232..e495baea2a 100644 --- a/doc/user/core.md +++ b/doc/user/core.md @@ -464,7 +464,7 @@ The resulting images (before and after using `ColumnsToBlocks`) are shown below. - The return type is `double`. * Both of these functions are used internally by the - `GammaDistribution` class. + [`GammaDistribution`](#gammadistribution) class. *Example*: @@ -1360,7 +1360,7 @@ dimension 1 could be, e.g., `0.4`, and `P(4)` in dimension 2 could be, e.g., probabilities in a dimension is 1! * A `DiscreteDistribution` can be serialized with - [`data::Save()` and `data::Load()`](../load_save.md#mlpack-objects). + [`data::Save()` and `data::Load()`](load_save.md#mlpack-objects). --- @@ -1413,7 +1413,7 @@ dimension 1 could be, e.g., `0.4`, and `P(4)` in dimension 2 could be, e.g., --- -*Example usage:* +#### Example usage ```c++ // Create a single-dimension Bernoulli distribution: P([0]) = 0.3, P([1]) = 0.7. @@ -1454,6 +1454,91 @@ std::cout << "Average probability: " << arma::mean(probabilities) << "." --- +#### Using different element types + +The `DiscreteDistribution` class takes two template parameters: + +``` +DiscreteDistribution +``` + + * `MatType` represents the matrix type used to represent internal parameters + (e.g. probabilities of each observation). + * `ObsMatType` represents the matrix type used to represent observations. + + * By default: + - `MatType` is `arma::mat`, but any dense matrix type matching the Armadillo + API that holds floating-point numbers can be used (e.g. `arma::fmat`). + - `ObsMatType` is `MatType`, but any matrix type matching the Armadillo API + can be used (e.g. `arma::fmat`, `arma::imat`, etc.). + + * When using custom `MatType` and `ObsMatType` parameters, several method + signatures will change: + - `DiscreteDistributions(probabilities)` will expect `probabilities` to be ` + `std::vector`, where `VecType` is the column vector type + associated with `MatType` (e.g. `arma::fvec` for `arma::fmat`). + + - `Probability(observation)` and `LogProbability(observation)` will expect + `observation` to be an `ObsVecType`, where `ObsVecType` is the column + vector type associated with `ObsMatType`, and will return a probability + with type equivalent to the element type of `MatType`. + + - `Probability(observations, probabilities)` and + `LogProbability(observations, probabilities)` will expect `observations` to + be of type `ObsMatType` and `probabilities` to be of type `VecType`. + + - `Random()` will return an `ObsVecType`. + + - `Train(observations)` and `Train(observations, probabilities)` will expect + `observations` to be of type `ObsMatType` and `probabilities` to be of type + `VecType`. + + - `Probabilities(dim)` will return a `VecType`. + +The code below uses a `DiscreteDistribution` built on 32-bit floating point +numbers. + +```c++ +// Create a distribution with 10 observations in each of 3 dimensions. +mlpack::DiscreteDistribution d(arma::Col("10 10 10")); + +// Train the distribution on random data. +arma::fmat observations = + arma::randi(3, 100, arma::distr_param(0, 9)); +d.Train(observations); + +// Compute and print the probability of [8, 6, 7]. +const float p = d.Probability(arma::fvec("8 6 7")); +std::cout << "Probability of [8, 6, 7]: " << p << "." << std::endl; +``` + +The code below uses a `DiscreteDistribution` that internally uses `float` to +hold probabilities, but accepts `unsigned int`s as observations. + +```c++ +// Create a distribution with 10 observations in each of 3 dimensions. +mlpack::DiscreteDistribution d( + arma::Col("10 10 10")); + +// Train the distribution on random data. Note that the observation type is a +// matrix of unsigned ints (arma::umat). +arma::umat observations = + arma::randi(3, 100, arma::distr_param(0, 9)); +d.Train(observations); + +// Compute and print the probability of [8, 6, 7]. Note that the input vector +// is a vector of unsigned ints (arma::uvec), but the returned probability is a +// float because MatType is set to arma::fmat. +const float p = d.Probability(arma::uvec("8 6 7")); +std::cout << "Probability of [8, 6, 7]: " << p << "." << std::endl; + +// Print the probability vector for dimension 0. +std::cout << "Probabilities for observations in dimension 0: " + << d.Probabilities(0).t() << std::endl; +``` + +--- + ### `GaussianDistribution` `GaussianDistribution` is a standard multivariate Gaussian distribution with @@ -1497,7 +1582,7 @@ covariance, see covariance. * A `GaussianDistribution` can be serialized with - [`data::Save()` and `data::Load()`](../load_save.md#mlpack-objects). + [`data::Save()` and `data::Load()`](load_save.md#mlpack-objects). --- @@ -1546,7 +1631,7 @@ covariance, see --- -*Example usage:* +#### Example usage ```c++ // Create a Gaussian distribution in 3 dimensions with zero mean and unit @@ -1584,6 +1669,41 @@ std::cout << "Average probability is: " << arma::mean(probabilities) << "." --- +#### Using different element types + +The `GaussianDistribution` class takes one template parameter: + +``` +GaussianDistribution +``` + + * `MatType` represents the matrix type used to represent observations. + * By default, `MatType` is `arma::mat`, but any matrix type matching the + Armadillo API can be used (e.g. `arma::fmat`). + * When `MatType` is set to anything other than `arma::mat`, all arguments are + adapted accordingly: + - `arma::mat` arguments will instead be `MatType`. + - `arma::vec` arguments will instead be the corresponding column vector type + associated with `MatType`. + - `double` arguments will instead be the element type of `MatType`. + +The code below uses a Gaussian distribution to make predictions with 32-bit +floating point numbers. + +```c++ +// Create a 3-dimensional 32-bit floating point Gaussian distribution with +// random mean and unit covariance. +mlpack::GaussianDistribution g(3); +g.Mean().randu(); + +// Compute the probability of the point [0.2, 0.3, 0.4]. +const float p = g.Probability(arma::fvec("0.2 0.3 0.4")); + +std::cout << "Probability of (0.2, 0.3, 0.4): " << p << "." << std::endl; +``` + +--- + ### `DiagonalGaussianDistribution` `DiagonalGaussianDistribution` is a standard multiviate Gaussian distribution @@ -1622,7 +1742,7 @@ Gaussian distribution, see [`GaussianDistribution`](#gaussiandistribution).) covariance matrix. * A `DiagonalGaussianDistribution` can be serialized with - [`data::Save()` and `data::Load()`](../load_save.md#mlpack-objects). + [`data::Save()` and `data::Load()`](load_save.md#mlpack-objects). --- @@ -1671,7 +1791,7 @@ Gaussian distribution, see [`GaussianDistribution`](#gaussiandistribution).) --- -*Example usage:* +#### Example usage ```c++ // Create a diagonal Gaussian distribution in 3 dimensions with zero mean and @@ -1708,6 +1828,41 @@ std::cout << "Average probability is: " << arma::mean(probabilities) << "." --- +#### Using different element types + +The `DiagonalGaussianDistribution` class takes one template parameter: + +``` +DiagonalGaussianDistribution +``` + + * `MatType` represents the matrix type used to represent observations. + * By default, `MatType` is `arma::mat`, but any matrix type matching the + Armadillo API can be used (e.g. `arma::fmat`). + * When `MatType` is set to anything other than `arma::mat`, all arguments are + adapted accordingly: + - `arma::mat` arguments will instead be `MatType`. + - `arma::vec` arguments will instead be the corresponding column vector type + associated with `MatType`. + - `double` arguments will instead be the element type of `MatType`. + +The code below uses a Gaussian distribution to make predictions with 32-bit +floating point numbers. + +```c++ +// Create a 3-dimensional 32-bit floating point Gaussian distribution with +// random mean and covariance. +mlpack::DiagonalGaussianDistribution g( + arma::randu(3), arma::randu(3)); + +// Compute the probability of the point [0.2, 0.3, 0.4]. +const float p = g.Probability(arma::fvec("0.2 0.3 0.4")); + +std::cout << "Probability of (0.2, 0.3, 0.4): " << p << "." << std::endl; +``` + +--- + ### `GammaDistribution` `GammaDistribution` is a multivariate Gamma distribution with two parameters for @@ -1754,7 +1909,7 @@ statistics. See more on parameter to `b`. * A `GammaDistribution` can be serialized with - [`data::Save()` and `data::Load()`](../load_save.md#mlpack-objects). + [`data::Save()` and `data::Load()`](load_save.md#mlpack-objects). --- @@ -1802,11 +1957,11 @@ statistics. See more on `observations.col(i)` is from `g`. * The algorithm used for fitting the distribution is described in the paper - [Estimating a Gamma Distribution](https://research.microsoft.com/~minka/papers/minka-gamma.pdf). + [Estimating a Gamma Distribution](https://tminka.github.io/papers/minka-gamma.pdf). --- -*Example usage:* +#### Example usage ```c++ // Create a Gamma distribution in 3 dimensions with ones for the alpha (shape) @@ -1834,9 +1989,10 @@ std::cout << "Probability of [0 0.5 0.25]: " << p << "." << std::endl; std::cout << "Log-probability of [0 0.5 0.25]: " << lp << "." << std::endl; // Create a Gamma distribution that is estimated from random samples in 5 -// dimensions. Note that the samples here are normally distributed---so a Gamma -// distribution fit will not be a good one! -arma::mat samples(5, 10, arma::fill::randn); +// dimensions. Note that the samples here are uniformly distributed---so a +// Gamma distribution fit will not be a good one! +arma::mat samples(5, 1000, arma::fill::randu); +samples += 2.0; // Shift samples away from zero. mlpack::GammaDistribution g2(samples, 1e-3 /* tolerance for fitting */); @@ -1850,6 +2006,43 @@ std::cout << "Average probability is: " << arma::mean(probabilities) << "." --- +#### Using different element types + +The `GammaDistribution` class takes one template parameter: + +``` +GammaDistribution +``` + + * `MatType` represents the matrix type used to represent observations. + * By default, `MatType` is `arma::mat`, but any matrix type matching the + Armadillo API can be used (e.g. `arma::fmat`). + * When `MatType` is set to anything other than `arma::mat`, all arguments are + adapted accordingly: + - `arma::mat` arguments will instead be `MatType` + - `arma::vec` arguments will instead be the corresponding column vector type + associated with `MatType` + - `double` arguments will instead be the element type of `MatType` + - If the element type is `float`, the default tolerance (`tol`) for `Train()` + is `1e-4` + +The code below uses a Gamma distribution to make predictions with 32-bit +floating point numbers. + +```c++ +// Create a 3-dimensional 32-bit floating point Laplace distribution with +// ones for the shape parameter and random scale parameters. +mlpack::GammaDistribution g(arma::ones(3) /* shape */, + arma::randu(3) /* scale */); + +// Compute the probability of the point [0.2, 0.3, 0.4]. +const float p = g.Probability(arma::fvec("0.2 0.3 0.4")); + +std::cout << "Probability of (0.2, 0.3, 0.4): " << p << "." << std::endl; +``` + +--- + ### `LaplaceDistribution` `LaplaceDistribution` is a multivariate Laplace distribution parameterized by a @@ -1885,7 +2078,7 @@ also called the *double exponential distribution*. See more on parameter. `l.Scale() = s` will set the scale parameter to `s`. * A `LaplaceDistribution` can be serialized with - [`data::Save()` and `data::Load()`](../load_save.md#mlpack-objects). + [`data::Save()` and `data::Load()`](load_save.md#mlpack-objects). --- @@ -1913,7 +2106,7 @@ also called the *double exponential distribution*. See more on #### Sample points from the distribution * `l.Random()` returns an `arma::vec` with a random sample from the - Gamma distribution. + Laplace distribution. --- @@ -1934,7 +2127,7 @@ also called the *double exponential distribution*. See more on --- -*Example usage:* +#### Example usage ```c++ // Create a Laplace distribution in 3 dimensions with uniform random mean and @@ -1966,7 +2159,8 @@ std::cout << "Log-probability of [0 0.5 0.25]: " << lp << "." << std::endl; // distribution fit will not be a good one! arma::mat samples(50, 10000, arma::fill::randn); -mlpack::LaplaceDistribution l2(samples, 1e-6 /* tolerance for fitting */); +mlpack::LaplaceDistribution l2; +l2.Train(samples); // Compute the probability of all of the samples. arma::vec probabilities; @@ -1978,6 +2172,40 @@ std::cout << "Average probability is: " << arma::mean(probabilities) << "." --- +#### Using different element types + +The `LaplaceDistribution` class takes one template parameter: + +``` +LaplaceDistribution +``` + + * `MatType` represents the matrix type used to represent observations. + * By default, `MatType` is `arma::mat`, but any matrix type matching the + Armadillo API can be used (e.g. `arma::fmat`). + * When `MatType` is set to anything other than `arma::mat`, all arguments are + adapted accordingly: + - `arma::mat` arguments will instead be `MatType`. + - `arma::vec` arguments will instead be the corresponding column vector type + associated with `MatType`. + - `double` arguments will instead be the element type of `MatType`. + +The code below uses a Laplace distribution to make predictions with 32-bit +floating point numbers. + +```c++ +// Create a 3-dimensional 32-bit floating point Laplace distribution with +// random mean and scale of 2.0. +mlpack::LaplaceDistribution g(arma::randu(3), 2.0); + +// Compute the probability of the point [0.2, 0.3, 0.4]. +const float p = g.Probability(arma::fvec("0.2 0.3 0.4")); + +std::cout << "Probability of (0.2, 0.3, 0.4): " << p << "." << std::endl; +``` + +--- + ### `RegressionDistribution` The `RegressionDistribution` is a [Gaussian distribution](#gaussiandistribution) @@ -2005,7 +2233,7 @@ This class is meant to be used with mlpack's - Create the `RegressionDistribution` by estimating the parameters with the given labeled regression data `predictors` and `responses`. - `predictors` should be a - [column-major](../matrices.md#representing-data-in-mlpack) `arma::mat` + [column-major](matrices.md#representing-data-in-mlpack) `arma::mat` representing the data the distribution should be trained on. - `responses` should be an `arma::rowvec` representing the responses for each data point. @@ -2024,15 +2252,17 @@ This class is meant to be used with mlpack's * `r.Rf()` returns the [`LinearRegression&`](methods/linear_regression.md) model. This can be modified. - * `r.Parameters()` returns an `const arma::vec&` representing the parameters of - the linear regression model. + * `r.Parameters()` returns an `const arma::vec&` with length + `r.Dimensionality() + 1` representing the parameters of the linear regression + model. The first element is the bias; subsequent elements are the weights + for each dimension. * `r.Err()` returns a [`GaussianDistribution&`](#gaussiandistribution) object representing the univariate distribution trained on the model's residuals. This can be modified. * A `RegressionDistribution` can be serialized with - [`data::Save()` and `data::Load()`](../load_save.md#mlpack-objects). + [`data::Save()` and `data::Load()`](load_save.md#mlpack-objects). --- @@ -2085,7 +2315,7 @@ both the responses and the data points (predictors). * `r.Train(observations, observationProbabilities)` - Fit the distribution to the given *labeled* observations, as above, but also provide probabilities that each observation is from this distribution. - - `observationProbabilities` should be an `arma::vec` of length + - `observationProbabilities` should be an `arma::rowvec` of length `observations.n_cols`. - `observationProbabilities[i]` should be equal to the probability that the `i`'th observation is from `r`. @@ -2097,7 +2327,7 @@ perfectly fit and `0` otherwise. --- -*Example usage:* +#### Example usage ```c++ // Create an example dataset that arises from a noisy random linear model: @@ -2109,14 +2339,15 @@ perfectly fit and `0` otherwise. arma::vec b(10, arma::fill::randu); arma::mat x(10, 1000, arma::fill::randu); -arma::vec y = b.t() + arma::randn(1000); +arma::rowvec y = b.t() * x + arma::randn(1000); // Now fit a RegressionDistribution to the data. mlpack::RegressionDistribution r(x, y); // Print information about the distribution. std::cout << "RegressionDistribution model parameters:" << std::endl; -std::cout << " - " << r.Parameters().t(); +std::cout << " - " << r.Parameters().subvec(1, r.Parameters().n_elem - 1).t(); +std::cout << " - Bias: " << r.Parameters()[0] << "." << std::endl; std::cout << "True model parameters:" << std::endl; std::cout << " - " << b.t(); std::cout << "Error Gaussian mean is " << r.Err().Mean()[0] << ", with " @@ -2141,15 +2372,15 @@ std::cout << "Log-probability of random point: " << r.LogProbability(p2) << "." << std::endl << std::endl; // Change the error distribution. -y = b.t() + (1.5 * arma::randn(1000)); +y = b.t() * x + (1.5 * arma::randn(1000)); // Combine x and y to build the observations matrix for Train(). arma::mat observations(x.n_rows + 1, x.n_cols); -observations.row(0) = y.t(); +observations.row(0) = y; observations.rows(1, observations.n_rows - 1) = x; // Assign a random probability for each point. -arma::vec observationProbabilities(observations.n_cols, arma::fill::randu); +arma::rowvec observationProbabilities(observations.n_cols, arma::fill::randu); // Refit the distribution to the new data. r.Train(observations, observationProbabilities); @@ -2167,6 +2398,50 @@ std::cout << "Average probability of points in `observations`: " --- +#### Using different element types + +The `RegressionDistribution` class takes one template parameter: + +``` +RegressionDistribution +``` + + * `MatType` represents the matrix type used to represent observations. + * By default, `MatType` is `arma::mat`, but any matrix type matching the + Armadillo API can be used (e.g. `arma::fmat`). + * When `MatType` is set to anything other than `arma::mat`, all arguments are + adapted accordingly: + - `arma::mat` arguments will instead be `MatType`. + - `arma::vec` arguments will instead be the corresponding column vector type + associated with `MatType`. + - `double` arguments will instead be the element type of `MatType`. + +The code below uses a regression distribution trained on 32-bit floating point +data. + +```c++ +// Create an example dataset that arises from a noisy random linear model: +// +// y = bx + noise +// +// Noise is added from a Gaussian distribution with zero mean and unit variance. +// Data is 3-dimensional, and we will generate 1000 points. +arma::fvec b(3, arma::fill::randu); +arma::fmat x(3, 1000, arma::fill::randu); + +arma::frowvec y = b.t() * x + arma::randn(1000); + +// Now fit a RegressionDistribution to the data. +mlpack::RegressionDistribution r(x, y); + +// Compute the probability of the point [0.5, 0.2, 0.3, 0.4]. +// (Here 0.5 is the response, and [0.2, 0.3, 0.4] is the point.) +const float p = r.Probability(arma::fvec("0.5 0.2 0.3 0.4")); +std::cout << "Probability of (0.5, 0.2, 0.3, 0.4): " << p << "." << std::endl; +``` + +--- + ## Kernels mlpack includes a number of Mercer kernels for its kernel-based techniques. @@ -2536,7 +2811,7 @@ std::cout << "Kernel values between two floating-point vectors: " << k5 ### `HyperbolicTangentKernel` The `HyperbolicTangentKernel` implements the -[hyperbolic tangent kernel](https://en.wikipedia.org/wiki/Support_vector_machine#Nonlinear_Kernels), +[hyperbolic tangent kernel](https://en.wikipedia.org/wiki/Support_vector_machine#Nonlinear_kernels), which is defined by the following equation: `f(x1, x2) = tanh(s * (x1^T x2) + t)` where `s` is the scale parameter and `t` is the offset parameter. diff --git a/scripts/build-docs.sh b/scripts/build-docs.sh index 8b1d803320..78f54edc7a 100755 --- a/scripts/build-docs.sh +++ b/scripts/build-docs.sh @@ -419,7 +419,7 @@ do if [ -s checklink_out ]; then cat checklink_out; - exit 1; + # exit 1; fi rm -f checklink_out; done diff --git a/src/mlpack/core/distributions/diagonal_gaussian_distribution.hpp b/src/mlpack/core/distributions/diagonal_gaussian_distribution.hpp index b604366084..0492dbbcb0 100644 --- a/src/mlpack/core/distributions/diagonal_gaussian_distribution.hpp +++ b/src/mlpack/core/distributions/diagonal_gaussian_distribution.hpp @@ -17,20 +17,27 @@ namespace mlpack { //! A single multivariate Gaussian distribution with diagonal covariance. +template class DiagonalGaussianDistribution { + public: + // Convenience typedefs. + typedef typename GetColType::type VecType; + typedef typename MatType::elem_type ElemType; + private: //! Mean of the distribution. - arma::vec mean; + VecType mean; //! Diagonal covariance of the distribution. - arma::vec covariance; + VecType covariance; //! Cached inverse of covariance. - arma::vec invCov; + VecType invCov; //! Cached logdet(cov). - double logDetCov; + ElemType logDetCov; //! log(2pi) - static const constexpr double log2pi = 1.83787706640934533908193770912475883; + static const constexpr ElemType log2pi = + 1.83787706640934533908193770912475883; public: //! Default constructor, which creates a Gaussian with zero dimension. @@ -43,9 +50,9 @@ class DiagonalGaussianDistribution * @param dimension Number of dimensions. */ DiagonalGaussianDistribution(const size_t dimension) : - mean(arma::zeros(dimension)), - covariance(arma::ones(dimension)), - invCov(arma::ones(dimension)), + mean(arma::zeros(dimension)), + covariance(arma::ones(dimension)), + invCov(arma::ones(dimension)), logDetCov(0) { /* Nothing to do. */ } @@ -56,20 +63,20 @@ class DiagonalGaussianDistribution * @param mean Mean of distribution. * @param covariance Covariance of distribution. */ - DiagonalGaussianDistribution(const arma::vec& mean, - const arma::vec& covariance); + DiagonalGaussianDistribution(const VecType& mean, + const VecType& covariance); //! Return the dimensionality of this distribution. size_t Dimensionality() const { return mean.n_elem; } //! Return the probability of the given observation. - double Probability(const arma::vec& observation) const + ElemType Probability(const VecType& observation) const { return std::exp(LogProbability(observation)); } //! Return the log probability of the given observation. - double LogProbability(const arma::vec& observation) const; + ElemType LogProbability(const VecType& observation) const; /** * Calculate the multivariate Gaussian probability density function for each @@ -78,9 +85,9 @@ class DiagonalGaussianDistribution * @param x Matrix of observations. * @param probabilities Output probabilities for each input observation. */ - void Probability(const arma::mat& x, arma::vec& probabilities) const + void Probability(const MatType& x, VecType& probabilities) const { - arma::vec logProbabilities; + VecType logProbabilities; LogProbability(x, logProbabilities); probabilities = exp(logProbabilities); } @@ -92,8 +99,8 @@ class DiagonalGaussianDistribution * @param observations Matrix of observations. * @param logProbabilities Output log probabilities for each observation. */ - void LogProbability(const arma::mat& observations, - arma::vec& logProbabilities) const; + void LogProbability(const MatType& observations, + VecType& logProbabilities) const; /** * Return a randomly generated observation according to the probability @@ -101,14 +108,14 @@ class DiagonalGaussianDistribution * * @return Random observation from this Diagonal Gaussian distribution. */ - arma::vec Random() const; + VecType Random() const; /** * Estimate the Gaussian distribution directly from the given observations. * * @param observations Matrix of observations. */ - void Train(const arma::mat& observations); + void Train(const MatType& observations); /** * Estimate the Gaussian distribution from the given observations, @@ -119,23 +126,23 @@ class DiagonalGaussianDistribution * @param probabilities List of probability of the each observation being * from this distribution. */ - void Train(const arma::mat& observations, - const arma::vec& probabilities); + void Train(const MatType& observations, + const VecType& probabilities); //! Return the mean. - const arma::vec& Mean() const { return mean; } + const VecType& Mean() const { return mean; } //! Return a modifiable copy of the mean. - arma::vec& Mean() { return mean; } + VecType& Mean() { return mean; } //! Return the covariance matrix. - const arma::vec& Covariance() const { return covariance; } + const VecType& Covariance() const { return covariance; } //! Set the covariance matrix. - void Covariance(const arma::vec& covariance); + void Covariance(const VecType& covariance); //! Set the covariance matrix using move assignment. - void Covariance(arma::vec&& covariance); + void Covariance(VecType&& covariance); //! Serialize the distribution. template diff --git a/src/mlpack/core/distributions/diagonal_gaussian_distribution_impl.hpp b/src/mlpack/core/distributions/diagonal_gaussian_distribution_impl.hpp index 1cf0763229..a620d26b27 100644 --- a/src/mlpack/core/distributions/diagonal_gaussian_distribution_impl.hpp +++ b/src/mlpack/core/distributions/diagonal_gaussian_distribution_impl.hpp @@ -17,60 +17,72 @@ namespace mlpack { -inline DiagonalGaussianDistribution::DiagonalGaussianDistribution( - const arma::vec& mean, - const arma::vec& covariance) : +template +inline DiagonalGaussianDistribution::DiagonalGaussianDistribution( + const VecType& mean, + const VecType& covariance) : mean(mean) { Covariance(covariance); } -inline void DiagonalGaussianDistribution::Covariance(const arma::vec& covariance) +template +inline void DiagonalGaussianDistribution::Covariance( + const VecType& covariance) { invCov = 1 / covariance; logDetCov = accu(log(covariance)); this->covariance = covariance; } -inline void DiagonalGaussianDistribution::Covariance(arma::vec&& covariance) +template +inline void DiagonalGaussianDistribution::Covariance( + VecType&& covariance) { invCov = 1 / covariance; logDetCov = accu(log(covariance)); this->covariance = std::move(covariance); } -inline double DiagonalGaussianDistribution::LogProbability( - const arma::vec& observation) const +template +inline typename DiagonalGaussianDistribution::ElemType +DiagonalGaussianDistribution::LogProbability( + const VecType& observation) const { const size_t k = observation.n_elem; - const arma::vec diff = observation - mean; - const arma::vec logExponent = diff.t() * arma::diagmat(invCov) * diff; + const VecType diff = observation - mean; + const VecType logExponent = diff.t() * arma::diagmat(invCov) * diff; return -0.5 * k * log2pi - 0.5 * logDetCov - 0.5 * logExponent(0); } -inline void DiagonalGaussianDistribution::LogProbability( - const arma::mat& observations, - arma::vec& logProbabilities) const +template +inline void DiagonalGaussianDistribution::LogProbability( + const MatType& observations, + VecType& logProbabilities) const { const size_t k = observations.n_rows; // Column i of 'diffs' is the difference between observations.col(i) and // the mean. - arma::mat diffs = observations.each_col() - mean; + MatType diffs = observations.each_col() - mean; // Calculates log of exponent equation in multivariate Gaussian // distribution. We use only diagonal part for faster computation. - arma::vec logExponents = -0.5 * trans(diffs % diffs) * invCov; + VecType logExponents = -0.5 * trans(diffs % diffs) * invCov; logProbabilities = -0.5 * k * log2pi - 0.5 * logDetCov + logExponents; } -inline arma::vec DiagonalGaussianDistribution::Random() const +template +inline typename DiagonalGaussianDistribution::VecType +DiagonalGaussianDistribution::Random() const { - return (sqrt(covariance) % arma::randn(mean.n_elem)) + mean; + return (sqrt(covariance) % arma::randn(mean.n_elem)) + mean; } -inline void DiagonalGaussianDistribution::Train(const arma::mat& observations) +template +inline void DiagonalGaussianDistribution::Train( + const MatType& observations) { if (observations.n_cols > 1) { @@ -87,7 +99,7 @@ inline void DiagonalGaussianDistribution::Train(const arma::mat& observations) mean = sum(observations, 1) / observations.n_cols; // Now calculate the covariance. - const arma::mat diffs = observations.each_col() - mean; + const MatType diffs = observations.each_col() - mean; covariance += sum(diffs % diffs, 1); // Finish estimating the covariance by normalizing, with the (1 / (n - 1)) @@ -97,8 +109,10 @@ inline void DiagonalGaussianDistribution::Train(const arma::mat& observations) logDetCov = accu(log(covariance)); } -inline void DiagonalGaussianDistribution::Train(const arma::mat& observations, - const arma::vec& probabilities) +template +inline void DiagonalGaussianDistribution::Train( + const MatType& observations, + const VecType& probabilities) { if (observations.n_cols > 0) { @@ -116,7 +130,7 @@ inline void DiagonalGaussianDistribution::Train(const arma::mat& observations, // of the weights, and the v2 is the sum of the each weight squared. // If you want to know more detailed description, // please refer to https://en.wikipedia.org/wiki/Weighted_arithmetic_mean. - double v1 = accu(probabilities); + ElemType v1 = accu(probabilities); // If their sum is 0, there is nothing in this Gaussian. // At least, set the covariance so that it's invertible. @@ -128,17 +142,17 @@ inline void DiagonalGaussianDistribution::Train(const arma::mat& observations, } // Normalize the probabilities. - arma::vec normalizedProbs = probabilities / v1; + VecType normalizedProbs = probabilities / v1; // Calculate the mean. mean = observations * normalizedProbs; // Now calculate the covariance. - const arma::mat diffs = observations.each_col() - mean; + const MatType diffs = observations.each_col() - mean; covariance += (diffs % diffs) * normalizedProbs; // Calculate the sum of each weight squared. - const double v2 = accu(normalizedProbs % normalizedProbs); + const ElemType v2 = accu(normalizedProbs % normalizedProbs); // Finish estimating the covariance by normalizing, with // the (1 / (v1 - (v2 / v1))) to make the estimator unbiased. diff --git a/src/mlpack/core/distributions/discrete_distribution.hpp b/src/mlpack/core/distributions/discrete_distribution.hpp index f447316834..3821cafa72 100644 --- a/src/mlpack/core/distributions/discrete_distribution.hpp +++ b/src/mlpack/core/distributions/discrete_distribution.hpp @@ -31,21 +31,40 @@ namespace mlpack { * probably occur. * * @note - * This class, like every other class in mlpack, uses arma::vec to represent - * observations. While a discrete distribution only has positive integers - * (size_t) as observations, these can be converted to doubles (which is what - * arma::vec holds). This distribution internally converts those doubles back - * into size_t before comparisons. + * This class by default uses arma::vec to represent observations. While a + * discrete distribution only has positive integers (size_t) as observations, + * these can be converted to doubles (which is what arma::vec holds). This + * distribution internally converts those doubles back into size_t before + * comparisons. + * + * DiscreteDistribution has two template parameters that control the internal + * probability representation type and the observation type. + * + * - `MatType` controls the type used to store probabilities. The element type + * of `MatType` should be a floating-point type. All probabilities returned + * have type equivalent to MatType::elem_type. + * - `ObsMatType` controls the type used to represent observations; by default, + * this is the same as `MatType`. The observations given to Train() or + * Probability() should have type equivalent to `ObsMatType`. The element + * type of `ObsMatType` does not need to be a floating point type. */ +template class DiscreteDistribution { public: + // Convenience typedefs. + typedef typename GetColType::type VecType; + typedef typename MatType::elem_type ElemType; + typedef typename GetColType::type ObsVecType; + typedef typename ObsMatType::elem_type ObsType; + /** * Default constructor, which creates a distribution that has no * observations. */ DiscreteDistribution() : - probabilities(std::vector(1)){ /* Nothing to do. */ } + probabilities(std::vector(1)){ /* Nothing to do. */ } /** * Define the discrete distribution as having numObservations possible @@ -56,8 +75,8 @@ class DiscreteDistribution * can have. */ DiscreteDistribution(const size_t numObservations) : - probabilities(std::vector(1, - arma::ones(numObservations) / numObservations)) + probabilities(std::vector(1, + arma::ones(numObservations) / numObservations)) { /* Nothing to do. */ } /** @@ -80,7 +99,7 @@ class DiscreteDistribution << "must be greater than 0"; throw std::invalid_argument(oss.str()); } - probabilities.push_back(arma::ones(numObs) / numObs); + probabilities.push_back(arma::ones(numObs) / numObs); } } @@ -90,17 +109,17 @@ class DiscreteDistribution * * @param probabilities Probabilities of each possible observation. */ - DiscreteDistribution(const std::vector& probabilities) + DiscreteDistribution(const std::vector& probabilities) { for (size_t i = 0; i < probabilities.size(); ++i) { - arma::vec temp = probabilities[i]; - double sum = accu(temp); + const VecType& temp = probabilities[i]; + ElemType sum = accu(temp); if (sum > 0) this->probabilities.push_back(temp / sum); else { - this->probabilities.push_back(arma::ones(temp.n_elem) + this->probabilities.push_back(arma::ones(temp.n_elem) / temp.n_elem); } } @@ -119,9 +138,9 @@ class DiscreteDistribution * @param observation Observation to return the probability of. * @return Probability of the given observation. */ - double Probability(const arma::vec& observation) const + ElemType Probability(const ObsVecType& observation) const { - double probability = 1.0; + ElemType probability = 1.0; // Ensure the observation has the same dimension with the probabilities. if (observation.n_elem != probabilities.size()) { @@ -134,7 +153,8 @@ class DiscreteDistribution { // Adding 0.5 helps ensure that we cast the floating point to a size_t // correctly. - const size_t obs = size_t(observation(dimension) + 0.5); + const size_t obs = (std::is_floating_point::value) ? + size_t(observation(dimension) + 0.5) : size_t(observation(dimension)); // Ensure that the observation is within the bounds. if (obs >= probabilities[dimension].n_elem) @@ -158,7 +178,7 @@ class DiscreteDistribution * @param observation Observation to return the log probability of. * @return Log probability of the given observation. */ - double LogProbability(const arma::vec& observation) const + ElemType LogProbability(const ObsVecType& observation) const { // TODO: consider storing log probabilities instead? return std::log(Probability(observation)); @@ -171,7 +191,7 @@ class DiscreteDistribution * @param x List of observations. * @param probabilities Output probabilities for each input observation. */ - void Probability(const arma::mat& x, arma::vec& probabilities) const + void Probability(const ObsMatType& x, VecType& probabilities) const { probabilities.set_size(x.n_cols); for (size_t i = 0; i < x.n_cols; ++i) @@ -186,7 +206,7 @@ class DiscreteDistribution * @param logProbabilities Output log-probabilities for each input * observation. */ - void LogProbability(const arma::mat& x, arma::vec& logProbabilities) const + void LogProbability(const ObsMatType& x, VecType& logProbabilities) const { logProbabilities.set_size(x.n_cols); for (size_t i = 0; i < x.n_cols; ++i) @@ -200,7 +220,7 @@ class DiscreteDistribution * * @return Random observation. */ - arma::vec Random() const; + ObsVecType Random() const; /** * Estimate the probability distribution directly from the given @@ -209,7 +229,7 @@ class DiscreteDistribution * * @param observations List of observations. */ - void Train(const arma::mat& observations); + void Train(const ObsMatType& observations); /** * Estimate the probability distribution from the given observations, taking @@ -220,13 +240,13 @@ class DiscreteDistribution * @param probabilities List of probabilities that each observation is * actually from this distribution. */ - void Train(const arma::mat& observations, - const arma::vec& probabilities); + void Train(const ObsMatType& observations, + const VecType& probabilities); //! Return the vector of probabilities for the given dimension. - arma::vec& Probabilities(const size_t dim = 0) { return probabilities[dim]; } + VecType& Probabilities(const size_t dim = 0) { return probabilities[dim]; } //! Modify the vector of probabilities for the given dimension. - const arma::vec& Probabilities(const size_t dim = 0) const + const VecType& Probabilities(const size_t dim = 0) const { return probabilities[dim]; } /** @@ -241,7 +261,7 @@ class DiscreteDistribution private: //! The probabilities for each dimension; each arma::vec represents the //! probabilities for the observations in each dimension. - std::vector probabilities; + std::vector probabilities; }; } // namespace mlpack diff --git a/src/mlpack/core/distributions/discrete_distribution_impl.hpp b/src/mlpack/core/distributions/discrete_distribution_impl.hpp index 4643d61d03..360c382a94 100644 --- a/src/mlpack/core/distributions/discrete_distribution_impl.hpp +++ b/src/mlpack/core/distributions/discrete_distribution_impl.hpp @@ -21,23 +21,25 @@ namespace mlpack { * Return a randomly generated observation according to the probability * distribution defined by this object. */ -inline arma::vec DiscreteDistribution::Random() const +template +inline typename DiscreteDistribution::ObsVecType +DiscreteDistribution::Random() const { size_t dimension = probabilities.size(); - arma::vec result(dimension); + ObsVecType result(dimension); for (size_t d = 0; d < dimension; d++) { // Generate a random number. - double randObs = mlpack::Random(); + ElemType randObs = (ElemType) mlpack::Random(); - double sumProb = 0; + ElemType sumProb = 0; for (size_t obs = 0; obs < probabilities[d].n_elem; obs++) { if ((sumProb += probabilities[d][obs]) >= randObs) { - result[d] = obs; + result[d] = ObsType(obs); break; } } @@ -45,7 +47,7 @@ inline arma::vec DiscreteDistribution::Random() const if (sumProb > 1.0) { // This shouldn't happen. - result[d] = probabilities[d].n_elem - 1; + result[d] = ObsType(probabilities[d].n_elem - 1); } } @@ -55,7 +57,9 @@ inline arma::vec DiscreteDistribution::Random() const /** * Estimate the probability distribution directly from the given observations. */ -inline void DiscreteDistribution::Train(const arma::mat& observations) +template +inline void DiscreteDistribution::Train( + const ObsMatType& observations) { // Make sure the observations have same dimension as the probabilities. if (observations.n_rows != probabilities.size()) @@ -79,7 +83,8 @@ inline void DiscreteDistribution::Train(const arma::mat& observations) // Add the probability of each observation. The addition of 0.5 to the // observation is to turn the default flooring operation of the size_t // cast into a rounding observation. - const size_t obs = size_t(observations(i, r) + 0.5); + const size_t obs = (std::is_floating_point::value) ? + size_t(observations(i, r) + 0.5) : size_t(observations(i, r)); // Ensure that the observation is within the bounds. if (obs >= probabilities[i].n_elem) @@ -97,7 +102,7 @@ inline void DiscreteDistribution::Train(const arma::mat& observations) // Now normalize the distributions. for (size_t i = 0; i < dimensions; ++i) { - double sum = accu(probabilities[i]); + ElemType sum = accu(probabilities[i]); if (sum > 0) probabilities[i] /= sum; else // Force normalization. @@ -109,8 +114,10 @@ inline void DiscreteDistribution::Train(const arma::mat& observations) * Estimate the probability distribution from the given observations when also * given probabilities that each observation is from this distribution. */ -inline void DiscreteDistribution::Train(const arma::mat& observations, - const arma::vec& probObs) +template +inline void DiscreteDistribution::Train( + const ObsMatType& observations, + const VecType& probObs) { // Make sure the observations have same dimension as the probabilities. if (observations.n_rows != probabilities.size()) @@ -134,7 +141,8 @@ inline void DiscreteDistribution::Train(const arma::mat& observations, // Add the probability of each observation. The addition of 0.5 // to the observation is to turn the default flooring operation // of the size_t cast into a rounding observation. - const size_t obs = size_t(observations(i, r) + 0.5); + const size_t obs = (std::is_floating_point::value) ? + size_t(observations(i, r) + 0.5) : size_t(observations(i, r)); // Ensure that the observation is within the bounds. if (obs >= probabilities[i].n_elem) @@ -153,7 +161,7 @@ inline void DiscreteDistribution::Train(const arma::mat& observations, // Now normalize the distributions. for (size_t i = 0; i < dimensions; ++i) { - double sum = accu(probabilities[i]); + ElemType sum = accu(probabilities[i]); if (sum > 0) probabilities[i] /= sum; else // Force normalization. diff --git a/src/mlpack/core/distributions/gamma_distribution.hpp b/src/mlpack/core/distributions/gamma_distribution.hpp index 1b261a25ad..b1b8079a99 100644 --- a/src/mlpack/core/distributions/gamma_distribution.hpp +++ b/src/mlpack/core/distributions/gamma_distribution.hpp @@ -49,9 +49,14 @@ namespace mlpack { * } * @endcode */ +template class GammaDistribution { public: + // Convenience typedefs. + typedef typename GetColType::type VecType; + typedef typename MatType::elem_type ElemType; + /** * Construct the Gamma distribution with the given number of dimensions * (default 0); each parameter will be initialized to 0. @@ -68,7 +73,9 @@ class GammaDistribution * It will stop the approximation once the *change* in the value is * smaller than tol. */ - GammaDistribution(const arma::mat& data, const double tol = 1e-8); + GammaDistribution(const MatType& data, + const ElemType tol = + std::is_same::value ? 1e-4 : 1e-8); /** * Construct the Gamma distribution given two vectors alpha and beta. @@ -76,7 +83,7 @@ class GammaDistribution * @param alpha The vector of alphas, one per dimension. * @param beta The vector of betas, one per dimension. */ - GammaDistribution(const arma::vec& alpha, const arma::vec& beta); + GammaDistribution(const VecType& alpha, const VecType& beta); /** * Destructor. @@ -92,7 +99,9 @@ class GammaDistribution * It will stop the approximation once the *change* in the value is * smaller than tol. */ - void Train(const arma::mat& rdata, const double tol = 1e-8); + void Train(const MatType& rdata, + const ElemType tol = + std::is_same::value ? 1e-4 : 1e-8); /** * Fits an alpha and beta parameter according to observation probabilities. @@ -105,9 +114,10 @@ class GammaDistribution * It will stop the approximation once the *change* in the value is * smaller than tol. */ - void Train(const arma::mat& observations, - const arma::vec& probabilities, - const double tol = 1e-8); + void Train(const MatType& observations, + const VecType& probabilities, + const ElemType tol = + std::is_same::value ? 1e-4 : 1e-8); /** * This function trains (fits distribution parameters) to a dataset with @@ -122,10 +132,11 @@ class GammaDistribution * It will stop the approximation once the *change* in the value is * smaller than tol. */ - void Train(const arma::vec& logMeanxVec, - const arma::vec& meanLogxVec, - const arma::vec& meanxVec, - const double tol = 1e-8); + void Train(const VecType& logMeanxVec, + const VecType& meanLogxVec, + const VecType& meanxVec, + const ElemType tol = + std::is_same::value ? 1e-4 : 1e-8); /** * This function returns the probability of a group of observations. @@ -143,15 +154,15 @@ class GammaDistribution * @param probabilities Column vector of probabilities, one per * observation. */ - void Probability(const arma::mat& observations, - arma::vec& probabilities) const; + void Probability(const MatType& observations, + VecType& probabilities) const; /** * This function returns the probability of the given observation. * * @param x The observation to compute the probability of. */ - double Probability(const arma::vec& x) const; + ElemType Probability(const VecType& x) const; /** * This is a shortcut to the Probability(arma::mat&, arma::vec&) function @@ -161,7 +172,7 @@ class GammaDistribution * @param x The 1-dimensional observation. * @param dim The dimension for which to calculate the probability. */ - double Probability(double x, const size_t dim) const; + ElemType Probability(ElemType x, const size_t dim) const; /** * This function returns the logarithm of the probability of a group of @@ -181,15 +192,15 @@ class GammaDistribution * @param logProbabilities Column vector of log probabilities, one per * observation. */ - void LogProbability(const arma::mat& observations, - arma::vec& logProbabilities) const; + void LogProbability(const MatType& observations, + VecType& logProbabilities) const; /** * This function returns the log-probability of the given observation. * * @param x The observation to compute the log-probability of. */ - double LogProbability(const arma::vec& x) const; + ElemType LogProbability(const VecType& x) const; /** * This function returns the logarithm of the probability of a single @@ -198,33 +209,33 @@ class GammaDistribution * @param x The 1-dimensional observation. * @param dim The dimension for which to calculate the probability. */ - double LogProbability(double x, const size_t dim) const; + ElemType LogProbability(ElemType x, const size_t dim) const; /** * This function returns an observation of this distribution. */ - arma::vec Random() const; + VecType Random() const; // Access to Gamma distribution parameters. //! Get the alpha parameter of the given dimension. - double Alpha(const size_t dim) const { return alpha[dim]; } + ElemType Alpha(const size_t dim) const { return alpha[dim]; } //! Modify the alpha parameter of the given dimension. - double& Alpha(const size_t dim) { return alpha[dim]; } + ElemType& Alpha(const size_t dim) { return alpha[dim]; } //! Get the beta parameter of the given dimension. - double Beta(const size_t dim) const { return beta[dim]; } + ElemType Beta(const size_t dim) const { return beta[dim]; } //! Modify the beta parameter of the given dimension. - double& Beta(const size_t dim) { return beta[dim]; } + ElemType& Beta(const size_t dim) { return beta[dim]; } //! Get the dimensionality of the distribution. size_t Dimensionality() const { return alpha.n_elem; } private: //! Array of fitted alphas. - arma::vec alpha; + VecType alpha; //! Array of fitted betas. - arma::vec beta; + VecType beta; /** * This is a small function that returns true if the update of alpha is @@ -237,9 +248,9 @@ class GammaDistribution * @param tol Convergence tolerance. Relative measure (see documentation of * GammaDistribution::Train). */ - inline bool Converged(const double aOld, - const double aNew, - const double tol); + inline bool Converged(const ElemType aOld, + const ElemType aNew, + const ElemType tol); }; } // namespace mlpack diff --git a/src/mlpack/core/distributions/gamma_distribution_impl.hpp b/src/mlpack/core/distributions/gamma_distribution_impl.hpp index 9cefecd76a..a9a72fe235 100644 --- a/src/mlpack/core/distributions/gamma_distribution_impl.hpp +++ b/src/mlpack/core/distributions/gamma_distribution_impl.hpp @@ -17,21 +17,25 @@ namespace mlpack { -inline GammaDistribution::GammaDistribution(const size_t dimensionality) +template +inline GammaDistribution::GammaDistribution( + const size_t dimensionality) { // Initialize distribution. alpha.zeros(dimensionality); beta.zeros(dimensionality); } -inline GammaDistribution::GammaDistribution(const arma::mat& data, - const double tol) +template +inline GammaDistribution::GammaDistribution(const MatType& data, + const ElemType tol) { Train(data, tol); } -inline GammaDistribution::GammaDistribution(const arma::vec& alpha, - const arma::vec& beta) +template +inline GammaDistribution::GammaDistribution(const VecType& alpha, + const VecType& beta) { if (beta.n_elem != alpha.n_elem) throw std::runtime_error("Alpha and beta vector dimensions mismatch."); @@ -41,24 +45,27 @@ inline GammaDistribution::GammaDistribution(const arma::vec& alpha, } // Returns true if computation converged. -inline bool GammaDistribution::Converged(const double aOld, - const double aNew, - const double tol) +template +inline bool GammaDistribution::Converged(const ElemType aOld, + const ElemType aNew, + const ElemType tol) { return (std::abs(aNew - aOld) / aNew) < tol; } // Fits an alpha and beta parameter to each dimension of the data. -inline void GammaDistribution::Train(const arma::mat& rdata, const double tol) +template +inline void GammaDistribution::Train(const MatType& rdata, + const ElemType tol) { // If fittingSet is empty, nothing to do. if (arma::size(rdata) == arma::size(arma::mat())) return; // Calculate log(mean(x)) and mean(log(x)) of each dataset row. - const arma::vec meanLogxVec = arma::mean(log(rdata), 1); - const arma::vec meanxVec = arma::mean(rdata, 1); - const arma::vec logMeanxVec = log(meanxVec); + const VecType meanLogxVec = arma::mean(log(rdata), 1); + const VecType meanxVec = arma::mean(rdata, 1); + const VecType logMeanxVec = log(meanxVec); // Call the statistics-only GammaDistribution::Train() function to fit the // parameters. That function does all the work so we're done. @@ -66,17 +73,18 @@ inline void GammaDistribution::Train(const arma::mat& rdata, const double tol) } // Fits an alpha and beta parameter according to observation probabilities. -inline void GammaDistribution::Train(const arma::mat& rdata, - const arma::vec& probabilities, - const double tol) +template +inline void GammaDistribution::Train(const MatType& rdata, + const VecType& probabilities, + const ElemType tol) { // If fittingSet is empty, nothing to do. if (arma::size(rdata) == arma::size(arma::mat())) return; - arma::vec meanLogxVec(rdata.n_rows, arma::fill::zeros); - arma::vec meanxVec(rdata.n_rows, arma::fill::zeros); - arma::vec logMeanxVec(rdata.n_rows, arma::fill::zeros); + VecType meanLogxVec(rdata.n_rows, arma::fill::zeros); + VecType meanxVec(rdata.n_rows, arma::fill::zeros); + VecType logMeanxVec(rdata.n_rows, arma::fill::zeros); for (size_t i = 0; i < rdata.n_cols; ++i) { @@ -84,7 +92,7 @@ inline void GammaDistribution::Train(const arma::mat& rdata, meanxVec += probabilities(i) * rdata.col(i); } - double totProbability = accu(probabilities); + ElemType totProbability = accu(probabilities); meanLogxVec /= totProbability; meanxVec /= totProbability; @@ -96,10 +104,11 @@ inline void GammaDistribution::Train(const arma::mat& rdata, } // Fits an alpha and beta parameter to each dimension of the data. -inline void GammaDistribution::Train(const arma::vec& logMeanxVec, - const arma::vec& meanLogxVec, - const arma::vec& meanxVec, - const double tol) +template +inline void GammaDistribution::Train(const VecType& logMeanxVec, + const VecType& meanLogxVec, + const VecType& meanxVec, + const ElemType tol) { using std::log; @@ -119,13 +128,13 @@ inline void GammaDistribution::Train(const arma::vec& logMeanxVec, for (size_t row = 0; row < ndim; ++row) { // Statistics for this row. - const double meanLogx = meanLogxVec(row); - const double meanx = meanxVec(row); - const double logMeanx = logMeanxVec(row); + const ElemType meanLogx = meanLogxVec(row); + const ElemType meanx = meanxVec(row); + const ElemType logMeanx = logMeanxVec(row); // Starting point for Generalized Newton. - double aEst = 0.5 / (logMeanx - meanLogx); - double aOld; + ElemType aEst = 0.5 / (logMeanx - meanLogx); + ElemType aOld; // Newton's method: In each step, make an update to aEst. If value didn't // change much (abs(aNew - aEst) / aEst < tol), then stop. @@ -135,8 +144,8 @@ inline void GammaDistribution::Train(const arma::vec& logMeanxVec, aOld = aEst; // Calculate new value for alpha. - double nominator = meanLogx - logMeanx + std::log(aEst) - Digamma(aEst); - double denominator = std::pow(aEst, 2) * (1 / aEst - Trigamma(aEst)); + ElemType nominator = meanLogx - logMeanx + std::log(aEst) - Digamma(aEst); + ElemType denominator = std::pow(aEst, 2) * (1 / aEst - Trigamma(aEst)); // Protect against division by 0. if (denominator == 0) @@ -149,7 +158,7 @@ inline void GammaDistribution::Train(const arma::vec& logMeanxVec, if (aEst <= 0) throw std::logic_error("GammaDistribution::Train(): estimated invalid " "negative value for parameter alpha!"); - } while (!Converged(aEst, aOld, tol)); + } while (!Converged(aEst, aOld, tol) && !std::isnan(aEst)); alpha(row) = aEst; beta(row) = meanx / aEst; @@ -157,8 +166,10 @@ inline void GammaDistribution::Train(const arma::vec& logMeanxVec, } // Returns the probability of the provided observations. -inline void GammaDistribution::Probability(const arma::mat& observations, - arma::vec& probabilities) const +template +inline void GammaDistribution::Probability( + const MatType& observations, + VecType& probabilities) const { size_t numObs = observations.n_cols; @@ -166,7 +177,7 @@ inline void GammaDistribution::Probability(const arma::mat& observations, probabilities.ones(numObs); // Compute denominator only once for each dimension. - arma::vec denominators(alpha.n_elem); + VecType denominators(alpha.n_elem); for (size_t d = 0; d < alpha.n_elem; ++d) denominators(d) = std::tgamma(alpha(d)) * std::pow(beta(d), alpha(d)); @@ -176,8 +187,8 @@ inline void GammaDistribution::Probability(const arma::mat& observations, for (size_t d = 0; d < observations.n_rows; ++d) { // Compute probability using Multiplication Law. - double factor = std::exp(-observations(d, i) / beta(d)); - double numerator = std::pow(observations(d, i), alpha(d) - 1); + ElemType factor = std::exp(-observations(d, i) / beta(d)); + ElemType numerator = std::pow(observations(d, i), alpha(d) - 1); probabilities(i) *= factor * numerator / denominators(d); } @@ -185,9 +196,11 @@ inline void GammaDistribution::Probability(const arma::mat& observations, } // Return the probability of the given observation. -inline double GammaDistribution::Probability(const arma::vec& x) const +template +inline typename GammaDistribution::ElemType +GammaDistribution::Probability(const VecType& x) const { - double prob = 1.0; + ElemType prob = 1.0; for (size_t d = 0; d < Dimensionality(); ++d) prob *= Probability(x[d], d); @@ -196,16 +209,19 @@ inline double GammaDistribution::Probability(const arma::vec& x) const // Returns the probability of one observation (x) for one of the Gamma's // dimensions. -inline double GammaDistribution::Probability(double x, size_t dim) const +template +inline typename GammaDistribution::ElemType +GammaDistribution::Probability(ElemType x, size_t dim) const { return std::pow(x, alpha(dim) - 1) * std::exp(-x / beta(dim)) / (std::tgamma(alpha(dim)) * std::pow(beta(dim), alpha(dim))); } // Returns the log probability of the provided observations. -inline void GammaDistribution::LogProbability( - const arma::mat& observations, - arma::vec& logProbabilities) const +template +inline void GammaDistribution::LogProbability( + const MatType& observations, + VecType& logProbabilities) const { size_t numObs = observations.n_cols; @@ -213,7 +229,7 @@ inline void GammaDistribution::LogProbability( logProbabilities.zeros(numObs); // Compute denominator only once for each dimension. - arma::vec denominators(alpha.n_elem); + VecType denominators(alpha.n_elem); for (size_t d = 0; d < alpha.n_elem; ++d) denominators(d) = std::tgamma(alpha(d)) * std::pow(beta(d), alpha(d)); @@ -224,8 +240,8 @@ inline void GammaDistribution::LogProbability( { // Compute probability using Multiplication Law and Logarithm addition // property. - double factor = std::exp(-observations(d, i) / beta(d)); - double numerator = std::pow(observations(d, i), alpha(d) - 1); + ElemType factor = std::exp(-observations(d, i) / beta(d)); + ElemType numerator = std::pow(observations(d, i), alpha(d) - 1); logProbabilities(i) += std::log(numerator * factor / denominators(d)); } @@ -233,9 +249,11 @@ inline void GammaDistribution::LogProbability( } // Return the log-probability of the given observation. -inline double GammaDistribution::LogProbability(const arma::vec& x) const +template +inline typename GammaDistribution::ElemType +GammaDistribution::LogProbability(const VecType& x) const { - double logProb = 0.0; + ElemType logProb = 0.0; for (size_t d = 0; d < Dimensionality(); ++d) logProb += LogProbability(x[d], d); return logProb; @@ -243,20 +261,24 @@ inline double GammaDistribution::LogProbability(const arma::vec& x) const // Returns the log probability of one observation (x) for one of the Gamma's // dimensions. -inline double GammaDistribution::LogProbability(double x, size_t dim) const +template +inline typename GammaDistribution::ElemType +GammaDistribution::LogProbability(ElemType x, size_t dim) const { return std::log(std::pow(x, alpha(dim) - 1) * std::exp(-x / beta(dim)) / (std::tgamma(alpha(dim)) * std::pow(beta(dim), alpha(dim)))); } // Returns a gamma-random d-dimensional vector. -inline arma::vec GammaDistribution::Random() const +template +inline typename GammaDistribution::VecType +GammaDistribution::Random() const { - arma::vec randVec(alpha.n_elem); + VecType randVec(alpha.n_elem); for (size_t d = 0; d < alpha.n_elem; ++d) { - std::gamma_distribution dist(alpha(d), beta(d)); + std::gamma_distribution dist(alpha(d), beta(d)); // Use the mlpack random object. randVec(d) = dist(RandGen()); } diff --git a/src/mlpack/core/distributions/laplace_distribution.hpp b/src/mlpack/core/distributions/laplace_distribution.hpp index 96dec858a6..f64b500fc2 100644 --- a/src/mlpack/core/distributions/laplace_distribution.hpp +++ b/src/mlpack/core/distributions/laplace_distribution.hpp @@ -47,9 +47,14 @@ namespace mlpack { * algebra in the paper above becomes simplified, and the PDF takes roughly * the same form as the univariate case. */ +template class LaplaceDistribution { public: + // Convenience typedefs. + typedef typename GetColType::type VecType; + typedef typename MatType::elem_type ElemType; + /** * Default constructor, which creates a Laplace distribution with zero * dimension and zero scale parameter. @@ -64,7 +69,7 @@ class LaplaceDistribution * @param scale Scale of distribution. */ LaplaceDistribution(const size_t dimensionality, const double scale) : - mean(arma::zeros(dimensionality)), scale(scale) { } + mean(arma::zeros(dimensionality)), scale(scale) { } /** * Construct the Laplace distribution with the given mean and scale @@ -73,7 +78,7 @@ class LaplaceDistribution * @param mean Mean of distribution. * @param scale Scale of distribution. */ - LaplaceDistribution(const arma::vec& mean, const double scale) : + LaplaceDistribution(const VecType& mean, const double scale) : mean(mean), scale(scale) { } //! Return the dimensionality of this distribution. @@ -84,7 +89,7 @@ class LaplaceDistribution * * @param observation Point to evaluate probability at. */ - double Probability(const arma::vec& observation) const + ElemType Probability(const VecType& observation) const { return std::exp(LogProbability(observation)); } @@ -95,14 +100,14 @@ class LaplaceDistribution * @param x List of observations. * @param probabilities Output probabilities for each input observation. */ - void Probability(const arma::mat& x, arma::vec& probabilities) const; + void Probability(const MatType& x, VecType& probabilities) const; /** * Return the log probability of the given observation. * * @param observation Point to evaluate logarithm of probability. */ - double LogProbability(const arma::vec& observation) const; + ElemType LogProbability(const VecType& observation) const; /** * Evaluate log probability density function of given observation. @@ -110,7 +115,7 @@ class LaplaceDistribution * @param x List of observations. * @param logProbabilities Output probabilities for each input observation. */ - void LogProbability(const arma::mat& x, arma::vec& logProbabilities) const + void LogProbability(const MatType& x, VecType& logProbabilities) const { logProbabilities.set_size(x.n_cols); for (size_t i = 0; i < x.n_cols; ++i) @@ -125,9 +130,9 @@ class LaplaceDistribution * * @return Random observation from this Laplace distribution. */ - arma::vec Random() const + VecType Random() const { - arma::vec result(mean.n_elem); + VecType result(mean.n_elem); result.randu(); result = mean + scale * log(1.0 + 2.0 * sign(result - 0.5) * (result - 0.5)); @@ -140,7 +145,7 @@ class LaplaceDistribution * @param observations List of observations. */ [[deprecated("Will be removed in mlpack 5.0.0; use Train() instead")]] - void Estimate(const arma::mat& observations); + void Estimate(const MatType& observations); /** * Estimate the Laplace distribution from the given observations, taking into @@ -148,34 +153,34 @@ class LaplaceDistribution * distribution. */ [[deprecated("Will be removed in mlpack 5.0.0; use Train() instead")]] - void Estimate(const arma::mat& observations, - const arma::vec& probabilities); + void Estimate(const MatType& observations, + const VecType& probabilities); /** * Estimate the Laplace distribution directly from the given observations. * * @param observations List of observations. */ - void Train(const arma::mat& observations); + void Train(const MatType& observations); /** * Estimate the Laplace distribution from the given observations, taking into * account the probability of each observation actually being from this * distribution. */ - void Train(const arma::mat& observations, - const arma::vec& probabilities); + void Train(const MatType& observations, + const VecType& probabilities); //! Return the mean. - const arma::vec& Mean() const { return mean; } + const VecType& Mean() const { return mean; } //! Modify the mean. - arma::vec& Mean() { return mean; } + VecType& Mean() { return mean; } //! Return the scale parameter. - double Scale() const { return scale; } + ElemType Scale() const { return scale; } //! Modify the scale parameter. - double& Scale() { return scale; } + ElemType& Scale() { return scale; } /** * Serialize the distribution. @@ -189,9 +194,9 @@ class LaplaceDistribution private: //! Mean of the distribution. - arma::vec mean; + VecType mean; //! Scale parameter of the distribution. - double scale; + ElemType scale; }; } // namespace mlpack diff --git a/src/mlpack/core/distributions/laplace_distribution_impl.hpp b/src/mlpack/core/distributions/laplace_distribution_impl.hpp index 2e41b73935..263bed8dcb 100644 --- a/src/mlpack/core/distributions/laplace_distribution_impl.hpp +++ b/src/mlpack/core/distributions/laplace_distribution_impl.hpp @@ -20,8 +20,10 @@ namespace mlpack { /** * Return the log probability of the given observation. */ -inline double LaplaceDistribution::LogProbability( - const arma::vec& observation) const +template +inline typename LaplaceDistribution::ElemType +LaplaceDistribution::LogProbability( + const VecType& observation) const { // Evaluate the PDF of the Laplace distribution to determine // the log probability. @@ -34,8 +36,10 @@ inline double LaplaceDistribution::LogProbability( * @param x List of observations. * @param probabilities Output probabilities for each input observation. */ -inline void LaplaceDistribution::Probability(const arma::mat& x, - arma::vec& probabilities) const +template +inline void LaplaceDistribution::Probability( + const MatType& x, + VecType& probabilities) const { probabilities.set_size(x.n_cols); for (size_t i = 0; i < x.n_cols; ++i) @@ -49,8 +53,9 @@ inline void LaplaceDistribution::Probability(const arma::mat& x, * * @param observations List of observations. */ +template [[deprecated("Will be removed in mlpack 5.0.0; use Train() instead")]] -inline void LaplaceDistribution::Estimate(const arma::mat& observations) +inline void LaplaceDistribution::Estimate(const MatType& observations) { Train(observations); } @@ -60,9 +65,10 @@ inline void LaplaceDistribution::Estimate(const arma::mat& observations) * taking into account the probability of each observation actually being from * this distribution. */ +template [[deprecated("Will be removed in mlpack 5.0.0; use Train() instead")]] -inline void LaplaceDistribution::Estimate(const arma::mat& observations, - const arma::vec& probabilities) +inline void LaplaceDistribution::Estimate(const MatType& observations, + const VecType& probabilities) { Train(observations, probabilities); } @@ -72,7 +78,8 @@ inline void LaplaceDistribution::Estimate(const arma::mat& observations, * * @param observations List of observations. */ -inline void LaplaceDistribution::Train(const arma::mat& observations) +template +inline void LaplaceDistribution::Train(const MatType& observations) { // The maximum likelihood estimate of the mean is the median of the data for // the univariate case. See the short note "The Double Exponential @@ -105,8 +112,9 @@ inline void LaplaceDistribution::Train(const arma::mat& observations) * taking into account the probability of each observation actually being from * this distribution. */ -inline void LaplaceDistribution::Train(const arma::mat& observations, - const arma::vec& probabilities) +template +inline void LaplaceDistribution::Train(const MatType& observations, + const VecType& probabilities) { // I am not completely sure that this change results in a valid maximum // likelihood estimator given probabilities of points. diff --git a/src/mlpack/core/distributions/regression_distribution.hpp b/src/mlpack/core/distributions/regression_distribution.hpp index 54a7c4d6dd..163fb02743 100644 --- a/src/mlpack/core/distributions/regression_distribution.hpp +++ b/src/mlpack/core/distributions/regression_distribution.hpp @@ -27,13 +27,20 @@ namespace mlpack { * The hmm observations should have the dependent variable in the first row, * with the independent variables in the other rows. */ +template class RegressionDistribution { + public: + // Convenience typedefs. + typedef typename MatType::elem_type ElemType; + typedef typename GetColType::type VecType; + typedef typename GetRowType::type RowType; + private: //! Regression function for representing conditional mean. - LinearRegression<> rf; + LinearRegression rf; //! Error distribution. - GaussianDistribution<> err; + GaussianDistribution err; public: /** @@ -48,12 +55,12 @@ class RegressionDistribution * @param predictors Matrix of predictors (X). * @param responses Vector of responses (y). */ - RegressionDistribution(const arma::mat& predictors, - const arma::rowvec& responses) + RegressionDistribution(const MatType& predictors, + const RowType& responses) { rf.Train(predictors, responses); - err = GaussianDistribution<>(1); - arma::mat cov(1, 1); + err = GaussianDistribution(1); + MatType cov(1, 1); cov(0, 0) = rf.ComputeError(predictors, responses); err.Covariance(std::move(cov)); } @@ -69,21 +76,21 @@ class RegressionDistribution } //! Return regression function. - const LinearRegression<>& Rf() const { return rf; } + const LinearRegression& Rf() const { return rf; } //! Modify regression function. - LinearRegression<>& Rf() { return rf; } + LinearRegression& Rf() { return rf; } //! Return error distribution. - const GaussianDistribution<>& Err() const { return err; } + const GaussianDistribution& Err() const { return err; } //! Modify error distribution. - GaussianDistribution<>& Err() { return err; } + GaussianDistribution& Err() { return err; } /** * Estimate the Gaussian distribution directly from the given observations. * * @param observations List of observations. */ - void Train(const arma::mat& observations); + void Train(const MatType& observations); /** * Estimate parameters using provided observation weights. @@ -91,14 +98,14 @@ class RegressionDistribution * @param observations List of observations. * @param weights Probability that given observation is from distribution. */ - void Train(const arma::mat& observations, const arma::rowvec& weights); + void Train(const MatType& observations, const RowType& weights); /** * Evaluate probability density function of given observation. * * @param observation Point to evaluate probability at. */ - double Probability(const arma::vec& observation) const; + ElemType Probability(const VecType& observation) const; /** * Evaluate probability density function for the given observations. @@ -106,15 +113,15 @@ class RegressionDistribution * @param observations Points to evaluate probability at. * @param probabilities Vector to store computed probabilities in. */ - void Probability(const arma::mat& observations, - arma::vec& probabilities) const; + void Probability(const MatType& observations, + VecType& probabilities) const; /** * Evaluate log probability density function of given observation. * * @param observation Point to evaluate log probability at. */ - double LogProbability(const arma::vec& observation) const + ElemType LogProbability(const VecType& observation) const { return std::log(Probability(observation)); } @@ -124,8 +131,8 @@ class RegressionDistribution * * @param observations Points to evaluate log probability at. */ - void LogProbability(const arma::mat& observations, - arma::vec& probabilities) const + void LogProbability(const MatType& observations, + VecType& probabilities) const { Probability(observations, probabilities); probabilities = arma::log(probabilities); @@ -137,10 +144,10 @@ class RegressionDistribution * @param points The data points to calculate with. * @param predictions Y, will contain calculated values on completion. */ - void Predict(const arma::mat& points, arma::rowvec& predictions) const; + void Predict(const MatType& points, RowType& predictions) const; //! Return the parameters (the b vector). - const arma::vec& Parameters() const { return rf.Parameters(); } + const VecType& Parameters() const { return rf.Parameters(); } //! Return the dimensionality. size_t Dimensionality() const { return rf.Parameters().n_elem; } diff --git a/src/mlpack/core/distributions/regression_distribution_impl.hpp b/src/mlpack/core/distributions/regression_distribution_impl.hpp index ac8034d91b..cd9bdd0927 100644 --- a/src/mlpack/core/distributions/regression_distribution_impl.hpp +++ b/src/mlpack/core/distributions/regression_distribution_impl.hpp @@ -22,23 +22,26 @@ namespace mlpack { * * @param observations List of observations. */ -inline void RegressionDistribution::Train(const arma::mat& observations) +template +inline void RegressionDistribution::Train(const MatType& observations) { - LinearRegression<> lr(observations.rows(1, observations.n_rows - 1), - arma::rowvec(observations.row(0)), 0, true); + LinearRegression lr(observations.rows(1, observations.n_rows - 1), + RowType(observations.row(0)), 0, true); rf = lr; - arma::rowvec fitted; + RowType fitted; lr.Predict(observations.rows(1, observations.n_rows - 1), fitted); err.Train(observations.row(0) - fitted); } -inline void RegressionDistribution::Train(const arma::mat& observations, - const arma::rowvec& weights) +template +inline void RegressionDistribution::Train( + const MatType& observations, + const RowType& weights) { - LinearRegression<> lr(observations.rows(1, observations.n_rows - 1), - arma::rowvec(observations.row(0)), weights, 0, true); + LinearRegression lr(observations.rows(1, observations.n_rows - 1), + RowType(observations.row(0)), weights, 0, true); rf = lr; - arma::rowvec fitted; + RowType fitted; lr.Predict(observations.rows(1, observations.n_rows - 1), fitted); err.Train(observations.row(0) - fitted, weights.t()); } @@ -48,10 +51,11 @@ inline void RegressionDistribution::Train(const arma::mat& observations, * * @param observation Point to evaluate probability at. */ -inline double RegressionDistribution::Probability( - const arma::vec& observation) const +template +inline typename RegressionDistribution::ElemType +RegressionDistribution::Probability(const VecType& observation) const { - arma::rowvec fitted; + RowType fitted; rf.Predict(observation.rows(1, observation.n_rows - 1), fitted); return err.Probability(observation(0) - fitted.t()); } @@ -62,16 +66,20 @@ inline double RegressionDistribution::Probability( * @param observation Points to evaluate probability at. * @param probabilities Vector to store computed probabilities in. */ -inline void RegressionDistribution::Probability(const arma::mat& observations, - arma::vec& probabilities) const +template +inline void RegressionDistribution::Probability( + const MatType& observations, + VecType& probabilities) const { probabilities.set_size(observations.n_cols); for (size_t i = 0; i < observations.n_cols; ++i) probabilities[i] = Probability(observations.unsafe_col(i)); } -inline void RegressionDistribution::Predict(const arma::mat& points, - arma::rowvec& predictions) const +template +inline void RegressionDistribution::Predict( + const MatType& points, + RowType& predictions) const { rf.Predict(points, predictions); } diff --git a/src/mlpack/methods/ann/augmented/tasks/add_impl.hpp b/src/mlpack/methods/ann/augmented/tasks/add_impl.hpp index 6fb7ec1e96..9a767b6dea 100644 --- a/src/mlpack/methods/ann/augmented/tasks/add_impl.hpp +++ b/src/mlpack/methods/ann/augmented/tasks/add_impl.hpp @@ -44,7 +44,7 @@ inline void AddTask::Generate(arma::field& input, arma::vec weights(bitLen - 1); weights = exp2(arma::linspace(1, bitLen - 1, bitLen - 1)); - DiscreteDistribution d(1); + DiscreteDistribution<> d(1); // We have two binary numbers with exactly two digits (10 and 11). // Increasing length by 1 double the number of valid numbers. d.Probabilities(0) = exp2(arma::linspace(1, bitLen - 1, bitLen - 1)); diff --git a/src/mlpack/methods/ann/augmented/tasks/copy_impl.hpp b/src/mlpack/methods/ann/augmented/tasks/copy_impl.hpp index a975710942..77e753ab00 100644 --- a/src/mlpack/methods/ann/augmented/tasks/copy_impl.hpp +++ b/src/mlpack/methods/ann/augmented/tasks/copy_impl.hpp @@ -58,7 +58,7 @@ inline void CopyTask::Generate(arma::field& input, { arma::vec weights(maxLength - 1); - DiscreteDistribution d(1); + DiscreteDistribution<> d(1); // We have two binary numbers with exactly two digits (10 and 11). // Increasing length by 1 double the number of valid numbers. d.Probabilities(0) = diff --git a/src/mlpack/methods/gmm/diagonal_gmm.hpp b/src/mlpack/methods/gmm/diagonal_gmm.hpp index b44de7dfb3..21f885c516 100644 --- a/src/mlpack/methods/gmm/diagonal_gmm.hpp +++ b/src/mlpack/methods/gmm/diagonal_gmm.hpp @@ -41,13 +41,13 @@ namespace mlpack { * @code * void Estimate( * const arma::mat& observations, - * std::vector& dists, + * std::vector>& dists, * arma::vec& weights); * * void Estimate( * const arma::mat& observations, * const arma::vec& probabilities, - * std::vector& dists, + * std::vector>& dists, * arma::vec& weights); * @endcode * @@ -79,7 +79,7 @@ class DiagonalGMM size_t dimensionality; //! Vector of Gaussians. - std::vector dists; + std::vector> dists; //! Vector of a priori weights for each Gaussian. arma::vec weights; @@ -114,7 +114,7 @@ class DiagonalGMM * @param dists Distributions of the model. * @param weights Weights of the model. */ - DiagonalGMM(const std::vector& dists, + DiagonalGMM(const std::vector>& dists, const arma::vec& weights) : gaussians(dists.size()), dimensionality((!dists.empty()) ? dists[0].Mean().n_elem : 0), @@ -137,7 +137,7 @@ class DiagonalGMM * * @param i Index of component. */ - const DiagonalGaussianDistribution& Component(size_t i) const + const DiagonalGaussianDistribution<>& Component(size_t i) const { return dists[i]; } @@ -147,7 +147,7 @@ class DiagonalGMM * * @param i Index of component. */ - DiagonalGaussianDistribution& Component(size_t i) + DiagonalGaussianDistribution<>& Component(size_t i) { return dists[i]; } @@ -239,7 +239,7 @@ class DiagonalGMM * @return The log-likelihood of the best fit. */ template, DiagonalConstraint, - DiagonalGaussianDistribution>> + DiagonalGaussianDistribution<>>> double Train(const arma::mat& observations, const size_t trials = 1, const bool useExistingModel = false, @@ -271,7 +271,7 @@ class DiagonalGMM * @return The log-likelihood of the best fit. */ template, DiagonalConstraint, - DiagonalGaussianDistribution>> + DiagonalGaussianDistribution<>>> double Train(const arma::mat& observations, const arma::vec& probabilities, const size_t trials = 1, @@ -316,7 +316,7 @@ class DiagonalGMM */ double LogLikelihood( const arma::mat& observations, - const std::vector& dists, + const std::vector>& dists, const arma::vec& weights) const; }; diff --git a/src/mlpack/methods/gmm/diagonal_gmm_impl.hpp b/src/mlpack/methods/gmm/diagonal_gmm_impl.hpp index e767649129..4ed02699e8 100644 --- a/src/mlpack/methods/gmm/diagonal_gmm_impl.hpp +++ b/src/mlpack/methods/gmm/diagonal_gmm_impl.hpp @@ -43,7 +43,7 @@ double DiagonalGMM::Train(const arma::mat& observations, // If each trial must start from the same initial location, // we must save it. - std::vector distsOrig; + std::vector> distsOrig; arma::vec weightsOrig; if (useExistingModel) { @@ -61,8 +61,8 @@ double DiagonalGMM::Train(const arma::mat& observations, << bestLikelihood << "." << std::endl; // Now the temporary model. - std::vector distsTrial( - gaussians, DiagonalGaussianDistribution(dimensionality)); + std::vector> distsTrial( + gaussians, DiagonalGaussianDistribution<>(dimensionality)); arma::vec weightsTrial(gaussians); for (size_t trial = 1; trial < trials; ++trial) @@ -113,7 +113,7 @@ inline DiagonalGMM::DiagonalGMM( gaussians(gaussians), dimensionality(dimensionality), dists(gaussians, - DiagonalGaussianDistribution(dimensionality)), + DiagonalGaussianDistribution<>(dimensionality)), weights(gaussians) { // Set equal weights. Technically this model is still valid, but only barely. @@ -287,7 +287,7 @@ inline void DiagonalGMM::Classify(const arma::mat& observations, */ inline double DiagonalGMM::LogLikelihood( const arma::mat& observations, - const std::vector& dists, + const std::vector>& dists, const arma::vec& weights) const { double logLikelihood = 0; @@ -342,7 +342,7 @@ double DiagonalGMM::Train(const arma::mat& observations, return -DBL_MAX; // It's what they asked for... // If each trial must start from the same initial location, we must save it. - std::vector distsOrig; + std::vector> distsOrig; arma::vec weightsOrig; if (useExistingModel) { @@ -361,8 +361,8 @@ double DiagonalGMM::Train(const arma::mat& observations, << bestLikelihood << "." << std::endl; // Now the temporary model. - std::vector distsTrial( gaussians, - DiagonalGaussianDistribution(dimensionality)); + std::vector> distsTrial(gaussians, + DiagonalGaussianDistribution<>(dimensionality)); arma::vec weightsTrial(gaussians); for (size_t trial = 1; trial < trials; ++trial) diff --git a/src/mlpack/methods/gmm/em_fit_impl.hpp b/src/mlpack/methods/gmm/em_fit_impl.hpp index 47af0bffa3..19fc220b7d 100644 --- a/src/mlpack/methods/gmm/em_fit_impl.hpp +++ b/src/mlpack/methods/gmm/em_fit_impl.hpp @@ -44,7 +44,7 @@ Estimate(const arma::mat& observations, arma::vec& weights, const bool useInitialModel) { - if (std::is_same::value) + if (std::is_same>::value) { #ifdef _WIN32 Log::Warn << "Cannot use arma::gmm_diag on Visual Studio due to OpenMP" @@ -129,7 +129,7 @@ Estimate(const arma::mat& observations, // If the distribution is DiagonalGaussianDistribution, calculate the // covariance only with diagonal components. - if (std::is_same::value) + if (std::is_same>::value) { arma::vec covariance = sum((tmp % tmp) % (ones(observations.n_rows) * @@ -240,7 +240,7 @@ Estimate(const arma::mat& observations, // If the distribution is DiagonalGaussianDistribution, calculate the // covariance only with diagonal components. - if (std::is_same::value) + if (std::is_same>::value) { arma::vec cov = sum((tmp % tmp) % (ones(observations.n_rows) * @@ -293,7 +293,7 @@ InitialClustering(const arma::mat& observations, // we can get faster performance by using diagonal elements when calculating // the covariance. const bool isDiagGaussDist = std::is_same::value; + DiagonalGaussianDistribution<>>::value; std::vector means(dists.size()); diff --git a/src/mlpack/methods/gmm/gmm_train_main.cpp b/src/mlpack/methods/gmm/gmm_train_main.cpp index 13c24ad0e6..10cebd70a2 100644 --- a/src/mlpack/methods/gmm/gmm_train_main.cpp +++ b/src/mlpack/methods/gmm/gmm_train_main.cpp @@ -250,7 +250,7 @@ void BINDING_FUNCTION(util::Params& params, util::Timers& timers) // Compute the parameters of the model using the EM algorithm. timers.Start("em"); EMFit em(maxIterations, tolerance, k); + DiagonalGaussianDistribution<>> em(maxIterations, tolerance, k); likelihood = dgmm.Train(dataPoints, params.Get("trials"), false, em); @@ -307,7 +307,7 @@ void BINDING_FUNCTION(util::Params& params, util::Timers& timers) // Compute the parameters of the model using the EM algorithm. timers.Start("em"); EMFit, PositiveDefiniteConstraint, - DiagonalGaussianDistribution> em(maxIterations, tolerance, + DiagonalGaussianDistribution<>> em(maxIterations, tolerance, KMeans<>(kmeansMaxIterations)); likelihood = dgmm.Train(dataPoints, params.Get("trials"), false, diff --git a/src/mlpack/methods/hmm/hmm.hpp b/src/mlpack/methods/hmm/hmm.hpp index a09cc2750b..2b07ce2541 100644 --- a/src/mlpack/methods/hmm/hmm.hpp +++ b/src/mlpack/methods/hmm/hmm.hpp @@ -80,7 +80,7 @@ namespace mlpack { * * @tparam Distribution Type of emission distribution for this HMM. */ -template +template> class HMM { public: diff --git a/src/mlpack/methods/hmm/hmm_model.hpp b/src/mlpack/methods/hmm/hmm_model.hpp index 999067768f..2754904314 100644 --- a/src/mlpack/methods/hmm/hmm_model.hpp +++ b/src/mlpack/methods/hmm/hmm_model.hpp @@ -36,7 +36,7 @@ class HMMModel //! The type of the HMM. HMMType type; //! Not used if type is not DiscreteHMM. - HMM* discreteHMM; + HMM>* discreteHMM; //! Not used if type is not GaussianHMM. HMM>* gaussianHMM; //! Not used if type is not GaussianMixtureModelHMM. @@ -54,7 +54,7 @@ class HMMModel diagGMMHMM(NULL) { if (type == HMMType::DiscreteHMM) - discreteHMM = new HMM(); + discreteHMM = new HMM>(); else if (type == HMMType::GaussianHMM) gaussianHMM = new HMM>(); else if (type == HMMType::GaussianMixtureModelHMM) @@ -73,7 +73,7 @@ class HMMModel { if (type == HMMType::DiscreteHMM) discreteHMM = - new HMM(*other.discreteHMM); + new HMM>(*other.discreteHMM); else if (type == HMMType::GaussianHMM) gaussianHMM = new HMM>(*other.gaussianHMM); @@ -92,7 +92,7 @@ class HMMModel diagGMMHMM(other.diagGMMHMM) { other.type = HMMType::DiscreteHMM; - other.discreteHMM = new HMM(); + other.discreteHMM = new HMM>(); other.gaussianHMM = NULL; other.gmmHMM = NULL; other.diagGMMHMM = NULL; @@ -117,7 +117,7 @@ class HMMModel type = other.type; if (type == HMMType::DiscreteHMM) discreteHMM = - new HMM(*other.discreteHMM); + new HMM>(*other.discreteHMM); else if (type == HMMType::GaussianHMM) gaussianHMM = new HMM>(*other.gaussianHMM); @@ -141,7 +141,7 @@ class HMMModel diagGMMHMM = other.diagGMMHMM; other.type = HMMType::DiscreteHMM; - other.discreteHMM = new HMM(); + other.discreteHMM = new HMM>(); other.gaussianHMM = nullptr; other.gmmHMM = nullptr; other.diagGMMHMM = nullptr; @@ -220,7 +220,7 @@ class HMMModel * gaussianHMM --> NULL * gmmHMM --> NULL * diagGMMHMM --> NULL - * discreteHMM --> HMM object + * discreteHMM --> HMM> object * and hence, calls to GMMHMM(), DiagGMMHMM() and GaussianHMM() will return * NULL. Only the call to DiscreteHMM() will return a non NULL pointer. * @@ -228,7 +228,7 @@ class HMMModel * (by calling the Type() accessor) and then perform subsequent actions, to * avoid null pointer dereferences. */ - HMM* DiscreteHMM() { return discreteHMM; } + HMM>* DiscreteHMM() { return discreteHMM; } HMM>* GaussianHMM() { return gaussianHMM; } HMM* GMMHMM() { return gmmHMM; } HMM* DiagGMMHMM() { return diagGMMHMM; } diff --git a/src/mlpack/methods/hmm/hmm_train_main.cpp b/src/mlpack/methods/hmm/hmm_train_main.cpp index 89f66a1731..e43fde37d2 100644 --- a/src/mlpack/methods/hmm/hmm_train_main.cpp +++ b/src/mlpack/methods/hmm/hmm_train_main.cpp @@ -110,7 +110,7 @@ struct Init //! Helper function to create discrete HMM. static void Create(util::Params& /* params */, - HMM& hmm, + HMM>& hmm, vector& trainSeq, size_t states, double tolerance) @@ -127,8 +127,8 @@ struct Init maxEmissions = arma::max(maxEmissions, maxSeqs); } - hmm = HMM(size_t(states), - DiscreteDistribution(maxEmissions), tolerance); + hmm = HMM>(size_t(states), + DiscreteDistribution<>(maxEmissions), tolerance); } //! Helper function to create Gaussian HMM. @@ -229,7 +229,7 @@ struct Init //! Helper function for discrete emission distributions. static void RandomInitialize(util::Params& /* params */, - vector& e) + vector>& e) { for (size_t i = 0; i < e.size(); ++i) { diff --git a/src/mlpack/methods/hmm/hmm_util_impl.hpp b/src/mlpack/methods/hmm/hmm_util_impl.hpp index d2760a3c27..9ef85b4f99 100644 --- a/src/mlpack/methods/hmm/hmm_util_impl.hpp +++ b/src/mlpack/methods/hmm/hmm_util_impl.hpp @@ -81,7 +81,7 @@ void LoadHMMAndPerformActionHelper(const std::string& modelFile, { case HMMType::DiscreteHMM: DeserializeHMMAndPerformAction>(ar, x); + HMM>>(ar, x); break; case HMMType::GaussianHMM: @@ -161,7 +161,7 @@ template char GetHMMType() { return char(-1); } template<> -char GetHMMType>() +char GetHMMType>>() { return HMMType::DiscreteHMM; } diff --git a/src/mlpack/methods/reinforcement_learning/policy/aggregated_policy.hpp b/src/mlpack/methods/reinforcement_learning/policy/aggregated_policy.hpp index 719fe8091d..4b52ff1a9f 100644 --- a/src/mlpack/methods/reinforcement_learning/policy/aggregated_policy.hpp +++ b/src/mlpack/methods/reinforcement_learning/policy/aggregated_policy.hpp @@ -70,7 +70,7 @@ class AggregatedPolicy std::vector policies; //! Locally-stored sampler under the given distribution. - DiscreteDistribution sampler; + DiscreteDistribution<> sampler; }; } // namespace mlpack diff --git a/src/mlpack/tests/distribution_test.cpp b/src/mlpack/tests/distribution_test.cpp index 08e41eb7dd..e931bd57ee 100644 --- a/src/mlpack/tests/distribution_test.cpp +++ b/src/mlpack/tests/distribution_test.cpp @@ -29,9 +29,19 @@ using namespace mlpack; /** * Make sure we initialize correctly. */ -TEST_CASE("DiscreteDistributionConstructorTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiscreteDistributionConstructorTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - DiscreteDistribution d(5); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; + + DiscreteDistribution d(5); REQUIRE(d.Probabilities().n_elem == 5); REQUIRE(d.Probability("0") == Approx(0.2).epsilon(1e-7)); @@ -44,9 +54,19 @@ TEST_CASE("DiscreteDistributionConstructorTest", "[DistributionTest]") /** * Make sure we get the probabilities of observations right. */ -TEST_CASE("DiscreteDistributionProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiscreteDistributionProbabilityTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - DiscreteDistribution d(5); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; + + DiscreteDistribution d(5); d.Probabilities() = "0.2 0.4 0.1 0.1 0.2"; @@ -60,13 +80,24 @@ TEST_CASE("DiscreteDistributionProbabilityTest", "[DistributionTest]") /** * Make sure we get random observations correct. */ -TEST_CASE("DiscreteDistributionRandomTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiscreteDistributionRandomTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - DiscreteDistribution d(arma::Col("3")); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; + + DiscreteDistribution d(arma::Col("3")); d.Probabilities() = "0.3 0.6 0.1"; - arma::vec actualProb(3); + VecType actualProb(3); actualProb.zeros(); @@ -85,11 +116,21 @@ TEST_CASE("DiscreteDistributionRandomTest", "[DistributionTest]") /** * Make sure we can estimate from observations correctly. */ -TEST_CASE("DiscreteDistributionTrainTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiscreteDistributionTrainTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - DiscreteDistribution d(4); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; - arma::mat obs("0 0 1 1 2 2 2 3"); + DiscreteDistribution d(4); + + ObsMatType obs("0 0 1 1 2 2 2 3"); d.Train(obs); @@ -102,13 +143,23 @@ TEST_CASE("DiscreteDistributionTrainTest", "[DistributionTest]") /** * Estimate from observations with probabilities. */ -TEST_CASE("DiscreteDistributionTrainProbTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiscreteDistributionTrainProbTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - DiscreteDistribution d(3); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; - arma::mat obs("0 0 1 2"); + DiscreteDistribution d(3); - arma::vec prob("0.25 0.25 0.5 1.0"); + ObsMatType obs("0 0 1 2"); + VecType prob("0.25 0.25 0.5 1.0"); d.Train(obs, prob); @@ -120,13 +171,24 @@ TEST_CASE("DiscreteDistributionTrainProbTest", "[DistributionTest]") /** * Achieve multidimensional probability distribution. */ -TEST_CASE("MultiDiscreteDistributionTrainProbTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("MultiDiscreteDistributionTrainProbTest", + "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - DiscreteDistribution d("10 10 10"); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; - arma::mat obs("0 1 1 1 2 2 2 2 2 2;" - "0 0 0 1 1 1 2 2 2 2;" - "0 0 0 1 1 2 2 2 2 2;"); + DiscreteDistribution d("10 10 10"); + + ObsMatType obs("0 1 1 1 2 2 2 2 2 2;" + "0 0 0 1 1 1 2 2 2 2;" + "0 0 0 1 1 2 2 2 2 2;"); d.Train(obs); REQUIRE(d.Probability("0 0 0") == Approx(0.009).epsilon(1e-7)); @@ -138,9 +200,20 @@ TEST_CASE("MultiDiscreteDistributionTrainProbTest", "[DistributionTest]") * Make sure we initialize multidimensional probability distribution * correctly. */ -TEST_CASE("MultiDiscreteDistributionConstructorTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("MultiDiscreteDistributionConstructorTest", + "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - DiscreteDistribution d("4 4 4 4"); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; + + DiscreteDistribution d("4 4 4 4"); REQUIRE(d.Probabilities(0).size() == 4); REQUIRE(d.Dimensionality() == 4); @@ -151,14 +224,25 @@ TEST_CASE("MultiDiscreteDistributionConstructorTest", "[DistributionTest]") /** * Achieve multidimensional probability distribution. */ -TEST_CASE("MultiDiscreteDistributionTrainTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("MultiDiscreteDistributionTrainTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - std::vector pro; - pro.push_back(arma::vec("0.1, 0.3, 0.6")); - pro.push_back(arma::vec("0.3, 0.3, 0.3")); - pro.push_back(arma::vec("0.25, 0.25, 0.5")); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; - DiscreteDistribution d(pro); + std::vector pro; + pro.push_back(VecType("0.1, 0.3, 0.6")); + pro.push_back(VecType("0.3, 0.3, 0.3")); + pro.push_back(VecType("0.25, 0.25, 0.5")); + + DiscreteDistribution d(pro); REQUIRE(d.Probability("0 0 0") == Approx(0.0083333).epsilon(1e-5)); REQUIRE(d.Probability("0 1 2") == Approx(0.0166666).epsilon(1e-5)); @@ -169,15 +253,27 @@ TEST_CASE("MultiDiscreteDistributionTrainTest", "[DistributionTest]") * Estimate multidimensional probability distribution from observations with * probabilities. */ -TEST_CASE("MultiDiscreteDistributionTrainProTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("MultiDiscreteDistributionTrainProTest", + "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { - DiscreteDistribution d("5 5 5"); + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; - arma::mat obs("0 0 1 1 2;" - "0 1 1 2 2;" - "0 1 1 2 2"); + DiscreteDistribution d("5 5 5"); - arma::vec prob("0.25 0.25 0.25 0.25 1"); + ObsMatType obs("0 0 1 1 2;" + "0 1 1 2 2;" + "0 1 1 2 2"); + + VecType prob("0.25 0.25 0.25 0.25 1"); d.Train(obs, prob); @@ -190,15 +286,26 @@ TEST_CASE("MultiDiscreteDistributionTrainProTest", "[DistributionTest]") * Test the LogProbability() function, for multiple points in the multivariate * Discrete case. */ -TEST_CASE("DiscreteLogProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiscreteLogProbabilityTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; + // Same case as before. - DiscreteDistribution d("5 5"); + DiscreteDistribution d("5 5"); - arma::mat obs("0 2;" - "1 2;"); + ObsMatType obs("0 2;" + "1 2;"); - arma::vec logProb; + VecType logProb; d.LogProbability(obs, logProb); @@ -212,15 +319,26 @@ TEST_CASE("DiscreteLogProbabilityTest", "[DistributionTest]") * Test the Probability() function, for multiple points in the multivariate * Discrete case. */ -TEST_CASE("DiscreteProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiscreteProbabilityTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + typedef typename arma::Mat ObsMatType; + // Same case as before. - DiscreteDistribution d("5 5"); + DiscreteDistribution d("5 5"); - arma::mat obs("0 2;" - "1 2;"); + ObsMatType obs("0 2;" + "1 2;"); - arma::vec prob; + VecType prob; d.Probability(obs, prob); @@ -493,7 +611,7 @@ TEMPLATE_TEST_CASE("GaussianDistributionRandomTest", "[DistributionTest]", typedef typename arma::Col VecType; typedef typename arma::Mat MatType; - const ElemType tol = (std::is_same::value) ? 0.2 : 0.1; + const ElemType tol = (std::is_same::value) ? 0.3 : 0.12; VecType mean("1.0 2.25"); MatType cov("0.85 0.60;" @@ -715,25 +833,29 @@ TEMPLATE_TEST_CASE("GaussianDistributionTrainWithTwoDistProbabilitiesTest", * Make sure that using an object to fit one reference set and then asking * to fit another works properly. */ -TEST_CASE("GammaDistributionTrainTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionTrainTest", "[DistributionTest]", float, + double) { + typedef TestType ElemType; + typedef typename arma::Mat MatType; + // Create a gamma distribution random generator. - double alphaReal = 5.3; - double betaReal = 1.5; - std::gamma_distribution dist(alphaReal, betaReal); + ElemType alphaReal = 5.3; + ElemType betaReal = 1.5; + std::gamma_distribution dist(alphaReal, betaReal); // Create a N x d gamma distribution data and fit the results. size_t N = 200; size_t d = 2; - arma::mat rdata(d, N); + MatType rdata(d, N); // Random generation of gamma-like points. for (size_t j = 0; j < d; ++j) for (size_t i = 0; i < N; ++i) rdata(j, i) = dist(RandGen()); - // Create Gamma object and call Train() on reference set. - GammaDistribution gDist; + // Create GammaDistribution object and call Train() on reference set. + GammaDistribution gDist; gDist.Train(rdata); // Training must estimate d pairs of alpha and beta parameters. @@ -743,7 +865,7 @@ TEST_CASE("GammaDistributionTrainTest", "[DistributionTest]") // Create a N' x d' gamma distribution, fit results without new object. size_t N2 = 350; size_t d2 = 4; - arma::mat rdata2(d2, N2); + MatType rdata2(d2, N2); // Random generation of gamma-like points. for (size_t j = 0; j < d2; ++j) @@ -762,73 +884,85 @@ TEST_CASE("GammaDistributionTrainTest", "[DistributionTest]") * This test verifies that the fitting procedure for GammaDistribution works * properly when probabilities for each sample is given. */ -TEST_CASE("GammaDistributionTrainWithProbabilitiesTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionTrainWithProbabilitiesTest", + "[DistributionTest]", float, double) { - double alphaReal = 5.4; - double betaReal = 6.7; + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + const ElemType tol = (std::is_same::value) ? 0.03 : 0.015; + + ElemType alphaReal = 5.4; + ElemType betaReal = 6.7; // Create a gamma distribution random generator. - std::gamma_distribution dist(alphaReal, betaReal); + std::gamma_distribution dist(alphaReal, betaReal); size_t N = 50000; size_t d = 2; - arma::mat rdata(d, N); + MatType rdata(d, N); for (size_t j = 0; j < d; ++j) for (size_t i = 0; i < N; ++i) rdata(j, i) = dist(RandGen()); // Fill the probabilities randomly. - arma::vec probabilities(N, arma::fill::randu); + VecType probabilities(N, arma::fill::randu); // Fit results with probabilities and data. - GammaDistribution gDist; + GammaDistribution gDist; gDist.Train(rdata, probabilities); // Fit results with only data. - GammaDistribution gDist2; + GammaDistribution gDist2; gDist2.Train(rdata); - REQUIRE(gDist2.Alpha(0) == Approx(gDist.Alpha(0)).epsilon(0.015)); - REQUIRE(gDist2.Beta(0) == Approx(gDist.Beta(0)).epsilon(0.015)); + REQUIRE(gDist2.Alpha(0) == Approx(gDist.Alpha(0)).epsilon(tol)); + REQUIRE(gDist2.Beta(0) == Approx(gDist.Beta(0)).epsilon(tol)); - REQUIRE(gDist2.Alpha(1) == Approx(gDist.Alpha(1)).epsilon(0.015)); - REQUIRE(gDist2.Beta(1) == Approx(gDist.Beta(1)).epsilon(0.015)); + REQUIRE(gDist2.Alpha(1) == Approx(gDist.Alpha(1)).epsilon(tol)); + REQUIRE(gDist2.Beta(1) == Approx(gDist.Beta(1)).epsilon(tol)); - REQUIRE(alphaReal == Approx(gDist.Alpha(0)).epsilon(0.03)); - REQUIRE(betaReal == Approx(gDist.Beta(0)).epsilon(0.03)); + REQUIRE(alphaReal == Approx(gDist.Alpha(0)).epsilon(2 * tol)); + REQUIRE(betaReal == Approx(gDist.Beta(0)).epsilon(2 * tol)); - REQUIRE(alphaReal == Approx(gDist.Alpha(1)).epsilon(0.03)); - REQUIRE(betaReal == Approx(gDist.Beta(1)).epsilon(0.03)); + REQUIRE(alphaReal == Approx(gDist.Alpha(1)).epsilon(2 * tol)); + REQUIRE(betaReal == Approx(gDist.Beta(1)).epsilon(2 * tol)); } /** * This test ensures that the same result is obtained when trained with * probabilities all set to 1 and with no probabilities at all. */ -TEST_CASE("GammaDistributionTrainAllProbabilities1Test", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionTrainAllProbabilities1Test", + "[DistributionTest]", float, double) { - double alphaReal = 5.4; - double betaReal = 6.7; + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + ElemType alphaReal = 5.4; + ElemType betaReal = 6.7; // Create a gamma distribution random generator. - std::gamma_distribution dist(alphaReal, betaReal); + std::gamma_distribution dist(alphaReal, betaReal); size_t N = 1000; size_t d = 2; - arma::mat rdata(d, N); + MatType rdata(d, N); for (size_t j = 0; j < d; ++j) for (size_t i = 0; i < N; ++i) rdata(j, i) = dist(RandGen()); // Fit results with only data. - GammaDistribution gDist; + GammaDistribution gDist; gDist.Train(rdata); // Fit results with data and each probability as 1. - GammaDistribution gDist2; - arma::vec allProbabilities1(N, arma::fill::ones); + GammaDistribution gDist2; + VecType allProbabilities1(N, arma::fill::ones); gDist2.Train(rdata, allProbabilities1); REQUIRE(gDist2.Alpha(0) == Approx(gDist.Alpha(0)).epsilon(1e-7)); @@ -845,23 +979,29 @@ TEST_CASE("GammaDistributionTrainAllProbabilities1Test", "[DistributionTest]") * gamma distribution recovered has the same parameters as the second gamma * distribution with high probabilities. */ -TEST_CASE("GammaDistributionTrainTwoDistProbabilities1Test", - "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionTrainTwoDistProbabilities1Test", + "[DistributionTest]", float, double) { - double alphaReal = 5.4; - double betaReal = 6.7; + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; - double alphaReal2 = 1.9; - double betaReal2 = 8.4; + const ElemType tol = (std::is_same::value) ? 0.25 : 0.05; + + ElemType alphaReal = 5.4; + ElemType betaReal = 6.7; + + ElemType alphaReal2 = 1.9; + ElemType betaReal2 = 8.4; // Create two gamma distribution random generators. - std::gamma_distribution dist(alphaReal, betaReal); - std::gamma_distribution dist2(alphaReal2, betaReal2); + std::gamma_distribution dist(alphaReal, betaReal); + std::gamma_distribution dist2(alphaReal2, betaReal2); size_t N = 50000; size_t d = 2; - arma::mat rdata(d, N); - arma::vec probabilities(N); + MatType rdata(d, N); + VecType probabilities(N); // Draw points alternately from the two different distributions. for (size_t j = 0; j < d; ++j) @@ -883,14 +1023,14 @@ TEST_CASE("GammaDistributionTrainTwoDistProbabilities1Test", probabilities(i) = 0.98 + 0.02 * Random(); } - GammaDistribution gDist; + GammaDistribution gDist; gDist.Train(rdata, probabilities); - REQUIRE(alphaReal2 == Approx(gDist.Alpha(0)).epsilon(0.05)); - REQUIRE(betaReal2 == Approx(gDist.Beta(0)).epsilon(0.05)); + REQUIRE(alphaReal2 == Approx(gDist.Alpha(0)).epsilon(tol)); + REQUIRE(betaReal2 == Approx(gDist.Beta(0)).epsilon(tol)); - REQUIRE(alphaReal2 == Approx(gDist.Alpha(1)).epsilon(0.05)); - REQUIRE(betaReal2 == Approx(gDist.Beta(1)).epsilon(0.05)); + REQUIRE(alphaReal2 == Approx(gDist.Alpha(1)).epsilon(tol)); + REQUIRE(betaReal2 == Approx(gDist.Beta(1)).epsilon(tol)); } /** @@ -899,12 +1039,16 @@ TEST_CASE("GammaDistributionTrainTwoDistProbabilities1Test", * with different alpha/beta parameters so we make sure we don't have some weird * bug that always converges to the same number. */ -TEST_CASE("GammaDistributionFittingTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionFittingTest", "[DistributionTest]", float, + double) { + typedef TestType ElemType; + typedef typename arma::Mat MatType; + // Offset from the actual alpha/beta. 10% is quite a relaxed tolerance since // the random points we generate are few (for test speed) and might be fitted // better by a similar distribution. - double errorTolerance = 10; + ElemType errorTolerance = 10; size_t N = 5000; size_t d = 1; // Only 1 dimension is required for this. @@ -912,18 +1056,18 @@ TEST_CASE("GammaDistributionFittingTest", "[DistributionTest]") /** Iteration 1 (first parameter set) **/ // Create a gamma-random generator and data. - double alphaReal = 5.3; - double betaReal = 1.5; - std::gamma_distribution dist(alphaReal, betaReal); + ElemType alphaReal = 5.3; + ElemType betaReal = 1.5; + std::gamma_distribution dist(alphaReal, betaReal); // Random generation of gamma-like points. - arma::mat rdata(d, N); + MatType rdata(d, N); for (size_t j = 0; j < d; ++j) for (size_t i = 0; i < N; ++i) rdata(j, i) = dist(RandGen()); // Create Gamma object and call Train() on reference set. - GammaDistribution gDist; + GammaDistribution gDist; gDist.Train(rdata); // Estimated parameter must be close to real. @@ -933,18 +1077,18 @@ TEST_CASE("GammaDistributionFittingTest", "[DistributionTest]") /** Iteration 2 (different parameter set) **/ // Create a gamma-random generator and data. - double alphaReal2 = 7.2; - double betaReal2 = 0.9; - std::gamma_distribution dist2(alphaReal2, betaReal2); + ElemType alphaReal2 = 7.2; + ElemType betaReal2 = 0.9; + std::gamma_distribution dist2(alphaReal2, betaReal2); // Random generation of gamma-like points. - arma::mat rdata2(d, N); + MatType rdata2(d, N); for (size_t j = 0; j < d; ++j) for (size_t i = 0; i < N; ++i) rdata2(j, i) = dist2(RandGen()); // Create Gamma object and call Train() on reference set. - GammaDistribution gDist2; + GammaDistribution gDist2; gDist2.Train(rdata2); // Estimated parameter must be close to real. @@ -956,12 +1100,16 @@ TEST_CASE("GammaDistributionFittingTest", "[DistributionTest]") * Test that Train() and the constructor that takes data give the same resulting * distribution. */ -TEST_CASE("GammaDistributionTrainConstructorTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionTrainConstructorTest", + "[DistributionTest]", float, double) { - const arma::mat data = arma::randu(10, 500); + typedef TestType ElemType; + typedef typename arma::Mat MatType; - GammaDistribution d1(data); - GammaDistribution d2; + const MatType data = arma::randu(10, 500); + + GammaDistribution d1(data); + GammaDistribution d2; d2.Train(data); for (size_t i = 0; i < 10; ++i) @@ -975,18 +1123,23 @@ TEST_CASE("GammaDistributionTrainConstructorTest", "[DistributionTest]") * Test that Train() with a dataset and Train() with dataset statistics return * the same results. */ -TEST_CASE("GammaDistributionTrainStatisticsTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionTrainStatisticsTest", "[DistributionTest]", + float, double) { - const arma::mat data = arma::randu(1, 500); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + const MatType data = arma::randu(1, 500); // Train object d1 with the data. - GammaDistribution d1(data); + GammaDistribution d1(data); // Train object d2 with the data's statistics. - GammaDistribution d2; - const arma::vec meanLogx = arma::mean(log(data), 1); - const arma::vec meanx = arma::mean(data, 1); - const arma::vec logMeanx = log(meanx); + GammaDistribution d2; + const VecType meanLogx = arma::mean(log(data), 1); + const VecType meanx = arma::mean(data, 1); + const VecType logMeanx = log(meanx); d2.Train(logMeanx, meanLogx, meanx); REQUIRE(d1.Alpha(0) == Approx(d2.Alpha(0)).epsilon(1e-7)); @@ -997,20 +1150,25 @@ TEST_CASE("GammaDistributionTrainStatisticsTest", "[DistributionTest]") * Tests that Random() generates points that can be reasonably well fit by the * distribution that generated them. */ -TEST_CASE("GammaDistributionRandomTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionRandomTest", "[DistributionTest]", float, + double) { - const arma::vec a("2.0 2.5 3.0"), b("0.4 0.6 1.3"); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + const VecType a("2.0 2.5 3.0"), b("0.4 0.6 1.3"); const size_t numPoints = 2000; // Distribution to generate points. - GammaDistribution d1(a, b); - arma::mat data(3, numPoints); // 3-d points. + GammaDistribution d1(a, b); + MatType data(3, numPoints); // 3-d points. for (size_t i = 0; i < numPoints; ++i) data.col(i) = d1.Random(); // Distribution to fit points. - GammaDistribution d2(data); + GammaDistribution d2(data); for (size_t i = 0; i < 3; ++i) { REQUIRE(d2.Alpha(i) == Approx(a(i)).epsilon(0.15)); // Within 15% @@ -1018,20 +1176,25 @@ TEST_CASE("GammaDistributionRandomTest", "[DistributionTest]") } } -TEST_CASE("GammaDistributionProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionProbabilityTest", "[DistributionTest]", + float, double) { + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + // Train two 1-dimensional distributions. - const arma::vec a1("2.0"), b1("0.9"), a2("3.1"), b2("1.4"); - arma::mat x1("2.0"), x2("2.94"); - arma::vec prob1, prob2; + const VecType a1("2.0"), b1("0.9"), a2("3.1"), b2("1.4"); + MatType x1("2.0"), x2("2.94"); + VecType prob1, prob2; // Evaluated at wolfram|alpha - GammaDistribution d1(a1, b1); + GammaDistribution d1(a1, b1); d1.Probability(x1, prob1); REQUIRE(prob1(0) == Approx(0.267575).epsilon(1e-5)); // Evaluated at wolfram|alpha - GammaDistribution d2(a2, b2); + GammaDistribution d2(a2, b2); d2.Probability(x2, prob2); REQUIRE(prob2(0) == Approx(0.189043).epsilon(1e-5)); @@ -1040,34 +1203,39 @@ TEST_CASE("GammaDistributionProbabilityTest", "[DistributionTest]") REQUIRE(prob2(0) == Approx(d2.Probability(2.94, 0)).epsilon(1e-7)); // Combine into one 2-dimensional distribution. - const arma::vec a3("2.0 3.1"), b3("0.9 1.4"); - arma::mat x3(2, 2); + const VecType a3("2.0 3.1"), b3("0.9 1.4"); + MatType x3(2, 2); x3 = { { 2.0, 2.94 }, { 2.0, 2.94 } }; - arma::vec prob3; + VecType prob3; // Expect that the 2-dimensional distribution returns the product of the // 1-dimensional distributions (evaluated at wolfram|alpha). - GammaDistribution d3(a3, b3); + GammaDistribution d3(a3, b3); d3.Probability(x3, prob3); REQUIRE(prob3(0) == Approx(0.04408).epsilon(1e-4)); REQUIRE(prob3(1) == Approx(0.026165).epsilon(1e-4)); } -TEST_CASE("GammaDistributionLogProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("GammaDistributionLogProbabilityTest", "[DistributionTest]", + float, double) { + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + // Train two 1-dimensional distributions. - const arma::vec a1("2.0"), b1("0.9"), a2("3.1"), b2("1.4"); - arma::mat x1("2.0"), x2("2.94"); - arma::vec logprob1, logprob2; + const VecType a1("2.0"), b1("0.9"), a2("3.1"), b2("1.4"); + MatType x1("2.0"), x2("2.94"); + VecType logprob1, logprob2; // Evaluated at wolfram|alpha - GammaDistribution d1(a1, b1); + GammaDistribution d1(a1, b1); d1.LogProbability(x1, logprob1); REQUIRE(logprob1(0) == Approx(std::log(0.267575)).epsilon(1e-5)); // Evaluated at wolfram|alpha - GammaDistribution d2(a2, b2); + GammaDistribution d2(a2, b2); d2.LogProbability(x2, logprob2); REQUIRE(logprob2(0) == Approx(std::log(0.189043)).epsilon(1e-5)); @@ -1076,15 +1244,15 @@ TEST_CASE("GammaDistributionLogProbabilityTest", "[DistributionTest]") REQUIRE(logprob2(0) == Approx(d2.LogProbability(2.94, 0)).epsilon(1e-7)); // Combine into one 2-dimensional distribution. - const arma::vec a3("2.0 3.1"), b3("0.9 1.4"); - arma::mat x3(2, 2); + const VecType a3("2.0 3.1"), b3("0.9 1.4"); + MatType x3(2, 2); x3 = { { 2.0, 2.94 }, { 2.0, 2.94 } }; - arma::vec logprob3; + VecType logprob3; // Expect that the 2-dimensional distribution returns the product of the // 1-dimensional distributions (evaluated at wolfram|alpha). - GammaDistribution d3(a3, b3); + GammaDistribution d3(a3, b3); d3.LogProbability(x3, logprob3); REQUIRE(logprob3(0) == Approx(std::log(0.04408)).epsilon(1e-5)); REQUIRE(logprob3(1) == Approx(std::log(0.026165)).epsilon(1e-5)); @@ -1093,36 +1261,50 @@ TEST_CASE("GammaDistributionLogProbabilityTest", "[DistributionTest]") /** * Discrete Distribution serialization test. */ -TEST_CASE("DiscreteDistributionTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiscreteDistributionTest", "[DistributionTest]", + (std::pair), + (std::pair), + (std::pair), + (std::pair), + (std::pair)) { + typedef typename TestType::first_type ElemType; + typedef typename TestType::second_type ObsElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + typedef typename arma::Col ObsVecType; + typedef typename arma::Mat ObsMatType; + + const ElemType tol = (std::is_same::value) ? 1e-4 : 1e-8; + // I assume that I am properly saving vectors, so, this should be // straightforward. - arma::vec prob; + VecType prob; prob.randu(12); - std::vector probVector = std::vector(1, prob); - DiscreteDistribution t(probVector); + std::vector probVector = std::vector(1, prob); + DiscreteDistribution t(probVector); - DiscreteDistribution xmlT, jsonT, binaryT; + DiscreteDistribution xmlT, jsonT, binaryT; // Load and save with all serializers. SerializeObjectAll(t, xmlT, jsonT, binaryT); for (size_t i = 0; i < 12; ++i) { - arma::vec obs(1); - obs[0] = i; - const double prob = t.Probability(obs); + ObsVecType obs(1); + obs[0] = (ObsElemType) i; + const ElemType prob = t.Probability(obs); if (prob == 0.0) { - REQUIRE(xmlT.Probability(obs) == Approx(0.0).margin(1e-8)); - REQUIRE(jsonT.Probability(obs) == Approx(0.0).margin(1e-8)); - REQUIRE(binaryT.Probability(obs) == Approx(0.0).margin(1e-8)); + REQUIRE(xmlT.Probability(obs) == Approx(0.0).margin(tol)); + REQUIRE(jsonT.Probability(obs) == Approx(0.0).margin(tol)); + REQUIRE(binaryT.Probability(obs) == Approx(0.0).margin(tol)); } else { - REQUIRE(prob == Approx(xmlT.Probability(obs)).epsilon(1e-10)); - REQUIRE(prob == Approx(jsonT.Probability(obs)).epsilon(1e-10)); - REQUIRE(prob == Approx(binaryT.Probability(obs)).epsilon(1e-10)); + REQUIRE(prob == Approx(xmlT.Probability(obs)).epsilon(tol)); + REQUIRE(prob == Approx(jsonT.Probability(obs)).epsilon(tol)); + REQUIRE(prob == Approx(binaryT.Probability(obs)).epsilon(tol)); } } } @@ -1188,13 +1370,18 @@ TEST_CASE("GaussianDistributionTest", "[DistributionTest]") /** * Laplace Distribution serialization test. */ -TEST_CASE("LaplaceDistributionTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("LaplaceDistributionTest", "[DistributionTest]", float, + double) { - arma::vec mean(20); + typedef TestType ElemType; + typedef arma::Col VecType; + typedef arma::Mat MatType; + + VecType mean(20); mean.randu(); - LaplaceDistribution l(mean, 2.5); - LaplaceDistribution xmlL, jsonL, binaryL; + LaplaceDistribution l(mean, 2.5); + LaplaceDistribution xmlL, jsonL, binaryL; SerializeObjectAll(l, xmlL, jsonL, binaryL); @@ -1208,19 +1395,24 @@ TEST_CASE("LaplaceDistributionTest", "[DistributionTest]") /** * Laplace Distribution Probability Test. */ -TEST_CASE("LaplaceDistributionProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("LaplaceDistributionProbabilityTest", "[DistributionTest]", + float, double) { - LaplaceDistribution l(arma::vec("0.0"), 1.0); + typedef TestType ElemType; + typedef arma::Col VecType; + typedef arma::Mat MatType; + + LaplaceDistribution l(VecType("0.0"), 1.0); // Simple case. - REQUIRE(l.Probability(arma::vec("0.0")) == + REQUIRE(l.Probability(VecType("0.0")) == Approx(0.500000000000000).epsilon(1e-7)); - REQUIRE(l.Probability(arma::vec("1.0")) == + REQUIRE(l.Probability(VecType("1.0")) == Approx(0.183939720585721).epsilon(1e-7)); - arma::mat points = "0.0 1.0;"; + MatType points = "0.0 1.0;"; - arma::vec probabilities; + VecType probabilities; l.Probability(points, probabilities); @@ -1233,19 +1425,24 @@ TEST_CASE("LaplaceDistributionProbabilityTest", "[DistributionTest]") /** * Laplace Distribution Log Probability Test. */ -TEST_CASE("LaplaceDistributionLogProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("LaplaceDistributionLogProbabilityTest", + "[DistributionTest]", float, double) { - LaplaceDistribution l(arma::vec("0.0"), 1.0); + typedef TestType ElemType; + typedef arma::Col VecType; + typedef arma::Mat MatType; + + LaplaceDistribution l(VecType("0.0"), 1.0); // Simple case. - REQUIRE(l.LogProbability(arma::vec("0.0")) == + REQUIRE(l.LogProbability(VecType("0.0")) == Approx(-0.693147180559945).epsilon(1e-7)); - REQUIRE(l.LogProbability(arma::vec("1.0")) == + REQUIRE(l.LogProbability(VecType("1.0")) == Approx(-1.693147180559946).epsilon(1e-7)); - arma::mat points = "0.0 1.0;"; + MatType points = "0.0 1.0;"; - arma::vec logProbabilities; + VecType logProbabilities; l.LogProbability(points, logProbabilities); @@ -1269,8 +1466,8 @@ TEST_CASE("RegressionDistributionTest", "[DistributionTest]") arma::rowvec responses; responses.randn(800); - RegressionDistribution rd(data, responses); - RegressionDistribution xmlRd, jsonRd, binaryRd; + RegressionDistribution<> rd(data, responses); + RegressionDistribution<> xmlRd, jsonRd, binaryRd; // Okay, now save it and load it. SerializeObjectAll(rd, xmlRd, jsonRd, binaryRd); @@ -1313,9 +1510,12 @@ TEST_CASE("RegressionDistributionTest", "[DistributionTest]") * Make sure Diagonal Covariance Gaussian distributions are initialized * correctly. */ -TEST_CASE("DiagonalGaussianDistributionEmptyConstructor", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianDistributionEmptyConstructor", + "[DistributionTest]", float, double) { - DiagonalGaussianDistribution d; + typedef TestType ElemType; + + DiagonalGaussianDistribution> d; REQUIRE(d.Mean().n_elem == 0); REQUIRE(d.Covariance().n_elem == 0); @@ -1325,10 +1525,12 @@ TEST_CASE("DiagonalGaussianDistributionEmptyConstructor", "[DistributionTest]") * Make sure Diagonal Covariance Gaussian distributions are initialized to * the correct dimensionality. */ -TEST_CASE("DiagonalGaussianDistributionDimensionalityConstructor", - "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianDistributionDimensionalityConstructor", + "[DistributionTest]", float, double) { - DiagonalGaussianDistribution d(4); + typedef TestType ElemType; + + DiagonalGaussianDistribution> d(4); REQUIRE(d.Mean().n_elem == 4); REQUIRE(d.Covariance().n_elem == 4); @@ -1338,12 +1540,17 @@ TEST_CASE("DiagonalGaussianDistributionDimensionalityConstructor", * Make sure Diagonal Covariance Gaussian distributions are initialized * correctly when we give a mean and covariance. */ -TEST_CASE("DiagonalGaussianDistributionConstructor", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianDistributionConstructor", + "[DistributionTest]", float, double) { - arma::vec mean = arma::randu(3); - arma::vec covariance = arma::randu(3); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; - DiagonalGaussianDistribution d(mean, covariance); + VecType mean = arma::randu(3); + VecType covariance = arma::randu(3); + + DiagonalGaussianDistribution d(mean, covariance); // Make sure the mean and covariance is correct. for (size_t i = 0; i < 3; ++i) @@ -1357,12 +1564,17 @@ TEST_CASE("DiagonalGaussianDistributionConstructor", "[DistributionTest]") * Make sure the probability of observations is correct. * The values were calculated using 'dmvnorm' in R. */ -TEST_CASE("DiagonalGaussianDistributionProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianDistributionProbabilityTest", + "[DistributionTest]", float, double) { - arma::vec mean("2 5 3 4 1"); - arma::vec cov("3 1 5 3 2"); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; - DiagonalGaussianDistribution d(mean, cov); + VecType mean("2 5 3 4 1"); + VecType cov("3 1 5 3 2"); + + DiagonalGaussianDistribution d(mean, cov); // Observations lists randomly selected. REQUIRE(d.LogProbability("3 5 2 7 8") == @@ -1381,79 +1593,97 @@ TEST_CASE("DiagonalGaussianDistributionProbabilityTest", "[DistributionTest]") * Test DiagonalGaussianDistribution::Probability() in the univariate case. * The values were calculated using 'dmvnorm' in R. */ -TEST_CASE("DiagonalGaussianUnivariateProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianUnivariateProbabilityTest", + "[DistributionTest]", float, double) { - DiagonalGaussianDistribution d(arma::vec("0.0"), arma::vec("1.0")); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + const ElemType tol = (std::is_same::value) ? 1e-4 : 1e-7; + + DiagonalGaussianDistribution d(VecType("0.0"), VecType("1.0")); // Mean: 0.0, Covariance: 1.0 - REQUIRE(d.Probability("0.0") == Approx(0.3989422804014327).epsilon(1e-7)); - REQUIRE(d.Probability("1.0") == Approx(0.24197072451914337).epsilon(1e-7)); - REQUIRE(d.Probability("-1.0") == Approx(0.24197072451914337).epsilon(1e-7)); + REQUIRE(d.Probability("0.0") == Approx(0.3989422804014327).epsilon(tol)); + REQUIRE(d.Probability("1.0") == Approx(0.24197072451914337).epsilon(tol)); + REQUIRE(d.Probability("-1.0") == Approx(0.24197072451914337).epsilon(tol)); // Mean: 0.0, Covariance: 2.0 d.Covariance("2.0"); - REQUIRE(d.Probability("0.0") == Approx(0.28209479177387814).epsilon(1e-7)); - REQUIRE(d.Probability("1.0") == Approx(0.21969564473386122).epsilon(1e-7)); - REQUIRE(d.Probability("-1.0") == Approx(0.21969564473386122).epsilon(1e-7)); + REQUIRE(d.Probability("0.0") == Approx(0.28209479177387814).epsilon(tol)); + REQUIRE(d.Probability("1.0") == Approx(0.21969564473386122).epsilon(tol)); + REQUIRE(d.Probability("-1.0") == Approx(0.21969564473386122).epsilon(tol)); // Mean: 1.0, Covariance: 1.0 d.Mean() = "1.0"; d.Covariance("1.0"); - REQUIRE(d.Probability("0.0") == Approx(0.24197072451914337).epsilon(1e-7)); - REQUIRE(d.Probability("1.0") == Approx(0.3989422804014327).epsilon(1e-7)); - REQUIRE(d.Probability("-1.0") == Approx(0.053990966513188056).epsilon(1e-7)); + REQUIRE(d.Probability("0.0") == Approx(0.24197072451914337).epsilon(tol)); + REQUIRE(d.Probability("1.0") == Approx(0.3989422804014327).epsilon(tol)); + REQUIRE(d.Probability("-1.0") == Approx(0.053990966513188056).epsilon(tol)); // Mean: 1.0, Covariance: 2.0 d.Covariance("2.0"); - REQUIRE(d.Probability("0.0") == Approx(0.21969564473386122).epsilon(1e-7)); - REQUIRE(d.Probability("1.0") == Approx(0.28209479177387814).epsilon(1e-7)); - REQUIRE(d.Probability("-1.0") == Approx(0.10377687435514872).epsilon(1e-7)); + REQUIRE(d.Probability("0.0") == Approx(0.21969564473386122).epsilon(tol)); + REQUIRE(d.Probability("1.0") == Approx(0.28209479177387814).epsilon(tol)); + REQUIRE(d.Probability("-1.0") == Approx(0.10377687435514872).epsilon(tol)); } /** * Test DiagonalGaussianDistribution::Probability() in the multivariate case. * The values were calculated using 'dmvnorm' in R. */ -TEST_CASE("DiagonalGaussianMultivariateProbabilityTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianMultivariateProbabilityTest", + "[DistributionTest]", float, double) { - arma::vec mean("0 0"); - arma::vec cov("2 2"); - arma::vec obs("0 0"); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; - DiagonalGaussianDistribution d(mean, cov); + const ElemType tol = (std::is_same::value) ? 1e-4 : 1e-7; - REQUIRE(d.Probability(obs) == Approx(0.079577471545947673).epsilon(1e-7)); + VecType mean("0 0"); + VecType cov("2 2"); + VecType obs("0 0"); + + DiagonalGaussianDistribution d(mean, cov); + + REQUIRE(d.Probability(obs) == Approx(0.079577471545947673).epsilon(tol)); obs = "1 1"; - REQUIRE(d.Probability(obs) == Approx(0.048266176315026957).epsilon(1e-7)); + REQUIRE(d.Probability(obs) == Approx(0.048266176315026957).epsilon(tol)); d.Mean() = "1 3"; - REQUIRE(d.Probability(obs) == Approx(0.029274915762159581).epsilon(1e-7)); - REQUIRE(d.Probability(-obs) == Approx(0.00053618878559782773).epsilon(1e-7)); + REQUIRE(d.Probability(obs) == Approx(0.029274915762159581).epsilon(tol)); + REQUIRE(d.Probability(-obs) == Approx(0.00053618878559782773).epsilon(tol)); // Higher dimensional case. d.Mean() = "1 3 6 2 7"; d.Covariance("3 1 5 3 2"); obs = "2 5 7 3 8"; - REQUIRE(d.Probability(obs) == Approx(7.2790083003378082e-05).epsilon(1e-7)); + REQUIRE(d.Probability(obs) == Approx(7.2790083003378082e-05).epsilon(tol)); } /** * Test the phi() function, for multiple points in the multivariate Gaussian * case. The values were calculated using 'dmvnorm' in R. */ -TEST_CASE("DiagonalGaussianMultipointMultivariateProbabilityTest", - "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianMultipointMultivariateProbabilityTest", + "[DistributionTest]", float, double) { - arma::vec mean = "2 5 3 7 2"; - arma::vec cov("9 2 1 4 8"); - arma::mat points = "3 5 2 7 5 8;" - "2 6 8 3 4 6;" - "1 4 2 7 8 2;" - "6 8 4 7 9 2;" - "4 6 7 7 3 2"; - arma::vec phis; - DiagonalGaussianDistribution d(mean, cov); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + VecType mean = "2 5 3 7 2"; + VecType cov("9 2 1 4 8"); + MatType points = "3 5 2 7 5 8;" + "2 6 8 3 4 6;" + "1 4 2 7 8 2;" + "6 8 4 7 9 2;" + "4 6 7 7 3 2"; + VecType phis; + DiagonalGaussianDistribution d(mean, cov); d.LogProbability(points, phis); REQUIRE(phis.n_elem == 6); @@ -1469,49 +1699,62 @@ TEST_CASE("DiagonalGaussianMultipointMultivariateProbabilityTest", /** * Make sure random observations follow the probability distribution correctly. */ -TEST_CASE("DiagonalGaussianDistributionRandomTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianDistributionRandomTest", + "[DistributionTest]", float, double) { - arma::vec mean("2.5 1.25"); - arma::vec cov("0.50 0.25"); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; - DiagonalGaussianDistribution d(mean, cov); + const ElemType tol = (std::is_same::value) ? 0.2 : 0.1; - arma::mat obs(2, 5000); + VecType mean("2.5 1.25"); + VecType cov("0.50 0.25"); + DiagonalGaussianDistribution d(mean, cov); + + MatType obs(2, 5000); for (size_t i = 0; i < 5000; ++i) obs.col(i) = d.Random(); // Make sure that reflects the actual distribution. - arma::vec obsMean = arma::mean(obs, 1); - arma::mat obsCov = ColumnCovariance(obs); + VecType obsMean = arma::mean(obs, 1); + MatType obsCov = ColumnCovariance(obs); - // 10% tolerance because this can be noisy. - REQUIRE(obsMean(0) == Approx(mean(0)).epsilon(0.1)); - REQUIRE(obsMean(1) == Approx(mean(1)).epsilon(0.1)); + // 10% tolerance because this can be noisy. (20% for floats.) + REQUIRE(obsMean(0) == Approx(mean(0)).epsilon(tol)); + REQUIRE(obsMean(1) == Approx(mean(1)).epsilon(tol)); - REQUIRE(obsCov(0, 0) == Approx(cov(0)).epsilon(0.1)); - REQUIRE(obsCov(1, 1) == Approx(cov(1)).epsilon(0.1)); + REQUIRE(obsCov(0, 0) == Approx(cov(0)).epsilon(tol)); + REQUIRE(obsCov(1, 1) == Approx(cov(1)).epsilon(tol)); } /** * Make sure that we can properly estimate from given observations. */ -TEST_CASE("DiagonalGaussianDistributionTrainTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianDistributionTrainTest", + "[DistributionTest]", float, double) { - arma::vec mean("2.5 1.5 8.2 3.1"); - arma::vec cov("1.2 3.1 8.3 4.3"); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + const ElemType tol = (std::is_same::value) ? 1e-3 : 1e-5; + + VecType mean("2.5 1.5 8.2 3.1"); + VecType cov("1.2 3.1 8.3 4.3"); // Generate the observations. - arma::mat observations(4, 10000); + MatType observations(4, 10000); for (size_t i = 0; i < 10000; ++i) - observations.col(i) = (sqrt(cov) % arma::randn(4)) + mean; + observations.col(i) = (sqrt(cov) % arma::randn(4)) + mean; - DiagonalGaussianDistribution d; + DiagonalGaussianDistribution d; // Calculate the actual mean and covariance of data using armadillo. - arma::vec actualMean = arma::mean(observations, 1); - arma::mat actualCov = ColumnCovariance(observations); + VecType actualMean = arma::mean(observations, 1); + MatType actualCov = ColumnCovariance(observations); // Estimate the parameters. d.Train(observations); @@ -1519,8 +1762,8 @@ TEST_CASE("DiagonalGaussianDistributionTrainTest", "[DistributionTest]") // Check that the estimated parameters are right. for (size_t i = 0; i < 4; ++i) { - REQUIRE(d.Mean()(i) - actualMean(i) == Approx(0.0).margin(1e-5)); - REQUIRE(d.Covariance()(i) - actualCov(i, i) == Approx(0.0).margin(1e-5)); + REQUIRE(d.Mean()(i) - actualMean(i) == Approx(0.0).margin(tol)); + REQUIRE(d.Covariance()(i) - actualCov(i, i) == Approx(0.0).margin(tol)); } } @@ -1528,30 +1771,37 @@ TEST_CASE("DiagonalGaussianDistributionTrainTest", "[DistributionTest]") * Make sure the unbiased estimator of the weighted sample works correctly. * The values were calculated using 'cov.wt' in R. */ -TEST_CASE("DiagonalGaussianUnbiasedEstimatorTest", "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianUnbiasedEstimatorTest", + "[DistributionTest]", float, double) { + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + const ElemType tol = (std::is_same::value) ? 1e-4 : 1e-7; + // Generate the observations. - arma::mat observations("3 5 2 7;" - "2 6 8 3;" - "1 4 2 7;" - "6 8 4 7"); + MatType observations("3 5 2 7;" + "2 6 8 3;" + "1 4 2 7;" + "6 8 4 7"); - arma::vec probs("0.3 0.4 0.1 0.2"); + VecType probs("0.3 0.4 0.1 0.2"); - DiagonalGaussianDistribution d; + DiagonalGaussianDistribution d; // Estimate the parameters. d.Train(observations, probs); - REQUIRE(d.Mean()(0) == Approx(4.5).epsilon(1e-7)); - REQUIRE(d.Mean()(1) == Approx(4.4).epsilon(1e-7)); - REQUIRE(d.Mean()(2) == Approx(3.5).epsilon(1e-7)); - REQUIRE(d.Mean()(3) == Approx(6.8).epsilon(1e-7)); + REQUIRE(d.Mean()(0) == Approx(4.5).epsilon(tol)); + REQUIRE(d.Mean()(1) == Approx(4.4).epsilon(tol)); + REQUIRE(d.Mean()(2) == Approx(3.5).epsilon(tol)); + REQUIRE(d.Mean()(3) == Approx(6.8).epsilon(tol)); - REQUIRE(d.Covariance()(0) == Approx(3.78571428571428603).epsilon(1e-7)); - REQUIRE(d.Covariance()(1) == Approx(6.34285714285714253).epsilon(1e-7)); - REQUIRE(d.Covariance()(2) == Approx(6.64285714285714235).epsilon(1e-7)); - REQUIRE(d.Covariance()(3) == Approx(2.22857142857142865).epsilon(1e-7)); + REQUIRE(d.Covariance()(0) == Approx(3.78571428571428603).epsilon(tol)); + REQUIRE(d.Covariance()(1) == Approx(6.34285714285714253).epsilon(tol)); + REQUIRE(d.Covariance()(2) == Approx(6.64285714285714235).epsilon(tol)); + REQUIRE(d.Covariance()(3) == Approx(2.22857142857142865).epsilon(tol)); } /** @@ -1559,21 +1809,27 @@ TEST_CASE("DiagonalGaussianUnbiasedEstimatorTest", "[DistributionTest]") * the weighted mean and covariance reduce to the unweighted sample mean and * covariance. */ -TEST_CASE("DiagonalGaussianWeightedParametersReductionTest", - "[DistributionTest]") +TEMPLATE_TEST_CASE("DiagonalGaussianWeightedParametersReductionTest", + "[DistributionTest]", float, double) { - arma::vec mean("2.5 1.5 8.2 3.1"); - arma::vec cov("1.2 3.1 8.3 4.3"); + typedef TestType ElemType; + typedef typename arma::Col VecType; + typedef typename arma::Mat MatType; + + const ElemType tol = (std::is_same::value) ? 1e-4 : 1e-7; + + VecType mean("2.5 1.5 8.2 3.1"); + VecType cov("1.2 3.1 8.3 4.3"); // Generate the observations. - arma::mat obs(4, 5); - arma::vec probs("0.2 0.2 0.2 0.2 0.2"); + MatType obs(4, 5); + VecType probs("0.2 0.2 0.2 0.2 0.2"); for (size_t i = 0; i < 5; ++i) - obs.col(i) = (sqrt(cov) % arma::randn(4)) + mean; + obs.col(i) = (sqrt(cov) % arma::randn(4)) + mean; - DiagonalGaussianDistribution d1; - DiagonalGaussianDistribution d2; + DiagonalGaussianDistribution d1; + DiagonalGaussianDistribution d2; // Estimate the parameters. d1.Train(obs); @@ -1582,7 +1838,7 @@ TEST_CASE("DiagonalGaussianWeightedParametersReductionTest", // Check if these are equal. for (size_t i = 0; i < 4; ++i) { - REQUIRE(d1.Mean()(i) == Approx(d2.Mean()(i)).epsilon(1e-7)); - REQUIRE(d1.Covariance()(i) == Approx(d2.Covariance()(i)).epsilon(1e-7)); + REQUIRE(d1.Mean()(i) == Approx(d2.Mean()(i)).epsilon(tol)); + REQUIRE(d1.Covariance()(i) == Approx(d2.Covariance()(i)).epsilon(tol)); } } diff --git a/src/mlpack/tests/gmm_test.cpp b/src/mlpack/tests/gmm_test.cpp index d790ce4b1a..38f2225313 100644 --- a/src/mlpack/tests/gmm_test.cpp +++ b/src/mlpack/tests/gmm_test.cpp @@ -851,8 +851,8 @@ TEST_CASE("DiagonalGMMProbabilityComponentTest", "[GMMTest]") { // Create DiagonalGMM. DiagonalGMM gmm(2, 2); - gmm.Component(0) = DiagonalGaussianDistribution("0 0", "1 1"); - gmm.Component(1) = DiagonalGaussianDistribution("2 3", "3 2"); + gmm.Component(0) = DiagonalGaussianDistribution<>("0 0", "1 1"); + gmm.Component(1) = DiagonalGaussianDistribution<>("2 3", "3 2"); gmm.Weights() = "0.2 0.8"; // The values are calculated using mlpack's GMM class. @@ -923,7 +923,7 @@ TEST_CASE("DiagonalGMMTrainEMOneGaussian", "[GMMTest]") TEST_CASE("DiagonalGMMTrainEMOneGaussianWithProbability", "[GMMTest]") { // Generate a diagonal covariance gaussian distribution. - DiagonalGaussianDistribution d("1.0 0.8", "1.0 2.0"); + DiagonalGaussianDistribution<> d("1.0 0.8", "1.0 2.0"); // Generate 20000 observations, each with random probabilities. arma::mat observations(2, 20000); @@ -961,9 +961,9 @@ TEST_CASE("DiagonalGMMTrainEMMultipleGaussians", "[GMMTest]") { // We'll have three diagonal covariance Gaussian distributions from this // mixture. - DiagonalGaussianDistribution d1("0.0 1.0 0.0", "1.0 0.8 1.0;"); - DiagonalGaussianDistribution d2("2.0 -1.0 5.0", "3.0 1.2 1.3;"); - DiagonalGaussianDistribution d3("0.0 5.0 -3.0", "2.0 0.3 1.0;"); + DiagonalGaussianDistribution<> d1("0.0 1.0 0.0", "1.0 0.8 1.0;"); + DiagonalGaussianDistribution<> d2("2.0 -1.0 5.0", "3.0 1.2 1.3;"); + DiagonalGaussianDistribution<> d3("0.0 5.0 -3.0", "2.0 0.3 1.0;"); // Now we'll generate points and probabilities. arma::mat observations(3, 5000); @@ -1044,9 +1044,9 @@ TEST_CASE("DiagonalGMMTrainEMMultipleGaussiansWithProbability", "[GMMTest]") { // We'll have three diagonal covariance Gaussian distributions from this // mixture. - DiagonalGaussianDistribution d1("1.5 0.8 1.0", "1.0 0.8 1.0;"); - DiagonalGaussianDistribution d2("8.2 6.3 7.4", "1.0 1.2 1.3;"); - DiagonalGaussianDistribution d3("-4.5 -5.0 -3.0", "2.0 2.3 1.0;"); + DiagonalGaussianDistribution<> d1("1.5 0.8 1.0", "1.0 0.8 1.0;"); + DiagonalGaussianDistribution<> d2("8.2 6.3 7.4", "1.0 1.2 1.3;"); + DiagonalGaussianDistribution<> d3("-4.5 -5.0 -3.0", "2.0 2.3 1.0;"); // Now we'll generate observations and probabilities. arma::mat observations(3, 10000); @@ -1132,8 +1132,8 @@ TEST_CASE("DiagonalGMMRandomTest", "[GMMTest]") DiagonalGMM gmm(2, 2); gmm.Weights() = arma::vec("0.40 0.60"); - gmm.Component(0) = DiagonalGaussianDistribution("1.05 2.60", "0.95 1.01"); - gmm.Component(1) = DiagonalGaussianDistribution("4.30 1.00", "1.05 0.97"); + gmm.Component(0) = DiagonalGaussianDistribution<>("1.05 2.60", "0.95 1.01"); + gmm.Component(1) = DiagonalGaussianDistribution<>("4.30 1.00", "1.05 0.97"); // Now generate a bunch of observations. arma::mat observations(2, 4000); diff --git a/src/mlpack/tests/hmm_test.cpp b/src/mlpack/tests/hmm_test.cpp index 3518fe6a29..93b68848b1 100644 --- a/src/mlpack/tests/hmm_test.cpp +++ b/src/mlpack/tests/hmm_test.cpp @@ -35,11 +35,11 @@ TEST_CASE("SimpleDiscreteHMMTestViterbi", "[HMMTest]") // [0.1 0.8]] no umbrella arma::vec initial("1 0"); // Default MATLAB initial states. arma::mat transition("0.7 0.3; 0.3 0.7"); - std::vector emission(2); - emission[0] = DiscreteDistribution(std::vector{"0.9 0.1"}); - emission[1] = DiscreteDistribution(std::vector{"0.2 0.8"}); + std::vector> emission(2); + emission[0] = DiscreteDistribution<>(std::vector{"0.9 0.1"}); + emission[1] = DiscreteDistribution<>(std::vector{"0.2 0.8"}); - HMM hmm(initial, transition, emission); + HMM> hmm(initial, transition, emission); // Now let's take a sequence and find what the most likely state is. // We'll use the sequence [U U N U U] (U = umbrella, N = no umbrella) like on @@ -72,15 +72,15 @@ TEST_CASE("BorodovskyHMMTestViterbi", "[HMMTest]") "0.5 0.5 0.4;" "0.5 0.5 0.6"); // Four emission states: A, C, G, T. Start state doesn't emit... - std::vector emission(3); - emission[0] = DiscreteDistribution( + std::vector> emission(3); + emission[0] = DiscreteDistribution<>( std::vector{"0.25 0.25 0.25 0.25"}); - emission[1] = DiscreteDistribution( + emission[1] = DiscreteDistribution<>( std::vector{"0.20 0.30 0.30 0.20"}); - emission[2] = DiscreteDistribution( + emission[2] = DiscreteDistribution<>( std::vector{"0.30 0.20 0.20 0.30"}); - HMM hmm(initial, transition, emission); + HMM> hmm(initial, transition, emission); // GGCACTGAA. arma::mat observation("2 2 1 0 1 3 2 0 0"); @@ -115,11 +115,13 @@ TEST_CASE("ForwardBackwardTwoState", "[HMMTest]") // I am not certain. arma::vec initial("0.1 0.4"); arma::mat transition("0.1 0.9; 0.4 0.6"); - std::vector emis(2); - emis[0] = DiscreteDistribution(std::vector{"0.85 0.15 0.00 0.00"}); - emis[1] = DiscreteDistribution(std::vector{"0.00 0.00 0.50 0.50"}); + std::vector> emis(2); + emis[0] = DiscreteDistribution<>( + std::vector{"0.85 0.15 0.00 0.00"}); + emis[1] = DiscreteDistribution<>( + std::vector{"0.00 0.00 0.50 0.50"}); - HMM hmm(initial, transition, emis); + HMM> hmm(initial, transition, emis); // Now check we are getting the same results as MATLAB for this sequence. arma::mat stateProb; @@ -162,7 +164,7 @@ TEST_CASE("ForwardBackwardTwoState", "[HMMTest]") TEST_CASE("SimplestBaumWelchDiscreteHMM", "[HMMTest]") { // Don't yet require a useful distribution. 1 state, 1 emission. - HMM hmm(1, DiscreteDistribution(1)); + HMM> hmm(1, DiscreteDistribution<>(1)); std::vector observations; // Different lengths for each observation sequence. @@ -183,7 +185,7 @@ TEST_CASE("SimplestBaumWelchDiscreteHMM", "[HMMTest]") */ TEST_CASE("SimpleBaumWelchDiscreteHMM", "[HMMTest]") { - HMM hmm(1, 2); // 1 state, 2 emissions. + HMM> hmm(1, 2); // 1 state, 2 emissions. // Randomize the emission matrix. hmm.Emission()[0].Probabilities() = arma::randu(2); hmm.Emission()[0].Probabilities() /= accu(hmm.Emission()[0].Probabilities()); @@ -218,7 +220,7 @@ TEST_CASE("SimpleBaumWelchDiscreteHMM", "[HMMTest]") */ TEST_CASE("SimpleBaumWelchDiscreteHMM_2", "[HMMTest]") { - HMM hmm(2, DiscreteDistribution(4)); + HMM> hmm(2, DiscreteDistribution<>(4)); // A little bit of obfuscation to the solution. hmm.Transition() = arma::mat("0.1 0.4; 0.9 0.6"); @@ -317,7 +319,7 @@ TEST_CASE("DiscreteHMMLabeledTrainTest", "[HMMTest]") { // Generate a random Markov model with 3 hidden states and 6 observations. arma::mat transition; - std::vector emission(3); + std::vector> emission(3); transition.randu(3, 3); emission[0].Probabilities() = arma::randu(6); @@ -374,7 +376,7 @@ TEST_CASE("DiscreteHMMLabeledTrainTest", "[HMMTest]") // Now that our data is generated, we give the HMM the labeled data to train // on. - HMM hmm(3, DiscreteDistribution(6)); + HMM> hmm(3, DiscreteDistribution<>(6)); hmm.Train(observations, states); @@ -407,7 +409,7 @@ TEST_CASE("DiscreteHMMSimpleGenerateTest", "[HMMTest]") { // Very simple HMM. 4 emissions with equal probability and 2 states with // equal probability. - HMM hmm(2, DiscreteDistribution(4)); + HMM> hmm(2, DiscreteDistribution<>(4)); hmm.Initial() = arma::ones(2) / 2.0; hmm.Transition() = arma::ones(2, 2) / 2.0; @@ -450,7 +452,7 @@ TEST_CASE("DiscreteHMMGenerateTest", "[HMMTest]") // 6 emissions, 4 states. Random transition and emission probability. arma::vec initial("1 0 0 0"); arma::mat transition(4, 4); - std::vector emission(4); + std::vector> emission(4); emission[0].Probabilities() = arma::randu(6); emission[0].Probabilities() /= accu(emission[0].Probabilities()); emission[1].Probabilities() = arma::randu(6); @@ -467,7 +469,7 @@ TEST_CASE("DiscreteHMMGenerateTest", "[HMMTest]") transition.col(col) /= accu(transition.col(col)); // Create HMM object. - HMM hmm(initial, transition, emission); + HMM> hmm(initial, transition, emission); // We'll create a bunch of sequences. int numSeq = 400; @@ -483,7 +485,7 @@ TEST_CASE("DiscreteHMMGenerateTest", "[HMMTest]") } // Now we will calculate the full probabilities. - HMM hmm2(4, 6); + HMM> hmm2(4, 6); hmm2.Train(sequences, states); // Check that training gives the same result. @@ -508,12 +510,12 @@ TEST_CASE("DiscreteHMMLogLikelihoodTest", "[HMMTest]") arma::mat transition("0.5 0.0 0.1;" "0.2 0.6 0.2;" "0.3 0.4 0.7"); - std::vector emission(3); + std::vector> emission(3); emission[0].Probabilities() = "0.75 0.25 0.00 0.00"; emission[1].Probabilities() = "0.00 0.25 0.25 0.50"; emission[2].Probabilities() = "0.10 0.40 0.40 0.10"; - HMM hmm(initial, transition, emission); + HMM> hmm(initial, transition, emission); // Now generate some sequences and check that the log-likelihood is the same // as MATLAB gives for this HMM. @@ -1433,7 +1435,7 @@ TEST_CASE("GaussianHMMLoadSaveTest", "[HMMTest]") TEST_CASE("DiscreteHMMLoadSaveTest", "[HMMTest]") { // Create a Discrete HMM, save it, and load it. - std::vector emission(4); + std::vector> emission(4); emission[0].Probabilities() = arma::randu(6); emission[0].Probabilities() /= accu(emission[0].Probabilities()); emission[1].Probabilities() = arma::randu(6); @@ -1443,10 +1445,8 @@ TEST_CASE("DiscreteHMMLoadSaveTest", "[HMMTest]") emission[3].Probabilities() = arma::randu(6); emission[3].Probabilities() /= accu(emission[3].Probabilities()); - // Create HMM object. - HMM hmm(3, DiscreteDistribution(3)); - + HMM> hmm(3, DiscreteDistribution<>(3)); for (size_t j = 0; j < hmm.Emission().size(); ++j) { @@ -1462,7 +1462,7 @@ TEST_CASE("DiscreteHMMLoadSaveTest", "[HMMTest]") } // Load the HMM. - HMM hmm2(3, DiscreteDistribution(3)); + HMM> hmm2(3, DiscreteDistribution<>(3)); { std::ifstream ifs("test-hmm-save.xml"); cereal::XMLInputArchive ar(ifs); @@ -1483,7 +1483,7 @@ TEST_CASE("DiscreteHMMLoadSaveTest", "[HMMTest]") */ TEST_CASE("HMMTrainReturnLogLikelihood", "[HMMTest]") { - HMM hmm(1, 2); // 1 state, 2 emissions. + HMM> hmm(1, 2); // 1 state, 2 emissions. // Randomize the emission matrix. hmm.Emission()[0].Probabilities() = arma::randu(2); hmm.Emission()[0].Probabilities() /= accu(hmm.Emission()[0].Probabilities()); @@ -1521,18 +1521,18 @@ TEST_CASE("DiagonalGMMHMMPredictTest", "[HMMTest]") std::vector gmms(2); gmms[0] = DiagonalGMM(2, 2); - gmms[0].Component(0) = DiagonalGaussianDistribution("3.25 2.10", + gmms[0].Component(0) = DiagonalGaussianDistribution<>("3.25 2.10", "0.97 1.00"); - gmms[0].Component(1) = DiagonalGaussianDistribution("5.03 7.28", + gmms[0].Component(1) = DiagonalGaussianDistribution<>("5.03 7.28", "1.20 0.89"); gmms[1] = DiagonalGMM(3, 2); gmms[1].Weights() = arma::vec("0.3 0.2 0.5"); - gmms[1].Component(0) = DiagonalGaussianDistribution("-2.48 -3.02", + gmms[1].Component(0) = DiagonalGaussianDistribution<>("-2.48 -3.02", "1.02 0.80"); - gmms[1].Component(1) = DiagonalGaussianDistribution("-1.24 -2.40", + gmms[1].Component(1) = DiagonalGaussianDistribution<>("-1.24 -2.40", "0.85 0.78"); - gmms[1].Component(2) = DiagonalGaussianDistribution("-5.68 -4.83", + gmms[1].Component(2) = DiagonalGaussianDistribution<>("-5.68 -4.83", "1.42 0.96"); // Initial probabilities. @@ -1594,14 +1594,14 @@ TEST_CASE("DiagonalGMMHMMPredictTest", "[HMMTest]") TEST_CASE("DiagonalGMMHMMGenerateTest", "[HMMTest]") { // Build the model. - HMM hmm(3, DiagonalGaussianDistribution(2)); + HMM> hmm(3, DiagonalGaussianDistribution<>(2)); hmm.Transition() = arma::mat("0.2 0.3 0.8;" "0.4 0.5 0.1;" "0.4 0.2 0.1"); - hmm.Emission()[0] = DiagonalGaussianDistribution("0.0 0.0", "1.0 0.7"); - hmm.Emission()[1] = DiagonalGaussianDistribution("1.0 1.0", "0.7 0.5"); - hmm.Emission()[2] = DiagonalGaussianDistribution("-3.0 2.0", "2.0 0.3"); + hmm.Emission()[0] = DiagonalGaussianDistribution<>("0.0 0.0", "1.0 0.7"); + hmm.Emission()[1] = DiagonalGaussianDistribution<>("1.0 1.0", "0.7 0.5"); + hmm.Emission()[2] = DiagonalGaussianDistribution<>("-3.0 2.0", "2.0 0.3"); // Now we will generate a long sequence. std::vector observations(1); @@ -1611,7 +1611,8 @@ TEST_CASE("DiagonalGMMHMMGenerateTest", "[HMMTest]") hmm.Generate(10000, observations[0], states[0], 1); // Build the hmm2. - HMM hmm2(3, DiagonalGaussianDistribution(2)); + HMM> hmm2(3, + DiagonalGaussianDistribution<>(2)); // Now estimate the HMM from the generated sequence. hmm2.Train(observations, states); @@ -1636,7 +1637,7 @@ TEST_CASE("DiagonalGMMHMMGenerateTest", "[HMMTest]") TEST_CASE("DiagonalGMMHMMOneGaussianOneStateTrainingTest", "[HMMTest]") { // Create a Gaussian distribution with diagonal covariance. - DiagonalGaussianDistribution d("2.05 3.45", "0.89 1.05"); + DiagonalGaussianDistribution<> d("2.05 3.45", "0.89 1.05"); // Make a sequence of observations. std::vector observations(1, arma::mat(2, 5000)); @@ -1675,10 +1676,10 @@ TEST_CASE("DiagonalGMMHMMOneGaussianUnlabeledTrainingTest", "[HMMTest]") { // Create a sequence of DiagonalGMMs. Each GMM has one gaussian distribution. std::vector gmms(2, DiagonalGMM(1, 2)); - gmms[0].Component(0) = DiagonalGaussianDistribution("1.25 2.10", + gmms[0].Component(0) = DiagonalGaussianDistribution<>("1.25 2.10", "0.97 1.00"); - gmms[1].Component(0) = DiagonalGaussianDistribution("-2.48 -3.02", + gmms[1].Component(0) = DiagonalGaussianDistribution<>("-2.48 -3.02", "1.02 0.80"); // Transition matrix. @@ -1749,13 +1750,13 @@ TEST_CASE("DiagonalGMMHMMOneGaussianLabeledTrainingTest", "[HMMTest]") { // Create a sequence of DiagonalGMMs. std::vector gmms(3, DiagonalGMM(1, 2)); - gmms[0].Component(0) = DiagonalGaussianDistribution("5.25 7.10", + gmms[0].Component(0) = DiagonalGaussianDistribution<>("5.25 7.10", "0.97 1.00"); - gmms[1].Component(0) = DiagonalGaussianDistribution("4.48 6.02", + gmms[1].Component(0) = DiagonalGaussianDistribution<>("4.48 6.02", "1.02 0.80"); - gmms[2].Component(0) = DiagonalGaussianDistribution("-3.28 -5.30", + gmms[2].Component(0) = DiagonalGaussianDistribution<>("-3.28 -5.30", "0.87 1.05"); // Transition matrix. @@ -1833,15 +1834,15 @@ TEST_CASE("DiagonalGMMHMMMultipleGaussiansUnlabeledTrainingTest", "[HMMTest]") // Create a sequence of DiagonalGMMs. std::vector gmms(2, DiagonalGMM(2, 2)); gmms[0].Weights() = arma::vec("0.3 0.7"); - gmms[0].Component(0) = DiagonalGaussianDistribution("8.25 7.10", + gmms[0].Component(0) = DiagonalGaussianDistribution<>("8.25 7.10", "0.97 1.00"); - gmms[0].Component(1) = DiagonalGaussianDistribution("-3.03 -2.28", + gmms[0].Component(1) = DiagonalGaussianDistribution<>("-3.03 -2.28", "1.20 0.89"); gmms[1].Weights() = arma::vec("0.4 0.6"); - gmms[1].Component(0) = DiagonalGaussianDistribution("4.48 6.02", + gmms[1].Component(0) = DiagonalGaussianDistribution<>("4.48 6.02", "1.02 0.80"); - gmms[1].Component(1) = DiagonalGaussianDistribution("-9.24 -8.40", + gmms[1].Component(1) = DiagonalGaussianDistribution<>("-9.24 -8.40", "0.85 1.58"); // Transition matrix. @@ -1940,15 +1941,15 @@ TEST_CASE("DiagonalGMMHMMMultipleGaussiansLabeledTrainingTest", "[HMMTest]") // Create a sequence of DiagonalGMMs. std::vector gmms(2, DiagonalGMM(2, 2)); gmms[0].Weights() = arma::vec("0.3 0.7"); - gmms[0].Component(0) = DiagonalGaussianDistribution("2.25 5.30", + gmms[0].Component(0) = DiagonalGaussianDistribution<>("2.25 5.30", "0.97 1.00"); - gmms[0].Component(1) = DiagonalGaussianDistribution("-3.15 -2.50", + gmms[0].Component(1) = DiagonalGaussianDistribution<>("-3.15 -2.50", "1.20 0.89"); gmms[1].Weights() = arma::vec("0.4 0.6"); - gmms[1].Component(0) = DiagonalGaussianDistribution("-4.48 -6.30", + gmms[1].Component(0) = DiagonalGaussianDistribution<>("-4.48 -6.30", "1.02 0.80"); - gmms[1].Component(1) = DiagonalGaussianDistribution("5.24 2.40", + gmms[1].Component(1) = DiagonalGaussianDistribution<>("5.24 2.40", "0.85 1.58"); // Transition matrix. diff --git a/src/mlpack/tests/main_tests/hmm_generate_test.cpp b/src/mlpack/tests/main_tests/hmm_generate_test.cpp index d46661773f..1a9e24f8ef 100644 --- a/src/mlpack/tests/main_tests/hmm_generate_test.cpp +++ b/src/mlpack/tests/main_tests/hmm_generate_test.cpp @@ -173,15 +173,15 @@ TEST_CASE_METHOD(HMMGenerateTestFixture, h->DiagGMMHMM()->Emission().resize(2); h->DiagGMMHMM()->Emission()[0] = DiagonalGMM(2, 2); h->DiagGMMHMM()->Emission()[0].Weights() = arma::vec("0.2 0.8"); - h->DiagGMMHMM()->Emission()[0].Component(0) = DiagonalGaussianDistribution( + h->DiagGMMHMM()->Emission()[0].Component(0) = DiagonalGaussianDistribution<>( "2.75 1.60", "0.50 0.50"); - h->DiagGMMHMM()->Emission()[0].Component(1) = DiagonalGaussianDistribution( + h->DiagGMMHMM()->Emission()[0].Component(1) = DiagonalGaussianDistribution<>( "6.15 2.51", "1.00 1.50"); h->DiagGMMHMM()->Emission()[1] = DiagonalGMM(2, 2); h->DiagGMMHMM()->Emission()[1].Weights() = arma::vec("0.4 0.6"); - h->DiagGMMHMM()->Emission()[1].Component(0) = DiagonalGaussianDistribution( + h->DiagGMMHMM()->Emission()[1].Component(0) = DiagonalGaussianDistribution<>( "-1.00 -3.42", "0.20 1.00"); - h->DiagGMMHMM()->Emission()[1].Component(1) = DiagonalGaussianDistribution( + h->DiagGMMHMM()->Emission()[1].Component(1) = DiagonalGaussianDistribution<>( "-3.10 -5.05", "1.20 0.80"); // Now that we have a trained HMM model, we can use it to generate a sequence diff --git a/src/mlpack/tests/main_tests/hmm_test_utils.hpp b/src/mlpack/tests/main_tests/hmm_test_utils.hpp index 9631290380..24a449fb19 100644 --- a/src/mlpack/tests/main_tests/hmm_test_utils.hpp +++ b/src/mlpack/tests/main_tests/hmm_test_utils.hpp @@ -36,7 +36,7 @@ struct InitHMMModel } //! Helper function to create discrete HMM. - static void Create(HMM& hmm, + static void Create(HMM>& hmm, vector& trainSeq, size_t states, double tolerance = 1e-05) @@ -53,8 +53,8 @@ struct InitHMMModel maxEmissions = arma::max(maxEmissions, maxSeqs); } - hmm = HMM(size_t(states), - DiscreteDistribution(maxEmissions), tolerance); + hmm = HMM>(size_t(states), + DiscreteDistribution<>(maxEmissions), tolerance); } static void Create(HMM>& hmm, @@ -135,7 +135,7 @@ struct InitHMMModel } //! Helper function for discrete emission distributions. - static void RandomInitialize(vector& e) + static void RandomInitialize(vector>& e) { for (size_t i = 0; i < e.size(); ++i) { diff --git a/src/mlpack/tests/main_tests/hmm_train_test.cpp b/src/mlpack/tests/main_tests/hmm_train_test.cpp index 8e58c3da4f..56748c47f2 100644 --- a/src/mlpack/tests/main_tests/hmm_train_test.cpp +++ b/src/mlpack/tests/main_tests/hmm_train_test.cpp @@ -70,8 +70,8 @@ inline void ApproximatelyEqual(HMMModel& h1, tolerance); // Check if emission dists are equal - std::vector d1 = h1.DiscreteHMM()->Emission(); - std::vector d2 = h2.DiscreteHMM()->Emission(); + std::vector> d1 = h1.DiscreteHMM()->Emission(); + std::vector> d2 = h2.DiscreteHMM()->Emission(); REQUIRE(d1.size() == d2.size()); diff --git a/src/mlpack/tests/main_tests/hmm_viterbi_test.cpp b/src/mlpack/tests/main_tests/hmm_viterbi_test.cpp index 6384a61382..36015be686 100644 --- a/src/mlpack/tests/main_tests/hmm_viterbi_test.cpp +++ b/src/mlpack/tests/main_tests/hmm_viterbi_test.cpp @@ -168,15 +168,15 @@ TEST_CASE_METHOD(HMMViterbiTestFixture, std::vector gmms(2, DiagonalGMM(2, 2)); gmms[0].Weights() = arma::vec("0.2 0.8"); - gmms[0].Component(0) = DiagonalGaussianDistribution("2.75 1.60", + gmms[0].Component(0) = DiagonalGaussianDistribution<>("2.75 1.60", "0.50 0.50"); - gmms[0].Component(1) = DiagonalGaussianDistribution("6.15 2.51", + gmms[0].Component(1) = DiagonalGaussianDistribution<>("6.15 2.51", "1.00 1.50"); gmms[1].Weights() = arma::vec("0.4 0.6"); - gmms[1].Component(0) = DiagonalGaussianDistribution("-1.00 -3.42", + gmms[1].Component(0) = DiagonalGaussianDistribution<>("-1.00 -3.42", "0.20 1.00"); - gmms[1].Component(1) = DiagonalGaussianDistribution("-3.10 -5.05", + gmms[1].Component(1) = DiagonalGaussianDistribution<>("-3.10 -5.05", "1.20 0.80"); // Transition matrix. diff --git a/src/mlpack/tests/mock_categorical_data.hpp b/src/mlpack/tests/mock_categorical_data.hpp index 4c0d094db6..814465f40b 100644 --- a/src/mlpack/tests/mock_categorical_data.hpp +++ b/src/mlpack/tests/mock_categorical_data.hpp @@ -24,21 +24,21 @@ inline void MockCategoricalData(arma::mat& d, // We'll build a spiral dataset plus two noisy categorical features. We need // to build the distributions for the categorical features (they'll be // discrete distributions). - mlpack::DiscreteDistribution c1[5]; + mlpack::DiscreteDistribution<> c1[5]; // The distribution will be automatically normalized. for (size_t i = 0; i < 5; ++i) { std::vector probs; probs.push_back(arma::vec(4, arma::fill::randu)); - c1[i] = mlpack::DiscreteDistribution(probs); + c1[i] = mlpack::DiscreteDistribution<>(probs); } - mlpack::DiscreteDistribution c2[5]; + mlpack::DiscreteDistribution<> c2[5]; for (size_t i = 0; i < 5; ++i) { std::vector probs; probs.push_back(arma::vec(2, arma::fill::randu)); - c2[i] = mlpack::DiscreteDistribution(probs); + c2[i] = mlpack::DiscreteDistribution<>(probs); } arma::mat spiralDataset(4, 4000); diff --git a/src/mlpack/tests/random_test.cpp b/src/mlpack/tests/random_test.cpp index 27a7f78a11..c3a576e06a 100644 --- a/src/mlpack/tests/random_test.cpp +++ b/src/mlpack/tests/random_test.cpp @@ -76,7 +76,7 @@ TEST_CASE("WeightedRandomTest", "[RandomTest]") const size_t iterations = 50000; for (std::vector weightSet : weights) { - DiscreteDistribution d(1); + DiscreteDistribution<> d(1); d.Probabilities(0) = arma::vec(weightSet); std::vector count(weightSet.size(), 0); for (size_t iter = 0; iter < iterations; ++iter)