linalg.h

Go to the documentation of this file.

//=============================================================================================================
 
#ifndef LINALG_H
#define LINALG_H
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "math_global.h"
 
#include <utility>
#include <vector>
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Eigen/SVD>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QString>
 
//=============================================================================================================
// DEFINE NAMESPACE UTILSLIB
//=============================================================================================================
 
namespace UTILSLIB
{
 
//=============================================================================================================
class MATHSHARED_EXPORT Linalg
{
public:
    typedef std::pair<int,int> IdxIntValue;         
 
    //=========================================================================================================
    ~Linalg() = default;
 
    //=========================================================================================================
    static Eigen::VectorXd combine_xyz(const Eigen::VectorXd& vec);
 
    //=========================================================================================================
    static double getConditionNumber(const Eigen::MatrixXd& A,
                                     Eigen::VectorXd &s);
 
    //=========================================================================================================
    static double getConditionSlope(const Eigen::MatrixXd& A,
                                    Eigen::VectorXd &s);
 
    //=========================================================================================================
    static void get_whitener(Eigen::MatrixXd& A,
                             bool pca,
                             QString ch_type,
                             Eigen::VectorXd& eig,
                             Eigen::MatrixXd& eigvec);
 
    //=========================================================================================================
    static void get_whitener(Eigen::MatrixXd& A,
                             bool pca,
                             const std::string& ch_type,
                             Eigen::VectorXd& eig,
                             Eigen::MatrixXd& eigvec);
 
    //=========================================================================================================
    static Eigen::VectorXi intersect(const Eigen::VectorXi &v1,
                                     const Eigen::VectorXi &v2,
                                     Eigen::VectorXi &idx_sel);
 
    //=========================================================================================================
    static Eigen::SparseMatrix<double> make_block_diag(const Eigen::MatrixXd &A,
                                                       qint32 n);
 
    //=========================================================================================================
    template<typename T>
    static Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> pinv(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& a);
 
    //=========================================================================================================
    static qint32 rank(const Eigen::MatrixXd& A,
                       double tol = 1e-8);
 
    //=========================================================================================================
    template<typename T>
    static Eigen::VectorXi sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v,
                                bool desc = true);
 
    //=========================================================================================================
    template<typename T>
    static Eigen::VectorXi sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v_prime,
                                Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &mat,
                                bool desc = true);
 
    //=========================================================================================================
    template<typename T>
    static std::vector<Eigen::Triplet<T> > sortrows(const std::vector<Eigen::Triplet<T> > &A,
                                                    qint32 column = 0);
 
    //=========================================================================================================
    template<typename T>
    static inline bool compareIdxValuePairBiggerThan(const std::pair<int,T>& lhs,
                                                     const std::pair<int,T>& rhs);
 
    //=========================================================================================================
    template<typename T>
    static inline bool compareIdxValuePairSmallerThan(const std::pair<int,T>& lhs,
                                                      const std::pair<int,T>& rhs);
 
    //=========================================================================================================
    template<typename T>
    static inline bool compareTripletFirstEntry(const Eigen::Triplet<T>& lhs,
                                                const Eigen::Triplet<T> & rhs);
 
    //=========================================================================================================
    template<typename T>
    static inline bool compareTripletSecondEntry(const Eigen::Triplet<T>& lhs,
                                                 const Eigen::Triplet<T> & rhs);
};
 
//=============================================================================================================
// INLINE & TEMPLATE DEFINITIONS
//=============================================================================================================
 
template<typename T>
Eigen::VectorXi Linalg::sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v,
                             bool desc)
{
    std::vector< std::pair<int,T> > t_vecIdxValue;
    Eigen::VectorXi idx(v.size());
 
    if(v.size() > 0)
    {
        for(qint32 i = 0; i < v.size(); ++i)
            t_vecIdxValue.push_back(std::pair<int,T>(i, v[i]));
 
        if(desc)
            std::sort(t_vecIdxValue.begin(), t_vecIdxValue.end(), Linalg::compareIdxValuePairBiggerThan<T>);
        else
            std::sort(t_vecIdxValue.begin(), t_vecIdxValue.end(), Linalg::compareIdxValuePairSmallerThan<T>);
 
        for(qint32 i = 0; i < v.size(); ++i)
        {
            idx[i] = t_vecIdxValue[i].first;
            v[i] = t_vecIdxValue[i].second;
        }
    }
 
    return idx;
}
 
//=============================================================================================================
 
template<typename T>
Eigen::VectorXi Linalg::sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v_prime,
                             Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &mat,
                             bool desc)
{
    Eigen::VectorXi idx = Linalg::sort<T>(v_prime, desc);
 
    if(v_prime.size() > 0)
    {
        Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> newMat(mat.rows(), mat.cols());
        for(qint32 i = 0; i < idx.size(); ++i)
            newMat.col(i) = mat.col(idx[i]);
        mat = newMat;
    }
 
    return idx;
}
 
//=============================================================================================================
 
template<typename T>
std::vector<Eigen::Triplet<T> > Linalg::sortrows(const std::vector<Eigen::Triplet<T> > &A,
                                                  qint32 column)
{
    std::vector<Eigen::Triplet<T> > p_ASorted;
 
    for(quint32 i = 0; i < A.size(); ++i)
        p_ASorted.push_back(A[i]);
 
    if(column == 0)
        std::sort(p_ASorted.begin(), p_ASorted.end(), Linalg::compareTripletFirstEntry<T>);
    if(column == 1)
        std::sort(p_ASorted.begin(), p_ASorted.end(), Linalg::compareTripletSecondEntry<T>);
 
    return p_ASorted;
}
 
//=============================================================================================================
 
template<typename T>
inline bool Linalg::compareIdxValuePairBiggerThan(const std::pair<int,T>& lhs,
                                                  const std::pair<int,T>& rhs)
{
    return lhs.second > rhs.second;
}
 
//=============================================================================================================
 
template<typename T>
inline bool Linalg::compareIdxValuePairSmallerThan(const std::pair<int,T>& lhs,
                                                   const std::pair<int,T>& rhs)
{
    return lhs.second < rhs.second;
}
 
//=============================================================================================================
 
template<typename T>
inline bool Linalg::compareTripletFirstEntry(const Eigen::Triplet<T>& lhs,
                                             const Eigen::Triplet<T> & rhs)
{
    return lhs.row() < rhs.row();
}
 
//=============================================================================================================
 
template<typename T>
inline bool Linalg::compareTripletSecondEntry(const Eigen::Triplet<T>& lhs,
                                              const Eigen::Triplet<T> & rhs)
{
    return lhs.col() < rhs.col();
}
 
//=============================================================================================================
 
template<typename T>
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> Linalg::pinv(const Eigen::Matrix<T,
                                                               Eigen::Dynamic,
                                                               Eigen::Dynamic>& a)
{
    double epsilon = std::numeric_limits<double>::epsilon();
    Eigen::JacobiSVD< Eigen::MatrixXd > svd(a, Eigen::ComputeThinU | Eigen::ComputeThinV);
    double tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().array().abs()(0);
    return svd.matrixV() * (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0).matrix().asDiagonal() * svd.matrixU().adjoint();
}
 
//=============================================================================================================
} // NAMESPACE
 
#endif // LINALG_H