{
 "metadata": {
  "name": "",
  "signature": "sha256:1ceb7f4a35d18593571f810414beb60bbee16b8291d496980dfaafdf217ceea4"
 },
 "nbformat": 3,
 "nbformat_minor": 0,
 "worksheets": [
  {
   "cells": [
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.stats import poisson, nbinom\n",
      "from scipy import __version__ as scipy_ver\n",
      "print \"scipy:\", scipy_ver\n",
      "np.set_printoptions(linewidth=80)\n",
      "np.set_printoptions(precision=3)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "# Muller's ratchet\n",
      "\n",
      "## Simulation"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy import stats"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 2
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def plot_freq(x, color='k', alpha=1):\n",
      "    bar(range(G), x, align=\"center\", color=color, alpha=alpha)\n",
      "    xlabel(\"# mutations\")\n",
      "    ylabel(\"frequency\")"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 3
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "U = 0.003\n",
      "s = 0.001\n",
      "N = 1000\n",
      "G = 25"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 50
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def mutation_matrix(U):\n",
      "    mutation_rv = stats.poisson(U)\n",
      "    M = zeros((G,G))\n",
      "    for i in range(G):\n",
      "        for j in range(i+1):\n",
      "            M[i,j] = mutation_rv.pmf(i-j)\n",
      "    M[G-1,:] += 1-M.sum(axis=0)\n",
      "    \n",
      "    assert (M.sum(axis=0)==1).all()\n",
      "    return M\n",
      "\n",
      "M = mutation_matrix(U)\n",
      "print M[:5,:5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  9.970e-01   0.000e+00   0.000e+00   0.000e+00   0.000e+00]\n",
        " [  2.991e-03   9.970e-01   0.000e+00   0.000e+00   0.000e+00]\n",
        " [  4.487e-06   2.991e-03   9.970e-01   0.000e+00   0.000e+00]\n",
        " [  4.487e-09   4.487e-06   2.991e-03   9.970e-01   0.000e+00]\n",
        " [  3.365e-12   4.487e-09   4.487e-06   2.991e-03   9.970e-01]]\n"
       ]
      }
     ],
     "prompt_number": 51
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "gjfhjdgf"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def selection_matrix(s):\n",
      "    S = diag((1-s)**arange(G))\n",
      "    return S\n",
      "\n",
      "S = selection_matrix(s)\n",
      "print S[:5,:5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[ 1.     0.     0.     0.     0.   ]\n",
        " [ 0.     0.999  0.     0.     0.   ]\n",
        " [ 0.     0.     0.998  0.     0.   ]\n",
        " [ 0.     0.     0.     0.997  0.   ]\n",
        " [ 0.     0.     0.     0.     0.996]]\n"
       ]
      }
     ],
     "prompt_number": 53
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "E = M.dot(S)\n",
      "print E[:5,:5]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[[  9.970e-01   0.000e+00   0.000e+00   0.000e+00   0.000e+00]\n",
        " [  2.991e-03   9.960e-01   0.000e+00   0.000e+00   0.000e+00]\n",
        " [  4.487e-06   2.988e-03   9.950e-01   0.000e+00   0.000e+00]\n",
        " [  4.487e-09   4.482e-06   2.985e-03   9.940e-01   0.000e+00]\n",
        " [  3.365e-12   4.482e-09   4.478e-06   2.982e-03   9.930e-01]]\n"
       ]
      }
     ],
     "prompt_number": 54
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x_msb = stats.poisson(U/s).pmf(range(G))"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 55
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x_free = zeros(G)\n",
      "x_free[0] = 1"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 56
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Evolve to MSB:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = x_free.copy()\n",
      "for _ in range(10000):\n",
      "    x = E.dot(x)\n",
      "    x /= x.sum()\n",
      "plot_freq(x, color='w')\n",
      "plot_freq(x_msb, color='b', alpha=0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAEPCAYAAABcA4N7AAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG0NJREFUeJzt3X1slfX9//HXwdZ2g4qArPes9AZasC2FlsKQrRPapkQ7\ngS3rRIVRtbIw1G1qljjbqiE0m8mQuoiozCl2RA2WQa0M4/GmgF20wJaCKTfVUqkQBgMKFFo+vz/8\neb4U6OecA1y9fT6Sk5zrOu/P6fs6V05fva6r13W5jDFGAAB0YVBPNwAA6N0ICgCAFUEBALAiKAAA\nVgQFAMCKoAAAWDkaFNXV1UpMTFRCQoLKysoueX3NmjVKTU1VSkqKpk2bpp07d3pei4mJUUpKitLS\n0jR58mQn2wQAWLicOo+io6NDY8eO1ebNmxUZGamMjAxVVFQoKSnJU7N161aNGzdOQ4cOVXV1tUpK\nSrRt2zZJ0ujRo/Xpp59q+PDhTrQHAPCRY1sUtbW1io+PV0xMjAIDA1VQUKDKyspONVOnTtXQoUMl\nSZmZmTpw4ECn1zkXEAB6nmNB0dzcrOjoaM90VFSUmpubu6x/6aWXNGvWLM+0y+XSzJkzlZ6erlWr\nVjnVJgDAiwCn3tjlcvlc+/777+vll19WTU2NZ15NTY3Cw8N1+PBhZWdnKzExUdOnT3eiVQCAhWNB\nERkZqaamJs90U1OToqKiLqnbuXOn7rvvPlVXV2vYsGGe+eHh4ZKkkSNHavbs2aqtrb0kKOLj47V3\n716HlgAA+qe4uDjt2bPH9wHGIefOnTOxsbFm//79pq2tzaSmppr6+vpONV988YWJi4szW7du7TS/\ntbXVHD9+3BhjzMmTJ80PfvAD8+67717yMxxsv1coLi7u6RYcxfL1Xf152Yzp/8vn7+9Ox7YoAgIC\nVF5ertzcXHV0dKiwsFBJSUlauXKlJKmoqEhPPvmkjh49qkWLFkmSAgMDVVtbq5aWFs2ZM0eS1N7e\nrnnz5iknJ8epVgEAFo4FhSTl5eUpLy+v07yioiLP8xdffFEvvvjiJeNiY2O1fft2J1sDAPiIM7N7\nsaysrJ5uwVEsX9/Vn5dN6v/L5y/HTrjrDi6Xi3MtAMBP/v7uZIsCAGBFUAAArAgKAIAVQQEAsCIo\nAABWBAUAwIqgAABYERQAACuCAgBgRVAAAKwICgCAFUEBALAiKAAAVgQFAMCKoAAAWBEUAAArggIA\nYOXoPbPROwQHh6utLchaExTUpjNnDvpVC2BgICgGgLa2IBUXN1prSktj/K4FMDCw6wkAYEVQAACs\nCAoAgBVBAQCwIigAAFYEBQDAiqAAAFgRFAAAK4ICAGBFUAAArAgKAIAVQQEAsCIoAABWBAUAwIqg\nAABYORoU1dXVSkxMVEJCgsrKyi55fc2aNUpNTVVKSoqmTZumnTt3+jwWANA9HAuKjo4OLV68WNXV\n1aqvr1dFRYV27drVqSY2NlYffvihdu7cqT/84Q+6//77fR4LAOgejgVFbW2t4uPjFRMTo8DAQBUU\nFKiysrJTzdSpUzV06FBJUmZmpg4cOODzWABA93AsKJqbmxUdHe2ZjoqKUnNzc5f1L730kmbNmnVF\nYwEAznHsntkul8vn2vfff18vv/yyampq/B5bUlLieZ6VlaWsrCyfxwLAQOB2u+V2u694vGNBERkZ\nqaamJs90U1OToqKiLqnbuXOn7rvvPlVXV2vYsGF+jZU6BwUA4FIX/xFdWlrq13jHdj2lp6eroaFB\njY2NOnv2rNauXav8/PxONV9++aXmzJmj1157TfHx8X6NBQB0D8e2KAICAlReXq7c3Fx1dHSosLBQ\nSUlJWrlypSSpqKhITz75pI4ePapFixZJkgIDA1VbW9vlWABA93MZY0xPN3GlXC6X+nD73cblilFx\ncaO1prQ0RsY0+lULoG/y93cnZ2YDAKwICgCAFUEBALAiKAAAVgQFAMCKoAAAWBEUAAArggIAYEVQ\nAACsCAoAgBVBAQCwIigAAFYEBQDAiqAAAFgRFAAAK4ICAGBFUAAArAgKAIAVQQEAsCIoAABWBAUA\nwIqgAABYERQAACuCAgBgFdDTDaDvCg4OV1tbkLUmKKhNZ84c7KaOADiBoMAVa2sLUnFxo7WmtDSm\nO1oB4CB2PQEArAgKAIAVQQEAsCIoAABWBAUAwIqgAABYERQAACuCAgBgRVAAAKwcDYrq6molJiYq\nISFBZWVll7y+e/duTZ06VcHBwXrmmWc6vRYTE6OUlBSlpaVp8uTJTrYJALBw7BIeHR0dWrx4sTZv\n3qzIyEhlZGQoPz9fSUlJnpoRI0ZoxYoVevvtty8Z73K55Ha7NXz4cKdaBAD4wLEtitraWsXHxysm\nJkaBgYEqKChQZWVlp5qRI0cqPT1dgYGBl30PY4xT7QEAfORYUDQ3Nys6OtozHRUVpebmZp/Hu1wu\nzZw5U+np6Vq1apUTLQIAfODYrieXy3VV42tqahQeHq7Dhw8rOztbiYmJmj59+jXqDgDgK69BMWnS\nJC1cuFB33nmnhg0b5vMbR0ZGqqmpyTPd1NSkqKgon8eHh4dL+mb31OzZs1VbW3vZoCgpKfE8z8rK\nUlZWls8/oy/jXhAAfOV2u+V2u694vNeg+Pvf/67Vq1crIyND6enp+uUvf6mcnByvWwzp6elqaGhQ\nY2OjIiIitHbtWlVUVFy29uJjEadOnVJHR4dCQkLU2tqqTZs2qbi4+LJjLwyKgYR7QQDw1cV/RJeW\nlvo13mtQJCQkaOnSpXr66ae1YcMGLVy4UIMGDdLChQv14IMPdvlfSQEBASovL1dubq46OjpUWFio\npKQkrVy5UpJUVFSklpYWZWRk6Pjx4xo0aJCWL1+u+vp6HTp0SHPmzJEktbe3a968ecrJyfFrwQAA\n14ZPxyh27Nih1atX65133tHcuXN155136uOPP9att96q7du3dzkuLy9PeXl5neYVFRV5noeFhXXa\nPfWtIUOGWN8XANB9fDpGMXToUN17770qKytTUNA3+8WnTJmimpoaxxsEAPQsr0HxxhtvKDY29rKv\nrVu37po3BADoXbyeR/Hiiy/q2LFjnumjR4/q8ccfd7QpAEDv4TUoqqqqdOONN3qmhw0bpo0bNzra\nFACg9/AaFOfPn9eZM2c806dPn9bZs2cdbQoA0Ht4PUYxb948zZgxQwsXLpQxRqtXr9Y999zTHb0B\nAHoBr0Hx2GOPKSUlRZs3b5bL5dITTzyh3Nzc7ugNANAL+HQexeXOhwAADAxej1G89dZbSkhI0A03\n3KCQkBCFhITohhtu6I7eAAC9gNctikcffVQbNmzodMMhAMDA4XWLIiwsjJAAgAHM6xZFenq6fv7z\nn+uOO+7Q9ddfL+mbe018e9E+AED/5jUo/ve//+k73/mONm3a1Gk+QQEAA4PXoPjrX//aDW0AAHor\nr8coPv/8c82YMUPjx4+XJO3cuVNPP/20440BAHoHr0Fx3333aenSpZ7jE8nJyV3eqQ4A0P94DYpT\np04pMzPTM+1yuRQYGOhoUwCA3sNrUIwcOVJ79uzxTL/55psKDw93tCkAQO/h9WB2eXm57r//fu3e\nvVsREREaPXq01qxZ0x29AQB6Aa9BERcXp/fee0+tra06f/68QkJCuqMvAEAv4TUoSktL5XK5ZIyR\ny+XyzH/iiSccbQwA0Dt4DYrBgwd7AuL06dPasGGDxo0b53hjAIDewWtQ/O53v+s0/cgjjygnJ8ex\nhgAAvYvX/3q6WGtrq5qbm53oBQDQC3ndokhOTvY8P3/+vA4dOsTxCQAYQLwGxT/+8Y//Kw4IUGho\nKCfcAcAA4jUoLr6b3YkTJzpNDx8+/Np2BADoVbwGxcSJE/Xll19q2LBhkqSjR49q1KhRcrlccrlc\n2rdvn+NNAgB6jteD2dnZ2dqwYYOOHDmiI0eOaOPGjcrJydH+/fsJCQAYALwGxdatWzVr1izPdF5e\nnrZs2eJoUwCA3sPrrqeIiAg9/fTTuuuuu2SM0euvv67IyMju6A0A0At43aKoqKjQoUOHNHv2bM2Z\nM0eHDh3ifhQAMIB43aIYMWKEnn32WbW2tmrw4MHd0RMAoBfxukWxZcsWjRs3TomJiZKkHTt26Fe/\n+pXjjQEAegevQfHQQw+purpaN910kyQpNTVVH3zwgeONAQB6B5+u9TRq1KhO0wEBXvdYAQD6Ca9B\nMWrUKNXU1EiSzp49qz/96U9KSkry6c2rq6uVmJiohIQElZWVXfL67t27NXXqVAUHB+uZZ57xaywA\noHt4DYrnn39ezz33nJqbmxUZGam6ujo999xzXt+4o6NDixcvVnV1terr61VRUaFdu3Z1qhkxYoRW\nrFhxyaXMfRkLAOge1n1I7e3tevDBB/X666/7/ca1tbWKj49XTEyMJKmgoECVlZWdtkZGjhypkSNH\nauPGjX6PBQB0D+sWRUBAgL744gu1tbX5/cbNzc2Kjo72TEdFRfl8H4urGQsAuLa8HpWOjY3VLbfc\novz8fH33u9+VJLlcLv3mN7+xjrvw/tr+upqx6J2Cg8PV1hZkrQkKatOZMwe7qSMAvuoyKO6++269\n+uqrWr9+vR5++GGdP39eJ0+e9PmNIyMj1dTU5JluampSVFTUNR9bUlLieZ6VlaWsrCyfe0T3aWsL\nUnFxo7WmtDSmO1oBBhy32y23233F47sMik8//VRfffWVRo0apV//+tcyxvj1xunp6WpoaFBjY6Mi\nIiK0du3aLi/9cfF7+zP2wqAAAFzq4j+iS0tL/RrfZVA88MADmjFjhvbt26dJkyZ1es2X+1AEBASo\nvLxcubm56ujoUGFhoZKSkrRy5UpJUlFRkVpaWpSRkaHjx49r0KBBWr58uerr6zVkyJDLjgUAdL8u\ng2LJkiVasmSJHnjgAT3//PNX9OZ5eXnKy8vrNK+oqMjzPCwsrNMuJm9jAQDdz6fzKAAAA5dPl/AA\nAAxcBAUAwIqgAABYERQAACuCAgBgRVAAAKwICgCAFUEBALAiKAAAVgQFAMCKoAAAWBEUAAArggIA\nYEVQAACsCAoAgBVBAQCwIigAAFYEBQDAiqAAAFgRFAAAK4ICAGBFUAAArAgKAIAVQQEAsCIoAABW\nBAUAwIqgAABYERQAACuCAgBgRVAAAKwICgCAFUEBALAiKAAAVgQFAMCKoAAAWDkaFNXV1UpMTFRC\nQoLKysouW7NkyRIlJCQoNTVVdXV1nvkxMTFKSUlRWlqaJk+e7GSbAACLAKfeuKOjQ4sXL9bmzZsV\nGRmpjIwM5efnKykpyVNTVVWlPXv2qKGhQZ988okWLVqkbdu2SZJcLpfcbreGDx/uVIsAAB84tkVR\nW1ur+Ph4xcTEKDAwUAUFBaqsrOxUs379es2fP1+SlJmZqWPHjunrr7/2vG6Mcao9AICPHAuK5uZm\nRUdHe6ajoqLU3Nzsc43L5dLMmTOVnp6uVatWOdUmAMALx3Y9uVwun+q62mr4+OOPFRERocOHDys7\nO1uJiYmaPn36JXUlJSWe51lZWcrKyrqSdnuF4OBwtbUFdfl6UFCbzpw52I0d9Qxvn4M0cD4L4Fpw\nu91yu91XPN6xoIiMjFRTU5NnuqmpSVFRUdaaAwcOKDIyUpIUEREhSRo5cqRmz56t2tpar0HR17W1\nBam4uLHL10tLY7qrlR7l7XOQBs5nAVwLF/8RXVpa6td4x3Y9paenq6GhQY2NjTp79qzWrl2r/Pz8\nTjX5+fn629/+Jknatm2bbrzxRoWGhurUqVM6ceKEJKm1tVWbNm1ScnKyU60CACwc26IICAhQeXm5\ncnNz1dHRocLCQiUlJWnlypWSpKKiIs2aNUtVVVWKj4/X4MGDtXr1aklSS0uL5syZI0lqb2/XvHnz\nlJOT41SrAAALx4JCkvLy8pSXl9dpXlFRUafp8vLyS8bFxsZq+/btTrYGAPARZ2YDAKwICgCAFUEB\nALAiKAAAVgQFAMCKoAAAWBEUAAArggIAYEVQAACsCAoAgBVBAQCwIigAAFYEBQDAiqAAAFgRFAAA\nK4ICAGBFUAAArAgKAIAVQQEAsCIoAABWBAUAwCqgpxsArlZwcLja2oKsNUFBbTpz5mA3dQT0LwQF\n+ry2tiAVFzdaa0pLY7qjFaBfYtcTAMCKoAAAWBEUAAArggIAYEVQAACsCAoAgBVBAQCwIigAAFYE\nBQDAiqAAAFhxCQ8MKFwXCvAfQYEBhetCAf5zdNdTdXW1EhMTlZCQoLKyssvWLFmyRAkJCUpNTVVd\nXZ1fYwEAznMsKDo6OrR48WJVV1ervr5eFRUV2rVrV6eaqqoq7dmzRw0NDXrhhRe0aNEin8cOBI2N\n7p5uwVH9ffncbndPt+CY/rxsUv9fPn85FhS1tbWKj49XTEyMAgMDVVBQoMrKyk4169ev1/z58yVJ\nmZmZOnbsmFpaWnwaOxD091+kvX35goPD5XLFWB/BweFdju/Pv2z687JJ/X/5/OXYMYrm5mZFR0d7\npqOiovTJJ594rWlubtZXX33ldSzgNI5nAN9wLChcLpdPdcYYp1oAus3l/5vqmEpL/+qZ4r+p0GcZ\nh2zdutXk5uZ6ppcuXWqWLVvWqaaoqMhUVFR4pseOHWtaWlp8GmuMMXFxcUYSDx48ePDw4xEXF+fX\n73PHtijS09PV0NCgxsZGRUREaO3ataqoqOhUk5+fr/LychUUFGjbtm268cYbFRoaqhEjRngdK0l7\n9uxxqn0AwP/nWFAEBASovLxcubm56ujoUGFhoZKSkrRy5UpJUlFRkWbNmqWqqirFx8dr8ODBWr16\ntXUsAKD7uYzhIAEAoGt9/lpPJSUlioqKUlpamtLS0lRdXd3TLV21/n6yYUxMjFJSUpSWlqbJkyf3\ndDtXbeHChQoNDVVycrJn3n//+19lZ2drzJgxysnJ0bFjx3qww6tzueXrT9+7pqYm/fjHP9b48eN1\n880369lnn5XUP9ZhV8vm9/rz64hGL1RSUmKeeeaZnm7jmmlvbzdxcXFm//795uzZsyY1NdXU19f3\ndFvXVExMjDly5EhPt3HNfPjhh+azzz4zN998s2feI488YsrKyowxxixbtsw89thjPdXeVbvc8vWn\n793BgwdNXV2dMcaYEydOmDFjxpj6+vp+sQ67WjZ/11+f36KQ1K/+xXagnGzYn9bZ9OnTNWzYsE7z\nLjyZdP78+Xr77bd7orVr4nLLJ/WfdRgWFqYJEyZIkoYMGaKkpCQ1Nzf3i3XY1bJJ/q2/fhEUK1as\nUGpqqgoLC/vk5uGFujoJsT9xuVyaOXOm0tPTtWrVqp5uxxFff/21QkNDJUmhoaH6+uuve7ija68/\nfe++1djYqLq6OmVmZva7dfjtsk2ZMkWSf+uvTwRFdna2kpOTL3msX79eixYt0v79+7V9+3aFh4fr\nt7/9bU+3e1V8PVGxL6upqVFdXZ3eeecdPffcc/roo496uiVHuVyufrde+9v3TpJOnjypuXPnavny\n5QoJCen0Wl9fhydPntRPf/pTLV++XEOGDPF7/fWJy4z/85//9Knu3nvv1e233+5wN86KjIxUU1OT\nZ7qpqUlRUVE92NG1Fx7+zfWRRo4cqdmzZ6u2tlbTp0/v4a6urdDQULW0tCgsLEwHDx7U9773vZ5u\n6Zq6cHn6w/fu3Llzmjt3ru6++27dcccdkvrPOvx22e666y7Psvm7/vrEFoXNwYP/d0mEdevWdfrP\njL7owhMVz549q7Vr1yo/P7+n27pmTp06pRMnTkiSWltbtWnTpj6/zi4nPz9fr7zyiiTplVde8XxB\n+4v+9L0zxqiwsFDjxo3TQw895JnfH9ZhV8vm9/q79sfZu9fdd99tkpOTTUpKivnJT35iWlpaerql\nq1ZVVWXGjBlj4uLizNKlS3u6nWtq3759JjU11aSmpprx48f3i+UrKCgw4eHhJjAw0ERFRZmXX37Z\nHDlyxMyYMcMkJCSY7Oxsc/To0Z5u84pdvHwvvfRSv/reffTRR8blcpnU1FQzYcIEM2HCBPPOO+/0\ni3V4uWWrqqrye/1xwh0AwKrP73oCADiLoAAAWBEUAAArggIAYEVQAACsCAoAgBVBgQHj97//vdxu\nt95++20tW7bMkZ9RWVmpXbt2+V1XXFys9957z5GegKtFUGDAqK2t1ZQpU/TBBx/ohz/8oSM/Y926\ndaqvr/e7rrS0VDNmzHCkJ+CqdcvpgUAPeuSRR0xKSooJCQkxEyZMMCEhISYlJcU89dRTl9TOnz/f\nLFq0yEyZMsXExsaa999/39xzzz0mKSnJLFiwwFM3ePBgz/M33njDLFiwwGzZssUMHz7cjB492qSl\npZm9e/eaF154wWRkZJjU1FQzd+5cc+rUKVNTU3NJ3fz5882bb75pjDFm8+bNJi0tzSQnJ5uFCxea\ntrY2Y4wx3//+901xcbGZOHGiSU5ONrt37zbGGON2uz1n3aalpZkTJ044+XFiACIoMCD861//MkuW\nLDHnzp0z06ZN67JuwYIF5he/+IUxxpjKykoTEhJi/vOf/5jz58+bSZMmmR07dhhjjBkyZIhnzJtv\nvukJkQULFpi33nrL89qFN2h6/PHHzYoVKy5b9+306dOnTXR0tGloaDDGGHPPPfeYP//5z8aYb274\nVF5ebowx5i9/+Yu59957jTHG3H777WbLli3GGGNaW1tNe3v7lXxEQJfY9YQB4dNPP1VKSop27dql\npKQka+23V9K8+eabFRYWpvHjx8vlcmn8+PFqbGz0+rPMBVfF+fe//63p06crJSVFa9as6bS7yVx0\n9RxjjD7//HONHj1a8fHxkr65Yc6HH37oqZkzZ44kaeLEiZ5epk2bpocfflgrVqzQ0aNHdd1113nt\nEfBHn7jMOHClduzYoQULFujAgQO66aabdOrUKRljNHHiRG3ZskXBwcGXjLn++uslSYMGDVJQUJBn\n/qBBg9Te3i6p831DTp8+3Wn8ha8tWLBA69evV3Jysl555RW53e7L1nU1zxjTad63/Vx33XWeXh57\n7DHddttt2rhxo6ZNm6Z3331XY8eOtX8wgB/YokC/lpqaqrq6Oo0ZM0a7du3Srbfeqk2bNumzzz67\nbEj4KjQ0VLt379b58+e1bt06zy/zkJAQHT9+3FN38uRJhYWF6dy5c3rttde6rJO+CYmxY8eqsbFR\ne/fulSS9+uqr+tGPfmTtZe/evRo/frweffRRZWRk6PPPP7/i5QIuh6BAv3f48GENHz5ckrR7924l\nJiZa6y/8C76ru5otW7ZMt912m6ZNm6aIiAjP/IKCAv3xj3/UpEmTtG/fPj311FPKzMzULbfc0mmX\n18V13woKCtLq1av1s5/9TCkpKQoICNADDzxw2V6+nV6+fLmSk5OVmpqq66+/Xnl5eb58LIDPuMw4\nAMCKLQoAgBVBAQCwIigAAFYEBQDAiqAAAFgRFAAAK4ICAGBFUAAArP4fYCgqLmNaoVsAAAAASUVO\nRK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x12065630>"
       ]
      }
     ],
     "prompt_number": 57
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Evolve to first click of the ratchet:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "x = x_free.copy()\n",
      "t = 0\n",
      "\n",
      "print x[:5]\n",
      "\n",
      "while x[0] > 0:\n",
      "    t += 1\n",
      "    if t % N == 0: print t, x[:5]\n",
      "    x = E.dot(x)    \n",
      "    x /= x.sum()\n",
      "    assert allclose(x.sum(),1)\n",
      "    \n",
      "    n = multinomial(N, x)\n",
      "    assert len(n) == G\n",
      "    x = n/float(N)\n",
      "\n",
      "print \"CLICK\"\n",
      "print t, x[:5]\n",
      "\n",
      "plot_freq(x, color='b', alpha=0.5)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "[ 1.  0.  0.  0.  0.]\n",
        "CLICK"
       ]
      },
      {
       "output_type": "stream",
       "stream": "stdout",
       "text": [
        "\n",
        "902 [ 0.     0.086  0.173  0.543  0.178]\n"
       ]
      },
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEPCAYAAACk43iMAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAG+ZJREFUeJzt3X9Q1HXix/HXGlT+INI0CRYPFZJVARE4bewHZRfhJWNq\nc9gPRdHIxuusqayZJqEaT6Zp5lK8O7rr9PqFTHaFmVJHtelIHV16dneg4Y+91i3qMi1/kMjy/v7R\ntV8RdMH8sMDn+ZjZmf3svj8fXqJ+Xvv5fPbz+TiMMUYAAFvqE+oAAIDQoQQAwMYoAQCwMUoAAGyM\nEgAAG6MEAMDGLC2ByspKJSYmKiEhQcXFxe2OcbvdSk1N1dixY5WZmWllHADAKRxWnSfg9/s1atQo\nVVVVKSYmRhkZGSorK5PL5QqMOXTokCZNmqQ333xTTqdTX331lQYPHmxFHABAOyzbEqipqVF8fLzi\n4uIUHh6u3NxcVVRUtBrz0ksvacaMGXI6nZJEAQBAF7OsBHw+n2JjYwPTTqdTPp+v1Zj6+np9/fXX\nuvbaa5Wenq7nn3/eqjgAgHaEWbVgh8MRdMyJEye0bds2vf322zp27JiuuOIKTZw4UQkJCVbFAgCc\nxLISiImJkdfrDUx7vd7Abp8fxMbGavDgwerbt6/69u2rq6++Wjt27GhTAvHx8dqzZ49VUQGgVxo5\ncqR279595kHGIidOnDAjRoww+/btM8ePHzcpKSmmtra21Zi6ujozefJk09zcbI4ePWrGjh1r/v3v\nf7dZloUxu8TSpUtDHeFH6cn5e3J2Y8gfaj09f0fWnZZtCYSFhamkpERZWVny+/3Kz8+Xy+VSaWmp\nJKmgoECJiYm68cYblZycrD59+mjBggUaPXq0VZEAAKewrAQkKTs7W9nZ2a1eKygoaDV9//336/77\n77cyBgDgNDhjuAv09JPgenL+npxdIn+o9fT8HWHZyWLnksPhUA+ICQDdSkfWnWwJAICNUQIAYGOU\nAADYGCUAADZGCQCAjVECAGBjlAAA2BglAAA2RgkAgI1RAgBgY5QAANgYJQAANkYJAICNUQIAYGOU\nAADYmKV3FgPOxkMPFauhoTHUMdoVFdVXy5cvCXUM4JyhBNDtNDQ0Ki6uMNQx2uXxFIY6AnBOsTsI\nAGyMEgAAG6MEAMDGKAEAsDFKAABsjBIAABujBADAxigBALAxSgAAbMzSEqisrFRiYqISEhJUXFzc\n5n23263IyEilpqYqNTVVTzzxhJVxAACnsOyyEX6/X4sWLVJVVZViYmKUkZGhnJwcuVyuVuOuueYa\nrV+/3qoYAIAzsGxLoKamRvHx8YqLi1N4eLhyc3NVUVHRZpwxxqoIAIAgLCsBn8+n2NjYwLTT6ZTP\n52s1xuFwqLq6WikpKZoyZYpqa2utigMAaIdlu4McDkfQMePHj5fX61W/fv20adMmTZs2TZ988olV\nkQAAp7CsBGJiYuT1egPTXq9XTqez1ZiIiIjA8+zsbN199936+uuvNWjQoDbLKywsDDzPzMxUZmbm\nOc8MAD2Z2+2W2+3u1DyWlUB6errq6+vl8XgUHR2t8vJylZWVtRrzxRdf6NJLL5XD4VBNTY2MMe0W\ngNS6BAAAbZ36AbmoqCjoPJaVQFhYmEpKSpSVlSW/36/8/Hy5XC6VlpZKkgoKCrRu3Tr97ne/U1hY\nmPr166e1a9daFQcA0A6H6QFfz3E4HHyLyEby8gq79Z3F1qwpDHUMoEM6su7kjGEAsDFKAABsjBIA\nABujBADAxigBALAxSgAAbIwSAAAbowQAwMYoAQCwMUoAAGyMEgAAG6MEAMDGKAEAsDFKAABsjBIA\nABujBADAxigBALAxSgAAbIwSAAAbowQAwMYoAQCwMUoAAGyMEgAAG6MEAMDGKAEAsDFKAABsjBIA\nABujBADAxigBALAxS0ugsrJSiYmJSkhIUHFx8WnHffjhhwoLC9Nf/vIXK+MAAE5hWQn4/X4tWrRI\nlZWVqq2tVVlZmerq6todt2TJEt14440yxlgVBwDQDstKoKamRvHx8YqLi1N4eLhyc3NVUVHRZtzK\nlSs1c+ZMDRkyxKooAIDTsKwEfD6fYmNjA9NOp1M+n6/NmIqKCi1cuFCS5HA4rIoDAGhHmFUL7sgK\nffHixVq+fLkcDoeMMWfcHVRYWBh4npmZqczMzHOQEgB6D7fbLbfb3al5LCuBmJgYeb3ewLTX65XT\n6Ww15qOPPlJubq4k6auvvtKmTZsUHh6unJycNss7uQQAAG2d+gG5qKgo6DyWlUB6errq6+vl8XgU\nHR2t8vJylZWVtRqzd+/ewPO5c+dq6tSp7RYAAMAalpVAWFiYSkpKlJWVJb/fr/z8fLlcLpWWlkqS\nCgoKrPrRAIAOsqwEJCk7O1vZ2dmtXjvdyn/16tVWRgEAtIMzhgHAxigBALAxSgAAbIwSAAAbowQA\nwMYoAQCwMUoAAGyMEgAAG6MEAMDGKAEAsLGgJZCWlqZVq1bp4MGDXZEHANCFgpbA2rVr5fP5lJGR\nodzcXL355pvcBhIAeomgJZCQkKBly5bpk08+0a233qp58+Zp2LBhWrp0qb7++uuuyAgAsEiHjgns\n2LFD9913nx544AHNmDFDL7/8siIiInTddddZnQ8AYKGgl5JOS0tTZGSk5s+fr+LiYl1wwQWSpIkT\nJ2rr1q2WBwQAWCdoCbz88ssaMWJEu++9+uqr5zwQAKDrBN0d9Mc//lGHDh0KTB88eFCPPPKIpaEA\nAF0jaAls3LhRF198cWB64MCBeuONNywNBQDoGkFLoKWlRd99911gurGxUU1NTZaGAgB0jaDHBG67\n7TZNnjxZ8+bNkzFGq1ev1uzZs7siGwDAYkFLYMmSJUpOTlZVVZUcDoceffRRZWVldUU2AIDFgpaA\nJGVnZys7O9vqLACALhb0mMArr7yihIQEXXTRRYqIiFBERIQuuuiirsgGALBY0C2BBx98UBs2bJDL\n5eqKPACALhR0SyAqKooCAIBeKuiWQHp6un7xi19o2rRpOv/88yVJDodD06dPtzwcAMBaQUvgm2++\nUd++ffXWW2+1ep0SAICeL2gJrFmzpgtiAABCIegxgV27dmny5MkaM2aMJOnjjz/WE0880aGFV1ZW\nKjExUQkJCSouLm7zfkVFhVJSUpSamqq0tDS98847nYwPAPgxgpbAggULtGzZssDxgKSkJJWVlQVd\nsN/v16JFi1RZWana2lqVlZWprq6u1Zjrr79eO3bs0Pbt27VmzRrdeeedZ/nHAACcjaAlcOzYMU2Y\nMCEw7XA4FB4eHnTBNTU1io+PV1xcnMLDw5Wbm6uKiopWY/r37x94fuTIEQ0ePLgz2QEAP1LQEhgy\nZIh2794dmF63bp0uu+yyoAv2+XyKjY0NTDudTvl8vjbjXnvtNblcLmVnZ2vFihUdzQ0AOAeCHhgu\nKSnRnXfeqZ07dyo6OlrDhw/Xiy++GHTBDoejQwGmTZumadOmacuWLbrjjju0a9eudscVFhYGnmdm\nZiozM7NDywcAu3C73XK73Z2aJ2gJjBw5Um+//baOHj2qlpYWRUREdGjBMTEx8nq9gWmv1yun03na\n8VdddZWam5t14MABXXLJJW3eP7kEAABtnfoBuaioKOg8QUugqKhIDodDxphWn+4fffTRM86Xnp6u\n+vp6eTweRUdHq7y8vM0B5T179mjEiBFyOBzatm2bJLVbAAAAawQtgf79+wdW/o2NjdqwYYNGjx4d\nfMFhYSopKVFWVpb8fr/y8/PlcrlUWloqSSooKNArr7yi5557TuHh4RowYIDWrl37I/84AIDOcBhj\nTGdmOH78uG644Qa99957VmVq44ctEdhDXl6h4uIKQx2jXR5PodasKQx1DKBDOrLuDPrtoFMdPXq0\n3W/5AAB6nqC7g5KSkgLPW1pa9OWXXwY9HgAA6BmClsDrr7/+/4PDwjR06NAOnSwGAOj+gpbAqXcR\nO3z4cKvpQYMGndtEAIAuE7QExo8fr08//VQDBw6UJB08eFDDhg2Tw+GQw+HQ3r17LQ8JALBG0APD\nP/vZz7RhwwYdOHBABw4c0BtvvKEbbrhB+/btowAAoIcLWgLvv/++pkyZEpjOzs5WdXW1paEAAF0j\n6O6g6OhoPfHEE7r99ttljNFLL72kmJiYrsgGALBY0C2BsrIyffnll7r55ps1ffp0ffnllx26nwAA\noPsLuiVwySWXaMWKFTp69Gir6/8DAHq+oFsC1dXVGj16tBITEyVJO3bs0N133215MACA9YKWwOLF\ni1VZWRm461dKSkqXXjcIAGCdDl07aNiwYa2mw8KC7kUCAPQAQdfmw4YN09atWyVJTU1NWrFihVwu\nl+XBAADWC7ol8Pvf/16rVq2Sz+dTTEyMtm/frlWrVnVFNgCAxc64JdDc3Kxf/epXeumll7oqDwCg\nC51xSyAsLEz/+c9/dPz48a7KAwDoQkGPCYwYMUJXXnmlcnJy1K9fP0nf363mvvvuszwcAMBap90S\nuOOOOyRJ69ev10033aSWlhYdOXJER44caXM5aQBAz3TaLYGPPvpIn332mYYNG6Zf/vKX3OMXAHqh\n05bAXXfdpcmTJ2vv3r1KS0tr9R73EQCA3uG0u4Puuece1dXVae7cudq3b1+rBwUAAL1Dh84TAAD0\nTh26bAQAoHeiBADAxigBALAxSgAAbIwSAAAbs7wEKisrlZiYqISEBBUXF7d5/8UXX1RKSoqSk5M1\nadIkffzxx1ZHAgD8j6V3h/H7/Vq0aJGqqqoUExOjjIwM5eTktLofwYgRI7R582ZFRkaqsrJSd955\npz744AMrYwEA/sfSLYGamhrFx8crLi5O4eHhys3NVUVFRasxV1xxhSIjIyVJEyZM0P79+62MBAA4\niaUl4PP5FBsbG5h2Op3y+XynHf/ss89qypQpVkYCAJzE0t1BDoejw2Pfffdd/elPfwrcyvJUhYWF\ngeeZmZnKzMz8kekAoHdxu91yu92dmsfSEoiJiZHX6w1Me71eOZ3ONuM+/vhjLViwQJWVlRo4cGC7\nyzq5BAAAbZ36AbmoqCjoPJbuDkpPT1d9fb08Ho+amppUXl6unJycVmM+/fRTTZ8+XS+88ILi4+Ot\njAMAOIWlWwJhYWEqKSlRVlaW/H6/8vPz5XK5VFpaKkkqKCjQY489poMHD2rhwoWSpPDwcNXU1FgZ\nCwDwPw7TA+4W43A4uKmNjeTlFSourjDUMdrl8RRqzZrCUMcAOqQj607OGAYAG6MEAMDGKAEAsDFK\nAABszNJvBwF29NBDxWpoaAx1jDaiovpq+fIloY6BboYS6IW660pIsseKqKGhsVt+u8njKQx1BHRD\nlEAv1F1XQhIrIqC74ZgAANgYJQAANkYJAICNUQIAYGOUAADYGCUAADZGCQCAjVECAGBjlAAA2Bgl\nAAA2RgkAgI1RAgBgY5QAANgYJQAANkYJAICNUQIAYGOUAADYGCUAADZGCQCAjVECAGBjlAAA2Jjl\nJVBZWanExEQlJCSouLi4zfs7d+7UFVdcoQsvvFBPPfWU1XEAACcJs3Lhfr9fixYtUlVVlWJiYpSR\nkaGcnBy5XK7AmEsuuUQrV67Ua6+9ZmUUAEA7LN0SqKmpUXx8vOLi4hQeHq7c3FxVVFS0GjNkyBCl\np6crPDzcyigAgHZYWgI+n0+xsbGBaafTKZ/PZ+WPBAB0gqUl4HA4rFw8AOBHsvSYQExMjLxeb2Da\n6/XK6XSe1bIKCwsDzzMzM5WZmfkj0wFA7+J2u+V2uzs1j6UlkJ6ervr6enk8HkVHR6u8vFxlZWXt\njjXGnHFZJ5cAAKCtUz8gFxUVBZ3H0hIICwtTSUmJsrKy5Pf7lZ+fL5fLpdLSUklSQUGBGhoalJGR\noW+//VZ9+vTR008/rdraWg0YMMDKaAAAWVwCkpSdna3s7OxWrxUUFASeR0VFtdplBADoOpwxDAA2\nZvmWQE/10EPFamhoDHWMNqKi+mr58iWhjgGgl6AETqOhoVFxcYWhjtGGx1MY6ggAehF2BwGAjVEC\nAGBjlAAA2BglAAA2RgkAgI1RAgBgY5QAANgYJQAANkYJAICNUQIAYGOUAADYGCUAADZGCQCAjVEC\nAGBjlAAA2BglAAA2RgkAgI1RAgBgY5QAANgYJQAANkYJAICNUQIAYGOUAADYGCUAADYWFuoAALqX\nhx4qVkNDY6hjtBEV1VfLly8JdYxex9ISqKys1OLFi+X3+zV//nwtWdL2L/Cee+7Rpk2b1K9fP61Z\ns0apqalWRgIQRENDo+LiCkMdow2PpzDUEXoly3YH+f1+LVq0SJWVlaqtrVVZWZnq6upajdm4caN2\n796t+vp6PfPMM1q4cKFVcULK43GHOsKP0pPz9+TsEvlDze12hzqC5SwrgZqaGsXHxysuLk7h4eHK\nzc1VRUVFqzHr16/XnDlzJEkTJkzQoUOH9MUXX1gVKWR6+n+Enpy/J2eXyB9qdigBy3YH+Xw+xcbG\nBqadTqf+9re/BR2zf/9+DR061KpYAHq5c3lM4x//cJ+z3VDd9ZiGZSXgcDg6NM4Yc1bzAUB7zuUx\nDY+n8Jwuq1syFnn//fdNVlZWYHrZsmVm+fLlrcYUFBSYsrKywPSoUaNMQ0NDm2WNHDnSSOLBgwcP\nHp14jBw5Mui62rItgfT0dNXX18vj8Sg6Olrl5eUqKytrNSYnJ0clJSXKzc3VBx98oIsvvrjdXUG7\nd++2KiYA2JplJRAWFqaSkhJlZWXJ7/crPz9fLpdLpaWlkqSCggJNmTJFGzduVHx8vPr376/Vq1db\nFQcA0A6HMafslAcA2Ea3vmxEZWWlEhMTlZCQoOLi4lDH6bR58+Zp6NChSkpKCnWUTvN6vbr22ms1\nZswYjR07VitWrAh1pE757rvvNGHCBI0bN06jR4/Www8/HOpIZ8Xv9ys1NVVTp04NdZROi4uLU3Jy\nslJTU/XTn/401HE65dChQ5o5c6ZcLpdGjx6tDz74INSROmzXrl1KTU0NPCIjI8/8//esj/xarLm5\n2YwcOdLs27fPNDU1mZSUFFNbWxvqWJ2yefNms23bNjN27NhQR+m0zz//3Gzfvt0YY8zhw4fN5Zdf\n3uN+/0ePHjXGGHPixAkzYcIEs2XLlhAn6rynnnrK3HrrrWbq1KmhjtJpcXFx5sCBA6GOcVZmz55t\nnn32WWPM9/9+Dh06FOJEZ8fv95uoqCjz6aefnnZMt90S6MjJZt3dVVddpYEDB4Y6xlmJiorSuHHj\nJEkDBgyQy+XSZ599FuJUndOvXz9JUlNTk/x+vwYNGhTiRJ2zf/9+bdy4UfPnz2/zVeqeoifm/uab\nb7RlyxbNmzdP0vfHNyMjI0Oc6uxUVVVp5MiRrc7HOlW3LYH2TiTz+XwhTGRfHo9H27dv14QJE0Id\npVNaWlo0btw4DR06VNdee61Gjx4d6kidcu+99+rJJ59Unz7d9r/pGTkcDl1//fVKT0/XH/7wh1DH\n6bB9+/ZpyJAhmjt3rsaPH68FCxbo2LFjoY51VtauXatbb731jGO67b8uThrrHo4cOaKZM2fq6aef\n1oABA0Idp1P69Omjf/zjH9q/f782b97coy4BsGHDBl166aVKTU3tkZ+mJWnr1q3avn27Nm3apFWr\nVmnLli2hjtQhzc3N2rZtm+6++25t27ZN/fv31/Lly0Mdq9Oampr0+uuv65ZbbjnjuG5bAjExMfJ6\nvYFpr9crp9MZwkT2c+LECc2YMUO33367pk2bFuo4Zy0yMlI///nP9fe//z3UUTqsurpa69ev1/Dh\nwzVr1iy98847mj17dqhjdcpll10mSRoyZIhuvvlm1dTUhDhRxzidTjmdTmVkZEiSZs6cqW3btoU4\nVedt2rRJaWlpGjJkyBnHddsSOPlks6amJpWXlysnJyfUsWzDGKP8/HyNHj1aixcvDnWcTvvqq690\n6NAhSVJjY6P++te/9qjLlC9btkxer1f79u3T2rVrdd111+m5554LdawOO3bsmA4fPixJOnr0qN56\n660e8y25qKgoxcbG6pNPPpH0/X71MWPGhDhV55WVlWnWrFlBx3Xbm8qc7mSznmTWrFl67733dODA\nAcXGxuqxxx7T3LlzQx2rQ7Zu3aoXXngh8BU/Sfr1r3+tG2+8McTJOubzzz/XnDlz1NLSopaWFt1x\nxx2aPHlyqGOdtZ62e/SLL77QzTffLOn73Su33XabbrjhhhCn6riVK1fqtttuU1NTk0aOHNnjTmQ9\nevSoqqqqOnQshpPFAMDGuu3uIACA9SgBALAxSgAAbIwSAAAbowQAwMYoAQCwMUoAvc7DDz8st9ut\n1157zbLT/SsqKlRXV9fpcUuXLtXbb79tSSbgbFAC6HVqamo0ceJEvffee7r66qst+Rmvvvqqamtr\nOz2uqKioR5+0hl6oSy5qDXSBBx54wCQnJ5uIiAgzbtw4ExERYZKTk83jjz/eZuycOXPMwoULzcSJ\nE82IESPMu+++a2bPnm1cLpfJy8sLjOvfv3/g+csvv2zy8vJMdXW1GTRokBk+fLhJTU01e/bsMc88\n84zJyMgwKSkpZsaMGebYsWNm69atbcbNmTPHrFu3zhhjTFVVlUlNTTVJSUlm3rx55vjx48YYY37y\nk5+YpUuXmvHjx5ukpCSzc+dOY4wxbrfbjBs3zowbN86kpqaaw4cPW/nrhE1QAuhVPvzwQ3PPPfeY\nEydOmEmTJp12XF5enpk1a5YxxpiKigoTERFh/vWvf5mWlhaTlpZmduzYYYwxZsCAAYF51q1bFyiI\nvLw888orrwTeO/nmKY888ohZuXJlu+N+mG5sbDSxsbGmvr7eGPP9TUx+85vfGGO+vxlLSUmJMcaY\n3/72t2b+/PnGGGOmTp1qqqurjTHf3zCnubn5bH5FQCvsDkKv8tFHHyk5OVl1dXVBrzX1wy0bx44d\nq6ioKI0ZM0YOh0NjxoyRx+MJ+rPMSVdc+ec//6mrrrpKycnJevHFF1vtAjKnXJnFGKNdu3Zp+PDh\nio+PlyTNmTNHmzdvDoyZPn26JGn8+PGBLJMmTdK9996rlStX6uDBgzrvvPOCZgSC6bYXkAM6Y8eO\nHcrLy9P+/fs1ePBgHTt2TMYYjR8/XtXV1brwwgvbzHP++edL+v6+AxdccEHg9T59+qi5uVlS6wu3\nNTY2tpr/5Pfy8vK0fv16JSUl6c9//nOrexe0d/G3U18zxrR67Yc85513XiDLkiVLdNNNN+mNN97Q\npEmT9Oabb2rUqFFn/sUAQbAlgF4hJSVF27dv1+WXX666ujpdd911euutt7Rt27Z2C6Cjhg4dqp07\nd6qlpUWvvvpqYEUdERGhb7/9NjDuyJEjioqK0okTJ/TCCy+cdpz0fQGMGjVKHo9He/bskSQ9//zz\nuuaaa86YZc+ePRozZowefPBBZWRkaNeuXWf95wJ+QAmg1/jvf/8buI/wzp07lZiYeMbxJ3/yPt2l\nmpcvX66bbrpJkyZNUnR0dOD13NxcPfnkk0pLS9PevXv1+OOPa8KECbryyitb7YY6ddwPLrjgAq1e\nvVq33HKLkpOTFRYWprvuuqvdLD9MP/3000pKSlJKSorOP/98ZWdnd+TXApwRl5IGABtjSwAAbIwS\nAAAbowQAwMYoAQCwMUoAAGyMEgAAG6MEAMDGKAEAsLH/A3nMVJQ2pttfAAAAAElFTkSuQmCC\n",
       "text": [
        "<matplotlib.figure.Figure at 0x859b198>"
       ]
      }
     ],
     "prompt_number": 58
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "def ratchet_click(s,U,N):\n",
      "    x = x_free.copy()\n",
      "    t = 0\n",
      "    \n",
      "    M = mutation_matrix(U)\n",
      "    S = selection_matrix(s)\n",
      "    E = M.dot(S)\n",
      "    \n",
      "    while x[0] > 0:\n",
      "        t += 1\n",
      "        x = E.dot(x)    \n",
      "        x /= x.sum()\n",
      "        assert allclose(x.sum(),1)\n",
      "        \n",
      "        n = multinomial(N, x)\n",
      "        assert len(n) == G\n",
      "        x = n/float(N)\n",
      "    \n",
      "    return t\n",
      "ratchet_click(s,U,N)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "pyout",
       "prompt_number": 60,
       "text": [
        "901"
       ]
      }
     ],
     "prompt_number": 60
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "clicks = array([ratchet_click(s,U,N) for i in range(100)])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 14
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "hist(clicks, normed=True)\n",
      "xlabel(\"time to click\")\n",
      "ylabel(\"frequency\");"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [