[mlpack-svn] r10305 - mlpack/trunk/src/mlpack/methods/gmm
fastlab-svn at coffeetalk-1.cc.gatech.edu
fastlab-svn at coffeetalk-1.cc.gatech.edu
Wed Nov 16 17:53:40 EST 2011
Author: rcurtin
Date: 2011-11-16 17:53:40 -0500 (Wed, 16 Nov 2011)
New Revision: 10305
Added:
mlpack/trunk/src/mlpack/methods/gmm/gmm.cpp
mlpack/trunk/src/mlpack/methods/gmm/gmm.hpp
mlpack/trunk/src/mlpack/methods/gmm/gmm_main.cpp
Removed:
mlpack/trunk/src/mlpack/methods/gmm/mog_em.cpp
mlpack/trunk/src/mlpack/methods/gmm/mog_em.hpp
mlpack/trunk/src/mlpack/methods/gmm/mog_em_main.cpp
mlpack/trunk/src/mlpack/methods/gmm/mog_l2e.cpp
mlpack/trunk/src/mlpack/methods/gmm/mog_l2e.hpp
mlpack/trunk/src/mlpack/methods/gmm/mog_l2e_main.cpp
mlpack/trunk/src/mlpack/methods/gmm/optimizers.cpp
mlpack/trunk/src/mlpack/methods/gmm/optimizers.hpp
Modified:
mlpack/trunk/src/mlpack/methods/gmm/CMakeLists.txt
Log:
Remove MoGL2E. That can be implemented some other day.
Modified: mlpack/trunk/src/mlpack/methods/gmm/CMakeLists.txt
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/CMakeLists.txt 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/CMakeLists.txt 2011-11-16 22:53:40 UTC (rev 10305)
@@ -5,16 +5,9 @@
set(SOURCES
kmeans.hpp
kmeans.cpp
- # for mog_em
- mog_em.hpp
- mog_em.cpp
+ gmm.hpp
+ gmm.cpp
phi.hpp
- # for mog_l2e
- mog_l2e.hpp
- mog_l2e.cpp
- # optimizers
- optimizers.hpp
- optimizers.cpp
)
# Add directory name to sources.
@@ -27,27 +20,10 @@
set(MLPACK_SRCS ${MLPACK_SRCS} ${DIR_SRCS} PARENT_SCOPE)
# main executable, em
-add_executable(mog_em
- mog_em_main.cpp
+add_executable(gmm
+ gmm_main.cpp
)
# link dependencies of mog_em
-target_link_libraries(mog_em
+target_link_libraries(gmm
mlpack
)
-
-# main executable, l2e
-add_executable(mog_l2e
- mog_l2e_main.cpp
-)
-target_link_libraries(mog_l2e
- mlpack
-)
-
-# Test executable.
-add_executable(gmm_test
- gmm_test.cpp
-)
-target_link_libraries(gmm_test
- mlpack
- boost_unit_test_framework
-)
Copied: mlpack/trunk/src/mlpack/methods/gmm/gmm.cpp (from rev 10304, mlpack/trunk/src/mlpack/methods/gmm/mog_em.cpp)
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/gmm.cpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/gmm/gmm.cpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -0,0 +1,160 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_em.cpp
+ *
+ * Implementation for the loglikelihood function, the EM algorithm
+ * and also computes the K-means for getting an initial point
+ *
+ */
+#include "mog_em.hpp"
+#include "phi.hpp"
+#include "kmeans.hpp"
+
+using namespace mlpack;
+using namespace gmm;
+
+void MoGEM::ExpectationMaximization(const arma::mat& data)
+{
+ // Create temporary models and set to the right size.
+ std::vector<arma::vec> means_trial(gaussians, arma::vec(dimension));
+ std::vector<arma::mat> covariances_trial(gaussians,
+ arma::mat(dimension, dimension));
+ arma::vec weights_trial(gaussians);
+
+ arma::mat cond_prob(gaussians, data.n_cols);
+
+ long double l, l_old, best_l, TINY = 1.0e-10;
+
+ best_l = -DBL_MAX;
+
+ // We will perform five trials, and then save the trial with the best result
+ // as our trained model.
+ for (size_t iteration = 0; iteration < 5; iteration++)
+ {
+ // Use k-means to find initial values for the parameters.
+ KMeans(data, gaussians, means_trial, covariances_trial, weights_trial);
+
+ Log::Warn << "K-Means results:" << std::endl;
+ for (size_t i = 0; i < gaussians; i++)
+ {
+ Log::Warn << "Mean " << i << ":" << std::endl;
+ Log::Warn << means_trial[i] << std::endl;
+ Log::Warn << "Covariance " << i << ":" << std::endl;
+ Log::Warn << covariances_trial[i] << std::endl;
+ }
+ Log::Warn << "Weights: " << std::endl << weights_trial << std::endl;
+
+ // Calculate the log likelihood of the model.
+ l = Loglikelihood(data, means_trial, covariances_trial, weights_trial);
+
+ l_old = -DBL_MAX;
+
+ // Iterate to update the model until no more improvement is found.
+ size_t max_iterations = 1000;
+ size_t iteration = 0;
+ while (std::abs(l - l_old) > TINY && iteration < max_iterations)
+ {
+ Log::Warn << "Iteration " << iteration << std::endl;
+ // Calculate the conditional probabilities of choosing a particular
+ // Gaussian given the data and the present theta value.
+ for (size_t j = 0; j < data.n_cols; j++)
+ {
+ for (size_t i = 0; i < gaussians; i++)
+ {
+ cond_prob(i, j) = phi(data.unsafe_col(j), means_trial[i],
+ covariances_trial[i]) * weights_trial[i];
+ }
+
+ // Normalize column to have sum probability of one.
+ cond_prob.col(j) /= arma::sum(cond_prob.col(j));
+ }
+
+ // Store the sums of each row because they are used multiple times.
+ arma::vec prob_row_sums = arma::sum(cond_prob, 1 /* row-wise */);
+
+ // Calculate the new value of the means using the updated conditional
+ // probabilities.
+ for (size_t i = 0; i < gaussians; i++)
+ {
+ means_trial[i].zeros();
+ for (size_t j = 0; j < data.n_cols; j++)
+ means_trial[i] += cond_prob(i, j) * data.col(j);
+
+ means_trial[i] /= prob_row_sums[i];
+ }
+
+ // Calculate the new value of the covariances using the updated
+ // conditional probabilities and the updated means.
+ for (size_t i = 0; i < gaussians; i++)
+ {
+ covariances_trial[i].zeros();
+ for (size_t j = 0; j < data.n_cols; j++)
+ {
+ arma::vec tmp = data.col(j) - means_trial[i];
+ covariances_trial[i] += cond_prob(i, j) * (tmp * trans(tmp));
+ }
+
+ covariances_trial[i] /= prob_row_sums[i];
+ }
+
+ // Calculate the new values for omega using the updated conditional
+ // probabilities.
+ weights_trial = prob_row_sums / data.n_cols;
+/*
+ Log::Warn << "Estimated weights:" << std::endl << weights_trial
+ << std::endl;
+
+ for (size_t i = 0; i < gaussians; i++)
+ {
+// Log::Warn << "Estimated mean " << i << ":" << std::endl;
+// Log::Warn << means_trial[i] << std::endl;
+ Log::Warn << "Estimated covariance " << i << ":" << std::endl;
+ Log::Warn << covariances_trial[i] << std::endl;
+ }
+*/
+
+ // Update values of l; calculate new log-likelihood.
+ l_old = l;
+ l = Loglikelihood(data, means_trial, covariances_trial, weights_trial);
+
+ Log::Warn << "Improved log likelihood to " << l << std::endl;
+
+ iteration++;
+ }
+
+ // The trial model is trained. Is it better than our existing model?
+ if (l > best_l)
+ {
+ best_l = l;
+
+ means = means_trial;
+ covariances = covariances_trial;
+ weights = weights_trial;
+ }
+ }
+
+ Log::Info << "Log likelihood value of the estimated model: " << best_l << "."
+ << std::endl;
+ return;
+}
+
+long double MoGEM::Loglikelihood(const arma::mat& data,
+ const std::vector<arma::vec>& means_l,
+ const std::vector<arma::mat>& covariances_l,
+ const arma::vec& weights_l) const
+{
+ long double loglikelihood = 0;
+ long double likelihood;
+
+ for (size_t j = 0; j < data.n_cols; j++)
+ {
+ likelihood = 0;
+ for(size_t i = 0; i < gaussians; i++)
+ likelihood += weights_l(i) * phi(data.unsafe_col(j), means_l[i],
+ covariances_l[i]);
+
+ loglikelihood += log(likelihood);
+ }
+
+ return loglikelihood;
+}
Copied: mlpack/trunk/src/mlpack/methods/gmm/gmm.hpp (from rev 10304, mlpack/trunk/src/mlpack/methods/gmm/mog_em.hpp)
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/gmm.hpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/gmm/gmm.hpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -0,0 +1,149 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_em.h
+ *
+ * Defines a Gaussian Mixture model and
+ * estimates the parameters of the model
+ */
+#ifndef __MLPACK_METHODS_MOG_MOG_EM_HPP
+#define __MLPACK_METHODS_MOG_MOG_EM_HPP
+
+#include <mlpack/core.h>
+
+PARAM_MODULE("mog", "Parameters for the Gaussian mixture model.");
+
+PARAM_INT("k", "The number of Gaussians in the mixture model (defaults to 1).",
+ "mog", 1);
+PARAM_INT("d", "The number of dimensions of the data on which the mixture "
+ "model is to be fit.", "mog", 0);
+
+namespace mlpack {
+namespace gmm {
+
+/**
+ * A Gaussian mixture model class.
+ *
+ * This class uses maximum likelihood loss functions to
+ * estimate the parameters of a gaussian mixture
+ * model on a given data via the EM algorithm.
+ *
+ *
+ * Example use:
+ *
+ * @code
+ * MoGEM mog;
+ * ArrayList<double> results;
+ *
+ * mog.Init(number_of_gaussians, dimension);
+ * mog.ExpectationMaximization(data, &results, optim_flag);
+ * @endcode
+ */
+class MoGEM {
+ private:
+ //! The number of Gaussians in the model.
+ size_t gaussians;
+ //! The dimensionality of the model.
+ size_t dimension;
+ //! Vector of means; one for each Gaussian.
+ std::vector<arma::vec> means;
+ //! Vector of covariances; one for each Gaussian.
+ std::vector<arma::mat> covariances;
+ //! Vector of a priori weights for each Gaussian.
+ arma::vec weights;
+
+ public:
+ /**
+ * Create a GMM with the given number of Gaussians, each of which have the
+ * specified dimensionality.
+ *
+ * @param gaussians Number of Gaussians in this GMM.
+ * @param dimension Dimensionality of each Gaussian.
+ */
+ MoGEM(size_t gaussians, size_t dimension) :
+ gaussians(gaussians),
+ dimension(dimension),
+ means(gaussians),
+ covariances(gaussians) { /* nothing to do */ }
+
+ //! Return the number of gaussians in the model.
+ const size_t Gaussians() { return gaussians; }
+
+ //! Return the dimensionality of the model.
+ const size_t Dimension() { return dimension; }
+
+ //! Return a const reference to the vector of means (mu).
+ const std::vector<arma::vec>& Means() const { return means; }
+ //! Return a reference to the vector of means (mu).
+ std::vector<arma::vec>& Means() { return means; }
+
+ //! Return a const reference to the vector of covariance matrices (sigma).
+ const std::vector<arma::mat>& Covariances() const { return covariances; }
+ //! Return a reference to the vector of covariance matrices (sigma).
+ std::vector<arma::mat>& Covariances() { return covariances; }
+
+ //! Return a const reference to the a priori weights of each Gaussian.
+ const arma::vec& Weights() const { return weights; }
+ //! Return a reference to the a priori weights of each Gaussian.
+ arma::vec& Weights() { return weights; }
+
+ /**
+ * This function outputs the parameters of the model
+ * to an arraylist of doubles
+ *
+ * @code
+ * ArrayList<double> results;
+ * mog.OutputResults(&results);
+ * @endcode
+ */
+ void OutputResults(std::vector<double>& results)
+ {
+ // Initialize the size of the output array
+ results.resize(gaussians * (1 + dimension * (1 + dimension)));
+
+ // Copy values to the array from the private variables of the class
+ for (size_t i = 0; i < gaussians; i++)
+ {
+ results[i] = weights[i];
+ for (size_t j = 0; j < dimension; j++)
+ {
+ results[gaussians + (i * dimension) + j] = (means[i])[j];
+ for (size_t k = 0; k < dimension; k++)
+ {
+ results[gaussians * (1 + dimension) +
+ (i * dimension * dimension) + (j * dimension) + k] =
+ (covariances[i])(j, k);
+ }
+ }
+ }
+ }
+
+ /**
+ * This function calculates the parameters of the model
+ * using the Maximum Likelihood function via the
+ * Expectation Maximization (EM) Algorithm.
+ *
+ * @code
+ * MoG mog;
+ * Matrix data = "the data on which you want to fit the model";
+ * ArrayList<double> results;
+ * mog.ExpectationMaximization(data, &results);
+ * @endcode
+ */
+ void ExpectationMaximization(const arma::mat& data_points);
+
+ /**
+ * This function computes the loglikelihood of model.
+ * This function is used by the 'ExpectationMaximization'
+ * function.
+ *
+ */
+ long double Loglikelihood(const arma::mat& data_points,
+ const std::vector<arma::vec>& means,
+ const std::vector<arma::mat>& covars,
+ const arma::vec& weights) const;
+};
+
+}; // namespace gmm
+}; // namespace mlpack
+
+#endif
Copied: mlpack/trunk/src/mlpack/methods/gmm/gmm_main.cpp (from rev 10304, mlpack/trunk/src/mlpack/methods/gmm/mog_em_main.cpp)
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/gmm_main.cpp (rev 0)
+++ mlpack/trunk/src/mlpack/methods/gmm/gmm_main.cpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -0,0 +1,68 @@
+/**
+ * @author Parikshit Ram (pram at cc.gatech.edu)
+ * @file mog_em_main.cpp
+ *
+ * This program test drives the parametric estimation
+ * of a Gaussian Mixture model using maximum likelihood.
+ *
+ * PARAMETERS TO BE INPUT:
+ *
+ * --data
+ * This is the file that contains the data on which
+ * the model is to be fit
+ *
+ * --mog_em/K
+ * This is the number of gaussians we want to fit
+ * on the data, defaults to '1'
+ *
+ * --output
+ * This file will contain the parameters estimated,
+ * defaults to 'ouotput.csv'
+ *
+ */
+#include "mog_em.hpp"
+
+PROGRAM_INFO("Mixture of Gaussians",
+ "This program takes a parametric estimate of a Gaussian mixture model (GMM)"
+ " using maximum likelihood.", "mog");
+
+PARAM_STRING_REQ("data", "A file containing the data on which the model has to "
+ "be fit.", "mog");
+PARAM_STRING("output", "The file into which the output is to be written into.",
+ "mog", "output.csv");
+
+using namespace mlpack;
+using namespace mlpack::gmm;
+
+int main(int argc, char* argv[]) {
+ CLI::ParseCommandLine(argc, argv);
+
+ ////// READING PARAMETERS AND LOADING DATA //////
+ arma::mat data_points;
+ data::Load(CLI::GetParam<std::string>("mog/data").c_str(), data_points, true);
+
+ ////// MIXTURE OF GAUSSIANS USING EM //////
+ MoGEM mog;
+
+ CLI::GetParam<int>("mog/k") = 1;
+ CLI::GetParam<int>("mog/d") = data_points.n_rows;
+
+ ////// Timing the initialization of the mixture model //////
+ Timers::StartTimer("mog/model_init");
+ mog.Init(1, data_points.n_rows);
+ Timers::StopTimer("mog/model_init");
+
+ ////// Computing the parameters of the model using the EM algorithm //////
+ std::vector<double> results;
+
+ Timers::StartTimer("mog/EM");
+ mog.ExpectationMaximization(data_points);
+ Timers::StopTimer("mog/EM");
+
+ mog.Display();
+ mog.OutputResults(results);
+
+ ////// OUTPUT RESULTS //////
+ // We need a better solution for this. So, currently, we do nothing.
+ // XML is probably the right tool for the job.
+}
Deleted: mlpack/trunk/src/mlpack/methods/gmm/mog_em.cpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/mog_em.cpp 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/mog_em.cpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -1,160 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_em.cpp
- *
- * Implementation for the loglikelihood function, the EM algorithm
- * and also computes the K-means for getting an initial point
- *
- */
-#include "mog_em.hpp"
-#include "phi.hpp"
-#include "kmeans.hpp"
-
-using namespace mlpack;
-using namespace gmm;
-
-void MoGEM::ExpectationMaximization(const arma::mat& data)
-{
- // Create temporary models and set to the right size.
- std::vector<arma::vec> means_trial(gaussians, arma::vec(dimension));
- std::vector<arma::mat> covariances_trial(gaussians,
- arma::mat(dimension, dimension));
- arma::vec weights_trial(gaussians);
-
- arma::mat cond_prob(gaussians, data.n_cols);
-
- long double l, l_old, best_l, TINY = 1.0e-10;
-
- best_l = -DBL_MAX;
-
- // We will perform five trials, and then save the trial with the best result
- // as our trained model.
- for (size_t iteration = 0; iteration < 5; iteration++)
- {
- // Use k-means to find initial values for the parameters.
- KMeans(data, gaussians, means_trial, covariances_trial, weights_trial);
-
- Log::Warn << "K-Means results:" << std::endl;
- for (size_t i = 0; i < gaussians; i++)
- {
- Log::Warn << "Mean " << i << ":" << std::endl;
- Log::Warn << means_trial[i] << std::endl;
- Log::Warn << "Covariance " << i << ":" << std::endl;
- Log::Warn << covariances_trial[i] << std::endl;
- }
- Log::Warn << "Weights: " << std::endl << weights_trial << std::endl;
-
- // Calculate the log likelihood of the model.
- l = Loglikelihood(data, means_trial, covariances_trial, weights_trial);
-
- l_old = -DBL_MAX;
-
- // Iterate to update the model until no more improvement is found.
- size_t max_iterations = 1000;
- size_t iteration = 0;
- while (std::abs(l - l_old) > TINY && iteration < max_iterations)
- {
- Log::Warn << "Iteration " << iteration << std::endl;
- // Calculate the conditional probabilities of choosing a particular
- // Gaussian given the data and the present theta value.
- for (size_t j = 0; j < data.n_cols; j++)
- {
- for (size_t i = 0; i < gaussians; i++)
- {
- cond_prob(i, j) = phi(data.unsafe_col(j), means_trial[i],
- covariances_trial[i]) * weights_trial[i];
- }
-
- // Normalize column to have sum probability of one.
- cond_prob.col(j) /= arma::sum(cond_prob.col(j));
- }
-
- // Store the sums of each row because they are used multiple times.
- arma::vec prob_row_sums = arma::sum(cond_prob, 1 /* row-wise */);
-
- // Calculate the new value of the means using the updated conditional
- // probabilities.
- for (size_t i = 0; i < gaussians; i++)
- {
- means_trial[i].zeros();
- for (size_t j = 0; j < data.n_cols; j++)
- means_trial[i] += cond_prob(i, j) * data.col(j);
-
- means_trial[i] /= prob_row_sums[i];
- }
-
- // Calculate the new value of the covariances using the updated
- // conditional probabilities and the updated means.
- for (size_t i = 0; i < gaussians; i++)
- {
- covariances_trial[i].zeros();
- for (size_t j = 0; j < data.n_cols; j++)
- {
- arma::vec tmp = data.col(j) - means_trial[i];
- covariances_trial[i] += cond_prob(i, j) * (tmp * trans(tmp));
- }
-
- covariances_trial[i] /= prob_row_sums[i];
- }
-
- // Calculate the new values for omega using the updated conditional
- // probabilities.
- weights_trial = prob_row_sums / data.n_cols;
-/*
- Log::Warn << "Estimated weights:" << std::endl << weights_trial
- << std::endl;
-
- for (size_t i = 0; i < gaussians; i++)
- {
-// Log::Warn << "Estimated mean " << i << ":" << std::endl;
-// Log::Warn << means_trial[i] << std::endl;
- Log::Warn << "Estimated covariance " << i << ":" << std::endl;
- Log::Warn << covariances_trial[i] << std::endl;
- }
-*/
-
- // Update values of l; calculate new log-likelihood.
- l_old = l;
- l = Loglikelihood(data, means_trial, covariances_trial, weights_trial);
-
- Log::Warn << "Improved log likelihood to " << l << std::endl;
-
- iteration++;
- }
-
- // The trial model is trained. Is it better than our existing model?
- if (l > best_l)
- {
- best_l = l;
-
- means = means_trial;
- covariances = covariances_trial;
- weights = weights_trial;
- }
- }
-
- Log::Info << "Log likelihood value of the estimated model: " << best_l << "."
- << std::endl;
- return;
-}
-
-long double MoGEM::Loglikelihood(const arma::mat& data,
- const std::vector<arma::vec>& means_l,
- const std::vector<arma::mat>& covariances_l,
- const arma::vec& weights_l) const
-{
- long double loglikelihood = 0;
- long double likelihood;
-
- for (size_t j = 0; j < data.n_cols; j++)
- {
- likelihood = 0;
- for(size_t i = 0; i < gaussians; i++)
- likelihood += weights_l(i) * phi(data.unsafe_col(j), means_l[i],
- covariances_l[i]);
-
- loglikelihood += log(likelihood);
- }
-
- return loglikelihood;
-}
Deleted: mlpack/trunk/src/mlpack/methods/gmm/mog_em.hpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/mog_em.hpp 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/mog_em.hpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -1,149 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_em.h
- *
- * Defines a Gaussian Mixture model and
- * estimates the parameters of the model
- */
-#ifndef __MLPACK_METHODS_MOG_MOG_EM_HPP
-#define __MLPACK_METHODS_MOG_MOG_EM_HPP
-
-#include <mlpack/core.h>
-
-PARAM_MODULE("mog", "Parameters for the Gaussian mixture model.");
-
-PARAM_INT("k", "The number of Gaussians in the mixture model (defaults to 1).",
- "mog", 1);
-PARAM_INT("d", "The number of dimensions of the data on which the mixture "
- "model is to be fit.", "mog", 0);
-
-namespace mlpack {
-namespace gmm {
-
-/**
- * A Gaussian mixture model class.
- *
- * This class uses maximum likelihood loss functions to
- * estimate the parameters of a gaussian mixture
- * model on a given data via the EM algorithm.
- *
- *
- * Example use:
- *
- * @code
- * MoGEM mog;
- * ArrayList<double> results;
- *
- * mog.Init(number_of_gaussians, dimension);
- * mog.ExpectationMaximization(data, &results, optim_flag);
- * @endcode
- */
-class MoGEM {
- private:
- //! The number of Gaussians in the model.
- size_t gaussians;
- //! The dimensionality of the model.
- size_t dimension;
- //! Vector of means; one for each Gaussian.
- std::vector<arma::vec> means;
- //! Vector of covariances; one for each Gaussian.
- std::vector<arma::mat> covariances;
- //! Vector of a priori weights for each Gaussian.
- arma::vec weights;
-
- public:
- /**
- * Create a GMM with the given number of Gaussians, each of which have the
- * specified dimensionality.
- *
- * @param gaussians Number of Gaussians in this GMM.
- * @param dimension Dimensionality of each Gaussian.
- */
- MoGEM(size_t gaussians, size_t dimension) :
- gaussians(gaussians),
- dimension(dimension),
- means(gaussians),
- covariances(gaussians) { /* nothing to do */ }
-
- //! Return the number of gaussians in the model.
- const size_t Gaussians() { return gaussians; }
-
- //! Return the dimensionality of the model.
- const size_t Dimension() { return dimension; }
-
- //! Return a const reference to the vector of means (mu).
- const std::vector<arma::vec>& Means() const { return means; }
- //! Return a reference to the vector of means (mu).
- std::vector<arma::vec>& Means() { return means; }
-
- //! Return a const reference to the vector of covariance matrices (sigma).
- const std::vector<arma::mat>& Covariances() const { return covariances; }
- //! Return a reference to the vector of covariance matrices (sigma).
- std::vector<arma::mat>& Covariances() { return covariances; }
-
- //! Return a const reference to the a priori weights of each Gaussian.
- const arma::vec& Weights() const { return weights; }
- //! Return a reference to the a priori weights of each Gaussian.
- arma::vec& Weights() { return weights; }
-
- /**
- * This function outputs the parameters of the model
- * to an arraylist of doubles
- *
- * @code
- * ArrayList<double> results;
- * mog.OutputResults(&results);
- * @endcode
- */
- void OutputResults(std::vector<double>& results)
- {
- // Initialize the size of the output array
- results.resize(gaussians * (1 + dimension * (1 + dimension)));
-
- // Copy values to the array from the private variables of the class
- for (size_t i = 0; i < gaussians; i++)
- {
- results[i] = weights[i];
- for (size_t j = 0; j < dimension; j++)
- {
- results[gaussians + (i * dimension) + j] = (means[i])[j];
- for (size_t k = 0; k < dimension; k++)
- {
- results[gaussians * (1 + dimension) +
- (i * dimension * dimension) + (j * dimension) + k] =
- (covariances[i])(j, k);
- }
- }
- }
- }
-
- /**
- * This function calculates the parameters of the model
- * using the Maximum Likelihood function via the
- * Expectation Maximization (EM) Algorithm.
- *
- * @code
- * MoG mog;
- * Matrix data = "the data on which you want to fit the model";
- * ArrayList<double> results;
- * mog.ExpectationMaximization(data, &results);
- * @endcode
- */
- void ExpectationMaximization(const arma::mat& data_points);
-
- /**
- * This function computes the loglikelihood of model.
- * This function is used by the 'ExpectationMaximization'
- * function.
- *
- */
- long double Loglikelihood(const arma::mat& data_points,
- const std::vector<arma::vec>& means,
- const std::vector<arma::mat>& covars,
- const arma::vec& weights) const;
-};
-
-}; // namespace gmm
-}; // namespace mlpack
-
-#endif
Deleted: mlpack/trunk/src/mlpack/methods/gmm/mog_em_main.cpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/mog_em_main.cpp 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/mog_em_main.cpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -1,68 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_em_main.cpp
- *
- * This program test drives the parametric estimation
- * of a Gaussian Mixture model using maximum likelihood.
- *
- * PARAMETERS TO BE INPUT:
- *
- * --data
- * This is the file that contains the data on which
- * the model is to be fit
- *
- * --mog_em/K
- * This is the number of gaussians we want to fit
- * on the data, defaults to '1'
- *
- * --output
- * This file will contain the parameters estimated,
- * defaults to 'ouotput.csv'
- *
- */
-#include "mog_em.hpp"
-
-PROGRAM_INFO("Mixture of Gaussians",
- "This program takes a parametric estimate of a Gaussian mixture model (GMM)"
- " using maximum likelihood.", "mog");
-
-PARAM_STRING_REQ("data", "A file containing the data on which the model has to "
- "be fit.", "mog");
-PARAM_STRING("output", "The file into which the output is to be written into.",
- "mog", "output.csv");
-
-using namespace mlpack;
-using namespace mlpack::gmm;
-
-int main(int argc, char* argv[]) {
- CLI::ParseCommandLine(argc, argv);
-
- ////// READING PARAMETERS AND LOADING DATA //////
- arma::mat data_points;
- data::Load(CLI::GetParam<std::string>("mog/data").c_str(), data_points, true);
-
- ////// MIXTURE OF GAUSSIANS USING EM //////
- MoGEM mog;
-
- CLI::GetParam<int>("mog/k") = 1;
- CLI::GetParam<int>("mog/d") = data_points.n_rows;
-
- ////// Timing the initialization of the mixture model //////
- Timers::StartTimer("mog/model_init");
- mog.Init(1, data_points.n_rows);
- Timers::StopTimer("mog/model_init");
-
- ////// Computing the parameters of the model using the EM algorithm //////
- std::vector<double> results;
-
- Timers::StartTimer("mog/EM");
- mog.ExpectationMaximization(data_points);
- Timers::StopTimer("mog/EM");
-
- mog.Display();
- mog.OutputResults(results);
-
- ////// OUTPUT RESULTS //////
- // We need a better solution for this. So, currently, we do nothing.
- // XML is probably the right tool for the job.
-}
Deleted: mlpack/trunk/src/mlpack/methods/gmm/mog_l2e.cpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/mog_l2e.cpp 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/mog_l2e.cpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -1,319 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_l2e.cpp
- *
- * Implementation for L2 loss function, and also some initial points generator.
- */
-#include "mog_l2e.hpp"
-#include "phi.hpp"
-#include "kmeans.hpp"
-
-using namespace mlpack;
-using namespace gmm;
-
-long double MoGL2E::L2Error(const arma::mat& data)
-{
- return RegularizationTerm_() - (2 * GoodnessOfFitTerm_(data)) / data.n_colse;
-}
-
-long double MoGL2E::L2Error(const arma::mat& data, arma::vec& gradients)
-{
- arma::vec g_reg, g_fit;
-
- long double l2e = RegularizationTerm_(g_reg) -
- (2 * GoodnessOfFitTerm_(data, g_fit)) / data.n_cols;
-
- gradients = g_reg - (2 * g_fit) / data.n_cols;
-
- return l2e;
-}
-
-long double MoGL2E::RegularizationTerm_()
-{
- arma::mat phi_mu(gaussians, gaussians)
-
- // Fill the phi_mu matrix (which is symmetric). Each entry of the matrix is
- // the phi() function evaluated on the difference between the means of each
- // Gaussian using the sum of the covariances of the two mixtures.
- for (size_t k = 1; k < gaussians; k++)
- {
- for (size_t j = 0; j < k; j++)
- {
- long double tmpVal = phi(means[k], means[j], covariances[k] +
- covariances[j]);
-
- phi_mu(j, k) = tmpVal;
- phi_mu(k, j) = tmpVal;
- }
- }
-
- // Because the difference between the means is 0, save a little time by only
- // doing the part of the calculation we need to (instead of calling phi()).
- for(size_t k = 0; k < gaussians; k++)
- phi_mu(k, k) = pow(2 * M_PI, (double) means[k].n_elem / -2.0)
- * pow(det(2 * covariances[k]), -0.5);
-
- return dot(weights, weights * phi_mu);
-}
-
-long double MoGL2E::RegularizationTerm_(arma::vec& g_reg)
-{
- arma::mat phi_mu(gaussians, gaussians);
- arma::vec x, y;
- long double reg, tmpVal;
-
- arma::vec df_dw, g_omega;
-
- std::vector<arma::vec> g_mu(gaussians, arma::vec(dimension));
- std::vector<arma::vec> g_sigma(gaussians, arma::vec((dimension * (dimension
- + 1)) / 2));
-
- std::vector<std::vector<arma::vec> > dp_d_mu(gaussians,
- std::vector<arma::vec>(gaussians));
- std::vector<std::vector<arma::vec> > dp_d_sigma(gaussians,
- std::vector<arma::vec>(gaussians));
-
- x = weights;
-
- // Fill the phi_mu matrix (which is symmetric). Each entry of the matrix is
- // the phi() function evaluated on the difference between the means of each
- // Gaussian using the sum of the covariances of the two mixtures.
- for(size_t k = 1; k < gaussians; k++)
- {
- for(size_t j = 0; j < k; j++)
- {
- std::vector<arma::mat> tmp_d_cov(dimension * (dimension + 1));
- arma::vec tmp_dp_d_sigma;
-
- // We should find a way to avoid all this copying to set up for the call
- // to phi().
- for(size_t i = 0; i < (dimension * (dimension + 1) / 2); i++)
- {
- tmp_d_cov[i] = (covariancesGradients[k])[i];
- tmp_d_cov[(dimension * (dimension + 1) / 2) + i] =
- (covariancesGradients[j])[i];
- }
-
- tmpVal = phi(means[k], means[j], covariances[k] + covariances[j],
- tmp_d_cov, dp_d_mu[j][k], tmp_dp_d_sigma);
-
- phi_mu(j, k) = tmpVal;
- phi_mu(k, j) = tmpVal;
-
- dp_d_mu[k][j] = -dp_d_mu[j][k];
-
- dp_d_sigma[j][k] = tmp_dp_d_sigma.rows(0,
- (dimension * (dimension + 1) / 2) - 1);
- dp_d_sigma[k][j] = tmp_dp_d_sigma.rows((dimension * (dimension + 1) / 2),
- tmp_dp_d_sigma.n_rows - 1);
- }
- }
-
- // Fill the diagonal elements of the phi_mu matrix.
- for (size_t k = 0; k < gaussians; k++)
- {
- arma::vec junk; // This result is not needed.
- phi_mu(k, k) = phi(means[k], means[k], 2 * covariances[k],
- covariancesGradients[k], junk, dp_d_sigma[k][k]);
-
- dp_d_mu[k][k].zeros(dimension);
- }
-
- // Calculate the regularization term value.
- arma::vec y = weights * phi_mu;
- long double reg = dot(weights, y);
-
- // Calculate the g_omega value; a vector of size K - 1
- df_dw = 2.0 * y;
- g_omega = weightsGradients * df_dw;
-
- // Calculate the g_mu values; K vectors of size D
- for (size_t k = 0; k < gaussians; k++)
- {
- for (size_t j = 0; j < gaussians; j++)
- g_mu[k] = 2.0 * weights[k] * weights[j] * dp_d_mu[j][k];
-
- // Calculating the g_sigma values - K vectors of size D(D+1)/2
- for (size_t k = 0; k < gaussians; k++)
- {
- for (size_t j = 0; j < gaussians; j++)
- g_sigma[k] += x[k] * dp_d_sigma[j][k];
- g_sigma[k] *= 2.0 * x[k];
- }
-
- // Making the single gradient vector of size K*(D+1)*(D+2)/2 - 1
- arma::vec tmp_g_reg((gaussians * (dimension + 1) *
- (dimension * 2) / 2) - 1);
- size_t j = 0;
- for (size_t k = 0; k < g_omega.n_elem; k++)
- tmp_g_reg[k] = g_omega[k];
- j = g_omega.n_elem;
-
- for (size_t k = 0; k < gaussians; k++) {
- for (size_t i = 0; i < dimension; i++)
- tmp_g_reg[j + (k * dimension) + i] = (g_mu[k])[i];
-
- for(size_t i = 0; i < (dimension * (dimension + 1) / 2); i++) {
- tmp_g_reg[j + (gaussians * dimension)
- + k * (dimension * (dimension + 1) / 2)
- + i] = (g_sigma[k])[i];
- }
- }
-
- g_reg = tmp_g_reg;
- }
-
- return reg;
-}
-
-long double MoGL2E::GoodnessOfFitTerm_(const arma::mat& data)
-{
- long double fit;
- arma::mat phi_x(gaussians, data.n_cols);
- arma::vec identity_vector;
-
- identity_vector.ones(data.n_cols);
-
- for (size_t k = 0; k < gaussians; k++)
- for (size_t i = 0; i < data.n_cols; i++)
- phi_x(k, i) = phi(data.unsafe_col(i), means[k], covariances[k]);
-
- fit = dot(weights * phi_x, identity_vector);
-
- return fit;
-}
-
-long double MoGL2E::GoodnessOfFitTerm_(const arma::mat& data,
- arma::vec& g_fit)
-{
- long double fit;
- arma::mat phi_x(gaussians, data.n_cols);
- arma::vec weights_l, x, y, identity_vector;
- arma::vec g_omega, tmp_g_omega;
- std::vector<arma::vec> g_mu, g_sigma;
-
- weights_l = weights;
- x.set_size(data.n_rows);
- identity_vector.ones(data.n_cols);
-
- g_mu.resize(gaussians);
- g_sigma.resize(gaussians);
-
- for(size_t k = 0; k < gaussians; k++)
- {
- g_mu[k].zeros(dimension);
- g_sigma[k].zeros(dimension * (dimension + 1) / 2);
-
- for (size_t i = 0; i < data.n_cols; i++) {
- arma::vec tmp_g_mu, tmp_g_sigma;
- phi_x(k, i) = phi(data.unsafe_col(i), means[k], covariances[k],
- d_sigma_[k], tmp_g_mu, tmp_g_sigma);
-
- g_mu[k] += tmp_g_mu;
- g_sigma[k] = tmp_g_sigma;
- }
-
- g_mu[k] *= weights_l[k];
- g_sigma[k] *= weights_l[k];
- }
-
- fit = dot(weights_l * phi_x, identity_vector);
-
- // Calculating the g_omega
- tmp_g_omega = phi_x * identity_vector;
- g_omega = d_omega_ * tmp_g_omega;
-
- // Making the single gradient vector of size K*(D+1)*(D+2)/2
- arma::vec tmp_g_fit((gaussians * (dimension + 1) *
- (dimension * 2) / 2) - 1);
- size_t j = 0;
- for (size_t k = 0; k < g_omega.n_elem; k++)
- tmp_g_fit[k] = g_omega[k];
- j = g_omega.n_elem;
-
- for (size_t k = 0; k < gaussians; k++)
- {
- for (size_t i = 0; i < dimension; i++)
- tmp_g_fit[j + (k * dimension) + i] = (g_mu[k])[i];
-
- for (size_t i = 0; i < (dimension * (dimension + 1) / 2); i++)
- tmp_g_fit[j + gaussians * dimension
- + k * (dimension * (dimension + 1) / 2) + i] = (g_sigma[k])[i];
- }
-
- g_fit = tmp_g_fit;
-
- return fit;
-}
-
-void MoGL2E::MultiplePointsGenerator(arma::mat& points,
- const arma::mat& d,
- size_t number_of_components)
-{
- size_t i, j, x;
-
- for (i = 0; i < points.n_rows; i++)
- for (j = 0; j < points.n_cols - 1; j++)
- points(i, j) = (rand() % 20001) / 1000 - 10;
-
- for (i = 0; i < points.n_rows; i++)
- {
- for (j = 0; j < points.n_cols; j++)
- {
- arma::vec tmp_mu = d.col(rand() % d.n_cols);
- for (x = 0; x < d.n_rows; x++)
- points(i, number_of_components - 1 + (j * d.n_rows) + x) = tmp_mu[x];
- }
- }
-
- for (i = 0; i < points.n_rows; i++)
- for (j = 0; j < points.n_cols; j++)
- for (x = 0; x < (d.n_rows * (d.n_rows + 1) / 2); x++)
- points(i, (number_of_components * (d.n_rows + 1) - 1)
- + (j * (d.n_rows * (d.n_rows + 1) / 2)) + x) = (rand() % 501) / 100;
-
- return;
-}
-
-void MoGL2E::InitialPointGenerator(arma::vec& theta,
- const arma::mat& data,
- size_t k_comp)
-{
- std::vector<arma::vec> means_l;
- std::vector<arma::mat> covars;
- arma::vec weights_l;
- double noise;
-
- weights_l.set_size(k_comp);
- means_l.resize(k_comp);
- covars.resize(k_comp);
-
- theta.set_size(k_comp);
-
- for (size_t i = 0; i < k_comp; i++)
- {
- means_l[i].set_size(data.n_rows);
- covars[i].set_size(data.n_rows, data.n_rows);
- }
-
- KMeans(data, k_comp, means_l, covars, weights_l);
-
- for (size_t k = 0; k < k_comp - 1; k++)
- {
- noise = (double) (rand() % 10000) / (double) 1000;
- theta[k] = noise - 5;
- }
-
- for (size_t k = 0; k < k_comp; k++)
- {
- for (size_t j = 0; j < data.n_rows; j++)
- theta[k_comp - 1 + k * data.n_rows + j] = (means_l[k])[j];
-
- arma::mat u = chol(covars[k]);
- for(size_t j = 0; j < data.n_rows; j++)
- for(size_t i = 0; i < j + 1; i++)
- theta[k_comp - 1 + (k_comp * data.n_rows)
- + (k * data.n_rows * (data.n_rows + 1) / 2)
- + (j * (j + 1) / 2 + i)] = u(i, j) + ((rand() % 501) / 100);
- }
-}
Deleted: mlpack/trunk/src/mlpack/methods/gmm/mog_l2e.hpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/mog_l2e.hpp 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/mog_l2e.hpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -1,397 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_l2e.hpp
- *
- * Defines a Gaussian Mixture model and estimates the parameters of the model.
- */
-#ifndef __MLPACK_METHODS_MOG_MOG_L2E_HPP
-#define __MLPACK_METHODS_MOG_MOG_L2E_HPP
-
-#include <mlpack/core.h>
-
-PARAM_MODULE("mog_l2e", "Parameters for the Gaussian mixture model.");
-
-PARAM_INT("k", "The number of Gaussians in the mixture model (defaults to 1).",
- "mog_l2e", 1);
-PARAM_INT("d", "The number of dimensions of the data on which the mixture "
- "model is to be fit.", "mog_l2e", 0);
-
-namespace mlpack {
-namespace gmm {
-
-/**
- * A Gaussian mixture model class.
- *
- * This class uses L2 loss function to estimate the parameters of a gaussian
- * mixture model on a given data.
- *
- * The parameters are converted for optimization to maintain the following
- * facts:
- * - the weights sum to one; for this, the weights were parameterized using
- * the logistic function
- * - the covariance matrix is always positive definite; for this, the Cholesky
- * decomposition is used
- *
- * Example use:
- *
- * @code
- * MoGL2E mog;
- * ArrayList<double> results;
- * double *params;
- *
- * mog.MakeModel(number_of_gaussians, dimension, params);
- * mog.L2Error(data);
- * mog.OutputResults(&results);
- * @endcode
- */
-class MoGL2E {
- private:
- // The parameters of the Mixture model
- size_t gaussians;
- size_t dimension;
- std::vector<arma::vec> means;
- std::vector<arma::mat> covariances;
- arma::vec weights;
-
- // The differential for the paramterization
- // for optimization
- arma::mat weightsGradients;
- std::vector<std::vector<arma::mat> > covariancesGradients;
-
- public:
-
- MoGL2E(size_t gaussians, size_t dimension) :
- gaussians(gaussians),
- dimension(dimension),
- means(gaussians),
- covariances(gaussians) { /* nothing to do */ }
-
- void Resize_d_sigma_()
- {
- d_sigma_.resize(gaussians);
- for(size_t i = 0; i < gaussians; i++)
- d_sigma_[i].resize(dimension * (dimension + 1) / 2);
- }
-
- /**
- * This function uses the parameters used for optimization
- * and converts it into athe parameters of a Gaussian
- * mixture model. This is to be used when you do not want
- * the gradient values.
- *
- * Example use:
- *
- * @code
- * MoGL2E mog;
- * mog.MakeModel(number_of_gaussians, dimension,
- * parameters_for_optimization);
- * @endcode
- */
-
- void MakeModel(size_t num_mods, size_t dimension, const arma::vec& theta) {
- arma::vec temp_mu(dimension);
- arma::mat lower_triangle_matrix;
- double sum, s_min = 0.01;
-
- gaussians = num_mods;
- this->dimension = dimension;
- lower_triangle_matrix.set_size(dimension, dimension);
-
- // calculating the omega values
- arma::vec temp_array = exp(theta);
- temp_array[num_mods - 1] = 1;
- sum = accu(temp_array);
-
- temp_array /= sum;
- weights = temp_array;
-
- // calculating the mu values
- for(size_t k = 0; k < gaussians; k++) {
- for(size_t j = 0; j < dimension; j++) {
- temp_mu[j] = theta[num_mods + (k * dimension) + j - 1];
- }
- means[k] = temp_mu;
- }
-
- // calculating the sigma values
- // using a lower triangular matrix and its transpose
- // to obtain a positive definite symmetric matrix
- arma::mat sigma_temp(dimension, dimension);
- for(size_t k = 0; k < gaussians; k++) {
- lower_triangle_matrix.zeros();
- for(size_t j = 0; j < dimension; j++) {
- for(size_t i = 0; i < j; i++) {
- lower_triangle_matrix(j, i) = theta[(gaussians - 1)
- + (gaussians * dimension) + k * (dimension * (dimension + 1) / 2)
- + (j * (j + 1) / 2) + i];
- }
- lower_triangle_matrix(j, j) = theta[(gaussians - 1)
- + (gaussians * dimension) + k * (dimension * (dimension + 1) / 2)
- + (j * (j + 1) / 2) + j] + s_min;
- }
- sigma_temp = lower_triangle_matrix * trans(lower_triangle_matrix);
- covariances[k] = sigma_temp;
- }
- }
-
- /**
- * This function uses the parameters used for optimization
- * and converts it into the parameters of a Gaussian
- * mixture model. This is to be used when you want
- * the gradient values.
- *
- * Example use:
- *
- * @code
- * MoGL2E mog;
- * mog.MakeModelWithGradients(number_of_gaussians, dimension,
- * parameters_for_optimization);
- * @endcode
- */
- void MakeModelWithGradients(size_t num_mods,
- size_t dimension,
- const arma::vec& theta) {
- arma::vec temp_mu(dimension);
- arma::mat lower_triangle_matrix;
- double sum, s_min = 0.01;
-
- gaussians = num_mods;
- this->dimension = dimension;
- lower_triangle_matrix.set_size(dimension, dimension);
-
- // calculating the omega values
- arma::vec temp_array = exp(theta);
- temp_array[num_mods - 1] = 1;
- sum = accu(temp_array);
-
- temp_array /= sum;
- weights = temp_array;
-
- // calculating the d_omega values
- arma::mat d_omega_temp(num_mods - 1, num_mods);
- d_omega_temp.zeros();
- for (size_t i = 0; i < num_mods - 1; i++) {
- for (size_t j = 0; j < i; j++) {
- d_omega_temp(i, j) = -(weights[i] * weights[j]);
- d_omega_temp(j, i) = -(weights[i] * weights[j]);
- }
- d_omega_temp(i, i) = weights[i] * (1 - weights[i]);
- }
-
- for (size_t i = 0; i < num_mods - 1; i++)
- d_omega_temp(i, num_mods - 1) = -(weights[i] * weights[num_mods - 1]);
-
- set_d_omega(d_omega_temp);
-
- // calculating the mu values
- for (size_t k = 0; k < num_mods; k++) {
- for (size_t j = 0; j < dimension; j++)
- temp_mu[j] = theta[num_mods + (k * dimension) + j - 1];
- means[k] = temp_mu;
- }
- // d_mu is not computed because it is implicitly known
- // since no parameterization is applied on them
-
- // using a lower triangular matrix and its transpose
- // to obtain a positive definite symmetric matrix
-
- // initializing the d_sigma values
-
- arma::mat d_sigma_temp(dimension, dimension);
- Resize_d_sigma_();
-
- // calculating the sigma values
- arma::mat sigma_temp(dimension, dimension);
- for (size_t k = 0; k < gaussians; k++) {
- lower_triangle_matrix.zeros();
- for (size_t j = 0; j < dimension; j++) {
- for (size_t i = 0; i < j; i++) {
- lower_triangle_matrix(j, i) = theta[(gaussians - 1)
- + (gaussians * dimension) + k * (dimension * (dimension + 1) / 2)
- + (j * (j + 1) / 2) + i];
- }
- lower_triangle_matrix(j, j) = theta[(gaussians - 1)
- + (gaussians * dimension) + k * (dimension * (dimension + 1) / 2)
- + (j * (j + 1) / 2) + j] + s_min;
- }
- sigma_temp = lower_triangle_matrix * trans(lower_triangle_matrix);
- covariances[k] = sigma_temp;
-
- // calculating the d_sigma values
- for (size_t i = 0; i < dimension; i++) {
- for (size_t in = 0; in < i + 1; in++) {
- arma::mat d_sigma_d_r(dimension, dimension);
- d_sigma_d_r.zeros();
- d_sigma_d_r(i, in) = 1.0;
-
- d_sigma_temp = d_sigma_d_r * trans(lower_triangle_matrix) +
- lower_triangle_matrix * trans(d_sigma_d_r);
- set_d_sigma(k, (i * (i + 1) / 2) + in, d_sigma_temp);
- }
- }
- }
- }
-
- arma::mat& d_omega() {
- return d_omega_;
- }
-
- std::vector<std::vector<arma::mat> >& d_sigma() {
- return d_sigma_;
- }
-
- std::vector<arma::mat>& d_sigma(size_t i) {
- return d_sigma_[i];
- }
-
- void set_d_omega(const arma::mat& d_omega) {
- d_omega_ = d_omega;
- }
-
- void set_d_sigma(size_t i, size_t j, const arma::mat& d_sigma_i_j) {
- d_sigma_[i][j] = d_sigma_i_j;
- }
-
- /**
- * This function outputs the parameters of the model
- * to an arraylist of doubles
- *
- * @code
- * ArrayList<double> results;
- * mog.OutputResults(&results);
- * @endcode
- */
- void OutputResults(std::vector<double>& results) {
-
- // Initialize the size of the output array
- results.resize(gaussians * (1 + dimension * (1 + dimension)));
-
- // Copy values to the array from the private variables of the class
- for (size_t i = 0; i < gaussians; i++) {
- results[i] = weights[i];
- for (size_t j = 0; j < dimension; j++) {
- results[gaussians + (i * dimension) + j] = (means[i])[j];
- for (size_t k = 0; k < dimension; k++) {
- results[gaussians * (1 + dimension)
- + (i * dimension * dimension) + (j * dimension)
- + k] = (covariances[i])(j, k);
- }
- }
- }
- }
-
- /**
- * This function calculates the L2 error and
- * the gradient of the error with respect to the
- * parameters given the data and the parameterized
- * mixture
- *
- * Example use:
- *
- * @code
- * const Matrix data;
- * MoGL2E mog;
- * size_t num_gauss, dimension;
- * double *params; // get the parameters
- *
- * mog.MakeModel(num_gauss, dimension, params);
- * mog.L2Error(data);
- * @endcode
- */
- long double L2Error(const arma::mat& data);
- long double L2Error(const arma::mat& data, arma::vec& gradients);
-
- /**
- * Calculates the regularization value for a
- * Gaussian mixture and its gradient with
- * respect to the parameters
- *
- * Used by the 'L2Error' function to calculate
- * the regularization part of the error
- */
- long double RegularizationTerm_();
- long double RegularizationTerm_(arma::vec& g_reg);
-
- /**
- * Calculates the goodness-of-fit value for a
- * Gaussian mixture and its gradient with
- * respect to the parameters
- *
- * Used by the 'L2Error' function to calculate
- * the goodness-of-fit part of the error
- */
- long double GoodnessOfFitTerm_(const arma::mat& data);
- long double GoodnessOfFitTerm_(const arma::mat& data, arma::vec& g_fit);
-
- /**
- * This function computes multiple number of starting points
- * required for the Nelder Mead method
- *
- * Example use:
- * @code
- * double **p;
- * size_t n, num_gauss;
- * const Matrix data;
- *
- * MoGL2E::MultiplePointsGeneratot(p, n, data, num_gauss);
- * @endcode
- */
- static void MultiplePointsGenerator(arma::mat& points,
- const arma::mat& d,
- size_t number_of_components);
-
- /**
- * This function parameterizes the starting point obtained
- * from the 'k_means" for optimization purposes using the
- * Quasi Newton method
- *
- * Example use:
- * @code
- * double *p;
- * size_t num_gauss;
- * const Matrix data;
- *
- * MoGL2E::InitialPointGeneratot(p, data, num_gauss);
- * @endcode
- */
- static void InitialPointGenerator(arma::vec& theta,
- const arma::mat& data,
- size_t k_comp);
-
- /**
- * This is the function which would be used for
- * optimization. It creates its own object of
- * class MoGL2E and returns the L2 error
- * and the gradient which are computed by
- * the functions of the class
- *
- */
- static long double L2ErrorForOpt(const arma::vec& params,
- const arma::mat& data) {
- size_t num_gauss = (params.n_elem + 1) * 2 /
- ((data.n_rows + 1) * (data.n_rows + 2));
-
- MoGL2E model(num_gauss, data.n_rows);
- model.MakeModel(num_gauss, data.n_rows, params);
-
- return model.L2Error(data);
- }
-
- static long double L2ErrorForOpt(const arma::vec& params,
- const arma::mat& data,
- arma::vec& gradient) {
-
- size_t num_gauss = (params.n_elem + 1) * 2 /
- ((data.n_rows + 1) * (data.n_rows + 2));
-
- MoGL2E model(num_gauss, data.n_rows);
- model.MakeModelWithGradients(num_gauss, data.n_rows, params);
-
- return model.L2Error(data, gradient);
- }
-};
-
-}; // namespace gmm
-}; // namespace mlpack
-
-#endif // __MLPACK_METHODS_MOG_MOG_L2E_HPP
Deleted: mlpack/trunk/src/mlpack/methods/gmm/mog_l2e_main.cpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/mog_l2e_main.cpp 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/mog_l2e_main.cpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -1,122 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file mog_l2e_main.cpp
- *
- * This program test drives the L2 estimation
- * of a Gaussian Mixture model.
- *
- * PARAMETERS TO BE INPUT:
- *
- * --data
- * This is the file that contains the data on which
- * the model is to be fit
- *
- * --mog_l2e/K
- * This is the number of gaussians we want to fit
- * on the data, defaults to '1'
- *
- * --output
- * This file will contain the parameters estimated,
- * defaults to 'output.csv'
- *
- */
-#include "mog_l2e.hpp"
-#include "optimizers.hpp"
-
-PROGRAM_INFO("Mixture of Gaussians",
- "This program takes a parametric estimate of a Gaussian mixture model (GMM)"
- " using the L2 loss function.", "mog_l2e");
-
-PARAM_STRING_REQ("data", "A file containing the data on which the model has to "
- "be fit.", "mog_l2e");
-PARAM_STRING("output", "The file into which the output is to be written into",
- "mog_l2e", "output.csv");
-
-using namespace mlpack;
-using namespace mlpack::gmm;
-
-int main(int argc, char* argv[]) {
- CLI::ParseCommandLine(argc, argv);
-
- ////// READING PARAMETERS AND LOADING DATA //////
- arma::mat data_points;
- data::Load(CLI::GetParam<std::string>("mog_l2e/data").c_str(), data_points,
- true);
-
- ////// MIXTURE OF GAUSSIANS USING L2 ESTIMATCLIN //////
- size_t number_of_gaussians = CLI::GetParam<int>("mog_l2e/k");
- CLI::GetParam<int>("mog_l2e/d") = data_points.n_rows;
- size_t dimension = CLI::GetParam<int>("mog_l2e/d");
-
- ////// RUNNING AN OPTIMIZER TO MINIMIZE THE L2 ERROR //////
- const char *opt_method = CLI::GetParam<std::string>("opt/method").c_str();
- size_t param_dim = (number_of_gaussians * (dimension + 1) * (dimension + 2)
- / 2 - 1);
- CLI::GetParam<int>("opt/param_space_dim") = param_dim;
-
- size_t optim_flag = (strcmp(opt_method, "NelderMead") == 0 ? 1 : 0);
- MoGL2E mog;
-
- if (optim_flag == 1) {
- ////// OPTIMIZER USING NELDER MEAD METHOD //////
- NelderMead opt;
-
- ////// Initializing the optimizer //////
- Timers::StartTimer("opt/init_opt");
- opt.Init(MoGL2E::L2ErrorForOpt, data_points);
- Timers::StopTimer("opt/init_opt");
-
- ////// Getting starting points for the optimization //////
- arma::mat pts(param_dim, param_dim + 1);
-
- Timers::StartTimer("opt/get_init_pts");
- MoGL2E::MultiplePointsGenerator(pts, data_points, number_of_gaussians);
- Timers::StopTimer("opt/get_init_pts");
-
- ////// The optimization //////
- Timers::StartTimer("opt/optimizing");
- opt.Eval(pts);
- Timers::StopTimer("opt/optimizing");
-
- ////// Making model with the optimal parameters //////
- // This is a stupid way to do it and putting the 0s there ensures it will
- // fail. Do it better!
- mog.MakeModel(0, 0, pts.col(0));
- } else {
- ////// OPTIMIZER USING QUASI NEWTON METHOD //////
- QuasiNewton opt;
-
- ////// Initializing the optimizer //////
- Timers::StartTimer("opt/init_opt");
- opt.Init(MoGL2E::L2ErrorForOpt, data_points);
- Timers::StopTimer("opt/init_opt");
-
- ////// Getting starting point for the optimization //////
- arma::vec pt(param_dim);
-
- Timers::StartTimer("opt/get_init_pt");
- MoGL2E::InitialPointGenerator(pt, data_points, number_of_gaussians);
- Timers::StopTimer("opt/get_init_pt");
-
- ////// The optimization //////
- Timers::StartTimer("opt/optimizing");
- opt.Eval(pt);
- Timers::StopTimer("opt/optimizing");
-
- ////// Making model with optimal parameters //////
- // This is a stupid way to do it and putting the 0s there ensures it will
- // fail. Do it better!
- mog.MakeModel(0, 0, pt);
-
- }
-
- long double error = mog.L2Error(data_points);
- Log::Info << "Minimum L2 error achieved: " << error << "." << std::endl;
- mog.Display();
-
- std::vector<double> results;
- mog.OutputResults(results);
-
- ////// OUTPUT RESULTS //////
- // We need a better way to do this (like XML). For now we just do nothing.
-}
Deleted: mlpack/trunk/src/mlpack/methods/gmm/optimizers.cpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/optimizers.cpp 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/optimizers.cpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -1,352 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file optimizers.cpp
- *
- * Implementation of the optimizers
- */
-#include <mlpack/core.h>
-#include "optimizers.hpp"
-
-using namespace mlpack;
-using namespace gmm;
-
-void NelderMead::Eval(arma::mat& pts) {
- size_t dim = dimension();
- size_t num_func_eval;
- size_t i, j, ihi, ilo, inhi;
- size_t mpts = dim + 1;
- double sum, swap;
- arma::vec psum;
- long double swap_y, rtol, ytry, TINY = 1.0e-10;
- arma::vec y;
- arma::vec param_passed;
-
- // Default value is 1e-5.
- long double tol = CLI::GetParam<double>("opt/tolerance");
-
- // Default value is 50000.
- size_t NMAX = CLI::GetParam<int>("opt/MAX_FUNC_EVAL");
-
- param_passed.set_size(dim);
- psum.set_size(dim);
- num_func_eval = 0;
- y.set_size(mpts);
- for(i = 0; i < mpts; i++) {
- param_passed = pts.row(i);
- y[i] = (*func_ptr_)(param_passed, data_);
- }
-
- for(;;) {
- ilo = 0;
- ihi = y[0] > y[1] ? (inhi = 1, 0) : (inhi = 0, 1);
- for (i = 0; i < mpts; i++) {
- if(y[i] <= y[ilo])
- ilo = i;
- if(y[i] > y[ihi]) {
- inhi = ihi;
- ihi = i;
- }
- else if((y[i] > y[inhi]) && (i != ihi))
- inhi = i;
- }
-
- rtol = 2.0 * fabs(y[ihi] - y[ilo]) / (fabs(y[ihi]) + fabs(y[ilo]) + TINY);
- if (rtol < tol) {
- swap_y = y[0];
- y[0] = y[ilo];
- y[ilo] = swap_y;
- for (i = 0; i < dim; i++) {
- swap = pts(0, i);
- pts(0, i) = pts(ilo, i);
- pts(ilo, i) = swap;
- }
-
- break;
- }
-
- if (num_func_eval > NMAX) {
- Log::Warn << "Nelder-Mead: Maximum number of function evaluations "
- "exceeded." << std::endl;
- break;
- }
-
- num_func_eval += 2;
-
- // Beginning a new iteration.
- // Extrapolating by a factor of -1.0 through the face of the simplex
- // across from the high point, i.e, reflect the simplex from the high point
- for (j = 0 ; j < dim; j++) {
- sum = 0.0;
- for(i = 0; i < mpts; i++)
- if (i != ihi)
- sum += pts(i, j);
-
- psum[j] = sum / dim;
- }
-
- ytry = ModSimplex_(pts, y, psum, ihi, -1.0);
-
- if(ytry <= y[ilo]) {
- // result better than best point
- // so additional extrapolation by a factor of 2
- ytry = ModSimplex_(pts, y, psum, ihi, 2.0);
- } else if(ytry >= y[ihi]) {
- // result worse than the worst point
- // so there is a lower intermediate point,
- // i.e., do a one dimensional contraction
- ytry = ModSimplex_(pts, y, psum, ihi, 0.5);
- if(ytry > y[ihi]) {
- // Can't get rid of the high point,
- // try to contract around the best point
- for (i = 0; i < mpts; i++) {
- if (i != ilo) {
- for (j = 0; j < dim; j++)
- pts(i, j) = psum[j] = 0.5 * (pts(i, j) + pts(ilo, j));
-
- param_passed = psum;
- y[i] = (*func_ptr_)(param_passed, data());
- }
- }
- num_func_eval += dim;
-
- for (j = 0 ; j < dim ; j++) {
- sum = 0.0;
- for (i = 0 ; i < mpts ; i++)
- if (i != ihi)
- sum += pts(i, j);
-
- psum[j] = sum / dim;
- }
- }
- }
- else
- --num_func_eval;
- }
-
- CLI::GetParam<int>("opt/func_evals") = num_func_eval;
-}
-
-long double NelderMead::ModSimplex_(arma::mat& pts, arma::vec& y,
- arma::vec& psum, size_t ihi,
- float fac) {
-
- size_t j, dim = dimension();
- long double ytry;
- arma::vec ptry(dim);
- arma::vec param_passed(dim);
-
- for (j = 0; j < dim; j++)
- ptry[j] = psum[j] * (1 - fac) + pts(ihi, j) * fac;
-
- param_passed = ptry;
- ytry = (*func_ptr_)(param_passed, data());
- if (ytry < y[ihi]) {
- y[ihi] = ytry;
- for (j = 0; j < dim; j++)
- pts(ihi, j) = ptry[j];
- }
-
- return ytry;
-}
-
-void QuasiNewton::Eval(arma::vec& pt) {
- size_t n = dimension(), iters;
- size_t i, its;
- // Default value is 200.
- size_t MAXIMUM_ITERATCLINS = CLI::GetParam<int>("opt/MAX_ITERS");
-
- long double temp_1, temp_2, temp_3, temp_4, f_previous, f_min,
- maximum_step_length, sum = 0.0, sumdg, sumxi, temp, test;
- arma::vec dgrad, grad, hdgrad, xi;
- arma::vec pold, pnew;
- arma::mat hessian;
-
- // Default value is 3.0e-8.
- double EPSILON = CLI::GetParam<double>("opt/EPSILON");
-
- // Default value is 1.0e-5.
- double TOLERANCE = CLI::GetParam<double>("optTOLERANCE");
-
- // Default value is 100.
- double MAX_STEP_SIZE = CLI::GetParam<double>("opt/MAX_STEP_SIZE");
-
- // Default value is 1.0e-7.
- double g_tol = CLI::GetParam<double>("opt/gtol");
-
- dgrad.set_size(n);
- grad.set_size(n);
- hdgrad.set_size(n);
- hessian.set_size(n, n);
- pnew.set_size(n);
- xi.set_size(n);
- pold = pt;
- f_previous = (*func_ptr_)(pold, data(), grad);
- arma::vec tmp;
- tmp.ones(n);
- hessian.diag() = tmp;
-
- xi = -1 * grad;
- sum = dot(pold, pold);
-
- double fmax;
- if(sqrt(sum) > (float) n)
- fmax = sqrt(sum);
- else
- fmax = (float) n;
-
- maximum_step_length = MAX_STEP_SIZE * fmax;
-
- for (its = 0; its < MAXIMUM_ITERATCLINS; its++) {
- dgrad = grad;
- LineSearch_(pold, f_previous, grad, xi,
- pnew, f_min, maximum_step_length);
- f_previous = f_min;
- xi = pold - pnew;
- pold = pnew;
- pt = pold;
-
- test = 0.0;
- for (i = 0; i < n; i++) {
- if (fabs(pold[i]) > 1.0)
- fmax = fabs(pold[i]);
- else
- fmax = 1.0;
-
- temp = fabs(xi[i]) / fmax;
- if (temp > test)
- test = temp;
- }
-
- if (test < TOLERANCE) {
- iters = its;
- CLI::GetParam<int>("opt/iters") = iters;
- return;
- }
-
- test = 0.0;
- if (f_min > 1.0)
- temp_1 = f_min;
- else
- temp_1 = 1.0;
-
- for (i = 0; i < n; i++) {
- if (fabs(pold[i]) > 1.0)
- fmax = pold[i];
- else
- fmax = 1.0;
-
- temp = fabs(grad[i]) * fmax / temp_1;
- if (temp > test)
- test = temp;
- }
-
- if (test < g_tol) {
- iters = its;
- CLI::GetParam<int>("opt/iters") = iters;
- return;
- }
-
- dgrad -= grad;
- dgrad *= -1.0;
- hdgrad = hessian * dgrad;
-
- temp_2 = dot(dgrad, xi);
- temp_4 = dot(dgrad, hdgrad);
- sumdg = dot(dgrad, dgrad);
- sumxi = dot(xi, xi);
-
- if (temp_2 > sqrt(EPSILON * sumdg * sumxi)) {
- temp_2 = 1.0 / temp_2;
- temp_3 = 1.0 / temp_4;
-
- dgrad = temp_2 * xi;
- dgrad -= temp_3 * hdgrad;
-
- hessian += temp_2 * (xi * trans(xi));
- hessian -= temp_3 * (hdgrad * trans(hdgrad));
- hessian += temp_4 * (dgrad * trans(dgrad));
- }
-
- xi = -hessian * grad;
- }
-
- Log::Warn << "Too many iterations in Quasi-Newton: giving up." << std::endl;
-}
-
-void QuasiNewton::LineSearch_(arma::vec& pold, long double fold,
- arma::vec& grad, arma::vec& xi,
- arma::vec& pnew, long double& f_min,
- long double maximum_step_length) {
-
- size_t i, n = dimension();
- long double a, step_length, previous_step_length = 0.0,
- minimum_step_length, b, disc, previous_f_value = 0.0,
- rhs1, rhs2, slope, sum, temp, test, temp_step_length,
- MIN_DECREASE = 1.0e-4, TOLERANCE = 1.0e-7;
-
- sum = sqrt(dot(xi, xi));
- if(sum > maximum_step_length)
- xi *= (maximum_step_length / sum);
-
- slope = dot(grad, xi);
- if (slope >= 0.0)
- return;
-
- test = 0.0;
- for (i = 0; i < n; i++) {
- double fmax;
- fmax = (fabs(pold[i]) > 1.0 ? fabs(pold[i]) : 1.0);
- temp = fabs(xi[i]) / fmax;
- if (temp > test)
- test = temp;
- }
-
- minimum_step_length = TOLERANCE / test;
- step_length = 1.0;
-
- for (;;) {
- pnew = pold + step_length * xi;
-
- f_min = (*func_ptr_)(pnew, data(), grad);
-
- if (step_length < minimum_step_length) {
- pnew = pold;
- return;
- } else if(f_min <= fold + MIN_DECREASE * step_length * slope) {
- return;
- } else {
- if (step_length == 1.0)
- temp_step_length = -slope / (2.0 * (f_min - fold - slope));
- else {
- rhs1 = f_min - fold - step_length * slope;
- rhs2 = previous_f_value - fold - previous_step_length * slope;
- a = (rhs1 / (step_length * step_length)
- - rhs2 / (previous_step_length * previous_step_length))
- / (step_length - previous_step_length);
- b = (-previous_step_length * rhs1 / (step_length * step_length) +
- step_length * rhs2 / (previous_step_length * previous_step_length))
- / (step_length - previous_step_length);
-
- if (a == 0.0)
- temp_step_length = -slope / (2.0*b);
- else {
- disc = b * b - 3.0 * a * slope;
- if(disc < 0.0)
- temp_step_length = 0.5 * step_length;
- else if (b <= 0.0)
- temp_step_length = (-b + sqrt(disc)) / (3.0 * a);
- else
- temp_step_length = -slope / (b + sqrt(disc));
- }
-
- if(temp_step_length > 0.5 * step_length)
- temp_step_length = 0.5 * step_length;
- }
- }
-
- previous_step_length = step_length;
- previous_f_value = f_min;
- step_length = (temp_step_length > 0.1 * step_length) ? temp_step_length
- : 0.1 * step_length;
- }
-}
Deleted: mlpack/trunk/src/mlpack/methods/gmm/optimizers.hpp
===================================================================
--- mlpack/trunk/src/mlpack/methods/gmm/optimizers.hpp 2011-11-16 22:52:16 UTC (rev 10304)
+++ mlpack/trunk/src/mlpack/methods/gmm/optimizers.hpp 2011-11-16 22:53:40 UTC (rev 10305)
@@ -1,153 +0,0 @@
-/**
- * @author Parikshit Ram (pram at cc.gatech.edu)
- * @file optimizers.hpp
- *
- * Declares classes for two types of optimizer
- *
- */
-#ifndef __MLPACK_METHODS_MOG_OPTIMIZERS_HPP
-#define __MLPACK_METHODS_MOG_OPTIMIZERS_HPP
-
-#include <mlpack/core.h>
-
-PARAM_STRING("method", "The method used to optimize", "opt", "");
-
-PARAM_INT_REQ("param_space_dim", "The dimension of the parameter space.", "opt");
-PARAM_INT("MAX_FUNC_EVAL", "The maximum number of function evaluations\
- allowed to the NelderMead optimizer (defaults to 50000)", "opt", 50000);
-
-PARAM_INT("func_evals", "The number of function evaluations taken by the algorithm", "opt", 0);
-PARAM_INT("MAX_ITERS", "The maximum number of iterations allowed to the function", "opt", 200);
-PARAM_INT("iters", "The number of iterations the algorithm actually went through", "opt", 0);
-
-PARAM_DOUBLE("EPSILON", "Value of epsilon.", "opt", 3.0e-8);
-PARAM_DOUBLE("TOLERANCE", "Tolerance for the minimum movement for the parameter value.", "opt", 1.0e-5);
-PARAM_DOUBLE("gtol", "Tolerance value for the gradient of the function", "opt", 1.0e-7);
-PARAM_DOUBLE("MAX_STEP_SIZE", "The maximum step size in the direction of the gradient.", "opt", 100.0);
-PARAM_DOUBLE("tolerance", "Undocumented parameter", "opt", 1.0e-5);
-
-PARAM_MODULE("opt", "This file contains two optimizers.");
-
-namespace mlpack {
-namespace gmm {
-
-/**
- * An optimizer using the Nelder Mead method,
- * also known as the polytope or the simplex
- * method.
- *
- * It does multivariate minimization of an
- * objective function. If it is optimizing in
- * 'd' dimensions, it would require 'd+1'
- * starting points.
- *
- * Example use:
- *
- * @code
- * double init_pts[d+1][d];
- * size_t number_of_function_evaluations;
- * struct datanode *opt_module = fx_submodule(NULL,"NelderMead","opt_module");
- * Matrix data;
- * size_t dim_param_space;
- *
- * ...
- * NelderMead opt;
- * opt.Init(obj_function, data, dim_param_space, opt_module);
- * ...
- * opt.Eval(init_pts);
- * // init_pts[0] contaings the optimal point found
- * @endcode
- *
- */
-class NelderMead {
- private:
- size_t dimension_;
- arma::mat data_;
- long double (*func_ptr_)(const arma::vec&, const arma::mat&);
-
- public:
- NelderMead() { }
- ~NelderMead() { }
-
- void Init(long double (*fun)(const arma::vec&, const arma::mat&),
- arma::mat& data) {
- data_ = data;
- func_ptr_ = fun;
- dimension_ = mlpack::CLI::GetParam<int>("opt/param_space_dim");
- }
-
- const arma::mat& data() {
- return data_;
- }
-
- size_t dimension() {
- return dimension_;
- }
-
- void Eval(arma::mat& pts);
- long double ModSimplex_(arma::mat& pts, arma::vec& y,
- arma::vec& psum, size_t ihi, float fac);
-};
-
-/**
- * An optimizer using the Quasi Newton method,
- * also known as the variable metrics
- * method.
- *
- * It does multivariate minimization of an
- * objective function using only the function
- * value and the gradients.
- *
- * Example use:
- *
- * @code
- * double init_pt[d];
- * size_t number_of_iters;
- * struct datanode *opt_module = fx_submodule(NULL,"QuasiNewton","opt_module");
- * Matrix data;
- * size_t dim_param_space;
- *
- * ...
- * QuasiNewton opt;
- * opt.Init(obj_function, data, dim_param_space, opt_module);
- * ...
- * opt.Eval(init_pt);
- * // init_pt contains the optimal point found
- * @endcode
- *
- */
-class QuasiNewton {
- private:
- size_t dimension_;
- arma::mat data_;
- long double (*func_ptr_)(const arma::vec&, const arma::mat&, arma::vec&);
-
- public:
- QuasiNewton(){ }
- ~QuasiNewton(){ }
-
- void Init(long double (*fun)(const arma::vec&, const arma::mat&, arma::vec&),
- arma::mat& data) {
- data_ = data;
- func_ptr_ = fun;
- dimension_ = mlpack::CLI::GetParam<int>("opt/param_space_dim");
- }
-
- const arma::mat& data() {
- return data_;
- }
-
- size_t dimension() {
- return dimension_;
- }
-
- void Eval(arma::vec& pt);
- void LineSearch_(arma::vec& pold, long double fold, arma::vec& grad,
- arma::vec& xi, arma::vec& pnew, long double& f_min,
- long double maximum_step_length);
-};
-
-}; // namespace gmm
-}; // namespace mlpack
-
-#endif // __MLPACK_METHODS_MOG_OPTIMIZERS_HPP
More information about the mlpack-svn
mailing list