{
  "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_distributed.FirstDerivative`,\n:py:class:`pylops_distributed.SecondDerivative`, and\n:py:class:`pylops_distributed.Laplacian`.\n\nThe derivative operators can be applied along any dimension of a N-dimensional\ninput. When the input is chuncked along any other direction(s) than the one\nthe derivative is applied, the derivative is efficiently performed without\nneither communication between workers nor duplication of part of the input\narray. On the other hand when the input is chuncked along the direction where\nthe derivative is applied, the chunks are partially overlapping such that no\ncommunication is required between the workers when applying the derivative.\n\nIn some applications, the user cannot avoid this second scenario to happen\n(e.g, when the derivative should be applied to all the dimensions of the\ndataset). Nevertheless our implementation makes this possible in a fully\ntransparent way and with very little additional overhead.\n\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import numpy as np\nimport dask.array as da\nimport matplotlib.pyplot as plt\n\nimport pylops\nimport pylops_distributed\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. We\nchunck the vector in 3 chunks.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nx = 100\nnchunks = 3\n\nx = np.zeros(nx)\nx[int(nx/2)] = 1\nxd = da.from_array(x, chunks=nx // nchunks + 1)\nprint('x:', xd)\n\nDop = pylops.FirstDerivative(nx, dtype='float32')\ndDop = pylops_distributed.FirstDerivative(nx, compute=(True, True),\n                                          dtype='float32')\n\ny = Dop * x\nxadj = Dop.H * y\n\nyd = Dop * xd\nxadjd = Dop.H * yd\n\n\nfig, axs = plt.subplots(3, 1, figsize=(13, 8))\naxs[0].stem(np.arange(nx), x, linefmt='k', markerfmt='ko',\n            use_line_collection=True)\naxs[0].set_title('Input', size=20, fontweight='bold')\naxs[1].stem(np.arange(nx), y, linefmt='--r', markerfmt='ro',\n            label='Numpy', use_line_collection=True)\naxs[1].stem(np.arange(nx), yd, linefmt='--r', markerfmt='ro',\n            label='Dask', use_line_collection=True)\naxs[1].set_title('Forward', size=20, fontweight='bold')\naxs[1].legend()\naxs[2].stem(np.arange(nx), xadj, linefmt='k', markerfmt='ko',\n            label='Numpy', use_line_collection=True)\naxs[2].stem(np.arange(nx), xadjd, linefmt='--r', markerfmt='ro',\n            label='Dask', use_line_collection=True)\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. Now we consider two cases, first when the data is\nchunked along the first direction and second when its chunked along the\nsecond direction.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nx, ny = 11, 21\nnchunks = 2\n\nA = np.zeros((nx, ny))\nA[nx//2, ny//2] = 1.\nA1d = da.from_array(A, chunks= (nx // nchunks + 1, ny))\nA2d = da.from_array(A, chunks= (nx , ny // nchunks + 1))\nprint('A1d:', A1d)\nprint('A2d:', A2d)\n\nDop = pylops_distributed.FirstDerivative(nx * ny, dims=(nx, ny),\n                                         compute=(True, True),\n                                         dir=0, dtype='float64')\n\nB1d = np.reshape(Dop * A1d.flatten(), (nx, ny))\nB2d = np.reshape(Dop * A2d.flatten(), (nx, ny))\n\nfig, axs = plt.subplots(1, 3, figsize=(12, 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(B1d, interpolation='nearest', cmap='rainbow')\naxs[1].axis('tight')\naxs[1].set_title('y (1st dir chunks)')\nplt.colorbar(im, ax=axs[1])\nim = axs[2].imshow(B2d, interpolation='nearest', cmap='rainbow')\naxs[2].axis('tight')\naxs[2].set_title('y (2nd dir chunks)')\nplt.colorbar(im, ax=axs[2])\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.\nA1d = da.from_array(A, chunks= (nx // nchunks + 1, ny))\nA2d = da.from_array(A, chunks= (nx , ny // nchunks + 1))\nprint('A1d:', A1d)\nprint('A2d:', A2d)\n\nDop = pylops_distributed.SecondDerivative(nx * ny, dims=(nx, ny),\n                                          compute=(True, True),\n                                          dir=0, dtype='float64')\n\nB1d = np.reshape(Dop * A1d.flatten(), (nx, ny))\nB2d = np.reshape(Dop * A2d.flatten(), (nx, ny))\n\nfig, axs = plt.subplots(1, 3, figsize=(12, 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(B1d, interpolation='nearest', cmap='rainbow')\naxs[1].axis('tight')\naxs[1].set_title('y (1st dir chunks)')\nplt.colorbar(im, ax=axs[1])\nim = axs[2].imshow(B2d, interpolation='nearest', cmap='rainbow')\naxs[2].axis('tight')\naxs[2].set_title('y (2nd dir chunks)')\nplt.colorbar(im, ax=axs[2])\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_distributed.Laplacian(dims=(nx, ny), weights=(1, 1),\n                                       compute=(True, True), dtype='float64')\n\n# asymmetrical\nL2asymop = pylops_distributed.Laplacian(dims=(nx, ny), weights=(3, 1),\n                                        compute=(True, True), dtype='float64')\n\nBsym = np.reshape(L2symop * A1d.flatten(), (nx, ny))\nBasym = np.reshape(L2asymop * A2d.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.8"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}