diff --git a/doc/guide/cv.hpp b/doc/guide/cv.hpp new file mode 100644 index 0000000000..4c6f8f9478 --- /dev/null +++ b/doc/guide/cv.hpp @@ -0,0 +1,372 @@ +namespace mlpack { +namespace cv { + +/*! @page cv Cross-Validation + +@section cvintro Introduction + +@b mlpack implements cross-validation support for its learning algorithms, for a +variety of performance measures. Cross-validation is useful for determining an +estimate of how well the learner will generalize to un-seen test data. It is a +commonly used part of the data science pipeline. + +In short, given some learner and some performance measure, we wish to get an +average of the performance measure given different splits of the dataset into +training data and validation data. The learner is trained on the training data, +and the performance measure is evaluated on the validation data. + +mlpack currently implements two easy-to-use forms of cross-validation: + + - @b simple @b cross-validation, where we simply desire the performance measure + on a single split of the data into a training set and validation set + + - @b k-fold @b cross-validation, where we split the data k ways and desire the + average performance measure on each of the k splits of the data + +In this tutorial we will see the usage examples and details of the +cross-validation module. Because the cross-validation code is generic and can +be used with any learner and performance measure, any use of the +cross-validation code in mlpack has to be in C++. + +This tutorial is split into the following sections: + + - @ref cvbasic Simple cross-validation examples + - @ref cvbasic_ex_1 10-fold cross-validation on softmax regression + - @ref cvbasic_ex_2 10-fold cross-validation on weighted decision trees + - @ref cvbasic_ex_3 10-fold cross-validation with categorical decision trees + - @ref cvbasic_ex_4 Simple cross-validation for linear regression + - @ref cvbasic_metrics Performance measures + - @ref cvbasic_api The \c KFoldCV and \c SimpleCV classes + - @ref cvbasic_further Further reference + +@section cvbasic Simple cross-validation examples + +@subsection cvbasic_ex_1 10-fold cross-validation on softmax regression + +Suppose we have some data to train and validate on, as defined below: + +@code + // 100-point 6-dimensional random dataset. + arma::mat data = arma::randu(6, 100); + // Random labels in the [0, 4] interval. + arma::Row labels = + arma::randi>(100, arma::distr_param(0, 4)); + size_t numClasses = 5; +@endcode + +The code above generates an 100-point random 6-dimensional dataset with 5 +classes. + +To run 10-fold cross-validation for softmax regression with accuracy as a +performance measure, we can write the following piece of code. + +@code + KFoldCV cv(10, data, labels, numClasses); + double lambda = 0.1; + double softmaxAccuracy = cv.Evaluate(lambda); +@endcode + +Note that the \c Evaluate method of \c KFoldCV takes any hyperparameters of an +algorithm---that is, anything that is not \c data, \c labels, \c numClasses, +\c datasetInfo, or \c weights (those last three may not be present for every +algorithm type). To be more specific, in this example the \c Evaluate method +relies on the following \ref regression::SoftmaxRegression "SoftmaxRegression" +constructor: + +@code + template + SoftmaxRegression(const arma::mat& data, + const arma::Row& labels, + const size_t numClasses, + const double lambda = 0.0001, + const bool fitIntercept = false, + OptimizerType optimizer = OptimizerType()); +@endcode + +which has the parameter \c lambda after three conventional arguments (\c data, +\c labels and \c numClasses). We can skip passing \c fitIntercept and \c +optimizer since there are the default values. (Technically, we don't even need +to pass \c lambda since there is a default value.) + +In general to cross-validate you need to specify what machine learning algorithm +and metric you are going to use, and then to pass some conventional data-related +parameters into one of the cross-validation constructors and all other +parameters (which are generally hyperparameters) into the \c Evaluate method. + +@subsection cvbasic_ex_2 10-fold cross-validation on weighted decision trees + +In the following example we will cross-validate +\ref tree::DecisionTree "DecisionTree" with weights. This is very similar to +the previous example, except that we also have instance weights for each point +in the dataset. We can generate weights for the dataset from the previous +example with the code below: + +@code + // Random weights for every point from the code snippet above. + arma::rowvec weights = arma::randu(1, 100); +@endcode + +Given those weights for each point, we can now perform cross-validation by also +passing the weights to the constructor of \c KFoldCV: + +@code + KFoldCV, Accuracy> cv2(10, data, labels, numClasses, weights); + size_t minimumLeafSize = 8; + double weightedDecisionTreeAccuracy = cv2.Evaluate(minimumLeafSize); +@endcode + +As with the previous example, internally this call to \c cv2.Evaluate() relies +on the following \ref tree::DecisionTree "DecisionTree" constructor: + +@code + template + DecisionTree(MatType&& data, + LabelsType&& labels, + const size_t numClasses, + WeightsType&& weights, + const size_t minimumLeafSize = 10, + const std::enable_if_t::type>::value>* + = 0); +@endcode + +@subsection cvbasic_ex_3 10-fold cross-validation with categorical decision trees + +\ref tree::DecisionTree "DecisionTree" models can be constructed in multiple +other ways. For example, if we have a dataset with both categorical and +numerical features, we can also perform cross-validation by using the associated +\c data::DatasetInfo object. Thus, given some \c data::DatasetInfo object +called \c datasetInfo (that perhaps was produced by a call to \c data::Load() ), +we can perform k-fold cross-validation in a similar manner to the other +examples: + +@code + KFoldCV, Accuracy> cv3(10, data, datasetInfo, labels, + numClasses); + double decisionTreeWithDIAccuracy = cv3.Evaluate(minimumLeafSize); +@endcode + +This particular call to \c cv3.Evaluate() relies on the following +\ref tree::DecisionTree "DecisionTree" constructor: + +@code + template + DecisionTree(MatType&& data, + const data::DatasetInfo& datasetInfo, + LabelsType&& labels, + const size_t numClasses, + const size_t minimumLeafSize = 10); +@endcode + +@subsection cvbasic_ex_4 Simple cross-validation for linear regression + +\c SimpleCV has the same interface as \c KFoldCV, except it takes as one of its +arguments a proportion (from 0 to 1) of data used as a validation set. For +example, to validate \ref regression::LinearRegression "LinearRegression" with +20\% of the data used in the validation set we can write the following code. + +@code + // Random responses for every point from the code snippet in the beginning of + // the tutorial. + arma::rowvec responses = arma::randu(100); + + SimpleCV cv4(0.2, data, responses); + double lrLambda = 0.05; + double lrMSE = cv4.Evaluate(lrLambda); +@endcode + +@section cvbasic_metrics Performance measures + +The cross-validation classes require a performance measure to be specified. +\b mlpack has a number of performance measures implemented; below is a list: + + - mlpack::cv::Accuracy: a simple measure of accuracy + - mlpack::cv::F1: the F1 score; depends on an averaging strategy + - mlpack::cv::MSE: minimum squared error (for regression problems) + - mlpack::cv::Precision: the precision, for classification problems + - mlpack::cv::Recall: the recall, for classification problems + +In addition, it is not difficult to implement a custom performance measure. A +class following the structure below can be used: + +@code +class CustomMeasure +{ + // + // This evaluates the metric given a trained model and a set of data (with + // labels or responses) to evaluate on. The data parameter will be a type of + // Armadillo matrix, and the labels will be the labels that go with the model. + // + // If you know that your model is a classification model (and thus that + // ResponsesType will be arma::Row), it is ok to replace the + // ResponsesType template parameter with arma::Row. + // + template + static double Evaluate(MLAlgorithm& model, + const DataType& data, + const ResponsesType& labels) + { + // Inside the method you should call model.Predict() and compare the + // values with the labels, in order to get the desired performance measure + // and return it. + } +}; +@endcode + +Once this is implemented, then \c CustomMeasure (or whatever the class is +called) is easy to use as a custom performance measure with \c KFoldCV or +\c SimpleCV. + +@section cvbasic_api The KFoldCV and SimpleCV classes + +This section provides details about the \c KFoldCV and \c SimpleCV classes. +The cross-validation infrastructure is based on heavy amounts of template +metaprogramming, so that any \b mlpack learner and any performance measure can +be used. Both classes have two required template parameters and one optional +parameter: + + - \c MLAlgorithm: the type of learner to be used + - \c Metric: the performance measure to be evaluated + - \c MatType: the type of matrix used to store the data + +In addition, there are two more template parameters, but these are automatically +extracted from the given \c MLAlgorithm class, and users should not need to +specify these parameters except when using an unconventional type like +\c arma::fmat for data points. + +The general structure of the \c KFoldCV and \c SimpleCV classes is split into +two parts: + + - The constructor: create the object, and store the data for the \c MLAlgorithm + training. + - The \c Evaluate() method: take any non-data parameters for the + \c MLAlgorithm and calculate the desired performance measure. + +This split is important because it defines the API: all data-related parameters +are passed to the constructor, whereas algorithm hyperparameters are passed to +the \c Evaluate() method. + +@subsection cvbasic_api_constructor The KFoldCV and SimpleCV constructors + +There are six constructors available for \c KFoldCV and \c SimpleCV, each +tailored for a different learning situation. Each is given below for the +\c KFoldCV class, but the same constructors are also available for the +\c SimpleCV class, with the exception that instead of specifying \c k, the +number of folds, the \c SimpleCV class takes a parameter between 0 and 1 +specifying the percentage of the dataset to use as a validation set. + + - `KFoldCV(k, xs, ys)`: this is for unweighted regression applications and + two-class classification applications; \c xs is the dataset and \c ys + are the responses or labels for each point in the dataset. + + - `KFoldCV(k, xs, ys, numClasses)`: this is for unweighted classification + applications; \c xs is the dataset, \c ys are the class labels for each + data point, and \c numClasses is the number of classes in the dataset. + + - `KFoldCV(k, xs, datasetInfo, ys, numClasses)`: this is for unweighted + categorical/numeric classification applications; \c xs is the dataset, + \c datasetInfo is a data::DatasetInfo object that holds the types of + each dimension in the dataset, \c ys are the class labels for each data + point, and \c numClasses is the number of classes in the dataset. + + - `KFoldCV(k, xs, ys, weights)`: this is for weighted regression or + two-class classification applications; \c xs is the dataset, \c ys are + the responses or labels for each point in the dataset, and \c weights + are the weights for each point in the dataset. + + - `KFoldCV(k, xs, ys, numClasses, weights)`: this is for weighted + classification applications; \c xs is the dataset, \c ys are the class + labels for each point in the dataset; \c numClasses is the number of + classes in the dataset, and \c weights holds the weights for each point + in the dataset. + + - `KFoldCV(k, xs, datasetInfo, ys, numClasses, weights)`: this is for + weighted cateogrical/numeric classification applications; \c xs is the + dataset, \c datasetInfo is a data::DatasetInfo object that holds the + types of each dimension in the dataset, \c ys are the class labels for + each data point, \c numClasses is the number of classes in each dataset, + and \c weights holds the weights for each point in the dataset. + +Note that the constructor you should use is the constructor that most closely +matches the constructor of the machine learning algorithm you would like +performance measures of. So, for instance, if you are doing multi-class softmax +regression, you could call the constructor +\c "SoftmaxRegression(xs, ys, numClasses)". Therefore, for \c KFoldCV you would +call the constructor \c "KFoldCV(k, xs, ys, numClasses)" and for \c SimpleCV you +would call the constructor \c "SimpleCV(pct, xs, ys, numClasses)". + +@subsection cvbasic_api_evaluate The Evaluate() method + +The other method that \c KFoldCV and \c SimpleCV have is the method to +actually calculate the performance measure: \c Evaluate(). The \c Evaluate() +method takes any hyperparameters that would follow the data arguments to the +constructor or \c Train() method of the given \c MLAlgorithm. The +\c Evaluate() method takes no more arguments than that, and returns the +desired performance measure on the dataset. + +Therefore, let us suppose that we are interested in cross-validating the +performance of a softmax regression model, and that we have constructed +the appropriate \c KFoldCV object using the code below: + +@code +KFoldCV cv(k, data, labels, numClasses); +@endcode + +The \ref regression::SoftmaxRegression "SoftmaxRegression" class has the +constructor + +@code + template + SoftmaxRegression(const arma::mat& data, + const arma::Row& labels, + const size_t numClasses, + const double lambda = 0.0001, + const bool fitIntercept = false, + OptimizerType optimizer = OptimizerType()); +@endcode + +Note that all parameters after are \c numClasses are optional. This means that +we can specify none or any of them in our call to \c Evaluate(). Below is some +example code showing three different ways we can call \c Evaluate() with the +\c cv object from the code snippet above. + +@code +// First, call with all defaults. +double result1 = cv.Evaluate(); + +// Next, call with lambda set to 0.1 and fitIntercept set to true. +double result2 = cv.Evaluate(0.1, true); + +// Lastly, create a custom optimizer to use for optimization, and use a lambda +// value of 0.5 and fit no intercept. +optimization::SGD<> sgd(0.05, 50000); // Step size of 0.05, 50k max iterations. +double result3 = cv.Evaluate(0.5, false, sgd); +@endcode + +The same general idea applies to any \c MLAlgorithm: all hyperparameters must be +passed to the \c Evaluate() method of \c KFoldCV or \c SimpleCV. + +@section cvbasic_further Further references + +For further documentation, please see the associated Doxygen documentation for +each of the relevant classes: + + - mlpack::cv::SimpleCV + - mlpack::cv::KFoldCV + - mlpack::cv::Accuracy + - mlpack::cv::F1 + - mlpack::cv::MSE + - mlpack::cv::Precision + - mlpack::cv::Recall + +If you are interested in implementing a different cross-validation strategy than +k-fold cross-validation or simple cross-validation, take a look at the +implementations of each of those classes to guide your implementation. + +In addition, the @ref hpt "hyperparameter tuner" documentation may also be +relevant. + +*/ + +} // namespace cv +} // namespace mlpack diff --git a/doc/guide/hpt.hpp b/doc/guide/hpt.hpp new file mode 100644 index 0000000000..ff6bc8ffa1 --- /dev/null +++ b/doc/guide/hpt.hpp @@ -0,0 +1,238 @@ +namespace mlpack { +namespace hpt { + +/*! @page hpt Hyper-Parameter Tuning + +@section hptintro Introduction + +\b mlpack implements a generic hyperparameter tuner that is able to tune both +continuous and discrete parameters of various different algorithms. This is an +important task---the performance of many machine learning algorithms can be +highly dependent on the hyperparameters that are chosen for that algorithm. +(One example: the choice of \f$k\f$ for a \f$k\f$-nearest-neighbors classifier.) + +This hyper-parameter tuner is built on the same general concept as the +cross-validation classes (see the @ref cv "cross-validation tutorial"): given +some machine learning algorithm, some data, some performance measure, and a set +of hyperparameters, attempt to find the hyperparameter set that best optimizes +the performance measure on the given data with the given algorithm. + +\b mlpack's implementation of hyperparameter tuning is flexible, and is built in +a way that supports many algorithms and many optimizers. At the time of this +writing, complex hyperparameter optimization techniques are not available, but +the hyperparameter tuner does support these, should they be implemented in the +future. + +In this tutorial we will see the usage examples of the hyper-parameter tuning +module, and also more details about the \c HyperParameterTuner class. + +@section hptbasic Basic Usage + +The interface of the hyper-parameter tuning module is quite similar to the +interface of the @ref cv "cross-validation module". To construct a \c +HyperParameterTuner object you need to specify as template parameters what +machine learning algorithm, cross-validation strategy, performance measure, and +optimization strategy (\ref optimization::GridSearch "GridSearch" will be used by +default) you are going to use. Then, you must pass the same arguments as for +the cross-validation classes: the data and labels (or responses) to use are +given to the constructor, and the possible hyperparameter values are given to +the \c HyperParameterTuner::Optimize() method, which returns the best +algorithm configuration as a \c std::tuple<>. + +Let's see some examples. + +Suppose we have the following data to train and validate on. +@code + // 100-point 5-dimensional random dataset. + arma::mat data = arma::randu(5, 100); + // Noisy responses retrieved by a random linear transformation of data. + arma::rowvec responses = arma::randu(5) * data + + 0.1 * arma::randn(100); +@endcode + +Given the dataset above, we can use the following code to try to find a good \c +lambda value for \ref regression::LinearRegression "LinearRegression". Here we +use \ref cv::SimpleCV "SimpleCV" instead of k-fold cross-validation to save +computation time. + +@code + // Using 80% of data for training and remaining 20% for assessing MSE. + double validationSize = 0.2; + HyperParameterTuner hpt(validationSize, + data, responses); + + // Finding a good value for lambda from the discrete set of values 0.0, 0.001, + // 0.01, 0.1, and 1.0. + arma::vec lambdas{0.0, 0.001, 0.01, 0.1, 1.0}; + double bestLambda; + std::tie(bestLambda) = hpt.Optimize(lambdas); +@endcode + +In this example we have used \ref optimization::GridSearch "GridSearch" (the +default optimizer) to find a good value for the \c lambda hyper-parameter. For +that we have specified what values should be tried. + +@section hptfixed Fixed Arguments + +When some hyper-parameters should not be optimized, you can specify values +for them with the \c Fixed() method as in the following example of trying to +find good \c lambda1 and \c lambda2 values for \ref regression::LARS "LARS" +(least-angle regression). + +@code + HyperParameterTuner hpt2(validationSize, data, + responses); + + // The hyper-parameter tuner should not try to change the transposeData or + // useCholesky parameters. + bool transposeData = true; + bool useCholesky = false; + + // We wish only to search for the best lambda1 and lambda2 values. + arma::vec lambda1Set{0.0, 0.001, 0.01, 0.1, 1.0}; + arma::vec lambda2Set{0.0, 0.002, 0.02, 0.2, 2.0}; + + double bestLambda1, bestLambda2; + std::tie(bestLambda1, bestLambda2) = hpt2.Optimize(Fixed(transposeData), + Fixed(useCholesky), lambda1Set, lambda2Set); +@endcode + +Note that for the call to \c hpt2.Optimize(), we have used the same order of +arguments as they appear in the corresponding \ref regression::LARS "LARS" +constructor: + +@code + LARS(const arma::mat& data, + const arma::rowvec& responses, + const bool transposeData = true, + const bool useCholesky = false, + const double lambda1 = 0.0, + const double lambda2 = 0.0, + const double tolerance = 1e-16); +@endcode + +@section hptgradient Gradient-Based Optimization + +In some cases we may wish to optimize a hyperparameter over the space of all +possible real values, instead of providing a grid in which to search. +Alternately, we may know approximately optimal values from a grid search for +real-valued hyperparameters, but wish to further tune those values. + +In this case, we can use a gradient-based optimizer for hyperparameter search. +In the following example, we try to optimize the \c lambda1 and \c lambda2 +hyper-parameters for \ref regression::LARS "LARS" with the +\ref optimization::GradientDescent "GradientDescent" optimizer. + +@code + HyperParameterTuner hpt3(validationSize, + data, responses); + + // GradientDescent can be adjusted in the following way. + hpt3.Optimizer().StepSize() = 0.1; + hpt3.Optimizer().Tolerance() = 1e-15; + + // We can set up values used for calculating gradients. + hpt3.RelativeDelta() = 0.01; + hpt3.MinDelta() = 1e-10; + + double initialLambda1 = 0.001; + double initialLambda2 = 0.002; + + double bestGDLambda1, bestGDLambda2; + std::tie(bestGDLambda1, bestGDLambda2) = hpt3.Optimize(Fixed(transposeData), + Fixed(useCholesky), initialLambda1, initialLambda2); +@endcode + +@section hpt_class The HyperParameterTuner class + +The \c HyperParameterTuner class is very similar to the +\ref cv::KFoldCV "KFoldCV" and \ref cv::SimpleCV "SimpleCV" classes (see the +@ref "cross-validation tutorial" for more information on those two classes), but +there are a few important differences. + +First, the \c HyperParameterTuner accepts five different hyperparameters; only +the first three of these are required: + + - \c MLAlgorithm This is the algorithm to be used. + - \c Metric This is the performance measure to be used; see + @ref cvbasic_metrics for more information. + - \c CVType This is the type of cross-validation to be used for evaluating the + performance measure; this should be \ref cv::KFoldCV "KFoldCV" or + \ref cv::SimpleCV "SimpleCV". + - \c OptimizerType This is the type of optimizer to use; it can be + \c GridSearch or a gradient-based optimizer. + - \c MatType This is the type of data matrix to use. The default is + \c arma::mat. This only needs to be changed if you are specifically + using sparse data, or if you want to use a numeric type other than + \c double. + +The last two template parameters are automatically inferred by the +\c HyperParameterTuner and should not need to be manually specified, unless an +unconventional data type like \c arma::fmat is being used for data points. + +Typically, \ref cv::SimpleCV "SimpleCV" is a good choice for \c CVType because +it takes so much less time to compute than full \ref cv::KFoldCV "KFoldCV"; +however, the disadvantage is that \ref cv::SimpleCV "SimpleCV" might give a +somewhat more noisy estimate of the performance measure on unseen test data. + +The constructor for the \c HyperParameterTuner is called with exactly the same +arguments as the corresponding \c CVType that has been chosen. For more +information on that, please see the +@ref cvbasic_api "cross-validation constructor tutorial". As an example, if we +are using \ref cv::SimpleCV "SimpleCV" and wish to hold out 20\% of the dataset +as a validation set, we might construct a \c HyperParameterTuner like this: + +@code +// We will use LinearRegression as the MLAlgorithm, and MSE as the performance +// measure. Our dataset is 'dataset' and the responses are 'responses'. +HyperParameterTuner hpt(0.2, dataset, + responses); +@endcode + +Next, we must set up the hyperparameters to be optimized. If we are doing a +grid search with the \ref optimization::GridSearch "GridSearch" optimizer (the +default), then we only need to pass a `std::vector` (for non-numeric +hyperparameters) or an `arma::vec` (for numeric hyperparameters) containing all +of the possible choices that we wish to search over. + +For instance, a set of numeric values might be chosen like this, for the +\c lambda parameter (of type \c double): + +@code +arma::vec lambdaSet = arma::vec("0.0 0.1 0.5 1.0"); +@endcode + +Similarly, a set of non-numeric values might be chosen like this, for the +\c intercept parameter: + +@code +std::vector interceptSet = { false, true }; +@endcode + +Once all of these are set up, the \c HyperParameterTuner::Optimize() method may +be called to find the best set of hyperparameters: + +@code +bool intercept; +double lambda; +std::tie(lambda, intercept) = hpt.Optimize(lambdaSet, interceptSet); +@endcode + +Alternately, the \c Fixed() method (detailed in the @ref hptfixed +"Fixed arguments" section) can be used to fix the values of some parameters. + +For continuous optimizers like +\ref optimization::GradientDescent "GradientDescent", a range does not need to +be specified but instead only a single value. See the +\ref hptgradient "Gradient-Based Optimization" section for more details. + +@section hptfurther Further documentation + +For more information on the \c HyperParameterTuner class, see the +mlpack::hpt::HyperParameterTuner class documentation and the +@ref cv "cross-validation tutorial". + +*/ + +} // namespace hpt +} // namespace mlpack diff --git a/doc/tutorials/tutorials.txt b/doc/tutorials/tutorials.txt index 204bfd1477..e364da6bd9 100644 --- a/doc/tutorials/tutorials.txt +++ b/doc/tutorials/tutorials.txt @@ -18,6 +18,8 @@ start. - \ref iodoc - \ref timer - \ref sample + - \ref cv + - \ref hpt @section method_tut Method-specific Tutorials diff --git a/src/mlpack/core.hpp b/src/mlpack/core.hpp index 30c65cf54d..03798eeba6 100644 --- a/src/mlpack/core.hpp +++ b/src/mlpack/core.hpp @@ -100,6 +100,8 @@ * - @ref iodoc * - @ref timer * - @ref sample + * - @ref cv + * - @ref hpt * - @ref verinfo * * Tutorials on specific methods are also available. diff --git a/src/mlpack/core/tree/CMakeLists.txt b/src/mlpack/core/tree/CMakeLists.txt index adc5bae0b2..0b2f01a5f3 100644 --- a/src/mlpack/core/tree/CMakeLists.txt +++ b/src/mlpack/core/tree/CMakeLists.txt @@ -120,6 +120,7 @@ set(SOURCES statistic.hpp traversal_info.hpp tree_traits.hpp + enumerate_tree.hpp ) # add directory name to sources diff --git a/src/mlpack/core/tree/cosine_tree/cosine_tree.cpp b/src/mlpack/core/tree/cosine_tree/cosine_tree.cpp index 929ca82b4e..61159a6001 100644 --- a/src/mlpack/core/tree/cosine_tree/cosine_tree.cpp +++ b/src/mlpack/core/tree/cosine_tree/cosine_tree.cpp @@ -94,7 +94,8 @@ CosineTree::CosineTree(const arma::mat& dataset, // Initialize Monte Carlo error estimate for comparison. double monteCarloError = root.FrobNormSquared(); - while (monteCarloError > epsilon * root.FrobNormSquared()) + while (treeQueue.top() && + (monteCarloError > epsilon * root.FrobNormSquared())) { // Pop node from queue with highest projection error. CosineTree* currentNode; diff --git a/src/mlpack/core/tree/enumerate_tree.hpp b/src/mlpack/core/tree/enumerate_tree.hpp new file mode 100644 index 0000000000..8811c8f4c2 --- /dev/null +++ b/src/mlpack/core/tree/enumerate_tree.hpp @@ -0,0 +1,64 @@ +/** + * @file enumerate_tree.hpp + * @author Ivan (Jonan) Georgiev + * + * This file contains function that performs a simple depth-first walk on the tree + * calling `Enter` and `Leave` methods of a provided walker. + * + * 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. + */ +#ifndef MLPACK_CORE_TREE_ENUMERATE_TREE_HPP +#define MLPACK_CORE_TREE_ENUMERATE_TREE_HPP + +namespace mlpack { +namespace tree /** Trees and tree-building procedures. */ { +namespace enumerate { + +// Actual implementation of the enumeration. The problem is the unified +// detection if we're on the root, because Enter and Leave expect the +// parent being passed. +template +void EnumerateTreeImpl(TreeType* tree, Walker& walker, bool root) +{ + if (root) + walker.Enter(tree, (const TreeType*)nullptr); + + const size_t numChildren = tree->NumChildren(); + for (size_t i = 0; i < numChildren; ++i) + { + TreeType* child = tree->ChildPtr(i); + walker.Enter(child, tree); + EnumerateTreeImpl(child, walker, false); + walker.Leave(child, tree); + } + + if (root) + walker.Leave(tree, (const TreeType*)nullptr); +} + +} // namespace enumerate + + +/** + * Traverses all nodes of the tree, including the inner ones. On each node + * two methods of the `enumer` are called: + * + * Enter(TreeType* node, TreeType* parent); + * Leave(TreeType* node, TreeType* parent); + * + * @param walker An instance of custom class, receiver of the enumeration. + */ +template +inline void EnumerateTree(TreeType* tree, Walker& walker) +{ + enumerate::EnumerateTreeImpl(tree, walker, true); +} + +} // namespace tree +} // namespace mlpack + + +#endif // MLPACK_CORE_TREE_ENUMERATE_TREE_HPP diff --git a/src/mlpack/core/util/timers.cpp b/src/mlpack/core/util/timers.cpp index 974b3eba40..d05312a791 100644 --- a/src/mlpack/core/util/timers.cpp +++ b/src/mlpack/core/util/timers.cpp @@ -111,7 +111,6 @@ void Timers::PrintTimer(const std::string& timerName) Log::Info << ", "; Log::Info << s.count() << "." << std::setw(1) << (totalDurationMicroSec.count() / 100000) << " secs"; - output = true; } Log::Info << ")"; diff --git a/src/mlpack/methods/det/det_main.cpp b/src/mlpack/methods/det/det_main.cpp index e9755fe415..41568dfbaa 100644 --- a/src/mlpack/methods/det/det_main.cpp +++ b/src/mlpack/methods/det/det_main.cpp @@ -33,6 +33,15 @@ PROGRAM_INFO("Density Estimation With Density Estimation Trees", " with the " + PRINT_PARAM_STRING("training_set_estimates") + " output " "parameter." "\n\n" + "Enabling path printing for each node outputs the path from the root node " + "to a leaf for each entry in the test set, or training set (if a test set " + "is not provided). Strings like 'LRLRLR' (indicating that traversal went " + "to the left child, then the right child, then the left child, and so " + "forth) will be output. If 'lr-id' or 'id-lr' are given as the " + + PRINT_PARAM_STRING("path_format") + " parameter, then the ID (tag) of " + "every node along the path will be printed after or before the L or R " + "character indicating the direction of traversal, respectively." + "\n\n" "This program also can provide density estimates for a set of test points, " "specified in the " + PRINT_PARAM_STRING("test") + " parameter. The " "density estimation tree used for this task will be the tree that was " @@ -61,6 +70,19 @@ PARAM_MATRIX_OUT("test_set_estimates", "The output estimates on the test set " PARAM_MATRIX_OUT("vi", "The output variable importance values for each " "feature.", "i"); +// Tagging and path printing options +PARAM_STRING_IN("path_format", "The format of path printing: 'lr', 'id-lr', or " + "'lr-id'.", "p", "lr"); + +PARAM_STRING_OUT("tag_counters_file", "The file to output the number of points " + "that went to each leaf.", "c"); + +PARAM_STRING_OUT("tag_file", "The file to output the tags (and possibly paths)" + " for each sample in the test set.", "g"); + +PARAM_FLAG("skip_pruning", "Whether to bypass the pruning process and output " + "the unpruned tree only.", "s"); + // Parameters for the training algorithm. PARAM_INT_IN("folds", "The number of folds of cross-validation to perform for " "the estimation (0 is LOOCV)", "f", 10); @@ -76,71 +98,80 @@ PARAM_FLAG("volume_regularization", "This flag gives the used the option to use" "penalize low volume leaves.", "R"); */ + void mlpackMain() { // Validate input parameters. if (CLI::HasParam("training") && CLI::HasParam("input_model")) - Log::Fatal << "Only one of --training_file (-t) or --input_model_file (-m) " - << "may be specified!" << endl; + Log::Fatal << "Only one of " << PRINT_PARAM_STRING("training") << " or " << + PRINT_PARAM_STRING("input_model") << " may be specified!" << endl; if (!CLI::HasParam("training") && !CLI::HasParam("input_model")) - Log::Fatal << "Neither --training_file (-t) nor --input_model_file (-m) " - << "are specified!" << endl; + Log::Fatal << "Neither " << PRINT_PARAM_STRING("training") << " nor " << + PRINT_PARAM_STRING("input_model") << " are specified!" << endl; + + if (CLI::HasParam("tag_file") && + !CLI::HasParam("training") && !CLI::HasParam("test")) + { + Log::Fatal << "Neither " << PRINT_PARAM_STRING("training") << " nor " << + PRINT_PARAM_STRING("test") << " are specified, but needed when " << + PRINT_PARAM_STRING("tag_file") << " is asked." << endl; + } if (!CLI::HasParam("training")) { if (CLI::HasParam("training_set_estimates")) - Log::Warn << "--training_set_estimates_file (-e) ignored because " - << "--training_file (-t) is not specified." << endl; + Log::Warn << PRINT_PARAM_STRING("training_set_estimates") << + " ignored because " << PRINT_PARAM_STRING("training") << + " is not specified." << endl; if (CLI::HasParam("folds")) - Log::Warn << "--folds (-f) ignored because --training_file (-t) is not " - << "specified." << endl; + Log::Warn << PRINT_PARAM_STRING("folds") << " ignored because " << + PRINT_PARAM_STRING("training") << " is not specified." << endl; if (CLI::HasParam("min_leaf_size")) - Log::Warn << "--min_leaf_size (-l) ignored because --training_file (-t) " - << "is not specified." << endl; + Log::Warn << PRINT_PARAM_STRING("min_leaf_size") << " ignored because " << + PRINT_PARAM_STRING("training") << " is not specified." << endl; if (CLI::HasParam("max_leaf_size")) - Log::Warn << "--max_leaf_size (-L) ignored because --training_file (-t) " - << "is not specified." << endl; + Log::Warn << PRINT_PARAM_STRING("max_leaf_size") << " ignored because " << + PRINT_PARAM_STRING("training") << " is not specified." << endl; } else if (!CLI::HasParam("output_model") && !CLI::HasParam("training_set_estimates") && !CLI::HasParam("vi")) { - Log::Warn << "None of --output_model_file (-M), --training_set_estimates " - << "(-e), or --vi (-i) are specified; no output will be saved!" << endl; + Log::Warn << "None of " << PRINT_PARAM_STRING("output_model") << ", " << + PRINT_PARAM_STRING("training_set_estimates") << ", or " << + PRINT_PARAM_STRING("vi") << " are specified; no output will be saved!" << + endl; } if (!CLI::HasParam("test") && CLI::HasParam("test_set_estimates")) - Log::Warn << "--test_set_estimates_file (-E) ignored because --test_file " - << "(-T) is not specified." << endl; + Log::Warn << PRINT_PARAM_STRING("test_set_estimates") << " ignored " << + "because " << PRINT_PARAM_STRING("test") << " is not specified." << endl; // Are we training a DET or loading from file? DTree* tree; + arma::mat trainingData; + arma::mat testData; + if (CLI::HasParam("training")) { - arma::mat trainingData = std::move(CLI::GetParam("training")); - - // Cross-validation here. - size_t folds = CLI::GetParam("folds"); - if (folds == 0) - { - folds = trainingData.n_cols; - Log::Info << "Performing leave-one-out cross validation." << endl; - } - else - { - Log::Info << "Performing " << folds << "-fold cross validation." << endl; - } + trainingData = std::move(CLI::GetParam("training")); const bool regularization = false; // const bool regularization = CLI::HasParam("volume_regularization"); const int maxLeafSize = CLI::GetParam("max_leaf_size"); const int minLeafSize = CLI::GetParam("min_leaf_size"); + const bool skipPruning = CLI::HasParam("skip_pruning"); + size_t folds = CLI::GetParam("folds"); + + if (folds == 0) + folds = trainingData.n_cols; // Obtain the optimal tree. Timer::Start("det_training"); tree = Trainer(trainingData, folds, regularization, - maxLeafSize, minLeafSize, ""); + maxLeafSize, minLeafSize, + skipPruning); Timer::Stop("det_training"); // Compute training set estimates, if desired. @@ -166,25 +197,112 @@ void mlpackMain() // the given file. if (CLI::HasParam("test")) { - arma::mat testData = std::move(CLI::GetParam("test")); - - // Compute test set densities. - Timer::Start("det_test_set_estimation"); - arma::rowvec testDensities(testData.n_cols); - for (size_t i = 0; i < testData.n_cols; i++) - testDensities[i] = tree->ComputeValue(testData.unsafe_col(i)); - Timer::Stop("det_test_set_estimation"); - + testData = std::move(CLI::GetParam("test")); if (CLI::HasParam("test_set_estimates")) + { + // Compute test set densities. + Timer::Start("det_test_set_estimation"); + arma::rowvec testDensities(testData.n_cols); + + for (size_t i = 0; i < testData.n_cols; i++) + testDensities[i] = tree->ComputeValue(testData.unsafe_col(i)); + + Timer::Stop("det_test_set_estimation"); + CLI::GetParam("test_set_estimates") = std::move(testDensities); + } + + // Print variable importance. + if (CLI::HasParam("vi")) + { + arma::vec importances; + tree->ComputeVariableImportance(importances); + CLI::GetParam("vi") = importances.t(); + } } - // Print variable importance. - if (CLI::HasParam("vi")) + if (CLI::HasParam("tag_file")) { - arma::vec importances; - tree->ComputeVariableImportance(importances); - CLI::GetParam("vi") = std::move(importances.t()); + const arma::mat& estimationData = + CLI::HasParam("test") ? testData : trainingData; + const string tagFile = CLI::GetParam("tag_file"); + std::ofstream ofs; + ofs.open(tagFile, std::ofstream::out); + + arma::Row counters; + + Timer::Start("det_test_set_tagging"); + if (!ofs.is_open()) + { + Log::Warn << "Unable to open file '" << tagFile + << "' to save tag membership info." + << std::endl; + } + else if (CLI::HasParam("path_format")) + { + const bool reqCounters = CLI::HasParam("tag_counters_file"); + const string pathFormat = CLI::GetParam("path_format"); + + PathCacher::PathFormat theFormat; + if (pathFormat == "lr" || pathFormat == "LR") + theFormat = PathCacher::FormatLR; + else if (pathFormat == "lr-id" || pathFormat == "LR-ID") + theFormat = PathCacher::FormatLR_ID; + else if (pathFormat == "id-lr" || pathFormat == "ID-LR") + theFormat = PathCacher::FormatID_LR; + else + { + Log::Warn << "Unknown path format specified: '" << pathFormat + << "'. Valid are: lr | lr-id | id-lr. Defaults to 'lr'." << endl; + theFormat = PathCacher::FormatLR; + } + + PathCacher path(theFormat, tree); + counters.zeros(path.NumNodes()); + + for (size_t i = 0; i < estimationData.n_cols; i++) + { + int tag = tree->FindBucket(estimationData.unsafe_col(i)); + + ofs << tag << " " << path.PathFor(tag) << std::endl; + for (; tag >= 0 && reqCounters; tag = path.ParentOf(tag)) + counters(tag) += 1; + } + + ofs.close(); + + if (reqCounters) + { + ofs.open(CLI::GetParam("tag_counters_file"), + std::ofstream::out); + + for (size_t j = 0; j < counters.n_elem; ++j) + ofs << j << " " + << counters(j) << " " + << path.PathFor(j) << endl; + + ofs.close(); + } + } + else + { + int numLeaves = tree->TagTree(); + counters.zeros(numLeaves); + + for (size_t i = 0; i < estimationData.n_cols; i++) + { + const int tag = tree->FindBucket(estimationData.unsafe_col(i)); + + ofs << tag << std::endl; + counters(tag) += 1; + } + + if (CLI::HasParam("tag_counters_file")) + data::Save(CLI::GetParam("tag_counters_file"), counters); + } + + Timer::Stop("det_test_set_tagging"); + ofs.close(); } // Save the model, if desired. diff --git a/src/mlpack/methods/det/dt_utils.hpp b/src/mlpack/methods/det/dt_utils.hpp index fd66e2e6ce..3e2c9ec7d4 100644 --- a/src/mlpack/methods/det/dt_utils.hpp +++ b/src/mlpack/methods/det/dt_utils.hpp @@ -35,7 +35,7 @@ void PrintLeafMembership(DTree* dtree, const MatType& data, const arma::Mat& labels, const size_t numClasses, - const std::string leafClassMembershipFile = ""); + const std::string& leafClassMembershipFile = ""); /** * Print the variable importance of each dimension of a density estimation tree. @@ -67,7 +67,79 @@ DTree* Trainer(MatType& dataset, const bool useVolumeReg = false, const size_t maxLeafSize = 10, const size_t minLeafSize = 5, - const std::string unprunedTreeOutput = ""); + const std::string unprunedTreeOutput = "", + const bool skipPruning = false); + +/** + * This class is responsible for caching the path to each node of the tree. Its + * instance is provided to EnumerateTree() utility ONCE and it caches the paths + * to all the leafs and then easily (and quickly) retrieves these paths for each + * test entry. + */ +class PathCacher +{ + public: + /** + * Possible formats to use for output. + */ + enum PathFormat + { + //! Print only whether we went left or right. + FormatLR, + //! Print the direction, then the tag of the node. + FormatLR_ID, + //! Print the tag of the node, then the direction. + FormatID_LR + }; + + /** + * Construct a PathCacher object on the given tree with the given format. + * + * @param fmt Format to use for output. + * @param tree Tree to cache paths in. + */ + template + PathCacher(PathFormat fmt, DTree* tree); + + /** + * Enter a given node. + */ + template + void Enter(const DTree* node, + const DTree* parent); + + /** + * Leave the given node. + */ + template + void Leave(const DTree* node, + const DTree* parent); + + /** + * Return the constructed path for a given tag. + */ + const std::string& PathFor(int tag) const; + + /** + * Get the parent tag of a given tag. + */ + int ParentOf(int tag) const; + + /** + * Get the number of nodes in the path cache. + */ + size_t NumNodes() const { return pathCache.size(); } + + protected: + typedef std::list> PathType; + typedef std::vector> PathCacheType; + + PathType path; + PathFormat format; + PathCacheType pathCache; + + std::string BuildString(); +}; } // namespace det } // namespace mlpack diff --git a/src/mlpack/methods/det/dt_utils_impl.hpp b/src/mlpack/methods/det/dt_utils_impl.hpp index d4dba48467..866413245b 100644 --- a/src/mlpack/methods/det/dt_utils_impl.hpp +++ b/src/mlpack/methods/det/dt_utils_impl.hpp @@ -14,19 +14,20 @@ #define MLPACK_METHODS_DET_DT_UTILS_IMPL_HPP #include "dt_utils.hpp" +#include namespace mlpack { namespace det { -template -void PrintLeafMembership(DTree* dtree, +template +void PrintLeafMembership(DTree* dtree, const MatType& data, const arma::Mat& labels, const size_t numClasses, - const std::string leafClassMembershipFile) + const std::string& leafClassMembershipFile) { // Tag the leaves with numbers. - TagType numLeaves = dtree->TagTree(); + int numLeaves = dtree->TagTree(); arma::Mat table(numLeaves, (numClasses + 1)); table.zeros(); @@ -34,7 +35,7 @@ void PrintLeafMembership(DTree* dtree, for (size_t i = 0; i < data.n_cols; i++) { const typename MatType::vec_type testPoint = data.unsafe_col(i); - const TagType leafTag = dtree->FindBucket(testPoint); + const int leafTag = dtree->FindBucket(testPoint); const size_t label = labels[i]; table(leafTag, label) += 1; } @@ -111,11 +112,12 @@ DTree* Trainer(MatType& dataset, const bool useVolumeReg, const size_t maxLeafSize, const size_t minLeafSize, - const std::string unprunedTreeOutput) + const bool skipPruning) { // Initialize the tree. - DTree dtree(dataset); + DTree* dtree = new DTree(dataset); + Timer::Start("tree_growing"); // Prepare to grow the tree... arma::Col oldFromNew(dataset.n_cols); for (size_t i = 0; i < oldFromNew.n_elem; i++) @@ -126,60 +128,52 @@ DTree* Trainer(MatType& dataset, // Growing the tree double oldAlpha = 0.0; - double alpha = dtree.Grow(newDataset, oldFromNew, useVolumeReg, maxLeafSize, + double alpha = dtree->Grow(newDataset, oldFromNew, useVolumeReg, maxLeafSize, minLeafSize); - Log::Info << dtree.SubtreeLeaves() << " leaf nodes in the tree using full " + Timer::Stop("tree_growing"); + Log::Info << dtree->SubtreeLeaves() << " leaf nodes in the tree using full " << "dataset; minimum alpha: " << alpha << "." << std::endl; - // Compute densities for the training points in the full tree, if we were - // asked for this. - if (unprunedTreeOutput != "") - { - std::ofstream outfile(unprunedTreeOutput.c_str()); - if (outfile.good()) - { - for (size_t i = 0; i < dataset.n_cols; ++i) - { - arma::vec testPoint = dataset.unsafe_col(i); - outfile << dtree.ComputeValue(testPoint) << std::endl; - } - } - else - { - Log::Warn << "Can't open '" << unprunedTreeOutput << "' to write computed" - << " densities to." << std::endl; - } + if (skipPruning) + return dtree; - outfile.close(); - } + if (folds == dataset.n_cols) + Log::Info << "Performing leave-one-out cross validation." << std::endl; + else + Log::Info << "Performing " << folds << "-fold cross validation." << + std::endl; + + Timer::Start("pruning_sequence"); // Sequentially prune and save the alpha values and the values of c_t^2 * r_t. std::vector > prunedSequence; - while (dtree.SubtreeLeaves() > 1) + while (dtree->SubtreeLeaves() > 1) { std::pair treeSeq(oldAlpha, - dtree.SubtreeLeavesLogNegError()); + dtree->SubtreeLeavesLogNegError()); prunedSequence.push_back(treeSeq); oldAlpha = alpha; - alpha = dtree.PruneAndUpdate(oldAlpha, dataset.n_cols, useVolumeReg); + alpha = dtree->PruneAndUpdate(oldAlpha, dataset.n_cols, useVolumeReg); // Some sanity checks. It seems that on some datasets, the error does not // increase as the tree is pruned but instead stays the same---hence the // "<=" in the final assert. Log::Assert((alpha < std::numeric_limits::max()) - || (dtree.SubtreeLeaves() == 1)); + || (dtree->SubtreeLeaves() == 1)); Log::Assert(alpha > oldAlpha); - Log::Assert(dtree.SubtreeLeavesLogNegError() <= treeSeq.second); + Log::Assert(dtree->SubtreeLeavesLogNegError() <= treeSeq.second); } - std::pair treeSeq(oldAlpha, dtree.SubtreeLeavesLogNegError()); + std::pair treeSeq(oldAlpha, + dtree->SubtreeLeavesLogNegError()); prunedSequence.push_back(treeSeq); + Timer::Stop("pruning_sequence"); Log::Info << prunedSequence.size() << " trees in the sequence; maximum alpha:" << " " << oldAlpha << "." << std::endl; - MatType cvData(dataset); + const MatType cvData(dataset); const size_t testSize = dataset.n_cols / folds; arma::vec regularizationConstants(prunedSequence.size()); @@ -191,7 +185,7 @@ DTree* Trainer(MatType& dataset, // implementation. omp_size_t is the appropriate type according to the // platform. #pragma omp parallel for default(none) \ - shared(cvData, prunedSequence, regularizationConstants) + shared(prunedSequence, regularizationConstants) for (omp_size_t fold = 0; fold < (omp_size_t) folds; fold++) { // Break up data into train and test sets. @@ -289,8 +283,9 @@ DTree* Trainer(MatType& dataset, Log::Info << "Optimal alpha: " << optimalAlpha << "." << std::endl; - // Initialize the tree. - DTree* dtreeOpt = new DTree(dataset); + // Re-Initialize the tree. + delete dtree; + dtree = new DTree(dataset); // Getting ready to grow the tree... for (size_t i = 0; i < oldFromNew.n_elem; i++) @@ -301,31 +296,97 @@ DTree* Trainer(MatType& dataset, // Grow the tree. oldAlpha = -DBL_MAX; - alpha = dtreeOpt->Grow(newDataset, + alpha = dtree->Grow(newDataset, oldFromNew, useVolumeReg, maxLeafSize, minLeafSize); // Prune with optimal alpha. - while ((oldAlpha < optimalAlpha) && (dtreeOpt->SubtreeLeaves() > 1)) + while ((oldAlpha < optimalAlpha) && (dtree->SubtreeLeaves() > 1)) { oldAlpha = alpha; - alpha = dtreeOpt->PruneAndUpdate(oldAlpha, newDataset.n_cols, useVolumeReg); + alpha = dtree->PruneAndUpdate(oldAlpha, newDataset.n_cols, useVolumeReg); // Some sanity checks. Log::Assert((alpha < std::numeric_limits::max()) || - (dtreeOpt->SubtreeLeaves() == 1)); + (dtree->SubtreeLeaves() == 1)); Log::Assert(alpha > oldAlpha); } - Log::Info << dtreeOpt->SubtreeLeaves() << " leaf nodes in the optimally " + Log::Info << dtree->SubtreeLeaves() << " leaf nodes in the optimally " << "pruned tree; optimal alpha: " << oldAlpha << "." << std::endl; - return dtreeOpt; + return dtree; +} + +template +PathCacher::PathCacher(PathCacher::PathFormat fmt, DTree* dtree) : + format(fmt) +{ + // Here we use TagTree()'s output to determine the + // number of _nodes_ in the tree. + pathCache.resize(dtree->TagTree(0, true)); + pathCache[0] = PathCacheType::value_type(-1, ""); + tree::EnumerateTree(dtree, *this); +} + +template +void PathCacher::Enter(const DTree* node, + const DTree* parent) +{ + if (parent == nullptr) + return; + + int tag = node->BucketTag(); + + path.push_back(PathType::value_type(parent->Left() == node, tag)); + pathCache[tag] = PathCacheType::value_type(parent->BucketTag(), + (node->SubtreeLeaves() > 1) ? + "" : BuildString()); +} + +template +void PathCacher::Leave(const DTree* /* node */, + const DTree* parent) +{ + if (parent != nullptr) + path.pop_back(); +} + +std::string PathCacher::BuildString() +{ + std::string str(""); + for (PathType::iterator it = path.begin(); it != path.end(); it++) + { + switch (format) + { + case FormatLR: + str += it->first ? "L" : "R"; + break; + case FormatLR_ID: + str += (it->first ? "L" : "R") + std::to_string(it->second); + break; + case FormatID_LR: + str += std::to_string(it->second) + (it->first ? "L" : "R"); + break; + } + } + + return str; +} + +int PathCacher::ParentOf(int tag) const +{ + return pathCache[tag].first; +} + +const std::string& PathCacher::PathFor(int tag) const +{ + return pathCache[tag].second; } } // namespace det } // namespace mlpack -#endif +#endif // MLPACK_METHODS_DET_DT_UTILS_IMPL_HPP diff --git a/src/mlpack/methods/det/dtree.hpp b/src/mlpack/methods/det/dtree.hpp index 04ce0b0743..88acb45990 100644 --- a/src/mlpack/methods/det/dtree.hpp +++ b/src/mlpack/methods/det/dtree.hpp @@ -46,12 +46,12 @@ template StatType; + //! The actual, underlying type we're working with. + typedef typename MatType::elem_type ElemType; + //! The type of vector we are using. + typedef typename MatType::vec_type VecType; + //! The statistic type we are holding. + typedef typename arma::Col StatType; /** * Create an empty density estimation tree. @@ -184,11 +184,14 @@ class DTree /** * Index the buckets for possible usage later; this results in every leaf in * the tree having a specific tag (accessible with BucketTag()). This - * function calls itself recursively. + * function calls itself recursively. The tag is incremented with + * `operator++()`, so any `TagType` overriding it will do. * * @param tag Tag for the next leaf; leave at 0 for the initial call. + * @param everyNodde Whether to increment on every node, not just leaves. */ - TagType TagTree(const TagType& tag = 0); + TagType TagTree(const TagType& tag = 0, bool everyNode = false); + /** * Return the tag of the leaf containing the query. This is useful for @@ -198,6 +201,7 @@ class DTree */ TagType FindBucket(const VecType& query) const; + /** * Compute the variable importance of each dimension in the learned tree. * @@ -301,7 +305,19 @@ class DTree //! Return the upper part of the alpha sum. double AlphaUpper() const { return alphaUpper; } //! Return the current bucket's ID, if leaf, or -1 otherwise - TagType BucketTag() const { return subtreeLeaves == 1 ? bucketTag : -1; } + TagType BucketTag() const { return bucketTag; } + //! Return the number of children in this node. + size_t NumChildren() const { return !left ? 0 : 2; } + + /** + * Return the specified child (0 will be left, 1 will be right). If the index + * is greater than 1, this will return the right child. + * + * @param child Index of child to return. + */ + DTree& Child(const size_t child) const { return !child ? *left : *right; } + + DTree*& ChildPtr(const size_t child) { return (!child) ? left : right; } //! Return the maximum values. const StatType& MaxVals() const { return maxVals; } @@ -335,6 +351,9 @@ class DTree const size_t splitDim, const ElemType splitValue, arma::Col& oldFromNew) const; + + void FillMinMax(const StatType& mins, + const StatType& maxs); }; } // namespace det diff --git a/src/mlpack/methods/det/dtree_impl.hpp b/src/mlpack/methods/det/dtree_impl.hpp index 6372b19c0f..5dbdad7254 100644 --- a/src/mlpack/methods/det/dtree_impl.hpp +++ b/src/mlpack/methods/det/dtree_impl.hpp @@ -20,134 +20,136 @@ using namespace det; namespace details { - /** - * This one sorts and scand the given per-dimension extract and puts all splits - * in a vector, that can easily be iterated afterwards. General implementation. - */ - template - void ExtractSplits(std::vector>& splitVec, - const MatType& data, - size_t dim, - const size_t start, - const size_t end, - const size_t minLeafSize) - { - static_assert( - std::is_same::value == true, - "The ElemType does not correspond to the matrix's element type."); - typedef std::pair SplitItem; - const typename MatType::row_type dimVec = +/** + * This one sorts and scand the given per-dimension extract and puts all splits + * in a vector, that can easily be iterated afterwards. General implementation. + */ +template +void ExtractSplits(std::vector>& splitVec, + const MatType& data, + size_t dim, + const size_t start, + const size_t end, + const size_t minLeafSize) +{ + static_assert( + std::is_same::value == true, + "The ElemType does not correspond to the matrix's element type."); + + typedef std::pair SplitItem; + const typename MatType::row_type dimVec = arma::sort(data(dim, arma::span(start, end - 1))); - // Ensure the minimum leaf size on both sides. We need to figure out why - // there are spikes if this minLeafSize is enforced here... - for (size_t i = minLeafSize - 1; i < dimVec.n_elem - minLeafSize; ++i) + // Ensure the minimum leaf size on both sides. We need to figure out why there + // are spikes if this minLeafSize is enforced here... + for (size_t i = minLeafSize - 1; i < dimVec.n_elem - minLeafSize; ++i) + { + // This makes sense for real continuous data. This kinda corrupts the data + // and estimation if the data is ordinal. Potentially we can fix that by + // taking into account ordinality later in the min/max update, but then we + // can end-up with a zero-volumed dimension. No good. + const ElemType split = (dimVec[i] + dimVec[i + 1]) / 2.0; + + // Check if we can split here (two points are different) + if (split != dimVec[i]) + splitVec.push_back(SplitItem(split, i + 1)); + } +} + +// Now the custom arma::Mat implementation. +template +void ExtractSplits(std::vector>& splitVec, + const arma::Mat& data, + size_t dim, + const size_t start, + const size_t end, + const size_t minLeafSize) +{ + typedef std::pair SplitItem; + arma::rowvec dimVec = data(dim, arma::span(start, end - 1)); + + // We sort these, in-place (it's a copy of the data, anyways). + std::sort(dimVec.begin(), dimVec.end()); + + for (size_t i = minLeafSize - 1; i < dimVec.n_elem - minLeafSize; ++i) + { + // This makes sense for real continuous data. This kinda corrupts the data + // and estimation if the data is ordinal. Potentially we can fix that by + // taking into account ordinality later in the min/max update, but then we + // can end-up with a zero-volumed dimension. No good. + const ElemType split = (dimVec[i] + dimVec[i + 1]) / 2.0; + + if (split != dimVec[i]) + splitVec.push_back(SplitItem(split, i + 1)); + } +} + +// This the custom, sparse optimized implementation of the same routine. +template +void ExtractSplits(std::vector>& splitVec, + const arma::SpMat& data, + size_t dim, + const size_t start, + const size_t end, + const size_t minLeafSize) +{ + // It's common sense, but we also use it in a check later. + Log::Assert(minLeafSize > 0); + + typedef std::pair SplitItem; + const size_t n_elem = end - start; + + // Construct a vector of values. + const arma::SpRow row = data(dim, arma::span(start, end - 1)); + std::vector valsVec(row.begin(), row.end()); + + // ... and sort it! + std::sort(valsVec.begin(), valsVec.end()); + + // Now iterate over the values, taking account for the over-the-zeroes jump + // and construct the splits vector. + const size_t zeroes = n_elem - valsVec.size(); + ElemType lastVal = -std::numeric_limits::max(); + size_t padding = 0; + + for (size_t i = 0; i < valsVec.size(); ++i) + { + const ElemType newVal = valsVec[i]; + if (lastVal < ElemType(0) && newVal > ElemType(0) && zeroes > 0) { - // This makes sense for real continuous data. This kinda corrupts the - // data and estimation if the data is ordinal. Potentially we can fix - // that by taking into account ordinality later in the min/max update, - // but then we can end-up with a zero-volumed dimension. No good. - const ElemType split = (dimVec[i] + dimVec[i + 1]) / 2.0; + Log::Assert(padding == 0); // We should arrive here once! + + // The minLeafSize > 0 also guarantees we're not entering right at the + // start. + if (i >= minLeafSize && i <= n_elem - minLeafSize) + splitVec.push_back(SplitItem(lastVal / 2.0, i)); + + padding = zeroes; + lastVal = ElemType(0); + } + + // This is the normal case. + if (i + padding >= minLeafSize && i + padding <= n_elem - minLeafSize) + { + // This makes sense for real continuous data. This kinda corrupts the + // data and estimation if the data is ordinal. Potentially we can fix that + // by taking into account ordinality later in the min/max update, but then + // we can end-up with a zero-volumed dimension. No good. + const ElemType split = (lastVal + newVal) / 2.0; // Check if we can split here (two points are different) - if (split != dimVec[i]) - splitVec.push_back(SplitItem(split, i + 1)); + if (split != newVal) + splitVec.push_back(SplitItem(split, i + padding)); } + + lastVal = newVal; } +} - // Now the custom arma::Mat implementation - template - void ExtractSplits(std::vector>& splitVec, - const arma::Mat& data, - size_t dim, - const size_t start, - const size_t end, - const size_t minLeafSize) - { - typedef std::pair SplitItem; - arma::vec dimVec = data(dim, arma::span(start, end - 1)).t(); +} // namespace details - // We sort these, in-place (it's a copy of the data, anyways). - std::sort(dimVec.begin(), dimVec.end()); - - for (size_t i = minLeafSize - 1; i < dimVec.n_elem - minLeafSize; ++i) - { - // This makes sense for real continuous data. This kinda corrupts the - // data and estimation if the data is ordinal. Potentially we can fix - // that by taking into account ordinality later in the min/max update, - // but then we can end-up with a zero-volumed dimension. No good. - const ElemType split = (dimVec[i] + dimVec[i + 1]) / 2.0; - - if (split != dimVec[i]) - splitVec.push_back(SplitItem(split, i + 1)); - } - } - - // This the custom, sparse optimized implementation of the same routine. - template - void ExtractSplits(std::vector>& splitVec, - const arma::SpMat& data, - size_t dim, - const size_t start, - const size_t end, - const size_t minLeafSize) - { - // It's common sense, but we also use it in a check later. - Log::Assert(minLeafSize > 0); - - typedef std::pair SplitItem; - const size_t n_elem = end - start; - - // Construct a vector of values. - const arma::SpRow row = data(dim, arma::span(start, end - 1)); - std::vector valsVec(row.begin(), row.end()); - - // ... and sort it! - std::sort(valsVec.begin(), valsVec.end()); - - // Now iterate over the values, taking account for the over-the-zeroes - // jump and construct the splits vector. - const size_t zeroes = n_elem - valsVec.size(); - ElemType lastVal = -std::numeric_limits::max(); - size_t padding = 0; - - for (size_t i = 0; i < valsVec.size(); ++i) - { - const ElemType newVal = valsVec[i]; - if (lastVal < ElemType(0) && newVal > ElemType(0) && zeroes > 0) - { - Log::Assert(padding == 0); // We should arrive here once! - - // The minLeafSize > 0 also guarantees we're not entering right at the - // start. - if (i >= minLeafSize && i <= n_elem - minLeafSize) - splitVec.push_back(SplitItem(lastVal / 2.0, i)); - - padding = zeroes; - lastVal = ElemType(0); - } - - // the normal case - if (i + padding >= minLeafSize && i + padding <= n_elem - minLeafSize) - { - // This makes sense for real continuous data. This kinda corrupts the - // data and estimation if the data is ordinal. Potentially we can fix - // that by taking into account ordinality later in the min/max update, - // but then we can end-up with a zero-volumed dimension. No good. - const ElemType split = (lastVal + newVal) / 2.0; - - // Check if we can split here (two points are different) - if (split != newVal) - splitVec.push_back(SplitItem(split, i + padding)); - } - - lastVal = newVal; - } - } -}; // namespace details - -template +template DTree::DTree() : start(0), end(0), @@ -165,7 +167,7 @@ DTree::DTree() : right(NULL) { /* Nothing to do. */ } -template +template DTree::DTree(const DTree& obj) : start(obj.start), end(obj.end), @@ -187,10 +189,13 @@ DTree::DTree(const DTree& obj) : /* Nothing to do. */ } -template +template DTree& DTree::operator=( const DTree& obj) { + if (this == &obj) + return *this; + // Copy the values from the other tree. start = obj.start; end = obj.end; @@ -213,12 +218,12 @@ DTree& DTree::operator=( // Copy the children. left = ((obj.left == NULL) ? NULL : new DTree(*obj.left)); - left = ((obj.right == NULL) ? NULL : new DTree(*obj.right)); + right = ((obj.right == NULL) ? NULL : new DTree(*obj.right)); return *this; } -template +template DTree::DTree(DTree&& obj): start(obj.start), end(obj.end), @@ -254,10 +259,13 @@ DTree::DTree(DTree&& obj): obj.right = NULL; } -template +template DTree& DTree::operator=( DTree&& obj) { + if (this == &obj) + return *this; + // Move the values from the other tree. start = obj.start; end = obj.end; @@ -302,8 +310,8 @@ DTree& DTree::operator=( } -// Root node initializers -template +// Root node initializers. +template DTree::DTree(const StatType& maxVals, const StatType& minVals, const size_t totalPoints) : @@ -325,7 +333,7 @@ DTree::DTree(const StatType& maxVals, right(NULL) { /* Nothing to do. */ } -template +template DTree::DTree(MatType & data) : start(0), end(data.n_cols), @@ -346,8 +354,8 @@ DTree::DTree(MatType & data) : logNegError = LogNegativeError(data.n_cols); } -// Non-root node initializers -template +// Non-root node initializers. +template DTree::DTree(const StatType& maxVals, const StatType& minVals, const size_t start, @@ -371,7 +379,7 @@ DTree::DTree(const StatType& maxVals, right(NULL) { /* Nothing to do. */ } -template +template DTree::DTree(const StatType& maxVals, const StatType& minVals, const size_t totalPoints, @@ -395,7 +403,7 @@ DTree::DTree(const StatType& maxVals, right(NULL) { /* Nothing to do. */ } -template +template DTree::~DTree() { delete left; @@ -404,7 +412,7 @@ DTree::~DTree() // This function computes the log-l2-negative-error of a given node from the // formula R(t) = log(|t|^2 / (N^2 V_t)). -template +template double DTree::LogNegativeError(const size_t totalPoints) const { // log(-|t|^2 / (N^2 V_t)) = log(-1) + 2 log(|t|) - 2 log(N) - log(V_t). @@ -425,7 +433,7 @@ double DTree::LogNegativeError(const size_t totalPoints) const // This function finds the best split with respect to the L2-error, by trying // all possible splits. The dataset is the full data set but the start and // end are used to obtain the point in this node. -template +template bool DTree::FindSplit(const MatType& data, size_t& splitDim, ElemType& splitValue, @@ -433,7 +441,7 @@ bool DTree::FindSplit(const MatType& data, double& rightError, const size_t minLeafSize) const { - typedef std::pair SplitItem; + typedef std::pair SplitItem; // Ensure the dimensionality of the data is the same as the dimensionality of // the bounding rectangle. @@ -543,7 +551,7 @@ bool DTree::FindSplit(const MatType& data, return splitFound; } -template +template size_t DTree::SplitData(MatType& data, const size_t splitDim, const ElemType splitValue, @@ -577,8 +585,8 @@ size_t DTree::SplitData(MatType& data, return left; } -// Greedily expand the tree -template +// Greedily expand the tree. +template double DTree::Grow(MatType& data, arma::Col& oldFromNew, const bool useVolReg, @@ -727,7 +735,7 @@ double DTree::Grow(MatType& data, } -template +template double DTree::PruneAndUpdate(const double oldAlpha, const size_t points, const bool useVolReg) @@ -842,7 +850,7 @@ double DTree::PruneAndUpdate(const double oldAlpha, // // Future improvement: Open up the range with epsilons on both sides where // epsilon depends on the density near the boundary. -template +template bool DTree::WithinRange(const VecType& query) const { for (size_t i = 0; i < query.n_elem; ++i) @@ -853,7 +861,7 @@ bool DTree::WithinRange(const VecType& query) const } -template +template double DTree::ComputeValue(const VecType& query) const { Log::Assert(query.n_elem == maxVals.n_elem); @@ -882,8 +890,8 @@ double DTree::ComputeValue(const VecType& query) const } // Index the buckets for possible usage later. -template -TagType DTree::TagTree(const TagType& tag) +template +TagType DTree::TagTree(const TagType& tag, bool every) { if (subtreeLeaves == 1) { @@ -891,18 +899,33 @@ TagType DTree::TagTree(const TagType& tag) bucketTag = tag; return (tag + 1); } - else + + TagType nextTag; + if (every) { - return right->TagTree(left->TagTree(tag)); + bucketTag = tag; + nextTag = (tag + 1); } + else + nextTag = tag; + + return right->TagTree(left->TagTree(nextTag, every), every); } -template +template TagType DTree::FindBucket(const VecType& query) const { Log::Assert(query.n_elem == maxVals.n_elem); - if (subtreeLeaves == 1) // If we are a leaf... + if (root == 1) // If we are the root... + { + // Check if the query is within range. + if (!WithinRange(query)) + return -1; + } + + // If we are a leaf... + if (subtreeLeaves == 1) { return bucketTag; } @@ -915,9 +938,9 @@ TagType DTree::FindBucket(const VecType& query) const } } -template -void -DTree::ComputeVariableImportance(arma::vec& importances) const +template +void DTree::ComputeVariableImportance(arma::vec& importances) + const { // Clear and set to right size. importances.zeros(maxVals.n_elem); @@ -944,8 +967,31 @@ DTree::ComputeVariableImportance(arma::vec& importances) const } } -template -template +template +void DTree::FillMinMax(const StatType& mins, + const StatType& maxs) +{ + if (!root) + { + minVals = mins; + maxVals = maxs; + } + + if (left && right) + { + StatType maxValsL(maxs); + StatType maxValsR(maxs); + StatType minValsL(mins); + StatType minValsR(mins); + + maxValsL[splitDim] = minValsR[splitDim] = splitValue; + left->FillMinMax(minValsL, maxValsL); + right->FillMinMax(minValsR, maxValsR); + } +} + +template +template void DTree::Serialize(Archive& ar, const unsigned int /* version */) { @@ -953,8 +999,6 @@ void DTree::Serialize(Archive& ar, ar & CreateNVP(start, "start"); ar & CreateNVP(end, "end"); - ar & CreateNVP(maxVals, "maxVals"); - ar & CreateNVP(minVals, "minVals"); ar & CreateNVP(splitDim, "splitDim"); ar & CreateNVP(splitValue, "splitValue"); ar & CreateNVP(logNegError, "logNegError"); @@ -976,5 +1020,14 @@ void DTree::Serialize(Archive& ar, ar & CreateNVP(left, "left"); ar & CreateNVP(right, "right"); -} + if (root) + { + ar & CreateNVP(maxVals, "maxVals"); + ar & CreateNVP(minVals, "minVals"); + + // This is added in order to reduce (dramatically!) the model file size. + if (Archive::is_loading::value && left && right) + FillMinMax(minVals, maxVals); + } +} diff --git a/src/mlpack/methods/gmm/diagonal_constraint.hpp b/src/mlpack/methods/gmm/diagonal_constraint.hpp index e0eafc02b7..864e2f5525 100644 --- a/src/mlpack/methods/gmm/diagonal_constraint.hpp +++ b/src/mlpack/methods/gmm/diagonal_constraint.hpp @@ -27,8 +27,7 @@ class DiagonalConstraint static void ApplyConstraint(arma::mat& covariance) { // Save the diagonal only. - arma::vec diagonal = covariance.diag(); - covariance = arma::diagmat(diagonal); + covariance = arma::diagmat(arma::clamp(covariance.diag(), 1e-10, DBL_MAX)); } //! Serialize the constraint (which holds nothing, so, nothing to do). diff --git a/src/mlpack/methods/gmm/em_fit.hpp b/src/mlpack/methods/gmm/em_fit.hpp index 71817f8cab..82a646db70 100644 --- a/src/mlpack/methods/gmm/em_fit.hpp +++ b/src/mlpack/methods/gmm/em_fit.hpp @@ -162,6 +162,26 @@ class EMFit dists, const arma::vec& weights) const; + // Armadillo uses uword internally as an OpenMP index type, which crashes + // Visual Studio. + #ifndef _WIN32 + /** + * Use the Armadillo gmm_diag clusterer to train a GMM with diagonal + * covariance. If InitialClusteringType == kmeans::KMeans<>, this will use + * Armadillo's initialization also. + * + * @param observations Data to train on. + * @param dists Distributions to store model in. + * @param weights Prior weights. + * @param useInitialModel If true, the existing model will be used. + */ + void ArmadilloGMMWrapper( + const arma::mat& observations, + std::vector& dists, + arma::vec& weights, + const bool useInitialModel); + #endif + //! Maximum iterations of EM algorithm. size_t maxIterations; //! Tolerance for convergence of EM. diff --git a/src/mlpack/methods/gmm/em_fit_impl.hpp b/src/mlpack/methods/gmm/em_fit_impl.hpp index feaa9d4ae1..40d7013fbf 100644 --- a/src/mlpack/methods/gmm/em_fit_impl.hpp +++ b/src/mlpack/methods/gmm/em_fit_impl.hpp @@ -15,6 +15,7 @@ // In case it hasn't been included yet. #include "em_fit.hpp" +#include "diagonal_constraint.hpp" namespace mlpack { namespace gmm { @@ -39,6 +40,17 @@ void EMFit::Estimate( arma::vec& weights, const bool useInitialModel) { + // Shortcut: if the user is using the DiagonalConstraint, then we will call + // out to Armadillo. But Armadillo uses uword internally as an OpenMP index + // type, which crashes Visual Studio, so don't do this on Windows. + #ifndef _WIN32 + if (std::is_same::value) + { + ArmadilloGMMWrapper(observations, dists, weights, useInitialModel); + return; + } + #endif + // Only perform initial clustering if the user wanted it. if (!useInitialModel) InitialClustering(observations, dists, weights); @@ -316,6 +328,68 @@ void EMFit::Serialize( ar & CreateNVP(constraint, "constraint"); } +// Armadillo uses uword internally as an OpenMP index type, which crashes Visual +// Studio. +#ifndef _WIN32 +template +void EMFit:: +ArmadilloGMMWrapper(const arma::mat& observations, + std::vector& dists, + arma::vec& weights, + const bool useInitialModel) +{ + arma::gmm_diag g; + + // Warn the user that tolerance isn't used for convergence here if they've + // specified a non-default value. + if (tolerance != EMFit().Tolerance()) + Log::Warn << "GMM::Train(): tolerance ignored when training GMMs with " + << "DiagonalConstraint." << std::endl; + + // If the initial clustering is the default k-means, we'll just use + // Armadillo's implementation. If mlpack ever changes k-means defaults to use + // something that is reliably quicker than the Lloyd iteration k-means update, + // then this code maybe should be revisited. + if (!std::is_same>::value || + useInitialModel) + { + // Use clusterer to get initial values. + if (!useInitialModel) + InitialClustering(observations, dists, weights); + + // Assemble matrix of means. + arma::mat means(observations.n_rows, dists.size()); + arma::mat covs(observations.n_rows, dists.size()); + for (size_t i = 0; i < dists.size(); ++i) + { + means.col(i) = dists[i].Mean(); + covs.col(i) = dists[i].Covariance().diag(); + } + + g.reset(observations.n_rows, dists.size()); + g.set_params(std::move(means), std::move(covs), weights.t()); + + g.learn(observations, dists.size(), arma::eucl_dist, arma::keep_existing, 0, + maxIterations, 1e-10, false /* no printing */); + } + else + { + // Use Armadillo for the initial clustering. We'll try and match mlpack + // defaults. + g.learn(observations, dists.size(), arma::eucl_dist, arma::static_subset, + 1000, maxIterations, 1e-10, false /* no printing */); + } + + // Extract means, covariances, and weights. + weights = g.hefts.t(); + for (size_t i = 0; i < dists.size(); ++i) + { + dists[i].Mean() = g.means.col(i); + dists[i].Covariance(std::move(arma::diagmat(g.dcovs.col(i)))); + } +} +#endif + } // namespace gmm } // namespace mlpack diff --git a/src/mlpack/methods/gmm/gmm_train_main.cpp b/src/mlpack/methods/gmm/gmm_train_main.cpp index 6f81b3e639..d8a5885ad2 100644 --- a/src/mlpack/methods/gmm/gmm_train_main.cpp +++ b/src/mlpack/methods/gmm/gmm_train_main.cpp @@ -15,6 +15,7 @@ #include "gmm.hpp" #include "no_constraint.hpp" +#include "diagonal_constraint.hpp" #include @@ -52,6 +53,11 @@ PROGRAM_INFO("Gaussian Mixture Model (GMM) Training", "Bradley-Fayyad refined start initialization will be used. This can often " "lead to better clustering results." "\n\n" + "The 'diagonal_covariance' flag will cause the learned covariances to be " + "diagonal matrices. This significantly simplifies the model itself and " + "causes training to be faster, but restricts the ability to fit more " + "complex GMMs." + "\n\n" "If GMM training fails with an error indicating that a covariance matrix " "could not be inverted, make sure that the " + PRINT_PARAM_STRING("no_force_positive") + " parameter is not " @@ -94,6 +100,8 @@ PARAM_FLAG("no_force_positive", "Do not force the covariance matrices to be " "positive definite.", "P"); PARAM_INT_IN("max_iterations", "Maximum number of iterations of EM algorithm " "(passing 0 will run until convergence).", "n", 250); +PARAM_FLAG("diagonal_covariance", "Force the covariance of the Gaussians to " + "be diagonal. This can accelerate training time significantly.", "d"); // Parameters for dataset modification. PARAM_DOUBLE_IN("noise", "Variance of zero-mean Gaussian noise to add to data.", @@ -128,6 +136,11 @@ void mlpackMain() "be greater than or equal to 1." << std::endl; } + if (CLI::HasParam("diagonal_covariance") && + CLI::HasParam("no_force_positive")) + Log::Warn << "--no_force_positive ignored because --diagonal_covariance is " + << "specified!" << endl; + if (!CLI::HasParam("output_model")) Log::Warn << "--output_model_file is not specified, so no model will be " << "saved!" << endl; @@ -163,6 +176,7 @@ void mlpackMain() const size_t maxIterations = (size_t) CLI::GetParam("max_iterations"); const double tolerance = CLI::GetParam("tolerance"); const bool forcePositive = !CLI::HasParam("no_force_positive"); + const bool diagonalCovariance = CLI::HasParam("diagonal_covariance"); // This gets a bit weird because we need different types depending on whether // --refined_start is specified. @@ -186,9 +200,18 @@ void mlpackMain() KMeansType k(1000, metric::SquaredEuclideanDistance(), RefinedStart(samplings, percentage)); - // Depending on the value of 'forcePositive', we have to use different - // types. - if (forcePositive) + // Depending on the value of forcePositive and diagonalCovariance, we have + // to use different types. + if (diagonalCovariance) + { + // Compute the parameters of the model using the EM algorithm. + Timer::Start("em"); + EMFit em(maxIterations, tolerance, k); + likelihood = gmm.Train(dataPoints, CLI::GetParam("trials"), false, + em); + Timer::Stop("em"); + } + else if (forcePositive) { // Compute the parameters of the model using the EM algorithm. Timer::Start("em"); @@ -209,8 +232,18 @@ void mlpackMain() } else { - // Depending on the value of forcePositive, we have to use different types. - if (forcePositive) + // Depending on the value of forcePositive and diagonalCovariance, we have + // to use different types. + if (diagonalCovariance) + { + // Compute the parameters of the model using the EM algorithm. + Timer::Start("em"); + EMFit, DiagonalConstraint> em(maxIterations, tolerance); + likelihood = gmm.Train(dataPoints, CLI::GetParam("trials"), false, + em); + Timer::Stop("em"); + } + else if (forcePositive) { // Compute the parameters of the model using the EM algorithm. Timer::Start("em"); diff --git a/src/mlpack/methods/rann/ra_util.cpp b/src/mlpack/methods/rann/ra_util.cpp index 867d5908fe..3788917261 100644 --- a/src/mlpack/methods/rann/ra_util.cpp +++ b/src/mlpack/methods/rann/ra_util.cpp @@ -30,11 +30,6 @@ size_t mlpack::neighbor::RAUtil::MinimumSamplesReqd(const size_t n, double prob; Log::Assert(alpha <= 1.0); - // going through all values of sample sizes - // to find the minimum samples required to satisfy the - // desired bound - bool done = false; - // This performs a binary search on the integer values between 'lb = k' // and 'ub = n' to find the minimum number of samples 'm' required to obtain // the desired success probability 'alpha'. @@ -46,7 +41,6 @@ size_t mlpack::neighbor::RAUtil::MinimumSamplesReqd(const size_t n, { if (prob - alpha < 0.001 || ub < lb + 2) { - done = true; break; } else @@ -66,12 +60,11 @@ size_t mlpack::neighbor::RAUtil::MinimumSamplesReqd(const size_t n, } else { - done = true; break; } } m = (ub + lb) / 2; - } while (!done); + } while (true); return (std::min(m + 1, n)); } diff --git a/src/mlpack/tests/gmm_test.cpp b/src/mlpack/tests/gmm_test.cpp index 1d38323e00..ae294ca253 100644 --- a/src/mlpack/tests/gmm_test.cpp +++ b/src/mlpack/tests/gmm_test.cpp @@ -756,5 +756,105 @@ BOOST_AUTO_TEST_CASE(UseExistingModelTest) } } +/** + * Make sure we can fit a diagonal GMM reasonably. + */ +BOOST_AUTO_TEST_CASE(DiagonalGMMTrainTest) +{ + Log::Warn.ignoreInput = false; + // We'll have three diagonal-covariance Gaussian distributions from this + // mixture. + distribution::GaussianDistribution d1("0.0 1.0 0.0", "1.0 0.0 0.0;" + "0.0 0.8 0.0;" + "0.0 0.0 1.0"); + distribution::GaussianDistribution d2("2.0 -1.0 5.0", "3.0 0.0 0.0;" + "0.0 1.2 0.0;" + "0.0 0.0 1.3"); + distribution::GaussianDistribution d3("0.0 5.0 -3.0", "2.0 0.0 0.0;" + "0.0 0.3 0.0;" + "0.0 0.0 1.0"); + + // Now we'll generate points and probabilities. 1500 points. Slower than I + // would like... + arma::mat points(3, 5000); + + for (size_t i = 0; i < 5000; i++) + { + double randValue = math::Random(); + + if (randValue <= 0.20) // p(d1) = 0.20 + points.col(i) = d1.Random(); + else if (randValue <= 0.50) // p(d2) = 0.30 + points.col(i) = d2.Random(); + else // p(d3) = 0.50 + points.col(i) = d3.Random(); + } + + // Now train the model. 3 dimensions, 3 components. + GMM g(3, 3); + + g.Train, DiagonalConstraint>>(points, 5); + + // Now check the results. We need to order by weights so that when we do the + // checking, things will be correct. + arma::uvec sortedIndices = sort_index(g.Weights()); + + // First Gaussian (d1). + BOOST_REQUIRE_SMALL(g.Weights()[sortedIndices[0]] - 0.2, 0.1); + + for (size_t i = 0; i < 3; i++) + BOOST_REQUIRE_SMALL((g.Component(sortedIndices[0]).Mean()[i] + - d1.Mean()[i]), 0.4); + + for (size_t row = 0; row < 3; ++row) + { + for (size_t col = 0; col < 3; ++col) + { + const double v = g.Component(sortedIndices[0]).Covariance()(row, col); + if (row == col) + BOOST_REQUIRE_SMALL(v - d1.Covariance()(row, col), 0.5); + else + BOOST_REQUIRE_SMALL(v, 1e-5); + } + } + + // Second Gaussian (d2). + BOOST_REQUIRE_SMALL(g.Weights()[sortedIndices[1]] - 0.3, 0.1); + + for (size_t i = 0; i < 3; i++) + BOOST_REQUIRE_SMALL((g.Component(sortedIndices[1]).Mean()[i] + - d2.Mean()[i]), 0.4); + + for (size_t row = 0; row < 3; ++row) + { + for (size_t col = 0; col < 3; ++col) + { + const double v = g.Component(sortedIndices[1]).Covariance()(row, col); + if (row == col) + BOOST_REQUIRE_SMALL(v - d2.Covariance()(row, col), 0.5); + else + BOOST_REQUIRE_SMALL(v, 1e-5); + } + } + + // Third Gaussian (d3). + BOOST_REQUIRE_SMALL(g.Weights()[sortedIndices[2]] - 0.5, 0.1); + + for (size_t i = 0; i < 3; ++i) + BOOST_REQUIRE_SMALL((g.Component(sortedIndices[2]).Mean()[i] + - d3.Mean()[i]), 0.4); + + for (size_t row = 0; row < 3; ++row) + { + for (size_t col = 0; col < 3; ++col) + { + const double v = g.Component(sortedIndices[2]).Covariance()(row, col); + if (row == col) + BOOST_REQUIRE_SMALL(v - d3.Covariance()(row, col), 0.5); + else + BOOST_REQUIRE_SMALL(v, 1e-5); + } + } +} BOOST_AUTO_TEST_SUITE_END();