EigenMatrixAlgebra.h

Go to the documentation of this file.

#ifndef EIGENMATRIXALGEBRA_H_
#define EIGENMATRIXALGEBRA_H_
 
#include <iostream>
 
#include <Eigen/Core>
#include <Eigen/Dense>
#include <unsupported/Eigen/MatrixFunctions>
 
#include <boost/none.hpp>
#include <boost/any.hpp>
 
namespace muq {
  namespace Modeling {
    class EigenMatrixAlgebra {
    public:
 
      EigenMatrixAlgebra();
 
      ~EigenMatrixAlgebra();
 
 
      static bool IsEigenMatrix(std::type_info const& obj_type);
 
 
      static bool IsZero(boost::any const& obj);
 
 
      static double Norm(boost::any const& mat);
 
 
      static unsigned int Size(boost::any const& mat, int const dim);
 
 
      static boost::any AccessElement(boost::any const& mat, unsigned int const i, unsigned int const j);
 
 
      static boost::any Identity(std::type_info const& type, unsigned int const rows, unsigned int const cols);
 
 
      static boost::any Add(boost::any const& in0, boost::any const& in1);
 
 
      static boost::any Subtract(boost::any const& in0, boost::any const& in1);
 
 
      static boost::any Multiply(boost::any const& in0, boost::any const& in1);
 
 
      static boost::any ScalarMultiply(boost::any const& in0, boost::any const& in1);
 
 
      static boost::any Zero(std::type_info const& type, unsigned int const rows, unsigned int const cols);
 
 
      static boost::any ApplyInverse(boost::any const& A, boost::any const& x);
 
 
      static boost::any Apply(boost::any const& A, boost::any const& x);
 
 
      static boost::any SquareRoot(boost::any const& obj);
 
 
      static double LogDeterminate(boost::any const& obj);
 
    private:
 
 
      template<typename type0, typename type1>
    static inline boost::any Add(boost::any const& in0, boost::any const& in1) {
    const type0& x0 = boost::any_cast<type0 const&>(in0);
    const type1& x1 = boost::any_cast<type1 const&>(in1);
    assert(x0.rows()==x1.rows());
    assert(x0.cols()==x1.cols());
 
    return (type0)(x0+x1);
      }
 
 
      template<typename mattype, typename vectype>
    static inline boost::any AddVector(boost::any const& in0, boost::any const& in1) {
    mattype mat = boost::any_cast<mattype>(in0);
    const vectype& vec = boost::any_cast<vectype const&>(in1);
    assert(mat.rows()==vec.size());
    assert(mat.cols()==vec.size());
    
    //Eigen::VectorXd& diag = mat.diagonal();
    mat.diagonal() += vec;
 
    return mat;
 
    //Eigen:: result = mat+vec.asDiagonal();
    //return result;
      }
 
 
      template<typename type>
    static inline unsigned int Size(boost::any const& mat, int const dim) {
    const type& eig = boost::any_cast<type const&>(mat);
    switch( dim ) {
    case 0:
      return eig.rows();
    case 1:
      return eig.cols();
    default:
      return eig.size();
    }
      }
 
 
      template<typename EigenType>
    static inline bool IsZero(boost::any const& obj) {
    const EigenType& v = boost::any_cast<EigenType const&>(obj);
 
    return (v.array()==EigenType::Zero(v.rows(), v.cols()).array()).all();
      }
 
 
      template<typename mattype>
    static inline boost::any AccessElement(boost::any const& mat, unsigned int const i, unsigned int const j) {
    // get a constant reference to the matrix
    const mattype& matref = boost::any_cast<const mattype&>(mat);
    
    // check the size
    assert(i<matref.rows());
    assert(j<matref.cols());
    
    // return ith element
    return matref(i,j);
      }
      
 
      template<typename type0, typename type1>
    static inline boost::any Subtract(boost::any const& in0, boost::any const& in1) {
    const type0& x0 = boost::any_cast<type0 const&>(in0);
    const type1& x1 = boost::any_cast<type1 const&>(in1);
    assert(x0.rows()==x1.rows());
    assert(x0.cols()==x1.cols());
 
    return (type0)(x0-x1);
      }
 
 
      template<typename type0, typename type1>
    static inline boost::any Multiply(boost::any const& in0, boost::any const& in1) {
    const type0& x0 = boost::any_cast<type0 const&>(in0);
    const type1& x1 = boost::any_cast<type1 const&>(in1);
    assert(x0.cols()==x1.rows());
 
    return (type0)(x0*x1);
      }
 
 
      template<typename type0, typename type1> 
    static inline boost::any ScalarMultiply(boost::any const& in0, boost::any const& in1) {
    const type0& x0 = boost::any_cast<type0 const&>(in0);
    const type1& x1 = boost::any_cast<type1 const&>(in1);
    
    return (type1)(x0 * x1);
      }
 
 
      template<typename mattype, typename vectype>
    inline static boost::any ApplyInverse(boost::any const& A, boost::any const& x) {
    const mattype& mat = boost::any_cast<mattype const&>(A);
    const vectype& vec = boost::any_cast<vectype const&>(x);
    assert(mat.rows()==vec.size());
    assert(mat.cols()==vec.size());
 
    // solve the system
    const vectype soln = mat.colPivHouseholderQr().solve(vec);
 
    return soln;
      }
 
 
      template<typename mattype, typename vectype>
    inline static boost::any ApplyCholeskyInverse(boost::any const& A, boost::any const& x) {
    const Eigen::LLT<mattype>& chol = boost::any_cast<Eigen::LLT<mattype> const&>(A);
    const mattype& L = chol.matrixL();
  
    const vectype& vec = boost::any_cast<vectype const&>(x);
    assert(L.rows()==vec.size());
    assert(L.cols()==vec.size());
    
    // solve the system
    vectype soln = L.template triangularView<Eigen::Lower>().solve(vec);
    L.template triangularView<Eigen::Lower>().transpose().solveInPlace(soln);
 
    return soln;
      }
 
 
      template<typename mattype, typename vectype>
    inline static boost::any Apply(boost::any const& A, boost::any const& x) {
    const mattype& mat = boost::any_cast<mattype const&>(A);
    const vectype& vec = boost::any_cast<vectype const&>(x);
    assert(mat.cols()==vec.size());
 
    return (vectype)(mat*vec);
      }
 
 
      template<typename mattype, typename vectype>
    inline static boost::any ApplyCholesky(boost::any const& A, boost::any const& x) {
    const Eigen::LLT<mattype>& chol = boost::any_cast<Eigen::LLT<mattype> const&>(A);
    const mattype& mat = chol.matrixL();
    const vectype& vec = boost::any_cast<vectype const&>(x);
    assert(mat.cols()==vec.size());
 
    vectype soln = mat.template triangularView<Eigen::Lower>().transpose()*vec;
    soln = mat.template triangularView<Eigen::Lower>()*soln;
 
    return soln;
      }
 
 
      template<typename mattype>
    inline static boost::any SquareRoot(boost::any const& obj) {
    const Eigen::LLT<mattype>& mat = boost::any_cast<Eigen::LLT<mattype> const&>(obj);
 
    return (mattype)(mat.matrixL());
      }
 
 
      template<typename mattype>
    inline static double LogDeterminate(boost::any const& obj) {
    const Eigen::LLT<mattype>& chol = boost::any_cast<Eigen::LLT<mattype> const&>(obj);
    const mattype& mat = chol.matrixL();
 
    return 2.0*mat.diagonal().array().log().sum();
      }
    };
  } // namespace Modeling
} // namespace muq
 
#endif