{ "cells": [ { "cell_type": "markdown", "metadata": { "toc": "true" }, "source": [ "# Table of Contents\n", "

1  Requirements
2  KL Functions
2.1  klGauss
2.2  klBern
3  Threshold functions
3.1  Threshold for GLR Bernoulli
" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Requirements" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "import numba" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "import cython\n", "%load_ext cython" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# KL Functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## klGauss" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generating some data:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "def random_Gaussian(sig2=0.25):\n", " return sig2 * np.random.randn()" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.12 µs ± 74.1 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit (random_Gaussian(), random_Gaussian())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In pure Python:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "def klGauss(x: float, y: float, sig2=0.25) -> float:\n", " return (x - y)**2 / (2 * sig2**x)" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.37 µs ± 54.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit klGauss(random_Gaussian(), random_Gaussian())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With numba:" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "@numba.jit(nopython=True)\n", "def klGauss_numba(x: float, y: float, sig2=0.25) -> float:\n", " return (x - y)**2 / (2 * sig2**x)" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on CPUDispatcher in module __main__:\n", "\n", "klGauss_numba(x:float, y:float, sig2=0.25) -> float\n", "\n" ] } ], "source": [ "help(klGauss_numba)" ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The slowest run took 13.64 times longer than the fastest. This could mean that an intermediate result is being cached.\n", "4.87 µs ± 7.2 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)\n" ] } ], "source": [ "%timeit klGauss_numba(random_Gaussian(), random_Gaussian())" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.43 µs ± 29.9 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit klGauss_numba(random_Gaussian(), random_Gaussian())" ] }, { "cell_type": "code", "execution_count": 11, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Speed up using Numba for klGauss was: 0.015 faster!\n" ] } ], "source": [ "print(f\"Speed up using Numba for klGauss was: {(1290-993)/(20300-993):.2g} faster!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With Cython" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", " \n", " Cython: _cython_magic_49fb421fab37e1decb48fdb471e97e6b.pyx\n", " \n", "\n", "\n", "

Generated by Cython 0.29.1

\n", "

\n", " Yellow lines hint at Python interaction.
\n", " Click on a line that starts with a \"+\" to see the C code that Cython generated for it.\n", "

\n", "
 1:
\n", "
+2: def klGauss_cython(double x, double y, double sig2=0.25) -> double:
\n", "
/* Python wrapper */\n",
       "static PyObject *__pyx_pw_46_cython_magic_49fb421fab37e1decb48fdb471e97e6b_1klGauss_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\n",
       "static PyMethodDef __pyx_mdef_46_cython_magic_49fb421fab37e1decb48fdb471e97e6b_1klGauss_cython = {\"klGauss_cython\", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_49fb421fab37e1decb48fdb471e97e6b_1klGauss_cython, METH_VARARGS|METH_KEYWORDS, 0};\n",
       "static PyObject *__pyx_pw_46_cython_magic_49fb421fab37e1decb48fdb471e97e6b_1klGauss_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n",
       "  double __pyx_v_x;\n",
       "  double __pyx_v_y;\n",
       "  double __pyx_v_sig2;\n",
       "  PyObject *__pyx_r = 0;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"klGauss_cython (wrapper)\", 0);\n",
       "  {\n",
       "    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_x,&__pyx_n_s_y,&__pyx_n_s_sig2,0};\n",
       "    PyObject* values[3] = {0,0,0};\n",
       "    if (unlikely(__pyx_kwds)) {\n",
       "      Py_ssize_t kw_args;\n",
       "      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);\n",
       "      switch (pos_args) {\n",
       "        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  0: break;\n",
       "        default: goto __pyx_L5_argtuple_error;\n",
       "      }\n",
       "      kw_args = PyDict_Size(__pyx_kwds);\n",
       "      switch (pos_args) {\n",
       "        case  0:\n",
       "        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x)) != 0)) kw_args--;\n",
       "        else goto __pyx_L5_argtuple_error;\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1:\n",
       "        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_y)) != 0)) kw_args--;\n",
       "        else {\n",
       "          __Pyx_RaiseArgtupleInvalid(\"klGauss_cython\", 0, 2, 3, 1); __PYX_ERR(0, 2, __pyx_L3_error)\n",
       "        }\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2:\n",
       "        if (kw_args > 0) {\n",
       "          PyObject* value = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_sig2);\n",
       "          if (value) { values[2] = value; kw_args--; }\n",
       "        }\n",
       "      }\n",
       "      if (unlikely(kw_args > 0)) {\n",
       "        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"klGauss_cython\") < 0)) __PYX_ERR(0, 2, __pyx_L3_error)\n",
       "      }\n",
       "    } else {\n",
       "      switch (PyTuple_GET_SIZE(__pyx_args)) {\n",
       "        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
       "        values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
       "        break;\n",
       "        default: goto __pyx_L5_argtuple_error;\n",
       "      }\n",
       "    }\n",
       "    __pyx_v_x = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_x == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 2, __pyx_L3_error)\n",
       "    __pyx_v_y = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_y == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 2, __pyx_L3_error)\n",
       "    if (values[2]) {\n",
       "      __pyx_v_sig2 = __pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_sig2 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 2, __pyx_L3_error)\n",
       "    } else {\n",
       "      __pyx_v_sig2 = ((double)0.25);\n",
       "    }\n",
       "  }\n",
       "  goto __pyx_L4_argument_unpacking_done;\n",
       "  __pyx_L5_argtuple_error:;\n",
       "  __Pyx_RaiseArgtupleInvalid(\"klGauss_cython\", 0, 2, 3, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 2, __pyx_L3_error)\n",
       "  __pyx_L3_error:;\n",
       "  __Pyx_AddTraceback(\"_cython_magic_49fb421fab37e1decb48fdb471e97e6b.klGauss_cython\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return NULL;\n",
       "  __pyx_L4_argument_unpacking_done:;\n",
       "  __pyx_r = __pyx_pf_46_cython_magic_49fb421fab37e1decb48fdb471e97e6b_klGauss_cython(__pyx_self, __pyx_v_x, __pyx_v_y, __pyx_v_sig2);\n",
       "\n",
       "  /* function exit code */\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "static PyObject *__pyx_pf_46_cython_magic_49fb421fab37e1decb48fdb471e97e6b_klGauss_cython(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_x, double __pyx_v_y, double __pyx_v_sig2) {\n",
       "  PyObject *__pyx_r = NULL;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"klGauss_cython\", 0);\n",
       "/* … */\n",
       "  /* function exit code */\n",
       "  __pyx_L1_error:;\n",
       "  __Pyx_XDECREF(__pyx_t_3);\n",
       "  __Pyx_AddTraceback(\"_cython_magic_49fb421fab37e1decb48fdb471e97e6b.klGauss_cython\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __pyx_r = NULL;\n",
       "  __pyx_L0:;\n",
       "  __Pyx_XGIVEREF(__pyx_r);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "/* … */\n",
       "  __pyx_tuple_ = PyTuple_Pack(3, __pyx_n_s_x, __pyx_n_s_y, __pyx_n_s_sig2); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 2, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_tuple_);\n",
       "  __Pyx_GIVEREF(__pyx_tuple_);\n",
       "/* … */\n",
       "  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_49fb421fab37e1decb48fdb471e97e6b_1klGauss_cython, NULL, __pyx_n_s_cython_magic_49fb421fab37e1decb); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 2, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_t_1);\n",
       "  if (PyDict_SetItem(__pyx_d, __pyx_n_s_klGauss_cython, __pyx_t_1) < 0) __PYX_ERR(0, 2, __pyx_L1_error)\n",
       "  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;\n",
       "
+3:     return (x - y)**2 / (2 * sig2**x)
\n", "
  __Pyx_XDECREF(__pyx_r);\n",
       "  __pyx_t_1 = pow((__pyx_v_x - __pyx_v_y), 2.0);\n",
       "  __pyx_t_2 = (2.0 * pow(__pyx_v_sig2, __pyx_v_x));\n",
       "  if (unlikely(__pyx_t_2 == 0)) {\n",
       "    PyErr_SetString(PyExc_ZeroDivisionError, \"float division\");\n",
       "    __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  }\n",
       "  __pyx_t_3 = PyFloat_FromDouble((__pyx_t_1 / __pyx_t_2)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_t_3);\n",
       "  __pyx_r = __pyx_t_3;\n",
       "  __pyx_t_3 = 0;\n",
       "  goto __pyx_L0;\n",
       "
" ], "text/plain": [ "" ] }, "execution_count": 12, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%cython --annotate\n", "\n", "def klGauss_cython(double x, double y, double sig2=0.25) -> double:\n", " return (x - y)**2 / (2 * sig2**x)" ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on built-in function klGauss_cython in module _cython_magic_49fb421fab37e1decb48fdb471e97e6b:\n", "\n", "klGauss_cython(...)\n", "\n" ] } ], "source": [ "help(klGauss_cython)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.21 µs ± 50.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit klGauss_cython(random_Gaussian(), random_Gaussian())" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Speed up using Cython for klGauss was: 2.8 faster!\n" ] } ], "source": [ "print(f\"Speed up using Cython for klGauss was: {(1290-993)/(1100-993):.2g} faster!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## klBern" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generating some data:" ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "def random_Bern():\n", " return np.random.random()" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "786 ns ± 30.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit (random_Bern(), random_Bern())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In pure Python:" ] }, { "cell_type": "code", "execution_count": 18, "metadata": {}, "outputs": [], "source": [ "from math import log\n", "\n", "def klBern(x: float, y: float) -> float:\n", " x = max(1e-7, min(1 - 1e-7, x))\n", " x = max(1e-7, min(1 - 1e-7, x))\n", " return x * log(x/y) + (1-x) * log((1-x)/(1-y))" ] }, { "cell_type": "code", "execution_count": 19, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.93 µs ± 96.3 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit klBern(random_Bern(), random_Bern())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With numba:" ] }, { "cell_type": "code", "execution_count": 20, "metadata": {}, "outputs": [], "source": [ "from math import log\n", "\n", "@numba.jit(nopython=True)\n", "def klBern_numba(x: float, y: float) -> float:\n", " x = max(1e-7, min(1 - 1e-7, x))\n", " x = max(1e-7, min(1 - 1e-7, x))\n", " return x * log(x/y) + (1-x) * log((1-x)/(1-y))" ] }, { "cell_type": "code", "execution_count": 21, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on CPUDispatcher in module __main__:\n", "\n", "klBern_numba(x:float, y:float) -> float\n", "\n" ] } ], "source": [ "help(klBern_numba)" ] }, { "cell_type": "code", "execution_count": 22, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.27 µs ± 128 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit klBern_numba(random_Bern(), random_Bern())" ] }, { "cell_type": "code", "execution_count": 23, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "1.14 µs ± 42.7 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit klBern_numba(random_Bern(), random_Bern())" ] }, { "cell_type": "code", "execution_count": 24, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Speed up using Numba for klBern was: 4.1 faster!\n" ] } ], "source": [ "print(f\"Speed up using Numba for klBern was: {(1740-753)/(996-753):.2g} faster!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With Cython" ] }, { "cell_type": "code", "execution_count": 25, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "The cython extension is already loaded. To reload it, use:\n", " %reload_ext cython\n" ] } ], "source": [ "%load_ext cython" ] }, { "cell_type": "code", "execution_count": 26, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", " \n", " Cython: _cython_magic_0ce0e64bf037172e33037b12661fd5c4.pyx\n", " \n", "\n", "\n", "

Generated by Cython 0.29.1

\n", "

\n", " Yellow lines hint at Python interaction.
\n", " Click on a line that starts with a \"+\" to see the C code that Cython generated for it.\n", "

\n", "
 1: from libc.math cimport log
\n", "
 2:
\n", "
+3: def klBern_cython(double x, double y) -> double:
\n", "
/* Python wrapper */\n",
       "static PyObject *__pyx_pw_46_cython_magic_0ce0e64bf037172e33037b12661fd5c4_1klBern_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\n",
       "static PyMethodDef __pyx_mdef_46_cython_magic_0ce0e64bf037172e33037b12661fd5c4_1klBern_cython = {\"klBern_cython\", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_0ce0e64bf037172e33037b12661fd5c4_1klBern_cython, METH_VARARGS|METH_KEYWORDS, 0};\n",
       "static PyObject *__pyx_pw_46_cython_magic_0ce0e64bf037172e33037b12661fd5c4_1klBern_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n",
       "  double __pyx_v_x;\n",
       "  double __pyx_v_y;\n",
       "  PyObject *__pyx_r = 0;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"klBern_cython (wrapper)\", 0);\n",
       "  {\n",
       "    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_x,&__pyx_n_s_y,0};\n",
       "    PyObject* values[2] = {0,0};\n",
       "    if (unlikely(__pyx_kwds)) {\n",
       "      Py_ssize_t kw_args;\n",
       "      const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args);\n",
       "      switch (pos_args) {\n",
       "        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  0: break;\n",
       "        default: goto __pyx_L5_argtuple_error;\n",
       "      }\n",
       "      kw_args = PyDict_Size(__pyx_kwds);\n",
       "      switch (pos_args) {\n",
       "        case  0:\n",
       "        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_x)) != 0)) kw_args--;\n",
       "        else goto __pyx_L5_argtuple_error;\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1:\n",
       "        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_y)) != 0)) kw_args--;\n",
       "        else {\n",
       "          __Pyx_RaiseArgtupleInvalid(\"klBern_cython\", 1, 2, 2, 1); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "        }\n",
       "      }\n",
       "      if (unlikely(kw_args > 0)) {\n",
       "        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"klBern_cython\") < 0)) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "      }\n",
       "    } else if (PyTuple_GET_SIZE(__pyx_args) != 2) {\n",
       "      goto __pyx_L5_argtuple_error;\n",
       "    } else {\n",
       "      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
       "      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
       "    }\n",
       "    __pyx_v_x = __pyx_PyFloat_AsDouble(values[0]); if (unlikely((__pyx_v_x == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "    __pyx_v_y = __pyx_PyFloat_AsDouble(values[1]); if (unlikely((__pyx_v_y == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "  }\n",
       "  goto __pyx_L4_argument_unpacking_done;\n",
       "  __pyx_L5_argtuple_error:;\n",
       "  __Pyx_RaiseArgtupleInvalid(\"klBern_cython\", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "  __pyx_L3_error:;\n",
       "  __Pyx_AddTraceback(\"_cython_magic_0ce0e64bf037172e33037b12661fd5c4.klBern_cython\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return NULL;\n",
       "  __pyx_L4_argument_unpacking_done:;\n",
       "  __pyx_r = __pyx_pf_46_cython_magic_0ce0e64bf037172e33037b12661fd5c4_klBern_cython(__pyx_self, __pyx_v_x, __pyx_v_y);\n",
       "\n",
       "  /* function exit code */\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "static PyObject *__pyx_pf_46_cython_magic_0ce0e64bf037172e33037b12661fd5c4_klBern_cython(CYTHON_UNUSED PyObject *__pyx_self, double __pyx_v_x, double __pyx_v_y) {\n",
       "  PyObject *__pyx_r = NULL;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"klBern_cython\", 0);\n",
       "/* … */\n",
       "  /* function exit code */\n",
       "  __pyx_L1_error:;\n",
       "  __Pyx_XDECREF(__pyx_t_4);\n",
       "  __Pyx_AddTraceback(\"_cython_magic_0ce0e64bf037172e33037b12661fd5c4.klBern_cython\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __pyx_r = NULL;\n",
       "  __pyx_L0:;\n",
       "  __Pyx_XGIVEREF(__pyx_r);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "/* … */\n",
       "  __pyx_tuple_ = PyTuple_Pack(2, __pyx_n_s_x, __pyx_n_s_y); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_tuple_);\n",
       "  __Pyx_GIVEREF(__pyx_tuple_);\n",
       "/* … */\n",
       "  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_0ce0e64bf037172e33037b12661fd5c4_1klBern_cython, NULL, __pyx_n_s_cython_magic_0ce0e64bf037172e33); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_t_1);\n",
       "  if (PyDict_SetItem(__pyx_d, __pyx_n_s_klBern_cython, __pyx_t_1) < 0) __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;\n",
       "
+4:     x = max(1e-7, min(1 - 1e-7, x))
\n", "
  __pyx_t_1 = __pyx_v_x;\n",
       "  __pyx_t_2 = (1.0 - 1e-7);\n",
       "  if (((__pyx_t_1 < __pyx_t_2) != 0)) {\n",
       "    __pyx_t_3 = __pyx_t_1;\n",
       "  } else {\n",
       "    __pyx_t_3 = __pyx_t_2;\n",
       "  }\n",
       "  __pyx_t_1 = __pyx_t_3;\n",
       "  __pyx_t_3 = 1e-7;\n",
       "  if (((__pyx_t_1 > __pyx_t_3) != 0)) {\n",
       "    __pyx_t_2 = __pyx_t_1;\n",
       "  } else {\n",
       "    __pyx_t_2 = __pyx_t_3;\n",
       "  }\n",
       "  __pyx_v_x = __pyx_t_2;\n",
       "
+5:     x = max(1e-7, min(1 - 1e-7, x))
\n", "
  __pyx_t_2 = __pyx_v_x;\n",
       "  __pyx_t_1 = (1.0 - 1e-7);\n",
       "  if (((__pyx_t_2 < __pyx_t_1) != 0)) {\n",
       "    __pyx_t_3 = __pyx_t_2;\n",
       "  } else {\n",
       "    __pyx_t_3 = __pyx_t_1;\n",
       "  }\n",
       "  __pyx_t_2 = __pyx_t_3;\n",
       "  __pyx_t_3 = 1e-7;\n",
       "  if (((__pyx_t_2 > __pyx_t_3) != 0)) {\n",
       "    __pyx_t_1 = __pyx_t_2;\n",
       "  } else {\n",
       "    __pyx_t_1 = __pyx_t_3;\n",
       "  }\n",
       "  __pyx_v_x = __pyx_t_1;\n",
       "
+6:     return x * log(x/y) + (1-x) * log((1-x)/(1-y))
\n", "
  __Pyx_XDECREF(__pyx_r);\n",
       "  if (unlikely(__pyx_v_y == 0)) {\n",
       "    PyErr_SetString(PyExc_ZeroDivisionError, \"float division\");\n",
       "    __PYX_ERR(0, 6, __pyx_L1_error)\n",
       "  }\n",
       "  __pyx_t_1 = (1.0 - __pyx_v_x);\n",
       "  __pyx_t_2 = (1.0 - __pyx_v_y);\n",
       "  if (unlikely(__pyx_t_2 == 0)) {\n",
       "    PyErr_SetString(PyExc_ZeroDivisionError, \"float division\");\n",
       "    __PYX_ERR(0, 6, __pyx_L1_error)\n",
       "  }\n",
       "  __pyx_t_4 = PyFloat_FromDouble(((__pyx_v_x * log((__pyx_v_x / __pyx_v_y))) + ((1.0 - __pyx_v_x) * log((__pyx_t_1 / __pyx_t_2))))); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 6, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_t_4);\n",
       "  __pyx_r = __pyx_t_4;\n",
       "  __pyx_t_4 = 0;\n",
       "  goto __pyx_L0;\n",
       "
" ], "text/plain": [ "" ] }, "execution_count": 26, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%cython --annotate\n", "from libc.math cimport log\n", "\n", "def klBern_cython(double x, double y) -> double:\n", " x = max(1e-7, min(1 - 1e-7, x))\n", " x = max(1e-7, min(1 - 1e-7, x))\n", " return x * log(x/y) + (1-x) * log((1-x)/(1-y))" ] }, { "cell_type": "code", "execution_count": 27, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on built-in function klBern_cython in module _cython_magic_0ce0e64bf037172e33037b12661fd5c4:\n", "\n", "klBern_cython(...)\n", "\n" ] } ], "source": [ "help(klBern_cython)" ] }, { "cell_type": "code", "execution_count": 28, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "922 ns ± 36.6 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)\n" ] } ], "source": [ "%timeit klBern_cython(random_Bern(), random_Bern())" ] }, { "cell_type": "code", "execution_count": 29, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Speed up using Cython for klBern was: 9.1 faster!\n" ] } ], "source": [ "print(f\"Speed up using Cython for klBern was: {(1740-753)/(861-753):.2g} faster!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Threshold functions" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Threshold for GLR Bernoulli" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Generating some data:" ] }, { "cell_type": "code", "execution_count": 30, "metadata": {}, "outputs": [], "source": [ "def random_t0_s_t_delta(min_t: int=100, max_t: int=1000) -> (int, int, int, float):\n", " t0 = 0\n", " t = np.random.randint(min_t, max_t + 1)\n", " s = np.random.randint(t0, t)\n", " delta = np.random.choice([0.1, 0.05, 0.01, 0.005, 0.001, max(0.0005, 1/t)])\n", " return (t0, s, t, delta)" ] }, { "cell_type": "code", "execution_count": 31, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.04 µs ± 148 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit random_t0_s_t_delta()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In pure Python:" ] }, { "cell_type": "code", "execution_count": 32, "metadata": {}, "outputs": [], "source": [ "def threshold(t0: int, s: int, t: int, delta: float) -> float:\n", " return np.log((s - t0 + 1) * (t - s) / delta)" ] }, { "cell_type": "code", "execution_count": 33, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "10.2 µs ± 175 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit threshold(*random_t0_s_t_delta())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "It's *way* faster to use `math.log` instead of `numpy.log` (of course)!" ] }, { "cell_type": "code", "execution_count": 34, "metadata": {}, "outputs": [], "source": [ "from math import log\n", "\n", "def threshold2(t0: int, s: int, t: int, delta: float) -> float:\n", " return log((s - t0 + 1) * (t - s) / delta)" ] }, { "cell_type": "code", "execution_count": 35, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.93 µs ± 132 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit threshold2(*random_t0_s_t_delta())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In numba:" ] }, { "cell_type": "code", "execution_count": 36, "metadata": {}, "outputs": [], "source": [ "from math import log\n", "\n", "@numba.jit(nopython=True)\n", "def threshold_numba(t0: int, s: int, t: int, delta: float) -> float:\n", " return log((s - t0 + 1) * (t - s) / delta)" ] }, { "cell_type": "code", "execution_count": 37, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on CPUDispatcher in module __main__:\n", "\n", "threshold_numba(t0:int, s:int, t:int, delta:float) -> float\n", "\n" ] } ], "source": [ "help(threshold_numba)" ] }, { "cell_type": "code", "execution_count": 38, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.57 µs ± 105 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit threshold_numba(*random_t0_s_t_delta())" ] }, { "cell_type": "code", "execution_count": 39, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Speed up using Cython for thresold was: 0.56 faster!\n" ] } ], "source": [ "print(f\"Speed up using Cython for thresold was: {(7510-7200)/(7750-7200):.2g} faster!\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In Cython:" ] }, { "cell_type": "code", "execution_count": 40, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", " \n", " Cython: _cython_magic_4e48b0e25538ff2883ba5a92832943be.pyx\n", " \n", "\n", "\n", "

Generated by Cython 0.29.1

\n", "

\n", " Yellow lines hint at Python interaction.
\n", " Click on a line that starts with a \"+\" to see the C code that Cython generated for it.\n", "

\n", "
 1: from libc.math cimport log
\n", "
 2:
\n", "
+3: cpdef double threshold_cython(int t0, int s, int t, double delta):
\n", "
static PyObject *__pyx_pw_46_cython_magic_4e48b0e25538ff2883ba5a92832943be_1threshold_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\n",
       "static double __pyx_f_46_cython_magic_4e48b0e25538ff2883ba5a92832943be_threshold_cython(int __pyx_v_t0, int __pyx_v_s, int __pyx_v_t, double __pyx_v_delta, CYTHON_UNUSED int __pyx_skip_dispatch) {\n",
       "  double __pyx_r;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"threshold_cython\", 0);\n",
       "/* … */\n",
       "  /* function exit code */\n",
       "  __pyx_L1_error:;\n",
       "  __Pyx_WriteUnraisable(\"_cython_magic_4e48b0e25538ff2883ba5a92832943be.threshold_cython\", __pyx_clineno, __pyx_lineno, __pyx_filename, 1, 0);\n",
       "  __pyx_r = 0;\n",
       "  __pyx_L0:;\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "/* Python wrapper */\n",
       "static PyObject *__pyx_pw_46_cython_magic_4e48b0e25538ff2883ba5a92832943be_1threshold_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\n",
       "static PyObject *__pyx_pw_46_cython_magic_4e48b0e25538ff2883ba5a92832943be_1threshold_cython(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n",
       "  int __pyx_v_t0;\n",
       "  int __pyx_v_s;\n",
       "  int __pyx_v_t;\n",
       "  double __pyx_v_delta;\n",
       "  PyObject *__pyx_r = 0;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"threshold_cython (wrapper)\", 0);\n",
       "  {\n",
       "    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_t0,&__pyx_n_s_s,&__pyx_n_s_t,&__pyx_n_s_delta,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 = PyTuple_GET_SIZE(__pyx_args);\n",
       "      switch (pos_args) {\n",
       "        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  0: break;\n",
       "        default: goto __pyx_L5_argtuple_error;\n",
       "      }\n",
       "      kw_args = PyDict_Size(__pyx_kwds);\n",
       "      switch (pos_args) {\n",
       "        case  0:\n",
       "        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_t0)) != 0)) kw_args--;\n",
       "        else goto __pyx_L5_argtuple_error;\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1:\n",
       "        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_s)) != 0)) kw_args--;\n",
       "        else {\n",
       "          __Pyx_RaiseArgtupleInvalid(\"threshold_cython\", 1, 4, 4, 1); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "        }\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2:\n",
       "        if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_t)) != 0)) kw_args--;\n",
       "        else {\n",
       "          __Pyx_RaiseArgtupleInvalid(\"threshold_cython\", 1, 4, 4, 2); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "        }\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  3:\n",
       "        if (likely((values[3] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_delta)) != 0)) kw_args--;\n",
       "        else {\n",
       "          __Pyx_RaiseArgtupleInvalid(\"threshold_cython\", 1, 4, 4, 3); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "        }\n",
       "      }\n",
       "      if (unlikely(kw_args > 0)) {\n",
       "        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"threshold_cython\") < 0)) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "      }\n",
       "    } else if (PyTuple_GET_SIZE(__pyx_args) != 4) {\n",
       "      goto __pyx_L5_argtuple_error;\n",
       "    } else {\n",
       "      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
       "      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
       "      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);\n",
       "      values[3] = PyTuple_GET_ITEM(__pyx_args, 3);\n",
       "    }\n",
       "    __pyx_v_t0 = __Pyx_PyInt_As_int(values[0]); if (unlikely((__pyx_v_t0 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "    __pyx_v_s = __Pyx_PyInt_As_int(values[1]); if (unlikely((__pyx_v_s == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "    __pyx_v_t = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_t == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "    __pyx_v_delta = __pyx_PyFloat_AsDouble(values[3]); if (unlikely((__pyx_v_delta == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "  }\n",
       "  goto __pyx_L4_argument_unpacking_done;\n",
       "  __pyx_L5_argtuple_error:;\n",
       "  __Pyx_RaiseArgtupleInvalid(\"threshold_cython\", 1, 4, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "  __pyx_L3_error:;\n",
       "  __Pyx_AddTraceback(\"_cython_magic_4e48b0e25538ff2883ba5a92832943be.threshold_cython\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return NULL;\n",
       "  __pyx_L4_argument_unpacking_done:;\n",
       "  __pyx_r = __pyx_pf_46_cython_magic_4e48b0e25538ff2883ba5a92832943be_threshold_cython(__pyx_self, __pyx_v_t0, __pyx_v_s, __pyx_v_t, __pyx_v_delta);\n",
       "\n",
       "  /* function exit code */\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "static PyObject *__pyx_pf_46_cython_magic_4e48b0e25538ff2883ba5a92832943be_threshold_cython(CYTHON_UNUSED PyObject *__pyx_self, int __pyx_v_t0, int __pyx_v_s, int __pyx_v_t, double __pyx_v_delta) {\n",
       "  PyObject *__pyx_r = NULL;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"threshold_cython\", 0);\n",
       "  __Pyx_XDECREF(__pyx_r);\n",
       "  __pyx_t_1 = PyFloat_FromDouble(__pyx_f_46_cython_magic_4e48b0e25538ff2883ba5a92832943be_threshold_cython(__pyx_v_t0, __pyx_v_s, __pyx_v_t, __pyx_v_delta, 0)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_t_1);\n",
       "  __pyx_r = __pyx_t_1;\n",
       "  __pyx_t_1 = 0;\n",
       "  goto __pyx_L0;\n",
       "\n",
       "  /* function exit code */\n",
       "  __pyx_L1_error:;\n",
       "  __Pyx_XDECREF(__pyx_t_1);\n",
       "  __Pyx_AddTraceback(\"_cython_magic_4e48b0e25538ff2883ba5a92832943be.threshold_cython\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __pyx_r = NULL;\n",
       "  __pyx_L0:;\n",
       "  __Pyx_XGIVEREF(__pyx_r);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "
+4:     return log((s - t0 + 1) * (t - s) / delta)
\n", "
  __pyx_t_1 = (((__pyx_v_s - __pyx_v_t0) + 1) * (__pyx_v_t - __pyx_v_s));\n",
       "  if (unlikely(__pyx_v_delta == 0)) {\n",
       "    PyErr_SetString(PyExc_ZeroDivisionError, \"float division\");\n",
       "    __PYX_ERR(0, 4, __pyx_L1_error)\n",
       "  }\n",
       "  __pyx_r = log((((double)__pyx_t_1) / __pyx_v_delta));\n",
       "  goto __pyx_L0;\n",
       "
" ], "text/plain": [ "" ] }, "execution_count": 40, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%cython --annotate\n", "from libc.math cimport log\n", "\n", "cpdef double threshold_cython(int t0, int s, int t, double delta):\n", " return log((s - t0 + 1) * (t - s) / delta)" ] }, { "cell_type": "code", "execution_count": 41, "metadata": {}, "outputs": [ { "data": { "text/html": [ "\n", "\n", "\n", "\n", " \n", " Cython: _cython_magic_8a5a8f9968b73b8b7b920f92f30405bc.pyx\n", " \n", "\n", "\n", "

Generated by Cython 0.29.1

\n", "

\n", " Yellow lines hint at Python interaction.
\n", " Click on a line that starts with a \"+\" to see the C code that Cython generated for it.\n", "

\n", "
 1: from libc.math cimport log
\n", "
 2:
\n", "
+3: def threshold_cython2(int t0, int s, int t, double delta) -> double:
\n", "
/* Python wrapper */\n",
       "static PyObject *__pyx_pw_46_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc_1threshold_cython2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/\n",
       "static PyMethodDef __pyx_mdef_46_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc_1threshold_cython2 = {\"threshold_cython2\", (PyCFunction)(void*)(PyCFunctionWithKeywords)__pyx_pw_46_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc_1threshold_cython2, METH_VARARGS|METH_KEYWORDS, 0};\n",
       "static PyObject *__pyx_pw_46_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc_1threshold_cython2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {\n",
       "  int __pyx_v_t0;\n",
       "  int __pyx_v_s;\n",
       "  int __pyx_v_t;\n",
       "  double __pyx_v_delta;\n",
       "  PyObject *__pyx_r = 0;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"threshold_cython2 (wrapper)\", 0);\n",
       "  {\n",
       "    static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_t0,&__pyx_n_s_s,&__pyx_n_s_t,&__pyx_n_s_delta,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 = PyTuple_GET_SIZE(__pyx_args);\n",
       "      switch (pos_args) {\n",
       "        case  4: values[3] = PyTuple_GET_ITEM(__pyx_args, 3);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  3: values[2] = PyTuple_GET_ITEM(__pyx_args, 2);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  0: break;\n",
       "        default: goto __pyx_L5_argtuple_error;\n",
       "      }\n",
       "      kw_args = PyDict_Size(__pyx_kwds);\n",
       "      switch (pos_args) {\n",
       "        case  0:\n",
       "        if (likely((values[0] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_t0)) != 0)) kw_args--;\n",
       "        else goto __pyx_L5_argtuple_error;\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  1:\n",
       "        if (likely((values[1] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_s)) != 0)) kw_args--;\n",
       "        else {\n",
       "          __Pyx_RaiseArgtupleInvalid(\"threshold_cython2\", 1, 4, 4, 1); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "        }\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  2:\n",
       "        if (likely((values[2] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_t)) != 0)) kw_args--;\n",
       "        else {\n",
       "          __Pyx_RaiseArgtupleInvalid(\"threshold_cython2\", 1, 4, 4, 2); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "        }\n",
       "        CYTHON_FALLTHROUGH;\n",
       "        case  3:\n",
       "        if (likely((values[3] = __Pyx_PyDict_GetItemStr(__pyx_kwds, __pyx_n_s_delta)) != 0)) kw_args--;\n",
       "        else {\n",
       "          __Pyx_RaiseArgtupleInvalid(\"threshold_cython2\", 1, 4, 4, 3); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "        }\n",
       "      }\n",
       "      if (unlikely(kw_args > 0)) {\n",
       "        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, \"threshold_cython2\") < 0)) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "      }\n",
       "    } else if (PyTuple_GET_SIZE(__pyx_args) != 4) {\n",
       "      goto __pyx_L5_argtuple_error;\n",
       "    } else {\n",
       "      values[0] = PyTuple_GET_ITEM(__pyx_args, 0);\n",
       "      values[1] = PyTuple_GET_ITEM(__pyx_args, 1);\n",
       "      values[2] = PyTuple_GET_ITEM(__pyx_args, 2);\n",
       "      values[3] = PyTuple_GET_ITEM(__pyx_args, 3);\n",
       "    }\n",
       "    __pyx_v_t0 = __Pyx_PyInt_As_int(values[0]); if (unlikely((__pyx_v_t0 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "    __pyx_v_s = __Pyx_PyInt_As_int(values[1]); if (unlikely((__pyx_v_s == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "    __pyx_v_t = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_t == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "    __pyx_v_delta = __pyx_PyFloat_AsDouble(values[3]); if (unlikely((__pyx_v_delta == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "  }\n",
       "  goto __pyx_L4_argument_unpacking_done;\n",
       "  __pyx_L5_argtuple_error:;\n",
       "  __Pyx_RaiseArgtupleInvalid(\"threshold_cython2\", 1, 4, 4, PyTuple_GET_SIZE(__pyx_args)); __PYX_ERR(0, 3, __pyx_L3_error)\n",
       "  __pyx_L3_error:;\n",
       "  __Pyx_AddTraceback(\"_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc.threshold_cython2\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return NULL;\n",
       "  __pyx_L4_argument_unpacking_done:;\n",
       "  __pyx_r = __pyx_pf_46_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc_threshold_cython2(__pyx_self, __pyx_v_t0, __pyx_v_s, __pyx_v_t, __pyx_v_delta);\n",
       "\n",
       "  /* function exit code */\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "\n",
       "static PyObject *__pyx_pf_46_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc_threshold_cython2(CYTHON_UNUSED PyObject *__pyx_self, int __pyx_v_t0, int __pyx_v_s, int __pyx_v_t, double __pyx_v_delta) {\n",
       "  PyObject *__pyx_r = NULL;\n",
       "  __Pyx_RefNannyDeclarations\n",
       "  __Pyx_RefNannySetupContext(\"threshold_cython2\", 0);\n",
       "/* … */\n",
       "  /* function exit code */\n",
       "  __pyx_L1_error:;\n",
       "  __Pyx_XDECREF(__pyx_t_2);\n",
       "  __Pyx_AddTraceback(\"_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc.threshold_cython2\", __pyx_clineno, __pyx_lineno, __pyx_filename);\n",
       "  __pyx_r = NULL;\n",
       "  __pyx_L0:;\n",
       "  __Pyx_XGIVEREF(__pyx_r);\n",
       "  __Pyx_RefNannyFinishContext();\n",
       "  return __pyx_r;\n",
       "}\n",
       "/* … */\n",
       "  __pyx_tuple_ = PyTuple_Pack(4, __pyx_n_s_t0, __pyx_n_s_s, __pyx_n_s_t, __pyx_n_s_delta); if (unlikely(!__pyx_tuple_)) __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_tuple_);\n",
       "  __Pyx_GIVEREF(__pyx_tuple_);\n",
       "/* … */\n",
       "  __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_46_cython_magic_8a5a8f9968b73b8b7b920f92f30405bc_1threshold_cython2, NULL, __pyx_n_s_cython_magic_8a5a8f9968b73b8b7b); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_t_1);\n",
       "  if (PyDict_SetItem(__pyx_d, __pyx_n_s_threshold_cython2, __pyx_t_1) < 0) __PYX_ERR(0, 3, __pyx_L1_error)\n",
       "  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;\n",
       "
+4:     return log((s - t0 + 1) * (t - s) / delta)
\n", "
  __Pyx_XDECREF(__pyx_r);\n",
       "  __pyx_t_1 = (((__pyx_v_s - __pyx_v_t0) + 1) * (__pyx_v_t - __pyx_v_s));\n",
       "  if (unlikely(__pyx_v_delta == 0)) {\n",
       "    PyErr_SetString(PyExc_ZeroDivisionError, \"float division\");\n",
       "    __PYX_ERR(0, 4, __pyx_L1_error)\n",
       "  }\n",
       "  __pyx_t_2 = PyFloat_FromDouble(log((((double)__pyx_t_1) / __pyx_v_delta))); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 4, __pyx_L1_error)\n",
       "  __Pyx_GOTREF(__pyx_t_2);\n",
       "  __pyx_r = __pyx_t_2;\n",
       "  __pyx_t_2 = 0;\n",
       "  goto __pyx_L0;\n",
       "
" ], "text/plain": [ "" ] }, "execution_count": 41, "metadata": {}, "output_type": "execute_result" } ], "source": [ "%%cython --annotate\n", "from libc.math cimport log\n", "\n", "def threshold_cython2(int t0, int s, int t, double delta) -> double:\n", " return log((s - t0 + 1) * (t - s) / delta)" ] }, { "cell_type": "code", "execution_count": 42, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Help on built-in function threshold_cython in module _cython_magic_4e48b0e25538ff2883ba5a92832943be:\n", "\n", "threshold_cython(...)\n", "\n" ] } ], "source": [ "help(threshold_cython)" ] }, { "cell_type": "code", "execution_count": 43, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "7.24 µs ± 103 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit threshold_cython(*random_t0_s_t_delta())" ] }, { "cell_type": "code", "execution_count": 44, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "8.8 µs ± 994 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)\n" ] } ], "source": [ "%timeit threshold_cython2(*random_t0_s_t_delta())" ] }, { "cell_type": "code", "execution_count": 45, "metadata": { "scrolled": true }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Speed up using Cython for thresold was: 2.4 faster!\n" ] } ], "source": [ "print(f\"Speed up using Cython for thresold was: {abs(7510-7200)/abs(7070-7200):.2g} faster!\")" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.7" }, "toc": { "colors": { "hover_highlight": "#DAA520", "running_highlight": "#FF0000", "selected_highlight": "#FFD700" }, "moveMenuLeft": true, "nav_menu": { "height": "162.683px", "width": "251.833px" }, "navigate_menu": true, "number_sections": true, "sideBar": false, "threshold": 4, "toc_cell": true, "toc_section_display": "block", "toc_window_display": true } }, "nbformat": 4, "nbformat_minor": 2 }