mnemath.h

Go to the documentation of this file.

//=============================================================================================================
#ifndef MNEMATH_H
#define MNEMATH_H
 
//ToDo move this to the new MNE math library
 
//=============================================================================================================
// INCLUDES
//=============================================================================================================
 
#include "utils_global.h"
 
#include <string>
#include <utility>
#include <vector>
 
//=============================================================================================================
// EIGEN INCLUDES
//=============================================================================================================
 
#include <Eigen/Core>
#include <Eigen/SparseCore>
#include <Eigen/SVD>
 
//=============================================================================================================
// QT INCLUDES
//=============================================================================================================
 
#include <QVariant>
 
//=============================================================================================================
// FORWARD DECLARATIONS
//=============================================================================================================
 
//=============================================================================================================
// DEFINE NAMESPACE UTILSLIB
//=============================================================================================================
 
namespace UTILSLIB
{
 
//=============================================================================================================
// UTILSLIB FORWARD DECLARATIONS
//=============================================================================================================
 
//=============================================================================================================
class UTILSSHARED_EXPORT MNEMath
{
public:
    typedef std::pair<int,int> IdxIntValue;         
    //=========================================================================================================
    virtual ~MNEMath()
    { }
 
    //=========================================================================================================
    static int gcd(int iA, int iB);
 
    //=========================================================================================================
    static Eigen::VectorXd* combine_xyz(const Eigen::VectorXd& vec);
 
//    //=========================================================================================================
//    /**
//    * ### MNE toolbox root function ###: Definition of the mne_block_diag function - decoding part
//    */
//    static inline Eigen::MatrixXd extract_block_diag(MatrixXd& A, qint32 n);
 
    //=========================================================================================================
    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 bool issparse(Eigen::VectorXd &v);
 
    //=========================================================================================================
    static Eigen::MatrixXd legendre(qint32 n,
                                    const Eigen::VectorXd &X,
                                    QString normalize = QString("unnorm"));
 
    //=========================================================================================================
    static Eigen::MatrixXd legendre(qint32 n,
                                    const Eigen::VectorXd &X,
                                    std::string normalize = "unnorm");
 
    //=========================================================================================================
    static Eigen::SparseMatrix<double>* make_block_diag(const Eigen::MatrixXd &A,
                                                        qint32 n);
 
    //=========================================================================================================
    static int nchoose2(int n);
 
    //=========================================================================================================
    static qint32 rank(const Eigen::MatrixXd& A,
                       double tol = 1e-8);
 
    //=========================================================================================================
    static Eigen::MatrixXd rescale(const Eigen::MatrixXd &data,
                                   const Eigen::RowVectorXf &times,
                                   const QPair<float, float> &baseline,
                                   QString mode);
 
    //=========================================================================================================
    static Eigen::MatrixXd rescale(const Eigen::MatrixXd& data,
                                   const Eigen::RowVectorXf& times,
                                   const std::pair<float, float>& baseline,
                                   const std::string& mode);
 
    //=========================================================================================================
    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);
 
    //=========================================================================================================
    template<typename T>
    static inline double log2(const T d);
 
    //=========================================================================================================
    template<typename T>
    static void histcounts(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matRawData,
                           bool bMakeSymmetrical,
                           int iClassAmount,
                           Eigen::VectorXd& vecResultClassLimits,
                           Eigen::VectorXi& vecResultFrequency,
                           double dGlobalMin = 0.0,
                           double dGlobalMax= 0.0);
    template<typename T>
    static void histcounts(const Eigen::Matrix<T, Eigen::Dynamic, 1>& matRawData,
                           bool bMakeSymmetrical,
                           int iClassAmount,
                           Eigen::VectorXd& vecResultClassLimits,
                           Eigen::VectorXi& vecResultFrequency,
                           double dGlobalMin = 0.0,
                           double dGlobalMax= 0.0);
    template<typename T>
    static void histcounts(const Eigen::Matrix<T, 1, Eigen::Dynamic>& matRawData,
                           bool bMakeSymmetrical,
                           int iClassAmount,
                           Eigen::VectorXd& vecResultClassLimits,
                           Eigen::VectorXi& vecResultFrequency,
                           double dGlobalMin = 0.0,
                           double dGlobalMax= 0.0);
 
    //=========================================================================================================
    template<typename T>
    static Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> pinv(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& a);
 
