{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Inverse potential problem \n",
    "\n",
    "For the inverse potential the goal is to find a homogeneous source given measurements at a circle enclosing the source. The Operator that maps the shape of a homogeneous heat source to the heat flux measured at some circle outside of the object is defined by the following partial differential equation. The heat distributions satisfies\n",
    "\n",
    "$$\n",
    "        \\begin{cases}\n",
    "            \\Delta u = 1_K & \\text{ in } \\Omega \\\\\n",
    "            u = 0          & \\text{ on } \\partial\\Omega\n",
    "        \\end{cases}\n",
    "$$\n",
    "\n",
    "where $\\partial\\Omega$ is the measurement circle and $K$ is the heat source. The forward operator\n",
    "maps the shape of the heat source to the Neumann data:\n",
    "\n",
    "$$\n",
    "        \\partial K \\mapsto \\frac{\\partial u}{\\partial\\nu}|_{\\partial\\Omega}.\n",
    "$$\n",
    "\n",
    "References:\n",
    "- F. Hettlich & W. Rundell \"Iterative methods for the reconstruction of an inverse potential problem\", Inverse Problems, 12 (1996) 251–266.\n",
    "- T. Hohage \"Logarithmic convergence rates of the iteratively regularized Gauss–Newton method for an inverse potential and an inverse scattering problem\", Inverse Problems, 13 (1997) 1279–1299.\n",
    "\n",
    "## Import auxiliary libraries and define logging\n",
    "\n",
    "We use `numpy` arrays as the main data structure for the domain and rely on `matplotlib.pyplot` as the tool for plots. \n",
    "\n",
    "To observe the information of the solver while computing we define the logging string and set the log level to info."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import logging\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "logging.basicConfig(\n",
    "    level=logging.INFO,\n",
    "    format='%(asctime)s %(levelname)s %(name)-40s :: %(message)s'\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Defining the Operator class\n",
    "\n",
    "We define a class for the potential operator that operators on a uniform grid to discretize the measurements on the circle and relies on the `StarTrigDiscr` class from `regpy` as a discretization of star shaped sources. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from regpy.operators import Operator\n",
    "from regpy.vecsps.curve import StarTrigDiscr\n",
    "from  regpy.vecsps import UniformGridFcts\n",
    "\n",
    "class Potential(Operator):\n",
    "    r\"\"\"Operator that maps the shape of a homogeneous heat source to the heat flux measured at some\n",
    "    circle outside of the object. The heat distributions satisfies\n",
    "\n",
    "\t.. math::\n",
    "        \\begin{cases}\n",
    "            \\Delta u = 1_K & \\text{ in } \\Omega \\\\\n",
    "            u = 0          & \\text{ on } \\partial\\Omega\n",
    "        \\end{cases}\n",
    "\n",
    "\n",
    "    where \\(\\partial\\Omega\\) is the measurement circle and \\(K\\) is the heat source. The operator\n",
    "    maps the shape of the heat source to the Neumann data:\n",
    "\n",
    "\t.. math::\n",
    "        \\partial K \\mapsto \\left.\\frac{\\partial u}{\\partial\\nu}\\right|_{\\partial\\Omega}.\n",
    "\n",
    "\n",
    "    Attributes\n",
    "    ----------\n",
    "    radius : float\n",
    "        The radius of the measurement circle.\n",
    "    N_means : int\n",
    "        The number of equispaced measurement points on the circle.\n",
    "    nforward : int, optional\n",
    "        The order of the Fourier expansion in the forward solver.\n",
    "\n",
    "    Raises\n",
    "    ------\n",
    "    ValueError\n",
    "        Will be raised on evaluating the operator if the object radius is negative or penetrates\n",
    "        the measurement circle.\n",
    "\n",
    "    References\n",
    "    ----------\n",
    "    - F. Hettlich & W. Rundell \"Iterative methods for the reconstruction of an inverse potential\n",
    "      problem\", Inverse Problems, 12 (1996) 251–266.\n",
    "    - T. Hohage \"Logarithmic convergence rates of the iteratively regularized\n",
    "      Gauss–Newton method for an inverse potential and an inverse scattering problem\", Inverse\n",
    "      Problems, 13 (1997) 1279–1299.\n",
    "    \"\"\"\n",
    "\n",
    "    def __init__(self, radius, nforward=128, N_means=128, N_ieq=128):\n",
    "        \n",
    "        self.radius = radius\n",
    "        \"\"\"The measurement radius.\"\"\"\n",
    "        self.nforward = nforward\n",
    "        \"\"\"The Fourier order of the forward solver.\"\"\"\n",
    "         \n",
    "        super().__init__(\n",
    "            domain=StarTrigDiscr(2*N_ieq),\n",
    "            codomain=UniformGridFcts(np.linspace(0, 2*np.pi, N_means, endpoint=False), dtype=complex)\n",
    "        )\n",
    "\n",
    "        k = 1 + np.arange(self.nforward)\n",
    "        k_t = np.outer(k, np.linspace(0, 2*np.pi, self.nforward, endpoint=False))\n",
    "        k_tfl = np.outer(k, self.codomain.coords[0])\n",
    "        self.cosin = np.cos(k_t)\n",
    "        self.sinus = np.sin(k_t)\n",
    "        self.cos_fl = np.cos(k_tfl)\n",
    "        self.sin_fl = np.sin(k_tfl)\n",
    "\n",
    "    def _eval(self, x, differentiate=False, adjoint_derivative=False):\n",
    "        nfwd = self.nforward\n",
    "        self._bd = self.domain.eval_curve(x, nvals=nfwd)\n",
    "\n",
    "        q = self._bd.radius[0]\n",
    "        if q.max() >= self.radius:\n",
    "            raise ValueError('Object penetrates measurement circle')\n",
    "        if q.min() <= 0:\n",
    "            raise ValueError('Radial function negative')\n",
    "\n",
    "        qq = q**2\n",
    "        flux = 1 / (2 * self.radius * nfwd) * np.sum(qq) * self.codomain.ones()\n",
    "        fac = 2 / (nfwd * self.radius)\n",
    "        for j in range(0, (nfwd - 1) // 2):\n",
    "            fac /= self.radius\n",
    "            qq *= q\n",
    "            flux += (\n",
    "                (fac / (j + 3)) * self.cos_fl[j, :] * np.sum(qq * self.cosin[j, :]) +\n",
    "                (fac / (j + 3)) * self.sin_fl[j, :] * np.sum(qq * self.sinus[j, :])\n",
    "            )\n",
    "\n",
    "        if nfwd % 2 == 0:\n",
    "            fac /= self.radius\n",
    "            qq *= q\n",
    "            flux += fac * self.cos_fl[:, nfwd // 2] * np.sum(qq * self.cosin[nfwd // 2, :])\n",
    "        return flux\n",
    "\n",
    "    def _derivative(self, h):\n",
    "        nfwd = self.nforward\n",
    "        q = self._bd.radius[0]\n",
    "        qqh = q * self._bd.derivative(h)\n",
    "\n",
    "        der = 1 / (self.radius * nfwd) * np.sum(qqh) * self.codomain.ones()\n",
    "        fac = 2 / (nfwd * self.radius)\n",
    "        for j in range((nfwd - 1) // 2):\n",
    "            fac /= self.radius\n",
    "            qqh *= q\n",
    "            der += fac * (\n",
    "                self.cos_fl[j, :] * np.sum(qqh * self.cosin[j, :]) +\n",
    "                self.sin_fl[j, :] * np.sum(qqh * self.sinus[j, :])\n",
    "            )\n",
    "\n",
    "        if nfwd % 2 == 0:\n",
    "            fac /= self.radius\n",
    "            qqh *= q\n",
    "            der += fac * self.cos_fl[nfwd // 2, :] * np.sum(qqh * self.cosin[nfwd // 2, :])\n",
    "        return der\n",
    "\n",
    "    def _adjoint(self, g):\n",
    "        nfwd = self.nforward\n",
    "        q = self._bd.radius[0]\n",
    "        qq = q.copy()\n",
    "\n",
    "        adj = 1 / (self.radius * nfwd) * np.sum(g) * qq\n",
    "        fac = 2 / (nfwd * self.radius)\n",
    "        for j in range((nfwd - 1) // 2):\n",
    "            fac /= self.radius\n",
    "            qq *= q\n",
    "            adj += fac * (\n",
    "                np.sum(g * self.cos_fl[j, :]) * (self.cosin[j, :] * qq) +\n",
    "                np.sum(g * self.sin_fl[j, :]) * (self.sinus[j, :] * qq)\n",
    "            )\n",
    "\n",
    "        if nfwd % 2 == 0:\n",
    "            fac /= self.radius\n",
    "            qq *= q\n",
    "            adj += fac * np.sum(g * self.cos_fl[nfwd // 2, :]) * (self.cosin[nfwd // 2, :] * qq)\n",
    "\n",
    "        return self._bd.adjoint(adj.real)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Define Operator and regularization setting\n",
    "\n",
    "Using the above class create the forward operator. \n",
    "\n",
    "The regularization setting is defined using a Sobolev penalty and $L^2$ data-fidelity. Note, that the construction here uses the abstract Sobolev which by default is the Sobolev space of with smoothness index 1. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from regpy.solvers import RegularizationSetting\n",
    "from regpy.hilbert import L2, Sobolev\n",
    "\n",
    "\n",
    "#Forward operator\n",
    "op = Potential(\n",
    "    radius=1.3\n",
    ")\n",
    "\n",
    "setting = RegularizationSetting(op=op, penalty=Sobolev, data_fid=L2)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The exact solution and data\n",
    "\n",
    "For the exact solution we define an star shaped domain with boundary given by $ \\sqrt{3\\cos{t}^2+1}/2$. We then define the exact data and add noise with 3% relative noise level. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Exact data and Poission data\n",
    "exact_solution = op.domain.sample(lambda t: np.sqrt(6*np.cos(1.5*t)**2+1)/3)\n",
    "exact_data = op(exact_solution)\n",
    "noise = op.codomain.randn()\n",
    "noise = 0.03*setting.h_codomain.norm(exact_data)/setting.h_codomain.norm(noise) * noise\n",
    "data = exact_data + noise\n",
    "\n",
    "fig = plt.figure(figsize=(7,7))\n",
    "plt.plot(*op.domain.eval_curve(exact_solution).curve[0])\n",
    "ax = plt.gca()\n",
    "ax.set_xlim([-1.1, 1.1])\n",
    "ax.set_ylim([-1.1, 1.1])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solving the problem\n",
    "\n",
    "With an initial guess that is the circle of radius 1 we use a Newton CG solver to get the solution. For a reference we plot the initial guess with the true objective. \n",
    "\n",
    "As stopping rules we rely on the discrepancy principle and to prevent an extensively long run we combine it with an iteration maximal count of 100.  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Initial guess\n",
    "init = op.domain.sample(lambda t: 0.2)\n",
    "\n",
    "fig = plt.figure(figsize=(7,7))\n",
    "plt.plot(*op.domain.eval_curve(exact_solution).curve[0])\n",
    "plt.plot(*op.domain.eval_curve(init).curve[0])\n",
    "ax = plt.gca()\n",
    "ax.set_xlim([-1.1, 1.1])\n",
    "ax.set_ylim([-1.1, 1.1])\n",
    "plt.show()\n",
    "\n",
    "from regpy.stoprules import CountIterations, Discrepancy\n",
    "\n",
    "from regpy.solvers.nonlinear.newton import NewtonCG\n",
    "solver = NewtonCG(\n",
    "    setting, \n",
    "    data, \n",
    "    init = init,\n",
    "    cgmaxit=50, \n",
    "    rho=0.6\n",
    ")\n",
    "\n",
    "stoprule = (\n",
    "    CountIterations(100) +\n",
    "    Discrepancy(\n",
    "        setting.h_codomain.norm, data,\n",
    "        noiselevel = setting.h_codomain.norm(noise),\n",
    "        tau=1.1\n",
    "    )\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot function\n",
    "for reco,reco_data in solver.until(stoprule=stoprule):\n",
    "    fig, axs = plt.subplots(1, 2,figsize=(15,7))\n",
    "    axs[0].set_title('Obstacle')\n",
    "    axs[0].set_xlim([-1.1, 1.1])\n",
    "    axs[0].set_ylim([-1.1, 1.1])\n",
    "    axs[1].set_title('Heat flux')\n",
    "\n",
    "    axs[0].plot(*op.domain.eval_curve(exact_solution).curve[0])\n",
    "    axs[0].plot(*op.domain.eval_curve(reco).curve[0])\n",
    "\n",
    "    axs[1].plot(exact_data.real, label='exact')\n",
    "    axs[1].plot(reco_data.real, label='reco')\n",
    "    axs[1].plot(data.real, label='measured')\n",
    "    axs[1].legend()\n",
    "    axs[1].set_ylim(ymin=0)\n",
    "\n",
    "    plt.show()"
   ]
  }
 ],
 "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
}