monomorphic-gamma-g-allele-only-Copy1.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import sympy as sym\n",
    "\n",
    "import models\n",
    "import payoffs\n",
    "import plotting\n",
    "import selection_functions\n",
    "import symbolics"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "plt.style.use(\"ggplot\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# One-locus model\n",
    "\n",
    "Start from full model. Set $x_3=x_4=0, x_1, x_2=1 - x1$. Monomorphic with $\\gamma=G$ only."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1 = sym.symbols('x1', real=True, positive=True)\n",
    "T, R, P, S = sym.symbols('T, R, P, S', real=True, positive=True)\n",
    "M = sym.symbols(\"M\", real=True, positive=True)\n",
    "epsilon = sym.symbols(\"epsilon\", real=True, positive=True)\n",
    "\n",
    "UGA = symbolics.UGA\n",
    "UgA = symbolics.UgA"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Total offspring (fitness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[x1], [1-x1], [0], [0]])\n",
    "payoff_kernel = np.array([[R, S], [T, P]])\n",
    "W = models.generalized_sexual_selection(x, UGA, UgA, payoff_kernel, M=M, m=0, mutation_rate=epsilon)\n",
    "N, = models.total_offspring(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "N = sym.simplify(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 4.0 \\left(- 0.5 M + 0.5 P + \\left(- 1.0 P + 0.5 S + 0.5 T\\right) \\operatorname{UGA}{\\left(x_{1} \\right)} + \\left(0.5 P + 0.5 R - 0.5 S - 0.5 T\\right) \\operatorname{UGA}^{2}{\\left(x_{1} \\right)}\\right)$"
      ],
      "text/plain": [
       "4.0*(-0.5*M + 0.5*P + (-1.0*P + 0.5*S + 0.5*T)*UGA(x1) + (0.5*P + 0.5*R - 0.5*S - 0.5*T)*UGA(x1)**2)"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(N, UGA(x1))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After a bit of algebra we arrive at the following expression for total offspring.\n",
    "\n",
    "\\begin{align}\n",
    "   % N(x_1) =& 2.0*P*UGA(x1)**2 - 4.0*P*UGA(x1) + 2.0*P + 2.0*R*UGA(x1)**2 - 2.0*S*UGA(x1)**2 + 2.0*S*UGA(x1) - 2.0*T*UGA(x1)**2 + 2.0*T*UGA(x1)\\\\\n",
    "   % N(x_1) =& 2\\bigg(PU_{GA}(x_1)^2 - 2PU_{GA}(x_1) + P + RU_{GA}(x_1)^2 - SU_{GA}(x_1)^2 + SU_{GA}(x_1) - TU_{GA}(x_1)^2 + TU_{GA}(x_1)\\bigg) \\\\\n",
    "   % N(x_1) =& 2\\bigg(RU_{GA}(x_1)^2 + PU_{GA}(x_1)^2 - TU_{GA}(x_1)^2 - SU_{GA}(x_1)^2 + SU_{GA}(x_1) + TU_{GA}(x_1) - 2PU_{GA}(x_1) + P\\bigg) \\\\\n",
    "    N(x_1) =& 2\\bigg(\\big((R + P) - (T + S)\\big)U_{GA}(x_1)^2 + \\big((T + S) - 2P\\big)U_{GA}(x_1) + P\\bigg) \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# The locus of potential equilibria (LPE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "(f,), _, _, _ = models.offspring_genotypes_evolution(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solving for the LPE using Python\n",
    "x1_star, = sym.solve(f, x1, implicit=True)\n",
    "x1_star =  sym.cancel(x1_star)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ \\frac{2\\big(U_G(S - M) + (R - S)U_G^2\\big)}{F(U_G)} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{1.0 \\left(4.0 M \\epsilon \\operatorname{UGA}{\\left(x_{1} \\right)} - 4.0 M \\epsilon - 2.0 M \\operatorname{UGA}{\\left(x_{1} \\right)} + 4.0 P \\epsilon \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} - 8.0 P \\epsilon \\operatorname{UGA}{\\left(x_{1} \\right)} + 4.0 P \\epsilon + 2.0 R \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} - 2.0 S \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} + 2.0 S \\operatorname{UGA}{\\left(x_{1} \\right)} - 4.0 T \\epsilon \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} + 4.0 T \\epsilon \\operatorname{UGA}{\\left(x_{1} \\right)}\\right)}{4.0 M \\epsilon + 2.0 M - 4.0 P \\epsilon \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} + 8.0 P \\epsilon \\operatorname{UGA}{\\left(x_{1} \\right)} - 4.0 P \\epsilon - 2.0 P \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} + 4.0 P \\operatorname{UGA}{\\left(x_{1} \\right)} - 2.0 P - 4.0 R \\epsilon \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} - 2.0 R \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} + 4.0 S \\epsilon \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} - 4.0 S \\epsilon \\operatorname{UGA}{\\left(x_{1} \\right)} + 2.0 S \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} - 2.0 S \\operatorname{UGA}{\\left(x_{1} \\right)} + 4.0 T \\epsilon \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} - 4.0 T \\epsilon \\operatorname{UGA}{\\left(x_{1} \\right)} + 2.0 T \\operatorname{UGA}^{2}{\\left(x_{1} \\right)} - 2.0 T \\operatorname{UGA}{\\left(x_{1} \\right)}}$"
      ],
      "text/plain": [
       "-1.0*(4.0*M*epsilon*UGA(x1) - 4.0*M*epsilon - 2.0*M*UGA(x1) + 4.0*P*epsilon*UGA(x1)**2 - 8.0*P*epsilon*UGA(x1) + 4.0*P*epsilon + 2.0*R*UGA(x1)**2 - 2.0*S*UGA(x1)**2 + 2.0*S*UGA(x1) - 4.0*T*epsilon*UGA(x1)**2 + 4.0*T*epsilon*UGA(x1))/(4.0*M*epsilon + 2.0*M - 4.0*P*epsilon*UGA(x1)**2 + 8.0*P*epsilon*UGA(x1) - 4.0*P*epsilon - 2.0*P*UGA(x1)**2 + 4.0*P*UGA(x1) - 2.0*P - 4.0*R*epsilon*UGA(x1)**2 - 2.0*R*UGA(x1)**2 + 4.0*S*epsilon*UGA(x1)**2 - 4.0*S*epsilon*UGA(x1) + 2.0*S*UGA(x1)**2 - 2.0*S*UGA(x1) + 4.0*T*epsilon*UGA(x1)**2 - 4.0*T*epsilon*UGA(x1) + 2.0*T*UGA(x1)**2 - 2.0*T*UGA(x1))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1_star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1_star = sym.factor(x1_star.subs({UGA(x1): 1}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{0.5}{1.0 \\epsilon + 0.5}$"
      ],
      "text/plain": [
       "0.5/(1.0*epsilon + 0.5)"
      ]
     },
     "execution_count": 157,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1_star.subs({UGA(x1): 1, M: 0})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Invadability\n",
    "\n",
    "Restrict the full model to obtain a monomorphic $\\gamma$ equilibrium, add a small amount of individuals carrying the $g$ allele of the $\\gamma$ gene and simulate."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2, x3 = sym.symbols(\"x1, x2, x3\", real=True, positive=True)\n",
    "x4 = 1 - x1 - x2 - x3\n",
    "xA, xG = sym.symbols(\"x_A, x_G\", real=True, positive=True)\n",
    "x = np.array([[x1], [x2], [x3], [x4]])\n",
    "W = models.generalized_sexual_selection(x, UGA, UgA, payoff_kernel, M=M, mutation_rate=epsilon)\n",
    "f1, f2, f3, _ = models.offspring_genotypes_evolution(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "F = sym.Matrix([f1, f2, f3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_jac = F.jacobian([x1, x2, x3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 206,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{4.0 M UG \\epsilon - 2.0 M UG - 4.0 M \\epsilon + 4.0 P UG^{2} \\epsilon - 8.0 P UG \\epsilon + 4.0 P \\epsilon + 2.0 R UG^{2} - 2.0 S UG^{2} + 2.0 S UG - 4.0 T UG^{2} \\epsilon + 4.0 T UG \\epsilon}{- 4.0 M \\epsilon - 2.0 M + 4.0 P UG^{2} \\epsilon + 2.0 P UG^{2} - 8.0 P UG \\epsilon - 4.0 P UG + 4.0 P \\epsilon + 2.0 P + 4.0 R UG^{2} \\epsilon + 2.0 R UG^{2} - 4.0 S UG^{2} \\epsilon - 2.0 S UG^{2} + 4.0 S UG \\epsilon + 2.0 S UG - 4.0 T UG^{2} \\epsilon - 2.0 T UG^{2} + 4.0 T UG \\epsilon + 2.0 T UG}$"
      ],
      "text/plain": [
       "(4.0*M*UG*epsilon - 2.0*M*UG - 4.0*M*epsilon + 4.0*P*UG**2*epsilon - 8.0*P*UG*epsilon + 4.0*P*epsilon + 2.0*R*UG**2 - 2.0*S*UG**2 + 2.0*S*UG - 4.0*T*UG**2*epsilon + 4.0*T*UG*epsilon)/(-4.0*M*epsilon - 2.0*M + 4.0*P*UG**2*epsilon + 2.0*P*UG**2 - 8.0*P*UG*epsilon - 4.0*P*UG + 4.0*P*epsilon + 2.0*P + 4.0*R*UG**2*epsilon + 2.0*R*UG**2 - 4.0*S*UG**2*epsilon - 2.0*S*UG**2 + 4.0*S*UG*epsilon + 2.0*S*UG - 4.0*T*UG**2*epsilon - 2.0*T*UG**2 + 4.0*T*UG*epsilon + 2.0*T*UG)"
      ]
     },
     "execution_count": 206,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(x1_star.subs({UGA(x1): UG}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "UG = sym.symbols(\"UG\", positive=True, real=True)\n",
    "Ug = sym.symbols(\"Ug\", positive=True, real=True)\n",
    "UG_prime = sym.symbols(\"UG_prime\", positive=True, real=True)\n",
    "\n",
    "evaluated_F_jac = sym.simplify(F_jac.subs({x2: 1 - x1, x3:0})\n",
    "                                    .doit()\n",
    "                                    .subs({UGA(x1): UG, UgA(x1): Ug, sym.Derivative(UGA(x1), x1): UG_prime}))\n",
    " "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "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-16-74d037378b80>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0meigenvals\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mevaluated_F_jac\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0meigenvals\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[0m",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/matrices/matrices.py\u001b[0m in \u001b[0;36meigenvals\u001b[0;34m(self, error_when_incomplete, **flags)\u001b[0m\n\u001b[1;32m   1276\u001b[0m                 \u001b[0meigs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mroots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcharpoly\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mDummy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'x'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msimplify\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0msimplify\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1277\u001b[0m             \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1278\u001b[0;31m                 \u001b[0meigs\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mroots\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mmat\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcharpoly\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mDummy\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'x'\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mflags\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   1279\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1280\u001b[0m         \u001b[0;31m# make sure the algebraic multiplicity sums to the\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/matrices/matrices.py\u001b[0m in \u001b[0;36mcharpoly\u001b[0;34m(self, x, simplify)\u001b[0m\n\u001b[1;32m    346\u001b[0m         \u001b[0mberk_vector\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_eval_berkowitz_vector\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    347\u001b[0m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_uniquely_named_symbol\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mberk_vector\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 348\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mPurePoly\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msimplify\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mberk_vector\u001b[0m\u001b[0;34m]\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[0m\n\u001b[0m\u001b[1;32m    349\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    350\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mcofactor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"berkowitz\"\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[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/matrices/matrices.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m    346\u001b[0m         \u001b[0mberk_vector\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_eval_berkowitz_vector\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    347\u001b[0m         \u001b[0mx\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_uniquely_named_symbol\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mberk_vector\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 348\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0mPurePoly\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0msimplify\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mberk_vector\u001b[0m\u001b[0;34m]\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[0m\n\u001b[0m\u001b[1;32m    349\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    350\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0mcofactor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mi\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mj\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mmethod\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;34m\"berkowitz\"\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[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/simplify/simplify.py\u001b[0m in \u001b[0;36msimplify\u001b[0;34m(expr, ratio, measure, rational, inverse, doit, **kwargs)\u001b[0m\n\u001b[1;32m    598\u001b[0m     \u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbottom_up\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mw\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mw\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'normal'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mw\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\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    599\u001b[0m     \u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpowsimp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_content_primitive\u001b[0m\u001b[0;34m(\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[0;32m--> 600\u001b[0;31m     \u001b[0m_e\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcancel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\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    601\u001b[0m     \u001b[0mexpr1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mshorter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_e\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_mexpand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_e\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcancel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# issue 6829\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    602\u001b[0m     \u001b[0mexpr2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mshorter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtogether\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdeep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtogether\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdeep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\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[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/polys/polytools.py\u001b[0m in \u001b[0;36mcancel\u001b[0;34m(f, *gens, **args)\u001b[0m\n\u001b[1;32m   6618\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   6619\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 6620\u001b[0;31m         \u001b[0;34m(\u001b[0m\u001b[0mF\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mG\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mopt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparallel_poly_from_expr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mgens\u001b[0m\u001b[0;34m,\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   6621\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mPolificationFailed\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   6622\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mtuple\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTuple\u001b[0m\u001b[0;34m)\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[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/polys/polytools.py\u001b[0m in \u001b[0;36mparallel_poly_from_expr\u001b[0;34m(exprs, *gens, **args)\u001b[0m\n\u001b[1;32m   4307\u001b[0m     \u001b[0;34m\"\"\"Construct polynomials from expressions. \"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4308\u001b[0m     \u001b[0mopt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moptions\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_options\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgens\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[0;32m-> 4309\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_parallel_poly_from_expr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexprs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mopt\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   4310\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4311\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/polys/polytools.py\u001b[0m in \u001b[0;36m_parallel_poly_from_expr\u001b[0;34m(exprs, opt)\u001b[0m\n\u001b[1;32m   4346\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4347\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mopt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpand\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4348\u001b[0;31m                     \u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpand\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[0m\u001b[1;32m   4349\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4350\u001b[0m             \u001b[0mfailed\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[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36mexpand\u001b[0;34m(self, deep, modulus, power_base, power_exp, mul, log, multinomial, basic, **hints)\u001b[0m\n\u001b[1;32m   3455\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mhints\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'mul'\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3456\u001b[0m                 expr, _ = Expr._expand_hint(\n\u001b[0;32m-> 3457\u001b[0;31m                     expr, '_eval_expand_mul', deep=deep, **hints)\n\u001b[0m\u001b[1;32m   3458\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mhints\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'log'\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3459\u001b[0m                 expr, _ = Expr._expand_hint(\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   3379\u001b[0m             \u001b[0msargs\u001b[0m \u001b[0;34m=\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   3380\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\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[0;32m-> 3381\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\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   3382\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3383\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   3387\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3388\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mhasattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\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[0;32m-> 3389\u001b[0;31m             \u001b[0mnewexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m**\u001b[0m\u001b[0mhints\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   3390\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mnewexpr\u001b[0m \u001b[0;34m!=\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3391\u001b[0m                 \u001b[0;32mreturn\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mnewexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/mul.py\u001b[0m in \u001b[0;36m_eval_expand_mul\u001b[0;34m(self, **hints)\u001b[0m\n\u001b[1;32m    888\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0msums\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    889\u001b[0m                 \u001b[0mdeep\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mhints\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\"deep\"\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--> 890\u001b[0;31m                 \u001b[0mterms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfunc\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expandsums\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msums\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    891\u001b[0m                 \u001b[0margs\u001b[0m \u001b[0;34m=\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    892\u001b[0m                 \u001b[0;32mfor\u001b[0m \u001b[0mterm\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mterms\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/mul.py\u001b[0m in \u001b[0;36m_expandsums\u001b[0;34m(sums)\u001b[0m\n\u001b[1;32m    853\u001b[0m         \u001b[0mright\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMul\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expandsums\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msums\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mL\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[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    854\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 855\u001b[0;31m         \u001b[0mterms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mMul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mleft\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mb\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mright\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    856\u001b[0m         \u001b[0madded\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mAdd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mterms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    857\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mAdd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmake_args\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0madded\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# it may have collapsed down to one term\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/mul.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m    853\u001b[0m         \u001b[0mright\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMul\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expandsums\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msums\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mL\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[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    854\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 855\u001b[0;31m         \u001b[0mterms\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m[\u001b[0m\u001b[0mMul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ma\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mb\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0ma\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mleft\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mb\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mright\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    856\u001b[0m         \u001b[0madded\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mAdd\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mterms\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    857\u001b[0m         \u001b[0;32mreturn\u001b[0m \u001b[0mAdd\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmake_args\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0madded\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# it may have collapsed down to one term\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/operations.py\u001b[0m in \u001b[0;36m__new__\u001b[0;34m(cls, *args, **options)\u001b[0m\n\u001b[1;32m     48\u001b[0m         \u001b[0mis_commutative\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mnc_part\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     49\u001b[0m         \u001b[0mobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_from_args\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc_part\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mnc_part\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mis_commutative\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 50\u001b[0;31m         \u001b[0mobj\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_exec_constructor_postprocessors\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mobj\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     51\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     52\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0morder_symbols\u001b[0m \u001b[0;32mis\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/basic.py\u001b[0m in \u001b[0;36m_exec_constructor_postprocessors\u001b[0;34m(cls, obj)\u001b[0m\n\u001b[1;32m   1832\u001b[0m                 postprocessor_mappings = (\n\u001b[1;32m   1833\u001b[0m                     \u001b[0mBasic\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_constructor_postprocessor_mapping\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mcls\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mitems\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[0;32m-> 1834\u001b[0;31m                     \u001b[0;32mfor\u001b[0m \u001b[0mcls\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mtype\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mi\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mmro\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[0m\u001b[1;32m   1835\u001b[0m                     \u001b[0;32mif\u001b[0m \u001b[0mcls\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mBasic\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_constructor_postprocessor_mapping\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1836\u001b[0m                 )\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "eigenvals = evaluated_F_jac.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "e1, e2, e3 = eigenvals.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 187,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{0.5 \\left(M - R\\right)}{\\left(1.0 M - 1.0 R\\right) \\left(1.0 \\epsilon + 0.5\\right)}$"
      ],
      "text/plain": [
       "0.5*(M - R)/((1.0*M - 1.0*R)*(1.0*epsilon + 0.5))"
      ]
     },
     "execution_count": 187,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1_star.subs({UGA(x1): UG})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [],
   "source": [
    "UG = sym.symbols(\"UG\", positive=True, real=True)\n",
    "Ug = sym.symbols(\"Ug\", positive=True, real=True)\n",
    "UG_prime = sym.symbols(\"UG_prime\", positive=True, real=True)\n",
    "\n",
    "evaluated_F_jac = sym.simplify(F_jac.subs({x2: 1 - x1, x3:0})\n",
    "           .doit()\n",
    "           .subs({UGA(x1): UG, UgA(x1): Ug, sym.Derivative(UGA(x1), x1): UG_prime})\n",
    "           .subs({x1: x1_star.subs({UGA(x1): UG})})\n",
    "           .subs({Ug: x1_star.subs({UGA(x1): UG})})\n",
    "           .subs({UG: 1}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 143,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvals = evaluated_F_jac.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [],
   "source": [
    "e1, e2, e3 = eigenvals.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{2 M UG_{prime} \\epsilon - M UG_{prime} + 2 M \\epsilon + M - 2 R \\epsilon - R - 2 T UG_{prime} \\epsilon + T UG_{prime}}{2 M - 2 R}$"
      ],
      "text/plain": [
       "-(2*M*UG_prime*epsilon - M*UG_prime + 2*M*epsilon + M - 2*R*epsilon - R - 2*T*UG_prime*epsilon + T*UG_prime)/(2*M - 2*R)"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{- 8 P \\epsilon^{3} + 12 P \\epsilon^{2} - 24 R \\epsilon^{2} - 26 R \\epsilon - 3 R - 4 S \\epsilon^{2} + 6 S \\epsilon - 4 T \\epsilon^{2} + 6 T \\epsilon + 2 \\epsilon \\sqrt{16 P^{2} \\epsilon^{4} + 128 P R \\epsilon^{4} + 32 P R \\epsilon^{3} - 8 P R \\epsilon^{2} + 16 P S \\epsilon^{3} + 16 P T \\epsilon^{3} + 16 R^{2} \\epsilon^{2} + 24 R^{2} \\epsilon + 9 R^{2} + 16 R S \\epsilon^{2} + 12 R S \\epsilon + 64 R T \\epsilon^{3} + 16 R T \\epsilon^{2} - 4 R T \\epsilon + 4 S^{2} \\epsilon^{2} + 8 S T \\epsilon^{2} + 4 T^{2} \\epsilon^{2}} + \\sqrt{16 P^{2} \\epsilon^{4} + 128 P R \\epsilon^{4} + 32 P R \\epsilon^{3} - 8 P R \\epsilon^{2} + 16 P S \\epsilon^{3} + 16 P T \\epsilon^{3} + 16 R^{2} \\epsilon^{2} + 24 R^{2} \\epsilon + 9 R^{2} + 16 R S \\epsilon^{2} + 12 R S \\epsilon + 64 R T \\epsilon^{3} + 16 R T \\epsilon^{2} - 4 R T \\epsilon + 4 S^{2} \\epsilon^{2} + 8 S T \\epsilon^{2} + 4 T^{2} \\epsilon^{2}}}{- 32 R \\epsilon^{2} - 32 R \\epsilon - 8 R}$"
      ],
      "text/plain": [
       "-(-8*P*epsilon**3 + 12*P*epsilon**2 - 24*R*epsilon**2 - 26*R*epsilon - 3*R - 4*S*epsilon**2 + 6*S*epsilon - 4*T*epsilon**2 + 6*T*epsilon + 2*epsilon*sqrt(16*P**2*epsilon**4 + 128*P*R*epsilon**4 + 32*P*R*epsilon**3 - 8*P*R*epsilon**2 + 16*P*S*epsilon**3 + 16*P*T*epsilon**3 + 16*R**2*epsilon**2 + 24*R**2*epsilon + 9*R**2 + 16*R*S*epsilon**2 + 12*R*S*epsilon + 64*R*T*epsilon**3 + 16*R*T*epsilon**2 - 4*R*T*epsilon + 4*S**2*epsilon**2 + 8*S*T*epsilon**2 + 4*T**2*epsilon**2) + sqrt(16*P**2*epsilon**4 + 128*P*R*epsilon**4 + 32*P*R*epsilon**3 - 8*P*R*epsilon**2 + 16*P*S*epsilon**3 + 16*P*T*epsilon**3 + 16*R**2*epsilon**2 + 24*R**2*epsilon + 9*R**2 + 16*R*S*epsilon**2 + 12*R*S*epsilon + 64*R*T*epsilon**3 + 16*R*T*epsilon**2 - 4*R*T*epsilon + 4*S**2*epsilon**2 + 8*S*T*epsilon**2 + 4*T**2*epsilon**2))/(-32*R*epsilon**2 - 32*R*epsilon - 8*R)"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e2.subs({M: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 175,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'-(-8*P*epsilon**3 + 12*P*epsilon**2 - 24*R*epsilon**2 - 26*R*epsilon - 3*R - 4*S*epsilon**2 + 6*S*epsilon - 4*T*epsilon**2 + 6*T*epsilon - 2*epsilon*sqrt(16*P**2*epsilon**4 + 128*P*R*epsilon**4 + 32*P*R*epsilon**3 - 8*P*R*epsilon**2 + 16*P*S*epsilon**3 + 16*P*T*epsilon**3 + 16*R**2*epsilon**2 + 24*R**2*epsilon + 9*R**2 + 16*R*S*epsilon**2 + 12*R*S*epsilon + 64*R*T*epsilon**3 + 16*R*T*epsilon**2 - 4*R*T*epsilon + 4*S**2*epsilon**2 + 8*S*T*epsilon**2 + 4*T**2*epsilon**2) - sqrt(16*P**2*epsilon**4 + 128*P*R*epsilon**4 + 32*P*R*epsilon**3 - 8*P*R*epsilon**2 + 16*P*S*epsilon**3 + 16*P*T*epsilon**3 + 16*R**2*epsilon**2 + 24*R**2*epsilon + 9*R**2 + 16*R*S*epsilon**2 + 12*R*S*epsilon + 64*R*T*epsilon**3 + 16*R*T*epsilon**2 - 4*R*T*epsilon + 4*S**2*epsilon**2 + 8*S*T*epsilon**2 + 4*T**2*epsilon**2))/(-32*R*epsilon**2 - 32*R*epsilon - 8*R)'"
      ]
     },
     "execution_count": 175,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "str(e3.subs({M: 0}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 159,
   "metadata": {},
   "outputs": [],
   "source": [
    "delta = (64*P**2*epsilon**4 + \n",
    "         32*P*R*epsilon**2 +\n",
    "         64*P*S*epsilon**3 +\n",
    "         64*P*T*epsilon**3 +\n",
    "         4*R**2 +\n",
    "         16*R*S*epsilon + \n",
    "         16*R*T*epsilon +\n",
    "         16*S**2*epsilon**2 +\n",
    "         32*S*T*epsilon**2 +\n",
    "         16*T**2*epsilon**2 +\n",
    "         \n",
    "         128*P*R*epsilon**4 + \n",
    "         32*P*R*epsilon**3 -  \n",
    "         16*R**2*epsilon**2 + \n",
    "         24*R**2*epsilon + \n",
    "         8*R**2 + \n",
    "         16*R*S*epsilon**2 +  \n",
    "         64*R*T*epsilon**3 + \n",
    "         16*R*T*epsilon**2 +\n",
    "         8*R*S*epsilon -    \n",
    "         48*P**2*epsilon**4 - \n",
    "         24*P*R*epsilon**2 - \n",
    "         48*P*S*epsilon**3 - \n",
    "         48*P*T*epsilon**3 - \n",
    "         3*R**2 -\n",
    "         12*R*S*epsilon - \n",
    "         12*R*T*epsilon - \n",
    "         12*S**2*epsilon**2 - \n",
    "         24*S*T*epsilon**2 - \n",
    "         12*T**2*epsilon**2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 176,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = -8*P*epsilon**3 + 12*P*epsilon**2 - 24*R*epsilon**2 - 26*R*epsilon - 3*R - 4*S*epsilon**2 + 6*S*epsilon - 4*T*epsilon**2 + 6*T*epsilon"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 178,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{- 8 P \\epsilon^{3} + 12 P \\epsilon^{2} - 24 R \\epsilon^{2} - 26 R \\epsilon - 3 R - 4 S \\epsilon^{2} + 6 S \\epsilon - 4 T \\epsilon^{2} + 6 T \\epsilon}{\\left(2 \\epsilon + 1\\right)^{2}}$"
      ],
      "text/plain": [
       "(-8*P*epsilon**3 + 12*P*epsilon**2 - 24*R*epsilon**2 - 26*R*epsilon - 3*R - 4*S*epsilon**2 + 6*S*epsilon - 4*T*epsilon**2 + 6*T*epsilon)/(2*epsilon + 1)**2"
      ]
     },
     "execution_count": 178,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A / (2*epsilon + 1)**2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 174,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{4 R \\left(32 P \\epsilon^{4} + 8 P \\epsilon^{3} - 2 R \\epsilon + R + 4 S \\epsilon^{2} + 2 S \\epsilon + 16 T \\epsilon^{3} + 4 T \\epsilon^{2}\\right)}{\\left(2 \\epsilon + 1\\right)^{4}}$"
      ],
      "text/plain": [
       "4*R*(32*P*epsilon**4 + 8*P*epsilon**3 - 2*R*epsilon + R + 4*S*epsilon**2 + 2*S*epsilon + 16*T*epsilon**3 + 4*T*epsilon**2)/(2*epsilon + 1)**4"
      ]
     },
     "execution_count": 174,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "n = (128*P*R*epsilon**4 + \n",
    "         32*P*R*epsilon**3 -  \n",
    "         #16*R**2*epsilon**2 + \n",
    "         8*R**2*epsilon + \n",
    "         4*R**2 + \n",
    "         16*R*S*epsilon**2 +  \n",
    "         64*R*T*epsilon**3 + \n",
    "         16*R*T*epsilon**2 +\n",
    "         8*R*S*epsilon)\n",
    "d = (2*epsilon + 1)**4\n",
    "sym.simplify(n / d)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 179,
   "metadata": {},
   "outputs": [],
   "source": [
    "K = (128*P*R*epsilon**4 + \n",
    "     32*P*R*epsilon**3 -  \n",
    "     8*R**2*epsilon + \n",
    "     4*R**2 + \n",
    "     16*R*S*epsilon**2 +  \n",
    "     64*R*T*epsilon**3 + \n",
    "     16*R*T*epsilon**2 +\n",
    "     8*R*S*epsilon)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 181,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 128 P R \\epsilon^{4} - 32 P R \\epsilon^{3} + 8 P \\epsilon^{3} - 12 P \\epsilon^{2} + 8 R^{2} \\epsilon - 4 R^{2} - 16 R S \\epsilon^{2} - 8 R S \\epsilon - 64 R T \\epsilon^{3} - 16 R T \\epsilon^{2} + 24 R \\epsilon^{2} + 26 R \\epsilon + 3 R + 4 S \\epsilon^{2} - 6 S \\epsilon + 4 T \\epsilon^{2} - 6 T \\epsilon$"
      ],
      "text/plain": [
       "-128*P*R*epsilon**4 - 32*P*R*epsilon**3 + 8*P*epsilon**3 - 12*P*epsilon**2 + 8*R**2*epsilon - 4*R**2 - 16*R*S*epsilon**2 - 8*R*S*epsilon - 64*R*T*epsilon**3 - 16*R*T*epsilon**2 + 24*R*epsilon**2 + 26*R*epsilon + 3*R + 4*S*epsilon**2 - 6*S*epsilon + 4*T*epsilon**2 - 6*T*epsilon"
      ]
     },
     "execution_count": 181,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(-K - A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\text{True}$"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 140,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "_epsilon = 4e-1\n",
    "e2.subs({T: 4, R: 2, P: 1, S: 0.5, M: 0, epsilon: _epsilon}) < 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\text{True}$"
      ],
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 141,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e3.subs({T: 4, R: 2, P: 1, S: 0.5, M: 0, epsilon: _epsilon}) < 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 162,
   "metadata": {},
   "outputs": [],
   "source": [
    "N, = models.total_offspring(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 167,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{4.0 P^{2} \\epsilon^{4}}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{2.0 P R \\epsilon^{2}}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{4.0 P S \\epsilon^{3}}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{4.0 P T \\epsilon^{3}}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{0.25 R^{2}}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{1.0 R S \\epsilon}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{1.0 R T \\epsilon}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{1.0 S^{2} \\epsilon^{2}}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{2.0 S T \\epsilon^{2}}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625} + \\frac{1.0 T^{2} \\epsilon^{2}}{1.0 \\epsilon^{4} + 2.0 \\epsilon^{3} + 1.5 \\epsilon^{2} + 0.5 \\epsilon + 0.0625}$"
      ],
      "text/plain": [
       "4.0*P**2*epsilon**4/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 2.0*P*R*epsilon**2/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 4.0*P*S*epsilon**3/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 4.0*P*T*epsilon**3/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 0.25*R**2/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 1.0*R*S*epsilon/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 1.0*R*T*epsilon/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 1.0*S**2*epsilon**2/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 2.0*S*T*epsilon**2/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625) + 1.0*T**2*epsilon**2/(1.0*epsilon**4 + 2.0*epsilon**3 + 1.5*epsilon**2 + 0.5*epsilon + 0.0625)"
      ]
     },
     "execution_count": 167,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.expand(sym.simplify(N.subs({x2: 1 - x1, x3: 0}).subs({UGA(x1): x1_star.subs({UGA(x1): 1, M: 0})}).subs({M: 0}))**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 3.85185185185185$"
      ],
      "text/plain": [
       "3.85185185185185"
      ]
     },
     "execution_count": 139,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(N.subs({x2: 1 - x1, x3: 0}).subs({UGA(x1): x1})).subs({x1: x1_star}).subs({M: 0, P: 1, S: 1 / R, T: 2 * R}).subs({R: 2}).subs({epsilon: 4e-1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle x_{1} \\left(\\frac{0.005 x_{3} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.25 \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{1} \\left(\\frac{0.495 x_{3} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.25 \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{1} \\left(\\frac{0.01 x_{1} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.5 x_{2} \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1} + \\frac{0.005 x_{3} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.25 \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{1} \\left(\\frac{0.99 x_{1} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.5 x_{2} \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1} + \\frac{0.495 x_{3} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.25 \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{2} \\left(\\frac{0.25 x_{3} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.005 \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{2} \\left(\\frac{0.25 x_{3} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.495 \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{2} \\left(\\frac{0.5 x_{1} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.01 x_{2} \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1} + \\frac{0.25 x_{3} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.005 \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{2} \\left(\\frac{0.5 x_{1} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.99 x_{2} \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1} + \\frac{0.25 x_{3} \\left(3.0 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.495 \\left(1 - \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UGA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{3} \\left(\\frac{0.005 x_{1} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.25 x_{2} \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{3} \\left(\\frac{0.495 x_{1} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.25 x_{2} \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{3} \\left(\\frac{0.005 x_{1} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.25 x_{2} \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1} + \\frac{0.01 x_{3} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.5 \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + x_{3} \\left(\\frac{0.495 x_{1} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.25 x_{2} \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1} + \\frac{0.99 x_{3} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.5 \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + \\left(\\frac{0.25 x_{1} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.005 x_{2} \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1}\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right) + \\left(\\frac{0.25 x_{1} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.495 x_{2} \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1}\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right) + \\left(- x_{1} - x_{2} - x_{3} + 1\\right) \\left(\\frac{0.25 x_{1} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.005 x_{2} \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1} + \\frac{0.5 x_{3} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.01 \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right) + \\left(- x_{1} - x_{2} - x_{3} + 1\\right) \\left(\\frac{0.25 x_{1} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.495 x_{2} \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right)}{- x_{1} - x_{3} + 1} + \\frac{0.5 x_{3} \\left(3.0 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 1.0\\right) \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}}{x_{1} + x_{3}} + \\frac{0.99 \\left(1 - \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)}\\right) \\left(6 \\operatorname{UgA}{\\left(x_{1} + x_{3} \\right)} + 2\\right) \\left(- x_{1} - x_{2} - x_{3} + 1\\right)}{- x_{1} - x_{3} + 1}\\right)$"
      ],
      "text/plain": [
       "x1*(0.005*x3*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.25*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x1*(0.495*x3*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.25*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x1*(0.01*x1*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.5*x2*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)/(-x1 - x3 + 1) + 0.005*x3*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.25*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x1*(0.99*x1*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.5*x2*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)/(-x1 - x3 + 1) + 0.495*x3*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.25*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x2*(0.25*x3*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.005*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x2*(0.25*x3*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.495*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x2*(0.5*x1*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.01*x2*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)/(-x1 - x3 + 1) + 0.25*x3*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.005*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x2*(0.5*x1*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.99*x2*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)/(-x1 - x3 + 1) + 0.25*x3*(3.0*UGA(x1 + x3) + 1.0)*UGA(x1 + x3)/(x1 + x3) + 0.495*(1 - UGA(x1 + x3))*(6*UGA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x3*(0.005*x1*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.25*x2*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)/(-x1 - x3 + 1)) + x3*(0.495*x1*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.25*x2*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)/(-x1 - x3 + 1)) + x3*(0.005*x1*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.25*x2*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)/(-x1 - x3 + 1) + 0.01*x3*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.5*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + x3*(0.495*x1*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.25*x2*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)/(-x1 - x3 + 1) + 0.99*x3*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.5*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + (0.25*x1*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.005*x2*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)/(-x1 - x3 + 1))*(-x1 - x2 - x3 + 1) + (0.25*x1*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.495*x2*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)/(-x1 - x3 + 1))*(-x1 - x2 - x3 + 1) + (-x1 - x2 - x3 + 1)*(0.25*x1*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.005*x2*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)/(-x1 - x3 + 1) + 0.5*x3*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.01*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1)) + (-x1 - x2 - x3 + 1)*(0.25*x1*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.495*x2*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)/(-x1 - x3 + 1) + 0.5*x3*(3.0*UgA(x1 + x3) + 1.0)*UgA(x1 + x3)/(x1 + x3) + 0.99*(1 - UgA(x1 + x3))*(6*UgA(x1 + x3) + 2)*(-x1 - x2 - x3 + 1)/(-x1 - x3 + 1))"
      ]
     },
     "execution_count": 107,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N.subs({T: 4, R: 2, P: 1, S: 0.5, M: 0, epsilon: 1e-2, UGA(x1): x1_star.subs({UGA(x1): 1})})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 114,
   "metadata": {},
   "outputs": [
    {
     "ename": "NameError",
     "evalue": "name 'eigenvals' is not defined",
     "output_type": "error",
     "traceback": [
      "\u001b[0;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m                                 Traceback (most recent call last)",
      "\u001b[0;32m<ipython-input-114-335ca9b280b0>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[1;32m      2\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m      3\u001b[0m \u001b[0;32mwith\u001b[0m \u001b[0mopen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'non-invadable-eigenvalues.pickle'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'wb'\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mas\u001b[0m \u001b[0mhandle\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m----> 4\u001b[0;31m     \u001b[0mpickle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mdump\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0meigenvals\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhandle\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mprotocol\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mpickle\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mHIGHEST_PROTOCOL\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;31mNameError\u001b[0m: name 'eigenvals' is not defined"
     ]
    }
   ],
   "source": [
    "import pickle\n",
    "\n",
    "with open('non-invadable-eigenvalues.pickle', 'wb') as handle:\n",
    "    pickle.dump(eigenvals, handle, protocol=pickle.HIGHEST_PROTOCOL)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "e1, e2, e3 = eigenvals.keys()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "numerator, denominator = sym.fraction(sym.factor(e1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2*(P*UGA_star**2 - 2*P*UGA_star + P + R*UGA_star**2 - S*UGA_star**2 + S*UGA_star - T*UGA_star**2 + T*UGA_star)**2"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "denominator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "upper_bound, = sym.solve(numerator, UGA_prime_star)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2*epsilon + 1)*(P*UGA_star**2 - 2*P*UGA_star + P + R*UGA_star**2 - S*UGA_star**2 + S*UGA_star - T*UGA_star**2 + T*UGA_star)**2/((2*epsilon - 1)*(2*P*R*UGA_star**2 - 2*P*R*UGA_star - P*S*UGA_star**2 + 2*P*S*UGA_star - P*S - R*T*UGA_star**2))"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(upper_bound)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-(2*R**2*epsilon + R**2)/(2*R**2)"
      ]
     },
     "execution_count": 95,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e1.subs({UGA_star: 1, UGA_prime_star: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 73,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# vectorized numeric repr for the eigenvalue\n",
    "_numeric_e2 = sym.lambdify((UGA_star, UgA_star, T, R, P, S, epsilon),\n",
    "                           e2.subs({x1: x1_star}),\n",
    "                           modules=\"numpy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 75,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def plot_second_eigenvalue(T, R, P, S, epsilon):\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10, 10))\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\"Eigenvalue, $e_2$\", fontsize=25)\n",
    "    ax.grid(\"off\")\n",
    "\n",
    "    equilibrium_selection_probs = np.linspace(0, 1, 100)\n",
    "    UgAs, UGAs = np.meshgrid(equilibrium_selection_probs, equilibrium_selection_probs)\n",
    "    Z = _numeric_e2(UGAs, UgAs, T, R, P, S, epsilon)\n",
    "    cax = ax.imshow(Z, origin=\"lower\")\n",
    "    contours = ax.contour(Z, colors='w', origin='lower')\n",
    "    ax.clabel(contours, contours.levels, inline=True, fontsize=10)\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": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "f12f73414a9c48f985f492a9b5cb74a6",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mpld3.disable_notebook()  # don't need interactive plotting for this!\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=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",
    "# slider used to control the mutation rate\n",
    "e_slider = widgets.FloatSlider(value=0.0, min=0.0, max=1.0, step=1e-3, description=r\"$\\epsilon$\", readout_format=\".3f\")\n",
    "\n",
    "w = widgets.interactive(plot_second_eigenvalue, T=T_slider, R=R_slider, P=P_slider, S=S_slider, epsilon=e_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [],
   "source": [
    "evaluated = (e2.subs({x1: x1_star, UgA_star: x1_star})\n",
    "               .subs({UGA_star: 1, UGA_prime_star: 0}))\n",
    "together = sym.together(evaluated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 107,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "expanded = sym.expand(together)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 108,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4.0*P*epsilon**4/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) + 6.0*P*epsilon**3/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) - 12.0*R*epsilon**3/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) - 13.0*R*epsilon**2/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) - 1.5*R*epsilon/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) - 2.0*S*epsilon**3/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) + 3.0*S*epsilon**2/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) - 2.0*T*epsilon**3/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) + 3.0*T*epsilon**2/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon) + 0.015625*sqrt(65536.0*P**2*epsilon**8 + 65536.0*P**2*epsilon**7 + 16384.0*P**2*epsilon**6 + 524288.0*P*R*epsilon**8 + 655360.0*P*R*epsilon**7 + 229376.0*P*R*epsilon**6 - 8192.0*P*R*epsilon**4 + 65536.0*P*S*epsilon**7 + 65536.0*P*S*epsilon**6 + 16384.0*P*S*epsilon**5 + 65536.0*P*T*epsilon**7 + 65536.0*P*T*epsilon**6 + 16384.0*P*T*epsilon**5 + 65536.0*R**2*epsilon**6 + 163840.0*R**2*epsilon**5 + 151552.0*R**2*epsilon**4 + 61440.0*R**2*epsilon**3 + 9216.0*R**2*epsilon**2 + 65536.0*R*S*epsilon**6 + 114688.0*R*S*epsilon**5 + 65536.0*R*S*epsilon**4 + 12288.0*R*S*epsilon**3 + 262144.0*R*T*epsilon**7 + 327680.0*R*T*epsilon**6 + 114688.0*R*T*epsilon**5 - 4096.0*R*T*epsilon**3 + 16384.0*S**2*epsilon**6 + 16384.0*S**2*epsilon**5 + 4096.0*S**2*epsilon**4 + 32768.0*S*T*epsilon**6 + 32768.0*S*T*epsilon**5 + 8192.0*S*T*epsilon**4 + 16384.0*T**2*epsilon**6 + 16384.0*T**2*epsilon**5 + 4096.0*T**2*epsilon**4)/(16.0*R*epsilon**3 + 16.0*R*epsilon**2 + 4.0*R*epsilon)"
      ]
     },
     "execution_count": 108,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "expanded"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 110,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp_expr = sym.together(expanded)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 111,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-4.0*P*epsilon**3 + 6.0*P*epsilon**2 - 12.0*R*epsilon**2 - 13.0*R*epsilon - 1.5*R - 2.0*S*epsilon**2 + 3.0*S*epsilon - 2.0*T*epsilon**2 + 3.0*T*epsilon + 0.015625*sqrt(65536.0*P**2*epsilon**6 + 65536.0*P**2*epsilon**5 + 16384.0*P**2*epsilon**4 + 524288.0*P*R*epsilon**6 + 655360.0*P*R*epsilon**5 + 229376.0*P*R*epsilon**4 - 8192.0*P*R*epsilon**2 + 65536.0*P*S*epsilon**5 + 65536.0*P*S*epsilon**4 + 16384.0*P*S*epsilon**3 + 65536.0*P*T*epsilon**5 + 65536.0*P*T*epsilon**4 + 16384.0*P*T*epsilon**3 + 65536.0*R**2*epsilon**4 + 163840.0*R**2*epsilon**3 + 151552.0*R**2*epsilon**2 + 61440.0*R**2*epsilon + 9216.0*R**2 + 65536.0*R*S*epsilon**4 + 114688.0*R*S*epsilon**3 + 65536.0*R*S*epsilon**2 + 12288.0*R*S*epsilon + 262144.0*R*T*epsilon**5 + 327680.0*R*T*epsilon**4 + 114688.0*R*T*epsilon**3 - 4096.0*R*T*epsilon + 16384.0*S**2*epsilon**4 + 16384.0*S**2*epsilon**3 + 4096.0*S**2*epsilon**2 + 32768.0*S*T*epsilon**4 + 32768.0*S*T*epsilon**3 + 8192.0*S*T*epsilon**2 + 16384.0*T**2*epsilon**4 + 16384.0*T**2*epsilon**3 + 4096.0*T**2*epsilon**2))/(R*(16.0*epsilon**2 + 16.0*epsilon + 4.0))"
      ]
     },
     "execution_count": 111,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_expr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-13.0*R*epsilon - 1.5*R + 3.0*S*epsilon + 3.0*T*epsilon + 0.015625*sqrt(61440.0*R**2*epsilon + 9216.0*R**2 + 12288.0*R*S*epsilon - 4096.0*R*T*epsilon))/(R*(16.0*epsilon + 4.0))"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_expr.subs({epsilon**2: 0})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "   % (-13.0*R*epsilon - 1.5*R + 3.0*S*epsilon + 3.0*T*epsilon + 0.015625*sqrt(61440.0*R**2*epsilon + 9216.0*R**2 + 12288.0*R*S*epsilon - 4096.0*R*T*epsilon))/(R*(16.0*epsilon + 4.0)) < 0 \\\\\n",
    "    \\frac{(-13R\\epsilon - 1.5R + 3S\\epsilon + 3T\\epsilon + 0.015625\\sqrt{61440R^2\\epsilon + 9216R^2 + 12288RS\\epsilon - 4096RT\\epsilon}}{R(16\\epsilon + 4.0)} < 0 \\\\\n",
    "    \\frac{(-13R\\epsilon - \\frac{3}{2}R + 3S\\epsilon + 3T\\epsilon + \\frac{1}{64}\\sqrt{61440R^2\\epsilon + 9216R^2 + 12288RS\\epsilon - 4096RT\\epsilon}}{R(16\\epsilon + 4.0)} < 0 \\\\\n",
    "    \\frac{(-13R\\epsilon - \\frac{3}{2}R + 3S\\epsilon + 3T\\epsilon + \\frac{1}{2}\\sqrt{60R^2\\epsilon + 9R^2 + 12RS\\epsilon - 4RT\\epsilon}}{R(16\\epsilon + 4)} < 0\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "numerator, denominator = sym.fraction(sym.factor(e2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8*x1*(x1 - 1)*(P*UGA_star**2 - 2*P*UGA_star + P + R*UGA_star**2 - S*UGA_star**2 + S*UGA_star - T*UGA_star**2 + T*UGA_star)"
      ]
     },
     "execution_count": 61,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "denominator"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# vectorized numeric repr for the eigenvalue\n",
    "_numeric_e3 = sym.lambdify((UGA_star, UgA_star, T, R, P, S, epsilon),\n",
    "                           e3.subs({x1: x1_star}),\n",
    "                           modules=\"numpy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def plot_third_eigenvalue(T, R, P, S, epsilon):\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10, 10))\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\"Eigenvalue, $e_3$\", fontsize=25)\n",
    "    ax.grid(\"off\")\n",
    "\n",
    "    equilibrium_selection_probs = np.linspace(0, 1, 100)\n",
    "    UgAs, UGAs = np.meshgrid(equilibrium_selection_probs, equilibrium_selection_probs)\n",
    "    Z = _numeric_e3(UGAs, UgAs, T, R, P, S, epsilon)\n",
    "    cax = ax.imshow(Z, origin=\"lower\")\n",
    "    contours = ax.contour(Z, colors='w', origin='lower')\n",
    "    ax.clabel(contours, contours.levels, inline=True, fontsize=10)\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": 79,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "32df2ee80b4f4d9d8b03c96b215f1e76",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mpld3.disable_notebook()  # don't need interactive plotting for this!\n",
    "\n",
    "# sliders used to control the Prisoner's Dilemma Payoffs\n",
    "T_slider = widgets.FloatSlider(value=25, 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",
    "# slider used to control the mutation rate\n",
    "e_slider = widgets.FloatSlider(value=0.0, min=0.0, max=1.0, step=1e-3, description=r\"$\\epsilon$\", readout_format=\".3f\")\n",
    "\n",
    "w = widgets.interactive(plot_third_eigenvalue, T=T_slider, R=R_slider, P=P_slider, S=S_slider, epsilon=e_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "False"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e2 == e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-13.0*R*epsilon - 1.5*R + 3.0*S*epsilon + 3.0*T*epsilon + 0.015625*sqrt(61440.0*R**2*epsilon + 9216.0*R**2 + 12288.0*R*S*epsilon - 4096.0*R*T*epsilon))/(R*(16.0*epsilon + 4.0))"
      ]
     },
     "execution_count": 164,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(e2.subs({x1: x1_star, UgA_star: x1_star}).subs({UGA_star: 1, UGA_prime_star: 0})).subs({epsilon**2: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 165,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-13.0*R*epsilon - 1.5*R + 3.0*S*epsilon + 3.0*T*epsilon - 0.015625*sqrt(61440.0*R**2*epsilon + 9216.0*R**2 + 12288.0*R*S*epsilon - 4096.0*R*T*epsilon))/(R*(16.0*epsilon + 4.0))"
      ]
     },
     "execution_count": 165,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(e3.subs({x1: x1_star, UgA_star: x1_star}).subs({UGA_star: 1, UGA_prime_star: 0})).subs({epsilon**2: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 153,
   "metadata": {},
   "outputs": [],
   "source": [
    "evaluated = (e3.subs({x1: x1_star, UgA_star: x1_star})\n",
    "               .subs({UGA_star: 1, UGA_prime_star: 0}))\n",
    "together = sym.together(evaluated)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 154,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "expanded = sym.expand(together)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "tmp_expr = sym.together(expanded)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 155,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-13.0*R*epsilon - 1.5*R + 3.0*S*epsilon + 3.0*T*epsilon + 0.015625*sqrt(61440.0*R**2*epsilon + 9216.0*R**2 + 12288.0*R*S*epsilon - 4096.0*R*T*epsilon))/(R*(16.0*epsilon + 4.0))"
      ]
     },
     "execution_count": 155,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tmp_expr.subs({epsilon**2: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "(-13.0*R*epsilon - 1.5*R + 3.0*S*epsilon + 3.0*T*epsilon + 0.015625*sqrt(61440.0*R**2*epsilon + 9216.0*R**2 + 12288.0*R*S*epsilon - 4096.0*R*T*epsilon))/(R*(16.0*epsilon + 4.0))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "First eigenvalue can be written as a hyperbola in $U_{gA}^*,U_{GA}^*$ plane.\n",
    "\n",
    "$$ e_1 = \\left(\\frac{\\big((R + P) - (T + S)\\big)U_{gA}^{*2} + \\big((T + S) - 2P\\big)U_{gA}^*}{N^*}\\right) - \\left(\\frac{\\big((R + P) - (T + S)\\big)U_{GA}^{*2} + \\big((T + S) - 2P\\big)U_{GA}^*}{N^*}\\right) $$\n",
    "\n",
    "### Completing the square\n",
    "\n",
    "To write this hyperbola in standard form we need to [complete the square](https://en.wikipedia.org/wiki/Completing_the_square). Completing the square for the quadratic polynomial in $U_{GA}(x^*)$ yields the following.\n",
    "\n",
    "\\begin{align}\n",
    "    \\big((R + P) - (T + S)\\big)\\left(U_{GA}^* - \\bar{U}_{GA}^*\\right)^2 - \\frac{1}{2}\\big(\\bar{N}^* - 2P\\big) \\\\\n",
    "\\end{align}\n",
    "\n",
    "where $\\bar{U}_{GA}^*$ is the value of $U_{GA}^*$ that maximizes total offspring $N$ derived above. Completing the square for the quadratic polynomial in $U_{gA}(x^*)$ yields a similar expression.\n",
    "\n",
    "\\begin{align}\n",
    "    \\big((R + P) - (T + S)\\big)\\left(U_{gA}^* - \\bar{U}_{GA}^*\\right)^2 - \\frac{1}{2}\\big(\\bar{N}^* - 2P\\big) \\\\\n",
    "\\end{align}\n",
    "\n",
    "Substituting these results into the expression for the eigenvalue yields the following.\n",
    "\n",
    "$$ e_1 = \\left(\\frac{\\big(U_{gA}^* - \\bar{U}_{GA}^*\\big)^{2}}{\\frac{N^*}{\\big((R + P) - (T + S)\\big)}}\\right) - \\left(\\frac{\\big(U_{GA}^* - \\bar{U}_{GA}^*\\big)^{2}}{\\frac{N^*}{\\big((R + P) - (T + S)\\big)}}\\right) $$\n",
    "\n",
    "\n",
    "### Conditions for negative eigenvalue\n",
    "\n",
    "Non-invadability requires that this eigenvalue is strictly negative.\n",
    "\n",
    "$$ \\left(\\frac{\\big(U_{gA}^* - \\bar{U}_{GA}^*\\big)^{2}}{\\frac{N^*}{\\big((R + P) - (T + S)\\big)}}\\right) - \\left(\\frac{\\big(U_{GA}^* - \\bar{U}_{GA}^*\\big)^{2}}{\\frac{N^*}{\\big((R + P) - (T + S)\\big)}}\\right) < 0$$\n",
    "\n",
    "#### Case: (R + P) < (T + S)\n",
    "In this case the inequality simplifies to the following.\n",
    "\n",
    "$$ \\big(U_{gA}^* - \\bar{U}_{GA}^*\\big)^{2} - \\big(U_{GA}^* - \\bar{U}_{GA}^*\\big)^{2} > 0 $$\n",
    "\n",
    "This implies that we have four sub-cases to consider.\n",
    "\n",
    "\\begin{align}\n",
    "    \\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) >& \\big(U_{GA}^* - \\bar{U}_{GA}^*\\big) \\\\\n",
    "    \\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) >& -\\big(U_{GA}^* - \\bar{U}_{GA}^*\\big) \\\\\n",
    "    -\\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) >& \\big(U_{GA}^* - \\bar{U}_{GA}^*\\big) \\\\\n",
    "    -\\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) >& -\\big(U_{GA}^* - \\bar{U}_{GA}^*\\big)\n",
    "\\end{align}\n",
    "\n",
    "These four sub-cases define the following two regions of the $U_{gA}^*,U_{GA}^*$ plane where the eigenvalue is negative.\n",
    "\n",
    "\\begin{align}\n",
    "   & 2\\big(\\bar{U}_{GA}^* - U_{gA}^*\\big) < U_{GA}^* - U_{gA}^* < 0\\\\\n",
    "   & 0 < U_{GA}^* - U_{gA}^* < 2\\big(\\bar{U}_{GA}^* - U_{gA}\\big)\n",
    "\\end{align}\n",
    "\n",
    "#### Case: (R + P) > (T + S)\n",
    "In this case the inequality simplifies to the following.\n",
    "\n",
    "$$ \\big(U_{gA}^* - \\bar{U}_{GA}^*\\big)^{2} - \\big(U_{GA}^* - \\bar{U}_{GA}^*\\big)^{2} < 0 $$\n",
    "\n",
    "This implies that we again have four sub-cases to consider.\n",
    "\n",
    "\\begin{align}\n",
    "    \\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) <& \\big(U_{GA}^* - \\bar{U}_{GA}^*\\big) \\\\\n",
    "    \\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) <& -\\big(U_{GA}^* - \\bar{U}_{GA}^*\\big) \\\\\n",
    "    -\\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) <& \\big(U_{GA}^* - \\bar{U}_{GA}^*\\big) \\\\\n",
    "    -\\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) <& -\\big(U_{GA}^* - \\bar{U}_{GA}^*\\big)\n",
    "\\end{align}\n",
    "\n",
    "These four sub-cases define the following two regions of the $U_{gA}^*,U_{GA}^*$ plane where the eigenvalue is negative.\n",
    "\n",
    "\\begin{align}\n",
    "   & 2\\big(\\bar{U}_{GA}^* - U_{gA}^*\\big) < U_{GA}^* - U_{gA}^* < 0\\\\\n",
    "   & 0 < U_{GA}^* - U_{gA}^* < 2\\big(\\bar{U}_{GA}^* - U_{gA}\\big)\n",
    "\\end{align}\n",
    "\n",
    "\n",
    "### Asymptotes \n",
    "\n",
    "Asymptotes of this hyperbola can be defined as follows.\n",
    "\n",
    "$$ U_{GA}^* - \\bar{U}_{GA}^* = \\pm\\big(U_{gA}^* - \\bar{U}_{GA}^*\\big) $$\n",
    "\n",
    "The above can be written as linear equations in $U_{gA}^*,U_{GA}^*$ plane.\n",
    "\n",
    "\\begin{align}\n",
    "    U_{GA}^* =& U_{gA}^* \\\\\n",
    "    U_{GA}^* =& 2\\bar{U}_{GA}^* - U_{gA}^*\n",
    "\\end{align}\n",
    "\n",
    "Note that the asymptotes of this hyperbola correspond to the locus of points for which the value of the first eigenvalue is exactly zero.\n",
    "\n",
    "#### Discussion\n",
    "\n",
    "Suppose $U_{GA}^* = \\bar{U}_{GA}^*$, then both of the sets of inequalities will hold for any value of $U_{gA}^*$. This implies that the eigenvalue will be negative for any value of $U_{gA}^*$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# vectorized numeric repr for the eigenvalue\n",
    "_numeric_e1 = sym.lambdify((UGA_star, UGA_prime_star, T, R, P, S, M), e1, modules=\"numpy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def plot_first_eigenvalue(T, R, P, S, M):\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10, 10))\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\"Eigenvalue, $e_1$\", fontsize=25)\n",
    "    ax.grid(\"off\")\n",
    "\n",
    "    equilibrium_selection_probs = np.linspace(0, 1, 100).reshape(-1, 1)\n",
    "    equilibrium_selection_derivs = np.linspace(0, 10, 1000).reshape(1, -1) # sensible UGA'(0) >= 1\n",
    "    Z = _numeric_e1(equilibrium_selection_probs, equilibrium_selection_derivs, T, R, P, S, M)\n",
    "    cax = ax.imshow(Z, origin=\"lower\", aspect=\"auto\", vmin=-0.5, vmax=1.0)\n",
    "\n",
    "    levels = np.arange(-0.25, 1.25, 0.25)\n",
    "    contours = ax.contour(Z, levels=levels, colors='w', origin='lower')\n",
    "    ax.clabel(contours, contours.levels, inline=True, fontsize=10)\n",
    "    \n",
    "    # adjust the tick labels\n",
    "    locs, _ = plt.xticks()\n",
    "    plt.xticks(locs[1:], np.linspace(0, 10, 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": 85,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "31c7e73c14f147b5b4a90a12dbfa23cf",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mpld3.disable_notebook()  # don't need interactive plotting for this!\n",
    "\n",
    "# sliders used to control the Prisoner's Dilemma Payoffs\n",
    "T_slider = widgets.FloatSlider(value=25, 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=1, min=0, max=100, step=0.1, description=r\"$M$\")\n",
    "\n",
    "w = widgets.interactive(plot_first_eigenvalue, T=T_slider, R=R_slider, P=P_slider, S=S_slider, M=M_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "UGA_prime_star, = sym.solve(numerator, UGA_prime_star)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-(P*UGA_star**2 - 2*P*UGA_star + P + R*UGA_star**2 - S*UGA_star**2 + S*UGA_star - T*UGA_star**2 + T*UGA_star)**2/(2*P*R*UGA_star**2 - 2*P*R*UGA_star - P*S*UGA_star**2 + 2*P*S*UGA_star - P*S - R*T*UGA_star**2)"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(UGA_prime_star)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After a bit of fiddling, we arrive at the same stability condition derived above."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "   % U'_{GA}(x^*) < -\\frac{(PU_{GA}(x^*)^2 - 2PU_{GA}(x^*) + P + RU_{GA}(x^*)^2 - SU_{GA}(x^*)^2 + SU_{GA}(x^*) - TU_{GA}(x^*)^2 + TU_{GA}(x^*))^2}{(2PRU_{GA}(x^*)^2 - 2PRU_{GA}(x^*) - PSU_{GA}(x^*)^2 + 2PSU_{GA}(x^*) - PS - RTU_{GA}(x^*)^2)} \\\\\n",
    "    U'_{GA}(x^*) < \\frac{\\bigg(\\big((R + P) - (T + S)\\big)U_{GA}(x^*)^2 + \\big((T + S) - 2P\\big)U_{GA}(x^*) + P\\bigg)^2}{\\big(R(T-2P) + PS\\big)U_{GA}(x^*)^2 + 2P(R - S)U_{GA}(x^*) + PS}\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-3*M + 3*P*UGA_star**2 - 6*P*UGA_star - P*UgA_star**2 + 2*P*UgA_star + 2*P + 3*R*UGA_star**2 - R*UgA_star**2 - 3*S*UGA_star**2 + 3*S*UGA_star + S*UgA_star**2 - S*UgA_star - 3*T*UGA_star**2 + 3*T*UGA_star + T*UgA_star**2 - T*UgA_star + m)/(4*M - 4*P*UGA_star**2 + 8*P*UGA_star - 4*P - 4*R*UGA_star**2 + 4*S*UGA_star**2 - 4*S*UGA_star + 4*T*UGA_star**2 - 4*T*UGA_star)"
      ]
     },
     "execution_count": 72,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$$ e_2 = \\left(\\frac{\\big((R + P) - (T + S)\\big)U_{gA}^{*2} + \\big((T + S) - 2P\\big)U_{gA}^* + P}{4N^*}\\right) - \\frac{3}{4} $$\n",
    "\n",
    "### Completing the square\n",
    "\n",
    "Completing the square for the quadratic polynomial in $U_{gA}(x^*)$ yields the following expression.\n",
    "\n",
    "\\begin{align}\n",
    "    \\big((R + P) - (T + S)\\big)\\left(U_{gA}^* - \\bar{U}_{GA}^*\\right)^2 + \\frac{1}{2}\\bar{N}^* \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we need to find conditions under which the third eigenvalue is strictly negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 78,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# vectorized numeric repr for the eigenvalue\n",
    "_numeric_e2 = sym.lambdify((UGA_star, UgA_star, T, R, P, S, M, m), e2, modules=\"numpy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 79,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def plot_second_eigenvalue(T, R, P, S, M, m):\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10, 10))\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\"Eigenvalue, $e_2$\", fontsize=25)\n",
    "    ax.grid(\"off\")\n",
    "\n",
    "    equilibrium_selection_probs = np.linspace(0, 1, 100)\n",
    "    UgAs, UGAs = np.meshgrid(equilibrium_selection_probs, equilibrium_selection_probs)\n",
    "    Z = _numeric_e2(UGAs, UgAs, T, R, P, S, M, m)\n",
    "    cax = ax.imshow(Z, origin=\"lower\")\n",
    "    contours = ax.contour(Z, colors='w', origin='lower')\n",
    "    ax.clabel(contours, contours.levels, inline=True, fontsize=10)\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": 80,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "6b8b65bf9d8e48fa90b4e6e3adc2c68f",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mpld3.disable_notebook()  # don't need interactive plotting for this!\n",
    "\n",
    "# sliders used to control the Prisoner's Dilemma Payoffs\n",
    "T_slider = widgets.FloatSlider(value=25, 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=1, min=0, max=100, step=0.1, description=r\"$M$\")\n",
    "m_slider = widgets.FloatSlider(value=0, min=0, max=100, step=0.1, description=r\"$m$\")\n",
    "\n",
    "w = widgets.interactive(plot_second_eigenvalue, T=T_slider, R=R_slider, P=P_slider, S=S_slider, M=M_slider, m=m_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(-M + P*UGA_star**2 - 2*P*UGA_star - P*UgA_star**2 + 2*P*UgA_star + R*UGA_star**2 - R*UgA_star**2 - S*UGA_star**2 + S*UGA_star + S*UgA_star**2 - S*UgA_star - T*UGA_star**2 + T*UGA_star + T*UgA_star**2 - T*UgA_star + m)/(2*M - 2*P*UGA_star**2 + 4*P*UGA_star - 2*P - 2*R*UGA_star**2 + 2*S*UGA_star**2 - 2*S*UGA_star + 2*T*UGA_star**2 - 2*T*UGA_star)"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "# vectorized numeric repr for the eigenvalue\n",
    "_numeric_e3 = sym.lambdify((UGA_star, UgA_star, T, R, P, S, M, m), e3, modules=\"numpy\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "def plot_third_eigenvalue(T, R, P, S, M, m):\n",
    "\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10, 10))\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\"Eigenvalue, $e_3$\", fontsize=25)\n",
    "    ax.grid(\"off\")\n",
    "\n",
    "    equilibrium_selection_probs = np.linspace(0, 1, 100)\n",
    "    UgAs, UGAs = np.meshgrid(equilibrium_selection_probs, equilibrium_selection_probs)\n",
    "    Z = _numeric_e3(UGAs, UgAs, T, R, P, S, M, m)\n",
    "    cax = ax.imshow(Z, origin=\"lower\")\n",
    "    contours = ax.contour(Z, colors='w', origin='lower')\n",
    "    ax.clabel(contours, contours.levels, inline=True, fontsize=10)\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": 90,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "40a14cf559fa4ca89dec92de7b9202f9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "A Jupyter Widget"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "mpld3.disable_notebook()  # don't need interactive plotting for this!\n",
    "\n",
    "# sliders used to control the Prisoner's Dilemma Payoffs\n",
    "T_slider = widgets.FloatSlider(value=25, 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=1, min=0, max=100, step=0.1, description=r\"$M$\")\n",
    "m_slider = widgets.FloatSlider(value=0, min=0, max=100, step=0.1, description=r\"$m$\")\n",
    "\n",
    "w = widgets.interactive(plot_third_eigenvalue, T=T_slider, R=R_slider, P=P_slider, S=S_slider, M=M_slider, m=m_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## When is the fitness-maximizing interior equilibrium un-invadable by randomista?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 147,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "evaluated_e3 = e3.subs({UGA_star: optimal_UGA_star, UgA_star: optimal_x1_star, m: 0})\n",
    "simplified_e3 = sym.factor(sym.cancel(sym.together(evaluated_e3)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 148,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "numerator, denominator = sym.fraction(simplified_e3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 151,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-4.0*M**3*P**3 - 12.0*M**3*P**2*R + 12.0*M**3*P**2*S + 12.0*M**3*P**2*T - 12.0*M**3*P*R**2 + 24.0*M**3*P*R*S + 24.0*M**3*P*R*T - 12.0*M**3*P*S**2 - 24.0*M**3*P*S*T - 12.0*M**3*P*T**2 - 4.0*M**3*R**3 + 12.0*M**3*R**2*S + 12.0*M**3*R**2*T - 12.0*M**3*R*S**2 - 24.0*M**3*R*S*T - 12.0*M**3*R*T**2 + 4.0*M**3*S**3 + 12.0*M**3*S**2*T + 12.0*M**3*S*T**2 + 4.0*M**3*T**3 + 8.0*M**2*P**3*R + 16.0*M**2*P**2*R**2 - 16.0*M**2*P**2*R*S - 16.0*M**2*P**2*R*T - 2.0*M**2*P**2*S**2 - 4.0*M**2*P**2*S*T - 2.0*M**2*P**2*T**2 + 8.0*M**2*P*R**3 - 16.0*M**2*P*R**2*S - 16.0*M**2*P*R**2*T + 4.0*M**2*P*R*S**2 + 8.0*M**2*P*R*S*T + 4.0*M**2*P*R*T**2 + 4.0*M**2*P*S**3 + 12.0*M**2*P*S**2*T + 12.0*M**2*P*S*T**2 + 4.0*M**2*P*T**3 - 2.0*M**2*R**2*S**2 - 4.0*M**2*R**2*S*T - 2.0*M**2*R**2*T**2 + 4.0*M**2*R*S**3 + 12.0*M**2*R*S**2*T + 12.0*M**2*R*S*T**2 + 4.0*M**2*R*T**3 - 2.0*M**2*S**4 - 8.0*M**2*S**3*T - 12.0*M**2*S**2*T**2 - 8.0*M**2*S*T**3 - 2.0*M**2*T**4 - 4.0*M*P**3*R**2 - 4.0*M*P**2*R**3 + 4.0*M*P**2*R**2*S + 4.0*M*P**2*R**2*T + 2.0*M*P**2*R*S**2 + 4.0*M*P**2*R*S*T + 2.0*M*P**2*R*T**2 + 2.0*M*P*R**2*S**2 + 4.0*M*P*R**2*S*T + 2.0*M*P*R**2*T**2 - 2.0*M*P*R*S**3 - 6.0*M*P*R*S**2*T - 6.0*M*P*R*S*T**2 - 2.0*M*P*R*T**3 - 0.25*M*P*S**4 - 1.0*M*P*S**3*T - 1.5*M*P*S**2*T**2 - 1.0*M*P*S*T**3 - 0.25*M*P*T**4 - 0.25*M*R*S**4 - 1.0*M*R*S**3*T - 1.5*M*R*S**2*T**2 - 1.0*M*R*S*T**3 - 0.25*M*R*T**4 + 0.25*M*S**5 + 1.25*M*S**4*T + 2.5*M*S**3*T**2 + 2.5*M*S**2*T**3 + 1.25*M*S*T**4 + 0.25*M*T**5 - 1.0*P**2*R**2*S**2 + 2.0*P**2*R**2*S*T - 1.0*P**2*R**2*T**2 + 1.0*P**2*R*S**3 - 1.0*P**2*R*S**2*T - 1.0*P**2*R*S*T**2 + 1.0*P**2*R*T**3 - 0.25*P**2*S**4 + 0.5*P**2*S**2*T**2 - 0.25*P**2*T**4 + 1.0*P*R**2*S**3 - 1.0*P*R**2*S**2*T - 1.0*P*R**2*S*T**2 + 1.0*P*R**2*T**3 - 1.0*P*R*S**4 + 2.0*P*R*S**2*T**2 - 1.0*P*R*T**4 + 0.25*P*S**5 + 0.25*P*S**4*T - 0.5*P*S**3*T**2 - 0.5*P*S**2*T**3 + 0.25*P*S*T**4 + 0.25*P*T**5 - 0.25*R**2*S**4 + 0.5*R**2*S**2*T**2 - 0.25*R**2*T**4 + 0.25*R*S**5 + 0.25*R*S**4*T - 0.5*R*S**3*T**2 - 0.5*R*S**2*T**3 + 0.25*R*S*T**4 + 0.25*R*T**5 - 0.0625*S**6 - 0.125*S**5*T + 0.0625*S**4*T**2 + 0.25*S**3*T**3 + 0.0625*S**2*T**4 - 0.125*S*T**5 - 0.0625*T**6"
      ]
     },
     "execution_count": 151,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "numerator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{align}\n",
    "    4\\bigg(-\\big(P^3 + 3P^2R + 3PR^2 + R^3\\big) - \\big(P^3 - 3P^2S + 3PS^2 + S^3\\big) + \\big(-P^3 + 3P^2T - 3PT^2 + T^3\\big) - \\big(R^3 - 3R^2S + 3RS^2 - S^3\\big) + \\big(T^3 - 3RT^2 + 3R^2T - R^3\\big) + 2P^3 + 2R^3 + 8PRS + 8PRT - 8PST - 8RST - 2T^3 + \\big(S^3 + 3S^2T + 3ST^2 + T^3\\big)\\bigg)M^3 + 2\\big(4P^3R + 8P^2R^2 - 8P^2RS - 8P^2RT - P^2S^2 - 2P^2ST - P^2T^2 + 4PR^3 - 8PR^2S - 8PR^2T + 2PRS^2 + 4PRST + 2PRT^2 + 2PS^3 + 6PS^2T + 6PST^2 + 2PT^3 - R^2S^2 - 2R^2ST - R^2T^2 + 2RS^3 + 6RS^2T + 6RST^2 + 2RT^3 - S^4 - 4S^3T - 6S^2T^2 - 4ST^3 - T^4\\big)M^2 + \\big(- 4P^3R^2 - 4P^2R^3 + 4P^2R^2S + 4P^2R^2T + 2P^2RS^2 + 4P^2RST + 2P^2RT^2 + 2PR^2S^2 + 4PR^2ST + 2PR^2T^2 - 2PRS^3 - 6PRS^2T - 6PRST^2 - 2PRT^3 - 0.25PS^4 -PS^3T - 1.5PS^2T^2 - PST^3 - 0.25PT^4 - 0.25RS^4 - RS^3T - 1.5RS^2T^2 - RST^3 - 0.25RT^4 + 0.25S^5 + 1.25S^4T + 2.5S^3T^2 + 2.5S^2T^3 + 1.25ST^4 + 0.25T^5\\big)M + \\big(- P^2R^2S^2 + 2P^2R^2ST - P^2R^2T^2 + P^2RS^3 - P^2RS^2T - P^2RST^2 + P^2RT^3 - 0.25P^2S^4 + 0.5P^2S^2T^2 - 0.25P^2T^4 + PR^2S^3 - PR^2S^2T - PR^2ST^2 + PR^2T^3 - PRS^4 + 2PRS^2T^2 - PRT^4 + 0.25PS^5 + 0.25PS^4T - 0.5PS^3T^2 - 0.5PS^2T^3 + 0.25PST^4 + 0.25PT^5 - 0.25R^2S^4 + 0.5R^2S^2T^2 - 0.25R^2T^4 + 0.25RS^5 + 0.25RS^4T - 0.5RS^3T^2 - 0.5RS^2T^3 + 0.25RST^4 + 0.25RT^5 - 0.0625S^6 - 0.125S^5T + 0.0625S^4T^2 + 0.25S^3T^3 + 0.0625S^2T^4 - 0.125ST^5 - 0.0625T^6\\big) \\\\\n",
    "    4\\bigg(-\\big(R + P\\big)^3 - \\big(P - S\\big)^3 + \\big(T - P\\big)^3 - \\big(R - S\\big)^3 + \\big(T - R\\big)^3 + 2P^3 + 2R^3 + 8PRS + 8PRT - 8PST - 8RST - 2T^3 + \\big(T + S\\big)^3\\bigg)M^3 + 2\\big(4P^3R + 8P^2R^2 - 8P^2RS - 8P^2RT - P^2S^2 - 2P^2ST - P^2T^2 + 4PR^3 - 8PR^2S - 8PR^2T + 2PRS^2 + 4PRST + 2PRT^2 + 2PS^3 + 6PS^2T + 6PST^2 + 2PT^3 - R^2S^2 - 2R^2ST - R^2T^2 + 2RS^3 + 6RS^2T + 6RST^2 + 2RT^3 - S^4 - 4S^3T - 6S^2T^2 - 4ST^3 - T^4\\big)M^2 + \\big(- 4P^3R^2 - 4P^2R^3 + 4P^2R^2S + 4P^2R^2T + 2P^2RS^2 + 4P^2RST + 2P^2RT^2 + 2PR^2S^2 + 4PR^2ST + 2PR^2T^2 - 2PRS^3 - 6PRS^2T - 6PRST^2 - 2PRT^3 - 0.25PS^4 -PS^3T - 1.5PS^2T^2 - PST^3 - 0.25PT^4 - 0.25RS^4 - RS^3T - 1.5RS^2T^2 - RST^3 - 0.25RT^4 + 0.25S^5 + 1.25S^4T + 2.5S^3T^2 + 2.5S^2T^3 + 1.25ST^4 + 0.25T^5\\big)M + \\big(- P^2R^2S^2 + 2P^2R^2ST - P^2R^2T^2 + P^2RS^3 - P^2RS^2T - P^2RST^2 + P^2RT^3 - 0.25P^2S^4 + 0.5P^2S^2T^2 - 0.25P^2T^4 + PR^2S^3 - PR^2S^2T - PR^2ST^2 + PR^2T^3 - PRS^4 + 2PRS^2T^2 - PRT^4 + 0.25PS^5 + 0.25PS^4T - 0.5PS^3T^2 - 0.5PS^2T^3 + 0.25PST^4 + 0.25PT^5 - 0.25R^2S^4 + 0.5R^2S^2T^2 - 0.25R^2T^4 + 0.25RS^5 + 0.25RS^4T - 0.5RS^3T^2 - 0.5RS^2T^3 + 0.25RST^4 + 0.25RT^5 - 0.0625S^6 - 0.125S^5T + 0.0625S^4T^2 + 0.25S^3T^3 + 0.0625S^2T^4 - 0.125ST^5 - 0.0625T^6\\big) \\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "(S - R)(S^2 - 2SR + R^2) = S^3 - 3RS^2 + 3SR^2 - R^3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 158,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-P**3 - 3*P**2*R + 3*P**2*S + 3*P**2*T - 3*P*R**2 + 8*P*R*S + 8*P*R*T - 3*P*S**2 - 8*P*S*T - 3*P*T**2 - R**3 + 3*R**2*S + 3*R**2*T - 3*R*S**2 - 8*R*S*T - 3*R*T**2 + S**3 + 3*S**2*T + 3*S*T**2 + T**3"
      ]
     },
     "execution_count": 158,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(-P**3 - 3*P**2*R + 3*P**2*S + 3*P**2*T - 3*P*R**2 + 8*P*R*S + 8*P*R*T - 3*P*S**2 - 8*P*S*T - 3*P*T**2 - R**3 + 3*R**2*S + 3*R**2*T - 3*R*S**2 - 8*R*S*T - 3*R*T**2 + S**3 + 3*S**2*T + 3*S*T**2 + T**3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 164,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "r1, r2, r3 = sym.solve(numerator, M)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "sym.cancel(sym.together(r1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 149,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.5*M*P + 0.5*M*R - 0.5*M*S - 0.5*M*T - 0.5*P*R + 0.125*S**2 + 0.25*S*T + 0.125*T**2)**3"
      ]
     },
     "execution_count": 149,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "denominator"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Denominator will be negative if and only if...\n",
    "\n",
    "\\begin{align}\n",
    "    % 4MP + 4MR - 4MS - 4MT - 4PR + S^2 + 2ST + T^2 < 0 \\\\\n",
    "    % 4M\\big((R + P) - (T + S)\\big) - 4PR + (T + S)^2< 0 \\\\\n",
    "    % 4M\\big((R + P) - (T + S)\\big) < 4RP - (T + S)^2 \\\\\n",
    "    M < \\frac{1}{4}\\frac{4RP - (T + S)^2}{(R + P) - (T + S)} \\\\\n",
    "    M < \\frac{RP - \\left(\\frac{T + S}{2}\\right)^2}{(R + P) - (T + S)} \\\\\n",
    "\\end{align}\n",
    "\n",
    "...since we are looking at interior equilibrium the RHS of the above inequality will be positive."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 139,
   "metadata": {
    "collapsed": true,
    "jupyter": {
     "outputs_hidden": true
    }
   },
   "outputs": [],
   "source": [
    "_f = sym.lambdify((T, R, P, S, M), evaluated_e3, modules = \"numpy\")\n",
    "\n",
    "def plot_non_invadable_by_randomista(T, R, P, S):\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(10,6))\n",
    "    Ms = np.linspace(0, min(0.5 * (T + S), P), 100)\n",
    "    ax.plot(Ms, _f(T, R, P, S, Ms))\n",
    "    ax.set_xlabel(r\"$M$\", fontsize=15)\n",
    "    ax.set_ylabel(r\"$e_3$\", rotation=\"horizontal\", fontsize=15)\n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 140,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "a6aaf85998844ff79ddb1c3e570a2981",
       "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",
    "w = widgets.interactive(plot_non_invadable_by_randomista, T=T_slider, R=R_slider, P=P_slider, S=S_slider)\n",
    "display(w)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 157,
   "metadata": {},
   "outputs": [
    {
     "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-157-fcea7e261b22>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m()\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msolutions\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msym\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msimplified_e3\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mM\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36msolve\u001b[0;34m(f, *symbols, **flags)\u001b[0m\n\u001b[1;32m   1063\u001b[0m     \u001b[0;31m###########################################################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1064\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mbare_f\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1065\u001b[0;31m         \u001b[0msolution\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_solve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0msymbols\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1066\u001b[0m     \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1067\u001b[0m         \u001b[0msolution\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_solve_system\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msymbols\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mflags\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36m_solve\u001b[0;34m(f, *symbols, **flags)\u001b[0m\n\u001b[1;32m   1354\u001b[0m             \u001b[0;31m# in any factor to zero\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1355\u001b[0m             \u001b[0mdens\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_simple_dens\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msymbols\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1356\u001b[0;31m             result = [s for s in result if\n\u001b[0m\u001b[1;32m   1357\u001b[0m                 all(not checksol(den, {symbol: s}, **flags) for den in\n\u001b[1;32m   1358\u001b[0m                 dens)]\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36m<listcomp>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m   1356\u001b[0m             result = [s for s in result if\n\u001b[1;32m   1357\u001b[0m                 all(not checksol(den, {symbol: s}, **flags) for den in\n\u001b[0;32m-> 1358\u001b[0;31m                 dens)]\n\u001b[0m\u001b[1;32m   1359\u001b[0m         \u001b[0;31m# set flags for quick exit at end\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1360\u001b[0m         \u001b[0mcheck\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mFalse\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36m<genexpr>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m   1355\u001b[0m             \u001b[0mdens\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_simple_dens\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0msymbols\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1356\u001b[0m             result = [s for s in result if\n\u001b[0;32m-> 1357\u001b[0;31m                 all(not checksol(den, {symbol: s}, **flags) for den in\n\u001b[0m\u001b[1;32m   1358\u001b[0m                 dens)]\n\u001b[1;32m   1359\u001b[0m         \u001b[0;31m# set flags for quick exit at end\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/solvers/solvers.py\u001b[0m in \u001b[0;36mchecksol\u001b[0;34m(f, symbol, sol, **flags)\u001b[0m\n\u001b[1;32m    245\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mflags\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'simplify'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    246\u001b[0m                 \u001b[0;32mfor\u001b[0m \u001b[0mk\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msol\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 247\u001b[0;31m                     \u001b[0msol\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mk\u001b[0m\u001b[0;34m]\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0msimplify\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msol\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[0m\n\u001b[0m\u001b[1;32m    248\u001b[0m             \u001b[0;31m# start over without the failed expanded form, possibly\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    249\u001b[0m             \u001b[0;31m# with a simplified solution\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/simplify/simplify.py\u001b[0m in \u001b[0;36msimplify\u001b[0;34m(expr, ratio, measure, fu)\u001b[0m\n\u001b[1;32m    544\u001b[0m     \u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mbottom_up\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;32mlambda\u001b[0m \u001b[0mw\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mw\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mnormal\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    545\u001b[0m     \u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mMul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mpowsimp\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_content_primitive\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[0;32m--> 546\u001b[0;31m     \u001b[0m_e\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcancel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    547\u001b[0m     \u001b[0mexpr1\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mshorter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_e\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0m_mexpand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0m_e\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcancel\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m  \u001b[0;31m# issue 6829\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    548\u001b[0m     \u001b[0mexpr2\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mshorter\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtogether\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdeep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mtogether\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr1\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdeep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mTrue\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/sympy/polys/polytools.py\u001b[0m in \u001b[0;36mcancel\u001b[0;34m(f, *gens, **args)\u001b[0m\n\u001b[1;32m   6432\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   6433\u001b[0m     \u001b[0;32mtry\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 6434\u001b[0;31m         \u001b[0;34m(\u001b[0m\u001b[0mF\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mG\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mopt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mparallel_poly_from_expr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mp\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mq\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m*\u001b[0m\u001b[0mgens\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   6435\u001b[0m     \u001b[0;32mexcept\u001b[0m \u001b[0mPolificationFailed\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   6436\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mf\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mtuple\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mTuple\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[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/polys/polytools.py\u001b[0m in \u001b[0;36mparallel_poly_from_expr\u001b[0;34m(exprs, *gens, **args)\u001b[0m\n\u001b[1;32m   4264\u001b[0m     \u001b[0;34m\"\"\"Construct polynomials from expressions. \"\"\"\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4265\u001b[0m     \u001b[0mopt\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0moptions\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mbuild_options\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mgens\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0margs\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4266\u001b[0;31m     \u001b[0;32mreturn\u001b[0m \u001b[0m_parallel_poly_from_expr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexprs\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mopt\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   4267\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4268\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/polys/polytools.py\u001b[0m in \u001b[0;36m_parallel_poly_from_expr\u001b[0;34m(exprs, opt)\u001b[0m\n\u001b[1;32m   4303\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4304\u001b[0m                 \u001b[0;32mif\u001b[0m \u001b[0mopt\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpand\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 4305\u001b[0;31m                     \u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mexpand\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   4306\u001b[0m         \u001b[0;32melse\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   4307\u001b[0m             \u001b[0mfailed\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mTrue\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36mexpand\u001b[0;34m(self, deep, modulus, power_base, power_exp, mul, log, multinomial, basic, **hints)\u001b[0m\n\u001b[1;32m   3049\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mhints\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'multinomial'\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[1;32m   3050\u001b[0m                 expr, _ = Expr._expand_hint(\n\u001b[0;32m-> 3051\u001b[0;31m                     expr, '_eval_expand_multinomial', deep=deep, **hints)\n\u001b[0m\u001b[1;32m   3052\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mhints\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mget\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m'mul'\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[1;32m   3053\u001b[0m                 expr, _ = Expr._expand_hint(\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2978\u001b[0;31m                 \u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0marghit\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mExpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_expand_hint\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m**\u001b[0m\u001b[0mhints\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2979\u001b[0m                 \u001b[0mhit\u001b[0m \u001b[0;34m|=\u001b[0m \u001b[0marghit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2980\u001b[0m                 \u001b[0msargs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mappend\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0marg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/anaconda3/lib/python3.6/site-packages/sympy/core/expr.py\u001b[0m in \u001b[0;36m_expand_hint\u001b[0;34m(expr, hint, deep, **hints)\u001b[0m\n\u001b[1;32m   2973\u001b[0m         \u001b[0;31m#              |\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2974\u001b[0m         \u001b[0;31m#              V\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 2975\u001b[0;31m         \u001b[0;32mif\u001b[0m \u001b[0mdeep\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0mgetattr\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m'args'\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mis_Atom\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2976\u001b[0m             \u001b[0msargs\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   2977\u001b[0m             \u001b[0;32mfor\u001b[0m \u001b[0marg\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0margs\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "solutions = sym.solve(simplified_e3, M)"
   ]
  },
  {
   "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
}