    //=========================================================================================================
    static bool compareTransformation(const Eigen::MatrixX4f& mDevHeadT,
                                      const Eigen::MatrixX4f& mDevHeadTNew,
                                      const float& fThreshRot,
                                      const float& fThreshTrans);
 
};
 
//=============================================================================================================
// INLINE & TEMPLATE DEFINITIONS
//=============================================================================================================
 
template< typename T>
Eigen::VectorXi MNEMath::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)
    {
        //fill temporal vector
        for(qint32 i = 0; i < v.size(); ++i)
            t_vecIdxValue.push_back(std::pair<int,T>(i, v[i]));
 
        //sort temporal vector
        if(desc)
            std::sort(t_vecIdxValue.begin(), t_vecIdxValue.end(), MNEMath::compareIdxValuePairBiggerThan<T>);
        else
            std::sort(t_vecIdxValue.begin(), t_vecIdxValue.end(), MNEMath::compareIdxValuePairSmallerThan<T>);
 
        //store results
        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 MNEMath::sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v_prime,
                              Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &mat,
                              bool desc)
{
    Eigen::VectorXi idx = MNEMath::sort<T>(v_prime, desc);
 
    if(v_prime.size() > 0)
    {
        //sort Matrix
        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> > MNEMath::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(), MNEMath::compareTripletFirstEntry<T>);
    if(column == 1)
        std::sort(p_ASorted.begin(), p_ASorted.end(), MNEMath::compareTripletSecondEntry<T>);
 
    return p_ASorted;
}
 
//=============================================================================================================
 
template<typename T>
inline bool MNEMath::compareIdxValuePairBiggerThan(const std::pair<int,T>& lhs,
                                                   const std::pair<int,T>& rhs)
{
    return lhs.second > rhs.second;
}
 
//=============================================================================================================
 
template<typename T>
inline bool MNEMath::compareIdxValuePairSmallerThan( const std::pair<int,T>& lhs,
                                                     const std::pair<int,T>& rhs)
{
    return lhs.second < rhs.second;
}
 
//=============================================================================================================
 
template<typename T>
inline bool MNEMath::compareTripletFirstEntry( const Eigen::Triplet<T>& lhs,
                                               const Eigen::Triplet<T> & rhs)
{
    return lhs.row() < rhs.row();
}
 
//=============================================================================================================
 
template<typename T>
inline bool MNEMath::compareTripletSecondEntry( const Eigen::Triplet<T>& lhs,
                                                const Eigen::Triplet<T> & rhs)
{
    return lhs.col() < rhs.col();
}
 
//=============================================================================================================
 
template<typename T>
inline double MNEMath::log2( const T d)
{
    return log(d)/log(2);
}
 
//=============================================================================================================
 
template<typename T>
void MNEMath::histcounts(const Eigen::Matrix<T, Eigen::Dynamic, 1>& matRawData,
                         bool bMakeSymmetrical,
                         int iClassAmount,
                         Eigen::VectorXd& vecResultClassLimits,
                         Eigen::VectorXi& vecResultFrequency,
                         double dGlobalMin,
                         double dGlobalMax)
{
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrixName(matRawData.rows(),1);
    matrixName.col(0)= matRawData;
    MNEMath::histcounts(matrixName, bMakeSymmetrical, iClassAmount, vecResultClassLimits, vecResultFrequency, dGlobalMin, dGlobalMax);
}
 
//=============================================================================================================
 
template<typename T>
void MNEMath::histcounts(const Eigen::Matrix<T, 1, Eigen::Dynamic>& matRawData,
                         bool bMakeSymmetrical,
                         int iClassAmount,
                         Eigen::VectorXd& vecResultClassLimits,
                         Eigen::VectorXi& vecResultFrequency,
                         double dGlobalMin,
                         double dGlobalMax)
{
    Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrixName(1,matRawData.cols());
    matrixName.row(0)= matRawData;
    MNEMath::histcounts(matrixName, bMakeSymmetrical, iClassAmount, vecResultClassLimits, vecResultFrequency, dGlobalMin, dGlobalMax);
}
 
//=============================================================================================================
 
