MNE-CPP  0.1.9
A Framework for Electrophysiology
mne_cov_matrix.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "mne_cov_matrix.h"
42 #include "mne_sss_data.h"
43 #include "mne_proj_item.h"
44 #include "mne_proj_op.h"
45 
46 #include <Eigen/Core>
47 #include <Eigen/Eigenvalues>
49 
50 //=============================================================================================================
51 // USED NAMESPACES
52 //=============================================================================================================
53 
54 using namespace Eigen;
55 using namespace FIFFLIB;
56 using namespace MNELIB;
57 
58 #ifndef FAIL
59 #define FAIL -1
60 #endif
61 
62 #ifndef OK
63 #define OK 0
64 #endif
65 
66 #define MALLOC_30(x,t) (t *)malloc((x)*sizeof(t))
67 #define REALLOC_30(x,y,t) (t *)((x == NULL) ? malloc((y)*sizeof(t)) : realloc((x),(y)*sizeof(t)))
68 
69 #define FREE_30(x) if ((char *)(x) != NULL) free((char *)(x))
70 
71 #define FREE_CMATRIX_30(m) mne_free_cmatrix_30((m))
72 #define FREE_DCMATRIX_30(m) mne_free_dcmatrix_30((m))
73 
74 void mne_free_cmatrix_30 (float **m)
75 {
76  if (m) {
77  FREE_30(*m);
78  FREE_30(m);
79  }
80 }
81 
82 void mne_free_dcmatrix_30 (double **m)
83 
84 {
85  if (m) {
86  FREE_30(*m);
87  FREE_30(m);
88  }
89 }
90 
91 #define ALLOC_CMATRIX_30(x,y) mne_cmatrix_30((x),(y))
92 
93 #define ALLOC_DCMATRIX_30(x,y) mne_dmatrix_30((x),(y))
94 
95 static void matrix_error_30(int kind, int nr, int nc)
96 
97 {
98  if (kind == 1)
99  printf("Failed to allocate memory pointers for a %d x %d matrix\n",nr,nc);
100  else if (kind == 2)
101  printf("Failed to allocate memory for a %d x %d matrix\n",nr,nc);
102  else
103  printf("Allocation error for a %d x %d matrix\n",nr,nc);
104  if (sizeof(void *) == 4) {
105  printf("This is probably because you seem to be using a computer with 32-bit architecture.\n");
106  printf("Please consider moving to a 64-bit platform.");
107  }
108  printf("Cannot continue. Sorry.\n");
109  exit(1);
110 }
111 
112 float **mne_cmatrix_30(int nr,int nc)
113 
114 {
115  int i;
116  float **m;
117  float *whole;
118 
119  m = MALLOC_30(nr,float *);
120  if (!m)
121  matrix_error_30(1,nr,nc);
122 
123  whole = MALLOC_30(nr*nc,float);
124  if (!whole)
125  matrix_error_30(2,nr,nc);
126 
127  for(i=0;i<nr;i++)
128  m[i] = whole + i*nc;
129 
130  return m;
131 }
132 
133 double **mne_dmatrix_30(int nr, int nc)
134 
135 {
136  int i;
137  double **m;
138  double *whole;
139 
140  m = MALLOC_30(nr,double *);
141  if (!m)
142  matrix_error_30(1,nr,nc);
143 
144  whole = MALLOC_30(nr*nc,double);
145  if (!whole)
146  matrix_error_30(2,nr,nc);
147 
148  for(i=0;i<nr;i++)
149  m[i] = whole + i*nc;
150 
151  return m;
152 }
153 
154 //============================= mne_decompose.c =============================
155 
156 int mne_decompose_eigen (double *mat,
157  double *lambda,
158  float **vectors, /* Eigenvectors fit into floats easily */
159  int dim)
160 /*
161  * Compute the eigenvalue decomposition of
162  * a symmetric matrix using the LAPACK routines
163  *
164  * 'mat' contains the lower triangle of the matrix
165  */
166 {
167  int np = dim*(dim+1)/2;
168  double *w = MALLOC_30(dim,double);
169  double *z = MALLOC_30(dim*dim,double);
170  double *work = MALLOC_30(3*dim,double);
171  double *dmat = MALLOC_30(np,double);
172  float *vecp = vectors[0];
173 
174  int info,k;
175  int maxi;
176  double scale;
177 
178 // maxi = idamax(&np,mat,&one);
179 // idamax workaround begin
180  maxi = 0;
181  for(int i = 0; i < np; ++i)
182  if (std::fabs(mat[i]) > std::fabs(mat[maxi]))
183  maxi = i;
184 // idamax workaround end
185 
186  scale = 1.0/mat[maxi];//scale = 1.0/mat[maxi-1];
187 
188  for (k = 0; k < np; k++)
189  dmat[k] = mat[k]*scale;
190 // dspev(compz,uplo,&dim,dmat,w,z,&dim,work,&info);
191 
192 // dspev workaround begin
193  MatrixXd dmat_tmp = MatrixXd::Zero(dim,dim);
194  int idx = 0;
195  for (int i = 0; i < dim; ++i) {
196  for(int j = 0; j <= i; ++j) {
197  dmat_tmp(i,j) = dmat[idx];
198  dmat_tmp(j,i) = dmat[idx];
199  ++idx;
200  }
201  }
202  SelfAdjointEigenSolver<MatrixXd> es;
203  es.compute(dmat_tmp);
204  for ( int i = 0; i < dim; ++i )
205  w[i] = es.eigenvalues()[i];
206 
207  idx = 0;
208  for ( int j = 0; j < dim; ++j ) {
209  for( int i = 0; i < dim; ++i ) {
210  z[idx] = es.eigenvectors()(i,j);// Column Major
211  ++idx;
212  }
213  }
214 // dspev workaround end
215 
216  info = 0;
217 
218  qDebug() << "!!!DEBUG ToDo: dspev(compz,uplo,&dim,dmat,w,z,&dim,work,&info);";
219 
220  FREE_30(work);
221  if (info != 0)
222  printf("Eigenvalue decomposition failed (LAPACK info = %d)",info);
223  else {
224  scale = 1.0/scale;
225  for (k = 0; k < dim; k++)
226  lambda[k] = scale*w[k];
227  for (k = 0; k < dim*dim; k++)
228  vecp[k] = z[k];
229  }
230  FREE_30(w);
231  FREE_30(z);
232  FREE_30(dmat);
233  if (info == 0)
234  return 0;
235  else
236  return -1;
237 }
238 
239 double **mne_dmatt_dmat_mult2 (double **m1,double **m2, int d1,int d2,int d3)
240 /* Matrix multiplication
241  * result(d1 x d3) = m1(d2 x d1)^T * m2(d2 x d3) */
242 
243 {
244  double **result = ALLOC_DCMATRIX_30(d1,d3);
245 #ifdef BLAS
246  char *transa = "N";
247  char *transb = "T";
248  double zero = 0.0;
249  double one = 1.0;
250 
251  dgemm (transa,transb,&d3,&d1,&d2,
252  &one,m2[0],&d3,m1[0],&d1,&zero,result[0],&d3);
253 
254  return result;
255 #else
256  int j,k,p;
257  double sum;
258 
259  for (j = 0; j < d1; j++)
260  for (k = 0; k < d3; k++) {
261  sum = 0.0;
262  for (p = 0; p < d2; p++)
263  sum = sum + m1[p][j]*m2[p][k];
264  result[j][k] = sum;
265  }
266  return result;
267 #endif
268 }
269 
270 //=============================================================================================================
271 // DEFINE MEMBER METHODS
272 //=============================================================================================================
273 
274 MneCovMatrix::MneCovMatrix(int p_kind,
275  int p_ncov,
276  const QStringList& p_names,
277  double *p_cov,
278  double *p_cov_diag,
279  FiffSparseMatrix* p_cov_sparse)
280 :kind(p_kind)
281 ,ncov(p_ncov)
282 ,nfree(1)
283 ,nproj(0)
284 ,nzero(0)
285 ,names(p_names)
286 ,cov(p_cov)
287 ,cov_diag(p_cov_diag)
288 ,cov_sparse(p_cov_sparse)
289 ,lambda(NULL)
290 ,inv_lambda(NULL)
291 ,eigen(NULL)
292 ,chol(NULL)
293 ,proj(NULL)
294 ,sss(NULL)
295 ,ch_class(NULL)
296 ,nbad(0)
297 {
298 }
299 
300 //=============================================================================================================
301 
303 {
304  FREE_30(cov);
305  FREE_30(cov_diag);
306  if(cov_sparse)
307  delete cov_sparse;
308  names.clear();
309  FREE_CMATRIX_30(eigen);
310  FREE_30(lambda);
311  FREE_30(inv_lambda);
312  FREE_30(chol);
313  FREE_30(ch_class);
314  if(proj)
315  delete proj;
316  if (sss)
317  delete sss;
318  bads.clear();
319 }
320 
321 //=============================================================================================================
322 
323 MneCovMatrix *MneCovMatrix::mne_dup_cov(MneCovMatrix *c)
324 {
325  double *vals;
326  int nval;
327  int k;
328  MneCovMatrix* res;
329 
330  if (c->cov_diag)
331  nval = c->ncov;
332  else
333  nval = (c->ncov*(c->ncov+1))/2;
334 
335  vals = MALLOC_30(nval,double);
336  if (c->cov_diag) {
337  for (k = 0; k < nval; k++)
338  vals[k] = c->cov_diag[k];
339  res = mne_new_cov(c->kind,c->ncov,c->names,NULL,vals);
340  }
341  else {
342  for (k = 0; k < nval; k++)
343  vals[k] = c->cov[k];
344  res = mne_new_cov(c->kind,c->ncov,c->names,vals,NULL);
345  }
346  /*
347  * Duplicate additional items
348  */
349  if (c->ch_class) {
350  res->ch_class = MALLOC_30(c->ncov,int);
351  for (k = 0; k < c->ncov; k++)
352  res->ch_class[k] = c->ch_class[k];
353  }
354  res->bads = c->bads;
355  res->nbad = c->nbad;
356  res->proj = MneProjOp::mne_dup_proj_op(c->proj);
357  res->sss = new MneSssData(*(c->sss));
358 
359  return res;
360 }
361 
362 //=============================================================================================================
363 
364 int MneCovMatrix::mne_is_diag_cov(MneCovMatrix *c)
365 {
366  return c->cov_diag != NULL;
367 }
368 
369 //=============================================================================================================
370 
371 int MneCovMatrix::mne_add_inv_cov(MneCovMatrix *c)
372 /*
373  * Calculate the inverse square roots for whitening
374  */
375 {
376  double *src = c->lambda ? c->lambda : c->cov_diag;
377  int k;
378 
379  if (src == NULL) {
380  qCritical("Covariance matrix is not diagonal or not decomposed.");
381  return FAIL;
382  }
383  c->inv_lambda = REALLOC_30(c->inv_lambda,c->ncov,double);
384  for (k = 0; k < c->ncov; k++) {
385  if (src[k] <= 0.0)
386  c->inv_lambda[k] = 0.0;
387  else
388  c->inv_lambda[k] = 1.0/sqrt(src[k]);
389  }
390  return OK;
391 }
392 
393 //=============================================================================================================
394 
395 int MneCovMatrix::condition_cov(MneCovMatrix *c, float rank_threshold, int use_rank)
396 {
397  double *scale = NULL;
398  double *cov = NULL;
399  double *lambda = NULL;
400  float **eigen = NULL;
401  double **data1 = NULL;
402  double **data2 = NULL;
403  double magscale,gradscale,eegscale;
404  int nmag,ngrad,neeg,nok;
405  int j,k;
406  int res = FAIL;
407 
408  if (c->cov_diag)
409  return OK;
410  if (!c->ch_class) {
411  qCritical("Channels not classified. Rank cannot be determined.");
412  return FAIL;
413  }
414  magscale = gradscale = eegscale = 0.0;
415  nmag = ngrad = neeg = 0;
416  for (k = 0; k < c->ncov; k++) {
417  if (c->ch_class[k] == MNE_COV_CH_MEG_MAG) {
418  magscale += c->cov[mne_lt_packed_index(k,k)]; nmag++;
419  }
420  else if (c->ch_class[k] == MNE_COV_CH_MEG_GRAD) {
421  gradscale += c->cov[mne_lt_packed_index(k,k)]; ngrad++;
422  }
423  else if (c->ch_class[k] == MNE_COV_CH_EEG) {
424  eegscale += c->cov[mne_lt_packed_index(k,k)]; neeg++;
425  }
426 #ifdef DEBUG
427  fprintf(stdout,"%d ",c->ch_class[k]);
428 #endif
429  }
430 #ifdef DEBUG
431  fprintf(stdout,"\n");
432 #endif
433  if (nmag > 0)
434  magscale = magscale > 0.0 ? sqrt(nmag/magscale) : 0.0;
435  if (ngrad > 0)
436  gradscale = gradscale > 0.0 ? sqrt(ngrad/gradscale) : 0.0;
437  if (neeg > 0)
438  eegscale = eegscale > 0.0 ? sqrt(neeg/eegscale) : 0.0;
439 #ifdef DEBUG
440  fprintf(stdout,"%d %g\n",nmag,magscale);
441  fprintf(stdout,"%d %g\n",ngrad,gradscale);
442  fprintf(stdout,"%d %g\n",neeg,eegscale);
443 #endif
444  scale = MALLOC_30(c->ncov,double);
445  for (k = 0; k < c->ncov; k++) {
446  if (c->ch_class[k] == MNE_COV_CH_MEG_MAG)
447  scale[k] = magscale;
448  else if (c->ch_class[k] == MNE_COV_CH_MEG_GRAD)
449  scale[k] = gradscale;
450  else if (c->ch_class[k] == MNE_COV_CH_EEG)
451  scale[k] = eegscale;
452  else
453  scale[k] = 1.0;
454  }
455  cov = MALLOC_30(c->ncov*(c->ncov+1)/2.0,double);
456  lambda = MALLOC_30(c->ncov,double);
457  eigen = ALLOC_CMATRIX_30(c->ncov,c->ncov);
458  for (j = 0; j < c->ncov; j++)
459  for (k = 0; k <= j; k++)
460  cov[mne_lt_packed_index(j,k)] = c->cov[mne_lt_packed_index(j,k)]*scale[j]*scale[k];
461  if (mne_decompose_eigen(cov,lambda,eigen,c->ncov) == 0) {
462 #ifdef DEBUG
463  for (k = 0; k < c->ncov; k++)
464  fprintf(stdout,"%g ",lambda[k]/lambda[c->ncov-1]);
465  fprintf(stdout,"\n");
466 #endif
467  nok = 0;
468  for (k = c->ncov-1; k >= 0; k--) {
469  if (lambda[k] >= rank_threshold*lambda[c->ncov-1])
470  nok++;
471  else
472  break;
473  }
474  printf("\n\tEstimated covariance matrix rank = %d (%g)\n",nok,lambda[c->ncov-nok]/lambda[c->ncov-1]);
475  if (use_rank > 0 && use_rank < nok) {
476  nok = use_rank;
477  printf("\tUser-selected covariance matrix rank = %d (%g)\n",nok,lambda[c->ncov-nok]/lambda[c->ncov-1]);
478  }
479  /*
480  * Put it back together
481  */
482  for (j = 0; j < c->ncov-nok; j++)
483  lambda[j] = 0.0;
484  data1 = ALLOC_DCMATRIX_30(c->ncov,c->ncov);
485  for (j = 0; j < c->ncov; j++) {
486 #ifdef DEBUG
487  mne_print_vector(stdout,NULL,eigen[j],c->ncov);
488 #endif
489  for (k = 0; k < c->ncov; k++)
490  data1[j][k] = sqrt(lambda[j])*eigen[j][k];
491  }
492  data2 = mne_dmatt_dmat_mult2(data1,data1,c->ncov,c->ncov,c->ncov);
493 #ifdef DEBUG
494  printf(">>>\n");
495  for (j = 0; j < c->ncov; j++)
496  mne_print_dvector(stdout,NULL,data2[j],c->ncov);
497  printf(">>>\n");
498 #endif
499  /*
500  * Scale back
501  */
502  for (k = 0; k < c->ncov; k++)
503  if (scale[k] > 0.0)
504  scale[k] = 1.0/scale[k];
505  for (j = 0; j < c->ncov; j++)
506  for (k = 0; k <= j; k++)
507  if (c->cov[mne_lt_packed_index(j,k)] != 0.0)
508  c->cov[mne_lt_packed_index(j,k)] = scale[j]*scale[k]*data2[j][k];
509  res = nok;
510  }
511  FREE_30(cov);
512  FREE_30(lambda);
513  FREE_CMATRIX_30(eigen);
514  FREE_DCMATRIX_30(data1);
515  FREE_DCMATRIX_30(data2);
516  return res;
517 }
518 
519 //=============================================================================================================
520 
521 int MneCovMatrix::mne_decompose_eigen_cov_small(MneCovMatrix *c, float p_small, int use_rank)
522 /*
523  * Do the eigenvalue decomposition
524  */
525 {
526  int np,k,p,rank;
527  float rank_threshold = 1e-6;
528 
529  if (p_small < 0)
530  p_small = 1.0;
531 
532  if (!c)
533  return OK;
534  if (c->cov_diag)
535  return mne_add_inv_cov(c);
536  if (c->lambda && c->eigen) {
537  printf("\n\tEigenvalue decomposition had been precomputed.\n");
538  c->nzero = 0;
539  for (k = 0; k < c->ncov; k++, c->nzero++)
540  if (c->lambda[k] > 0)
541  break;
542  }
543  else {
544  FREE_30(c->lambda); c->lambda = NULL;
545  FREE_CMATRIX_30(c->eigen); c->eigen = NULL;
546 
547  if ((rank = condition_cov(c,rank_threshold,use_rank)) < 0)
548  return FAIL;
549 
550  np = c->ncov*(c->ncov+1)/2;
551  c->lambda = MALLOC_30(c->ncov,double);
552  c->eigen = ALLOC_CMATRIX_30(c->ncov,c->ncov);
553  if (mne_decompose_eigen (c->cov,c->lambda,c->eigen,c->ncov) != 0)
554  goto bad;
555  c->nzero = c->ncov - rank;
556  for (k = 0; k < c->nzero; k++)
557  c->lambda[k] = 0.0;
558  /*
559  * Find which eigenvectors correspond to EEG/MEG
560  */
561  {
562  float meglike,eeglike;
563  int nmeg,neeg;
564 
565  nmeg = neeg = 0;
566  for (k = c->nzero; k < c->ncov; k++) {
567  meglike = eeglike = 0.0;
568  for (p = 0; p < c->ncov; p++) {
569  if (c->ch_class[p] == MNE_COV_CH_EEG)
570  eeglike += std::fabs(c->eigen[k][p]);
571  else if (c->ch_class[p] == MNE_COV_CH_MEG_MAG || c->ch_class[p] == MNE_COV_CH_MEG_GRAD)
572  meglike += std::fabs(c->eigen[k][p]);
573  }
574  if (meglike > eeglike)
575  nmeg++;
576  else
577  neeg++;
578  }
579  printf("\t%d MEG and %d EEG-like channels remain in the whitened data\n",nmeg,neeg);
580  }
581  }
582  return mne_add_inv_cov(c);
583 
584 bad : {
585  FREE_30(c->lambda); c->lambda = NULL;
586  FREE_CMATRIX_30(c->eigen); c->eigen = NULL;
587  return FAIL;
588  }
589 }
590 
591 //=============================================================================================================
592 
593 int MneCovMatrix::mne_decompose_eigen_cov(MneCovMatrix *c)
594 
595 {
596  return mne_decompose_eigen_cov_small(c,-1.0,-1);
597 }
598 
599 //=============================================================================================================
600 
601 int MneCovMatrix::mne_lt_packed_index(int j, int k)
602 
603 {
604  if (j >= k)
605  return k + j*(j+1)/2;
606  else
607  return j + k*(k+1)/2;
608 }
MNE SSS Data (MneSssData) class declaration.
MNEProjItem class declaration.
MNEProjOp class declaration.
MNE SSS Data description.
Definition: mne_sss_data.h:82
FiffSparseMatrix class declaration.
MneCovMatrix class declaration.
Data associated with MNE computations for each mneMeasDataSet.
Covariance matrix storage.