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 //=============================================================================================================
48 // EIGEN INCLUDES
49 //=============================================================================================================
50 
51 #include <Eigen/Core>
52 #include <Eigen/SparseCore>
53 #include <Eigen/SVD>
54 
55 //=============================================================================================================
56 // QT INCLUDES
57 //=============================================================================================================
58 
59 #include <QVariant>
60 
61 //=============================================================================================================
62 // FORWARD DECLARATIONS
63 //=============================================================================================================
64 
65 //=============================================================================================================
66 // DEFINE NAMESPACE UTILSLIB
67 //=============================================================================================================
68 
69 namespace UTILSLIB
70 {
71 
72 //=============================================================================================================
73 // UTILSLIB FORWARD DECLARATIONS
74 //=============================================================================================================
75 
76 //=============================================================================================================
84 {
85 public:
86  typedef std::pair<int,int> IdxIntValue;
88  //=========================================================================================================
92  virtual ~MNEMath()
93  { }
94 
95  //=========================================================================================================
104  static int gcd(int iA, int iB);
105 
106  //=========================================================================================================
120  static Eigen::VectorXd* combine_xyz(const Eigen::VectorXd& vec);
121 
122 // //=========================================================================================================
123 // /**
124 // * ### MNE toolbox root function ###: Definition of the mne_block_diag function - decoding part
125 // */
126 // static inline Eigen::MatrixXd extract_block_diag(MatrixXd& A, qint32 n);
127 
128  //=========================================================================================================
136  static double getConditionNumber(const Eigen::MatrixXd& A,
137  Eigen::VectorXd &s);
138 
139  //=========================================================================================================
147  static double getConditionSlope(const Eigen::MatrixXd& A,
148  Eigen::VectorXd &s);
149 
150  //=========================================================================================================
159  static void get_whitener(Eigen::MatrixXd& A,
160  bool pca,
161  QString ch_type,
162  Eigen::VectorXd& eig,
163  Eigen::MatrixXd& eigvec);
164 
165  //=========================================================================================================
175  static Eigen::VectorXi intersect(const Eigen::VectorXi &v1,
176  const Eigen::VectorXi &v2,
177  Eigen::VectorXi &idx_sel);
178 
179  //=========================================================================================================
188  static bool issparse(Eigen::VectorXd &v);
189 
190  //=========================================================================================================
201  static Eigen::MatrixXd legendre(qint32 n,
202  const Eigen::VectorXd &X,
203  QString normalize = QString("unnorm"));
204 
205  //=========================================================================================================
222  static Eigen::SparseMatrix<double>* make_block_diag(const Eigen::MatrixXd &A,
223  qint32 n);
224 
225  //=========================================================================================================
232  static int nchoose2(int n);
233 
234  //=========================================================================================================
245  static qint32 rank(const Eigen::MatrixXd& A,
246  double tol = 1e-8);
247 
248  //=========================================================================================================
264  static Eigen::MatrixXd rescale(const Eigen::MatrixXd &data,
265  const Eigen::RowVectorXf &times,
266  const QPair<float, float> &baseline,
267  QString mode);
268 
269  //=========================================================================================================
278  template<typename T>
279  static Eigen::VectorXi sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v,
280  bool desc = true);
281 
282  //=========================================================================================================
293  template<typename T>
294  static Eigen::VectorXi sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v_prime,
295  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &mat,
296  bool desc = true);
297 
298  //=========================================================================================================
307  template<typename T>
308  static std::vector<Eigen::Triplet<T> > sortrows(const std::vector<Eigen::Triplet<T> > &A,
309  qint32 column = 0);
310 
311  //=========================================================================================================
320  template<typename T>
321  static inline bool compareIdxValuePairBiggerThan(const std::pair<int,T>& lhs,
322  const std::pair<int,T>& rhs);
323 
324  //=========================================================================================================
333  template<typename T>
334  static inline bool compareIdxValuePairSmallerThan(const std::pair<int,T>& lhs,
335  const std::pair<int,T>& rhs);
336 
337  //=========================================================================================================
346  template<typename T>
347  static inline bool compareTripletFirstEntry(const Eigen::Triplet<T>& lhs,
348  const Eigen::Triplet<T> & rhs);
349 
350  //=========================================================================================================
359  template<typename T>
360  static inline bool compareTripletSecondEntry(const Eigen::Triplet<T>& lhs,
361  const Eigen::Triplet<T> & rhs);
362 
363  //=========================================================================================================
371  template<typename T>
372  static inline double log2(const T d);
373 
374  //=========================================================================================================
386  template<typename T>
387  static void histcounts(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matRawData,
388  bool bMakeSymmetrical,
389  int iClassAmount,
390  Eigen::VectorXd& vecResultClassLimits,
391  Eigen::VectorXi& vecResultFrequency,
392  double dGlobalMin = 0.0,
393  double dGlobalMax= 0.0);
394  template<typename T>
395  static void histcounts(const Eigen::Matrix<T, Eigen::Dynamic, 1>& matRawData,
396  bool bMakeSymmetrical,
397  int iClassAmount,
398  Eigen::VectorXd& vecResultClassLimits,
399  Eigen::VectorXi& vecResultFrequency,
400  double dGlobalMin = 0.0,
401  double dGlobalMax= 0.0);
402  template<typename T>
403  static void histcounts(const Eigen::Matrix<T, 1, Eigen::Dynamic>& matRawData,
404  bool bMakeSymmetrical,
405  int iClassAmount,
406  Eigen::VectorXd& vecResultClassLimits,
407  Eigen::VectorXi& vecResultFrequency,
408  double dGlobalMin = 0.0,
409  double dGlobalMax= 0.0);
410 
411  //=========================================================================================================
417  template<typename T>
418  static Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> pinv(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& a);
419 
420  //=========================================================================================================
432  static bool compareTransformation(const Eigen::MatrixX4f& mDevHeadT,
433  const Eigen::MatrixX4f& mDevHeadTNew,
434  const float& fThreshRot,
435  const float& fThreshTrans);
436 
437 };
438 
439 //=============================================================================================================
440 // INLINE & TEMPLATE DEFINITIONS
441 //=============================================================================================================
442 
443 template< typename T>
444 Eigen::VectorXi MNEMath::sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v,
445  bool desc)
446 {
447  std::vector< std::pair<int,T> > t_vecIdxValue;
448  Eigen::VectorXi idx(v.size());
449 
450  if(v.size() > 0)
451  {
452  //fill temporal vector
453  for(qint32 i = 0; i < v.size(); ++i)
454  t_vecIdxValue.push_back(std::pair<int,T>(i, v[i]));
455 
456  //sort temporal vector
457  if(desc)
458  std::sort(t_vecIdxValue.begin(), t_vecIdxValue.end(), MNEMath::compareIdxValuePairBiggerThan<T>);
459  else
460  std::sort(t_vecIdxValue.begin(), t_vecIdxValue.end(), MNEMath::compareIdxValuePairSmallerThan<T>);
461 
462  //store results
463  for(qint32 i = 0; i < v.size(); ++i)
464  {
465  idx[i] = t_vecIdxValue[i].first;
466  v[i] = t_vecIdxValue[i].second;
467  }
468  }
469 
470  return idx;
471 }
472 
473 //=============================================================================================================
474 
475 template<typename T>
476 Eigen::VectorXi MNEMath::sort(Eigen::Matrix<T, Eigen::Dynamic, 1> &v_prime,
477  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> &mat,
478  bool desc)
479 {
480  Eigen::VectorXi idx = MNEMath::sort<T>(v_prime, desc);
481 
482  if(v_prime.size() > 0)
483  {
484  //sort Matrix
485  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> newMat(mat.rows(), mat.cols());
486  for(qint32 i = 0; i < idx.size(); ++i)
487  newMat.col(i) = mat.col(idx[i]);
488  mat = newMat;
489  }
490 
491  return idx;
492 }
493 
494 //=============================================================================================================
495 
496 template<typename T>
497 std::vector<Eigen::Triplet<T> > MNEMath::sortrows(const std::vector<Eigen::Triplet<T> > &A,
498  qint32 column)
499 {
500  std::vector<Eigen::Triplet<T> > p_ASorted;
501 
502  for(quint32 i = 0; i < A.size(); ++i)
503  p_ASorted.push_back(A[i]);
504 
505  if(column == 0)
506  std::sort(p_ASorted.begin(), p_ASorted.end(), MNEMath::compareTripletFirstEntry<T>);
507  if(column == 1)
508  std::sort(p_ASorted.begin(), p_ASorted.end(), MNEMath::compareTripletSecondEntry<T>);
509 
510  return p_ASorted;
511 }
512 
513 //=============================================================================================================
514 
515 template<typename T>
516 inline bool MNEMath::compareIdxValuePairBiggerThan(const std::pair<int,T>& lhs,
517  const std::pair<int,T>& rhs)
518 {
519  return lhs.second > rhs.second;
520 }
521 
522 //=============================================================================================================
523 
524 template<typename T>
525 inline bool MNEMath::compareIdxValuePairSmallerThan( const std::pair<int,T>& lhs,
526  const std::pair<int,T>& rhs)
527 {
528  return lhs.second < rhs.second;
529 }
530 
531 //=============================================================================================================
532 
533 template<typename T>
534 inline bool MNEMath::compareTripletFirstEntry( const Eigen::Triplet<T>& lhs,
535  const Eigen::Triplet<T> & rhs)
536 {
537  return lhs.row() < rhs.row();
538 }
539 
540 //=============================================================================================================
541 
542 template<typename T>
543 inline bool MNEMath::compareTripletSecondEntry( const Eigen::Triplet<T>& lhs,
544  const Eigen::Triplet<T> & rhs)
545 {
546  return lhs.col() < rhs.col();
547 }
548 
549 //=============================================================================================================
550 
551 template<typename T>
552 inline double MNEMath::log2( const T d)
553 {
554  return log(d)/log(2);
555 }
556 
557 //=============================================================================================================
558 
559 template<typename T>
560 void MNEMath::histcounts(const Eigen::Matrix<T, Eigen::Dynamic, 1>& matRawData,
561  bool bMakeSymmetrical,
562  int iClassAmount,
563  Eigen::VectorXd& vecResultClassLimits,
564  Eigen::VectorXi& vecResultFrequency,
565  double dGlobalMin,
566  double dGlobalMax)
567 {
568  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrixName(matRawData.rows(),1);
569  matrixName.col(0)= matRawData;
570  MNEMath::histcounts(matrixName, bMakeSymmetrical, iClassAmount, vecResultClassLimits, vecResultFrequency, dGlobalMin, dGlobalMax);
571 }
572 
573 //=============================================================================================================
574 
575 template<typename T>
576 void MNEMath::histcounts(const Eigen::Matrix<T, 1, Eigen::Dynamic>& matRawData,
577  bool bMakeSymmetrical,
578  int iClassAmount,
579  Eigen::VectorXd& vecResultClassLimits,
580  Eigen::VectorXi& vecResultFrequency,
581  double dGlobalMin,
582  double dGlobalMax)
583 {
584  Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> matrixName(1,matRawData.cols());
585  matrixName.row(0)= matRawData;
586  MNEMath::histcounts(matrixName, bMakeSymmetrical, iClassAmount, vecResultClassLimits, vecResultFrequency, dGlobalMin, dGlobalMax);
587 }
588 
589 //=============================================================================================================
590 
591 template<typename T>
592 void MNEMath::histcounts(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& matRawData,
593  bool bMakeSymmetrical,
594  int iClassAmount,
595  Eigen::VectorXd& vecResultClassLimits,
596  Eigen::VectorXi& vecResultFrequency,
597  double dGlobalMin,
598  double dGlobalMax)
599 {
600  if(matRawData.rows() == 0 || matRawData.cols() == 0) {
601  return;
602  }
603 
604  vecResultClassLimits.resize(iClassAmount + 1);
605  vecResultFrequency.resize(iClassAmount);
606 
607  for (int count = 0; count < iClassAmount; ++count) //initialize the vector with zero values
608  {
609  vecResultFrequency(count) = 0;
610  }
611 
612  double desiredMin,
613  desiredMax;
614  double rawMin(0.0),
615  rawMax(0.0),
616  localMin(0.0),
617  localMax(0.0);
618 
619  rawMin = matRawData.minCoeff(); //finds the raw matrix minimum value
620  rawMax = matRawData.maxCoeff(); //finds the raw matrix maximum value
621 
622  if (bMakeSymmetrical == true) //in case the user wants the histogram to have symmetrical class ranges
623  {
624  if (std::fabs(rawMin) > rawMax) //in case the negative side is larger than the positive side
625  {
626  localMax = std::fabs(rawMin); //positive side is "stretched" to the exact length as negative side
627  localMin = rawMin;
628  }
629  else if (rawMax > std::fabs(rawMin)) //in case the positive side is larger than the negative side
630  {
631  localMin = -(rawMax); //negative side is "stretched" to the exact length as positive side
632  localMax = rawMax;
633  }
634  else //in case both sides are exactly the same
635  {
636  localMin = rawMin;
637  localMax = rawMax;
638  }
639  }
640  else //in case bMakeSymmetrical == false
641  {
642  localMin = rawMin;
643  localMax = rawMax;
644  }
645  //selects either local or global range (according to user preference and input)
646  if (dGlobalMin == 0.0 && dGlobalMax == 0.0) //if global range is NOT given by the user, use local ranges
647  {
648  desiredMin = localMin;
649  desiredMax = localMax;
650  vecResultClassLimits[0] = desiredMin; //replace default value with local minimum at position 0
651  vecResultClassLimits[iClassAmount] = desiredMax; //replace default value with local maximum at position n
652  }
653  else
654  {
655  desiredMin = dGlobalMin;
656  desiredMax = dGlobalMax;
657  vecResultClassLimits(0)= desiredMin; //replace default value with global minimum at position 0
658  vecResultClassLimits(iClassAmount)= desiredMax; //replace default value with global maximum at position n
659  }
660 
661  double range = (vecResultClassLimits(iClassAmount) - vecResultClassLimits(0)), //calculates the length from maximum positive value to zero
662  dynamicUpperClassLimit;
663 
664  for (int kr = 0; kr < iClassAmount; ++kr) //dynamically initialize the upper class limit values
665  {
666  dynamicUpperClassLimit = (vecResultClassLimits(0) + (kr*(range/iClassAmount))); //generic formula to determine the upper class limit with respect to range and number of class
667  vecResultClassLimits(kr) = dynamicUpperClassLimit; //places the appropriate upper class limit value to the right position in the QVector
668  }
669 
670  for (int ir = 0; ir < matRawData.rows(); ++ir) //iterates through all columns of the data matrix
671  {
672  for (int jr = 0; jr<matRawData.cols(); ++jr) //iterates through all rows of the data matrix
673  {
674  if(matRawData(ir,jr) != 0.0) {
675  for (int kr = 0; kr < iClassAmount; ++kr) //starts iteration from 1 to iClassAmount
676  {
677  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
678  {
679  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
680  {
681  vecResultFrequency(kr) = vecResultFrequency(kr) + 1 ; //if the value fits both arguments, the appropriate class frequency is increased by 1
682  }
683  }
684  else
685  {
686  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
687  {
688  vecResultFrequency(kr) = vecResultFrequency(kr) + 1 ; //if the value fits both arguments, the appropriate class frequency is increased by 1
689  }
690  }
691  }
692  }
693  }
694  }
695 }
696 
697 //=============================================================================================================
698 
699 template<typename T>
700 Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> MNEMath::pinv(const Eigen::Matrix<T,
701  Eigen::Dynamic,
702  Eigen::Dynamic>& a)
703 {
704  double epsilon = std::numeric_limits<double>::epsilon();
705  Eigen::JacobiSVD< Eigen::MatrixXd > svd(a ,Eigen::ComputeThinU | Eigen::ComputeThinV);
706  double tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().array().abs()(0);
707  return svd.matrixV() * (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(),0).matrix().asDiagonal() * svd.matrixU().adjoint();
708 }
709 
710 //=============================================================================================================
711 } // NAMESPACE
712 
713 #endif // MNEMATH_H
utils library export/import macros.
virtual ~MNEMath()
Definition: mnemath.h:92
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:592
static std::vector< Eigen::Triplet< T > > sortrows(const std::vector< Eigen::Triplet< T > > &A, qint32 column=0)
Definition: mnemath.h:497
static Eigen::VectorXi sort(Eigen::Matrix< T, Eigen::Dynamic, 1 > &v, bool desc=true)
Definition: mnemath.h:444
static bool compareIdxValuePairBiggerThan(const std::pair< int, T > &lhs, const std::pair< int, T > &rhs)
Definition: mnemath.h:516
static double log2(const T d)
Definition: mnemath.h:552
std::pair< int, int > IdxIntValue
Definition: mnemath.h:86
Math methods.
Definition: mnemath.h:83
static Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > pinv(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
Definition: mnemath.h:700
#define UTILSSHARED_EXPORT
Definition: utils_global.h:58
static bool compareTripletSecondEntry(const Eigen::Triplet< T > &lhs, const Eigen::Triplet< T > &rhs)
Definition: mnemath.h:543
static bool compareTripletFirstEntry(const Eigen::Triplet< T > &lhs, const Eigen::Triplet< T > &rhs)
Definition: mnemath.h:534
static bool compareIdxValuePairSmallerThan(const std::pair< int, T > &lhs, const std::pair< int, T > &rhs)
Definition: mnemath.h:525