{
  "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\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')\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))\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',\n            markerfmt='ko', label='direct')\naxs[2].stem(np.arange(nx), xadj_lop, linefmt='--r',\n            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.\n\nD1op = pylops.FirstDerivative(nx * ny, dims=(nx, ny), dir=0, dtype='float64')\nB = np.reshape(D1op * A.flatten(), (nx, ny))\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 3))\nfig.suptitle('First Derivative in 1st direction', fontsize=12,\n             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(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.\n\nD2op = pylops.SecondDerivative(nx * ny, dims=(nx, ny), dir=0, dtype='float64')\nB = np.reshape(D2op * A.flatten(), (nx, ny))\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 3))\nfig.suptitle('Second Derivative in 1st direction', fontsize=12,\n             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(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 (``dir=1``)\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "D2op = pylops.SecondDerivative(nx * ny, dims=(nx, ny),\n                               dir=1, dtype='float64')\nB = np.reshape(D2op*np.ndarray.flatten(A), (nx, ny))\n\nfig, axs = plt.subplots(1, 2, figsize=(10, 3))\nfig.suptitle('Second Derivative in 2nd direction', fontsize=12,\n             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(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 = np.reshape(L2symop * A.flatten(), (nx, ny))\nBasym = np.reshape(L2asymop * A.flatten(), (nx, ny))\n\nfig, axs = plt.subplots(1, 3, figsize=(10, 3))\nfig.suptitle('Laplacian', fontsize=12,\n             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 = np.reshape(Gop * A.flatten(), (2*nx, ny))\nC = np.reshape(Gop.H * B.flatten(), (nx, ny))\n\nfig, axs = plt.subplots(1, 3, figsize=(10, 3))\nfig.suptitle('Gradient', fontsize=12,\n             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(B, interpolation='nearest', cmap='rainbow')\naxs[1].axis('tight')\naxs[1].set_title('y')\nplt.colorbar(im, ax=axs[1])\nim = axs[2].imshow(C, interpolation='nearest', cmap='rainbow')\naxs[2].axis('tight')\naxs[2].set_title('xadj')\nplt.colorbar(im, ax=axs[2])\nplt.tight_layout()\nplt.subplots_adjust(top=0.8)"
      ]
    },
    {
      "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(np.int))\nhorlayers = horlayers[horlayers < nz//2]\nnhorlayers = len(horlayers)\n\nvertlayers = np.cumsum(np.random.uniform(2, 20, 10).astype(np.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.ravel()\ndirder = dirder.reshape(nz, nx)\ndir2der = D2dop * A.ravel()\ndir2der = dir2der.reshape(nz, nx)\n\njump = 4\nfig, axs = plt.subplots(3, 1, figsize=(4, 9))\nim = axs[0].imshow(A, cmap='gist_rainbow', extent=(0, nx//jump, nz//jump, 0))\nq = axs[0].quiver(np.arange(nx//jump)+.5, np.arange(nz//jump)+.5,\n                  np.flipud(v[1, ::jump, ::jump]),\n                  np.flipud(v[0, ::jump, ::jump]),\n                  color='w', linewidths=20)\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.6.12"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}