MNE-CPP  beta 1.0
fiff_proj.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //*************************************************************************************************************
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "fiff_proj.h"
42 #include <utils/mnemath.h>
43 
44 
45 //*************************************************************************************************************
46 //=============================================================================================================
47 // Eigen INCLUDES
48 //=============================================================================================================
49 
50 #include <Eigen/SVD>
51 
52 
53 //*************************************************************************************************************
54 //=============================================================================================================
55 // USED NAMESPACES
56 //=============================================================================================================
57 
58 using namespace UTILSLIB;
59 using namespace FIFFLIB;
60 using namespace Eigen;
61 
62 
63 //*************************************************************************************************************
64 //=============================================================================================================
65 // DEFINE MEMBER METHODS
66 //=============================================================================================================
67 
68 FiffProj::FiffProj()
69 : kind(-1)
70 , active(false)
71 , desc("")
72 , data(new FiffNamedMatrix)
73 {
74 
75 }
76 
77 
78 //*************************************************************************************************************
79 
80 FiffProj::FiffProj(const FiffProj& p_FiffProj)
81 : kind(p_FiffProj.kind)
82 , active(p_FiffProj.active)
83 , desc(p_FiffProj.desc)
84 , data(p_FiffProj.data)
85 {
86 
87 }
88 
89 
90 //*************************************************************************************************************
91 
92 FiffProj::FiffProj( fiff_int_t p_kind, bool p_active, QString p_desc, FiffNamedMatrix& p_data)
93 : kind(p_kind)
94 , active(p_active)
95 , desc(p_desc)
96 , data(new FiffNamedMatrix(p_data))
97 {
98 
99 }
100 
101 
102 //*************************************************************************************************************
103 
105 {
106 
107 }
108 
109 
110 //*************************************************************************************************************
111 
112 void FiffProj::activate_projs(QList<FiffProj> &p_qListFiffProj)
113 {
114  // Activate the projection items
115  QList<FiffProj>::Iterator it;
116  for(it = p_qListFiffProj.begin(); it != p_qListFiffProj.end(); ++it)
117  it->active = true;
118 
119  printf("\t%d projection items activated.\n", p_qListFiffProj.size());
120 }
121 
122 
123 //*************************************************************************************************************
124 
125 fiff_int_t FiffProj::make_projector(const QList<FiffProj>& projs, const QStringList& ch_names, MatrixXd& proj, const QStringList& bads, MatrixXd& U, bool include_active)
126 {
127  fiff_int_t nchan = ch_names.size();
128  if (nchan == 0)
129  {
130  printf("No channel names specified\n");//ToDo throw here
131  return 0;
132  }
133 
134 // if(proj)
135 // delete proj;
136  proj = MatrixXd::Identity(nchan,nchan);
137  fiff_int_t nproj = 0;
138  U = MatrixXd();
139 
140  //
141  // Check trivial cases first
142  //
143  if (projs.size() == 0)
144  return 0;
145 
146  fiff_int_t nvec = 0;
147  fiff_int_t k, l;
148  for (k = 0; k < projs.size(); ++k)
149  {
150  if (!projs[k].active || include_active)
151  {
152  ++nproj;
153  nvec += projs[k].data->nrow;
154  }
155  }
156 
157  if (nproj == 0)
158  return 0;
159 
160  //
161  // Pick the appropriate entries
162  //
163  MatrixXd vecs = MatrixXd::Zero(nchan,nvec);
164  nvec = 0;
165  fiff_int_t nonzero = 0;
166  qint32 p, c, i, j, v;
167  double onesize;
168  bool isBad = false;
169  RowVectorXi sel(nchan);
170  RowVectorXi vecSel(nchan);
171  sel.setConstant(-1);
172  vecSel.setConstant(-1);
173  for (k = 0; k < projs.size(); ++k)
174  {
175  if (!projs[k].active || include_active)
176  {
177  FiffProj one = projs[k];
178 
179  QMap<QString, int> uniqueMap;
180  for(l = 0; l < one.data->col_names.size(); ++l)
181  uniqueMap[one.data->col_names[l] ] = 0;
182 
183  if (one.data->col_names.size() != uniqueMap.keys().size())
184  {
185  printf("Channel name list in projection item %d contains duplicate items",k);
186  return 0;
187  }
188 
189  //
190  // Get the two selection vectors to pick correct elements from
191  // the projection vectors omitting bad channels
192  //
193  sel.resize(nchan);
194  vecSel.resize(nchan);
195  sel.setConstant(-1);
196  vecSel.setConstant(-1);
197  p = 0;
198  for (c = 0; c < nchan; ++c)
199  {
200  for (i = 0; i < one.data->col_names.size(); ++i)
201  {
202  if (QString::compare(ch_names.at(c),one.data->col_names[i]) == 0)
203  {
204  isBad = false;
205  for (j = 0; j < bads.size(); ++j)
206  {
207  if (QString::compare(ch_names.at(c),bads.at(j)) == 0)
208  {
209  isBad = true;
210  }
211  }
212 
213  if (!isBad && sel[p] != c)
214  {
215  sel[p] = c;
216  vecSel[p] = i;
217  ++p;
218  }
219 
220  }
221  }
222  }
223  sel.conservativeResize(p);
224  vecSel.conservativeResize(p);
225  //
226  // If there is something to pick, pickit
227  //
228  if (sel.cols() > 0)
229  for (v = 0; v < one.data->nrow; ++v)
230  for (i = 0; i < p; ++i)
231  vecs(sel[i],nvec+v) = one.data->data(v,vecSel[i]);
232 
233  //
234  // Rescale for more straightforward detection of small singular values
235  //
236  for (v = 0; v < one.data->nrow; ++v)
237  {
238  onesize = sqrt((vecs.col(nvec+v).transpose()*vecs.col(nvec+v))(0,0));
239  if (onesize > 0.0)
240  {
241  vecs.col(nvec+v) = vecs.col(nvec+v)/onesize;
242  ++nonzero;
243  }
244  }
245  nvec += one.data->nrow;
246  }
247  }
248  //
249  // Check whether all of the vectors are exactly zero
250  //
251  if (nonzero == 0)
252  return 0;
253 
254  //
255  // Reorthogonalize the vectors
256  //
257  JacobiSVD<MatrixXd> svd(vecs.block(0,0,vecs.rows(),nvec), ComputeFullU);
258  //Sort singular values and singular vectors
259  VectorXd S = svd.singularValues();
260  MatrixXd t_U = svd.matrixU();
261  MNEMath::sort<double>(S, t_U);
262 
263  //
264  // Throw away the linearly dependent guys
265  //
266  nproj = 0;
267  for(k = 0; k < S.size(); ++k)
268  if (S[k]/S[0] > 1e-2)
269  ++nproj;
270 
271  U = t_U.block(0, 0, t_U.rows(), nproj);
272 
273  //
274  // Here is the celebrated result
275  //
276  proj -= U*U.transpose();
277 
278  return nproj;
279 }
FiffProj class declaration.
static void activate_projs(QList< FiffProj > &p_qListFiffProj)
Definition: fiff_proj.cpp:112
FiffNamedMatrix::SDPtr data
Definition: fiff_proj.h:168
SSP projector data.
Definition: fiff_proj.h:89
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
Definition: fiff.h:98
MNEMath class declaration.