Untitled.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 2,
   "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": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "x1, x2, x3 = sym.symbols('x1, x2, x3', real=True, positive=True)\n",
    "x4 = 1 - x1 - x2 - x3\n",
    "T, R, P, S = sym.symbols('T, R, P, S', real=True, positive=True)\n",
    "M, epsilon = sym.symbols(\"M, epsilon\", real=True, positive=True)\n",
    "\n",
    "UGA = symbolics.UGA\n",
    "UgA = symbolics.UgA\n",
    "\n",
    "UG, Ug = sym.symbols(\"UG, Ug\", real=True, nonnegative=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# $M=0$, $\\epsilon=0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monomorphic G"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "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)\n",
    "N, = models.total_offspring(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 4.0 \\left(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*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": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(N, UGA(x1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "(f,), _, _, _ = models.offspring_genotypes_evolution(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solving for the LPE using Python\n",
    "x_star, = sym.solve(f, x1, implicit=True)\n",
    "x_star =  sym.cancel(x_star.subs({UGA(x1): UG}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{0.5 UG \\left(S + UG \\left(R - S\\right)\\right)}{0.5 P + UG^{2} \\left(0.5 P + 0.5 R - 0.5 S - 0.5 T\\right) + UG \\left(- 1.0 P + 0.5 S + 0.5 T\\right)}$"
      ],
      "text/plain": [
       "0.5*UG*(S + UG*(R - S))/(0.5*P + UG**2*(0.5*P + 0.5*R - 0.5*S - 0.5*T) + UG*(-1.0*P + 0.5*S + 0.5*T))"
      ]
     },
     "execution_count": 31,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, UG)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 32,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, UG).subs({UG: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.0$"
      ],
      "text/plain": [
       "1.00000000000000"
      ]
     },
     "execution_count": 33,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, UG).subs({UG: 1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.00000000000000, (P - S)/(P + R - S - T)]"
      ]
     },
     "execution_count": 35,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.solve(x1 - x_star.subs({UG: x1}), x1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monomorphic A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[1-x3-epsilon], [epsilon/2], [x3], [epsilon/2]])\n",
    "payoff_kernel = np.array([[R, S], [T, P]])\n",
    "W = models.generalized_sexual_selection(x, UGA, UgA, payoff_kernel)\n",
    "N, = models.total_offspring(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.0 R$"
      ],
      "text/plain": [
       "2.0*R"
      ]
     },
     "execution_count": 44,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.limit(N, epsilon, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, _, (f,), _ = models.offspring_genotypes_evolution(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.limit(f, epsilon, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[x1], [x2], [x3], [1 - x1 - x2 - x3]])\n",
    "payoff_kernel = np.array([[R, S], [T, P]])\n",
    "W = models.generalized_sexual_selection(x, UGA, UgA, payoff_kernel, M)\n",
    "\n",
    "(f1,), (f2,), (f3,), _ = models.offspring_genotypes_evolution(W, x)\n",
    "F = sym.Matrix([f1, f2, f3])\n",
    "F_jac = F.jacobian([x1, x2, x3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 112,
   "metadata": {},
   "outputs": [],
   "source": [
    "dUGdxA = sym.symbols(\"dUGdxA\", positive=True, real=True)\n",
    "dUgdxA = sym.symbols(\"dUgdxA\", positive=True, real=True)\n",
    "xA = sym.symbols(\"xA\", nonnegative=True, real=True)\n",
    "\n",
    "evaluated_F_jac = (F_jac.subs({x1 + x3: xA})\n",
    "                        .doit()\n",
    "                        .subs({sym.Derivative(UGA(xA), xA): dUGdxA, sym.Derivative(UgA(xA), xA): dUgdxA})\n",
    "                        .subs({x2: 0, xA: x1 + x3})\n",
    "                        .subs({UGA(x1 + x3): 1, UgA(x1 + x3): 1})\n",
    "                        .subs({x3: 1 - x1}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{1.0 \\left(1.0 R x_{1} \\left(1 - x_{1}\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right)\\right) \\left(2.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right) + 1.0 T dUgdxA \\left(1 - x_{1}\\right) + 0.5 dUGdxA x_{1} \\left(- 2 M + 2 T\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) - x_{1} \\left(0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) - x_{1} \\left(- 2.0 M + 2.0 R + 1.0 dUGdxA x_{1} \\left(- 2 M + 2 R\\right) + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) + 1.0 x_{1} \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) - 1.0 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} - 1.0 R x_{1} + 1.0 R + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right) - \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} + 2.0 R dUgdxA \\left(1 - x_{1}\\right) - 1.0 R x_{1} - 2.0 R \\left(1 - x_{1}\\right) + 1.0 R - 1.0 T dUgdxA + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right) + 1.0 \\left(1 - x_{1}\\right) \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right)\\right)}{\\left(R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right)^{2}} - 1 + \\frac{- 0.5 R x_{1} + 1.0 x_{1} \\left(- 2 M + 2 R\\right) + x_{1} \\left(- 2.0 M + 2.0 R + 1.0 dUGdxA x_{1} \\left(- 2 M + 2 R\\right) + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) + 1.0 x_{1} \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} - 1.0 R x_{1} + 1.0 R + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)} & \\frac{1.0 \\left(1.0 R x_{1} \\left(1 - x_{1}\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right)\\right) \\left(2.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) - 1.0 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right)}{\\left(R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right)^{2}} + \\frac{- 0.5 R x_{1} + 0.5 x_{1} \\left(- 2 M + 2 R\\right) + 0.25 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)} & \\frac{1.0 \\left(1.0 R x_{1} \\left(1 - x_{1}\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right)\\right) \\left(1.0 T dUgdxA \\left(1 - x_{1}\\right) + 0.5 dUGdxA x_{1} \\left(- 2 M + 2 T\\right) - x_{1} \\left(- 1.0 M + 1.0 R + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) - x_{1} \\left(- 1.0 M + 1.0 R + 1.0 dUGdxA x_{1} \\left(- 2 M + 2 R\\right) + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) + 1.0 x_{1} \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) - \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} - 1.0 R x_{1} + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right) - \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} + 2.0 R dUgdxA \\left(1 - x_{1}\\right) - 1.0 R x_{1} - 2.0 R \\left(1 - x_{1}\\right) + 2.0 R - 1.0 T dUgdxA + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right) + 1.0 \\left(1 - x_{1}\\right) \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right)\\right)}{\\left(R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right)^{2}} + \\frac{0.5 R x_{1} + x_{1} \\left(- 1.0 M + 1.0 R + 1.0 dUGdxA x_{1} \\left(- 2 M + 2 R\\right) + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) + 1.0 x_{1} \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} - 1.0 R x_{1} + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)}\\\\\\frac{- 0.5 R x_{1} - 0.25 dUGdxA x_{1} \\left(- 2 M + 2 T\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)} & \\frac{- 0.5 R x_{1} + 0.5 x_{1} \\left(- 2 M + 2 R\\right) + 0.25 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)} - 1 & \\frac{- 0.5 R x_{1} - 0.25 dUGdxA x_{1} \\left(- 2 M + 2 T\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)}\\\\\\frac{1.0 \\left(0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right) \\left(2.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right) + 1.0 T dUgdxA \\left(1 - x_{1}\\right) + 0.5 dUGdxA x_{1} \\left(- 2 M + 2 T\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) - x_{1} \\left(0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) - x_{1} \\left(- 2.0 M + 2.0 R + 1.0 dUGdxA x_{1} \\left(- 2 M + 2 R\\right) + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) + 1.0 x_{1} \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) - 1.0 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} - 1.0 R x_{1} + 1.0 R + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right) - \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} + 2.0 R dUgdxA \\left(1 - x_{1}\\right) - 1.0 R x_{1} - 2.0 R \\left(1 - x_{1}\\right) + 1.0 R - 1.0 T dUgdxA + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right) + 1.0 \\left(1 - x_{1}\\right) \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right)\\right)}{\\left(R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right)^{2}} + \\frac{- 0.5 R x_{1} - 1.0 R \\left(1 - x_{1}\\right) + x_{1} \\left(0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} + 2.0 R dUgdxA \\left(1 - x_{1}\\right) - 1.0 R x_{1} - 2.0 R \\left(1 - x_{1}\\right) + 1.0 R - 1.0 T dUgdxA + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right) + 1.0 \\left(1 - x_{1}\\right) \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)} & \\frac{1.0 \\left(0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right) \\left(2.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) - 1.0 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right)}{\\left(R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right)^{2}} + \\frac{- 0.5 R x_{1} - 1.0 R \\left(1 - x_{1}\\right) + 0.25 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)} & \\frac{1.0 \\left(0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right) \\left(1.0 T dUgdxA \\left(1 - x_{1}\\right) + 0.5 dUGdxA x_{1} \\left(- 2 M + 2 T\\right) - x_{1} \\left(- 1.0 M + 1.0 R + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) - x_{1} \\left(- 1.0 M + 1.0 R + 1.0 dUGdxA x_{1} \\left(- 2 M + 2 R\\right) + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 1.0 x_{1} \\left(- 2 M + 2 R\\right) + 1.0 x_{1} \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) - \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} - 1.0 R x_{1} + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right) - \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} + 2.0 R dUgdxA \\left(1 - x_{1}\\right) - 1.0 R x_{1} - 2.0 R \\left(1 - x_{1}\\right) + 2.0 R - 1.0 T dUgdxA + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right) + 1.0 \\left(1 - x_{1}\\right) \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right)\\right)}{\\left(R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)\\right)^{2}} + \\frac{0.5 R x_{1} + 1.0 R \\left(1 - x_{1}\\right) + x_{1} \\left(- 1.0 M + 1.0 R + 0.5 dUGdxA \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) - 0.25 dUGdxA \\left(- 2 M + 2 T\\right) - 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(dUGdxA \\left(- 2 M + 2 R\\right) - 2 dUGdxA \\left(- M + S\\right)\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R dUgdxA x_{1} + 2.0 R dUgdxA \\left(1 - x_{1}\\right) - 1.0 R x_{1} - 2.0 R \\left(1 - x_{1}\\right) + 2.0 R - 1.0 T dUgdxA + 0.5 x_{1} \\left(2 R dUgdxA - 2 S dUgdxA\\right) + 1.0 \\left(1 - x_{1}\\right) \\left(2 R dUgdxA - 2 S dUgdxA\\right)\\right)}{1.0 R x_{1} \\left(1 - x_{1}\\right) + 0.5 x_{1} \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right) + x_{1} \\left(1.0 x_{1} \\left(- 2 M + 2 R\\right) + 0.5 \\left(1 - x_{1}\\right) \\left(- 2 M + 2 R\\right)\\right) + \\left(1 - x_{1}\\right) \\left(1.0 R x_{1} + 2.0 R \\left(1 - x_{1}\\right)\\right)} - 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1.0*(1.0*R*x1*(1 - x1) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)))*(2.0*R*x1 + 2.0*R*(1 - x1) + 1.0*T*dUgdxA*(1 - x1) + 0.5*dUGdxA*x1*(-2*M + 2*T) - 1.0*x1*(-2*M + 2*R) - x1*(0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) - x1*(-2.0*M + 2.0*R + 1.0*dUGdxA*x1*(-2*M + 2*R) + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 1.0*x1*(-2*M + 2*R) + 1.0*x1*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S)) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) - 1.0*(1 - x1)*(-2*M + 2*R) - (1 - x1)*(1.0*R*dUgdxA*x1 - 1.0*R*x1 + 1.0*R + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA)) - (1 - x1)*(1.0*R*dUgdxA*x1 + 2.0*R*dUgdxA*(1 - x1) - 1.0*R*x1 - 2.0*R*(1 - x1) + 1.0*R - 1.0*T*dUgdxA + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA) + 1.0*(1 - x1)*(2*R*dUgdxA - 2*S*dUgdxA)))/(R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))**2 - 1 + (-0.5*R*x1 + 1.0*x1*(-2*M + 2*R) + x1*(-2.0*M + 2.0*R + 1.0*dUGdxA*x1*(-2*M + 2*R) + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 1.0*x1*(-2*M + 2*R) + 1.0*x1*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S)) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) + 0.5*(1 - x1)*(-2*M + 2*R) + (1 - x1)*(1.0*R*dUgdxA*x1 - 1.0*R*x1 + 1.0*R + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA)))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1))), 1.0*(1.0*R*x1*(1 - x1) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)))*(2.0*R*x1 + 2.0*R*(1 - x1) - 1.0*x1*(-2*M + 2*R) - 1.0*(1 - x1)*(-2*M + 2*R))/(R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))**2 + (-0.5*R*x1 + 0.5*x1*(-2*M + 2*R) + 0.25*(1 - x1)*(-2*M + 2*R))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1))),                    1.0*(1.0*R*x1*(1 - x1) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)))*(1.0*T*dUgdxA*(1 - x1) + 0.5*dUGdxA*x1*(-2*M + 2*T) - x1*(-1.0*M + 1.0*R + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) - x1*(-1.0*M + 1.0*R + 1.0*dUGdxA*x1*(-2*M + 2*R) + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 1.0*x1*(-2*M + 2*R) + 1.0*x1*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S)) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) - (1 - x1)*(1.0*R*dUgdxA*x1 - 1.0*R*x1 + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA)) - (1 - x1)*(1.0*R*dUgdxA*x1 + 2.0*R*dUgdxA*(1 - x1) - 1.0*R*x1 - 2.0*R*(1 - x1) + 2.0*R - 1.0*T*dUgdxA + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA) + 1.0*(1 - x1)*(2*R*dUgdxA - 2*S*dUgdxA)))/(R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))**2 + (0.5*R*x1 + x1*(-1.0*M + 1.0*R + 1.0*dUGdxA*x1*(-2*M + 2*R) + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 1.0*x1*(-2*M + 2*R) + 1.0*x1*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S)) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) + (1 - x1)*(1.0*R*dUgdxA*x1 - 1.0*R*x1 + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA)))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))],\n",
       "[                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            (-0.5*R*x1 - 0.25*dUGdxA*x1*(-2*M + 2*T))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1))),                                                                                                                                                                                                                                                                                                            (-0.5*R*x1 + 0.5*x1*(-2*M + 2*R) + 0.25*(1 - x1)*(-2*M + 2*R))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1))) - 1,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                            (-0.5*R*x1 - 0.25*dUGdxA*x1*(-2*M + 2*T))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))],\n",
       "[                                    1.0*(0.5*x1*(1 - x1)*(-2*M + 2*R) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))*(2.0*R*x1 + 2.0*R*(1 - x1) + 1.0*T*dUgdxA*(1 - x1) + 0.5*dUGdxA*x1*(-2*M + 2*T) - 1.0*x1*(-2*M + 2*R) - x1*(0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) - x1*(-2.0*M + 2.0*R + 1.0*dUGdxA*x1*(-2*M + 2*R) + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 1.0*x1*(-2*M + 2*R) + 1.0*x1*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S)) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) - 1.0*(1 - x1)*(-2*M + 2*R) - (1 - x1)*(1.0*R*dUgdxA*x1 - 1.0*R*x1 + 1.0*R + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA)) - (1 - x1)*(1.0*R*dUgdxA*x1 + 2.0*R*dUgdxA*(1 - x1) - 1.0*R*x1 - 2.0*R*(1 - x1) + 1.0*R - 1.0*T*dUgdxA + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA) + 1.0*(1 - x1)*(2*R*dUgdxA - 2*S*dUgdxA)))/(R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))**2 + (-0.5*R*x1 - 1.0*R*(1 - x1) + x1*(0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) + 0.5*(1 - x1)*(-2*M + 2*R) + (1 - x1)*(1.0*R*dUgdxA*x1 + 2.0*R*dUgdxA*(1 - x1) - 1.0*R*x1 - 2.0*R*(1 - x1) + 1.0*R - 1.0*T*dUgdxA + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA) + 1.0*(1 - x1)*(2*R*dUgdxA - 2*S*dUgdxA)))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1))),           1.0*(0.5*x1*(1 - x1)*(-2*M + 2*R) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))*(2.0*R*x1 + 2.0*R*(1 - x1) - 1.0*x1*(-2*M + 2*R) - 1.0*(1 - x1)*(-2*M + 2*R))/(R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))**2 + (-0.5*R*x1 - 1.0*R*(1 - x1) + 0.25*(1 - x1)*(-2*M + 2*R))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1))), 1.0*(0.5*x1*(1 - x1)*(-2*M + 2*R) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))*(1.0*T*dUgdxA*(1 - x1) + 0.5*dUGdxA*x1*(-2*M + 2*T) - x1*(-1.0*M + 1.0*R + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) - x1*(-1.0*M + 1.0*R + 1.0*dUGdxA*x1*(-2*M + 2*R) + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 1.0*x1*(-2*M + 2*R) + 1.0*x1*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S)) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) - (1 - x1)*(1.0*R*dUgdxA*x1 - 1.0*R*x1 + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA)) - (1 - x1)*(1.0*R*dUgdxA*x1 + 2.0*R*dUgdxA*(1 - x1) - 1.0*R*x1 - 2.0*R*(1 - x1) + 2.0*R - 1.0*T*dUgdxA + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA) + 1.0*(1 - x1)*(2*R*dUgdxA - 2*S*dUgdxA)))/(R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1)))**2 + (0.5*R*x1 + 1.0*R*(1 - x1) + x1*(-1.0*M + 1.0*R + 0.5*dUGdxA*(1 - x1)*(-2*M + 2*R) - 0.25*dUGdxA*(-2*M + 2*T) - 0.5*(1 - x1)*(-2*M + 2*R) + 0.5*(1 - x1)*(dUGdxA*(-2*M + 2*R) - 2*dUGdxA*(-M + S))) + (1 - x1)*(1.0*R*dUgdxA*x1 + 2.0*R*dUgdxA*(1 - x1) - 1.0*R*x1 - 2.0*R*(1 - x1) + 2.0*R - 1.0*T*dUgdxA + 0.5*x1*(2*R*dUgdxA - 2*S*dUgdxA) + 1.0*(1 - x1)*(2*R*dUgdxA - 2*S*dUgdxA)))/(1.0*R*x1*(1 - x1) + 0.5*x1*(1 - x1)*(-2*M + 2*R) + x1*(1.0*x1*(-2*M + 2*R) + 0.5*(1 - x1)*(-2*M + 2*R)) + (1 - x1)*(1.0*R*x1 + 2.0*R*(1 - x1))) - 1]])"
      ]
     },
     "execution_count": 113,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluated_F_jac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "simplified_F_jac = sym.simplify(evaluated_F_jac)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvals = simplified_F_jac.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [],
   "source": [
    "e1, e2, e3 = eigenvals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e1.subs({M: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{- 5 R + 2 T dUGdxA x_{1} - 2 T dUgdxA x_{1} + 2 T dUgdxA + \\left|{R + 2 T dUGdxA x_{1} - 2 T dUgdxA x_{1} + 2 T dUgdxA}\\right|}{8 R}$"
      ],
      "text/plain": [
       "(-5*R + 2*T*dUGdxA*x1 - 2*T*dUgdxA*x1 + 2*T*dUgdxA + Abs(R + 2*T*dUGdxA*x1 - 2*T*dUgdxA*x1 + 2*T*dUgdxA))/(8*R)"
      ]
     },
     "execution_count": 131,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(e2.subs({M: 0}), deep=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{- 5 R + 2 T dUGdxA x_{1} - 2 T dUgdxA x_{1} + 2 T dUgdxA - \\left|{R + 2 T dUGdxA x_{1} - 2 T dUgdxA x_{1} + 2 T dUgdxA}\\right|}{8 R}$"
      ],
      "text/plain": [
       "(-5*R + 2*T*dUGdxA*x1 - 2*T*dUgdxA*x1 + 2*T*dUgdxA - Abs(R + 2*T*dUGdxA*x1 - 2*T*dUgdxA*x1 + 2*T*dUgdxA))/(8*R)"
      ]
     },
     "execution_count": 130,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(e3.subs({M: 0}), deep=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 126,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(R + 2 T dUGdxA x_{1} - 2 T dUgdxA x_{1} + 2 T dUgdxA\\right)^{2}$"
      ],
      "text/plain": [
       "(R + 2*T*dUGdxA*x1 - 2*T*dUgdxA*x1 + 2*T*dUgdxA)**2"
      ]
     },
     "execution_count": 126,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(R**2 + 4*R*T*dUGdxA*x1 - 4*R*T*dUgdxA*x1 + 4*R*T*dUgdxA + 4*T**2*dUGdxA**2*x1**2 - 8*T**2*dUGdxA*dUgdxA*x1**2 + 8*T**2*dUGdxA*dUgdxA*x1 + 4*T**2*dUgdxA**2*x1**2 - 8*T**2*dUgdxA**2*x1 + 4*T**2*dUgdxA**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.array([[x1], [x2], [x3], [1 - x1 - x2 - x3]])\n",
    "payoff_kernel = np.array([[R, S], [T, P]])\n",
    "W = models.generalized_sexual_selection(x, UGA, UgA, payoff_kernel)\n",
    "N, = models.total_offspring(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.0 P$"
      ],
      "text/plain": [
       "2.0*P"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.limit(N, epsilon, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, (f,), _, _ = models.offspring_genotypes_evolution(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.limit(f, epsilon, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [],
   "source": [
    "(f1,), (f2,), (f3,), _ = models.offspring_genotypes_evolution(W, x)\n",
    "F = sym.Matrix([f1, f2, f3])\n",
    "F_jac = F.jacobian([x1, x2, x3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "dUGdxA = sym.symbols(\"dUGdxA\", positive=True, real=True)\n",
    "dUgdxA = sym.symbols(\"dUgdxA\", positive=True, real=True)\n",
    "xA = sym.symbols(\"xA\", nonnegative=True, real=True)\n",
    "\n",
    "evaluated_F_jac = (F_jac.subs({x1: 0})\n",
    "      .doit()\n",
    "      .subs({sym.Derivative(UGA(x3), x3): dUGdxA, sym.Derivative(UgA(x3), x3): dUgdxA})\n",
    "      .subs({UGA(x3): 0, UgA(x3): 0})\n",
    "      .subs({x3: 0}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{0.25 dUGdxA x_{2} \\left(- 2 M + 2 S\\right) + 0.5 x_{2} \\left(- 2 M + 2 P\\right) + 0.25 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)}{1.0 P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)} - 1 & 0 & \\frac{0.5 P x_{2} + 0.25 dUGdxA x_{2} \\left(- 2 M + 2 S\\right)}{1.0 P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)}\\\\\\frac{1.0 \\left(1.0 P x_{2} \\left(1 - x_{2}\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right)\\right) \\left(2.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right) - 1.0 S dUgdxA \\left(1 - x_{2}\\right) - 0.5 dUGdxA x_{2} \\left(- 2 M + 2 S\\right) - 1.0 x_{2} \\left(- 2 M + 2 P\\right) - x_{2} \\left(1.0 M - 1.0 P - 0.5 dUGdxA \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.25 dUGdxA \\left(- 2 M + 2 S\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right)\\right) - x_{2} \\left(1.0 M - 1.0 P - 1.0 dUGdxA x_{2} \\left(- 2 M + 2 P\\right) - 0.5 dUGdxA \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.25 dUGdxA \\left(- 2 M + 2 S\\right) + 1.0 x_{2} \\left(- 2 M + 2 P\\right) + 1.0 x_{2} \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right)\\right) - 1.0 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) - \\left(1 - x_{2}\\right) \\left(- 1.0 P dUgdxA x_{2} + 1.0 P x_{2} + 0.5 x_{2} \\left(- 2 P dUgdxA + 2 T dUgdxA\\right)\\right) - \\left(1 - x_{2}\\right) \\left(- 1.0 P dUgdxA x_{2} - 2.0 P dUgdxA \\left(1 - x_{2}\\right) + 1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right) - 2.0 P + 1.0 S dUgdxA + 0.5 x_{2} \\left(- 2 P dUgdxA + 2 T dUgdxA\\right) + 1.0 \\left(1 - x_{2}\\right) \\left(- 2 P dUgdxA + 2 T dUgdxA\\right)\\right)\\right)}{\\left(P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)\\right)^{2}} + \\frac{- 1.0 P x_{2} + 0.5 x_{2} \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 M - 1.0 P - 1.0 dUGdxA x_{2} \\left(- 2 M + 2 P\\right) - 0.5 dUGdxA \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.25 dUGdxA \\left(- 2 M + 2 S\\right) + 1.0 x_{2} \\left(- 2 M + 2 P\\right) + 1.0 x_{2} \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right)\\right) + 0.25 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + \\left(1 - x_{2}\\right) \\left(- 1.0 P dUgdxA x_{2} + 1.0 P x_{2} + 0.5 x_{2} \\left(- 2 P dUgdxA + 2 T dUgdxA\\right)\\right)}{1.0 P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)} & \\frac{1.0 \\left(1.0 P x_{2} \\left(1 - x_{2}\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right)\\right) \\left(2.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right) - 1.0 x_{2} \\left(- 2 M + 2 P\\right) - 1.0 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right)}{\\left(P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)\\right)^{2}} + \\frac{- 1.0 P x_{2} + 1.0 P \\left(1 - x_{2}\\right) + 1.5 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)}{1.0 P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)} - 1 & \\frac{1.0 \\left(1.0 P x_{2} \\left(1 - x_{2}\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right)\\right) \\left(- 1.0 S dUgdxA \\left(1 - x_{2}\\right) - 0.5 dUGdxA x_{2} \\left(- 2 M + 2 S\\right) - x_{2} \\left(1.0 M - 1.0 P - 0.5 dUGdxA \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.25 dUGdxA \\left(- 2 M + 2 S\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right)\\right) - x_{2} \\left(1.0 M - 1.0 P - 1.0 dUGdxA x_{2} \\left(- 2 M + 2 P\\right) - 0.5 dUGdxA \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.25 dUGdxA \\left(- 2 M + 2 S\\right) + 1.0 x_{2} \\left(- 2 M + 2 P\\right) + 1.0 x_{2} \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right)\\right) - \\left(1 - x_{2}\\right) \\left(- 1.0 P dUgdxA x_{2} + 1.0 P x_{2} + 0.5 x_{2} \\left(- 2 P dUgdxA + 2 T dUgdxA\\right)\\right) - \\left(1 - x_{2}\\right) \\left(- 1.0 P dUgdxA x_{2} - 2.0 P dUgdxA \\left(1 - x_{2}\\right) + 1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right) - 2.0 P + 1.0 S dUgdxA + 0.5 x_{2} \\left(- 2 P dUgdxA + 2 T dUgdxA\\right) + 1.0 \\left(1 - x_{2}\\right) \\left(- 2 P dUgdxA + 2 T dUgdxA\\right)\\right)\\right)}{\\left(P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)\\right)^{2}} + \\frac{- 0.5 P x_{2} + x_{2} \\left(1.0 M - 1.0 P - 1.0 dUGdxA x_{2} \\left(- 2 M + 2 P\\right) - 0.5 dUGdxA \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.25 dUGdxA \\left(- 2 M + 2 S\\right) + 1.0 x_{2} \\left(- 2 M + 2 P\\right) + 1.0 x_{2} \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(dUGdxA \\left(- 2 M + 2 T\\right) - 2 dUGdxA \\left(- M + P\\right)\\right)\\right) + \\left(1 - x_{2}\\right) \\left(- 1.0 P dUgdxA x_{2} + 1.0 P x_{2} + 0.5 x_{2} \\left(- 2 P dUgdxA + 2 T dUgdxA\\right)\\right)}{1.0 P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)}\\\\\\frac{1.0 S dUgdxA \\left(1 - x_{2}\\right) + 0.25 dUGdxA x_{2} \\left(- 2 M + 2 S\\right) + 0.25 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)}{1.0 P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)} & 0 & \\frac{0.5 P x_{2} + 1.0 P \\left(1 - x_{2}\\right) + 1.0 S dUgdxA \\left(1 - x_{2}\\right) + 0.25 dUGdxA x_{2} \\left(- 2 M + 2 S\\right)}{1.0 P x_{2} \\left(1 - x_{2}\\right) + 0.5 x_{2} \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right) + x_{2} \\left(1.0 x_{2} \\left(- 2 M + 2 P\\right) + 0.5 \\left(1 - x_{2}\\right) \\left(- 2 M + 2 P\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)} - 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   (0.25*dUGdxA*x2*(-2*M + 2*S) + 0.5*x2*(-2*M + 2*P) + 0.25*(1 - x2)*(-2*M + 2*P))/(1.0*P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2))) - 1,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  0,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                (0.5*P*x2 + 0.25*dUGdxA*x2*(-2*M + 2*S))/(1.0*P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2)))],\n",
       "[1.0*(1.0*P*x2*(1 - x2) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)))*(2.0*P*x2 + 2.0*P*(1 - x2) - 1.0*S*dUgdxA*(1 - x2) - 0.5*dUGdxA*x2*(-2*M + 2*S) - 1.0*x2*(-2*M + 2*P) - x2*(1.0*M - 1.0*P - 0.5*dUGdxA*(1 - x2)*(-2*M + 2*P) + 0.25*dUGdxA*(-2*M + 2*S) + 0.5*(1 - x2)*(-2*M + 2*P) + 0.5*(1 - x2)*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P))) - x2*(1.0*M - 1.0*P - 1.0*dUGdxA*x2*(-2*M + 2*P) - 0.5*dUGdxA*(1 - x2)*(-2*M + 2*P) + 0.25*dUGdxA*(-2*M + 2*S) + 1.0*x2*(-2*M + 2*P) + 1.0*x2*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P)) + 0.5*(1 - x2)*(-2*M + 2*P) + 0.5*(1 - x2)*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P))) - 1.0*(1 - x2)*(-2*M + 2*P) - (1 - x2)*(-1.0*P*dUgdxA*x2 + 1.0*P*x2 + 0.5*x2*(-2*P*dUgdxA + 2*T*dUgdxA)) - (1 - x2)*(-1.0*P*dUgdxA*x2 - 2.0*P*dUgdxA*(1 - x2) + 1.0*P*x2 + 2.0*P*(1 - x2) - 2.0*P + 1.0*S*dUgdxA + 0.5*x2*(-2*P*dUgdxA + 2*T*dUgdxA) + 1.0*(1 - x2)*(-2*P*dUgdxA + 2*T*dUgdxA)))/(P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2)))**2 + (-1.0*P*x2 + 0.5*x2*(-2*M + 2*P) + x2*(1.0*M - 1.0*P - 1.0*dUGdxA*x2*(-2*M + 2*P) - 0.5*dUGdxA*(1 - x2)*(-2*M + 2*P) + 0.25*dUGdxA*(-2*M + 2*S) + 1.0*x2*(-2*M + 2*P) + 1.0*x2*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P)) + 0.5*(1 - x2)*(-2*M + 2*P) + 0.5*(1 - x2)*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P))) + 0.25*(1 - x2)*(-2*M + 2*P) + (1 - x2)*(-1.0*P*dUgdxA*x2 + 1.0*P*x2 + 0.5*x2*(-2*P*dUgdxA + 2*T*dUgdxA)))/(1.0*P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2))), 1.0*(1.0*P*x2*(1 - x2) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)))*(2.0*P*x2 + 2.0*P*(1 - x2) - 1.0*x2*(-2*M + 2*P) - 1.0*(1 - x2)*(-2*M + 2*P))/(P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2)))**2 + (-1.0*P*x2 + 1.0*P*(1 - x2) + 1.5*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P))/(1.0*P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2))) - 1, 1.0*(1.0*P*x2*(1 - x2) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)))*(-1.0*S*dUgdxA*(1 - x2) - 0.5*dUGdxA*x2*(-2*M + 2*S) - x2*(1.0*M - 1.0*P - 0.5*dUGdxA*(1 - x2)*(-2*M + 2*P) + 0.25*dUGdxA*(-2*M + 2*S) + 0.5*(1 - x2)*(-2*M + 2*P) + 0.5*(1 - x2)*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P))) - x2*(1.0*M - 1.0*P - 1.0*dUGdxA*x2*(-2*M + 2*P) - 0.5*dUGdxA*(1 - x2)*(-2*M + 2*P) + 0.25*dUGdxA*(-2*M + 2*S) + 1.0*x2*(-2*M + 2*P) + 1.0*x2*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P)) + 0.5*(1 - x2)*(-2*M + 2*P) + 0.5*(1 - x2)*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P))) - (1 - x2)*(-1.0*P*dUgdxA*x2 + 1.0*P*x2 + 0.5*x2*(-2*P*dUgdxA + 2*T*dUgdxA)) - (1 - x2)*(-1.0*P*dUgdxA*x2 - 2.0*P*dUgdxA*(1 - x2) + 1.0*P*x2 + 2.0*P*(1 - x2) - 2.0*P + 1.0*S*dUgdxA + 0.5*x2*(-2*P*dUgdxA + 2*T*dUgdxA) + 1.0*(1 - x2)*(-2*P*dUgdxA + 2*T*dUgdxA)))/(P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2)))**2 + (-0.5*P*x2 + x2*(1.0*M - 1.0*P - 1.0*dUGdxA*x2*(-2*M + 2*P) - 0.5*dUGdxA*(1 - x2)*(-2*M + 2*P) + 0.25*dUGdxA*(-2*M + 2*S) + 1.0*x2*(-2*M + 2*P) + 1.0*x2*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P)) + 0.5*(1 - x2)*(-2*M + 2*P) + 0.5*(1 - x2)*(dUGdxA*(-2*M + 2*T) - 2*dUGdxA*(-M + P))) + (1 - x2)*(-1.0*P*dUgdxA*x2 + 1.0*P*x2 + 0.5*x2*(-2*P*dUgdxA + 2*T*dUgdxA)))/(1.0*P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2)))],\n",
       "[                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                     (1.0*S*dUgdxA*(1 - x2) + 0.25*dUGdxA*x2*(-2*M + 2*S) + 0.25*(1 - x2)*(-2*M + 2*P))/(1.0*P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2))),                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                  0,                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                                   (0.5*P*x2 + 1.0*P*(1 - x2) + 1.0*S*dUgdxA*(1 - x2) + 0.25*dUGdxA*x2*(-2*M + 2*S))/(1.0*P*x2*(1 - x2) + 0.5*x2*(1 - x2)*(-2*M + 2*P) + x2*(1.0*x2*(-2*M + 2*P) + 0.5*(1 - x2)*(-2*M + 2*P)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2))) - 1]])"
      ]
     },
     "execution_count": 88,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "evaluated_F_jac"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "eigenvals = evaluated_F_jac.eigenvals()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "e1, e2, e3 = eigenvals"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{3}{4}$"
      ],
      "text/plain": [
       "-3/4"
      ]
     },
     "execution_count": 47,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{- P + S dUGdxA x_{2} - S dUgdxA x_{2} + S dUgdxA}{2 P}$"
      ],
      "text/plain": [
       "(-P + S*dUGdxA*x2 - S*dUgdxA*x2 + S*dUgdxA)/(2*P)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(e2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 49,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "e3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Monomorphic g"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 124,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = sym.symbols(\"alpha\")\n",
    "x = np.array([[0], [0], [x3], [1 - x3]])\n",
    "payoff_kernel = np.array([[R, S], [T, P]])\n",
    "W = models.generalized_sexual_selection(x, lambda x: x**alpha, lambda x: x, payoff_kernel)\n",
    "N, = models.total_offspring(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 125,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0.5 x_{3} \\left(1 - x_{3}\\right) \\left(2 P \\left(1 - x_{3}\\right) + 2 T x_{3}\\right) + 0.5 x_{3} \\left(1 - x_{3}\\right) \\left(2 R x_{3} + 2 S \\left(1 - x_{3}\\right)\\right) + x_{3} \\left(1.0 x_{3} \\left(2 R x_{3} + 2 S \\left(1 - x_{3}\\right)\\right) + 0.5 \\left(1 - x_{3}\\right) \\left(2 P \\left(1 - x_{3}\\right) + 2 T x_{3}\\right)\\right) + \\left(1 - x_{3}\\right) \\left(0.5 x_{3} \\left(2 R x_{3} + 2 S \\left(1 - x_{3}\\right)\\right) + 1.0 \\left(1 - x_{3}\\right) \\left(2 P \\left(1 - x_{3}\\right) + 2 T x_{3}\\right)\\right)$"
      ],
      "text/plain": [
       "0.5*x3*(1 - x3)*(2*P*(1 - x3) + 2*T*x3) + 0.5*x3*(1 - x3)*(2*R*x3 + 2*S*(1 - x3)) + x3*(1.0*x3*(2*R*x3 + 2*S*(1 - x3)) + 0.5*(1 - x3)*(2*P*(1 - x3) + 2*T*x3)) + (1 - x3)*(0.5*x3*(2*R*x3 + 2*S*(1 - x3)) + 1.0*(1 - x3)*(2*P*(1 - x3) + 2*T*x3))"
      ]
     },
     "execution_count": 125,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "_, _, (f,), _ = models.offspring_genotypes_evolution(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# solving for the LPE using Python\n",
    "x_star, = sym.solve(f, x3, implicit=True)\n",
    "x_star =  sym.cancel(x_star.subs({UgA(x3): Ug}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{0.5 Ug \\left(S + Ug \\left(R - S\\right)\\right)}{0.5 P + Ug^{2} \\left(0.5 P + 0.5 R - 0.5 S - 0.5 T\\right) + Ug \\left(- 1.0 P + 0.5 S + 0.5 T\\right)}$"
      ],
      "text/plain": [
       "0.5*Ug*(S + Ug*(R - S))/(0.5*P + Ug**2*(0.5*P + 0.5*R - 0.5*S - 0.5*T) + Ug*(-1.0*P + 0.5*S + 0.5*T))"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, Ug)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 54,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, Ug).subs({Ug: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1.0$"
      ],
      "text/plain": [
       "1.00000000000000"
      ]
     },
     "execution_count": 55,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, Ug).subs({Ug: 1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[1.00000000000000, (P - S)/(P + R - S - T)]"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.solve(x3 - x_star.subs({Ug: x3}), x3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 128,
   "metadata": {},
   "outputs": [],
   "source": [
    "alpha = sym.symbols(\"alpha\", positive=True, real=True)\n",
    "x = np.array([[x1], [x2], [x3], [1 - x2 - x3]])\n",
    "payoff_kernel = np.array([[R, S], [T, P]])\n",
    "W = models.generalized_sexual_selection(x, lambda x: x**alpha, lambda x: x, payoff_kernel)\n",
    "N, = models.total_offspring(W, x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 129,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 2.0 P$"
      ],
      "text/plain": [
       "2.0*P"
      ]
     },
     "execution_count": 129,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(N.subs({x1: 0, x3: 0}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 130,
   "metadata": {},
   "outputs": [],
   "source": [
    "(f1,), (f2,), (f3,), _ = models.offspring_genotypes_evolution(W, x)\n",
    "F = sym.Matrix([f1, f2, f3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 141,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - x_{2} + \\frac{1.0 P x_{2} \\left(1 - x_{2}\\right) + x_{2} \\left(2.0 P x_{2} + 1.0 P \\left(1 - x_{2}\\right)\\right)}{2.0 P x_{2} \\left(1 - x_{2}\\right) + x_{2} \\left(2.0 P x_{2} + 1.0 P \\left(1 - x_{2}\\right)\\right) + \\left(1 - x_{2}\\right) \\left(1.0 P x_{2} + 2.0 P \\left(1 - x_{2}\\right)\\right)}$"
      ],
      "text/plain": [
       "-x2 + (1.0*P*x2*(1 - x2) + x2*(2.0*P*x2 + 1.0*P*(1 - x2)))/(2.0*P*x2*(1 - x2) + x2*(2.0*P*x2 + 1.0*P*(1 - x2)) + (1 - x2)*(1.0*P*x2 + 2.0*P*(1 - x2)))"
      ]
     },
     "execution_count": 141,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f2.subs({x1: 0, x3: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 142,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[]"
      ]
     },
     "execution_count": 142,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.solve(f2.subs({x1: 0, x3: 0}), x2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 144,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 0$"
      ],
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 144,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(f2.subs({x1: 0, x3: 0}).diff(x2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 131,
   "metadata": {},
   "outputs": [],
   "source": [
    "F_jac = F.jacobian([x1, x2, x3])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 145,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{\\infty}{\\operatorname{sign}{\\left(4.0 P^{2} x_{2} - 2.0 P S x_{2} - 2.0 P T x_{2} \\right)}}$"
      ],
      "text/plain": [
       "-oo/sign(4.0*P**2*x2 - 2.0*P*S*x2 - 2.0*P*T*x2)"
      ]
     },
     "execution_count": 145,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.limit(F_jac[0,0].subs({x3: 0}), x1, 0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\frac{1.0 R x_{2}}{xA} + \\frac{0.25 \\left(1 - x_{2}\\right) \\left(2 R \\operatorname{UgA}{\\left(xA \\right)} + 2 S \\left(1 - \\operatorname{UgA}{\\left(xA \\right)}\\right)\\right) \\operatorname{UgA}{\\left(xA \\right)}}{xA}}{0.5 x_{2} \\left(1 - x_{2}\\right) \\left(1 - \\operatorname{UgA}{\\left(xA \\right)}\\right) \\left(2 P \\left(1 - \\operatorname{UgA}{\\left(xA \\right)}\\right) + 2 T \\operatorname{UgA}{\\left(xA \\right)}\\right) + \\left(1 - x_{2}\\right) \\left(0.5 x_{2} \\left(1 - \\operatorname{UgA}{\\left(xA \\right)}\\right) \\left(2 P \\left(1 - \\operatorname{UgA}{\\left(xA \\right)}\\right) + 2 T \\operatorname{UgA}{\\left(xA \\right)}\\right) + 1.0 \\left(1 - x_{2}\\right) \\left(1 - \\operatorname{UgA}{\\left(xA \\right)}\\right) \\left(2 P \\left(1 - \\operatorname{UgA}{\\left(xA \\right)}\\right) + 2 T \\operatorname{UgA}{\\left(xA \\right)}\\right)\\right)} - 1$"
      ],
      "text/plain": [
       "(1.0*R*x2/xA + 0.25*(1 - x2)*(2*R*UgA(xA) + 2*S*(1 - UgA(xA)))*UgA(xA)/xA)/(0.5*x2*(1 - x2)*(1 - UgA(xA))*(2*P*(1 - UgA(xA)) + 2*T*UgA(xA)) + (1 - x2)*(0.5*x2*(1 - UgA(xA))*(2*P*(1 - UgA(xA)) + 2*T*UgA(xA)) + 1.0*(1 - x2)*(1 - UgA(xA))*(2*P*(1 - UgA(xA)) + 2*T*UgA(xA)))) - 1"
      ]
     },
     "execution_count": 118,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "dUGdxA = sym.symbols(\"dUGdxA\", positive=True, real=True)\n",
    "dUgdxA = sym.symbols(\"dUgdxA\", positive=True, real=True)\n",
    "xA = sym.symbols(\"xA\", nonnegative=True, real=True)\n",
    "\n",
    "(F_jac[0,0].subs({x1 + x3: xA})\n",
    "      .doit()\n",
    "      .subs({UGA(xA): 1, sym.Derivative(UGA(xA), xA): 0, sym.Derivative(UgA(xA), xA): 1})\n",
    "      .subs({x1: 0, x3: 0}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 105,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\text{NaN} & \\frac{0.5 P R}{2.0 P^{2} x_{2}^{2} - 4.0 P^{2} x_{2} + 2.0 P^{2} - 4.0 P R x_{2}^{2} + 4.0 P R x_{2} + 2.0 R^{2} x_{2}^{2}} & - \\frac{1.0 x_{2} \\left(0.5 P^{2} x_{2} - 0.5 P^{2} + 0.5 P R x_{2} - 1.0 P R - 0.5 R S x_{2} + 0.5 R S - 0.5 R T x_{2} + 0.5 R T\\right)}{2.0 P^{2} x_{2}^{2} - 4.0 P^{2} x_{2} + 2.0 P^{2} - 4.0 P R x_{2}^{2} + 4.0 P R x_{2} + 2.0 R^{2} x_{2}^{2}}\\\\\\text{NaN} & \\frac{- 1.0 P^{2} x_{2}^{2} + 2.0 P^{2} x_{2} - 1.0 P^{2} + 3.0 P R x_{2}^{2} - 4.0 P R x_{2} + 0.5 P R - 2.0 R^{2} x_{2}^{2}}{2.0 P^{2} x_{2}^{2} - 4.0 P^{2} x_{2} + 2.0 P^{2} - 4.0 P R x_{2}^{2} + 4.0 P R x_{2} + 2.0 R^{2} x_{2}^{2}} & \\frac{1.0 x_{2} \\left(1.0 P^{2} x_{2}^{2} - 1.5 P^{2} x_{2} + 0.5 P^{2} + 1.0 P R x_{2}^{2} - 2.5 P R x_{2} + 1.0 P R - 1.0 P S x_{2}^{2} + 2.0 P S x_{2} - 1.0 P S + 0.5 R S x_{2} - 0.5 R S - 1.0 R T x_{2}^{2} + 1.5 R T x_{2} - 0.5 R T\\right)}{2.0 P^{2} x_{2}^{2} - 4.0 P^{2} x_{2} + 2.0 P^{2} - 4.0 P R x_{2}^{2} + 4.0 P R x_{2} + 2.0 R^{2} x_{2}^{2}}\\\\\\text{NaN} & \\frac{0.5 P R}{2.0 P^{2} x_{2}^{2} - 4.0 P^{2} x_{2} + 2.0 P^{2} - 4.0 P R x_{2}^{2} + 4.0 P R x_{2} + 2.0 R^{2} x_{2}^{2}} & \\frac{- 1.5 P^{2} x_{2}^{2} + 2.5 P^{2} x_{2} - 1.0 P^{2} + 2.5 P R x_{2}^{2} - 2.0 P R x_{2} + 1.0 P S x_{2}^{2} - 2.0 P S x_{2} + 1.0 P S - 2.0 R^{2} x_{2}^{2} - 0.5 R S x_{2}^{2} + 0.5 R S x_{2} + 0.5 R T x_{2}^{2} - 0.5 R T x_{2}}{2.0 P^{2} x_{2}^{2} - 4.0 P^{2} x_{2} + 2.0 P^{2} - 4.0 P R x_{2}^{2} + 4.0 P R x_{2} + 2.0 R^{2} x_{2}^{2}}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[nan,                                                                                            0.5*P*R/(2.0*P**2*x2**2 - 4.0*P**2*x2 + 2.0*P**2 - 4.0*P*R*x2**2 + 4.0*P*R*x2 + 2.0*R**2*x2**2),                                                                                        -1.0*x2*(0.5*P**2*x2 - 0.5*P**2 + 0.5*P*R*x2 - 1.0*P*R - 0.5*R*S*x2 + 0.5*R*S - 0.5*R*T*x2 + 0.5*R*T)/(2.0*P**2*x2**2 - 4.0*P**2*x2 + 2.0*P**2 - 4.0*P*R*x2**2 + 4.0*P*R*x2 + 2.0*R**2*x2**2)],\n",
       "[nan, (-1.0*P**2*x2**2 + 2.0*P**2*x2 - 1.0*P**2 + 3.0*P*R*x2**2 - 4.0*P*R*x2 + 0.5*P*R - 2.0*R**2*x2**2)/(2.0*P**2*x2**2 - 4.0*P**2*x2 + 2.0*P**2 - 4.0*P*R*x2**2 + 4.0*P*R*x2 + 2.0*R**2*x2**2), 1.0*x2*(1.0*P**2*x2**2 - 1.5*P**2*x2 + 0.5*P**2 + 1.0*P*R*x2**2 - 2.5*P*R*x2 + 1.0*P*R - 1.0*P*S*x2**2 + 2.0*P*S*x2 - 1.0*P*S + 0.5*R*S*x2 - 0.5*R*S - 1.0*R*T*x2**2 + 1.5*R*T*x2 - 0.5*R*T)/(2.0*P**2*x2**2 - 4.0*P**2*x2 + 2.0*P**2 - 4.0*P*R*x2**2 + 4.0*P*R*x2 + 2.0*R**2*x2**2)],\n",
       "[nan,                                                                                            0.5*P*R/(2.0*P**2*x2**2 - 4.0*P**2*x2 + 2.0*P**2 - 4.0*P*R*x2**2 + 4.0*P*R*x2 + 2.0*R**2*x2**2),    (-1.5*P**2*x2**2 + 2.5*P**2*x2 - 1.0*P**2 + 2.5*P*R*x2**2 - 2.0*P*R*x2 + 1.0*P*S*x2**2 - 2.0*P*S*x2 + 1.0*P*S - 2.0*R**2*x2**2 - 0.5*R*S*x2**2 + 0.5*R*S*x2 + 0.5*R*T*x2**2 - 0.5*R*T*x2)/(2.0*P**2*x2**2 - 4.0*P**2*x2 + 2.0*P**2 - 4.0*P*R*x2**2 + 4.0*P*R*x2 + 2.0*R**2*x2**2)]])"
      ]
     },
     "execution_count": 105,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.simplify(F_jac.subs({x1: 0})\n",
    "            .doit()\n",
    "            .subs({UGA(x3): 1, sym.Derivative(UGA(x3), x3): 0, sym.Derivative(UgA(x3), x3): 1})\n",
    "            .subs({x3: 0}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{8.0 \\left(0.5 P \\epsilon + Ug^{2} \\left(0.5 P \\epsilon + 0.25 R - 0.25 S - 0.5 T \\epsilon\\right) + Ug \\left(- 1.0 P \\epsilon + 0.25 S + 0.5 T \\epsilon\\right)\\right)}{\\left(8.0 \\epsilon + 4.0\\right) \\left(0.5 P + Ug^{2} \\left(0.5 P + 0.5 R - 0.5 S - 0.5 T\\right) + Ug \\left(- 1.0 P + 0.5 S + 0.5 T\\right)\\right)}$"
      ],
      "text/plain": [
       "8.0*(0.5*P*epsilon + Ug**2*(0.5*P*epsilon + 0.25*R - 0.25*S - 0.5*T*epsilon) + Ug*(-1.0*P*epsilon + 0.25*S + 0.5*T*epsilon))/((8.0*epsilon + 4.0)*(0.5*P + Ug**2*(0.5*P + 0.5*R - 0.5*S - 0.5*T) + Ug*(-1.0*P + 0.5*S + 0.5*T)))"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, Ug)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{8.0 \\epsilon}{8.0 \\epsilon + 4.0}$"
      ],
      "text/plain": [
       "8.0*epsilon/(8.0*epsilon + 4.0)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, Ug).subs({Ug: 0})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{4.0}{8.0 \\epsilon + 4.0}$"
      ],
      "text/plain": [
       "4.0/(8.0*epsilon + 4.0)"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sym.factor(x_star, Ug).subs({Ug: 1})"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "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-17-93d632dfcaa2>\u001b[0m in \u001b[0;36m<module>\u001b[0;34m\u001b[0m\n\u001b[0;32m----> 1\u001b[0;31m \u001b[0msym\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msolve\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx3\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0msym\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfactor\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mx_star\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mUg\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msubs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m{\u001b[0m\u001b[0mUg\u001b[0m\u001b[0;34m:\u001b[0m \u001b[0mx3\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[0m",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/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   1172\u001b[0m     \u001b[0;31m###########################################################################\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1173\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mbare_f\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1174\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   1175\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   1176\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[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/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   1428\u001b[0m         \u001b[0;32mfor\u001b[0m \u001b[0ms\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mfailed\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1429\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-> 1430\u001b[0;31m                 \u001b[0msoln\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[0ms\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   1431\u001b[0m                 \u001b[0;32mfor\u001b[0m \u001b[0msol\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msoln\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1432\u001b[0m                     \u001b[0;32mif\u001b[0m \u001b[0mgot_s\u001b[0m \u001b[0;32mand\u001b[0m \u001b[0many\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0mss\u001b[0m \u001b[0;32min\u001b[0m \u001b[0msol\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mfree_symbols\u001b[0m \u001b[0;32mfor\u001b[0m \u001b[0mss\u001b[0m \u001b[0;32min\u001b[0m \u001b[0mgot_s\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/solvers/solvers.py\u001b[0m in \u001b[0;36m_solve\u001b[0;34m(f, *symbols, **flags)\u001b[0m\n\u001b[1;32m   1758\u001b[0m         \u001b[0;31m# if in doubt, keep it\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1759\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[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 1760\u001b[0;31m         result = [s for s in result if\n\u001b[0m\u001b[1;32m   1761\u001b[0m                   all(not checksol(d, {symbol: s}, **flags)\n\u001b[1;32m   1762\u001b[0m                     for d in dens)]\n",
      "\u001b[0;32m~/Research/evolution-of-cooperation/env/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   1760\u001b[0m         result = [s for s in result if\n\u001b[1;32m   1761\u001b[0m                   all(not checksol(d, {symbol: s}, **flags)\n\u001b[0;32m-> 1762\u001b[0;31m                     for d in dens)]\n\u001b[0m\u001b[1;32m   1763\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mcheck\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1764\u001b[0m         \u001b[0;31m# keep only results if the check is not False\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/solvers/solvers.py\u001b[0m in \u001b[0;36m<genexpr>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m   1760\u001b[0m         result = [s for s in result if\n\u001b[1;32m   1761\u001b[0m                   all(not checksol(d, {symbol: s}, **flags)\n\u001b[0;32m-> 1762\u001b[0;31m                     for d in dens)]\n\u001b[0m\u001b[1;32m   1763\u001b[0m     \u001b[0;32mif\u001b[0m \u001b[0mcheck\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   1764\u001b[0m         \u001b[0;31m# keep only results if the check is not False\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/solvers/solvers.py\u001b[0m in \u001b[0;36mchecksol\u001b[0;34m(f, symbol, sol, **flags)\u001b[0m\n\u001b[1;32m    294\u001b[0m                 \u001b[0;31m# there are free symbols -- simple expansion might work\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    295\u001b[0m                 \u001b[0m_\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mval\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--> 296\u001b[0;31m                 \u001b[0mval\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0m_mexpand\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mval\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mas_numer_denom\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m[\u001b[0m\u001b[0;36m0\u001b[0m\u001b[0;34m]\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mrecursive\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[0m\u001b[1;32m    297\u001b[0m         \u001b[0;32melif\u001b[0m \u001b[0mattempt\u001b[0m \u001b[0;34m==\u001b[0m \u001b[0;36m2\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    298\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0mminimal\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/function.py\u001b[0m in \u001b[0;36m_mexpand\u001b[0;34m(expr, recursive)\u001b[0m\n\u001b[1;32m   2788\u001b[0m     \u001b[0mwas\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;32mNone\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2789\u001b[0m     \u001b[0;32mwhile\u001b[0m \u001b[0mwas\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[0;32m-> 2790\u001b[0;31m         \u001b[0mwas\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexpr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mexpand_mul\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mexpand_multinomial\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[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m   2791\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0;32mnot\u001b[0m \u001b[0mrecursive\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2792\u001b[0m             \u001b[0;32mbreak\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/function.py\u001b[0m in \u001b[0;36mexpand_mul\u001b[0;34m(expr, deep)\u001b[0m\n\u001b[1;32m   2812\u001b[0m     \"\"\"\n\u001b[1;32m   2813\u001b[0m     return sympify(expr).expand(deep=deep, mul=True, power_exp=False,\n\u001b[0;32m-> 2814\u001b[0;31m     power_base=False, basic=False, multinomial=False, log=False)\n\u001b[0m\u001b[1;32m   2815\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   2816\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   3446\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0muse_hint\u001b[0m\u001b[0;34m:\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3447\u001b[0m                 \u001b[0mhint\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m'_eval_expand_'\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m-> 3448\u001b[0;31m                 \u001b[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhit\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[0mexpr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mhint\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mdeep\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0mdeep\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   3449\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m   3450\u001b[0m         \u001b[0;32mwhile\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/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     45\u001b[0m             \u001b[0;32mreturn\u001b[0m \u001b[0margs\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[0m\n\u001b[1;32m     46\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 47\u001b[0;31m         \u001b[0mc_part\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnc_part\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0morder_symbols\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mcls\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mflatten\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     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~/Research/evolution-of-cooperation/env/lib/python3.6/site-packages/sympy/core/mul.py\u001b[0m in \u001b[0;36mflatten\u001b[0;34m(cls, seq)\u001b[0m\n\u001b[1;32m    620\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    621\u001b[0m         \u001b[0;31m# order commutative part canonically\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 622\u001b[0;31m         \u001b[0m_mulsort\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mc_part\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    623\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    624\u001b[0m         \u001b[0;31m# current code expects coeff to be always in slot-0\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_mulsort\u001b[0;34m(args)\u001b[0m\n\u001b[1;32m     32\u001b[0m \u001b[0;32mdef\u001b[0m \u001b[0m_mulsort\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[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     33\u001b[0m     \u001b[0;31m# in-place sorting of args\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m---> 34\u001b[0;31m     \u001b[0margs\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msort\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mkey\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0m_args_sortkey\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     35\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m     36\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;36mcompare\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    209\u001b[0m             \u001b[0mr\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mBasic\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m*\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mfrozenset\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;32melse\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    210\u001b[0m             \u001b[0;32mif\u001b[0m \u001b[0misinstance\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0ml\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mBasic\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--> 211\u001b[0;31m                 \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0ml\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mcompare\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mr\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    212\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    213\u001b[0m                 \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0ml\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mr\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0ml\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mr\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;36mcompare\u001b[0;34m(self, other)\u001b[0m\n\u001b[1;32m    201\u001b[0m         \u001b[0;31m#\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    202\u001b[0m         \u001b[0mst\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_hashable_content\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--> 203\u001b[0;31m         \u001b[0mot\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0mother\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0m_hashable_content\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    204\u001b[0m         \u001b[0mc\u001b[0m \u001b[0;34m=\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mst\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m>\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mot\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m-\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mst\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m<\u001b[0m \u001b[0mlen\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mot\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    205\u001b[0m         \u001b[0;32mif\u001b[0m \u001b[0mc\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/symbol.py\u001b[0m in \u001b[0;36m_hashable_content\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    272\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_hashable_content\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\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    273\u001b[0m         \u001b[0;31m# Note: user-specified assumptions not hashed, just derived ones\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0;32m--> 274\u001b[0;31m         \u001b[0;32mreturn\u001b[0m \u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0mname\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m)\u001b[0m \u001b[0;34m+\u001b[0m \u001b[0mtuple\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0msorted\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0massumptions0\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[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n\u001b[0m\u001b[1;32m    275\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    276\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0m_eval_subs\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mold\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mnew\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/core/symbol.py\u001b[0m in \u001b[0;36massumptions0\u001b[0;34m(self)\u001b[0m\n\u001b[1;32m    282\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0massumptions0\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\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    283\u001b[0m         return dict((key, value) for key, value\n\u001b[0;32m--> 284\u001b[0;31m                 in self._assumptions.items() if value is not None)\n\u001b[0m\u001b[1;32m    285\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    286\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mcacheit\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/symbol.py\u001b[0m in \u001b[0;36m<genexpr>\u001b[0;34m(.0)\u001b[0m\n\u001b[1;32m    282\u001b[0m     \u001b[0;32mdef\u001b[0m \u001b[0massumptions0\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mself\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    283\u001b[0m         return dict((key, value) for key, value\n\u001b[0;32m--> 284\u001b[0;31m                 in self._assumptions.items() if value is not None)\n\u001b[0m\u001b[1;32m    285\u001b[0m \u001b[0;34m\u001b[0m\u001b[0m\n\u001b[1;32m    286\u001b[0m     \u001b[0;34m@\u001b[0m\u001b[0mcacheit\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n",
      "\u001b[0;31mKeyboardInterrupt\u001b[0m: "
     ]
    }
   ],
   "source": [
    "sym.solve(x3 - sym.factor(x_star, Ug).subs({Ug: x3}))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "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
}