
      {
       "metadata": {},
       "output_type": "display_data",
       "png": "iVBORw0KGgoAAAANSUhEUgAAAZ4AAAEPCAYAAAByRqLpAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3X10U/X9B/D3xVYnpQg6SCHBE5qmawttWumM4tRIbbuC\nVtANiwO7ARs6HYpbAT2bK2cOWpVtKjhww4qK+HikHZQOHwgqUHvAKmfWCUjRNn1Q+yBQsQ/p5/cH\nIz9CHxJK7jelfb/OyTm96ff7vZ/7bZJ3b+7NjSYiAiIiIkWGBLsAIiIaXBg8RESkFIOHiIiUYvAQ\nEZFSDB4iIlKKwUNERErpGjwlJSWIiYmB1WpFfn5+t20WLlwIq9UKm82G8vJyn31feeUVTJgwAeed\ndx4++OADr7FWrFgBq9WKmJgYbNu2TZ+NIiKisyM66ejoEIvFIpWVldLW1iY2m00qKiq82mzZskUy\nMjJERKS0tFTsdrvPvp988ol8+umn4nA4ZO/evZ6xPv74Y7HZbNLW1iaVlZVisVjE7XbrtXlERNRH\nuu3xlJWVISoqCmazGaGhocjKykJhYaFXm6KiImRnZwMA7HY7mpubUVdX12vfmJgYREdHd1lfYWEh\nZs2ahdDQUJjNZkRFRaGsrEyvzSMioj7SLXhcLhfGjRvnWTaZTHC5XH61qamp8dn3dDU1NTCZTGfU\nh4iI1NMteDRN86ud6HjFHn9rICIidUL0GthoNKKqqsqzXFVV5bVH0l2b6upqmEwmtLe3++zra33V\n1dUwGo1d2kVFReGzzz474+0hIhrMLBYLDh48GJjB9Dp41N7eLpGRkVJZWSmtra0+Ty7YvXu35+QC\nf/o6HA7Zs2ePZ/nkyQWtra1y6NAhiYyMlM7Ozi516bjJffbHP/4x2CV00R9rEumfdbEm/7Am//XH\nugL52qnbHk9ISAhWrVqF9PR0uN1uzJs3D7GxsVi7di0AYMGCBZg6dSqKi4sRFRWFsLAwFBQU9NoX\nAF5//XUsXLgQX3/9NaZNm4akpCRs3boVcXFxmDlzJuLi4hASEoInn3ySb7UREfVDugUPAGRkZCAj\nI8PrvgULFngtr1q1yu++ADBjxgzMmDGj2z4PPPAAHnjggT5WS0REKvDKBf2Aw+EIdgld9MeagP5Z\nF2vyD2vyX3+tK1C0/713N2homqbrmXRERANRIF87ucdDRERKMXiIiEgpBg8RESnF4CEiIqUYPERE\npBSDh4iIlGLwEBGRUgweIiJSisFzjhs+/GJomqb0Nnz4xcHebCI6h/HKBee4ExdCVb09A2sOicg3\nXrmAiIjOWQweIiJSisFDRERKMXiIiEgpBg8RESnF4CEiIqUYPEREpBSDh4iIlGLwEBGRUgweIiJS\nisFDRERKMXiIiEgpBg8RESnF4CEiIqUYPEREpBSDh4iIlGLwEBGRUgweIiJSisFDRERKMXiIiEgp\nBg8RESnF4CEiIqUYPEREpBSDh4iIlNI1eEpKShATEwOr1Yr8/Pxu2yxcuBBWqxU2mw3l5eU++zY2\nNiI1NRXR0dFIS0tDc3MzAOC7777DrFmzkJCQgLi4OOTl5em5aURE1Ee6BY/b7cbdd9+NkpISVFRU\nYOPGjfjkk0+82hQXF+PgwYM4cOAAnnrqKdx5550+++bl5SE1NRX79+9HSkqKJ2BefPFFAMC+ffuw\nd+9erF27Fl988YVem0dERH2kW/CUlZUhKioKZrMZoaGhyMrKQmFhoVeboqIiZGdnAwDsdjuam5tR\nV1fXa99T+2RnZ2PTpk0AgDFjxqClpQVutxstLS04//zzMXz4cL02j4iI+ki34HG5XBg3bpxn2WQy\nweVy+dWmpqamx7719fUwGAwAAIPBgPr6egBAeno6hg8fjjFjxsBsNiMnJwcjRozQa/OIiKiPQvQa\nWNM0v9qJiF9tuhtP0zTP/c8//zyOHz+O2tpaNDY24uqrr0ZKSgrGjx/fpV9ubq7nZ4fDAYfD4Vet\nRESDhdPphNPp1GVs3YLHaDSiqqrKs1xVVQWTydRrm+rqaphMJrS3t3e532g0Ajixl1NXV4eIiAjU\n1tZi9OjRAIBdu3ZhxowZOO+88zBq1ChcddVV2LNnj8/gISKirk7/p3zZsmUBG1u3t9qSk5Nx4MAB\nHD58GG1tbXjppZeQmZnp1SYzMxPPPvssAKC0tBQjRoyAwWDotW9mZibWr18PAFi/fj2mT58OAIiJ\nicHbb78NAGhpaUFpaSliY2P12jwiIuoj3fZ4QkJCsGrVKqSnp8PtdmPevHmIjY3F2rVrAQALFizA\n1KlTUVxcjKioKISFhaGgoKDXvgCwdOlSzJw5E+vWrYPZbMbLL7/sGW/evHmIj49HZ2cn5s6di4kT\nJ+q1eURE1Eea+HOQZQDRNM2v40rnihPHuFRvz8CaQyLyLZCvnbxyARERKcXgISIipRg8RESkFIOH\niIiUYvAQEZFSDB4iIlKKwUNEREoxeIiISCkGDxERKcXgISIipRg8RESkFIOHiIiUYvAQEZFSDB4i\nIlKKwUNEREoxeIiISCkGDxERKcXgISIipRg8RESkFIOHiIiUYvAQEZFSDB4iIlKKwUNEREoxeIiI\nSCkGDxERKcXgISIipRg8RESkFIOHiIiUYvAQEZFSDB4iIlKKwUNEREoxeIiISCkGDxERKcXgISIi\npRg8RESklK7BU1JSgpiYGFitVuTn53fbZuHChbBarbDZbCgvL/fZt7GxEampqYiOjkZaWhqam5s9\nv9u3bx+uvPJKTJw4EQkJCWhtbdVv44iIqG9EJx0dHWKxWKSyslLa2trEZrNJRUWFV5stW7ZIRkaG\niIiUlpaK3W732TcnJ0fy8/NFRCQvL0+WLFkiIiLt7e2SkJAg+/btExGRxsZGcbvdXerScZODAoAA\novg2sOaQiHwL5PNetz2esrIyREVFwWw2IzQ0FFlZWSgsLPRqU1RUhOzsbACA3W5Hc3Mz6urqeu17\nap/s7Gxs2rQJALBt2zYkJCQgPj4eADBy5EgMGcJ3EomI+hvdXpldLhfGjRvnWTaZTHC5XH61qamp\n6bFvfX09DAYDAMBgMKC+vh4AsH//fmiahh//+MeYNGkSHnnkEb02jYiIzkKIXgNrmuZXuxN7cL7b\ndDeepmme+zs6OvDee+9hz549uPDCC5GSkoJJkyZhypQpZ1Y4ERHpymfwTJo0CXPnzsVtt92GkSNH\n+j2w0WhEVVWVZ7mqqgomk6nXNtXV1TCZTGhvb+9yv9FoBHBiL6eurg4RERGora3F6NGjAQDjxo3D\nNddcg4svvhgAMHXqVHzwwQfdBk9ubq7nZ4fDAYfD4fd2ERENBk6nE06nU5/BfR0E2r9/v9x///1i\nsVjk1ltvlZKSEuns7PR58Ki9vV0iIyOlsrJSWltbfZ5csHv3bs/JBb31zcnJkby8PBERWbFihefk\ngsbGRrnsssvk22+/lfb2drn++uuluLi4S11+bPI5BTy5gIgUCOTz3u+R3G63FBYWytixY8VkMsmD\nDz4oDQ0NvfYpLi6W6OhosVgssnz5chERWbNmjaxZs8bT5q677hKLxSIJCQmyd+/eXvuKiDQ0NEhK\nSopYrVZJTU2VpqYmz++ef/55mTBhgkycONETSF02eIC9aDJ4iEiFQD7vtf8N2KuPPvoIBQUF2Lp1\nK9LT03Hbbbfhvffew/PPP48PP/xQn10xnWia5tdxpXPFiWNcqrdnYM0hEfkWyNdOv47xXHTRRZg/\nfz7y8/NxwQUXAACuuOIK7Ny5MyBFEBHR4OFzj+fQoUOIjIxUVY/u9NzjaWtrw5EjR3QZuyejRo0C\n93iISG9K93j++c9/YvHixRgxYgQAoKmpCStXrsRDDz0UkAIGkp/+NBtbt27GeeddoGR9fPEnonOR\nzz2exMTELsdxkpKSvK6rdi7Rc4/H4bgJO3bMBXCTLuN35caJ/x24x0NE+grka6fPKxd0dnbiu+++\n8ywfP34cbW1tAVk5ERENPj7favvZz36GlJQUzJ07FyKCgoIC3H777SpqIyKiAchn8CxZsgQJCQl4\n8803oWkaHnzwQaSnp6uojYiIBiC/rtWWkZGBjIwMvWshIqJBwOcxntdeew1WqxXDhw9HeHg4wsPD\nMXz4cBW1ERHRAORzj2fx4sXYvHkzYmNjVdRDREQDnM89noiICIYOEREFjM89nuTkZNx6662YPn06\nzj//fAAnzue++eabdS+OiIgGHp/B88033+DCCy/Etm3bvO5n8BARUV/4DJ5nnnlGQRlERDRY+DzG\n8+mnnyIlJQUTJkwAAOzbt4/XaSMioj7zGTy//OUvsXz5cs/xnfj4eGzcuFH3woiIaGDyGTzffvst\n7Ha7Z1nTNISGhupaFBERDVw+g2fUqFE4ePCgZ/nVV1/FmDFjdC2KiIgGLp8nF6xatQq/+tWv8N//\n/hdjx47F+PHjsWHDBhW1ERHRAOQzeCwWC9566y20tLSgs7MT4eHhKuoiIqIBymfwLFu2zPMFQJqm\nee5/8MEHdS2MiIgGJp/BExYW5gmc48ePY/PmzYiLi9O9MCIiGph8Bs/vfvc7r+WcnBykpaXpVhAR\nEQ1sPs9qO11LSwtcLpcetRAR0SDgc48nPj7e83NnZye+/PJLHt8hIqI+8xk8//rXv/6/cUgIDAYD\nP0BKRER95jN4Tv+20aNHj3otX3zxxYGtiM4BIV5nOKoQHj4SR440Kl0nEenDZ/Bcdtll+OKLLzBy\n5EgAQFNTEy699FJomgZN03Do0CHdi6T+pgOAKF3j0aNqg46I9OPz5ILU1FRs3rwZDQ0NaGhowJYt\nW5CWlobKykqGDhERnTGfwbN7925MnTrVs5yRkYFdu3bpWhQREQ1cPt9qGzt2LB566CHMnj0bIoIX\nXngBRqNRRW1ERDQA+dzj2bhxI7788kvMmDEDN998M7788kt+Hw8REfWZzz2eSy65BI8//jhaWloQ\nFhamoiYiIhrAfO7x7Nq1C3FxcYiJiQEAfPTRR/j1r3+te2FERDQw+Qyee++9FyUlJfj+978PALDZ\nbNixY4fuhRER0cDk17XaLr30Uq/lkBCf79ARERF1y2fwXHrppdi5cycAoK2tDY8++ihiY2P9Gryk\npAQxMTGwWq3Iz8/vts3ChQthtVphs9lQXl7us29jYyNSU1MRHR2NtLQ0NDc3e433xRdfYNiwYVi5\ncqVfNRIRkVo+g2fNmjVYvXo1XC4XjEYjysvLsXr1ap8Du91u3H333SgpKUFFRQU2btyITz75xKtN\ncXExDh48iAMHDuCpp57CnXfe6bNvXl4eUlNTsX//fqSkpCAvL89rzPvuuw/Tpk3zewKIiEitXt8z\n6+jowD333IMXXnjhjAcuKytDVFQUzGYzACArKwuFhYVee0tFRUXIzs4GANjtdjQ3N6Ourg6VlZU9\n9i0qKvIcY8rOzobD4fCEz6ZNmxAZGcmz74iI+rFe93hCQkLw+eefo7W19YwHdrlcGDdunGfZZDJ1\n+R6fntrU1NT02Le+vh4GgwEAYDAYUF9fDwA4duwYHn74YeTm5p5xrUREpI7PswQiIyPxox/9CJmZ\nmRg6dCgAQNM03Hfffb328/fqxSK+LzYpIt2Od/JCpQCQm5uLRYsWYejQoX6NSUREwdFj8MyZMwfP\nPfccioqKsGjRInR2duLYsWN+D2w0GlFVVeVZrqqqgslk6rVNdXU1TCYT2tvbu9x/8jI9BoMBdXV1\niIiIQG1tLUaPHg3gxFt7r732GhYvXozm5mYMGTIEF154YbefOTp1r8jhcMDhcPi9XUREg4HT6YTT\n6dRncOlBbGysuFwuiY+Pl4aGBvn666+9br60t7dLZGSkVFZWSmtrq9hsNqmoqPBqs2XLFsnIyBAR\nkd27d4vdbvfZNycnR/Ly8kREZMWKFbJkyZIu687NzZWVK1d2W1cvm3zWrr02U4BNAoiiW4cAULi+\nk7fgrJOIgieQz8Ee93juuOMOpKSk4NChQ5g0aZLX7/z5Hp6QkBCsWrUK6enpcLvdmDdvHmJjY7F2\n7VoAwIIFCzB16lQUFxcjKioKYWFhKCgo6LUvACxduhQzZ87EunXrYDab8fLLL/c1c4mIKAi0/yVZ\nj+644w6sWbNGVT260zRNt2NADsdN2LFjLoCbdBm/KzdOvFuq+piWFpR16vV3IyLfAvna6dfneIiI\niALFr0vmEBERBQqDh4iIlGLwEBGRUgweIiJSisFDRERKMXiIiEgpBg8RESnF4CEiIqUYPEREpBSD\nh4iIlGLwEBGRUgweIiJSisFDRERKMXiIiEipHr8Ijqh/CYGmaUrXGB4+EkeONCpdJ9FgwOChc0QH\nVH/53NGjaoOOaLDgW21ERKQUg4eIiJRi8BARkVIMHiIiUorBQ0RESjF4iIhIKQYPEREpxeAhIiKl\nGDxERKQUg4eIiJRi8BARkVIMHiIiUorBQ0RESjF4iIhIKQYPEREpxeAhIiKlGDxERKQUg4eIiJRi\n8BARkVK6B09JSQliYmJgtVqRn5/fbZuFCxfCarXCZrOhvLzcZ9/GxkakpqYiOjoaaWlpaG5uBgC8\n8cYbSE5ORkJCApKTk7F9+3Z9N46IiM6c6Kijo0MsFotUVlZKW1ub2Gw2qaio8GqzZcsWycjIEBGR\n0tJSsdvtPvvm5ORIfn6+iIjk5eXJkiVLRESkvLxcamtrRUTkP//5jxiNxi416bnJ116bKcAmAUTR\nrUMAKFzfydvgWScRnRDI54OuezxlZWWIioqC2WxGaGgosrKyUFhY6NWmqKgI2dnZAAC73Y7m5mbU\n1dX12vfUPtnZ2di0aRMAIDExEREREQCAuLg4HD9+HO3t7XpuIhERnSFdg8flcmHcuHGeZZPJBJfL\n5VebmpqaHvvW19fDYDAAAAwGA+rr67us+7XXXsOkSZMQGhoa0G0iIqKzE6Ln4Jqm+dXuxF6c7zbd\njadpWpf7P/74YyxduhRvvPFGt2Pl5uZ6fnY4HHA4HH7VSUQ0WDidTjidTl3G1jV4jEYjqqqqPMtV\nVVUwmUy9tqmurobJZEJ7e3uX+41GI4ATezl1dXWIiIhAbW0tRo8e7dXu5ptvxnPPPYfx48d3W9ep\nwUNERF2d/k/5smXLAja2rm+1JScn48CBAzh8+DDa2trw0ksvITMz06tNZmYmnn32WQBAaWkpRowY\nAYPB0GvfzMxMrF+/HgCwfv16TJ8+HQDQ3NyMadOmIT8/H1deeaWem0ZERH0VsNMUelBcXCzR0dFi\nsVhk+fLlIiKyZs0aWbNmjafNXXfdJRaLRRISEmTv3r299hURaWhokJSUFLFarZKamipNTU0iIvKn\nP/1JwsLCJDEx0XP76quvvOrRc5N5VtvAWycRnRDI54P2vwEHDU3T/Dqm1BcOx03YsWMugJt0Gb8r\nN068W6r6T6gNmnUOsqcHUY8C+drJKxcQEZFSDB4iIlKKwUNEREoxeIiISCkGDxERKaXrB0iJzm0h\nfl99IxDCw0fiyJFGZesjChYGD1GPOqDyFO6jR9WFHFEw8a02IiJSisFDRERKMXiIiEgpBg8RESnF\n4CEiIqUYPEREpBSDh4iIlGLwEBGRUgweIiJSisFDRERK8ZI5RP2G2mvDAbw+HAUHg4eo31B7bTiA\n14ej4OBbbUREpBSDh4iIlGLwEBGRUgweIiJSisFDRERKMXiIiEgpBg8RESnF4CEiIqUYPEREpBSD\nh4iIlGLwEBGRUgweIiJSisFDRERKMXiIiEgpBg8RESnF7+MhogFt+PCLcfRok+K1hgJoV7rGc+lL\n/XTd4ykpKUFMTAysVivy8/O7bbNw4UJYrVbYbDaUl5f77NvY2IjU1FRER0cjLS0Nzc3Nnt+tWLEC\nVqsVMTEx2LZtm34bRkTnjBOhI4pv7crXqT5cz4LopKOjQywWi1RWVkpbW5vYbDapqKjwarNlyxbJ\nyMgQEZHS0lKx2+0+++bk5Eh+fr6IiOTl5cmSJUtEROTjjz8Wm80mbW1tUllZKRaLRdxud5e6dNxk\nufbaTAE2CSBneNvehz4iQMf/HnV96Xs2Nem1zt5uJ9fZ17k6m3Xq9ffTc179rQm6PR9Ot337dmXr\nOlXv86rX4+ls/5Z9qUvfv2Ugx9dtj6esrAxRUVEwm80IDQ1FVlYWCgsLvdoUFRUhOzsbAGC329Hc\n3Iy6urpe+57aJzs7G5s2bQIAFBYWYtasWQgNDYXZbEZUVBTKysr02rwAcwa7gG44g11AD5zBLqAb\nzmAX0A1nsAvowul0BruEbjiDXUAPnMEuQFe6BY/L5cK4ceM8yyaTCS6Xy682NTU1Pfatr6+HwWAA\nABgMBtTX1wMAampqYDKZel0fEREFn24nF2ia5le7E3twvtt0N56mab2ux98aAiU0dAiGDn0IISH/\nPKN+3333Kb73vb19WKPgyJE+dCMiCiLdgsdoNKKqqsqzXFVV5bVH0l2b6upqmEwmtLe3d7nfaDQC\nOLGXU1dXh4iICNTW1mL06NE9jnWyz6ksFovyQPJHW9uBs+it1/YsC8I6e3Nynb3Vpdc6fQlUTYGc\nV/9qUvl8WLZM5d/uVL1to141ne28nnldev4tLRZL4AYL2NGi07S3t0tkZKRUVlZKa2urz5MLdu/e\n7Tm5oLe+OTk5kpeXJyIiK1as6HJyQWtrqxw6dEgiIyOls7NTr80jIqI+0m2PJyQkBKtWrUJ6ejrc\nbjfmzZuH2NhYrF27FgCwYMECTJ06FcXFxYiKikJYWBgKCgp67QsAS5cuxcyZM7Fu3TqYzWa8/PLL\nAIC4uDjMnDkTcXFxCAkJwZNPPtkv92yIiAY7TcSPgyxEREQBMqgumePPB1r1YDabkZCQgKSkJFx+\n+eUAgvNB2Llz58JgMCA+Pt5zX1/q2Lt3L+Lj42G1WnHPPfcEvKbc3FyYTCYkJSUhKSkJW7duVVpT\nVVUVrrvuOkyYMAETJ07E448/DiC4c9VTTcGcq++++w52ux2JiYmIi4vD/fffDyC489RTTcF+TJ3k\ndruRlJSEG2+8EUDwn3/d1aRkroL9Xp8q/nygVS9ms1kaGhq87jvbD8L2xTvvvCMffPCBTJw4sU91\nnDxm9sMf/lDef/99ERHJyMiQrVu3BrSm3NxcWblyZZe2qmqqra2V8vJyERE5evSoREdHS0VFRVDn\nqqeagj1XLS0tInLiuKzdbpd333036I+p7moK9jydtHLlSrntttvkxhtvFJHgP/+6q0nFXA2aPR5/\nPtCqJzntHc1gfBD26quvxsiRI/tcx/vvv4/a2locPXrUs+d2++23e/oEqiag+9PsVdUUERGBxMRE\nAMCwYcMQGxsLl8sV1LnqqSYguHM1dOhQAEBbWxvcbjdGjhwZ9MdUdzUBwZ0n4MSZtsXFxZg/f76n\nlmDPVXc1iYjuczVogsefD7TqRdM0XH/99UhOTsY//vEPAP3ng7BnWsfp9xuNRl3qe+KJJ2Cz2TBv\n3jzP2w/BqOnw4cMoLy+H3W7vN3N1sqYrrrgCQHDnqrOzE4mJiTAYDJ63AoM9T93VBAT/MbVo0SI8\n8sgjGDLk/192gz1X3dWkaZruczVogieYZ7jt3LkT5eXl2Lp1K1avXo13333X6/f95YOwvupQ5c47\n70RlZSU+/PBDjBkzBr/97W+DUsexY8dwyy234LHHHkN4eLjX74I1V8eOHcNPfvITPPbYYxg2bFjQ\n52rIkCH48MMPUV1djXfeeQfbt2/3+n0w5un0mpxOZ9DnafPmzRg9ejSSkpJ6/NC86rnqqSYVczVo\ngsefD7TqZcyYMQCAUaNGYcaMGSgrK/N8EBZAnz4IGyhnUofJZILRaER1dbWu9Y0ePdrzJJw/f77n\nrUaVNbW3t+OWW27BnDlzMH36dADBn6uTNc2ePdtTU3+YKwC46KKLMG3aNOzduzfo83R6TXv27An6\nPO3atQtFRUUYP348Zs2ahbfffhtz5swJ6lx1V9Ptt9+uZq7O6qjUOcSfD7TqoaWlRY4cOSIiIseO\nHZPJkyfLv//976B9ELaysrLLyQVnWsfll18upaWl0tnZGZCDm6fXVFNT4/n5L3/5i8yaNUtpTZ2d\nnTJnzhy59957ve4P5lz1VFMw5+qrr76SpqYmERH59ttv5eqrr5Y333wzqPPUU021tbWeNsF4TJ3K\n6XTKDTfcICL94/l3ek0qHlODJnhERIqLiyU6OlosFossX75cyToPHTokNptNbDabTJgwwbPehoYG\nSUlJEavVKqmpqZ4ni4jIn//8Z7FYLPKDH/xASkpKAlZLVlaWjBkzRkJDQ8VkMsnTTz/dpzr27Nkj\nEydOFIvFIr/5zW8CWtO6detkzpw5Eh8fLwkJCXLTTTdJXV2d0preffdd0TRNbDabJCYmSmJiomzd\nujWoc9VdTcXFxUGdq3379klSUpLYbDaJj4+Xhx9+WET69tjWu6ZgP6ZO5XQ6PWeQBfv5d9L27ds9\nNc2ePVv3ueIHSImISKlBc4yHiIj6BwYPEREpxeAhIiKlGDxERKQUg4eIiJRi8BARkVIMHqLTfPPN\nN/j73//uWa6pqcFPf/rTgK9nx44d2L17d0DGMpvNaGxsBABcddVVvbZ1OBzYu3dvQNZL1BcMHqLT\nNDU14cknn/Qsjx07Fq+88krA17N9+3bs2rUrIGOdeo2vnTt3+mzbH67JR4MXg4foNEuXLsVnn32G\npKQkLFmyBJ9//rnni+qeeeYZTJ8+HWlpaRg/fjxWrVqFRx99FJdddhmuvPJKNDU1AQA+++wzZGRk\nIDk5Gddccw0+/fRTr3UcPnwYa9euxV//+lckJSVh586dOHz4MKZMmQKbzYbrr7/e67pYJx07dgy/\n+MUvkJCQAJvNhtdff71Lm2HDhnl+zs/PR0JCAhITE/HAAw94tevs7MTPf/5z/OEPfzjrOSM6IwG5\n3gLRAHL48GGva8edei25goICiYqKkmPHjslXX30lw4cPl7Vr14qIyKJFi+Rvf/ubiIhMmTJFDhw4\nICIipaWlMmXKlC7rOf0Lt2644QZ59tlnRUTk6aeflunTp3fps3jxYlm0aJFn+eQlVk79ssFhw4aJ\nyIlLRE1of23vAAAB90lEQVSePFmOHz/u1dbhcEhpaalkZWUpu3QU0alCgh18RP2N+LiK1HXXXYew\nsDCEhYVhxIgRnq8Mjo+Px759+9DS0oJdu3Z5HRdqa2vzua7S0lLPF2jNnj0bixcv7tL+rbfewksv\nveRZHjFiRI91vvnmm5g7dy6+973vebUVESxYsAC33nqr56uhiVRi8BCdoQsuuMDz85AhQzzLQ4YM\nQUdHBzo7OzFy5EiUl5ef8di+Qs/fNsCJYzndtdU0DZMnT8bbb7+N++67z2t7iFTgMR6i04SHh+Po\n0aNn3O/ki3x4eDjGjx+PV1991XP/vn37fK5n8uTJePHFFwEAGzZswDXXXNOlT2pqKlavXu1ZPvnt\nkN1JTU1FQUEBjh8/DgCe408AMH/+fEydOhUzZ86E2+0+k80kOmsMHqLTXHLJJbjqqqsQHx+PJUuW\neJ0FdvoZYaf/fHJ5w4YNWLduHRITEzFx4kQUFRV1Wc+NN96I119/3XNywRNPPIGCggLYbDZs2LAB\njz32WJc+v//979HU1IT4+HgkJibC6XR2aXOyhvT0dGRmZiI5ORlJSUlYuXKlV7tFixYhKSkJc+bM\n8XsviigQ+LUIRESkFPd4iIhIKQYPEREpxeAhIiKlGDxERKQUg4eIiJRi8BARkVIMHiIiUorBQ0RE\nSv0f06cKmjKS7kIAAAAASUVORK5CYII=\n",
       "text": [
        "<matplotlib.figure.Figure at 0x11b276d8>"
       ]
      }
     ],
     "prompt_number": 20
    },