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