{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Derivatives\nThis example shows how to use the suite of derivative operators, namely\n:py:class:`pylops.FirstDerivative`, :py:class:`pylops.SecondDerivative`,\n:py:class:`pylops.Laplacian` and :py:class:`pylops.Gradient`,\n:py:class:`pylops.FirstDirectionalDerivative` and\n:py:class:`pylops.SecondDirectionalDerivative`.\n\nThe derivative operators are very useful when the model to be inverted for\nis expect to be smooth in one or more directions. As shown in the\n*Optimization* tutorial, these operators will be used as part of the\nregularization term to obtain a smooth solution.\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\")\nnp.random.seed(0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by looking at a simple first-order centered derivative and how\ncould implement it naively by creating a dense matrix. Note that we will not\napply the derivative where the stencil is partially outside of the range of\nthe input signal (i.e., at the edge of the signal)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nx = 10\n\nD = np.diag(0.5 * np.ones(nx - 1), k=1) - np.diag(0.5 * np.ones(nx - 1), -1)\nD[0] = D[-1] = 0\n\nfig, ax = plt.subplots(1, 1, figsize=(6, 4))\nim = plt.imshow(D, cmap=\"rainbow\", vmin=-0.5, vmax=0.5)\nax.set_title(\"First derivative\", size=14, fontweight=\"bold\")\nax.set_xticks(np.arange(nx - 1) + 0.5)\nax.set_yticks(np.arange(nx - 1) + 0.5)\nax.grid(linewidth=3, color=\"white\")\nax.xaxis.set_ticklabels([])\nax.yaxis.set_ticklabels([])\nfig.colorbar(im, ax=ax, ticks=[-0.5, 0.5], shrink=0.7)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We now create a signal filled with zero and a single one at its center and\napply the derivative matrix by means of a dot product\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "x = np.zeros(nx)\nx[int(nx / 2)] = 1\n\ny_dir = np.dot(D, x)\nxadj_dir = np.dot(D.T, y_dir)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's now do the same using the :py:class:`pylops.FirstDerivative` operator\nand compare its outputs after applying the forward and adjoint operators\nto those from the dense matrix.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "D1op = pylops.FirstDerivative(nx, dtype=\"float32\")\n\ny_lop = D1op * x\nxadj_lop = D1op.H * y_lop\n\nfig, axs = plt.subplots(3, 1, figsize=(13, 8), sharex=True)\naxs[0].stem(np.arange(nx), x, linefmt=\"k\", markerfmt=\"ko\")\naxs[0].set_title(\"Input\", size=20, fontweight=\"bold\")\naxs[1].stem(np.arange(nx), y_dir, linefmt=\"k\", markerfmt=\"ko\", label=\"direct\")\naxs[1].stem(np.arange(nx), y_lop, linefmt=\"--r\", markerfmt=\"ro\", label=\"lop\")\naxs[1].set_title(\"Forward\", size=20, fontweight=\"bold\")\naxs[1].legend()\naxs[2].stem(np.arange(nx), xadj_dir, linefmt=\"k\", markerfmt=\"ko\", label=\"direct\")\naxs[2].stem(np.arange(nx), xadj_lop, linefmt=\"--r\", markerfmt=\"ro\", label=\"lop\")\naxs[2].set_title(\"Adjoint\", size=20, fontweight=\"bold\")\naxs[2].legend()\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As expected we obtain the same result, with the only difference that\nin the second case we did not need to explicitly create a matrix,\nsaving memory and computational time.\n\nLet's move onto applying the same first derivative to a 2d array in\nthe first direction\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nx, ny = 11, 21\nA = np.zeros((nx, ny))\nA[nx // 2, ny // 2] = 1.0\n\nD1op = pylops.FirstDerivative((nx, ny), axis=0, dtype=\"float64\")\nB = D1op * A\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 3), sharey=True)\nfig.suptitle(\n    \"First Derivative in 1st direction\", fontsize=12, fontweight=\"bold\", y=0.95\n)\nim = axs[0].imshow(A, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[0].axis(\"tight\")\naxs[0].set_title(\"x\")\nplt.colorbar(im, ax=axs[0])\nim = axs[1].imshow(B, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[1].axis(\"tight\")\naxs[1].set_title(\"y\")\nplt.colorbar(im, ax=axs[1])\nplt.tight_layout()\nplt.subplots_adjust(top=0.8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now do the same for the second derivative\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "A = np.zeros((nx, ny))\nA[nx // 2, ny // 2] = 1.0\n\nD2op = pylops.SecondDerivative(dims=(nx, ny), axis=0, dtype=\"float64\")\nB = D2op * A\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 3), sharey=True)\nfig.suptitle(\n    \"Second Derivative in 1st direction\", fontsize=12, fontweight=\"bold\", y=0.95\n)\nim = axs[0].imshow(A, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[0].axis(\"tight\")\naxs[0].set_title(\"x\")\nplt.colorbar(im, ax=axs[0])\nim = axs[1].imshow(B, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[1].axis(\"tight\")\naxs[1].set_title(\"y\")\nplt.colorbar(im, ax=axs[1])\nplt.tight_layout()\nplt.subplots_adjust(top=0.8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can also apply the second derivative to the second direction of\nour data (``axis=1``)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "D2op = pylops.SecondDerivative(dims=(nx, ny), axis=1, dtype=\"float64\")\nB = D2op * A\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 3), sharey=True)\nfig.suptitle(\n    \"Second Derivative in 2nd direction\", fontsize=12, fontweight=\"bold\", y=0.95\n)\nim = axs[0].imshow(A, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[0].axis(\"tight\")\naxs[0].set_title(\"x\")\nplt.colorbar(im, ax=axs[0])\nim = axs[1].imshow(B, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[1].axis(\"tight\")\naxs[1].set_title(\"y\")\nplt.colorbar(im, ax=axs[1])\nplt.tight_layout()\nplt.subplots_adjust(top=0.8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We use the symmetrical Laplacian operator as well\nas a asymmetrical version of it (by adding more weight to the\nderivative along one direction)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# symmetrical\nL2symop = pylops.Laplacian(dims=(nx, ny), weights=(1, 1), dtype=\"float64\")\n\n# asymmetrical\nL2asymop = pylops.Laplacian(dims=(nx, ny), weights=(3, 1), dtype=\"float64\")\n\nBsym = L2symop * A\nBasym = L2asymop * A\n\nfig, axs = plt.subplots(1, 3, figsize=(10, 3), sharey=True)\nfig.suptitle(\"Laplacian\", fontsize=12, fontweight=\"bold\", y=0.95)\nim = axs[0].imshow(A, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[0].axis(\"tight\")\naxs[0].set_title(\"x\")\nplt.colorbar(im, ax=axs[0])\nim = axs[1].imshow(Bsym, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[1].axis(\"tight\")\naxs[1].set_title(\"y sym\")\nplt.colorbar(im, ax=axs[1])\nim = axs[2].imshow(Basym, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[2].axis(\"tight\")\naxs[2].set_title(\"y asym\")\nplt.colorbar(im, ax=axs[2])\nplt.tight_layout()\nplt.subplots_adjust(top=0.8)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We consider now the gradient operator. Given a 2-dimensional array,\nthis operator applies first-order derivatives on both dimensions and\nconcatenates them.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Gop = pylops.Gradient(dims=(nx, ny), dtype=\"float64\")\n\nB = Gop * A\nC = Gop.H * B\n\nfig, axs = plt.subplots(2, 2, figsize=(10, 6), sharex=True, sharey=True)\nfig.suptitle(\"Gradient\", fontsize=12, fontweight=\"bold\", y=0.95)\nim = axs[0, 0].imshow(A, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[0, 0].axis(\"tight\")\naxs[0, 0].set_title(\"x\")\nplt.colorbar(im, ax=axs[0, 0])\nim = axs[0, 1].imshow(B[0, ...], interpolation=\"nearest\", cmap=\"rainbow\")\naxs[0, 1].axis(\"tight\")\naxs[0, 1].set_title(\"y - 1st direction\")\nplt.colorbar(im, ax=axs[0, 1])\nim = axs[1, 1].imshow(B[1, ...], interpolation=\"nearest\", cmap=\"rainbow\")\naxs[1, 1].axis(\"tight\")\naxs[1, 1].set_title(\"y - 2nd direction\")\nplt.colorbar(im, ax=axs[1, 1])\nim = axs[1, 0].imshow(C, interpolation=\"nearest\", cmap=\"rainbow\")\naxs[1, 0].axis(\"tight\")\naxs[1, 0].set_title(\"xadj\")\nplt.colorbar(im, ax=axs[1, 0])\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we use the Gradient operator to compute directional derivatives.\nWe create a model which has some layering in the horizontal and vertical\ndirections and show how the direction derivatives differs from standard\nderivatives\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nx, nz = 60, 40\n\nhorlayers = np.cumsum(np.random.uniform(2, 10, 20).astype(int))\nhorlayers = horlayers[horlayers < nz // 2]\nnhorlayers = len(horlayers)\n\nvertlayers = np.cumsum(np.random.uniform(2, 20, 10).astype(int))\nvertlayers = vertlayers[vertlayers < nx]\nnvertlayers = len(vertlayers)\n\nA = 1500 * np.ones((nz, nx))\nfor top, base in zip(horlayers[:-1], horlayers[1:]):\n    A[top:base] = np.random.normal(2000, 200)\nfor top, base in zip(vertlayers[:-1], vertlayers[1:]):\n    A[horlayers[-1] :, top:base] = np.random.normal(2000, 200)\n\nv = np.zeros((2, nz, nx))\nv[0, : horlayers[-1]] = 1\nv[1, horlayers[-1] :] = 1\n\nDdop = pylops.FirstDirectionalDerivative((nz, nx), v=v, sampling=(nz, nx))\nD2dop = pylops.SecondDirectionalDerivative((nz, nx), v=v, sampling=(nz, nx))\n\ndirder = Ddop * A\ndir2der = D2dop * A\n\njump = 4\nfig, axs = plt.subplots(3, 1, figsize=(4, 9), sharex=True)\nim = axs[0].imshow(A, cmap=\"gist_rainbow\", extent=(0, nx // jump, nz // jump, 0))\nq = axs[0].quiver(\n    np.arange(nx // jump) + 0.5,\n    np.arange(nz // jump) + 0.5,\n    np.flipud(v[1, ::jump, ::jump]),\n    np.flipud(v[0, ::jump, ::jump]),\n    color=\"w\",\n    linewidths=20,\n)\naxs[0].set_title(\"x\")\naxs[0].axis(\"tight\")\naxs[1].imshow(dirder, cmap=\"gray\", extent=(0, nx // jump, nz // jump, 0))\naxs[1].set_title(\"y = D * x\")\naxs[1].axis(\"tight\")\naxs[2].imshow(dir2der, cmap=\"gray\", extent=(0, nx // jump, nz // jump, 0))\naxs[2].set_title(\"y = D2 * x\")\naxs[2].axis(\"tight\")\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
}