{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Multi-Dimensional Convolution\nThis example shows how to use the :py:class:`pylops.waveeqprocessing.MDC` operator\nto convolve a 3D kernel with an input seismic data. The resulting data is\na blurred version of the input data and the problem of removing such blurring\nis reffered to as *Multi-dimensional Deconvolution (MDD)* and its implementation\nis discussed in more details in the **MDD** tutorial.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\n\nimport pylops\nfrom pylops.utils.seismicevents import hyperbolic2d, makeaxis\nfrom pylops.utils.tapers import taper3d\nfrom pylops.utils.wavelets import ricker\n\nplt.close(\"all\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start by creating a set of hyperbolic events to be used as our MDC kernel\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Input parameters\npar = {\n    \"ox\": -300,\n    \"dx\": 10,\n    \"nx\": 61,\n    \"oy\": -500,\n    \"dy\": 10,\n    \"ny\": 101,\n    \"ot\": 0,\n    \"dt\": 0.004,\n    \"nt\": 400,\n    \"f0\": 20,\n    \"nfmax\": 200,\n}\n\nt0_m = 0.2\nvrms_m = 1100.0\namp_m = 1.0\n\nt0_G = (0.2, 0.5, 0.7)\nvrms_G = (1200.0, 1500.0, 2000.0)\namp_G = (1.0, 0.6, 0.5)\n\n# Taper\ntap = taper3d(par[\"nt\"], (par[\"ny\"], par[\"nx\"]), (5, 5), tapertype=\"hanning\")\n\n# Create axis\nt, t2, x, y = makeaxis(par)\n\n# Create wavelet\nwav = ricker(t[:41], f0=par[\"f0\"])[0]\n\n# Generate model\nm, mwav = hyperbolic2d(x, t, t0_m, vrms_m, amp_m, wav)\n\n# Generate operator\nG, Gwav = np.zeros((par[\"ny\"], par[\"nx\"], par[\"nt\"])), np.zeros(\n    (par[\"ny\"], par[\"nx\"], par[\"nt\"])\n)\nfor iy, y0 in enumerate(y):\n    G[iy], Gwav[iy] = hyperbolic2d(x - y0, t, t0_G, vrms_G, amp_G, wav)\nG, Gwav = G * tap, Gwav * tap\n\n# Add negative part to data and model\nm = np.concatenate((np.zeros((par[\"nx\"], par[\"nt\"] - 1)), m), axis=-1)\nmwav = np.concatenate((np.zeros((par[\"nx\"], par[\"nt\"] - 1)), mwav), axis=-1)\nGwav2 = np.concatenate((np.zeros((par[\"ny\"], par[\"nx\"], par[\"nt\"] - 1)), Gwav), axis=-1)\n\n# Define MDC linear operator\nGwav_fft = np.fft.rfft(Gwav2, 2 * par[\"nt\"] - 1, axis=-1)\nGwav_fft = Gwav_fft[..., : par[\"nfmax\"]]\n\n# Move frequency/time to first axis\nm, mwav = m.T, mwav.T\nGwav_fft = Gwav_fft.transpose(2, 0, 1)\n\n# Create operator\nMDCop = pylops.waveeqprocessing.MDC(\n    Gwav_fft,\n    nt=2 * par[\"nt\"] - 1,\n    nv=1,\n    dt=0.004,\n    dr=1.0,\n)\n\n# Create data\nd = MDCop * m.ravel()\nd = d.reshape(2 * par[\"nt\"] - 1, par[\"ny\"])\n\n# Apply adjoint operator to data\nmadj = MDCop.H * d.ravel()\nmadj = madj.reshape(2 * par[\"nt\"] - 1, par[\"nx\"])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally let's display the operator, input model, data and adjoint model\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(1, 2, figsize=(9, 6))\naxs[0].imshow(\n    Gwav2[int(par[\"ny\"] / 2)].T,\n    aspect=\"auto\",\n    interpolation=\"nearest\",\n    cmap=\"gray\",\n    vmin=-Gwav2.max(),\n    vmax=Gwav2.max(),\n    extent=(x.min(), x.max(), t2.max(), t2.min()),\n)\naxs[0].set_title(\"G - inline view\", fontsize=15)\naxs[0].set_xlabel(\"r\")\naxs[1].set_ylabel(\"t\")\naxs[1].imshow(\n    Gwav2[:, int(par[\"nx\"] / 2)].T,\n    aspect=\"auto\",\n    interpolation=\"nearest\",\n    cmap=\"gray\",\n    vmin=-Gwav2.max(),\n    vmax=Gwav2.max(),\n    extent=(y.min(), y.max(), t2.max(), t2.min()),\n)\naxs[1].set_title(\"G - inline view\", fontsize=15)\naxs[1].set_xlabel(\"s\")\naxs[1].set_ylabel(\"t\")\nfig.tight_layout()\n\nfig, axs = plt.subplots(1, 3, figsize=(9, 6))\naxs[0].imshow(\n    mwav,\n    aspect=\"auto\",\n    interpolation=\"nearest\",\n    cmap=\"gray\",\n    vmin=-mwav.max(),\n    vmax=mwav.max(),\n    extent=(x.min(), x.max(), t2.max(), t2.min()),\n)\naxs[0].set_title(r\"$m$\", fontsize=15)\naxs[0].set_xlabel(\"r\")\naxs[0].set_ylabel(\"t\")\naxs[1].imshow(\n    d,\n    aspect=\"auto\",\n    interpolation=\"nearest\",\n    cmap=\"gray\",\n    vmin=-d.max(),\n    vmax=d.max(),\n    extent=(x.min(), x.max(), t2.max(), t2.min()),\n)\naxs[1].set_title(r\"$d$\", fontsize=15)\naxs[1].set_xlabel(\"s\")\naxs[1].set_ylabel(\"t\")\naxs[2].imshow(\n    madj,\n    aspect=\"auto\",\n    interpolation=\"nearest\",\n    cmap=\"gray\",\n    vmin=-madj.max(),\n    vmax=madj.max(),\n    extent=(x.min(), x.max(), t2.max(), t2.min()),\n)\naxs[2].set_title(r\"$m_{adj}$\", fontsize=15)\naxs[2].set_xlabel(\"s\")\naxs[2].set_ylabel(\"t\")\nfig.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
}