mamba mutation cython.ipynb

{
 "metadata": {
  "name": "mamba mutation cython"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Mutation function optimization\n",
      "This is an attempt to optimize the mutation function of `mamba`.\n",
      "\n",
      "First I need to change dir:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "cd d:\\workspace\\mamba"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "d:\\workspace\\mamba\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now use `simulation` to create a `population` and some `genomes`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import simulation"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "p,g = simulation.run(0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Starting simulation with  0 ticks\n",
        "Tick"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0\n",
        "Simulation finished, 0 ticks, time elapsed 0.213297071155 seconds\n"
       ]
      }
     ],
     "prompt_number": 3
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "load the models where the python and cython functions are written:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import cython_load\n",
      "import model\n",
      "import model_c\n",
      "reload(model);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 28
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "mutation_numpy = model.mutation_explicit_genomes\n",
      "mutation_numpy2 = model.mutation_explicit_genomes2\n",
      "# cython functions removed"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Time the functions:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit mutation_numpy(p, g, array([0.003]*len(p)), g.shape[1], g[0,:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 3.31 s per loop\n"
       ]
      }
     ],
     "prompt_number": 11
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit mutation_numpy2(p, g, array([0.003]*len(p)), g.shape[1], g[0,:]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 1.76 s per loop\n"
       ]
      }
     ],
     "prompt_number": 35
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit mutation_cython(p, g, array([0.003]*len(p)), g.shape[1], g[0,:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 2.9 s per loop\n"
       ]
      }
     ],
     "prompt_number": 7
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit mutation_cython_cdef(p, g, array([0.003]*len(p)), g.shape[1], g[0,:])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1 loops, best of 3: 2.98 s per loop\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So the second `numpy` function wins - `cython` has no advantage here. \n",
      "\n",
      "The difference between the winner and the other `numpy` function is that `find_row` is at the end, at the `update` stage, so it gets called less times. \n",
      "\n",
      "Let's verify this using profiling:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%prun mutation_numpy(p, g, array([0.003]*len(p)), g.shape[1], g[0,:]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 51
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "13515 function calls in 2.314 seconds\n",
      "\n",
      "ncalls  tottime  percall  cumtime  percall filename:lineno(function)\n",
      "\n",
      "2598    2.180    0.001    2.180    0.001 {model_c.find_row}"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%prun mutation_numpy2(p, g, array([0.003]*len(p)), g.shape[1], g[0,:]);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 52
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "12562 function calls in 1.229 seconds\n",
      "\n",
      "ncalls  tottime  percall  cumtime  percall filename:lineno(function)\n",
      "\n",
      "1295    1.123    0.001    1.123    0.001 {model_c.find_row}"
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So I got an improvement of roughly half the time because I prevented roughly half of the `find_row` calls.\n",
      "\n",
      "## Optimize `find_row`\n",
      "The next step is to optimize the `cython` code of `find_row`.\n",
      "\n",
      "The naive `numpy` implelemnation is:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def find_row_naive(matrix, target):\n",
      "\tfor i, row in enumerate(matrix):\n",
      "\t\tfound = True\n",
      "\t\tfor j in range(target.shape[0]):\n",
      "\t\t\tif target[j] != row[j]:\n",
      "\t\t\t\tfound = False\n",
      "\t\t\t\tbreak\n",
      "\t\tif found:\n",
      "\t\t\treturn i\n",
      "\treturn -1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 41
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "But I'll use a faster `cython` implementation. First I load the `cython` magic extensions and the genome that is in the middle of the `genomes` array (average case):"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%load_ext cythonmagic\n",
      "ge = g[g.shape[0]/2,:]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 13
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This is the function from `model_c.pyx` used above:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "cimport cython\n",
      "from cpython cimport bool\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False) # TODO int -> DTYPE_...\n",
      "def find_row(np.ndarray[int, ndim=2] matrix, np.ndarray[int, ndim=1] target):\n",
      "\tcdef np.ndarray[int, ndim=1] row \n",
      "\tcdef Py_ssize_t i, j\n",
      "\tcdef bool found\n",
      "\tfor i, row in enumerate(matrix):\n",
      "\t\tfound = True\n",
      "\t\tfor j in range(target.shape[0]):\n",
      "\t\t\tif target[j] != row[j]:\n",
      "\t\t\t\tfound = False\n",
      "\t\t\t\tbreak\n",
      "\t\tif found:\n",
      "\t\t\treturn i\n",
      "\treturn -1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 42
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit find_row(g, ge)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 530 us per loop\n"
       ]
      }
     ],
     "prompt_number": 18
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So the target is to **improve below 500 us**.\n",
      "\n",
      "The next function tries to use `numpy`'s array comparison to avoid the inner loop. \n",
      "\n",
      "I think it is less effective, as we compare *all* the elements, whereas in the previous implementation we stop comparing when we find a mismatch. \n",
      "\n",
      "But still worth the try:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "cimport cython\n",
      "from cpython cimport bool\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False) # TODO int -> DTYPE_...\n",
      "def find_row_np_compare(np.ndarray[int, ndim=2] matrix, np.ndarray[int, ndim=1] target):\n",
      "    cdef np.ndarray[int, ndim=1] row \n",
      "    cdef Py_ssize_t i, j\n",
      "    cdef bool found\n",
      "    for i, row in enumerate(matrix):\n",
      "        if (row==target).all():\n",
      "            return i\n",
      "    return -1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit find_row_np_compare(g, ge)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 3.5 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Not as good as the previous `find_row`.\n",
      "\n",
      "The next step is to try to use some other `numpy` features:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "cimport cython\n",
      "from cpython cimport bool\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False) # TODO int -> DTYPE_...\n",
      "def find_row_np_tostring(np.ndarray[int, ndim=2] matrix, np.ndarray[int, ndim=1] target):\n",
      "    cdef np.ndarray[int, ndim=1] row \n",
      "    cdef Py_ssize_t i, j\n",
      "    cdef bool found\n",
      "    for i, row in enumerate(matrix):\n",
      "        if row.tostring()==target.tostring():\n",
      "            return i\n",
      "    return -1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit find_row_np_tostring(g, ge)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 996 us per loop\n"
       ]
      }
     ],
     "prompt_number": 22
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Another approach is to have an array of `nums` ready, and to check on that. It is costly to create but can save time when searching:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit model_c.genomes_to_nums(g)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 103 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 23
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "nums = model_c.genomes_to_nums(g)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 31
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "cimport cython\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "def find_row_nums(np.ndarray[double, ndim=1] nums, np.ndarray[int, ndim=1] target):\n",
      "    cdef target_num\n",
      "    cdef Py_ssize_t i\n",
      "    \n",
      "    target_num =  (2. ** np.arange(target.shape[0]) * target).sum()\n",
      "    for i,n in enumerate(nums):\n",
      "        if n == target_num:\n",
      "            return i\n",
      "    return -1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 32
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now benchmark against the original function:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit find_row(g, ge)\n",
      "%timeit find_row_nums(nums, ge)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 539 us per loop\n",
        "10000 loops, best of 3: 165 us per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 27
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So this approach is **3-5-fold faster on searching but slow on inserting**.\n",
      "\n",
      "## Optimize `genome_to_num`\n",
      "\n",
      "Let's see if we can optimize inserting. The first function is the one used above, the others are new:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def genome_to_num(genome):\n",
      "\ti = np.arange(genome.shape[0])\n",
      "\treturn (2. ** i * genome).sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 29
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "cimport cython\n",
      "DTYPE_INT = np.int\n",
      "ctypedef np.int_t DTYPE_INT_t\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "def genome_to_num_cython(np.ndarray[DTYPE_INT_t, ndim=1] genome):\n",
      "    cdef np.ndarray[DTYPE_INT_t, ndim=1] i\n",
      "    i = np.arange(genome.shape[0])\n",
      "    return (2. ** i * genome).sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 30
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import numexpr\n",
      "def genome_to_num_numexpr(genome):\n",
      "\ti = np.arange(genome.shape[0])\n",
      "\treturn numexpr.evaluate(\"sum(2. ** i * genome)\").item()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 105
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "print genome_to_num(ge)\n",
      "#genome_to_num_cython(ge) # this keeps throwing errors, and I think it isn't faster anyway\n",
      "print genome_to_num_numexpr(ge)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1.7014118346e+38\n",
        "1.7014118346e+38\n"
       ]
      }
     ],
     "prompt_number": 32
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit [genome_to_num(ge) for ge in g];\n",
      "%timeit [genome_to_num_numexpr(ge) for ge in g];"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 105 ms per loop\n",
        "10 loops, best of 3: 51.3 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 33
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So here the `numexpr` approach is ~2-fold faster for the entire `genomes` table!\n",
      "\n",
      "Now to check if list comprehension is better than `for`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit [genome_to_num_numexpr(ge) for ge in g];"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 50.6 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 34
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%timeit \n",
      "a = [None]*g.shape[0]\n",
      "for i,ge in enumerate(g): \n",
      "    a[i]=genome_to_num_numexpr(ge)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 49.2 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 35
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Not much of a difference - I think this is because the *list comprehension* allocates memory, just like the `for` loop I used (in contrast to an approach that uses `lst.append(item)`).\n",
      "\n",
      "Can I improve on **50 ms for the `genomes` table**?\n",
      "\n",
      "Another `cython`, non-`numexpr` attempt:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "import numpy.random as npr\n",
      "cimport cython\n",
      "from cpython cimport bool\n",
      "\n",
      "DTYPE_INT = np.int\n",
      "ctypedef np.int_t DTYPE_INT_t\n",
      "DTYPE_FLOAT = np.float\n",
      "ctypedef np.float_t DTYPE_FLOAT_t\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "def genomes_to_nums(np.ndarray[DTYPE_INT_t, ndim=2] genomes):\n",
      "    cdef np.ndarray[DTYPE_FLOAT_t, ndim=1] nums\n",
      "    cdef np.ndarray[DTYPE_FLOAT_t, ndim=1] i\n",
      "    cdef DTYPE_FLOAT_t num_loci_f \n",
      "    num_loci_f = np.float(genomes.shape[1])\n",
      "    nums = np.zeros(genomes.shape[0], dtype=DTYPE_FLOAT)\n",
      "    cdef Py_ssize_t j\n",
      "    i = np.arange(num_loci_f)\n",
      "    for j in range(genomes.shape[0]):\n",
      "        nums[j] = (2. ** i * genomes[j,:]).sum()\n",
      "    return nums"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit genomes_to_nums(g);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 100 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 37
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "And one with `numexpr`, if possible:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "import numpy.random as npr\n",
      "cimport cython\n",
      "import model\n",
      "reload(model)\n",
      "\n",
      "DTYPE_INT = np.int\n",
      "ctypedef np.int_t DTYPE_INT_t\n",
      "DTYPE_FLOAT = np.float\n",
      "ctypedef np.float_t DTYPE_FLOAT_t\n",
      "\n",
      "@cython.boundscheck(False)\n",
      "@cython.wraparound(False)\n",
      "def genomes_to_nums_numexpr(np.ndarray[DTYPE_INT_t, ndim=2] genomes):\n",
      "    cdef np.ndarray[DTYPE_FLOAT_t, ndim=1] nums\n",
      "    cdef np.ndarray[DTYPE_INT_t, ndim=1] gen\n",
      "    cdef np.ndarray[DTYPE_FLOAT_t, ndim=1] i\n",
      "    cdef DTYPE_FLOAT_t num_loci_f \n",
      "    num_loci_f = np.float(genomes.shape[1])\n",
      "    nums = np.zeros(genomes.shape[0], dtype=DTYPE_FLOAT)\n",
      "    cdef Py_ssize_t j\n",
      "    i = np.arange(num_loci_f)\n",
      "    for j in range(genomes.shape[0]):\n",
      "        gen = genomes[j, :]\n",
      "        nums[j] = model.genome_to_num_numexpr(gen)\n",
      "    return nums"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 38
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This keeps throwing errors, though."
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "#%timeit genomes_to_nums_numexpr(g)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 40
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "The next approach is to use the fact that `genome` is a sparse matrix:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def genome_to_num_numexpr_nonzero(genome):\n",
      "    i = genome.nonzero()[0]\n",
      "    return numexpr.evaluate(\"sum(2. ** i )\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 112
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def genome_to_num_nonzero(genome):\n",
      "    i = genome.nonzero()[0]\n",
      "    return (2. ** i).sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 103
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit [genome_to_num_numexpr(ge) for ge in g]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 49.3 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 114
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit [genome_to_num_numexpr_nonzero(ge) for ge in g]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "10 loops, best of 3: 30.1 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 115
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit [genome_to_num_nonzero(ge) for ge in g]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 13.9 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 116
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "This approach works nice, we **cut the running time by 60%** compared to the dense `numexpr` approach, and by 1/8 compared to the `cython`+`numpy` which we started from. But now the `numexpr` approach isn'y better than a simple `numpy` approach!\n",
      "\n",
      "The **sparse** approach is the definite winner, cutting the time by another factor of ~50.\n",
      "\n",
      "Now I want to check if `list comprehension` is better than allocating an array and filling it, either with `python` or with `cython`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def genomes_to_nums_compr(genomes):\n",
      "    return [genome_to_num_nonzero(ge) for ge in g]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 119
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def genomes_to_nums_alloc(genomes):\n",
      "    nums = [None]*genomes.shape[0]\n",
      "    for i,ge in enumerate(genomes):\n",
      "        nums[i] = genome_to_num_nonzero(ge)\n",
      "    return nums"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 120
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%cython\n",
      "import numpy as np\n",
      "cimport numpy as np\n",
      "import numpy.random as npr\n",
      "cimport cython\n",
      "\n",
      "DTYPE_INT = np.int\n",
      "ctypedef np.int_t DTYPE_INT_t\n",
      "DTYPE_FLOAT = np.float\n",
      "ctypedef np.float_t DTYPE_FLOAT_t\n",
      "\n",
      "def genome_to_num_nonzero(genome):\n",
      "    i = genome.nonzero()[0]\n",
      "    return (2. ** i).sum()\n",
      "\n",
      "def genomes_to_nums_cython(np.ndarray[DTYPE_INT_t, ndim=2] genomes):\n",
      "    cdef np.ndarray[DTYPE_FLOAT_t, ndim=1] nums\n",
      "    cdef Py_ssize_t i\n",
      "    cdef np.ndarray[DTYPE_INT_t, ndim=1] ge \n",
      "    \n",
      "    nums = np.empty(genomes.shape[0])\n",
      "    for i,ge in enumerate(genomes):\n",
      "        nums[i] = genome_to_num_nonzero(ge)\n",
      "    return nums"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 126
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit genomes_to_nums_compr(g)\n",
      "%timeit genomes_to_nums_alloc(g)\n",
      "%timeit genomes_to_nums_cython(g)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 13.2 ms per loop\n",
        "100 loops, best of 3: 13.6 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100 loops, best of 3: 23.2 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 127
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So here the `cython` approach isn't worth it, and *list comprehension* is as good as pre-allocating a list. I also tris using a `list` for `nums` in the `cython` version, it didn't matter."
     ]
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Sparse matrix \n",
      "\n",
      "But now that I think of it, this *sparse* approach with calculating the num might be improved if instead of calculating the nums I simply compare non-zero elements.\n",
      "So let's convert `genomes` to a sparse matrix and then search on it in various ways.\n",
      "\n",
      "Start with the naive approach:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now let's try to use the sparse matrix of `scipy`. \n",
      "\n",
      "Import `scipy` and create different kinds of sparse matrix versions of `genomes`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import scipy.sparse as spr"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 4
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gcsr = spr.csr_matrix(g)\n",
      "gcsc = spr.csc_matrix(g)\n",
      "gdok = spr.dok_matrix(g)\n",
      "glil = spr.lil_matrix(g)\n",
      "gcoo = spr.coo_matrix(g)\n",
      "gbsr = spr.bsr_matrix(g)\n",
      "gdia = spr.dia_matrix(g)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 63
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Check which sparse matrix type gives `nonzero()` faster:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit -n 100 g.nonzero();\n",
      "%timeit -n 100 gcsr.nonzero();\n",
      "%timeit -n 100 gcsc.nonzero();\n",
      "%timeit -n 100 gdok.nonzero();\n",
      "%timeit -n 100 glil.nonzero();\n",
      "%timeit -n 100 gcoo.nonzero();\n",
      "%timeit -n 100 gbsr.nonzero();\n",
      "%timeit -n 100 gdia.nonzero();"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 14.1 ms per loop\n",
        "100 loops, best of 3: 124 us per loop\n",
        "100 loops, best of 3: 125 us per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100 loops, best of 3: 7.87 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100 loops, best of 3: 1.42 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100 loops, best of 3: 28.8 us per loop\n",
        "100 loops, best of 3: 305 us per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100 loops, best of 3: 28.1 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 106
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "**`coo_matrix` is the winner** - fastest `nonzero()`. \n",
      "Let's verify on random matrix:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "r = spr.rand(g.shape[0],g.shape[1],g.sum()/float(len(g.flatten())),format='csr')\n",
      "%timeit -n 100 r.nonzero();\n",
      "r = spr.rand(g.shape[0],g.shape[1],g.sum()/float(len(g.flatten())),format='csc')\n",
      "%timeit -n 100 r.nonzero();\n",
      "r = spr.rand(g.shape[0],g.shape[1],g.sum()/float(len(g.flatten())),format='dok')\n",
      "%timeit -n 100 r.nonzero();\n",
      "r = spr.rand(g.shape[0],g.shape[1],g.sum()/float(len(g.flatten())),format='lil')\n",
      "%timeit -n 100 r.nonzero();\n",
      "r = spr.rand(g.shape[0],g.shape[1],g.sum()/float(len(g.flatten())),format='coo')\n",
      "%timeit -n 100 r.nonzero();\n",
      "r = spr.rand(g.shape[0],g.shape[1],g.sum()/float(len(g.flatten())),format='bsr')\n",
      "%timeit -n 100 r.nonzero();\n",
      "r = spr.rand(g.shape[0],g.shape[1],g.sum()/float(len(g.flatten())),format='dia')\n",
      "%timeit -n 100 r.nonzero();"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 131 us per loop\n",
        "100 loops, best of 3: 132 us per loop\n",
        "100 loops, best of 3: 7.78 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100 loops, best of 3: 1.38 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100 loops, best of 3: 32 us per loop\n",
        "100 loops, best of 3: 295 us per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "100 loops, best of 3: 34.6 ms per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 8
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Indeed, `coo_matrix` is the winner. \n",
      "\n",
      "Let's use it. Assume that `genomes` is `coo_matrix`...\n",
      "\n",
      "I didn't find a nice way to do it.\n",
      "\n",
      "I did find a way to use `nonzero()` with the `lil_matrix`:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def find_row_sparse_naive(genomes, target):    \n",
      "    #genomes_lil = spr.lil_matrix(genomes)\n",
      "    genz = target.nonzero()[0]\n",
      "    for i,row in enumerate(glil.rows):\n",
      "        x = (row==genz)\n",
      "        if x.all() and x.any():\n",
      "            return i\n",
      "    return -1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 97
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit find_row(g, ge)\n",
      "%timeit find_row_nums(nums, ge)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "1000 loops, best of 3: 979 us per loop\n",
        "1000 loops, best of 3: 273 us per loop"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n"
       ]
      }
     ],
     "prompt_number": 98
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit find_row_sparse_naive(glil, ge)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "100 loops, best of 3: 11.4 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 99
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "But it is not as fast as the `nums` approach.\n",
      "\n",
      "## Test the `genome_to_num` function"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a=g[0,:]\n",
      "genome_to_num_nonzero(a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 132,
       "text": [
        "0.0"
       ]
      }
     ],
     "prompt_number": 132
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a[0]=1\n",
      "genome_to_num_nonzero(a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 133,
       "text": [
        "1.0"
       ]
      }
     ],
     "prompt_number": 133
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a[0]=0\n",
      "a[-1]=1\n",
      "genome_to_num_nonzero(a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 134,
       "text": [
        "5.3575430359313366e+300"
       ]
      }
     ],
     "prompt_number": 134
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "2**999==genome_to_num_nonzero(a)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 136,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 136
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a[-1]=0\n",
      "a.sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 138,
       "text": [
        "0"
       ]
      }
     ],
     "prompt_number": 138
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a.fill(1)\n",
      "a.sum()"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 140,
       "text": [
        "1000"
       ]
      }
     ],
     "prompt_number": 140
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "genome_to_num_nonzero(a)==(2**1000)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "pyout",
       "prompt_number": 143,
       "text": [
        "True"
       ]
      }
     ],
     "prompt_number": 143
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "## Summary\n",
      "\n",
      "I started with 3 secs per 1 generation of simulation. Where am I now?"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "import simulation\n",
      "reload(simulation);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 21
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%timeit simulation.run(0,0);"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Starting simulation with  0 ticks\n",
        "Simulation finished,"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0 ticks, time elapsed 0.319544448837 seconds\n",
        "Starting simulation with  0 ticks\n",
        "Simulation finished,"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0 ticks, time elapsed 0.285992677134 seconds\n",
        "Starting simulation with  0 ticks\n",
        "Simulation finished,"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0 ticks, time elapsed 0.283181251529 seconds\n",
        "Starting simulation with  0 ticks\n",
        "Simulation finished,"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 0 ticks, time elapsed 0.290546981349 seconds\n",
        "1 loops, best of 3: 284 ms per loop\n"
       ]
      }
     ],
     "prompt_number": 26
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "So I got it down from 3 seconds to 0.3 seconds, **~10 times faster**.\n",
      "That means that a 500 generations simulations (like the ones I used in my previous work) should take ~100 seconds,\n",
      "in comparison to several hours that they took in my Java implementation.\n",
      "\n",
      "**This is progress.**\n",
      "\n",
      "The profiler for 5 generations:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "%%prun\n",
      "p,g = simulation.run(5,0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "Starting simulation with  5 ticks\n",
        "Simulation finished,"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        " 5 ticks, time elapsed 2.63202320618 seconds\n",
        "\n"
       ]
      }
     ],
     "prompt_number": 19
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "11789 function calls in 0.193 seconds\n",
      "\n",
      "ncalls  tottime  percall  cumtime  percall filename:lineno(function)\n",
      "\n",
      "5696    1.792    0.000    1.792    0.000 model.py:59(find_row_nums)\n",
      "\n",
      "     4823    0.160    0.000    0.160    0.000 {numpy.core.multiarray.array}\n",
      "\n",
      "        6    0.088    0.015    0.089    0.015 model.py:95(clear_empty_classes)\n",
      "\n",
      "    11641    0.085    0.000    0.085    0.000 {method 'nonzero' of 'numpy.ndarray' objects}\n",
      "\n",
      "     5945    0.074    0.000    0.074    0.000 {method 'put' of 'numpy.ndarray' objects}\n",
      "\n",
      "    11641    0.072    0.000    0.184    0.000 model.py:49(genome_to_num)\n",
      "\n",
      "        6    0.071    0.012    0.128    0.021 shape_base.py:348(_replace_zero_by_x_arrays)\n",
      "\n",
      "        6    0.063    0.010    2.228    0.371 model.py:106(mutation_explicit_genomes)\n",
      "\n",
      "    17610    0.061    0.000    0.061    0.000 {method 'sum' of 'numpy.ndarray' objects}\n",
      "\n",
      "     5945    0.060    0.000    0.093    0.000 model.py:38(hamming_fitness_genome)\n",
      "\n",
      "        7    0.043    0.006    0.213    0.030 shape_base.py:11(apply_along_axis)\n",
      "\n",
      "     4773    0.014    0.000    0.014    0.000 {method 'any' of 'numpy.ndarray' objects}\n",
      "\n",
      "Only functions with >0.009 total time are shown.\n",
      "Not much I can do here - maybe some playing around with `numpy` `dtype`s will help,\n",
      "or moving some stuff to `cython`, but based on my experiments I doubt it.\n",
      "\n",
      "At any case I'm not sure it is worth it, and apart from `find_row_nums` there is no function that is a clear case for improvement.\n",
      "\n",
      "## FIN"
     ]
    }
   ],
   "metadata": {}
  }
 ]
}