MNE-CPP  0.1.9
A Framework for Electrophysiology
mnemath.h
Go to the documentation of this file.
1 //=============================================================================================================
36 #ifndef MNEMATH_H
37 #define MNEMATH_H
38 
39 //ToDo move this to the new MNE math library
40 
41 //=============================================================================================================
42 // INCLUDES
43 //=============================================================================================================
44 
45 #include "utils_global.h"
46 
47 #include <string>
48 #include <utility>
49 #include <vector>
50 
51 //=============================================================================================================
52 // EIGEN INCLUDES
53 //=============================================================================================================
54 
55 #include <Eigen/Core>
56 #include <Eigen/SparseCore>
57 #include <Eigen/SVD>
58 
59 //=============================================================================================================
60 // QT INCLUDES
61 //=============================================================================================================
62 
63 #include <QVariant>
64 
65 //=============================================================================================================
66 // FORWARD DECLARATIONS
67 //=============================================================================================================
68 
69 //=============================================================================================================
70 // DEFINE NAMESPACE UTILSLIB
71 //=============================================================================================================
72 
73 namespace UTILSLIB
74 {
75 
76 //=============================================================================================================
77 // UTILSLIB FORWARD DECLARATIONS
78 //=============================================================================================================
79 
80 //=============================================================================================================
88 {
89 public:
90  typedef std::pair<int,int> IdxIntValue;
92  //=========================================================================================================
96  virtual ~MNEMath()
97  { }
98 
99  //=========================================================================================================
108  static int gcd(int iA, int iB);
109 
110  //=========================================================================================================
124  static Eigen::VectorXd* combine_xyz(const Eigen::VectorXd& vec);
125 
126 // //=========================================================================================================
127 // /**
128 // * ### MNE toolbox root function ###: Definition of the mne_block_diag function - decoding part
129 // */
130 // static inline Eigen::MatrixXd extract_block_diag(MatrixXd& A, qint32 n);
131 
132  //=========================================================================================================
140  static double getConditionNumber(const Eigen::MatrixXd& A,
141  Eigen::VectorXd &s);
142 
143  //=========================================================================================================
151  static double getConditionSlope(const Eigen::MatrixXd& A,
152  Eigen::VectorXd &s);
153 
154  //=========================================================================================================
163  static void get_whitener(Eigen::MatrixXd& A,
164  bool pca,
165  QString ch_type,
166  Eigen::VectorXd& eig,
167  Eigen::MatrixXd& eigvec);
168 
169  //=========================================================================================================
178  static void get_whitener(Eigen::MatrixXd& A,
179  bool pca,
180  const std::string& ch_type,
181  Eigen::VectorXd& eig,
182  Eigen::MatrixXd& eigvec);
183 
184  //=========================================================================================================
194  static Eigen::VectorXi intersect(const Eigen::VectorXi &v1,
195  const Eigen::VectorXi &v2,
196  Eigen::VectorXi &idx_sel);
197 
198  //=========================================================================================================
207  static bool issparse(Eigen::VectorXd &v);
208 
209  //=========================================================================================================
220  static Eigen::MatrixXd legendre(qint32 n,
221  const Eigen::VectorXd &X,
222  QString normalize = QString("unnorm"));
223 
224  //=========================================================================================================
235  static Eigen::MatrixXd legendre(qint32 n,
236  const Eigen::VectorXd &X,
237  std::string normalize = "unnorm");
238 
239  //=========================================================================================================
256  static Eigen::SparseMatrix<double>* make_block_diag(const Eigen::MatrixXd &A,
257  qint32 n);
258 
259  //=========================================================================================================
266  static int nchoose2(int n);
267 
268  //=========================================================================================================
279  static qint32 rank(const Eigen::MatrixXd& A,
280  double tol = 1e-8);
281 
282  //=========================================================================================================
298  static Eigen::MatrixXd rescale(const Eigen::MatrixXd &data,
299  const Eigen::RowVectorXf &times,
300  const QPair<float, float> &baseline,
301  QString mode);
302 
303  //=========================================================================================================
319  static Eigen::MatrixXd rescale(const Eigen::MatrixXd& data,
320  const Eigen::RowVectorXf& times,
321  const std::pair<float, float>& baseline,
322  const std::string& mode);
323 
324  //=========================================================================================================
333  template<typename T>
334  static Eigen::VectorXi sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v,
335  bool desc = true);
336 
337  //=========================================================================================================
348  template<typename T>
349  static Eigen::VectorXi sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v_prime,
350  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &mat,
351  bool desc = true);
352 
353  //=========================================================================================================
362  template<typename T>
363  static std::vector<Eigen::Triplet<T> > sortrows(const std::vector<Eigen::Triplet<T> > &A,
364  qint32 column = 0);
365 
366  //=========================================================================================================
375  template<typename T>
376  static inline bool compareIdxValuePairBiggerThan(const std::pair<int,T>& lhs,
377  const std::pair<int,T>& rhs);
378 
379  //=========================================================================================================
388  template<typename T>
389  static inline bool compareIdxValuePairSmallerThan(const std::pair<int,T>& lhs,
390  const std::pair<int,T>& rhs);
391 
392  //=========================================================================================================
401  template<typename T>
402  static inline bool compareTripletFirstEntry(const Eigen::Triplet<T>& lhs,
403  const Eigen::Triplet<T> & rhs);
404 
405  //=========================================================================================================
414  template<typename T>
415  static inline bool compareTripletSecondEntry(const Eigen::Triplet<T>& lhs,
416  const Eigen::Triplet<T> & rhs);
417 
418  //=========================================================================================================
426  template<typename T>
427  static inline double log2(const T d);
428 
429  //=========================================================================================================
441  template<typename T>
442  static void histcounts(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matRawData,
443  bool bMakeSymmetrical,
444  int iClassAmount,
445  Eigen::VectorXd& vecResultClassLimits,
446  Eigen::VectorXi& vecResultFrequency,
447  double dGlobalMin = 0.0,
448  double dGlobalMax= 0.0);
449  template<typename T>
450  static void histcounts(const Eigen::Matrix<T, Eigen::Dynamic, 1>& matRawData,
451  bool bMakeSymmetrical,
452  int iClassAmount,
453  Eigen::VectorXd& vecResultClassLimits,
454  Eigen::VectorXi& vecResultFrequency,
455  double dGlobalMin = 0.0,
456  double dGlobalMax= 0.0);
457  template<typename T>
458  static void histcounts(const Eigen::Matrix<T, 1, Eigen::Dynamic>& matRawData,
459  bool bMakeSymmetrical,
460  int iClassAmount,
461  Eigen::VectorXd& vecResultClassLimits,
462  Eigen::VectorXi& vecResultFrequency,
463  double dGlobalMin = 0.0,
464  double dGlobalMax= 0.0);
465 
466  //=========================================================================================================
472  template<typename T>
473  static Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> pinv(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& a);
474 
475  //=========================================================================================================
487  static bool compareTransformation(const Eigen::MatrixX4f& mDevHeadT,
488  const Eigen::MatrixX4f& mDevHeadTNew,
489  const float& fThreshRot,
490  const float& fThreshTrans);
491 
492 };
493 
494 //=============================================================================================================
495 // INLINE & TEMPLATE DEFINITIONS
496 //=============================================================================================================
497 
498 template< typename T>
499 Eigen::VectorXi MNEMath::sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v,
500  bool desc)
501 {
502  std::vector< std::pair<int,T> > t_vecIdxValue;
503  Eigen::VectorXi idx(v.size());
504 
505  if(v.size() > 0)
506  {
507  //fill temporal vector
508  for(qint32 i = 0; i < v.size(); ++i)
509  t_vecIdxValue.push_back(std::pair<int,T>(i, v[i]));
510 
511  //sort temporal vector
512  if(desc)
513  std::sort(t_vecIdxValue.begin(), t_vecIdxValue.end(), MNEMath::compareIdxValuePairBiggerThan<T>);
514  else
515  std::sort(t_vecIdxValue.begin(), t_vecIdxValue.end(), MNEMath::compareIdxValuePairSmallerThan<T>);
516 
517  //store results
518  for(qint32 i = 0; i < v.size(); ++i)
519  {
520  idx[i] = t_vecIdxValue[i].first;
521  v[i] = t_vecIdxValue[i].second;
522  }
523  }
524 
525  return idx;
526 }
527 
528 //=============================================================================================================
529 
530 template<typename T>
531 Eigen::VectorXi MNEMath::sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v_prime,
532  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &mat,
533  bool desc)
534 {
535  Eigen::VectorXi idx = MNEMath::sort<T>(v_prime, desc);
536 
537  if(v_prime.size() > 0)
538  {
539  //sort Matrix
540  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> newMat(mat.rows(), mat.cols());
541  for(qint32 i = 0; i < idx.size(); ++i)
542  newMat.col(i) = mat.col(idx[i]);
543  mat = newMat;
544  }
545 
546  return idx;
547 }
548 
549 //=============================================================================================================
550 
551 template<typename T>
552 std::vector<Eigen::Triplet<T> > MNEMath::sortrows(const std::vector<Eigen::Triplet<T> > &A,
553  qint32 column)
554 {
555  std::vector<Eigen::Triplet<T> > p_ASorted;
556 
557  for(quint32 i = 0; i < A.size(); ++i)
558  p_ASorted.push_back(A[i]);
559 
560  if(column == 0)
561  std::sort(p_ASorted.begin(), p_ASorted.end(), MNEMath::compareTripletFirstEntry<T>);
562  if(column == 1)
563  std::sort(p_ASorted.begin(), p_ASorted.end(), MNEMath::compareTripletSecondEntry<T>);
564 
565  return p_ASorted;
566 }
567 
568 //=============================================================================================================
569 
570 template<typename T>
571 inline bool MNEMath::compareIdxValuePairBiggerThan(const std::pair<int,T>& lhs,
572  const std::pair<int,T>& rhs)
573 {
574  return lhs.second > rhs.second;
575 }
576 
577 //=============================================================================================================
578 
579 template<typename T>
580 inline bool MNEMath::compareIdxValuePairSmallerThan( const std::pair<int,T>& lhs,
581  const std::pair<int,T>& rhs)
582 {
583  return lhs.second < rhs.second;
584 }
585 
586 //=============================================================================================================
587 
588 template<typename T>
589 inline bool MNEMath::compareTripletFirstEntry( const Eigen::Triplet<T>& lhs,
590  const Eigen::Triplet<T> & rhs)
591 {
592  return lhs.row() < rhs.row();
593 }
594 
595 //=============================================================================================================
596 
597 template<typename T>
598 inline bool MNEMath::compareTripletSecondEntry( const Eigen::Triplet<T>& lhs,
599  const Eigen::Triplet<T> & rhs)
600 {
601  return lhs.col() < rhs.col();
602 }
603 
604 //=============================================================================================================
605 
606 template<typename T>
607 inline double MNEMath::log2( const T d)
608 {
609  return log(d)/log(2);
610 }
611 
612 //=============================================================================================================
613 
614 template<typename T>
615 void MNEMath::histcounts(const Eigen::Matrix<T, Eigen::Dynamic, 1>& matRawData,
616  bool bMakeSymmetrical,
617  int iClassAmount,
618  Eigen::VectorXd& vecResultClassLimits,
619  Eigen::VectorXi& vecResultFrequency,
620  double dGlobalMin,
621  double dGlobalMax)
622 {
623  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrixName(matRawData.rows(),1);
624  matrixName.col(0)= matRawData;
625  MNEMath::histcounts(matrixName, bMakeSymmetrical, iClassAmount, vecResultClassLimits, vecResultFrequency, dGlobalMin, dGlobalMax);
626 }
627 
628 //=============================================================================================================
629 
630 template<typename T>
631 void MNEMath::histcounts(const Eigen::Matrix<T, 1, Eigen::Dynamic>& matRawData,
632  bool bMakeSymmetrical,
633  int iClassAmount,
634  Eigen::VectorXd& vecResultClassLimits,
635  Eigen::VectorXi& vecResultFrequency,
636  double dGlobalMin,
637  double dGlobalMax)
638 {
639  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrixName(1,matRawData.cols());
640  matrixName.row(0)= matRawData;
641  MNEMath::histcounts(matrixName, bMakeSymmetrical, iClassAmount, vecResultClassLimits, vecResultFrequency, dGlobalMin, dGlobalMax);
642 }
643 
644 //=============================================================================================================
645 
646 template<typename T>
647 void MNEMath::histcounts(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matRawData,
648  bool bMakeSymmetrical,
649  int iClassAmount,
650  Eigen::VectorXd& vecResultClassLimits,
651  Eigen::VectorXi& vecResultFrequency,
652  double dGlobalMin,
653  double dGlobalMax)
654 {
655  if(matRawData.rows() == 0 || matRawData.cols() == 0) {
656  return;
657  }
658 
659  vecResultClassLimits.resize(iClassAmount + 1);
660  vecResultFrequency.resize(iClassAmount);
661 
662  for (int count = 0; count < iClassAmount; ++count) //initialize the vector with zero values
663  {
664  vecResultFrequency(count) = 0;
665  }
666 
667  double desiredMin,
668  desiredMax;
669  double rawMin(0.0),
670  rawMax(0.0),
671  localMin(0.0),
672  localMax(0.0);
673 
674  rawMin = matRawData.minCoeff(); //finds the raw matrix minimum value
675  rawMax = matRawData.maxCoeff(); //finds the raw matrix maximum value
676 
677  if (bMakeSymmetrical == true) //in case the user wants the histogram to have symmetrical class ranges
678  {
679  if (std::fabs(rawMin) > rawMax) //in case the negative side is larger than the positive side
680  {
681  localMax = std::fabs(rawMin); //positive side is "stretched" to the exact length as negative side
682  localMin = rawMin;
683  }
684  else if (rawMax > std::fabs(rawMin)) //in case the positive side is larger than the negative side
685  {
686  localMin = -(rawMax); //negative side is "stretched" to the exact length as positive side
687  localMax = rawMax;
688  }
689  else //in case both sides are exactly the same
690  {
691  localMin = rawMin;
692  localMax = rawMax;
693  }
694  }
695  else //in case bMakeSymmetrical == false
696  {
697  localMin = rawMin;
698  localMax = rawMax;
699  }
700  //selects either local or global range (according to user preference and input)
701  if (dGlobalMin == 0.0 && dGlobalMax == 0.0) //if global range is NOT given by the user, use local ranges
702  {
703  desiredMin = localMin;
704  desiredMax = localMax;
705  vecResultClassLimits[0] = desiredMin; //replace default value with local minimum at position 0
706  vecResultClassLimits[iClassAmount] = desiredMax; //replace default value with local maximum at position n
707  }
708  else
709  {
710  desiredMin = dGlobalMin;
711  desiredMax = dGlobalMax;
712  vecResultClassLimits(0)= desiredMin; //replace default value with global minimum at position 0
713  vecResultClassLimits(iClassAmount)= desiredMax; //replace default value with global maximum at position n
714  }
715 
716  double range = (vecResultClassLimits(iClassAmount) - vecResultClassLimits(0)), //calculates the length from maximum positive value to zero
717  dynamicUpperClassLimit;
718 
719  for (int kr = 0; kr < iClassAmount; ++kr) //dynamically initialize the upper class limit values
720  {
721  dynamicUpperClassLimit = (vecResultClassLimits(0) + (kr*(range/iClassAmount))); //generic formula to determine the upper class limit with respect to range and number of class
722  vecResultClassLimits(kr) = dynamicUpperClassLimit; //places the appropriate upper class limit value to the right position in the QVector
723  }
724 
725  for (int ir = 0; ir < matRawData.rows(); ++ir) //iterates through all columns of the data matrix
726  {
727  for (int jr = 0; jr<matRawData.cols(); ++jr) //iterates through all rows of the data matrix
728  {
729  if(matRawData(ir,jr) != 0.0) {
730  for (int kr = 0; kr < iClassAmount; ++kr) //starts iteration from 1 to iClassAmount
731  {
732  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
733  {
734  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
735  {
736  vecResultFrequency(kr) = vecResultFrequency(kr) + 1 ; //if the value fits both arguments, the appropriate class frequency is increased by 1
737  }
738  }
739  else
740  {
741  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
742  {
743  vecResultFrequency(kr) = vecResultFrequency(kr) + 1 ; //if the value fits both arguments, the appropriate class frequency is increased by 1
744  }
745  }
746  }
747  }
748  }
749  }
750 }
751 
752 //=============================================================================================================
753 
754 template<typename T>
755 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> MNEMath::pinv(const Eigen::Matrix<T,
756  Eigen::Dynamic,
757  Eigen::Dynamic>& a)
758 {
759  double epsilon = std::numeric_limits<double>::epsilon();
760  Eigen::JacobiSVD< Eigen::MatrixXd > svd(a ,Eigen::ComputeThinU | Eigen::ComputeThinV);
761  double tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().array().abs()(0);
762  return svd.matrixV() * (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(),0).matrix().asDiagonal() * svd.matrixU().adjoint();
763 }
764 
765 //=============================================================================================================
766 } // NAMESPACE
767 
768 #endif // MNEMATH_H
virtual ~MNEMath()
Definition: mnemath.h:96
utils library export/import macros.
Math methods.
Definition: mnemath.h:87
static Eigen::VectorXi sort(Eigen::Matrix< T, Eigen::Dynamic, 1 > &v, bool desc=true)
Definition: mnemath.h:499
static double log2(const T d)
Definition: mnemath.h:607
static bool compareIdxValuePairSmallerThan(const std::pair< int, T > &lhs, const std::pair< int, T > &rhs)
Definition: mnemath.h:580
static std::vector< Eigen::Triplet< T > > sortrows(const std::vector< Eigen::Triplet< T > > &A, qint32 column=0)
Definition: mnemath.h:552
#define UTILSSHARED_EXPORT
Definition: utils_global.h:58
std::pair< int, int > IdxIntValue
Definition: mnemath.h:90
static bool compareTripletSecondEntry(const Eigen::Triplet< T > &lhs, const Eigen::Triplet< T > &rhs)
Definition: mnemath.h:598
static Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > pinv(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
Definition: mnemath.h:755
static bool compareTripletFirstEntry(const Eigen::Triplet< T > &lhs, const Eigen::Triplet< T > &rhs)
Definition: mnemath.h:589
static bool compareIdxValuePairBiggerThan(const std::pair< int, T > &lhs, const std::pair< int, T > &rhs)
Definition: mnemath.h:571
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)
Definition: mnemath.h:647