{
  "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_gpu.FirstDerivative`, :py:class:`pylops_gpu.SecondDerivative`\nand :py:class:`pylops_gpu.Laplacian`.\n\nThe derivative operators are very useful when the model to be inverted for\nis expect to be smooth in one or more directions. These operators\ncan in fact be used as part of the regularization term to obtain a smooth\nsolution.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import torch\nimport numpy as np\nimport matplotlib.pyplot as plt\nimport pylops_gpu\n\nfrom pylops_gpu.utils.backend import device\n\ndev = device()\nprint('PyLops-gpu working on %s...' % dev)\nplt.close('all')"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by looking at a simple first-order centered derivative. We\ncompute it by means of the :py:class:`pylops_gpu.FirstDerivative` operator.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nx = 10\nx = torch.zeros(nx, dtype=torch.float32)\nx[int(nx/2)] = 1\n\nD1op = pylops_gpu.FirstDerivative(nx, dtype=torch.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, basefmt='k', linefmt='k',\n            markerfmt='ko', use_line_collection=True)\naxs[0].set_title('Input', size=20, fontweight='bold')\naxs[1].stem(np.arange(nx), y_lop, basefmt='k', linefmt='k',\n            markerfmt='ko', use_line_collection=True)\naxs[1].set_title('Forward', size=20, fontweight='bold')\naxs[2].stem(np.arange(nx), xadj_lop, basefmt='k', linefmt='k',\n            markerfmt='ko', use_line_collection=True)\naxs[2].set_title('Adjoint', size=20, fontweight='bold')\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let'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 = torch.zeros((nx, ny), dtype=torch.float32)\nA[nx//2, ny//2] = 1.\n\nD1op = pylops_gpu.FirstDerivative(nx * ny, dims=(nx, ny),\n                                  dir=0, dtype=torch.float32)\nB = torch.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 = torch.zeros((nx, ny), dtype=torch.float32)\nA[nx//2, ny//2] = 1.\n\nD2op = pylops_gpu.SecondDerivative(nx * ny, dims=(nx, ny),\n                                   dir=0, dtype=torch.float32)\nB = torch.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": [
        "And finally 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_gpu.Laplacian(dims=(nx, ny), weights=(1, 1),\n                               dtype=torch.float32)\n\n# asymmetrical\nL2asymop = pylops_gpu.Laplacian(dims=(nx, ny), weights=(3, 1),\n                                dtype=torch.float32)\n\nBsym = torch.reshape(L2symop * A.flatten(), (nx, ny))\nBasym = torch.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)"
      ]
    }
  ],
  "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
}