OpenMEEG
Loading...
Searching...
No Matches
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
22namespace 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 ThreadException e;
33 #pragma omp parallel for
34 #ifdef OPENMP_UNSIGNED
35 for (unsigned i=0; i<S.nlin(); ++i) {
36 #else
37 for (int i=0; i<static_cast<int>(S.nlin()); ++i) {
38 #endif
39 e.Run([&](){
40 Vector vtemp(H.nlin());
41 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)
42 res.setlin(i,vtemp);
43 #pragma omp critical
44 PROGRESSBAR(i,S.nlin());
45 });
46 }
47 e.Rethrow();
48 return res
49 }
50#else
51 template <typename SelectionMatrix>
52 Matrix linsolve(const SymMatrix& H,const SelectionMatrix& S) {
53 Matrix res(S.transpose());
54 H.solveLin(res); // solving the system AX=B with LAPACK
55 return res.transpose();
56 }
57#endif
58
59 class GainMEG: public Matrix {
60 public:
61 using Matrix::operator=;
62 GainMEG (const Matrix& GainMat): Matrix(GainMat) {}
63 GainMEG(const SymMatrix& HeadMatInv,const Matrix& SourceMat,const Matrix& Head2MEGMat,const Matrix& Source2MEGMat):
64 Matrix(Source2MEGMat+(Head2MEGMat*HeadMatInv)*SourceMat)
65 { }
66 };
67
68 class GainEEG: public Matrix {
69 public:
70 using Matrix::operator=;
71 GainEEG (const Matrix& GainMat): Matrix(GainMat) {}
72 GainEEG (const SymMatrix& HeadMatInv,const Matrix& SourceMat,const SparseMatrix& Head2EEGMat):
73 Matrix((Head2EEGMat*HeadMatInv)*SourceMat)
74 { }
75 };
76
77 class GainEEGadjoint: public Matrix {
78 public:
79
80 using Matrix::operator=;
81
82 GainEEGadjoint(const Geometry& geo,const Matrix& dipoles,const SymMatrix& HeadMat,const SparseMatrix& Head2EEGMat): Matrix(Head2EEGMat.nlin(),dipoles.nlin()) {
83 const Matrix& Hinv = linsolve(HeadMat,Head2EEGMat);
84 ProgressBar pb(ncol());
85 for (unsigned i=0; i<ncol(); ++i,++pb)
86 setcol(i,Hinv*DipSourceMat(geo,dipoles.submat(i,1,0,dipoles.ncol()),"").getcol(0)); // TODO ugly
87 }
88 };
89
90 class GainMEGadjoint: public Matrix {
91 public:
92
93 using Matrix::operator=;
94
95 GainMEGadjoint(const Geometry& geo,const Matrix& dipoles,const SymMatrix& HeadMat,const Matrix& Head2MEGMat,const Matrix& Source2MEGMat):
96 Matrix(Head2MEGMat.nlin(),dipoles.nlin())
97 {
98 const Matrix& Hinv = linsolve(HeadMat,Head2MEGMat);
99 ProgressBar pb(ncol());
100 for (unsigned i=0; i<ncol(); ++i,++pb)
101 setcol(i,Hinv*DipSourceMat(geo,dipoles.submat(i,1,0,dipoles.ncol()),"").getcol(0)+Source2MEGMat.getcol(i)); // TODO ugly
102 }
103 };
104
106 public:
107 GainEEGMEGadjoint(const Geometry& geo,const Matrix& dipoles,const SymMatrix& HeadMat,const SparseMatrix& Head2EEGMat,const Matrix& Head2MEGMat,const Matrix& Source2MEGMat):
108 EEGleadfield(Head2EEGMat.nlin(),dipoles.nlin()),MEGleadfield(Head2MEGMat.nlin(),dipoles.nlin())
109 {
111 for (unsigned i=0; i<Head2EEGMat.nlin(); ++i)
112 RHS.setlin(i,Head2EEGMat.getlin(i));
113 for (unsigned i=0; i<Head2MEGMat.nlin(); ++i)
115
116 const Matrix& Hinv = linsolve(HeadMat,RHS);
117
118 ProgressBar pb(dipoles.nlin());
119 for (unsigned i=0; i<dipoles.nlin(); ++i,++pb) {
120 const Vector& dsm = DipSourceMat(geo,dipoles.submat(i,1,0,dipoles.ncol()),"").getcol(0); // TODO ugly
121 EEGleadfield.setcol(i,Hinv.submat(0,Head2EEGMat.nlin(),0,HeadMat.nlin())*dsm);
122 MEGleadfield.setcol(i,Hinv.submat(Head2EEGMat.nlin(),Head2MEGMat.nlin(),0,HeadMat.nlin())*dsm+Source2MEGMat.getcol(i));
123 }
124 }
125
126 void saveEEG( const std::string filename ) const { EEGleadfield.save(filename); }
127 void saveMEG( const std::string filename ) const { MEGleadfield.save(filename); }
128
129 size_t nlin() const { return MEGleadfield.nlin() + EEGleadfield.nlin(); }
130
131 private:
132
133 Matrix EEGleadfield;
134 Matrix MEGleadfield;
135 };
136
137 class GainInternalPot : public Matrix {
138 public:
139 using Matrix::operator=;
140 GainInternalPot (const SymMatrix& HeadMatInv,const Matrix& SourceMat,const Matrix& Head2IPMat,const Matrix& Source2IPMat):
141 Matrix(Source2IPMat+(Head2IPMat*HeadMatInv)*SourceMat)
142 { }
143 };
144
145 class GainEITInternalPot : public Matrix {
146 public:
147 using Matrix::operator=;
148 GainEITInternalPot (const SymMatrix& HeadMatInv,const Matrix& SourceMat,const Matrix& Head2IPMat):
149 Matrix((Head2IPMat*HeadMatInv)*SourceMat)
150 { }
151 };
152}
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:107
void saveMEG(const std::string filename) const
Definition gain.h:127
void saveEEG(const std::string filename) const
Definition gain.h:126
size_t nlin() const
Definition gain.h:129
GainEEG(const Matrix &GainMat)
Definition gain.h:71
GainEEG(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const SparseMatrix &Head2EEGMat)
Definition gain.h:72
GainEEGadjoint(const Geometry &geo, const Matrix &dipoles, const SymMatrix &HeadMat, const SparseMatrix &Head2EEGMat)
Definition gain.h:82
GainEITInternalPot(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2IPMat)
Definition gain.h:148
GainInternalPot(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2IPMat, const Matrix &Source2IPMat)
Definition gain.h:140
GainMEG(const SymMatrix &HeadMatInv, const Matrix &SourceMat, const Matrix &Head2MEGMat, const Matrix &Source2MEGMat)
Definition gain.h:63
GainMEG(const Matrix &GainMat)
Definition gain.h:62
GainMEGadjoint(const Geometry &geo, const Matrix &dipoles, const SymMatrix &HeadMat, const Matrix &Head2MEGMat, const Matrix &Source2MEGMat)
Definition gain.h:95
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:233
Vector getlin(const Index i) const
Definition matrix.h:246
Matrix submat(const Index istart, const Index isize, const Index jstart, const Index jsize) const
Definition matrix.h:208
void setcol(const Index j, const Vector &v)
Definition matrix.h:260
void setlin(const Index i, const Vector &v)
Definition matrix.h:271
Matrix transpose() const
Vector getlin(const size_t i) const
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:52
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)