71 using namespace Eigen;
99 QList<FiffDirEntry> t_Dir;
101 if(!t_pStream->open(t_Tree, t_Dir))
103 printf(
"\tNot able to open IODevice.\n");
108 printf(
"\tFiff covariance not found.\n");
115 : QSharedData(p_FiffCov)
116 , kind(p_FiffCov.kind)
117 , diag(p_FiffCov.diag)
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)
125 , eigvec(p_FiffCov.eigvec)
164 res.
dim = sel.size();
166 for(qint32 k = 0; k < sel.size(); ++k)
167 res.
names << this->names[sel(k)];
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));
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;
190 VectorXi C_ch_idx = VectorXi::Zero(p_NoiseCov.
names.size());
192 for(qint32 i = 0; i < p_ChNames.size(); ++i)
194 qint32 idx = p_NoiseCov.
names.indexOf(p_ChNames[i]);
197 C_ch_idx[count] = idx;
201 C_ch_idx.conservativeResize(count);
203 MatrixXd C(count, count);
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));
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);
223 printf(
"Created an SSP operator (subspace dimension = %d)\n", ncomp);
224 C = proj * (C * proj.transpose());
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);
230 QStringList meg_names, eeg_names;
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());
236 for(qint32 k = 0; k < C.rows(); ++k)
238 if(meg_names.indexOf(p_ChNames[k]) > -1)
240 C_meg_idx[count] = k;
245 C_meg_idx.conservativeResize(count);
247 C_meg_idx = VectorXi();
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());
254 for(qint32 k = 0; k < C.rows(); ++k)
256 if(eeg_names.indexOf(p_ChNames[k]) > -1)
258 C_eeg_idx[count] = k;
264 C_eeg_idx.conservativeResize(count);
266 C_eeg_idx = VectorXi();
268 bool has_meg = C_meg_idx.size() > 0;
269 bool has_eeg = C_eeg_idx.size() > 0;
271 MatrixXd C_meg, C_eeg;
272 VectorXd C_meg_eig, C_eeg_eig;
273 MatrixXd C_meg_eigvec, C_eeg_eigvec;
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);
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);
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);
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];
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];
315 if (C_meg_idx.size() + C_eeg_idx.size() != n_chan)
317 printf(
"Error in FiffCov::prepare_noise_cov: channel sizes do no match!\n");
322 p_NoiseCov.
dim = p_ChNames.size();
323 p_NoiseCov.
diag =
false;
324 p_NoiseCov.
names = p_ChNames;
336 if(p_exclude.size() == 0)
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];
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);
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)];
360 QStringList ch_names = cov_good.
names;
362 std::vector<qint32> idx_eeg, idx_mag, idx_grad;
363 for(qint32 i = 0; i < ch_names.size(); ++i)
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);
373 MatrixXd C(cov_good.
data);
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");
378 QList<FiffProj> t_listProjs;
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));
394 QMap<QString, QPair<double, std::vector<qint32> > >::Iterator it;
395 for(it = regData.begin(); it != regData.end(); ++it)
397 QString desc(it.key());
398 double reg = it.value().first;
399 std::vector<qint32> idx = it.value().second;
401 if(idx.size() == 0 || reg == 0.0)
402 printf(
"\tNothing to regularize within %s data.\n", desc.toLatin1().constData());
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]);
415 QStringList this_ch_names;
416 for(quint32 k = 0; k < idx.size(); ++k)
417 this_ch_names << ch_names[idx[k]];
422 JacobiSVD<MatrixXd> svd(P, ComputeFullU);
424 VectorXd t_s = svd.singularValues();
425 MatrixXd t_U = svd.matrixU();
426 MNEMath::sort<double>(t_s, t_U);
428 U = t_U.block(0,0, t_U.rows(), t_U.cols()-ncomp);
432 printf(
"\tCreated an SSP operator for %s (dimension = %d).\n", desc.toLatin1().constData(), ncomp);
433 this_C = U.transpose() * (this_C * U);
437 double sigma = this_C.diagonal().mean();
438 this_C.diagonal() = this_C.diagonal().array() + reg * sigma;
439 if(p_bProj && ncomp > 0)
440 this_C = U * (this_C * U.transpose());
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);
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);
#define FIFFV_MNE_NOISE_COV
FIFF measurement file information.
static void activate_projs(QList< FiffProj > &p_qListFiffProj)
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
QSharedPointer< FiffStream > SPtr
FiffCov & operator=(const FiffCov &rhs)
FiffInfoBase class declaration.
FiffStream class declaration.
Directory tree structure.
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)
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
qint32 make_projector(MatrixXd &proj) const
FiffCov pick_channels(const QStringList &p_include=defaultQStringList, const QStringList &p_exclude=defaultQStringList)
MNEMath class declaration.
static RowVectorXi pick_channels(const QStringList &ch_names, const QStringList &include=defaultQStringList, const QStringList &exclude=defaultQStringList)
FiffCov class declaration.