Merge pull request #3195 from shubham1206agra/header-only

Header Only #2
This commit is contained in:
Omar Shrit
2022-06-09 14:45:48 +01:00
committed by GitHub
107 changed files with 3874 additions and 3912 deletions
+6 -6
View File
@@ -2,17 +2,17 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
discrete_distribution.hpp
discrete_distribution.cpp
discrete_distribution_impl.hpp
gaussian_distribution.hpp
gaussian_distribution.cpp
gaussian_distribution_impl.hpp
laplace_distribution.hpp
laplace_distribution.cpp
laplace_distribution_impl.hpp
regression_distribution.hpp
regression_distribution.cpp
regression_distribution_impl.hpp
gamma_distribution.hpp
gamma_distribution.cpp
gamma_distribution_impl.hpp
diagonal_gaussian_distribution.hpp
diagonal_gaussian_distribution.cpp
diagonal_gaussian_distribution_impl.hpp
)
# add directory name to sources
@@ -15,7 +15,7 @@
#include <mlpack/prereqs.hpp>
namespace mlpack {
namespace distribution {
namespace distribution /** Probability distributions. */ {
//! A single multivariate Gaussian distribution with diagonal covariance.
class DiagonalGaussianDistribution
@@ -153,4 +153,7 @@ class DiagonalGaussianDistribution
} // namespace distribution
} // namespace mlpack
// Include implementation.
#include "diagonal_gaussian_distribution_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file core/dists/diagonal_gaussian_distribution.cpp
* @file core/dists/diagonal_gaussian_distribution_impl.hpp
* @author Kim SangYeon
*
* Implementation of Gaussian distribution class with diagonal covariance.
@@ -9,13 +9,16 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_CORE_DISTRIBUTIONS_DIAGONAL_GAUSSIAN_DISTRIBUTION_IMPL_HPP
#define MLPACK_CORE_DISTRIBUTIONS_DIAGONAL_GAUSSIAN_DISTRIBUTION_IMPL_HPP
#include "diagonal_gaussian_distribution.hpp"
#include <mlpack/methods/gmm/diagonal_constraint.hpp>
using namespace mlpack;
using namespace mlpack::distribution;
namespace mlpack {
namespace distribution /** Probability distributions. */ {
DiagonalGaussianDistribution::DiagonalGaussianDistribution(
inline DiagonalGaussianDistribution::DiagonalGaussianDistribution(
const arma::vec& mean,
const arma::vec& covariance) :
mean(mean)
@@ -23,21 +26,21 @@ DiagonalGaussianDistribution::DiagonalGaussianDistribution(
Covariance(covariance);
}
void DiagonalGaussianDistribution::Covariance(const arma::vec& covariance)
inline void DiagonalGaussianDistribution::Covariance(const arma::vec& covariance)
{
this->invCov = 1 / covariance;
this->logDetCov = arma::accu(log(covariance));
invCov = 1 / covariance;
logDetCov = arma::accu(log(covariance));
this->covariance = covariance;
}
void DiagonalGaussianDistribution::Covariance(arma::vec&& covariance)
inline void DiagonalGaussianDistribution::Covariance(arma::vec&& covariance)
{
this->invCov = 1 / covariance;
this->logDetCov = arma::accu(log(covariance));
invCov = 1 / covariance;
logDetCov = arma::accu(log(covariance));
this->covariance = std::move(covariance);
}
double DiagonalGaussianDistribution::LogProbability(
inline double DiagonalGaussianDistribution::LogProbability(
const arma::vec& observation) const
{
const size_t k = observation.n_elem;
@@ -46,7 +49,7 @@ double DiagonalGaussianDistribution::LogProbability(
return -0.5 * k * log2pi - 0.5 * logDetCov - 0.5 * logExponent(0);
}
void DiagonalGaussianDistribution::LogProbability(
inline void DiagonalGaussianDistribution::LogProbability(
const arma::mat& observations,
arma::vec& logProbabilities) const
{
@@ -63,12 +66,12 @@ void DiagonalGaussianDistribution::LogProbability(
logProbabilities = -0.5 * k * log2pi - 0.5 * logDetCov + logExponents;
}
arma::vec DiagonalGaussianDistribution::Random() const
inline arma::vec DiagonalGaussianDistribution::Random() const
{
return (arma::sqrt(covariance) % arma::randn<arma::vec>(mean.n_elem)) + mean;
}
void DiagonalGaussianDistribution::Train(const arma::mat& observations)
inline void DiagonalGaussianDistribution::Train(const arma::mat& observations)
{
if (observations.n_cols > 1)
{
@@ -95,8 +98,8 @@ void DiagonalGaussianDistribution::Train(const arma::mat& observations)
logDetCov = arma::accu(log(covariance));
}
void DiagonalGaussianDistribution::Train(const arma::mat& observations,
const arma::vec& probabilities)
inline void DiagonalGaussianDistribution::Train(const arma::mat& observations,
const arma::vec& probabilities)
{
if (observations.n_cols > 0)
{
@@ -146,3 +149,8 @@ void DiagonalGaussianDistribution::Train(const arma::mat& observations,
invCov = 1 / covariance;
logDetCov = arma::accu(log(covariance));
}
} // namespace distribution
} // namespace mlpack
#endif
@@ -252,4 +252,7 @@ class DiscreteDistribution
} // namespace distribution
} // namespace mlpack
// Include implementation.
#include "discrete_distribution_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file core/dists/discrete_distribution.cpp
* @file core/dists/discrete_distribution_impl.hpp
* @author Ryan Curtin
* @author Rohan Raj
*
@@ -10,16 +10,19 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_CORE_DISTRIBUTIONS_DISCRETE_DISTRIBUTION_IMPL_HPP
#define MLPACK_CORE_DISTRIBUTIONS_DISCRETE_DISTRIBUTION_IMPL_HPP
#include "discrete_distribution.hpp"
using namespace mlpack;
using namespace mlpack::distribution;
namespace mlpack {
namespace distribution /** Probability distributions. */ {
/**
* Return a randomly generated observation according to the probability
* distribution defined by this object.
*/
arma::vec DiscreteDistribution::Random() const
inline arma::vec DiscreteDistribution::Random() const
{
size_t dimension = probabilities.size();
arma::vec result(dimension);
@@ -53,7 +56,7 @@ arma::vec DiscreteDistribution::Random() const
/**
* Estimate the probability distribution directly from the given observations.
*/
void DiscreteDistribution::Train(const arma::mat& observations)
inline void DiscreteDistribution::Train(const arma::mat& observations)
{
// Make sure the observations have same dimension as the probabilities.
if (observations.n_rows != probabilities.size())
@@ -107,8 +110,8 @@ 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.
*/
void DiscreteDistribution::Train(const arma::mat& observations,
const arma::vec& probObs)
inline void DiscreteDistribution::Train(const arma::mat& observations,
const arma::vec& probObs)
{
// Make sure the observations have same dimension as the probabilities.
if (observations.n_rows != probabilities.size())
@@ -158,3 +161,8 @@ void DiscreteDistribution::Train(const arma::mat& observations,
probabilities[i].fill(1.0 / probabilities[i].n_elem);
}
}
} // namespace distribution
} // namespace mlpack
#endif
+6 -3
View File
@@ -16,8 +16,8 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef _MLPACK_CORE_DISTRIBUTIONS_GAMMA_DISTRIBUTION_HPP
#define _MLPACK_CORE_DISTRIBUTIONS_GAMMA_DISTRIBUTION_HPP
#ifndef MLPACK_CORE_DISTRIBUTIONS_GAMMA_DISTRIBUTION_HPP
#define MLPACK_CORE_DISTRIBUTIONS_GAMMA_DISTRIBUTION_HPP
#include <mlpack/prereqs.hpp>
#include <mlpack/core/math/random.hpp>
@@ -25,7 +25,7 @@
#include <mlpack/core/math/trigamma.hpp>
namespace mlpack {
namespace distribution {
namespace distribution /** Probability distributions. */ {
/**
* This class represents the Gamma distribution. It supports training a Gamma
@@ -232,4 +232,7 @@ class GammaDistribution
} // namespace distribution
} // namespace mlpack
// Include implementation.
#include "gamma_distribution_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file core/dists/gamma_distribution.cpp
* @file core/dists/gamma_distribution_impl.hpp
* @author Yannis Mentekidis
* @author Rohan Raj
*
@@ -10,26 +10,29 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_CORE_DISTRIBUTIONS_GAMMA_DISTRIBUTION_IMPL_HPP
#define MLPACK_CORE_DISTRIBUTIONS_GAMMA_DISTRIBUTION_IMPL_HPP
#include "gamma_distribution.hpp"
using namespace mlpack;
using namespace mlpack::distribution;
namespace mlpack {
namespace distribution /** Probability distributions. */ {
GammaDistribution::GammaDistribution(const size_t dimensionality)
inline GammaDistribution::GammaDistribution(const size_t dimensionality)
{
// Initialize distribution.
alpha.zeros(dimensionality);
beta.zeros(dimensionality);
}
GammaDistribution::GammaDistribution(const arma::mat& data,
const double tol)
inline GammaDistribution::GammaDistribution(const arma::mat& data,
const double tol)
{
Train(data, tol);
}
GammaDistribution::GammaDistribution(const arma::vec& alpha,
const arma::vec& beta)
inline GammaDistribution::GammaDistribution(const arma::vec& alpha,
const arma::vec& beta)
{
if (beta.n_elem != alpha.n_elem)
throw std::runtime_error("Alpha and beta vector dimensions mismatch.");
@@ -47,7 +50,7 @@ inline bool GammaDistribution::Converged(const double aOld,
}
// Fits an alpha and beta parameter to each dimension of the data.
void GammaDistribution::Train(const arma::mat& rdata, const double tol)
inline void GammaDistribution::Train(const arma::mat& rdata, const double tol)
{
// If fittingSet is empty, nothing to do.
if (arma::size(rdata) == arma::size(arma::mat()))
@@ -64,9 +67,9 @@ void GammaDistribution::Train(const arma::mat& rdata, const double tol)
}
// Fits an alpha and beta parameter according to observation probabilities.
void GammaDistribution::Train(const arma::mat& rdata,
const arma::vec& probabilities,
const double tol)
inline void GammaDistribution::Train(const arma::mat& rdata,
const arma::vec& probabilities,
const double tol)
{
// If fittingSet is empty, nothing to do.
if (arma::size(rdata) == arma::size(arma::mat()))
@@ -94,10 +97,10 @@ void GammaDistribution::Train(const arma::mat& rdata,
}
// Fits an alpha and beta parameter to each dimension of the data.
void GammaDistribution::Train(const arma::vec& logMeanxVec,
const arma::vec& meanLogxVec,
const arma::vec& meanxVec,
const double tol)
inline void GammaDistribution::Train(const arma::vec& logMeanxVec,
const arma::vec& meanLogxVec,
const arma::vec& meanxVec,
const double tol)
{
using std::log;
@@ -155,8 +158,8 @@ void GammaDistribution::Train(const arma::vec& logMeanxVec,
}
// Returns the probability of the provided observations.
void GammaDistribution::Probability(const arma::mat& observations,
arma::vec& probabilities) const
inline void GammaDistribution::Probability(const arma::mat& observations,
arma::vec& probabilities) const
{
size_t numObs = observations.n_cols;
@@ -184,15 +187,16 @@ void GammaDistribution::Probability(const arma::mat& observations,
// Returns the probability of one observation (x) for one of the Gamma's
// dimensions.
double GammaDistribution::Probability(double x, size_t dim) const
inline double GammaDistribution::Probability(double 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.
void GammaDistribution::LogProbability(const arma::mat& observations,
arma::vec& logProbabilities) const
inline void GammaDistribution::LogProbability(
const arma::mat& observations,
arma::vec& logProbabilities) const
{
size_t numObs = observations.n_cols;
@@ -221,14 +225,14 @@ void GammaDistribution::LogProbability(const arma::mat& observations,
// Returns the log probability of one observation (x) for one of the Gamma's
// dimensions.
double GammaDistribution::LogProbability(double x, size_t dim) const
inline double GammaDistribution::LogProbability(double 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.
arma::vec GammaDistribution::Random() const
inline arma::vec GammaDistribution::Random() const
{
arma::vec randVec(alpha.n_elem);
@@ -236,8 +240,13 @@ arma::vec GammaDistribution::Random() const
{
std::gamma_distribution<double> dist(alpha(d), beta(d));
// Use the mlpack random object.
randVec(d) = dist(mlpack::math::randGen);
randVec(d) = dist(mlpack::math::RandGen());
}
return randVec;
}
} // namespace distribution
} // namespace mlpack
#endif
@@ -16,7 +16,7 @@
#include <mlpack/prereqs.hpp>
namespace mlpack {
namespace distribution {
namespace distribution /** Probability distributions. */ {
/**
* A single multivariate Gaussian distribution.
@@ -194,4 +194,7 @@ class GaussianDistribution
} // namespace distribution
} // namespace mlpack
// Include implementation.
#include "gaussian_distribution_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file core/dists/gaussian_distribution.cpp
* @file core/dists/gaussian_distribution_impl.hpp
* @author Ryan Curtin
* @author Michael Fox
*
@@ -10,33 +10,37 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_CORE_DISTRIBUTIONS_GAUSSIAN_DISTRIBUTION_IMPL_HPP
#define MLPACK_CORE_DISTRIBUTIONS_GAUSSIAN_DISTRIBUTION_IMPL_HPP
#include "gaussian_distribution.hpp"
#include <mlpack/methods/gmm/positive_definite_constraint.hpp>
using namespace mlpack;
using namespace mlpack::distribution;
namespace mlpack {
namespace distribution /** Probability distributions. */ {
GaussianDistribution::GaussianDistribution(const arma::vec& mean,
const arma::mat& covariance)
: mean(mean), logDetCov(0.0)
inline GaussianDistribution::GaussianDistribution(
const arma::vec& mean,
const arma::mat& covariance) :
mean(mean),
logDetCov(0.0)
{
Covariance(covariance);
}
void GaussianDistribution::Covariance(const arma::mat& covariance)
inline void GaussianDistribution::Covariance(const arma::mat& covariance)
{
this->covariance = covariance;
FactorCovariance();
}
void GaussianDistribution::Covariance(arma::mat&& covariance)
inline void GaussianDistribution::Covariance(arma::mat&& covariance)
{
this->covariance = std::move(covariance);
FactorCovariance();
}
void GaussianDistribution::FactorCovariance()
inline void GaussianDistribution::FactorCovariance()
{
// On Armadillo < 4.500, the "lower" option isn't available.
@@ -68,7 +72,8 @@ void GaussianDistribution::FactorCovariance()
logDetCov *= 2;
}
double GaussianDistribution::LogProbability(const arma::vec& observation) const
inline double GaussianDistribution::LogProbability(
const arma::vec& observation) const
{
const size_t k = observation.n_elem;
const arma::vec diff = mean - observation;
@@ -76,7 +81,7 @@ double GaussianDistribution::LogProbability(const arma::vec& observation) const
return -0.5 * k * log2pi - 0.5 * logDetCov - 0.5 * v(0);
}
arma::vec GaussianDistribution::Random() const
inline arma::vec GaussianDistribution::Random() const
{
return covLower * arma::randn<arma::vec>(mean.n_elem) + mean;
}
@@ -86,7 +91,7 @@ arma::vec GaussianDistribution::Random() const
*
* @param observations List of observations.
*/
void GaussianDistribution::Train(const arma::mat& observations)
inline void GaussianDistribution::Train(const arma::mat& observations)
{
if (observations.n_cols > 0)
{
@@ -95,10 +100,7 @@ void GaussianDistribution::Train(const arma::mat& observations)
}
else // This will end up just being empty.
{
// TODO(stephentu): why do we allow this case? why not throw an error?
mean.zeros(0);
covariance.zeros(0);
return;
Log::Fatal << "Observation columns equal to 0." << std::endl;
}
// Calculate the mean.
@@ -130,8 +132,8 @@ void GaussianDistribution::Train(const arma::mat& observations)
* account the probability of each observation actually being from this
* distribution.
*/
void GaussianDistribution::Train(const arma::mat& observations,
const arma::vec& probabilities)
inline void GaussianDistribution::Train(const arma::mat& observations,
const arma::vec& probabilities)
{
if (observations.n_cols > 0)
{
@@ -140,10 +142,7 @@ void GaussianDistribution::Train(const arma::mat& observations,
}
else // This will end up just being empty.
{
// TODO(stephentu): same as above
mean.zeros(0);
covariance.zeros(0);
return;
Log::Fatal << "Observation columns equal to 0." << std::endl;
}
double sumProb = 0;
@@ -185,3 +184,8 @@ void GaussianDistribution::Train(const arma::mat& observations,
FactorCovariance();
}
} // namespace distribution
} // namespace mlpack
#endif
@@ -14,8 +14,10 @@
#ifndef MLPACK_CORE_DISTRIBUTIONS_LAPLACE_DISTRIBUTION_HPP
#define MLPACK_CORE_DISTRIBUTIONS_LAPLACE_DISTRIBUTION_HPP
#include <mlpack/prereqs.hpp>
namespace mlpack {
namespace distribution {
namespace distribution /** Probability distributions. */ {
/**
* The multivariate Laplace distribution centered at 0 has pdf
@@ -189,4 +191,7 @@ class LaplaceDistribution
} // namespace distribution
} // namespace mlpack
// Include implementation.
#include "laplace_distribution_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/*
* @file core/dists/laplace_distribution.cpp
* @file core/dists/laplace_distribution_impl.hpp
* @author Zhihao Lou
* @author Rohan Raj
*
@@ -10,17 +10,19 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include <mlpack/prereqs.hpp>
#ifndef MLPACK_CORE_DISTRIBUTIONS_LAPLACE_DISTRIBUTION_IMPL_HPP
#define MLPACK_CORE_DISTRIBUTIONS_LAPLACE_DISTRIBUTION_IMPL_HPP
#include "laplace_distribution.hpp"
using namespace mlpack;
using namespace mlpack::distribution;
namespace mlpack {
namespace distribution /** Probability distributions. */ {
/**
* Return the log probability of the given observation.
*/
double LaplaceDistribution::LogProbability(const arma::vec& observation) const
inline double LaplaceDistribution::LogProbability(
const arma::vec& observation) const
{
// Evaluate the PDF of the Laplace distribution to determine
// the log probability.
@@ -33,8 +35,8 @@ double LaplaceDistribution::LogProbability(const arma::vec& observation) const
* @param x List of observations.
* @param probabilities Output probabilities for each input observation.
*/
void LaplaceDistribution::Probability(const arma::mat& x,
arma::vec& probabilities) const
inline void LaplaceDistribution::Probability(const arma::mat& x,
arma::vec& probabilities) const
{
probabilities.set_size(x.n_cols);
for (size_t i = 0; i < x.n_cols; ++i)
@@ -48,7 +50,7 @@ void LaplaceDistribution::Probability(const arma::mat& x,
*
* @param observations List of observations.
*/
void LaplaceDistribution::Estimate(const arma::mat& observations)
inline void LaplaceDistribution::Estimate(const arma::mat& 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
@@ -81,8 +83,8 @@ void LaplaceDistribution::Estimate(const arma::mat& observations)
* taking into account the probability of each observation actually being from
* this distribution.
*/
void LaplaceDistribution::Estimate(const arma::mat& observations,
const arma::vec& probabilities)
inline void LaplaceDistribution::Estimate(const arma::mat& observations,
const arma::vec& probabilities)
{
// I am not completely sure that this change results in a valid maximum
// likelihood estimator given probabilities of points.
@@ -99,3 +101,8 @@ void LaplaceDistribution::Estimate(const arma::mat& observations,
scale += probabilities(i) * arma::norm(observations.col(i) - mean, 2);
scale /= arma::accu(probabilities);
}
} // namespace distribution
} // namespace mlpack
#endif
@@ -18,7 +18,7 @@
#include <mlpack/methods/linear_regression/linear_regression.hpp>
namespace mlpack {
namespace distribution {
namespace distribution /** Probability distributions. */ {
/**
* A class that represents a univariate conditionally Gaussian distribution.
@@ -160,4 +160,7 @@ class RegressionDistribution
} // namespace distribution
} // namespace mlpack
// Include implementation.
#include "regression_distribution_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file core/dists/regression_distribution.cpp
* @file core/dists/regression_distribution_impl.hpp
* @author Michael Fox
*
* Implementation of conditional Gaussian distribution for HMM regression
@@ -11,17 +11,20 @@
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_CORE_DISTRIBUTIONS_REGRESSION_DISTRIBUTION_IMPL_HPP
#define MLPACK_CORE_DISTRIBUTIONS_REGRESSION_DISTRIBUTION_IMPL_HPP
#include "regression_distribution.hpp"
using namespace mlpack;
using namespace mlpack::distribution;
namespace mlpack {
namespace distribution /** Probability distributions. */ {
/**
* Estimate parameters using provided observation weights.
*
* @param observations List of observations.
*/
void RegressionDistribution::Train(const arma::mat& observations)
inline void RegressionDistribution::Train(const arma::mat& observations)
{
regression::LinearRegression lr(observations.rows(1, observations.n_rows - 1),
arma::rowvec(observations.row(0)), 0, true);
@@ -36,14 +39,14 @@ void RegressionDistribution::Train(const arma::mat& observations)
*
* @param weights Probability that given observation is from distribution.
*/
void RegressionDistribution::Train(const arma::mat& observations,
const arma::vec& weights)
inline void RegressionDistribution::Train(const arma::mat& observations,
const arma::vec& weights)
{
Train(observations, arma::rowvec(weights.t()));
}
void RegressionDistribution::Train(const arma::mat& observations,
const arma::rowvec& weights)
inline void RegressionDistribution::Train(const arma::mat& observations,
const arma::rowvec& weights)
{
regression::LinearRegression lr(observations.rows(1, observations.n_rows - 1),
arma::rowvec(observations.row(0)), weights, 0, true);
@@ -58,23 +61,29 @@ void RegressionDistribution::Train(const arma::mat& observations,
*
* @param observation Point to evaluate probability at.
*/
double RegressionDistribution::Probability(const arma::vec& observation) const
inline double RegressionDistribution::Probability(
const arma::vec& observation) const
{
arma::rowvec fitted;
rf.Predict(observation.rows(1, observation.n_rows-1), fitted);
return err.Probability(observation(0)-fitted.t());
}
void RegressionDistribution::Predict(const arma::mat& points,
arma::vec& predictions) const
inline void RegressionDistribution::Predict(const arma::mat& points,
arma::vec& predictions) const
{
arma::rowvec rowPredictions;
Predict(points, rowPredictions);
predictions = rowPredictions.t();
}
void RegressionDistribution::Predict(const arma::mat& points,
arma::rowvec& predictions) const
inline void RegressionDistribution::Predict(const arma::mat& points,
arma::rowvec& predictions) const
{
rf.Predict(points, predictions);
}
} // namespace distribution
} // namespace mlpack
#endif
+1 -2
View File
@@ -3,7 +3,7 @@
set(SOURCES
clamp.hpp
columns_to_blocks.hpp
columns_to_blocks.cpp
columns_to_blocks_impl.hpp
digamma.hpp
lin_alg.hpp
lin_alg_impl.hpp
@@ -14,7 +14,6 @@ set(SOURCES
multiply_slices.hpp
quantile.hpp
random.hpp
random.cpp
random_basis.hpp
random_basis_impl.hpp
range.hpp
@@ -224,4 +224,7 @@ class ColumnsToBlocks
} // namespace math
} // namespace mlpack
// Include implementation.
#include "columns_to_blocks_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file core/math/columns_to_blocks.cpp
* @file core/math/columns_to_blocks_impl.hpp
* @author Tham Ngap Wei
*
* Implementation of the ColumnsToBlocks class.
@@ -9,15 +9,18 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_NN_COLUMNS_TO_BLOCKS_IMPL_HPP
#define MLPACK_METHODS_NN_COLUMNS_TO_BLOCKS_IMPL_HPP
#include "columns_to_blocks.hpp"
namespace mlpack {
namespace math {
ColumnsToBlocks::ColumnsToBlocks(const size_t rows,
const size_t cols,
const size_t blockHeight,
const size_t blockWidth) :
inline ColumnsToBlocks::ColumnsToBlocks(const size_t rows,
const size_t cols,
const size_t blockHeight,
const size_t blockWidth) :
blockHeight(blockHeight),
blockWidth(blockWidth),
bufSize(1),
@@ -30,15 +33,16 @@ ColumnsToBlocks::ColumnsToBlocks(const size_t rows,
{
}
bool ColumnsToBlocks::IsPerfectSquare(const size_t value) const
inline bool ColumnsToBlocks::IsPerfectSquare(const size_t value) const
{
const size_t root = (size_t) std::round(std::sqrt(value));
return (value == root * root);
}
void ColumnsToBlocks::Transform(const arma::mat& maximalInputs,
arma::mat& output)
inline void ColumnsToBlocks::Transform(const arma::mat& maximalInputs,
arma::mat& output)
{
//! TODO: Maybe replace std::runtime_error with Log::Fatal.
if (!IsPerfectSquare(maximalInputs.n_rows))
{
throw std::runtime_error("maximalInputs.n_rows should be perfect square");
@@ -97,3 +101,5 @@ void ColumnsToBlocks::Transform(const arma::mat& maximalInputs,
} // namespace math
} // namespace mlpack
#endif
-25
View File
@@ -1,25 +0,0 @@
/**
* @file core/math/random.cpp
*
* Declarations of global random number generators.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include <random>
#include <mlpack/mlpack_export.hpp>
namespace mlpack {
namespace math {
// Global random object.
MLPACK_EXPORT std::mt19937 randGen;
// Global uniform distribution.
MLPACK_EXPORT std::uniform_real_distribution<> randUniformDist(0.0, 1.0);
// Global normal distribution.
MLPACK_EXPORT std::normal_distribution<> randNormalDist(0.0, 1.0);
} // namespace math
} // namespace mlpack
+29 -15
View File
@@ -23,12 +23,26 @@ namespace math /** Miscellaneous math routines. */ {
* correctly on Windows.
*/
// Global random object.
extern MLPACK_EXPORT std::mt19937 randGen;
// Global uniform distribution.
extern MLPACK_EXPORT std::uniform_real_distribution<> randUniformDist;
// Global normal distribution.
extern MLPACK_EXPORT std::normal_distribution<> randNormalDist;
//! Global random object.
inline std::mt19937& RandGen()
{
static thread_local std::mt19937 randGen;
return randGen;
}
//! Global uniform distribution.
inline std::uniform_real_distribution<>& RandUniformDist()
{
static thread_local std::uniform_real_distribution<> randUniformDist(0.0, 1.0);
return randUniformDist;
}
//! Global normal distribution.
inline std::normal_distribution<>& RandNormalDist()
{
static thread_local std::normal_distribution<> randNormalDist(0.0, 1.0);
return randNormalDist;
}
/**
* Set the random seed used by the random functions (Random() and RandInt()).
@@ -40,7 +54,7 @@ extern MLPACK_EXPORT std::normal_distribution<> randNormalDist;
inline void RandomSeed(const size_t seed)
{
#if (!defined(BINDING_TYPE) || BINDING_TYPE != BINDING_TYPE_TEST)
randGen.seed((uint32_t) seed);
RandGen().seed((uint32_t) seed);
#if (BINDING_TYPE == BINDING_TYPE_R)
// To suppress Found 'srand', possibly from 'srand' (C).
(void) seed;
@@ -64,14 +78,14 @@ inline void RandomSeed(const size_t seed)
inline void FixedRandomSeed()
{
const static size_t seed = rand();
randGen.seed((uint32_t) seed);
RandGen().seed((uint32_t) seed);
srand((unsigned int) seed);
arma::arma_rng::set_seed(seed);
}
inline void CustomRandomSeed(const size_t seed)
{
randGen.seed((uint32_t) seed);
RandGen().seed((uint32_t) seed);
srand((unsigned int) seed);
arma::arma_rng::set_seed(seed);
}
@@ -82,7 +96,7 @@ inline void CustomRandomSeed(const size_t seed)
*/
inline double Random()
{
return randUniformDist(randGen);
return RandUniformDist()(RandGen());
}
/**
@@ -90,7 +104,7 @@ inline double Random()
*/
inline double Random(const double lo, const double hi)
{
return lo + (hi - lo) * randUniformDist(randGen);
return lo + (hi - lo) * RandUniformDist()(RandGen());
}
/**
@@ -109,7 +123,7 @@ inline double RandBernoulli(const double input)
*/
inline int RandInt(const int hiExclusive)
{
return (int) std::floor((double) hiExclusive * randUniformDist(randGen));
return (int) std::floor((double) hiExclusive * RandUniformDist()(RandGen()));
}
/**
@@ -118,7 +132,7 @@ inline int RandInt(const int hiExclusive)
inline int RandInt(const int lo, const int hiExclusive)
{
return lo + (int) std::floor((double) (hiExclusive - lo)
* randUniformDist(randGen));
* RandUniformDist()(RandGen()));
}
/**
@@ -126,7 +140,7 @@ inline int RandInt(const int lo, const int hiExclusive)
*/
inline double RandNormal()
{
return randNormalDist(randGen);
return RandNormalDist()(RandGen());
}
/**
@@ -138,7 +152,7 @@ inline double RandNormal()
*/
inline double RandNormal(const double mean, const double variance)
{
return variance * randNormalDist(randGen) + mean;
return variance * RandNormalDist()(RandGen()) + mean;
}
/**
+1 -1
View File
@@ -32,7 +32,7 @@ set(SOURCES
cellbound.hpp
cellbound_impl.hpp
cosine_tree/cosine_tree.hpp
cosine_tree/cosine_tree.cpp
cosine_tree/cosine_tree_impl.hpp
cover_tree.hpp
cover_tree/cover_tree.hpp
cover_tree/cover_tree_impl.hpp
@@ -13,6 +13,8 @@
#define MLPACK_CORE_TREE_COSINE_TREE_COSINE_TREE_HPP
#include <mlpack/prereqs.hpp>
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/math/quantile.hpp>
namespace mlpack {
namespace tree {
@@ -285,4 +287,7 @@ class CompareCosineNode
} // namespace tree
} // namespace mlpack
// Include implementation.
#include "cosine_tree_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file core/tree/cosine_tree/cosine_tree.cpp
* @file core/tree/cosine_tree/cosine_tree_impl.hpp
* @author Siddharth Agrawal
*
* Implementation of cosine tree.
@@ -9,15 +9,15 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "cosine_tree.hpp"
#include <mlpack/core/util/log.hpp>
#ifndef MLPACK_CORE_TREE_COSINE_TREE_COSINE_TREE_IMPL_HPP
#define MLPACK_CORE_TREE_COSINE_TREE_COSINE_TREE_IMPL_HPP
#include <mlpack/core/math/quantile.hpp>
#include "cosine_tree.hpp"
namespace mlpack {
namespace tree {
CosineTree::CosineTree(const arma::mat& dataset) :
inline CosineTree::CosineTree(const arma::mat& dataset) :
dataset(&dataset),
parent(NULL),
left(NULL),
@@ -46,8 +46,8 @@ CosineTree::CosineTree(const arma::mat& dataset) :
splitPointIndex = ColumnSampleLS();
}
CosineTree::CosineTree(CosineTree& parentNode,
const std::vector<size_t>& subIndices) :
inline CosineTree::CosineTree(CosineTree& parentNode,
const std::vector<size_t>& subIndices) :
dataset(&parentNode.GetDataset()),
parent(&parentNode),
left(NULL),
@@ -75,9 +75,9 @@ CosineTree::CosineTree(CosineTree& parentNode,
splitPointIndex = ColumnSampleLS();
}
CosineTree::CosineTree(const arma::mat& dataset,
const double epsilon,
const double delta) :
inline CosineTree::CosineTree(const arma::mat& dataset,
const double epsilon,
const double delta) :
dataset(&dataset),
delta(delta),
left(NULL),
@@ -158,7 +158,7 @@ CosineTree::CosineTree(const arma::mat& dataset,
}
//! Copy the given tree.
CosineTree::CosineTree(const CosineTree& other) :
inline CosineTree::CosineTree(const CosineTree& other) :
// Copy matrix, but only if we are the root.
dataset((other.parent == NULL) ? new arma::mat(*other.dataset) : NULL),
delta(other.delta),
@@ -211,7 +211,7 @@ CosineTree::CosineTree(const CosineTree& other) :
}
//! Copy assignment operator: copy the given other tree.
CosineTree& CosineTree::operator=(const CosineTree& other)
inline CosineTree& CosineTree::operator=(const CosineTree& other)
{
// Return if it's the same tree.
if (this == &other)
@@ -279,7 +279,7 @@ CosineTree& CosineTree::operator=(const CosineTree& other)
}
//! Move the given tree.
CosineTree::CosineTree(CosineTree&& other) :
inline CosineTree::CosineTree(CosineTree&& other) :
dataset(other.dataset),
delta(std::move(other.delta)),
parent(other.parent),
@@ -314,7 +314,7 @@ CosineTree::CosineTree(CosineTree&& other) :
}
//! Move assignment operator: take ownership of the given tree.
CosineTree& CosineTree::operator=(CosineTree&& other)
inline CosineTree& CosineTree::operator=(CosineTree&& other)
{
// Return if it's the same tree.
if (this == &other)
@@ -361,7 +361,7 @@ CosineTree& CosineTree::operator=(CosineTree&& other)
return *this;
}
CosineTree::~CosineTree()
inline CosineTree::~CosineTree()
{
if (localDataset)
delete dataset;
@@ -371,10 +371,10 @@ CosineTree::~CosineTree()
delete right;
}
void CosineTree::ModifiedGramSchmidt(CosineNodeQueue& treeQueue,
arma::vec& centroid,
arma::vec& newBasisVector,
arma::vec* addBasisVector)
inline void CosineTree::ModifiedGramSchmidt(CosineNodeQueue& treeQueue,
arma::vec& centroid,
arma::vec& newBasisVector,
arma::vec* addBasisVector)
{
// Set new basis vector to centroid.
newBasisVector = centroid;
@@ -405,10 +405,10 @@ void CosineTree::ModifiedGramSchmidt(CosineNodeQueue& treeQueue,
newBasisVector /= arma::norm(newBasisVector, 2);
}
double CosineTree::MonteCarloError(CosineTree* node,
CosineNodeQueue& treeQueue,
arma::vec* addBasisVector1,
arma::vec* addBasisVector2)
inline double CosineTree::MonteCarloError(CosineTree* node,
CosineNodeQueue& treeQueue,
arma::vec* addBasisVector1,
arma::vec* addBasisVector2)
{
std::vector<size_t> sampledIndices;
arma::vec probabilities;
@@ -489,7 +489,7 @@ double CosineTree::MonteCarloError(CosineTree* node,
return (node->FrobNormSquared() - lowerBound);
}
void CosineTree::ConstructBasis(CosineNodeQueue& treeQueue)
inline void CosineTree::ConstructBasis(CosineNodeQueue& treeQueue)
{
// Initialize basis as matrix of zeros.
basis.zeros(dataset->n_rows, treeQueue.size());
@@ -507,7 +507,7 @@ void CosineTree::ConstructBasis(CosineNodeQueue& treeQueue)
}
}
void CosineTree::CosineNodeSplit()
inline void CosineTree::CosineNodeSplit()
{
// If less than two points, splitting does not make sense---there is nothing
// to split.
@@ -544,9 +544,9 @@ void CosineTree::CosineNodeSplit()
right = new CosineTree(*this, rightIndices);
}
void CosineTree::ColumnSamplesLS(std::vector<size_t>& sampledIndices,
arma::vec& probabilities,
size_t numSamples)
inline void CosineTree::ColumnSamplesLS(std::vector<size_t>& sampledIndices,
arma::vec& probabilities,
size_t numSamples)
{
// Initialize the cumulative distribution vector size.
arma::vec cDistribution;
@@ -576,7 +576,7 @@ void CosineTree::ColumnSamplesLS(std::vector<size_t>& sampledIndices,
}
}
size_t CosineTree::ColumnSampleLS()
inline size_t CosineTree::ColumnSampleLS()
{
// If only one element is present, there can only be one sample.
if (numColumns < 2)
@@ -603,10 +603,10 @@ size_t CosineTree::ColumnSampleLS()
return BinarySearch(cDistribution, randValue, start, end);
}
size_t CosineTree::BinarySearch(arma::vec& cDistribution,
double value,
size_t start,
size_t end)
inline size_t CosineTree::BinarySearch(arma::vec& cDistribution,
double value,
size_t start,
size_t end)
{
size_t pivot = (start + end) / 2;
@@ -631,7 +631,7 @@ size_t CosineTree::BinarySearch(arma::vec& cDistribution,
}
}
void CosineTree::CalculateCosines(arma::vec& cosines)
inline void CosineTree::CalculateCosines(arma::vec& cosines)
{
// Initialize cosine vector as a vector of zeros.
cosines.zeros(numColumns);
@@ -653,7 +653,7 @@ void CosineTree::CalculateCosines(arma::vec& cosines)
}
}
void CosineTree::CalculateCentroid()
inline void CosineTree::CalculateCentroid()
{
// Initialize centroid as vector of zeros.
centroid.zeros(dataset->n_rows);
@@ -668,3 +668,5 @@ void CosineTree::CalculateCentroid()
} // namespace tree
} // namespace mlpack
#endif
+1 -1
View File
@@ -4,7 +4,7 @@ set(SOURCES
adaboost.hpp
adaboost_impl.hpp
adaboost_model.hpp
adaboost_model.cpp
adaboost_model_impl.hpp
)
# Add directory name to sources.
@@ -126,4 +126,7 @@ class AdaBoostModel
} // namespace adaboost
} // namespace mlpack
// Include implementation.
#include "adaboost_model_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/adaboost/adaboost_model.cpp
* @file methods/adaboost/adaboost_model_impl.hpp
* @author Ryan Curtin
*
* A serializable AdaBoost model, used by the main program.
@@ -9,18 +9,17 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_ADABOOST_ADABOOST_MODEL_IMPL_HPP
#define MLPACK_METHODS_ADABOOST_ADABOOST_MODEL_IMPL_HPP
#include "adaboost.hpp"
#include "adaboost_model.hpp"
using namespace mlpack;
using namespace std;
using namespace arma;
using namespace mlpack::adaboost;
using namespace mlpack::tree;
using namespace mlpack::perceptron;
namespace mlpack {
namespace adaboost {
//! Create an empty AdaBoost model.
AdaBoostModel::AdaBoostModel() :
inline AdaBoostModel::AdaBoostModel() :
weakLearnerType(0),
dsBoost(NULL),
pBoost(NULL),
@@ -30,8 +29,8 @@ AdaBoostModel::AdaBoostModel() :
}
//! Create the AdaBoost model with the given mappings and type.
AdaBoostModel::AdaBoostModel(
const Col<size_t>& mappings,
inline AdaBoostModel::AdaBoostModel(
const arma::Col<size_t>& mappings,
const size_t weakLearnerType) :
mappings(mappings),
weakLearnerType(weakLearnerType),
@@ -43,20 +42,20 @@ AdaBoostModel::AdaBoostModel(
}
//! Copy constructor.
AdaBoostModel::AdaBoostModel(const AdaBoostModel& other) :
inline AdaBoostModel::AdaBoostModel(const AdaBoostModel& other) :
mappings(other.mappings),
weakLearnerType(other.weakLearnerType),
dsBoost(other.dsBoost == NULL ? NULL :
new AdaBoost<ID3DecisionStump>(*other.dsBoost)),
new AdaBoost<tree::ID3DecisionStump>(*other.dsBoost)),
pBoost(other.pBoost == NULL ? NULL :
new AdaBoost<Perceptron<>>(*other.pBoost)),
new AdaBoost<perceptron::Perceptron<>>(*other.pBoost)),
dimensionality(other.dimensionality)
{
// Nothing to do.
}
//! Move constructor.
AdaBoostModel::AdaBoostModel(AdaBoostModel&& other) :
inline AdaBoostModel::AdaBoostModel(AdaBoostModel&& other) :
mappings(std::move(other.mappings)),
weakLearnerType(other.weakLearnerType),
dsBoost(other.dsBoost),
@@ -70,7 +69,7 @@ AdaBoostModel::AdaBoostModel(AdaBoostModel&& other) :
}
//! Copy assignment operator.
AdaBoostModel& AdaBoostModel::operator=(const AdaBoostModel& other)
inline AdaBoostModel& AdaBoostModel::operator=(const AdaBoostModel& other)
{
if (this != &other)
{
@@ -79,11 +78,11 @@ AdaBoostModel& AdaBoostModel::operator=(const AdaBoostModel& other)
delete dsBoost;
dsBoost = (other.dsBoost == NULL) ? NULL :
new AdaBoost<ID3DecisionStump>(*other.dsBoost);
new AdaBoost<tree::ID3DecisionStump>(*other.dsBoost);
delete pBoost;
pBoost = (other.pBoost == NULL) ? NULL :
new AdaBoost<Perceptron<>>(*other.pBoost);
new AdaBoost<perceptron::Perceptron<>>(*other.pBoost);
dimensionality = other.dimensionality;
}
@@ -91,7 +90,7 @@ AdaBoostModel& AdaBoostModel::operator=(const AdaBoostModel& other)
}
//! Move assignment operator.
AdaBoostModel& AdaBoostModel::operator=(AdaBoostModel&& other)
inline AdaBoostModel& AdaBoostModel::operator=(AdaBoostModel&& other)
{
if (this != &other)
{
@@ -109,40 +108,40 @@ AdaBoostModel& AdaBoostModel::operator=(AdaBoostModel&& other)
return *this;
}
AdaBoostModel::~AdaBoostModel()
inline AdaBoostModel::~AdaBoostModel()
{
delete dsBoost;
delete pBoost;
}
//! Train the model.
void AdaBoostModel::Train(const mat& data,
const Row<size_t>& labels,
const size_t numClasses,
const size_t iterations,
const double tolerance)
inline void AdaBoostModel::Train(const arma::mat& data,
const arma::Row<size_t>& labels,
const size_t numClasses,
const size_t iterations,
const double tolerance)
{
dimensionality = data.n_rows;
if (weakLearnerType == WeakLearnerTypes::DECISION_STUMP)
{
delete dsBoost;
ID3DecisionStump ds(data, labels, max(labels) + 1);
dsBoost = new AdaBoost<ID3DecisionStump>(data, labels, numClasses, ds,
iterations, tolerance);
tree::ID3DecisionStump ds(data, labels, max(labels) + 1);
dsBoost = new AdaBoost<tree::ID3DecisionStump>(data, labels, numClasses,
ds, iterations, tolerance);
}
else if (weakLearnerType == WeakLearnerTypes::PERCEPTRON)
{
delete pBoost;
Perceptron<> p(data, labels, max(labels) + 1);
pBoost = new AdaBoost<Perceptron<>>(data, labels, numClasses, p, iterations,
tolerance);
perceptron::Perceptron<> p(data, labels, max(labels) + 1);
pBoost = new AdaBoost<perceptron::Perceptron<>>(data, labels, numClasses,
p, iterations, tolerance);
}
}
//! Classify test points.
void AdaBoostModel::Classify(const mat& testData,
Row<size_t>& predictions,
mat& probabilities)
inline void AdaBoostModel::Classify(const arma::mat& testData,
arma::Row<size_t>& predictions,
arma::mat& probabilities)
{
if (weakLearnerType == WeakLearnerTypes::DECISION_STUMP)
dsBoost->Classify(testData, predictions, probabilities);
@@ -151,11 +150,16 @@ void AdaBoostModel::Classify(const mat& testData,
}
//! Classify test points.
void AdaBoostModel::Classify(const mat& testData,
Row<size_t>& predictions)
inline void AdaBoostModel::Classify(const arma::mat& testData,
arma::Row<size_t>& predictions)
{
if (weakLearnerType == WeakLearnerTypes::DECISION_STUMP)
dsBoost->Classify(testData, predictions);
else if (weakLearnerType == WeakLearnerTypes::PERCEPTRON)
pBoost->Classify(testData, predictions);
}
} // namespace adaboost
} // namespace mlpack
#endif
@@ -3,7 +3,6 @@
set(SOURCES
bayesian_linear_regression.hpp
bayesian_linear_regression_impl.hpp
bayesian_linear_regression.cpp
)
# add directory name to sources
@@ -1,184 +0,0 @@
/**
* @file methods/bayesian_linear_regression/bayesian_linear_regression.cpp
* @author Clement Mercier
*
* Implementation of Bayesian linear regression.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "bayesian_linear_regression.hpp"
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/util/timers.hpp>
using namespace mlpack;
using namespace mlpack::regression;
BayesianLinearRegression::BayesianLinearRegression(const bool centerData,
const bool scaleData,
const size_t maxIterations,
const double tolerance) :
centerData(centerData),
scaleData(scaleData),
maxIterations(maxIterations),
tolerance(tolerance),
responsesOffset(0.0),
alpha(0.0),
beta(0.0),
gamma(0.0)
{/* Nothing to do */}
double BayesianLinearRegression::Train(const arma::mat& data,
const arma::rowvec& responses)
{
arma::mat phi;
arma::rowvec t;
arma::colvec eigVal;
arma::mat eigVec;
// Preprocess the data. Center and scale.
responsesOffset = CenterScaleData(data, responses, phi, t);
if (!arma::eig_sym(eigVal, eigVec, arma::symmatu(phi * phi.t())))
{
Log::Fatal << "BayesianLinearRegression::Train(): Eigendecomposition "
<< "of covariance failed!" << std::endl;
}
// Compute this quantities once and for all.
const arma::mat eigVecInv = inv(eigVec);
const arma::colvec eigVecInvPhitT = eigVecInv * phi * t.t();
// Initialize the hyperparameters and begin with an infinitely broad prior.
alpha = 1e-6;
beta = 1 / (var(t, 1) * 0.1);
unsigned short i = 0;
double deltaAlpha = 1.0, crit = 1.0;
while ((crit > tolerance) && (i < maxIterations))
{
deltaAlpha = -alpha;
double deltaBeta = -beta;
// Update the solution.
omega = eigVec * diagmat(1 / (eigVal + (alpha / beta))) * eigVecInvPhitT;
// Update alpha.
gamma = sum(eigVal / (alpha / beta + eigVal));
alpha = gamma / dot(omega, omega);
// Update beta.
const arma::rowvec temp = t - omega.t() * phi;
beta = (data.n_cols - gamma) / dot(temp, temp);
// Compute the stopping criterion.
deltaAlpha += alpha;
deltaBeta += beta;
crit = std::abs(deltaAlpha / alpha + deltaBeta / beta);
i++;
}
// Compute the covariance matrix for the uncertainties later.
matCovariance = eigVec * diagmat(1 / (beta * eigVal + alpha)) * eigVecInv;
return RMSE(data, responses);
}
void BayesianLinearRegression::Predict(const arma::mat& points,
arma::rowvec& predictions) const
{
// Center and scale the points before applying the model.
arma::mat matX;
CenterScaleDataPred(points, matX);
predictions = omega.t() * matX + responsesOffset;
}
void BayesianLinearRegression::Predict(const arma::mat& points,
arma::rowvec& predictions,
arma::rowvec& std) const
{
// Center and scale the points before applying the model.
arma::mat matX;
CenterScaleDataPred(points, matX);
predictions = omega.t() * matX + responsesOffset;
// Compute the standard deviation for each point.
std = sqrt(Variance() + sum(matX % (matCovariance * matX), 0));
}
double BayesianLinearRegression::RMSE(const arma::mat& data,
const arma::rowvec& responses) const
{
arma::rowvec predictions;
Predict(data, predictions);
return sqrt(mean(square(responses - predictions)));
}
double BayesianLinearRegression::CenterScaleData(const arma::mat& data,
const arma::rowvec& responses,
arma::mat& dataProc,
arma::rowvec& responsesProc)
{
if (!centerData && !scaleData)
{
dataProc = arma::mat(const_cast<double*>(data.memptr()), data.n_rows,
data.n_cols, false, true);
responsesProc = arma::rowvec(const_cast<double*>(responses.memptr()),
responses.n_elem, false,
true);
}
else if (centerData && !scaleData)
{
dataOffset = mean(data, 1);
responsesOffset = mean(responses);
dataProc = data.each_col() - dataOffset;
responsesProc = responses - responsesOffset;
}
else if (!centerData && scaleData)
{
dataScale = stddev(data, 0, 1);
dataProc = data.each_col() / dataScale;
responsesProc = arma::rowvec(const_cast<double*>(responses.memptr()),
responses.n_elem, false,
true);
}
else
{
dataOffset = mean(data, 1);
dataScale = stddev(data, 0, 1);
responsesOffset = mean(responses);
dataProc = (data.each_col() - dataOffset).each_col() / dataScale;
responsesProc = responses - responsesOffset;
}
return responsesOffset;
}
void BayesianLinearRegression::CenterScaleDataPred(
const arma::mat& data,
arma::mat& dataProc) const
{
if (!centerData && !scaleData)
{
dataProc = arma::mat(const_cast<double*>(data.memptr()), data.n_rows,
data.n_cols, false, true);
}
else if (centerData && !scaleData)
{
dataProc = data.each_col() - dataOffset;
}
else if (!centerData && scaleData)
{
dataProc = data.each_col() / dataScale;
}
else
{
dataProc = (data.each_col() - dataOffset).each_col() / dataScale;
}
}
@@ -16,6 +16,8 @@
#define MLPACK_METHODS_BAYESIAN_LINEAR_REGRESSION_HPP
#include <mlpack/prereqs.hpp>
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/util/timers.hpp>
namespace mlpack {
namespace regression {
@@ -17,6 +17,176 @@
namespace mlpack {
namespace regression {
inline BayesianLinearRegression::BayesianLinearRegression(
const bool centerData,
const bool scaleData,
const size_t maxIterations,
const double tolerance) :
centerData(centerData),
scaleData(scaleData),
maxIterations(maxIterations),
tolerance(tolerance),
responsesOffset(0.0),
alpha(0.0),
beta(0.0),
gamma(0.0)
{/* Nothing to do */}
inline double BayesianLinearRegression::Train(const arma::mat& data,
const arma::rowvec& responses)
{
arma::mat phi;
arma::rowvec t;
arma::colvec eigVal;
arma::mat eigVec;
// Preprocess the data. Center and scale.
responsesOffset = CenterScaleData(data, responses, phi, t);
if (!arma::eig_sym(eigVal, eigVec, arma::symmatu(phi * phi.t())))
{
Log::Fatal << "BayesianLinearRegression::Train(): Eigendecomposition "
<< "of covariance failed!" << std::endl;
}
// Compute this quantities once and for all.
const arma::mat eigVecInv = inv(eigVec);
const arma::colvec eigVecInvPhitT = eigVecInv * phi * t.t();
// Initialize the hyperparameters and begin with an infinitely broad prior.
alpha = 1e-6;
beta = 1 / (var(t, 1) * 0.1);
unsigned short i = 0;
double crit = 1.0;
while ((crit > tolerance) && (i < maxIterations))
{
double deltaAlpha = -alpha;
double deltaBeta = -beta;
// Update the solution.
omega = eigVec * diagmat(1 / (eigVal + (alpha / beta))) * eigVecInvPhitT;
// Update alpha.
gamma = sum(eigVal / (alpha / beta + eigVal));
alpha = gamma / dot(omega, omega);
// Update beta.
const arma::rowvec temp = t - omega.t() * phi;
beta = (data.n_cols - gamma) / dot(temp, temp);
// Compute the stopping criterion.
deltaAlpha += alpha;
deltaBeta += beta;
crit = std::abs(deltaAlpha / alpha + deltaBeta / beta);
i++;
}
// Compute the covariance matrix for the uncertainties later.
matCovariance = eigVec * diagmat(1 / (beta * eigVal + alpha)) * eigVecInv;
return RMSE(data, responses);
}
inline void BayesianLinearRegression::Predict(const arma::mat& points,
arma::rowvec& predictions) const
{
// Center and scale the points before applying the model.
arma::mat matX;
CenterScaleDataPred(points, matX);
predictions = omega.t() * matX + responsesOffset;
}
inline void BayesianLinearRegression::Predict(const arma::mat& points,
arma::rowvec& predictions,
arma::rowvec& std) const
{
// Center and scale the points before applying the model.
arma::mat matX;
CenterScaleDataPred(points, matX);
predictions = omega.t() * matX + responsesOffset;
// Compute the standard deviation for each point.
std = sqrt(Variance() + sum(matX % (matCovariance * matX), 0));
}
inline double BayesianLinearRegression::RMSE(
const arma::mat& data,
const arma::rowvec& responses) const
{
arma::rowvec predictions;
Predict(data, predictions);
return sqrt(mean(square(responses - predictions)));
}
inline double BayesianLinearRegression::CenterScaleData(
const arma::mat& data,
const arma::rowvec& responses,
arma::mat& dataProc,
arma::rowvec& responsesProc)
{
if (!centerData && !scaleData)
{
dataProc = arma::mat(const_cast<double*>(data.memptr()), data.n_rows,
data.n_cols, false, true);
responsesProc = arma::rowvec(const_cast<double*>(responses.memptr()),
responses.n_elem, false,
true);
}
else if (centerData && !scaleData)
{
dataOffset = mean(data, 1);
responsesOffset = mean(responses);
dataProc = data.each_col() - dataOffset;
responsesProc = responses - responsesOffset;
}
else if (!centerData && scaleData)
{
dataScale = stddev(data, 0, 1);
dataProc = data.each_col() / dataScale;
responsesProc = arma::rowvec(const_cast<double*>(responses.memptr()),
responses.n_elem, false,
true);
}
else
{
dataOffset = mean(data, 1);
dataScale = stddev(data, 0, 1);
responsesOffset = mean(responses);
dataProc = (data.each_col() - dataOffset).each_col() / dataScale;
responsesProc = responses - responsesOffset;
}
return responsesOffset;
}
inline void BayesianLinearRegression::CenterScaleDataPred(
const arma::mat& data,
arma::mat& dataProc) const
{
if (!centerData && !scaleData)
{
dataProc = arma::mat(const_cast<double*>(data.memptr()), data.n_rows,
data.n_cols, false, true);
}
else if (centerData && !scaleData)
{
dataProc = data.each_col() - dataOffset;
}
else if (!centerData && scaleData)
{
dataProc = data.each_col() / dataScale;
}
else
{
dataProc = (data.each_col() - dataOffset).each_col() / dataScale;
}
}
/**
* Serialize the Bayesian linear regression model.
*/
@@ -317,7 +317,7 @@ inline double ParallelSGD<ExponentialBackoff>::Optimize(
if (shuffle) // Determine order of visitation.
std::shuffle(visitationOrder.begin(), visitationOrder.end(),
mlpack::math::randGen);
mlpack::math::RandGen());
#pragma omp parallel
{
@@ -2,7 +2,7 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
randomized_block_krylov_svd.hpp
randomized_block_krylov_svd.cpp
randomized_block_krylov_svd_impl.hpp
)
# Add directory name to sources.
@@ -125,4 +125,7 @@ class RandomizedBlockKrylovSVD
} // namespace svd
} // namespace mlpack
// Include implementation.
#include "randomized_block_krylov_svd_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/block_krylov_svd/randomized_block_krylov_svd.cpp
* @file methods/block_krylov_svd/randomized_block_krylov_svd_impl.hpp
* @author Marcus Edel
*
* Implementation of the randomized block krylov SVD method.
@@ -9,19 +9,22 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_BLOCK_KRYLOV_SVD_RANDOMIZED_BLOCK_KRYLOV_SVD_IMPL_HPP
#define MLPACK_METHODS_BLOCK_KRYLOV_SVD_RANDOMIZED_BLOCK_KRYLOV_SVD_IMPL_HPP
#include "randomized_block_krylov_svd.hpp"
namespace mlpack {
namespace svd {
RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t maxIterations,
const size_t rank,
const size_t blockSize) :
inline RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(
const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t maxIterations,
const size_t rank,
const size_t blockSize) :
maxIterations(maxIterations),
blockSize(blockSize)
{
@@ -35,19 +38,20 @@ RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(const arma::mat& data,
}
}
RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(const size_t maxIterations,
const size_t blockSize) :
inline RandomizedBlockKrylovSVD::RandomizedBlockKrylovSVD(
const size_t maxIterations,
const size_t blockSize) :
maxIterations(maxIterations),
blockSize(blockSize)
{
/* Nothing to do here */
}
void RandomizedBlockKrylovSVD::Apply(const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t rank)
inline void RandomizedBlockKrylovSVD::Apply(const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t rank)
{
arma::mat Q, R, block, blockIteration;
@@ -94,3 +98,5 @@ void RandomizedBlockKrylovSVD::Apply(const arma::mat& data,
} // namespace svd
} // namespace mlpack
#endif
-1
View File
@@ -5,7 +5,6 @@ set(SOURCES
cf_impl.hpp
cf_model.hpp
cf_model_impl.hpp
cf_model.cpp
svd_wrapper.hpp
svd_wrapper_impl.hpp
)
-202
View File
@@ -1,202 +0,0 @@
/**
* @file methods/cf/cf_model_impl.hpp
* @author Wenhao Huang
*
* A serializable CF model, used by the main program.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "cf_model.hpp"
namespace mlpack {
namespace cf {
CFModel::CFModel() :
decompositionType(NMF),
normalizationType(NO_NORMALIZATION),
cf(NULL)
{
// Nothing else to do.
}
CFModel::CFModel(const CFModel& other) :
decompositionType(other.decompositionType),
normalizationType(other.normalizationType),
cf(other.cf->Clone())
{
// Nothing else to do.
}
CFModel::CFModel(CFModel&& other) :
decompositionType(other.decompositionType),
normalizationType(other.normalizationType),
cf(std::move(other.cf))
{
// Reset properties of the other one.
other.decompositionType = NMF;
other.normalizationType = NO_NORMALIZATION;
}
CFModel& CFModel::operator=(const CFModel& other)
{
if (this != &other)
{
decompositionType = other.decompositionType;
normalizationType = other.normalizationType;
cf = other.cf->Clone();
}
return *this;
}
CFModel& CFModel::operator=(CFModel&& other)
{
if (this != &other)
{
decompositionType = other.decompositionType;
normalizationType = other.normalizationType;
cf = std::move(other.cf);
// Reset the other object.
other.decompositionType = NMF;
other.normalizationType = NO_NORMALIZATION;
}
return *this;
}
CFModel::~CFModel()
{
delete cf;
}
template<typename DecompositionPolicy>
CFWrapperBase* TrainHelper(const DecompositionPolicy& decomposition,
const CFModel::NormalizationTypes normalizationType,
const arma::mat& data,
const size_t numUsersForSimilarity,
const size_t rank,
const size_t maxIterations,
const double minResidue,
const bool mit)
{
switch (normalizationType)
{
case CFModel::NO_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, NoNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
case CFModel::ITEM_MEAN_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, ItemMeanNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
case CFModel::USER_MEAN_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, UserMeanNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
case CFModel::OVERALL_MEAN_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, OverallMeanNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
case CFModel::Z_SCORE_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, ZScoreNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
}
// This shouldn't ever happen.
return NULL;
}
void CFModel::Train(const arma::mat& data,
const size_t numUsersForSimilarity,
const size_t rank,
const size_t maxIterations,
const double minResidue,
const bool mit)
{
// Delete the current CFType object, if there is one.
delete cf;
switch (decompositionType)
{
case NMF:
cf = TrainHelper(NMFPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case BATCH_SVD:
cf = TrainHelper(BatchSVDPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case RANDOMIZED_SVD:
cf = TrainHelper(RandomizedSVDPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case REG_SVD:
cf = TrainHelper(RegSVDPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case SVD_COMPLETE:
cf = TrainHelper(SVDCompletePolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case SVD_INCOMPLETE:
cf = TrainHelper(SVDIncompletePolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case BIAS_SVD:
cf = TrainHelper(BiasSVDPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case SVD_PLUS_PLUS:
cf = TrainHelper(SVDPlusPlusPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
}
}
//! Make predictions.
void CFModel::Predict(const NeighborSearchTypes nsType,
const InterpolationTypes interpolationType,
const arma::Mat<size_t>& combinations,
arma::vec& predictions)
{
cf->Predict(nsType, interpolationType, combinations, predictions);
}
//! Compute recommendations for queried users.
void CFModel::GetRecommendations(const NeighborSearchTypes nsType,
const InterpolationTypes interpolationType,
const size_t numRecs,
arma::Mat<size_t>& recommendations,
const arma::Col<size_t>& users)
{
cf->GetRecommendations(nsType, interpolationType, numRecs, recommendations,
users);
}
//! Compute recommendations for all users.
void CFModel::GetRecommendations(const NeighborSearchTypes nsType,
const InterpolationTypes interpolationType,
const size_t numRecs,
arma::Mat<size_t>& recommendations)
{
cf->GetRecommendations(nsType, interpolationType, numRecs, recommendations);
}
} // namespace cf
} // namespace mlpack
+188
View File
@@ -40,6 +40,107 @@
namespace mlpack {
namespace cf {
inline CFModel::CFModel() :
decompositionType(NMF),
normalizationType(NO_NORMALIZATION),
cf(NULL)
{
// Nothing else to do.
}
inline CFModel::CFModel(const CFModel& other) :
decompositionType(other.decompositionType),
normalizationType(other.normalizationType),
cf(other.cf->Clone())
{
// Nothing else to do.
}
inline CFModel::CFModel(CFModel&& other) :
decompositionType(other.decompositionType),
normalizationType(other.normalizationType),
cf(std::move(other.cf))
{
// Reset properties of the other one.
other.decompositionType = NMF;
other.normalizationType = NO_NORMALIZATION;
}
inline CFModel& CFModel::operator=(const CFModel& other)
{
if (this != &other)
{
decompositionType = other.decompositionType;
normalizationType = other.normalizationType;
cf = other.cf->Clone();
}
return *this;
}
inline CFModel& CFModel::operator=(CFModel&& other)
{
if (this != &other)
{
decompositionType = other.decompositionType;
normalizationType = other.normalizationType;
cf = std::move(other.cf);
// Reset the other object.
other.decompositionType = NMF;
other.normalizationType = NO_NORMALIZATION;
}
return *this;
}
inline CFModel::~CFModel()
{
delete cf;
}
template<typename DecompositionPolicy>
CFWrapperBase* TrainHelper(const DecompositionPolicy& decomposition,
const CFModel::NormalizationTypes normalizationType,
const arma::mat& data,
const size_t numUsersForSimilarity,
const size_t rank,
const size_t maxIterations,
const double minResidue,
const bool mit)
{
switch (normalizationType)
{
case CFModel::NO_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, NoNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
case CFModel::ITEM_MEAN_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, ItemMeanNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
case CFModel::USER_MEAN_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, UserMeanNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
case CFModel::OVERALL_MEAN_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, OverallMeanNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
case CFModel::Z_SCORE_NORMALIZATION:
return new CFWrapper<DecompositionPolicy, ZScoreNormalization>(data,
decomposition, numUsersForSimilarity, rank, maxIterations, minResidue,
mit);
}
// This shouldn't ever happen.
return NULL;
}
template<typename NeighborSearchPolicy, typename CFType>
void PredictHelper(CFType& cf,
const InterpolationTypes interpolationType,
@@ -321,6 +422,93 @@ void SerializeHelper(Archive& ar,
}
}
inline void CFModel::Train(
const arma::mat& data,
const size_t numUsersForSimilarity,
const size_t rank,
const size_t maxIterations,
const double minResidue,
const bool mit)
{
// Delete the current CFType object, if there is one.
delete cf;
switch (decompositionType)
{
case NMF:
cf = TrainHelper(NMFPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case BATCH_SVD:
cf = TrainHelper(BatchSVDPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case RANDOMIZED_SVD:
cf = TrainHelper(RandomizedSVDPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case REG_SVD:
cf = TrainHelper(RegSVDPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case SVD_COMPLETE:
cf = TrainHelper(SVDCompletePolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case SVD_INCOMPLETE:
cf = TrainHelper(SVDIncompletePolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case BIAS_SVD:
cf = TrainHelper(BiasSVDPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
case SVD_PLUS_PLUS:
cf = TrainHelper(SVDPlusPlusPolicy(), normalizationType, data,
numUsersForSimilarity, rank, maxIterations, minResidue, mit);
break;
}
}
//! Make predictions.
inline void CFModel::Predict(
const NeighborSearchTypes nsType,
const InterpolationTypes interpolationType,
const arma::Mat<size_t>& combinations,
arma::vec& predictions)
{
cf->Predict(nsType, interpolationType, combinations, predictions);
}
//! Compute recommendations for queried users.
inline void CFModel::GetRecommendations(
const NeighborSearchTypes nsType,
const InterpolationTypes interpolationType,
const size_t numRecs,
arma::Mat<size_t>& recommendations,
const arma::Col<size_t>& users)
{
cf->GetRecommendations(nsType, interpolationType, numRecs, recommendations,
users);
}
//! Compute recommendations for all users.
inline void CFModel::GetRecommendations(
const NeighborSearchTypes nsType,
const InterpolationTypes interpolationType,
const size_t numRecs,
arma::Mat<size_t>& recommendations)
{
cf->GetRecommendations(nsType, interpolationType, numRecs, recommendations);
}
template<typename Archive>
void CFModel::serialize(Archive& ar, const uint32_t /* version */)
{
@@ -5,7 +5,6 @@ set(SOURCES
fastmks_impl.hpp
fastmks_model.hpp
fastmks_model_impl.hpp
fastmks_model.cpp
fastmks_rules.hpp
fastmks_rules_impl.hpp
fastmks_stat.hpp
@@ -1,318 +0,0 @@
/**
* @file methods/fastmks/fastmks_model.cpp
* @author Ryan Curtin
*
* Implementation of non-templatized functions of FastMKSModel.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "fastmks_model.hpp"
using namespace mlpack;
using namespace mlpack::fastmks;
using namespace mlpack::kernel;
FastMKSModel::FastMKSModel(const int kernelType) :
kernelType(kernelType),
linear(NULL),
polynomial(NULL),
cosine(NULL),
gaussian(NULL),
epan(NULL),
triangular(NULL),
hyptan(NULL)
{
// Nothing to do.
}
FastMKSModel::FastMKSModel(const FastMKSModel& other) :
kernelType(other.kernelType),
linear(other.linear == NULL ? NULL :
new FastMKS<LinearKernel>(*other.linear)),
polynomial(other.polynomial == NULL ? NULL :
new FastMKS<PolynomialKernel>(*other.polynomial)),
cosine(other.cosine == NULL ? NULL :
new FastMKS<CosineDistance>(*other.cosine)),
gaussian(other.gaussian == NULL ? NULL :
new FastMKS<GaussianKernel>(*other.gaussian)),
epan(other.epan == NULL ? NULL :
new FastMKS<EpanechnikovKernel>(*other.epan)),
triangular(other.triangular == NULL ? NULL :
new FastMKS<TriangularKernel>(*other.triangular)),
hyptan(other.hyptan == NULL ? NULL :
new FastMKS<HyperbolicTangentKernel>(*other.hyptan))
{
// Nothing to do.
}
FastMKSModel::FastMKSModel(FastMKSModel&& other) :
kernelType(other.kernelType),
linear(other.linear),
polynomial(other.polynomial),
cosine(other.cosine),
gaussian(other.gaussian),
epan(other.epan),
triangular(other.triangular),
hyptan(other.hyptan)
{
// Clear other object.
other.kernelType = KernelTypes::LINEAR_KERNEL;
other.linear = NULL;
other.polynomial = NULL;
other.cosine = NULL;
other.gaussian = NULL;
other.epan = NULL;
other.triangular = NULL;
other.hyptan = NULL;
}
FastMKSModel& FastMKSModel::operator=(const FastMKSModel& other)
{
if (this != &other)
{
// Clear memory.
delete linear;
delete polynomial;
delete cosine;
delete gaussian;
delete epan;
delete triangular;
delete hyptan;
// Set pointers to null.
linear = NULL;
polynomial = NULL;
cosine = NULL;
gaussian = NULL;
epan = NULL;
triangular = NULL;
hyptan = NULL;
kernelType = other.kernelType;
if (other.linear)
linear = new FastMKS<LinearKernel>(*other.linear);
if (other.polynomial)
polynomial = new FastMKS<PolynomialKernel>(*other.polynomial);
if (other.cosine)
cosine = new FastMKS<CosineDistance>(*other.cosine);
if (other.gaussian)
gaussian = new FastMKS<GaussianKernel>(*other.gaussian);
if (other.epan)
epan = new FastMKS<EpanechnikovKernel>(*other.epan);
if (other.triangular)
triangular = new FastMKS<TriangularKernel>(*other.triangular);
if (other.hyptan)
hyptan = new FastMKS<HyperbolicTangentKernel>(*other.hyptan);
}
return *this;
}
FastMKSModel& FastMKSModel::operator=(FastMKSModel&& other)
{
if (this != &other)
{
kernelType = other.kernelType;
linear = other.linear;
polynomial = other.polynomial;
cosine = other.cosine;
gaussian = other.gaussian;
epan = other.epan;
triangular = other.triangular;
hyptan = other.hyptan;
// Clear other object.
other.kernelType = KernelTypes::LINEAR_KERNEL;
other.linear = nullptr;
other.polynomial = nullptr;
other.cosine = nullptr;
other.gaussian = nullptr;
other.epan = nullptr;
other.triangular = nullptr;
other.hyptan = nullptr;
}
return *this;
}
FastMKSModel::~FastMKSModel()
{
// Clean memory.
if (linear)
delete linear;
if (polynomial)
delete polynomial;
if (cosine)
delete cosine;
if (gaussian)
delete gaussian;
if (epan)
delete epan;
if (triangular)
delete triangular;
if (hyptan)
delete hyptan;
}
bool FastMKSModel::Naive() const
{
switch (kernelType)
{
case LINEAR_KERNEL:
return linear->Naive();
case POLYNOMIAL_KERNEL:
return polynomial->Naive();
case COSINE_DISTANCE:
return cosine->Naive();
case GAUSSIAN_KERNEL:
return gaussian->Naive();
case EPANECHNIKOV_KERNEL:
return epan->Naive();
case TRIANGULAR_KERNEL:
return triangular->Naive();
case HYPTAN_KERNEL:
return hyptan->Naive();
}
throw std::runtime_error("invalid model type");
}
bool& FastMKSModel::Naive()
{
switch (kernelType)
{
case LINEAR_KERNEL:
return linear->Naive();
case POLYNOMIAL_KERNEL:
return polynomial->Naive();
case COSINE_DISTANCE:
return cosine->Naive();
case GAUSSIAN_KERNEL:
return gaussian->Naive();
case EPANECHNIKOV_KERNEL:
return epan->Naive();
case TRIANGULAR_KERNEL:
return triangular->Naive();
case HYPTAN_KERNEL:
return hyptan->Naive();
}
throw std::runtime_error("invalid model type");
}
bool FastMKSModel::SingleMode() const
{
switch (kernelType)
{
case LINEAR_KERNEL:
return linear->SingleMode();
case POLYNOMIAL_KERNEL:
return polynomial->SingleMode();
case COSINE_DISTANCE:
return cosine->SingleMode();
case GAUSSIAN_KERNEL:
return gaussian->SingleMode();
case EPANECHNIKOV_KERNEL:
return epan->SingleMode();
case TRIANGULAR_KERNEL:
return triangular->SingleMode();
case HYPTAN_KERNEL:
return hyptan->SingleMode();
}
throw std::runtime_error("invalid model type");
}
bool& FastMKSModel::SingleMode()
{
switch (kernelType)
{
case LINEAR_KERNEL:
return linear->SingleMode();
case POLYNOMIAL_KERNEL:
return polynomial->SingleMode();
case COSINE_DISTANCE:
return cosine->SingleMode();
case GAUSSIAN_KERNEL:
return gaussian->SingleMode();
case EPANECHNIKOV_KERNEL:
return epan->SingleMode();
case TRIANGULAR_KERNEL:
return triangular->SingleMode();
case HYPTAN_KERNEL:
return hyptan->SingleMode();
}
throw std::runtime_error("invalid model type");
}
void FastMKSModel::Search(util::Timers& timers,
const arma::mat& querySet,
const size_t k,
arma::Mat<size_t>& indices,
arma::mat& kernels,
const double base)
{
switch (kernelType)
{
case LINEAR_KERNEL:
Search(timers, *linear, querySet, k, indices, kernels, base);
break;
case POLYNOMIAL_KERNEL:
Search(timers, *polynomial, querySet, k, indices, kernels, base);
break;
case COSINE_DISTANCE:
Search(timers, *cosine, querySet, k, indices, kernels, base);
break;
case GAUSSIAN_KERNEL:
Search(timers, *gaussian, querySet, k, indices, kernels, base);
break;
case EPANECHNIKOV_KERNEL:
Search(timers, *epan, querySet, k, indices, kernels, base);
break;
case TRIANGULAR_KERNEL:
Search(timers, *triangular, querySet, k, indices, kernels, base);
break;
case HYPTAN_KERNEL:
Search(timers, *hyptan, querySet, k, indices, kernels, base);
break;
default:
throw std::runtime_error("invalid model type");
}
}
void FastMKSModel::Search(util::Timers& timers,
const size_t k,
arma::Mat<size_t>& indices,
arma::mat& kernels)
{
timers.Start("computing_products");
switch (kernelType)
{
case LINEAR_KERNEL:
linear->Search(k, indices, kernels);
break;
case POLYNOMIAL_KERNEL:
polynomial->Search(k, indices, kernels);
break;
case COSINE_DISTANCE:
cosine->Search(k, indices, kernels);
break;
case GAUSSIAN_KERNEL:
gaussian->Search(k, indices, kernels);
break;
case EPANECHNIKOV_KERNEL:
epan->Search(k, indices, kernels);
break;
case TRIANGULAR_KERNEL:
triangular->Search(k, indices, kernels);
break;
case HYPTAN_KERNEL:
hyptan->Search(k, indices, kernels);
break;
default:
throw std::invalid_argument("invalid model type");
}
timers.Stop("computing_products");
}
@@ -17,6 +17,238 @@
namespace mlpack {
namespace fastmks {
inline FastMKSModel::FastMKSModel(const int kernelType) :
kernelType(kernelType),
linear(NULL),
polynomial(NULL),
cosine(NULL),
gaussian(NULL),
epan(NULL),
triangular(NULL),
hyptan(NULL)
{
// Nothing to do.
}
inline FastMKSModel::FastMKSModel(const FastMKSModel& other) :
kernelType(other.kernelType),
linear(other.linear == NULL ? NULL :
new FastMKS<kernel::LinearKernel>(*other.linear)),
polynomial(other.polynomial == NULL ? NULL :
new FastMKS<kernel::PolynomialKernel>(*other.polynomial)),
cosine(other.cosine == NULL ? NULL :
new FastMKS<kernel::CosineDistance>(*other.cosine)),
gaussian(other.gaussian == NULL ? NULL :
new FastMKS<kernel::GaussianKernel>(*other.gaussian)),
epan(other.epan == NULL ? NULL :
new FastMKS<kernel::EpanechnikovKernel>(*other.epan)),
triangular(other.triangular == NULL ? NULL :
new FastMKS<kernel::TriangularKernel>(*other.triangular)),
hyptan(other.hyptan == NULL ? NULL :
new FastMKS<kernel::HyperbolicTangentKernel>(*other.hyptan))
{
// Nothing to do.
}
inline FastMKSModel::FastMKSModel(FastMKSModel&& other) :
kernelType(other.kernelType),
linear(other.linear),
polynomial(other.polynomial),
cosine(other.cosine),
gaussian(other.gaussian),
epan(other.epan),
triangular(other.triangular),
hyptan(other.hyptan)
{
// Clear other object.
other.kernelType = KernelTypes::LINEAR_KERNEL;
other.linear = NULL;
other.polynomial = NULL;
other.cosine = NULL;
other.gaussian = NULL;
other.epan = NULL;
other.triangular = NULL;
other.hyptan = NULL;
}
inline FastMKSModel& FastMKSModel::operator=(const FastMKSModel& other)
{
if (this != &other)
{
// Clear memory.
delete linear;
delete polynomial;
delete cosine;
delete gaussian;
delete epan;
delete triangular;
delete hyptan;
// Set pointers to null.
linear = NULL;
polynomial = NULL;
cosine = NULL;
gaussian = NULL;
epan = NULL;
triangular = NULL;
hyptan = NULL;
kernelType = other.kernelType;
if (other.linear)
linear = new FastMKS<kernel::LinearKernel>(*other.linear);
if (other.polynomial)
polynomial = new FastMKS<kernel::PolynomialKernel>(*other.polynomial);
if (other.cosine)
cosine = new FastMKS<kernel::CosineDistance>(*other.cosine);
if (other.gaussian)
gaussian = new FastMKS<kernel::GaussianKernel>(*other.gaussian);
if (other.epan)
epan = new FastMKS<kernel::EpanechnikovKernel>(*other.epan);
if (other.triangular)
triangular = new FastMKS<kernel::TriangularKernel>(*other.triangular);
if (other.hyptan)
hyptan = new FastMKS<kernel::HyperbolicTangentKernel>(*other.hyptan);
}
return *this;
}
inline FastMKSModel& FastMKSModel::operator=(FastMKSModel&& other)
{
if (this != &other)
{
kernelType = other.kernelType;
linear = other.linear;
polynomial = other.polynomial;
cosine = other.cosine;
gaussian = other.gaussian;
epan = other.epan;
triangular = other.triangular;
hyptan = other.hyptan;
// Clear other object.
other.kernelType = KernelTypes::LINEAR_KERNEL;
other.linear = nullptr;
other.polynomial = nullptr;
other.cosine = nullptr;
other.gaussian = nullptr;
other.epan = nullptr;
other.triangular = nullptr;
other.hyptan = nullptr;
}
return *this;
}
inline FastMKSModel::~FastMKSModel()
{
// Clean memory.
if (linear)
delete linear;
if (polynomial)
delete polynomial;
if (cosine)
delete cosine;
if (gaussian)
delete gaussian;
if (epan)
delete epan;
if (triangular)
delete triangular;
if (hyptan)
delete hyptan;
}
inline bool FastMKSModel::Naive() const
{
switch (kernelType)
{
case LINEAR_KERNEL:
return linear->Naive();
case POLYNOMIAL_KERNEL:
return polynomial->Naive();
case COSINE_DISTANCE:
return cosine->Naive();
case GAUSSIAN_KERNEL:
return gaussian->Naive();
case EPANECHNIKOV_KERNEL:
return epan->Naive();
case TRIANGULAR_KERNEL:
return triangular->Naive();
case HYPTAN_KERNEL:
return hyptan->Naive();
}
throw std::runtime_error("invalid model type");
}
inline bool& FastMKSModel::Naive()
{
switch (kernelType)
{
case LINEAR_KERNEL:
return linear->Naive();
case POLYNOMIAL_KERNEL:
return polynomial->Naive();
case COSINE_DISTANCE:
return cosine->Naive();
case GAUSSIAN_KERNEL:
return gaussian->Naive();
case EPANECHNIKOV_KERNEL:
return epan->Naive();
case TRIANGULAR_KERNEL:
return triangular->Naive();
case HYPTAN_KERNEL:
return hyptan->Naive();
}
throw std::runtime_error("invalid model type");
}
inline bool FastMKSModel::SingleMode() const
{
switch (kernelType)
{
case LINEAR_KERNEL:
return linear->SingleMode();
case POLYNOMIAL_KERNEL:
return polynomial->SingleMode();
case COSINE_DISTANCE:
return cosine->SingleMode();
case GAUSSIAN_KERNEL:
return gaussian->SingleMode();
case EPANECHNIKOV_KERNEL:
return epan->SingleMode();
case TRIANGULAR_KERNEL:
return triangular->SingleMode();
case HYPTAN_KERNEL:
return hyptan->SingleMode();
}
throw std::runtime_error("invalid model type");
}
inline bool& FastMKSModel::SingleMode()
{
switch (kernelType)
{
case LINEAR_KERNEL:
return linear->SingleMode();
case POLYNOMIAL_KERNEL:
return polynomial->SingleMode();
case COSINE_DISTANCE:
return cosine->SingleMode();
case GAUSSIAN_KERNEL:
return gaussian->SingleMode();
case EPANECHNIKOV_KERNEL:
return epan->SingleMode();
case TRIANGULAR_KERNEL:
return triangular->SingleMode();
case HYPTAN_KERNEL:
return hyptan->SingleMode();
}
throw std::runtime_error("invalid model type");
}
//! This is called when the KernelType is the same as the model.
template<typename KernelType>
void BuildFastMKSModel(util::Timers& timers,
@@ -229,6 +461,78 @@ void FastMKSModel::Search(util::Timers& timers,
}
}
inline void FastMKSModel::Search(
util::Timers& timers,
const arma::mat& querySet,
const size_t k,
arma::Mat<size_t>& indices,
arma::mat& kernels,
const double base)
{
switch (kernelType)
{
case LINEAR_KERNEL:
Search(timers, *linear, querySet, k, indices, kernels, base);
break;
case POLYNOMIAL_KERNEL:
Search(timers, *polynomial, querySet, k, indices, kernels, base);
break;
case COSINE_DISTANCE:
Search(timers, *cosine, querySet, k, indices, kernels, base);
break;
case GAUSSIAN_KERNEL:
Search(timers, *gaussian, querySet, k, indices, kernels, base);
break;
case EPANECHNIKOV_KERNEL:
Search(timers, *epan, querySet, k, indices, kernels, base);
break;
case TRIANGULAR_KERNEL:
Search(timers, *triangular, querySet, k, indices, kernels, base);
break;
case HYPTAN_KERNEL:
Search(timers, *hyptan, querySet, k, indices, kernels, base);
break;
default:
throw std::runtime_error("invalid model type");
}
}
inline void FastMKSModel::Search(
util::Timers& timers,
const size_t k,
arma::Mat<size_t>& indices,
arma::mat& kernels)
{
timers.Start("computing_products");
switch (kernelType)
{
case LINEAR_KERNEL:
linear->Search(k, indices, kernels);
break;
case POLYNOMIAL_KERNEL:
polynomial->Search(k, indices, kernels);
break;
case COSINE_DISTANCE:
cosine->Search(k, indices, kernels);
break;
case GAUSSIAN_KERNEL:
gaussian->Search(k, indices, kernels);
break;
case EPANECHNIKOV_KERNEL:
epan->Search(k, indices, kernels);
break;
case TRIANGULAR_KERNEL:
triangular->Search(k, indices, kernels);
break;
case HYPTAN_KERNEL:
hyptan->Search(k, indices, kernels);
break;
default:
throw std::invalid_argument("invalid model type");
}
timers.Stop("computing_products");
}
} // namespace fastmks
} // namespace mlpack
-2
View File
@@ -2,10 +2,8 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
gmm.hpp
gmm.cpp
gmm_impl.hpp
diagonal_gmm.hpp
diagonal_gmm.cpp
diagonal_gmm_impl.hpp
em_fit.hpp
em_fit_impl.hpp
-230
View File
@@ -1,230 +0,0 @@
/**
* @file methods/gmm/diagonal_gmm.cpp
* @author Kim SangYeon
*
* Implementation of template-based GMM methods.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "diagonal_gmm.hpp"
#include <mlpack/core/math/log_add.hpp>
namespace mlpack {
namespace gmm {
/**
* Create a DiagonalGMM with the given number of Gaussians, each of which have
* the specified dimensionality. The means and covariances will be set to 0.
*
* @param gaussians Number of Gaussians in this GMM.
* @param dimensionality Dimensionality of each Gaussian.
*/
DiagonalGMM::DiagonalGMM(const size_t gaussians, const size_t dimensionality) :
gaussians(gaussians),
dimensionality(dimensionality),
dists(gaussians,
distribution::DiagonalGaussianDistribution(dimensionality)),
weights(gaussians)
{
// Set equal weights. Technically this model is still valid, but only barely.
weights.fill(1.0 / gaussians);
}
// Copy constructor for when the other GMM uses the same fitting type.
DiagonalGMM::DiagonalGMM(const DiagonalGMM& other) :
gaussians(other.Gaussians()),
dimensionality(other.dimensionality),
dists(other.dists),
weights(other.weights) { /* Nothing to do. */ }
DiagonalGMM& DiagonalGMM::operator=(const DiagonalGMM& other)
{
gaussians = other.gaussians;
dimensionality = other.dimensionality;
dists = other.dists;
weights = other.weights;
return *this;
}
/**
* Return the log probability of the given observation being from this GMM.
*/
double DiagonalGMM::LogProbability(const arma::vec& observation) const
{
// Sum the probability for each Gaussian in our mixture (and we have to
// multiply by the prior for each Gaussian too).
double sum = -std::numeric_limits<double>::infinity();
for (size_t i = 0; i < gaussians; ++i)
{
sum = math::LogAdd(sum, log(weights[i]) +
dists[i].LogProbability(observation));
}
return sum;
}
/**
* Return the log probability of the given observation GMM matrix.
*
* @param observation Observation matrix to compute log-probabilty.
* @param logProbs Stores the value of log-probability for input.
*/
void DiagonalGMM::LogProbability(const arma::mat& observation,
arma::vec& logProbs) const
{
// Sum the probability for each Gaussian in our mixture (and we have to
// multiply by the prior for each Gaussian too).
logProbs.set_size(observation.n_cols);
// Store log-probability value in a matrix.
arma::mat logProb(observation.n_cols, gaussians);
// Assign value to the matrix.
for (size_t i = 0; i < gaussians; i++)
{
arma::vec temp(logProb.colptr(i), observation.n_cols, false, true);
dists[i].LogProbability(observation, temp);
}
// Save log(weights) as a vector.
arma::vec logWeights = arma::log(weights);
// Compute log-probability.
logProb += repmat(logWeights.t(), logProb.n_rows, 1);
math::LogSumExp(logProb, logProbs);
}
/**
* Return the probability of the given observation being from this GMM.
*/
double DiagonalGMM::Probability(const arma::vec& observation) const
{
return exp(LogProbability(observation));
}
/**
* Return the probability of the given observation GMM matrix.
*
* @param observation Observation matrix to compute probabilty.
* @param probs Stores the value of probability for observation.
*/
void DiagonalGMM::Probability(const arma::mat& observation,
arma::vec& probs) const
{
LogProbability(observation, probs);
probs = exp(probs);
}
/**
* Return the log probability of the given observation being from the given
* component in the mixture.
*/
double DiagonalGMM::LogProbability(const arma::vec& observation,
const size_t component) const
{
// We are only considering one Gaussian component -- so we only need to call
// Probability() once. We do consider the prior probability!
return log(weights[component]) +
dists[component].LogProbability(observation);
}
/**
* Return the probability of the given observation being from the given
* component in the mixture.
*/
double DiagonalGMM::Probability(const arma::vec& observation,
const size_t component) const
{
return exp(LogProbability(observation, component));
}
/**
* Return a randomly generated observation according to the probability
* distribution defined by this object.
*/
arma::vec DiagonalGMM::Random() const
{
// Determine which Gaussian it will be coming from.
double gaussRand = math::Random();
size_t gaussian = 0;
double sumProb = 0;
for (size_t g = 0; g < gaussians; g++)
{
sumProb += weights(g);
if (gaussRand <= sumProb)
{
gaussian = g;
break;
}
}
return arma::sqrt(dists[gaussian].Covariance()) %
arma::randn<arma::vec>(dimensionality) + dists[gaussian].Mean();
}
/**
* Classify the given observations as being from an individual component in
* this GMM.
*/
void DiagonalGMM::Classify(const arma::mat& observations,
arma::Row<size_t>& labels) const
{
// This is not the best way to do this!
// We should not have to fill this with values, because each one should be
// overwritten.
labels.set_size(observations.n_cols);
for (size_t i = 0; i < observations.n_cols; ++i)
{
// Find maximum probability component.
double probability = 0;
for (size_t j = 0; j < gaussians; ++j)
{
double newProb = Probability(observations.unsafe_col(i), j);
if (newProb >= probability)
{
probability = newProb;
labels[i] = j;
}
}
}
}
/**
* Get the log-likelihood of this data's fit to the model.
*/
double DiagonalGMM::LogLikelihood(
const arma::mat& observations,
const std::vector<distribution::DiagonalGaussianDistribution>& dists,
const arma::vec& weights) const
{
double logLikelihood = 0;
arma::vec phis;
arma::mat likelihoods(gaussians, observations.n_cols);
for (size_t i = 0; i < gaussians; ++i)
{
dists[i].Probability(observations, phis);
likelihoods.row(i) = weights(i) * trans(phis);
}
// Now sum over every point.
for (size_t j = 0; j < observations.n_cols; ++j)
{
if (accu(likelihoods.col(j)) == 0)
Log::Info << "Likelihood of point " << j << " is 0! It is probably an "
<< "outlier." << std::endl;
logLikelihood += log(accu(likelihoods.col(j)));
}
return logLikelihood;
}
} // namespace gmm
} // namespace mlpack
+1
View File
@@ -16,6 +16,7 @@
#include <mlpack/prereqs.hpp>
#include <mlpack/core/dists/diagonal_gaussian_distribution.hpp>
#include <mlpack/core/math/log_add.hpp>
// This is the default fitting method class.
#include "em_fit.hpp"
@@ -101,6 +101,218 @@ double DiagonalGMM::Train(const arma::mat& observations,
return bestLikelihood;
}
/**
* Create a DiagonalGMM with the given number of Gaussians, each of which have
* the specified dimensionality. The means and covariances will be set to 0.
*
* @param gaussians Number of Gaussians in this GMM.
* @param dimensionality Dimensionality of each Gaussian.
*/
inline DiagonalGMM::DiagonalGMM(
const size_t gaussians,
const size_t dimensionality) :
gaussians(gaussians),
dimensionality(dimensionality),
dists(gaussians,
distribution::DiagonalGaussianDistribution(dimensionality)),
weights(gaussians)
{
// Set equal weights. Technically this model is still valid, but only barely.
weights.fill(1.0 / gaussians);
}
// Copy constructor for when the other GMM uses the same fitting type.
inline DiagonalGMM::DiagonalGMM(const DiagonalGMM& other) :
gaussians(other.Gaussians()),
dimensionality(other.dimensionality),
dists(other.dists),
weights(other.weights) { /* Nothing to do. */ }
inline DiagonalGMM& DiagonalGMM::operator=(const DiagonalGMM& other)
{
gaussians = other.gaussians;
dimensionality = other.dimensionality;
dists = other.dists;
weights = other.weights;
return *this;
}
/**
* Return the log probability of the given observation being from this GMM.
*/
inline double DiagonalGMM::LogProbability(const arma::vec& observation) const
{
// Sum the probability for each Gaussian in our mixture (and we have to
// multiply by the prior for each Gaussian too).
double sum = -std::numeric_limits<double>::infinity();
for (size_t i = 0; i < gaussians; ++i)
{
sum = math::LogAdd(sum, log(weights[i]) +
dists[i].LogProbability(observation));
}
return sum;
}
/**
* Return the log probability of the given observation GMM matrix.
*
* @param observation Observation matrix to compute log-probabilty.
* @param logProbs Stores the value of log-probability for input.
*/
inline void DiagonalGMM::LogProbability(const arma::mat& observation,
arma::vec& logProbs) const
{
// Sum the probability for each Gaussian in our mixture (and we have to
// multiply by the prior for each Gaussian too).
logProbs.set_size(observation.n_cols);
// Store log-probability value in a matrix.
arma::mat logProb(observation.n_cols, gaussians);
// Assign value to the matrix.
for (size_t i = 0; i < gaussians; i++)
{
arma::vec temp(logProb.colptr(i), observation.n_cols, false, true);
dists[i].LogProbability(observation, temp);
}
// Save log(weights) as a vector.
arma::vec logWeights = arma::log(weights);
// Compute log-probability.
logProb += repmat(logWeights.t(), logProb.n_rows, 1);
math::LogSumExp(logProb, logProbs);
}
/**
* Return the probability of the given observation being from this GMM.
*/
inline double DiagonalGMM::Probability(const arma::vec& observation) const
{
return exp(LogProbability(observation));
}
/**
* Return the probability of the given observation GMM matrix.
*
* @param observation Observation matrix to compute probabilty.
* @param probs Stores the value of probability for observation.
*/
inline void DiagonalGMM::Probability(const arma::mat& observation,
arma::vec& probs) const
{
LogProbability(observation, probs);
probs = exp(probs);
}
/**
* Return the log probability of the given observation being from the given
* component in the mixture.
*/
inline double DiagonalGMM::LogProbability(const arma::vec& observation,
const size_t component) const
{
// We are only considering one Gaussian component -- so we only need to call
// Probability() once. We do consider the prior probability!
return log(weights[component]) +
dists[component].LogProbability(observation);
}
/**
* Return the probability of the given observation being from the given
* component in the mixture.
*/
inline double DiagonalGMM::Probability(const arma::vec& observation,
const size_t component) const
{
return exp(LogProbability(observation, component));
}
/**
* Return a randomly generated observation according to the probability
* distribution defined by this object.
*/
inline arma::vec DiagonalGMM::Random() const
{
// Determine which Gaussian it will be coming from.
double gaussRand = math::Random();
size_t gaussian = 0;
double sumProb = 0;
for (size_t g = 0; g < gaussians; g++)
{
sumProb += weights(g);
if (gaussRand <= sumProb)
{
gaussian = g;
break;
}
}
return arma::sqrt(dists[gaussian].Covariance()) %
arma::randn<arma::vec>(dimensionality) + dists[gaussian].Mean();
}
/**
* Classify the given observations as being from an individual component in
* this GMM.
*/
inline void DiagonalGMM::Classify(const arma::mat& observations,
arma::Row<size_t>& labels) const
{
// This is not the best way to do this!
// We should not have to fill this with values, because each one should be
// overwritten.
labels.set_size(observations.n_cols);
for (size_t i = 0; i < observations.n_cols; ++i)
{
// Find maximum probability component.
double probability = 0;
for (size_t j = 0; j < gaussians; ++j)
{
double newProb = Probability(observations.unsafe_col(i), j);
if (newProb >= probability)
{
probability = newProb;
labels[i] = j;
}
}
}
}
/**
* Get the log-likelihood of this data's fit to the model.
*/
inline double DiagonalGMM::LogLikelihood(
const arma::mat& observations,
const std::vector<distribution::DiagonalGaussianDistribution>& dists,
const arma::vec& weights) const
{
double logLikelihood = 0;
arma::vec phis;
arma::mat likelihoods(gaussians, observations.n_cols);
for (size_t i = 0; i < gaussians; ++i)
{
dists[i].Probability(observations, phis);
likelihoods.row(i) = weights(i) * trans(phis);
}
// Now sum over every point.
for (size_t j = 0; j < observations.n_cols; ++j)
{
if (accu(likelihoods.col(j)) == 0)
Log::Info << "Likelihood of point " << j << " is 0! It is probably an "
<< "outlier." << std::endl;
logLikelihood += log(accu(likelihoods.col(j)));
}
return logLikelihood;
}
/**
* Fit the DiagonalGMM to the given observations, each of which has a certain
* probability of being from this distribution.
-247
View File
@@ -1,247 +0,0 @@
/**
* @file methods/gmm/gmm.cpp
* @author Parikshit Ram (pram@cc.gatech.edu)
* @author Ryan Curtin
* @author Michael Fox
*
* Implementation of template-based GMM methods.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "gmm.hpp"
#include <mlpack/core/math/log_add.hpp>
namespace mlpack {
namespace gmm {
/**
* Create a GMM with the given number of Gaussians, each of which have the
* specified dimensionality. The means and covariances will be set to 0.
*
* @param gaussians Number of Gaussians in this GMM.
* @param dimensionality Dimensionality of each Gaussian.
*/
GMM::GMM(const size_t gaussians, const size_t dimensionality) :
gaussians(gaussians),
dimensionality(dimensionality),
dists(gaussians, distribution::GaussianDistribution(dimensionality)),
weights(gaussians)
{
// Set equal weights. Technically this model is still valid, but only barely.
weights.fill(1.0 / gaussians);
}
// Copy constructor for when the other GMM uses the same fitting type.
GMM::GMM(const GMM& other) :
gaussians(other.Gaussians()),
dimensionality(other.dimensionality),
dists(other.dists),
weights(other.weights) { /* Nothing to do. */ }
GMM& GMM::operator=(const GMM& other)
{
gaussians = other.gaussians;
dimensionality = other.dimensionality;
dists = other.dists;
weights = other.weights;
return *this;
}
/**
* Return the log probability of the given observation being from this GMM.
*
* @param observation Observation vector to compute log-probabilty.
*/
double GMM::LogProbability(const arma::vec& observation) const
{
// Sum the probability for each Gaussian in our mixture (and we have to
// multiply by the prior for each Gaussian too).
double sum = -std::numeric_limits<double>::infinity();
for (size_t i = 0; i < gaussians; ++i)
sum = math::LogAdd(sum, log(weights[i]) +
dists[i].LogProbability(observation));
return sum;
}
/**
* Return the log probability of the given observation GMM matrix.
*
* @param observation Observation matrix to compute log-probabilty.
* @param logProbs Stores the value of log-probability for Observation.
*/
void GMM::LogProbability(const arma::mat& observation,
arma::vec& logProbs) const
{
// Sum the probability for each Gaussian in our mixture (and we have to
// multiply by the prior for each Gaussian too).
logProbs.set_size(observation.n_cols);
// Store log-probability value in a matrix.
arma::mat logProb(observation.n_cols, gaussians);
// Assign value to the matrix.
for (size_t i = 0; i < gaussians; i++)
{
arma::vec temp(logProb.colptr(i), observation.n_cols, false, true);
dists[i].LogProbability(observation, temp);
}
// Save log(weights) as a vector.
arma::vec logWeights = arma::log(weights);
// Compute log-probability.
logProb += repmat(logWeights.t(), logProb.n_rows, 1);
math::LogSumExp(logProb, logProbs);
}
/**
* Return the probability of the given observation being from this GMM.
*
* @param observation Observation vector to compute probabilty.
*/
double GMM::Probability(const arma::vec& observation) const
{
return exp(LogProbability(observation));
}
/**
* Return the probability of the given observation GMM matrix.
*
* @param observation Observation matrix to compute probabilty.
* @param probs Stores the value of probability for x.
*/
void GMM::Probability(const arma::mat& observation,
arma::vec& probs) const
{
LogProbability(observation, probs);
probs = exp(probs);
}
/**
* Return the log probability of the given observation being from the given
* component in the mixture.
*
* @param observation Observation vector to compute log-probabilty.
* @param component Calculate the log-probability for given observation vector.
*/
double GMM::LogProbability(const arma::vec& observation,
const size_t component) const
{
// We are only considering one Gaussian component -- so we only need to call
// Probability() once. We do consider the prior probability!
return log(weights[component]) + dists[component].LogProbability(observation);
}
/**
* Return the probability of the given observation being from the given
* component in the mixture.
*
* @param observation Observation matrix to compute probabilty.
* @param component Calculate the probability for given component.
*/
double GMM::Probability(const arma::vec& observation,
const size_t component) const
{
return exp(LogProbability(observation, component));
}
/**
* Return a randomly generated observation according to the probability
* distribution defined by this object.
*/
arma::vec GMM::Random() const
{
// Determine which Gaussian it will be coming from.
double gaussRand = math::Random();
size_t gaussian = 0;
double sumProb = 0;
for (size_t g = 0; g < gaussians; g++)
{
sumProb += weights(g);
if (gaussRand <= sumProb)
{
gaussian = g;
break;
}
}
arma::mat cholDecomp;
if (!arma::chol(cholDecomp, dists[gaussian].Covariance()))
{
Log::Fatal << "Cholesky decomposition failed." << std::endl;
}
return trans(cholDecomp) *
arma::randn<arma::vec>(dimensionality) + dists[gaussian].Mean();
}
/**
* Classify the given observations as being from an individual component in this
* GMM.
*
* @param observation Observation matrix for classification.
* @param labels Save the labels for the given observation matrix.
*/
void GMM::Classify(const arma::mat& observations,
arma::Row<size_t>& labels) const
{
// This is not the best way to do this!
// We should not have to fill this with values, because each one should be
// overwritten.
labels.set_size(observations.n_cols);
for (size_t i = 0; i < observations.n_cols; ++i)
{
// Find maximum probability component.
double probability = -std::numeric_limits<double>::infinity();
for (size_t j = 0; j < gaussians; ++j)
{
// We have to use LogProbability() otherwise Probability() would overflow
// easily.
double newProb = LogProbability(observations.unsafe_col(i), j);
if (newProb >= probability)
{
probability = newProb;
labels[i] = j;
}
}
}
}
/**
* Get the log-likelihood of this data's fit to the model.
*
* @param data Data matrix to compute log-likelihood.
* @parma distsL Vector of Gaussian distribution.
* @param weightsL Vector of weights for computing likelihoods.
*/
double GMM::LogLikelihood(
const arma::mat& data,
const std::vector<distribution::GaussianDistribution>& distsL,
const arma::vec& weightsL) const
{
double loglikelihood = 0;
arma::vec logPhis;
arma::mat logLikelihoods(gaussians, data.n_cols);
// It has to be LogProbability() otherwise Probability() would overflow easily
for (size_t i = 0; i < gaussians; ++i)
{
distsL[i].LogProbability(data, logPhis);
logLikelihoods.row(i) = log(weightsL(i)) + trans(logPhis);
}
// Now sum over every point.
for (size_t j = 0; j < data.n_cols; ++j)
loglikelihood += mlpack::math::AccuLog(logLikelihoods.col(j));
return loglikelihood;
}
} // namespace gmm
} // namespace mlpack
+4 -2
View File
@@ -10,14 +10,16 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_MOG_MOG_EM_HPP
#define MLPACK_METHODS_MOG_MOG_EM_HPP
#ifndef MLPACK_METHODS_GMM_GMM_HPP
#define MLPACK_METHODS_GMM_GMM_HPP
#include <mlpack/prereqs.hpp>
// This is the default fitting method class.
#include "em_fit.hpp"
#include <mlpack/core/math/log_add.hpp>
namespace mlpack {
namespace gmm /** Gaussian Mixture Models. */ {
+226
View File
@@ -20,6 +20,232 @@
namespace mlpack {
namespace gmm {
/**
* Create a GMM with the given number of Gaussians, each of which have the
* specified dimensionality. The means and covariances will be set to 0.
*
* @param gaussians Number of Gaussians in this GMM.
* @param dimensionality Dimensionality of each Gaussian.
*/
inline GMM::GMM(const size_t gaussians, const size_t dimensionality) :
gaussians(gaussians),
dimensionality(dimensionality),
dists(gaussians, distribution::GaussianDistribution(dimensionality)),
weights(gaussians)
{
// Set equal weights. Technically this model is still valid, but only barely.
weights.fill(1.0 / gaussians);
}
// Copy constructor for when the other GMM uses the same fitting type.
inline GMM::GMM(const GMM& other) :
gaussians(other.Gaussians()),
dimensionality(other.dimensionality),
dists(other.dists),
weights(other.weights) { /* Nothing to do. */ }
inline GMM& GMM::operator=(const GMM& other)
{
gaussians = other.gaussians;
dimensionality = other.dimensionality;
dists = other.dists;
weights = other.weights;
return *this;
}
/**
* Return the log probability of the given observation being from this GMM.
*
* @param observation Observation vector to compute log-probabilty.
*/
inline double GMM::LogProbability(const arma::vec& observation) const
{
// Sum the probability for each Gaussian in our mixture (and we have to
// multiply by the prior for each Gaussian too).
double sum = -std::numeric_limits<double>::infinity();
for (size_t i = 0; i < gaussians; ++i)
sum = math::LogAdd(sum, log(weights[i]) +
dists[i].LogProbability(observation));
return sum;
}
/**
* Return the log probability of the given observation GMM matrix.
*
* @param observation Observation matrix to compute log-probabilty.
* @param logProbs Stores the value of log-probability for Observation.
*/
inline void GMM::LogProbability(const arma::mat& observation,
arma::vec& logProbs) const
{
// Sum the probability for each Gaussian in our mixture (and we have to
// multiply by the prior for each Gaussian too).
logProbs.set_size(observation.n_cols);
// Store log-probability value in a matrix.
arma::mat logProb(observation.n_cols, gaussians);
// Assign value to the matrix.
for (size_t i = 0; i < gaussians; i++)
{
arma::vec temp(logProb.colptr(i), observation.n_cols, false, true);
dists[i].LogProbability(observation, temp);
}
// Save log(weights) as a vector.
arma::vec logWeights = arma::log(weights);
// Compute log-probability.
logProb += repmat(logWeights.t(), logProb.n_rows, 1);
math::LogSumExp(logProb, logProbs);
}
/**
* Return the probability of the given observation being from this GMM.
*
* @param observation Observation vector to compute probabilty.
*/
inline double GMM::Probability(const arma::vec& observation) const
{
return exp(LogProbability(observation));
}
/**
* Return the probability of the given observation GMM matrix.
*
* @param observation Observation matrix to compute probabilty.
* @param probs Stores the value of probability for x.
*/
inline void GMM::Probability(const arma::mat& observation,
arma::vec& probs) const
{
LogProbability(observation, probs);
probs = exp(probs);
}
/**
* Return the log probability of the given observation being from the given
* component in the mixture.
*
* @param observation Observation vector to compute log-probabilty.
* @param component Calculate the log-probability for given observation vector.
*/
inline double GMM::LogProbability(const arma::vec& observation,
const size_t component) const
{
// We are only considering one Gaussian component -- so we only need to call
// Probability() once. We do consider the prior probability!
return log(weights[component]) + dists[component].LogProbability(observation);
}
/**
* Return the probability of the given observation being from the given
* component in the mixture.
*
* @param observation Observation matrix to compute probabilty.
* @param component Calculate the probability for given component.
*/
inline double GMM::Probability(const arma::vec& observation,
const size_t component) const
{
return exp(LogProbability(observation, component));
}
/**
* Return a randomly generated observation according to the probability
* distribution defined by this object.
*/
inline arma::vec GMM::Random() const
{
// Determine which Gaussian it will be coming from.
double gaussRand = math::Random();
size_t gaussian = 0;
double sumProb = 0;
for (size_t g = 0; g < gaussians; g++)
{
sumProb += weights(g);
if (gaussRand <= sumProb)
{
gaussian = g;
break;
}
}
arma::mat cholDecomp;
if (!arma::chol(cholDecomp, dists[gaussian].Covariance()))
{
Log::Fatal << "Cholesky decomposition failed." << std::endl;
}
return trans(cholDecomp) *
arma::randn<arma::vec>(dimensionality) + dists[gaussian].Mean();
}
/**
* Classify the given observations as being from an individual component in this
* GMM.
*
* @param observation Observation matrix for classification.
* @param labels Save the labels for the given observation matrix.
*/
inline void GMM::Classify(const arma::mat& observations,
arma::Row<size_t>& labels) const
{
// This is not the best way to do this!
// We should not have to fill this with values, because each one should be
// overwritten.
labels.set_size(observations.n_cols);
for (size_t i = 0; i < observations.n_cols; ++i)
{
// Find maximum probability component.
double probability = -std::numeric_limits<double>::infinity();
for (size_t j = 0; j < gaussians; ++j)
{
// We have to use LogProbability() otherwise Probability() would overflow
// easily.
double newProb = LogProbability(observations.unsafe_col(i), j);
if (newProb >= probability)
{
probability = newProb;
labels[i] = j;
}
}
}
}
/**
* Get the log-likelihood of this data's fit to the model.
*
* @param data Data matrix to compute log-likelihood.
* @parma distsL Vector of Gaussian distribution.
* @param weightsL Vector of weights for computing likelihoods.
*/
inline double GMM::LogLikelihood(
const arma::mat& data,
const std::vector<distribution::GaussianDistribution>& distsL,
const arma::vec& weightsL) const
{
double loglikelihood = 0;
arma::vec logPhis;
arma::mat logLikelihoods(gaussians, data.n_cols);
// It has to be LogProbability() otherwise Probability() would overflow easily
for (size_t i = 0; i < gaussians; ++i)
{
distsL[i].LogProbability(data, logPhis);
logLikelihoods.row(i) = log(weightsL(i)) + trans(logPhis);
}
// Now sum over every point.
for (size_t j = 0; j < data.n_cols; ++j)
loglikelihood += mlpack::math::AccuLog(logLikelihoods.col(j));
return loglikelihood;
}
/**
* Fit the GMM to the given observations.
*/
@@ -13,7 +13,7 @@ set(SOURCES
hoeffding_tree.hpp
hoeffding_tree_impl.hpp
hoeffding_tree_model.hpp
hoeffding_tree_model.cpp
hoeffding_tree_model_impl.hpp
information_gain.hpp
numeric_split_info.hpp
typedef.hpp
@@ -15,6 +15,7 @@
#include "hoeffding_tree.hpp"
#include "binary_numeric_split.hpp"
#include "information_gain.hpp"
#include <queue>
namespace mlpack {
namespace tree {
@@ -220,4 +221,7 @@ class HoeffdingTreeModel
} // namespace tree
} // namespace mlpack
// Include implementation.
#include "hoeffding_tree_model_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @param hoeffding_tree_model.cpp
* @file methods/hoeffding_trees/hoeffding_tree_model_impl.hpp
* @author Ryan Curtin
*
* Implementation of the HoeffdingTreeModel class.
@@ -9,15 +9,16 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_HOEFFDING_TREE_HOEFFDING_TREE_MODEL_IMPL_HPP
#define MLPACK_METHODS_HOEFFDING_TREE_HOEFFDING_TREE_MODEL_IMPL_HPP
#include "hoeffding_tree_model.hpp"
#include <queue>
using namespace mlpack;
using namespace mlpack::tree;
namespace mlpack {
namespace tree {
// Constructor.
HoeffdingTreeModel::HoeffdingTreeModel(const TreeType& type) :
inline HoeffdingTreeModel::HoeffdingTreeModel(const TreeType& type) :
type(type),
giniHoeffdingTree(NULL),
giniBinaryTree(NULL),
@@ -28,7 +29,7 @@ HoeffdingTreeModel::HoeffdingTreeModel(const TreeType& type) :
}
// Copy constructor.
HoeffdingTreeModel::HoeffdingTreeModel(const HoeffdingTreeModel& other) :
inline HoeffdingTreeModel::HoeffdingTreeModel(const HoeffdingTreeModel& other) :
type(other.type),
giniHoeffdingTree(other.giniHoeffdingTree ? new GiniHoeffdingTreeType(
*other.giniHoeffdingTree) : NULL),
@@ -43,7 +44,7 @@ HoeffdingTreeModel::HoeffdingTreeModel(const HoeffdingTreeModel& other) :
}
// Move constructor.
HoeffdingTreeModel::HoeffdingTreeModel(HoeffdingTreeModel&& other) :
inline HoeffdingTreeModel::HoeffdingTreeModel(HoeffdingTreeModel&& other) :
type(other.type),
giniHoeffdingTree(other.giniHoeffdingTree),
giniBinaryTree(other.giniBinaryTree),
@@ -59,7 +60,7 @@ HoeffdingTreeModel::HoeffdingTreeModel(HoeffdingTreeModel&& other) :
}
// Copy operator.
HoeffdingTreeModel& HoeffdingTreeModel::operator=(
inline HoeffdingTreeModel& HoeffdingTreeModel::operator=(
const HoeffdingTreeModel& other)
{
if (this != &other)
@@ -90,7 +91,8 @@ HoeffdingTreeModel& HoeffdingTreeModel::operator=(
}
// Move operator.
HoeffdingTreeModel& HoeffdingTreeModel::operator=(HoeffdingTreeModel&& other)
inline HoeffdingTreeModel& HoeffdingTreeModel::operator=(
HoeffdingTreeModel&& other)
{
if (this != &other)
{
@@ -117,7 +119,7 @@ HoeffdingTreeModel& HoeffdingTreeModel::operator=(HoeffdingTreeModel&& other)
}
// Destructor.
HoeffdingTreeModel::~HoeffdingTreeModel()
inline HoeffdingTreeModel::~HoeffdingTreeModel()
{
delete giniHoeffdingTree;
delete giniBinaryTree;
@@ -126,7 +128,7 @@ HoeffdingTreeModel::~HoeffdingTreeModel()
}
// Create the model.
void HoeffdingTreeModel::BuildModel(
inline void HoeffdingTreeModel::BuildModel(
const arma::mat& dataset,
const data::DatasetInfo& datasetInfo,
const arma::Row<size_t>& labels,
@@ -189,9 +191,9 @@ void HoeffdingTreeModel::BuildModel(
}
// Train the model on one pass of the dataset.
void HoeffdingTreeModel::Train(const arma::mat& dataset,
const arma::Row<size_t>& labels,
const bool batchTraining)
inline void HoeffdingTreeModel::Train(const arma::mat& dataset,
const arma::Row<size_t>& labels,
const bool batchTraining)
{
// Depending on the type, pass through once.
switch (type)
@@ -215,8 +217,9 @@ void HoeffdingTreeModel::Train(const arma::mat& dataset,
}
// Classify the given points.
void HoeffdingTreeModel::Classify(const arma::mat& dataset,
arma::Row<size_t>& predictions) const
inline void HoeffdingTreeModel::Classify(const arma::mat& dataset,
arma::Row<size_t>& predictions)
const
{
// Call Classify() with the right model.
switch (type)
@@ -240,9 +243,10 @@ void HoeffdingTreeModel::Classify(const arma::mat& dataset,
}
// Classify the given points.
void HoeffdingTreeModel::Classify(const arma::mat& dataset,
arma::Row<size_t>& predictions,
arma::rowvec& probabilities) const
inline void HoeffdingTreeModel::Classify(const arma::mat& dataset,
arma::Row<size_t>& predictions,
arma::rowvec& probabilities)
const
{
// Call Classify() with the right model.
switch (type)
@@ -286,7 +290,7 @@ size_t CountNodes(TreeType& tree)
}
// Get the number of nodes in the tree.
size_t HoeffdingTreeModel::NumNodes() const
inline size_t HoeffdingTreeModel::NumNodes() const
{
// Call CountNodes() with the right type of tree.
switch (type)
@@ -303,3 +307,8 @@ size_t HoeffdingTreeModel::NumNodes() const
return 0; // This should never happen!
}
} // namespace tree
} // namespace mlpack
#endif
-1
View File
@@ -8,7 +8,6 @@ set(SOURCES
kde_stat.hpp
kde_model.hpp
kde_model_impl.hpp
kde_model.cpp
)
# Add directory name to sources.
-319
View File
@@ -1,319 +0,0 @@
/**
* @file methods/kde/kde_model.cpp
* @author Roberto Hueso
*
* Implementation of KDE Model.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "kde_model.hpp"
namespace mlpack {
namespace kde {
//! Initialize the KDEModel with the given parameters.
KDEModel::KDEModel(const double bandwidth,
const double relError,
const double absError,
const KernelTypes kernelType,
const TreeTypes treeType,
const bool monteCarlo,
const double mcProb,
const size_t initialSampleSize,
const double mcEntryCoef,
const double mcBreakCoef) :
bandwidth(bandwidth),
relError(relError),
absError(absError),
kernelType(kernelType),
treeType(treeType),
monteCarlo(monteCarlo),
mcProb(mcProb),
initialSampleSize(initialSampleSize),
mcEntryCoef(mcEntryCoef),
mcBreakCoef(mcBreakCoef),
kdeModel(NULL)
{
// Nothing to do.
}
// Copy constructor.
KDEModel::KDEModel(const KDEModel& other) :
bandwidth(other.bandwidth),
relError(other.relError),
absError(other.absError),
kernelType(other.kernelType),
treeType(other.treeType),
monteCarlo(other.monteCarlo),
mcProb(other.mcProb),
initialSampleSize(other.initialSampleSize),
mcEntryCoef(other.mcEntryCoef),
mcBreakCoef(other.mcBreakCoef),
kdeModel(other.kdeModel->Clone())
{
// Nothing to do.
}
// Move constructor.
KDEModel::KDEModel(KDEModel&& other) :
bandwidth(other.bandwidth),
relError(other.relError),
absError(other.absError),
kernelType(other.kernelType),
treeType(other.treeType),
monteCarlo(other.monteCarlo),
mcProb(other.mcProb),
initialSampleSize(other.initialSampleSize),
mcEntryCoef(other.mcEntryCoef),
mcBreakCoef(other.mcBreakCoef),
kdeModel(std::move(other.kdeModel))
{
// Reset other model.
other.bandwidth = 1.0;
other.relError = KDEDefaultParams::relError;
other.absError = KDEDefaultParams::absError;
other.kernelType = KernelTypes::GAUSSIAN_KERNEL;
other.treeType = TreeTypes::KD_TREE;
other.monteCarlo = KDEDefaultParams::monteCarlo;
other.mcProb = KDEDefaultParams::mcProb;
other.initialSampleSize = KDEDefaultParams::initialSampleSize;
other.mcEntryCoef = KDEDefaultParams::mcEntryCoef;
other.mcBreakCoef = KDEDefaultParams::mcBreakCoef;
}
KDEModel& KDEModel::operator=(const KDEModel& other)
{
if (this != &other)
{
delete kdeModel;
bandwidth = other.bandwidth;
relError = other.relError;
absError = other.absError;
kernelType = other.kernelType;
treeType = other.treeType;
monteCarlo = other.monteCarlo;
mcProb = other.mcProb;
initialSampleSize = other.initialSampleSize;
mcEntryCoef = other.mcEntryCoef;
mcBreakCoef = other.mcBreakCoef;
kdeModel = other.kdeModel->Clone();
}
return *this;
}
KDEModel& KDEModel::operator=(KDEModel&& other)
{
if (this != &other)
{
delete kdeModel;
bandwidth = other.bandwidth;
relError = other.relError;
absError = other.absError;
kernelType = other.kernelType;
treeType = other.treeType;
monteCarlo = other.monteCarlo;
mcProb = other.mcProb;
initialSampleSize = other.initialSampleSize;
mcEntryCoef = other.mcEntryCoef;
mcBreakCoef = other.mcBreakCoef;
kdeModel = std::move(other.kdeModel);
// Reset other model.
other.bandwidth = 1.0;
other.relError = KDEDefaultParams::relError;
other.absError = KDEDefaultParams::absError;
other.kernelType = KernelTypes::GAUSSIAN_KERNEL;
other.treeType = TreeTypes::KD_TREE;
other.monteCarlo = KDEDefaultParams::monteCarlo;
other.mcProb = KDEDefaultParams::mcProb;
other.initialSampleSize = KDEDefaultParams::initialSampleSize;
other.mcEntryCoef = KDEDefaultParams::mcEntryCoef;
other.mcBreakCoef = KDEDefaultParams::mcBreakCoef;
}
return *this;
}
// Clean memory.
KDEModel::~KDEModel()
{
delete kdeModel;
}
template<template<typename TreeMetricType,
typename TreeMatType,
typename TreeStatType> class TreeType>
KDEWrapperBase* InitializeModelHelper(const KDEModel::KernelTypes kernelType,
const double relError,
const double absError,
const double bandwidth)
{
switch (kernelType)
{
case KDEModel::GAUSSIAN_KERNEL:
return new KDEWrapper<kernel::GaussianKernel, TreeType>(
relError, absError, kernel::GaussianKernel(bandwidth));
case KDEModel::EPANECHNIKOV_KERNEL:
return new KDEWrapper<kernel::EpanechnikovKernel, TreeType>(
relError, absError, kernel::EpanechnikovKernel(bandwidth));
case KDEModel::LAPLACIAN_KERNEL:
return new KDEWrapper<kernel::LaplacianKernel, TreeType>(
relError, absError, kernel::LaplacianKernel(bandwidth));
case KDEModel::SPHERICAL_KERNEL:
return new KDEWrapper<kernel::SphericalKernel, TreeType>(
relError, absError, kernel::SphericalKernel(bandwidth));
case KDEModel::TRIANGULAR_KERNEL:
return new KDEWrapper<kernel::TriangularKernel, TreeType>(
relError, absError, kernel::TriangularKernel(bandwidth));
}
// This should never happen.
return NULL;
}
void KDEModel::InitializeModel()
{
// Clean memory, if necessary.
delete kdeModel;
// Build the actual model.
switch (treeType)
{
case KD_TREE:
kdeModel = InitializeModelHelper<tree::KDTree>(kernelType, relError,
absError, bandwidth);
break;
case BALL_TREE:
kdeModel = InitializeModelHelper<tree::BallTree>(kernelType, relError,
absError, bandwidth);
break;
case COVER_TREE:
kdeModel = InitializeModelHelper<tree::StandardCoverTree>(kernelType,
relError, absError, bandwidth);
break;
case OCTREE:
kdeModel = InitializeModelHelper<tree::Octree>(kernelType, relError,
absError, bandwidth);
break;
case R_TREE:
kdeModel = InitializeModelHelper<tree::RTree>(kernelType, relError,
absError, bandwidth);
break;
}
}
void KDEModel::BuildModel(util::Timers& timers, arma::mat&& referenceSet)
{
InitializeModel();
// Set whether to use Monte Carlo estimations or not.
kdeModel->MonteCarlo() = monteCarlo;
// Set Monte Carlo probability.
kdeModel->MCProb(mcProb);
// Set Monte Carlo initial sample size.
kdeModel->MCInitialSampleSize() = initialSampleSize;
// Set Monte Carlo entry coefficient.
kdeModel->MCEntryCoef(mcEntryCoef);
// Set Monte Carlo break coefficient.
kdeModel->MCBreakCoef(mcBreakCoef);
// Train the model.
kdeModel->Train(timers, std::move(referenceSet));
}
// Perform bichromatic evaluation.
void KDEModel::Evaluate(util::Timers& timers,
arma::mat&& querySet,
arma::vec& estimates)
{
kdeModel->Evaluate(timers, std::move(querySet), estimates);
}
// Perform monochromatic evaluation.
void KDEModel::Evaluate(util::Timers& timers, arma::vec& estimates)
{
kdeModel->Evaluate(timers, estimates);
}
// Clean memory.
void KDEModel::CleanMemory()
{
delete kdeModel;
}
// Modify model kernel bandwidth.
void KDEModel::Bandwidth(const double newBandwidth)
{
bandwidth = newBandwidth;
kdeModel->Bandwidth(bandwidth);
}
// Modify model relative error tolerance.
void KDEModel::RelativeError(const double newRelError)
{
relError = newRelError;
kdeModel->RelativeError(relError);
}
// Modify model absolute error tolerance.
void KDEModel::AbsoluteError(const double newAbsError)
{
absError = newAbsError;
kdeModel->AbsoluteError(absError);
}
// Modify whether Monte Carlo estimations will be used.
void KDEModel::MonteCarlo(const bool newMonteCarlo)
{
monteCarlo = newMonteCarlo;
kdeModel->MonteCarlo() = monteCarlo;
}
// Modify model Monte Carlo probability.
void KDEModel::MCProbability(const double newMCProb)
{
mcProb = newMCProb;
kdeModel->MCProb(mcProb);
}
// Modify model Monte Carlo initial sample size.
void KDEModel::MCInitialSampleSize(const size_t newSampleSize)
{
initialSampleSize = newSampleSize;
kdeModel->MCInitialSampleSize() = initialSampleSize;
}
// Modify model Monte Carlo entry coefficient.
void KDEModel::MCEntryCoefficient(const double newEntryCoef)
{
mcEntryCoef = newEntryCoef;
kdeModel->MCEntryCoef(mcEntryCoef);
}
// Modify model Monte Carlo break coefficient.
void KDEModel::MCBreakCoefficient(const double newBreakCoef)
{
mcBreakCoef = newBreakCoef;
kdeModel->MCBreakCoef(mcBreakCoef);
}
} // namespace kde
} // namespace mlpack
+304
View File
@@ -18,6 +18,310 @@
namespace mlpack {
namespace kde {
//! Initialize the KDEModel with the given parameters.
inline KDEModel::KDEModel(
const double bandwidth,
const double relError,
const double absError,
const KernelTypes kernelType,
const TreeTypes treeType,
const bool monteCarlo,
const double mcProb,
const size_t initialSampleSize,
const double mcEntryCoef,
const double mcBreakCoef) :
bandwidth(bandwidth),
relError(relError),
absError(absError),
kernelType(kernelType),
treeType(treeType),
monteCarlo(monteCarlo),
mcProb(mcProb),
initialSampleSize(initialSampleSize),
mcEntryCoef(mcEntryCoef),
mcBreakCoef(mcBreakCoef),
kdeModel(NULL)
{
// Nothing to do.
}
// Copy constructor.
inline KDEModel::KDEModel(const KDEModel& other) :
bandwidth(other.bandwidth),
relError(other.relError),
absError(other.absError),
kernelType(other.kernelType),
treeType(other.treeType),
monteCarlo(other.monteCarlo),
mcProb(other.mcProb),
initialSampleSize(other.initialSampleSize),
mcEntryCoef(other.mcEntryCoef),
mcBreakCoef(other.mcBreakCoef),
kdeModel(other.kdeModel->Clone())
{
// Nothing to do.
}
// Move constructor.
inline KDEModel::KDEModel(KDEModel&& other) :
bandwidth(other.bandwidth),
relError(other.relError),
absError(other.absError),
kernelType(other.kernelType),
treeType(other.treeType),
monteCarlo(other.monteCarlo),
mcProb(other.mcProb),
initialSampleSize(other.initialSampleSize),
mcEntryCoef(other.mcEntryCoef),
mcBreakCoef(other.mcBreakCoef),
kdeModel(std::move(other.kdeModel))
{
// Reset other model.
other.bandwidth = 1.0;
other.relError = KDEDefaultParams::relError;
other.absError = KDEDefaultParams::absError;
other.kernelType = KernelTypes::GAUSSIAN_KERNEL;
other.treeType = TreeTypes::KD_TREE;
other.monteCarlo = KDEDefaultParams::monteCarlo;
other.mcProb = KDEDefaultParams::mcProb;
other.initialSampleSize = KDEDefaultParams::initialSampleSize;
other.mcEntryCoef = KDEDefaultParams::mcEntryCoef;
other.mcBreakCoef = KDEDefaultParams::mcBreakCoef;
}
inline KDEModel& KDEModel::operator=(const KDEModel& other)
{
if (this != &other)
{
delete kdeModel;
bandwidth = other.bandwidth;
relError = other.relError;
absError = other.absError;
kernelType = other.kernelType;
treeType = other.treeType;
monteCarlo = other.monteCarlo;
mcProb = other.mcProb;
initialSampleSize = other.initialSampleSize;
mcEntryCoef = other.mcEntryCoef;
mcBreakCoef = other.mcBreakCoef;
kdeModel = other.kdeModel->Clone();
}
return *this;
}
inline KDEModel& KDEModel::operator=(KDEModel&& other)
{
if (this != &other)
{
delete kdeModel;
bandwidth = other.bandwidth;
relError = other.relError;
absError = other.absError;
kernelType = other.kernelType;
treeType = other.treeType;
monteCarlo = other.monteCarlo;
mcProb = other.mcProb;
initialSampleSize = other.initialSampleSize;
mcEntryCoef = other.mcEntryCoef;
mcBreakCoef = other.mcBreakCoef;
kdeModel = std::move(other.kdeModel);
// Reset other model.
other.bandwidth = 1.0;
other.relError = KDEDefaultParams::relError;
other.absError = KDEDefaultParams::absError;
other.kernelType = KernelTypes::GAUSSIAN_KERNEL;
other.treeType = TreeTypes::KD_TREE;
other.monteCarlo = KDEDefaultParams::monteCarlo;
other.mcProb = KDEDefaultParams::mcProb;
other.initialSampleSize = KDEDefaultParams::initialSampleSize;
other.mcEntryCoef = KDEDefaultParams::mcEntryCoef;
other.mcBreakCoef = KDEDefaultParams::mcBreakCoef;
}
return *this;
}
// Clean memory.
inline KDEModel::~KDEModel()
{
delete kdeModel;
}
template<template<typename TreeMetricType,
typename TreeMatType,
typename TreeStatType> class TreeType>
KDEWrapperBase* InitializeModelHelper(const KDEModel::KernelTypes kernelType,
const double relError,
const double absError,
const double bandwidth)
{
switch (kernelType)
{
case KDEModel::GAUSSIAN_KERNEL:
return new KDEWrapper<kernel::GaussianKernel, TreeType>(
relError, absError, kernel::GaussianKernel(bandwidth));
case KDEModel::EPANECHNIKOV_KERNEL:
return new KDEWrapper<kernel::EpanechnikovKernel, TreeType>(
relError, absError, kernel::EpanechnikovKernel(bandwidth));
case KDEModel::LAPLACIAN_KERNEL:
return new KDEWrapper<kernel::LaplacianKernel, TreeType>(
relError, absError, kernel::LaplacianKernel(bandwidth));
case KDEModel::SPHERICAL_KERNEL:
return new KDEWrapper<kernel::SphericalKernel, TreeType>(
relError, absError, kernel::SphericalKernel(bandwidth));
case KDEModel::TRIANGULAR_KERNEL:
return new KDEWrapper<kernel::TriangularKernel, TreeType>(
relError, absError, kernel::TriangularKernel(bandwidth));
}
// This should never happen.
return NULL;
}
inline void KDEModel::InitializeModel()
{
// Clean memory, if necessary.
delete kdeModel;
// Build the actual model.
switch (treeType)
{
case KD_TREE:
kdeModel = InitializeModelHelper<tree::KDTree>(kernelType, relError,
absError, bandwidth);
break;
case BALL_TREE:
kdeModel = InitializeModelHelper<tree::BallTree>(kernelType, relError,
absError, bandwidth);
break;
case COVER_TREE:
kdeModel = InitializeModelHelper<tree::StandardCoverTree>(kernelType,
relError, absError, bandwidth);
break;
case OCTREE:
kdeModel = InitializeModelHelper<tree::Octree>(kernelType, relError,
absError, bandwidth);
break;
case R_TREE:
kdeModel = InitializeModelHelper<tree::RTree>(kernelType, relError,
absError, bandwidth);
break;
}
}
inline void KDEModel::BuildModel(util::Timers& timers,
arma::mat&& referenceSet)
{
InitializeModel();
// Set whether to use Monte Carlo estimations or not.
kdeModel->MonteCarlo() = monteCarlo;
// Set Monte Carlo probability.
kdeModel->MCProb(mcProb);
// Set Monte Carlo initial sample size.
kdeModel->MCInitialSampleSize() = initialSampleSize;
// Set Monte Carlo entry coefficient.
kdeModel->MCEntryCoef(mcEntryCoef);
// Set Monte Carlo break coefficient.
kdeModel->MCBreakCoef(mcBreakCoef);
// Train the model.
kdeModel->Train(timers, std::move(referenceSet));
}
// Perform bichromatic evaluation.
inline void KDEModel::Evaluate(util::Timers& timers,
arma::mat&& querySet,
arma::vec& estimates)
{
kdeModel->Evaluate(timers, std::move(querySet), estimates);
}
// Perform monochromatic evaluation.
inline void KDEModel::Evaluate(util::Timers& timers,
arma::vec& estimates)
{
kdeModel->Evaluate(timers, estimates);
}
// Clean memory.
inline void KDEModel::CleanMemory()
{
delete kdeModel;
}
// Modify model kernel bandwidth.
inline void KDEModel::Bandwidth(const double newBandwidth)
{
bandwidth = newBandwidth;
kdeModel->Bandwidth(bandwidth);
}
// Modify model relative error tolerance.
inline void KDEModel::RelativeError(const double newRelError)
{
relError = newRelError;
kdeModel->RelativeError(relError);
}
// Modify model absolute error tolerance.
inline void KDEModel::AbsoluteError(const double newAbsError)
{
absError = newAbsError;
kdeModel->AbsoluteError(absError);
}
// Modify whether Monte Carlo estimations will be used.
inline void KDEModel::MonteCarlo(const bool newMonteCarlo)
{
monteCarlo = newMonteCarlo;
kdeModel->MonteCarlo() = monteCarlo;
}
// Modify model Monte Carlo probability.
inline void KDEModel::MCProbability(const double newMCProb)
{
mcProb = newMCProb;
kdeModel->MCProb(mcProb);
}
// Modify model Monte Carlo initial sample size.
inline void KDEModel::MCInitialSampleSize(const size_t newSampleSize)
{
initialSampleSize = newSampleSize;
kdeModel->MCInitialSampleSize() = initialSampleSize;
}
// Modify model Monte Carlo entry coefficient.
inline void KDEModel::MCEntryCoefficient(const double newEntryCoef)
{
mcEntryCoef = newEntryCoef;
kdeModel->MCEntryCoef(mcEntryCoef);
}
// Modify model Monte Carlo break coefficient.
inline void KDEModel::MCBreakCoefficient(const double newBreakCoef)
{
mcBreakCoef = newBreakCoef;
kdeModel->MCBreakCoef(mcBreakCoef);
}
//! Train the model (build the tree).
template<typename KernelType,
template<typename TreeMetricType,
-1
View File
@@ -3,7 +3,6 @@
set(SOURCES
lars.hpp
lars_impl.hpp
lars.cpp
)
# add directory name to sources
-652
View File
@@ -1,652 +0,0 @@
/**
* @file methods/lars/lars.cpp
* @author Nishant Mehta (niche)
*
* Implementation of LARS and LASSO.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "lars.hpp"
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/util/timers.hpp>
using namespace mlpack;
using namespace mlpack::regression;
LARS::LARS(const bool useCholesky,
const double lambda1,
const double lambda2,
const double tolerance) :
matGram(&matGramInternal),
useCholesky(useCholesky),
lasso((lambda1 != 0)),
lambda1(lambda1),
elasticNet((lambda1 != 0) && (lambda2 != 0)),
lambda2(lambda2),
tolerance(tolerance)
{ /* Nothing left to do. */ }
LARS::LARS(const bool useCholesky,
const arma::mat& gramMatrix,
const double lambda1,
const double lambda2,
const double tolerance) :
matGram(&gramMatrix),
useCholesky(useCholesky),
lasso((lambda1 != 0)),
lambda1(lambda1),
elasticNet((lambda1 != 0) && (lambda2 != 0)),
lambda2(lambda2),
tolerance(tolerance)
{ /* Nothing left to do */ }
LARS::LARS(const arma::mat& data,
const arma::rowvec& responses,
const bool transposeData,
const bool useCholesky,
const double lambda1,
const double lambda2,
const double tolerance) :
LARS(useCholesky, lambda1, lambda2, tolerance)
{
Train(data, responses, transposeData);
}
LARS::LARS(const arma::mat& data,
const arma::rowvec& responses,
const bool transposeData,
const bool useCholesky,
const arma::mat& gramMatrix,
const double lambda1,
const double lambda2,
const double tolerance) :
LARS(useCholesky, gramMatrix, lambda1, lambda2, tolerance)
{
Train(data, responses, transposeData);
}
// Copy Constructor.
LARS::LARS(const LARS& other) :
matGramInternal(other.matGramInternal),
matGram(other.matGram != &other.matGramInternal ?
other.matGram : &matGramInternal),
matUtriCholFactor(other.matUtriCholFactor),
useCholesky(other.useCholesky),
lasso(other.lasso),
lambda1(other.lambda1),
elasticNet(other.elasticNet),
lambda2(other.lambda2),
tolerance(other.tolerance),
betaPath(other.betaPath),
lambdaPath(other.lambdaPath),
activeSet(other.activeSet),
isActive(other.isActive),
ignoreSet(other.ignoreSet),
isIgnored(other.isIgnored)
{
// Nothing to do here.
}
// Move constructor.
LARS::LARS(LARS&& other) :
matGramInternal(std::move(other.matGramInternal)),
matGram(other.matGram != &other.matGramInternal ?
other.matGram : &matGramInternal),
matUtriCholFactor(std::move(other.matUtriCholFactor)),
useCholesky(other.useCholesky),
lasso(other.lasso),
lambda1(other.lambda1),
elasticNet(other.elasticNet),
lambda2(other.lambda2),
tolerance(other.tolerance),
betaPath(std::move(other.betaPath)),
lambdaPath(std::move(other.lambdaPath)),
activeSet(std::move(other.activeSet)),
isActive(std::move(other.isActive)),
ignoreSet(std::move(other.ignoreSet)),
isIgnored(std::move(other.isIgnored))
{
// Nothing to do here.
}
// Copy operator.
LARS& LARS::operator=(const LARS& other)
{
if (&other == this)
return *this;
matGramInternal = other.matGramInternal;
matGram = other.matGram != &other.matGramInternal ?
other.matGram : &matGramInternal;
matUtriCholFactor = other.matUtriCholFactor;
useCholesky = other.useCholesky;
lasso = other.lasso;
lambda1 = other.lambda1;
elasticNet = other.elasticNet;
lambda2 = other.lambda2;
tolerance = other.tolerance;
betaPath = other.betaPath;
lambdaPath = other.lambdaPath;
activeSet = other.activeSet;
isActive = other.isActive;
ignoreSet = other.ignoreSet;
isIgnored = other.isIgnored;
return *this;
}
// Move Operator.
LARS& LARS::operator=(LARS&& other)
{
if (&other == this)
return *this;
matGramInternal = std::move(other.matGramInternal);
matGram = other.matGram != &other.matGramInternal ?
other.matGram : &matGramInternal;
matUtriCholFactor = std::move(other.matUtriCholFactor);
useCholesky = other.useCholesky;
lasso = other.lasso;
lambda1 = other.lambda1;
elasticNet = other.elasticNet;
lambda2 = other.lambda2;
tolerance = other.tolerance;
betaPath = std::move(other.betaPath);
lambdaPath = std::move(other.lambdaPath);
activeSet = std::move(other.activeSet);
isActive = std::move(other.isActive);
ignoreSet = std::move(other.ignoreSet);
isIgnored = std::move(other.isIgnored);
return *this;
}
double LARS::Train(const arma::mat& matX,
const arma::rowvec& y,
arma::vec& beta,
const bool transposeData)
{
// Clear any previous solution information.
betaPath.clear();
lambdaPath.clear();
activeSet.clear();
isActive.clear();
ignoreSet.clear();
isIgnored.clear();
matUtriCholFactor.reset();
// Update values in case lambda1 or lambda2 changed.
lasso = (lambda1 != 0);
elasticNet = (lambda1 != 0 && lambda2 != 0);
// This matrix may end up holding the transpose -- if necessary.
arma::mat dataTrans;
// dataRef is row-major.
const arma::mat& dataRef = (transposeData ? dataTrans : matX);
if (transposeData)
dataTrans = trans(matX);
// Compute X' * y.
arma::vec vecXTy = trans(y * dataRef);
// Set up active set variables. In the beginning, the active set has size 0
// (all dimensions are inactive).
isActive.resize(dataRef.n_cols, false);
// Set up ignores set variables. Initialized empty.
isIgnored.resize(dataRef.n_cols, false);
// Initialize yHat and beta.
beta = arma::zeros(dataRef.n_cols);
arma::vec yHat = arma::zeros(dataRef.n_rows);
arma::vec yHatDirection(dataRef.n_rows);
bool lassocond = false;
// Compute the initial maximum correlation among all dimensions.
arma::vec corr = vecXTy;
double maxCorr = 0;
size_t changeInd = 0;
for (size_t i = 0; i < vecXTy.n_elem; ++i)
{
if (fabs(corr(i)) > maxCorr)
{
maxCorr = fabs(corr(i));
changeInd = i;
}
}
betaPath.push_back(beta);
lambdaPath.push_back(maxCorr);
// If the maximum correlation is too small, there is no reason to continue.
if (maxCorr < lambda1)
{
lambdaPath[0] = lambda1;
return maxCorr;
}
// Compute the Gram matrix. If this is the elastic net problem, we will add
// lambda2 * I_n to the matrix.
if (matGram->n_elem != dataRef.n_cols * dataRef.n_cols)
{
// In this case, matGram should reference matGramInternal.
matGramInternal = trans(dataRef) * dataRef;
if (elasticNet && !useCholesky)
matGramInternal += lambda2 * arma::eye(dataRef.n_cols, dataRef.n_cols);
}
// Main loop.
while (((activeSet.size() + ignoreSet.size()) < dataRef.n_cols) &&
(maxCorr > tolerance))
{
// Compute the maximum correlation among inactive dimensions.
maxCorr = 0;
for (size_t i = 0; i < dataRef.n_cols; ++i)
{
if ((!isActive[i]) && (!isIgnored[i]) && (fabs(corr(i)) > maxCorr))
{
maxCorr = fabs(corr(i));
changeInd = i;
}
}
if (!lassocond)
{
if (useCholesky)
{
// vec newGramCol = vec(activeSet.size());
// for (size_t i = 0; i < activeSet.size(); ++i)
// {
// newGramCol[i] = dot(matX.col(activeSet[i]), matX.col(changeInd));
// }
// This is equivalent to the above 5 lines.
arma::vec newGramCol = matGram->elem(changeInd * dataRef.n_cols +
arma::conv_to<arma::uvec>::from(activeSet));
CholeskyInsert((*matGram)(changeInd, changeInd), newGramCol);
}
// Add variable to active set.
Activate(changeInd);
}
// Compute signs of correlations.
arma::vec s = arma::vec(activeSet.size());
for (size_t i = 0; i < activeSet.size(); ++i)
s(i) = corr(activeSet[i]) / fabs(corr(activeSet[i]));
// Compute the "equiangular" direction in parameter space (betaDirection).
// We use quotes because in the case of non-unit norm variables, this need
// not be equiangular.
arma::vec unnormalizedBetaDirection;
double normalization;
arma::vec betaDirection;
if (useCholesky)
{
// Check for singularity.
const double lastUtriElement = matUtriCholFactor(
matUtriCholFactor.n_cols - 1, matUtriCholFactor.n_rows - 1);
if (std::abs(lastUtriElement) > tolerance)
{
// Ok, no singularity.
/**
* Note that:
* R^T R % S^T % S = (R % S)^T (R % S)
* Now, for 1 the ones vector:
* inv( (R % S)^T (R % S) ) 1
* = inv(R % S) inv((R % S)^T) 1
* = inv(R % S) Solve((R % S)^T, 1)
* = inv(R % S) Solve(R^T, s)
* = Solve(R % S, Solve(R^T, s)
* = s % Solve(R, Solve(R^T, s))
*/
unnormalizedBetaDirection = solve(trimatu(matUtriCholFactor),
solve(trimatl(trans(matUtriCholFactor)), s));
normalization = 1.0 / sqrt(dot(s, unnormalizedBetaDirection));
betaDirection = normalization * unnormalizedBetaDirection;
}
else
{
// Singularity, so remove variable from active set, add to ignores set,
// and look for new variable to add.
Log::Warn << "Encountered singularity when adding variable "
<< changeInd << " to active set; permanently removing."
<< std::endl;
Deactivate(activeSet.size() - 1);
Ignore(changeInd);
CholeskyDelete(matUtriCholFactor.n_rows - 1);
continue;
}
}
else
{
arma::mat matGramActive = arma::mat(activeSet.size(), activeSet.size());
for (size_t i = 0; i < activeSet.size(); ++i)
for (size_t j = 0; j < activeSet.size(); ++j)
matGramActive(i, j) = (*matGram)(activeSet[i], activeSet[j]);
// Check for singularity.
arma::mat matS = s * arma::ones<arma::mat>(1, activeSet.size());
const bool solvedOk = solve(unnormalizedBetaDirection,
matGramActive % trans(matS) % matS,
arma::ones<arma::mat>(activeSet.size(), 1));
if (solvedOk)
{
// Ok, no singularity.
normalization = 1.0 / sqrt(sum(unnormalizedBetaDirection));
betaDirection = normalization * unnormalizedBetaDirection % s;
}
else
{
// Singularity, so remove variable from active set, add to ignores set,
// and look for new variable to add.
Deactivate(activeSet.size() - 1);
Ignore(changeInd);
Log::Warn << "Encountered singularity when adding variable "
<< changeInd << " to active set; permanently removing."
<< std::endl;
continue;
}
}
// compute "equiangular" direction in output space
ComputeYHatDirection(dataRef, betaDirection, yHatDirection);
double gamma = maxCorr / normalization;
// If not all variables are active.
if ((activeSet.size() + ignoreSet.size()) < dataRef.n_cols)
{
// Compute correlations with direction.
for (size_t ind = 0; ind < dataRef.n_cols; ind++)
{
if (isActive[ind] || isIgnored[ind])
continue;
const double dirCorr = dot(dataRef.col(ind), yHatDirection);
const double val1 = (maxCorr - corr(ind)) / (normalization - dirCorr);
const double val2 = (maxCorr + corr(ind)) / (normalization + dirCorr);
if ((val1 > 0.0) && (val1 < gamma))
gamma = val1;
if ((val2 > 0.0) && (val2 < gamma))
gamma = val2;
// Handle edge case where the largest actually is equal to 0.
if (std::max(val1, val2) == 0.0)
gamma = 0.0;
}
}
// Bound gamma according to LASSO.
if (lasso)
{
lassocond = false;
double lassoboundOnGamma = DBL_MAX;
size_t activeIndToKickOut = -1;
for (size_t i = 0; i < activeSet.size(); ++i)
{
double val = -beta(activeSet[i]) / betaDirection(i);
if ((val > 0) && (val < lassoboundOnGamma))
{
lassoboundOnGamma = val;
activeIndToKickOut = i;
}
}
if (lassoboundOnGamma < gamma)
{
gamma = lassoboundOnGamma;
lassocond = true;
changeInd = activeIndToKickOut;
}
}
// Update the prediction.
yHat += gamma * yHatDirection;
// Update the estimator.
for (size_t i = 0; i < activeSet.size(); ++i)
{
beta(activeSet[i]) += gamma * betaDirection(i);
}
// Sanity check to make sure the kicked out dimension is actually zero.
if (lassocond)
{
if (beta(activeSet[changeInd]) != 0)
beta(activeSet[changeInd]) = 0;
}
betaPath.push_back(beta);
if (lassocond)
{
// Index is in position changeInd in activeSet.
if (useCholesky)
CholeskyDelete(changeInd);
Deactivate(changeInd);
}
corr = vecXTy - trans(dataRef) * yHat;
if (elasticNet)
corr -= lambda2 * beta;
double curLambda = 0;
for (size_t i = 0; i < activeSet.size(); ++i)
curLambda += fabs(corr(activeSet[i]));
curLambda /= ((double) activeSet.size());
lambdaPath.push_back(curLambda);
// Time to stop for LASSO?
if (lasso)
{
if (curLambda <= lambda1)
{
InterpolateBeta();
break;
}
}
}
// Unfortunate copy...
beta = betaPath.back();
return ComputeError(matX, y, !transposeData);
}
double LARS::Train(const arma::mat& data,
const arma::rowvec& responses,
const bool transposeData)
{
arma::vec beta;
return Train(data, responses, beta, transposeData);
}
void LARS::Predict(const arma::mat& points,
arma::rowvec& predictions,
const bool rowMajor) const
{
// We really only need to store beta internally...
if (rowMajor)
predictions = trans(points * betaPath.back());
else
predictions = betaPath.back().t() * points;
}
// Private functions.
void LARS::Deactivate(const size_t activeVarInd)
{
isActive[activeSet[activeVarInd]] = false;
activeSet.erase(activeSet.begin() + activeVarInd);
}
void LARS::Activate(const size_t varInd)
{
isActive[varInd] = true;
activeSet.push_back(varInd);
}
void LARS::Ignore(const size_t varInd)
{
isIgnored[varInd] = true;
ignoreSet.push_back(varInd);
}
void LARS::ComputeYHatDirection(const arma::mat& matX,
const arma::vec& betaDirection,
arma::vec& yHatDirection)
{
yHatDirection.fill(0);
for (size_t i = 0; i < activeSet.size(); ++i)
yHatDirection += betaDirection(i) * matX.col(activeSet[i]);
}
void LARS::InterpolateBeta()
{
int pathLength = betaPath.size();
// interpolate beta and stop
double ultimateLambda = lambdaPath[pathLength - 1];
double penultimateLambda = lambdaPath[pathLength - 2];
double interp = (penultimateLambda - lambda1)
/ (penultimateLambda - ultimateLambda);
betaPath[pathLength - 1] = (1 - interp) * (betaPath[pathLength - 2])
+ interp * betaPath[pathLength - 1];
lambdaPath[pathLength - 1] = lambda1;
}
void LARS::CholeskyInsert(const arma::vec& newX, const arma::mat& X)
{
if (matUtriCholFactor.n_rows == 0)
{
matUtriCholFactor = arma::mat(1, 1);
if (elasticNet)
matUtriCholFactor(0, 0) = sqrt(dot(newX, newX) + lambda2);
else
matUtriCholFactor(0, 0) = norm(newX, 2);
}
else
{
arma::vec newGramCol = trans(X) * newX;
CholeskyInsert(dot(newX, newX), newGramCol);
}
}
void LARS::CholeskyInsert(double sqNormNewX, const arma::vec& newGramCol)
{
int n = matUtriCholFactor.n_rows;
if (n == 0)
{
matUtriCholFactor = arma::mat(1, 1);
if (elasticNet)
matUtriCholFactor(0, 0) = sqrt(sqNormNewX + lambda2);
else
matUtriCholFactor(0, 0) = sqrt(sqNormNewX);
}
else
{
arma::mat matNewR = arma::mat(n + 1, n + 1);
if (elasticNet)
sqNormNewX += lambda2;
arma::vec matUtriCholFactork = solve(trimatl(trans(matUtriCholFactor)),
newGramCol);
matNewR(arma::span(0, n - 1), arma::span(0, n - 1)) = matUtriCholFactor;
matNewR(arma::span(0, n - 1), n) = matUtriCholFactork;
matNewR(n, arma::span(0, n - 1)).fill(0.0);
matNewR(n, n) = sqrt(sqNormNewX - dot(matUtriCholFactork,
matUtriCholFactork));
matUtriCholFactor = matNewR;
}
}
void LARS::GivensRotate(const arma::vec::fixed<2>& x,
arma::vec::fixed<2>& rotatedX,
arma::mat& matG)
{
if (x(1) == 0)
{
matG = arma::eye(2, 2);
rotatedX = x;
}
else
{
double r = norm(x, 2);
matG = arma::mat(2, 2);
double scaledX1 = x(0) / r;
double scaledX2 = x(1) / r;
matG(0, 0) = scaledX1;
matG(1, 0) = -scaledX2;
matG(0, 1) = scaledX2;
matG(1, 1) = scaledX1;
rotatedX = arma::vec(2);
rotatedX(0) = r;
rotatedX(1) = 0;
}
}
void LARS::CholeskyDelete(const size_t colToKill)
{
size_t n = matUtriCholFactor.n_rows;
if (colToKill == (n - 1))
{
matUtriCholFactor = matUtriCholFactor(arma::span(0, n - 2),
arma::span(0, n - 2));
}
else
{
matUtriCholFactor.shed_col(colToKill); // remove column colToKill
n--;
for (size_t k = colToKill; k < n; ++k)
{
arma::mat matG;
arma::vec::fixed<2> rotatedVec;
GivensRotate(matUtriCholFactor(arma::span(k, k + 1), k), rotatedVec,
matG);
matUtriCholFactor(arma::span(k, k + 1), k) = rotatedVec;
if (k < n - 1)
{
matUtriCholFactor(arma::span(k, k + 1), arma::span(k + 1, n - 1)) =
matG * matUtriCholFactor(arma::span(k, k + 1),
arma::span(k + 1, n - 1));
}
}
matUtriCholFactor.shed_row(n);
}
}
double LARS::ComputeError(const arma::mat& matX,
const arma::rowvec& y,
const bool rowMajor)
{
if (rowMajor)
{
return arma::accu(arma::pow(y - trans(matX * betaPath.back()), 2.0));
}
else
{
return arma::accu(arma::pow(y - betaPath.back().t() * matX, 2.0));
}
}
+2
View File
@@ -25,6 +25,8 @@
#define MLPACK_METHODS_LARS_LARS_HPP
#include <mlpack/prereqs.hpp>
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/util/timers.hpp>
namespace mlpack {
namespace regression {
+641
View File
@@ -18,6 +18,647 @@
namespace mlpack {
namespace regression {
inline LARS::LARS(
const bool useCholesky,
const double lambda1,
const double lambda2,
const double tolerance) :
matGram(&matGramInternal),
useCholesky(useCholesky),
lasso((lambda1 != 0)),
lambda1(lambda1),
elasticNet((lambda1 != 0) && (lambda2 != 0)),
lambda2(lambda2),
tolerance(tolerance)
{ /* Nothing left to do. */ }
inline LARS::LARS(
const bool useCholesky,
const arma::mat& gramMatrix,
const double lambda1,
const double lambda2,
const double tolerance) :
matGram(&gramMatrix),
useCholesky(useCholesky),
lasso((lambda1 != 0)),
lambda1(lambda1),
elasticNet((lambda1 != 0) && (lambda2 != 0)),
lambda2(lambda2),
tolerance(tolerance)
{ /* Nothing left to do */ }
inline LARS::LARS(
const arma::mat& data,
const arma::rowvec& responses,
const bool transposeData,
const bool useCholesky,
const double lambda1,
const double lambda2,
const double tolerance) :
LARS(useCholesky, lambda1, lambda2, tolerance)
{
Train(data, responses, transposeData);
}
inline LARS::LARS(
const arma::mat& data,
const arma::rowvec& responses,
const bool transposeData,
const bool useCholesky,
const arma::mat& gramMatrix,
const double lambda1,
const double lambda2,
const double tolerance) :
LARS(useCholesky, gramMatrix, lambda1, lambda2, tolerance)
{
Train(data, responses, transposeData);
}
// Copy Constructor.
inline LARS::LARS(const LARS& other) :
matGramInternal(other.matGramInternal),
matGram(other.matGram != &other.matGramInternal ?
other.matGram : &matGramInternal),
matUtriCholFactor(other.matUtriCholFactor),
useCholesky(other.useCholesky),
lasso(other.lasso),
lambda1(other.lambda1),
elasticNet(other.elasticNet),
lambda2(other.lambda2),
tolerance(other.tolerance),
betaPath(other.betaPath),
lambdaPath(other.lambdaPath),
activeSet(other.activeSet),
isActive(other.isActive),
ignoreSet(other.ignoreSet),
isIgnored(other.isIgnored)
{
// Nothing to do here.
}
// Move constructor.
inline LARS::LARS(LARS&& other) :
matGramInternal(std::move(other.matGramInternal)),
matGram(other.matGram != &other.matGramInternal ?
other.matGram : &matGramInternal),
matUtriCholFactor(std::move(other.matUtriCholFactor)),
useCholesky(other.useCholesky),
lasso(other.lasso),
lambda1(other.lambda1),
elasticNet(other.elasticNet),
lambda2(other.lambda2),
tolerance(other.tolerance),
betaPath(std::move(other.betaPath)),
lambdaPath(std::move(other.lambdaPath)),
activeSet(std::move(other.activeSet)),
isActive(std::move(other.isActive)),
ignoreSet(std::move(other.ignoreSet)),
isIgnored(std::move(other.isIgnored))
{
// Nothing to do here.
}
// Copy operator.
inline LARS& LARS::operator=(const LARS& other)
{
if (&other == this)
return *this;
matGramInternal = other.matGramInternal;
matGram = other.matGram != &other.matGramInternal ?
other.matGram : &matGramInternal;
matUtriCholFactor = other.matUtriCholFactor;
useCholesky = other.useCholesky;
lasso = other.lasso;
lambda1 = other.lambda1;
elasticNet = other.elasticNet;
lambda2 = other.lambda2;
tolerance = other.tolerance;
betaPath = other.betaPath;
lambdaPath = other.lambdaPath;
activeSet = other.activeSet;
isActive = other.isActive;
ignoreSet = other.ignoreSet;
isIgnored = other.isIgnored;
return *this;
}
// Move Operator.
inline LARS& LARS::operator=(LARS&& other)
{
if (&other == this)
return *this;
matGramInternal = std::move(other.matGramInternal);
matGram = other.matGram != &other.matGramInternal ?
other.matGram : &matGramInternal;
matUtriCholFactor = std::move(other.matUtriCholFactor);
useCholesky = other.useCholesky;
lasso = other.lasso;
lambda1 = other.lambda1;
elasticNet = other.elasticNet;
lambda2 = other.lambda2;
tolerance = other.tolerance;
betaPath = std::move(other.betaPath);
lambdaPath = std::move(other.lambdaPath);
activeSet = std::move(other.activeSet);
isActive = std::move(other.isActive);
ignoreSet = std::move(other.ignoreSet);
isIgnored = std::move(other.isIgnored);
return *this;
}
inline double LARS::Train(const arma::mat& matX,
const arma::rowvec& y,
arma::vec& beta,
const bool transposeData)
{
// Clear any previous solution information.
betaPath.clear();
lambdaPath.clear();
activeSet.clear();
isActive.clear();
ignoreSet.clear();
isIgnored.clear();
matUtriCholFactor.reset();
// Update values in case lambda1 or lambda2 changed.
lasso = (lambda1 != 0);
elasticNet = (lambda1 != 0 && lambda2 != 0);
// This matrix may end up holding the transpose -- if necessary.
arma::mat dataTrans;
// dataRef is row-major.
const arma::mat& dataRef = (transposeData ? dataTrans : matX);
if (transposeData)
dataTrans = trans(matX);
// Compute X' * y.
arma::vec vecXTy = trans(y * dataRef);
// Set up active set variables. In the beginning, the active set has size 0
// (all dimensions are inactive).
isActive.resize(dataRef.n_cols, false);
// Set up ignores set variables. Initialized empty.
isIgnored.resize(dataRef.n_cols, false);
// Initialize yHat and beta.
beta = arma::zeros(dataRef.n_cols);
arma::vec yHat = arma::zeros(dataRef.n_rows);
arma::vec yHatDirection(dataRef.n_rows);
bool lassocond = false;
// Compute the initial maximum correlation among all dimensions.
arma::vec corr = vecXTy;
double maxCorr = 0;
size_t changeInd = 0;
for (size_t i = 0; i < vecXTy.n_elem; ++i)
{
if (fabs(corr(i)) > maxCorr)
{
maxCorr = fabs(corr(i));
changeInd = i;
}
}
betaPath.push_back(beta);
lambdaPath.push_back(maxCorr);
// If the maximum correlation is too small, there is no reason to continue.
if (maxCorr < lambda1)
{
lambdaPath[0] = lambda1;
return maxCorr;
}
// Compute the Gram matrix. If this is the elastic net problem, we will add
// lambda2 * I_n to the matrix.
if (matGram->n_elem != dataRef.n_cols * dataRef.n_cols)
{
// In this case, matGram should reference matGramInternal.
matGramInternal = trans(dataRef) * dataRef;
if (elasticNet && !useCholesky)
matGramInternal += lambda2 * arma::eye(dataRef.n_cols, dataRef.n_cols);
}
// Main loop.
while (((activeSet.size() + ignoreSet.size()) < dataRef.n_cols) &&
(maxCorr > tolerance))
{
// Compute the maximum correlation among inactive dimensions.
maxCorr = 0;
for (size_t i = 0; i < dataRef.n_cols; ++i)
{
if ((!isActive[i]) && (!isIgnored[i]) && (fabs(corr(i)) > maxCorr))
{
maxCorr = fabs(corr(i));
changeInd = i;
}
}
if (!lassocond)
{
if (useCholesky)
{
// vec newGramCol = vec(activeSet.size());
// for (size_t i = 0; i < activeSet.size(); ++i)
// {
// newGramCol[i] = dot(matX.col(activeSet[i]), matX.col(changeInd));
// }
// This is equivalent to the above 5 lines.
arma::vec newGramCol = matGram->elem(changeInd * dataRef.n_cols +
arma::conv_to<arma::uvec>::from(activeSet));
CholeskyInsert((*matGram)(changeInd, changeInd), newGramCol);
}
// Add variable to active set.
Activate(changeInd);
}
// Compute signs of correlations.
arma::vec s = arma::vec(activeSet.size());
for (size_t i = 0; i < activeSet.size(); ++i)
s(i) = corr(activeSet[i]) / fabs(corr(activeSet[i]));
// Compute the "equiangular" direction in parameter space (betaDirection).
// We use quotes because in the case of non-unit norm variables, this need
// not be equiangular.
arma::vec unnormalizedBetaDirection;
double normalization;
arma::vec betaDirection;
if (useCholesky)
{
// Check for singularity.
const double lastUtriElement = matUtriCholFactor(
matUtriCholFactor.n_cols - 1, matUtriCholFactor.n_rows - 1);
if (std::abs(lastUtriElement) > tolerance)
{
// Ok, no singularity.
/**
* Note that:
* R^T R % S^T % S = (R % S)^T (R % S)
* Now, for 1 the ones vector:
* inv( (R % S)^T (R % S) ) 1
* = inv(R % S) inv((R % S)^T) 1
* = inv(R % S) Solve((R % S)^T, 1)
* = inv(R % S) Solve(R^T, s)
* = Solve(R % S, Solve(R^T, s)
* = s % Solve(R, Solve(R^T, s))
*/
unnormalizedBetaDirection = solve(trimatu(matUtriCholFactor),
solve(trimatl(trans(matUtriCholFactor)), s));
normalization = 1.0 / sqrt(dot(s, unnormalizedBetaDirection));
betaDirection = normalization * unnormalizedBetaDirection;
}
else
{
// Singularity, so remove variable from active set, add to ignores set,
// and look for new variable to add.
Log::Warn << "Encountered singularity when adding variable "
<< changeInd << " to active set; permanently removing."
<< std::endl;
Deactivate(activeSet.size() - 1);
Ignore(changeInd);
CholeskyDelete(matUtriCholFactor.n_rows - 1);
continue;
}
}
else
{
arma::mat matGramActive = arma::mat(activeSet.size(), activeSet.size());
for (size_t i = 0; i < activeSet.size(); ++i)
for (size_t j = 0; j < activeSet.size(); ++j)
matGramActive(i, j) = (*matGram)(activeSet[i], activeSet[j]);
// Check for singularity.
arma::mat matS = s * arma::ones<arma::mat>(1, activeSet.size());
const bool solvedOk = solve(unnormalizedBetaDirection,
matGramActive % trans(matS) % matS,
arma::ones<arma::mat>(activeSet.size(), 1));
if (solvedOk)
{
// Ok, no singularity.
normalization = 1.0 / sqrt(sum(unnormalizedBetaDirection));
betaDirection = normalization * unnormalizedBetaDirection % s;
}
else
{
// Singularity, so remove variable from active set, add to ignores set,
// and look for new variable to add.
Deactivate(activeSet.size() - 1);
Ignore(changeInd);
Log::Warn << "Encountered singularity when adding variable "
<< changeInd << " to active set; permanently removing."
<< std::endl;
continue;
}
}
// compute "equiangular" direction in output space
ComputeYHatDirection(dataRef, betaDirection, yHatDirection);
double gamma = maxCorr / normalization;
// If not all variables are active.
if ((activeSet.size() + ignoreSet.size()) < dataRef.n_cols)
{
// Compute correlations with direction.
for (size_t ind = 0; ind < dataRef.n_cols; ind++)
{
if (isActive[ind] || isIgnored[ind])
continue;
const double dirCorr = dot(dataRef.col(ind), yHatDirection);
const double val1 = (maxCorr - corr(ind)) / (normalization - dirCorr);
const double val2 = (maxCorr + corr(ind)) / (normalization + dirCorr);
if ((val1 > 0.0) && (val1 < gamma))
gamma = val1;
if ((val2 > 0.0) && (val2 < gamma))
gamma = val2;
// Handle edge case where the largest actually is equal to 0.
if (std::max(val1, val2) == 0.0)
gamma = 0.0;
}
}
// Bound gamma according to LASSO.
if (lasso)
{
lassocond = false;
double lassoboundOnGamma = DBL_MAX;
size_t activeIndToKickOut = -1;
for (size_t i = 0; i < activeSet.size(); ++i)
{
double val = -beta(activeSet[i]) / betaDirection(i);
if ((val > 0) && (val < lassoboundOnGamma))
{
lassoboundOnGamma = val;
activeIndToKickOut = i;
}
}
if (lassoboundOnGamma < gamma)
{
gamma = lassoboundOnGamma;
lassocond = true;
changeInd = activeIndToKickOut;
}
}
// Update the prediction.
yHat += gamma * yHatDirection;
// Update the estimator.
for (size_t i = 0; i < activeSet.size(); ++i)
{
beta(activeSet[i]) += gamma * betaDirection(i);
}
// Sanity check to make sure the kicked out dimension is actually zero.
if (lassocond)
{
if (beta(activeSet[changeInd]) != 0)
beta(activeSet[changeInd]) = 0;
}
betaPath.push_back(beta);
if (lassocond)
{
// Index is in position changeInd in activeSet.
if (useCholesky)
CholeskyDelete(changeInd);
Deactivate(changeInd);
}
corr = vecXTy - trans(dataRef) * yHat;
if (elasticNet)
corr -= lambda2 * beta;
double curLambda = 0;
for (size_t i = 0; i < activeSet.size(); ++i)
curLambda += fabs(corr(activeSet[i]));
curLambda /= ((double) activeSet.size());
lambdaPath.push_back(curLambda);
// Time to stop for LASSO?
if (lasso)
{
if (curLambda <= lambda1)
{
InterpolateBeta();
break;
}
}
}
// Unfortunate copy...
beta = betaPath.back();
return ComputeError(matX, y, !transposeData);
}
inline double LARS::Train(const arma::mat& data,
const arma::rowvec& responses,
const bool transposeData)
{
arma::vec beta;
return Train(data, responses, beta, transposeData);
}
inline void LARS::Predict(const arma::mat& points,
arma::rowvec& predictions,
const bool rowMajor) const
{
// We really only need to store beta internally...
if (rowMajor)
predictions = trans(points * betaPath.back());
else
predictions = betaPath.back().t() * points;
}
// Private functions.
inline void LARS::Deactivate(const size_t activeVarInd)
{
isActive[activeSet[activeVarInd]] = false;
activeSet.erase(activeSet.begin() + activeVarInd);
}
inline void LARS::Activate(const size_t varInd)
{
isActive[varInd] = true;
activeSet.push_back(varInd);
}
inline void LARS::Ignore(const size_t varInd)
{
isIgnored[varInd] = true;
ignoreSet.push_back(varInd);
}
inline void LARS::ComputeYHatDirection(const arma::mat& matX,
const arma::vec& betaDirection,
arma::vec& yHatDirection)
{
yHatDirection.fill(0);
for (size_t i = 0; i < activeSet.size(); ++i)
yHatDirection += betaDirection(i) * matX.col(activeSet[i]);
}
inline void LARS::InterpolateBeta()
{
int pathLength = betaPath.size();
// interpolate beta and stop
double ultimateLambda = lambdaPath[pathLength - 1];
double penultimateLambda = lambdaPath[pathLength - 2];
double interp = (penultimateLambda - lambda1)
/ (penultimateLambda - ultimateLambda);
betaPath[pathLength - 1] = (1 - interp) * (betaPath[pathLength - 2])
+ interp * betaPath[pathLength - 1];
lambdaPath[pathLength - 1] = lambda1;
}
inline void LARS::CholeskyInsert(const arma::vec& newX,
const arma::mat& X)
{
if (matUtriCholFactor.n_rows == 0)
{
matUtriCholFactor = arma::mat(1, 1);
if (elasticNet)
matUtriCholFactor(0, 0) = sqrt(dot(newX, newX) + lambda2);
else
matUtriCholFactor(0, 0) = norm(newX, 2);
}
else
{
arma::vec newGramCol = trans(X) * newX;
CholeskyInsert(dot(newX, newX), newGramCol);
}
}
inline void LARS::CholeskyInsert(double sqNormNewX,
const arma::vec& newGramCol)
{
int n = matUtriCholFactor.n_rows;
if (n == 0)
{
matUtriCholFactor = arma::mat(1, 1);
if (elasticNet)
matUtriCholFactor(0, 0) = sqrt(sqNormNewX + lambda2);
else
matUtriCholFactor(0, 0) = sqrt(sqNormNewX);
}
else
{
arma::mat matNewR = arma::mat(n + 1, n + 1);
if (elasticNet)
sqNormNewX += lambda2;
arma::vec matUtriCholFactork = solve(trimatl(trans(matUtriCholFactor)),
newGramCol);
matNewR(arma::span(0, n - 1), arma::span(0, n - 1)) = matUtriCholFactor;
matNewR(arma::span(0, n - 1), n) = matUtriCholFactork;
matNewR(n, arma::span(0, n - 1)).fill(0.0);
matNewR(n, n) = sqrt(sqNormNewX - dot(matUtriCholFactork,
matUtriCholFactork));
matUtriCholFactor = matNewR;
}
}
inline void LARS::GivensRotate(const arma::vec::fixed<2>& x,
arma::vec::fixed<2>& rotatedX,
arma::mat& matG)
{
if (x(1) == 0)
{
matG = arma::eye(2, 2);
rotatedX = x;
}
else
{
double r = norm(x, 2);
matG = arma::mat(2, 2);
double scaledX1 = x(0) / r;
double scaledX2 = x(1) / r;
matG(0, 0) = scaledX1;
matG(1, 0) = -scaledX2;
matG(0, 1) = scaledX2;
matG(1, 1) = scaledX1;
rotatedX = arma::vec(2);
rotatedX(0) = r;
rotatedX(1) = 0;
}
}
inline void LARS::CholeskyDelete(const size_t colToKill)
{
size_t n = matUtriCholFactor.n_rows;
if (colToKill == (n - 1))
{
matUtriCholFactor = matUtriCholFactor(arma::span(0, n - 2),
arma::span(0, n - 2));
}
else
{
matUtriCholFactor.shed_col(colToKill); // remove column colToKill
n--;
for (size_t k = colToKill; k < n; ++k)
{
arma::mat matG;
arma::vec::fixed<2> rotatedVec;
GivensRotate(matUtriCholFactor(arma::span(k, k + 1), k), rotatedVec,
matG);
matUtriCholFactor(arma::span(k, k + 1), k) = rotatedVec;
if (k < n - 1)
{
matUtriCholFactor(arma::span(k, k + 1), arma::span(k + 1, n - 1)) =
matG * matUtriCholFactor(arma::span(k, k + 1),
arma::span(k + 1, n - 1));
}
}
matUtriCholFactor.shed_row(n);
}
}
inline double LARS::ComputeError(const arma::mat& matX,
const arma::rowvec& y,
const bool rowMajor)
{
if (rowMajor)
{
return arma::accu(arma::pow(y - trans(matX * betaPath.back()), 2.0));
}
else
{
return arma::accu(arma::pow(y - betaPath.back().t() * matX, 2.0));
}
}
/**
* Serialize the LARS model.
*/
@@ -3,7 +3,7 @@
# Do not include test programs here
set(SOURCES
linear_regression.hpp
linear_regression.cpp
linear_regression_impl.hpp
)
# add directory name to sources
@@ -14,6 +14,8 @@
#define MLPACK_METHODS_LINEAR_REGRESSION_LINEAR_REGRESSION_HPP
#include <mlpack/prereqs.hpp>
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/util/size_checks.hpp>
namespace mlpack {
namespace regression /** Regression methods. */ {
@@ -167,4 +169,7 @@ class LinearRegression
} // namespace regression
} // namespace mlpack
// Include implementation.
#include "linear_regression_impl.hpp"
#endif // MLPACK_METHODS_LINEAR_REGRESSION_HPP
@@ -1,5 +1,5 @@
/**
* @file methods/linear_regression/linear_regression.cpp
* @file methods/linear_regression/linear_regression_impl.hpp
* @author James Cline
* @author Michael Fox
*
@@ -10,42 +10,45 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_LINEAR_REGRESSION_LINEAR_REGRESSION_IMPL_HPP
#define MLPACK_METHODS_LINEAR_REGRESSION_LINEAR_REGRESSION_IMPL_HPP
#include "linear_regression.hpp"
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/util/size_checks.hpp>
using namespace mlpack;
using namespace mlpack::regression;
namespace mlpack {
namespace regression {
LinearRegression::LinearRegression(const arma::mat& predictors,
const arma::rowvec& responses,
const double lambda,
const bool intercept) :
inline LinearRegression::LinearRegression(
const arma::mat& predictors,
const arma::rowvec& responses,
const double lambda,
const bool intercept) :
LinearRegression(predictors, responses, arma::rowvec(), lambda, intercept)
{}
{ /* Nothing to do. */ }
LinearRegression::LinearRegression(const arma::mat& predictors,
const arma::rowvec& responses,
const arma::rowvec& weights,
const double lambda,
const bool intercept) :
inline LinearRegression::LinearRegression(
const arma::mat& predictors,
const arma::rowvec& responses,
const arma::rowvec& weights,
const double lambda,
const bool intercept) :
lambda(lambda),
intercept(intercept)
{
Train(predictors, responses, weights, intercept);
}
double LinearRegression::Train(const arma::mat& predictors,
const arma::rowvec& responses,
const bool intercept)
inline double LinearRegression::Train(const arma::mat& predictors,
const arma::rowvec& responses,
const bool intercept)
{
return Train(predictors, responses, arma::rowvec(), intercept);
}
double LinearRegression::Train(const arma::mat& predictors,
const arma::rowvec& responses,
const arma::rowvec& weights,
const bool intercept)
inline double LinearRegression::Train(const arma::mat& predictors,
const arma::rowvec& responses,
const arma::rowvec& weights,
const bool intercept)
{
this->intercept = intercept;
@@ -93,7 +96,8 @@ double LinearRegression::Train(const arma::mat& predictors,
return ComputeError(predictors, responses);
}
void LinearRegression::Predict(const arma::mat& points,
inline void LinearRegression::Predict(
const arma::mat& points,
arma::rowvec& predictions) const
{
if (intercept)
@@ -122,8 +126,9 @@ void LinearRegression::Predict(const arma::mat& points,
}
}
double LinearRegression::ComputeError(const arma::mat& predictors,
const arma::rowvec& responses) const
inline double LinearRegression::ComputeError(
const arma::mat& predictors,
const arma::rowvec& responses) const
{
// Sanity check on data.
util::CheckSameSizes(predictors, responses, "LinearRegression::Train()");
@@ -160,3 +165,8 @@ double LinearRegression::ComputeError(const arma::mat& predictors,
return cost;
}
} // namespace regression
} // namespace mlpack
#endif
@@ -5,7 +5,6 @@
# that you have files in both sections
set(SOURCES
lcc.hpp
lcc.cpp
lcc_impl.hpp
)
@@ -1,247 +0,0 @@
/**
* @file methods/local_coordinate_coding/lcc.cpp
* @author Nishant Mehta
*
* Implementation of Local Coordinate Coding.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "lcc.hpp"
#include <mlpack/core/math/lin_alg.hpp>
namespace mlpack {
namespace lcc {
LocalCoordinateCoding::LocalCoordinateCoding(
const size_t atoms,
const double lambda,
const size_t maxIterations,
const double tolerance) :
atoms(atoms),
lambda(lambda),
maxIterations(maxIterations),
tolerance(tolerance)
{
// Nothing to do.
}
void LocalCoordinateCoding::Encode(const arma::mat& data, arma::mat& codes)
{
arma::mat invSqDists = 1.0 / (repmat(trans(sum(square(dictionary))), 1,
data.n_cols) + repmat(sum(square(data)), atoms, 1) - 2 * trans(dictionary)
* data);
arma::mat dictGram = trans(dictionary) * dictionary;
arma::mat dictGramTD(dictGram.n_rows, dictGram.n_cols);
codes.set_size(atoms, data.n_cols);
for (size_t i = 0; i < data.n_cols; ++i)
{
// Report progress.
if ((i % 100) == 0)
{
Log::Debug << "Optimization at point " << i << "." << std::endl;
}
arma::vec invW = invSqDists.unsafe_col(i);
arma::mat dictPrime = dictionary * diagmat(invW);
arma::mat dictGramTD = diagmat(invW) * dictGram * diagmat(invW);
bool useCholesky = false;
regression::LARS lars(useCholesky, dictGramTD, 0.5 * lambda);
// Run LARS for this point, by making an alias of the point and passing
// that.
arma::vec beta = codes.unsafe_col(i);
arma::rowvec responses = data.unsafe_col(i).t();
lars.Train(dictPrime, responses, beta, false);
beta %= invW; // Remember, beta is an alias of codes.col(i).
}
}
void LocalCoordinateCoding::OptimizeDictionary(const arma::mat& data,
const arma::mat& codes,
const arma::uvec& adjacencies)
{
// Count number of atomic neighbors for each point x^i.
arma::uvec neighborCounts = arma::zeros<arma::uvec>(data.n_cols, 1);
if (adjacencies.n_elem > 0)
{
// This gets the column index. Intentional integer division.
size_t curPointInd = (size_t) (adjacencies(0) / atoms);
++neighborCounts(curPointInd);
size_t nextColIndex = (curPointInd + 1) * atoms;
for (size_t l = 1; l < adjacencies.n_elem; l++)
{
// If l no longer refers to an element in this column, advance the column
// number accordingly.
if (adjacencies(l) >= nextColIndex)
{
curPointInd = (size_t) (adjacencies(l) / atoms);
nextColIndex = (curPointInd + 1) * atoms;
}
++neighborCounts(curPointInd);
}
}
// Build dataPrime := [X x^1 ... x^1 ... x^n ... x^n]
// where each x^i is repeated for the number of neighbors x^i has.
arma::mat dataPrime = arma::zeros(data.n_rows,
data.n_cols + adjacencies.n_elem);
dataPrime(arma::span::all, arma::span(0, data.n_cols - 1)) = data;
size_t curCol = data.n_cols;
for (size_t i = 0; i < data.n_cols; ++i)
{
if (neighborCounts(i) > 0)
{
dataPrime(arma::span::all, arma::span(curCol, curCol + neighborCounts(i)
- 1)) = repmat(data.col(i), 1, neighborCounts(i));
}
curCol += neighborCounts(i);
}
// Handle the case of inactive atoms (atoms not used in the given coding).
std::vector<size_t> inactiveAtoms;
for (size_t j = 0; j < atoms; ++j)
if (accu(codes.row(j) != 0) == 0)
inactiveAtoms.push_back(j);
const size_t nInactiveAtoms = inactiveAtoms.size();
const size_t nActiveAtoms = atoms - nInactiveAtoms;
// Efficient construction of codes restricted to active atoms.
arma::mat codesPrime = arma::zeros(nActiveAtoms, data.n_cols +
adjacencies.n_elem);
arma::vec wSquared = arma::ones(data.n_cols + adjacencies.n_elem, 1);
if (nInactiveAtoms > 0)
{
Log::Warn << "There are " << nInactiveAtoms
<< " inactive atoms. They will be re-initialized randomly.\n";
// Create matrix holding only active codes.
arma::mat activeCodes;
math::RemoveRows(codes, inactiveAtoms, activeCodes);
// Create reverse atom lookup for active atoms.
arma::uvec atomReverseLookup(atoms);
size_t inactiveOffset = 0;
for (size_t i = 0; i < atoms; ++i)
{
if (inactiveAtoms[inactiveOffset] == i)
++inactiveOffset;
else
atomReverseLookup(i - inactiveOffset) = i;
}
codesPrime(arma::span::all, arma::span(0, data.n_cols - 1)) = activeCodes;
// Fill the rest of codesPrime.
for (size_t l = 0; l < adjacencies.n_elem; ++l)
{
// Recover the location in the codes matrix that this adjacency refers to.
size_t atomInd = adjacencies(l) % atoms;
size_t pointInd = (size_t) (adjacencies(l) / atoms);
// Fill matrix.
codesPrime(atomReverseLookup(atomInd), data.n_cols + l) = 1.0;
wSquared(data.n_cols + l) = codes(atomInd, pointInd);
}
}
else
{
// All atoms are active.
codesPrime(arma::span::all, arma::span(0, data.n_cols - 1)) = codes;
for (size_t l = 0; l < adjacencies.n_elem; ++l)
{
// Recover the location in the codes matrix that this adjacency refers to.
size_t atomInd = adjacencies(l) % atoms;
size_t pointInd = (size_t) (adjacencies(l) / atoms);
// Fill matrix.
codesPrime(atomInd, data.n_cols + l) = 1.0;
wSquared(data.n_cols + l) = codes(atomInd, pointInd);
}
}
wSquared.subvec(data.n_cols, wSquared.n_elem - 1) = lambda *
abs(wSquared.subvec(data.n_cols, wSquared.n_elem - 1));
// Solve system.
if (nInactiveAtoms == 0)
{
// No inactive atoms. We can solve directly.
arma::mat A = codesPrime * diagmat(wSquared) * trans(codesPrime);
arma::mat B = codesPrime * diagmat(wSquared) * trans(dataPrime);
dictionary = trans(solve(A, B));
/*
dictionary = trans(solve(codesPrime * diagmat(wSquared) * trans(codesPrime),
codesPrime * diagmat(wSquared) * trans(dataPrime)));
*/
}
else
{
// Inactive atoms must be reinitialized randomly, so we cannot solve
// directly for the entire dictionary estimate.
arma::mat dictionaryActive =
trans(solve(codesPrime * diagmat(wSquared) * trans(codesPrime),
codesPrime * diagmat(wSquared) * trans(dataPrime)));
// Update all atoms.
size_t currentInactiveIndex = 0;
for (size_t i = 0; i < atoms; ++i)
{
if (inactiveAtoms[currentInactiveIndex] == i)
{
// This atom is inactive. Reinitialize it randomly.
dictionary.col(i) = (data.col(math::RandInt(data.n_cols)) +
data.col(math::RandInt(data.n_cols)) +
data.col(math::RandInt(data.n_cols)));
// Now normalize the atom.
dictionary.col(i) /= norm(dictionary.col(i), 2);
// Increment inactive atom counter.
++currentInactiveIndex;
}
else
{
// Update estimate.
dictionary.col(i) = dictionaryActive.col(i - currentInactiveIndex);
}
}
}
}
double LocalCoordinateCoding::Objective(const arma::mat& data,
const arma::mat& codes,
const arma::uvec& adjacencies) const
{
double weightedL1NormZ = 0;
for (size_t l = 0; l < adjacencies.n_elem; l++)
{
// Map adjacency back to its location in the codes matrix.
const size_t atomInd = adjacencies(l) % atoms;
const size_t pointInd = (size_t) (adjacencies(l) / atoms);
weightedL1NormZ += fabs(codes(atomInd, pointInd)) * arma::as_scalar(
arma::sum(arma::square(dictionary.col(atomInd) - data.col(pointInd))));
}
double froNormResidual = norm(data - dictionary * codes, "fro");
return std::pow(froNormResidual, 2.0) + lambda * weightedL1NormZ;
}
} // namespace lcc
} // namespace mlpack
@@ -15,11 +15,12 @@
#include <mlpack/prereqs.hpp>
#include <mlpack/methods/lars/lars.hpp>
#include <mlpack/core/math/lin_alg.hpp>
// Include three simple dictionary initializers from sparse coding.
#include "../sparse_coding/nothing_initializer.hpp"
#include "../sparse_coding/data_dependent_random_initializer.hpp"
#include "../sparse_coding/random_initializer.hpp"
#include <mlpack/methods/sparse_coding/nothing_initializer.hpp>
#include <mlpack/methods/sparse_coding/data_dependent_random_initializer.hpp>
#include <mlpack/methods/sparse_coding/random_initializer.hpp>
namespace mlpack {
namespace lcc {
@@ -35,6 +35,19 @@ LocalCoordinateCoding::LocalCoordinateCoding(
Train(data, initializer);
}
inline LocalCoordinateCoding::LocalCoordinateCoding(
const size_t atoms,
const double lambda,
const size_t maxIterations,
const double tolerance) :
atoms(atoms),
lambda(lambda),
maxIterations(maxIterations),
tolerance(tolerance)
{
// Nothing to do.
}
template<typename DictionaryInitializer>
double LocalCoordinateCoding::Train(
const arma::mat& data,
@@ -103,6 +116,224 @@ double LocalCoordinateCoding::Train(
return lastObjVal;
}
inline void LocalCoordinateCoding::Encode(const arma::mat& data,
arma::mat& codes)
{
arma::mat invSqDists = 1.0 / (repmat(trans(sum(square(dictionary))), 1,
data.n_cols) + repmat(sum(square(data)), atoms, 1) - 2 * trans(dictionary)
* data);
arma::mat dictGram = trans(dictionary) * dictionary;
arma::mat dictGramTD(dictGram.n_rows, dictGram.n_cols);
codes.set_size(atoms, data.n_cols);
for (size_t i = 0; i < data.n_cols; ++i)
{
// Report progress.
if ((i % 100) == 0)
{
Log::Debug << "Optimization at point " << i << "." << std::endl;
}
arma::vec invW = invSqDists.unsafe_col(i);
arma::mat dictPrime = dictionary * diagmat(invW);
arma::mat dictGramTD = diagmat(invW) * dictGram * diagmat(invW);
bool useCholesky = false;
regression::LARS lars(useCholesky, dictGramTD, 0.5 * lambda);
// Run LARS for this point, by making an alias of the point and passing
// that.
arma::vec beta = codes.unsafe_col(i);
arma::rowvec responses = data.unsafe_col(i).t();
lars.Train(dictPrime, responses, beta, false);
beta %= invW; // Remember, beta is an alias of codes.col(i).
}
}
inline void LocalCoordinateCoding::OptimizeDictionary(
const arma::mat& data,
const arma::mat& codes,
const arma::uvec& adjacencies)
{
// Count number of atomic neighbors for each point x^i.
arma::uvec neighborCounts = arma::zeros<arma::uvec>(data.n_cols, 1);
if (adjacencies.n_elem > 0)
{
// This gets the column index. Intentional integer division.
size_t curPointInd = (size_t) (adjacencies(0) / atoms);
++neighborCounts(curPointInd);
size_t nextColIndex = (curPointInd + 1) * atoms;
for (size_t l = 1; l < adjacencies.n_elem; l++)
{
// If l no longer refers to an element in this column, advance the column
// number accordingly.
if (adjacencies(l) >= nextColIndex)
{
curPointInd = (size_t) (adjacencies(l) / atoms);
nextColIndex = (curPointInd + 1) * atoms;
}
++neighborCounts(curPointInd);
}
}
// Build dataPrime := [X x^1 ... x^1 ... x^n ... x^n]
// where each x^i is repeated for the number of neighbors x^i has.
arma::mat dataPrime = arma::zeros(data.n_rows,
data.n_cols + adjacencies.n_elem);
dataPrime(arma::span::all, arma::span(0, data.n_cols - 1)) = data;
size_t curCol = data.n_cols;
for (size_t i = 0; i < data.n_cols; ++i)
{
if (neighborCounts(i) > 0)
{
dataPrime(arma::span::all, arma::span(curCol, curCol + neighborCounts(i)
- 1)) = repmat(data.col(i), 1, neighborCounts(i));
}
curCol += neighborCounts(i);
}
// Handle the case of inactive atoms (atoms not used in the given coding).
std::vector<size_t> inactiveAtoms;
for (size_t j = 0; j < atoms; ++j)
if (accu(codes.row(j) != 0) == 0)
inactiveAtoms.push_back(j);
const size_t nInactiveAtoms = inactiveAtoms.size();
const size_t nActiveAtoms = atoms - nInactiveAtoms;
// Efficient construction of codes restricted to active atoms.
arma::mat codesPrime = arma::zeros(nActiveAtoms, data.n_cols +
adjacencies.n_elem);
arma::vec wSquared = arma::ones(data.n_cols + adjacencies.n_elem, 1);
if (nInactiveAtoms > 0)
{
Log::Warn << "There are " << nInactiveAtoms
<< " inactive atoms. They will be re-initialized randomly.\n";
// Create matrix holding only active codes.
arma::mat activeCodes;
math::RemoveRows(codes, inactiveAtoms, activeCodes);
// Create reverse atom lookup for active atoms.
arma::uvec atomReverseLookup(atoms);
size_t inactiveOffset = 0;
for (size_t i = 0; i < atoms; ++i)
{
if (inactiveAtoms[inactiveOffset] == i)
++inactiveOffset;
else
atomReverseLookup(i - inactiveOffset) = i;
}
codesPrime(arma::span::all, arma::span(0, data.n_cols - 1)) = activeCodes;
// Fill the rest of codesPrime.
for (size_t l = 0; l < adjacencies.n_elem; ++l)
{
// Recover the location in the codes matrix that this adjacency refers to.
size_t atomInd = adjacencies(l) % atoms;
size_t pointInd = (size_t) (adjacencies(l) / atoms);
// Fill matrix.
codesPrime(atomReverseLookup(atomInd), data.n_cols + l) = 1.0;
wSquared(data.n_cols + l) = codes(atomInd, pointInd);
}
}
else
{
// All atoms are active.
codesPrime(arma::span::all, arma::span(0, data.n_cols - 1)) = codes;
for (size_t l = 0; l < adjacencies.n_elem; ++l)
{
// Recover the location in the codes matrix that this adjacency refers to.
size_t atomInd = adjacencies(l) % atoms;
size_t pointInd = (size_t) (adjacencies(l) / atoms);
// Fill matrix.
codesPrime(atomInd, data.n_cols + l) = 1.0;
wSquared(data.n_cols + l) = codes(atomInd, pointInd);
}
}
wSquared.subvec(data.n_cols, wSquared.n_elem - 1) = lambda *
abs(wSquared.subvec(data.n_cols, wSquared.n_elem - 1));
// Solve system.
if (nInactiveAtoms == 0)
{
// No inactive atoms. We can solve directly.
arma::mat A = codesPrime * diagmat(wSquared) * trans(codesPrime);
arma::mat B = codesPrime * diagmat(wSquared) * trans(dataPrime);
dictionary = trans(solve(A, B));
/*
dictionary = trans(solve(codesPrime * diagmat(wSquared) * trans(codesPrime),
codesPrime * diagmat(wSquared) * trans(dataPrime)));
*/
}
else
{
// Inactive atoms must be reinitialized randomly, so we cannot solve
// directly for the entire dictionary estimate.
arma::mat dictionaryActive =
trans(solve(codesPrime * diagmat(wSquared) * trans(codesPrime),
codesPrime * diagmat(wSquared) * trans(dataPrime)));
// Update all atoms.
size_t currentInactiveIndex = 0;
for (size_t i = 0; i < atoms; ++i)
{
if (inactiveAtoms[currentInactiveIndex] == i)
{
// This atom is inactive. Reinitialize it randomly.
dictionary.col(i) = (data.col(math::RandInt(data.n_cols)) +
data.col(math::RandInt(data.n_cols)) +
data.col(math::RandInt(data.n_cols)));
// Now normalize the atom.
dictionary.col(i) /= norm(dictionary.col(i), 2);
// Increment inactive atom counter.
++currentInactiveIndex;
}
else
{
// Update estimate.
dictionary.col(i) = dictionaryActive.col(i - currentInactiveIndex);
}
}
}
}
inline double LocalCoordinateCoding::Objective(
const arma::mat& data,
const arma::mat& codes,
const arma::uvec& adjacencies) const
{
double weightedL1NormZ = 0;
for (size_t l = 0; l < adjacencies.n_elem; l++)
{
// Map adjacency back to its location in the codes matrix.
const size_t atomInd = adjacencies(l) % atoms;
const size_t pointInd = (size_t) (adjacencies(l) / atoms);
weightedL1NormZ += fabs(codes(atomInd, pointInd)) * arma::as_scalar(
arma::sum(arma::square(dictionary.col(atomInd) - data.col(pointInd))));
}
double froNormResidual = norm(data - dictionary * codes, "fro");
return std::pow(froNormResidual, 2.0) + lambda * weightedL1NormZ;
}
template<typename Archive>
void LocalCoordinateCoding::serialize(Archive& ar,
const uint32_t /* version */)
@@ -2,7 +2,7 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
matrix_completion.hpp
matrix_completion.cpp
matrix_completion_impl.hpp
)
# Add directory name to sources.
@@ -145,4 +145,7 @@ class MatrixCompletion
} // namespace matrix_completion
} // namespace mlpack
// Include implementation.
#include "matrix_completion_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/matrix_completion/matrix_completion.cpp
* @file methods/matrix_completion/matrix_completion_impl.hpp
* @author Stephen Tu
*
* Implementation of MatrixCompletion class.
@@ -9,6 +9,8 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_MATRIX_COMPLETION_MATRIX_COMPLETION_IMPL_HPP
#define MLPACK_METHODS_MATRIX_COMPLETION_MATRIX_COMPLETION_IMPL_HPP
#include "matrix_completion.hpp"
#include <mlpack/core/util/size_checks.hpp>
@@ -16,35 +18,47 @@
namespace mlpack {
namespace matrix_completion {
MatrixCompletion::MatrixCompletion(const size_t m,
const size_t n,
const arma::umat& indices,
const arma::vec& values,
const size_t r) :
m(m), n(n), indices(indices), values(values),
inline MatrixCompletion::MatrixCompletion(
const size_t m,
const size_t n,
const arma::umat& indices,
const arma::vec& values,
const size_t r) :
m(m),
n(n),
indices(indices),
values(values),
sdp(indices.n_cols, 0, arma::randu<arma::mat>(m + n, r))
{
CheckValues();
InitSDP();
}
MatrixCompletion::MatrixCompletion(const size_t m,
const size_t n,
const arma::umat& indices,
const arma::vec& values,
const arma::mat& initialPoint) :
m(m), n(n), indices(indices), values(values),
inline MatrixCompletion::MatrixCompletion(
const size_t m,
const size_t n,
const arma::umat& indices,
const arma::vec& values,
const arma::mat& initialPoint) :
m(m),
n(n),
indices(indices),
values(values),
sdp(indices.n_cols, 0, initialPoint)
{
CheckValues();
InitSDP();
}
MatrixCompletion::MatrixCompletion(const size_t m,
const size_t n,
const arma::umat& indices,
const arma::vec& values) :
m(m), n(n), indices(indices), values(values),
inline MatrixCompletion::MatrixCompletion(
const size_t m,
const size_t n,
const arma::umat& indices,
const arma::vec& values) :
m(m),
n(n),
indices(indices),
values(values),
sdp(indices.n_cols, 0,
arma::randu<arma::mat>(m + n, DefaultRank(m, n, indices.n_cols)))
{
@@ -52,7 +66,7 @@ MatrixCompletion::MatrixCompletion(const size_t m,
InitSDP();
}
void MatrixCompletion::CheckValues()
inline void MatrixCompletion::CheckValues()
{
if (indices.n_rows != 2)
{
@@ -73,7 +87,7 @@ void MatrixCompletion::CheckValues()
}
}
void MatrixCompletion::InitSDP()
inline void MatrixCompletion::InitSDP()
{
sdp.SDP().C().eye(m + n, m + n);
sdp.SDP().SparseB() = 2. * values;
@@ -86,7 +100,7 @@ void MatrixCompletion::InitSDP()
}
}
void MatrixCompletion::Recover(arma::mat& recovered)
inline void MatrixCompletion::Recover(arma::mat& recovered)
{
recovered = sdp.Function().GetInitialPoint();
sdp.Optimize(recovered);
@@ -94,9 +108,9 @@ void MatrixCompletion::Recover(arma::mat& recovered)
recovered = recovered(arma::span(0, m - 1), arma::span(m, m + n - 1));
}
size_t MatrixCompletion::DefaultRank(const size_t m,
const size_t n,
const size_t p)
inline size_t MatrixCompletion::DefaultRank(const size_t m,
const size_t n,
const size_t p)
{
// If r = O(sqrt(p)), then we are guaranteed an exact solution.
// For more details, see
@@ -114,3 +128,5 @@ size_t MatrixCompletion::DefaultRank(const size_t m,
} // namespace matrix_completion
} // namespace mlpack
#endif
@@ -14,7 +14,7 @@ set(SOURCES
sort_policies/furthest_neighbor_sort_impl.hpp
typedef.hpp
unmap.hpp
unmap.cpp
unmap_impl.hpp
)
# Add directory name to sources.
@@ -63,4 +63,7 @@ void Unmap(const arma::Mat<size_t>& neighbors,
} // namespace neighbor
} // namespace mlpack
// Include implementation.
#include "unmap_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/neighbor_search/unmap.cpp
* @file methods/neighbor_search/unmap_impl.hpp
* @author Ryan Curtin
*
* Auxiliary function to unmap neighbor search results.
@@ -9,19 +9,22 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_NEIGHBOR_SEARCH_UNMAP_IMPL_HPP
#define MLPACK_METHODS_NEIGHBOR_SEARCH_UNMAP_IMPL_HPP
#include "unmap.hpp"
namespace mlpack {
namespace neighbor {
// Useful in the dual-tree setting.
void Unmap(const arma::Mat<size_t>& neighbors,
const arma::mat& distances,
const std::vector<size_t>& referenceMap,
const std::vector<size_t>& queryMap,
arma::Mat<size_t>& neighborsOut,
arma::mat& distancesOut,
const bool squareRoot)
inline void Unmap(const arma::Mat<size_t>& neighbors,
const arma::mat& distances,
const std::vector<size_t>& referenceMap,
const std::vector<size_t>& queryMap,
arma::Mat<size_t>& neighborsOut,
arma::mat& distancesOut,
const bool squareRoot)
{
// Set matrices to correct size.
neighborsOut.set_size(neighbors.n_rows, neighbors.n_cols);
@@ -44,12 +47,12 @@ void Unmap(const arma::Mat<size_t>& neighbors,
}
// Useful in the single-tree setting.
void Unmap(const arma::Mat<size_t>& neighbors,
const arma::mat& distances,
const std::vector<size_t>& referenceMap,
arma::Mat<size_t>& neighborsOut,
arma::mat& distancesOut,
const bool squareRoot)
inline void Unmap(const arma::Mat<size_t>& neighbors,
const arma::mat& distances,
const std::vector<size_t>& referenceMap,
arma::Mat<size_t>& neighborsOut,
arma::mat& distancesOut,
const bool squareRoot)
{
// Set matrices to correct size.
neighborsOut.set_size(neighbors.n_rows, neighbors.n_cols);
@@ -67,3 +70,5 @@ void Unmap(const arma::Mat<size_t>& neighbors,
} // namespace neighbor
} // namespace mlpack
#endif
+1 -1
View File
@@ -2,7 +2,7 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
quic_svd.hpp
quic_svd.cpp
quic_svd_impl.hpp
)
# Add directory name to sources.
+3
View File
@@ -94,4 +94,7 @@ class QUIC_SVD
} // namespace svd
} // namespace mlpack
// Include implementation.
#include "quic_svd_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/quic_svd/quic_svd.cpp
* @file methods/quic_svd/quic_svd_impl.hpp
* @author Siddharth Agrawal
*
* An implementation of QUIC-SVD.
@@ -9,30 +9,31 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_QUIC_SVD_QUIC_SVD_IMPL_HPP
#define MLPACK_METHODS_QUIC_SVD_QUIC_SVD_IMPL_HPP
// In case it hasn't been included yet.
#include "quic_svd.hpp"
using namespace mlpack::tree;
namespace mlpack {
namespace svd {
QUIC_SVD::QUIC_SVD(const arma::mat& dataset,
arma::mat& u,
arma::mat& v,
arma::mat& sigma,
const double epsilon,
const double delta) :
inline QUIC_SVD::QUIC_SVD(
const arma::mat& dataset,
arma::mat& u,
arma::mat& v,
arma::mat& sigma,
const double epsilon,
const double delta) :
dataset(dataset)
{
// Since columns are sample in the implementation, the matrix is transposed if
// necessary for maximum speedup.
CosineTree* ctree;
tree::CosineTree* ctree;
if (dataset.n_cols > dataset.n_rows)
ctree = new CosineTree(dataset, epsilon, delta);
ctree = new tree::CosineTree(dataset, epsilon, delta);
else
ctree = new CosineTree(dataset.t(), epsilon, delta);
ctree = new tree::CosineTree(dataset.t(), epsilon, delta);
// Get subspace basis by creating the cosine tree.
ctree->GetFinalBasis(basis);
@@ -45,9 +46,9 @@ QUIC_SVD::QUIC_SVD(const arma::mat& dataset,
ExtractSVD(u, v, sigma);
}
void QUIC_SVD::ExtractSVD(arma::mat& u,
arma::mat& v,
arma::mat& sigma)
inline void QUIC_SVD::ExtractSVD(arma::mat& u,
arma::mat& v,
arma::mat& sigma)
{
// Calculate A * V_hat, necessary for further calculations.
arma::mat projectedMat;
@@ -82,3 +83,5 @@ void QUIC_SVD::ExtractSVD(arma::mat& u,
} // namespace svd
} // namespace mlpack
#endif
+1 -1
View File
@@ -2,7 +2,7 @@
# Anything not in this list will not be compiled into the output library
set(SOURCES
radical.hpp
radical.cpp
radical_impl.hpp
)
# add directory name to sources
+5
View File
@@ -16,6 +16,8 @@
#include <mlpack/prereqs.hpp>
#include <mlpack/core/util/io.hpp>
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/util/timers.hpp>
namespace mlpack {
namespace radical {
@@ -148,4 +150,7 @@ void WhitenFeatureMajorMatrix(const arma::mat& matX,
} // namespace radical
} // namespace mlpack
// Include implementation.
#include "radical_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/radical/radical.cpp
* @file methods/radical/radical_impl.hpp
* @author Nishant Mehta
*
* Implementation of Radical class
@@ -9,22 +9,21 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_RADICAL_RADICAL_IMPL_HPP
#define MLPACK_METHODS_RADICAL_RADICAL_IMPL_HPP
#include "radical.hpp"
#include <mlpack/core/util/log.hpp>
#include <mlpack/core/util/timers.hpp>
using namespace std;
using namespace arma;
using namespace mlpack;
using namespace mlpack::radical;
namespace mlpack {
namespace radical {
// Set the parameters to RADICAL.
Radical::Radical(const double noiseStdDev,
const size_t replicates,
const size_t angles,
const size_t sweeps,
const size_t m) :
inline Radical::Radical(
const double noiseStdDev,
const size_t replicates,
const size_t angles,
const size_t sweeps,
const size_t m) :
noiseStdDev(noiseStdDev),
replicates(replicates),
angles(angles),
@@ -34,14 +33,15 @@ Radical::Radical(const double noiseStdDev,
// Nothing to do here.
}
void Radical::CopyAndPerturb(mat& xNew, const mat& x) const
inline void Radical::CopyAndPerturb(arma::mat& xNew,
const arma::mat& x) const
{
xNew = repmat(x, replicates, 1) + noiseStdDev * randn(replicates * x.n_rows,
xNew = arma::repmat(x, replicates, 1) + noiseStdDev * arma::randn(replicates * x.n_rows,
x.n_cols);
}
double Radical::Vasicek(vec& z) const
inline double Radical::Vasicek(arma::vec& z) const
{
z = sort(z);
@@ -54,25 +54,26 @@ double Radical::Vasicek(vec& z) const
// Apparently faster.
double sum = 0;
uword range = z.n_elem - m;
for (uword i = 0; i < range; ++i)
arma::uword range = z.n_elem - m;
for (arma::uword i = 0; i < range; ++i)
{
sum += log(max(z(i + m) - z(i), DBL_MIN));
sum += log(std::max(z(i + m) - z(i), DBL_MIN));
}
return sum;
}
double Radical::DoRadical2D(const mat& matX, util::Timers& timers)
inline double Radical::DoRadical2D(const arma::mat& matX,
util::Timers& timers)
{
timers.Start("radical_copy_and_perturb");
CopyAndPerturb(perturbed, matX);
timers.Stop("radical_copy_and_perturb");
mat::fixed<2, 2> matJacobi;
arma::mat::fixed<2, 2> matJacobi;
vec values(angles);
arma::vec values(angles);
for (size_t i = 0; i < angles; ++i)
{
@@ -86,29 +87,29 @@ double Radical::DoRadical2D(const mat& matX, util::Timers& timers)
matJacobi(1, 1) = cosTheta;
candidate = perturbed * matJacobi;
vec candidateY1 = candidate.unsafe_col(0);
vec candidateY2 = candidate.unsafe_col(1);
arma::vec candidateY1 = candidate.unsafe_col(0);
arma::vec candidateY2 = candidate.unsafe_col(1);
values(i) = Vasicek(candidateY1) + Vasicek(candidateY2);
}
uword indOpt = 0;
arma::uword indOpt = 0;
values.min(indOpt); // we ignore the return value; we don't care about it
return (indOpt / (double) angles) * M_PI / 2.0;
}
void Radical::DoRadical(const mat& matXT,
mat& matY,
mat& matW,
util::Timers& timers)
inline void Radical::DoRadical(const arma::mat& matXT,
arma::mat& matY,
arma::mat& matW,
util::Timers& timers)
{
// matX is nPoints by nDims (although less intuitive than columns being
// points, and although this is the transpose of the ICA literature, this
// choice is for computational efficiency when repeatedly generating
// two-dimensional coordinate projections for Radical2D).
timers.Start("radical_transpose_data");
mat matX = trans(matXT);
arma::mat matX = trans(matXT);
timers.Stop("radical_transpose_data");
// If m was not specified, initialize m as recommended in
@@ -120,8 +121,8 @@ void Radical::DoRadical(const mat& matXT,
const size_t nPoints = matX.n_rows;
timers.Start("radical_whiten_data");
mat matXWhitened;
mat matWhitening;
arma::mat matXWhitened;
arma::mat matWhitening;
WhitenFeatureMajorMatrix(matX, matY, matWhitening);
timers.Stop("radical_whiten_data");
// matY is now the whitened form of matX.
@@ -135,9 +136,9 @@ void Radical::DoRadical(const mat& matXT,
timers.Start("radical_do_radical");
matW = matWhitening;
mat matYSubspace(nPoints, 2);
arma::mat matYSubspace(nPoints, 2);
mat matJ = eye(nDims, nDims);
arma::mat matJ = arma::eye(nDims, nDims);
for (size_t sweepNum = 0; sweepNum < sweeps; sweepNum++)
{
@@ -184,13 +185,18 @@ void Radical::DoRadical(const mat& matXT,
timers.Stop("radical_transpose_data");
}
void mlpack::radical::WhitenFeatureMajorMatrix(const mat& matX,
mat& matXWhitened,
mat& matWhitening)
inline void WhitenFeatureMajorMatrix(const arma::mat& matX,
arma::mat& matXWhitened,
arma::mat& matWhitening)
{
mat matU, matV;
vec s;
arma::mat matU, matV;
arma::vec s;
arma::svd(matU, s, matV, cov(matX));
matWhitening = matU * diagmat(1 / sqrt(s)) * trans(matV);
matXWhitened = matX * matWhitening;
}
} // namespace radical
} // namespace mlpack
#endif
@@ -2,7 +2,7 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
randomized_svd.hpp
randomized_svd.cpp
randomized_svd_impl.hpp
)
# Add directory name to sources.
@@ -260,4 +260,7 @@ class RandomizedSVD
} // namespace svd
} // namespace mlpack
// Include implementation.
#include "randomized_svd_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/randomized_svd/randomized_svd.cpp
* @file methods/randomized_svd/randomized_svd_impl.hpp
* @author Marcus Edel
*
* Implementation of the randomized SVD method.
@@ -10,19 +10,23 @@
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_RANDOMIZED_SVD_RANDOMIZED_SVD_IMPL_HPP
#define MLPACK_METHODS_RANDOMIZED_SVD_RANDOMIZED_SVD_IMPL_HPP
#include "randomized_svd.hpp"
namespace mlpack {
namespace svd {
RandomizedSVD::RandomizedSVD(const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t iteratedPower,
const size_t maxIterations,
const size_t rank,
const double eps) :
inline RandomizedSVD::RandomizedSVD(
const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t iteratedPower,
const size_t maxIterations,
const size_t rank,
const double eps) :
iteratedPower(iteratedPower),
maxIterations(maxIterations),
eps(eps)
@@ -37,9 +41,10 @@ RandomizedSVD::RandomizedSVD(const arma::mat& data,
}
}
RandomizedSVD::RandomizedSVD(const size_t iteratedPower,
const size_t maxIterations,
const double eps) :
inline RandomizedSVD::RandomizedSVD(
const size_t iteratedPower,
const size_t maxIterations,
const double eps) :
iteratedPower(iteratedPower),
maxIterations(maxIterations),
eps(eps)
@@ -48,11 +53,11 @@ RandomizedSVD::RandomizedSVD(const size_t iteratedPower,
}
void RandomizedSVD::Apply(const arma::sp_mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t rank)
inline void RandomizedSVD::Apply(const arma::sp_mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t rank)
{
// Center the data into a temporary matrix for sparse matrix.
arma::sp_mat rowMean = arma::sum(data, 1) / data.n_cols;
@@ -60,11 +65,11 @@ void RandomizedSVD::Apply(const arma::sp_mat& data,
Apply(data, u, s, v, rank, rowMean);
}
void RandomizedSVD::Apply(const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t rank)
inline void RandomizedSVD::Apply(const arma::mat& data,
arma::mat& u,
arma::vec& s,
arma::mat& v,
const size_t rank)
{
// Center the data into a temporary matrix.
arma::mat rowMean = arma::sum(data, 1) / data.n_cols + eps;
@@ -74,3 +79,5 @@ void RandomizedSVD::Apply(const arma::mat& data,
} // namespace svd
} // namespace mlpack
#endif
@@ -8,7 +8,6 @@ set(SOURCES
range_search_stat.hpp
rs_model.hpp
rs_model_impl.hpp
rs_model.cpp
)
# Add directory name to sources.
@@ -1,289 +0,0 @@
/**
* @file methods/range_search/rs_model.cpp
* @author Ryan Curtin
*
* Implementation of serialize() and inline functions for RSModel.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "rs_model.hpp"
#include <mlpack/core/math/random_basis.hpp>
namespace mlpack {
namespace range {
/**
* Initialize the RSModel with the given tree type and whether or not a random
* basis should be used.
*/
RSModel::RSModel(TreeTypes treeType, bool randomBasis) :
treeType(treeType),
leafSize(0),
randomBasis(randomBasis),
rSearch(NULL)
{
// Nothing to do.
}
// Copy constructor.
RSModel::RSModel(const RSModel& other) :
treeType(other.treeType),
leafSize(other.leafSize),
randomBasis(other.randomBasis),
q(other.q),
rSearch(other.rSearch->Clone())
{
// Nothing to do.
}
// Move constructor.
RSModel::RSModel(RSModel&& other) :
treeType(other.treeType),
leafSize(other.leafSize),
randomBasis(other.randomBasis),
q(std::move(other.q)),
rSearch(std::move(other.rSearch))
{
// Reset other model.
other.treeType = TreeTypes::KD_TREE;
other.leafSize = 0;
other.randomBasis = false;
}
// Copy operator.
RSModel& RSModel::operator=(const RSModel& other)
{
if (this != &other)
{
delete rSearch;
treeType = other.treeType;
leafSize = other.leafSize;
randomBasis = other.randomBasis;
q = other.q;
rSearch = other.rSearch->Clone();
}
return *this;
}
// Move operator.
RSModel& RSModel::operator=(RSModel&& other)
{
if (this != &other)
{
delete rSearch;
treeType = other.treeType;
leafSize = other.leafSize;
randomBasis = other.randomBasis;
q = std::move(other.q);
rSearch = std::move(other.rSearch);
other.treeType = TreeTypes::KD_TREE;
other.leafSize = 0;
other.randomBasis = false;
}
return *this;
}
// Clean memory, if necessary.
RSModel::~RSModel()
{
delete rSearch;
}
void RSModel::InitializeModel(const bool naive, const bool singleMode)
{
// Clean memory, if necessary.
delete rSearch;
switch (treeType)
{
case KD_TREE:
rSearch = new LeafSizeRSWrapper<tree::KDTree>(naive, singleMode);
break;
case COVER_TREE:
rSearch = new RSWrapper<tree::StandardCoverTree>(naive, singleMode);
break;
case R_TREE:
rSearch = new RSWrapper<tree::RTree>(naive, singleMode);
break;
case R_STAR_TREE:
rSearch = new RSWrapper<tree::RStarTree>(naive, singleMode);
break;
case BALL_TREE:
rSearch = new LeafSizeRSWrapper<tree::BallTree>(naive, singleMode);
break;
case X_TREE:
rSearch = new RSWrapper<tree::XTree>(naive, singleMode);
break;
case HILBERT_R_TREE:
rSearch = new RSWrapper<tree::HilbertRTree>(naive, singleMode);
break;
case R_PLUS_TREE:
rSearch = new RSWrapper<tree::RPlusTree>(naive, singleMode);
break;
case R_PLUS_PLUS_TREE:
rSearch = new RSWrapper<tree::RPlusPlusTree>(naive, singleMode);
break;
case VP_TREE:
rSearch = new LeafSizeRSWrapper<tree::VPTree>(naive, singleMode);
break;
case RP_TREE:
rSearch = new LeafSizeRSWrapper<tree::RPTree>(naive, singleMode);
break;
case MAX_RP_TREE:
rSearch = new LeafSizeRSWrapper<tree::MaxRPTree>(naive, singleMode);
break;
case UB_TREE:
rSearch = new LeafSizeRSWrapper<tree::UBTree>(naive, singleMode);
break;
case OCTREE:
rSearch = new LeafSizeRSWrapper<tree::Octree>(naive, singleMode);
break;
}
}
void RSModel::BuildModel(util::Timers& timers,
arma::mat&& referenceSet,
const size_t leafSize,
const bool naive,
const bool singleMode)
{
// Initialize random basis if necessary.
if (randomBasis)
{
timers.Start("computing_random_basis");
Log::Info << "Creating random basis..." << std::endl;
math::RandomBasis(q, referenceSet.n_rows);
// Do we need to modify the reference set?
if (randomBasis)
referenceSet = q * referenceSet;
timers.Stop("computing_random_basis");
}
this->leafSize = leafSize;
if (!naive)
Log::Info << "Building reference tree..." << std::endl;
InitializeModel(naive, singleMode);
rSearch->Train(timers, std::move(referenceSet), leafSize);
if (!naive)
Log::Info << "Tree built." << std::endl;
}
// Perform range search.
void RSModel::Search(util::Timers& timers,
arma::mat&& querySet,
const math::Range& range,
std::vector<std::vector<size_t>>& neighbors,
std::vector<std::vector<double>>& distances)
{
// We may need to map the query set randomly.
if (randomBasis)
{
timers.Start("applying_random_basis");
querySet = q * querySet;
timers.Stop("applying_random_basis");
}
Log::Info << "Search for points in the range [" << range.Lo() << ", "
<< range.Hi() << "] with ";
if (!Naive() && !SingleMode())
Log::Info << "dual-tree " << TreeName() << " search..." << std::endl;
else if (!Naive())
Log::Info << "single-tree " << TreeName() << " search..." << std::endl;
else
Log::Info << "brute-force (naive) search..." << std::endl;
rSearch->Search(timers, std::move(querySet), range, neighbors, distances,
leafSize);
}
// Perform range search (monochromatic case).
void RSModel::Search(util::Timers& timers,
const math::Range& range,
std::vector<std::vector<size_t>>& neighbors,
std::vector<std::vector<double>>& distances)
{
Log::Info << "Search for points in the range [" << range.Lo() << ", "
<< range.Hi() << "] with ";
if (!Naive() && !SingleMode())
Log::Info << "dual-tree " << TreeName() << " search..." << std::endl;
else if (!Naive())
Log::Info << "single-tree " << TreeName() << " search..." << std::endl;
else
Log::Info << "brute-force (naive) search..." << std::endl;
rSearch->Search(timers, range, neighbors, distances);
}
// Get the name of the tree type.
std::string RSModel::TreeName() const
{
switch (treeType)
{
case KD_TREE:
return "kd-tree";
case COVER_TREE:
return "cover tree";
case R_TREE:
return "R tree";
case R_STAR_TREE:
return "R* tree";
case BALL_TREE:
return "ball tree";
case X_TREE:
return "X tree";
case HILBERT_R_TREE:
return "Hilbert R tree";
case R_PLUS_TREE:
return "R+ tree";
case R_PLUS_PLUS_TREE:
return "R++ tree";
case VP_TREE:
return "vantage point tree";
case RP_TREE:
return "random projection tree (mean split)";
case MAX_RP_TREE:
return "random projection tree (max split)";
case UB_TREE:
return "UB tree";
case OCTREE:
return "octree";
default:
return "unknown tree";
}
}
// Clean memory.
void RSModel::CleanMemory()
{
delete rSearch;
}
} // namespace range
} // namespace mlpack
@@ -20,6 +20,276 @@
namespace mlpack {
namespace range {
/**
* Initialize the RSModel with the given tree type and whether or not a random
* basis should be used.
*/
inline RSModel::RSModel(TreeTypes treeType, bool randomBasis) :
treeType(treeType),
leafSize(0),
randomBasis(randomBasis),
rSearch(NULL)
{
// Nothing to do.
}
// Copy constructor.
inline RSModel::RSModel(const RSModel& other) :
treeType(other.treeType),
leafSize(other.leafSize),
randomBasis(other.randomBasis),
q(other.q),
rSearch(other.rSearch->Clone())
{
// Nothing to do.
}
// Move constructor.
inline RSModel::RSModel(RSModel&& other) :
treeType(other.treeType),
leafSize(other.leafSize),
randomBasis(other.randomBasis),
q(std::move(other.q)),
rSearch(std::move(other.rSearch))
{
// Reset other model.
other.treeType = TreeTypes::KD_TREE;
other.leafSize = 0;
other.randomBasis = false;
}
// Copy operator.
inline RSModel& RSModel::operator=(const RSModel& other)
{
if (this != &other)
{
delete rSearch;
treeType = other.treeType;
leafSize = other.leafSize;
randomBasis = other.randomBasis;
q = other.q;
rSearch = other.rSearch->Clone();
}
return *this;
}
// Move operator.
inline RSModel& RSModel::operator=(RSModel&& other)
{
if (this != &other)
{
delete rSearch;
treeType = other.treeType;
leafSize = other.leafSize;
randomBasis = other.randomBasis;
q = std::move(other.q);
rSearch = std::move(other.rSearch);
other.treeType = TreeTypes::KD_TREE;
other.leafSize = 0;
other.randomBasis = false;
}
return *this;
}
// Clean memory, if necessary.
inline RSModel::~RSModel()
{
delete rSearch;
}
inline void RSModel::InitializeModel(const bool naive,
const bool singleMode)
{
// Clean memory, if necessary.
delete rSearch;
switch (treeType)
{
case KD_TREE:
rSearch = new LeafSizeRSWrapper<tree::KDTree>(naive, singleMode);
break;
case COVER_TREE:
rSearch = new RSWrapper<tree::StandardCoverTree>(naive, singleMode);
break;
case R_TREE:
rSearch = new RSWrapper<tree::RTree>(naive, singleMode);
break;
case R_STAR_TREE:
rSearch = new RSWrapper<tree::RStarTree>(naive, singleMode);
break;
case BALL_TREE:
rSearch = new LeafSizeRSWrapper<tree::BallTree>(naive, singleMode);
break;
case X_TREE:
rSearch = new RSWrapper<tree::XTree>(naive, singleMode);
break;
case HILBERT_R_TREE:
rSearch = new RSWrapper<tree::HilbertRTree>(naive, singleMode);
break;
case R_PLUS_TREE:
rSearch = new RSWrapper<tree::RPlusTree>(naive, singleMode);
break;
case R_PLUS_PLUS_TREE:
rSearch = new RSWrapper<tree::RPlusPlusTree>(naive, singleMode);
break;
case VP_TREE:
rSearch = new LeafSizeRSWrapper<tree::VPTree>(naive, singleMode);
break;
case RP_TREE:
rSearch = new LeafSizeRSWrapper<tree::RPTree>(naive, singleMode);
break;
case MAX_RP_TREE:
rSearch = new LeafSizeRSWrapper<tree::MaxRPTree>(naive, singleMode);
break;
case UB_TREE:
rSearch = new LeafSizeRSWrapper<tree::UBTree>(naive, singleMode);
break;
case OCTREE:
rSearch = new LeafSizeRSWrapper<tree::Octree>(naive, singleMode);
break;
}
}
inline void RSModel::BuildModel(util::Timers& timers,
arma::mat&& referenceSet,
const size_t leafSize,
const bool naive,
const bool singleMode)
{
// Initialize random basis if necessary.
if (randomBasis)
{
timers.Start("computing_random_basis");
Log::Info << "Creating random basis..." << std::endl;
math::RandomBasis(q, referenceSet.n_rows);
// Do we need to modify the reference set?
if (randomBasis)
referenceSet = q * referenceSet;
timers.Stop("computing_random_basis");
}
this->leafSize = leafSize;
if (!naive)
Log::Info << "Building reference tree..." << std::endl;
InitializeModel(naive, singleMode);
rSearch->Train(timers, std::move(referenceSet), leafSize);
if (!naive)
Log::Info << "Tree built." << std::endl;
}
// Perform range search.
inline void RSModel::Search(util::Timers& timers,
arma::mat&& querySet,
const math::Range& range,
std::vector<std::vector<size_t>>& neighbors,
std::vector<std::vector<double>>& distances)
{
// We may need to map the query set randomly.
if (randomBasis)
{
timers.Start("applying_random_basis");
querySet = q * querySet;
timers.Stop("applying_random_basis");
}
Log::Info << "Search for points in the range [" << range.Lo() << ", "
<< range.Hi() << "] with ";
if (!Naive() && !SingleMode())
Log::Info << "dual-tree " << TreeName() << " search..." << std::endl;
else if (!Naive())
Log::Info << "single-tree " << TreeName() << " search..." << std::endl;
else
Log::Info << "brute-force (naive) search..." << std::endl;
rSearch->Search(timers, std::move(querySet), range, neighbors, distances,
leafSize);
}
// Perform range search (monochromatic case).
inline void RSModel::Search(util::Timers& timers,
const math::Range& range,
std::vector<std::vector<size_t>>& neighbors,
std::vector<std::vector<double>>& distances)
{
Log::Info << "Search for points in the range [" << range.Lo() << ", "
<< range.Hi() << "] with ";
if (!Naive() && !SingleMode())
Log::Info << "dual-tree " << TreeName() << " search..." << std::endl;
else if (!Naive())
Log::Info << "single-tree " << TreeName() << " search..." << std::endl;
else
Log::Info << "brute-force (naive) search..." << std::endl;
rSearch->Search(timers, range, neighbors, distances);
}
// Get the name of the tree type.
inline std::string RSModel::TreeName() const
{
switch (treeType)
{
case KD_TREE:
return "kd-tree";
case COVER_TREE:
return "cover tree";
case R_TREE:
return "R tree";
case R_STAR_TREE:
return "R* tree";
case BALL_TREE:
return "ball tree";
case X_TREE:
return "X tree";
case HILBERT_R_TREE:
return "Hilbert R tree";
case R_PLUS_TREE:
return "R+ tree";
case R_PLUS_PLUS_TREE:
return "R++ tree";
case VP_TREE:
return "vantage point tree";
case RP_TREE:
return "random projection tree (mean split)";
case MAX_RP_TREE:
return "random projection tree (max split)";
case UB_TREE:
return "UB tree";
case OCTREE:
return "octree";
default:
return "unknown tree";
}
}
// Clean memory.
inline void RSModel::CleanMemory()
{
delete rSearch;
}
template<template<typename TreeMetricType,
typename TreeStatType,
typename TreeMatType> class TreeType>
+1 -2
View File
@@ -18,12 +18,11 @@ set(SOURCES
# utilities
ra_util.hpp
ra_util.cpp
ra_util_impl.hpp
# model
ra_model.hpp
ra_model_impl.hpp
ra_model.cpp
)
# add directory name to sources
-238
View File
@@ -1,238 +0,0 @@
/**
* @file methods/rann/ra_model.cpp
* @author Ryan Curtin
*
* Implementation of the RAModel class.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "ra_model.hpp"
#include <mlpack/core/math/random_basis.hpp>
namespace mlpack {
namespace neighbor {
RAModel::RAModel(const TreeTypes treeType, const bool randomBasis) :
treeType(treeType),
leafSize(20),
randomBasis(randomBasis),
raSearch(NULL)
{
// Nothing to do.
}
// Copy constructor.
RAModel::RAModel(const RAModel& other) :
treeType(other.treeType),
leafSize(other.leafSize),
randomBasis(other.randomBasis),
q(other.q),
raSearch(other.raSearch->Clone())
{
// Nothing to do.
}
// Move constructor.
RAModel::RAModel(RAModel&& other) :
treeType(other.treeType),
leafSize(other.leafSize),
randomBasis(other.randomBasis),
q(std::move(other.q)),
raSearch(std::move(other.raSearch))
{
// Clear other model.
other.treeType = TreeTypes::KD_TREE;
other.leafSize = 20;
other.randomBasis = false;
}
// Copy operator.
RAModel& RAModel::operator=(const RAModel& other)
{
if (this != &other)
{
// Clear current model.
delete raSearch;
treeType = other.treeType;
leafSize = other.leafSize;
randomBasis = other.randomBasis;
q = other.q;
raSearch = other.raSearch->Clone();
}
return *this;
}
RAModel& RAModel::operator=(RAModel&& other)
{
if (this != &other)
{
// Clear current model.
delete raSearch;
treeType = other.treeType;
leafSize = other.leafSize;
randomBasis = other.randomBasis;
q = std::move(other.q);
raSearch = std::move(other.raSearch);
// Reset other model.
other.treeType = TreeTypes::KD_TREE;
other.leafSize = 20;
other.randomBasis = false;
}
return *this;
}
// Clean memory, if necessary
RAModel::~RAModel()
{
delete raSearch;
}
void RAModel::InitializeModel(const bool naive, const bool singleMode)
{
// Clean memory, if necessary.
delete raSearch;
switch (treeType)
{
case KD_TREE:
raSearch = new LeafSizeRAWrapper<tree::KDTree>(naive, singleMode);
break;
case COVER_TREE:
raSearch = new RAWrapper<tree::StandardCoverTree>(naive, singleMode);
break;
case R_TREE:
raSearch = new RAWrapper<tree::RTree>(naive, singleMode);
break;
case R_STAR_TREE:
raSearch = new RAWrapper<tree::RStarTree>(naive, singleMode);
break;
case X_TREE:
raSearch = new RAWrapper<tree::XTree>(naive, singleMode);
break;
case HILBERT_R_TREE:
raSearch = new RAWrapper<tree::HilbertRTree>(naive, singleMode);
break;
case R_PLUS_TREE:
raSearch = new RAWrapper<tree::RPlusTree>(naive, singleMode);
break;
case R_PLUS_PLUS_TREE:
raSearch = new RAWrapper<tree::RPlusPlusTree>(naive, singleMode);
break;
case UB_TREE:
raSearch = new LeafSizeRAWrapper<tree::UBTree>(naive, singleMode);
break;
case OCTREE:
raSearch = new LeafSizeRAWrapper<tree::Octree>(naive, singleMode);
break;
}
}
void RAModel::BuildModel(util::Timers& timers,
arma::mat&& referenceSet,
const size_t leafSize,
const bool naive,
const bool singleMode)
{
// Initialize random basis, if necessary.
if (randomBasis)
{
timers.Start("computing_random_basis");
Log::Info << "Creating random basis..." << std::endl;
math::RandomBasis(q, referenceSet.n_rows);
referenceSet = q * referenceSet;
timers.Stop("computing_random_basis");
}
this->leafSize = leafSize;
if (!naive)
Log::Info << "Building reference tree..." << std::endl;
InitializeModel(naive, singleMode);
raSearch->Train(timers, std::move(referenceSet), leafSize);
if (!naive)
Log::Info << "Tree built." << std::endl;
}
void RAModel::Search(util::Timers& timers,
arma::mat&& querySet,
const size_t k,
arma::Mat<size_t>& neighbors,
arma::mat& distances)
{
// Apply the random basis if necessary.
if (randomBasis)
querySet = q * querySet;
Log::Info << "Searching for " << k << " approximate nearest neighbors with ";
if (!Naive() && !SingleMode())
Log::Info << "dual-tree rank-approximate " << TreeName() << " search...";
else if (!Naive())
Log::Info << "single-tree rank-approximate " << TreeName() << " search...";
else
Log::Info << "brute-force (naive) rank-approximate search...";
Log::Info << std::endl;
raSearch->Search(timers, std::move(querySet), k, neighbors, distances,
leafSize);
}
void RAModel::Search(util::Timers& timers,
const size_t k,
arma::Mat<size_t>& neighbors,
arma::mat& distances)
{
Log::Info << "Searching for " << k << " approximate nearest neighbors with ";
if (!Naive() && !SingleMode())
Log::Info << "dual-tree rank-approximate " << TreeName() << " search...";
else if (!Naive())
Log::Info << "single-tree rank-approximate " << TreeName() << " search...";
else
Log::Info << "brute-force (naive) rank-approximate search...";
Log::Info << std::endl;
raSearch->Search(timers, k, neighbors, distances);
}
std::string RAModel::TreeName() const
{
switch (treeType)
{
case KD_TREE:
return "kd-tree";
case COVER_TREE:
return "cover tree";
case R_TREE:
return "R tree";
case R_STAR_TREE:
return "R* tree";
case X_TREE:
return "X tree";
case HILBERT_R_TREE:
return "Hilbert R tree";
case R_PLUS_TREE:
return "R+ tree";
case R_PLUS_PLUS_TREE:
return "R++ tree";
case UB_TREE:
return "UB tree";
case OCTREE:
return "octree";
default:
return "unknown tree";
}
}
} // namespace neighbor
} // namespace mlpack
+220
View File
@@ -19,6 +19,226 @@
namespace mlpack {
namespace neighbor {
inline RAModel::RAModel(const TreeTypes treeType, const bool randomBasis) :
treeType(treeType),
leafSize(20),
randomBasis(randomBasis),
raSearch(NULL)
{
// Nothing to do.
}
// Copy constructor.
inline RAModel::RAModel(const RAModel& other) :
treeType(other.treeType),
leafSize(other.leafSize),
randomBasis(other.randomBasis),
q(other.q),
raSearch(other.raSearch->Clone())
{
// Nothing to do.
}
// Move constructor.
inline RAModel::RAModel(RAModel&& other) :
treeType(other.treeType),
leafSize(other.leafSize),
randomBasis(other.randomBasis),
q(std::move(other.q)),
raSearch(std::move(other.raSearch))
{
// Clear other model.
other.treeType = TreeTypes::KD_TREE;
other.leafSize = 20;
other.randomBasis = false;
}
// Copy operator.
inline RAModel& RAModel::operator=(const RAModel& other)
{
if (this != &other)
{
// Clear current model.
delete raSearch;
treeType = other.treeType;
leafSize = other.leafSize;
randomBasis = other.randomBasis;
q = other.q;
raSearch = other.raSearch->Clone();
}
return *this;
}
inline RAModel& RAModel::operator=(RAModel&& other)
{
if (this != &other)
{
// Clear current model.
delete raSearch;
treeType = other.treeType;
leafSize = other.leafSize;
randomBasis = other.randomBasis;
q = std::move(other.q);
raSearch = std::move(other.raSearch);
// Reset other model.
other.treeType = TreeTypes::KD_TREE;
other.leafSize = 20;
other.randomBasis = false;
}
return *this;
}
// Clean memory, if necessary
inline RAModel::~RAModel()
{
delete raSearch;
}
inline void RAModel::InitializeModel(const bool naive,
const bool singleMode)
{
// Clean memory, if necessary.
delete raSearch;
switch (treeType)
{
case KD_TREE:
raSearch = new LeafSizeRAWrapper<tree::KDTree>(naive, singleMode);
break;
case COVER_TREE:
raSearch = new RAWrapper<tree::StandardCoverTree>(naive, singleMode);
break;
case R_TREE:
raSearch = new RAWrapper<tree::RTree>(naive, singleMode);
break;
case R_STAR_TREE:
raSearch = new RAWrapper<tree::RStarTree>(naive, singleMode);
break;
case X_TREE:
raSearch = new RAWrapper<tree::XTree>(naive, singleMode);
break;
case HILBERT_R_TREE:
raSearch = new RAWrapper<tree::HilbertRTree>(naive, singleMode);
break;
case R_PLUS_TREE:
raSearch = new RAWrapper<tree::RPlusTree>(naive, singleMode);
break;
case R_PLUS_PLUS_TREE:
raSearch = new RAWrapper<tree::RPlusPlusTree>(naive, singleMode);
break;
case UB_TREE:
raSearch = new LeafSizeRAWrapper<tree::UBTree>(naive, singleMode);
break;
case OCTREE:
raSearch = new LeafSizeRAWrapper<tree::Octree>(naive, singleMode);
break;
}
}
inline void RAModel::BuildModel(util::Timers& timers,
arma::mat&& referenceSet,
const size_t leafSize,
const bool naive,
const bool singleMode)
{
// Initialize random basis, if necessary.
if (randomBasis)
{
timers.Start("computing_random_basis");
Log::Info << "Creating random basis..." << std::endl;
math::RandomBasis(q, referenceSet.n_rows);
referenceSet = q * referenceSet;
timers.Stop("computing_random_basis");
}
this->leafSize = leafSize;
if (!naive)
Log::Info << "Building reference tree..." << std::endl;
InitializeModel(naive, singleMode);
raSearch->Train(timers, std::move(referenceSet), leafSize);
if (!naive)
Log::Info << "Tree built." << std::endl;
}
inline void RAModel::Search(util::Timers& timers,
arma::mat&& querySet,
const size_t k,
arma::Mat<size_t>& neighbors,
arma::mat& distances)
{
// Apply the random basis if necessary.
if (randomBasis)
querySet = q * querySet;
Log::Info << "Searching for " << k << " approximate nearest neighbors with ";
if (!Naive() && !SingleMode())
Log::Info << "dual-tree rank-approximate " << TreeName() << " search...";
else if (!Naive())
Log::Info << "single-tree rank-approximate " << TreeName() << " search...";
else
Log::Info << "brute-force (naive) rank-approximate search...";
Log::Info << std::endl;
raSearch->Search(timers, std::move(querySet), k, neighbors, distances,
leafSize);
}
inline void RAModel::Search(util::Timers& timers,
const size_t k,
arma::Mat<size_t>& neighbors,
arma::mat& distances)
{
Log::Info << "Searching for " << k << " approximate nearest neighbors with ";
if (!Naive() && !SingleMode())
Log::Info << "dual-tree rank-approximate " << TreeName() << " search...";
else if (!Naive())
Log::Info << "single-tree rank-approximate " << TreeName() << " search...";
else
Log::Info << "brute-force (naive) rank-approximate search...";
Log::Info << std::endl;
raSearch->Search(timers, k, neighbors, distances);
}
inline std::string RAModel::TreeName() const
{
switch (treeType)
{
case KD_TREE:
return "kd-tree";
case COVER_TREE:
return "cover tree";
case R_TREE:
return "R tree";
case R_STAR_TREE:
return "R* tree";
case X_TREE:
return "X tree";
case HILBERT_R_TREE:
return "Hilbert R tree";
case R_PLUS_TREE:
return "R+ tree";
case R_PLUS_PLUS_TREE:
return "R++ tree";
case UB_TREE:
return "UB tree";
case OCTREE:
return "octree";
default:
return "unknown tree";
}
}
template<template<typename TreeMetricType,
typename TreeStatType,
typename TreeMatType> class TreeType>
+3 -13
View File
@@ -48,22 +48,12 @@ class RAUtil
const size_t k,
const size_t m,
const size_t t);
/**
* Pick up desired number of samples (with replacement) from a given range
* of integers so that only the distinct samples are returned from
* the range [0 - specified upper bound)
*
* @param numSamples Number of random samples.
* @param rangeUpperBound The upper bound on the range of integers.
* @param distinctSamples The list of the distinct samples.
*/
static void ObtainDistinctSamples(const size_t numSamples,
const size_t rangeUpperBound,
arma::uvec& distinctSamples);
};
} // namespace neighbor
} // namespace mlpack
// Include implementation.
#include "ra_util_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/rann/ra_util.cpp
* @file methods/rann/ra_util_impl.hpp
* @author Parikshit Ram
* @author Ryan Curtin
*
@@ -10,15 +10,18 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_RANN_RA_UTIL_IMPL_HPP
#define MLPACK_METHODS_RANN_RA_UTIL_IMPL_HPP
#include "ra_util.hpp"
using namespace mlpack;
using namespace mlpack::neighbor;
namespace mlpack {
namespace neighbor {
size_t mlpack::neighbor::RAUtil::MinimumSamplesReqd(const size_t n,
const size_t k,
const double tau,
const double alpha)
inline size_t RAUtil::MinimumSamplesReqd(const size_t n,
const size_t k,
const double tau,
const double alpha)
{
size_t ub = n; // The upper bound on the binary search.
size_t lb = k; // The lower bound on the binary search.
@@ -69,10 +72,10 @@ size_t mlpack::neighbor::RAUtil::MinimumSamplesReqd(const size_t n,
return (std::min(m + 1, n));
}
double mlpack::neighbor::RAUtil::SuccessProbability(const size_t n,
const size_t k,
const size_t m,
const size_t t)
inline double RAUtil::SuccessProbability(const size_t n,
const size_t k,
const size_t m,
const size_t t)
{
if (k == 1)
{
@@ -161,3 +164,8 @@ double mlpack::neighbor::RAUtil::SuccessProbability(const size_t n,
return sum;
} // For k > 1.
}
} // namespace neighbor
} // namespace mlpack
#endif
@@ -272,7 +272,7 @@ inline double ParallelSGD<ExponentialBackoff>::Optimize(
if (shuffle) // Determine order of visitation.
std::shuffle(visitationOrder.begin(), visitationOrder.end(),
mlpack::math::randGen);
mlpack::math::RandGen());
#pragma omp parallel
{
@@ -2,7 +2,6 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
env_type.hpp
env_type.cpp
mountain_car.hpp
cart_pole.hpp
continuous_mountain_car.hpp
@@ -1,29 +0,0 @@
/**
* @file methods/reinforcement_learning/environment/env_type.cpp
* @author Nishant Kumar
*
* This file defines the static variables used by the discrete and continuous
* environments.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "env_type.hpp"
namespace mlpack {
namespace rl {
// Instantiate static members.
size_t DiscreteActionEnv::State::dimension = 0;
size_t DiscreteActionEnv::Action::size = 0;
size_t DiscreteActionEnv::rewardSize = 0;
size_t ContinuousActionEnv::State::dimension = 0;
size_t ContinuousActionEnv::Action::size = 0;
size_t ContinuousActionEnv::rewardSize = 0;
} // namespace rl
} // namespace mlpack
@@ -206,6 +206,19 @@ class ContinuousActionEnv
static size_t rewardSize;
};
#ifndef MLPACK_METHODS_RL_ENVIRONMENT_ENV_TYPE_VARIABLES
#define MLPACK_METHODS_RL_ENVIRONMENT_ENV_TYPE_VARIABLES
size_t DiscreteActionEnv::State::dimension = 0;
size_t DiscreteActionEnv::Action::size = 0;
size_t DiscreteActionEnv::rewardSize = 0;
size_t ContinuousActionEnv::State::dimension = 0;
size_t ContinuousActionEnv::Action::size = 0;
size_t ContinuousActionEnv::rewardSize = 0;
#endif
} // namespace rl
} // namespace mlpack
@@ -2,10 +2,9 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
softmax_regression.hpp
softmax_regression.cpp
softmax_regression_impl.hpp
softmax_regression_function.hpp
softmax_regression_function.cpp
softmax_regression_function_impl.hpp
)
# Add directory name to sources.
@@ -1,140 +0,0 @@
/**
* @file methods/softmax_regression/softmax_regression.cpp
* @author Siddharth Agrawal
*
* Implementation of softmax regression.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "softmax_regression.hpp"
namespace mlpack {
namespace regression {
SoftmaxRegression::
SoftmaxRegression(const size_t inputSize,
const size_t numClasses,
const bool fitIntercept) :
numClasses(numClasses),
lambda(0.0001),
fitIntercept(fitIntercept)
{
SoftmaxRegressionFunction::InitializeWeights(
parameters, inputSize, numClasses, fitIntercept);
}
void SoftmaxRegression::Classify(const arma::mat& dataset,
arma::Row<size_t>& labels)
const
{
arma::mat probabilities;
Classify(dataset, probabilities);
// Prepare necessary data.
labels.zeros(dataset.n_cols);
double maxProbability = 0;
// For each test input.
for (size_t i = 0; i < dataset.n_cols; ++i)
{
// For each class.
for (size_t j = 0; j < numClasses; ++j)
{
// If a higher class probability is encountered, change prediction.
if (probabilities(j, i) > maxProbability)
{
maxProbability = probabilities(j, i);
labels(i) = j;
}
}
// Set maximum probability to zero for the next input.
maxProbability = 0;
}
}
void SoftmaxRegression::Classify(const arma::mat& dataset,
arma::Row<size_t>& labels,
arma::mat& probabilities)
const
{
Classify(dataset, probabilities);
// Prepare necessary data.
labels.zeros(dataset.n_cols);
double maxProbability = 0;
// For each test input.
for (size_t i = 0; i < dataset.n_cols; ++i)
{
// For each class.
for (size_t j = 0; j < numClasses; ++j)
{
// If a higher class probability is encountered, change prediction.
if (probabilities(j, i) > maxProbability)
{
maxProbability = probabilities(j, i);
labels(i) = j;
}
}
// Set maximum probability to zero for the next input.
maxProbability = 0;
}
}
void SoftmaxRegression::Classify(const arma::mat& dataset,
arma::mat& probabilities)
const
{
util::CheckSameDimensionality(dataset, FeatureSize(),
"SoftmaxRegression::Classify()");
// Calculate the probabilities for each test input.
arma::mat hypothesis;
if (fitIntercept)
{
// In order to add the intercept term, we should compute following matrix:
// [1; data] = arma::join_cols(ones(1, data.n_cols), data)
// hypothesis = arma::exp(parameters * [1; data]).
//
// Since the cost of join maybe high due to the copy of original data,
// split the hypothesis computation to two components.
hypothesis = arma::exp(
arma::repmat(parameters.col(0), 1, dataset.n_cols) +
parameters.cols(1, parameters.n_cols - 1) * dataset);
}
else
{
hypothesis = arma::exp(parameters * dataset);
}
probabilities = hypothesis / arma::repmat(arma::sum(hypothesis, 0),
numClasses, 1);
}
double SoftmaxRegression::ComputeAccuracy(
const arma::mat& testData,
const arma::Row<size_t>& labels) const
{
arma::Row<size_t> predictions;
// Get predictions for the provided data.
Classify(testData, predictions);
// Increment count for every correctly predicted label.
size_t count = 0;
for (size_t i = 0; i < predictions.n_elem; ++i)
if (predictions(i) == labels(i))
count++;
// Return percentage accuracy.
return (count * 100.0) / predictions.n_elem;
}
} // namespace regression
} // namespace mlpack
@@ -14,6 +14,7 @@
#define MLPACK_METHODS_SOFTMAX_REGRESSION_SOFTMAX_REGRESSION_FUNCTION_HPP
#include <mlpack/prereqs.hpp>
#include <mlpack/core/math/make_alias.hpp>
namespace mlpack {
namespace regression {
@@ -205,4 +206,7 @@ class SoftmaxRegressionFunction
} // namespace regression
} // namespace mlpack
// Include implementation.
#include "softmax_regression_function_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/softmax_regression/softmax_regression_function.cpp
* @file methods/softmax_regression/softmax_regression_function_impl.hpp
* @author Siddharth Agrawal
*
* Implementation of function to be optimized for softmax regression.
@@ -9,13 +9,15 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_SOFTMAX_REGRESSION_SOFTMAX_REGRESSION_FUNCTION_IMPL_HPP
#define MLPACK_METHODS_SOFTMAX_REGRESSION_SOFTMAX_REGRESSION_FUNCTION_IMPL_HPP
#include "softmax_regression_function.hpp"
#include <mlpack/core/math/make_alias.hpp>
using namespace mlpack;
using namespace mlpack::regression;
namespace mlpack {
namespace regression {
SoftmaxRegressionFunction::SoftmaxRegressionFunction(
inline SoftmaxRegressionFunction::SoftmaxRegressionFunction(
const arma::mat& data,
const arma::Row<size_t>& labels,
const size_t numClasses,
@@ -36,7 +38,7 @@ SoftmaxRegressionFunction::SoftmaxRegressionFunction(
/**
* Shuffle the data.
*/
void SoftmaxRegressionFunction::Shuffle()
inline void SoftmaxRegressionFunction::Shuffle()
{
// Determine new ordering.
arma::uvec ordering = arma::shuffle(arma::linspace<arma::uvec>(0,
@@ -75,12 +77,12 @@ void SoftmaxRegressionFunction::Shuffle()
* normal distribution. The weights cannot be initialized to zero, as that will
* lead to each class output being the same.
*/
const arma::mat SoftmaxRegressionFunction::InitializeWeights()
inline const arma::mat SoftmaxRegressionFunction::InitializeWeights()
{
return InitializeWeights(data.n_rows, numClasses, fitIntercept);
}
const arma::mat SoftmaxRegressionFunction::InitializeWeights(
inline const arma::mat SoftmaxRegressionFunction::InitializeWeights(
const size_t featureSize,
const size_t numClasses,
const bool fitIntercept)
@@ -90,7 +92,7 @@ const arma::mat SoftmaxRegressionFunction::InitializeWeights(
return parameters;
}
void SoftmaxRegressionFunction::InitializeWeights(
inline void SoftmaxRegressionFunction::InitializeWeights(
arma::mat &weights,
const size_t featureSize,
const size_t numClasses,
@@ -111,7 +113,7 @@ void SoftmaxRegressionFunction::InitializeWeights(
* labels. The output is in the form of a matrix, which leads to simpler
* calculations in the Evaluate() and Gradient() methods.
*/
void SoftmaxRegressionFunction::GetGroundTruthMatrix(
inline void SoftmaxRegressionFunction::GetGroundTruthMatrix(
const arma::Row<size_t>& labels, arma::sp_mat& groundTruth)
{
// Calculate the ground truth matrix according to the labels passed. The
@@ -147,7 +149,7 @@ void SoftmaxRegressionFunction::GetGroundTruthMatrix(
* Evaluate the probabilities matrix. If fitIntercept flag is true,
* it should consider the parameters.cols(0) intercept term.
*/
void SoftmaxRegressionFunction::GetProbabilitiesMatrix(
inline void SoftmaxRegressionFunction::GetProbabilitiesMatrix(
const arma::mat& parameters,
arma::mat& probabilities,
const size_t start,
@@ -181,7 +183,8 @@ void SoftmaxRegressionFunction::GetProbabilitiesMatrix(
/**
* Evaluates the objective function given the parameters.
*/
double SoftmaxRegressionFunction::Evaluate(const arma::mat& parameters) const
inline double SoftmaxRegressionFunction::Evaluate(const arma::mat& parameters)
const
{
// The objective function is the negative log likelihood of the model
// calculated over all the training examples. Mathematically it is as follows:
@@ -219,9 +222,9 @@ double SoftmaxRegressionFunction::Evaluate(const arma::mat& parameters) const
/**
* Evaluate the objective function for the given points given the parameters.
*/
double SoftmaxRegressionFunction::Evaluate(const arma::mat& parameters,
const size_t start,
const size_t batchSize) const
inline double SoftmaxRegressionFunction::Evaluate(const arma::mat& parameters,
const size_t start,
const size_t batchSize) const
{
arma::mat probabilities;
GetProbabilitiesMatrix(parameters, probabilities, start, batchSize);
@@ -239,8 +242,8 @@ double SoftmaxRegressionFunction::Evaluate(const arma::mat& parameters,
/**
* Calculates and stores the gradient values given a set of parameters.
*/
void SoftmaxRegressionFunction::Gradient(const arma::mat& parameters,
arma::mat& gradient) const
inline void SoftmaxRegressionFunction::Gradient(const arma::mat& parameters,
arma::mat& gradient) const
{
// Calculate the class probabilities for each training example. The
// probabilities for each of the classes are given by:
@@ -272,10 +275,11 @@ void SoftmaxRegressionFunction::Gradient(const arma::mat& parameters,
}
}
void SoftmaxRegressionFunction::Gradient(const arma::mat& parameters,
const size_t start,
arma::mat& gradient,
const size_t batchSize) const
inline void SoftmaxRegressionFunction::Gradient(const arma::mat& parameters,
const size_t start,
arma::mat& gradient,
const size_t batchSize)
const
{
arma::mat probabilities;
GetProbabilitiesMatrix(parameters, probabilities, start, batchSize);
@@ -301,9 +305,10 @@ void SoftmaxRegressionFunction::Gradient(const arma::mat& parameters,
}
}
void SoftmaxRegressionFunction::PartialGradient(const arma::mat& parameters,
const size_t j,
arma::sp_mat& gradient) const
inline void SoftmaxRegressionFunction::PartialGradient(const arma::mat& parameters,
const size_t j,
arma::sp_mat& gradient)
const
{
gradient.zeros(arma::size(parameters));
@@ -332,3 +337,8 @@ void SoftmaxRegressionFunction::PartialGradient(const arma::mat& parameters,
parameters.col(j);
}
}
} // namespace regression
} // namespace mlpack
#endif
@@ -49,6 +49,108 @@ SoftmaxRegression::SoftmaxRegression(
Train(data, labels, numClasses, optimizer, callbacks...);
}
inline SoftmaxRegression::SoftmaxRegression(
const size_t inputSize,
const size_t numClasses,
const bool fitIntercept) :
numClasses(numClasses),
lambda(0.0001),
fitIntercept(fitIntercept)
{
SoftmaxRegressionFunction::InitializeWeights(
parameters, inputSize, numClasses, fitIntercept);
}
inline void SoftmaxRegression::Classify(const arma::mat& dataset,
arma::Row<size_t>& labels)
const
{
arma::mat probabilities;
Classify(dataset, probabilities);
// Prepare necessary data.
labels.zeros(dataset.n_cols);
double maxProbability = 0;
// For each test input.
for (size_t i = 0; i < dataset.n_cols; ++i)
{
// For each class.
for (size_t j = 0; j < numClasses; ++j)
{
// If a higher class probability is encountered, change prediction.
if (probabilities(j, i) > maxProbability)
{
maxProbability = probabilities(j, i);
labels(i) = j;
}
}
// Set maximum probability to zero for the next input.
maxProbability = 0;
}
}
inline void SoftmaxRegression::Classify(const arma::mat& dataset,
arma::Row<size_t>& labels,
arma::mat& probabilities)
const
{
Classify(dataset, probabilities);
// Prepare necessary data.
labels.zeros(dataset.n_cols);
double maxProbability = 0;
// For each test input.
for (size_t i = 0; i < dataset.n_cols; ++i)
{
// For each class.
for (size_t j = 0; j < numClasses; ++j)
{
// If a higher class probability is encountered, change prediction.
if (probabilities(j, i) > maxProbability)
{
maxProbability = probabilities(j, i);
labels(i) = j;
}
}
// Set maximum probability to zero for the next input.
maxProbability = 0;
}
}
inline void SoftmaxRegression::Classify(const arma::mat& dataset,
arma::mat& probabilities)
const
{
util::CheckSameDimensionality(dataset, FeatureSize(),
"SoftmaxRegression::Classify()");
// Calculate the probabilities for each test input.
arma::mat hypothesis;
if (fitIntercept)
{
// In order to add the intercept term, we should compute following matrix:
// [1; data] = arma::join_cols(ones(1, data.n_cols), data)
// hypothesis = arma::exp(parameters * [1; data]).
//
// Since the cost of join maybe high due to the copy of original data,
// split the hypothesis computation to two components.
hypothesis = arma::exp(
arma::repmat(parameters.col(0), 1, dataset.n_cols) +
parameters.cols(1, parameters.n_cols - 1) * dataset);
}
else
{
hypothesis = arma::exp(parameters * dataset);
}
probabilities = hypothesis / arma::repmat(arma::sum(hypothesis, 0),
numClasses, 1);
}
template<typename VecType>
size_t SoftmaxRegression::Classify(const VecType& point) const
{
@@ -57,6 +159,25 @@ size_t SoftmaxRegression::Classify(const VecType& point) const
return size_t(label(0));
}
inline double SoftmaxRegression::ComputeAccuracy(
const arma::mat& testData,
const arma::Row<size_t>& labels) const
{
arma::Row<size_t> predictions;
// Get predictions for the provided data.
Classify(testData, predictions);
// Increment count for every correctly predicted label.
size_t count = 0;
for (size_t i = 0; i < predictions.n_elem; ++i)
if (predictions(i) == labels(i))
count++;
// Return percentage accuracy.
return (count * 100.0) / predictions.n_elem;
}
template<typename OptimizerType>
double SoftmaxRegression::Train(const arma::mat& data,
const arma::Row<size_t>& labels,
@@ -2,12 +2,11 @@
# Anything not in this list will not be compiled into mlpack.
set(SOURCES
sparse_autoencoder.hpp
sparse_autoencoder.cpp
sparse_autoencoder_impl.hpp
sparse_autoencoder_function.hpp
sparse_autoencoder_function.cpp
sparse_autoencoder_function_impl.hpp
maximal_inputs.hpp
maximal_inputs.cpp
maximal_inputs_impl.hpp
)
# Add directory name to sources.
@@ -93,4 +93,7 @@ void NormalizeColByMax(const arma::mat& input, arma::mat& output);
} // namespace nn
} // namespace mlpack
// Include implementation.
#include "maximal_inputs_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/sparse_autoencoder/maximal_inputs.cpp
* @file methods/sparse_autoencoder/maximal_inputs_impl.hpp
* @author Tham Ngap Wei
*
* Implementation of MaximalInputs().
@@ -9,12 +9,15 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_NN_MAXIMAL_INPUTS_IMPL_HPP
#define MLPACK_METHODS_NN_MAXIMAL_INPUTS_IMPL_HPP
#include "maximal_inputs.hpp"
namespace mlpack {
namespace nn {
void MaximalInputs(const arma::mat& parameters, arma::mat& output)
inline void MaximalInputs(const arma::mat& parameters, arma::mat& output)
{
arma::mat paramTemp(parameters.submat(0, 0, (parameters.n_rows - 1) / 2 - 1,
parameters.n_cols - 2).t());
@@ -24,8 +27,8 @@ void MaximalInputs(const arma::mat& parameters, arma::mat& output)
NormalizeColByMax(paramTemp, output);
}
void NormalizeColByMax(const arma::mat &input,
arma::mat &output)
inline void NormalizeColByMax(const arma::mat &input,
arma::mat &output)
{
output.set_size(input.n_rows, input.n_cols);
for (arma::uword i = 0; i != input.n_cols; ++i)
@@ -44,3 +47,5 @@ void NormalizeColByMax(const arma::mat &input,
} // namespace nn
} // namespace mlpack
#endif
@@ -1,30 +0,0 @@
/**
* @file methods/sparse_autoencoder/sparse_autoencoder.cpp
* @author Siddharth Agrawal
*
* Implementation of sparse autoencoders.
*
* mlpack is free software; you may redistribute it and/or modify it under the
* terms of the 3-clause BSD license. You should have received a copy of the
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "sparse_autoencoder.hpp"
namespace mlpack {
namespace nn {
void SparseAutoencoder::GetNewFeatures(arma::mat& data,
arma::mat& features)
{
const size_t l1 = hiddenSize;
const size_t l2 = visibleSize;
Sigmoid(parameters.submat(0, 0, l1 - 1, l2 - 1) * data +
arma::repmat(parameters.submat(0, l2, l1 - 1, l2), 1, data.n_cols),
features);
}
} // namespace nn
} // namespace mlpack
@@ -165,4 +165,7 @@ class SparseAutoencoderFunction
} // namespace nn
} // namespace mlpack
// Include implementation.
#include "sparse_autoencoder_function_impl.hpp"
#endif
@@ -1,5 +1,5 @@
/**
* @file methods/sparse_autoencoder/sparse_autoencoder_function.cpp
* @file methods/sparse_autoencoder/sparse_autoencoder_function_impl.hpp
* @author Siddharth Agrawal
*
* Implementation of function to be optimized for sparse autoencoders.
@@ -9,18 +9,21 @@
* 3-clause BSD license along with mlpack. If not, see
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#ifndef MLPACK_METHODS_SPARSE_AUTOENCODER_SPARSE_AUTOENCODER_FUNCTION_IMPL_HPP
#define MLPACK_METHODS_SPARSE_AUTOENCODER_SPARSE_AUTOENCODER_FUNCTION_IMPL_HPP
#include "sparse_autoencoder_function.hpp"
using namespace mlpack;
using namespace mlpack::nn;
using namespace std;
namespace mlpack {
namespace nn {
SparseAutoencoderFunction::SparseAutoencoderFunction(const arma::mat& data,
const size_t visibleSize,
const size_t hiddenSize,
const double lambda,
const double beta,
const double rho) :
inline SparseAutoencoderFunction::SparseAutoencoderFunction(
const arma::mat& data,
const size_t visibleSize,
const size_t hiddenSize,
const double lambda,
const double beta,
const double rho) :
data(data),
visibleSize(visibleSize),
hiddenSize(hiddenSize),
@@ -37,7 +40,7 @@ SparseAutoencoderFunction::SparseAutoencoderFunction(const arma::mat& data,
* [-r, r] where 'r' is decided using the sizes of the visible and hidden
* layers. The biases b1, b2 are initialized to 0.
*/
const arma::mat SparseAutoencoderFunction::InitializeWeights()
inline const arma::mat SparseAutoencoderFunction::InitializeWeights()
{
// The module uses a matrix to store the parameters, its structure looks like:
// vSize 1
@@ -73,7 +76,8 @@ const arma::mat SparseAutoencoderFunction::InitializeWeights()
/** Evaluates the objective function given the parameters.
*/
double SparseAutoencoderFunction::Evaluate(const arma::mat& parameters) const
inline double SparseAutoencoderFunction::Evaluate(const arma::mat& parameters)
const
{
// The objective function is the average squared reconstruction error of the
// network. w1 and b1 are the weights and biases associated with the hidden
@@ -141,8 +145,8 @@ double SparseAutoencoderFunction::Evaluate(const arma::mat& parameters) const
/** Calculates and stores the gradient values given a set of parameters.
*/
void SparseAutoencoderFunction::Gradient(const arma::mat& parameters,
arma::mat& gradient) const
inline void SparseAutoencoderFunction::Gradient(const arma::mat& parameters,
arma::mat& gradient) const
{
// Performs a feedforward pass of the neural network, and computes the
// activations of the output layer as in the Evaluate() method. It uses the
@@ -208,3 +212,8 @@ void SparseAutoencoderFunction::Gradient(const arma::mat& parameters,
gradient.submat(0, l2, l1 - 1, l2) = arma::sum(delHid, 1) / data.n_cols;
gradient.submat(l3, 0, l3, l2 - 1) = (arma::sum(delOut, 1) / data.n_cols).t();
}
} // namespace nn
} // namespace mlpack
#endif

Some files were not shown because too many files have changed in this diff Show More