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