MNE-CPP 0.1.9
A Framework for Electrophysiology
Loading...
Searching...
No Matches
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>
48#include <fiff/c/fiff_sparse_matrix.h>
49
50//=============================================================================================================
51// USED NAMESPACES
52//=============================================================================================================
53
54using namespace Eigen;
55using namespace FIFFLIB;
56using 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
74void mne_free_cmatrix_30 (float **m)
75{
76 if (m) {
77 FREE_30(*m);
78 FREE_30(m);
79 }
80}
81
82void 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
95static 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
112float **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
133double **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
156int 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
239double **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
274MneCovMatrix::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
302MneCovMatrix::~MneCovMatrix()
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
323MneCovMatrix *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
364int MneCovMatrix::mne_is_diag_cov(MneCovMatrix *c)
365{
366 return c->cov_diag != NULL;
367}
368
369//=============================================================================================================
370
371int 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
395int 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
521int 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 (void)np; // squash unused variable warning, what is this for?
552 c->lambda = MALLOC_30(c->ncov,double);
553 c->eigen = ALLOC_CMATRIX_30(c->ncov,c->ncov);
554 if (mne_decompose_eigen (c->cov,c->lambda,c->eigen,c->ncov) != 0)
555 goto bad;
556 c->nzero = c->ncov - rank;
557 for (k = 0; k < c->nzero; k++)
558 c->lambda[k] = 0.0;
559 /*
560 * Find which eigenvectors correspond to EEG/MEG
561 */
562 {
563 float meglike,eeglike;
564 int nmeg,neeg;
565
566 nmeg = neeg = 0;
567 for (k = c->nzero; k < c->ncov; k++) {
568 meglike = eeglike = 0.0;
569 for (p = 0; p < c->ncov; p++) {
570 if (c->ch_class[p] == MNE_COV_CH_EEG)
571 eeglike += std::fabs(c->eigen[k][p]);
572 else if (c->ch_class[p] == MNE_COV_CH_MEG_MAG || c->ch_class[p] == MNE_COV_CH_MEG_GRAD)
573 meglike += std::fabs(c->eigen[k][p]);
574 }
575 if (meglike > eeglike)
576 nmeg++;
577 else
578 neeg++;
579 }
580 printf("\t%d MEG and %d EEG-like channels remain in the whitened data\n",nmeg,neeg);
581 }
582 }
583 return mne_add_inv_cov(c);
584
585bad : {
586 FREE_30(c->lambda); c->lambda = NULL;
587 FREE_CMATRIX_30(c->eigen); c->eigen = NULL;
588 return FAIL;
589 }
590}
591
592//=============================================================================================================
593
594int MneCovMatrix::mne_decompose_eigen_cov(MneCovMatrix *c)
595
596{
597 return mne_decompose_eigen_cov_small(c,-1.0,-1);
598}
599
600//=============================================================================================================
601
602int MneCovMatrix::mne_lt_packed_index(int j, int k)
603
604{
605 if (j >= k)
606 return k + j*(j+1)/2;
607 else
608 return j + k*(k+1)/2;
609}
fiff_sparse_matrix.h
FiffSparseMatrix class declaration.
k
int k
Definition fiff_tag.cpp:324
mne_cov_matrix.h
MneCovMatrix class declaration.
mne_proj_op.h
MNEProjOp class declaration.
mne_proj_item.h
MNEProjItem class declaration.
mne_sss_data.h
MNE SSS Data (MneSssData) class declaration.
FIFFLIB::FiffSparseMatrix
Data associated with MNE computations for each mneMeasDataSet.
Definition fiff_sparse_matrix.h:75
MNELIB::MneCovMatrix
Covariance matrix storage.
Definition mne_cov_matrix.h:97
MNELIB::MneCovMatrix::MneCovMatrix
MneCovMatrix(int p_kind, int p_ncov, const QStringList &p_names, double *p_cov, double *p_cov_diag, FIFFLIB::FiffSparseMatrix *p_cov_sparse)
Definition mne_cov_matrix.cpp:274
MNELIB::MneSssData
MNE SSS Data description.
Definition mne_sss_data.h:83
Generated on Sat Feb 14 2026 16:30:24 for MNE-CPP by doxygen 1.9.8