diff --git a/src/mlpack/core.hpp b/src/mlpack/core.hpp index 132ed5c19d..101f17d5f3 100644 --- a/src/mlpack/core.hpp +++ b/src/mlpack/core.hpp @@ -81,6 +81,7 @@ #include #include #include +#include #include #include #include diff --git a/src/mlpack/core/arma_extend/CMakeLists.txt b/src/mlpack/core/arma_extend/CMakeLists.txt index fc218442f4..30195a9b76 100644 --- a/src/mlpack/core/arma_extend/CMakeLists.txt +++ b/src/mlpack/core/arma_extend/CMakeLists.txt @@ -2,13 +2,8 @@ # Anything not in this list will not be compiled into mlpack. set(SOURCES arma_extend.hpp - fn_ccov.hpp fn_inplace_reshape.hpp - glue_ccov_meat.hpp - glue_ccov_proto.hpp hdf5_misc.hpp - op_ccov_meat.hpp - op_ccov_proto.hpp SpMat_extra_bones.hpp SpMat_extra_meat.hpp Mat_extra_bones.hpp diff --git a/src/mlpack/core/arma_extend/arma_extend.hpp b/src/mlpack/core/arma_extend/arma_extend.hpp index bcf6088a4b..c2aa4f5db2 100644 --- a/src/mlpack/core/arma_extend/arma_extend.hpp +++ b/src/mlpack/core/arma_extend/arma_extend.hpp @@ -5,9 +5,6 @@ * Include Armadillo extensions which currently are not part of the main * Armadillo codebase. * - * This will allow the use of the ccov() function (which performs the same - * function as cov(trans(X)) but without the cost of computing trans(X)). This - * also gives sparse matrix support, if it is necessary. */ #ifndef MLPACK_CORE_ARMA_EXTEND_ARMA_EXTEND_HPP #define MLPACK_CORE_ARMA_EXTEND_ARMA_EXTEND_HPP @@ -55,12 +52,6 @@ namespace arma { // u64/s64 #include "hdf5_misc.hpp" - // ccov() - #include "op_ccov_proto.hpp" - #include "op_ccov_meat.hpp" - #include "glue_ccov_proto.hpp" - #include "glue_ccov_meat.hpp" - #include "fn_ccov.hpp" // inplace_reshape() #include "fn_inplace_reshape.hpp" diff --git a/src/mlpack/core/arma_extend/fn_ccov.hpp b/src/mlpack/core/arma_extend/fn_ccov.hpp deleted file mode 100644 index 86f3ecb81d..0000000000 --- a/src/mlpack/core/arma_extend/fn_ccov.hpp +++ /dev/null @@ -1,34 +0,0 @@ -//! \addtogroup fn_ccov -//! @{ - - - -template -inline -const Op -ccov(const Base& X, const uword norm_type = 0) - { - arma_extra_debug_sigprint(); - - arma_debug_check( (norm_type > 1), "ccov(): norm_type must be 0 or 1"); - - return Op(X.get_ref(), norm_type, 0); - } - - - -template -inline -const Glue -cov(const Base& A, const Base& B, const uword norm_type = 0) - { - arma_extra_debug_sigprint(); - - arma_debug_check( (norm_type > 1), "ccov(): norm_type must be 0 or 1"); - - return Glue(A.get_ref(), B.get_ref(), norm_type); - } - - - -//! @} diff --git a/src/mlpack/core/arma_extend/glue_ccov_meat.hpp b/src/mlpack/core/arma_extend/glue_ccov_meat.hpp deleted file mode 100644 index c3589c368b..0000000000 --- a/src/mlpack/core/arma_extend/glue_ccov_meat.hpp +++ /dev/null @@ -1,144 +0,0 @@ -//! \addtogroup glue_cov -//! @{ - - - -template -inline -void -glue_ccov::direct_ccov(Mat& out, const Mat& A, const Mat& B, const uword norm_type) - { - arma_extra_debug_sigprint(); - - if(A.is_vec() && B.is_vec()) - { - arma_debug_check( (A.n_elem != B.n_elem), "ccov(): the number of elements in A and B must match" ); - - const eT* A_ptr = A.memptr(); - const eT* B_ptr = B.memptr(); - - eT A_acc = eT(0); - eT B_acc = eT(0); - eT out_acc = eT(0); - - const uword N = A.n_elem; - - for(uword i=0; i 1) ? eT(N-1) : eT(1) ) : eT(N); - - out.set_size(1,1); - out[0] = out_acc/norm_val; - } - else - { - arma_debug_assert_same_size(A, B, "ccov()"); - - const uword N = A.n_cols; - const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); - - out = A * trans(B); - out -= (sum(A) * trans(sum(B))) / eT(N); - out /= norm_val; - } - } - - - -template -inline -void -glue_ccov::direct_ccov(Mat< std::complex >& out, const Mat< std::complex >& A, const Mat< std::complex >& B, const uword norm_type) - { - arma_extra_debug_sigprint(); - - typedef typename std::complex eT; - - if(A.is_vec() && B.is_vec()) - { - arma_debug_check( (A.n_elem != B.n_elem), "cov(): the number of elements in A and B must match" ); - - const eT* A_ptr = A.memptr(); - const eT* B_ptr = B.memptr(); - - eT A_acc = eT(0); - eT B_acc = eT(0); - eT out_acc = eT(0); - - const uword N = A.n_elem; - - for(uword i=0; i 1) ? eT(N-1) : eT(1) ) : eT(N); - - out.set_size(1,1); - out[0] = out_acc/norm_val; - } - else - { - arma_debug_assert_same_size(A, B, "ccov()"); - - const uword N = A.n_cols; - const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); - - out = A * trans(conj(B)); - out -= (sum(A) * trans(conj(sum(B)))) / eT(N); - out /= norm_val; - } - } - - - -template -inline -void -glue_ccov::apply(Mat& out, const Glue& X) - { - arma_extra_debug_sigprint(); - - typedef typename T1::elem_type eT; - - const unwrap_check A_tmp(X.A, out); - const unwrap_check B_tmp(X.B, out); - - const Mat& A = A_tmp.M; - const Mat& B = B_tmp.M; - - const uword norm_type = X.aux_uword; - - if(&A != &B) - { - glue_ccov::direct_ccov(out, A, B, norm_type); - } - else - { - op_ccov::direct_ccov(out, A, norm_type); - } - - } - - - -//! @} diff --git a/src/mlpack/core/arma_extend/glue_ccov_proto.hpp b/src/mlpack/core/arma_extend/glue_ccov_proto.hpp deleted file mode 100644 index f5531175de..0000000000 --- a/src/mlpack/core/arma_extend/glue_ccov_proto.hpp +++ /dev/null @@ -1,15 +0,0 @@ -//! \addtogroup glue_ccov -//! @{ - -class glue_ccov - { - public: - - template inline static void direct_ccov(Mat& out, const Mat& A, const Mat& B, const uword norm_type); - template inline static void direct_ccov(Mat< std::complex >& out, const Mat< std::complex >& A, const Mat< std::complex >& B, const uword norm_type); - - template inline static void apply(Mat& out, const Glue& X); - }; - -//! @} - diff --git a/src/mlpack/core/arma_extend/op_ccov_meat.hpp b/src/mlpack/core/arma_extend/op_ccov_meat.hpp deleted file mode 100644 index 93c09f280a..0000000000 --- a/src/mlpack/core/arma_extend/op_ccov_meat.hpp +++ /dev/null @@ -1,97 +0,0 @@ -//! \addtogroup op_cov -//! @{ - - - -template -inline -void -op_ccov::direct_ccov(Mat& out, const Mat& A, const uword norm_type) - { - arma_extra_debug_sigprint(); - - if(A.is_vec()) - { - if(A.n_rows == 1) - { - out = var(trans(A), norm_type); - } - else - { - out = var(A, norm_type); - } - } - else - { - const uword N = A.n_cols; - const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); - - const Col acc = sum(A, 1); - - out = A * trans(A); - out -= (acc * trans(acc)) / eT(N); - out /= norm_val; - } - } - - - -template -inline -void -op_ccov::direct_ccov(Mat< std::complex >& out, const Mat< std::complex >& A, const uword norm_type) - { - arma_extra_debug_sigprint(); - - typedef typename std::complex eT; - - if(A.is_vec()) - { - if(A.n_rows == 1) - { - const Mat tmp_mat = var(trans(A), norm_type); - out.set_size(1,1); - out[0] = tmp_mat[0]; - } - else - { - const Mat tmp_mat = var(A, norm_type); - out.set_size(1,1); - out[0] = tmp_mat[0]; - } - } - else - { - const uword N = A.n_cols; - const eT norm_val = (norm_type == 0) ? ( (N > 1) ? eT(N-1) : eT(1) ) : eT(N); - - const Col acc = sum(A, 1); - - out = A * trans(conj(A)); - out -= (acc * trans(conj(acc))) / eT(N); - out /= norm_val; - } - } - - - -template -inline -void -op_ccov::apply(Mat& out, const Op& in) - { - arma_extra_debug_sigprint(); - - typedef typename T1::elem_type eT; - - const unwrap_check tmp(in.m, out); - const Mat& A = tmp.M; - - const uword norm_type = in.aux_uword_a; - - op_ccov::direct_ccov(out, A, norm_type); - } - - - -//! @} diff --git a/src/mlpack/core/arma_extend/op_ccov_proto.hpp b/src/mlpack/core/arma_extend/op_ccov_proto.hpp deleted file mode 100644 index 4fb49eb65a..0000000000 --- a/src/mlpack/core/arma_extend/op_ccov_proto.hpp +++ /dev/null @@ -1,18 +0,0 @@ -//! \addtogroup op_cov -//! @{ - - - -class op_ccov - { - public: - - template inline static void direct_ccov(Mat& out, const Mat& X, const uword norm_type); - template inline static void direct_ccov(Mat< std::complex >& out, const Mat< std::complex >& X, const uword norm_type); - - template inline static void apply(Mat& out, const Op& in); - }; - - - -//! @} diff --git a/src/mlpack/core/math/CMakeLists.txt b/src/mlpack/core/math/CMakeLists.txt index 5188956838..6bacd597f8 100644 --- a/src/mlpack/core/math/CMakeLists.txt +++ b/src/mlpack/core/math/CMakeLists.txt @@ -18,6 +18,8 @@ set(SOURCES range_impl.hpp round.hpp shuffle_data.hpp + ccov.hpp + ccov_impl.hpp ) # add directory name to sources diff --git a/src/mlpack/core/math/ccov.hpp b/src/mlpack/core/math/ccov.hpp new file mode 100644 index 0000000000..d2e24f8ed1 --- /dev/null +++ b/src/mlpack/core/math/ccov.hpp @@ -0,0 +1,39 @@ +/** + * @file ccov.hpp + * @author Ryan Curtin + * @author Conrad Sanderson + * + * ColumnCovariance(X) is same as cov(trans(X)) but without the cost + * of computing trans(X) + * + * 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_MATH_CCOV_HPP +#define MLPACK_CORE_MATH_CCOV_HPP + +#include + +namespace mlpack { +namespace math /** Miscellaneous math routines. */ { + +template +inline +arma::Mat +ColumnCovariance(const arma::Mat& A, const size_t norm_type = 0); + +template +inline +arma::Mat< std::complex > +ColumnCovariance(const arma::Mat< std::complex >& A, + const size_t norm_type = 0); + +} // namespace math +} // namespace mlpack + +// Include implementation +#include "ccov_impl.hpp" + +#endif // MLPACK_CORE_MATH_CCOV_HPP diff --git a/src/mlpack/core/math/ccov_impl.hpp b/src/mlpack/core/math/ccov_impl.hpp new file mode 100644 index 0000000000..1a9103f995 --- /dev/null +++ b/src/mlpack/core/math/ccov_impl.hpp @@ -0,0 +1,100 @@ +/** + * @file ccov_impl.hpp + * @author Ryan Curtin + * @author Conrad Sanderson + * + * ColumnCovariance(X) is same as cov(trans(X)) but without the cost of computing trans(X) + * + * 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_MATH_CCOV_IMPL_HPP +#define MLPACK_CORE_MATH_CCOV_IMPL_HPP + +#include "ccov.hpp" + +namespace mlpack { +namespace math /** Miscellaneous math routines. */ { + +template +inline arma::Mat ColumnCovariance(const arma::Mat& x, + const size_t normType) +{ + if (normType > 1) + { + Log::Fatal << "ColumnCovariance(): norm_type must be 0 or 1!" << std::endl; + } + + arma::Mat out; + + if (x.n_elem > 0) + { + const arma::Mat& xAlias = (x.n_cols == 1) ? + arma::Mat(const_cast(x.memptr()), x.n_cols, x.n_rows, false, + false) : + arma::Mat(const_cast(x.memptr()), x.n_rows, x.n_cols, false, + false); + + const size_t n = xAlias.n_cols; + const eT normVal = (normType == 0) ? ((n > 1) ? eT(n - 1) : eT(1)) : eT(n); + + const arma::Mat tmp = xAlias.each_col() - arma::mean(xAlias, 1); + + out = tmp * tmp.t(); + out /= normVal; + } + + return out; +} + +template +inline arma::Mat> ColumnCovariance( + const arma::Mat>& x, + const size_t normType) +{ + if (normType > 1) + { + Log::Fatal << "ColumnCovariance(): norm_type must be 0 or 1" << std::endl; + } + + typedef typename std::complex eT; + + arma::Mat out; + + if (x.is_vec()) + { + if (x.n_rows == 1) + { + const arma::Mat tmpMat = arma::var(arma::trans(x), normType); + out.set_size(1, 1); + out[0] = tmpMat[0]; + } + else + { + const arma::Mat tmpMat = arma::var(x, normType); + out.set_size(1, 1); + out[0] = tmpMat[0]; + } + } + else + { + const size_t n = x.n_cols; + const eT normVal = (normType == 0) ? + ((n > 1) ? eT(n - 1) : eT(1)) : eT(n); + + const arma::Col acc = arma::sum(x, 1); + + out = x * arma::trans(arma::conj(x)); + out -= (acc * arma::trans(arma::conj(acc))) / eT(n); + out /= normVal; + } + + return out; +} + +} // namespace math +} // namespace mlpack + +#endif // MLPACK_CORE_MATH_CCOV_IMPL_HPP diff --git a/src/mlpack/core/math/lin_alg.cpp b/src/mlpack/core/math/lin_alg.cpp index 915d1a36a9..822625a59b 100644 --- a/src/mlpack/core/math/lin_alg.cpp +++ b/src/mlpack/core/math/lin_alg.cpp @@ -10,7 +10,7 @@ * http://www.opensource.org/licenses/BSD-3-Clause for more information. */ #include "lin_alg.hpp" -#include +#include #include using namespace mlpack; @@ -60,7 +60,7 @@ void mlpack::math::WhitenUsingSVD(const arma::mat& x, arma::mat covX, u, v, invSMatrix, temp1; arma::vec sVector; - covX = ccov(x); + covX = mlpack::math::ColumnCovariance(x); svd(u, sVector, v, covX); @@ -85,7 +85,7 @@ void mlpack::math::WhitenUsingEig(const arma::mat& x, arma::vec eigenvalues; // Get eigenvectors of covariance of input matrix. - eig_sym(eigenvalues, eigenvectors, ccov(x)); + eig_sym(eigenvalues, eigenvectors, mlpack::math::ColumnCovariance(x)); // Generate diagonal matrix using 1 / sqrt(eigenvalues) for each value. VectorPower(eigenvalues, -0.5); @@ -135,7 +135,7 @@ void mlpack::math::Orthogonalize(const arma::mat& x, arma::mat& W) // eigendecomposition of the matrix A. arma::mat eigenvalues, eigenvectors; arma::vec egval; - eig_sym(egval, eigenvectors, ccov(x)); + eig_sym(egval, eigenvectors, mlpack::math::ColumnCovariance(x)); VectorPower(egval, -0.5); eigenvalues.zeros(egval.n_elem, egval.n_elem); diff --git a/src/mlpack/tests/distribution_test.cpp b/src/mlpack/tests/distribution_test.cpp index 7add802d9e..83ecda846e 100644 --- a/src/mlpack/tests/distribution_test.cpp +++ b/src/mlpack/tests/distribution_test.cpp @@ -461,7 +461,7 @@ BOOST_AUTO_TEST_CASE(GaussianDistributionRandomTest) // Now make sure that reflects the actual distribution. arma::vec obsMean = arma::mean(obs, 1); - arma::mat obsCov = ccov(obs); + arma::mat obsCov = mlpack::math::ColumnCovariance(obs); // 10% tolerance because this can be noisy. BOOST_REQUIRE_CLOSE(obsMean[0], mean[0], 10.0); @@ -496,7 +496,7 @@ BOOST_AUTO_TEST_CASE(GaussianDistributionTrainTest) // Find actual mean and covariance of data. arma::vec actualMean = arma::mean(observations, 1); - arma::mat actualCov = ccov(observations); + arma::mat actualCov = mlpack::math::ColumnCovariance(observations); d.Train(observations); @@ -1418,7 +1418,7 @@ BOOST_AUTO_TEST_CASE(DiagonalGaussianDistributionRandomTest) // Make sure that reflects the actual distribution. arma::vec obsMean = arma::mean(obs, 1); - arma::mat obsCov = arma::ccov(obs); + arma::mat obsCov = mlpack::math::ColumnCovariance(obs); // 10% tolerance because this can be noisy. BOOST_REQUIRE_CLOSE(obsMean(0), mean(0), 10.0); @@ -1446,7 +1446,7 @@ BOOST_AUTO_TEST_CASE(DiagonalGaussianDistributionTrainTest) // Calculate the actual mean and covariance of data using armadillo. arma::vec actualMean = arma::mean(observations, 1); - arma::mat actualCov = arma::ccov(observations); + arma::mat actualCov = mlpack::math::ColumnCovariance(observations); // Estimate the parameters. d.Train(observations); diff --git a/src/mlpack/tests/gmm_test.cpp b/src/mlpack/tests/gmm_test.cpp index fd8edd8d4e..f0a523bdff 100644 --- a/src/mlpack/tests/gmm_test.cpp +++ b/src/mlpack/tests/gmm_test.cpp @@ -109,7 +109,8 @@ BOOST_AUTO_TEST_CASE(GMMTrainEMOneGaussian) gmm.Train(data, 10); arma::vec actualMean = arma::mean(data, 1); - arma::mat actualCovar = ccov(data, 1 /* biased estimator */); + arma::mat actualCovar = mlpack::math::ColumnCovariance(data, + 1 /* biased estimator */); // Check the model to see that it is correct. BOOST_REQUIRE_LT(arma::norm(gmm.Component(0).Mean() - actualMean), 1e-5); @@ -198,7 +199,8 @@ BOOST_AUTO_TEST_CASE(GMMTrainEMMultipleGaussians) // Calculate the actual means and covariances because they will probably // be different (this is easier to do before we shuffle the points). means[i] = arma::mean(data.cols(point, point + counts[i] - 1), 1); - covars[i] = ccov(data.cols(point, point + counts[i] - 1), 1 /* biased */); + covars[i] = mlpack::math::ColumnCovariance(arma::mat(data.cols(point, + point + counts[i] - 1)), 1 /* biased */); point += counts[i]; } @@ -694,7 +696,8 @@ BOOST_AUTO_TEST_CASE(UseExistingModelTest) // Calculate the actual means and covariances because they will probably // be different (this is easier to do before we shuffle the points). means[i] = arma::mean(data.cols(point, point + counts[i] - 1), 1); - covars[i] = ccov(data.cols(point, point + counts[i] - 1), 1 /* biased */); + covars[i] = mlpack::math::ColumnCovariance(arma::mat(data.cols(point, + point + counts[i] - 1)), 1 /* biased */); point += counts[i]; } @@ -855,7 +858,8 @@ BOOST_AUTO_TEST_CASE(DiagonalGMMTrainEMOneGaussian) arma::vec actualMean = arma::mean(data, 1); arma::vec actualCovar = arma::diagvec( - arma::ccov(data, 1 /* biased estimator */)); + mlpack::math::ColumnCovariance(data, + 1 /* biased estimator */)); // Check the model to see that it is correct. CheckMatrices(gmm.Component(0).Mean(), actualMean); diff --git a/src/mlpack/tests/hmm_test.cpp b/src/mlpack/tests/hmm_test.cpp index aa26425512..903f7374ab 100644 --- a/src/mlpack/tests/hmm_test.cpp +++ b/src/mlpack/tests/hmm_test.cpp @@ -1410,7 +1410,8 @@ BOOST_AUTO_TEST_CASE(DiagonalGMMHMMOneGaussianOneStateTrainingTest) // Generate the ground truth values. arma::vec actualMean = arma::mean(observations[0], 1); arma::vec actualCovar = arma::diagvec( - arma::ccov(observations[0], 1 /* biased estimator */)); + mlpack::math::ColumnCovariance(observations[0], + 1 /* biased estimator */)); // Check the model to see that it is correct. CheckMatrices(hmm.Emission()[0].Component(0).Mean(), actualMean); diff --git a/src/mlpack/tests/lin_alg_test.cpp b/src/mlpack/tests/lin_alg_test.cpp index 0e335049c9..cd1f3ced85 100644 --- a/src/mlpack/tests/lin_alg_test.cpp +++ b/src/mlpack/tests/lin_alg_test.cpp @@ -89,7 +89,7 @@ BOOST_AUTO_TEST_CASE(TestWhitenUsingEig) Center(tmp, tmp_centered); WhitenUsingEig(tmp_centered, whitened, whitening_matrix); - mat newcov = ccov(whitened); + mat newcov = mlpack::math::ColumnCovariance(whitened); for (int row = 0; row < 5; row++) { for (int col = 0; col < 5; col++) @@ -118,7 +118,7 @@ BOOST_AUTO_TEST_CASE(TestOrthogonalize) Orthogonalize(tmp, orth); // test orthogonality - mat test = ccov(orth); + mat test = mlpack::math::ColumnCovariance(orth); double ival = test(0, 0); for (size_t row = 0; row < test.n_rows; row++) {