
    {
     "cell_type": "markdown",
     "metadata": {},
     "source": [
      "Now using Diffusion equation [@Gordo2000]:"
     ]
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "a = lambda x,x0,s: 0.6 * s * (1 - x/x0) * x\n",
      "b = lambda x,N: x * (1-x) / N"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "G = lambda x,x0,s,N: exp(2 * N * 0.6 * s * x * (x/2 - x0) / x0)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "from scipy.integrate import quad as integrate"
     ],
     "language": "python",
     "metadata": {},
     "outputs": [],
     "prompt_number": 36
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "intG = lambda z,w,x0,s,N: integrate(G, z, w, args=(x0,s,N))[0]"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "T0x0 = lambda x0: integrate( lambda x,x0,s,N: 2 * N / ( x * G(x,x0) ) * intG(0, x, x0, s, N), 0, x0, args=(x0, s, N))[0]\n",
      "Tx01 = lambda x0: integrate( lambda x,x0,s,N: 2 * N / ( x * G(x,x0) ) * intG(0, x0, x0, s, N), x0, 1, args=(x0, s, N))[0]\n",
      "T = lambda x0,s,N: T0x0(x0,s,N) + Tx01(x0,s,N)\n",
      "\n",
      "print T(1)"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "clicks = array([[[ratchet_click(s,U,N) for i in range(100)] for s in [0.001, 0.01, 0.1]] for U in [0.0003, 0.003, 0.03]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [
      "Ts = array([[[T(x_msb[0],s,N) for i in range(100)] for s in [0.001, 0.01, 0.1]] for U in [0.0003, 0.003, 0.03]])"
     ],
     "language": "python",
     "metadata": {},
     "outputs": []
    },
    {
     "cell_type": "code",
     "collapsed": false,
     "input": [],
     "language": "python",
     "metadata": {},
     "outputs": []
    }
   ],
   "metadata": {}
  }
 ]
}