chaotic-dynamics-of-social-behavior.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import ipywidgets as widgets\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import integrate\n",
    "\n",
    "import models\n",
    "import payoffs\n",
    "import plotting\n",
    "import selection_functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.style.use(\"ggplot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Basic simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "integrate.solve_ivp?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# random initial condition\n",
    "prng = np.random.RandomState(42)\n",
    "number_of_genotypes = 4\n",
    "initial_offspring_share, = prng.dirichlet(np.ones(number_of_genotypes), 1)\n",
    "y0 = initial_offspring_share\n",
    "\n",
    "# define the selection functions\n",
    "d1, d3 = 1.5, 1\n",
    "UGA = lambda x_A: selection_functions.kirkpatrick_selection(x_A, d1)\n",
    "UgA = lambda x_A: selection_functions.kirkpatrick_selection(x_A, d3)\n",
    "\n",
    "# define the payoffs\n",
    "payoff_kernel = payoffs.prisoners_dilemma_payoffs(prng)\n",
    "\n",
    "def f(t, y):\n",
    "    W = models.generalized_sexual_selection(y, UGA, UgA, payoff_kernel)\n",
    "    y_dot = models.offspring_genotypes_evolution(W, y)\n",
    "    return y_dot\n",
    "\n",
    "result = integrate.solve_ivp(f, t_span=(0, 10), y0=y0, rtol=1e-9, atol=1e-12, vectorized=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "  message: 'The solver successfully reached the interval end.'\n",
       "     nfev: 530\n",
       "     njev: 0\n",
       "      nlu: 0\n",
       "      sol: None\n",
       "   status: 0\n",
       "  success: True\n",
       "        t: array([  0.00000000e+00,   7.73035243e-03,   8.50338767e-02,\n",
       "         1.80589388e-01,   2.77094581e-01,   3.74725993e-01,\n",
       "         4.73415427e-01,   5.73101566e-01,   6.73726179e-01,\n",
       "         7.75234035e-01,   8.77572804e-01,   9.80692986e-01,\n",
       "         1.08454925e+00,   1.18912246e+00,   1.29436792e+00,\n",
       "         1.40024417e+00,   1.50671259e+00,   1.61373734e+00,\n",
       "         1.72128527e+00,   1.82932584e+00,   1.93783106e+00,\n",
       "         2.04677537e+00,   2.15613552e+00,   2.26589048e+00,\n",
       "         2.37602129e+00,   2.48651099e+00,   2.59734446e+00,\n",
       "         2.70850831e+00,   2.81999079e+00,   2.93178165e+00,\n",
       "         3.04387207e+00,   3.15625453e+00,   3.26892274e+00,\n",
       "         3.38187156e+00,   3.49509690e+00,   3.60859565e+00,\n",
       "         3.72236565e+00,   3.83640559e+00,   3.95071495e+00,\n",
       "         4.06529401e+00,   4.18014374e+00,   4.29526578e+00,\n",
       "         4.41066243e+00,   4.52633656e+00,   4.64229163e+00,\n",
       "         4.75853165e+00,   4.87506115e+00,   4.99188511e+00,\n",
       "         5.10900903e+00,   5.22643881e+00,   5.34418077e+00,\n",
       "         5.46224161e+00,   5.58062841e+00,   5.69934856e+00,\n",
       "         5.81840975e+00,   5.93781994e+00,   6.05758734e+00,\n",
       "         6.17772033e+00,   6.29822750e+00,   6.41911753e+00,\n",
       "         6.54039920e+00,   6.66208137e+00,   6.78417290e+00,\n",
       "         6.90668261e+00,   7.02961931e+00,   7.15299170e+00,\n",
       "         7.27680837e+00,   7.40107777e+00,   7.52580820e+00,\n",
       "         7.65100776e+00,   7.77668439e+00,   7.90284584e+00,\n",
       "         8.02949965e+00,   8.15665322e+00,   8.28431377e+00,\n",
       "         8.41248839e+00,   8.54118407e+00,   8.67040773e+00,\n",
       "         8.80016625e+00,   8.93046655e+00,   9.06131561e+00,\n",
       "         9.19272051e+00,   9.32468854e+00,   9.45722720e+00,\n",
       "         9.59034429e+00,   9.72404796e+00,   9.85834679e+00,\n",
       "         9.99324982e+00,   1.00000000e+01])\n",
       " t_events: None\n",
       "        y: array([[ 0.08219683,  0.08240963,  0.08440244,  0.08654122,  0.08836118,\n",
       "         0.08987964,  0.09110892,  0.0920622 ,  0.09275319,  0.093196  ,\n",
       "         0.09340494,  0.09339443,  0.0931788 ,  0.09277216,  0.09218835,\n",
       "         0.09144098,  0.09054328,  0.0895081 ,  0.08834782,  0.08707433,\n",
       "         0.08569903,  0.08423277,  0.08268586,  0.0810681 ,  0.07938869,\n",
       "         0.07765636,  0.07587927,  0.07406509,  0.07222098,  0.07035363,\n",
       "         0.06846928,  0.0665737 ,  0.06467227,  0.06276995,  0.06087132,\n",
       "         0.05898058,  0.05710162,  0.05523799,  0.0533929 ,  0.05156932,\n",
       "         0.04976991,  0.04799708,  0.046253  ,  0.04453959,  0.04285859,\n",
       "         0.04121149,  0.03959962,  0.03802411,  0.03648595,  0.03498592,\n",
       "         0.0335247 ,  0.0321028 ,  0.0307206 ,  0.02937838,  0.02807627,\n",
       "         0.02681431,  0.02559244,  0.02441052,  0.02326829,  0.02216544,\n",
       "         0.02110158,  0.02007623,  0.01908888,  0.01813895,  0.01722579,\n",
       "         0.01634874,  0.01550706,  0.01470001,  0.01392679,  0.01318658,\n",
       "         0.01247854,  0.0118018 ,  0.01115549,  0.01053871,  0.00995055,\n",
       "         0.00939012,  0.00885648,  0.00834874,  0.00786597,  0.00740727,\n",
       "         0.00697174,  0.00655848,  0.00616662,  0.00579529,  0.00544364,\n",
       "         0.00511083,  0.00479604,  0.00449848,  0.00448406],\n",
       "       [ 0.52725178,  0.52746882,  0.52969426,  0.53257623,  0.53562215,\n",
       "         0.53882962,  0.5421885 ,  0.54568852,  0.54931917,  0.55306988,\n",
       "         0.55693002,  0.56088904,  0.56493655,  0.56906328,  0.57325906,\n",
       "         0.57751402,  0.58181859,  0.58616349,  0.59053983,  0.59493909,\n",
       "         0.5993531 ,  0.60377411,  0.60819475,  0.61260806,  0.61700745,\n",
       "         0.62138672,  0.62574004,  0.63006196,  0.63434736,  0.63859149,\n",
       "         0.64278992,  0.64693852,  0.6510335 ,  0.65507134,  0.65904881,\n",
       "         0.66296295,  0.66681105,  0.67059068,  0.6742996 ,  0.67793584,\n",
       "         0.68149764,  0.68498343,  0.68839187,  0.69172179,  0.69497224,\n",
       "         0.69814241,  0.70123169,  0.70423963,  0.70716593,  0.71001045,\n",
       "         0.71277321,  0.71545434,  0.71805413,  0.720573  ,  0.72301146,\n",
       "         0.72537018,  0.72764991,  0.72985152,  0.73197597,  0.73402432,\n",
       "         0.7359977 ,  0.73789735,  0.73972456,  0.7414807 ,  0.74316719,\n",
       "         0.74478553,  0.74633725,  0.74782393,  0.74924719,  0.7506087 ,\n",
       "         0.75191014,  0.7531532 ,  0.75433964,  0.75547117,  0.75654954,\n",
       "         0.75757652,  0.75855384,  0.75948325,  0.76036648,  0.76120525,\n",
       "         0.76200125,  0.76275618,  0.76347168,  0.76414937,  0.76479085,\n",
       "         0.76539769,  0.76597141,  0.7665135 ,  0.76653975],\n",
       "       [ 0.2306407 ,  0.22972875,  0.22080735,  0.21026191,  0.20012819,\n",
       "         0.19037847,  0.181011  ,  0.17202198,  0.16340608,  0.15515673,\n",
       "         0.14726626,  0.13972614,  0.13252698,  0.12565733,  0.11910735,\n",
       "         0.11286673,  0.10692489,  0.10127102,  0.09589425,  0.09078367,\n",
       "         0.08592841,  0.08131772,  0.07694102,  0.07278794,  0.06884832,\n",
       "         0.0651123 ,  0.06157031,  0.05821308,  0.05503167,  0.05201747,\n",
       "         0.04916221,  0.04645794,  0.04389706,  0.04147228,  0.03917667,\n",
       "         0.03700359,  0.03494672,  0.03300003,  0.03115779,  0.02941457,\n",
       "         0.02776516,  0.02620467,  0.02472841,  0.02333195,  0.0220111 ,\n",
       "         0.02076185,  0.01958042,  0.01846325,  0.01740691,  0.0164082 ,\n",
       "         0.01546407,  0.01457161,  0.01372811,  0.01293095,  0.01217769,\n",
       "         0.01146599,  0.01079366,  0.0101586 ,  0.00955883,  0.00899248,\n",
       "         0.00845777,  0.00795302,  0.00747664,  0.0070271 ,  0.00660299,\n",
       "         0.00620293,  0.00582566,  0.00546993,  0.0051346 ,  0.00481857,\n",
       "         0.00452079,  0.00424029,  0.00397612,  0.00372739,  0.00349327,\n",
       "         0.00327295,  0.00306567,  0.00287071,  0.0026874 ,  0.00251507,\n",
       "         0.00235313,  0.00220098,  0.00205807,  0.00192388,  0.0017979 ,\n",
       "         0.00167968,  0.00156877,  0.00146473,  0.00145971],\n",
       "       [ 0.15991069,  0.1603928 ,  0.16509596,  0.17062064,  0.17588847,\n",
       "         0.18091228,  0.18569157,  0.19022731,  0.19452156,  0.1985774 ,\n",
       "         0.20239878,  0.20599039,  0.20935767,  0.21250723,  0.21544524,\n",
       "         0.21817826,  0.22071324,  0.22305739,  0.2252181 ,  0.22720292,\n",
       "         0.22901947,  0.2306754 ,  0.23217836,  0.23353591,  0.23475554,\n",
       "         0.23584462,  0.23681038,  0.23765988,  0.23839999,  0.2390374 ,\n",
       "         0.23957859,  0.24002983,  0.24039717,  0.24068642,  0.2409032 ,\n",
       "         0.24105288,  0.2411406 ,  0.24117131,  0.2411497 ,  0.24108027,\n",
       "         0.24096729,  0.24081482,  0.24062673,  0.24040666,  0.24015808,\n",
       "         0.23988426,  0.23958827,  0.23927301,  0.23894121,  0.23859542,\n",
       "         0.23823802,  0.23787125,  0.23749715,  0.23711767,  0.23673458,\n",
       "         0.23634952,  0.23596398,  0.23557936,  0.23519691,  0.23481776,\n",
       "         0.23444295,  0.23407339,  0.23370992,  0.23335325,  0.23300403,\n",
       "         0.2326628 ,  0.23233004,  0.23200613,  0.23169142,  0.23138615,\n",
       "         0.23109053,  0.2308047 ,  0.23052876,  0.23026273,  0.23000663,\n",
       "         0.22976042,  0.22952401,  0.2292973 ,  0.22908016,  0.22887241,\n",
       "         0.22867388,  0.22848436,  0.22830363,  0.22813146,  0.2279676 ,\n",
       "         0.22781179,  0.22766378,  0.2275233 ,  0.22751647]])"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_differences(selection_function, dG, dg, T, R, P, S, M, m, epsilon):\n",
    "\n",
    "    # random initial condition for first simulation\n",
    "    prng = np.random.RandomState(42)\n",
    "    \n",
    "    # create the selection functions\n",
    "    if selection_function == \"kirkpatrick\":\n",
    "        UGA = lambda x_A: selection_functions.kirkpatrick_selection(x_A, dG)\n",
    "        UgA = lambda x_A: selection_functions.kirkpatrick_selection(x_A, dg)\n",
    "    elif selection_function == \"seger\":\n",
    "        UGA = lambda x_A: selection_functions.seger_selection(x_A, dG)\n",
    "        UgA = lambda x_A: selection_functions.seger_selection(x_A, dg)\n",
    "    elif selection_function == \"wright\":\n",
    "        UGA = lambda x_A: selection_functions.wright_selection(x_A, dG)\n",
    "        UgA = lambda x_A: selection_functions.wright_selection(x_A, dg)\n",
    "    else:\n",
    "        valid_funcs = (\"kirkpatrick\", \"seger\", \"wright\")\n",
    "        msg = \"Selection_function must be one of {}, {}, or {}.\".format(*valid_funcs)\n",
    "        raise ValueError(msg)\n",
    "    # define the payoffs\n",
    "    payoff_kernel = np.array([[R, S], [T, P]])\n",
    "\n",
    "    def f(t, y):\n",
    "        W = models.generalized_sexual_selection(y, UGA, UgA, payoff_kernel, M, m, epsilon)\n",
    "        y_dot = models.offspring_genotypes_evolution(W, y)\n",
    "        return y_dot\n",
    "\n",
    "    # common time points at which solution should be evaluated\n",
    "    tmin, tmax = (0, 10000)\n",
    "    ts = np.linspace(tmin, tmax, 100)\n",
    "    \n",
    "    # random initial condition\n",
    "    number_of_genotypes = 4\n",
    "    initial_offspring_share, = prng.dirichlet(np.ones(number_of_genotypes), 1)\n",
    "    y0 = initial_offspring_share\n",
    "    solution1 = integrate.solve_ivp(f, t_span=(tmin, tmax), y0=y0, rtol=1e-9, atol=1e-12, t_eval=ts, vectorized=True)\n",
    "    \n",
    "    # generate new random, nearby initial condition for second simulation\n",
    "    initial_offspring_share, = prng.dirichlet(np.ones(number_of_genotypes), 1)\n",
    "    y0 = (1e-6 * initial_offspring_share) + (1 - (1e-6 * initial_offspring_share)) * y0\n",
    "    solution2 = integrate.solve_ivp(f, t_span=(tmin, tmax), y0=y0, rtol=1e-9, atol=1e-12, t_eval=ts, vectorized=True)\n",
    "    \n",
    "    diff = np.abs(solution1.y - solution2.y)\n",
    "    \n",
    "    return ts, diff\n",
    "    \n",
    "\n",
    "def plot_lyapunov_exponents(selection_function, dG, dg, T, R, P, S, M, m, epsilon):\n",
    "\n",
    "    ts, diff = compute_differences(selection_function, dG, dg, T, R, P, S, M, m, epsilon)\n",
    "    \n",
    "    fig, axes = plt.subplots(1, 4, sharex=True, sharey=True, figsize=(12, 8))\n",
    "    for i in range(4):\n",
    "        axes[i].plot(ts, diff[i])\n",
    "        axes[i].set_yscale(\"log\")\n",
    "        \n",
    "    return (ts, diff)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Interactive simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "e31047e1b54e4185be1698341572c6c6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# sliders used to control the Prisoner's Dilemma Payoffs\n",
    "T_slider = widgets.FloatSlider(value=10, min=0, max=100, step=0.1, description=r\"$T$\")\n",
    "R_slider = widgets.FloatSlider(value=8, min=0, max=100, step=0.1, description=r\"$R$\")\n",
    "P_slider = widgets.FloatSlider(value=6, min=0, max=100, step=0.1, description=r\"$P$\")\n",
    "S_slider = widgets.FloatSlider(value=4, min=0, max=100, step=0.1, description=r\"$S$\")\n",
    "\n",
    "# sliders used to control the metabolic costs\n",
    "M_slider = widgets.FloatSlider(value=1, min=0, max=100, step=0.1, description=r\"$M_G$\")\n",
    "m_slider = widgets.FloatSlider(value=0, min=0, max=100, step=0.1, description=r\"$m_g$\")\n",
    "\n",
    "# slider used to control which selection function is being used\n",
    "U_slider = widgets.Dropdown(options=[\"kirkpatrick\", \"seger\", \"wright\"], index=0, description=r\"$U_{\\gamma(j)A}$\")\n",
    "\n",
    "# slider that controls the parameters of the selection function\n",
    "dG_slider = widgets.FloatSlider(value=1, min=0.0, max=10, step=0.05, description=r\"$d_G$\")\n",
    "dg_slider = widgets.FloatSlider(value=1, min=0.0, max=10, step=0.05, description=r\"$d_g$\")\n",
    "\n",
    "# slider used to control the mutation rate\n",
    "e_slider = widgets.FloatSlider(value=0.0, min=0.0, max=1.0, step=1e-4, description=r\"$\\epsilon$\", readout_format=\".4f\")\n",
    "\n",
    "w = widgets.interactive(plot_lyapunov_exponents,\n",
    "                        selection_function=U_slider, dG=dG_slider, dg=dg_slider, \n",
    "                        T=T_slider, R=R_slider, P=P_slider, S=S_slider,\n",
    "                        M=M_slider, m=m_slider, epsilon=e_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# can get access to the solution!\n",
    "ts, diff = w.result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from scipy import optimize"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "def max_exponent(x, selection_function, epsilon):\n",
    "    dG, dg, T, R, P, S, M, m = x\n",
    "    ts, diff = compute_differences(selection_function, dG, dg, T, R, P, S, M, m, epsilon)\n",
    "    exponent = np.max(np.corrcoef(np.log(diff), ts)[-1,:-1])\n",
    "    return exponent"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/drpugh/anaconda3/lib/python3.6/site-packages/ipykernel_launcher.py:4: RuntimeWarning: divide by zero encountered in log\n",
      "  after removing the cwd from sys.path.\n",
      "/Users/drpugh/anaconda3/lib/python3.6/site-packages/numpy/lib/function_base.py:3087: RuntimeWarning: invalid value encountered in subtract\n",
      "  X -= avg[:, None]\n",
      "/Users/drpugh/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/rk.py:140: RuntimeWarning: invalid value encountered in maximum\n",
      "  scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol\n"
     ]
    },
    {
     "ename": "KeyboardInterrupt",
     "evalue": "",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m                         Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-50-1d2f8807f23b>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0mx0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m8\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     13\u001b[0m result = optimize.minimize(lambda x, selection_function, epsilon: -max_exponent(x, selection_function, epsilon),\n\u001b[0;32m---> 14\u001b[0;31m                            x0, args=(\"kirkpatrick\", 1e-3), constraints=constraints, options={\"disp\": True, \"maxiter\": 10})\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/optimize/_minimize.py\u001b[0m in \u001b[0;36mminimize\u001b[0;34m(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)\u001b[0m\n\u001b[1;32m    493\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'slsqp'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    494\u001b[0m         return _minimize_slsqp(fun, x0, args, jac, bounds,\n\u001b[0;32m--> 495\u001b[0;31m                                constraints, callback=callback, **options)\n\u001b[0m\u001b[1;32m    496\u001b[0m     \u001b[0;32melif\u001b[0m \u001b[0mmeth\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'dogleg'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    497\u001b[0m         return _minimize_dogleg(fun, x0, args, jac, hess,\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/optimize/slsqp.py\u001b[0m in \u001b[0;36m_minimize_slsqp\u001b[0;34m(func, x0, args, jac, bounds, constraints, maxiter, ftol, iprint, disp, eps, callback, **unknown_options)\u001b[0m\n\u001b[1;32m    376\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    377\u001b[0m             \u001b[0;31m# Compute objective function\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 378\u001b[0;31m             \u001b[0mfx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfunc\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    379\u001b[0m             \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    380\u001b[0m                 \u001b[0mfx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfloat\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfx\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/optimize/optimize.py\u001b[0m in \u001b[0;36mfunction_wrapper\u001b[0;34m(*wrapper_args)\u001b[0m\n\u001b[1;32m    290\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mwrapper_args\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    291\u001b[0m         \u001b[0mncalls\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 292\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mfunction\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mwrapper_args\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    293\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    294\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mncalls\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfunction_wrapper\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-50-1d2f8807f23b>\u001b[0m in \u001b[0;36m<lambda>\u001b[0;34m(x, selection_function, epsilon)\u001b[0m\n\u001b[1;32m     11\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     12\u001b[0m \u001b[0mx0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0marray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m10\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m8\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m6\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m4\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m0.5\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 13\u001b[0;31m result = optimize.minimize(lambda x, selection_function, epsilon: -max_exponent(x, selection_function, epsilon),\n\u001b[0m\u001b[1;32m     14\u001b[0m                            x0, args=(\"kirkpatrick\", 1e-3), constraints=constraints, options={\"disp\": True, \"maxiter\": 10})\n",
      "\u001b[0;32m<ipython-input-36-5e1484b7ac69>\u001b[0m in \u001b[0;36mmax_exponent\u001b[0;34m(x, selection_function, epsilon)\u001b[0m\n\u001b[1;32m      1\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mmax_exponent\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mselection_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepsilon\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      2\u001b[0m     \u001b[0mdG\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mR\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mS\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 3\u001b[0;31m     \u001b[0mts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdiff\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcompute_differences\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mselection_function\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdG\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mR\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mP\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mS\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepsilon\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m      4\u001b[0m     \u001b[0mexponent\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmax\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcorrcoef\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mlog\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mdiff\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mts\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      5\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mexponent\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-29-4a75619ac6dd>\u001b[0m in \u001b[0;36mcompute_differences\u001b[0;34m(selection_function, dG, dg, T, R, P, S, M, m, epsilon)\u001b[0m\n\u001b[1;32m     34\u001b[0m     \u001b[0minitial_offspring_share\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mprng\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdirichlet\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mones\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnumber_of_genotypes\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     35\u001b[0m     \u001b[0my0\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0minitial_offspring_share\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 36\u001b[0;31m     \u001b[0msolution1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mintegrate\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve_ivp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_span\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtmin\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtmax\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0my0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrtol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-9\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0matol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m1e-12\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mt_eval\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mts\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mvectorized\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     37\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     38\u001b[0m     \u001b[0;31m# generate new random, nearby initial condition for second simulation\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/ivp.py\u001b[0m in \u001b[0;36msolve_ivp\u001b[0;34m(fun, t_span, y0, method, t_eval, dense_output, events, vectorized, **options)\u001b[0m\n\u001b[1;32m    474\u001b[0m     \u001b[0mstatus\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    475\u001b[0m     \u001b[0;32mwhile\u001b[0m \u001b[0mstatus\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 476\u001b[0;31m         \u001b[0mmessage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstep\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    477\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    478\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0msolver\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mstatus\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;34m'finished'\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36mstep\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    180\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    181\u001b[0m             \u001b[0mt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 182\u001b[0;31m             \u001b[0msuccess\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmessage\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_step_impl\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    183\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    184\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0msuccess\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/rk.py\u001b[0m in \u001b[0;36m_step_impl\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    137\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    138\u001b[0m             y_new, f_new, error = rk_step(self.fun, t, y, self.f, h, self.A,\n\u001b[0;32m--> 139\u001b[0;31m                                           self.B, self.C, self.E, self.K)\n\u001b[0m\u001b[1;32m    140\u001b[0m             \u001b[0mscale\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0matol\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmaximum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mabs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my_new\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mrtol\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m             \u001b[0merror_norm\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnorm\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0merror\u001b[0m \u001b[0;34m/\u001b[0m \u001b[0mscale\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/rk.py\u001b[0m in \u001b[0;36mrk_step\u001b[0;34m(fun, t, y, f, h, A, B, C, E, K)\u001b[0m\n\u001b[1;32m     68\u001b[0m     \u001b[0;32mfor\u001b[0m \u001b[0ms\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mc\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32min\u001b[0m \u001b[0menumerate\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mzip\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mC\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     69\u001b[0m         \u001b[0mdy\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0ms\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 70\u001b[0;31m         \u001b[0mK\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0ms\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mc\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mh\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mdy\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     71\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     72\u001b[0m     \u001b[0my_new\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0my\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mh\u001b[0m \u001b[0;34m*\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdot\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mK\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m-\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mT\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mB\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36mfun\u001b[0;34m(t, y)\u001b[0m\n\u001b[1;32m    137\u001b[0m         \u001b[0;32mdef\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    138\u001b[0m             \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnfev\u001b[0m \u001b[0;34m+=\u001b[0m \u001b[0;36m1\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 139\u001b[0;31m             \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfun_single\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    140\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    141\u001b[0m         \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfun\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mfun\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36mfun_single\u001b[0;34m(t, y)\u001b[0m\n\u001b[1;32m    124\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mvectorized\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    125\u001b[0m             \u001b[0;32mdef\u001b[0m \u001b[0mfun_single\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 126\u001b[0;31m                 \u001b[0;32mreturn\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_fun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mravel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    127\u001b[0m             \u001b[0mfun_vectorized\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_fun\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    128\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/scipy/integrate/_ivp/base.py\u001b[0m in \u001b[0;36mfun_wrapped\u001b[0;34m(t, y)\u001b[0m\n\u001b[1;32m     19\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     20\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mfun_wrapped\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 21\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0masarray\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mfun\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     22\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m     \u001b[0;32mreturn\u001b[0m \u001b[0mfun_wrapped\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my0\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m<ipython-input-29-4a75619ac6dd>\u001b[0m in \u001b[0;36mf\u001b[0;34m(t, y)\u001b[0m\n\u001b[1;32m     22\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     23\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mf\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mt\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 24\u001b[0;31m         \u001b[0mW\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mgeneralized_sexual_selection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0my\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mUGA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mUgA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpayoff_kernel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepsilon\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     25\u001b[0m         \u001b[0my_dot\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mmodels\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0moffspring_genotypes_evolution\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mW\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0my\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     26\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0my_dot\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/models.py\u001b[0m in \u001b[0;36mgeneralized_sexual_selection\u001b[0;34m(x, UGA, UgA, payoff_kernel, M, m, epsilon)\u001b[0m\n\u001b[1;32m     67\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0mgeneralized_sexual_selection\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mUGA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mUgA\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mpayoff_kernel\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mm\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepsilon\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;36m0.0\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     68\u001b[0m     \u001b[0mnumber_of_genotypes\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mx\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mshape\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 69\u001b[0;31m     \u001b[0mx_A\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mx_a\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnp\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m1\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;36m2\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     70\u001b[0m     phenotype_selection_kernel = np.vstack((np.tile(np.array([UGA(x_A), 1 - UGA(x_A)]), (2,2)),\n\u001b[1;32m     71\u001b[0m                                             np.tile(np.array([UgA(x_A), 1 - UgA(x_A)]), (2,2))))\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/numpy/core/fromnumeric.py\u001b[0m in \u001b[0;36msum\u001b[0;34m(a, axis, dtype, out, keepdims)\u001b[0m\n\u001b[1;32m   1832\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0msum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mkwargs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1833\u001b[0m     return _methods._sum(a, axis=axis, dtype=dtype,\n\u001b[0;32m-> 1834\u001b[0;31m                          out=out, **kwargs)\n\u001b[0m\u001b[1;32m   1835\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1836\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/numpy/core/_methods.py\u001b[0m in \u001b[0;36m_sum\u001b[0;34m(a, axis, dtype, out, keepdims)\u001b[0m\n\u001b[1;32m     30\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     31\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 32\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0mumr_sum\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m     33\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     34\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_prod\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0maxis\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdtype\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mout\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mkeepdims\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mFalse\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "\n",
    "constraints = [{\"type\": \"ineq\", \"fun\": lambda x: x[2] - x[3]},\n",
    "               {\"type\": \"ineq\", \"fun\": lambda x: x[3] - x[4]},\n",
    "               {\"type\": \"ineq\", \"fun\": lambda x: x[4] - x[5]},\n",
    "               {\"type\": \"ineq\", \"fun\": lambda x: x[6] - x[7]},\n",
    "               {\"type\": \"ineq\", \"fun\": lambda x: x[5]},\n",
    "               {\"type\": \"ineq\", \"fun\": lambda x: x[7]},\n",
    "               {\"type\": \"ineq\", \"fun\": lambda x: min(0.5 * (x[2] + x[5]), x[4]) - x[6]},\n",
    "               {\"type\": \"ineq\", \"fun\": lambda x: x[0] - x[1]}]\n",
    "               \n",
    "\n",
    "x0 = np.array([2, 1, 10, 8, 6, 4, 1, 0.5])\n",
    "result = optimize.minimize(lambda x, selection_function, epsilon: -max_exponent(x, selection_function, epsilon),\n",
    "                           x0, args=(\"kirkpatrick\", 1e-3), constraints=constraints, options={\"disp\": True, \"maxiter\": 10})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "     fun: -0.6630840114590145\n",
       "     jac: array([   43153.08966418,   496300.09653103,  1330683.69134033,\n",
       "        5510700.25715718,  4114374.28457757,  5349557.40669724,\n",
       "        4907095.94448072,  1580258.93468051])\n",
       " message: 'Positive directional derivative for linesearch'\n",
       "    nfev: 10\n",
       "     nit: 5\n",
       "    njev: 1\n",
       "  status: 8\n",
       " success: False\n",
       "       x: array([  2. ,   1. ,  10. ,   8. ,   6. ,   4. ,   1. ,   0.5])"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}