<commit_before>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Easy \"Hard\" Way: Cythonizing\n",
    "\n",
    "[Cython](http://cython.org/) is a **compiler** and a **programming language** used to generate C extensions for Python.\n",
    "\n",
    "The Cython language is a Python/C [creole](https://en.wikipedia.org/wiki/Creole_language) which is essentially Python with some additional keywords for specifying static data types. It looks something like this:\n",
    "\n",
    "```cython\n",
    "def cython_sum(int n):\n",
    "    cdef float s = 0.0\n",
    "    cdef int i\n",
    "    for i in range(n):\n",
    "        s += i\n",
    "    return s\n",
    "```\n",
    "\n",
    "You can write in this language, then use the Cython compiler to *transpile* to efficient C code which can in turn be compiled into a Python extension module. This extension module can be imported like a normal Python module, but it is backed by compiled C, offering potentially drastic speedups over equivalent functions written in pure Python.\n",
    "\n",
    "The other major use case for Cython -- the one we will focus on here -- is writing wrappers around existing C code so that the functions therein can be made available in an extension module as described above. We will use this technique to make the SymPy-generated C code accessible to Python for use in SciPy's `odeint`.\n",
    "\n",
    "**Learning Objectives**\n",
    "\n",
    "After this lesson, you will be able to:\n",
    "\n",
    "- write a simple Cython function and run it a Jupyter notebook using the `%%cython` magic command\n",
    "- use Cython to wrap a function implemented in plain C and make it callable from Python\n",
    "- use the SymPy `codegen` function to output compilable C code\n",
    "- wrap `codegen`-generated code with Cython, compile it into an extension module, and call it from Python\n",
    "- use SymPy's `autowrap` function to do all of this behind the scenes\n",
    "- pass a custom code printer to `autowrap` to make use of an efficient C library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sympy as sm\n",
    "sm.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction to Cython\n",
    "\n",
    "Cython is used in two scenarios:\n",
    "\n",
    "1. writing Cython computation code (i.e. Python with statically typed variables) to generate fast extension modules\n",
    "2. writing Cython wrappers around existing C code\n",
    "\n",
    "In both cases, efficient C code can be generated and made available via import from Python.\n",
    "\n",
    "### Writing Cython Code\n",
    "\n",
    "In the first use case, you write in the Cython language, compile this into CPython C code, then compile it into a Python extension module. There are several ways to ways to go about the compilation process, and in many cases, Cython's tooling makes it fairly simple. For example, Jupyter notebooks can make use of a `%%cython` magic command that will do all of compilation in the background for us. To make use of it, we need to load the `cython` extension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can write a Cython function. Note that the `--annotate` (or `-a`) flag of the `%%cython` magic command will produce an interactive annotated printout of the Cython code, allowing us to see the C code that is generated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython --annotate\n",
    "def cython_fib(int n):\n",
    "    cdef int i\n",
    "    cdef double a = 0.0, b = 1.0\n",
    "    for i in range(n):\n",
    "        a, b = a + b, a\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit cython_fib(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: remove the type declarations to turn the Cython code above into a plain Python function and use `%timeit` to compare performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def python_fib(n):\n",
    "    a, b = 0.0, 1.0\n",
    "    for i in range(n):\n",
    "        a, b = a + b, a\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit python_fib(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see a bit more about writing Cython and its potential performance benefits, see [this Cython examples notebook](cython-examples.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Wrapping C Code using Cython\n",
    "\n",
    "Our main goal in using Cython is to make numerical code generated by SymPy available to Python. To achieve this, we'll focus on the second use case for Cython -- wrapping existing C code. Wrapping a set of C functions involves wrapping existing C code with a fairly simple Cython script that specifies the Python interface to the C functions. This script must do two things:\n",
    "\n",
    "1. specify the function signatures as found in the C header file\n",
    "2. implement the Python interface to the functions by wrapping them\n",
    "\n",
    "### Example\n",
    "\n",
    "In this example, we'll walk through wrapping a simple C library that generates an array of the values in the [Fibonacci sequence](https://en.wikipedia.org/wiki/Fibonacci_number):\n",
    "\n",
    "$$\n",
    "F_n = F_{n-1} + F_{n-2}\n",
    "$$\n",
    "\n",
    "initialized with $F_0 = 0$ and $F_1 = 1$. Here is the C code that does this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cat ../scipy2017codegen/cfib/cfib.c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A function is provided in our tutorial package to simplify the build configuration so that our `cfib.{c,h}` files are found, built, and linked correctly with our wrapper. It takes the name of the module we're creating, then some keyword arguments:\n",
    "\n",
    "- `sources`: C source files that need to be built\n",
    "- `include_dirs`: locations of header files. In our case, both the `cfib.h` file and the NumPy headers need to be available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a path to the cfib source\n",
    "import os\n",
    "import scipy2017codegen\n",
    "cfib_path = os.path.join(os.path.dirname(scipy2017codegen.__file__), 'cfib')\n",
    "\n",
    "# generate a \"pyxbld\" file that specifies how to build everything\n",
    "from scipy2017codegen.templates import render_pyxbld\n",
    "render_pyxbld('fib',\n",
    "              sources=[os.path.join(cfib_path, 'cfib.c')],\n",
    "              include_dirs=[np.get_include(), cfib_path])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can write our wrapper code.\n",
    "\n",
    "We make use of the `%%cython_pyximport` magic function, which is able to read the Cython code in the cell, compile it using our `pyxbld` file, then import its functions as though we had just written a Python function in a normal code cell.\n",
    "\n",
    "To write the wrapper, you first write the function signature as specified by the C library. Then, you create a *wrapper function* that makes use of the implementation and returns the result. This wrapper function becomes the interface to the compiled code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython_pyximport fib\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "\n",
    "# here we just replicate the function signature from the header\n",
    "cdef extern from \"cfib.h\":\n",
    "    void cfib(int n, double *x)\n",
    "    \n",
    "# here is the \"wrapper\"\n",
    "def fib(int n):\n",
    "    # preallocate an array of doubles\n",
    "    cdef np.ndarray[np.float64_t, ndim=1, mode='c'] x\n",
    "    x = np.empty(n, dtype=np.float64)\n",
    "    \n",
    "    # call the c function and return the result\n",
    "    cfib(n, <double *> x.data)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can call the function `fib` to check that everything worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fib(10)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `codegen()` to Generate C Code\n",
    "\n",
    "In the previous notebook focusing on generating C code using SymPy's codegen capabilities, we saw how to form a Jacobian from a system of ODEs, then output C code that can numerically evaluate the right hand side of the ODEs themselves as well as the Jacobian, given the current state vector. Here we'll do that again so that we can wrap it in Cython for use in `scipy.integrate.odeint`.\n",
    "\n",
    "Start by importing the system of ODEs and compute the Jacobian analytically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy2017codegen.chem import load_large_ode\n",
    "\n",
    "rhs_of_odes, states = load_large_ode()\n",
    "jac_of_odes = rhs_of_odes.jacobian(states)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll use `sympy.utilities.codegen.codegen` to output C source and header files which can compute the right hand side (RHS) of the ODEs and the Jacobian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy.utilities.codegen import codegen\n",
    "\n",
    "# form equations so we can put the OutputArgument in the argument_sequence\n",
    "dY = sm.MatrixSymbol('dY', len(states), 1)\n",
    "dY_eq = sm.Eq(dY, rhs_of_odes)\n",
    "t = sm.symbols('t')\n",
    "\n",
    "codegen([('c_odes', dY_eq)], 'c', argument_sequence=[*list(states), t, dY], to_files=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the generated C code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cat c_odes.c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wrapping the Generated Code with Cython\n",
    "\n",
    "Now we want to wrap the function that was generated `c_odes` with a Cython function. Let's output just the header `c_odes.h` to see just the function signature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cat c_odes.h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this function signature is what we need to wrap using Cython. Just as with the Fibonacci sequence example, we start by specifying how the extension module will be built."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "render_pyxbld('cy_odes',\n",
    "              sources=['c_odes.c'],\n",
    "              include_dirs=[np.get_include()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Fill in the wrapper implementation below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython_pyximport cy_odes\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "\n",
    "cdef extern from \"c_odes.h\":\n",
    "    void c_odes(double y0, double y1, double y2, double y3, double y4,\n",
    "                double y5, double y6, double y7, double y8, double y9,\n",
    "                double y10, double y11, double y12, double y13, double t,\n",
    "                double *dY)\n",
    "    \n",
    "def cy_odes(y, t):\n",
    "    cdef np.ndarray[np.float64_t, ndim=1, mode='c'] dY\n",
    "    dY = np.empty(y.size, dtype=np.float64)\n",
    "    c_odes(y[0], y[1], y[2], y[3], y[4], y[5], y[6], y[7], y[8], y[9],\n",
    "           y[10], y[11], y[12], y[13], t, <double *> dY.data)\n",
    "    return dY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use `odeint` to integrate the equations and plot the results to check that it worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "from scipy2017codegen.chem import watrad_init, watrad_plot\n",
    "\n",
    "y_init, t_vals = watrad_init()\n",
    "y_vals = odeint(cy_odes, y_init, t_vals)\n",
    "watrad_plot(t_vals, y_vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Generating and Compiling a C Extension Module Automatically\n",
    "\n",
    "Here we'll use SymPy's `autowrap` function to automatically generate the C code for an expression, generate a Cython wrapper, compile the wrapper into an extension module, and make the function available as a callable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy.utilities.autowrap import autowrap\n",
    "\n",
    "# need to specify input args so they are in correct order\n",
    "args = list(states)\n",
    "args.append(t)\n",
    "\n",
    "f_rhs = autowrap(rhs_of_odes, args=args, backend='cython', tempdir='./autowraptmp')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`autowrap` essentially does exactly what we just did manually: generate C code to evaluate the right hand side of the system of ODEs, write a Cython wrapper, compile the C code and wrapper into an extension module, and import it for use in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: just like the `codegen`, `autowrap` generates a C function that takes each of the state variables as separate arguments, as seen in the generated Cython wrapper file seen above. In order to use the callable returned by `autowrap` within `odeint`, write a Python wrapper that takes the state vector `y` and the current time `t` and uses `f_rhs` to compute the derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rhs_eval(y, t):\n",
    "    yout = f_rhs(*y, t)\n",
    "    return yout[:, 0]\n",
    "\n",
    "y_vals = odeint(rhs_eval, y_init, t_vals)\n",
    "watrad_plot(t_vals, y_vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using a Custom Library\n",
    "\n",
    "In our set of ODEs, there are quite a few instances of $y_i^2$, where $y_i$ is one of the 14 state variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy.utilities.codegen import C99CodeGen\n",
    "from sympy.printing.ccode import C99CodePrinter\n",
    "\n",
    "class CustomPrinter(C99CodePrinter):\n",
    "    def _print_Pow(self, expr):\n",
    "        return 'fastpow({}, {})'.format(expr.base, expr.exp)\n",
    "\n",
    "printer = CustomPrinter()\n",
    "gen = C99CodeGen(printer=printer)\n",
    "gen.preprocessor_statements.append('#include \"fastpow.h\"')\n",
    "\n",
    "import scipy2017codegen\n",
    "import os\n",
    "fastapprox_dir = os.path.join(os.path.dirname(scipy2017codegen.__file__), 'fastapprox')\n",
    "\n",
    "f_rhs = autowrap(rhs_of_odes, args=args, code_gen=gen, backend='cython',\n",
    "                 tempdir='./autowraptmp2', include_dirs=[fastapprox_dir])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_vals, info = odeint(rhs_eval, y_init, t_vals, full_output=True)\n",
    "watrad_plot(t_vals, y_vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "info"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,<commit_msg>Minor fix.<commit_after>{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The Easy \"Hard\" Way: Cythonizing\n",
    "\n",
    "[Cython](http://cython.org/) is a **compiler** and a **programming language** used to generate C extensions for Python.\n",
    "\n",
    "The Cython language is a Python/C [creole](https://en.wikipedia.org/wiki/Creole_language) which is essentially Python with some additional keywords for specifying static data types. It looks something like this:\n",
    "\n",
    "```cython\n",
    "def cython_sum(int n):\n",
    "    cdef float s = 0.0\n",
    "    cdef int i\n",
    "    for i in range(n):\n",
    "        s += i\n",
    "    return s\n",
    "```\n",
    "\n",
    "You can write in this language, then use the Cython compiler to *transpile* to efficient C code which can in turn be compiled into a Python extension module. This extension module can be imported like a normal Python module, but it is backed by compiled C, offering potentially drastic speedups over equivalent functions written in pure Python.\n",
    "\n",
    "The other major use case for Cython -- the one we will focus on here -- is writing wrappers around existing C code so that the functions therein can be made available in an extension module as described above. We will use this technique to make the SymPy-generated C code accessible to Python for use in SciPy's `odeint`.\n",
    "\n",
    "**Learning Objectives**\n",
    "\n",
    "After this lesson, you will be able to:\n",
    "\n",
    "- write a simple Cython function and run it a Jupyter notebook using the `%%cython` magic command\n",
    "- use Cython to wrap a function implemented in plain C and make it callable from Python\n",
    "- use the SymPy `codegen` function to output compilable C code\n",
    "- wrap `codegen`-generated code with Cython, compile it into an extension module, and call it from Python\n",
    "- use SymPy's `autowrap` function to do all of this behind the scenes\n",
    "- pass a custom code printer to `autowrap` to make use of an efficient C library"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sympy as sm\n",
    "sm.init_printing()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Introduction to Cython\n",
    "\n",
    "Cython is used in two scenarios:\n",
    "\n",
    "1. writing Cython computation code (i.e. Python with statically typed variables) to generate fast extension modules\n",
    "2. writing Cython wrappers around existing C code\n",
    "\n",
    "In both cases, efficient C code can be generated and made available via import from Python.\n",
    "\n",
    "### Writing Cython Code\n",
    "\n",
    "In the first use case, you write in the Cython language, compile this into CPython C code, then compile it into a Python extension module. There are several ways to ways to go about the compilation process, and in many cases, Cython's tooling makes it fairly simple. For example, Jupyter notebooks can make use of a `%%cython` magic command that will do all of compilation in the background for us. To make use of it, we need to load the `cython` extension."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%load_ext cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can write a Cython function. Note that the `--annotate` (or `-a`) flag of the `%%cython` magic command will produce an interactive annotated printout of the Cython code, allowing us to see the C code that is generated."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%cython --annotate\n",
    "def cython_fib(int n):\n",
    "    cdef int i\n",
    "    cdef double a = 0.0, b = 1.0\n",
    "    for i in range(n):\n",
    "        a, b = a + b, a\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit cython_fib(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: remove the type declarations to turn the Cython code above into a plain Python function and use `%timeit` to compare performance"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "def python_fib(n):\n",
    "    a, b = 0.0, 1.0\n",
    "    for i in range(n):\n",
    "        a, b = a + b, a\n",
    "    return a"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%timeit python_fib(100)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see a bit more about writing Cython and its potential performance benefits, see [this Cython examples notebook](cython-examples.ipynb)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Wrapping C Code using Cython\n",
    "\n",
    "Our main goal in using Cython is to make numerical code generated by SymPy available to Python. To achieve this, we'll focus on the second use case for Cython -- wrapping existing C code. Wrapping a set of C functions involves wrapping existing C code with a fairly simple Cython script that specifies the Python interface to the C functions. This script must do two things:\n",
    "\n",
    "1. specify the function signatures as found in the C header file\n",
    "2. implement the Python interface to the functions by wrapping them\n",
    "\n",
    "### Example\n",
    "\n",
    "In this example, we'll walk through wrapping a simple C library that generates an array of the values in the [Fibonacci sequence](https://en.wikipedia.org/wiki/Fibonacci_number):\n",
    "\n",
    "$$\n",
    "F_i = F_{i-1} + F_{i-2}\n",
    "$$\n",
    "\n",
    "initialized with $F_0 = 0$ and $F_1 = 1$. Here is the C code that does this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cat ../scipy2017codegen/cfib/cfib.c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A function is provided in our tutorial package to simplify the build configuration so that our `cfib.{c,h}` files are found, built, and linked correctly with our wrapper. It takes the name of the module we're creating, then some keyword arguments:\n",
    "\n",
    "- `sources`: C source files that need to be built\n",
    "- `include_dirs`: locations of header files. In our case, both the `cfib.h` file and the NumPy headers need to be available."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# get a path to the cfib source\n",
    "import os\n",
    "import scipy2017codegen\n",
    "cfib_path = os.path.join(os.path.dirname(scipy2017codegen.__file__), 'cfib')\n",
    "\n",
    "# generate a \"pyxbld\" file that specifies how to build everything\n",
    "from scipy2017codegen.templates import render_pyxbld\n",
    "render_pyxbld('fib',\n",
    "              sources=[os.path.join(cfib_path, 'cfib.c')],\n",
    "              include_dirs=[np.get_include(), cfib_path])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can write our wrapper code.\n",
    "\n",
    "We make use of the `%%cython_pyximport` magic function, which is able to read the Cython code in the cell, compile it using our `pyxbld` file, then import its functions as though we had just written a Python function in a normal code cell.\n",
    "\n",
    "To write the wrapper, you first write the function signature as specified by the C library. Then, you create a *wrapper function* that makes use of the implementation and returns the result. This wrapper function becomes the interface to the compiled code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython_pyximport fib\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "\n",
    "# here we just replicate the function signature from the header\n",
    "cdef extern from \"cfib.h\":\n",
    "    void cfib(int n, double *x)\n",
    "    \n",
    "# here is the \"wrapper\"\n",
    "def fib(int n):\n",
    "    # preallocate an array of doubles\n",
    "    cdef np.ndarray[np.float64_t, ndim=1, mode='c'] x\n",
    "    x = np.empty(n, dtype=np.float64)\n",
    "    \n",
    "    # call the c function and return the result\n",
    "    cfib(n, <double *> x.data)\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can call the function `fib` to check that everything worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fib(15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using `codegen()` to Generate C Code\n",
    "\n",
    "In the previous notebook focusing on generating C code using SymPy's codegen capabilities, we saw how to form a Jacobian from a system of ODEs, then output C code that can numerically evaluate the right hand side of the ODEs themselves as well as the Jacobian, given the current state vector. Here we'll do that again so that we can wrap it in Cython for use in `scipy.integrate.odeint`.\n",
    "\n",
    "Start by importing the system of ODEs and compute the Jacobian analytically."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy2017codegen.chem import load_large_ode\n",
    "\n",
    "rhs_of_odes, states = load_large_ode()\n",
    "jac_of_odes = rhs_of_odes.jacobian(states)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we'll use `sympy.utilities.codegen.codegen` to output C source and header files which can compute the right hand side (RHS) of the ODEs and the Jacobian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy.utilities.codegen import codegen\n",
    "\n",
    "# form equations so we can put the OutputArgument in the argument_sequence\n",
    "dY = sm.MatrixSymbol('dY', len(states), 1)\n",
    "dY_eq = sm.Eq(dY, rhs_of_odes)\n",
    "t = sm.symbols('t')\n",
    "\n",
    "codegen([('c_odes', dY_eq)], 'c', argument_sequence=[*list(states), t, dY], to_files=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's take a look at the generated C code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cat c_odes.c"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Wrapping the Generated Code with Cython\n",
    "\n",
    "Now we want to wrap the function that was generated `c_odes` with a Cython function. Let's output just the header `c_odes.h` to see just the function signature."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%cat c_odes.h"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So this function signature is what we need to wrap using Cython. Just as with the Fibonacci sequence example, we start by specifying how the extension module will be built."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "render_pyxbld('cy_odes',\n",
    "              sources=['c_odes.c'],\n",
    "              include_dirs=[np.get_include()])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: Fill in the wrapper implementation below."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%%cython_pyximport cy_odes\n",
    "import numpy as np\n",
    "cimport numpy as np\n",
    "\n",
    "cdef extern from \"c_odes.h\":\n",
    "    void c_odes(double y0, double y1, double y2, double y3, double y4,\n",
    "                double y5, double y6, double y7, double y8, double y9,\n",
    "                double y10, double y11, double y12, double y13, double t,\n",
    "                double *dY)\n",
    "    \n",
    "def cy_odes(y, t):\n",
    "    cdef np.ndarray[np.float64_t, ndim=1, mode='c'] dY\n",
    "    dY = np.empty(y.size, dtype=np.float64)\n",
    "    c_odes(y[0], y[1], y[2], y[3], y[4], y[5], y[6], y[7], y[8], y[9],\n",
    "           y[10], y[11], y[12], y[13], t, <double *> dY.data)\n",
    "    return dY"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now use `odeint` to integrate the equations and plot the results to check that it worked."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.integrate import odeint\n",
    "from scipy2017codegen.chem import watrad_init, watrad_plot\n",
    "\n",
    "y_init, t_vals = watrad_init()\n",
    "y_vals = odeint(cy_odes, y_init, t_vals)\n",
    "watrad_plot(t_vals, y_vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Generating and Compiling a C Extension Module Automatically\n",
    "\n",
    "Here we'll use SymPy's `autowrap` function to automatically generate the C code for an expression, generate a Cython wrapper, compile the wrapper into an extension module, and make the function available as a callable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy.utilities.autowrap import autowrap\n",
    "\n",
    "# need to specify input args so they are in correct order\n",
    "args = list(states)\n",
    "args.append(t)\n",
    "\n",
    "f_rhs = autowrap(rhs_of_odes, args=args, backend='cython', tempdir='./autowraptmp')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`autowrap` essentially does exactly what we just did manually: generate C code to evaluate the right hand side of the system of ODEs, write a Cython wrapper, compile the C code and wrapper into an extension module, and import it for use in Python."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise**: just like the `codegen`, `autowrap` generates a C function that takes each of the state variables as separate arguments, as seen in the generated Cython wrapper file seen above. In order to use the callable returned by `autowrap` within `odeint`, write a Python wrapper that takes the state vector `y` and the current time `t` and uses `f_rhs` to compute the derivatives."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def rhs_eval(y, t):\n",
    "    yout = f_rhs(*y, t)\n",
    "    return yout[:, 0]\n",
    "\n",
    "y_vals = odeint(rhs_eval, y_init, t_vals)\n",
    "watrad_plot(t_vals, y_vals)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Using a Custom Library\n",
    "\n",
    "In our set of ODEs, there are quite a few instances of $y_i^2$, where $y_i$ is one of the 14 state variables."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from sympy.utilities.codegen import C99CodeGen\n",
    "from sympy.printing.ccode import C99CodePrinter\n",
    "\n",
    "class CustomPrinter(C99CodePrinter):\n",
    "    def _print_Pow(self, expr):\n",
    "        return 'fastpow({}, {})'.format(expr.base, expr.exp)\n",
    "\n",
    "printer = CustomPrinter()\n",
    "gen = C99CodeGen(printer=printer)\n",
    "gen.preprocessor_statements.append('#include \"fastpow.h\"')\n",
    "\n",
    "import scipy2017codegen\n",
    "import os\n",
    "fastapprox_dir = os.path.join(os.path.dirname(scipy2017codegen.__file__), 'fastapprox')\n",
    "\n",
    "f_rhs = autowrap(rhs_of_odes, args=args, code_gen=gen, backend='cython',\n",
    "                 tempdir='./autowraptmp2', include_dirs=[fastapprox_dir])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "y_vals, info = odeint(rhs_eval, y_init, t_vals, full_output=True)\n",
    "watrad_plot(t_vals, y_vals)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "info"
   ]
  }
 ],
 "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.1"
  }
 },
 "nbformat": 4,
{'fraction_non_alphanumeric': 0.15864770501575406, 'fraction_numerical': 0.007437054699253456, 'mean_word_length': 2.807413811736734, 'pattern_counts': {'":': 398, '<': 7, '<?xml version=': 0, '>': 7, 'https://': 4, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6935183', 'n_tokens_mistral': 12297, 'n_tokens_neox': 12018, 'n_words': 3790}