{
  "cells": [
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "%matplotlib inline"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "\n# Seislet transform\nThis example shows how to use the :py:class:`pylops.signalprocessing.Seislet`\noperator. This operator the forward, adjoint and inverse Seislet transform\nthat is a modification of the well-know Wavelet transform where local slopes\nare used in the prediction and update steps to further improve the prediction\nof a trace from its previous (or subsequent) one and reduce the amount of\ninformation passed to the subsequent scale. While this transform was initially\ndeveloped in the context of processing and compression of seismic data, it is\nalso suitable to any other oscillatory dataset such as GPR or Acoustic\nrecordings.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "import matplotlib.pyplot as plt\nimport numpy as np\nfrom matplotlib.ticker import MaxNLocator\nfrom mpl_toolkits.axes_grid1 import make_axes_locatable\n\nimport pylops\n\nplt.close(\"all\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "In this example we use the same benchmark\n[dataset](http://ahay.org/blog/2014/10/08/program-of-the-month-sfsigmoid/)\nthat was used in the original paper describing the Seislet transform. First,\nlocal slopes are estimated using\n:py:func:`pylops.utils.signalprocessing.slope_estimate`.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "inputfile = \"../testdata/sigmoid.npz\"\n\nd = np.load(inputfile)\nd = d[\"sigmoid\"]\nnx, nt = d.shape\ndx, dt = 0.008, 0.004\nx, t = np.arange(nx) * dx, np.arange(nt) * dt\n\n# slope estimation\nslope, _ = pylops.utils.signalprocessing.slope_estimate(d.T, dt, dx, smooth=2.5)\nslope = -slope.T  # t-axis points down, reshape\n# clip slopes above 80\u00b0\npmax = np.arctan(80 * np.pi / 180)\nslope[slope > pmax] = pmax\nslope[slope < -pmax] = -pmax"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "clip = 0.5 * np.max(np.abs(d))\nclip_s = min(pmax, np.max(np.abs(slope)))\nopts = dict(aspect=2, extent=(x[0], x[-1], t[-1], t[0]))\n\nfig, axs = plt.subplots(1, 2, figsize=(14, 7), sharey=True, sharex=True)\naxs[0].imshow(d.T, cmap=\"gray\", vmin=-clip, vmax=clip, **opts)\naxs[0].set(xlabel=\"Position [km]\", ylabel=\"Time [s]\", title=\"Data\")\n\nim = axs[1].imshow(slope.T, cmap=\"RdBu_r\", vmin=-clip_s, vmax=clip_s, **opts)\naxs[1].set(xlabel=\"Position [km]\", title=\"Slopes\")\nfig.tight_layout()\n\npos = axs[1].get_position()\ncbpos = [\n    pos.x0 + 0.1 * pos.width,\n    pos.y0 + 0.9 * pos.height,\n    0.8 * pos.width,\n    0.05 * pos.height,\n]\ncax = fig.add_axes(cbpos)\ncb = fig.colorbar(im, cax=cax, orientation=\"horizontal\")\ncb.set_label(\"[s/km]\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Next the Seislet transform is computed.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Sop = pylops.signalprocessing.Seislet(slope, sampling=(dx, dt))\n\nseis = Sop * d\n\nnlevels_max = int(np.log2(nx))\nlevels_size = np.flip(np.array([2**i for i in range(nlevels_max)]))\nlevels_cum = np.cumsum(levels_size)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(14, 6))\nim = ax.imshow(\n    seis.T,\n    cmap=\"gray\",\n    vmin=-clip,\n    vmax=clip,\n    aspect=\"auto\",\n    interpolation=\"none\",\n    extent=(1, seis.shape[0], t[-1], t[0]),\n)\nax.xaxis.set_major_locator(MaxNLocator(nbins=20, integer=True))\nfor level in levels_cum:\n    ax.axvline(level + 0.5, color=\"w\")\nax.set(xlabel=\"Scale\", ylabel=\"Time [s]\", title=\"Seislet transform\")\ncax = make_axes_locatable(ax).append_axes(\"right\", size=\"2%\", pad=0.1)\ncb = fig.colorbar(im, cax=cax, orientation=\"vertical\")\ncb.formatter.set_powerlimits((0, 0))\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "We may also stretch the finer scales to be the width of the image\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(2, nlevels_max // 2, figsize=(14, 7), sharex=True, sharey=True)\nfor i, ax in enumerate(axs.ravel()[:-1]):\n    curdata = seis[levels_cum[i] : levels_cum[i + 1], :].T\n    vmax = np.max(np.abs(curdata))\n    ax.imshow(curdata, vmin=-vmax, vmax=vmax, cmap=\"gray\", interpolation=\"none\", **opts)\n    ax.set(title=f\"Scale {i+1}\")\n    if i + 1 > nlevels_max // 2:\n        ax.set(xlabel=\"Position [km]\")\ncurdata = seis[levels_cum[-1] :, :].T\nvmax = np.max(np.abs(curdata))\naxs[-1, -1].imshow(\n    curdata, vmin=-vmax, vmax=vmax, cmap=\"gray\", interpolation=\"none\", **opts\n)\naxs[0, 0].set(ylabel=\"Time [s]\")\naxs[1, 0].set(ylabel=\"Time [s]\")\naxs[-1, -1].set(xlabel=\"Position [km]\", title=f\"Scale {nlevels_max}\")\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "As a comparison we also compute the Seislet transform fixing slopes to zero.\nThis way we turn the Seislet tranform into a basic 1D Wavelet transform\nperformed over the spatial axis.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "Wop = pylops.signalprocessing.Seislet(np.zeros_like(slope), sampling=(dx, dt))\ndwt = Wop * d"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots(figsize=(14, 6))\nim = ax.imshow(\n    dwt.T,\n    cmap=\"gray\",\n    vmin=-clip,\n    vmax=clip,\n    aspect=\"auto\",\n    interpolation=\"none\",\n    extent=(1, dwt.shape[0], t[-1], t[0]),\n)\nax.xaxis.set_major_locator(MaxNLocator(nbins=20, integer=True))\nfor level in levels_cum:\n    ax.axvline(level + 0.5, color=\"w\")\nax.set(xlabel=\"Scale\", ylabel=\"Time [s]\", title=\"Wavelet transform\")\ncax = make_axes_locatable(ax).append_axes(\"right\", size=\"2%\", pad=0.1)\ncb = fig.colorbar(im, cax=cax, orientation=\"vertical\")\ncb.formatter.set_powerlimits((0, 0))\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Again, we may decompress the finer scales\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, axs = plt.subplots(2, nlevels_max // 2, figsize=(14, 7), sharex=True, sharey=True)\nfor i, ax in enumerate(axs.ravel()[:-1]):\n    curdata = dwt[levels_cum[i] : levels_cum[i + 1], :].T\n    vmax = np.max(np.abs(curdata))\n    ax.imshow(curdata, vmin=-vmax, vmax=vmax, cmap=\"gray\", interpolation=\"none\", **opts)\n    ax.set(title=f\"Scale {i+1}\")\n    if i + 1 > nlevels_max // 2:\n        ax.set(xlabel=\"Position [km]\")\ncurdata = dwt[levels_cum[-1] :, :].T\nvmax = np.max(np.abs(curdata))\naxs[-1, -1].imshow(\n    curdata, vmin=-vmax, vmax=vmax, cmap=\"gray\", interpolation=\"none\", **opts\n)\naxs[0, 0].set(ylabel=\"Time [s]\")\naxs[1, 0].set(ylabel=\"Time [s]\")\naxs[-1, -1].set(xlabel=\"Position [km]\", title=f\"Scale {nlevels_max}\")\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "Finally we evaluate the compression capabilities of the Seislet transform\ncompared to the 1D Wavelet transform. We zero-out all but the strongest 25%\nof the components. We perform the inverse transforms and assess the\ncompression error.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "perc = 0.25\nseis_strong_idx = np.argsort(-np.abs(seis.ravel()))\ndwt_strong_idx = np.argsort(-np.abs(dwt.ravel()))\nseis_strong = np.abs(seis.ravel())[seis_strong_idx]\ndwt_strong = np.abs(dwt.ravel())[dwt_strong_idx]"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "fig, ax = plt.subplots()\nax.plot(range(1, len(seis_strong) + 1), seis_strong / seis_strong[0], label=\"Seislet\")\nax.plot(\n    range(1, len(dwt_strong) + 1), dwt_strong / dwt_strong[0], \"--\", label=\"Wavelet\"\n)\nax.set(xlabel=\"n\", ylabel=\"Coefficient strength [%]\", title=\"Transform Coefficients\")\nax.axvline(np.rint(len(seis_strong) * perc), color=\"k\", label=f\"{100*perc:.0f}%\")\nax.legend()\nfig.tight_layout()"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "seis1 = np.zeros_like(seis.ravel())\nseis_strong_idx = seis_strong_idx[: int(np.rint(len(seis_strong) * perc))]\nseis1[seis_strong_idx] = seis.ravel()[seis_strong_idx]\nd_seis = Sop.inverse(seis1).reshape(Sop.dims)\n\ndwt1 = np.zeros_like(dwt.ravel())\ndwt_strong_idx = dwt_strong_idx[: int(np.rint(len(dwt_strong) * perc))]\ndwt1[dwt_strong_idx] = dwt.ravel()[dwt_strong_idx]\nd_dwt = Wop.inverse(dwt1).reshape(Wop.dims)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "opts.update(dict(cmap=\"gray\", vmin=-clip, vmax=clip))\nfig, axs = plt.subplots(2, 3, figsize=(14, 7), sharex=True, sharey=True)\naxs[0, 0].imshow(d.T, **opts)\naxs[0, 0].set(title=\"Data\")\naxs[0, 1].imshow(d_seis.T, **opts)\naxs[0, 1].set(title=f\"Rec. from Seislet ({100*perc:.0f}% of coeffs.)\")\naxs[0, 2].imshow((d - d_seis).T, **opts)\naxs[0, 2].set(title=\"Error from Seislet Rec.\")\naxs[1, 0].imshow(d.T, **opts)\naxs[1, 0].set(ylabel=\"Time [s]\", title=\"Data [Repeat]\")\naxs[1, 1].imshow(d_dwt.T, **opts)\naxs[1, 1].set(title=f\"Rec. from Wavelet ({100*perc:.0f}% of coeffs.)\")\naxs[1, 2].imshow((d - d_dwt).T, **opts)\naxs[1, 2].set(title=\"Error from Wavelet Rec.\")\nfor i in range(3):\n    axs[1, i].set(xlabel=\"Position [km]\")\nplt.tight_layout()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {},
      "source": [
        "To conclude it is worth noting that the Seislet transform, differently to the\nWavelet transform, is not orthogonal: in other words, its adjoint and\ninverse are not equivalent. While we have used the forward and inverse\ntransformations, when used as linear operator in composition with other\noperators, the Seislet transform requires the adjoint be defined and that it\nalso passes the dot-test pair that is. As shown below, this is the case\nwhen using the implementation in the PyLops package.\n\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "collapsed": false
      },
      "outputs": [],
      "source": [
        "pylops.utils.dottest(Sop, verb=True)"
      ]
    }
  ],
  "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
}