MNE-CPP  beta 1.0
mnemath.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //*************************************************************************************************************
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "mnemath.h"
42 
43 
44 //*************************************************************************************************************
45 //=============================================================================================================
46 // STL INCLUDES
47 //=============================================================================================================
48 
49 #include <iostream>
50 #include <algorithm> // std::sort
51 #include <vector> // std::vector
52 
53 //DEBUG fstream
54 //#include <fstream>
55 
56 
57 //*************************************************************************************************************
58 //=============================================================================================================
59 // QT INCLUDES
60 //=============================================================================================================
61 
62 #include <QFile>
63 #include <QDebug>
64 
65 
66 //*************************************************************************************************************
67 //=============================================================================================================
68 // USED NAMESPACES
69 //=============================================================================================================
70 
71 using namespace UTILSLIB;
72 
73 
74 //*************************************************************************************************************
75 //=============================================================================================================
76 // DEFINE MEMBER METHODS
77 //=============================================================================================================
78 
79 VectorXd* MNEMath::combine_xyz(const VectorXd& vec)
80 {
81  if (vec.size() % 3 != 0)
82  {
83  printf("Input must be a row or a column vector with 3N components");
84  return NULL;
85  }
86 
87  MatrixXd tmp = MatrixXd(vec.transpose());
88  SparseMatrix<double>* s = make_block_diag(tmp,3);
89 
90  SparseMatrix<double> sC = *s*s->transpose();
91  VectorXd* comb = new VectorXd(sC.rows());
92 
93  for(qint32 i = 0; i < sC.rows(); ++i)
94  (*comb)[i] = sC.coeff(i,i);
95 
96  delete s;
97  return comb;
98 }
99 
100 
101 //*************************************************************************************************************
102 
103 double MNEMath::getConditionNumber(const MatrixXd& A, VectorXd &s)
104 {
105  JacobiSVD<MatrixXd> svd(A);
106  s = svd.singularValues();
107 
108  double c = s.maxCoeff()/s.minCoeff();
109 
110  return c;
111 }
112 
113 
114 //*************************************************************************************************************
115 
116 double MNEMath::getConditionSlope(const MatrixXd& A, VectorXd &s)
117 {
118  JacobiSVD<MatrixXd> svd(A);
119  s = svd.singularValues();
120 
121  double c = s.maxCoeff()/s.mean();
122 
123  return c;
124 }
125 
126 
127 //*************************************************************************************************************
128 
129 void MNEMath::get_whitener(MatrixXd &A, bool pca, QString ch_type, VectorXd &eig, MatrixXd &eigvec)
130 {
131  // whitening operator
132  SelfAdjointEigenSolver<MatrixXd> t_eigenSolver(A);//Can be used because, covariance matrices are self-adjoint matrices.
133 
134  eig = t_eigenSolver.eigenvalues();
135  eigvec = t_eigenSolver.eigenvectors().transpose();
136 
137  MNEMath::sort<double>(eig, eigvec, false);
138  qint32 rnk = MNEMath::rank(A);
139 
140  for(qint32 i = 0; i < eig.size()-rnk; ++i)
141  eig(i) = 0;
142 
143  printf("Setting small %s eigenvalues to zero.\n", ch_type.toLatin1().constData());
144  if (!pca) // No PCA case.
145  printf("Not doing PCA for %s\n", ch_type.toLatin1().constData());
146  else
147  {
148  printf("Doing PCA for %s.",ch_type.toLatin1().constData());
149  // This line will reduce the actual number of variables in data
150  // and leadfield to the true rank.
151  eigvec = eigvec.block(eigvec.rows()-rnk, 0, rnk, eigvec.cols());
152  }
153 }
154 
155 
156 //*************************************************************************************************************
157 
158 VectorXi MNEMath::intersect(const VectorXi &v1, const VectorXi &v2, VectorXi &idx_sel)
159 {
160  std::vector<int> tmp;
161 
162  std::vector< std::pair<int,int> > t_vecIntIdxValue;
163 
164  //ToDo:Slow; map VectorXi to stl container
165  for(qint32 i = 0; i < v1.size(); ++i)
166  tmp.push_back(v1[i]);
167 
168  std::vector<int>::iterator it;
169  for(qint32 i = 0; i < v2.size(); ++i)
170  {
171  it = std::search(tmp.begin(), tmp.end(), &v2[i], &v2[i]+1);
172  if(it != tmp.end())
173  t_vecIntIdxValue.push_back(std::pair<int,int>(v2[i], it-tmp.begin()));//Index and int value are swapped // to sort using the idx
174  }
175 
176  std::sort(t_vecIntIdxValue.begin(), t_vecIntIdxValue.end(), MNEMath::compareIdxValuePairSmallerThan<int>);
177 
178  VectorXi p_res(t_vecIntIdxValue.size());
179  idx_sel = VectorXi(t_vecIntIdxValue.size());
180 
181  for(quint32 i = 0; i < t_vecIntIdxValue.size(); ++i)
182  {
183  p_res[i] = t_vecIntIdxValue[i].first;
184  idx_sel[i] = t_vecIntIdxValue[i].second;
185  }
186 
187  return p_res;
188 }
189 
190 
191 //*************************************************************************************************************
192 
193 // static inline MatrixXd extract_block_diag(MatrixXd& A, qint32 n)
194 // {
195 
196 
197 // //
198 // // Principal Investigators and Developers:
199 // // ** Richard M. Leahy, PhD, Signal & Image Processing Institute,
200 // // University of Southern California, Los Angeles, CA
201 // // ** John C. Mosher, PhD, Biophysics Group,
202 // // Los Alamos National Laboratory, Los Alamos, NM
203 // // ** Sylvain Baillet, PhD, Cognitive Neuroscience & Brain Imaging Laboratory,
204 // // CNRS, Hopital de la Salpetriere, Paris, France
205 // //
206 // // Copyright (c) 2005 BrainStorm by the University of Southern California
207 // // This software distributed under the terms of the GNU General Public License
208 // // as published by the Free Software Foundation. Further details on the GPL
209 // // license can be found at http://www.gnu.org/copyleft/gpl.html .
210 // //
211 // //FOR RESEARCH PURPOSES ONLY. THE SOFTWARE IS PROVIDED "AS IS," AND THE
212 // // UNIVERSITY OF SOUTHERN CALIFORNIA AND ITS COLLABORATORS DO NOT MAKE ANY
213 // // WARRANTY, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO WARRANTIES OF
214 // // MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE, NOR DO THEY ASSUME ANY
215 // // LIABILITY OR RESPONSIBILITY FOR THE USE OF THIS SOFTWARE.
216 // //
217 // // Author: John C. Mosher 1993 - 2004
218 // //
219 // //
220 // // Modifications for mne Matlab toolbox
221 // //
222 // // Matti Hamalainen
223 // // 2006
224 
225 
226 // [mA,na] = size(A); % matrix always has na columns
227 // % how many entries in the first column?
228 // bdn = na/n; % number of blocks
229 // ma = mA/bdn; % rows in first block
230 
231 // % blocks may themselves contain zero entries. Build indexing as above
232 // tmp = reshape([1:(ma*bdn)]',ma,bdn);
233 // i = zeros(ma*n,bdn);
234 // for iblock = 1:n,
235 // i((iblock-1)*ma+[1:ma],:) = tmp;
236 // end
237 
238 // i = i(:); % row indices foreach sparse bd
239 
240 
241 // j = [0:mA:(mA*(na-1))];
242 // j = j(ones(ma,1),:);
243 // j = j(:);
244 
245 // i = i + j;
246 
247 // bd = full(A(i)); % column vector
248 // bd = reshape(bd,ma,na); % full matrix
249 
250 // }
251 
252 
253 //*************************************************************************************************************
254 
255 bool MNEMath::issparse(VectorXd &v)
256 {
257  qDebug() << "ToDo: Figure out how to accelerate MNEMath::issparse(VectorXd &v).";
258 
259  qint32 c = 0;
260  qint32 n = v.rows();
261  qint32 t = n/2;
262 
263  for(qint32 i = 0; i < n; ++i)
264  {
265  if(v(i) == 0)
266  ++c;
267  if(c > t)
268  return true;
269  }
270 
271  return false;
272 }
273 
274 
275 //*************************************************************************************************************
276 
277 MatrixXd MNEMath::legendre(qint32 n, const VectorXd &X, QString normalize)
278 {
279  MatrixXd y;
280 
281  Q_UNUSED(y);
282 
283  Q_UNUSED(n);
284  Q_UNUSED(X);
285  Q_UNUSED(normalize);
286 
287  //ToDo
288 
289  return y;
290 }
291 
292 
293 //*************************************************************************************************************
294 
295 SparseMatrix<double>* MNEMath::make_block_diag(const MatrixXd &A, qint32 n)
296 {
297 
298  qint32 ma = A.rows();
299  qint32 na = A.cols();
300  float bdn = ((float)na)/n; // number of submatrices
301 
302 // std::cout << std::endl << "ma " << ma << " na " << na << " bdn " << bdn << std::endl;
303 
304  if(bdn - floor(bdn))
305  {
306  printf("Width of matrix must be even multiple of n\n");
307  return NULL;
308  }
309 
310  typedef Eigen::Triplet<double> T;
311  std::vector<T> tripletList;
312  tripletList.reserve(bdn*ma*n);
313 
314  qint32 current_col, current_row, i, r, c;
315  for(i = 0; i < bdn; ++i)
316  {
317  current_col = i * n;
318  current_row = i * ma;
319 
320  for(r = 0; r < ma; ++r)
321  for(c = 0; c < n; ++c)
322  tripletList.push_back(T(r+current_row, c+current_col, A(r, c+current_col)));
323  }
324 
325  SparseMatrix<double>* bd = new SparseMatrix<double>((int)floor((float)ma*bdn+0.5),na);
326 // SparseMatrix<double> p_Matrix(nrow, ncol);
327  bd->setFromTriplets(tripletList.begin(), tripletList.end());
328 
329  return bd;
330 }
331 
332 
333 //*************************************************************************************************************
334 
336 {
337 
338  //nchoosek(n, k) with k = 2, equals n*(n-1)*0.5
339 
340  int t_iNumOfCombination = (int)(n*(n-1)*0.5);
341 
342  return t_iNumOfCombination;
343 }
344 
345 
346 //*************************************************************************************************************
347 
348 qint32 MNEMath::rank(const MatrixXd& A, double tol)
349 {
350  JacobiSVD<MatrixXd> t_svdA(A);//U and V are not computed
351  VectorXd s = t_svdA.singularValues();
352  double t_dMax = s.maxCoeff();
353  t_dMax *= tol;
354  qint32 sum = 0;
355  for(qint32 i = 0; i < s.size(); ++i)
356  sum += s[i] > t_dMax ? 1 : 0;
357  return sum;
358 }
359 
360 
361 //*************************************************************************************************************
362 
363 MatrixXd MNEMath::rescale(const MatrixXd &data, const RowVectorXf &times, QPair<QVariant,QVariant> baseline, QString mode)
364 {
365  MatrixXd data_out = data;
366  QStringList valid_modes;
367  valid_modes << "logratio" << "ratio" << "zscore" << "mean" << "percent";
368  if(!valid_modes.contains(mode))
369  {
370  qWarning() << "\tWarning: mode should be any of : " << valid_modes;
371  return data_out;
372  }
373  printf("\tApplying baseline correction ... (mode: %s)\n", mode.toLatin1().constData());
374 
375  qint32 imin, imax;
376  float bmin, bmax;
377 
378  if(!baseline.first.isValid())
379  imin = 0;
380  else
381  {
382  bmin = baseline.first.toFloat();
383  for(qint32 i = 0; i < times.size(); ++i)
384  {
385  if(times[i] >= bmin)
386  {
387  imin = i;
388  break;
389  }
390  }
391  }
392  if (!baseline.second.isValid())
393  imax = times.size();
394  else
395  {
396  bmax = baseline.second.toFloat();
397  for(qint32 i = times.size()-1; i >= 0; --i)
398  {
399  if(times[i] <= bmax)
400  {
401  imax = i+1;
402  break;
403  }
404  }
405  }
406 
407  VectorXd mean = data_out.block(0, imin,data_out.rows(),imax-imin).rowwise().mean();
408  if(mode.compare("mean") == 0)
409  {
410  data_out -= mean.rowwise().replicate(data.cols());
411  }
412  else if(mode.compare("logratio") == 0)
413  {
414  for(qint32 i = 0; i < data_out.rows(); ++i)
415  for(qint32 j = 0; j < data_out.cols(); ++j)
416  data_out(i,j) = log10(data_out(i,j)/mean[i]); // a value of 1 means 10 times bigger
417  }
418  else if(mode.compare("ratio") == 0)
419  {
420  data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
421  }
422  else if(mode.compare("zscore") == 0)
423  {
424  MatrixXd std_mat = data.block(0, imin, data.rows(), imax-imin) - mean.rowwise().replicate(imax-imin);
425  std_mat = std_mat.cwiseProduct(std_mat);
426  VectorXd std_v = std_mat.rowwise().mean();
427  for(qint32 i = 0; i < std_v.size(); ++i)
428  std_v[i] = sqrt(std_v[i] / (float)(imax-imin));
429 
430  data_out -= mean.rowwise().replicate(data_out.cols());
431  data_out = data_out.cwiseQuotient(std_v.rowwise().replicate(data_out.cols()));
432  }
433  else if(mode.compare("percent") == 0)
434  {
435  data_out -= mean.rowwise().replicate(data_out.cols());
436  data_out = data_out.cwiseQuotient(mean.rowwise().replicate(data_out.cols()));
437  }
438 
439  return data_out;
440 }
static void get_whitener(MatrixXd &A, bool pca, QString ch_type, VectorXd &eig, MatrixXd &eigvec)
Definition: mnemath.cpp:129
static VectorXd * combine_xyz(const VectorXd &vec)
Definition: mnemath.cpp:79
static bool issparse(VectorXd &v)
Definition: mnemath.cpp:255
static SparseMatrix< double > * make_block_diag(const MatrixXd &A, qint32 n)
Definition: mnemath.cpp:295
static double getConditionNumber(const MatrixXd &A, VectorXd &s)
Definition: mnemath.cpp:103
static VectorXi intersect(const VectorXi &v1, const VectorXi &v2, VectorXi &idx_sel)
Definition: mnemath.cpp:158
static qint32 rank(const MatrixXd &A, double tol=1e-8)
Definition: mnemath.cpp:348
MNEMath class declaration.
static MatrixXd rescale(const MatrixXd &data, const RowVectorXf &times, QPair< QVariant, QVariant > baseline, QString mode)
Definition: mnemath.cpp:363
static MatrixXd legendre(qint32 n, const VectorXd &X, QString normalize=QString("unnorm"))
Definition: mnemath.cpp:277
static int nchoose2(int n)
Definition: mnemath.cpp:335
static double getConditionSlope(const MatrixXd &A, VectorXd &s)
Definition: mnemath.cpp:116