Mercurial > repos > pointsource_pde / file comparison
comparison: src/laser_sampling.py
src/laser_sampling.py
- changeset 3
- c3a4f4bb87f7
- parent 1
- a4137aedcb3a
equal
deleted
inserted
replaced
| 1:a4137aedcb3a | 3:c3a4f4bb87f7 |
|---|---|
| 1 import numpy as np | 1 import numpy as np |
| 2 from dolfinx import fem, geometry | 2 from dolfinx import fem, geometry |
| 3 from dolfinx_access import cell_FunctionSpace_f64 | |
| 3 from mpi4py import MPI | 4 from mpi4py import MPI |
| 4 from petsc4py import PETSc | 5 from petsc4py import PETSc |
| 5 from pointsource_pde.dolfinx_extras import get_mass_matrix | 6 from pointsource_pde.dolfinx_extras import get_mass_matrix |
| 6 from scipy.sparse.linalg import svds | 7 from scipy.sparse.linalg import svds |
| 7 from slepc4py import SLEPc | 8 from slepc4py import SLEPc |
| 8 | |
| 9 from dolfinx_access import cell_FunctionSpace_f64 | |
| 10 | 9 |
| 11 | 10 |
| 12 class LaserSampling: | 11 class LaserSampling: |
| 13 """ | 12 """ |
| 14 A linear operator class for laser sampling of a concentration (fem.Function) | 13 A linear operator class for laser sampling of a concentration (fem.Function) |
| 50 E.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_MAGNITUDE) | 49 E.setWhichEigenpairs(SLEPc.EPS.Which.SMALLEST_MAGNITUDE) |
| 51 # E.setWhichEigenpairs(SLEPc.EPS.Which.LARGEST_MAGNITUDE) | 50 # E.setWhichEigenpairs(SLEPc.EPS.Which.LARGEST_MAGNITUDE) |
| 52 E.setFromOptions() | 51 E.setFromOptions() |
| 53 E.solve() | 52 E.solve() |
| 54 | 53 |
| 55 sigma = abs(E.getEigenvalue(0)) | 54 # sigma = abs(E.getEigenvalue(0)) |
| 56 # self._opnorm = 1.0 / sigma | 55 # self._opnorm = 1.0 / sigma |
| 57 # print("Full sampling opnorm = %f" % self._opnorm) | 56 # print("Full sampling opnorm = %f" % self._opnorm) |
| 58 | 57 |
| 59 # SVD operator norm | 58 # SVD operator norm |
| 60 _u, s, _v = svds(self.R, k=1) | 59 _u, s, _v = svds(self.R, k=1) |
| 61 self._opnorm = s[0] # / sigma | 60 self._opnorm = s[0] # / sigma |
| 62 print(f"Laser sampling opnorm = {self._opnorm:.6f}") | 61 # print(f"Laser sampling opnorm = {self._opnorm:.6f}") |
| 63 | 62 |
| 64 def matrix(self, beams): | 63 def matrix(self, beams): |
| 65 """Combined: generate beams + segment + cell collision + dof weighting""" | 64 """Combined: generate beams + segment + cell collision + dof weighting""" |
| 66 dofmap = self.V.dofmap | 65 dofmap = self.V.dofmap |
| 67 tdim = self.domain.topology.dim | 66 tdim = self.domain.topology.dim |
| 159 return self._opnorm | 158 return self._opnorm |
| 160 | 159 |
| 161 def lipschitz_factor(self, *args): | 160 def lipschitz_factor(self, *args): |
| 162 return self.opnorm(*args) | 161 return self.opnorm(*args) |
| 163 | 162 |
| 164 def diff_adj_lipschitz_factor(self, *args): | 163 def diff_chain_lipschitz_factor(self, *args): |
| 165 return self.opnorm(*args) | 164 return self.opnorm(*args) |
| 166 | 165 |
| 167 def diff_bound(self, *args): | 166 def diff_bound(self, *args): |
| 168 return self.opnorm(*args) | 167 return self.opnorm(*args) |
| 169 | 168 |
| 170 def codomain_bound(self, xbound=None): | |
| 171 if xbound is None: | |
| 172 raise Exception("Linear operators have unbounded range") | |
| 173 else: | |
| 174 return self.opnorm() * xbound | |
| 175 | |
| 176 # Construct observation noise | 169 # Construct observation noise |
| 177 def noise(self, noise_level): | 170 def noise(self, noise_level, rng=None): |
| 171 if rng is None: | |
| 172 rng = np.random.default_rng() | |
| 178 M = self.M | 173 M = self.M |
| 179 # Generate noise with shape (M, 1) using scalar standard deviation | 174 # Generate noise with shape (M, 1) using scalar standard deviation |
| 180 noise = np.random.normal(0, noise_level, size=(M,)) | 175 noise = rng.normal(0, noise_level, size=(M,)) |
| 181 return noise | 176 return noise |
