Inverse medium scattering problem¶
In the following example we use the iteratively regularized Gauss-Newton solver from regpy.solvers.nonlinear.irgnm.IrgnmCG for the example of medium scattering.
In medium scattering we try to determine a perturbation \(f\) or its corresponding refractive index \(1+f\) of a medium from measurements of far field patterns \(u_{\infty}\) of scattered time-harmonic acoustic waves \(u_{sc}:=u - u_{inc}\) in this medium given some incident field \(u_{inc}\). The total field \(u\) satisfies
The Iteratively Regularized Gauss-Newton Method (IRGNM) minimizes in each iteration
where \(F\) is a Fréchet-differentiable operator. The minimum is determined using an implementation of the Tikhonov regularization using a CG method that is implemented as regpy.solvers.linear.tikhonov.TikhonovCG. The regularization parameter \(\alpha_n\) in this structure is a decreasing geometric sequence.
Basic imports and definition of logging level¶
[1]:
import numpy as np
import logging
import matplotlib.pyplot as plt
import matplotlib.colorbar as cbar
logging.basicConfig(
level=logging.INFO,
format='%(asctime)s %(levelname)s %(name)-20s :: %(message)s'
)
Defining the Operator¶
First define the medium scattering operator for fixed measurement directions using the general medium scattering operator supplied in the mediumscattering.py.
[2]:
from mediumscattering import MediumScatteringBase
from regpy.vecsps import UniformGridFcts
class MediumScatteringFixed(MediumScatteringBase):
"""Acoustic medium scattering with fixed measurement directions.
Parameters
----------
farfield_directions : array-like
Array of measurement directions of the farfield, shape `(n, 2)` or `(n, 3)` depending on
the problem dimension. All directions must be normalized.
**kwargs
All other (keyword-only) arguments are passed to the base class, which
see.
"""
def __init__(self, *, farfield_directions, **kwargs):
super().__init__(**kwargs)
farfield_directions = np.asarray(farfield_directions)
assert farfield_directions.ndim == 2
assert farfield_directions.shape[1] == self.domain.ndim
assert np.allclose(np.linalg.norm(farfield_directions, axis=-1), 1)
self.farfield_directions = farfield_directions
"""The farfield directions."""
self.farfield_matrix = self.normalization_factor * np.exp(
-1j * self.wave_number * (farfield_directions @ self.domain.coords[:, self.support])
)
"""The farfield matrix."""
self.codomain = UniformGridFcts(
axisdata=(self.farfield_directions, self.inc_directions),
dtype=complex
)
def _compute_farfield(self, farfield, inc_idx, v):
farfield[:, inc_idx] = self.farfield_matrix @ v[self.support]
def _compute_farfield_adjoint(self, farfield, inc_idx, v):
v[self.support] = farfield[:, inc_idx] @ self.farfield_matrix.conj()
---------------------------------------------------------------------------
ModuleNotFoundError Traceback (most recent call last)
Cell In[2], line 1
----> 1 from mediumscattering import MediumScatteringBase
2 from regpy.vecsps import UniformGridFcts
4 class MediumScatteringFixed(MediumScatteringBase):
ModuleNotFoundError: No module named 'mediumscattering'
Explicitly setting up a scattering operator using the above defined class.
Moreover, using the defined support and domain of the scattering we can define an embedding operator as the adjoint to a coordinate projection operator.
Finally we define the full forward operator as the composition of the embedding with the scattering operator.
[3]:
import regpy.util as util
# building the forward operator
radius = 1
scattering = MediumScatteringFixed(
gridshape=(64, 64),
radius=radius,
wave_number=1,
inc_directions=util.linspace_circle(16),
farfield_directions=util.linspace_circle(16),
)
from regpy.operators import CoordinateProjection
projection = CoordinateProjection(
scattering.domain,
scattering.support
)
embedding = projection.adjoint
op = scattering * embedding
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[3], line 4
2 # building the forward operator
3 radius = 1
----> 4 scattering = MediumScatteringFixed(
5 gridshape=(64, 64),
6 radius=radius,
7 wave_number=1,
8 inc_directions=util.linspace_circle(16),
9 farfield_directions=util.linspace_circle(16),
10 )
12 from regpy.operators import CoordinateProjection
14 projection = CoordinateProjection(
15 scattering.domain,
16 scattering.support
17 )
NameError: name 'MediumScatteringFixed' is not defined
Creating data¶
First we define a true solution as a contrast. Using this we can compute the exact data which we then perturb by some added Gaussian noise.
[4]:
#creating data
contrast = scattering.domain.zeros()
r = np.linalg.norm(scattering.domain.coords, axis=0)
contrast[r < radius] = np.exp(-1/(radius - r[r < radius]**2))
exact_solution = projection(contrast)
exact_data = op(exact_solution)
# create and add noise
noise = 0.005 * op.codomain.randn()
data = exact_data + noise
#plotting
fig,axs = plt.subplots(2,3,figsize=(8,4))
fig.tight_layout()
ax = axs[0,0]
im = ax.imshow(contrast.real)
fig.colorbar(im,ax=ax)
ax.set_ylabel('real')
ax.set_title('exact solution')
ax = axs[0,1]
im = ax.imshow(exact_data.real)
fig.colorbar(im,ax=ax)
ax.set_title('exact data')
ax = axs[0,2]
im = ax.imshow(data.real)
fig.colorbar(im,ax=ax)
ax.set_title('noisy data')
ax = axs[1,0]
ax.set_ylabel('imaginary')
im = ax.imshow(contrast.imag)
fig.colorbar(im,ax=ax)
ax.set_title('exact solution')
ax = axs[1,1]
im = ax.imshow(exact_data.imag)
fig.colorbar(im,ax=ax)
ax.set_title('exact data')
ax = axs[1,2]
im = ax.imshow(data.imag)
fig.colorbar(im,ax=ax)
ax.set_title('noisy data')
plt.show();
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[4], line 2
1 #creating data
----> 2 contrast = scattering.domain.zeros()
3 r = np.linalg.norm(scattering.domain.coords, axis=0)
4 contrast[r < radius] = np.exp(-1/(radius - r[r < radius]**2))
NameError: name 'scattering' is not defined
Define setting¶
In the regularization we want to have a Sobolev type smoothing penalty. Thus as a penalty we use the \(H^2\) space on the domain where we can explicitly pass the support. For the data fidelity we simply use an \(L^2\) setting.
[5]:
from regpy.hilbert import L2, HmDomain
from regpy.solvers import RegularizationSetting
#create setting
myh_domain = HmDomain(scattering.domain,scattering.support,dtype=complex,index=2)
setting = RegularizationSetting(
op=op,
# Define Sobolev norm on support via embedding
penalty = myh_domain,
data_fid=L2
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[5], line 4
2 from regpy.solvers import RegularizationSetting
3 #create setting
----> 4 myh_domain = HmDomain(scattering.domain,scattering.support,dtype=complex,index=2)
5 setting = RegularizationSetting(
6 op=op,
7 # Define Sobolev norm on support via embedding
8 penalty = myh_domain,
9 data_fid=L2
10 )
NameError: name 'scattering' is not defined
Regularization¶
Now using the setting and data we can setup the IRGNM solver choosing some of the parameters and defining the initial guess as the constant zero function.
Additionally, we define a stopping rule that is composed of an maximum iteration count of 100 and a discrepancy principle with \(\tau=1.1\).
[6]:
from regpy.solvers.nonlinear.irgnm import IrgnmCG
import regpy.stoprules as rules
init = op.domain.zeros()
#set up solver
solver = IrgnmCG(
setting, data,
regpar=0.0001, regpar_step=0.8,
init=init,
cg_pars=dict(
tol=1e-8,
reltolx=1e-8,
reltoly=1e-8
)
)
#set up stopping creiteria
stoprule = (
rules.CountIterations(100) +
rules.Discrepancy(
setting.h_codomain.norm, data,
noiselevel=setting.h_codomain.norm(noise),
tau=2.1
)
)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[6], line 4
1 from regpy.solvers.nonlinear.irgnm import IrgnmCG
2 import regpy.stoprules as rules
----> 4 init = op.domain.zeros()
5 #set up solver
6 solver = IrgnmCG(
7 setting, data,
8 regpar=0.0001, regpar_step=0.8,
(...) 14 )
15 )
NameError: name 'op' is not defined
Solving¶
Now we can run the solver using the stopping rule:
[7]:
reco, reco_data = solver.run(stoprule)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[7], line 1
----> 1 reco, reco_data = solver.run(stoprule)
NameError: name 'solver' is not defined
Plotting the results¶
[8]:
fig, axes = plt.subplots(ncols=3, nrows=2, constrained_layout=True,figsize=(8,4))
bars = np.vectorize(lambda ax: cbar.make_axes(ax)[0], otypes=[object])(axes)
axes[0, 0].set_title('exact contrast')
axes[1, 0].set_title('exact data')
axes[0, 1].set_title('reco contrast')
axes[1, 1].set_title('reco data')
axes[0, 2].set_title('difference')
def show(i, j, x):
im = axes[i, j].imshow(x)
bars[i, j].clear()
fig.colorbar(im, cax=bars[i, j])
show(0, 0, np.abs(contrast))
show(1, 0, np.abs(exact_data))
solution = embedding(reco)
show(0, 1, np.abs(solution))
show(1, 1, np.abs(reco_data))
show(0, 2, np.abs(solution - contrast))
show(1, 2, np.abs(exact_data - reco_data))
plt.show();
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In[8], line 15
12 bars[i, j].clear()
13 fig.colorbar(im, cax=bars[i, j])
---> 15 show(0, 0, np.abs(contrast))
16 show(1, 0, np.abs(exact_data))
17 solution = embedding(reco)
NameError: name 'contrast' is not defined
