OpenMEEG
gain.h
Go to the documentation of this file.
1 /*
2 Project Name : OpenMEEG
3 
4 © INRIA and ENPC (contributors: Geoffray ADDE, Maureen CLERC, Alexandre
5 GRAMFORT, Renaud KERIVEN, Jan KYBIC, Perrine LANDREAU, Théodore PAPADOPOULO,
6 Emmanuel OLIVI
7 Maureen.Clerc.AT.inria.fr, keriven.AT.certis.enpc.fr,
8 kybic.AT.fel.cvut.cz, papadop.AT.inria.fr)
9 
10 The OpenMEEG software is a C++ package for solving the forward/inverse
11 problems of electroencephalography and magnetoencephalography.
12 
13 This software is governed by the CeCILL-B license under French law and
14 abiding by the rules of distribution of free software. You can use,
15 modify and/ or redistribute the software under the terms of the CeCILL-B
16 license as circulated by CEA, CNRS and INRIA at the following URL
17 "http://www.cecill.info".
18 
19 As a counterpart to the access to the source code and rights to copy,
20 modify and redistribute granted by the license, users are provided only
21 with a limited warranty and the software's authors, the holders of the
22 economic rights, and the successive licensors have only limited
23 liability.
24 
25 In this respect, the user's attention is drawn to the risks associated
26 with loading, using, modifying and/or developing or reproducing the
27 software by the user in light of its specific status of free software,
28 that may mean that it is complicated to manipulate, and that also
29 therefore means that it is reserved for developers and experienced
30 professionals having in-depth computer knowledge. Users are therefore
31 encouraged to load and test the software's suitability as regards their
32 requirements in conditions enabling the security of their systems and/or
33 data to be ensured and, more generally, to use and operate it in the
34 same conditions as regards security.
35 
36 The fact that you are presently reading this means that you have had
37 knowledge of the CeCILL-B license and that you accept its terms.
38 */
39 
40 #pragma once
41 
42 #define USE_GMRES 0
43 
44 #include "matrix.h"
45 #include "sparse_matrix.h"
46 #include "symmatrix.h"
47 #include "geometry.h"
48 #include "assemble.h"
49 #include "gmres.h"
50 
51 namespace OpenMEEG {
52 
53  class GainMEG : public Matrix {
54  public:
55  using Matrix::operator=;
56  GainMEG (const Matrix& GainMat): Matrix(GainMat) {}
57  GainMEG (const SymMatrix& HeadMatInv,const Matrix& SourceMat, const Matrix& Head2MEGMat, const Matrix& Source2MEGMat) {
58  *this = Source2MEGMat+(Head2MEGMat*HeadMatInv)*SourceMat;
59  }
60  ~GainMEG () {};
61  };
62 
63  class GainEEG : public Matrix {
64  public:
65  using Matrix::operator=;
66  GainEEG (const Matrix& GainMat): Matrix(GainMat) {}
67  GainEEG (const SymMatrix& HeadMatInv,const Matrix& SourceMat, const SparseMatrix& Head2EEGMat) {
68  *this = (Head2EEGMat*HeadMatInv)*SourceMat;
69  }
70  ~GainEEG () {};
71  };
72 
73  class GainEEGadjoint : public Matrix {
74  public:
75  using Matrix::operator=;
76  GainEEGadjoint (const Geometry& geo,const Matrix& dipoles,const SymMatrix& HeadMat, const SparseMatrix& Head2EEGMat) {
77  Matrix LeadField(Head2EEGMat.nlin(),dipoles.nlin());
78  int gauss_order = 3;
79  // Consider the GMRes solver for problem with dimension > 15,000 (3,000 vertices per interface) else use LAPACK solver
80  #if USE_GMRES
81  Matrix mtemp(Head2EEGMat.nlin(),HeadMat.nlin());
82  Jacobi<SymMatrix> M(HeadMat); // Jacobi preconditionner
83  #pragma omp parallel for
84  #ifndef OPENMP_3_0
85  for (int i=0;i<static_cast<int>(LeadField.nlin());++i) {
86  #else
87  for (unsigned i=0;i<LeadField.nlin();++i) {
88  #endif
89  Vector vtemp(HeadMat.nlin());
90  GMRes(HeadMat, M, vtemp, Head2EEGMat.getlin(i), 1e3, 1e-7, HeadMat.nlin()); // max number of iteration = 1000, and precision=1e-7 (1e-5 for faster resolution)
91  mtemp.setlin(i, vtemp);
92  #pragma omp critical
93  PROGRESSBAR(i, LeadField.nlin());
94  }
95  #else
96  Matrix mtemp(Head2EEGMat.transpose());
97  HeadMat.solveLin(mtemp); // solving the system AX=B with LAPACK
98  mtemp=mtemp.transpose();
99  #endif
100  for ( unsigned i = 0; i < LeadField.ncol(); ++i) {
101  LeadField.setcol(i,mtemp * DipSourceMat(geo, dipoles.submat(i, 1, 0, dipoles.ncol()), gauss_order, true, "").getcol(0)); // TODO ugly
102  PROGRESSBAR(i,LeadField.ncol());
103  }
104  *this = LeadField;
105  }
107  };
108 
109  class GainMEGadjoint : public Matrix {
110  public:
111  using Matrix::operator=;
112  GainMEGadjoint (const Geometry& geo, const Matrix& dipoles,
113  const SymMatrix& HeadMat,
114  const Matrix& Head2MEGMat,
115  const Matrix& Source2MEGMat) {
116  Matrix LeadField(Head2MEGMat.nlin(),dipoles.nlin());
117  int gauss_order = 3;
118  // Consider the GMRes solver for problem with dimension > 15,000 (3,000 vertices per interface) else use LAPACK solver
119  #if USE_GMRES
120  Matrix mtemp(Head2MEGMat.nlin(),HeadMat.nlin());
121  Jacobi<SymMatrix> M(HeadMat); // Jacobi preconditionner
122  #pragma omp parallel for
123  #ifndef OPENMP_3_0
124  for (int i=0;i<static_cast<int>(LeadField.nlin());++i) {
125  #else
126  for (unsigned i=0;i<LeadField.nlin();++i) {
127  #endif
128  Vector vtemp(HeadMat.nlin());
129  GMRes(HeadMat, M, vtemp, Head2MEGMat.getlin(i), 1e3, 1e-7, HeadMat.nlin()); // max number of iteration = 1000, and precision=1e-7 (1e-5 for faster resolution)
130  mtemp.setlin(i,vtemp);
131  #pragma omp critical
132  PROGRESSBAR(i,LeadField.nlin());
133  }
134  #else
135  Matrix mtemp(Head2MEGMat.transpose());
136  HeadMat.solveLin(mtemp); // solving the system AX=B with LAPACK
137  mtemp=mtemp.transpose();
138  #endif
139  for (unsigned i=0;i<LeadField.ncol();i++) {
140  LeadField.setcol(i, mtemp * DipSourceMat(geo, dipoles.submat(i, 1, 0, dipoles.ncol()), gauss_order, true, "").getcol(0)+Source2MEGMat.getcol(i)); // TODO ugly
141  PROGRESSBAR(i,LeadField.ncol());
142  }
143  *this = LeadField;
144  }
146  };
147 
149  public:
150  GainEEGMEGadjoint (const Geometry& geo,const Matrix& dipoles,const SymMatrix& HeadMat, const SparseMatrix& Head2EEGMat, const Matrix& Head2MEGMat, const Matrix& Source2MEGMat) {
151  unsigned gauss_order = 3;
152  this->EEGleadfield = Matrix(Head2EEGMat.nlin(), dipoles.nlin());
153  this->MEGleadfield = Matrix(Head2MEGMat.nlin(), dipoles.nlin());
155 
156  for ( unsigned i = 0; i < Head2EEGMat.nlin(); ++i) {
157  RHS.setlin(i, Head2EEGMat.getlin(i));
158  }
159  for ( unsigned i = 0; i < Head2MEGMat.nlin(); ++i) {
160  RHS.setlin(i + Head2EEGMat.nlin(), Head2MEGMat.getlin(i));
161  }
162  // Consider the GMRes solver for problem with dimension > 15,000 (3,000 vertices per interface) else use LAPACK solver
163  #if USE_GMRES
164  Matrix mtemp(RHS.nlin(), HeadMat.nlin());
165  Jacobi<SymMatrix> M(HeadMat); // Jacobi preconditionner
166  #pragma omp parallel for
167  #ifndef OPENMP_3_0
168  for (int i=0;i<static_cast<int>(RHS.nlin());++i) {
169  #else
170  for (unsigned i=0;i<RHS.nlin();++i) {
171  #endif
172  Vector vtemp(HeadMat.nlin());
173  GMRes(HeadMat, M, vtemp, RHS.getlin(i), 1e3, 1e-7, HeadMat.nlin()); // max number of iteration = 1000, and precision=1e-7 (1e-5 for faster resolution)
174  mtemp.setlin(i, vtemp);
175  #pragma omp critical
176  PROGRESSBAR(i, RHS.nlin());
177  }
178  #else
179  Matrix mtemp(RHS.transpose());
180  HeadMat.solveLin(mtemp); // solving the system AX=B with LAPACK
181  mtemp = mtemp.transpose();
182  #endif
183  for ( unsigned i = 0; i < dipoles.nlin(); ++i) {
184  Vector dsm = DipSourceMat(geo,dipoles.submat(i, 1, 0, dipoles.ncol()), gauss_order, true, "").getcol(0); // TODO ugly
185  EEGleadfield.setcol(i, mtemp.submat(0, Head2EEGMat.nlin(), 0, HeadMat.nlin()) * dsm);
186  MEGleadfield.setcol(i, mtemp.submat(Head2EEGMat.nlin(), Head2MEGMat.nlin(), 0, HeadMat.nlin()) * dsm + Source2MEGMat.getcol(i));
187  PROGRESSBAR(i, dipoles.nlin());
188  }
189  }
190 
191  void saveEEG( const std::string filename ) const { EEGleadfield.save(filename); }
192 
193  void saveMEG( const std::string filename ) const { MEGleadfield.save(filename); }
194 
196  private:
197  Matrix EEGleadfield;
199  };
200 
201  class GainInternalPot : public Matrix {
202  public:
203  using Matrix::operator=;
204  GainInternalPot (const SymMatrix& HeadMatInv, const Matrix& SourceMat, const Matrix& Head2IPMat, const Matrix& Source2IPMat) {
205  *this = Source2IPMat + (Head2IPMat * HeadMatInv) * SourceMat;
206  }
208  };
209 
210  class GainEITInternalPot : public Matrix {
211  public:
212  using Matrix::operator=;
213  GainEITInternalPot (const SymMatrix& HeadMatInv, const Matrix& SourceMat, const Matrix& Head2IPMat) {
214  *this = (Head2IPMat * HeadMatInv) * SourceMat;
215  }
217  };
218 }
GainEEG(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const SparseMatrix &Head2EEGMat)
Definition: gain.h:67
void saveMEG(const std::string filename) const
Definition: gain.h:193
Geometry contains the electrophysiological model Here are stored the vertices, meshes and domains...
Definition: geometry.h:61
GainMEG(const Matrix &GainMat)
Definition: gain.h:56
Vector getcol(size_t j) const
Definition: matrix.h:245
Matrix submat(size_t istart, size_t isize, size_t jstart, size_t jsize) const
Definition: matrix.h:219
SparseMatrix transpose() const
size_t nlin() const
Definition: linop.h:85
void saveEEG(const std::string filename) const
Definition: gain.h:191
Matrix transpose() const
GainEEGMEGadjoint(const Geometry &geo, const Matrix &dipoles, const SymMatrix &HeadMat, const SparseMatrix &Head2EEGMat, const Matrix &Head2MEGMat, const Matrix &Source2MEGMat)
Definition: gain.h:150
Vector getlin(size_t i) const
Definition: symmatrix.h:356
GainMEGadjoint(const Geometry &geo, const Matrix &dipoles, const SymMatrix &HeadMat, const Matrix &Head2MEGMat, const Matrix &Source2MEGMat)
Definition: gain.h:112
GainEEG(const Matrix &GainMat)
Definition: gain.h:66
Vector getlin(size_t i) const
Definition: matrix.h:256
GainEITInternalPot(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2IPMat)
Definition: gain.h:213
GainEEGadjoint(const Geometry &geo, const Matrix &dipoles, const SymMatrix &HeadMat, const SparseMatrix &Head2EEGMat)
Definition: gain.h:76
void save(const char *filename) const
Save Matrix to file (Format set using file name extension)
Vector getlin(size_t i) const
GainMEG(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2MEGMat, const Matrix &Source2MEGMat)
Definition: gain.h:57
GainInternalPot(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2IPMat, const Matrix &Source2IPMat)
Definition: gain.h:204
void setcol(size_t j, const Vector &v)
Definition: matrix.h:267
unsigned GMRes(const T &A, const P &M, Vector &x, const Vector &b, int max_iter, double tol, unsigned m)
Definition: gmres.h:115
Vector solveLin(const Vector &B) const
Definition: symmatrix.h:142
#define PROGRESSBAR(a, b)
Definition: om_utils.h:61
virtual size_t ncol() const
Definition: linop.h:88
Matrix class.
Definition: matrix.h:61
void setlin(size_t i, const Vector &v)
Definition: matrix.h:276