{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Volterra example with iterative regularised Gauss-Newton semismooth solver\n", "The discrete Volterra operator. The Volterra operator $V_n$ is defined as\n", "\n", "$$ (V_n f)(x) = \\int_0^x f(t)^n dt. $$\n", "\n", "Its discrete form, using a Riemann sum, is simply\n", "\n", "$$ (V_n x)_i = h \\sum_{j \\leq i} x_j^n, $$\n", "\n", "where $h$ is the grid spacing. $V_1$ is linear.\n", "\n", "Therefore, we create a subclass of `regpy.Operator` type in the next code block. In the following create new operators of that type can simply be done by initiating with that class. This class will have the following parameters. (Side remark you may also do the definition in a separate file for this example check out `volterra.py` in this example folder. )\n", "\n", "---\n", "\n", "### Parameters\n", "`domain` : `regpy.vecsps.UniformGridFcts`\n", " > The domain on which the operator is defined. Must be one-dimensional.\n", " \n", "`exponent` : `float`\n", " > The exponent \\(n\\). Default is 1.\n" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from regpy.operators import Operator\n", "from regpy.vecsps import UniformGridFcts\n", "\n", "\n", "class Volterra(Operator):\n", " def __init__(self, domain, exponent=1):\n", " assert isinstance(domain, UniformGridFcts)\n", " assert domain.ndim == 1\n", " self.exponent = exponent\n", " \"\"\"The exponent.\"\"\"\n", " super().__init__(domain, domain, linear=(exponent == 1))\n", "\n", " def _eval(self, x, differentiate=False):\n", " if differentiate:\n", " self._factor = self.exponent * x**(self.exponent - 1)\n", " return self.domain.volume_elem * np.cumsum(x**self.exponent)\n", "\n", " def _derivative(self, x):\n", " return self.domain.volume_elem * np.cumsum(self._factor * x)\n", "\n", " def _adjoint(self, y):\n", " x = self.domain.volume_elem * np.flipud(np.cumsum(np.flipud(y)))\n", " if self.linear:\n", " return x\n", " else:\n", " return self._factor * x" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Initialization \n", "\n", "Construction of a Volterra operator with exponent $n=2$ on the interval $[0,2\\pi]$ discretized by a `UniformGridFcts` with 200 points. Note, the operator will be a non-linear operator. Moreover, choose the exact solution \n", "$$f^\\dagger(x) = (1-\\cos(x))/2$$ \n", "for which the exact solution is given as \n", "$$g^\\dagger = (3x-4\\sin(x)+\\cos(x)\\sin(x))/8.$$ " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "grid = UniformGridFcts(np.linspace(0, 2 * np.pi, 200))\n", "op = Volterra(grid,exponent=2)\n", "\n", "exact_solution = (1-np.cos(grid.coords[0]))**2/4 \n", "plt.figure()\n", "plt.plot(grid.coords[0],exact_solution,color=\"blue\",label='exact solution')\n", "\n", "exact_data = (3*grid.coords[0] - 4*np.sin(grid.coords[0]) + np.cos(grid.coords[0])*np.sin(grid.coords[0]))/8\n", "noise = 0.3 * op.domain.randn()\n", "data = exact_data + noise\n", "plt.plot(grid.coords[0], exact_data, color=\"green\", label='exact data')\n", "plt.plot(grid.coords[0], data, color=\"lightgreen\", label='noisy data')\n", "plt.legend()\n", "\n", "print(\"Maximal error of operator evaluation on exact solution to exact data: {:.4E}\".format(np.max(np.abs(op(exact_solution)-exact_data))))\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Regularization\n", "\n", "First we define the regularization setting using the `RegularizationSetting` and choosing $L^2$ spaces defined on the `UniformGridFcts`.\n", "\n", "As a next step we define the we choose an initial guess, in this case we choose something that is almost zero. \n", "\n", "Use the Semi-Smooth Iterative Regularized Gauss Newton Method implemented in `regpy.solvers.nonlinear.irgnm_semismooth` to reconstruct the exact solution from the above constructed noisy data `data`. This method requires an a-priori lower and upper bound on the exact solution here chosen to be $0$ and $1.5$. In each step the initial regularization parameter `regpar` of the solver gets refined by `regpar_step`. We stop the iterative algorithm after at most $100$ steps and have as early stopping criterion the discrepancy rule implemented. This can be easily done by summing the two instances of the `regpy.stoprules`. \n", "\n", "After everything is defined run the solver with the specified stopping rule using the method `run()` of the solver. " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from regpy.hilbert import L2\n", "from regpy.solvers import RegularizationSetting\n", "from regpy.solvers.nonlinear.irgnm_semismooth import IrgnmSemiSmooth\n", "import regpy.stoprules as rules\n", "\n", "setting = RegularizationSetting(op, L2, L2)\n", "\n", "init = op.domain.ones()*0.05\n", "\n", "solver = IrgnmSemiSmooth(\n", " setting, # Regularization setting\n", " data, # noisy data constructed above\n", " psi_minus=0, # lower bound\n", " psi_plus=1.2, # upper bound\n", " regpar=1, # initial regularization parameter\n", " regpar_step=0.7, # regularization parameter step\n", " init=init, # initial guess\n", ")\n", "stoprule = (\n", " rules.CountIterations(100) +\n", " rules.Discrepancy(\n", " setting.h_codomain.norm, data,\n", " noiselevel=setting.h_codomain.norm(noise),\n", " tau=1.10\n", " )\n", ")\n", "\n", "recos = []\n", "recos_data = []\n", "for reco, reco_data in solver.while_(stoprule):\n", " recos.append(1*reco)\n", " recos_data.append(reco_data)\n", "\n", "import matplotlib.animation\n", "\n", "plt.rcParams[\"animation.html\"] = \"jshtml\"\n", "plt.rcParams['figure.dpi'] = 150 \n", "plt.ioff()\n", "fig, axs = plt.subplots(1,2)\n", "\n", "def animate(t):\n", " axs[0].cla()\n", " axs[0].plot(exact_solution,color=\"blue\")\n", " axs[0].plot(recos[t],color=\"lightblue\")\n", " axs[0].set_ylim(ymin=0,ymax=1.5)\n", "\n", " axs[1].cla()\n", " axs[1].plot(exact_data,color=\"green\", label='exact')\n", " axs[1].plot(recos_data[t],color=\"orange\", label='reco')\n", " axs[1].plot(data, color=\"lightgreen\", label='measured')\n", " axs[1].set_ylim(ymin=-0.2,ymax=2.5)\n", "\n", "matplotlib.animation.FuncAnimation(fig, animate, frames=len(recos))\n" ] } ], "metadata": { "kernelspec": { "display_name": "Python 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.10.12" } }, "nbformat": 4, "nbformat_minor": 2 }