{
  "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\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport matplotlib.pyplot as plt\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 = .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)dt=-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 dt = \\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, halfcurrent=True)\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')"
      ]
    },
    {
      "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.RegularizedInversion(Cop, [Rop], yn, epsRs=[1e0],\n                                   **dict(iter_lim=100, atol=1e-5))\n\n# Preconditioned derivative\nSop = pylops.Smoothing1D(41, nt)\nxp = pylops.PreconditionedInversion(Cop, Sop, yn,\n                                    **dict(iter_lim=10, atol=1e-3))\n\n# Visualize data and inversion\nfig, axs = plt.subplots(1, 2, figsize=(18, 5))\naxs[0].plot(t, y, 'k', LineWidth=3, label='data')\naxs[0].plot(t, yn, '--g', LineWidth=3, label='noisy data')\naxs[0].legend()\naxs[0].set_title('Causal integration')\naxs[1].plot(t, x, 'k', LineWidth=8, label='original')\naxs[1].plot(t[1:-1], xder[1:-1], 'r', LineWidth=3, label='numerical derivative')\naxs[1].plot(t, xreg, 'g', LineWidth=3, label='regularized')\naxs[1].plot(t, xp, 'm', LineWidth=3, label='preconditioned')\naxs[1].legend()\naxs[1].set_title('Inverse causal integration')"
      ]
    },
    {
      "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 = .3\not = 0\nt = np.arange(nt)*dt+ot\nx = np.outer(np.sin(t), np.ones(nx))\n\nCop = pylops.CausalIntegration(nt*nx, dims=(nt, nx),\n                               sampling=dt, dir=0, halfcurrent=True)\n\ny = Cop*x.flatten()\ny = y.reshape(nt, nx)\nyn = y + np.random.normal(0, 4e-1, y.shape)\n\n# Numerical derivative\nDop = pylops.FirstDerivative(nt*nx, dims=(nt, nx), dir=0, sampling=dt)\nxder = Dop*yn.flatten()\nxder = xder.reshape(nt, nx)\n\n# Regularized derivative\nRop = pylops.Laplacian(dims=(nt, nx))\nxreg = pylops.RegularizedInversion(Cop, [Rop], yn.flatten(), epsRs=[1e0],\n                                   **dict(iter_lim=100, atol=1e-5))\nxreg = xreg.reshape(nt, nx)\n\n# Preconditioned derivative\nSop = pylops.Smoothing2D((11, 21), dims=(nt, nx))\nxp = pylops.PreconditionedInversion(Cop, Sop, yn.flatten(),\n                                    **dict(iter_lim=10, atol=1e-2))\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')\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', LineWidth=3, label='data')\naxs[0].plot(t, yn[:, nx//2], '--g', LineWidth=3, label='noisy data')\naxs[0].legend()\naxs[0].set_title('Causal integration')\naxs[1].plot(t, x[:, nx//2], 'k', LineWidth=8, label='original')\naxs[1].plot(t, xder[:, nx//2], 'r', LineWidth=3, label='numerical derivative')\naxs[1].plot(t, xreg[:, nx//2], 'g', LineWidth=3, label='regularized')\naxs[1].plot(t, xp[:, nx//2], 'm', LineWidth=3, label='preconditioned')\naxs[1].legend()\naxs[1].set_title('Inverse causal integration')"
      ]
    }
  ],
  "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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}