{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# 08. Pre-stack (AVO) inversion\nPre-stack inversion represents one step beyond post-stack inversion in that\nnot only the profile of acoustic impedance can be inferred from seismic data,\nrather a set of elastic parameters is estimated from pre-stack data\n(i.e., angle gathers) using the information contained in the so-called\nAVO (amplitude versus offset) response. Such elastic parameters represent\nvital information for more sophisticated geophysical subsurface\ncharacterization than it would be possible to achieve working with\npost-stack seismic data.\n\nIn this tutorial, the :py:class:`pylops.avo.prestack.PrestackLinearModelling`\noperator is used for modelling of both 1d and 2d synthetic pre-stack seismic\ndata using 1d profiles or 2d models of different subsurface elastic parameters\n(P-wave velocity, S-wave velocity, and density) as input.\n\n\\begin{align}d(t, \\theta) = w(t) * \\sum_{i=1}^N G_i(t, \\theta) \\frac{\\mathrm{d}\\ln m_i(t)}{\\mathrm{d}t}\\end{align}\n\nwhere $\\mathbf{m}(t)=[V_P(t), V_S(t), \\rho(t)]$ is a vector containing\nthree elastic parameters at time $t$, $G_i(t, \\theta)$ are the\ncoefficients of the AVO parametrization used to model pre-stack data\nand $w(t)$ is the time domain seismic wavelet.\nIn compact form:\n\n\\begin{align}\\mathbf{d}= \\mathbf{W} \\mathbf{G} \\mathbf{D} \\mathbf{m}\\end{align}\n\nwhere $\\mathbf{W}$ is a convolution operator, $\\mathbf{G}$ is\nthe AVO modelling operator, $\\mathbf{D}$ is a block-diagonal\nderivative operator, and $\\mathbf{m}$ is the input model.\nSubsequently the elastic parameters are estimated via the\n:py:class:`pylops.avo.prestack.PrestackInversion` module.\nOnce again, a two-steps inversion strategy can also be used to deal\nwith the case of noisy data.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.signal import filtfilt\n\nimport pylops\nfrom pylops.utils.wavelets import ricker\n\nplt.close(\"all\")\nnp.random.seed(0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's start with a 1d example. A synthetic profile of acoustic impedance\nis created and data is modelled using both the dense and linear operator\nversion of :py:class:`pylops.avo.prestack.PrestackLinearModelling` operator\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# sphinx_gallery_thumbnail_number = 5\n\n# model\nnt0 = 301\ndt0 = 0.004\n\nt0 = np.arange(nt0) * dt0\nvp = 1200 + np.arange(nt0) + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 80, nt0))\nvs = 600 + vp / 2 + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 20, nt0))\nrho = 1000 + vp + filtfilt(np.ones(5) / 5.0, 1, np.random.normal(0, 30, nt0))\nvp[131:] += 500\nvs[131:] += 200\nrho[131:] += 100\nvsvp = 0.5\nm = np.stack((np.log(vp), np.log(vs), np.log(rho)), axis=1)\n\n# background model\nnsmooth = 50\nmback = filtfilt(np.ones(nsmooth) / float(nsmooth), 1, m, axis=0)\n\n# angles\nntheta = 21\nthetamin, thetamax = 0, 40\ntheta = np.linspace(thetamin, thetamax, ntheta)\n\n# wavelet\nntwav = 41\nwav = ricker(t0[: ntwav // 2 + 1], 15)[0]\n\n# lop\nPPop = pylops.avo.prestack.PrestackLinearModelling(\n    wav, theta, vsvp=vsvp, nt0=nt0, linearization=\"akirich\"\n)\n\n# dense\nPPop_dense = pylops.avo.prestack.PrestackLinearModelling(\n    wav, theta, vsvp=vsvp, nt0=nt0, linearization=\"akirich\", explicit=True\n)\n\n# data lop\ndPP = PPop * m.ravel()\ndPP = dPP.reshape(nt0, ntheta)\n\n# data dense\ndPP_dense = PPop_dense * m.T.ravel()\ndPP_dense = dPP_dense.reshape(ntheta, nt0).T\n\n# noisy data\ndPPn_dense = dPP_dense + np.random.normal(0, 1e-2, dPP_dense.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now invert our data and retrieve elastic profiles for both noise-free\nand noisy data using :py:class:`pylops.avo.prestack.PrestackInversion`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# dense\nminv_dense, dPP_dense_res = pylops.avo.prestack.PrestackInversion(\n    dPP_dense,\n    theta,\n    wav,\n    m0=mback,\n    linearization=\"akirich\",\n    explicit=True,\n    returnres=True,\n    **dict(cond=1e-10)\n)\n\n# lop\nminv, dPP_res = pylops.avo.prestack.PrestackInversion(\n    dPP,\n    theta,\n    wav,\n    m0=mback,\n    linearization=\"akirich\",\n    explicit=False,\n    returnres=True,\n    **dict(damp=1e-10, iter_lim=2000)\n)\n\n# dense noisy\nminv_dense_noise, dPPn_dense_res = pylops.avo.prestack.PrestackInversion(\n    dPPn_dense,\n    theta,\n    wav,\n    m0=mback,\n    linearization=\"akirich\",\n    explicit=True,\n    returnres=True,\n    **dict(cond=1e-1)\n)\n\n# lop noisy (with vertical smoothing)\nminv_noise, dPPn_res = pylops.avo.prestack.PrestackInversion(\n    dPPn_dense,\n    theta,\n    wav,\n    m0=mback,\n    linearization=\"akirich\",\n    explicit=False,\n    epsR=5e-1,\n    returnres=True,\n    **dict(damp=1e-1, iter_lim=100)\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "The data, inverted models and residuals are now displayed.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# data and model\nfig, (axd, axdn, axvp, axvs, axrho) = plt.subplots(1, 5, figsize=(8, 5), sharey=True)\naxd.imshow(\n    dPP_dense,\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-np.abs(dPP_dense).max(),\n    vmax=np.abs(dPP_dense).max(),\n)\naxd.set_title(\"Data\")\naxd.axis(\"tight\")\naxdn.imshow(\n    dPPn_dense,\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-np.abs(dPP_dense).max(),\n    vmax=np.abs(dPP_dense).max(),\n)\naxdn.set_title(\"Noisy Data\")\naxdn.axis(\"tight\")\naxvp.plot(vp, t0, \"k\", lw=4, label=\"True\")\naxvp.plot(np.exp(mback[:, 0]), t0, \"--r\", lw=4, label=\"Back\")\naxvp.plot(np.exp(minv_dense[:, 0]), t0, \"--b\", lw=2, label=\"Inv Dense\")\naxvp.plot(np.exp(minv[:, 0]), t0, \"--m\", lw=2, label=\"Inv Lop\")\naxvp.plot(np.exp(minv_dense_noise[:, 0]), t0, \"--c\", lw=2, label=\"Noisy Dense\")\naxvp.plot(np.exp(minv_noise[:, 0]), t0, \"--g\", lw=2, label=\"Noisy Lop\")\naxvp.set_title(r\"$V_P$\")\naxvs.plot(vs, t0, \"k\", lw=4, label=\"True\")\naxvs.plot(np.exp(mback[:, 1]), t0, \"--r\", lw=4, label=\"Back\")\naxvs.plot(np.exp(minv_dense[:, 1]), t0, \"--b\", lw=2, label=\"Inv Dense\")\naxvs.plot(np.exp(minv[:, 1]), t0, \"--m\", lw=2, label=\"Inv Lop\")\naxvs.plot(np.exp(minv_dense_noise[:, 1]), t0, \"--c\", lw=2, label=\"Noisy Dense\")\naxvs.plot(np.exp(minv_noise[:, 1]), t0, \"--g\", lw=2, label=\"Noisy Lop\")\naxvs.set_title(r\"$V_S$\")\naxrho.plot(rho, t0, \"k\", lw=4, label=\"True\")\naxrho.plot(np.exp(mback[:, 2]), t0, \"--r\", lw=4, label=\"Back\")\naxrho.plot(np.exp(minv_dense[:, 2]), t0, \"--b\", lw=2, label=\"Inv Dense\")\naxrho.plot(np.exp(minv[:, 2]), t0, \"--m\", lw=2, label=\"Inv Lop\")\naxrho.plot(np.exp(minv_dense_noise[:, 2]), t0, \"--c\", lw=2, label=\"Noisy Dense\")\naxrho.plot(np.exp(minv_noise[:, 2]), t0, \"--g\", lw=2, label=\"Noisy Lop\")\naxrho.set_title(r\"$\\rho$\")\naxrho.legend(loc=\"center left\", bbox_to_anchor=(1, 0.5))\naxd.axis(\"tight\")\nplt.tight_layout()\n\n# residuals\nfig, axs = plt.subplots(1, 4, figsize=(8, 5), sharey=True)\nfig.suptitle(\"Residuals\", fontsize=14, fontweight=\"bold\", y=0.95)\nim = axs[0].imshow(\n    dPP_dense_res,\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-0.1,\n    vmax=0.1,\n)\naxs[0].set_title(\"Dense\")\naxs[0].set_xlabel(r\"$\\theta$\")\naxs[0].set_ylabel(\"t[s]\")\naxs[0].axis(\"tight\")\naxs[1].imshow(\n    dPP_res,\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-0.1,\n    vmax=0.1,\n)\naxs[1].set_title(\"Lop\")\naxs[1].set_xlabel(r\"$\\theta$\")\naxs[1].axis(\"tight\")\naxs[2].imshow(\n    dPPn_dense_res,\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-0.1,\n    vmax=0.1,\n)\naxs[2].set_title(\"Noisy Dense\")\naxs[2].set_xlabel(r\"$\\theta$\")\naxs[2].axis(\"tight\")\naxs[3].imshow(\n    dPPn_res,\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-0.1,\n    vmax=0.1,\n)\naxs[3].set_title(\"Noisy Lop\")\naxs[3].set_xlabel(r\"$\\theta$\")\naxs[3].axis(\"tight\")\nplt.tight_layout()\nplt.subplots_adjust(top=0.85)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally before moving to the 2d example, we consider the case when both PP\nand PS data are available. A joint PP-PS inversion can be easily solved\nas follows.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "PSop = pylops.avo.prestack.PrestackLinearModelling(\n    2 * wav, theta, vsvp=vsvp, nt0=nt0, linearization=\"ps\"\n)\nPPPSop = pylops.VStack((PPop, PSop))\n\n# data\ndPPPS = PPPSop * m.ravel()\ndPPPS = dPPPS.reshape(2, nt0, ntheta)\n\ndPPPSn = dPPPS + np.random.normal(0, 1e-2, dPPPS.shape)\n\n# Invert\nminvPPSP, dPPPS_res = pylops.avo.prestack.PrestackInversion(\n    dPPPS,\n    theta,\n    [wav, 2 * wav],\n    m0=mback,\n    linearization=[\"fatti\", \"ps\"],\n    epsR=5e-1,\n    returnres=True,\n    **dict(damp=1e-1, iter_lim=100)\n)\n\n# Data and model\nfig, (axd, axdn, axvp, axvs, axrho) = plt.subplots(1, 5, figsize=(8, 5), sharey=True)\naxd.imshow(\n    dPPPSn[0],\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-np.abs(dPPPSn[0]).max(),\n    vmax=np.abs(dPPPSn[0]).max(),\n)\naxd.set_title(\"PP Data\")\naxd.axis(\"tight\")\naxdn.imshow(\n    dPPPSn[1],\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], t0[-1], t0[0]),\n    vmin=-np.abs(dPPPSn[1]).max(),\n    vmax=np.abs(dPPPSn[1]).max(),\n)\naxdn.set_title(\"PS Data\")\naxdn.axis(\"tight\")\naxvp.plot(vp, t0, \"k\", lw=4, label=\"True\")\naxvp.plot(np.exp(mback[:, 0]), t0, \"--r\", lw=4, label=\"Back\")\naxvp.plot(np.exp(minv_noise[:, 0]), t0, \"--g\", lw=2, label=\"PP\")\naxvp.plot(np.exp(minvPPSP[:, 0]), t0, \"--b\", lw=2, label=\"PP+PS\")\naxvp.set_title(r\"$V_P$\")\naxvs.plot(vs, t0, \"k\", lw=4, label=\"True\")\naxvs.plot(np.exp(mback[:, 1]), t0, \"--r\", lw=4, label=\"Back\")\naxvs.plot(np.exp(minv_noise[:, 1]), t0, \"--g\", lw=2, label=\"PP\")\naxvs.plot(np.exp(minvPPSP[:, 1]), t0, \"--b\", lw=2, label=\"PP+PS\")\naxvs.set_title(r\"$V_S$\")\naxrho.plot(rho, t0, \"k\", lw=4, label=\"True\")\naxrho.plot(np.exp(mback[:, 2]), t0, \"--r\", lw=4, label=\"Back\")\naxrho.plot(np.exp(minv_noise[:, 2]), t0, \"--g\", lw=2, label=\"PP\")\naxrho.plot(np.exp(minvPPSP[:, 2]), t0, \"--b\", lw=2, label=\"PP+PS\")\naxrho.set_title(r\"$\\rho$\")\naxrho.legend(loc=\"center left\", bbox_to_anchor=(1, 0.5))\naxd.axis(\"tight\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We move now to a 2d example. First of all the model is loaded and\ndata generated.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# model\ninputfile = \"../testdata/avo/poststack_model.npz\"\n\nmodel = np.load(inputfile)\nx, z = model[\"x\"][::6] / 1000.0, model[\"z\"][:300] / 1000.0\nnx, nz = len(x), len(z)\nm = 1000 * model[\"model\"][:300, ::6]\n\nmvp = m.copy()\nmvs = m / 2\nmrho = m / 3 + 300\nm = np.log(np.stack((mvp, mvs, mrho), axis=1))\n\n# smooth model\nnsmoothz, nsmoothx = 30, 25\nmback = filtfilt(np.ones(nsmoothz) / float(nsmoothz), 1, m, axis=0)\nmback = filtfilt(np.ones(nsmoothx) / float(nsmoothx), 1, mback, axis=2)\n\n# dense operator\nPPop_dense = pylops.avo.prestack.PrestackLinearModelling(\n    wav,\n    theta,\n    vsvp=vsvp,\n    nt0=nz,\n    spatdims=(nx,),\n    linearization=\"akirich\",\n    explicit=True,\n)\n\n# lop operator\nPPop = pylops.avo.prestack.PrestackLinearModelling(\n    wav, theta, vsvp=vsvp, nt0=nz, spatdims=(nx,), linearization=\"akirich\"\n)\n\n# data\ndPP = PPop_dense * m.swapaxes(0, 1).ravel()\ndPP = dPP.reshape(ntheta, nz, nx).swapaxes(0, 1)\ndPPn = dPP + np.random.normal(0, 5e-2, dPP.shape)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we perform the same 4 different inversions as in the post-stack\ntutorial (see `sphx_glr_tutorials_poststack.py` for more details).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# dense inversion with noise-free data\nminv_dense = pylops.avo.prestack.PrestackInversion(\n    dPP, theta, wav, m0=mback, explicit=True, simultaneous=False\n)\n\n# dense inversion with noisy data\nminv_dense_noisy = pylops.avo.prestack.PrestackInversion(\n    dPPn, theta, wav, m0=mback, explicit=True, epsI=4e-2, simultaneous=False\n)\n\n# spatially regularized lop inversion with noisy data\nminv_lop_reg = pylops.avo.prestack.PrestackInversion(\n    dPPn,\n    theta,\n    wav,\n    m0=minv_dense_noisy,\n    explicit=False,\n    epsR=1e1,\n    **dict(damp=np.sqrt(1e-4), iter_lim=20)\n)\n\n# blockiness promoting inversion with noisy data\nminv_blocky = pylops.avo.prestack.PrestackInversion(\n    dPPn,\n    theta,\n    wav,\n    m0=mback,\n    explicit=False,\n    epsR=0.4,\n    epsRL1=0.1,\n    **dict(mu=0.1, niter_outer=3, niter_inner=3, iter_lim=5, damp=1e-3)\n)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Let's now visualize the inverted elastic parameters for the different\nscenarios\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "def plotmodel(\n    axs,\n    m,\n    x,\n    z,\n    vmin,\n    vmax,\n    params=(\"VP\", \"VS\", \"Rho\"),\n    cmap=\"gist_rainbow\",\n    title=None,\n):\n    \"\"\"Quick visualization of model\"\"\"\n    for ip, param in enumerate(params):\n        axs[ip].imshow(\n            m[:, ip], extent=(x[0], x[-1], z[-1], z[0]), vmin=vmin, vmax=vmax, cmap=cmap\n        )\n        axs[ip].set_title(\"%s - %s\" % (param, title))\n        axs[ip].axis(\"tight\")\n    plt.setp(axs[1].get_yticklabels(), visible=False)\n    plt.setp(axs[2].get_yticklabels(), visible=False)\n\n\n# data\nfig = plt.figure(figsize=(8, 9))\nax1 = plt.subplot2grid((2, 3), (0, 0), colspan=3)\nax2 = plt.subplot2grid((2, 3), (1, 0))\nax3 = plt.subplot2grid((2, 3), (1, 1), sharey=ax2)\nax4 = plt.subplot2grid((2, 3), (1, 2), sharey=ax2)\nax1.imshow(\n    dPP[:, 0], cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-0.4, vmax=0.4\n)\nax1.vlines(\n    [x[nx // 5], x[nx // 2], x[4 * nx // 5]],\n    ymin=z[0],\n    ymax=z[-1],\n    colors=\"w\",\n    linestyles=\"--\",\n)\nax1.set_xlabel(\"x [km]\")\nax1.set_ylabel(\"z [km]\")\nax1.set_title(r\"Stack ($\\theta$=0)\")\nax1.axis(\"tight\")\nax2.imshow(\n    dPP[:, :, nx // 5],\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], z[-1], z[0]),\n    vmin=-0.4,\n    vmax=0.4,\n)\nax2.set_xlabel(r\"$\\theta$\")\nax2.set_ylabel(\"z [km]\")\nax2.set_title(r\"Gather (x=%.2f)\" % x[nx // 5])\nax2.axis(\"tight\")\nax3.imshow(\n    dPP[:, :, nx // 2],\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], z[-1], z[0]),\n    vmin=-0.4,\n    vmax=0.4,\n)\nax3.set_xlabel(r\"$\\theta$\")\nax3.set_title(r\"Gather (x=%.2f)\" % x[nx // 2])\nax3.axis(\"tight\")\nax4.imshow(\n    dPP[:, :, 4 * nx // 5],\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], z[-1], z[0]),\n    vmin=-0.4,\n    vmax=0.4,\n)\nax4.set_xlabel(r\"$\\theta$\")\nax4.set_title(r\"Gather (x=%.2f)\" % x[4 * nx // 5])\nax4.axis(\"tight\")\nplt.setp(ax3.get_yticklabels(), visible=False)\nplt.setp(ax4.get_yticklabels(), visible=False)\n\n# noisy data\nfig = plt.figure(figsize=(8, 9))\nax1 = plt.subplot2grid((2, 3), (0, 0), colspan=3)\nax2 = plt.subplot2grid((2, 3), (1, 0))\nax3 = plt.subplot2grid((2, 3), (1, 1), sharey=ax2)\nax4 = plt.subplot2grid((2, 3), (1, 2), sharey=ax2)\nax1.imshow(\n    dPPn[:, 0], cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]), vmin=-0.4, vmax=0.4\n)\nax1.vlines(\n    [x[nx // 5], x[nx // 2], x[4 * nx // 5]],\n    ymin=z[0],\n    ymax=z[-1],\n    colors=\"w\",\n    linestyles=\"--\",\n)\nax1.set_xlabel(\"x [km]\")\nax1.set_ylabel(\"z [km]\")\nax1.set_title(r\"Noisy Stack ($\\theta$=0)\")\nax1.axis(\"tight\")\nax2.imshow(\n    dPPn[:, :, nx // 5],\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], z[-1], z[0]),\n    vmin=-0.4,\n    vmax=0.4,\n)\nax2.set_xlabel(r\"$\\theta$\")\nax2.set_ylabel(\"z [km]\")\nax2.set_title(r\"Gather (x=%.2f)\" % x[nx // 5])\nax2.axis(\"tight\")\nax3.imshow(\n    dPPn[:, :, nx // 2],\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], z[-1], z[0]),\n    vmin=-0.4,\n    vmax=0.4,\n)\nax3.set_title(r\"Gather (x=%.2f)\" % x[nx // 2])\nax3.set_xlabel(r\"$\\theta$\")\nax3.axis(\"tight\")\nax4.imshow(\n    dPPn[:, :, 4 * nx // 5],\n    cmap=\"gray\",\n    extent=(theta[0], theta[-1], z[-1], z[0]),\n    vmin=-0.4,\n    vmax=0.4,\n)\nax4.set_xlabel(r\"$\\theta$\")\nax4.set_title(r\"Gather (x=%.2f)\" % x[4 * nx // 5])\nax4.axis(\"tight\")\nplt.setp(ax3.get_yticklabels(), visible=False)\nplt.setp(ax4.get_yticklabels(), visible=False)\n\n# inverted models\nfig, axs = plt.subplots(6, 3, figsize=(8, 19))\nfig.suptitle(\"Model\", fontsize=12, fontweight=\"bold\", y=0.95)\nplotmodel(axs[0], m, x, z, m.min(), m.max(), title=\"True\")\nplotmodel(axs[1], mback, x, z, m.min(), m.max(), title=\"Back\")\nplotmodel(axs[2], minv_dense, x, z, m.min(), m.max(), title=\"Dense\")\nplotmodel(axs[3], minv_dense_noisy, x, z, m.min(), m.max(), title=\"Dense noisy\")\nplotmodel(axs[4], minv_lop_reg, x, z, m.min(), m.max(), title=\"Lop regularized\")\nplotmodel(axs[5], minv_blocky, x, z, m.min(), m.max(), title=\"Lop blocky\")\nplt.tight_layout()\nplt.subplots_adjust(top=0.92)\n\nfig, axs = plt.subplots(1, 3, figsize=(8, 7))\nfor ip, param in enumerate([\"VP\", \"VS\", \"Rho\"]):\n    axs[ip].plot(m[:, ip, nx // 2], z, \"k\", lw=4, label=\"True\")\n    axs[ip].plot(mback[:, ip, nx // 2], z, \"--r\", lw=4, label=\"Back\")\n    axs[ip].plot(minv_dense[:, ip, nx // 2], z, \"--b\", lw=2, label=\"Inv Dense\")\n    axs[ip].plot(\n        minv_dense_noisy[:, ip, nx // 2], z, \"--m\", lw=2, label=\"Inv Dense noisy\"\n    )\n    axs[ip].plot(\n        minv_lop_reg[:, ip, nx // 2], z, \"--g\", lw=2, label=\"Inv Lop regularized\"\n    )\n    axs[ip].plot(minv_blocky[:, ip, nx // 2], z, \"--y\", lw=2, label=\"Inv Lop blocky\")\n    axs[ip].set_title(param)\n    axs[ip].invert_yaxis()\naxs[2].legend(loc=8, fontsize=\"small\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "While the background model ``m0`` has been provided in all the examples so\nfar, it is worth showing that the module\n:py:class:`pylops.avo.prestack.PrestackInversion` can also produce so-called\nrelative elastic parameters (i.e., variations from an average medium\nproperty) when the background model ``m0`` is not available.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "dminv = pylops.avo.prestack.PrestackInversion(\n    dPP, theta, wav, m0=None, explicit=True, simultaneous=False\n)\n\nfig, axs = plt.subplots(1, 3, figsize=(8, 3))\nplotmodel(axs, dminv, x, z, -dminv.max(), dminv.max(), cmap=\"seismic\", title=\"relative\")\n\nfig, axs = plt.subplots(1, 3, figsize=(8, 7))\nfor ip, param in enumerate([\"VP\", \"VS\", \"Rho\"]):\n    axs[ip].plot(dminv[:, ip, nx // 2], z, \"k\", lw=2)\n    axs[ip].set_title(param)\n    axs[ip].invert_yaxis()"
      ]
    }
  ],
  "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
}