template<typename T>
void MNEMath::histcounts(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matRawData,
                         bool bMakeSymmetrical,
                         int iClassAmount,
                         Eigen::VectorXd& vecResultClassLimits,
                         Eigen::VectorXi& vecResultFrequency,
                         double dGlobalMin,
                         double dGlobalMax)
{
   if(matRawData.rows() == 0 || matRawData.cols() == 0) {
       return;
   }
 
   vecResultClassLimits.resize(iClassAmount + 1);
   vecResultFrequency.resize(iClassAmount);
 
    for (int count = 0; count < iClassAmount; ++count) //initialize the vector with zero values
    {
        vecResultFrequency(count) = 0;
    }
 
    double desiredMin,
           desiredMax;
    double rawMin(0.0),
           rawMax(0.0),
           localMin(0.0),
           localMax(0.0);
 
    rawMin = matRawData.minCoeff();        //finds the raw matrix minimum value
    rawMax = matRawData.maxCoeff();        //finds the raw matrix maximum value
 
    if (bMakeSymmetrical == true)          //in case the user wants the histogram to have symmetrical class ranges
    {
        if (std::fabs(rawMin) > rawMax)          //in case the negative side is larger than the positive side
        {
            localMax = std::fabs(rawMin);        //positive side is "stretched" to the exact length as negative side
            localMin = rawMin;
        }
        else if (rawMax > std::fabs(rawMin))     //in case the positive side is larger than the negative side
        {
            localMin = -(rawMax);          //negative side is "stretched" to the exact length as positive side
            localMax = rawMax;
        }
        else                               //in case both sides are exactly the same
        {
            localMin = rawMin;
            localMax = rawMax;
        }
    }
    else                                   //in case bMakeSymmetrical == false
    {
        localMin = rawMin;
        localMax = rawMax;
    }
    //selects either local or global range (according to user preference and input)
    if (dGlobalMin == 0.0 && dGlobalMax == 0.0)               //if global range is NOT given by the user, use local ranges
    {
        desiredMin = localMin;
        desiredMax = localMax;
        vecResultClassLimits[0] = desiredMin;                 //replace default value with local minimum at position 0
        vecResultClassLimits[iClassAmount] = desiredMax;      //replace default value with local maximum at position n
    }
    else
    {
        desiredMin = dGlobalMin;
        desiredMax = dGlobalMax;
        vecResultClassLimits(0)= desiredMin;                 //replace default value with global minimum at position 0
        vecResultClassLimits(iClassAmount)= desiredMax;      //replace default value with global maximum at position n
    }
 
    double  range = (vecResultClassLimits(iClassAmount) - vecResultClassLimits(0)),      //calculates the length from maximum positive value to zero
            dynamicUpperClassLimit;
 
    for (int kr = 0; kr < iClassAmount; ++kr)                                            //dynamically initialize the upper class limit values
    {
        dynamicUpperClassLimit = (vecResultClassLimits(0) + (kr*(range/iClassAmount)));  //generic formula to determine the upper class limit with respect to range and number of class
        vecResultClassLimits(kr) = dynamicUpperClassLimit;                               //places the appropriate upper class limit value to the right position in the QVector
    }
 
    for (int ir = 0; ir < matRawData.rows(); ++ir)       //iterates through all columns of the data matrix
    {
        for (int jr = 0; jr<matRawData.cols(); ++jr)     //iterates through all rows of the data matrix
        {
            if(matRawData(ir,jr) != 0.0) {
                for (int kr = 0; kr < iClassAmount; ++kr)    //starts iteration from 1 to iClassAmount
                {
                    if (kr == iClassAmount-1)                //used for the final iteration; if the data value is exactly the same as the final upper class limit, it will be included in the histogram
                    {
                        if (matRawData(ir,jr) >= vecResultClassLimits(kr) && matRawData(ir,jr) <= vecResultClassLimits(kr + 1))    //compares value in the matrix with lower and upper limit of each class
                        {
                             vecResultFrequency(kr) = vecResultFrequency(kr) + 1 ;           //if the value fits both arguments, the appropriate class frequency is increased by 1
                        }
                    }
                    else
                    {
                        if (matRawData(ir,jr) >= vecResultClassLimits(kr) && matRawData(ir,jr) < vecResultClassLimits(kr + 1))    //compares value in the matrix with lower and upper limit of each class
                        {
                            vecResultFrequency(kr) = vecResultFrequency(kr) + 1 ;           //if the value fits both arguments, the appropriate class frequency is increased by 1
                        }
                    }
                }
            }
        }
    }
}
 
//=============================================================================================================
 
template<typename T>
Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> MNEMath::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 // MNEMATH_H