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 #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)