apaam/mlpack
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge branch 'jeffin143-ccov'

Browse Source
This commit is contained in:
Ryan Curtin
2019-04-23 13:25:52 -04:00
parent 354938177a c45afea4ed
commit 24619a7bfc
16 changed files with 162 additions and 337 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/mlpack/core.hpp
+1
View File
@@ -81,6 +81,7 @@
#include <mlpack/core/math/range.hpp>
#include <mlpack/core/math/round.hpp>
#include <mlpack/core/math/shuffle_data.hpp>
#include <mlpack/core/math/ccov.hpp>
#include <mlpack/core/math/make_alias.hpp>
#include <mlpack/core/dists/discrete_distribution.hpp>
#include <mlpack/core/dists/gaussian_distribution.hpp>
src/mlpack/core/arma_extend/CMakeLists.txt
-5
View File
@@ -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
src/mlpack/core/arma_extend/arma_extend.hpp
-9
View File
@@ -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"
src/mlpack/core/arma_extend/fn_ccov.hpp
-34
View File
@@ -1,34 +0,0 @@
//! \addtogroup fn_ccov
//! @{
template<typename T1>
inline
const Op<T1, op_ccov>
ccov(const Base<typename T1::elem_type,T1>& 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<T1, op_ccov>(X.get_ref(), norm_type, 0);
}
template<typename T1, typename T2>
inline
const Glue<T1,T2,glue_ccov>
cov(const Base<typename T1::elem_type, T1>& A, const Base<typename T1::elem_type,T2>& 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<T1, T2, glue_ccov>(A.get_ref(), B.get_ref(), norm_type);
}
//! @}
src/mlpack/core/arma_extend/glue_ccov_meat.hpp
-144
View File
@@ -1,144 +0,0 @@
//! \addtogroup glue_cov
//! @{
template<typename eT>
inline
void
glue_ccov::direct_ccov(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& 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<N; ++i)
{
const eT A_tmp = A_ptr[i];
const eT B_tmp = B_ptr[i];
A_acc += A_tmp;
B_acc += B_tmp;
out_acc += A_tmp * B_tmp;
}
out_acc -= (A_acc * B_acc)/eT(N);
const eT norm_val = (norm_type == 0) ? ( (N > 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<typename T>
inline
void
glue_ccov::direct_ccov(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const Mat< std::complex<T> >& B, const uword norm_type)
{
arma_extra_debug_sigprint();
typedef typename std::complex<T> 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<N; ++i)
{
const eT A_tmp = A_ptr[i];
const eT B_tmp = B_ptr[i];
A_acc += A_tmp;
B_acc += B_tmp;
out_acc += std::conj(A_tmp) * B_tmp;
}
out_acc -= (std::conj(A_acc) * B_acc)/eT(N);
const eT norm_val = (norm_type == 0) ? ( (N > 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<typename T1, typename T2>
inline
void
glue_ccov::apply(Mat<typename T1::elem_type>& out, const Glue<T1,T2,glue_ccov>& X)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
const unwrap_check<T1> A_tmp(X.A, out);
const unwrap_check<T2> B_tmp(X.B, out);
const Mat<eT>& A = A_tmp.M;
const Mat<eT>& 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);
}
}
//! @}
src/mlpack/core/arma_extend/glue_ccov_proto.hpp
-15
View File
@@ -1,15 +0,0 @@
//! \addtogroup glue_ccov
//! @{
class glue_ccov
{
public:
template<typename eT> inline static void direct_ccov(Mat<eT>& out, const Mat<eT>& A, const Mat<eT>& B, const uword norm_type);
template<typename T> inline static void direct_ccov(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const Mat< std::complex<T> >& B, const uword norm_type);
template<typename T1, typename T2> inline static void apply(Mat<typename T1::elem_type>& out, const Glue<T1, T2, glue_ccov>& X);
};
//! @}
src/mlpack/core/arma_extend/op_ccov_meat.hpp
-97
View File
@@ -1,97 +0,0 @@
//! \addtogroup op_cov
//! @{
template<typename eT>
inline
void
op_ccov::direct_ccov(Mat<eT>& out, const Mat<eT>& 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<eT> acc = sum(A, 1);
out = A * trans(A);
out -= (acc * trans(acc)) / eT(N);
out /= norm_val;
}
}
template<typename T>
inline
void
op_ccov::direct_ccov(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& A, const uword norm_type)
{
arma_extra_debug_sigprint();
typedef typename std::complex<T> eT;
if(A.is_vec())
{
if(A.n_rows == 1)
{
const Mat<T> tmp_mat = var(trans(A), norm_type);
out.set_size(1,1);
out[0] = tmp_mat[0];
}
else
{
const Mat<T> 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<eT> acc = sum(A, 1);
out = A * trans(conj(A));
out -= (acc * trans(conj(acc))) / eT(N);
out /= norm_val;
}
}
template<typename T1>
inline
void
op_ccov::apply(Mat<typename T1::elem_type>& out, const Op<T1,op_ccov>& in)
{
arma_extra_debug_sigprint();
typedef typename T1::elem_type eT;
const unwrap_check<T1> tmp(in.m, out);
const Mat<eT>& A = tmp.M;
const uword norm_type = in.aux_uword_a;
op_ccov::direct_ccov(out, A, norm_type);
}
//! @}
src/mlpack/core/arma_extend/op_ccov_proto.hpp
-18
View File
@@ -1,18 +0,0 @@
//! \addtogroup op_cov
//! @{
class op_ccov
{
public:
template<typename eT> inline static void direct_ccov(Mat<eT>& out, const Mat<eT>& X, const uword norm_type);
template<typename T> inline static void direct_ccov(Mat< std::complex<T> >& out, const Mat< std::complex<T> >& X, const uword norm_type);
template<typename T1> inline static void apply(Mat<typename T1::elem_type>& out, const Op<T1,op_ccov>& in);
};
//! @}
src/mlpack/core/math/CMakeLists.txt
+2
View File
@@ -18,6 +18,8 @@ set(SOURCES
range_impl.hpp
round.hpp
shuffle_data.hpp
ccov.hpp
ccov_impl.hpp
)
# add directory name to sources
src/mlpack/core/math/ccov.hpp
+39
View File
@@ -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 <mlpack/prereqs.hpp>
namespace mlpack {
namespace math /** Miscellaneous math routines. */ {
template<typename eT>
inline
arma::Mat<eT>
ColumnCovariance(const arma::Mat<eT>& A, const size_t norm_type = 0);
template<typename T>
inline
arma::Mat< std::complex<T> >
ColumnCovariance(const arma::Mat< std::complex<T> >& A,
const size_t norm_type = 0);
} // namespace math
} // namespace mlpack
// Include implementation
#include "ccov_impl.hpp"
#endif // MLPACK_CORE_MATH_CCOV_HPP
src/mlpack/core/math/ccov_impl.hpp
+100
View File
@@ -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<typename eT>
inline arma::Mat<eT> ColumnCovariance(const arma::Mat<eT>& x,
const size_t normType)
{
if (normType > 1)
{
Log::Fatal << "ColumnCovariance(): norm_type must be 0 or 1!" << std::endl;
}
arma::Mat<eT> out;
if (x.n_elem > 0)
{
const arma::Mat<eT>& xAlias = (x.n_cols == 1) ?
arma::Mat<eT>(const_cast<eT*>(x.memptr()), x.n_cols, x.n_rows, false,
false) :
arma::Mat<eT>(const_cast<eT*>(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<eT> tmp = xAlias.each_col() - arma::mean(xAlias, 1);
out = tmp * tmp.t();
out /= normVal;
}
return out;
}
template<typename T>
inline arma::Mat<std::complex<T>> ColumnCovariance(
const arma::Mat<std::complex<T>>& x,
const size_t normType)
{
if (normType > 1)
{
Log::Fatal << "ColumnCovariance(): norm_type must be 0 or 1" << std::endl;
}
typedef typename std::complex<T> eT;
arma::Mat<eT> out;
if (x.is_vec())
{
if (x.n_rows == 1)
{
const arma::Mat<T> tmpMat = arma::var(arma::trans(x), normType);
out.set_size(1, 1);
out[0] = tmpMat[0];
}
else
{
const arma::Mat<T> 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<eT> 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
src/mlpack/core/math/lin_alg.cpp
+4 -4
View File
@@ -10,7 +10,7 @@
* http://www.opensource.org/licenses/BSD-3-Clause for more information.
*/
#include "lin_alg.hpp"
#include <mlpack/prereqs.hpp>
#include <mlpack/core.hpp>
#include <mlpack/core/math/random.hpp>
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);
src/mlpack/tests/distribution_test.cpp
+4 -4
View File
@@ -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);
src/mlpack/tests/gmm_test.cpp
+8 -4
View File
@@ -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);
src/mlpack/tests/hmm_test.cpp
+2 -1
View File
@@ -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);
src/mlpack/tests/lin_alg_test.cpp
+2 -2
View File
@@ -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++)
{