maximizing-total-offspring.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import ipywidgets as widgets\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from scipy import optimize\n",
    "import sympy as sym\n",
    "\n",
    "import models\n",
    "import symbolics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Maximizing equilibrium total offspring"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2, x3, x4 = sym.symbols('x1, x2, x3, x4', real=True, nonnegative=True)\n",
    "T, R, P, S = sym.symbols('T, R, P, S', real=True, positive=True)\n",
    "M, m = sym.symbols(\"M, m\", real=True, nonnegative=True)\n",
    "epsilon = sym.symbols(\"epsilon\", real=True, nonnegative=True)\n",
    "\n",
    "UGA = symbolics.UGA\n",
    "UgA = symbolics.UgA"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[x1], [x2], [x3], [1 - x1 - x2 - x3]])\n",
    "payoff_kernel = np.array([[R, S], [T, P]])\n",
    "W = models.generalized_sexual_selection(x, UGA, UgA, payoff_kernel, M, m, epsilon)\n",
    "N, = models.total_offspring(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 4.0 \\left(0.5 P - 0.5 m + x_{1} \\left(- 0.5 M + 0.5 m\\right) + x_{1} \\left(- 1.0 P + 0.5 S + 0.5 T\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + x_{1} \\left(1.0 P - 0.5 S - 0.5 T\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + x_{1} \\left(- 0.5 P - 0.5 R + 0.5 S + 0.5 T\\right) \\operatorname{UgA}^{2}{\\left(x_{1} + x_{3} \\right)} + x_{1} \\left(0.5 P + 0.5 R - 0.5 S - 0.5 T\\right) \\operatorname{UGA}^{2}{\\left(x_{1} + x_{3} \\right)} + x_{2} \\left(- 0.5 M + 0.5 m\\right) + x_{2} \\left(- 1.0 P + 0.5 S + 0.5 T\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + x_{2} \\left(1.0 P - 0.5 S - 0.5 T\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + x_{2} \\left(- 0.5 P - 0.5 R + 0.5 S + 0.5 T\\right) \\operatorname{UgA}^{2}{\\left(x_{1} + x_{3} \\right)} + x_{2} \\left(0.5 P + 0.5 R - 0.5 S - 0.5 T\\right) \\operatorname{UGA}^{2}{\\left(x_{1} + x_{3} \\right)} + \\left(- 1.0 P + 0.5 S + 0.5 T\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + \\left(0.5 P + 0.5 R - 0.5 S - 0.5 T\\right) \\operatorname{UgA}^{2}{\\left(x_{1} + x_{3} \\right)}\\right)$"
      ],
      "text/plain": [
       "4.0*(0.5*P - 0.5*m + x1*(-0.5*M + 0.5*m) + x1*(-1.0*P + 0.5*S + 0.5*T)*UGA(x1 + x3) + x1*(1.0*P - 0.5*S - 0.5*T)*UgA(x1 + x3) + x1*(-0.5*P - 0.5*R + 0.5*S + 0.5*T)*UgA(x1 + x3)**2 + x1*(0.5*P + 0.5*R - 0.5*S - 0.5*T)*UGA(x1 + x3)**2 + x2*(-0.5*M + 0.5*m) + x2*(-1.0*P + 0.5*S + 0.5*T)*UGA(x1 + x3) + x2*(1.0*P - 0.5*S - 0.5*T)*UgA(x1 + x3) + x2*(-0.5*P - 0.5*R + 0.5*S + 0.5*T)*UgA(x1 + x3)**2 + x2*(0.5*P + 0.5*R - 0.5*S - 0.5*T)*UGA(x1 + x3)**2 + (-1.0*P + 0.5*S + 0.5*T)*UgA(x1 + x3) + (0.5*P + 0.5*R - 0.5*S - 0.5*T)*UgA(x1 + x3)**2)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(sym.cancel(sym.together(sym.expand(N))), UGA(x1 + x3), UgA(x1 + x3), x1, x2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Total offspring in the two-locus model is a convex combination of total offspring in the one-locus model where all females carry G-allele of the $\\gamma$ gene and the one-locus model where all females carry the g-allele of the $\\gamma$ gene.\n",
    "\n",
    "\\begin{align}\n",
    "N(x_G, x_A; U_G, G_g) =& x_G \\Bigg(2\\bigg(\\big((R + P) - (T + S)\\big)U_G(x_A)^2 + \\big((T + S) - 2P\\big)U_G(x_A) + (P - M)\\bigg)\\Bigg) + (1 - x_G)\\Bigg(2\\bigg(\\big((R + P) - (T + S)\\big)U_g(x_A)^2 + \\big((T + S) - 2P\\big)U_g(x_A) + (P - m)\\bigg)\\Bigg) \\\\\n",
    "=& x_G N(x_A; U_G) + (1 - x_G)N(x_A, U_g)\n",
    "\\end{align}\n",
    "\n",
    "Note that the function for total offspring is linear in $x_G$. This fact implies that equilibrium total offspring will be maximized at either $x_G*=0$ or $x_G^*=1$, depending on parameters.  Thus any stable, fully polymorphic equilibrium will *not* maximize total offspring in equilibrium."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Substitute equilibrium values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "UGA_star, UgA_star, xG_star = sym.symbols(\"UGA_star, UgA_star, xG_star\")\n",
    "equilibrium_total_offspring = N.subs({UGA(x1+x3): UGA_star, UgA(x1+x3): UgA_star, x1: xG_star - x2}).simplify()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 4.0 \\left(0.5 P + UgA_{star}^{2} \\left(0.5 P + 0.5 R - 0.5 S - 0.5 T\\right) + UgA_{star} \\left(- 1.0 P + 0.5 S + 0.5 T\\right)\\right)$"
      ],
      "text/plain": [
       "4.0*(0.5*P + UgA_star**2*(0.5*P + 0.5*R - 0.5*S - 0.5*T) + UgA_star*(-1.0*P + 0.5*S + 0.5*T))"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(sym.cancel(equilibrium_total_offspring.subs({xG_star: 0, m: 0})), UgA_star)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ N\\bigg(\\frac{1}{2\\epsilon + 1}, 1, Ug\\bigg) = \\frac{2(R - M) + 4\\epsilon\\bigg(Ug^2\\big((R + P) - (T + S)\\big) + Ug\\big((T + S) - 2P\\big) + P\\bigg)}{2\\epsilon + 1} $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Make an interactive plot"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "_equilibrium_total_offspring = sym.lambdify((xG_star, UGA_star, UgA_star, T, R, P, S, M, m),\n",
    "                                            equilibrium_total_offspring,\n",
    "                                            modules=\"numpy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_total_offspring(xG_star, T, R, P, S, M):\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(20, 10))\n",
    "    \n",
    "    equilibrium_selection_probs = np.linspace(0, 1, 100).reshape(-1, 1)\n",
    "    UGAs = equilibrium_selection_probs.reshape(-1, 1)\n",
    "    UgAs = equilibrium_selection_probs.reshape(1, -1)\n",
    "    Z = _equilibrium_total_offspring(xG_star, UGAs, UgAs, T, R, P, S, M, 0)\n",
    "    cax = ax.imshow(Z, origin=\"lower\")\n",
    "\n",
    "    contours = ax.contour(Z, colors='w', origin='lower')\n",
    "    ax.clabel(contours, contours.levels, inline=True, fontsize=10)\n",
    "    \n",
    "    ax.set_ylabel(r\"$U_{GA}^*$\", fontsize=20, rotation=\"horizontal\")\n",
    "    ax.set_xlabel(r\"$U_{gA}^*$\", fontsize=20)\n",
    "    ax.set_title(r\"Equilibrium max total offspring for $x_G^*$={} is {}\".format(xG_star, Z.max()), fontsize=25)\n",
    "    ax.grid(False)\n",
    "\n",
    "    # adjust the tick labels\n",
    "    locs, _ = plt.xticks()\n",
    "    plt.xticks(locs[1:], np.linspace(0, 1, locs.size-1))\n",
    "    locs, _ = plt.yticks()\n",
    "    plt.yticks(locs[1:], np.linspace(0, 1, locs.size-1))\n",
    "\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "bf9688247aa446c28b0b0327e7303ad7",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.5, description='$x_G^*$', max=1.0, step=0.01), FloatSlider(value=10.…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "xG_slider = widgets.FloatSlider(value=0.5, min=0, max=1, step=0.01, description=r\"$x_G^*$\")\n",
    "\n",
    "# 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=3, min=0, max=100, step=0.1, description=r\"$R$\")\n",
    "P_slider = widgets.FloatSlider(value=2, min=0, max=100, step=0.1, description=r\"$P$\")\n",
    "S_slider = widgets.FloatSlider(value=1, min=0, max=100, step=0.1, description=r\"$S$\")\n",
    "\n",
    "M_slider = widgets.FloatSlider(value=0, min=0, max=100, step=0.1, description=r\"$M$\")\n",
    "\n",
    "w = widgets.interactive(plot_total_offspring, xG_star=xG_slider, T=T_slider, R=R_slider, P=P_slider, S=S_slider, M=M_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Find the optimal values of $x_G^*, U_{GA}^*, U_{gA}^*$.\n",
    "\n",
    "The number of total offspring can be written as a function of the equilibrium selection probability.\n",
    "\n",
    "$$ N\\big(x_G^*, U_{GA}^*, U_{gA}^*\\big) = 2\\bigg(\\big((R + P) - (T + S)\\big)x_G^*U_{GA}^{*2} + \\big((T + S) - 2P\\big)x_G^*U_{GA}^* + P - \\big((R + P) - (T + S)\\big)x_G^*U_{gA}^{*2} - \\big((T + S) - 2P\\big)x_G^*U_{gA}^* + \\big((R + P) - (T + S)\\big)U_{gA}^{*2} + \\big((T + S) - 2P\\big)U_{gA}^* \\bigg)$$\n",
    "\n",
    "\n",
    "To find the equilibrium selection probability that maximizes the number of total offspring we need to solve the following constrained optimization problem.\n",
    "\n",
    "$$ \\max_{x_G^*, U_{GA}^*, U_{gA}^*}\\ N\\big(x_G^*, U_{GA}^*, U_{gA}^*\\big)  $$\n",
    "\n",
    "subject to the following inequality constraints.\n",
    "\n",
    "\\begin{align}\n",
    "    -x_G^* \\le& 0 \\\\\n",
    "    x_G^* - 1 \\le& 0 \\\\\n",
    "    -U_{GA}^* \\le& 0 \\\\\n",
    "    U_{GA}^* - 1 \\le& 0 \\\\\n",
    "    -U_{gA}^* \\le& 0 \\\\\n",
    "    U_{gA}^* - 1 \\le& 0\n",
    "\\end{align}\n",
    "\n",
    "First-order conditions are as follows.\n",
    "\\begin{align}\n",
    "    2\\bigg(\\big((R + P) - (T + S)\\big)U_{GA}^{*2} + \\big((T + S) - 2P\\big)U_{GA}^* + P - \\big((R + P) - (T + S)\\big)U_{gA}^{*2} - \\big((T + S) - 2P\\big)U_{gA}^*\\bigg) =& -\\mu_{x_G^*, 0} + \\mu_{x_G^*,1} \\\\\n",
    "    2\\bigg(2\\big((R + P) - (T + S)\\big)x_G^*U_{GA}^* + \\big((T + S) - 2P\\big)x_G^*\\bigg) =& -\\mu_{U_{GA}^*, 0} + \\mu_{U_{GA}^*,1} \\\\\n",
    "    2\\bigg(-2\\big((R + P) - (T + S)\\big)x_G^*U_{gA}^* - \\big((T + S) - 2P\\big)x_G^* + 2\\big((R + P) - (T + S)\\big)U_{gA}^* + \\big((T + S) - 2P\\big) \\bigg) =& -\\mu_{U_{gA}^*, 0} + \\mu_{U_{gA}^*,1}\n",
    "\\end{align}\n",
    "\n",
    "Complementary slackness conditions are\n",
    "\n",
    "\\begin{align}\n",
    "    -\\mu_{x_G^*,0}x_G^* =& 0 \\\\\n",
    "    \\mu_{x_G^*,1}\\big(x_G^* - 1\\big) =& 0 \\\\\n",
    "    -\\mu_{U_{GA}^*,0}U_{GA}^* =& 0 \\\\\n",
    "    \\mu_{U_{GA}^*,1}\\big(U_{GA}^* - 1\\big) =& 0\n",
    "    -\\mu_{U_{gA}^*,0}U_{gA}^* =& 0 \\\\\n",
    "    \\mu_{U_{gA}^*,1}\\big(U_{gA}^* - 1\\big) =& 0\n",
    "\\end{align}\n",
    "\n",
    "where $\\mu_0, \\mu_1$ are Lagrange multipliers.\n",
    "\n",
    "### Case 1: interior equilibrium $(0 < x_G^* < 1, 0 < U_{GA}^* < 1, 0 < U_{gA}^* < 1)$\n",
    "\n",
    "In an interior equilibrium, complementary slackness conditions imply that all Lagrange multipliers are zero (i.e., $\\mu_{x_G^*, 0} =\\mu_{x_G^*, 1} =\\mu_{U_{GA}^*, 0} =\\mu_{U_{GA}^*, 0} =\\mu_{U_{gA}^*, 0} =\\mu_{U_{gA}^*, 0}=0$).\n",
    "\n",
    "Our first order conditions reduce to the following.\n",
    "\n",
    "\\begin{align}\n",
    "    2\\bigg(\\big((R + P) - (T + S)\\big)U_{GA}^{*2} + \\big((T + S) - 2P\\big)U_{GA}^* + P - \\big((R + P) - (T + S)\\big)U_{gA}^{*2} - \\big((T + S) - 2P\\big)U_{gA}^*\\bigg) =& 0 \\\\\n",
    "    2\\bigg(2\\big((R + P) - (T + S)\\big)x_G^*U_{GA}^* + \\big((T + S) - 2P\\big)x_G^*\\bigg) =& 0 \\\\\n",
    "    2\\bigg(-2\\big((R + P) - (T + S)\\big)x_G^*U_{gA}^* - \\big((T + S) - 2P\\big)x_G^* + 2\\big((R + P) - (T + S)\\big)U_{gA}^* + \\big((T + S) - 2P\\big) \\bigg) =& 0\n",
    "\\end{align}\n",
    "\n",
    "Rearranging the second first-order condition yields an expression for the optimal value of $U_{GA}^*$.\n",
    "\n",
    "$$ \\bar{U}_{GA}^* = \\frac{1}{2}\\left(\\frac{2P - (T + S)}{(R + P) - (T + S)}\\right) $$\n",
    "\n",
    "Substituting this result into the first first-order condition and rearranging yields an identical exprssion for the optimal value of $U_{gA}^*$.\n",
    "\n",
    "$$ \\bar{U}_{gA}^* = \\frac{1}{2}\\left(\\frac{2P - (T + S)}{(R + P) - (T + S)}\\right) $$\n",
    "\n",
    "Substituting this result into the third first-order condition yields a result which implies that the optimal value for $x_G^*$ is indeterminate (i.e., the objective is flat when holding $U_{GA}^*$ and $U_{gA}^*$ fixed). \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "first_order_conditions = sym.Matrix([equilibrium_total_offspring.diff(xG_star, 1),\n",
    "                                     equilibrium_total_offspring.diff(UGA_star, 1),\n",
    "                                     equilibrium_total_offspring.diff(UgA_star, 1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "optimal_UGA_star, = sym.solve(first_order_conditions[1,0], UGA_star)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "optimal_UgA_star, = sym.solve(first_order_conditions[0,0].subs({UGA_star: optimal_UGA_star}), UgA_star)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# optimal value for xG_star is indeterminate!\n",
    "sym.simplify(first_order_conditions[2,0].subs({UGA_star: optimal_UGA_star, UgA_star: optimal_UgA_star}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "jacobian = first_order_conditions.jacobian([xG_star, UGA_star, UgA_star])\n",
    "simplified_jacobian = sym.simplify(jacobian.subs({UGA_star: optimal_UGA_star, UgA_star: optimal_UgA_star}))\n",
    "e1, e2, e3 = (simplified_jacobian.eigenvals()\n",
    "                                 .keys())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4*P*xG_star + 4*P - 4*R*xG_star + 4*R + 4*S*xG_star - 4*S + 4*T*xG_star - 4*T"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "4*P*xG_star + 4*R*xG_star - 4*S*xG_star - 4*T*xG_star"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Requirement for total offspring to optimal a local maximum at the above values derived above is for the Hessian to be negative semi-definite. This requirement will be satisfied if and only if \n",
    "\n",
    "$$ R + P < T + S. $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Case 2: equilibrium with $\\bar{x}_G^*=1$, $(0 < U_{GA}^* < 1, 0 < U_{gA}^* < 1)$\n",
    "\n",
    "In this equilibrium, complementary slackness conditions imply that all Lagrange multipliers are zero (i.e., $\\mu_{x_G^*, 0} =\\mu_{U_{GA}^*, 0} =\\mu_{U_{GA}^*, 0} =\\mu_{U_{gA}^*, 0} =\\mu_{U_{gA}^*, 0}=0$) except $\\mu_{x_G^*, 1}  > 0$.\n",
    "\n",
    "Our first order conditions reduce to the following.\n",
    "\n",
    "\\begin{align}\n",
    "    2\\bigg(\\big((R + P) - (T + S)\\big)U_{GA}^{*2} + \\big((T + S) - 2P\\big)U_{GA}^* + P - \\big((R + P) - (T + S)\\big)U_{gA}^{*2} - \\big((T + S) - 2P\\big)U_{gA}^*\\bigg) =& \\mu_{x_G^*, 1} \\\\\n",
    "    2\\bigg(2\\big((R + P) - (T + S)\\big)U_{GA}^* + \\big((T + S) - 2P\\big)\\bigg) =& 0 \\\\\n",
    "    2\\bigg(-2\\big((R + P) - (T + S)\\big)U_{gA}^* - \\big((T + S) - 2P\\big) + 2\\big((R + P) - (T + S)\\big)U_{gA}^* + \\big((T + S) - 2P\\big) \\bigg) =& 0\n",
    "\\end{align}\n",
    "\n",
    "Rearranging the second first-order condition yields an expression for the optimal value of $U_{GA}^*$.\n",
    "\n",
    "$$ \\bar{U}_{GA}^* = \\frac{1}{2}\\left(\\frac{2P - (T + S)}{(R + P) - (T + S)}\\right) $$\n",
    "\n",
    "Substituting this optimal value of $U_{GA}^*$ into the first first-order condition and rearranging we find that the inequality will hold so long as \n",
    "\n",
    "$$ \\big((R + P) - (T + S)\\big)\\big(U_{gA}^* - \\bar{U}_{GA}^*\\big)^2 + P > 0 $$\n",
    "\n",
    "which requires $R + P > T + S$. Finally, rearranging the third first-order condition implies that the optimal value for $U_{gA}^*$ is indeterminate: so long as $R + P > T + S$, then $\\bar{x}_G^*=1$ for any value of $U_{gA}^*$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [
    {
     "ename": "SyntaxError",
     "evalue": "unexpected EOF while parsing (<ipython-input-62-fadb5b1c3254>, line 12)",
     "output_type": "error",
     "traceback": [
      "\u001b[0;36m  File \u001b[0;32m\"<ipython-input-62-fadb5b1c3254>\"\u001b[0;36m, line \u001b[0;32m12\u001b[0m\n\u001b[0;31m    else:\u001b[0m\n\u001b[0m         ^\u001b[0m\n\u001b[0;31mSyntaxError\u001b[0m\u001b[0;31m:\u001b[0m unexpected EOF while parsing\n"
     ]
    }
   ],
   "source": [
    "def interior_optimal_UGA(T, R, P, S):\n",
    "    return 0.5 * ((2 * P - (T + S)) / ((R + P) - (T + S)))\n",
    "\n",
    "\n",
    "def interior_optimal_UgA(T, R, P, S):\n",
    "    return interior_optimal_UGA(T, R, P, S)\n",
    "\n",
    "\n",
    "def _mu_xG_1(UGA, UgA, T, R, P, S):\n",
    "    multiplier = 2 * (((R + P) - (T + S)) * UGA**2 + ((T + S) - 2 * P) * UGA + P - \n",
    "                      ((R + P) - (T + S)) * UgA**2 - ((T + S) - 2 * P) * UgA)\n",
    "    return multiplier\n",
    "\n",
    "\n",
    "def max_total_fitness(T, R, P, S):\n",
    "    if _mu_xG_1(UGA, UgA, T, R, P, S) > 0:\n",
    "        pass # max at xG=1\n",
    "    elif _mu_xG_1(UGA, UgA, T, R, P, S) < 0:\n",
    "        pass # max at xG=0\n",
    "    else:\n",
    "        \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "objective = lambda x: -_equilibrium_total_offspring(x[0], x[1], x[2], 25, 3, 2, 1)\n",
    "x0 = 0.5 * np.ones(3)\n",
    "res = optimize.minimize(objective, x0, bounds=[(0,1), (0,1), (0,1)])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "15.523809523809522"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "-res.fun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": 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.10"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}