MNE-CPP  beta 1.0
fiff_cov.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //*************************************************************************************************************
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "fiff_cov.h"
42 #include "fiff_stream.h"
43 #include "fiff_info_base.h"
44 
45 #include <utils/mnemath.h>
46 
47 
48 //*************************************************************************************************************
49 //=============================================================================================================
50 // Qt INCLUDES
51 //=============================================================================================================
52 
53 #include <QPair>
54 
55 
56 //*************************************************************************************************************
57 //=============================================================================================================
58 // Eigen INCLUDES
59 //=============================================================================================================
60 
61 #include <Eigen/SVD>
62 
63 
64 //*************************************************************************************************************
65 //=============================================================================================================
66 // USED NAMESPACES
67 //=============================================================================================================
68 
69 using namespace FIFFLIB;
70 using namespace UTILSLIB;
71 using namespace Eigen;
72 
73 
74 //*************************************************************************************************************
75 //=============================================================================================================
76 // DEFINE MEMBER METHODS
77 //=============================================================================================================
78 
80 : kind(-1)
81 , diag(false)
82 , dim(-1)
83 , nfree(-1)
84 {
85 
86 }
87 
88 
89 //*************************************************************************************************************
90 
91 FiffCov::FiffCov(QIODevice &p_IODevice)
92 : kind(-1)
93 , diag(false)
94 , dim(-1)
95 , nfree(-1)
96 {
97  FiffStream::SPtr t_pStream(new FiffStream(&p_IODevice));
98  FiffDirTree t_Tree;
99  QList<FiffDirEntry> t_Dir;
100 
101  if(!t_pStream->open(t_Tree, t_Dir))
102  {
103  printf("\tNot able to open IODevice.\n");//ToDo Throw here
104  return;
105  }
106 
107  if(!t_pStream->read_cov(t_Tree, FIFFV_MNE_NOISE_COV, *this))
108  printf("\tFiff covariance not found.\n");//ToDo Throw here
109 }
110 
111 
112 //*************************************************************************************************************
113 
114 FiffCov::FiffCov(const FiffCov &p_FiffCov)
115 : QSharedData(p_FiffCov)
116 , kind(p_FiffCov.kind)
117 , diag(p_FiffCov.diag)
118 , dim(p_FiffCov.dim)
119 , names(p_FiffCov.names)
120 , data(p_FiffCov.data)
121 , projs(p_FiffCov.projs)
122 , bads(p_FiffCov.bads)
123 , nfree(p_FiffCov.nfree)
124 , eig(p_FiffCov.eig)
125 , eigvec(p_FiffCov.eigvec)
126 {
127 
128 }
129 
130 
131 //*************************************************************************************************************
132 
134 {
135 }
136 
137 
138 //*************************************************************************************************************
139 
141 {
142  kind = -1;
143  diag = false;
144  dim = -1;
145  names.clear();
146  data = MatrixXd();
147  projs.clear();
148  bads.clear();
149  nfree = -1;
150  eig = VectorXd();
151  eigvec = MatrixXd();
152 }
153 
154 
155 //*************************************************************************************************************
156 
157 FiffCov FiffCov::pick_channels(const QStringList &p_include, const QStringList &p_exclude)
158 {
159  RowVectorXi sel = FiffInfoBase::pick_channels(this->names, p_include, p_exclude);
160  FiffCov res;//No deep copy here - since almost everything else is adapted anyway
161 
162  res.kind = this->kind;
163  res.diag = this->diag;
164  res.dim = sel.size();
165 
166  for(qint32 k = 0; k < sel.size(); ++k)
167  res.names << this->names[sel(k)];
168 
169  res.data.resize(res.dim, res.dim);
170  for(qint32 i = 0; i < res.dim; ++i)
171  for(qint32 j = 0; j < res.dim; ++j)
172  res.data(i, j) = this->data(sel(i), sel(j));
173  res.projs = this->projs;
174 
175  for(qint32 k = 0; k < this->bads.size(); ++k)
176  if(res.names.contains(this->bads[k]))
177  res.bads << this->bads[k];
178  res.nfree = this->nfree;
179 
180  return res;
181 }
182 
183 
184 //*************************************************************************************************************
185 
186 FiffCov FiffCov::prepare_noise_cov(const FiffInfo &p_Info, const QStringList &p_ChNames) const
187 {
188  FiffCov p_NoiseCov(*this);
189 
190  VectorXi C_ch_idx = VectorXi::Zero(p_NoiseCov.names.size());
191  qint32 count = 0;
192  for(qint32 i = 0; i < p_ChNames.size(); ++i)
193  {
194  qint32 idx = p_NoiseCov.names.indexOf(p_ChNames[i]);
195  if(idx > -1)
196  {
197  C_ch_idx[count] = idx;
198  ++count;
199  }
200  }
201  C_ch_idx.conservativeResize(count);
202 
203  MatrixXd C(count, count);
204 
205  if(!p_NoiseCov.diag)
206  for(qint32 i = 0; i < count; ++i)
207  for(qint32 j = 0; j < count; ++j)
208  C(i,j) = p_NoiseCov.data(C_ch_idx(i), C_ch_idx(j));
209  else
210  {
211  qWarning("Warning in FiffCov::prepare_noise_cov: This has to be debugged - not done before!");
212  C = MatrixXd::Zero(count, count);
213  for(qint32 i = 0; i < count; ++i)
214  C.diagonal()[i] = p_NoiseCov.data(C_ch_idx(i),0);
215  }
216 
217  MatrixXd proj;
218  qint32 ncomp = p_Info.make_projector(proj, p_ChNames);
219 
220  //Create the projection operator
221  if (ncomp > 0)
222  {
223  printf("Created an SSP operator (subspace dimension = %d)\n", ncomp);
224  C = proj * (C * proj.transpose());
225  }
226 
227  RowVectorXi pick_meg = p_Info.pick_types(true, false, false, defaultQStringList, p_Info.bads);
228  RowVectorXi pick_eeg = p_Info.pick_types(false, true, false, defaultQStringList, p_Info.bads);
229 
230  QStringList meg_names, eeg_names;
231 
232  for(qint32 i = 0; i < pick_meg.size(); ++i)
233  meg_names << p_Info.chs[pick_meg[i]].ch_name;
234  VectorXi C_meg_idx = VectorXi::Zero(p_NoiseCov.names.size());
235  count = 0;
236  for(qint32 k = 0; k < C.rows(); ++k)
237  {
238  if(meg_names.indexOf(p_ChNames[k]) > -1)
239  {
240  C_meg_idx[count] = k;
241  ++count;
242  }
243  }
244  if(count > 0)
245  C_meg_idx.conservativeResize(count);
246  else
247  C_meg_idx = VectorXi();
248 
249  //
250  for(qint32 i = 0; i < pick_eeg.size(); ++i)
251  eeg_names << p_Info.chs[pick_eeg(0,i)].ch_name;
252  VectorXi C_eeg_idx = VectorXi::Zero(p_NoiseCov.names.size());
253  count = 0;
254  for(qint32 k = 0; k < C.rows(); ++k)
255  {
256  if(eeg_names.indexOf(p_ChNames[k]) > -1)
257  {
258  C_eeg_idx[count] = k;
259  ++count;
260  }
261  }
262 
263  if(count > 0)
264  C_eeg_idx.conservativeResize(count);
265  else
266  C_eeg_idx = VectorXi();
267 
268  bool has_meg = C_meg_idx.size() > 0;
269  bool has_eeg = C_eeg_idx.size() > 0;
270 
271  MatrixXd C_meg, C_eeg;
272  VectorXd C_meg_eig, C_eeg_eig;
273  MatrixXd C_meg_eigvec, C_eeg_eigvec;
274  if (has_meg)
275  {
276  count = C_meg_idx.rows();
277  C_meg = MatrixXd(count,count);
278  for(qint32 i = 0; i < count; ++i)
279  for(qint32 j = 0; j < count; ++j)
280  C_meg(i,j) = C(C_meg_idx(i), C_meg_idx(j));
281  MNEMath::get_whitener(C_meg, false, QString("MEG"), C_meg_eig, C_meg_eigvec);
282  }
283 
284  if (has_eeg)
285  {
286  count = C_eeg_idx.rows();
287  C_eeg = MatrixXd(count,count);
288  for(qint32 i = 0; i < count; ++i)
289  for(qint32 j = 0; j < count; ++j)
290  C_eeg(i,j) = C(C_eeg_idx(i), C_eeg_idx(j));
291  MNEMath::get_whitener(C_eeg, false, QString("EEG"), C_eeg_eig, C_eeg_eigvec);
292  }
293 
294  qint32 n_chan = p_ChNames.size();
295  p_NoiseCov.eigvec = MatrixXd::Zero(n_chan, n_chan);
296  p_NoiseCov.eig = VectorXd::Zero(n_chan);
297 
298  if(has_meg)
299  {
300  for(qint32 i = 0; i < C_meg_idx.rows(); ++i)
301  for(qint32 j = 0; j < C_meg_idx.rows(); ++j)
302  p_NoiseCov.eigvec(C_meg_idx[i], C_meg_idx[j]) = C_meg_eigvec(i, j);
303  for(qint32 i = 0; i < C_meg_idx.rows(); ++i)
304  p_NoiseCov.eig(C_meg_idx[i]) = C_meg_eig[i];
305  }
306  if(has_eeg)
307  {
308  for(qint32 i = 0; i < C_eeg_idx.rows(); ++i)
309  for(qint32 j = 0; j < C_eeg_idx.rows(); ++j)
310  p_NoiseCov.eigvec(C_eeg_idx[i], C_eeg_idx[j]) = C_eeg_eigvec(i, j);
311  for(qint32 i = 0; i < C_eeg_idx.rows(); ++i)
312  p_NoiseCov.eig(C_eeg_idx[i]) = C_eeg_eig[i];
313  }
314 
315  if (C_meg_idx.size() + C_eeg_idx.size() != n_chan)
316  {
317  printf("Error in FiffCov::prepare_noise_cov: channel sizes do no match!\n");//ToDo Throw here
318  return FiffCov();
319  }
320 
321  p_NoiseCov.data = C;
322  p_NoiseCov.dim = p_ChNames.size();
323  p_NoiseCov.diag = false;
324  p_NoiseCov.names = p_ChNames;
325 
326  return p_NoiseCov;
327 }
328 
329 
330 //*************************************************************************************************************
331 
332 FiffCov FiffCov::regularize(const FiffInfo& p_info, double p_fRegMag, double p_fRegGrad, double p_fRegEeg, bool p_bProj, QStringList p_exclude) const
333 {
334  FiffCov cov(*this);
335 
336  if(p_exclude.size() == 0)
337  {
338  p_exclude = p_info.bads;
339  for(qint32 i = 0; i < cov.bads.size(); ++i)
340  if(!p_exclude.contains(cov.bads[i]))
341  p_exclude << cov.bads[i];
342  }
343 
344  RowVectorXi sel_eeg = p_info.pick_types(false, true, false, defaultQStringList, p_exclude);
345  RowVectorXi sel_mag = p_info.pick_types(QString("mag"), false, false, defaultQStringList, p_exclude);
346  RowVectorXi sel_grad = p_info.pick_types(QString("grad"), false, false, defaultQStringList, p_exclude);
347 
348  QStringList info_ch_names = p_info.ch_names;
349  QStringList ch_names_eeg, ch_names_mag, ch_names_grad;
350  for(qint32 i = 0; i < sel_eeg.size(); ++i)
351  ch_names_eeg << info_ch_names[sel_eeg(i)];
352  for(qint32 i = 0; i < sel_mag.size(); ++i)
353  ch_names_mag << info_ch_names[sel_mag(i)];
354  for(qint32 i = 0; i < sel_grad.size(); ++i)
355  ch_names_grad << info_ch_names[sel_grad(i)];
356 
357  // This actually removes bad channels from the cov, which is not backward
358  // compatible, so let's leave all channels in
359  FiffCov cov_good = cov.pick_channels(info_ch_names, p_exclude);
360  QStringList ch_names = cov_good.names;
361 
362  std::vector<qint32> idx_eeg, idx_mag, idx_grad;
363  for(qint32 i = 0; i < ch_names.size(); ++i)
364  {
365  if(ch_names_eeg.contains(ch_names[i]))
366  idx_eeg.push_back(i);
367  else if(ch_names_mag.contains(ch_names[i]))
368  idx_mag.push_back(i);
369  else if(ch_names_grad.contains(ch_names[i]))
370  idx_grad.push_back(i);
371  }
372 
373  MatrixXd C(cov_good.data);
374 
375  if((unsigned) C.rows() != idx_eeg.size() + idx_mag.size() + idx_grad.size())
376  printf("Error in FiffCov::regularize: Channel dimensions do not fit.\n");//ToDo Throw
377 
378  QList<FiffProj> t_listProjs;
379  if(p_bProj)
380  {
381  t_listProjs = p_info.projs + cov_good.projs;
382  FiffProj::activate_projs(t_listProjs);
383  }
384 
385  //Build regularization MAP
386  QMap<QString, QPair<double, std::vector<qint32> > > regData;
387  regData.insert("EEG", QPair<double, std::vector<qint32> >(p_fRegEeg, idx_eeg));
388  regData.insert("MAG", QPair<double, std::vector<qint32> >(p_fRegMag, idx_mag));
389  regData.insert("GRAD", QPair<double, std::vector<qint32> >(p_fRegGrad, idx_grad));
390 
391  //
392  //Regularize
393  //
394  QMap<QString, QPair<double, std::vector<qint32> > >::Iterator it;
395  for(it = regData.begin(); it != regData.end(); ++it)
396  {
397  QString desc(it.key());
398  double reg = it.value().first;
399  std::vector<qint32> idx = it.value().second;
400 
401  if(idx.size() == 0 || reg == 0.0)
402  printf("\tNothing to regularize within %s data.\n", desc.toLatin1().constData());
403  else
404  {
405  printf("\tRegularize %s: %f\n", desc.toLatin1().constData(), reg);
406  MatrixXd this_C(idx.size(), idx.size());
407  for(quint32 i = 0; i < idx.size(); ++i)
408  for(quint32 j = 0; j < idx.size(); ++j)
409  this_C(i,j) = cov_good.data(idx[i], idx[j]);
410 
411  MatrixXd U;
412  qint32 ncomp;
413  if(p_bProj)
414  {
415  QStringList this_ch_names;
416  for(quint32 k = 0; k < idx.size(); ++k)
417  this_ch_names << ch_names[idx[k]];
418 
419  MatrixXd P;
420  ncomp = FiffProj::make_projector(t_listProjs, this_ch_names, P); //ToDo: Synchronize with mne-python and debug
421 
422  JacobiSVD<MatrixXd> svd(P, ComputeFullU);
423  //Sort singular values and singular vectors
424  VectorXd t_s = svd.singularValues();
425  MatrixXd t_U = svd.matrixU();
426  MNEMath::sort<double>(t_s, t_U);
427 
428  U = t_U.block(0,0, t_U.rows(), t_U.cols()-ncomp);
429 
430  if (ncomp > 0)
431  {
432  printf("\tCreated an SSP operator for %s (dimension = %d).\n", desc.toLatin1().constData(), ncomp);
433  this_C = U.transpose() * (this_C * U);
434  }
435  }
436 
437  double sigma = this_C.diagonal().mean();
438  this_C.diagonal() = this_C.diagonal().array() + reg * sigma; // modify diag inplace
439  if(p_bProj && ncomp > 0)
440  this_C = U * (this_C * U.transpose());
441 
442  for(qint32 i = 0; i < this_C.rows(); ++i)
443  for(qint32 j = 0; j < this_C.cols(); ++j)
444  C(idx[i],idx[j]) = this_C(i,j);
445  }
446  }
447 
448  // Put data back in correct locations
449  RowVectorXi idx = FiffInfo::pick_channels(cov.names, info_ch_names, p_exclude);
450  for(qint32 i = 0; i < idx.size(); ++i)
451  for(qint32 j = 0; j < idx.size(); ++j)
452  cov.data(idx[i], idx[j]) = C(i, j);
453 
454  return cov;
455 }
456 
457 
458 //*************************************************************************************************************
459 
461 {
462  if (this != &rhs) // protect against invalid self-assignment
463  {
464  kind = rhs.kind;
465  diag = rhs.diag;
466  dim = rhs.dim;
467  names = rhs.names;
468  data = rhs.data;
469  projs = rhs.projs;
470  bads = rhs.bads;
471  nfree = rhs.nfree;
472  eig = rhs.eig;
473  eigvec = rhs.eigvec;
474  }
475  // to support chained assignment operators (a=b=c), always return *this
476  return *this;
477 }
#define FIFFV_MNE_NOISE_COV
MatrixXd data
Definition: fiff_cov.h:211
FIFF measurement file information.
Definition: fiff_info.h:96
fiff_int_t kind
Definition: fiff_cov.h:207
QList< FiffProj > projs
Definition: fiff_info.h:266
fiff_int_t nfree
Definition: fiff_cov.h:214
static void activate_projs(QList< FiffProj > &p_qListFiffProj)
Definition: fiff_proj.cpp:112
QStringList names
Definition: fiff_cov.h:210
QStringList bads
Definition: fiff_cov.h:213
FiffCov regularize(const FiffInfo &p_info, double p_fMag=0.1, double p_fGrad=0.1, double p_fEeg=0.1, bool p_bProj=true, QStringList p_exclude=defaultQStringList) const
Definition: fiff_cov.cpp:332
QSharedPointer< FiffStream > SPtr
Definition: fiff_stream.h:132
FiffCov & operator=(const FiffCov &rhs)
Definition: fiff_cov.cpp:460
FiffInfoBase class declaration.
FiffStream class declaration.
Directory tree structure.
Definition: fiff_dir_tree.h:80
QList< FiffChInfo > chs
fiff_int_t dim
Definition: fiff_cov.h:209
static fiff_int_t make_projector(const QList< FiffProj > &projs, const QStringList &ch_names, MatrixXd &proj, const QStringList &bads=defaultQStringList, MatrixXd &U=defaultMatrixXd, bool include_active=true)
Definition: fiff_proj.cpp:125
covariance data
Definition: fiff_cov.h:94
MatrixXd eigvec
Definition: fiff_cov.h:216
Definition: fiff.h:98
RowVectorXi pick_types(const QString meg, bool eeg=false, bool stim=false, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList) const
FiffCov prepare_noise_cov(const FiffInfo &p_info, const QStringList &p_chNames) const
Definition: fiff_cov.cpp:186
qint32 make_projector(MatrixXd &proj) const
Definition: fiff_info.h:277
FiffCov pick_channels(const QStringList &p_include=defaultQStringList, const QStringList &p_exclude=defaultQStringList)
Definition: fiff_cov.cpp:157
MNEMath class declaration.
QList< FiffProj > projs
Definition: fiff_cov.h:212
FIFF File I/O routines.
Definition: fiff_stream.h:129
VectorXd eig
Definition: fiff_cov.h:215
static RowVectorXi pick_channels(const QStringList &ch_names, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList)
FiffCov class declaration.