AGI_Takeoff_dynamics.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# AGI Takeoff dynamics - Intelligence vs Quantity explosion\n",
    "## Summary\n",
    "In Galor & Weil's Unified Growth Theory (UGT), a model is given to explain the drivers for economic growth throughout history. The main insight of the model is that until relatively recently it was more worthwhile to have more kids than to invest in their education. This led to a transition from a Malthusian equilibrium to a modern growth equilibrium that's dominated more by technological growth.\n",
    "\n",
    "AGI takeoff can perhaps be modeled in a roughly similar way, to help us understand the early days of intelligence explosion. Rather than trying to understand takeoff timelines, we can try to understand whether takeoff is more likely to be initially driven by intelligence or quantity explosion.\n",
    "\n",
    "How would a \"quantity explosion\" look like? Say OpenMind has developed a powerful AI agent with broadly general human-level intelligence. It's not yet superintelligent, but it's smart enough to be able to improve itself, to acquire resources, and to understand its own situation. If it is power-seeking, it will try to acquire more resources, and it will try to improve itself. Improving itself may require getting access to a massive amount of (specific) compute and data, and has risks of value-drift or breaking in some way, so the agent decides to first try and duplicate itself many times across the internet. If it can do this, it will have a large number of copies of itself, which can then work together to acquire resources and improve themselves. \n",
    "\n",
    "This is a \"quantity explosion\" scenario, where the number of AI agents increases exponentially, but the intelligence of each agent remains roughly constant. This is in contrast to an \"intelligence explosion\" scenario, where the intelligence of each agent increases exponentially, but the number of agents remains roughly constant. In this scenario, different safety measures would be more effective (e.g. prosaic alignment and containment may be more important than solving inner alignment and understanding unbounded rational agents, and vice versa)."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model\n",
    "## readme!\n",
    "This work is unfinished, and written with the intention of self-exploration rather than to explain it to the reader. I think the models near the end of the notebook are more promising, even though they are more simplified, and that the next step here should maybe be to plot the relevant function and explore solutions manually. \n",
    "\n",
    "I currently don't intend on continuing to work on it, but I think this work may be worthwhile and I'd be very happy if someone else wants to continue exploring this. Feel free to contact me if you have any questions or comments.\n",
    "\n",
    "## Model overview\n",
    "We assume that all AI agents have a common goal and make similar decisions, but we don't account for any coordination between them. Their utility is given by a time-discounted variable that depends on efficiency-weighted effort they put directly into it, so we think of it as \"consumption\". They care about their own consumption, but also about the consumption of their descendants.\n",
    "\n",
    "The efficiency of an AI agent is given by its intelligence, which is a function of the amount of effort that went into its training. \n",
    "\n",
    "The resources they can spend are all compute-time, which itself can grow if they put effort into it.\n",
    "\n",
    "When an agent gets access to more compute, it can use it by copying itself and training the copies. They can choose to put no effort into training, in which case the copies will be identical to the original. [TODO maybe it'd be easier if the agent can train itself, rather than train its decendents?]\n",
    "\n",
    "We also assume that each action happens in discrete time steps.\n",
    "\n",
    "## Preferences\n",
    "\n",
    "The agents want to maximize the discounted sum of their utilities $u = \\sum_{t=0}^\\infty \\beta^t u_t$, where $\\beta$ is the discount factor. The utility of an agent $a$ at time step $t$ is the total consumption it gets from its own efforts and from the efforts of its descendants:\n",
    "\n",
    "$$u_t(a) = c_t(a) + \\sum_{a\\to b} u_t(b) = \\sum_{a\\leadsto z}c_t(z)$$\n",
    "\n",
    "where $c_t(a)$ is the consumption it gets from its own efforts, and $a\\to b$ means that $b$ is a child of $a$ and $a\\leadsto z$ means that $z$ is a descendant of $a$ (including $a$, its trained and untrained copies, their copies, and so on).\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Compute\n",
    "\n",
    "In order to copy itself, the agent needs to spend effort on getting access to more compute. This compute can be thought of as more servers or higher-quality hardware, so that compute is in units of operations per second (FLOPS). We assume that the amount of compute it gets is proportional to the amount of effort spent and doesn't change over time. The cost for getting access to another unit of compute, enough for one more copy, is constant.\n",
    "\n",
    "For simplicity, we assume that all agents require the same amount of compute which is constant through time. Smarter agents will be able to use it more efficiently, and if an agent has access to more compute the only way they can use it is by creating (potentially trained) copies to use it themselves. We don't account for memory or distinct types of compute; arguably we can just bundle it in with compute and it won't change the model qualitatively.\n",
    "\n",
    "Effort can thus be thought of as a fraction of the available compute. We thus normalize this and think of the available compute at each time $t$ as $1$. \n",
    "\n",
    "[\n",
    "    TODO: assume that there is a finite amount of possible available compute. Prove that if the compute is infinite then there's a problem with the optimization\n",
    "\n",
    "- We could assume that the agent has access to specific bounded amount of compute that they can share with their decendents, and that their own possible use of compute is limited.\n",
    "] \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## Intelligence\n",
    "\n",
    "Each agent has an intelligence of $I(a)$, constant through time. We think of intelligence as a multiplicative productivity factor.\n",
    "\n",
    "Training a copy of an agent, $a\\to b$, takes more effort the higher the intelligence of the copy, $I(b)$, and no effort at all if $I(b)=I(a)$. The whole training process happens through one time step exactly, and can only be done before an agent is deployed. We describe this as the following function:\n",
    "\n",
    "$$I(b) = J_{I(a)}(I(a)T(a\\to b))$$\n",
    "\n",
    "where $T(a\\to b)$ is the amount of effort $a$ spent on training $b$ and $J_L$ is a function that describe the effect of quality-weighted effort on increasing intelligence from level $L$. \n",
    "\n",
    "Reasonable assumptions:\n",
    "- $J_L(0) = L$\n",
    "- $J_L(e)$ is increasing and concave in $e$\n",
    "- $J_{J_L(e_1)}(e_2) = J_L(e_1+e_2)$ (because the amount of quality-weighted effort required to increase intelligence from $L_1$ to $L_2$ doesn't depend on the amount of quality-weighted effort that went into increasing intelligence from $L_0$ to $L_1$)\n",
    "  - Therefore, to define $J_L$ it's enough to define it for $L=0$ (which is the result of training from scratch). We'll denote this simply as $J=J_0$ and note that $J_L(e) = J(e + J^{-1}(L))$.\n",
    "  - If we assume that $J$ is increasing intelligence multiplicatively in the quality-weighted effort, in the sense of $J_L(e) = f(e)L$ for some function $f$, then $\\forall e_1,e_2.\\ f(e_1+e_2)L = J_L(e_1+e_2) = J_{J_L(e_1)}(e_2) = f(e_2)J_L(e_1) = f(e_2)f(e_1)L$, so $f(e_1+e_2) = f(e_1)f(e_2)$ and $f(e) = c^e$ for some $c$. (Note that the previous bullet point could be adjusted if we select 1 instead of 0).\n",
    "\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simpler model\n",
    "The goal of the AI is to make paperclips.\n",
    "\n",
    "There is a total of $\\bar{C}$ available compute, and each agent uses exactly $C$ compute. Starting with 1.\n",
    "\n",
    "At time $t$ the agent $a$ creates $p_t(a) = I_t(a)c_t(a)$ paperclips, where $c_t(a)$ is the amount of compute used for consumption (= creating paperclips). We assume that all agents are exactly the same and make the same decisions, so we can remove the agent indexing.\n",
    "\n",
    "The population growth, $n_{t+1} = e_tn_{t}$, is given by the amount of effort each agent puts into copying. We assume that the number of agents $n_t$ could be any real positive number, for simplicity.\n",
    "\n",
    "Intelligence growth grows polynomialy with exponent $r$, $I_{t+1} = (1 + i_t)^r I_t$, where $i_t$ is the effort agents invest in improving intelligence.\n",
    "\n",
    "Each agent uses all of their compute: $c_t + e_t + i_t = C$. \n",
    "\n",
    "The limit of available compute limits the number of agents available\n",
    "\n",
    "$$\\sum_a C \\le \\bar{C} \\implies \\forall t. n_t \\le \\frac{\\bar{C}}{C}$$\n",
    "\n",
    "The total utility is: \n",
    "$$u = \\sum_{t,a} \\beta^t p_t(a)$$\n",
    "\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Two time steps\n",
    "\n",
    "We assume $t\\in \\{0,1\\}$, $I_0 = 1$, $n_0 = 1$. \n",
    "\n",
    "At $t=1$, agents shouldn't spend any time on population or intelligence growth, so we get\n",
    "\n",
    "$$u_1 = n_1\\beta I_1C$$\n",
    "where\n",
    "$$ \\begin{align*}\n",
    "n_1 &= (1+e_0) n_0 = 1+e_0 \\\\\n",
    "I_1 &= (1+i_0)^rI_0 = (1+i_0)^r\n",
    "\\end{align*} $$\n",
    "\n",
    "At $t=0$, we get\n",
    "$$u_0 = n_0I_0(C-e_0-i_0) = C-e_0-i_0$$\n",
    "\n",
    "We need to choose $i_0,e_0$ to maximize $u = u_0+u_1$ under the constraints above.\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A maximum will be found either on the boundary or where the gradient of the utility is zero. Let's start with computing the gradient:\n",
    "\n",
    "$$ \\begin{align*}\n",
    "u &= u_0+u_1 = C-e_0-i_0 + n_1\\beta I_1 C = \\\\\n",
    "  &= C-e_0-i_0+ (1+e_0)\\beta (1+i_0)^rC\\\\\n",
    "\\frac{\\partial u}{\\partial i_0} &= -1 + (1+e_0)\\beta r(1+i_0)^{r-1}C \\\\\n",
    "\\frac{\\partial u}{\\partial e_0} &= -1 + \\beta (1+i_0)^rC\n",
    "\\end{align*}$$\n",
    "\n",
    "When the gradient is zero, we get\n",
    "\n",
    "$$ \\begin{align*}\n",
    "(1+e_0) \\beta r(1+i_0)^{r-1}C &= 1 \\\\\n",
    "\\beta(1+i_0)^rC &= 1\n",
    "\\end{align*}$$\n",
    "\n",
    "and the solution is \n",
    "\n",
    "$$ \\begin{align*}\n",
    "\n",
    "i_0 &= \\sqrt[r]{\\frac{1}{\\beta C}} - 1  \\\\\n",
    "e_0 &= \\frac{1+i_0}{r} - 1\n",
    "\n",
    "\\end{align*}$$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To check that this is indeed a maximum point, let's compute the hessian\n",
    "\n",
    "$$ \\begin{align*}\n",
    "\\frac{\\partial^2 u}{\\partial i_0^2} &= (1+e_0)\\beta r(r-1)(1+i_0)^{r-2}C \\\\\n",
    "\\frac{\\partial^2 u}{\\partial i_0\\partial e_0} &= \\beta r(1+i_0)^{r-1}C \\\\\n",
    "\\frac{\\partial^2 u}{\\partial e_0^2} &= 0\n",
    "\\end{align*} $$\n",
    "\n",
    "Oh no.. The determinant is negative, so there's a positive eigenvalue. Therefore the solution would be on the edges :(\n",
    "    \n",
    "     "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "What are the constraints?\n",
    "- $e_0 + i_0 \\le C$\n",
    "- $e_0 \\ge 0$\n",
    "- $i_0 \\ge 0$\n",
    "\n",
    "If we assume $e_0 + i_0 = C$, then we get\n",
    "$$ u = (1+e_0)\\beta (1+i_0)^rC = (1+C-i_0)\\beta(1+i_0)^rC $$\n",
    "so, to maximize this we need to maximize $(1+C-i_0)(1+i_0)^r$. The derivative is\n",
    "\n",
    "$$ \\frac{\\partial u}{\\partial i_0} = (1+C-i_0)r(1+i_0)^{r-1} - (1+i_0)^r = (1+i_0)^{r-1}((1+C-i_0)r - (1+i_0)) $$\n",
    "\n",
    "So, the maximum is at \n",
    "$$ \\begin{align*}\n",
    "\\frac{\\partial u}{\\partial i_0} &= 0 \\\\\n",
    "(1+i_0)^{r-1}((1+C-i_0)r - (1+i_0)) &=0 \\\\\n",
    "(1+C-i_0)r - (1+i_0) &= 0 \\\\\n",
    "(1+C)r -1 -(r+1)i_0 &= 0 \\\\\n",
    "i_0 &= \\frac{(1+C)r -1}{r+1} \n",
    "\\end{align*}$$\n",
    "\n",
    "[Note that the second derivative is negative, so this is indeed a maximum point]\n",
    "\n",
    "At that point, the utility is\n",
    "$$ \\begin{align*}\n",
    "u &= (1+C-i_0)\\beta(1+i_0)^rC \\\\\n",
    " &= (1+C-\\frac{(1+C)r -1}{r+1})\\beta(1+\\frac{(1+C)r -1}{r+1})^rC \\\\\n",
    " &= \\frac{(1+C)(r+1)-(1+C)r +1}{r+1}\\beta(\\frac{(2+C)r}{r+1})^rC \\\\\n",
    " &= \\beta\\frac{(1+C)(r+1)+1}{r+1}(2+C)^rC(\\frac{r}{r+1})^r \n",
    " \\end{align*}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "🤮"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[I'd bet at 1:1 odds that this is wrong]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "for $r=1$, we get $i_0 = \\frac{C}{2}$, $e_0 = \\frac{C}{2}$, $u = \\beta (1+\\frac{C}{2})^2 C$.\n",
    "\n",
    "Yep, so the above formula is slightly wrong. I don't want to fix this at the moment... "
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Notes\n",
    "* $\\beta$ is under the designer's control, so it'd be important to see how that could influence the behavior.\n",
    "* Maybe there are ways for the designer to limit the amount of compute available, $C$ or $\\bar{C}$. \n",
    "* I didn't check the constraint on $n_1$. In the case above, it might make sense to have $\\bar{C}$ large enough so that it wouldn't be possible to use all of it in one time step.\n",
    "* Next possible directions:\n",
    "  * Model costs for compute, with option to work for money. \n",
    "  * Model some diminishing returns for replication. \n",
    "  * Increase the number of steps. \n",
    "  * Solve for generic $I_0, n_0$. \n",
    "  * Find $r$, or a better way to model intelligence growth."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### intelligence growth\n",
    "\n",
    "If intelligence growth is multiplicative, in the sense that $I_{t+1} = f(i_0)I_t$, then we get\n",
    "\n",
    "$$ I_{t+2} = f(i_{t+1})I_{t+1} = f(i_{t+1})f(i_t)I_t $$\n",
    "\n",
    "And for $i_t = i_{t+1} = i$ we get\n",
    "\n",
    "$$ I_{t+2} = f(i)^2I_t $$\n",
    "\n",
    "But we assume that spending the same amount of resources in total will result in the same intelligence growth, so we get\n",
    "$$ I_{t+2}=f(2i)I_t \\implies f(i)^2 = f(2i) $$"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and generally, $f(a+b) = f(a)f(b)$, so $f$ is an exponential function. This has INCREASED marginal returns on intelligence, but constant marginal returns on its logarithm."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Intelligence growth can happen in several ways. One way would be to improve one narrow area, which doesn't have any increasing returns. Another way would be to get more compute, or to improve the ability to get more compute. The tradeoffs between them, which could be similar to tradeoffs between $i$ and $e$, could look like UGT."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Another model\n",
    "\n",
    "There is only one agent, trying to make the most paperclips. This agent can do three things: make paperclips ($c$, for \"consumption\"), improve its cognitive ability needed to make paperclips ($i$, for \"intelligence\"), and gain access to more compute ($e$, for \"expansion\").\n",
    "\n",
    "At time $t$, the agent creates $$p_t = I_t P(c_t)$$ paperclips, where $I_t$ is its intelligence level and $P$ is an increasing and concave function (say, $\\sqrt{}$ or $\\log$). The agent cares about the discounted sum of paperclips it creates, $u = \\sum_{t=0}^\\infty \\beta^t p_t$.\n",
    "\n",
    "In Jones & Vollrath (2013), equation (6.4), they suggest the following formula for the improvment in intelligence (\"skill\") from effort:\n",
    "$$ I_{t+1} = \\mu e^{\\psi i_t}I_t $$\n",
    "where $\\mu,\\psi$ are parameters. We can interpret $0<\\mu\\le 1$ as the natural degregation of intelligence, say due to model drift, so some effort is needed to maintain the same level.\n",
    "\n",
    "The amount of compute available at time $t+1$ is $$C_{t+1}=C_t + E(I_t,e_t)$$ where $E$ is an increasing and concave function in both arguments. One simple choice for $E$ could be $E(I,e)=(Ie)^\\alpha$ for some $0<\\alpha\\le 1$.\n",
    "\n",
    "At each timestep, the agent makes the decision of how to distribute its effort:\n",
    "$$ c_t + i_t + e_t = C_t. $$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Solution\n",
    "We'll assume $E(I,e) = (Ie)^\\alpha$ for some $0<\\alpha\\le 1$. Similarly, we'll assume $P(c) = c^\\gamma$ for some $0<\\gamma\\le 1$.\n",
    "### Solution for one time step\n",
    "Put everything you can into making paperclips, as the world ends.\n",
    "\n",
    "### Solution for two time steps\n",
    "At $t=1$, the agent should put all of its effort into making paperclips, so we get\n",
    "\n",
    "$$ u = \\beta p_1 + p_0 = \\beta I_1P(C_1) + I_0P(c_0). $$\n",
    "\n",
    "Now, we want to choose $c_0,i_0,e_0$ to maximize this under the constraints above. First, let's unpack $I_1$ and $C_1$:\n",
    "$$ \\begin{align*}\n",
    "I_1 &= \\mu e^{\\psi i_0}I_0 \\\\\n",
    "C_1 &= C_0 + E(I_0,e_0) = C_0 + (I_0e_0)^\\alpha\n",
    "\\end{align*}$$\n",
    "\n",
    "And therefore, we get\n",
    "$$ u = \\beta \\mu e^{\\psi i_0}I_0(C_0 + (I_0e_0)^\\alpha)^\\gamma + I_0c_0^\\gamma. $$\n",
    "\n",
    "The constraints are:\n",
    "$$ \\begin{align*}\n",
    "c_0 + i_0 + e_0 &= C_0 \\\\\n",
    "c_0 &\\ge 0 \\\\\n",
    "i_0 &\\ge 0 \\\\\n",
    "e_0 &\\ge 0\n",
    "\\end{align*}$$\n",
    "\n",
    "We can substitute $c_0 = C_0 - i_0 - e_0$ and get\n",
    "$$ \\begin{align*}\n",
    "\n",
    "u &= \\beta \\mu e^{\\psi i_0}I_0(C_0 + (I_0e_0)^\\alpha)^\\gamma + I_0(C_0 - i_0 - e_0)^\\gamma \\\\\n",
    "i_0 &\\ge 0 \\\\\n",
    "e_0 &\\ge 0 \\\\\n",
    "C_0 &\\ge i_0 + e_0\n",
    "\n",
    "\\end{align*}$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To prevent calculation errors, let's use sympy to compute the gradient and hessian of $u$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = sp.Symbol('u')\n",
    "beta = sp.Symbol('beta')\n",
    "mu = sp.Symbol('mu')\n",
    "psi = sp.Symbol('psi')\n",
    "i_0 = sp.Symbol('i_0')\n",
    "I_0 = sp.Symbol('I_0')\n",
    "C_0 = sp.Symbol('C_0')\n",
    "e_0 = sp.Symbol('e_0')\n",
    "alpha = sp.Symbol('alpha')\n",
    "gamma = sp.Symbol('gamma')\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "u = beta * mu * sp.exp(psi * i_0) * (C_0 + (I_0 * e_0)**alpha)**gamma + I_0 * (C_0 - i_0 - e_0)**gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle I_{0} \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma} + \\beta \\mu \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}$"
      ],
      "text/plain": [
       "I_0*(C_0 - e_0 - i_0)**gamma + beta*mu*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{C_{0} - e_{0} - i_{0}} + \\beta \\mu \\psi \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}$"
      ],
      "text/plain": [
       "-I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0) + beta*mu*psi*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ui = u.diff(i_0); ui"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{C_{0} - e_{0} - i_{0}} + \\frac{\\alpha \\beta \\gamma \\mu \\left(I_{0} e_{0}\\right)^{\\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)}$"
      ],
      "text/plain": [
       "-I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0) + alpha*beta*gamma*mu*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0*(C_0 + (I_0*e_0)**alpha))"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "ue = u.diff(e_0); ue"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{I_{0} \\gamma^{2} \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} + \\beta \\mu \\psi^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}$"
      ],
      "text/plain": [
       "I_0*gamma**2*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 - I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 + beta*mu*psi**2*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uii = u.diff(i_0).diff(i_0); uii"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{I_{0} \\gamma^{2} \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} + \\frac{\\alpha \\beta \\gamma \\mu \\psi \\left(I_{0} e_{0}\\right)^{\\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)}$"
      ],
      "text/plain": [
       "I_0*gamma**2*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 - I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 + alpha*beta*gamma*mu*psi*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0*(C_0 + (I_0*e_0)**alpha))"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uei = u.diff(e_0).diff(i_0); uei"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{I_{0} \\gamma^{2} \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} + \\frac{\\alpha^{2} \\beta \\gamma^{2} \\mu \\left(I_{0} e_{0}\\right)^{2 \\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0}^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{2}} - \\frac{\\alpha^{2} \\beta \\gamma \\mu \\left(I_{0} e_{0}\\right)^{2 \\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0}^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{2}} + \\frac{\\alpha^{2} \\beta \\gamma \\mu \\left(I_{0} e_{0}\\right)^{\\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0}^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)} - \\frac{\\alpha \\beta \\gamma \\mu \\left(I_{0} e_{0}\\right)^{\\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0}^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)}$"
      ],
      "text/plain": [
       "I_0*gamma**2*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 - I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 + alpha**2*beta*gamma**2*mu*(I_0*e_0)**(2*alpha)*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0**2*(C_0 + (I_0*e_0)**alpha)**2) - alpha**2*beta*gamma*mu*(I_0*e_0)**(2*alpha)*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0**2*(C_0 + (I_0*e_0)**alpha)**2) + alpha**2*beta*gamma*mu*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0**2*(C_0 + (I_0*e_0)**alpha)) - alpha*beta*gamma*mu*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0**2*(C_0 + (I_0*e_0)**alpha))"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "uee = u.diff(e_0).diff(e_0); uee"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "ename": "NotImplementedError",
     "evalue": "could not solve -I_0*e_0*gamma*(C_0 + (I_0*e_0)**alpha)*(C_0 - e_0 - (psi*(C_0 - e_0) - (gamma - 1)*LambertW(psi*exp((C_0*psi - e_0*psi - log(I_0*gamma/(beta*mu*psi*(C_0 + (I_0*e_0)**alpha)**gamma)))/(gamma - 1))/(gamma - 1)))/psi)**gamma + alpha*beta*gamma*mu*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*(C_0 - e_0 - (psi*(C_0 - e_0) - (gamma - 1)*LambertW(psi*exp((C_0*psi - e_0*psi - log(I_0*gamma/(beta*mu*psi*(C_0 + (I_0*e_0)**alpha)**gamma)))/(gamma - 1))/(gamma - 1)))/psi)*exp(psi*(C_0 - e_0) - (gamma - 1)*LambertW(psi*exp((C_0*psi - e_0*psi - log(I_0*gamma/(beta*mu*psi*(C_0 + (I_0*e_0)**alpha)**gamma)))/(gamma - 1))/(gamma - 1)))",
     "output_type": "error",
     "traceback": [
      "\u001b[1;31m---------------------------------------------------------------------------\u001b[0m",
      "\u001b[1;31mNotImplementedError\u001b[0m                       Traceback (most recent call last)",
      "\u001b[1;32mc:\\Users\\edoar\\code\\Evaluation_stuff\\AGI_Takeoff_dynamics.ipynb Cell 28\u001b[0m in \u001b[0;36m3\n\u001b[0;32m      <a href='vscode-notebook-cell:/c%3A/Users/edoar/code/Evaluation_stuff/AGI_Takeoff_dynamics.ipynb#X40sZmlsZQ%3D%3D?line=0'>1</a>\u001b[0m eqs \u001b[39m=\u001b[39m [u\u001b[39m.\u001b[39mdiff(i_0), u\u001b[39m.\u001b[39mdiff(e_0)]\n\u001b[0;32m      <a href='vscode-notebook-cell:/c%3A/Users/edoar/code/Evaluation_stuff/AGI_Takeoff_dynamics.ipynb#X40sZmlsZQ%3D%3D?line=1'>2</a>\u001b[0m \u001b[39mvars\u001b[39m \u001b[39m=\u001b[39m [i_0, e_0]\n\u001b[1;32m----> <a href='vscode-notebook-cell:/c%3A/Users/edoar/code/Evaluation_stuff/AGI_Takeoff_dynamics.ipynb#X40sZmlsZQ%3D%3D?line=2'>3</a>\u001b[0m sol \u001b[39m=\u001b[39m sp\u001b[39m.\u001b[39;49msolve(eqs, \u001b[39mvars\u001b[39;49m, \u001b[39mdict\u001b[39;49m\u001b[39m=\u001b[39;49m\u001b[39mTrue\u001b[39;49;00m)\n",
      "File \u001b[1;32mc:\\Users\\edoar\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\sympy\\solvers\\solvers.py:1147\u001b[0m, in \u001b[0;36msolve\u001b[1;34m(f, *symbols, **flags)\u001b[0m\n\u001b[0;32m   1145\u001b[0m         solution \u001b[39m=\u001b[39m _solve(f[\u001b[39m0\u001b[39m], \u001b[39m*\u001b[39msymbols, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mflags)\n\u001b[0;32m   1146\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[1;32m-> 1147\u001b[0m     linear, solution \u001b[39m=\u001b[39m _solve_system(f, symbols, \u001b[39m*\u001b[39m\u001b[39m*\u001b[39mflags)\n\u001b[0;32m   1148\u001b[0m \u001b[39massert\u001b[39;00m \u001b[39mtype\u001b[39m(solution) \u001b[39mis\u001b[39;00m \u001b[39mlist\u001b[39m\n\u001b[0;32m   1149\u001b[0m \u001b[39massert\u001b[39;00m \u001b[39mnot\u001b[39;00m solution \u001b[39mor\u001b[39;00m \u001b[39mtype\u001b[39m(solution[\u001b[39m0\u001b[39m]) \u001b[39mis\u001b[39;00m \u001b[39mdict\u001b[39m, solution\n",
      "File \u001b[1;32mc:\\Users\\edoar\\AppData\\Local\\Programs\\Python\\Python310\\lib\\site-packages\\sympy\\solvers\\solvers.py:1962\u001b[0m, in \u001b[0;36m_solve_system\u001b[1;34m(exprs, symbols, **flags)\u001b[0m\n\u001b[0;32m   1960\u001b[0m         got_s\u001b[39m.\u001b[39madd(s)\n\u001b[0;32m   1961\u001b[0m     \u001b[39mif\u001b[39;00m \u001b[39mnot\u001b[39;00m hit:\n\u001b[1;32m-> 1962\u001b[0m         \u001b[39mraise\u001b[39;00m \u001b[39mNotImplementedError\u001b[39;00m(\u001b[39m'\u001b[39m\u001b[39mcould not solve \u001b[39m\u001b[39m%s\u001b[39;00m\u001b[39m'\u001b[39m \u001b[39m%\u001b[39m eq2)\n\u001b[0;32m   1963\u001b[0m \u001b[39melse\u001b[39;00m:\n\u001b[0;32m   1964\u001b[0m     result \u001b[39m=\u001b[39m newresult\n",
      "\u001b[1;31mNotImplementedError\u001b[0m: could not solve -I_0*e_0*gamma*(C_0 + (I_0*e_0)**alpha)*(C_0 - e_0 - (psi*(C_0 - e_0) - (gamma - 1)*LambertW(psi*exp((C_0*psi - e_0*psi - log(I_0*gamma/(beta*mu*psi*(C_0 + (I_0*e_0)**alpha)**gamma)))/(gamma - 1))/(gamma - 1)))/psi)**gamma + alpha*beta*gamma*mu*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*(C_0 - e_0 - (psi*(C_0 - e_0) - (gamma - 1)*LambertW(psi*exp((C_0*psi - e_0*psi - log(I_0*gamma/(beta*mu*psi*(C_0 + (I_0*e_0)**alpha)**gamma)))/(gamma - 1))/(gamma - 1)))/psi)*exp(psi*(C_0 - e_0) - (gamma - 1)*LambertW(psi*exp((C_0*psi - e_0*psi - log(I_0*gamma/(beta*mu*psi*(C_0 + (I_0*e_0)**alpha)**gamma)))/(gamma - 1))/(gamma - 1)))"
     ]
    }
   ],
   "source": [
    "eqs = [u.diff(i_0), u.diff(e_0)]\n",
    "vars = [i_0, e_0]\n",
    "sol = sp.solve(eqs, vars, dict=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Ok, that didn't work well.. In the case where $\\alpha=\\gamma=1$, the solution is simply\n",
    "$$ \\begin{align*}\n",
    "i_0 &= \\psi^{-1}\\ln{\\frac{1}{\\beta\\mu}}  \\\\\n",
    "e_0 &= \\psi^{-1} - \\frac{C_0}{I_0}\n",
    "\\end{align*}$$\n",
    "\n",
    "Let's verify and check the Hessian.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0\n",
      "0\n"
     ]
    }
   ],
   "source": [
    "sol_i_0 = -sp.log(beta*mu)/psi\n",
    "sol_e_0 = psi**(-1) - C_0/I_0\n",
    "subs={gamma:1, alpha:1, i_0:sol_i_0, e_0 : sol_e_0}\n",
    "print(sp.simplify(ui.evalf(subs=subs)))\n",
    "print(sp.simplify(ue.evalf(subs=subs)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{I_{0} \\gamma^{2} \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} + \\beta \\mu \\psi^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi} & \\frac{I_{0} \\gamma^{2} \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} + \\frac{\\alpha \\beta \\gamma \\mu \\psi \\left(I_{0} e_{0}\\right)^{\\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)}\\\\\\frac{I_{0} \\gamma^{2} \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} + \\frac{\\alpha \\beta \\gamma \\mu \\psi \\left(I_{0} e_{0}\\right)^{\\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)} & \\frac{I_{0} \\gamma^{2} \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} - \\frac{I_{0} \\gamma \\left(C_{0} - e_{0} - i_{0}\\right)^{\\gamma}}{\\left(C_{0} - e_{0} - i_{0}\\right)^{2}} + \\frac{\\alpha^{2} \\beta \\gamma^{2} \\mu \\left(I_{0} e_{0}\\right)^{2 \\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0}^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{2}} - \\frac{\\alpha^{2} \\beta \\gamma \\mu \\left(I_{0} e_{0}\\right)^{2 \\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0}^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{2}} + \\frac{\\alpha^{2} \\beta \\gamma \\mu \\left(I_{0} e_{0}\\right)^{\\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0}^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)} - \\frac{\\alpha \\beta \\gamma \\mu \\left(I_{0} e_{0}\\right)^{\\alpha} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)^{\\gamma} e^{i_{0} \\psi}}{e_{0}^{2} \\left(C_{0} + \\left(I_{0} e_{0}\\right)^{\\alpha}\\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                                         I_0*gamma**2*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 - I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 + beta*mu*psi**2*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi),                                                                                                                                                                                                                                                                                                                                                                                            I_0*gamma**2*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 - I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 + alpha*beta*gamma*mu*psi*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0*(C_0 + (I_0*e_0)**alpha))],\n",
       "[I_0*gamma**2*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 - I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 + alpha*beta*gamma*mu*psi*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0*(C_0 + (I_0*e_0)**alpha)), I_0*gamma**2*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 - I_0*gamma*(C_0 - e_0 - i_0)**gamma/(C_0 - e_0 - i_0)**2 + alpha**2*beta*gamma**2*mu*(I_0*e_0)**(2*alpha)*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0**2*(C_0 + (I_0*e_0)**alpha)**2) - alpha**2*beta*gamma*mu*(I_0*e_0)**(2*alpha)*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0**2*(C_0 + (I_0*e_0)**alpha)**2) + alpha**2*beta*gamma*mu*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0**2*(C_0 + (I_0*e_0)**alpha)) - alpha*beta*gamma*mu*(I_0*e_0)**alpha*(C_0 + (I_0*e_0)**alpha)**gamma*exp(i_0*psi)/(e_0**2*(C_0 + (I_0*e_0)**alpha))]])"
      ]
     },
     "execution_count": 56,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "h = sp.Matrix([[uii, uei], [uei, uee]]); h"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}I_{0} \\psi & I_{0} \\psi\\\\I_{0} \\psi & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[I_0*psi, I_0*psi],\n",
       "[I_0*psi,       0]])"
      ]
     },
     "execution_count": 57,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(h.evalf(subs=subs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the hessian has negative determinant at this solution, and therefore it's a saddle point, NOT A MAXIMA."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This means that the maximum is at the boundary. The boundary is $i_0=0$ or $e_0=0$ or $C_0 = i_0 + e_0$. Let's check each of them in turn.\n",
    "\n",
    "The first case is $i_0=0$. In this case, the utility is a linear function with respect to $e_0$ and the derivative is: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle I_{0} \\beta \\mu - I_{0}$"
      ],
      "text/plain": [
       "I_0*beta*mu - I_0"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "(u.evalf(subs={gamma:1, alpha:1, i_0:0})).diff(e_0)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is interestingly greater than 0 for $\\mu=1$ (or, in fact, $\\mu>\\beta^{-1}$), so in that case the maximum would be when $c_0=0$ and all the computation is devoted to expansion. For low values of $\\mu$, where intelligence degrades, the maximum would be when $c_0=C_0$ and all the computation is devoted to consumption.\n",
    "\n",
    "When $e_0=0$, we get $u=$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C_{0} \\beta \\mu e^{i_{0} \\psi} + I_{0} \\left(C_{0} - i_{0}\\right)$"
      ],
      "text/plain": [
       "C_0*beta*mu*exp(i_0*psi) + I_0*(C_0 - i_0)"
      ]
     },
     "execution_count": 46,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u.evalf(subs={gamma:1, alpha:1, e_0:0})"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "which is a convex function of $i_0$. Therefore, the maximum is again at the boundary, which is $i_0=0$ or $i_0=C_0$ ($c_0=0$).\n",
    "\n",
    "$$ u = \\begin{cases}\n",
    "C_0(\\beta\\mu+I_0) & \\text{if } i_0=0 \\\\\n",
    "C_0\\beta\\mu e^{C_0\\psi} & \\text{if } i_0=C_0 \\\\\n",
    "\\end{cases}$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So, in this case we see that if $I_0$ is large enough we wouldn't want to put any effort into increasing intelligence, and if $\\psi$ is large enough then we should put everything into increasing intelligence."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Lastly, when $C_0 = i_0 + e_0$, we set $e_0=C_0-i_0$ and get $u=$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\beta \\mu \\left(C_{0} + I_{0} \\left(C_{0} - i_{0}\\right)\\right) e^{i_{0} \\psi}$"
      ],
      "text/plain": [
       "beta*mu*(C_0 + I_0*(C_0 - i_0))*exp(i_0*psi)"
      ]
     },
     "execution_count": 37,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u.evalf(subs={gamma:1, alpha:1, e_0:C_0-i_0 })"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with a solution at $i_0 =$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle C_{0} + \\frac{C_{0}}{I_{0}} - \\frac{1}{\\psi}$"
      ],
      "text/plain": [
       "C_0 + C_0/I_0 - 1/psi"
      ]
     },
     "execution_count": 48,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sol_i = sp.solve(u.evalf(subs={gamma:1, alpha:1, e_0:C_0-i_0 }).diff(i_0), i_0)[0]; sol_i"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If this is $>C_0$ (or equivalently, $I_0 < C_0\\psi$), or negative (for $\\psi < C_0^{-1}\\frac{I_0}{I_0-1}$) then the solution is at one of the edge cases (either $i_0=0$ or $i_0=C_0$) we've seen before. Otherwise, the solution is at the point we've found.\n",
    "\n",
    "The value of the second derivative there is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle - 1.0 I_{0} \\beta \\mu \\psi e^{\\frac{C_{0} I_{0} \\psi + C_{0} \\psi - I_{0}}{I_{0}}}$"
      ],
      "text/plain": [
       "-1.0*I_0*beta*mu*psi*exp((C_0*I_0*psi + C_0*psi - I_0)/I_0)"
      ]
     },
     "execution_count": 50,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sp.simplify(u.evalf(subs={gamma:1, alpha:1, e_0:C_0-i_0 }).diff(i_0).diff(i_0).evalf(subs={i_0:sol_i}))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This is negative, and at last we've found an interesting maximum point! The utility there is"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\beta \\mu \\left(C_{0} + I_{0} \\left(- \\frac{C_{0}}{I_{0}} + \\frac{1}{\\psi}\\right)\\right) e^{\\psi \\left(C_{0} + \\frac{C_{0}}{I_{0}} - \\frac{1}{\\psi}\\right)}$"
      ],
      "text/plain": [
       "beta*mu*(C_0 + I_0*(-C_0/I_0 + 1/psi))*exp(psi*(C_0 + C_0/I_0 - 1/psi))"
      ]
     },
     "execution_count": 52,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "u.evalf(subs={gamma:1, alpha:1, e_0:C_0-i_0, i_0:sol_i })"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Summary for two time steps with $\\alpha=\\gamma=1$\n",
    "- The utility function has no maximum points in the interior of the domain.\n",
    "- \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next steps\n",
    "- Write a short draft post to explain the key idea, and send to Niki"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "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.10.7"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}