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