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