FastPatching2D.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "82506aa5",
   "metadata": {},
   "source": [
    "## Patching 2D"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "78f3b931",
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "%matplotlib inline\n",
    "\n",
    "import warnings\n",
    "warnings.filterwarnings('ignore')\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy as sp\n",
    "import pylops\n",
    "\n",
    "from pylops.signalprocessing.patch2d import patch2d_design\n",
    "from pylops.utils.describe import describe\n",
    "\n",
    "from pylops.signalprocessing import Patch2D\n",
    "from patch2dol import Patch2D as Patch2DOLD\n",
    "\n",
    "USE_CUPY = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "63b05131",
   "metadata": {},
   "outputs": [],
   "source": [
    "if USE_CUPY:\n",
    "    import cupy as np\n",
    "    from cupyx.profiler import benchmark\n",
    "    np_asarray = np.asarray\n",
    "    np_asnumpy = np.asnumpy\n",
    "    np_float = np.float32\n",
    "    np_floatc = np.complex64\n",
    "    mempool = np.get_default_memory_pool()\n",
    "    fftengine = 'numpy'\n",
    "else:\n",
    "    np_asarray = np.asarray\n",
    "    np_asnumpy = np.asarray\n",
    "    np_float = np.float64\n",
    "    np_floatc = np.complex128\n",
    "    fftengine = 'scipy'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "57ac53a4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def bench_Op(Op, x):\n",
    "    return Op @ x\n",
    "\n",
    "def bench_OpH(Op, x):\n",
    "    return Op.H @ x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "f797a610",
   "metadata": {},
   "outputs": [],
   "source": [
    "savetaper = True\n",
    "tapertype = 'cosine'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "b29494ab",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "WARNING: 30-16 windows required...\n",
      "WARNING: data wins - start:[  0  32  64  96 128 160 192 224 256 288 320 352 384 416 448 480 512 544\n",
      " 576 608 640 672 704 736 768 800 832 864 896 928], end:[ 42  74 106 138 170 202 234 266 298 330 362 394 426 458 490 522 554 586\n",
      " 618 650 682 714 746 778 810 842 874 906 938 970] / start:[  0  30  60  90 120 150 180 210 240 270 300 330 360 390 420 450], end:[ 34  64  94 124 154 184 214 244 274 304 334 364 394 424 454 484]\n",
      "WARNING: model wins - start:[   0   64  128  192  256  320  384  448  512  576  640  704  768  832\n",
      "  896  960 1024 1088 1152 1216 1280 1344 1408 1472 1536 1600 1664 1728\n",
      " 1792 1856], end:[  64  128  192  256  320  384  448  512  576  640  704  768  832  896\n",
      "  960 1024 1088 1152 1216 1280 1344 1408 1472 1536 1600 1664 1728 1792\n",
      " 1856 1920] / start:[  0  64 128 192 256 320 384 448 512 576 640 704 768 832 896 960], end:[  64  128  192  256  320  384  448  512  576  640  704  768  832  896\n",
      "  960 1024]\n"
     ]
    }
   ],
   "source": [
    "nwin = (42, 34)\n",
    "nover = (10, 4)\n",
    "nop = (64, 64)\n",
    "#dimsd = (200, 200) # small\n",
    "dimsd = (1000, 500) # large\n",
    "\n",
    "y = np.random.normal(0, 1, dimsd[0]*dimsd[1]).reshape(dimsd).astype(np_float)\n",
    "\n",
    "nwins, dims, mwin_inends, dwin_inends = patch2d_design(dimsd, nwin, nover, nop)\n",
    "\n",
    "# no operator broadcast\n",
    "Op = pylops.signalprocessing.FFT2D(nwin, nffts=nop, dtype=np_floatc)\n",
    "Slid = Patch2DOLD(Op.H, dims, dimsd, nwin, nover, nop, tapertype=tapertype)\n",
    "Slid1a = Patch2D(Op.H, dims, dimsd, nwin, nover, nop, tapertype=tapertype, savetaper=savetaper)\n",
    "\n",
    "# with operator broadcast\n",
    "Op = pylops.signalprocessing.FFT2D((*nwins, *nwin), nffts=nop, dtype=np_floatc)\n",
    "Slid1b = Patch2D(Op.H, dims, dimsd, nwin, nover, nop, tapertype=tapertype, savetaper=savetaper)\n",
    "\n",
    "x = Slid.H * y.ravel()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "dfcd25c4",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "True True\n",
      "True True\n"
     ]
    }
   ],
   "source": [
    "print(np.allclose(Slid @ x, Slid1a @ x), np.allclose(Slid.H @ y, Slid1a.H @ y))\n",
    "print(np.allclose(Slid @ x, Slid1b @ x), np.allclose(Slid.H @ y, Slid1b.H @ y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "7eff976a",
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "92 ms ± 2.11 ms per loop (mean ± std. dev. of 50 runs, 5 loops each)\n",
      "47.4 ms ± 3.35 ms per loop (mean ± std. dev. of 50 runs, 5 loops each)\n",
      "31.2 ms ± 1.89 ms per loop (mean ± std. dev. of 50 runs, 5 loops each)\n"
     ]
    }
   ],
   "source": [
    "if not USE_CUPY:\n",
    "    %timeit -n 5 -r 50 Slid * x # OLD\n",
    "    %timeit -n 5 -r 50 Slid1a * x # NEW\n",
    "    %timeit -n 5 -r 50 Slid1b * x # NEW with Op broadcasted\n",
    "else:\n",
    "    print(benchmark(bench_Op, (Slid, x,), n_repeat=500))\n",
    "    print(benchmark(bench_Op, (Slid1a, x,), n_repeat=500))\n",
    "    print(benchmark(bench_Op, (Slid1b, x,), n_repeat=500))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "05f63666",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "50.5 ms ± 2.17 ms per loop (mean ± std. dev. of 50 runs, 5 loops each)\n",
      "37.5 ms ± 2 ms per loop (mean ± std. dev. of 50 runs, 5 loops each)\n",
      "24.1 ms ± 1.03 ms per loop (mean ± std. dev. of 50 runs, 5 loops each)\n"
     ]
    }
   ],
   "source": [
    "if not USE_CUPY:\n",
    "    %timeit -n 5 -r 50 Slid.H * y # OLD\n",
    "    %timeit -n 5 -r 50 Slid1a.H * y # NEW\n",
    "    %timeit -n 5 -r 50 Slid1b.H * y # NEW with Op broadcasted\n",
    "else:\n",
    "    print(benchmark(bench_OpH, (Slid, y,), n_repeat=500))\n",
    "    print(benchmark(bench_OpH, (Slid1a, y,), n_repeat=500))\n",
    "    print(benchmark(bench_OpH, (Slid1b, y,), n_repeat=500))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.18"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}