{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# 15. Least-squares migration\nSeismic migration is the process by which seismic data are manipulated to create\nan image of the subsurface reflectivity.\n\nWhile traditionally solved as the adjont of the demigration operator,\nit is becoming more and more common to solve the underlying inverse problem\nin the quest for more accurate and detailed subsurface images.\n\nIndipendently of the choice of the modelling operator (i.e., ray-based or\nfull wavefield-based), the demigration/migration process can be expressed as\na linear operator of such a kind:\n\n\\begin{align}d(\\mathbf{x_r}, \\mathbf{x_s}, t) =\n        w(t) * \\int\\limits_V G(\\mathbf{x}, \\mathbf{x_s}, t)\n        G(\\mathbf{x_r}, \\mathbf{x}, t) m(\\mathbf{x})\\,\\mathrm{d}\\mathbf{x}\\end{align}\n\nwhere $m(\\mathbf{x})$ is the reflectivity\nat every location in the subsurface, $G(\\mathbf{x}, \\mathbf{x_s}, t)$\nand $G(\\mathbf{x_r}, \\mathbf{x}, t)$ are the Green's functions\nfrom source-to-subsurface-to-receiver and finally  $w(t)$ is the\nwavelet. Ultimately, while the Green's functions can be computed in many different\nways, solving this system of equations for the reflectivity model is what\nwe generally refer to as Least-squares migration (LSM).\n\nIn this tutorial we will consider the most simple scenario where we use an\neikonal solver to compute the Green's functions and show how we can use the\n:py:class:`pylops.waveeqprocessing.LSM` operator to perform LSM.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom scipy.sparse.linalg import lsqr\n\nimport pylops\n\nplt.close(\"all\")\nnp.random.seed(0)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To start we create a simple model with 2 interfaces\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "# Velocity Model\nnx, nz = 81, 60\ndx, dz = 4, 4\nx, z = np.arange(nx) * dx, np.arange(nz) * dz\nv0 = 1000  # initial velocity\nkv = 0.0  # gradient\nvel = np.outer(np.ones(nx), v0 + kv * z)\n\n# Reflectivity Model\nrefl = np.zeros((nx, nz))\nrefl[:, 30] = -1\nrefl[:, 50] = 0.5\n\n# Receivers\nnr = 11\nrx = np.linspace(10 * dx, (nx - 10) * dx, nr)\nrz = 20 * np.ones(nr)\nrecs = np.vstack((rx, rz))\ndr = recs[0, 1] - recs[0, 0]\n\n# Sources\nns = 10\nsx = np.linspace(dx * 10, (nx - 10) * dx, ns)\nsz = 10 * np.ones(ns)\nsources = np.vstack((sx, sz))\nds = sources[0, 1] - sources[0, 0]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "plt.figure(figsize=(10, 5))\nim = plt.imshow(vel.T, cmap=\"summer\", extent=(x[0], x[-1], z[-1], z[0]))\nplt.scatter(recs[0], recs[1], marker=\"v\", s=150, c=\"b\", edgecolors=\"k\")\nplt.scatter(sources[0], sources[1], marker=\"*\", s=150, c=\"r\", edgecolors=\"k\")\ncb = plt.colorbar(im)\ncb.set_label(\"[m/s]\")\nplt.axis(\"tight\")\nplt.xlabel(\"x [m]\"), plt.ylabel(\"z [m]\")\nplt.title(\"Velocity\")\nplt.xlim(x[0], x[-1])\nplt.tight_layout()\n\nplt.figure(figsize=(10, 5))\nim = plt.imshow(refl.T, cmap=\"gray\", extent=(x[0], x[-1], z[-1], z[0]))\nplt.scatter(recs[0], recs[1], marker=\"v\", s=150, c=\"b\", edgecolors=\"k\")\nplt.scatter(sources[0], sources[1], marker=\"*\", s=150, c=\"r\", edgecolors=\"k\")\nplt.colorbar(im)\nplt.axis(\"tight\")\nplt.xlabel(\"x [m]\"), plt.ylabel(\"z [m]\")\nplt.title(\"Reflectivity\")\nplt.xlim(x[0], x[-1])\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We can now create our LSM object and invert for the reflectivity using two\ndifferent solvers: :py:func:`scipy.sparse.linalg.lsqr` (LS solution) and\n:py:func:`pylops.optimization.sparsity.fista` (LS solution with sparse model).\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "nt = 651\ndt = 0.004\nt = np.arange(nt) * dt\nwav, wavt, wavc = pylops.utils.wavelets.ricker(t[:41], f0=20)\n\n\nlsm = pylops.waveeqprocessing.LSM(\n    z,\n    x,\n    t,\n    sources,\n    recs,\n    v0,\n    wav,\n    wavc,\n    mode=\"analytic\",\n    engine=\"numba\",\n)\n\nd = lsm.Demop * refl\n\nmadj = lsm.Demop.H * d\n\nminv = lsm.solve(d.ravel(), solver=lsqr, **dict(iter_lim=100))\nminv = minv.reshape(nx, nz)\n\nminv_sparse = lsm.solve(\n    d.ravel(), solver=pylops.optimization.sparsity.fista, **dict(eps=1e2, niter=100)\n)\nminv_sparse = minv_sparse.reshape(nx, nz)\n\n# demigration\nd = d.reshape(ns, nr, nt)\n\ndadj = lsm.Demop * madj  # (ns * nr, nt)\ndadj = dadj.reshape(ns, nr, nt)\n\ndinv = lsm.Demop * minv\ndinv = dinv.reshape(ns, nr, nt)\n\ndinv_sparse = lsm.Demop * minv_sparse\ndinv_sparse = dinv_sparse.reshape(ns, nr, nt)\n\n# sphinx_gallery_thumbnail_number = 2\nfig, axs = plt.subplots(2, 2, figsize=(10, 8))\naxs[0][0].imshow(refl.T, cmap=\"gray\", vmin=-1, vmax=1)\naxs[0][0].axis(\"tight\")\naxs[0][0].set_title(r\"$m$\")\naxs[0][1].imshow(madj.T, cmap=\"gray\", vmin=-madj.max(), vmax=madj.max())\naxs[0][1].set_title(r\"$m_{adj}$\")\naxs[0][1].axis(\"tight\")\naxs[1][0].imshow(minv.T, cmap=\"gray\", vmin=-1, vmax=1)\naxs[1][0].axis(\"tight\")\naxs[1][0].set_title(r\"$m_{inv}$\")\naxs[1][1].imshow(minv_sparse.T, cmap=\"gray\", vmin=-1, vmax=1)\naxs[1][1].axis(\"tight\")\naxs[1][1].set_title(r\"$m_{FISTA}$\")\nplt.tight_layout()\n\nfig, axs = plt.subplots(1, 4, figsize=(10, 4))\naxs[0].imshow(d[0, :, :300].T, cmap=\"gray\", vmin=-d.max(), vmax=d.max())\naxs[0].set_title(r\"$d$\")\naxs[0].axis(\"tight\")\naxs[1].imshow(dadj[0, :, :300].T, cmap=\"gray\", vmin=-dadj.max(), vmax=dadj.max())\naxs[1].set_title(r\"$d_{adj}$\")\naxs[1].axis(\"tight\")\naxs[2].imshow(dinv[0, :, :300].T, cmap=\"gray\", vmin=-d.max(), vmax=d.max())\naxs[2].set_title(r\"$d_{inv}$\")\naxs[2].axis(\"tight\")\naxs[3].imshow(dinv_sparse[0, :, :300].T, cmap=\"gray\", vmin=-d.max(), vmax=d.max())\naxs[3].set_title(r\"$d_{fista}$\")\naxs[3].axis(\"tight\")\nplt.tight_layout()\n\nfig, axs = plt.subplots(1, 4, figsize=(10, 4))\naxs[0].imshow(d[ns // 2, :, :300].T, cmap=\"gray\", vmin=-d.max(), vmax=d.max())\naxs[0].set_title(r\"$d$\")\naxs[0].axis(\"tight\")\naxs[1].imshow(dadj[ns // 2, :, :300].T, cmap=\"gray\", vmin=-dadj.max(), vmax=dadj.max())\naxs[1].set_title(r\"$d_{adj}$\")\naxs[1].axis(\"tight\")\naxs[2].imshow(dinv[ns // 2, :, :300].T, cmap=\"gray\", vmin=-d.max(), vmax=d.max())\naxs[2].set_title(r\"$d_{inv}$\")\naxs[2].axis(\"tight\")\naxs[3].imshow(dinv_sparse[ns // 2, :, :300].T, cmap=\"gray\", vmin=-d.max(), vmax=d.max())\naxs[3].set_title(r\"$d_{fista}$\")\naxs[3].axis(\"tight\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "This was just a short teaser, for a more advanced set of examples of 2D and\n3D traveltime-based LSM head over to this\n[notebook](https://github.com/mrava87/pylops_notebooks/blob/master/developement/LeastSquaresMigration.ipynb).\n\n"
      ]
    }
  ],
  "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
}