MNE-CPP  beta 1.0
pwlrapmusic.cpp
Go to the documentation of this file.
1 //=============================================================================================================
35 //*************************************************************************************************************
36 //=============================================================================================================
37 // INCLUDES
38 //=============================================================================================================
39 
40 #include "pwlrapmusic.h"
41 
42 #ifdef _OPENMP
43 #include <omp.h>
44 #endif
45 
46 
47 //*************************************************************************************************************
48 //=============================================================================================================
49 // USED NAMESPACES
50 //=============================================================================================================
51 
52 using namespace INVERSELIB;
53 
54 
55 //*************************************************************************************************************
56 //=============================================================================================================
57 // DEFINE MEMBER METHODS
58 //=============================================================================================================
59 
61 : RapMusic()
62 {
63 }
64 
65 
66 //*************************************************************************************************************
67 
68 PwlRapMusic::PwlRapMusic(MNEForwardSolution& p_pFwd, bool p_bSparsed, int p_iN, double p_dThr)
69 : RapMusic(p_pFwd, p_bSparsed, p_iN, p_dThr)
70 {
71  //Init
72  init(p_pFwd, p_bSparsed, p_iN, p_dThr);
73 }
74 
75 
76 //*************************************************************************************************************
77 
78 PwlRapMusic::~PwlRapMusic()
79 {
80 
81 }
82 
83 
84 //*************************************************************************************************************
85 
86 const char* PwlRapMusic::getName() const
87 {
88  return "Powell RAP MUSIC";
89 }
90 
91 
92 //*************************************************************************************************************
93 
94 MNESourceEstimate PwlRapMusic::calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal)
95 {
96  return RapMusic::calculateInverse(p_fiffEvoked, pick_normal);
97 }
98 
99 
100 //*************************************************************************************************************
101 
102 MNESourceEstimate PwlRapMusic::calculateInverse(const MatrixXd &data, float tmin, float tstep) const
103 {
104  return RapMusic::calculateInverse(data, tmin, tstep);
105 }
106 
107 
108 //*************************************************************************************************************
109 
110 MNESourceEstimate PwlRapMusic::calculateInverse(const MatrixXd& p_matMeasurement, QList< DipolePair<double> > &p_RapDipoles) const
111 {
112  MNESourceEstimate p_SourceEstimate;
113 
114  //if not initialized -> break
115  if(!m_bIsInit)
116  {
117  std::cout << "RAP MUSIC wasn't initialized!"; //ToDo: catch this earlier
118  return p_SourceEstimate; //false
119  }
120 
121  //Test if data are correct
122  if(p_matMeasurement.rows() != m_iNumChannels)
123  {
124  std::cout << "Lead Field channels do not fit to number of measurement channels!"; //ToDo: catch this earlier
125  return p_SourceEstimate;
126  }
127 
128  //Inits
129  //Stop the time for benchmark purpose
130  clock_t start, end;
131  start = clock();
132 
133 // //Map HPCMatrix to Eigen Matrix
134 // Eigen::Map<MatrixXT>
135 // t_MappedMatMeasurement( p_pMatMeasurement->data(),
136 // p_pMatMeasurement->rows(),
137 // p_pMatMeasurement->cols() );
138 
139  //Calculate the signal subspace (t_pMatPhi_s)
140  MatrixXT* t_pMatPhi_s = NULL;//(m_iNumChannels, m_iN < t_r ? m_iN : t_r);
141  int t_r = calcPhi_s(/*(MatrixXT)*/p_matMeasurement, t_pMatPhi_s);
142 
143  int t_iMaxSearch = m_iN < t_r ? m_iN : t_r; //The smallest of Rank and Iterations
144 
145  if (t_r < m_iN)
146  {
147  std::cout << "Warning: Rank " << t_r << " of the measurement data is smaller than the " << m_iN;
148  std::cout << " sources to find." << std::endl;
149  std::cout << " Searching now for " << t_iMaxSearch << " correlated sources.";
150  std::cout << std::endl << std::endl;
151  }
152 
153  //Create Orthogonal Projector
154  //OrthProj
155  MatrixXT t_matOrthProj(m_iNumChannels,m_iNumChannels);
156  t_matOrthProj.setIdentity();
157 
158  //A_k_1
159  MatrixXT t_matA_k_1(m_iNumChannels, t_iMaxSearch);
160  t_matA_k_1.setZero();
161 
162 // if (m_pMatGrid != NULL)
163 // {
164 // if(p_pRapDipoles != NULL)
165 // p_pRapDipoles->initRapDipoles(m_pMatGrid);
166 // else
167 // p_pRapDipoles = new RapDipoles<T>(m_pMatGrid);
168 // }
169 // else
170 // {
171 // if(p_pRapDipoles != NULL)
172 // delete p_pRapDipoles;
173 
174 // p_pRapDipoles = new RapDipoles<T>();
175 // }
176  p_RapDipoles.clear();
177 
178  std::cout << "##### Calculation of PWL RAP MUSIC started ######\n\n";
179 
180  MatrixXT t_matProj_Phi_s(t_matOrthProj.rows(), t_pMatPhi_s->cols());
181  //new Version: Calculate projection before
182  MatrixXT t_matProj_LeadField(m_ForwardSolution.sol->data.rows(), m_ForwardSolution.sol->data.cols());
183 
184  for(int r = 0; r < t_iMaxSearch ; ++r)
185  {
186  t_matProj_Phi_s = t_matOrthProj*(*t_pMatPhi_s);
187 
188  //new Version: Calculating Projection before
189  t_matProj_LeadField = t_matOrthProj * m_ForwardSolution.sol->data;//Subtract the found sources from the current found source
190 
191  //###First Option###
192  //Step 1: lt. Mosher 1998 -> Maybe tmp_Proj_Phi_S is already orthogonal -> so no SVD needed -> U_B = tmp_Proj_Phi_S;
193  Eigen::JacobiSVD< MatrixXT > t_svdProj_Phi_S(t_matProj_Phi_s, Eigen::ComputeThinU);
194  MatrixXT t_matU_B;
195  useFullRank(t_svdProj_Phi_S.matrixU(), t_svdProj_Phi_S.singularValues().asDiagonal(), t_matU_B);
196 
197  //Inits
199  t_vecRoh.setZero();
200 
201  //subcorr benchmark
202  //Stop the time
203  clock_t start_subcorr, end_subcorr;
204  start_subcorr = clock();
205 
206  double t_val_roh_k;
207 
208  //Powell
209  int t_iCurrentRow = 2;
210 
211  int t_iIdx1 = -1;
212  int t_iIdx2 = -1;
213 
214  int t_iMaxIdx_old = -1;
215 
216  int t_iMaxFound = 0;
217 
218  Eigen::VectorXi t_pVecIdxElements(m_iNumGridPoints);
219 
220  PowellIdxVec(t_iCurrentRow, m_iNumGridPoints, t_pVecIdxElements);
221 
222  int t_iNumVecElements = m_iNumGridPoints;
223 
224  while(t_iMaxFound == 0)
225  {
226 
227  //Multithreading correlation calculation
228  #ifdef _OPENMP
229  #pragma omp parallel num_threads(m_iMaxNumThreads)
230  #endif
231  {
232  #ifdef _OPENMP
233  #pragma omp for
234  #endif
235  for(int i = 0; i < t_iNumVecElements; i++)
236  {
237  int k = t_pVecIdxElements(i);
238  //new Version: calculate matrix multiplication before
239  //Create Lead Field combinations -> It would be better to use a pointer construction, to increase performance
240  MatrixX6T t_matProj_G(t_matProj_LeadField.rows(),6);
241 
242  int idx1 = m_ppPairIdxCombinations[k]->x1;
243  int idx2 = m_ppPairIdxCombinations[k]->x2;
244 
245  RapMusic::getGainMatrixPair(t_matProj_LeadField, t_matProj_G, idx1, idx2);
246 
247  t_vecRoh(k) = RapMusic::subcorr(t_matProj_G, t_matU_B);//t_vecRoh holds the correlations roh_k
248  }
249  }
250 
251  // if(r==0)
252  // {
253  // std::fstream filestr;
254  // std::stringstream filename;
255  // filename << "Roh_gold.txt";
256  //
257  // filestr.open ( filename.str().c_str(), std::fstream::out);
258  // for(int i = 0; i < m_iNumLeadFieldCombinations; ++i)
259  // {
260  // filestr << t_vecRoh(i) << "\n";
261  // }
262  // filestr.close();
263  //
264  // //exit(0);
265  // }
266 
267  //Find the maximum of correlation - can't put this in the for loop because it's running in different threads.
268 
269  VectorXT::Index t_iMaxIdx;
270 
271  t_val_roh_k = t_vecRoh.maxCoeff(&t_iMaxIdx);//p_vecCor = ^roh_k
272 
273  if((int)t_iMaxIdx == t_iMaxIdx_old)
274  {
275  t_iMaxFound = 1;
276  break;
277  }
278  else
279  {
280  t_iMaxIdx_old = t_iMaxIdx;
281  //get positions in sparsed leadfield from index combinations;
282  t_iIdx1 = m_ppPairIdxCombinations[t_iMaxIdx]->x1;
283  t_iIdx2 = m_ppPairIdxCombinations[t_iMaxIdx]->x2;
284  }
285 
286 
287  //set new index
288  if(t_iIdx1 == t_iCurrentRow)
289  t_iCurrentRow = t_iIdx2;
290  else
291  t_iCurrentRow = t_iIdx1;
292 
293  PowellIdxVec(t_iCurrentRow, m_iNumGridPoints, t_pVecIdxElements);
294  }
295 
296  //subcorr benchmark
297  end_subcorr = clock();
298 
299  float t_fSubcorrElapsedTime = ( (float)(end_subcorr-start_subcorr) / (float)CLOCKS_PER_SEC ) * 1000.0f;
300  std::cout << "Time Elapsed: " << t_fSubcorrElapsedTime << " ms" << std::endl;
301 
302 
303  // (Idx+1) because of MATLAB positions -> starting with 1 not with 0
304  std::cout << "Iteration: " << r+1 << " of " << t_iMaxSearch
305  << "; Correlation: " << t_val_roh_k<< "; Position (Idx+1): " << t_iIdx1+1 << " - " << t_iIdx2+1 <<"\n\n";
306 
307  //Calculations with the max correlated dipole pair G_k_1
308  MatrixX6T t_matG_k_1(m_ForwardSolution.sol->data.rows(),6);
309  RapMusic::getGainMatrixPair(m_ForwardSolution.sol->data, t_matG_k_1, t_iIdx1, t_iIdx2);
310 
311  MatrixX6T t_matProj_G_k_1(t_matOrthProj.rows(), t_matG_k_1.cols());
312  t_matProj_G_k_1 = t_matOrthProj * t_matG_k_1;//Subtract the found sources from the current found source
313 // MatrixX6T t_matProj_G_k_1(t_matProj_LeadField.rows(), 6);
314 // getLeadFieldPair(t_matProj_LeadField, t_matProj_G_k_1, t_iIdx1, t_iIdx2);
315 
316  //Calculate source direction
317  //source direction (p_pMatPhi) for current source r (phi_k_1)
318  Vector6T t_vec_phi_k_1(6, 1);
319  RapMusic::subcorr(t_matProj_G_k_1, t_matU_B, t_vec_phi_k_1);//Correlate the current source to calculate the direction
320 
321  //Set return values
322  RapMusic::insertSource(t_iIdx1, t_iIdx2, t_vec_phi_k_1, t_val_roh_k, p_RapDipoles);
323 
324  //Stop Searching when Correlation is smaller then the Threshold
325  if (t_val_roh_k < m_dThreshold)
326  {
327  std::cout << "Searching stopped, last correlation " << t_val_roh_k;
328  std::cout << " is smaller then the given threshold " << m_dThreshold << std::endl << std::endl;
329  break;
330  }
331 
332  //Calculate A_k_1 = [a_theta_1..a_theta_k_1] matrix for subtraction of found source
333  RapMusic::calcA_k_1(t_matG_k_1, t_vec_phi_k_1, r, t_matA_k_1);
334 
335  //Calculate new orthogonal Projector (Pi_k_1)
336  calcOrthProj(t_matA_k_1, t_matOrthProj);
337 
338  //garbage collecting
339  //ToDo
340  }
341 
342  std::cout << "##### Calculation of PWL RAP MUSIC completed ######"<< std::endl << std::endl << std::endl;
343 
344  end = clock();
345 
346  float t_fElapsedTime = ( (float)(end-start) / (float)CLOCKS_PER_SEC ) * 1000.0f;
347  std::cout << "Total Time Elapsed: " << t_fElapsedTime << " ms" << std::endl << std::endl;
348 
349  //garbage collecting
350  delete t_pMatPhi_s;
351 
352  return p_SourceEstimate;
353 }
354 
355 
356 //*************************************************************************************************************
357 
358 int PwlRapMusic::PowellOffset(int p_iRow, int p_iNumPoints)
359 {
360 
361  return p_iRow*p_iNumPoints - (( (p_iRow-1)*p_iRow) / 2); //triangular series 1 3 6 10 ... = (num_pairs*(num_pairs+1))/2
362 
363 }
364 
365 
366 //*************************************************************************************************************
367 
368 void PwlRapMusic::PowellIdxVec(int p_iRow, int p_iNumPoints, Eigen::VectorXi& p_pVecElements)
369 {
370 
371  // if(p_pVecElements != NULL)
372  // delete[] p_pVecElements;
373  //
374  // p_pVecElements = new int(p_iNumPoints);
375 
376 
377  //col combination index
378  for(int i = 0; i <= p_iRow; ++i)//=p_iNumPoints-1
379  p_pVecElements(i) = PwlRapMusic::PowellOffset(i+1,p_iNumPoints)-(p_iNumPoints-p_iRow);
380 
381 
382  //row combination index
383  int off = PwlRapMusic::PowellOffset(p_iRow,p_iNumPoints);
384  int length = p_iNumPoints - p_iRow;
385  int k=0;
386  for(int i = p_iRow; i < p_iRow+length; ++i)//=p_iNumPoints-1
387  {
388  p_pVecElements(i) = off+k;
389  k = k + 1;
390  }
391 }
Eigen::Matrix< double, 6, 1 > Vector6T
Definition: rapmusic.h:130
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
Eigen::Matrix< double, Eigen::Dynamic, 6 > MatrixX6T
Definition: rapmusic.h:122
int m_iNumLeadFieldCombinations
Definition: rapmusic.h:325
Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic > MatrixXT
Definition: rapmusic.h:120
evoked data
Definition: fiff_evoked.h:91
static void getGainMatrixPair(const MatrixXT &p_matGainMarix, MatrixX6T &p_matGainMarix_Pair, int p_iIdx1, int p_iIdx2)
Definition: rapmusic.cpp:859
bool init(MNEForwardSolution &p_pFwd, bool p_bSparsed=false, int p_iN=2, double p_dThr=0.5)
Definition: rapmusic.cpp:107
PwlRapMusic algorithm class declaration.
virtual MNESourceEstimate calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal=false)
Definition: rapmusic.cpp:208
virtual MNESourceEstimate calculateInverse(const FiffEvoked &p_fiffEvoked, bool pick_normal=false)
Definition: pwlrapmusic.cpp:94
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
Eigen::Matrix< double, Eigen::Dynamic, 1 > VectorXT
Definition: rapmusic.h:128
The RapMusic class provides the RAP MUSIC Algorithm CPU implementation. ToDo: Paper references...
Definition: rapmusic.h:109
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
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
virtual const char * getName() const
Definition: pwlrapmusic.cpp:86