{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Inverse medium scattering problem\n", "\n", "In the following example we use the iteratively regularized Gauss-Newton solver from `regpy.solvers.nonlinear.irgnm.IrgnmCG` for the example of __medium scattering__.\n", "\n", "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 \n", "\n", "$$\n", " \\Delta u + k^2 (1+f) u = 0 \\qquad \\text{in } \\mathbb{R}^2\n", "$$\n", "\n", "The __Iteratively Regularized Gauss-Newton Method (IRGNM)__ minimizes in each iteration \n", "$$\n", " f_{n+1} = f_n + \\argmin_{h} \\Vert F(f_{n}) + F'[f_n] h - u_{\\infty}^{obs}\\Vert^{2} + \\alpha_{n} \\Vert f_{n} + h - f_0\\Vert^{2}\n", "$$\n", "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`.\n", "The regularization parameter $\\alpha_n$ in this structure is a decreasing geometric sequence." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Basic imports and definition of logging level" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import logging\n", "\n", "import matplotlib.pyplot as plt\n", "import matplotlib.colorbar as cbar\n", "\n", "\n", "logging.basicConfig(\n", " level=logging.INFO,\n", " format='%(asctime)s %(levelname)s %(name)-20s :: %(message)s'\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Defining the Operator\n", "\n", "First define the medium scattering operator for fixed measurement directions using the general medium scattering operator supplied in the `mediumscattering.py`. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from mediumscattering import MediumScatteringBase\n", "from regpy.vecsps import UniformGridFcts\n", "\n", "class MediumScatteringFixed(MediumScatteringBase):\n", " \"\"\"Acoustic medium scattering with fixed measurement directions.\n", "\n", " Parameters\n", " ----------\n", " farfield_directions : array-like\n", " Array of measurement directions of the farfield, shape `(n, 2)` or `(n, 3)` depending on\n", " the problem dimension. All directions must be normalized.\n", " **kwargs\n", " All other (keyword-only) arguments are passed to the base class, which\n", " see.\n", " \"\"\"\n", "\n", " def __init__(self, *, farfield_directions, **kwargs):\n", " super().__init__(**kwargs)\n", "\n", " farfield_directions = np.asarray(farfield_directions)\n", " assert farfield_directions.ndim == 2\n", " assert farfield_directions.shape[1] == self.domain.ndim\n", " assert np.allclose(np.linalg.norm(farfield_directions, axis=-1), 1)\n", " self.farfield_directions = farfield_directions\n", " \"\"\"The farfield directions.\"\"\"\n", " self.farfield_matrix = self.normalization_factor * np.exp(\n", " -1j * self.wave_number * (farfield_directions @ self.domain.coords[:, self.support])\n", " )\n", " \"\"\"The farfield matrix.\"\"\"\n", "\n", " self.codomain = UniformGridFcts(\n", " axisdata=(self.farfield_directions, self.inc_directions),\n", " dtype=complex\n", " )\n", "\n", " def _compute_farfield(self, farfield, inc_idx, v):\n", " farfield[:, inc_idx] = self.farfield_matrix @ v[self.support]\n", "\n", " def _compute_farfield_adjoint(self, farfield, inc_idx, v):\n", " v[self.support] = farfield[:, inc_idx] @ self.farfield_matrix.conj()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Explicitly setting up a scattering operator using the above defined class. \n", "\n", "Moreover, using the defined support and domain of the scattering we can define an embedding operator as the adjoint to a coordinate projection operator.\n", "\n", "Finally we define the full forward operator as the composition of the embedding with the scattering operator. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import regpy.util as util\n", "# building the forward operator\n", "radius = 1\n", "scattering = MediumScatteringFixed(\n", " gridshape=(64, 64),\n", " radius=radius,\n", " wave_number=1,\n", " inc_directions=util.linspace_circle(16),\n", " farfield_directions=util.linspace_circle(16),\n", ")\n", "\n", "from regpy.operators import CoordinateProjection\n", "\n", "projection = CoordinateProjection(\n", " scattering.domain,\n", " scattering.support\n", ")\n", "embedding = projection.adjoint\n", "\n", "op = scattering * embedding\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Creating data\n", "\n", "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. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "#creating data\n", "contrast = scattering.domain.zeros()\n", "r = np.linalg.norm(scattering.domain.coords, axis=0)\n", "contrast[r < radius] = np.exp(-1/(radius - r[r < radius]**2))\n", "\n", "exact_solution = projection(contrast)\n", "exact_data = op(exact_solution)\n", "# create and add noise\n", "noise = 0.005 * op.codomain.randn()\n", "data = exact_data + noise\n", "\n", "#plotting \n", "fig,axs = plt.subplots(2,3,figsize=(8,4))\n", "fig.tight_layout()\n", "ax = axs[0,0]\n", "im = ax.imshow(contrast.real)\n", "fig.colorbar(im,ax=ax)\n", "ax.set_ylabel('real')\n", "ax.set_title('exact solution')\n", "ax = axs[0,1]\n", "im = ax.imshow(exact_data.real)\n", "fig.colorbar(im,ax=ax)\n", "ax.set_title('exact data')\n", "ax = axs[0,2]\n", "im = ax.imshow(data.real)\n", "fig.colorbar(im,ax=ax)\n", "ax.set_title('noisy data')\n", "ax = axs[1,0]\n", "ax.set_ylabel('imaginary')\n", "im = ax.imshow(contrast.imag)\n", "fig.colorbar(im,ax=ax)\n", "ax.set_title('exact solution')\n", "ax = axs[1,1]\n", "im = ax.imshow(exact_data.imag)\n", "fig.colorbar(im,ax=ax)\n", "ax.set_title('exact data')\n", "ax = axs[1,2]\n", "im = ax.imshow(data.imag)\n", "fig.colorbar(im,ax=ax)\n", "ax.set_title('noisy data')\n", "plt.show();" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Define setting\n", "\n", "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. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from regpy.hilbert import L2, HmDomain\n", "from regpy.solvers import RegularizationSetting\n", "#create setting\n", "myh_domain = HmDomain(scattering.domain,scattering.support,dtype=complex,index=2)\n", "setting = RegularizationSetting(\n", " op=op,\n", " # Define Sobolev norm on support via embedding\n", " penalty = myh_domain, \n", " data_fid=L2\n", ")\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Regularization\n", "\n", "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. \n", "\n", "Additionally, we define a stopping rule that is composed of an maximum iteration count of 100 and a discrepancy principle with $\\tau=1.1$." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from regpy.solvers.nonlinear.irgnm import IrgnmCG\n", "import regpy.stoprules as rules\n", "\n", "init = op.domain.zeros()\n", "#set up solver\n", "solver = IrgnmCG(\n", " setting, data,\n", " regpar=0.0001, regpar_step=0.8,\n", " init=init,\n", " cg_pars=dict(\n", " tol=1e-8,\n", " reltolx=1e-8,\n", " reltoly=1e-8\n", " )\n", ")\n", "#set up stopping creiteria\n", "stoprule = (\n", " rules.CountIterations(100) +\n", " rules.Discrepancy(\n", " setting.h_codomain.norm, data,\n", " noiselevel=setting.h_codomain.norm(noise),\n", " tau=2.1\n", " )\n", ")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Solving\n", "\n", "Now we can run the solver using the stopping rule:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "reco, reco_data = solver.run(stoprule)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Plotting the results" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fig, axes = plt.subplots(ncols=3, nrows=2, constrained_layout=True,figsize=(8,4))\n", "bars = np.vectorize(lambda ax: cbar.make_axes(ax)[0], otypes=[object])(axes)\n", "\n", "axes[0, 0].set_title('exact contrast')\n", "axes[1, 0].set_title('exact data')\n", "axes[0, 1].set_title('reco contrast')\n", "axes[1, 1].set_title('reco data')\n", "axes[0, 2].set_title('difference')\n", "\n", "def show(i, j, x):\n", " im = axes[i, j].imshow(x)\n", " bars[i, j].clear()\n", " fig.colorbar(im, cax=bars[i, j])\n", "\n", "show(0, 0, np.abs(contrast))\n", "show(1, 0, np.abs(exact_data))\n", "solution = embedding(reco)\n", "show(0, 1, np.abs(solution))\n", "show(1, 1, np.abs(reco_data))\n", "show(0, 2, np.abs(solution - contrast))\n", "show(1, 2, np.abs(exact_data - reco_data))\n", "plt.show();\n" ] } ], "metadata": { "kernelspec": { "display_name": "3.12.3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.12.3" } }, "nbformat": 4, "nbformat_minor": 2 }