OpenMEEG
gain.h
Go to the documentation of this file.
1 // Project Name: OpenMEEG (http://openmeeg.github.io)
2 // © INRIA and ENPC under the French open source license CeCILL-B.
3 // See full copyright notice in the file LICENSE.txt
4 // If you make a copy of this file, you must either:
5 // - provide also LICENSE.txt and modify this header to refer to it.
6 // - replace this header by the LICENSE.txt content.
7 
8 #pragma once
9 
10 #include "matrix.h"
11 #include "sparse_matrix.h"
12 #include "symmatrix.h"
13 #include "geometry.h"
14 #include "progressbar.h"
15 #include "assemble.h"
16 
17 #define USE_GMRES 0
18 #if USE_GMRES
19 #include "gmres.h"
20 #endif
21 
22 namespace OpenMEEG {
23 
24 #if USE_GMRES
25 
26  // Consider the GMRes solver for problems with dimension>15,000 (3,000 vertices per interface)
27 
28  template <typename MATRIX>
29  Matrix linsolve(const SymMatrix& H,const MATRIX& S) {
30  Matrix res(S.nlin(),H.nlin());
31  Jacobi<SymMatrix> M(H); // Jacobi preconditionner
32  #pragma omp parallel for
33  #ifdef OPENMP_UNSIGNED
34  for (unsigned i=0; i<S.nlin(); ++i) {
35  #else
36  for (int i=0; i<static_cast<int>(S.nlin()); ++i) {
37  #endif
38  Vector vtemp(H.nlin());
39  GMRes(H,M,vtemp,S.getlin(i),1000,1e-7,H.nlin()); // max number of iteration=1000, and precision=1e-7 (1e-5 for faster resolution)
40  res.setlin(i,vtemp);
41  #pragma omp critical
42  PROGRESSBAR(i,S.nlin());
43  }
44  return res
45  }
46 #else
47  template <typename SelectionMatrix>
48  Matrix linsolve(const SymMatrix& H,const SelectionMatrix& S) {
49  Matrix res(S.transpose());
50  H.solveLin(res); // solving the system AX=B with LAPACK
51  return res.transpose();
52  }
53 #endif
54 
55  class GainMEG: public Matrix {
56  public:
57  using Matrix::operator=;
58  GainMEG (const Matrix& GainMat): Matrix(GainMat) {}
59  GainMEG(const SymMatrix& HeadMatInv,const Matrix& SourceMat,const Matrix& Head2MEGMat,const Matrix& Source2MEGMat):
60  Matrix(Source2MEGMat+(Head2MEGMat*HeadMatInv)*SourceMat)
61  { }
62  };
63 
64  class GainEEG: public Matrix {
65  public:
66  using Matrix::operator=;
67  GainEEG (const Matrix& GainMat): Matrix(GainMat) {}
68  GainEEG (const SymMatrix& HeadMatInv,const Matrix& SourceMat,const SparseMatrix& Head2EEGMat):
69  Matrix((Head2EEGMat*HeadMatInv)*SourceMat)
70  { }
71  };
72 
73  class GainEEGadjoint: public Matrix {
74  public:
75 
76  using Matrix::operator=;
77 
78  GainEEGadjoint(const Geometry& geo,const Matrix& dipoles,const SymMatrix& HeadMat,const SparseMatrix& Head2EEGMat): Matrix(Head2EEGMat.nlin(),dipoles.nlin()) {
79  const Matrix& Hinv = linsolve(HeadMat,Head2EEGMat);
80  ProgressBar pb(ncol());
81  for (unsigned i=0; i<ncol(); ++i,++pb)
82  setcol(i,Hinv*DipSourceMat(geo,dipoles.submat(i,1,0,dipoles.ncol()),"").getcol(0)); // TODO ugly
83  }
84  };
85 
86  class GainMEGadjoint: public Matrix {
87  public:
88 
89  using Matrix::operator=;
90 
91  GainMEGadjoint(const Geometry& geo,const Matrix& dipoles,const SymMatrix& HeadMat,const Matrix& Head2MEGMat,const Matrix& Source2MEGMat):
92  Matrix(Head2MEGMat.nlin(),dipoles.nlin())
93  {
94  const Matrix& Hinv = linsolve(HeadMat,Head2MEGMat);
95  ProgressBar pb(ncol());
96  for (unsigned i=0; i<ncol(); ++i,++pb)
97  setcol(i,Hinv*DipSourceMat(geo,dipoles.submat(i,1,0,dipoles.ncol()),"").getcol(0)+Source2MEGMat.getcol(i)); // TODO ugly
98  }
99  };
100 
102  public:
103  GainEEGMEGadjoint(const Geometry& geo,const Matrix& dipoles,const SymMatrix& HeadMat,const SparseMatrix& Head2EEGMat,const Matrix& Head2MEGMat,const Matrix& Source2MEGMat):
104  EEGleadfield(Head2EEGMat.nlin(),dipoles.nlin()),MEGleadfield(Head2MEGMat.nlin(),dipoles.nlin())
105  {
107  for (unsigned i=0; i<Head2EEGMat.nlin(); ++i)
108  RHS.setlin(i,Head2EEGMat.getlin(i));
109  for (unsigned i=0; i<Head2MEGMat.nlin(); ++i)
111 
112  const Matrix& Hinv = linsolve(HeadMat,RHS);
113 
114  ProgressBar pb(dipoles.nlin());
115  for (unsigned i=0; i<dipoles.nlin(); ++i,++pb) {
116  const Vector& dsm = DipSourceMat(geo,dipoles.submat(i,1,0,dipoles.ncol()),"").getcol(0); // TODO ugly
117  EEGleadfield.setcol(i,Hinv.submat(0,Head2EEGMat.nlin(),0,HeadMat.nlin())*dsm);
118  MEGleadfield.setcol(i,Hinv.submat(Head2EEGMat.nlin(),Head2MEGMat.nlin(),0,HeadMat.nlin())*dsm+Source2MEGMat.getcol(i));
119  }
120  }
121 
122  void saveEEG( const std::string filename ) const { EEGleadfield.save(filename); }
123  void saveMEG( const std::string filename ) const { MEGleadfield.save(filename); }
124 
125  size_t nlin() const { return MEGleadfield.nlin() + EEGleadfield.nlin(); }
126 
127  private:
128 
129  Matrix EEGleadfield;
130  Matrix MEGleadfield;
131  };
132 
133  class GainInternalPot : public Matrix {
134  public:
135  using Matrix::operator=;
136  GainInternalPot (const SymMatrix& HeadMatInv,const Matrix& SourceMat,const Matrix& Head2IPMat,const Matrix& Source2IPMat):
137  Matrix(Source2IPMat+(Head2IPMat*HeadMatInv)*SourceMat)
138  { }
139  };
140 
141  class GainEITInternalPot : public Matrix {
142  public:
143  using Matrix::operator=;
144  GainEITInternalPot (const SymMatrix& HeadMatInv,const Matrix& SourceMat,const Matrix& Head2IPMat):
145  Matrix((Head2IPMat*HeadMatInv)*SourceMat)
146  { }
147  };
148 }
Various helper functions for assembling matrices.
GainEEGMEGadjoint(const Geometry &geo, const Matrix &dipoles, const SymMatrix &HeadMat, const SparseMatrix &Head2EEGMat, const Matrix &Head2MEGMat, const Matrix &Source2MEGMat)
Definition: gain.h:103
void saveMEG(const std::string filename) const
Definition: gain.h:123
void saveEEG(const std::string filename) const
Definition: gain.h:122
size_t nlin() const
Definition: gain.h:125
GainEEG(const Matrix &GainMat)
Definition: gain.h:67
GainEEG(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const SparseMatrix &Head2EEGMat)
Definition: gain.h:68
GainEEGadjoint(const Geometry &geo, const Matrix &dipoles, const SymMatrix &HeadMat, const SparseMatrix &Head2EEGMat)
Definition: gain.h:78
GainEITInternalPot(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2IPMat)
Definition: gain.h:144
GainInternalPot(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2IPMat, const Matrix &Source2IPMat)
Definition: gain.h:136
GainMEG(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2MEGMat, const Matrix &Source2MEGMat)
Definition: gain.h:59
GainMEG(const Matrix &GainMat)
Definition: gain.h:58
GainMEGadjoint(const Geometry &geo, const Matrix &dipoles, const SymMatrix &HeadMat, const Matrix &Head2MEGMat, const Matrix &Source2MEGMat)
Definition: gain.h:91
Geometry contains the electrophysiological model Vertices, meshes and domains are stored in this geom...
Definition: geometry.h:31
Dimension nlin() const
Definition: linop.h:48
virtual Dimension ncol() const
Definition: linop.h:51
Matrix class Matrix class.
Definition: matrix.h:28
void save(const char *filename) const
Save Matrix to file (Format set using file name extension)
Vector getcol(const Index j) const
Definition: matrix.h:223
Vector getlin(const Index i) const
Definition: matrix.h:236
Matrix submat(const Index istart, const Index isize, const Index jstart, const Index jsize) const
Definition: matrix.h:198
void setcol(const Index j, const Vector &v)
Definition: matrix.h:250
void setlin(const Index i, const Vector &v)
Definition: matrix.h:261
Matrix transpose() const
Vector getlin(const size_t i) const
Definition: sparse_matrix.h:67
Vector solveLin(const Vector &B) const
Definition: symmatrix.h:105
SymMatrix HeadMat(const Geometry &geo, const Integrator &integrator=Integrator(3, 0, 0.005))
Matrix linsolve(const SymMatrix &H, const SelectionMatrix &S)
Definition: gain.h:48
unsigned GMRes(const T &A, const P &M, Vector &x, const Vector &b, int max_iter, double tol, unsigned m)
Definition: gmres.h:83
Matrix DipSourceMat(const Geometry &geo, const Matrix &dipoles, const Integrator &integrator, const std::string &domain_name)
SparseMatrix Head2EEGMat(const Geometry &geo, const Sensors &electrodes)
Matrix Head2MEGMat(const Geometry &geo, const Sensors &sensors)