MNE-CPP  0.1.9
A Framework for Electrophysiology
hpifitdata.cpp
Go to the documentation of this file.
1 //=============================================================================================================
37 //=============================================================================================================
38 // INCLUDES
39 //=============================================================================================================
40 
41 #include "hpifitdata.h"
42 #include "hpifit.h"
43 #include "sensorset.h"
44 
45 #include <utils/mnemath.h>
46 
47 #include <iostream>
48 #include <algorithm>
49 
50 //=============================================================================================================
51 // EIGEN INCLUDES
52 //=============================================================================================================
53 
54 //=============================================================================================================
55 // QT INCLUDES
56 //=============================================================================================================
57 
58 #include <qmath.h>
59 
60 //=============================================================================================================
61 // USED NAMESPACES
62 //=============================================================================================================
63 
64 using namespace INVERSELIB;
65 
66 //=============================================================================================================
67 // DEFINE GLOBAL METHODS
68 //=============================================================================================================
69 
70 //=============================================================================================================
71 // DEFINE MEMBER METHODS
72 //=============================================================================================================
73 
75  : m_sensors(SensorSet())
76 {
77 }
78 
79 //=============================================================================================================
80 
82 {
83  // Initialize variables
84  Eigen::RowVectorXd vecCurrentCoil = this->m_coilPos;
85  Eigen::VectorXd vecCurrentData = this->m_sensorData;
86  SensorSet currentSensors = this->m_sensors;
87 
88  int iDisplay = 0;
89  int iMaxiter = m_iMaxIterations;
90  int iSimplexNumitr = 0;
91 
92  this->m_coilPos = fminsearch(vecCurrentCoil,
93  iMaxiter,
94  2 * iMaxiter * vecCurrentCoil.cols(),
95  iDisplay,
96  vecCurrentData,
97  this->m_matProjector,
98  currentSensors,
99  iSimplexNumitr);
100 
101  this->m_errorInfo = dipfitError(vecCurrentCoil,
102  vecCurrentData,
103  currentSensors,
104  this->m_matProjector);
105 
106  this->m_errorInfo.numIterations = iSimplexNumitr;
107 }
108 
109 //=============================================================================================================
110 
111 Eigen::MatrixXd HPIFitData::magnetic_dipole(Eigen::MatrixXd matPos,
112  Eigen::MatrixXd matPnt,
113  Eigen::MatrixXd matOri)
114 {
115  double u0 = 1e-7;
116  int iNchan;
117  Eigen::MatrixXd r, r2, r5, x, y, z, mx, my, mz, Tx, Ty, Tz, lf;
118 
119  iNchan = matPnt.rows();
120 
121  // Shift the magnetometers so that the dipole is in the origin
122  matPnt.array().col(0) -=matPos(0);
123  matPnt.array().col(1) -=matPos(1);
124  matPnt.array().col(2) -=matPos(2);
125 
126  r = matPnt.array().square().rowwise().sum().sqrt();
127 
128  r2 = r5 = x = y = z = mx = my = mz = Tx = Ty = Tz = lf = Eigen::MatrixXd::Zero(iNchan,3);
129 
130  for(int i = 0;i < iNchan;i++) {
131  r2.row(i).array().fill(pow(r(i),2));
132  r5.row(i).array().fill(pow(r(i),5));
133  }
134 
135  for(int i = 0;i < iNchan;i++) {
136  x.row(i).array().fill(matPnt(i,0));
137  y.row(i).array().fill(matPnt(i,1));
138  z.row(i).array().fill(matPnt(i,2));
139  }
140 
141  mx.col(0).array().fill(1);
142  my.col(1).array().fill(1);
143  mz.col(2).array().fill(1);
144 
145  Tx = 3 * x.cwiseProduct(matPnt) - mx.cwiseProduct(r2);
146  Ty = 3 * y.cwiseProduct(matPnt) - my.cwiseProduct(r2);
147  Tz = 3 * z.cwiseProduct(matPnt) - mz.cwiseProduct(r2);
148 
149  for(int i = 0;i < iNchan;i++) {
150  lf(i,0) = Tx.row(i).dot(matOri.row(i));
151  lf(i,1) = Ty.row(i).dot(matOri.row(i));
152  lf(i,2) = Tz.row(i).dot(matOri.row(i));
153  }
154 
155  for(int i = 0;i < iNchan;i++) {
156  for(int j = 0;j < 3;j++) {
157  lf(i,j) = u0 * lf(i,j)/(4 * M_PI * r5(i,j));
158  }
159  }
160 
161  return lf;
162 }
163 
164 //=============================================================================================================
165 
166 Eigen::MatrixXd HPIFitData::compute_leadfield(const Eigen::MatrixXd& matPos, const SensorSet& sensors)
167 {
168 
169  Eigen::MatrixXd matPnt, matOri, matLf;
170  matPnt = sensors.rmag(); // position of each integrationpoint
171  matOri = sensors.cosmag(); // mOrientation of each coil
172 
173  matLf = magnetic_dipole(matPos, matPnt, matOri);
174 
175  return matLf;
176 }
177 
178 //=============================================================================================================
179 
180 DipFitError HPIFitData::dipfitError(const Eigen::MatrixXd& matPos,
181  const Eigen::MatrixXd& matData,
182  const SensorSet& sensors,
183  const Eigen::MatrixXd& matProjectors)
184 {
185  // Variable Declaration
186  struct DipFitError e;
187  Eigen::MatrixXd matLfSensor, matDif;
188  Eigen::MatrixXd matLf(matData.size(),3);
189  int iNp = sensors.np();
190 
191  // calculate lf for all sensorpoints
192  matLfSensor = compute_leadfield(matPos, sensors);
193 
194  // apply averaging per coil
195  for(int i = 0; i < sensors.ncoils(); i++){
196  matLf.row(i) = sensors.w(i) * matLfSensor.block(i*iNp,0,iNp,matLfSensor.cols());
197  }
198  //matLf = sensors.tra * matLf;
199 
200  // Compute lead field for a magnetic dipole in infinite vacuum
201  e.moment = UTILSLIB::MNEMath::pinv(matLf) * matData;
202 
203  //matDif = matData - matLf * e.moment;
204  matDif = matData - matProjectors * matLf * e.moment;
205 
206  e.error = matDif.array().square().sum()/matData.array().square().sum();
207 
208  e.numIterations = 0;
209 
210  return e;
211 }
212 
213 //=============================================================================================================
214 
216 {
217  return (a.base_arr < b.base_arr);
218 }
219 
220 //=============================================================================================================
221 
222 Eigen::MatrixXd HPIFitData::fminsearch(const Eigen::MatrixXd& matPos,
223  int iMaxiter,
224  int iMaxfun,
225  int iDisplay,
226  const Eigen::MatrixXd& matData,
227  const Eigen::MatrixXd& matProjectors,
228  const SensorSet& sensors,
229  int &iSimplexNumitr)
230 {
231  double tolx, tolf, rho, chi, psi, sigma, func_evals, usual_delta, zero_term_delta, temp1, temp2;
232  std::string header, how;
233  int n, itercount, prnt;
234  Eigen::MatrixXd onesn, two2np1, one2n, v, y, v1, tempX1, tempX2, xbar, xr, x, xe, xc, xcc, xin,posCopy;
235  std::vector <double> fv, fv1;
236  std::vector <int> idx;
237 
238  DipFitError tempdip, fxr, fxe, fxc, fxcc;
239 
240  tolx = tolf = m_fAbortError;
241 
242  switch(iDisplay) {
243  case 0:
244  prnt = 0;
245  break;
246  default:
247  prnt = 1;
248  }
249 
250  header = " Iteration Func-count min f(x) Procedure";
251 
252  posCopy = matPos;
253 
254  n = posCopy.cols();
255 
256  // Initialize parameters
257  rho = 1; chi = 2; psi = 0.5; sigma = 0.5;
258  onesn = Eigen::MatrixXd::Ones(1,n);
259  two2np1 = one2n = Eigen::MatrixXd::Zero(1,n);
260 
261  for(int i = 0;i < n;i++) {
262  two2np1(i) = 1 + i;
263  one2n(i) = i;
264  }
265 
266  v = v1 = Eigen::MatrixXd::Zero(n, n+1);
267  fv.resize(n+1);
268  idx.resize(n+1);
269  fv1.resize(n+1);
270 
271  for(int i = 0;i < n; i++) {
272  v(i,0) = posCopy(i);
273  }
274 
275  tempdip = dipfitError(posCopy, matData, sensors, matProjectors);
276  fv[0] = tempdip.error;
277 
278  func_evals = 1;
279  itercount = 0;
280  how = "";
281 
282  // Continue setting up the initial simplex.
283  // Following improvement suggested by L.Pfeffer at Stanford
284  usual_delta = 0.05; // 5 percent deltas for non-zero terms
285  zero_term_delta = 0.00025; // Even smaller delta for zero elements of x
286  xin = posCopy.transpose();
287 
288  for(int j = 0;j < n;j++) {
289  y = xin;
290 
291  if(y(j) != 0) {
292  y(j) = (1 + usual_delta) * y(j);
293  } else {
294  y(j) = zero_term_delta;
295  }
296 
297  v.col(j+1).array() = y;
298  posCopy = y.transpose();
299  tempdip = dipfitError(posCopy, matData, sensors, matProjectors);
300  fv[j+1] = tempdip.error;
301  }
302 
303  // Sort elements of fv
304  std::vector<HPISortStruct> vecSortStruct;
305 
306  for (int i = 0; i < fv.size(); i++) {
307  HPISortStruct structTemp;
308  structTemp.base_arr = fv[i];
309  structTemp.idx = i;
310  vecSortStruct.push_back(structTemp);
311  }
312 
313  std::sort(vecSortStruct.begin(), vecSortStruct.end(), compare);
314 
315  for (int i = 0; i < vecSortStruct.size(); i++) {
316  idx[i] = vecSortStruct[i].idx;
317  }
318 
319  for (int i = 0;i < n+1;i++) {
320  v1.col(i) = v.col(idx[i]);
321  fv1[i] = fv[idx[i]];
322  }
323 
324  v = v1;fv = fv1;
325 
326  how = "initial simplex";
327  itercount = itercount + 1;
328  func_evals = n + 1;
329 
330  tempX1 = Eigen::MatrixXd::Zero(1,n);
331 
332  while ((func_evals < iMaxfun) && (itercount < iMaxiter)) {
333 
334  for (int i = 0;i < n;i++) {
335  tempX1(i) = std::fabs(fv[0] - fv[i+1]);
336  }
337 
338  temp1 = tempX1.maxCoeff();
339 
340  tempX2 = Eigen::MatrixXd::Zero(n,n);
341 
342  for(int i = 0;i < n;i++) {
343  tempX2.col(i) = v.col(i+1) - v.col(0);
344  }
345 
346  tempX2 = tempX2.array().abs();
347 
348  temp2 = tempX2.maxCoeff();
349 
350  if((temp1 <= tolf) && (temp2 <= tolx)) {
351  break;
352  }
353 
354  xbar = v.block(0,0,n,n).rowwise().sum();
355  xbar /= n;
356 
357  xr = (1+rho) * xbar - rho * v.block(0,n,v.rows(),1);
358 
359  x = xr.transpose();
360  //std::cout << "Iteration Count: " << itercount << ":" << x << std::endl;
361 
362  fxr = dipfitError(x, matData, sensors, matProjectors);
363 
364  func_evals = func_evals+1;
365 
366  if (fxr.error < fv[0]) {
367  // Calculate the expansion point
368  xe = (1 + rho * chi) * xbar - rho * chi * v.col(v.cols()-1);
369  x = xe.transpose();
370  fxe = dipfitError(x, matData, sensors, matProjectors);
371  func_evals = func_evals+1;
372 
373  if(fxe.error < fxr.error) {
374  v.col(v.cols()-1) = xe;
375  fv[n] = fxe.error;
376  how = "expand";
377  } else {
378  v.col(v.cols()-1) = xr;
379  fv[n] = fxr.error;
380  how = "reflect";
381  }
382  }
383  else {
384  if(fxr.error < fv[n-1]) {
385  v.col(v.cols()-1) = xr;
386  fv[n] = fxr.error;
387  how = "reflect";
388  } else { // fxr.error >= fv[:,n-1]
389  // Perform contraction
390  if(fxr.error < fv[n]) {
391  // Perform an outside contraction
392  xc = (1 + psi * rho) * xbar - psi * rho * v.col(v.cols()-1);
393  x = xc.transpose();
394  fxc = dipfitError(x, matData, sensors, matProjectors);
395  func_evals = func_evals + 1;
396 
397  if(fxc.error <= fxr.error) {
398  v.col(v.cols()-1) = xc;
399  fv[n] = fxc.error;
400  how = "contract outside";
401  } else {
402  // perform a shrink
403  how = "shrink";
404  }
405  } else {
406  xcc = (1 - psi) * xbar + psi * v.col(v.cols()-1);
407  x = xcc.transpose();
408  fxcc = dipfitError(x, matData, sensors, matProjectors);
409  func_evals = func_evals+1;
410  if(fxcc.error < fv[n]) {
411  v.col(v.cols()-1) = xcc;
412  fv[n] = fxcc.error;
413  how = "contract inside";
414  } else {
415  // perform a shrink
416  how = "shrink";
417  }
418  }
419 
420  if(how.compare("shrink") == 0) {
421  for(int j = 1;j < n+1;j++) {
422  v.col(j).array() = v.col(0).array() + sigma * (v.col(j).array() - v.col(0).array());
423  x = v.col(j).array().transpose();
424  tempdip = dipfitError(x,matData, sensors, matProjectors);
425  fv[j] = tempdip.error;
426  }
427  }
428  }
429  }
430 
431  // Sort elements of fv
432  vecSortStruct.clear();
433 
434  for (int i = 0; i < fv.size(); i++) {
435  HPISortStruct structTemp;
436  structTemp.base_arr = fv[i];
437  structTemp.idx = i;
438  vecSortStruct.push_back(structTemp);
439  }
440 
441  std::sort(vecSortStruct.begin(), vecSortStruct.end(), compare);
442  for (int i = 0; i < vecSortStruct.size(); i++) {
443  idx[i] = vecSortStruct[i].idx;
444  }
445 
446  for (int i = 0;i < n+1;i++) {
447  v1.col(i) = v.col(idx[i]);
448  fv1[i] = fv[idx[i]];
449  }
450 
451  v = v1;
452  fv = fv1;
453  itercount = itercount + 1;
454  }
455 
456  x = v.col(0).transpose();
457 
458  // Seok
459  iSimplexNumitr = itercount;
460 
461  return x;
462 }
463 
Eigen::MatrixXd fminsearch(const Eigen::MatrixXd &matPos, int iMaxiter, int iMaxfun, int iDisplay, const Eigen::MatrixXd &matData, const Eigen::MatrixXd &matProjectors, const SensorSet &sensors, int &iSimplexNumitr)
Definition: hpifitdata.cpp:222
Eigen::MatrixXd compute_leadfield(const Eigen::MatrixXd &matPos, const SensorSet &sensors)
Definition: hpifitdata.cpp:166
DipFitError dipfitError(const Eigen::MatrixXd &matPos, const Eigen::MatrixXd &matData, const SensorSet &sensors, const Eigen::MatrixXd &matProjectors)
Definition: hpifitdata.cpp:180
Eigen::MatrixXd magnetic_dipole(Eigen::MatrixXd matPos, Eigen::MatrixXd matPnt, Eigen::MatrixXd matOri)
Definition: hpifitdata.cpp:111
MNEMath class declaration.
HPIFit class declaration.
static bool compare(HPISortStruct a, HPISortStruct b)
Definition: hpifitdata.cpp:215
SensorSet class declaration.
HPIFitData class declaration.
static Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > pinv(const Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &a)
Definition: mnemath.h:755