SciPy Cython Optimize Zeros API

cython_optimize_brentq_benchmark.ipynb

cython_optimize_brentq_example.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "from scipy.optimize.cython_optimize cimport brentq\n",
    "\n",
    "# import math from Cython\n",
    "from libc cimport math\n",
    "\n",
    "ARGS = {'C0': 1.0, 'C1': 0.7}  # a dictionary of extra arguments\n",
    "XLO, XHI = 0.5, 1.0  # lower and upper search boundaries\n",
    "XTOL, RTOL, MITR = 1e-3, 1e-3, 10  # other solver parameters\n",
    "\n",
    "# user defined struct for extra parameters\n",
    "ctypedef struct test_params:\n",
    "    double C0\n",
    "    double C1\n",
    "\n",
    "\n",
    "# user defined callback\n",
    "cdef double f(double x, void *args):\n",
    "    cdef test_params *myargs = <test_params *> args\n",
    "    return myargs.C0 - math.exp(-(x - myargs.C1))\n",
    "\n",
    "\n",
    "# Cython wrapper function\n",
    "cdef double brentq_wrapper_example(dict args, double xa, double xb,\n",
    "                                   double xtol, double rtol, int mitr):\n",
    "    # Cython automatically casts dictionary to struct\n",
    "    cdef test_params myargs = args\n",
    "    return brentq(\n",
    "        f, xa, xb, <test_params *> &myargs, xtol, rtol, mitr, NULL)\n",
    "\n",
    "\n",
    "# Python function\n",
    "def brentq_example(args=ARGS, xa=XLO, xb=XHI, xtol=XTOL, rtol=RTOL,\n",
    "                   mitr=MITR):\n",
    "    '''Calls Cython wrapper from Python.'''\n",
    "    return brentq_wrapper_example(args, xa, xb, xtol, rtol, mitr)\n",
    "\n",
    "\n",
    "x = brentq_example()\n",
    "# 0.6999942848231314"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.6999942848231314"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "from scipy.optimize.cython_optimize cimport brentq\n",
    "\n",
    "# import math from Cython\n",
    "from libc cimport math\n",
    "\n",
    "ARGS = {'C0': 1.0, 'C1': 0.7}  # a dictionary of extra arguments\n",
    "XLO, XHI = 0.5, 1.0  # lower and upper search boundaries\n",
    "XTOL, RTOL, MITR = 1e-3, 1e-3, 10  # other solver parameters\n",
    "\n",
    "# user defined struct for extra parameters\n",
    "ctypedef struct test_params:\n",
    "    double C0\n",
    "    double C1\n",
    "\n",
    "\n",
    "# user defined callback\n",
    "cdef double f(double x, void *args):\n",
    "    cdef test_params *myargs = <test_params *> args\n",
    "    return myargs.C0 - math.exp(-(x - myargs.C1))\n",
    "\n",
    "\n",
    "# Cython wrapper function\n",
    "cdef double brentq_wrapper_example(dict args, double xa, double xb,\n",
    "                                   double xtol, double rtol, int mitr):\n",
    "    # Cython automatically casts dictionary to struct\n",
    "    cdef test_params myargs = args\n",
    "    return brentq(\n",
    "        f, xa, xb, <test_params *> &myargs, xtol, rtol, mitr, NULL)\n",
    "\n",
    "\n",
    "# Python function\n",
    "def brentq_example(args=ARGS, xa=XLO, xb=XHI, xtol=XTOL, rtol=RTOL,\n",
    "                   mitr=MITR):\n",
    "    '''Calls Cython wrapper from Python.'''\n",
    "    return brentq_wrapper_example(args, xa, xb, xtol, rtol, mitr)\n",
    "\n",
    "\n",
    "x = brentq_example()\n",
    "# 0.6999942848231314\n",
    "\n",
    "\n",
    "from scipy.optimize.cython_optimize cimport zeros_full_output\n",
    "\n",
    "\n",
    "# cython brentq solver with full output\n",
    "cdef zeros_full_output brentq_full_output_wrapper_example(\n",
    "        dict args, double xa, double xb, double xtol, double rtol,\n",
    "        int mitr):\n",
    "    cdef test_params myargs = args\n",
    "    cdef zeros_full_output my_full_output\n",
    "    # use my_full_output instead of NULL\n",
    "    brentq(f, xa, xb, &myargs, xtol, rtol, mitr, &my_full_output)\n",
    "    return my_full_output\n",
    "\n",
    "\n",
    "# Python function\n",
    "def brent_full_output_example(args=ARGS, xa=XLO, xb=XHI, xtol=XTOL,\n",
    "                              rtol=RTOL, mitr=MITR):\n",
    "    '''Returns full output'''\n",
    "    return brentq_full_output_wrapper_example(args, xa, xb, xtol, rtol,\n",
    "                                              mitr)\n",
    "\n",
    "result = brent_full_output_example()\n",
    "# {'error_num': 0,\n",
    "#  'funcalls': 6,\n",
    "#  'iterations': 5,\n",
    "#  'root': 0.6999942848231314}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{'error_num': 0, 'funcalls': 6, 'iterations': 5, 'root': 0.6999942848231314}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "import numpy as np\n",
    "from scipy.optimize.cython_optimize cimport brentq\n",
    "\n",
    "# import math from Cython\n",
    "from libc cimport math\n",
    "\n",
    "CNT = 1e5\n",
    "C0_ARG = 1.0 + np.linspace(-0.05, 0.05, int(CNT), dtype=np.dtype(\"f8\"))\n",
    "C1_ARG = 0.7\n",
    "XLO, XHI = 0.0, 1.0  # lower and upper search boundaries\n",
    "XTOL, RTOL = 1e-3, 1e-3\n",
    "MITR = 10  # other solver parameters\n",
    "\n",
    "# user defined struct for extra parameters\n",
    "ctypedef struct test_params:\n",
    "    double C0\n",
    "    double C1\n",
    "\n",
    "\n",
    "# user defined callback\n",
    "cdef double f(double x, void *args):\n",
    "    cdef test_params *myargs = <test_params *> args\n",
    "    return myargs.C0 - math.exp(-(x - myargs.C1))\n",
    "\n",
    "\n",
    "# Cython wrapper function\n",
    "cdef int brentq_wrapper_loop_example(double[:] c0, double c1, double xa, double xb,\n",
    "                                     double xtol, double rtol, int mitr, double[:] result):\n",
    "    cdef Py_ssize_t i = 0\n",
    "    cdef test_params myargs\n",
    "\n",
    "    myargs.C1 = c1\n",
    "\n",
    "    #for i in prange(c0.shape[0], nogil=True):\n",
    "    while i < CNT:\n",
    "        myargs.C0 = c0[i]\n",
    "        result[i] = brentq(\n",
    "            f, xa, xb, <test_params *> &myargs, xtol, rtol, mitr, NULL)\n",
    "        i += 1\n",
    "    return 0\n",
    "\n",
    "\n",
    "# Python function\n",
    "def brentq_loop_example(c0=C0_ARG, c1=C1_ARG, xa=XLO, xb=XHI, xtol=XTOL, rtol=RTOL,\n",
    "                        mitr=MITR):\n",
    "    '''Calls Cython wrapper from Python.'''\n",
    "    cdef double[:] result = np.empty_like(c0)\n",
    "    if not brentq_wrapper_loop_example(c0, c1, xa, xb, xtol, rtol, mitr, result):\n",
    "        return result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "27.9 ms ± 447 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit brentq_loop_example()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.75149248, 0.75149143, 0.75149038, 0.75148933, 0.75148828,\n",
       "       0.75148723, 0.75148618, 0.75148513, 0.75148408, 0.75148303])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x_loop = brentq_loop_example()\n",
    "np.asarray(x_loop)[:10]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython\n",
    "\n",
    "from scipy.optimize.cython_optimize cimport brentq\n",
    "import numpy as np\n",
    "\n",
    "# import math from Cython\n",
    "from libc cimport math\n",
    "\n",
    "CNT = 1e5\n",
    "C0_ARG = 1.0 + np.linspace(-0.05, 0.05, int(CNT), dtype=np.dtype(\"f8\"))\n",
    "C1_ARG = 0.7\n",
    "XLO, XHI = 0., 1.0  # lower and upper search boundaries\n",
    "XTOL, RTOL, MITR = 1e-3, 1e-3, 10  # other solver parameters\n",
    "\n",
    "# user defined struct for extra parameters\n",
    "ctypedef struct test_params:\n",
    "    double C0\n",
    "    double C1\n",
    "\n",
    "\n",
    "# user defined callback\n",
    "cdef double f(double x, void *args):\n",
    "    cdef test_params *myargs = <test_params *> args\n",
    "    return myargs.C0 - math.exp(-(x - myargs.C1))\n",
    "\n",
    "\n",
    "# Cython wrapper function\n",
    "cdef double brentq_wrapper_map_example(dict args, double xa, double xb,\n",
    "                                       double xtol, double rtol, int mitr):\n",
    "    # Cython automatically casts dictionary to struct\n",
    "    cdef test_params myargs = args\n",
    "    return brentq(\n",
    "        f, xa, xb, <test_params *> &myargs, xtol, rtol, mitr, NULL)\n",
    "\n",
    "\n",
    "# Python function\n",
    "def brentq_map_example(c0=C0_ARG, c1=C1_ARG, xa=XLO, xb=XHI, xtol=XTOL, rtol=RTOL,\n",
    "                       mitr=MITR):\n",
    "    '''Calls Cython wrapper from Python.'''\n",
    "    args = [{'C0': c0_, 'C1': c1} for c0_ in c0]\n",
    "    return map(lambda a: brentq_wrapper_map_example(a, xa, xb, xtol, rtol, mitr), args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "68.9 ms ± 1.83 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n"
     ]
    }
   ],
   "source": [
    "%timeit list(brentq_map_example())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_map = list(brentq_map_example())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.allclose(np.asarray(x_map), np.asarray(x_loop))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.optimize import brentq\n",
    "import numpy as np\n",
    "\n",
    "CNT = 1e5\n",
    "C0_ARG = 1.0 + np.linspace(-0.05, 0.05, int(CNT), np.float64)\n",
    "C1_ARG = 0.7\n",
    "XLO, XHI = 0., 1.0  # lower and upper search boundaries\n",
    "XTOL, RTOL, MITR = 1e-3, 1e-3, 10  # other solver parameters\n",
    "\n",
    "# user defined callback\n",
    "def f(x, args):\n",
    "    return args['C0'] - np.exp(-(x - args['C1']))\n",
    "\n",
    "\n",
    "# Python function\n",
    "def brentq_np_example(c0=C0_ARG, c1=C1_ARG, xa=XLO, xb=XHI, xtol=XTOL, rtol=RTOL,\n",
    "                      mitr=MITR):\n",
    "    '''Calls Cython wrapper from Python.'''\n",
    "    args = [{'C0': c0_, 'C1': c1} for c0_ in c0]\n",
    "    return map(lambda a: brentq(f, xa, xb, a, xtol, rtol, mitr), args)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "719 ms ± 16 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)\n"
     ]
    }
   ],
   "source": [
    "%timeit list(brentq_np_example())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_np = list(brentq_np_example())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "assert np.allclose(np.asarray(x_map), np.asarray(x_np))"
   ]
  },
  {
   "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.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}

cython_optimize_example.pyx

from scipy.optimize.cython_optimize cimport brentq

# import math from Cython
from libc cimport math


ARGS = {'C0': 1.0, 'C1': 0.7}  # a dictionary of extra arguments
XLO, XHI = 0.5, 1.0  # lower and upper search boundaries
XTOL, RTOL, MITR = 1e-3, 1e-3, 10  # other solver parameters

# user defined struct for extra parameters
ctypedef struct test_params:
   double C0
   double C1


# user defined callback
cdef double f(double x, void *args):
   cdef test_params *myargs = <test_params *> args
   return myargs.C0 - math.exp(-(x - myargs.C1))


# Cython wrapper function
cdef double brentq_wrapper_example(dict args, double xa, double xb,
                                  double xtol, double rtol, int mitr):
   # Cython automatically casts dictionary to struct
   cdef test_params myargs = args  
   return brentq(
       f, xa, xb, <test_params *> &myargs, xtol, rtol, mitr, NULL)


# Python function
def brentq_example(args=ARGS, xa=XLO, xb=XHI, xtol=XTOL, rtol=RTOL,
                  mitr=MITR):
   '''Calls Cython wrapper from Python.'''
   return brentq_wrapper_example(args, xa, xb, xtol, rtol, mitr)


x = brentq_example()
# 0.6999942848231314


from scipy.optimize.cython_optimize cimport zeros_full_output


# cython brentq solver with full output
cdef zeros_full_output brentq_full_output_wrapper_example(
        dict args, double xa, double xb, double xtol, double rtol,
        int mitr):
    cdef test_params myargs = args
    cdef zeros_full_output my_full_output  # simplified output
    # put result into full_output
    brentq(f, xa, xb, &myargs, xtol, rtol, mitr, &my_full_output)
    return my_full_output


# Python function
def brent_full_output_example(args=ARGS, xa=XLO, xb=XHI, xtol=XTOL,
                              rtol=RTOL, mitr=MITR):
    '''Returns full output'''
    return brentq_full_output_wrapper_example(args, xa, xb, xtol, rtol,
                                              mitr)

result = brent_full_output_example()
# {'error_num': 0,
#  'funcalls': 6,
#  'iterations': 5,
#  'root': 0.6999942848231314}