Time_RegGridProlongator_vs_scipyRegGridInterpolator.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Runtime of `RegGridProlongator` and `RegGridInterpolator`\n",
    "### ( `emg3d.solver.RegGridProlongator` and `scipy.interpolate.RegGridInterpolator`)\n",
    "\n",
    "For background and motivation see: https://mail.python.org/pipermail/scipy-dev/2019-May/023537.html.\n",
    "\n",
    "Requires [emg3d](https://empymod.github.io) v0.6.0 or higher:\n",
    "```\n",
    "conda install -c prisae emg3d\n",
    "```\n",
    "or\n",
    "```\n",
    "pip install emg3d\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import scipy.interpolate as si\n",
    "from emg3d import utils, solver"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def prolon_scipy(grid, cgrid, efield, cefield, x_points):\n",
    "    \"\"\"Calculate SciPy alternative.\"\"\"\n",
    "    \n",
    "    for ixc in range(cgrid.nCx):\n",
    "        \n",
    "        # Bilinear interpolation in the y-z plane\n",
    "        fn = si.RegularGridInterpolator(\n",
    "                (cgrid.vectorNy, cgrid.vectorNz), cefield.fx[ixc, :, :],\n",
    "                bounds_error=False, fill_value=None)\n",
    "        hh = fn(x_points).reshape(grid.vnEx[1:], order='F')\n",
    "\n",
    "        # Piecewise constant interpolation in x-direction\n",
    "        efield[2*ixc, :, :] += hh\n",
    "        efield[2*ixc+1, :, :] += hh\n",
    "        \n",
    "    return efield\n",
    "        \n",
    "\n",
    "def prolon_emg3d(grid, cgrid, efield, cefield, x_points):\n",
    "    \"\"\"Calculate emg3d alternative.\"\"\"\n",
    "    \n",
    "    fn = solver.RegularGridProlongator(cgrid.vectorNy, cgrid.vectorNz, x_points)\n",
    "    for ixc in range(cgrid.nCx):\n",
    "        \n",
    "        # Bilinear interpolation in the y-z plane\n",
    "        hh = fn(cefield.fx[ixc, :, :]).reshape(grid.vnEx[1:], order='F')\n",
    "        \n",
    "        # Piecewise constant interpolation in x-direction\n",
    "        efield[2*ixc, :, :] += hh\n",
    "        efield[2*ixc+1, :, :] += hh\n",
    "        \n",
    "    return efield"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "def test_prolong(n):\n",
    "    # Nr of cells of fine grid.\n",
    "    nx = 2**n\n",
    "\n",
    "    # Create fine grid.\n",
    "    hx = 50*np.ones(nx)\n",
    "    grid = utils.TensorMesh([hx, hx, hx], x0=np.array([0, 0, 0]))\n",
    "\n",
    "    # Create coarse grid.\n",
    "    chx = np.diff(grid.vectorNx[::2])\n",
    "    cgrid = utils.TensorMesh([chx, chx, chx], x0=np.array([0, 0, 0]))\n",
    "    print(f\"\\n  === n : {n} ===; fine: {grid.nC}, coarse: {cgrid.nC} (first SciPy, then emg3d)\")\n",
    "\n",
    "    # Create empty fine grid fields.\n",
    "    efield1 = utils.Field(grid)\n",
    "    efield2 = utils.Field(grid)\n",
    "\n",
    "    # Create coarse grid field with some values.\n",
    "    cefield = utils.Field(cgrid)\n",
    "    cefield.fx = np.arange(cefield.fx.size)\n",
    "    cefield.fx = 1j*np.arange(cefield.fx.size)/10\n",
    "\n",
    "    # Required interpolation points.\n",
    "    yz_points = solver._get_prolongation_coordinates(grid, 'y', 'z')\n",
    "\n",
    "    # Timeit\n",
    "    t1 = %timeit -o prolon_scipy(grid, cgrid, efield1.fx.copy(), cefield, yz_points)\n",
    "    t2 = %timeit -o prolon_emg3d(grid, cgrid, efield2.fx.copy(), cefield, yz_points)\n",
    "    print(f\"Speedup: {t1.best/t2.best:.2f}\")\n",
    "    \n",
    "    # Compare\n",
    "    out1 = prolon_scipy(grid, cgrid, efield1.fx.copy(), cefield, yz_points)\n",
    "    out2 = prolon_emg3d(grid, cgrid, efield2.fx.copy(), cefield, yz_points)\n",
    "    np.allclose(out1, out2, rtol=1e-7, atol=0)  # Same as assert_allclose"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "  === n : 2 ===; fine: 64, coarse: 8 (first SciPy, then emg3d)\n",
      "364 µs ± 22.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "192 µs ± 36.6 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      "Speedup: 2.14\n",
      "\n",
      "  === n : 3 ===; fine: 512, coarse: 64 (first SciPy, then emg3d)\n",
      "709 µs ± 18.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "289 µs ± 9.67 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "Speedup: 2.44\n",
      "\n",
      "  === n : 4 ===; fine: 4096, coarse: 512 (first SciPy, then emg3d)\n",
      "1.69 ms ± 30.2 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "578 µs ± 2.77 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "Speedup: 2.85\n",
      "\n",
      "  === n : 5 ===; fine: 32768, coarse: 4096 (first SciPy, then emg3d)\n",
      "4.91 ms ± 212 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "1.62 ms ± 29.5 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "Speedup: 2.99\n",
      "\n",
      "  === n : 6 ===; fine: 262144, coarse: 32768 (first SciPy, then emg3d)\n",
      "19.9 ms ± 735 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "8.37 ms ± 568 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)\n",
      "Speedup: 2.40\n",
      "\n",
      "  === n : 7 ===; fine: 2097152, coarse: 262144 (first SciPy, then emg3d)\n",
      "193 ms ± 10.3 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "106 ms ± 6.24 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Speedup: 1.81\n",
      "\n",
      "  === n : 8 ===; fine: 16777216, coarse: 2097152 (first SciPy, then emg3d)\n",
      "2.07 s ± 278 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "1.11 s ± 42.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "Speedup: 1.82\n"
     ]
    }
   ],
   "source": [
    "for n in range(2, 9):\n",
    "    test_prolong(n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table style='border: 3px solid #ddd;'>\n",
       "  <tr>\n",
       "     <td style='text-align: center; font-weight: bold; font-size: 1.2em; border: 2px solid #fff;' colspan='8'>Tue May 28 07:03:55 2019 CEST</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>Linux</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>OS</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>4</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>CPU(s)</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>1.16.2</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>numpy</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>1.2.1</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>scipy</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>0.43.1</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>numba</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>0.6.0</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>emg3d</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>7.4.0</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>IPython</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>3.0.3</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>matplotlib</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "     <td style='text-align: center; border: 2px solid #fff;' colspan='8'>3.7.3 (default, Mar 27 2019, 22:11:17) \n",
       "[GCC 7.3.0]</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "     <td style='text-align: center; background-color: #ddd;border: 2px solid #fff;' colspan='8'>Intel(R) Math Kernel Library Version 2019.0.3 Product Build 20190125 for Intel(R) 64 architecture applications</td>\n",
       "  </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "\n",
       "------------------------------------------------------\n",
       "  Tue May 28 07:03:55 2019 CEST\n",
       "\n",
       "          Linux : OS\n",
       "              4 : CPU(s)\n",
       "         1.16.2 : numpy\n",
       "          1.2.1 : scipy\n",
       "         0.43.1 : numba\n",
       "          0.6.0 : emg3d\n",
       "          7.4.0 : IPython\n",
       "          3.0.3 : matplotlib\n",
       "\n",
       "  3.7.3 (default, Mar 27 2019, 22:11:17)  [GCC\n",
       "  7.3.0]\n",
       "\n",
       "  Intel(R) Math Kernel Library Version 2019.0.3\n",
       "  Product Build 20190125 for Intel(R) 64\n",
       "  architecture applications\n",
       "------------------------------------------------------"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "utils.Versions()"
   ]
  }
 ],
 "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}