{
 "metadata": {
  "name": "",
  "signature": "sha256:78aeab2a700bfc36e29a4a74c032f90fbeb069bb9dc791e12ab75f712750d589"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "#Cython"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Source of that code: Vogel http://nbviewer.ipython.org/github/carljv/cython_testing/blob/master/cython_linalg.ipynb"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from sympy.interactive.printing import init_printing\n",
      "init_printing(use_latex=True)\n",
      "from sympy import Rational, pi\n",
      "Rational(3,2)*pi"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "$$\\frac{3 \\pi}{2}$$"
       ],
       "metadata": {},
       "output_type": "pyout",
       "png": "iVBORw0KGgoAAAANSUhEUgAAABcAAAAqBAMAAACq4N3gAAAAMFBMVEX///8AAAAAAAAAAAAAAAAA\nAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAv3aB7AAAAD3RSTlMAIom7VJlmdt1E7xDN\nMqsI8sYEAAAA40lEQVQoFWNgEFIyYYABxgAG/wQYh/0LA38DjMO5kmH/ARgHSCOUATlHGRjq////\n/xXI5NUIYGDSdRcrFAAr1rpwhUGR/QKYzcC1moHBgG8DkMMowMD8hYH3AQdIgv83A/M3BtYEfpAM\nlwEDx0eG/Rv4JwA5bAkM/gUMFQxcCkAOw9RQSwYGMwY+MAckQAcA9BUMfKCOdYyhxhfgJu1jYPoJ\n57xmYFgF56xgYHgPV3d+AxIHqKR/A1wdA+cnBJuBwwCJE4zEZn+AxHEFhT8U8AYwsMM5EkqqL2AS\nDOf///8B4gAAIqI8anv4CrQAAAAASUVORK5CYII=\n",
       "prompt_number": 1,
       "text": [
        "3\u22c5\u03c0\n",
        "\u2500\u2500\u2500\n",
        " 2 "
       ]
      }
     ],
     "prompt_number": 1
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%latex\n",
      "\\begin{align}\n",
      "a_{00}\\,x_0 + a_{01}\\,x_1 + \\ldots + a_{0,n-1}\\,x_{n-1}   &= b_0 \\\\\\\n",
      "a_{10}\\,x_0  + a_{11}\\,x_1 + \\ldots + a_{1,n-1}\\,x_{n-1}   &= b_1 \\\\\\\n",
      "\\vdots &  \\\\\\\n",
      "a_{n-1,0}\\,x_0 + a_{n-1,1}\\,x_1 + \\ldots + a_{n-1,n-1}\\,x_{n-1} &= b_{n-1}\\ .\n",
      "\\end{align}"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "latex": [
        "\\begin{align}\n",
        "a_{00}\\,x_0 + a_{01}\\,x_1 + \\ldots + a_{0,n-1}\\,x_{n-1}   &= b_0 \\\\\\\n",
        "a_{10}\\,x_0  + a_{11}\\,x_1 + \\ldots + a_{1,n-1}\\,x_{n-1}   &= b_1 \\\\\\\n",
        "\\vdots &  \\\\\\\n",
        "a_{n-1,0}\\,x_0 + a_{n-1,1}\\,x_1 + \\ldots + a_{n-1,n-1}\\,x_{n-1} &= b_{n-1}\\ .\n",
        "\\end{align}"
       ],
       "metadata": {},
       "output_type": "display_data",
       "text": [
        "<IPython.core.display.Latex at 0x1047ac190>"
       ]
      }
     ],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext cythonmagic\n",
      "import numpy as np\n",
      "import scipy.linalg as la"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "cimport cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "\n",
      "def cython_matprod(np.ndarray[double, ndim = 2] A, B):\n",
      "    '''\n",
      "    Matrix-by-vector or matrix-by-matrix multiplication.\n",
      "\n",
      "    The arguments are dispatched to one of two functions\n",
      "    depending on whether B is a vector or a matrix.\n",
      "    '''\n",
      "    if B.ndim == 1:\n",
      "       # B is a vector\n",
      "       return matvecprod(A, B)\n",
      "    else:\n",
      "        # B is a matrix\n",
      "        return matmatprod(A, B)\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "cdef np.ndarray[double, ndim=2] matmatprod(\n",
      "    np.ndarray[double, ndim=2] A,\n",
      "    np.ndarray[double, ndim=2] B):\n",
      "    '''\n",
      "    Matrix-matrix multiplication.\n",
      "    '''\n",
      "    cdef: \n",
      "        int i, j, k\n",
      "        int A_n = A.shape[0]\n",
      "        int A_m = A.shape[1]\n",
      "        int B_n = B.shape[0]\n",
      "        int B_m = B.shape[1]\n",
      "        np.ndarray[double, ndim=2] C\n",
      "    \n",
      "    # Are matrices conformable?\n",
      "    assert A_m == B_n, \\\n",
      "        'Non-conformable shapes.'\n",
      "    \n",
      "    # Initialize the results matrix.\n",
      "    C = np.zeros((A_n, B_m))\n",
      "    for i in xrange(A_n):\n",
      "        for j in xrange(B_m):\n",
      "            for k in xrange(A_m):\n",
      "                C[i, j] += A[i, k] * B[k, j]\n",
      "    return C\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "cdef np.ndarray[double, ndim=1] matvecprod(\n",
      "    np.ndarray[double, ndim=2] A,\n",
      "    np.ndarray[double, ndim=1] b):\n",
      "    '''\n",
      "    Matrix-vector multiplication.\n",
      "    '''\n",
      "    cdef: \n",
      "        Py_ssize_t i, j, k\n",
      "        Py_ssize_t A_n = A.shape[0]\n",
      "        Py_ssize_t A_m = A.shape[1]\n",
      "        Py_ssize_t b_n = b.shape[0]\n",
      "        np.ndarray[double, ndim=1] c\n",
      "    \n",
      "    # Are matrices conformable?\n",
      "    assert A_m == b_n, \\\n",
      "        'Non-conformable shapes.'\n",
      "    \n",
      "    # Initialize the results matrix.\n",
      "    c = np.zeros(A_n)\n",
      "    for i in xrange(A_n):\n",
      "            for k in xrange(b_n):\n",
      "                c[i] += A[i, k] * b[k]\n",
      "    return c"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def python_matmatprod(A, B):\n",
      "    '''\n",
      "    Matrix-matrix multiplication\n",
      "    '''\n",
      "    A_n, A_m = A.shape\n",
      "    B_n, B_m = B.shape\n",
      "    assert A_m == B_n, \"Non-conformable shapes.\"\n",
      "    C = np.zeros((A_n, B_m))\n",
      "    for i in xrange(A_n):\n",
      "        for j in xrange(B_m):\n",
      "            for k in xrange(A_m):\n",
      "                C[i, j] += A[i, k] * B[k, j]\n",
      "    return C"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 12
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "# A is 2x3\n",
      "A = np.array([[2.0, 0.25, -1.0], \n",
      "              [3.0, 0.0 ,  5.0]])\n",
      "# B is 3x2\n",
      "B = np.array([[-3.0,  0.5], \n",
      "              [ 2.0,  1.5], \n",
      "              [ 4.0, -4.0]])\n",
      "# C is 2x2\n",
      "C = np.array([[1.0,  1.5], \n",
      "              [2.5, -1.0]])\n",
      "# b is 3x1 (a vector)\n",
      "b = np.array([1.0, -2.0, 0.5])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print 'Cython:'\n",
      "print '-------'\n",
      "print \"A x B =\\n\", cython_matprod(A, B), \"\\n\"\n",
      "print \"A x b =\\n\", cython_matprod(A, b), \"\\n\"\n",
      "print 'Numpy dot:'\n",
      "print '----------'\n",
      "print \"A x B =\\n\", np.dot(A, B), \"\\n\"\n",
      "print \"A x b =\\n\", np.dot(A, b), \"\\n\"\n",
      "print 'Python loops:'\n",
      "print '-------------'\n",
      "print \"A x B =\\n\", python_matmatprod(A, B), \"\\n\""
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Cython:\n",
        "-------\n",
        "A x B =\n",
        "[[ -9.5     5.375]\n",
        " [ 11.    -18.5  ]] \n",
        "\n",
        "A x b =\n",
        "[ 1.   5.5] \n",
        "\n",
        "Numpy dot:\n",
        "----------\n",
        "A x B =\n",
        "[[ -9.5     5.375]\n",
        " [ 11.    -18.5  ]] \n",
        "\n",
        "A x b =\n",
        "[ 1.   5.5] \n",
        "\n",
        "Python loops:\n",
        "-------------\n",
        "A x B =\n",
        "[[ -9.5     5.375]\n",
        " [ 11.    -18.5  ]] \n",
        "\n"
       ]
      }
     ],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit np.dot(A, B)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000000 loops, best of 3: 951 ns per loop\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit cython_matprod(A, B)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 4.27 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 15
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit python_matmatprod(A, B)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100000 loops, best of 3: 12.7 \u00b5s per loop\n"
       ]
      }
     ],
     "prompt_number": 16
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}