--- a/src/laser_sampling.py Thu Feb 26 09:32:12 2026 -0500
+++ b/src/laser_sampling.py Wed Apr 22 22:32:00 2026 -0500
@@ -1,13 +1,12 @@
import numpy as np
from dolfinx import fem, geometry
+from dolfinx_access import cell_FunctionSpace_f64
from mpi4py import MPI
from petsc4py import PETSc
from pointsource_pde.dolfinx_extras import get_mass_matrix
from scipy.sparse.linalg import svds
from slepc4py import SLEPc
-from dolfinx_access import cell_FunctionSpace_f64
-
class LaserSampling:
"""
@@ -52,14 +51,14 @@
E.setFromOptions()
E.solve()
- sigma = abs(E.getEigenvalue(0))
+ # sigma = abs(E.getEigenvalue(0))
# self._opnorm = 1.0 / sigma
# print("Full sampling opnorm = %f" % self._opnorm)
# SVD operator norm
_u, s, _v = svds(self.R, k=1)
self._opnorm = s[0] # / sigma
- print(f"Laser sampling opnorm = {self._opnorm:.6f}")
+ # print(f"Laser sampling opnorm = {self._opnorm:.6f}")
def matrix(self, beams):
"""Combined: generate beams + segment + cell collision + dof weighting"""
@@ -161,21 +160,17 @@
def lipschitz_factor(self, *args):
return self.opnorm(*args)
- def diff_adj_lipschitz_factor(self, *args):
+ def diff_chain_lipschitz_factor(self, *args):
return self.opnorm(*args)
def diff_bound(self, *args):
return self.opnorm(*args)
- def codomain_bound(self, xbound=None):
- if xbound is None:
- raise Exception("Linear operators have unbounded range")
- else:
- return self.opnorm() * xbound
-
# Construct observation noise
- def noise(self, noise_level):
+ def noise(self, noise_level, rng=None):
+ if rng is None:
+ rng = np.random.default_rng()
M = self.M
# Generate noise with shape (M, 1) using scalar standard deviation
- noise = np.random.normal(0, noise_level, size=(M,))
+ noise = rng.normal(0, noise_level, size=(M,))
return noise