sci_libs.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Other Scientific Libraries"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Other Scientific Libraries](#Other-Scientific-Libraries)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Cython](#Cython)  \n",
    "  - [Joblib](#Joblib)  \n",
    "  - [Other Options](#Other-Options)  \n",
    "  - [Exercises](#Exercises)  \n",
    "  - [Solutions](#Solutions)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In this lecture we review some other scientific libraries that are useful for\n",
    "economic research and analysis\n",
    "\n",
    "We have, however, already picked most of the low hanging fruit in terms of\n",
    "economic research\n",
    "\n",
    "Hence you should feel free to skip this lecture on first pass"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cython\n",
    "\n",
    "\n",
    "<a id='index-1'></a>\n",
    "Like [Numba](https://lectures.quantecon.org/py/numba.html),  [Cython](http://cython.org/) provides an approach to generating fast compiled code that can be used from Python\n",
    "\n",
    "As was the case with Numba, a key problem is the fact that Python is dynamically typed\n",
    "\n",
    "As you’ll recall, Numba solves this problem (where possible) by inferring type\n",
    "\n",
    "Cython’s approach is different — programmers add type definitions directly to their “Python” code\n",
    "\n",
    "As such, the Cython language can be thought of as Python with type definitions\n",
    "\n",
    "In addition to a language specification, Cython is also a language translator, transforming Cython code into optimized C and C++ code\n",
    "\n",
    "Cython also takes care of building language extensions — the wrapper code that interfaces between the resulting compiled code and Python\n",
    "\n",
    "**Important Note:**\n",
    "\n",
    "In what follows code is executed in a Jupyter notebook\n",
    "\n",
    "This is to take advantage of a Cython [cell magic](http://ipython.readthedocs.org/en/stable/interactive/magics.html#cell-magics) that makes Cython particularly easy to use\n",
    "\n",
    "Some modifications are required to run the code outside a notebook\n",
    "\n",
    "- See the book [Cython](http://shop.oreilly.com/product/0636920033431.do) by Kurt Smith or [the online documentation](http://cython.org/)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A First Example\n",
    "\n",
    "Let’s start with a rather artificial example\n",
    "\n",
    "Suppose that we want to compute the sum $ \\sum_{i=0}^n \\alpha^i $ for given $ \\alpha, n $\n",
    "\n",
    "Suppose further that we’ve forgotten the basic formula\n",
    "\n",
    "$$\n",
    "\\sum_{i=0}^n \\alpha^i = \\frac{1 - \\alpha^{n+1}}{1 - \\alpha}\n",
    "$$\n",
    "\n",
    "for a geometric progression and hence have resolved to rely on a loop"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Python vs C\n",
    "\n",
    "Here’s a pure Python function that does the job"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def geo_prog(alpha, n):\n",
    "    current = 1.0\n",
    "    sum = current\n",
    "    for i in range(n):\n",
    "        current = current * alpha\n",
    "        sum = sum + current\n",
    "    return sum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This works fine but for large $ n $ it is slow\n",
    "\n",
    "Here’s a C function that will do the same thing"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "hide-output": false
   },
   "source": [
    "```c\n",
    "double geo_prog(double alpha, int n) {\n",
    "    double current = 1.0;\n",
    "    double sum = current;\n",
    "    int i;\n",
    "    for (i = 1; i <= n; i++) {\n",
    "        current = current * alpha;\n",
    "        sum = sum + current;\n",
    "    }\n",
    "    return sum;\n",
    "}\n",
    "```\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you’re not familiar with C, the main thing you should take notice of is the\n",
    "type definitions\n",
    "\n",
    "- `int` means integer  \n",
    "- `double` means double precision floating point number  \n",
    "- the `double` in `double geo_prog(...` indicates that the function will\n",
    "  return a double  \n",
    "\n",
    "\n",
    "Not surprisingly, the C code is faster than the Python code"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### A Cython Implementation\n",
    "\n",
    "Cython implementations look like a convex combination of Python and C\n",
    "\n",
    "We’re going to run our Cython code in the Jupyter notebook, so we’ll start by\n",
    "loading the Cython extension in a notebook cell"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%load_ext Cython"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the next cell, we execute the following"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "def geo_prog_cython(double alpha, int n):\n",
    "    cdef double current = 1.0\n",
    "    cdef double sum = current\n",
    "    cdef int i\n",
    "    for i in range(n):\n",
    "        current = current * alpha\n",
    "        sum = sum + current\n",
    "    return sum"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here `cdef` is a Cython keyword indicating a variable declaration, and is followed by a type\n",
    "\n",
    "The `%%cython` line at the top is not actually Cython code — it’s a Jupyter cell magic indicating the start of Cython code\n",
    "\n",
    "After executing the cell, you can now call the function `geo_prog_cython` from within Python\n",
    "\n",
    "What you are in fact calling is compiled C code with a Python call interface"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TOC: Elapsed: 0:00:0.11\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.11625146865844727"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import quantecon as qe\n",
    "qe.util.tic()\n",
    "geo_prog(0.99, int(10**6))\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TOC: Elapsed: 0:00:0.05\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.052478790283203125"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qe.util.tic()\n",
    "geo_prog_cython(0.99, int(10**6))\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Example 2: Cython with NumPy Arrays\n",
    "\n",
    "Let’s go back to the first problem that we worked with: generating the iterates of the quadratic map\n",
    "\n",
    "$$\n",
    "x_{t+1} = 4 x_t (1 - x_t)\n",
    "$$\n",
    "\n",
    "The problem of computing iterates and returning a time series requires us to work with arrays\n",
    "\n",
    "The natural array type to work with is NumPy arrays\n",
    "\n",
    "Here’s a Cython implementation that initializes, populates and returns a NumPy\n",
    "array"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "\n",
    "def qm_cython_first_pass(double x0, int n):\n",
    "    cdef int t\n",
    "    x = np.zeros(n+1, float)\n",
    "    x[0] = x0\n",
    "    for t in range(n):\n",
    "        x[t+1] = 4.0 * x[t] * (1 - x[t])\n",
    "    return np.asarray(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you run this code and time it, you will see that it’s performance is disappointing — nothing like the speed gain we got from Numba"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TOC: Elapsed: 0:00:0.03\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.03825521469116211"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qe.util.tic()\n",
    "qm_cython_first_pass(0.1, int(10**5))\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example was also computed in the [Numba lecture](https://lectures.quantecon.org/py/numba.html#qm-numba-result), and you can see Numba is around 90 times faster\n",
    "\n",
    "The reason is that working with NumPy arrays incurs substantial Python overheads\n",
    "\n",
    "We can do better by using Cython’s [typed memoryviews](http://docs.cython.org/src/userguide/memoryviews.html), which provide more direct access to arrays in memory\n",
    "\n",
    "When using them, the first step is to create a NumPy array\n",
    "\n",
    "Next, we declare a memoryview and bind it to the NumPy array\n",
    "\n",
    "Here’s an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "from numpy cimport float_t\n",
    "\n",
    "def qm_cython(double x0, int n):\n",
    "    cdef int t\n",
    "    x_np_array = np.zeros(n+1, dtype=float)\n",
    "    cdef float_t [:] x = x_np_array\n",
    "    x[0] = x0\n",
    "    for t in range(n):\n",
    "        x[t+1] = 4.0 * x[t] * (1 - x[t])\n",
    "    return np.asarray(x)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here\n",
    "\n",
    "- `cimport` pulls in some compile-time information from NumPy  \n",
    "- `cdef float_t [:] x = x_np_array` creates a memoryview on the NumPy array `x_np_array`  \n",
    "- the return statement uses `np.asarray(x)` to convert the memoryview back to a NumPy array  \n",
    "\n",
    "\n",
    "Let’s time it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TOC: Elapsed: 0:00:0.00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0008511543273925781"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qe.util.tic()\n",
    "qm_cython(0.1, int(10**5))\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is fast, although still slightly slower than `qm_numba`"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary\n",
    "\n",
    "Cython requires more expertise than Numba, and is a little more fiddly in terms of getting good performance\n",
    "\n",
    "In fact, it’s surprising how difficult it is to beat the speed improvements provided by Numba\n",
    "\n",
    "Nonetheless,\n",
    "\n",
    "- Cython is a very mature, stable and widely used tool  \n",
    "- Cython can be more useful than Numba when working with larger, more sophisticated applications  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Joblib\n",
    "\n",
    "[Joblib](https://pythonhosted.org/joblib/) is a popular Python library for\n",
    "caching and parallelization\n",
    "\n",
    "To install it, start Jupyter and type"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: joblib in /home/quantecon/anaconda3/lib/python3.7/site-packages (0.13.0)\r\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\u001b[33mYou are using pip version 18.1, however version 19.0.1 is available.\r\n",
      "You should consider upgrading via the 'pip install --upgrade pip' command.\u001b[0m\r\n"
     ]
    }
   ],
   "source": [
    "!pip install joblib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "from within a notebook\n",
    "\n",
    "Here we review just the basics"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Caching\n",
    "\n",
    "Perhaps, like us, you sometimes run a long computation that simulates a model\n",
    "at a given set of parameters — to generate a figure, say, or a table\n",
    "\n",
    "20 minutes later you realize that you want to tweak the figure and now you have to\n",
    "do it all again\n",
    "\n",
    "What caching will do is automatically store results at each parameterization\n",
    "\n",
    "With Joblib, results are compressed and stored on file, and automatically served\n",
    "back up to you when you repeat the calculation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### An Example\n",
    "\n",
    "Let’s look at a toy example, related to the quadratic map model discussed [above](https://lectures.quantecon.org/py/numba.html#quad-map-eg)\n",
    "\n",
    "Let’s say we want to generate a long trajectory from a certain initial\n",
    "condition $ x_0 $ and see what fraction of the sample is below 0.1\n",
    "\n",
    "(We’ll omit JIT compilation or other speed ups for simplicity)\n",
    "\n",
    "Here’s our code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/quantecon/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:3: DeprecationWarning: The 'cachedir' parameter has been deprecated in version 0.12 and will be removed in version 0.14.\n",
      "You provided \"cachedir='./joblib_cache'\", use \"location='./joblib_cache'\" instead.\n",
      "  This is separate from the ipykernel package so we can avoid doing imports until\n"
     ]
    }
   ],
   "source": [
    "from joblib import Memory\n",
    "\n",
    "memory = Memory(cachedir='./joblib_cache')\n",
    "\n",
    "@memory.cache\n",
    "def qm(x0, n):\n",
    "    x = np.empty(n+1)\n",
    "    x[0] = x0\n",
    "    for t in range(n):\n",
    "        x[t+1] = 4 * x[t] * (1 - x[t])\n",
    "    return np.mean(x < 0.1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are using [joblib](http://pythonhosted.org/joblib/) to cache the result of calling qm at a given set of parameters\n",
    "\n",
    "With the argument cachedir=’./joblib_cache’, any call to this function results in both the input values and output values being stored a subdirectory joblib_cache of the present working directory\n",
    "\n",
    "(In UNIX shells, . refers to the present working directory)\n",
    "\n",
    "The first time we call the function with a given set of parameters we see some\n",
    "extra output that notes information being cached"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "________________________________________________________________________________\n",
      "[Memory] Calling __main__--home-quantecon-repos-collab-quantecon.build.lectures-_build_jupyter-py-__ipython-input__.qm...\n",
      "qm(0.2, 10000000)\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "_______________________________________________________________qm - 6.8s, 0.1min\n",
      "TOC: Elapsed: 0:00:6.81\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "6.8197596073150635"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qe.util.tic()\n",
    "n = int(1e7)\n",
    "qm(0.2, n)\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next time we call the function with the same set of parameters, the result is returned almost instantaneously"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TOC: Elapsed: 0:00:0.00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0008389949798583984"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qe.util.tic()\n",
    "n = int(1e7)\n",
    "qm(0.2, n)\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Other Options\n",
    "\n",
    "There are in fact many other approaches to speeding up your Python code\n",
    "\n",
    "One is interfacing with Fortran\n",
    "\n",
    "\n",
    "<a id='index-2'></a>\n",
    "If you are comfortable writing Fortran you will find it very easy to create\n",
    "extension modules from Fortran code using [F2Py](https://docs.scipy.org/doc/numpy/f2py/)\n",
    "\n",
    "F2Py is a Fortran-to-Python interface generator that is particularly simple to\n",
    "use\n",
    "\n",
    "Robert Johansson provides a [very nice introduction](http://nbviewer.jupyter.org/github/jrjohansson/scientific-python-lectures/blob/master/Lecture-6A-Fortran-and-C.ipynb) to F2Py, among other things\n",
    "\n",
    "Recently, [a Jupyter cell magic for Fortran](http://nbviewer.jupyter.org/github/mgaitan/fortran_magic/blob/master/documentation.ipynb) has been developed — you might want to give it a try"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercises\n",
    "\n",
    "\n",
    "<a id='speed-ex1'></a>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "Later we’ll learn all about [finite state Markov chains](https://lectures.quantecon.org/py/finite_markov.html)\n",
    "\n",
    "For now, let’s just concentrate on simulating a very simple example of such a chain\n",
    "\n",
    "Suppose that the volatility of returns on an asset can be in one of two regimes — high or low\n",
    "\n",
    "The transition probabilities across states are as follows\n",
    "\n",
    "<img src=\"https://s3-ap-southeast-2.amazonaws.com/lectures.quantecon.org/py/_static/figures/nfs_ex1.png\" style=\"\">\n",
    "\n",
    "  \n",
    "For example, let the period length be one month, and suppose the current state is high\n",
    "\n",
    "We see from the graph that the state next month will be\n",
    "\n",
    "- high with probability 0.8  \n",
    "- low with probability 0.2  \n",
    "\n",
    "\n",
    "Your task is to simulate a sequence of monthly volatility states according to this rule\n",
    "\n",
    "Set the length of the sequence to `n = 100000` and start in the high state\n",
    "\n",
    "Implement a pure Python version, a Numba version and a Cython version, and compare speeds\n",
    "\n",
    "To test your code, evaluate the fraction of time that the chain spends in the low state\n",
    "\n",
    "If your code is correct, it should be about 2/3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solutions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise 1\n",
    "\n",
    "We let\n",
    "\n",
    "- 0 represent “low”  \n",
    "- 1 represent “high”  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "p, q = 0.1, 0.2  # Prob of leaving low and high state respectively"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here’s a pure Python version of the function"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "def compute_series(n):\n",
    "    x = np.empty(n, dtype=int)\n",
    "    x[0] = 1  # Start in state 1\n",
    "    U = np.random.uniform(0, 1, size=n)\n",
    "    for t in range(1, n):\n",
    "        current_x = x[t-1]\n",
    "        if current_x == 0:\n",
    "            x[t] = U[t] < p\n",
    "        else:\n",
    "            x[t] = U[t] > q\n",
    "    return x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s run this code and check that the fraction of time spent in the low\n",
    "state is about 0.666"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.66879\n"
     ]
    }
   ],
   "source": [
    "n = 100000\n",
    "x = compute_series(n)\n",
    "print(np.mean(x == 0))  # Fraction of time x is in state 0"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now let’s time it"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TOC: Elapsed: 0:00:0.07\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0782933235168457"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qe.util.tic()\n",
    "compute_series(n)\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next let’s implement a Numba version, which is easy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "from numba import jit\n",
    "\n",
    "compute_series_numba = jit(compute_series)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s check we still get the right numbers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.67038\n"
     ]
    }
   ],
   "source": [
    "x = compute_series_numba(n)\n",
    "print(np.mean(x == 0))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let’s see the time"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TOC: Elapsed: 0:00:0.00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.0016720294952392578"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qe.util.tic()\n",
    "compute_series_numba(n)\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is a nice speed improvement for one line of code\n",
    "\n",
    "Now let’s implement a Cython version"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "hide-output": false
   },
   "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": 22,
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "%%cython\n",
    "import numpy as np\n",
    "from numpy cimport int_t, float_t\n",
    "\n",
    "def compute_series_cy(int n):\n",
    "    # == Create NumPy arrays first == #\n",
    "    x_np = np.empty(n, dtype=int)\n",
    "    U_np = np.random.uniform(0, 1, size=n)\n",
    "    # == Now create memoryviews of the arrays == #\n",
    "    cdef int_t [:] x = x_np\n",
    "    cdef float_t [:] U = U_np\n",
    "    # == Other variable declarations == #\n",
    "    cdef float p = 0.1\n",
    "    cdef float q = 0.2\n",
    "    cdef int t\n",
    "    # == Main loop == #\n",
    "    x[0] = 1\n",
    "    for t in range(1, n):\n",
    "        current_x = x[t-1]\n",
    "        if current_x == 0:\n",
    "            x[t] = U[t] < p\n",
    "        else:\n",
    "            x[t] = U[t] > q\n",
    "    return np.asarray(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([1, 1, 1, 0, 0, 0, 0, 0, 0, 0])"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "compute_series_cy(10)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.66692\n"
     ]
    }
   ],
   "source": [
    "x = compute_series_cy(n)\n",
    "print(np.mean(x == 0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "hide-output": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "TOC: Elapsed: 0:00:0.00\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.003400087356567383"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "qe.util.tic()\n",
    "compute_series_cy(n)\n",
    "qe.util.toc()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The Cython implementation is fast, but not as fast as Numba"
   ]
  }
 ],
 "metadata": {
  "filename": "sci_libs.rst",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.0"
  },
  "title": "Other Scientific Libraries"
 },
 "nbformat": 4,
 "nbformat_minor": 2
}