MNE-CPP  0.1.9
A Framework for Electrophysiology
rapmusic.cpp
Go to the documentation of this file.
1 //=============================================================================================================
36 //=============================================================================================================
37 // INCLUDES
38 //=============================================================================================================
39 
40 #include "rapmusic.h"
41 
42 #include <utils/mnemath.h>
43 
44 #ifdef _OPENMP
45 #include <omp.h>
46 #endif
47 
48 //=============================================================================================================
49 // USED NAMESPACES
50 //=============================================================================================================
51 
52 using namespace INVERSELIB;
53 using namespace MNELIB;
54 using namespace FIFFLIB;
55 using namespace UTILSLIB;
56 
57 //=============================================================================================================
58 // DEFINE MEMBER METHODS
59 //=============================================================================================================
60 
62 : m_iN(0)
63 , m_dThreshold(0)
64 , m_iNumGridPoints(0)
65 , m_iNumChannels(0)
66 , m_iNumLeadFieldCombinations(0)
67 , m_ppPairIdxCombinations(NULL)
68 , m_iMaxNumThreads(1)
69 , m_bIsInit(false)
70 , m_iSamplesStcWindow(-1)
71 , m_fStcOverlap(-1)
72 {
73 }
74 
75 //=============================================================================================================
76 
77 RapMusic::RapMusic(MNEForwardSolution& p_pFwd, bool p_bSparsed, int p_iN, double p_dThr)
78 : m_iN(0)
79 , m_dThreshold(0)
81 , m_iNumChannels(0)
85 , m_bIsInit(false)
87 , m_fStcOverlap(-1)
88 {
89  //Init
90  init(p_pFwd, p_bSparsed, p_iN, p_dThr);
91 }
92 
93 //=============================================================================================================
94 
95 RapMusic::~RapMusic()
96 {
97  if(m_ppPairIdxCombinations != NULL)
99 }
100 
101 //=============================================================================================================
102 
103 bool RapMusic::init(MNEForwardSolution& p_pFwd, bool p_bSparsed, int p_iN, double p_dThr)
104 {
105  //Get available thread number
106  #ifdef _OPENMP
107  std::cout << "OpenMP enabled" << std::endl;
108  m_iMaxNumThreads = omp_get_max_threads();
109  #else
110  std::cout << "OpenMP disabled (to enable it: VS2010->Project Properties->C/C++->Language, then modify OpenMP Support)" << std::endl;
111  m_iMaxNumThreads = 1;
112  #endif
113  std::cout << "Available Threats: " << m_iMaxNumThreads << std::endl << std::endl;
114 
115  //Initialize RAP MUSIC
116  std::cout << "##### Initialization RAP MUSIC started ######\n\n";
117 
118  m_iN = p_iN;
119  m_dThreshold = p_dThr;
120 
121 // //Grid check
122 // if(p_pMatGrid != NULL)
123 // {
124 // if ( p_pMatGrid->rows() != p_pMatLeadField->cols() / 3 )
125 // {
126 // std::cout << "Grid does not fit to given Lead Field!\n";
127 // return false;
128 // }
129 // }
130 
131 // m_pMatGrid = p_pMatGrid;
132 
133  //Lead Field check
134  if ( p_pFwd.sol->data.cols() % 3 != 0 )
135  {
136  std::cout << "Gain matrix is not associated with a 3D grid!\n";
137  return false;
138  }
139 
140  m_iNumGridPoints = p_pFwd.sol->data.cols()/3;
141 
142  m_iNumChannels = p_pFwd.sol->data.rows();
143 
144 // m_pMappedMatLeadField = new Eigen::Map<MatrixXT>
145 // ( p_pMatLeadField->data(),
146 // p_pMatLeadField->rows(),
147 // p_pMatLeadField->cols() );
148 
149  m_ForwardSolution = p_pFwd;
150 
151  //##### Calc lead field combination #####
152 
153  std::cout << "Calculate gain matrix combinations. \n";
154 
155  m_iNumLeadFieldCombinations = MNEMath::nchoose2(m_iNumGridPoints+1);
156 
158 
160 
161  std::cout << "Gain matrix combinations calculated. \n\n";
162 
163  //##### Calc lead field combination end #####
164 
165  std::cout << "Number of grid points: " << m_iNumGridPoints << "\n\n";
166 
167  std::cout << "Number of combinated points: " << m_iNumLeadFieldCombinations << "\n\n";
168 
169  std::cout << "Number of sources to find: " << m_iN << "\n\n";
170 
171  std::cout << "Threshold: " << m_dThreshold << "\n\n";
172 
173  //Init end
174 
175  std::cout << "##### Initialization RAP MUSIC completed ######\n\n\n";
176 
177  Q_UNUSED(p_bSparsed);
178 
179  m_bIsInit = true;
180 
181  return m_bIsInit;
182 }
183 
184 //=============================================================================================================
185 
186 const char* RapMusic::getName() const
187 {
188  return "RAP MUSIC";
189 }
190 
191 //=============================================================================================================
192 
194 {
195  return m_ForwardSolution.src;
196 }
197 
198 //=============================================================================================================
199 
200 MNESourceEstimate RapMusic::calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal)
201 {
202  Q_UNUSED(pick_normal);
203 
204  MNESourceEstimate p_sourceEstimate;
205 
206  if(p_fiffEvoked.data.rows() != m_iNumChannels)
207  {
208  std::cout << "Number of FiffEvoked channels (" << p_fiffEvoked.data.rows() << ") doesn't match the number of channels (" << m_iNumChannels << ") of the forward solution." << std::endl;
209  return p_sourceEstimate;
210  }
211 // else
212 // std::cout << "Number of FiffEvoked channels (" << p_fiffEvoked.data.rows() << ") matchs the number of channels (" << m_iNumChannels << ") of the forward solution." << std::endl;
213 
214  //
215  // Rap MUSIC Source estimate
216  //
217  p_sourceEstimate.data = MatrixXd::Zero(m_ForwardSolution.nsource, p_fiffEvoked.data.cols());
218 
219  //Results
220  p_sourceEstimate.vertices = VectorXi(m_ForwardSolution.src[0].vertno.size() + m_ForwardSolution.src[1].vertno.size());
221  p_sourceEstimate.vertices << m_ForwardSolution.src[0].vertno, m_ForwardSolution.src[1].vertno;
222 
223  p_sourceEstimate.times = p_fiffEvoked.times;
224  p_sourceEstimate.tmin = p_fiffEvoked.times[0];
225  p_sourceEstimate.tstep = p_fiffEvoked.times[1] - p_fiffEvoked.times[0];
226 
227  if(m_iSamplesStcWindow <= 3) //if samples per stc aren't set -> use full window
228  {
229  QList< DipolePair<double> > t_RapDipoles;
230  calculateInverse(p_fiffEvoked.data, t_RapDipoles);
231 
232  for(qint32 i = 0; i < t_RapDipoles.size(); ++i)
233  {
234  double dip1 = sqrt( pow(t_RapDipoles[i].m_Dipole1.phi_x(),2) +
235  pow(t_RapDipoles[i].m_Dipole1.phi_y(),2) +
236  pow(t_RapDipoles[i].m_Dipole1.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
237 
238  double dip2 = sqrt( pow(t_RapDipoles[i].m_Dipole2.phi_x(),2) +
239  pow(t_RapDipoles[i].m_Dipole2.phi_y(),2) +
240  pow(t_RapDipoles[i].m_Dipole2.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
241 
242  RowVectorXd dip1Time = RowVectorXd::Constant(p_fiffEvoked.data.cols(), dip1);
243  RowVectorXd dip2Time = RowVectorXd::Constant(p_fiffEvoked.data.cols(), dip2);
244 
245  p_sourceEstimate.data.block(t_RapDipoles[i].m_iIdx1, 0, 1, p_fiffEvoked.data.cols()) = dip1Time;
246  p_sourceEstimate.data.block(t_RapDipoles[i].m_iIdx2, 0, 1, p_fiffEvoked.data.cols()) = dip2Time;
247  }
248  }
249  else
250  {
251  bool first = true;
252  bool last = false;
253 
254  qint32 t_iNumSensors = p_fiffEvoked.data.rows();
255  qint32 t_iNumSteps = p_fiffEvoked.data.cols();
256 
257  qint32 t_iSamplesOverlap = (qint32)floor(((float)m_iSamplesStcWindow)*m_fStcOverlap);
258  qint32 t_iSamplesDiscard = t_iSamplesOverlap/2;
259 
260  MatrixXd data = MatrixXd::Zero(t_iNumSensors, m_iSamplesStcWindow);
261 
262  qint32 curSample = 0;
263  qint32 curResultSample = 0;
264  qint32 stcWindowSize = m_iSamplesStcWindow - 2*t_iSamplesDiscard;
265 
266  while(!last)
267  {
268  QList< DipolePair<double> > t_RapDipoles;
269 
270  //Data
271  if(curSample + m_iSamplesStcWindow >= t_iNumSteps) //last
272  {
273  last = true;
274  data = p_fiffEvoked.data.block(0, p_fiffEvoked.data.cols()-m_iSamplesStcWindow, t_iNumSensors, m_iSamplesStcWindow);
275  }
276  else
277  data = p_fiffEvoked.data.block(0, curSample, t_iNumSensors, m_iSamplesStcWindow);
278 
279  curSample += (m_iSamplesStcWindow - t_iSamplesOverlap);
280  if(first)
281  curSample -= t_iSamplesDiscard; //shift on start t_iSamplesDiscard backwards
282 
283  //Calculate
284  calculateInverse(data, t_RapDipoles);
285 
286  //Assign Result
287  if(last)
288  stcWindowSize = p_sourceEstimate.data.cols() - curResultSample;
289 
290  for(qint32 i = 0; i < t_RapDipoles.size(); ++i)
291  {
292  double dip1 = sqrt( pow(t_RapDipoles[i].m_Dipole1.phi_x(),2) +
293  pow(t_RapDipoles[i].m_Dipole1.phi_y(),2) +
294  pow(t_RapDipoles[i].m_Dipole1.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
295 
296  double dip2 = sqrt( pow(t_RapDipoles[i].m_Dipole2.phi_x(),2) +
297  pow(t_RapDipoles[i].m_Dipole2.phi_y(),2) +
298  pow(t_RapDipoles[i].m_Dipole2.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
299 
300  RowVectorXd dip1Time = RowVectorXd::Constant(stcWindowSize, dip1);
301  RowVectorXd dip2Time = RowVectorXd::Constant(stcWindowSize, dip2);
302 
303  p_sourceEstimate.data.block(t_RapDipoles[i].m_iIdx1, curResultSample, 1, stcWindowSize) = dip1Time;
304  p_sourceEstimate.data.block(t_RapDipoles[i].m_iIdx2, curResultSample, 1, stcWindowSize) = dip2Time;
305 
306  }
307 
308  curResultSample += stcWindowSize;
309 
310  if(first)
311  first = false;
312  }
313  }
314 
315  return p_sourceEstimate;
316 }
317 
318 //=============================================================================================================
319 
320 MNESourceEstimate RapMusic::calculateInverse(const MatrixXd &data, float tmin, float tstep, bool pick_normal) const
321 {
322  Q_UNUSED(pick_normal);
323 
324  MNESourceEstimate p_sourceEstimate;
325 
326  if(data.rows() != m_iNumChannels)
327  {
328  std::cout << "Number of FiffEvoked channels (" << data.rows() << ") doesn't match the number of channels (" << m_iNumChannels << ") of the forward solution." << std::endl;
329  return p_sourceEstimate;
330  }
331 // else
332 // std::cout << "Number of FiffEvoked channels (" << data.rows() << ") matchs the number of channels (" << m_iNumChannels << ") of the forward solution." << std::endl;
333 
334  //
335  // Rap MUSIC Source estimate
336  //
337  p_sourceEstimate.data = MatrixXd::Zero(m_ForwardSolution.nsource, data.cols());
338 
339  //Results
340  p_sourceEstimate.vertices = VectorXi(m_ForwardSolution.src[0].vertno.size() + m_ForwardSolution.src[1].vertno.size());
341  p_sourceEstimate.vertices << m_ForwardSolution.src[0].vertno, m_ForwardSolution.src[1].vertno;
342 
343  p_sourceEstimate.times = RowVectorXf::Zero(data.cols());
344  p_sourceEstimate.times[0] = tmin;
345  for(qint32 i = 1; i < p_sourceEstimate.times.size(); ++i)
346  p_sourceEstimate.times[i] = p_sourceEstimate.times[i-1] + tstep;
347  p_sourceEstimate.tmin = tmin;
348  p_sourceEstimate.tstep = tstep;
349 
350  QList< DipolePair<double> > t_RapDipoles;
351  calculateInverse(data, t_RapDipoles);
352 
353  for(qint32 i = 0; i < t_RapDipoles.size(); ++i)
354  {
355  double dip1 = sqrt( pow(t_RapDipoles[i].m_Dipole1.phi_x(),2) +
356  pow(t_RapDipoles[i].m_Dipole1.phi_y(),2) +
357  pow(t_RapDipoles[i].m_Dipole1.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
358 
359  double dip2 = sqrt( pow(t_RapDipoles[i].m_Dipole2.phi_x(),2) +
360  pow(t_RapDipoles[i].m_Dipole2.phi_y(),2) +
361  pow(t_RapDipoles[i].m_Dipole2.phi_z(),2) ) * t_RapDipoles[i].m_vCorrelation;
362 
363  RowVectorXd dip1Time = RowVectorXd::Constant(data.cols(), dip1);
364  RowVectorXd dip2Time = RowVectorXd::Constant(data.cols(), dip2);
365 
366  p_sourceEstimate.data.block(t_RapDipoles[i].m_iIdx1, 0, 1, data.cols()) = dip1Time;
367  p_sourceEstimate.data.block(t_RapDipoles[i].m_iIdx2, 0, 1, data.cols()) = dip2Time;
368  }
369 
370  return p_sourceEstimate;
371 }
372 
373 //=============================================================================================================
374 
375 MNESourceEstimate RapMusic::calculateInverse(const MatrixXd& p_matMeasurement, QList< DipolePair<double> > &p_RapDipoles) const
376 {
377  MNESourceEstimate p_SourceEstimate;
378 
379  //if not initialized -> break
380  if(!m_bIsInit)
381  {
382  std::cout << "RAP MUSIC wasn't initialized!"; //ToDo: catch this earlier
383  return p_SourceEstimate; //false
384  }
385 
386  //Test if data are correct
387  if(p_matMeasurement.rows() != m_iNumChannels)
388  {
389  std::cout << "Lead Field channels do not fit to number of measurement channels!"; //ToDo: catch this earlier
390  return p_SourceEstimate;
391  }
392 
393  //Inits
394  //Stop the time for benchmark purpose
395  clock_t start, end;
396  start = clock();
397 
398 // //Map HPCMatrix to Eigen Matrix
399 // Eigen::Map<MatrixXT>
400 // t_MappedMatMeasurement( p_pMatMeasurement->data(),
401 // p_pMatMeasurement->rows(),
402 // p_pMatMeasurement->cols() );
403 
404  //Calculate the signal subspace (t_pMatPhi_s)
405  MatrixXT* t_pMatPhi_s = NULL;//(m_iNumChannels, m_iN < t_r ? m_iN : t_r);
406  int t_r = calcPhi_s(/*(MatrixXT)*/p_matMeasurement, t_pMatPhi_s);
407 
408  int t_iMaxSearch = m_iN < t_r ? m_iN : t_r; //The smallest of Rank and Iterations
409 
410  if (t_r < m_iN)
411  {
412  std::cout << "Warning: Rank " << t_r << " of the measurement data is smaller than the " << m_iN;
413  std::cout << " sources to find." << std::endl;
414  std::cout << " Searching now for " << t_iMaxSearch << " correlated sources.";
415  std::cout << std::endl << std::endl;
416  }
417 
418  //Create Orthogonal Projector
419  //OrthProj
420  MatrixXT t_matOrthProj(m_iNumChannels,m_iNumChannels);
421  t_matOrthProj.setIdentity();
422 
423  //A_k_1
424  MatrixXT t_matA_k_1(m_iNumChannels, t_iMaxSearch);
425  t_matA_k_1.setZero();
426 
427 // if (m_pMatGrid != NULL)
428 // {
429 // if(p_pRapDipoles != NULL)
430 // p_pRapDipoles->initRapDipoles(m_pMatGrid);
431 // else
432 // p_pRapDipoles = new RapDipoles<T>(m_pMatGrid);
433 // }
434 // else
435 // {
436 // if(p_pRapDipoles != NULL)
437 // delete p_pRapDipoles;
438 
439 // p_pRapDipoles = new RapDipoles<T>();
440 // }
441  p_RapDipoles.clear();
442 
443  std::cout << "##### Calculation of RAP MUSIC started ######\n\n";
444 
445  MatrixXT t_matProj_Phi_s(t_matOrthProj.rows(), t_pMatPhi_s->cols());
446  //new Version: Calculate projection before
447  MatrixXT t_matProj_LeadField(m_ForwardSolution.sol->data.rows(), m_ForwardSolution.sol->data.cols());
448 
449  for(int r = 0; r < t_iMaxSearch ; ++r)
450  {
451  t_matProj_Phi_s = t_matOrthProj*(*t_pMatPhi_s);
452 
453  //new Version: Calculating Projection before
454  t_matProj_LeadField = t_matOrthProj * m_ForwardSolution.sol->data;//Subtract the found sources from the current found source
455 
456  //###First Option###
457  //Step 1: lt. Mosher 1998 -> Maybe tmp_Proj_Phi_S is already orthogonal -> so no SVD needed -> U_B = tmp_Proj_Phi_S;
458  Eigen::JacobiSVD< MatrixXT > t_svdProj_Phi_S(t_matProj_Phi_s, Eigen::ComputeThinU);
459  MatrixXT t_matU_B;
460  useFullRank(t_svdProj_Phi_S.matrixU(), t_svdProj_Phi_S.singularValues().asDiagonal(), t_matU_B);
461 
462  //Inits
464  t_vecRoh.setZero();
465 
466  //subcorr benchmark
467  //Stop the time
468  clock_t start_subcorr, end_subcorr;
469  start_subcorr = clock();
470 
471  //Multithreading correlation calculation
472  #ifdef _OPENMP
473  #pragma omp parallel num_threads(m_iMaxNumThreads)
474  #endif
475  {
476  #ifdef _OPENMP
477  #pragma omp for
478  #endif
479  for(int i = 0; i < m_iNumLeadFieldCombinations; i++)
480  {
481  //new Version: calculate matrix multiplication before
482  //Create Lead Field combinations -> It would be better to use a pointer construction, to increase performance
483  MatrixX6T t_matProj_G(t_matProj_LeadField.rows(),6);
484 
485  int idx1 = m_ppPairIdxCombinations[i]->x1;
486  int idx2 = m_ppPairIdxCombinations[i]->x2;
487 
488  RapMusic::getGainMatrixPair(t_matProj_LeadField, t_matProj_G, idx1, idx2);
489 
490  t_vecRoh(i) = RapMusic::subcorr(t_matProj_G, t_matU_B);//t_vecRoh holds the correlations roh_k
491  }
492  }
493 
494 // if(r==0)
495 // {
496 // std::fstream filestr;
497 // std::stringstream filename;
498 // filename << "Roh_gold.txt";
499 //
500 // filestr.open ( filename.str().c_str(), std::fstream::out);
501 // for(int i = 0; i < m_iNumLeadFieldCombinations; ++i)
502 // {
503 // filestr << t_vecRoh(i) << "\n";
504 // }
505 // filestr.close();
506 //
507 // //exit(0);
508 // }
509 
510  //subcorr benchmark
511  end_subcorr = clock();
512 
513  float t_fSubcorrElapsedTime = ( (float)(end_subcorr-start_subcorr) / (float)CLOCKS_PER_SEC ) * 1000.0f;
514  std::cout << "Time Elapsed: " << t_fSubcorrElapsedTime << " ms" << std::endl;
515 
516  //Find the maximum of correlation - can't put this in the for loop because it's running in different threads.
517  double t_val_roh_k;
518 
519  VectorXT::Index t_iMaxIdx;
520 
521  t_val_roh_k = t_vecRoh.maxCoeff(&t_iMaxIdx);//p_vecCor = ^roh_k
522 
523  //get positions in sparsed leadfield from index combinations;
524  int t_iIdx1 = m_ppPairIdxCombinations[t_iMaxIdx]->x1;
525  int t_iIdx2 = m_ppPairIdxCombinations[t_iMaxIdx]->x2;
526 
527  // (Idx+1) because of MATLAB positions -> starting with 1 not with 0
528  std::cout << "Iteration: " << r+1 << " of " << t_iMaxSearch
529  << "; Correlation: " << t_val_roh_k<< "; Position (Idx+1): " << t_iIdx1+1 << " - " << t_iIdx2+1 <<"\n\n";
530 
531  //Calculations with the max correlated dipole pair G_k_1 -> ToDo Obsolet when taking direkt Projected Lead Field
532  MatrixX6T t_matG_k_1(m_ForwardSolution.sol->data.rows(),6);
533  RapMusic::getGainMatrixPair(m_ForwardSolution.sol->data, t_matG_k_1, t_iIdx1, t_iIdx2);
534 
535  MatrixX6T t_matProj_G_k_1(t_matOrthProj.rows(), t_matG_k_1.cols());
536  t_matProj_G_k_1 = t_matOrthProj * t_matG_k_1;//Subtract the found sources from the current found source
537 // MatrixX6T t_matProj_G_k_1(t_matProj_LeadField.rows(), 6);
538 // getLeadFieldPair(t_matProj_LeadField, t_matProj_G_k_1, t_iIdx1, t_iIdx2);
539 
540  //Calculate source direction
541  //source direction (p_pMatPhi) for current source r (phi_k_1)
542  Vector6T t_vec_phi_k_1(6);
543  RapMusic::subcorr(t_matProj_G_k_1, t_matU_B, t_vec_phi_k_1);//Correlate the current source to calculate the direction
544 
545  //Set return values
546  RapMusic::insertSource(t_iIdx1, t_iIdx2, t_vec_phi_k_1, t_val_roh_k, p_RapDipoles);
547 
548  //Stop Searching when Correlation is smaller then the Threshold
549  if (t_val_roh_k < m_dThreshold)
550  {
551  std::cout << "Searching stopped, last correlation " << t_val_roh_k;
552  std::cout << " is smaller then the given threshold " << m_dThreshold << std::endl << std::endl;
553  break;
554  }
555 
556  //Calculate A_k_1 = [a_theta_1..a_theta_k_1] matrix for subtraction of found source
557  RapMusic::calcA_k_1(t_matG_k_1, t_vec_phi_k_1, r, t_matA_k_1);
558 
559  //Calculate new orthogonal Projector (Pi_k_1)
560  calcOrthProj(t_matA_k_1, t_matOrthProj);
561 
562  //garbage collecting
563  //ToDo
564  }
565 
566  std::cout << "##### Calculation of RAP MUSIC completed ######"<< std::endl << std::endl << std::endl;
567 
568  end = clock();
569 
570  float t_fElapsedTime = ( (float)(end-start) / (float)CLOCKS_PER_SEC ) * 1000.0f;
571  std::cout << "Total Time Elapsed: " << t_fElapsedTime << " ms" << std::endl << std::endl;
572 
573  //garbage collecting
574  delete t_pMatPhi_s;
575 
576  return p_SourceEstimate;
577 }
578 
579 //=============================================================================================================
580 
581 int RapMusic::calcPhi_s(const MatrixXT& p_matMeasurement, MatrixXT* &p_pMatPhi_s) const
582 {
583  //Calculate p_pMatPhi_s
584  MatrixXT t_matF;//t_matF = makeSquareMat(p_pMatMeasurement); //FF^T -> ToDo Check this
585  if (p_matMeasurement.cols() > p_matMeasurement.rows())
586  t_matF = makeSquareMat(p_matMeasurement); //FF^T
587  else
588  t_matF = MatrixXT(p_matMeasurement);
589 
590  Eigen::JacobiSVD<MatrixXT> t_svdF(t_matF, Eigen::ComputeThinU);
591 
592  int t_r = getRank(t_svdF.singularValues().asDiagonal());
593 
594  int t_iCols = t_r;//t_r < m_iN ? m_iN : t_r;
595 
596  if (p_pMatPhi_s != NULL)
597  delete p_pMatPhi_s;
598 
599  //m_iNumChannels has to be equal to t_svdF.matrixU().rows()
600  p_pMatPhi_s = new MatrixXT(m_iNumChannels, t_iCols);
601 
602  //assign the signal subspace
603  memcpy(p_pMatPhi_s->data(), t_svdF.matrixU().data(), sizeof(double) * m_iNumChannels * t_iCols);
604 
605  return t_r;
606 }
607 
608 //=============================================================================================================
609 
610 double RapMusic::subcorr(MatrixX6T& p_matProj_G, const MatrixXT& p_matU_B)
611 {
612  //Orthogonalisierungstest wegen performance weggelassen -> ohne is es viel schneller
613 
614  Matrix6T t_matSigma_A(6, 6);
615  Matrix6XT t_matU_A_T(6, p_matProj_G.rows()); //rows and cols are changed, because of CV_SVD_U_T
616 
617  Eigen::JacobiSVD<MatrixXT> t_svdProj_G(p_matProj_G, Eigen::ComputeThinU);
618 
619  t_matSigma_A = t_svdProj_G.singularValues().asDiagonal();
620  t_matU_A_T = t_svdProj_G.matrixU().transpose();
621 
622  //lt. Mosher 1998 ToDo: Only Retain those Components of U_A and U_B that correspond to nonzero singular values
623  //for U_A and U_B the number of columns corresponds to their ranks
624  MatrixXT t_matU_A_T_full;
625  //reduce to rank only when directions aren't calculated, otherwise use the full t_matU_A_T
626 
627  useFullRank(t_matU_A_T, t_matSigma_A, t_matU_A_T_full, IS_TRANSPOSED);//lt. Mosher 1998: Only Retain those Components of U_A that correspond to nonzero singular values -> for U_A the number of columns corresponds to their ranks
628 
629  MatrixXT t_matCor(t_matU_A_T_full.rows(), p_matU_B.cols());
630 
631  //Step 2: compute the subspace correlation
632  t_matCor = t_matU_A_T_full*p_matU_B;//lt. Mosher 1998: C = U_A^T * U_B
633 
634  VectorXT t_vecSigma_C;
635 
636  if (t_matCor.cols() > t_matCor.rows())
637  {
638  MatrixXT t_matCor_H(t_matCor.cols(), t_matCor.rows());
639  t_matCor_H = t_matCor.adjoint(); //for complex it has to be adjunct
640  //ToDo -> use instead adjointInPlace
641 
642  Eigen::JacobiSVD<MatrixXT> t_svdCor_H(t_matCor_H);
643 
644  t_vecSigma_C = t_svdCor_H.singularValues();
645  }
646  else
647  {
648  Eigen::JacobiSVD<MatrixXT> t_svdCor(t_matCor);
649 
650  t_vecSigma_C = t_svdCor.singularValues();
651  }
652 
653  //Step 3
654  double t_dRetSigma_C;
655  t_dRetSigma_C = t_vecSigma_C(0); //Take only the correlation of the first principal components
656 
657  //garbage collecting
658  //ToDo
659 
660  return t_dRetSigma_C;
661 }
662 
663 //=============================================================================================================
664 
665 double RapMusic::subcorr(MatrixX6T& p_matProj_G, const MatrixXT& p_matU_B, Vector6T& p_vec_phi_k_1)
666 {
667  //Orthogonalisierungstest wegen performance weggelassen -> ohne is es viel schneller
668 
669  Matrix6T sigma_A(6, 6);
670  Matrix6XT U_A_T(6, p_matProj_G.rows()); //rows and cols are changed, because of CV_SVD_U_T
671  Matrix6T V_A(6, 6);
672 
673  Eigen::JacobiSVD<MatrixXT> svdOfProj_G(p_matProj_G, Eigen::ComputeThinU | Eigen::ComputeThinV);
674 
675  sigma_A = svdOfProj_G.singularValues().asDiagonal();
676  U_A_T = svdOfProj_G.matrixU().transpose();
677  V_A = svdOfProj_G.matrixV();
678 
679  //lt. Mosher 1998 ToDo: Only Retain those Components of U_A and U_B that correspond to nonzero singular values
680  //for U_A and U_B the number of columns corresponds to their ranks
681  //-> reduce to rank only when directions aren't calculated, otherwise use the full U_A_T
682 
683  Matrix6XT t_matCor(6, p_matU_B.cols());
684 
685  //Step 2: compute the subspace correlation
686  t_matCor = U_A_T*p_matU_B;//lt. Mosher 1998: C = U_A^T * U_B
687 
688  VectorXT sigma_C;
689 
690  //Step 4
691  Matrix6XT U_C;
692 
693  if (t_matCor.cols() > t_matCor.rows())
694  {
695  MatrixX6T Cor_H(t_matCor.cols(), 6);
696  Cor_H = t_matCor.adjoint(); //for complex it has to be adjunct
697 
698  Eigen::JacobiSVD<MatrixXT> svdOfCor_H(Cor_H, Eigen::ComputeThinV);
699 
700  U_C = svdOfCor_H.matrixV(); //because t_matCor Hermitesch U and V are exchanged
701  sigma_C = svdOfCor_H.singularValues();
702  }
703  else
704  {
705  Eigen::JacobiSVD<MatrixXT> svdOfCor(t_matCor, Eigen::ComputeThinU);
706 
707  U_C = svdOfCor.matrixU();
708  sigma_C = svdOfCor.singularValues();
709  }
710 
711  Matrix6T sigma_a_inv;
712  sigma_a_inv = sigma_A.inverse();
713 
714  Matrix6XT X;
715  X = (V_A*sigma_a_inv)*U_C;//X = V_A*Sigma_A^-1*U_C
716 
717  Vector6T X_max;//only for the maximum c - so instead of X->cols use 1
718  X_max = X.col(0);
719 
720  double norm_X = 1/(X_max.norm());
721 
722  //Multiply a scalar with an Array -> linear transform
723  p_vec_phi_k_1 = X_max*norm_X;//u1 = x1/||x1|| this is the orientation
724 
725  //garbage collecting
726  //ToDo
727 
728  //Step 3
729  double ret_sigma_C;
730  ret_sigma_C = sigma_C(0); //Take only the correlation of the first principal components
731 
732  //garbage collecting
733  //ToDo
734 
735  return ret_sigma_C;
736 }
737 
738 //=============================================================================================================
739 
740 void RapMusic::calcA_k_1( const MatrixX6T& p_matG_k_1,
741  const Vector6T& p_matPhi_k_1,
742  const int p_iIdxk_1,
743  MatrixXT& p_matA_k_1)
744 {
745  //Calculate A_k_1 = [a_theta_1..a_theta_k_1] matrix for subtraction of found source
746  VectorXT t_vec_a_theta_k_1(p_matG_k_1.rows(),1);
747 
748  t_vec_a_theta_k_1 = p_matG_k_1*p_matPhi_k_1; // a_theta_k_1 = G_k_1*phi_k_1 this corresponds to the normalized signal component in subspace r
749 
750  p_matA_k_1.block(0,p_iIdxk_1,p_matA_k_1.rows(),1) = t_vec_a_theta_k_1;
751 }
752 
753 //=============================================================================================================
754 
755 void RapMusic::calcOrthProj(const MatrixXT& p_matA_k_1, MatrixXT& p_matOrthProj) const
756 {
757  //Calculate OrthProj=I-A_k_1*(A_k_1'*A_k_1)^-1*A_k_1' //Wetterling -> A_k_1 = Gain
758 
759  MatrixXT t_matA_k_1_tmp(p_matA_k_1.cols(), p_matA_k_1.cols());
760  t_matA_k_1_tmp = p_matA_k_1.adjoint()*p_matA_k_1;//A_k_1'*A_k_1 = A_k_1_tmp -> A_k_1' has to be adjoint for complex
761 
762  int t_size = t_matA_k_1_tmp.cols();
763 
764  while (!t_matA_k_1_tmp.block(0,0,t_size,t_size).fullPivLu().isInvertible())
765  {
766  --t_size;
767  }
768 
769  MatrixXT t_matA_k_1_tmp_inv(t_matA_k_1_tmp.rows(), t_matA_k_1_tmp.cols());
770  t_matA_k_1_tmp_inv.setZero();
771 
772  t_matA_k_1_tmp_inv.block(0,0,t_size,t_size) = t_matA_k_1_tmp.block(0,0,t_size,t_size).inverse();//(A_k_1_tmp)^-1 = A_k_1_tmp_inv
773 
774  t_matA_k_1_tmp = MatrixXT::Zero(p_matA_k_1.rows(), p_matA_k_1.cols());
775 
776  t_matA_k_1_tmp = p_matA_k_1*t_matA_k_1_tmp_inv;//(A_k_1*A_k_1_tmp_inv) = A_k_1_tmp
777 
778  MatrixXT t_matA_k_1_tmp2(p_matA_k_1.rows(), p_matA_k_1.rows());
779 
780  t_matA_k_1_tmp2 = t_matA_k_1_tmp*p_matA_k_1.adjoint();//(A_k_1_tmp)*A_k_1' -> here A_k_1' is only transposed - it has to be adjoint
781 
783  I.setIdentity();
784 
785  p_matOrthProj = I-t_matA_k_1_tmp2; //OrthProj=I-A_k_1*(A_k_1'*A_k_1)^-1*A_k_1';
786 
787  //garbage collecting
788  //ToDo
789 }
790 
791 //=============================================================================================================
792 
793 void RapMusic::calcPairCombinations( const int p_iNumPoints,
794  const int p_iNumCombinations,
795  Pair** p_ppPairIdxCombinations) const
796 {
797  int idx1 = 0;
798  int idx2 = 0;
799 
800  //Process Code in {m_max_num_threads} threads -> When compile with Intel Compiler -> probably obsolete
801  #ifdef _OPENMP
802  #pragma omp parallel num_threads(m_iMaxNumThreads) private(idx1, idx2)
803  #endif
804  {
805  #ifdef _OPENMP
806  #pragma omp for
807  #endif
808  for (int i = 0; i < p_iNumCombinations; ++i)
809  {
810  RapMusic::getPointPair(p_iNumPoints, i, idx1, idx2);
811 
812  Pair* t_pairCombination = new Pair();
813  t_pairCombination->x1 = idx1;
814  t_pairCombination->x2 = idx2;
815 
816  p_ppPairIdxCombinations[i] = t_pairCombination;
817  }
818  }
819 }
820 
821 //=============================================================================================================
822 
823 void RapMusic::getPointPair(const int p_iPoints, const int p_iCurIdx, int &p_iIdx1, int &p_iIdx2)
824 {
825  int ii = p_iPoints*(p_iPoints+1)/2-1-p_iCurIdx;
826  int K = (int)floor((sqrt((double)(8*ii+1))-1)/2);
827 
828  p_iIdx1 = p_iPoints-1-K;
829 
830  p_iIdx2 = (p_iCurIdx-p_iPoints*(p_iPoints+1)/2 + (K+1)*(K+2)/2)+p_iIdx1;
831 }
832 
833 //=============================================================================================================
834 //ToDo don't make a real copy
835 void RapMusic::getGainMatrixPair( const MatrixXT& p_matGainMarix,
836  MatrixX6T& p_matGainMarix_Pair,
837  int p_iIdx1, int p_iIdx2)
838 {
839  p_matGainMarix_Pair.block(0,0,p_matGainMarix.rows(),3) = p_matGainMarix.block(0, p_iIdx1*3, p_matGainMarix.rows(), 3);
840 
841  p_matGainMarix_Pair.block(0,3,p_matGainMarix.rows(),3) = p_matGainMarix.block(0, p_iIdx2*3, p_matGainMarix.rows(), 3);
842 }
843 
844 //=============================================================================================================
845 
846 void RapMusic::insertSource( int p_iDipoleIdx1, int p_iDipoleIdx2,
847  const Vector6T &p_vec_phi_k_1,
848  double p_valCor,
849  QList< DipolePair<double> > &p_RapDipoles)
850 {
851  DipolePair<double> t_pRapDipolePair;
852 
853  t_pRapDipolePair.m_iIdx1 = p_iDipoleIdx1; //p_iDipoleIdx1+1 because of MATLAB index
854  t_pRapDipolePair.m_iIdx2 = p_iDipoleIdx2;
855 
856  t_pRapDipolePair.m_Dipole1.x() = 0; //m_bGridInitialized ? (*m_pMatGrid)(p_iDipoleIdx1, 0) : 0;
857  t_pRapDipolePair.m_Dipole1.y() = 0; //m_bGridInitialized ? (*m_pMatGrid)(p_iDipoleIdx1, 1) : 0;
858  t_pRapDipolePair.m_Dipole1.z() = 0; //m_bGridInitialized ? (*m_pMatGrid)(p_iDipoleIdx1, 2) : 0;
859 
860  t_pRapDipolePair.m_Dipole2.x() = 0; //m_bGridInitialized ? (*m_pMatGrid)(p_iDipoleIdx2, 0) : 0;
861  t_pRapDipolePair.m_Dipole2.y() = 0; //m_bGridInitialized ? (*m_pMatGrid)(p_iDipoleIdx2, 1) : 0;
862  t_pRapDipolePair.m_Dipole2.z() = 0; //m_bGridInitialized ? (*m_pMatGrid)(p_iDipoleIdx2, 2) : 0;
863 
864  t_pRapDipolePair.m_Dipole1.phi_x() = p_vec_phi_k_1[0];
865  t_pRapDipolePair.m_Dipole1.phi_y() = p_vec_phi_k_1[1];
866  t_pRapDipolePair.m_Dipole1.phi_z() = p_vec_phi_k_1[2];
867 
868  t_pRapDipolePair.m_Dipole2.phi_x() = p_vec_phi_k_1[3];
869  t_pRapDipolePair.m_Dipole2.phi_y() = p_vec_phi_k_1[4];
870  t_pRapDipolePair.m_Dipole2.phi_z() = p_vec_phi_k_1[5];
871 
872  t_pRapDipolePair.m_vCorrelation = p_valCor;
873 
874  p_RapDipoles.append(t_pRapDipolePair);
875 }
876 
877 //=============================================================================================================
878 
879 void RapMusic::setStcAttr(int p_iSampStcWin, float p_fStcOverlap)
880 {
881  m_iSamplesStcWindow = p_iSampStcWin;
882  m_fStcOverlap = p_fStcOverlap;
883 }
virtual const char * getName() const
Definition: rapmusic.cpp:186
Eigen::Matrix< double, 6, 6 > Matrix6T
Definition: rapmusic.h:109
virtual MNELIB::MNESourceEstimate calculateInverse(const FIFFLIB::FiffEvoked &p_fiffEvoked, bool pick_normal=false)
Definition: rapmusic.cpp:200
Eigen::RowVectorXf times
Definition: fiff_evoked.h:230
static int useFullRank(const MatrixXT &p_Mat, const MatrixXT &p_matSigma_src, MatrixXT &p_matFull_Rank, int type=NOT_TRANSPOSED)
Definition: rapmusic.h:374
RapMusic algorithm class declaration.
void setStcAttr(int p_iSampStcWin, float p_fStcOverlap)
Definition: rapmusic.cpp:879
Eigen::Matrix< double, 6, Eigen::Dynamic > Matrix6XT
Definition: rapmusic.h:107
int calcPhi_s(const MatrixXT &p_matMeasurement, MatrixXT *&p_pMatPhi_s) const
Definition: rapmusic.cpp:581
virtual const MNELIB::MNESourceSpace & getSourceSpace() const
Definition: rapmusic.cpp:193
static void getGainMatrixPair(const MatrixXT &p_matGainMarix, MatrixX6T &p_matGainMarix_Pair, int p_iIdx1, int p_iIdx2)
Definition: rapmusic.cpp:835
static void getPointPair(const int p_iPoints, const int p_iCurIdx, int &p_iIdx1, int &p_iIdx2)
Definition: rapmusic.cpp:823
Dipole< T > m_Dipole1
Definition: dipole.h:59
int m_iNumLeadFieldCombinations
Definition: rapmusic.h:307
#define IS_TRANSPOSED
Definition: pwlrapmusic.h:74
static int getRank(const MatrixXT &p_matSigma)
Definition: rapmusic.h:358
Source Space descritpion.
evoked data
Definition: fiff_evoked.h:77
static double subcorr(MatrixX6T &p_matProj_G, const MatrixXT &p_pMatU_B)
Definition: rapmusic.cpp:610
static void insertSource(int p_iDipoleIdx1, int p_iDipoleIdx2, const Vector6T &p_vec_phi_k_1, double p_valCor, QList< DipolePair< double > > &p_RapDipoles)
Definition: rapmusic.cpp:846
MNEMath class declaration.
static void calcA_k_1(const MatrixX6T &p_matG_k_1, const Vector6T &p_matPhi_k_1, const int p_iIdxk_1, MatrixXT &p_matA_k_1)
Definition: rapmusic.cpp:740
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixXT
Definition: rapmusic.h:103
FIFFLIB::FiffNamedMatrix::SDPtr sol
void calcPairCombinations(const int p_iNumPoints, const int p_iNumCombinations, Pair **p_ppPairIdxCombinations) const
Definition: rapmusic.cpp:793
static MatrixXT makeSquareMat(const MatrixXT &p_matF)
Definition: rapmusic.h:391
bool init(MNELIB::MNEForwardSolution &p_pFwd, bool p_bSparsed=false, int p_iN=2, double p_dThr=0.5)
Definition: rapmusic.cpp:103
void calcOrthProj(const MatrixXT &p_matA_k_1, MatrixXT &p_matOrthProj) const
Definition: rapmusic.cpp:755
Eigen::Matrix< double, Eigen::Dynamic, 6 > MatrixX6T
Definition: rapmusic.h:105
Eigen::Matrix< double, 6, 1 > Vector6T
Definition: rapmusic.h:113
Eigen::MatrixXd data
Definition: fiff_evoked.h:231
MNELIB::MNEForwardSolution m_ForwardSolution
Definition: rapmusic.h:298
Dipole< T > m_Dipole2
Definition: dipole.h:62
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorXT
Definition: rapmusic.h:111
Pair ** m_ppPairIdxCombinations
Definition: rapmusic.h:309