EigenVectorAlgebra.h

Go to the documentation of this file.

#ifndef EIGENVECTORALGEBRA_H_
#define EIGENVECTORALGEBRA_H_
 
#include <Eigen/Core>
 
#include <boost/none.hpp>
#include <boost/any.hpp>
 
namespace muq {
  namespace Modeling {
    class EigenVectorAlgebra {
    public:
 
      EigenVectorAlgebra();
 
      virtual ~EigenVectorAlgebra();
 
 
      static bool IsEigenVector(std::type_info const& obj_type);
 
 
      static bool IsZero(boost::any const& obj);
 
 
      static unsigned int Size(boost::any const& vec);
 
 
      static double Norm(boost::any const& vec);
 
 
      static double InnerProduct(boost::any const& vec1, boost::any const& vec2);
 
 
      static boost::any OuterProduct(boost::any const& vec1, boost::any const& vec2);
 
 
      static boost::any AccessElement(boost::any const& obj, unsigned int const i);
 
 
      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 ScalarMultiply(boost::any const& in0, boost::any const& in1);
 
 
      static boost::any Apply(boost::any const& A, boost::any const& x);
 
 
      static boost::any ApplyInverse(boost::any const& A, boost::any const& x);
 
 
      static boost::any Zero(std::type_info const& type, unsigned int const size);
 
 
      static boost::any SquareRoot(boost::any const& obj);
 
 
      static double LogDeterminate(boost::any const& obj);
 
    private:
      
 
      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.size()).array()).all();
      }
      
 
      template<typename EigenType1, typename EigenType2>
    static inline double InnerProduct(boost::any const& vec1, boost::any const& vec2) {
    const EigenType1& x1 = boost::any_cast<EigenType1 const&>(vec1);
    const EigenType2& x2 = boost::any_cast<EigenType2 const&>(vec2);
    assert(x1.size()==x2.size());
    
    return x1.dot(x2);    
      }
 
 
      template<typename mattype, typename EigenType1, typename EigenType2>
    static inline boost::any OuterProduct(boost::any const& vec1, boost::any const& vec2) {
    const EigenType1& x1 = boost::any_cast<EigenType1 const&>(vec1);
    const EigenType2& x2 = boost::any_cast<EigenType2 const&>(vec2);
    
    return (mattype)(x1*x2.transpose());    
      }
 
 
      template<typename vectype>
    static inline boost::any AccessElement(boost::any const& vec, unsigned int const i) {
    // get a constant reference to the vector
    const vectype& vecref = boost::any_cast<const vectype&>(vec);
    
    // check the size
    assert(i<vecref.size());
    
    // return ith element
    return vecref(i);
      }
 
 
      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.size()==x1.size());
 
    return (type0)(x0+x1);
      }
 
 
      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.size()==x1.size());
 
    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>
    static inline 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.size()==vec.size());
    
    return (vectype)(mat.asDiagonal()*vec);
      }
 
 
      template<typename mattype, typename vectype>
    static inline 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.size()==vec.size());
    
    return (vectype)((1/mat.array()).matrix().asDiagonal()*vec);
      }
      
 
      template<typename type>
    inline static double LogDeterminate(boost::any const& obj) {
    const type& mat = boost::any_cast<type const&>(obj);
 
    return mat.array().log().sum();
      }
    };
  } // namespace Modeling
} // namespace muq
 
#endif