{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Causal Integration\n\nThis example shows how to use the :py:class:`pylops.CausalIntegration`\noperator to integrate an input signal (in forward mode) and to apply a smooth,\nregularized derivative (in inverse mode). This is a very interesting\nby-product of this operator which may result very useful when the data\nto which you want to apply a numerical derivative is noisy.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\n\nplt.close(\"all\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start with a 1D example. Define the input parameters: number of samples\nof input signal (``nt``), sampling step (``dt``) as well as the input\nsignal which will be equal to $x(t)=\\sin(t)$:\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt = 81\ndt = 0.3\nt = np.arange(nt) * dt\nx = np.sin(t)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now create our causal integration operator and apply it to the input\nsignal. We can also compute the analytical integral\n$y(t)=\\int \\sin(t)\\,\\mathrm{d}t=-\\cos(t)$ and compare the results. We can also\ninvert the integration operator and by remembering that this is equivalent\nto a first order derivative, we will compare our inverted model with the\nresult obtained by simply applying the :py:class:`pylops.FirstDerivative`\nforward operator to the same data.\n\nNote that, as explained in details in :py:class:`pylops.CausalIntegration`,\nintegration has no unique solution, as any constant $c$ can be added\nto the integrated signal $y$, for example if $x(t)=t^2$ the\n$y(t) = \\int t^2 \\,\\mathrm{d}t = \\frac{t^3}{3} + c$. We thus subtract first\nsample from the analytical integral to obtain the same result as the\nnumerical one.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Cop = pylops.CausalIntegration(nt, sampling=dt, kind=\"half\")\n\nyana = -np.cos(t) + np.cos(t[0])\ny = Cop * x\nxinv = Cop / y\n\n# Numerical derivative\nDop = pylops.FirstDerivative(nt, sampling=dt)\nxder = Dop * y\n\n# Visualize data and inversion\nfig, axs = plt.subplots(1, 2, figsize=(18, 5))\naxs[0].plot(t, yana, \"r\", lw=5, label=\"analytic integration\")\naxs[0].plot(t, y, \"--g\", lw=3, label=\"numerical integration\")\naxs[0].legend()\naxs[0].set_title(\"Causal integration\")\n\naxs[1].plot(t, x, \"k\", lw=8, label=\"original\")\naxs[1].plot(t[1:-1], xder[1:-1], \"r\", lw=5, label=\"numerical\")\naxs[1].plot(t, xinv, \"--g\", lw=3, label=\"inverted\")\naxs[1].legend()\naxs[1].set_title(\"Inverse causal integration = Derivative\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As expected we obtain the same result. Let's see what happens if we now\nadd some random noise to our data.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Add noise\nyn = y + np.random.normal(0, 4e-1, y.shape)\n\n# Numerical derivative\nDop = pylops.FirstDerivative(nt, sampling=dt)\nxder = Dop * yn\n\n# Regularized derivative\nRop = pylops.SecondDerivative(nt)\nxreg = pylops.optimization.leastsquares.regularized_inversion(\n    Cop, yn, [Rop], epsRs=[1e0], **dict(iter_lim=100, atol=1e-5)\n)[0]\n\n# Preconditioned derivative\nSop = pylops.Smoothing1D(41, nt)\nxp = pylops.optimization.leastsquares.preconditioned_inversion(\n    Cop, yn, Sop, **dict(iter_lim=10, atol=1e-3)\n)[0]\n\n# Visualize data and inversion\nfig, axs = plt.subplots(1, 2, figsize=(18, 5))\naxs[0].plot(t, y, \"k\", lw=3, label=\"data\")\naxs[0].plot(t, yn, \"--g\", lw=3, label=\"noisy data\")\naxs[0].legend()\naxs[0].set_title(\"Causal integration\")\naxs[1].plot(t, x, \"k\", lw=8, label=\"original\")\naxs[1].plot(t[1:-1], xder[1:-1], \"r\", lw=3, label=\"numerical derivative\")\naxs[1].plot(t, xreg, \"g\", lw=3, label=\"regularized\")\naxs[1].plot(t, xp, \"m\", lw=3, label=\"preconditioned\")\naxs[1].legend()\naxs[1].set_title(\"Inverse causal integration\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can see here the great advantage of framing our numerical derivative\nas an inverse problem, and more specifically as the inverse of the\ncausal integration operator.\n\nLet's conclude with a 2d example where again the integration/derivative will\nbe performed along the first axis\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt, nx = 41, 11\ndt = 0.3\not = 0\nt = np.arange(nt) * dt + ot\nx = np.outer(np.sin(t), np.ones(nx))\n\nCop = pylops.CausalIntegration(dims=(nt, nx), sampling=dt, axis=0, kind=\"half\")\n\ny = Cop * x\nyn = y + np.random.normal(0, 4e-1, y.shape)\n\n# Numerical derivative\nDop = pylops.FirstDerivative(dims=(nt, nx), axis=0, sampling=dt)\nxder = Dop * yn\n\n# Regularized derivative\nRop = pylops.Laplacian(dims=(nt, nx))\nxreg = pylops.optimization.leastsquares.regularized_inversion(\n    Cop, yn.ravel(), [Rop], epsRs=[1e0], **dict(iter_lim=100, atol=1e-5)\n)[0]\nxreg = xreg.reshape(nt, nx)\n\n# Preconditioned derivative\nSop = pylops.Smoothing2D((11, 21), dims=(nt, nx))\nxp = pylops.optimization.leastsquares.preconditioned_inversion(\n    Cop, yn.ravel(), Sop, **dict(iter_lim=10, atol=1e-2)\n)[0]\nxp = xp.reshape(nt, nx)\n\n# Visualize data and inversion\nvmax = 2 * np.max(np.abs(x))\nfig, axs = plt.subplots(2, 3, figsize=(18, 12))\naxs[0][0].imshow(x, cmap=\"seismic\", vmin=-vmax, vmax=vmax)\naxs[0][0].set_title(\"Model\")\naxs[0][0].axis(\"tight\")\naxs[0][1].imshow(y, cmap=\"seismic\", vmin=-vmax, vmax=vmax)\naxs[0][1].set_title(\"Data\")\naxs[0][1].axis(\"tight\")\naxs[0][2].imshow(yn, cmap=\"seismic\", vmin=-vmax, vmax=vmax)\naxs[0][2].set_title(\"Noisy data\")\naxs[0][2].axis(\"tight\")\naxs[1][0].imshow(xder, cmap=\"seismic\", vmin=-vmax, vmax=vmax)\naxs[1][0].set_title(\"Numerical derivative\")\naxs[1][0].axis(\"tight\")\naxs[1][1].imshow(xreg, cmap=\"seismic\", vmin=-vmax, vmax=vmax)\naxs[1][1].set_title(\"Regularized\")\naxs[1][1].axis(\"tight\")\naxs[1][2].imshow(xp, cmap=\"seismic\", vmin=-vmax, vmax=vmax)\naxs[1][2].set_title(\"Preconditioned\")\naxs[1][2].axis(\"tight\")\nplt.tight_layout()\n\n# Visualize data and inversion at a chosen xlocation\nfig, axs = plt.subplots(1, 2, figsize=(18, 5))\naxs[0].plot(t, y[:, nx // 2], \"k\", lw=3, label=\"data\")\naxs[0].plot(t, yn[:, nx // 2], \"--g\", lw=3, label=\"noisy data\")\naxs[0].legend()\naxs[0].set_title(\"Causal integration\")\naxs[1].plot(t, x[:, nx // 2], \"k\", lw=8, label=\"original\")\naxs[1].plot(t, xder[:, nx // 2], \"r\", lw=3, label=\"numerical derivative\")\naxs[1].plot(t, xreg[:, nx // 2], \"g\", lw=3, label=\"regularized\")\naxs[1].plot(t, xp[:, nx // 2], \"m\", lw=3, label=\"preconditioned\")\naxs[1].legend()\naxs[1].set_title(\"Inverse causal integration\")\nplt.tight_layout()"
      ]
    }
  ],
  "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.9.15"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}