Complex-Sqrt_Numpy-Numexpr-Numba-Cython.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simple Test between NumPy, Numexpr, Numba, and Cython\n",
    "\n",
    "## for complex square root of big matrices\n",
    "\n",
    "The calculation is\n",
    "$$\\sqrt{A\\ L^2 + B}\\ ,$$\n",
    "where\n",
    "- $A$ and $B$ are complex matrices (complex128) with dimension (#frequencies, #layers), and\n",
    "- L is a real matrix (floats64) with dimension (#offsets, #wavenumbers).\n",
    "\n",
    "Their dimensions are broadcasted (`A[:, None, :, None]`; `B[:, None, :, None]`; and `L[None, :, None, :]`) so that the result is a complex matrix of dimension (#frequencies, #offsets, #layers, #wavenumbers), which can become quite big.\n",
    "\n",
    "Here I compare my go at `Numpy`, `Numexpr`, `Numba`, and `Cython` implementations. The `Numexpr` was very fast; however, it used 4 threads by default, whereas all others only used 1 thread. So I restricted `Numexpr` to 1 thread to, which then turns out to be the slowest of the four.\n",
    "\n",
    "In this comparison `Numba` and `Cython` are about the same, give or take (run the notebook a few times), where the `Numba` implementation is slightly easier.\n",
    "\n",
    "**Please let me know if you know how to improve any of these!**\n",
    "\n",
    "**Note:** GitHub does not render the Cython output nicely. Have a look at nbviewer to see a nicer rendering of the notebook:  \n",
    "http://nbviewer.jupyter.org/github/prisae/tmp-share/blob/master/Complex-Sqrt_Numpy-Numexpr-Numba-Cython.ipynb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numba\n",
    "import cython\n",
    "import numexpr\n",
    "import numpy as np\n",
    "\n",
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## NumPy\n",
    "\n",
    "Numpy version to check result and compare times"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "def test_numpy(a, b, l):\n",
    "    return np.sqrt(a*l**2 + b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numexpr\n",
    "\n",
    "Numexpr seems to be slower, but is really good if run in parallel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "numexpr.set_num_threads(1)  # Set nr of threads to 1 for a fair comparison\n",
    "def test_numexpr(a, b, l):\n",
    "    return numexpr.evaluate(\"sqrt(a*l**2 + b)\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba @vectorize\n",
    "\n",
    "This is exactly the same function as with NumPy, just added the @vectorize decorater."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "@numba.vectorize('c16(c16, c16, f8)')\n",
    "def test_numba_v(a, b, l):\n",
    "    return np.sqrt(a*l**2 + b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Numba @njit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "@numba.njit\n",
    "def test_numba_n(a, b, l, out):\n",
    "    for nf in range(out.shape[0]):\n",
    "        for no in range(out.shape[1]):\n",
    "            for nl in range(out.shape[2]):\n",
    "                for ni in range(out.shape[3]):\n",
    "                    out[nf, no, nl, ni] = np.sqrt(a[nf, 0, nl, 0] * l[0, no, 0, ni] ** 2 + b[nf, 0, nl, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<!DOCTYPE html>\n",
       "<!-- Generated by Cython 0.28.3 -->\n",
       "<html>\n",
       "<head>\n",
       "    <meta http-equiv=\"Content-Type\" content=\"text/html; charset=utf-8\" />\n",
       "    <title>Cython: _cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d.pyx</title>\n",
       "    <style type=\"text/css\">\n",
       "    \n",
       "body.cython { font-family: courier; font-size: 12; }\n",
       "\n",
       ".cython.tag  {  }\n",
       ".cython.line { margin: 0em }\n",
       ".cython.code { font-size: 9; color: #444444; display: none; margin: 0px 0px 0px 8px; border-left: 8px none; }\n",
       "\n",
       ".cython.line .run { background-color: #B0FFB0; }\n",
       ".cython.line .mis { background-color: #FFB0B0; }\n",
       ".cython.code.run  { border-left: 8px solid #B0FFB0; }\n",
       ".cython.code.mis  { border-left: 8px solid #FFB0B0; }\n",
       "\n",
       ".cython.code .py_c_api  { color: red; }\n",
       ".cython.code .py_macro_api  { color: #FF7000; }\n",
       ".cython.code .pyx_c_api  { color: #FF3000; }\n",
       ".cython.code .pyx_macro_api  { color: #FF7000; }\n",
       ".cython.code .refnanny  { color: #FFA000; }\n",
       ".cython.code .trace  { color: #FFA000; }\n",
       ".cython.code .error_goto  { color: #FFA000; }\n",
       "\n",
       ".cython.code .coerce  { color: #008000; border: 1px dotted #008000 }\n",
       ".cython.code .py_attr { color: #FF0000; font-weight: bold; }\n",
       ".cython.code .c_attr  { color: #0000FF; }\n",
       ".cython.code .py_call { color: #FF0000; font-weight: bold; }\n",
       ".cython.code .c_call  { color: #0000FF; }\n",
       "\n",
       ".cython.score-0 {background-color: #FFFFff;}\n",
       ".cython.score-1 {background-color: #FFFFe7;}\n",
       ".cython.score-2 {background-color: #FFFFd4;}\n",
       ".cython.score-3 {background-color: #FFFFc4;}\n",
       ".cython.score-4 {background-color: #FFFFb6;}\n",
       ".cython.score-5 {background-color: #FFFFaa;}\n",
       ".cython.score-6 {background-color: #FFFF9f;}\n",
       ".cython.score-7 {background-color: #FFFF96;}\n",
       ".cython.score-8 {background-color: #FFFF8d;}\n",
       ".cython.score-9 {background-color: #FFFF86;}\n",
       ".cython.score-10 {background-color: #FFFF7f;}\n",
       ".cython.score-11 {background-color: #FFFF79;}\n",
       ".cython.score-12 {background-color: #FFFF73;}\n",
       ".cython.score-13 {background-color: #FFFF6e;}\n",
       ".cython.score-14 {background-color: #FFFF6a;}\n",
       ".cython.score-15 {background-color: #FFFF66;}\n",
       ".cython.score-16 {background-color: #FFFF62;}\n",
       ".cython.score-17 {background-color: #FFFF5e;}\n",
       ".cython.score-18 {background-color: #FFFF5b;}\n",
       ".cython.score-19 {background-color: #FFFF57;}\n",
       ".cython.score-20 {background-color: #FFFF55;}\n",
       ".cython.score-21 {background-color: #FFFF52;}\n",
       ".cython.score-22 {background-color: #FFFF4f;}\n",
       ".cython.score-23 {background-color: #FFFF4d;}\n",
       ".cython.score-24 {background-color: #FFFF4b;}\n",
       ".cython.score-25 {background-color: #FFFF48;}\n",
       ".cython.score-26 {background-color: #FFFF46;}\n",
       ".cython.score-27 {background-color: #FFFF44;}\n",
       ".cython.score-28 {background-color: #FFFF43;}\n",
       ".cython.score-29 {background-color: #FFFF41;}\n",
       ".cython.score-30 {background-color: #FFFF3f;}\n",
       ".cython.score-31 {background-color: #FFFF3e;}\n",
       ".cython.score-32 {background-color: #FFFF3c;}\n",
       ".cython.score-33 {background-color: #FFFF3b;}\n",
       ".cython.score-34 {background-color: #FFFF39;}\n",
       ".cython.score-35 {background-color: #FFFF38;}\n",
       ".cython.score-36 {background-color: #FFFF37;}\n",
       ".cython.score-37 {background-color: #FFFF36;}\n",
       ".cython.score-38 {background-color: #FFFF35;}\n",
       ".cython.score-39 {background-color: #FFFF34;}\n",
       ".cython.score-40 {background-color: #FFFF33;}\n",
       ".cython.score-41 {background-color: #FFFF32;}\n",
       ".cython.score-42 {background-color: #FFFF31;}\n",
       ".cython.score-43 {background-color: #FFFF30;}\n",
       ".cython.score-44 {background-color: #FFFF2f;}\n",
       ".cython.score-45 {background-color: #FFFF2e;}\n",
       ".cython.score-46 {background-color: #FFFF2d;}\n",
       ".cython.score-47 {background-color: #FFFF2c;}\n",
       ".cython.score-48 {background-color: #FFFF2b;}\n",
       ".cython.score-49 {background-color: #FFFF2b;}\n",
       ".cython.score-50 {background-color: #FFFF2a;}\n",
       ".cython.score-51 {background-color: #FFFF29;}\n",
       ".cython.score-52 {background-color: #FFFF29;}\n",
       ".cython.score-53 {background-color: #FFFF28;}\n",
       ".cython.score-54 {background-color: #FFFF27;}\n",
       ".cython.score-55 {background-color: #FFFF27;}\n",
       ".cython.score-56 {background-color: #FFFF26;}\n",
       ".cython.score-57 {background-color: #FFFF26;}\n",
       ".cython.score-58 {background-color: #FFFF25;}\n",
       ".cython.score-59 {background-color: #FFFF24;}\n",
       ".cython.score-60 {background-color: #FFFF24;}\n",
       ".cython.score-61 {background-color: #FFFF23;}\n",
       ".cython.score-62 {background-color: #FFFF23;}\n",
       ".cython.score-63 {background-color: #FFFF22;}\n",
       ".cython.score-64 {background-color: #FFFF22;}\n",
       ".cython.score-65 {background-color: #FFFF22;}\n",
       ".cython.score-66 {background-color: #FFFF21;}\n",
       ".cython.score-67 {background-color: #FFFF21;}\n",
       ".cython.score-68 {background-color: #FFFF20;}\n",
       ".cython.score-69 {background-color: #FFFF20;}\n",
       ".cython.score-70 {background-color: #FFFF1f;}\n",
       ".cython.score-71 {background-color: #FFFF1f;}\n",
       ".cython.score-72 {background-color: #FFFF1f;}\n",
       ".cython.score-73 {background-color: #FFFF1e;}\n",
       ".cython.score-74 {background-color: #FFFF1e;}\n",
       ".cython.score-75 {background-color: #FFFF1e;}\n",
       ".cython.score-76 {background-color: #FFFF1d;}\n",
       ".cython.score-77 {background-color: #FFFF1d;}\n",
       ".cython.score-78 {background-color: #FFFF1c;}\n",
       ".cython.score-79 {background-color: #FFFF1c;}\n",
       ".cython.score-80 {background-color: #FFFF1c;}\n",
       ".cython.score-81 {background-color: #FFFF1c;}\n",
       ".cython.score-82 {background-color: #FFFF1b;}\n",
       ".cython.score-83 {background-color: #FFFF1b;}\n",
       ".cython.score-84 {background-color: #FFFF1b;}\n",
       ".cython.score-85 {background-color: #FFFF1a;}\n",
       ".cython.score-86 {background-color: #FFFF1a;}\n",
       ".cython.score-87 {background-color: #FFFF1a;}\n",
       ".cython.score-88 {background-color: #FFFF1a;}\n",
       ".cython.score-89 {background-color: #FFFF19;}\n",
       ".cython.score-90 {background-color: #FFFF19;}\n",
       ".cython.score-91 {background-color: #FFFF19;}\n",
       ".cython.score-92 {background-color: #FFFF19;}\n",
       ".cython.score-93 {background-color: #FFFF18;}\n",
       ".cython.score-94 {background-color: #FFFF18;}\n",
       ".cython.score-95 {background-color: #FFFF18;}\n",
       ".cython.score-96 {background-color: #FFFF18;}\n",
       ".cython.score-97 {background-color: #FFFF17;}\n",
       ".cython.score-98 {background-color: #FFFF17;}\n",
       ".cython.score-99 {background-color: #FFFF17;}\n",
       ".cython.score-100 {background-color: #FFFF17;}\n",
       ".cython.score-101 {background-color: #FFFF16;}\n",
       ".cython.score-102 {background-color: #FFFF16;}\n",
       ".cython.score-103 {background-color: #FFFF16;}\n",
       ".cython.score-104 {background-color: #FFFF16;}\n",
       ".cython.score-105 {background-color: #FFFF16;}\n",
       ".cython.score-106 {background-color: #FFFF15;}\n",
       ".cython.score-107 {background-color: #FFFF15;}\n",
       ".cython.score-108 {background-color: #FFFF15;}\n",
       ".cython.score-109 {background-color: #FFFF15;}\n",
       ".cython.score-110 {background-color: #FFFF15;}\n",
       ".cython.score-111 {background-color: #FFFF15;}\n",
       ".cython.score-112 {background-color: #FFFF14;}\n",
       ".cython.score-113 {background-color: #FFFF14;}\n",
       ".cython.score-114 {background-color: #FFFF14;}\n",
       ".cython.score-115 {background-color: #FFFF14;}\n",
       ".cython.score-116 {background-color: #FFFF14;}\n",
       ".cython.score-117 {background-color: #FFFF14;}\n",
       ".cython.score-118 {background-color: #FFFF13;}\n",
       ".cython.score-119 {background-color: #FFFF13;}\n",
       ".cython.score-120 {background-color: #FFFF13;}\n",
       ".cython.score-121 {background-color: #FFFF13;}\n",
       ".cython.score-122 {background-color: #FFFF13;}\n",
       ".cython.score-123 {background-color: #FFFF13;}\n",
       ".cython.score-124 {background-color: #FFFF13;}\n",
       ".cython.score-125 {background-color: #FFFF12;}\n",
       ".cython.score-126 {background-color: #FFFF12;}\n",
       ".cython.score-127 {background-color: #FFFF12;}\n",
       ".cython.score-128 {background-color: #FFFF12;}\n",
       ".cython.score-129 {background-color: #FFFF12;}\n",
       ".cython.score-130 {background-color: #FFFF12;}\n",
       ".cython.score-131 {background-color: #FFFF12;}\n",
       ".cython.score-132 {background-color: #FFFF11;}\n",
       ".cython.score-133 {background-color: #FFFF11;}\n",
       ".cython.score-134 {background-color: #FFFF11;}\n",
       ".cython.score-135 {background-color: #FFFF11;}\n",
       ".cython.score-136 {background-color: #FFFF11;}\n",
       ".cython.score-137 {background-color: #FFFF11;}\n",
       ".cython.score-138 {background-color: #FFFF11;}\n",
       ".cython.score-139 {background-color: #FFFF11;}\n",
       ".cython.score-140 {background-color: #FFFF11;}\n",
       ".cython.score-141 {background-color: #FFFF10;}\n",
       ".cython.score-142 {background-color: #FFFF10;}\n",
       ".cython.score-143 {background-color: #FFFF10;}\n",
       ".cython.score-144 {background-color: #FFFF10;}\n",
       ".cython.score-145 {background-color: #FFFF10;}\n",
       ".cython.score-146 {background-color: #FFFF10;}\n",
       ".cython.score-147 {background-color: #FFFF10;}\n",
       ".cython.score-148 {background-color: #FFFF10;}\n",
       ".cython.score-149 {background-color: #FFFF10;}\n",
       ".cython.score-150 {background-color: #FFFF0f;}\n",
       ".cython.score-151 {background-color: #FFFF0f;}\n",
       ".cython.score-152 {background-color: #FFFF0f;}\n",
       ".cython.score-153 {background-color: #FFFF0f;}\n",
       ".cython.score-154 {background-color: #FFFF0f;}\n",
       ".cython.score-155 {background-color: #FFFF0f;}\n",
       ".cython.score-156 {background-color: #FFFF0f;}\n",
       ".cython.score-157 {background-color: #FFFF0f;}\n",
       ".cython.score-158 {background-color: #FFFF0f;}\n",
       ".cython.score-159 {background-color: #FFFF0f;}\n",
       ".cython.score-160 {background-color: #FFFF0f;}\n",
       ".cython.score-161 {background-color: #FFFF0e;}\n",
       ".cython.score-162 {background-color: #FFFF0e;}\n",
       ".cython.score-163 {background-color: #FFFF0e;}\n",
       ".cython.score-164 {background-color: #FFFF0e;}\n",
       ".cython.score-165 {background-color: #FFFF0e;}\n",
       ".cython.score-166 {background-color: #FFFF0e;}\n",
       ".cython.score-167 {background-color: #FFFF0e;}\n",
       ".cython.score-168 {background-color: #FFFF0e;}\n",
       ".cython.score-169 {background-color: #FFFF0e;}\n",
       ".cython.score-170 {background-color: #FFFF0e;}\n",
       ".cython.score-171 {background-color: #FFFF0e;}\n",
       ".cython.score-172 {background-color: #FFFF0e;}\n",
       ".cython.score-173 {background-color: #FFFF0d;}\n",
       ".cython.score-174 {background-color: #FFFF0d;}\n",
       ".cython.score-175 {background-color: #FFFF0d;}\n",
       ".cython.score-176 {background-color: #FFFF0d;}\n",
       ".cython.score-177 {background-color: #FFFF0d;}\n",
       ".cython.score-178 {background-color: #FFFF0d;}\n",
       ".cython.score-179 {background-color: #FFFF0d;}\n",
       ".cython.score-180 {background-color: #FFFF0d;}\n",
       ".cython.score-181 {background-color: #FFFF0d;}\n",
       ".cython.score-182 {background-color: #FFFF0d;}\n",
       ".cython.score-183 {background-color: #FFFF0d;}\n",
       ".cython.score-184 {background-color: #FFFF0d;}\n",
       ".cython.score-185 {background-color: #FFFF0d;}\n",
       ".cython.score-186 {background-color: #FFFF0d;}\n",
       ".cython.score-187 {background-color: #FFFF0c;}\n",
       ".cython.score-188 {background-color: #FFFF0c;}\n",
       ".cython.score-189 {background-color: #FFFF0c;}\n",
       ".cython.score-190 {background-color: #FFFF0c;}\n",
       ".cython.score-191 {background-color: #FFFF0c;}\n",
       ".cython.score-192 {background-color: #FFFF0c;}\n",
       ".cython.score-193 {background-color: #FFFF0c;}\n",
       ".cython.score-194 {background-color: #FFFF0c;}\n",
       ".cython.score-195 {background-color: #FFFF0c;}\n",
       ".cython.score-196 {background-color: #FFFF0c;}\n",
       ".cython.score-197 {background-color: #FFFF0c;}\n",
       ".cython.score-198 {background-color: #FFFF0c;}\n",
       ".cython.score-199 {background-color: #FFFF0c;}\n",
       ".cython.score-200 {background-color: #FFFF0c;}\n",
       ".cython.score-201 {background-color: #FFFF0c;}\n",
       ".cython.score-202 {background-color: #FFFF0c;}\n",
       ".cython.score-203 {background-color: #FFFF0b;}\n",
       ".cython.score-204 {background-color: #FFFF0b;}\n",
       ".cython.score-205 {background-color: #FFFF0b;}\n",
       ".cython.score-206 {background-color: #FFFF0b;}\n",
       ".cython.score-207 {background-color: #FFFF0b;}\n",
       ".cython.score-208 {background-color: #FFFF0b;}\n",
       ".cython.score-209 {background-color: #FFFF0b;}\n",
       ".cython.score-210 {background-color: #FFFF0b;}\n",
       ".cython.score-211 {background-color: #FFFF0b;}\n",
       ".cython.score-212 {background-color: #FFFF0b;}\n",
       ".cython.score-213 {background-color: #FFFF0b;}\n",
       ".cython.score-214 {background-color: #FFFF0b;}\n",
       ".cython.score-215 {background-color: #FFFF0b;}\n",
       ".cython.score-216 {background-color: #FFFF0b;}\n",
       ".cython.score-217 {background-color: #FFFF0b;}\n",
       ".cython.score-218 {background-color: #FFFF0b;}\n",
       ".cython.score-219 {background-color: #FFFF0b;}\n",
       ".cython.score-220 {background-color: #FFFF0b;}\n",
       ".cython.score-221 {background-color: #FFFF0b;}\n",
       ".cython.score-222 {background-color: #FFFF0a;}\n",
       ".cython.score-223 {background-color: #FFFF0a;}\n",
       ".cython.score-224 {background-color: #FFFF0a;}\n",
       ".cython.score-225 {background-color: #FFFF0a;}\n",
       ".cython.score-226 {background-color: #FFFF0a;}\n",
       ".cython.score-227 {background-color: #FFFF0a;}\n",
       ".cython.score-228 {background-color: #FFFF0a;}\n",
       ".cython.score-229 {background-color: #FFFF0a;}\n",
       ".cython.score-230 {background-color: #FFFF0a;}\n",
       ".cython.score-231 {background-color: #FFFF0a;}\n",
       ".cython.score-232 {background-color: #FFFF0a;}\n",
       ".cython.score-233 {background-color: #FFFF0a;}\n",
       ".cython.score-234 {background-color: #FFFF0a;}\n",
       ".cython.score-235 {background-color: #FFFF0a;}\n",
       ".cython.score-236 {background-color: #FFFF0a;}\n",
       ".cython.score-237 {background-color: #FFFF0a;}\n",
       ".cython.score-238 {background-color: #FFFF0a;}\n",
       ".cython.score-239 {background-color: #FFFF0a;}\n",
       ".cython.score-240 {background-color: #FFFF0a;}\n",
       ".cython.score-241 {background-color: #FFFF0a;}\n",
       ".cython.score-242 {background-color: #FFFF0a;}\n",
       ".cython.score-243 {background-color: #FFFF0a;}\n",
       ".cython.score-244 {background-color: #FFFF0a;}\n",
       ".cython.score-245 {background-color: #FFFF0a;}\n",
       ".cython.score-246 {background-color: #FFFF09;}\n",
       ".cython.score-247 {background-color: #FFFF09;}\n",
       ".cython.score-248 {background-color: #FFFF09;}\n",
       ".cython.score-249 {background-color: #FFFF09;}\n",
       ".cython.score-250 {background-color: #FFFF09;}\n",
       ".cython.score-251 {background-color: #FFFF09;}\n",
       ".cython.score-252 {background-color: #FFFF09;}\n",
       ".cython.score-253 {background-color: #FFFF09;}\n",
       ".cython.score-254 {background-color: #FFFF09;}\n",
       ".cython .hll { background-color: #ffffcc }\n",
       ".cython  { background: #f8f8f8; }\n",
       ".cython .c { color: #408080; font-style: italic } /* Comment */\n",
       ".cython .err { border: 1px solid #FF0000 } /* Error */\n",
       ".cython .k { color: #008000; font-weight: bold } /* Keyword */\n",
       ".cython .o { color: #666666 } /* Operator */\n",
       ".cython .ch { color: #408080; font-style: italic } /* Comment.Hashbang */\n",
       ".cython .cm { color: #408080; font-style: italic } /* Comment.Multiline */\n",
       ".cython .cp { color: #BC7A00 } /* Comment.Preproc */\n",
       ".cython .cpf { color: #408080; font-style: italic } /* Comment.PreprocFile */\n",
       ".cython .c1 { color: #408080; font-style: italic } /* Comment.Single */\n",
       ".cython .cs { color: #408080; font-style: italic } /* Comment.Special */\n",
       ".cython .gd { color: #A00000 } /* Generic.Deleted */\n",
       ".cython .ge { font-style: italic } /* Generic.Emph */\n",
       ".cython .gr { color: #FF0000 } /* Generic.Error */\n",
       ".cython .gh { color: #000080; font-weight: bold } /* Generic.Heading */\n",
       ".cython .gi { color: #00A000 } /* Generic.Inserted */\n",
       ".cython .go { color: #888888 } /* Generic.Output */\n",
       ".cython .gp { color: #000080; font-weight: bold } /* Generic.Prompt */\n",
       ".cython .gs { font-weight: bold } /* Generic.Strong */\n",
       ".cython .gu { color: #800080; font-weight: bold } /* Generic.Subheading */\n",
       ".cython .gt { color: #0044DD } /* Generic.Traceback */\n",
       ".cython .kc { color: #008000; font-weight: bold } /* Keyword.Constant */\n",
       ".cython .kd { color: #008000; font-weight: bold } /* Keyword.Declaration */\n",
       ".cython .kn { color: #008000; font-weight: bold } /* Keyword.Namespace */\n",
       ".cython .kp { color: #008000 } /* Keyword.Pseudo */\n",
       ".cython .kr { color: #008000; font-weight: bold } /* Keyword.Reserved */\n",
       ".cython .kt { color: #B00040 } /* Keyword.Type */\n",
       ".cython .m { color: #666666 } /* Literal.Number */\n",
       ".cython .s { color: #BA2121 } /* Literal.String */\n",
       ".cython .na { color: #7D9029 } /* Name.Attribute */\n",
       ".cython .nb { color: #008000 } /* Name.Builtin */\n",
       ".cython .nc { color: #0000FF; font-weight: bold } /* Name.Class */\n",
       ".cython .no { color: #880000 } /* Name.Constant */\n",
       ".cython .nd { color: #AA22FF } /* Name.Decorator */\n",
       ".cython .ni { color: #999999; font-weight: bold } /* Name.Entity */\n",
       ".cython .ne { color: #D2413A; font-weight: bold } /* Name.Exception */\n",
       ".cython .nf { color: #0000FF } /* Name.Function */\n",
       ".cython .nl { color: #A0A000 } /* Name.Label */\n",
       ".cython .nn { color: #0000FF; font-weight: bold } /* Name.Namespace */\n",
       ".cython .nt { color: #008000; font-weight: bold } /* Name.Tag */\n",
       ".cython .nv { color: #19177C } /* Name.Variable */\n",
       ".cython .ow { color: #AA22FF; font-weight: bold } /* Operator.Word */\n",
       ".cython .w { color: #bbbbbb } /* Text.Whitespace */\n",
       ".cython .mb { color: #666666 } /* Literal.Number.Bin */\n",
       ".cython .mf { color: #666666 } /* Literal.Number.Float */\n",
       ".cython .mh { color: #666666 } /* Literal.Number.Hex */\n",
       ".cython .mi { color: #666666 } /* Literal.Number.Integer */\n",
       ".cython .mo { color: #666666 } /* Literal.Number.Oct */\n",
       ".cython .sa { color: #BA2121 } /* Literal.String.Affix */\n",
       ".cython .sb { color: #BA2121 } /* Literal.String.Backtick */\n",
       ".cython .sc { color: #BA2121 } /* Literal.String.Char */\n",
       ".cython .dl { color: #BA2121 } /* Literal.String.Delimiter */\n",
       ".cython .sd { color: #BA2121; font-style: italic } /* Literal.String.Doc */\n",
       ".cython .s2 { color: #BA2121 } /* Literal.String.Double */\n",
       ".cython .se { color: #BB6622; font-weight: bold } /* Literal.String.Escape */\n",
       ".cython .sh { color: #BA2121 } /* Literal.String.Heredoc */\n",
       ".cython .si { color: #BB6688; font-weight: bold } /* Literal.String.Interpol */\n",
       ".cython .sx { color: #008000 } /* Literal.String.Other */\n",
       ".cython .sr { color: #BB6688 } /* Literal.String.Regex */\n",
       ".cython .s1 { color: #BA2121 } /* Literal.String.Single */\n",
       ".cython .ss { color: #19177C } /* Literal.String.Symbol */\n",
       ".cython .bp { color: #008000 } /* Name.Builtin.Pseudo */\n",
       ".cython .fm { color: #0000FF } /* Name.Function.Magic */\n",
       ".cython .vc { color: #19177C } /* Name.Variable.Class */\n",
       ".cython .vg { color: #19177C } /* Name.Variable.Global */\n",
       ".cython .vi { color: #19177C } /* Name.Variable.Instance */\n",
       ".cython .vm { color: #19177C } /* Name.Variable.Magic */\n",
       ".cython .il { color: #666666 } /* Literal.Number.Integer.Long */\n",
       "    </style>\n",
       "</head>\n",
       "<body class=\"cython\">\n",
       "<p><span style=\"border-bottom: solid 1px grey;\">Generated by Cython 0.28.3</span></p>\n",
       "<p>\n",
       "    <span style=\"background-color: #FFFF00\">Yellow lines</span> hint at Python interaction.<br />\n",
       "    Click on a line that starts with a \"<code>+</code>\" to see the C code that Cython generated for it.\n",
       "</p>\n",
       "<div class=\"cython\"><pre class=\"cython line score-8\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">01</span>: <span class=\"k\">import</span> <span class=\"nn\">cython</span></pre>\n",
       "<pre class='cython code score-8 '>  __pyx_t_1 = <span class='pyx_c_api'>__Pyx_PyDict_NewPresized</span>(0);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 1, __pyx_L1_error)</span>\n",
       "  <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n",
       "  if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_test, __pyx_t_1) &lt; 0) <span class='error_goto'>__PYX_ERR(0, 1, __pyx_L1_error)</span>\n",
       "  <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n",
       "</pre><pre class=\"cython line score-0\">&#xA0;<span class=\"\">02</span>: <span class=\"k\">from</span> <span class=\"nn\">libc.math</span> <span class=\"k\">cimport</span> <span class=\"n\">sqrt</span></pre>\n",
       "<pre class=\"cython line score-0\">&#xA0;<span class=\"\">03</span>: </pre>\n",
       "<pre class=\"cython line score-0\">&#xA0;<span class=\"\">04</span>: <span class=\"k\">cdef</span> <span class=\"kr\">extern</span> <span class=\"k\">from</span> <span class=\"s\">&quot;complex.h&quot;</span><span class=\"p\">:</span></pre>\n",
       "<pre class=\"cython line score-0\">&#xA0;<span class=\"\">05</span>:     <span class=\"n\">double</span> <span class=\"nb\">complex</span> <span class=\"n\">csqrt</span><span class=\"p\">(</span><span class=\"n\">double</span> <span class=\"nb\">complex</span> <span class=\"n\">z</span><span class=\"p\">)</span></pre>\n",
       "<pre class=\"cython line score-0\">&#xA0;<span class=\"\">06</span>: </pre>\n",
       "<pre class=\"cython line score-0\">&#xA0;<span class=\"\">07</span>: <span class=\"nd\">@cython</span><span class=\"o\">.</span><span class=\"n\">boundscheck</span><span class=\"p\">(</span><span class=\"bp\">False</span><span class=\"p\">)</span></pre>\n",
       "<pre class=\"cython line score-0\">&#xA0;<span class=\"\">08</span>: <span class=\"nd\">@cython</span><span class=\"o\">.</span><span class=\"n\">wraparound</span><span class=\"p\">(</span><span class=\"bp\">False</span><span class=\"p\">)</span></pre>\n",
       "<pre class=\"cython line score-58\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">09</span>: <span class=\"k\">def</span> <span class=\"nf\">test_cython</span><span class=\"p\">(</span><span class=\"nb\">complex</span> <span class=\"p\">[:,:,:,:]</span> <span class=\"n\">a</span><span class=\"p\">,</span> <span class=\"nb\">complex</span> <span class=\"p\">[:,:,:,:]</span> <span class=\"n\">b</span><span class=\"p\">,</span> <span class=\"n\">double</span> <span class=\"p\">[:,:,:,:]</span> <span class=\"n\">l</span><span class=\"p\">,</span> <span class=\"n\">double</span> <span class=\"nb\">complex</span> <span class=\"p\">[:,:,:,:]</span> <span class=\"n\">out</span><span class=\"p\">):</span></pre>\n",
       "<pre class='cython code score-58 '>/* Python wrapper */\n",
       "static PyObject *__pyx_pw_46_cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d_1test_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\n",
       "static PyMethodDef __pyx_mdef_46_cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d_1test_cython = {\"test_cython\", (PyCFunction)__pyx_pw_46_cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d_1test_cython, METH_VARARGS|METH_KEYWORDS, 0};\n",
       "static PyObject *__pyx_pw_46_cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d_1test_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n",
       "  __Pyx_memviewslice __pyx_v_a = { 0, 0, { 0 }, { 0 }, { 0 } };\n",
       "  __Pyx_memviewslice __pyx_v_b = { 0, 0, { 0 }, { 0 }, { 0 } };\n",
       "  __Pyx_memviewslice __pyx_v_l = { 0, 0, { 0 }, { 0 }, { 0 } };\n",
       "  __Pyx_memviewslice __pyx_v_out = { 0, 0, { 0 }, { 0 }, { 0 } };\n",
       "  PyObject *__pyx_r = 0;\n",
       "  <span class='refnanny'>__Pyx_RefNannyDeclarations</span>\n",
       "  <span class='refnanny'>__Pyx_RefNannySetupContext</span>(\"test_cython (wrapper)\", 0);\n",
       "  {\n",
       "    static PyObject **__pyx_pyargnames[] = {&amp;__pyx_n_s_a,&amp;__pyx_n_s_b,&amp;__pyx_n_s_l,&amp;__pyx_n_s_out,0};\n",
       "    PyObject* values[4] = {0,0,0,0};\n",
       "    if (unlikely(__pyx_kwds)) {\n",
       "      Py_ssize_t kw_args;\n",
       "      const Py_ssize_t pos_args = <span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args);\n",
       "      switch (pos_args) {\n",
       "        case  4: values[3] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 3);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  3: values[2] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 2);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2: values[1] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 1);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1: values[0] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 0);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  0: break;\n",
       "        default: goto __pyx_L5_argtuple_error;\n",
       "      }\n",
       "      kw_args = <span class='py_c_api'>PyDict_Size</span>(__pyx_kwds);\n",
       "      switch (pos_args) {\n",
       "        case  0:\n",
       "        if (likely((values[0] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_a)) != 0)) kw_args--;\n",
       "        else goto __pyx_L5_argtuple_error;\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1:\n",
       "        if (likely((values[1] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_b)) != 0)) kw_args--;\n",
       "        else {\n",
       "          <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>(\"test_cython\", 1, 4, 4, 1); <span class='error_goto'>__PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "        }\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2:\n",
       "        if (likely((values[2] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_l)) != 0)) kw_args--;\n",
       "        else {\n",
       "          <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>(\"test_cython\", 1, 4, 4, 2); <span class='error_goto'>__PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "        }\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  3:\n",
       "        if (likely((values[3] = <span class='pyx_c_api'>__Pyx_PyDict_GetItemStr</span>(__pyx_kwds, __pyx_n_s_out)) != 0)) kw_args--;\n",
       "        else {\n",
       "          <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>(\"test_cython\", 1, 4, 4, 3); <span class='error_goto'>__PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "        }\n",
       "      }\n",
       "      if (unlikely(kw_args &gt; 0)) {\n",
       "        if (unlikely(<span class='pyx_c_api'>__Pyx_ParseOptionalKeywords</span>(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"test_cython\") &lt; 0)) <span class='error_goto'>__PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "      }\n",
       "    } else if (<span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args) != 4) {\n",
       "      goto __pyx_L5_argtuple_error;\n",
       "    } else {\n",
       "      values[0] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 0);\n",
       "      values[1] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 1);\n",
       "      values[2] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 2);\n",
       "      values[3] = <span class='py_macro_api'>PyTuple_GET_ITEM</span>(__pyx_args, 3);\n",
       "    }\n",
       "    __pyx_v_a = <span class='pyx_c_api'>__Pyx_PyObject_to_MemoryviewSlice_dsdsdsds___pyx_t_double_complex</span>(values[0], PyBUF_WRITABLE);<span class='error_goto'> if (unlikely(!__pyx_v_a.memview)) __PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "    __pyx_v_b = <span class='pyx_c_api'>__Pyx_PyObject_to_MemoryviewSlice_dsdsdsds___pyx_t_double_complex</span>(values[1], PyBUF_WRITABLE);<span class='error_goto'> if (unlikely(!__pyx_v_b.memview)) __PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "    __pyx_v_l = <span class='pyx_c_api'>__Pyx_PyObject_to_MemoryviewSlice_dsdsdsds_double</span>(values[2], PyBUF_WRITABLE);<span class='error_goto'> if (unlikely(!__pyx_v_l.memview)) __PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "    __pyx_v_out = <span class='pyx_c_api'>__Pyx_PyObject_to_MemoryviewSlice_dsdsdsds___pyx_t_double_complex</span>(values[3], PyBUF_WRITABLE);<span class='error_goto'> if (unlikely(!__pyx_v_out.memview)) __PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "  }\n",
       "  goto __pyx_L4_argument_unpacking_done;\n",
       "  __pyx_L5_argtuple_error:;\n",
       "  <span class='pyx_c_api'>__Pyx_RaiseArgtupleInvalid</span>(\"test_cython\", 1, 4, 4, <span class='py_macro_api'>PyTuple_GET_SIZE</span>(__pyx_args)); <span class='error_goto'>__PYX_ERR(0, 9, __pyx_L3_error)</span>\n",
       "  __pyx_L3_error:;\n",
       "  <span class='pyx_c_api'>__Pyx_AddTraceback</span>(\"_cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d.test_cython\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n",
       "  return NULL;\n",
       "  __pyx_L4_argument_unpacking_done:;\n",
       "  __pyx_r = __pyx_pf_46_cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d_test_cython(__pyx_self, __pyx_v_a, __pyx_v_b, __pyx_v_l, __pyx_v_out);\n",
       "\n",
       "  /* function exit code */\n",
       "  <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "static PyObject *__pyx_pf_46_cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d_test_cython(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_a, __Pyx_memviewslice __pyx_v_b, __Pyx_memviewslice __pyx_v_l, __Pyx_memviewslice __pyx_v_out) {\n",
       "  size_t __pyx_v_nf;\n",
       "  size_t __pyx_v_no;\n",
       "  size_t __pyx_v_nl;\n",
       "  size_t __pyx_v_nw;\n",
       "  PyObject *__pyx_r = NULL;\n",
       "  <span class='refnanny'>__Pyx_RefNannyDeclarations</span>\n",
       "  <span class='refnanny'>__Pyx_RefNannySetupContext</span>(\"test_cython\", 0);\n",
       "/* … */\n",
       "  /* function exit code */\n",
       "  __pyx_r = Py_None; <span class='pyx_macro_api'>__Pyx_INCREF</span>(Py_None);\n",
       "  __PYX_XDEC_MEMVIEW(&amp;__pyx_v_a, 1);\n",
       "  __PYX_XDEC_MEMVIEW(&amp;__pyx_v_b, 1);\n",
       "  __PYX_XDEC_MEMVIEW(&amp;__pyx_v_l, 1);\n",
       "  __PYX_XDEC_MEMVIEW(&amp;__pyx_v_out, 1);\n",
       "  <span class='refnanny'>__Pyx_XGIVEREF</span>(__pyx_r);\n",
       "  <span class='refnanny'>__Pyx_RefNannyFinishContext</span>();\n",
       "  return __pyx_r;\n",
       "}\n",
       "/* … */\n",
       "  __pyx_tuple__22 = <span class='py_c_api'>PyTuple_Pack</span>(8, __pyx_n_s_a, __pyx_n_s_b, __pyx_n_s_l, __pyx_n_s_out, __pyx_n_s_nf, __pyx_n_s_no, __pyx_n_s_nl, __pyx_n_s_nw);<span class='error_goto'> if (unlikely(!__pyx_tuple__22)) __PYX_ERR(0, 9, __pyx_L1_error)</span>\n",
       "  <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_tuple__22);\n",
       "  <span class='refnanny'>__Pyx_GIVEREF</span>(__pyx_tuple__22);\n",
       "/* … */\n",
       "  __pyx_t_1 = PyCFunction_NewEx(&amp;__pyx_mdef_46_cython_magic_8b5ebf652ec46e4303ed6c9de0b6e87d_1test_cython, NULL, __pyx_n_s_cython_magic_8b5ebf652ec46e4303);<span class='error_goto'> if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 9, __pyx_L1_error)</span>\n",
       "  <span class='refnanny'>__Pyx_GOTREF</span>(__pyx_t_1);\n",
       "  if (<span class='py_c_api'>PyDict_SetItem</span>(__pyx_d, __pyx_n_s_test_cython, __pyx_t_1) &lt; 0) <span class='error_goto'>__PYX_ERR(0, 9, __pyx_L1_error)</span>\n",
       "  <span class='pyx_macro_api'>__Pyx_DECREF</span>(__pyx_t_1); __pyx_t_1 = 0;\n",
       "  __pyx_codeobj__23 = (PyObject*)<span class='pyx_c_api'>__Pyx_PyCode_New</span>(4, 0, 8, 0, CO_OPTIMIZED|CO_NEWLOCALS, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__22, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_dtr_cache_ipython_cython, __pyx_n_s_test_cython, 9, __pyx_empty_bytes);<span class='error_goto'> if (unlikely(!__pyx_codeobj__23)) __PYX_ERR(0, 9, __pyx_L1_error)</span>\n",
       "</pre><pre class=\"cython line score-0\">&#xA0;<span class=\"\">10</span>:     <span class=\"k\">cdef</span> <span class=\"kt\">size_t</span> <span class=\"nf\">nf</span><span class=\"p\">,</span> <span class=\"nf\">no</span><span class=\"p\">,</span> <span class=\"nf\">nl</span><span class=\"p\">,</span> <span class=\"nf\">nw</span></pre>\n",
       "<pre class=\"cython line score-0\">&#xA0;<span class=\"\">11</span>: </pre>\n",
       "<pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">12</span>:     <span class=\"k\">for</span> <span class=\"n\">nf</span> <span class=\"ow\">in</span> <span class=\"nb\">xrange</span><span class=\"p\">(</span><span class=\"n\">out</span><span class=\"o\">.</span><span class=\"n\">shape</span><span class=\"p\">[</span><span class=\"mf\">0</span><span class=\"p\">]):</span></pre>\n",
       "<pre class='cython code score-0 '>  __pyx_t_1 = (__pyx_v_out.shape[0]);\n",
       "  __pyx_t_2 = __pyx_t_1;\n",
       "  for (__pyx_t_3 = 0; __pyx_t_3 &lt; __pyx_t_2; __pyx_t_3+=1) {\n",
       "    __pyx_v_nf = __pyx_t_3;\n",
       "</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">13</span>:         <span class=\"k\">for</span> <span class=\"n\">no</span> <span class=\"ow\">in</span> <span class=\"nb\">xrange</span><span class=\"p\">(</span><span class=\"n\">out</span><span class=\"o\">.</span><span class=\"n\">shape</span><span class=\"p\">[</span><span class=\"mf\">1</span><span class=\"p\">]):</span></pre>\n",
       "<pre class='cython code score-0 '>    __pyx_t_4 = (__pyx_v_out.shape[1]);\n",
       "    __pyx_t_5 = __pyx_t_4;\n",
       "    for (__pyx_t_6 = 0; __pyx_t_6 &lt; __pyx_t_5; __pyx_t_6+=1) {\n",
       "      __pyx_v_no = __pyx_t_6;\n",
       "</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">14</span>:             <span class=\"k\">for</span> <span class=\"n\">nl</span> <span class=\"ow\">in</span> <span class=\"nb\">xrange</span><span class=\"p\">(</span><span class=\"n\">out</span><span class=\"o\">.</span><span class=\"n\">shape</span><span class=\"p\">[</span><span class=\"mf\">2</span><span class=\"p\">]):</span></pre>\n",
       "<pre class='cython code score-0 '>      __pyx_t_7 = (__pyx_v_out.shape[2]);\n",
       "      __pyx_t_8 = __pyx_t_7;\n",
       "      for (__pyx_t_9 = 0; __pyx_t_9 &lt; __pyx_t_8; __pyx_t_9+=1) {\n",
       "        __pyx_v_nl = __pyx_t_9;\n",
       "</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">15</span>:                 <span class=\"k\">for</span> <span class=\"n\">nw</span> <span class=\"ow\">in</span> <span class=\"nb\">xrange</span><span class=\"p\">(</span><span class=\"n\">out</span><span class=\"o\">.</span><span class=\"n\">shape</span><span class=\"p\">[</span><span class=\"mf\">3</span><span class=\"p\">]):</span></pre>\n",
       "<pre class='cython code score-0 '>        __pyx_t_10 = (__pyx_v_out.shape[3]);\n",
       "        __pyx_t_11 = __pyx_t_10;\n",
       "        for (__pyx_t_12 = 0; __pyx_t_12 &lt; __pyx_t_11; __pyx_t_12+=1) {\n",
       "          __pyx_v_nw = __pyx_t_12;\n",
       "</pre><pre class=\"cython line score-0\" onclick=\"(function(s){s.display=s.display==='block'?'none':'block'})(this.nextElementSibling.style)\">+<span class=\"\">16</span>:                     <span class=\"n\">out</span><span class=\"p\">[</span><span class=\"n\">nf</span><span class=\"p\">,</span> <span class=\"n\">no</span><span class=\"p\">,</span> <span class=\"n\">nl</span><span class=\"p\">,</span> <span class=\"n\">nw</span><span class=\"p\">]</span> <span class=\"o\">=</span> <span class=\"n\">csqrt</span><span class=\"p\">(</span><span class=\"n\">a</span><span class=\"p\">[</span><span class=\"n\">nf</span><span class=\"p\">,</span> <span class=\"mf\">0</span><span class=\"p\">,</span> <span class=\"n\">nl</span><span class=\"p\">,</span> <span class=\"mf\">0</span><span class=\"p\">]</span> <span class=\"o\">*</span> <span class=\"n\">l</span><span class=\"p\">[</span><span class=\"mf\">0</span><span class=\"p\">,</span> <span class=\"n\">no</span><span class=\"p\">,</span> <span class=\"mf\">0</span><span class=\"p\">,</span> <span class=\"n\">nw</span><span class=\"p\">]</span> <span class=\"o\">**</span> <span class=\"mf\">2</span> <span class=\"o\">+</span> <span class=\"n\">b</span><span class=\"p\">[</span><span class=\"n\">nf</span><span class=\"p\">,</span> <span class=\"mf\">0</span><span class=\"p\">,</span> <span class=\"n\">nl</span><span class=\"p\">,</span> <span class=\"mf\">0</span><span class=\"p\">])</span></pre>\n",
       "<pre class='cython code score-0 '>          __pyx_t_13 = __pyx_v_nf;\n",
       "          __pyx_t_14 = 0;\n",
       "          __pyx_t_15 = __pyx_v_nl;\n",
       "          __pyx_t_16 = 0;\n",
       "          __pyx_t_17 = 0;\n",
       "          __pyx_t_18 = __pyx_v_no;\n",
       "          __pyx_t_19 = 0;\n",
       "          __pyx_t_20 = __pyx_v_nw;\n",
       "          __pyx_t_21 = __pyx_v_nf;\n",
       "          __pyx_t_22 = 0;\n",
       "          __pyx_t_23 = __pyx_v_nl;\n",
       "          __pyx_t_24 = 0;\n",
       "          __pyx_t_25 = __pyx_v_nf;\n",
       "          __pyx_t_26 = __pyx_v_no;\n",
       "          __pyx_t_27 = __pyx_v_nl;\n",
       "          __pyx_t_28 = __pyx_v_nw;\n",
       "          *((__pyx_t_double_complex *) ( /* dim=3 */ (( /* dim=2 */ (( /* dim=1 */ (( /* dim=0 */ (__pyx_v_out.data + __pyx_t_25 * __pyx_v_out.strides[0]) ) + __pyx_t_26 * __pyx_v_out.strides[1]) ) + __pyx_t_27 * __pyx_v_out.strides[2]) ) + __pyx_t_28 * __pyx_v_out.strides[3]) )) = csqrt(__Pyx_c_sum_double(__Pyx_c_prod_double((*((__pyx_t_double_complex *) ( /* dim=3 */ (( /* dim=2 */ (( /* dim=1 */ (( /* dim=0 */ (__pyx_v_a.data + __pyx_t_13 * __pyx_v_a.strides[0]) ) + __pyx_t_14 * __pyx_v_a.strides[1]) ) + __pyx_t_15 * __pyx_v_a.strides[2]) ) + __pyx_t_16 * __pyx_v_a.strides[3]) ))), __pyx_t_double_complex_from_parts(pow((*((double *) ( /* dim=3 */ (( /* dim=2 */ (( /* dim=1 */ (( /* dim=0 */ (__pyx_v_l.data + __pyx_t_17 * __pyx_v_l.strides[0]) ) + __pyx_t_18 * __pyx_v_l.strides[1]) ) + __pyx_t_19 * __pyx_v_l.strides[2]) ) + __pyx_t_20 * __pyx_v_l.strides[3]) ))), 2.0), 0)), (*((__pyx_t_double_complex *) ( /* dim=3 */ (( /* dim=2 */ (( /* dim=1 */ (( /* dim=0 */ (__pyx_v_b.data + __pyx_t_21 * __pyx_v_b.strides[0]) ) + __pyx_t_22 * __pyx_v_b.strides[1]) ) + __pyx_t_23 * __pyx_v_b.strides[2]) ) + __pyx_t_24 * __pyx_v_b.strides[3]) )))));\n",
       "        }\n",
       "      }\n",
       "    }\n",
       "  }\n",
       "</pre></div></body></html>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%%cython -a\n",
    "import cython\n",
    "from libc.math cimport sqrt\n",
    "\n",
    "cdef extern from \"complex.h\":\n",
    "    double complex csqrt(double complex z)\n",
    "\n",
    "@cython.boundscheck(False)\n",
    "@cython.wraparound(False)\n",
    "def test_cython(complex [:,:,:,:] a, complex [:,:,:,:] b, double [:,:,:,:] l, double complex [:,:,:,:] out):\n",
    "    cdef size_t nf, no, nl, nw\n",
    "    \n",
    "    for nf in xrange(out.shape[0]):\n",
    "        for no in xrange(out.shape[1]):\n",
    "            for nl in xrange(out.shape[2]):\n",
    "                for nw in xrange(out.shape[3]):\n",
    "                    out[nf, no, nl, nw] = csqrt(a[nf, 0, nl, 0] * l[0, no, 0, nw] ** 2 + b[nf, 0, nl, 0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run comparison for a small and a big matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " Shape Test Matrix :: (1, 1, 1, 1) ; total # elements:: 1\n",
      "------------------------------------------------------------------------------------------\n",
      " NumPy            ::  5.57 µs ± 447 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      " Numexpr          ::  18.6 µs ± 684 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)\n",
      " Numba @vectorize ::  3.05 µs ± 371 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      " Numba @njit      ::  3.34 µs ± 157 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      " Cython           ::  5.16 µs ± 765 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n",
      "\n",
      " Shape Test Matrix :: (1, 10, 20, 50) ; total # elements:: 10000\n",
      "------------------------------------------------------------------------------------------\n",
      " NumPy            ::  599 µs ± 81.9 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      " Numexpr          ::  837 µs ± 47.1 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      " Numba @vectorize ::  350 µs ± 36.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      " Numba @njit      ::  285 µs ± 3.44 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      " Cython           ::  306 µs ± 4.22 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)\n",
      "\n",
      " Shape Test Matrix :: (10, 50, 200, 500) ; total # elements:: 50000000\n",
      "------------------------------------------------------------------------------------------\n",
      " NumPy            ::  2.52 s ± 9.85 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      " Numexpr          ::  4.12 s ± 72.8 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      " Numba @vectorize ::  1.9 s ± 36.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      " Numba @njit      ::  1.78 s ± 21.6 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      " Cython           ::  1.92 s ± 44.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n",
      "\n"
     ]
    }
   ],
   "source": [
    "fre = [1, 1, 10]\n",
    "off = [1, 10, 50]\n",
    "lay = [1, 20, 200]\n",
    "wav = [1, 50, 500] \n",
    "\n",
    "for i in range(3):\n",
    "    nlay, nfre, noff, nwav = lay[i], fre[i], off[i], wav[i]\n",
    "    a = np.random.rand(nfre, nlay) + 1j*np.random.rand(nfre, nlay)\n",
    "    b = np.random.rand(nfre, nlay) + 1j*np.random.rand(nfre, nlay)\n",
    "    l = np.random.rand(noff, nwav)\n",
    "    \n",
    "    # Broadcast shapes\n",
    "    a_bc = a[:, None, :, None]\n",
    "    b_bc = b[:, None, :, None]\n",
    "    l_bc = l[None, :, None, :]\n",
    "    \n",
    "    # Output shape\n",
    "    out_shape = (nfre, noff, nlay, nwav)\n",
    "    \n",
    "    print(' Shape Test Matrix ::', out_shape, '; total # elements:: '+str(nfre*noff*nlay*nwav))\n",
    "    print('------------------------------------------------------------------------------------------')\n",
    "    \n",
    "    print(' NumPy            ::  ', end='')\n",
    "    # Get NumPy result for comparison\n",
    "    numpy_result = test_numpy(a_bc, b_bc, l_bc)\n",
    "    # Get runtime\n",
    "    %timeit test_numpy(a_bc, b_bc, l_bc)\n",
    "    \n",
    "    print(' Numexpr          ::  ', end='')\n",
    "    # Ensure it agrees with NumPy\n",
    "    numexpr_result = test_numexpr(a_bc, b_bc, l_bc)\n",
    "    if not np.allclose(numpy_result, numexpr_result, atol=0, rtol=1e-10):\n",
    "        print(' * FAIL, DOES NOT AGREE WITH NumPy RESULT!')\n",
    "    # Get runtime\n",
    "    %timeit test_numexpr(a_bc, b_bc, l_bc)\n",
    "    \n",
    "    print(' Numba @vectorize ::  ', end='')\n",
    "    # Ensure it agrees with NumPy\n",
    "    numba_v_result = test_numba_v(a_bc, b_bc, l_bc)\n",
    "    if not np.allclose(numpy_result, numba_v_result, atol=0, rtol=1e-10):\n",
    "        print(' * FAIL, DOES NOT AGREE WITH NumPy RESULT!')\n",
    "    # Get runtime\n",
    "    %timeit test_numba_v(a_bc, b_bc, l_bc)\n",
    "    \n",
    "    print(' Numba @njit      ::  ', end='')\n",
    "    # Ensure it agrees with NumPy\n",
    "    numba_n_result = np.empty(out_shape, dtype=complex)\n",
    "    test_numba_n(a_bc, b_bc, l_bc, numba_n_result)\n",
    "    if not np.allclose(numpy_result, numba_n_result, atol=0, rtol=1e-10):\n",
    "        print(' * FAIL, DOES NOT AGREE WITH NumPy RESULT!')\n",
    "    # Get runtime\n",
    "    %timeit test_numba_n(a_bc, b_bc, l_bc, np.empty(out_shape, dtype=complex))\n",
    "    \n",
    "    print(' Cython           ::  ', end='')\n",
    "    # Ensure it agrees with NumPy\n",
    "    cython_result = np.empty(out_shape, dtype=complex)\n",
    "    test_cython(a_bc, b_bc, l_bc, cython_result)\n",
    "    if not np.allclose(numpy_result, cython_result, atol=0, rtol=1e-10):\n",
    "        print(' * FAIL, DOES NOT AGREE WITH NumPy RESULT!')\n",
    "    # Get runtime\n",
    "    %timeit test_cython(a_bc, b_bc, l_bc, np.empty(out_shape, dtype=complex))\n",
    "        \n",
    "    print()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "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='10'>Fri Jun 29 11:23:12 2018 CDT</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.13.3</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.1.0</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>scipy</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>1.7.1</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>empymod</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>6.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;'>2.6.5</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>numexpr</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>2.2.2</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>matplotlib</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>0.28.3</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>cython</td>\n",
       "    <td style='text-align: right; background-color: #ccc; border: 2px solid #fff;'>0.38.1+1.gc42707d0f.dirty</td>\n",
       "    <td style='text-align: left; border: 2px solid #fff;'>numba</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "     <td style='text-align: center; border: 2px solid #fff;' colspan='10'>3.6.5 |Anaconda custom (64-bit)| (default, Apr 29 2018, 16:14:56) \n",
       "[GCC 7.2.0]</td>\n",
       "  </tr>\n",
       "  <tr>\n",
       "     <td style='text-align: center; background-color: #ddd;border: 2px solid #fff;' colspan='10'>Intel(R) Math Kernel Library Version 2018.0.3 Product Build 20180406 for Intel(R) 64 architecture applications</td>\n",
       "  </tr>\n",
       "</table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from empymod import versions\n",
    "versions('HTML', add_pckg=[cython, numba], ncol=5)"
   ]
  }
 ],
 "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.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}