Coase EGM.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Coase's Theory of the Firm with Endgoenous Grid"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy.optimize import fminbound\n",
    "from interpolation import interp"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Getting the Euler equation"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Recall the bellman equation for this problem\n",
    "\n",
    "$$\n",
    "p^*(s) = \\min_{t \\leq s}  \\,\n",
    "        \\{ c(s-t) + \\delta p^*(t) \\}\n",
    "       \\quad \\text{for all} \\quad\n",
    "       s \\in [0,1]. \\tag{1}\n",
    "$$\n",
    "\n",
    "We assume that $ c $ is continuously differentiable, with $ c(0)=0 $ and $ c'(0) > 0 $.\n",
    "\n",
    "we know that the range of task ($l$) is detemined by $l=s-t$.\n",
    "\n",
    "This means we can change the above bellman equation to.\n",
    "\n",
    "$$\n",
    "p^*(s) = \\min_{0\\leq l \\leq s}  \\,\n",
    "        \\{ c(l) + \\delta p^*(s-l) \\}\n",
    "       \\quad \\text{for all} \\quad\n",
    "       s \\in [0,1]. \\tag{2}\n",
    "$$\n",
    "\n",
    "This is change the firms choice from choosing $t$ to choosing $l$ but as $s$ is given when the firm chooses $t$ their effectivly choosing $l$. So this choice is interchangeable. \n",
    "\n",
    "\n",
    "let $\\sigma^*$ denote the optimal policy.\n",
    "\n",
    "When means that $\\sigma^*(s)$ is the minimizer of $p^*(s)$.\n",
    "\n",
    "If we take the deriviative of ($1$) naively with respect to $s$ we get.\n",
    "\n",
    "$$\n",
    "p'(s)=c'(s-t)\n",
    "$$\n",
    "evaluating this at the optimum of $\\sigma^*(s)$ gives\n",
    "\n",
    "$$\n",
    "(p^*)'(s)=c'(\\sigma^*(s)) = (c'\\circ \\sigma^*)(s) \\tag{3}\n",
    "$$\n",
    "\n",
    "Now taking the first order condistion of ($2$) gives.\n",
    "\n",
    "$$\n",
    "c'(\\sigma^*(s))=\\delta (p^*)'(s-\\sigma^*(s)) \n",
    "$$\n",
    "\n",
    "using the fact of ($3$) gives.\n",
    "\n",
    "$$\n",
    "c'(\\sigma^*(s))=\\delta (c'\\circ \\sigma^*)(s-\\sigma^*(s)) \n",
    "$$\n",
    "This gives us the Euler equation\n",
    "$$\n",
    "c'(\\sigma(s))=\\delta (c'\\circ \\sigma)(s-\\sigma(s)) \n",
    "$$\n",
    "\n",
    "Now New can turn this intpo a functional equation\n",
    "$$\n",
    "(c'\\circ\\sigma)(s)=\\delta (c'\\circ \\sigma)(s-\\sigma(s)) \n",
    "$$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## The Operator\n",
    "\n",
    "Our operator $K$ will act on the set of $\\sigma\\in\\Sigma$.\n",
    "\n",
    "THe set of policies will be denoted by $P$.\n",
    "\n",
    "This operator $K$ will take a policy $\\sigma\\in P$ and will return a new function $K\\sigma$\n",
    "\n",
    "This means that $K\\sigma(s)$ is the $l\\in[0,s]$ to solve the Euler equation. \n",
    "\n",
    "Now by construction K fixed point will conincided with the solution of the Euler functional equation.\n",
    "\n",
    "Indeed for a fixed $s$ $K\\sigma^*(s)$ is the optimal $l\\in[0,s]$ Euler equation.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Endogenous grid\n",
    "\n",
    "becuase we knwo that $l=s-t$ we can make our $s$ grid endogenous. Making use of the Euler equation we get.\n",
    "$$\n",
    "l = (c')^{-1}\\left(\\delta(c'\\circ \\sigma)(t)\\right)\n",
    "$$\n",
    "Then using  $l$ and $t$ we can get $s$, and now we can use the points $(s,l)$ to evaluate $K\\sigma$\n",
    "\n",
    "Now we have a path to the solution. Iterating through the process above til convergence we will get the optimal $l^*$ given $s$ for each firm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implementation\n",
    "Set up a calss ``ProductionChain`` to for our primitives. We require that our cost function be increase, continuously differentiable and strictly convex and $c(0)=0$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ProductionChain:\n",
    "\n",
    "    def __init__(self,\n",
    "            n=10000,      \n",
    "            delta=1.05 #Transaction cost muiltiplyer\n",
    "                ):\n",
    "\n",
    "        self.n, self.delta = n, delta\n",
    "        \n",
    "        # each firms t postion held in array t_grid\n",
    "        self.t_grid = np.linspace(0, 1, n)\n",
    "        \n",
    "    # The cost function of having range of tasks l\n",
    "    def c(self, x):\n",
    "        return np.exp(10*x)-1\n",
    "    # Derivative of cost function\n",
    "    def c_prime(self,x):\n",
    "        return 10*np.exp(10*x)\n",
    "    # Inverse of cost functions derivative.\n",
    "    def c_prime_inverse(self,dc):\n",
    "        return (1/10)*(np.log(dc)-np.log(10))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "pc = ProductionChain()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I need to create the operator ``K``. For a given function of the optimal policy (``σ``) determine what point ``l`` takes for each ``t``. Then determine ``s`` and update ``σ`` using the points from ``s`` and ``l``."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "def K(σ, pc):\n",
    "    \"\"\"\n",
    "    operater K to determine updated σ_new\n",
    "    σ is our current policy function\n",
    "    pc is an instance of the ProductionChain class.\n",
    "    \"\"\"\n",
    "    \n",
    "    \n",
    "    t_grid = pc.t_grid\n",
    "    \n",
    "    l = np.empty_like(t_grid)\n",
    "    ell = np.empty_like(l)\n",
    "    \n",
    "    for i,t in enumerate(t_grid):\n",
    "        x = pc.delta*pc.c_prime(σ(t))\n",
    "        l[i] = pc.c_prime_inverse(x)\n",
    "\n",
    "    \n",
    "    # update policy\n",
    "    l = np.minimum(l, t_grid)\n",
    "    s = l+  t_grid\n",
    "\n",
    "    \n",
    "    σ_new = lambda x: interp(s, l, x)\n",
    "    \n",
    "\n",
    "    \n",
    "    return σ_new"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now I'll create a function which iterates through the operator ``K`` until the updated ``σ`` function has converged."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "def solve_model(pc,tol=1e-05, max_iter=1000, verbose=True, print_skip=25):\n",
    "    \n",
    "    \n",
    "    error = tol+1\n",
    "    t_grid =pc.t_grid\n",
    "    i=0\n",
    "    σ= lambda x: x # intial guess\n",
    "    l_new = np.zeros(len(t_grid))\n",
    "    \n",
    "    while error > tol and i < max_iter and np.max(l_new)<=1:\n",
    "        \n",
    "        σ_new = K(σ, pc)\n",
    "        l = σ(t_grid)\n",
    "        l_new =  σ_new(t_grid)\n",
    "        \n",
    "        error = np.max(np.abs(l-l_new))\n",
    "        \n",
    "        if verbose and i % print_skip == 0:\n",
    "            print(f\"Error at iteration {i} is {error}.\")\n",
    "            \n",
    "        σ = σ_new\n",
    "        i+=1\n",
    "        \n",
    "    if i < max_iter:\n",
    "        print(f\"Iteration converged in {i} steps\")\n",
    "    else:\n",
    "        print(f\"Warning: iteration hit upper bound {max_iter}\")\n",
    "    \n",
    "    return σ"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Error at iteration 0 is 0.5.\n",
      "Iteration converged in 20 steps\n"
     ]
    }
   ],
   "source": [
    "σ_star = solve_model(pc)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAELCAYAAAAoUKpTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxU9b3/8dcnhF0IW8AQCIuA7AgExFZR61JcEK0bWrcW9d6qvb9fe+utvb39idrV3q5XuqBS1wpWq8a9Ciq4QlhlJyyBEAgJSyBAyDKf3x8zeiMNJEByTmbm/Xw8eDgz5zuZz9dM5j3f8z3nfM3dEREROZqUsAsQEZGmT2EhIiJ1UliIiEidFBYiIlInhYWIiNRJYSEiInUKLCzMbIKZrTGzPDO7p5bt481skZlVmdlVh2272czWxf7dHFTNIiISZUGcZ2FmzYC1wAVAAbAAuM7dV9Zo0xtoD3wPyHH352KPdwJygWzAgYXAaHfffaTX69Kli/fu3bsxuiIikrAWLlxY4u7ptW1LDaiGsUCeu28AMLOZwCTg87Bw902xbZHDnvtV4C133xXb/hYwAXjmSC/Wu3dvcnNzG7J+EZGEZ2b5R9oW1G6oTGBLjfsFscca7LlmdruZ5ZpZbnFx8XEXKiIi/yyosLBaHqvv/q96Pdfdp7t7trtnp6fXOooSEZHjFFRYFAA9a9zvARQG8FwREWkAQYXFAqC/mfUxsxbAZCCnns99E7jQzDqaWUfgwthjIiISkEAmuN29yszuIvoh3wyY4e4rzOx+INfdc8xsDPAC0BGYaGb3ufsQd99lZg8QDRyA+z+b7D4ee/fuZceOHVRWVp5gr5qetm3b0qNHD1JSdPqMiDSsQA6dDVp2drbXdjTU3r17KSoqIjMzk9atW2NW23RIfIpEImzdupWWLVvStWvXsMsRkThkZgvdPbu2bUn1FXTHjh1kZmbSpk2bhAoKgJSUFLp160ZpaWnYpYhIAgrqPIsmobKyktatW4ddRqNp3rw5VVVVYZchIiEo3neIN5Zvo1PbllwyPKPBf35ShQWQcCOKmhK5byLyz3btr+CN5dt5ZVkhH2/YScRh4ojuCgsRkWS350AFb67YzivLtvHh+p1UR5w+Xdpy57n9uHR4dwZ0O6lRXldhESd2795NaWkp7777LpMmTaJjx45hlyQiASk9WMlbK4t4ZVkh768roSriZHVqw7+M78slwzMYnNG+0fcsKCzixMKFC3n99dfJz8+nsrKS2267LeySRKQR7Suv5O1VRby6bBtz15ZQUR0hs0NrppzZh0uHd2doZuMHRE0KiziSiIc5i8j/Kq+sZs7qHby0ZCvvrCmmoipCRlorbjqjF5cMz+C0nh1Cm5tUWDRBd999N6+++irvvvvu5+dMjB49mn79+n2+G0pEEkNVdYQP1+/kpSWFvLliO2WHqkhv15Lrx2YxcUQGI3t2JCUl/INXFBZNzJo1a3jvvfcYPHgwJSUln4dFx44d6dixI7fccku4BYrICXN3Fm/ZQ86SQl5ZVkhJWQXtWqVy8bCTmXRaJuP6dqZZEwiImhQWTcyLL77Iddddx4wZM+jfv3/Y5YhIA1pXtI+XlhTy0tKtbNl1kBapKZw/qCuXjcjknFPTadW8WdglHpHCoomZM2cOX//61xkyZAjNmzcPuxwROUFb9xzk5aWFvLSkkFXb9pJi8OV+Xfg/5w3gq0O60a5VfPydKyyamGXLltG/f3/uuOOOsEsRkeO0a38Fr326jZwlhczfFL3u6cisDkydOJhLhncnvV3LkCs8dkkdFve9vIKVhXsb9TUGd2/PvROH1KttWVkZ7k5+fj7jx49v1LpEpGHtP1TF26uKeGlJIXPXFlMVcfp1PYnvXTiAy0ZkktW5TdglnpCkDoumZu/evZSUlPDjH/847FJEpB4qqiLMW1fMS0sKeWtlEQcrq+me1oopZ/Vh0ohMBmW0S5jL8CR1WNT3G39QcnNzycjIYPjw4bz11lucf/75CfNGE0kUkYizYNMuXlpayGufbmPPgUo6tmnOlaMzmXRaJqOzmsahrg0tqcOiKamsrGTWrFlMmTKFs88+m2uuuYYLLrgg7LJEhOihriu37SVnSSE5SwvZVlpOmxbNuHBwNyadlsmZ/bvQvFlir/igsGgimjdvztNPPw3A1KlTwy1GRAAo2H2AFxdv5cUlheTtKCM1xTh7QDo/uHgQ5w/qSpsWyfMRmjw9FRGph/2Hqnh9+XaeX1jARxt2AjC2Tyd+csVQLh6aQce2LUKuMBwKCxFJepGI8/GGnTy3qIA3lm/nQEU1vTq34bsXDOCKkZn07BTfRzI1BIWFiCStjSX7eX5hAS8s3srWPQdp1zKVSad158pRPRjdq6MOMKlBYSEiSaX0YCWvLCvk+YUFLNq8hxSDs/qn8/2LBnLh4G5N+pIbYVJYiEjCq6qOMG9dCc8tKuCtlUVUVEXo3/Uk7rloIFeMzKRb+1Zhl9jkJV1YuHvCDi213oXIF63evpfnFxbw4pJCivcdomOb5lw3pidXju7BsMy0hP0saAxJFRbNmzfn4MGDtGmTmJNVlZWVpKYm1a9U5J/sLDvES0sKeX5RASsK95KaYpw7sCtXjurBVwZ2pUVqYp8P0ViS6pOla9eubN26lczMTFq3bp1Q3yoikQhFRUWkpaWFXYpI4CqqIsxZXcRzC7fy7podVEWcoZntuXfiYC4b0Z3OJ8XfhfuamqQKi/bt2wNQWFhIZWVlyNU0vLZt29KlS5ewyxAJhLuzrKCU5xcVkLO0kD0HKklv15JvntmHK0f14NST24VdYkJJqrCAaGB8FhoiEn927Cvn74u28tzCAvJ2lNEiNYULB3fjytE9OKtfF1IT/LIbYUm6sBCR+FMdceauLWbmgs3MXhXdzTS6V0d+esUwLhmeQVrr+FhAKJ4pLESkySrYfYBncwv4W+4WtpWW07ltC6ac2YdrxvTklPSTwi4vqSgsRKRJqaiKMHtVETMXbGHuumIgetLcjy4dzPmDuuloppAoLESkSdhQXMasBVt4flEBJWUVZKS14ttf6c/Vo3vo2kxNgMJCREJTXlnN68u3MXP+Fj7ZuItmKcZ5A7syeWxPzh7QlWYJuIhQvAosLMxsAvA7oBnwiLv//LDtLYEngNHATuBad99kZs2BR4BRsXqfcPefBVW3iDS8Vdv2MnP+Zl5YvJW95VX06tyG/5hwKleN6kFXXXqjSQokLMysGTANuAAoABaYWY67r6zRbAqw2937mdlk4BfAtcDVQEt3H2ZmbYCVZvaMu28KonYRaRhlh6p4eWkhMxdsYemWPbRolsKEoSczeUxPxvXtnJBLkSaSoEYWY4E8d98AYGYzgUlAzbCYBEyN3X4OeMiip1g70NbMUoHWQAWwN6C6ReQEuDtLtuxh5vwtvLyskAMV1QzodhL/79LBXDEyM2kXEopHQYVFJrClxv0C4PQjtXH3KjMrBToTDY5JwDagDfAdd991+AuY2e3A7QBZWVkNXb+IHIM9Byp4YfFWZi3Ywurt+2jdvBkTR2QweWwWI3t2SKhL7SSLoMKitnfG4ZdIPVKbsUA10B3oCMwzs7c/G6V83tB9OjAdIDs7W5dfFQmYu/PJxl08M38zry/fTkVVhOE90vjpFcOYOCKDdq104lw8CyosCoCeNe73AAqP0KYgtsspDdgFXA+84e6VwA4z+wDIBjYgIqErO1TFC4sKePLjfNYWldGuVSqTx/Tk2jE9GdJdF7ZMFEGFxQKgv5n1AbYCk4mGQE05wM3AR8BVwBx3dzPbDHzFzJ4iuhtqHPDbgOoWkSNYV7SPJz/O5++LtlJ2qIphmWk8eNVwJg7vTusWWm0u0QQSFrE5iLuAN4keOjvD3VeY2f1ArrvnAI8CT5pZHtERxeTY06cBfwGWE91V9Rd3XxZE3SLyRZXVEd5aWcQTH23i4w27aNEshUuHZ3DjGb04TXMRCc0ScXW17Oxsz83NDbsMkYSxY285z8zfwl/n51O09xCZHVpzw7heXJPdQ2tFJBAzW+ju2bVt0xncIlIrd2f+xl08+XE+byzfTlXEGT8gnZ9c3otzB+rs6mSjsBCRL9h/qIoXFm/lqY/zWb19H+1bpXLzl3pzw7he9OnSNuzyJCQKCxEBIG9HGU99nM/zCwvYd6iKId3b84srh3HZiExNWIvCQiSZVVVHeHtVEU9+nM8HeTtp0SyFi4edzI1n9GZUlias5X8pLESSUPG+Q8ycv5m/zt/MttJyuqe14u6vnsq1Y3rSRRPWUguFhUiScHcW5u/miY/yeX35NiqrnbP6d2HqZUM4b2BXrV0tR6WwEElwldURXvt0G4++v5FlBaW0a5XKDeN6ccO4XlqaVOpNYSGSoEoPVjJz/mYe+3AT20rL6dulLQ9cPpQrR2XSpoX+9OXY6B0jkmA27zzAjA828mzuFg5UVHNG3878+PKhnHtqV60ZIcdNYSGSAD6bj3hk3kbeXLmdZmZcNqI73zyzD0MzdTE/OXEKC5E4VlUd4fXl23nk/Y0s3bKHtNbNueOcU7jpjN500/Kk0oAUFiJxaG95JbPmb+GxDzexdc9B+mg+QhqZ3lUicWTLrgP85YNNzFqwmf0V1Yzr24n7LhvCVwZqPkIal8JCJA4szN/No+9v4I3l20kxY+KI7kzRfIQESGEh0kRVVUd4Y8V2Hpm3kSWx+Yh/PTs6H3FymuYjJFgKC5EmZm95Jc8u2MJfPojOR/Tu3IYHJg3hytE9NB8hodE7T6SJKCk7xIz3N/LkR/nsO1TF6X06MTU2H6G1IyRsCguRkG3dc5CH527gmfmbqaiOcPGwDP51/CkM66H5CGk6FBYiIdlQXMYf313PC4u3AvC1UZn869mn0FfXa5ImSGEhErAVhaX84d31vPbpNlqmpnDDuF7cNr4vmR1ah12ayBEpLEQCkrtpF9PeyeOdNcW0a5nKHeecwje+3EfrR0hcUFiINCJ3Z+66Eqa9k8f8jbvo1LYFd3/1VG4Y14u01s3DLk+k3hQWIo0gEnHeXLGdae/msXzrXjLSWnHvxMFMHpOl9awlLiksRBpQZXWEnCWF/OHdPNYX76dPl7Y8eOVwLh+ZSYtUrUQn8UthIdIAyiur+VvuFv703ga27jnIoIz2PHT9SC4amqFzJCQhKCxETsC+8kqe/mQzj8zbSEnZIUZldeCBy4dw7qldMVNISOJQWIgch7JDVTz2wUamz93A3vIqzurfhTvPHcnpfTopJCQhKSxEjsHBimqe+GgTf3pvPbsPVHL+oG7823n9GN6jQ9iliTQqhYVIPZRXVvPM/M1Me2c9JWWHGD8gne9eMIDTeiokJDkoLESOoqIqwt8WbuGhOXlsKy1nXN9O/PGGUYzp3Sns0kQCpbAQqUVVdYQXFm/ld7PXUbD7IKN7deRXV4/gS/26hF2aSCgUFiI1VEecV5YV8ru317GhZD/DMtN44PKhnDMgXRPXktQCCwszmwD8DmgGPOLuPz9se0vgCWA0sBO41t03xbYNB/4MtAciwBh3Lw+qdkl8n51x/Zu317K2qIyBJ7fjzzeO5sLB3RQSIgQUFmbWDJgGXAAUAAvMLMfdV9ZoNgXY7e79zGwy8AvgWjNLBZ4CbnT3pWbWGagMom5JfO7OnNU7+PVba1lRuJe+6W35n+tGcsmwDFJ0Mp3I54IaWYwF8tx9A4CZzQQmATXDYhIwNXb7OeAhi36luxBY5u5LAdx9Z0A1SwJzd97PK+FX/1jLki17yOrUhl9dPYJJp3UntZkuyyFyuKDCIhPYUuN+AXD6kdq4e5WZlQKdgQGAm9mbQDow090fbPySJVF9smEnv/rHWuZv2kX3tFb8/GvDuHJ0D5orJESOKKiwqG087/VskwqcCYwBDgCzzWyhu8/+wpPNbgduB8jKyjrhgiXxrCgs5eevr2beuhK6tmvJfZcNYfLYnrRM1VVgReoSVFgUAD1r3O8BFB6hTUFsniIN2BV7/D13LwEws9eAUcAXwsLdpwPTAbKzsw8PIkli20vL+e9/rOH5RQWktW7ODy8exI1n9KJVc4WESH0FFRYLgP5m1gfYCkwGrj+sTQ5wM/ARcBUwx90/2/30H2bWBqgAzgZ+E1DdEsf2H6riz++tZ/q8DUQicNtZfbnz3H5adEjkOAQSFrE5iLuAN4keOjvD3VeY2f1ArrvnAI8CT5pZHtERxeTYc3eb2a+JBo4Dr7n7q0HULfGpOuI8m7uFX/1jLSVlh7h0eAbfnzCQnp3ahF2aSNwy98TbY5Odne25ublhlyEheHfNDn722mrWFO1jdK+O/PCSQYzK6hh2WSJxITYfnF3bNp3BLQlh1ba9/PS1VcxbV0Kvzm3449dHMWHoyTqhTqSBKCwkrhXtLedX/1jD3xYW0L5Vc3506WBuHNdLS5iKNDCFhcSl/YeqmD53A9PnbqAqEmHKl/vw7a/0J62NJq9FGoPCQuJKdcR5bmF08nrHvkNcMjyD7391IFmdNXkt0pgUFhI35q4t5qevrWL19n2MyurAH28YzehemrwWCYLCQpq89cVl3PfySuauLaZnp9ZMu34UFw/T5LVIkBQW0mQdrKjmf+as4+F5G2jVvBn/dUn0zGtdnkMkeAoLaXLcnbdWFnHfyyvZuucgXxuVyQ8uGkR6u5ZhlyaStBQW0qRs3nmAqS+vYM7qHQzodhKzbh/H6X07h12WSNJTWEiTUF5ZzfS5G5j2Th6pKcYPLx7ELV/urcuGizQRCgsJ3Xtri7n3peVs2nmAS4Zl8F+XDiIjrXXYZYlIDQoLCU3hnoM88MpKXl++nT5d2vLEN8cyfkB62GWJSC0UFhK4yuoIM97fyO9mr6M64vz7BQO4/ey+OspJpAlTWEigPt6wkx+9uJx1O8o4f1BX7p04RJcOF4kDCgsJxI595fz01VW8uKSQzA6tefimbC4Y3C3sskSknhQW0qgiEeepT/L55RtrKK+q5q5z+3Hnuf1o3UK7nETiicJCGk3+zv3c/dwy5m/cxZf7deb+SUM5Jf2ksMsSkeOgsJAGF4k4T3y0iV+8sYbUFOPBK4dzdXYPXctJJI4dNSzMrAvQBVjn7tXBlCTxrOZo4uwB6fzsa8Po3kHnTIjEu7pGFu8BK4BTzWwF8Ad3fx/AzM4HBrj7Hxq5RokDGk2IJLa6wuI1YDwwA9gFXGNmdwM/dPe3zexHgMIiyeXv3M9/PLeMTzSaEElYRw0Ld7/bzAYAlwNfB4YDBiw2s0Lg0cYvUZoqjSZEkkedE9zuvhZ4MPYPM2sLtAfGAFc3anXSZGk0IZJcjvloKHffD+wHcszstYYvSZoyjSZEktMJHTrr7lUNVYg0fZt3HuDu55ZqNCGShHSehdTJ3Xny43x+9tpqjSZEkpTCQo5q9/4Kvve3pcxevYPxA9L5uUYTIklJYSFHtGDTLv7tmcXsLKtg6sTB3Pyl3hpNiCQphYX8k0jE+eN76/n1W2vp0bE1z3/rSwzrkRZ2WSISIoWFfMGOfeV8d9ZS3s8rYeKI7vz0iqG0a9U87LJEJGQKC/nc++tK+L+zlrCvvJKff20Y147pqd1OIgIoLASoqo7w27fXMe3dPPqln8TTt57OqSe3C7ssEWlCFBZJblvpQf7tmcUs2LSba7J7MPWyIbRpobeFiHxRSlAvZGYTzGyNmeWZ2T21bG9pZrNi2z8xs96Hbc8yszIz+15QNSe62auKuOh381hZuJffXnsaD141QkEhIrUK5JPBzJoB04ALgAJggZnluPvKGs2mALvdvZ+ZTQZ+AVxbY/tvgNeDqDfRVVRFePCN1Tzy/kYGZ7TnoetH0lcr2InIUQT1NXIskOfuGwDMbCYwCagZFpOAqbHbzwEPmZm5u5vZ5cAGotekkhOwZdcB7vrrIpYWlHLTGb34z4sH0aq51sMWkaMLKiwygS017hcApx+pjbtXmVkp0NnMDgLfJzoqOeIuKDO7HbgdICsrq+EqTyAfri/hW08tIuLOn24YxYShGWGXJCJxIqg5i9qOv/R6trkP+I27lx3tBdx9urtnu3t2enr6cZaZuJ7+JJ+bHp1PeruWvPLtMxUUInJMghpZFAA9a9zvARQeoU2BmaUCaURX5zsduMrMHgQ6ABEzK3f3hxq/7PhXVR3hx6+u4rEPN3HOqen8/rqRtNdJdiJyjIIKiwVAfzPrA2wFJgPXH9YmB7gZ+Ai4Cpjj7g6c9VkDM5sKlCko6qf0YCXffmYxc9cWM+XMPvznxYNolqKT7ETk2AUSFrE5iLuAN4FmwAx3X2Fm9wO57p5DdInWJ80sj+iIYnIQtSWqTSX7mfL4AvJ3HuDnXxvG5LGaxxGR42fRL++JJTs723Nzc8MuIzSfTWSbwR+/PpozTukcdkkiEgfMbKG7Z9e2TWdgJZinP8nn3pdW0LtLWx69OZtenduGXZKIJACFRYLQRLaINCaFRQLQRLaINDaFRZzbVLKfbz6+gM2ayBaRRqSwiGPLt5Zy04z5uDtP3Xo64/pqIltEGofCIk4tzN/FLTMW0L51c5669XT6dNFEtog0HoVFHPogr4RbH8/l5LRWPHXr6WR2aB12SSKS4BQWcebtlUXc8ddF9OnclidvHUvXdq3CLklEkoDCIo68vLSQ78xawuDu7Xn8G2Pp2LZF2CWJSJJQWMSJZxds4Z6/LyO7VycevSWbdjqHQkQCpLCIA3/5YCP3vbySs/p3YfqN2bRuocWKRCRYCosmbto7efzyzTV8dUg3fn/dSFqmKihEJHgKiybK3fnlm2v4w7vrufy07vz31SNIbRbUWlUiIl+ksGiCIhHn/ldW8tiHm7hubBY/uXwoKbp8h4iESGHRxLg79728gsc/yufWM/vww0sGYaagEJFwKSyamGnv5PH4R/ncdlb0goAKChFpCrQTvAmZtWAz//2PtVwxMpMfXKSgEJGmQ2HRRLy9sogf/P1Txg9I58GrhmuOQkSaFIVFE7Awfxd3/nURwzLT+OPXR9FcRz2JSBOjT6WQrS3axzcfy6V7h9bMuGUMbVtqGklEmh6FRYgK9xzk5hnzaZGawhPfHEvnk1qGXZKISK0UFiHZc6CCm2fMZ195FY99Yww9O7UJuyQRkSNSWISgvLKaWx/PJX/nAabfNJoh3dPCLklE5Ki0gzxgVdUR7vrrYhZu3s1D143iS6d0CbskEZE6aWQRsAdeWcnbq4qYOnEIlwzPCLscEZF6UVgE6OWlhTz+UT5TzuzDzV/qHXY5IiL1prAIyMaS/fzg758yKqsD91w0MOxyRESOicIiAOWV1dzx9CJSmxkPXa+T7kQk/miCOwD3vbySVdv2MuOWbLp3aB12OSIix0xfcRvZS0u28sz8zfzL2X35ysBuYZcjInJcFBaNaH1xGf/590/J7tWR7114atjliIgcN4VFIzlYUc2dTy+iZfNm/M/1IzVPISJxLbBPMDObYGZrzCzPzO6pZXtLM5sV2/6JmfWOPX6BmS00s09j//1KUDWfiKk5K1i9fR+/vmYEGWmapxCR+BZIWJhZM2AacBEwGLjOzAYf1mwKsNvd+wG/AX4Re7wEmOjuw4CbgSeDqPlE/H1RAbNyt3Dnuadwzqldwy5HROSEBTWyGAvkufsGd68AZgKTDmszCXg8dvs54DwzM3df7O6FscdXAK3MrMlennVjyX5++MJyxvbpxHfOHxB2OSIiDSKosMgEttS4XxB7rNY27l4FlAKdD2tzJbDY3Q8d/gJmdruZ5ZpZbnFxcYMVfizcnf/30nJSU4zfTx5JquYpRCRBBPVpVtsaoX4sbcxsCNFdU/9S2wu4+3R3z3b37PT09OMu9ES8sXw789aV8N0LB3ByWqtQahARaQxBhUUB0LPG/R5A4ZHamFkqkAbsit3vAbwA3OTu6xu92uOw/1AV97+ykkEZ7blxXK+wyxERaVBBhcUCoL+Z9TGzFsBkIOewNjlEJ7ABrgLmuLubWQfgVeAH7v5BQPUes9/PWce20nJ+fPkQ7X4SkYQTyKdabA7iLuBNYBXwrLuvMLP7zeyyWLNHgc5mlgd8F/js8Nq7gH7Aj8xsSexfkzrEKG/HPh6dt5GrR/dgdK9OYZcjItLgzP3wqYP4l52d7bm5uYG8lrtz/cOfsKKwlHe+d47W0RaRuGVmC909u7Zt2l9ygnKWFvLRhp3cPWGggkJEEpbC4gTsK6/kJ6+uYniPNK4fmxV2OSIijUaXKD8Bv317HcVlh3j4pmyapdR25K+ISGLQyOI4rS3ax2MfbuK6sVmM6Nkh7HJERBqVwuI4/eGdPFqlpnC3Lj0uIklAYXEcCvcc5OVl25g8NouObVuEXY6ISKNTWByHGe9vBOAbX+4dbiEiIgFRWByj0oOVPDN/M5cOz6BHxzZhlyMiEgiFxTF6Zv5m9ldUc9tZfcMuRUQkMAqLY1BRFeEvH2zky/06MzQzLexyREQCo7A4BjlLCynae4jbx58SdikiIoFSWNSTu/Pw3A0MPLkd4/t3CbscEZFAKSzq6b21xawp2sdtZ/XFTGdri0hyUVjU08PzNnBy+1ZMHNE97FJERAKnsKiHor3lfJC3k+tPz6JFqv6XiUjy0SdfPcxZvQOAC4d0C7kSEZFwKCzqYfaqIjI7tObUbu3CLkVEJBQKizocrKhm3roSLhjcTRPbIpK0FBZ1+CCvhENVEc4b1KSW/RYRCZTCog6zVxdxUstUTu/TOexSRERCo7A4ikjEmb1qB+MHdNFRUCKS1PQJeBTLC0vZse8Q5w3UUVAiktwUFkfx9qodpBicO1DzFSKS3BQWRzF7VRGjsjrSSavhiUiSU1gcwbbSg6wo3Mt5g7QLSkREYXEEs1dFz9o+X4fMiogoLI5k9qoisjq1oV/Xk8IuRUQkdAqLWhyoqOKD9Ts5b1BXnbUtIoLColbz1pVQURXhfM1XiIgACotazV5VRLuWqYzp3SnsUkREmgSFxWEiEWfO6mLOPjVdZ22LiMTo0/AwSwv2UFJ2SLugRERqCCwszGyCma0xszwzu6eW7S3NbFZs+ydm1rvGth/EHl9jZl9tzDpnr9pBsxTjnFPTG/NlRETiSiBhYWbNgGnARcBg4DozG6nn5EcAAAYmSURBVHxYsynAbnfvB/wG+EXsuYOBycAQYALwh9jPaxRvrypidK+OdGijs7ZFRD4T1MhiLJDn7hvcvQKYCUw6rM0k4PHY7eeA8yx63OokYKa7H3L3jUBe7Oc1uILdB1i9fZ9OxBMROUxQYZEJbKlxvyD2WK1t3L0KKAU61/O5mNntZpZrZrnFxcXHVWR5ZTUXDO6m+QoRkcMEFRa1ndnm9WxTn+fi7tPdPdvds9PTj2++oV/Xdjx8UzZ903XWtohITUGFRQHQs8b9HkDhkdqYWSqQBuyq53NFRKQRBRUWC4D+ZtbHzFoQnbDOOaxNDnBz7PZVwBx399jjk2NHS/UB+gPzA6pbRESA1CBexN2rzOwu4E2gGTDD3VeY2f1ArrvnAI8CT5pZHtERxeTYc1eY2bPASqAKuNPdq4OoW0REoiz65T2xZGdne25ubthliIjEFTNb6O7ZtW3TGdwiIlInhYWIiNRJYSEiInVSWIiISJ0ScoLbzIqB/BP4EV2AkgYqJx4kW39BfU4W6vOx6eXutZ7VnJBhcaLMLPdIRwQkomTrL6jPyUJ9bjjaDSUiInVSWIiISJ0UFrWbHnYBAUu2/oL6nCzU5waiOQsREamTRhYiIlInhYWIiNQpacPCzCaY2RozyzOze2rZ3tLMZsW2f2JmvYOvsmHVo8/fNbOVZrbMzGabWa8w6mxIdfW5RrurzMzNLO4Ps6xPn83smtjveoWZ/TXoGhtaPd7bWWb2jpktjr2/Lw6jzoZiZjPMbIeZLT/CdjOz38f+fywzs1En/KLunnT/iF4mfT3QF2gBLAUGH9bmDuBPsduTgVlh1x1An88F2sRufysZ+hxr1w6YC3wMZIdddwC/5/7AYqBj7H7XsOsOoM/TgW/Fbg8GNoVd9wn2eTwwClh+hO0XA68TXWl0HPDJib5mso4sxgJ57r7B3SuAmcCkw9pMAh6P3X4OOM/MalviNV7U2Wd3f8fdD8Tufkx0VcJ4Vp/fM8ADwINAeZDFNZL69Pk2YJq77wZw9x0B19jQ6tNnB9rHbqcR56ttuvtcouv+HMkk4AmP+hjoYGYZJ/KayRoWmcCWGvcLYo/V2sbdq4BSoHMg1TWO+vS5pilEv5nEszr7bGYjgZ7u/kqQhTWi+vyeBwADzOwDM/vYzCYEVl3jqE+fpwI3mFkB8Brw7WBKC82x/r3XKZCV8pqg2kYIhx9DXJ828aTe/TGzG4Bs4OxGrajxHbXPZpYC/Aa4JaiCAlCf33Mq0V1R5xAdPc4zs6HuvqeRa2ss9enzdcBj7v4rMzuD6KqcQ9090vjlhaLBP7+SdWRRAPSscb8H/zws/byNmaUSHboebdjX1NWnz5jZ+cAPgcvc/VBAtTWWuvrcDhgKvGtmm4ju282J80nu+r63X3L3SnffCKwhGh7xqj59ngI8C+DuHwGtiF5wL1HV6+/9WCRrWCwA+ptZHzNrQXQCO+ewNjnAzbHbVwFzPDZzFKfq7HNsl8yfiQZFvO/Hhjr67O6l7t7F3Xu7e2+i8zSXuXs8r8lbn/f2i0QPZsDMuhDdLbUh0CobVn36vBk4D8DMBhENi+JAqwxWDnBT7KiocUCpu287kR+YlLuh3L3KzO4C3iR6JMUMd19hZvcDue6eAzxKdKiaR3REMTm8ik9cPfv8S+Ak4G+xufzN7n5ZaEWfoHr2OaHUs89vAhea2UqgGrjb3XeGV/WJqWef/x142My+Q3R3zC3x/OXPzJ4huhuxS2we5l6gOYC7/4novMzFQB5wAPjGCb9mHP//EhGRgCTrbigRETkGCgsREamTwkJEROqksBARkTopLEREpE4KCxERqZPCQiRAZrbVzE4Luw6RY6XzLEQCEjtbehtwUgJcSkWSjEYWIgEws35ErwKaAuw0s52xa46JxAWNLEQCYmZ3AuPd/dqwaxE5VhpZiARnBLAk7CJEjofCQiQ4pxFd8lMk7mg3lEgAYgstlQH93D2ul/SU5KSRhUgwWsf+6W9O4pLeuCIBcPf9wJ+AlbH1B0TiinZDiYhInTSyEBGROiksRESkTgoLERGpk8JCRETqpLAQEZE6KSxERKROCgsREamTwkJEROr0/wE1C+qur5HZ1AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "l_star = σ_star(pc.t_grid)\n",
    "\n",
    "fig,ax = plt.subplots()\n",
    "\n",
    "ax.plot(pc.t_grid, l_star, label=r'$\\ell^*$')\n",
    "\n",
    "\n",
    "ax.set_ylabel(r'$\\ell$', fontsize=12)\n",
    "ax.set_xlabel(r'$t$', fontsize=12)\n",
    "\n",
    "\n",
    "ax.legend(fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now determine how many firms there are and range of tasks each will perform"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def t(σ_star):    \n",
    "    s=1.0\n",
    "    t = [s]\n",
    "    while σ_star(s)!=0:\n",
    "        s = s - σ_star(s)\n",
    "        t.append(s)\n",
    "    return np.array(t)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "t = t(σ_star)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAuIAAAK8CAYAAACqZ9q4AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzdd3TV9eH/8dc7e4cQVggz7L2HCxyouIoVUNwogrOOfqu1tnWgbdU6sEVUFAURcNUBioONCDJkEwiEkJCEEUJIQha5yX3//gD9UYqKIbmfO56Pczzn5t7Pvfd1e0zuq2/fw1hrBQAAAMCzgpwOAAAAAAQiijgAAADgAIo4AAAA4ACKOAAAAOAAijgAAADggBCnA9SFqKgom5CQIGOMQkNDlZiYqIMHD0qSEhMTf7zuZPed6FSuOZ3ra/v5tf06nnpdT7+HN7ynN7z3ibwpy4m8OdsPfCHjiXwx8/F8Pf/x/OmznMifP9sPAuEznkygfu4fHDx4UJmZmfnW2oY1eb5fFvF69erp9ttvV1hYmJo0aaLRo0dr6tSpkqTRo0f/eN3J7jvRqVxzOtfX9vNr+3U89bqefg9veE9veO8TeVOWE3lzth/4QsYT+WLm4/l6/uP502c5kT9/th8Ewmc8mUD93D+YOnWqbrnllqyaPp+pKQAAAIADPFbEjTFDjTFpxph0Y8zDJ3l8kDFmrTGmyhgz4oTHbjbG7Dj2z82eygwAAADUFY8UcWNMsKSXJV0iqbOka40xnU+4bLek0ZJmnvDc+pIekzRAUn9JjxljEuo6MwAAAFCXPDUi3l9SurU2w1pbKeldScOOv8Bam2mt3SjJfcJzL5Y0z1pbYK09JGmepKGeCA0AAADUFU8t1kyWlH3czzk6OsJd0+cmn3iRMWacpHGSlJSUVLOUAAAAqJHi4mLl5eXJ5XI5HaXWRUdHq1mzZgoKqt0xbE8VcXOS+2xtPtdaO1nSZElq2rTpqb42AAAATlNxcbH279+v5ORkRUZGypiT1Tff5Ha7lZubq/z8fDVq1KhWX9tTU1NyJDU/7udmkvZ44LkAAACoY3l5eUpOTlZUVJRflXBJCgoKUuPGjVVUVFT7r13rr3hyqyW1M8a0NsaESRolafYpPvcrSRcZYxKOLdK86Nh9AAAA8AIul0uRkZFOx6gzoaGhqqqqqvXX9UgRt9ZWSbpHRwv0VknvW2u3GGPGG2N+I0nGmH7GmBxJIyW9ZozZcuy5BZKe1NEyv1rS+GP3AQAAwEv420j48erqs3nsZE1r7VxJc0+479Hjbq/W0WknJ3vum5LerNOAAAAAgAdxsiYAAAAC1qFDh5SZmampU6fq0KFDHn1vj42IAwAAAN7m+++/1xdffKGsrCy5XC6NHTvWY+/NiDgAAAACmrXO7HxNEQcAAEBAePDBB9W5c2fl5eX9eF+fPn1077336vLLL9eIESM8moepKQAAAPB7aWlpWrJkiTp37vxfh/MkJCQoISFBo0eP9ngmRsQBAADg9z755BNde+21SktLU7t27ZyOI4kiDgAAgACwcOFCJSYmqkuXLgoNDXU6jiSmpgAAACAAbNy4Ue3atdNdd93ldJQfUcQBAABQq56Ys0Wpe4rr9D06N43TY1d0OaVrS0pKZK1VVlaWBg0aVKe5fg2mpgAAAMCvFRcXKz8/X0899ZTTUf4LI+IAAACoVac6Uu0pa9asUVJSkrp376558+ZpyJAhMsY4HYsRcQAAAPgvl8ul9957T2PGjNHgwYOVlpbmFSVcYkQcAAAAfiw0NFQzZsyQJD3++OPOhjkBI+IAAACAAyjiAAAAgAMo4gAAAIADKOIAAACAAyjiAAAAwK90sOSIcg6Vn9ZrsGsKAAAATpu11mu2Baxt1tofb2cdLNUb3+zS+2uydV5Q2Wm9LkUcAAAApyU0NFTl5eWKiopyOkqdcLlccrmN7p65Vl9s2quQoCD9tleyWh3ac1qvSxEHAADAaWnUqJFyc3OVnJysyMhIvxkZt9aquLxS6bt2683v9mnp7iMaN6iNbjmrlRrHRWjq1LWn9foUcQAAAJyWuLg4SdKePXvkcrkcTnP6rLUqd1XrcEWVKqvcyiisUs92zfT09S0VGxFaa+9DEQcAAMBpi4uL+7GQ+6rSI1V6b3W2pizbpdzCcrVrFKNxg1J068XJCgup/T1OKOIAAAAIaPklRzRteabeXpGlonKX+reqr/HDuui8Do0UFFR302wo4gAAAAhImfmlev2bDH3wfY5c1W5d1Lmxxg1qoz4tEzzy/hRxAAAABJT12YWavHSnvti8T6FBQRreJ1m3nZOiNg1jPJqDIg4AAAC/Z63V4u0H9OrinVq5q0BxESG6c3AbjT6rlRrFRjiSiSIOAAAAv1VZ5dacDXs0eWmG0vYfVlJ8hP5yWSeN6t9CMeHOVmGKOAAAAPxOyZEqvbtqt6Ys26W9RRXq0DhWL1zdQ5d3b1onO6DUBEUcAAAAfiPvcIWmfpupd77LUnFFlQam1Nffr+qmc9s39LqDhijiAAAA8HkZB0r0+jcZ+s/3uXK53bqkaxONG9RGPZvXczraT6KIAwAAwGet3X1Iry3Zqa9T9ys0OEgj+zbTbeekqHWDaKej/SKKOAAAAHyK2221KC1Pry3J0KrMAsVHhuqe89rqpjNaqWFsuNPxThlFHAAAAD6hssqtT9fnavLSDO3IK1FyvUg9enlnXdOvuaId3gGlJnwvMQAAAALK4QqXZq3arTeXZWpfcYU6NonVhGt66rLuSQoN9o4dUGqCIg4AAACvlFdcoTe/zdSM77J0+EiVzmyTqGdGdNegdg28bgeUmqCIAwAAwKvsPlimV5bs1H++z1GV261LuiXp9kEp6t7Me3dAqQmKOAAAALxCet5hTVq0U59u2KNgYzSybzONG5SilonevwNKTVDEAQAA4Kgte4o0adFOzd28VxEhwbrlzFYaOyhFjeMinI5WpyjiAAAAcMTa3Yf08sJ0LdiWp9jwEN11bhvdelZrJcb4zhaEp4MiDgAAAI+x1uq7jAJNXLRD36YfVEJUqP7vwva66cxWio8MdTqeR1HEAQAAUOestVq8/YBeXpiuNVmH1DA2XH++tJOuG9DCJ/cArw2B+akBAADgEW631dep+zVx0Q5tzi1Wcr1IPTmsi0b2ba6I0GCn4zmKIg4AAIBaV1Xt1ueb9urlRenavr9ErRKj9Ozw7rqyV7LCQnz3EJ7aRBEHAABAramscuvjdTl6ZfFOZR4sU/vGMXppVE9d1i1JIT58CmZdoIgDAADgtFW4qvX+mmy9unin9hRVqFtyvF67sY8u7NRYQUG+fwpmXaCIAwAAoMaq3VaTl+7U5KW7lF9yRH1bJujvV3XT4PYN/eIY+rpEEQcAAMCvVlTuUu6hcu0trtBn27fp7LYNdM/5vTSgdX0K+CmiiAMAAOCUHSw5oinLdmn6iiydZcuUEBWmj285U71aJDgdzedQxAEAAPCL9hdXaPLSDM1cuVsVVdW6tFuSupTUU1RYMCW8hijiAAAA+EnZBWV6dclOfbAmR9XWaljPprrr3LZq2yhGU6dudDqeT6OIAwAA4H/sPFCiSYt26pP1uQo2RiP6NtOdg9uoef0op6P5DYo4AAAAfrR1b7FeXpSuzzftVXhIkG4+o5XGDUpRk/gIp6P5HYo4AAAAtD67UBMXpmv+1v2KCQ/RnYPb6NazW6tBTLjT0fwWRRwAACCArcw4qImL0vXNjnzViwrVA0Paa/SZrRQfFep0NL9HEQcAAAgw1lot3ZGvlxema1VmgRrEhOlPl3TU9QNbKiaceugp/C8NAAAQINxuq/lb92vionRtzClSUnyEHr+is0b1b6GI0GCn4wUcijgAAICfq3Zbfb5pr15emK60/YfVon6Unr6qm67q3UxhIUFOxwtYFHEAAAA/5ap26+N1uXpl8U7tyi9V20YxmnBNT13ePUkhwRRwp1HEAQAA/EyFq1offJ+jVxfvVG5hubo0jdMr1/fWxV2aKCjIOB0Px1DEAQAA/ERZZZVmrtytyUszlHf4iHq3qKenruyqczs0lDEUcG9DEQcAAPBxxRUuvb08U1OW7dKhMpfObJOoCaN66oyURAq4F6OIAwAA+KiC0kq99e0uTV2eqcMVVTq/YyPdfV5b9WmZ4HQ0nAKKOAAAgI/JK67Q699kaMbK3Sp3VWtolya6+7y26poc73Q0/AoUcQAAAB+RW1iu15bs1Lurs1VV7dawnsm669w2atc41uloqAGKOAAAgJfblV+qVxan66O1uTJGGtGnme4Y3EYtE6OdjobTQBEHAADwUjsPlOhfC3ZozoY9Cg0O0g0DW2rcoBQ1rRfpdDTUAoo4AACAl8k6WKqXFuzQJ+tyFR4SrLHnpOi2c1LUMDbc6WioRRRxAAAAL5FzqEz/XpCuD9fmKCTIaMzZrXX74DZqEEMB90cUcQAAAIftLSrXxIXpen9NtowxunFgS911bhs1iotwOhrqEEUcAADAIXnFFZq0eKdmrtwtK6tr+jXX3ee1VVI8c8ADAUUcAADAw/JLjujVxTs1/bssVbmtRvZppnvOb6tmCVFOR4MHUcQBAAA85FBppV5bmqFpyzN1pKpav+3VTPde0JZtCAMURRwAAKCOFZW59MayDL25bJfKXNX6TY+muu+CdkppGON0NDiIIg4AAFBHDle49OayTL2xLEOHK6p0Wbck3TekndpzEiZEEQcAAKh1pUeqNHV5pl7/JkOFZS5d1Lmx7h/SXp2bxjkdDV6EIg4AAFBLyiurNf27TL26JEMFpZU6v2MjPTCkvbo1i3c6GrwQRRwAAOA0ua1VXvERnfPsIuWXHNE57Rro9xe2V68WCU5HgxejiAMAANTQkapqvb86W+uzC1VZ5Va7ZjF65Ybe6teqvtPR4AMo4gAAAL+Sq9qtD7/P0cSF6cotLNfIuGC1bRije8cNdDoafAhFHAAA4BRVVbv18bpc/WvhDmUXlKtn83r6x1XdlLH8gNPR4IMo4gAAAL+g2m01Z8MevbRgh3bll6prcpyeGN1F53VoJGOMMpY7nRC+iCIOAADwE9xuqy8279OE+du1I69EHZvE6rUb++iizo1ljHE6HnwcRRwAAOAE1lp9nbpfL87brm37Dqtdoxi9fF1vXdK1iYKCKOCoHRRxAACAY6y1WpSWpxfmbdfm3GKlNIjWS6N66vLuTRVMAUcto4gDAICAZ63V0h35emHedm3ILlSL+lF6bmQPXdmzqUKCg5yOBz9FEQcAAAFtefrRAr4m65CS60Xq6au6aXifZgqlgKOOUcQBAEBAWp1ZoOe/TtN3GQVqEhehJ6/sqmv6NldYCAUcnkERBwAAAWXt7kN6cd52fbMjXw1iwvXYFZ11bf8WiggNdjoaAgxFHAAABIRNOUV6YV6aFqUdUP3oMP350k66YWBLRYZRwOEMijgAAPBrqXuK9eL87ZqXul/1okL10NAOuvmMVooOpwbBWfwbCAAA/NL2/Yc1Yf52zd20T7ERIfr9he11y1mtFBsR6nQ0QBJFHAAA+JmdB0r00vwdmrNxj6LDQnTv+W015pwUxUdSwOFdKOIAAMAvZB0s1UsLduiTdbkKDwnWHYPbaNw5KUqIDnM6GnBSFHEAAODTcg6V6d8L0vXh2hyFBBmNObu1bh/cRg1iwp2OBvwsijgAAPBJe4vKNXFhut5fky0joxsHttRd57ZRo7gIp6MBp4QiDgAAfEp+yRFNXJiumSt3y8rqmn7Ndfd5bZUUH+l0NOBXoYgDAACfUHKkSm98k6HXl2aoosqtkX2a6e7z2qp5/SinowE1QhEHAABezVXt1rurduulBTuUX1KpS7s10R8u6qCUhjFORwNOC0UcAAB4Jbfbau7mvXruqzRlHizTgNb19fpNHdWrRYLT0YBaQREHAABeZ3l6vp7+cps25hSpY5NYvTW6n87t0FDGGKejAbWGIg4AALxG6p5iPf3lNi3dfkDJ9SL1/MgeurJXsoKDKODwPxRxAADguOyCMr0wb7s+WZ+ruIhQ/fnSTrrxjJaKCA12OhpQZyjiAADAMQWllZq4MF3vfJclY6Q7BrfRHYPbcBw9AgJFHAAAeFxZZZXe+jZTry7eqdLKKl3dt7nuH9JeTeI5jAeBgyIOAAA8pqrarffX5GjC/O3KO3xEF3ZurD8O7aC2jWKdjgZ4HEUcAAB4xJeb9+rZr9KUcaBUfVomaNL1vdW3VX2nYwGOoYgDAIA6dbiiSrsLSvX4O2vVtlGMXr+pr4Z0asRWhAh4FHEAAFAn0vYd1rNfblPoniKFhQTp2eHddVXvZIUEBzkdDfAKFHEAAFCrcgvL9eK87frP2hzFhIfo2vpRahIfoav7NXc6GuBVKOIAAKBWFJZV6pXFO/XW8kzJSred3Vp3n9dWn7y/1+logFeiiAMAgNNS4arW1OWZmrQoXYePVOmqXs30wIXt1CwhyulogFejiAMAgBqpdlv9Z22OXpy3XXuLKnR+x0Z6aGgHdWwS53Q0wCdQxAEAwK9irdWCrXl69qtt2r6/RD2b19OL1/TUwJREp6MBPsVjRdwYM1TSS5KCJb1hrX36hMfDJb0tqY+kg5KusdZmGmNCJb0hqfexvG9ba//hqdwAAOD/+z6rQE9/sU2rMw8ppUG0Xrm+t4Z2bcJWhEANeKSIG2OCJb0s6UJJOZJWG2NmW2tTj7tsjKRD1tq2xphRkp6RdI2kkZLCrbXdjDFRklKNMbOstZmeyA4AAKT0vBI9++U2fZ26Xw1jw/W333bV1X2bK5StCIEa89SIeH9J6dbaDEkyxrwraZik44v4MEmPH7v9oaSJ5uj/vbaSoo0xIZIiJVVKKvZQbgAAAtr+4gpNmL9d763OVlRYiP5wUXvdenZrRYUxuxU4XZ76LUqWlH3czzmSBvzUNdbaKmNMkaREHS3lwyTtlRQl6QFrbUGdJwYAIIAVV7j02pKdmrJsl6rdVjef2Ur3nNdWiTHhTkcD/IanivjJJo7ZU7ymv6RqSU0lJUj6xhgz/4fR9R+fbMw4SeMkKSkp6bQDAwAQiI5UVWv6iixNXJSuwjKXruzZVP93UQc1r89WhEBt81QRz5F0/HFazSTt+Ylrco5NQ4mXVCDpOklfWmtdkvKMMd9K6ivpv4q4tXaypMmS1LRp0xNLPgAA+BnVbqtP1+fq+a+3K7ewXOe0a6A/Du2orsnxTkcD/JanivhqSe2MMa0l5UoapaMF+3izJd0saYWkEZIWWmutMWa3pPONMe/o6NSUgZImeCg3AAB+zVqrJdsP6OkvtmnbvsPqmhynZ4Z319ntGjgdDfB7Hinix+Z83yPpKx3dvvBNa+0WY8x4SWustbMlTZE03RiTrqMj4aOOPf1lSW9J2qyj01festZu9ERuAAD82YbsQj39xTatyDioFvWj9O9re+mybkkKCmIrQsATPLbk2Vo7V9LcE+579LjbFTq6VeGJzys52f0AAKBmduWX6rmv0vT5pr1KjA7TE7/pomv7t1BYCFsRAp7E3kMAAASIvMMV+veCdM1atVthIUG674J2GjsoRTHh1AHACfzmAQDg50qOVGny0gy98U2GKqvcurZ/C917QTs1jGUrQsBJFHEAAPxUZZVbM1dm6d8L03WwtFKXdU/SHy7qoNYNop2OBkAUcQAA/I61VnM27tVzX6Vpd0GZzkhJ1MOXdFSP5vWcjgbgOBRxAAD8yPrsQj0xZ4vW7S5Up6Q4Tbu1vwa1ayBj2AkF8DYUcQAA/MC+ogo9++U2fbQuVw1jw/XsiO4a0bsZWxECXowiDgCAD6twVev1pRmatHinqq3VXee20V3ntWUnFMAH8FsKAIAPstbq80179Y+525RbWK5LujbRI5d2UvP6UU5HA3CKKOIAAPiYTTlFGv/ZFq3OPKROSXF6bmQPndEm0elYAH4lijgAAD4i73CF/vllmj5cm6P6UWH6x1XddHXf5gpmHjjgkyjiAAB4uQpXtaYs26VJi9JVWe3WuHNSdPf5bRUXEep0NACngSIOAICXstbqy8379Pcvtiq7oFwXdm6sP1/aSa04kAfwCxRxAAC80JY9RRo/J1UrdxWoQ+NYzbhtgM5q28DpWABqEUUcAAAvcuDwEb0wL03vrs5WvchQPXVlV43q11whwUFORwNQyyjiAAB4Abe12ldUofOeW6wKV7VuPau17r2gneIjmQcO+CuKOAAADrLWal7qfm3MKVKFq1r9U+rrz5d1UpuGMU5HA1DHKOIAADhk275iPflZqr5NP6jfxkgdk+J0/+h+TscC4CEUcQAAPOxgyRG9MG+7Zq3ardiIUD3xmy5ypx2QYTtwIKBQxAEA8JDKKrfeXpGplxbsUFlltW46o5XuH9JO9aLCNHW70+kAeBpFHACAOmat1cJtefrb51uVkV+qQe0b6q+XdVK7xrFORwPgIIo4AAB1aMf+wxr/Waq+2ZGvlAbRenN0X53XoZEM81CAgEcRBwCgDhwqrdSE+dv1zsrdig4L1l8v76wbB7ZUWAj7gQM4iiIOAEAtclW7NeO7LL04f4cOV7h03YAW+v2FHVQ/OszpaAC8DEUcAIBasjgtT099vlXpeSU6q22i/np5Z3VsEud0LABeiiIOAMBp2nmgRE99lqpFaQfUKjFKr9/UV0M6MQ8cwM+jiAMAUENFZS69tGCH3l6RqcjQYD1yaUfdfGYrhYcEOx0NgA+giAMA8CtVVbs1a9VuvTBvuwrLXRrVr7n+76IOahAT7nQ0AD6EIg4AwK+wbEe+nvwsVWn7D2tgSn399fLO6tI03ulYAHwQRRwAgFOwK79Uf/s8VfO35ql5/Ui9ekNvXdylCfPAAdQYRRwAgJ9RXOHSvxfs0NTlmQoLDtJDQzvo1rNaKyKUeeAATg9FHACAk6h2W723OlvPf52mgrJKjezTTH+4qIMaxUU4HQ2An6CIAwBwguU78zV+Tqq27Tusfq0SNPXy/urWjHngAGoXRRwAgGN2HyzT3+am6qst+5VcL1ITr+uly7olMQ8cQJ2giAMAAt7hCpdeXrRTby7bpZBgoz9c1F63nZPCPHAAdYoiDgAIWNZazdm4V099lqq8w0d0Ve9k/XFoRzVmHjgAD6CIAwACUsaBEj366RYtS89X1+Q4vXZjH/VqkeB0LAABhCIOAAgoFa5qTVqUrleXZCg8JEhP/KaLbhjYUsFBzAMH4FkUcQBAwFiUlqfHPt2i3QVlGtazqf58WSc1imUaCgBnUMQBAH5vb1G5xs9J1Reb9ymlYbRm3jZAZ7Zt4HQsAAGOIg4A8FuuaremfpupF+dvV7Xb6g8XtdfYQSkKD2E3FADOo4gDAPzSmswC/eWTzdq277DO79hIT/ymi5rXj3I6FgD8iCIOAPArBaWVevqLrXp/TY6axkfotRv76KLOjTmUB4DXoYgDAPyC2231/ppsPf3lNpVUVOn2QSm694J2ig7nqw6Ad+KvEwDA56XuKdZfPtmktbsL1b9VfT15ZVd1aBLrdCwA+FkUcQCAz6p2W+UcKteTE5cpPjJUz43soeG9k5mGAsAnUMQBAD7HWqu5m/ZpQ06hXNVuXdOvuR66uIPqRYU5HQ0AThlFHADgUzLzS/Xo7C1auv2AhscEqX3jWN3z225OxwKAX40iDgDwCRWuar26ZKcmLd6psOAgPXZFZ2n7ATELBYCvoogDALze0u0H9Oinm5V5sEyXd0/SXy/vrMZxEZq6Y4nT0QCgxijiAACvta+oQk9+nqrPN+5V6wbRmj6mv85p19DpWABQKyjiAACvU1Xt1rQVWXpx3nZVVrv1+wvba9ygFEWEcjQ9AP9BEQcAeJXvsw7pL59s1ta9xRrcvqHGD+uilonRTscCgFpHEQcAeIVDpZV69qttmrUqW03iIvTK9b01tGsT9gQH4Lco4gAAR7ndVh+uzdHTX2xTUblLY89prfuGtFcMR9MD8HP8lQMAOGbbvmL99ZPNWp15SH1aJuipK7uqU1Kc07EAwCMo4gAAjys9UqWXFuzQlGW7FBcRomeHd9eIPs0UFMQ0FACBgyIOAPAYa62+2rJPT8xJ1d6iCo3q11x/HNpRCdEcTQ8g8FDEAQAesftgmR6bvVmL0g6oY5NYTbyul/q0rO90LABwDEUcAFCnjlRVa/KSDE1clK6QIKO/XNZJo89spZDgIKejAYCjKOIAgDrzbXq+/vrJZmXkl+rSbk3018s7Kyk+0ulYAOAVKOIAgFqXV1yhpz7fqtkb9qhlYpSm3tJP53Zo5HQsAPAqFHEAQK2pdltNX5Gp57/eriNVbt17QTvddW4bjqYHgJOgiAMAasX67EL9+eNN2rKnWOe0a6Dxw7qqdQOOpgeAn0IRBwCclqIyl579aptmrtqthjHhmnhdL13WLYmj6QHgF1DEAQA1Yq3V55v26vHZW1RQWqlbzmytBy5sp9iIUKejAYBPoIgDAH61vMMV+usnm/XVlv3q3ixe027try5N452OBQA+hSIOADhl1lp9tDZX4z9LVbmrWg9f0lG3nd2aPcEBoAYo4gCAU5JbWK5HPtqkJdsPqG/LBD0zorvaNIxxOhYA+CyKOADgZ7ndVrNW79Y/5m5Ttdvq8Ss666YzWikoiMWYAHA6KOIAgJ+UdbBUD/9nk1ZkHNSZbRL19FXd1SIxyulYAOAXKOIAgP9R7baatjxT//wqTSFBRv+4qptG9WvOloQAUIso4gCA/5KeV6KHPtygtbsLdV6Hhvr7Vd2UFB/pdCwA8DsUcQCAJKmq2q3J32RowvwdigoL1ovX9NCVPZMZBQeAOkIRBwBo695iPfThRm3KLdKl3Zroid90VcPYcKdjAYBfo4gDQACz9ui2hE/9e5nqRYVq0vW9dWm3JKdjAUBAoIgDQIDakF2oTblFKqus0hU9murRyzsrITrM6VgAEDAo4gAQYCpc1Xpx3na9/k2GLo90q2OTWN1/TU+nYwFAwKGIA0AAWZ1ZoIc+3Khd+aW6tn9ztS3MUzAH8wCAI4KcDgAAqHulR6r0+Owtuvq1FXJVuzXjtgH6x1XdKeEA4CBGxAHAzy3bka+HP9qo3MJy3XxGKz14cQdFh/PnHwCcxl9iAPBTxTso83YAACAASURBVBUu/f3zrXp3dbZSGkTr/dvPUL9W9Z2OBQA4hiIOAH5owdb9+vPHm5V3uEK3D07RA0PaKyI02OlYAIDjUMQBwI8cKq3U+M9S9fG6XHVoHKvXbuyjHs3rOR0LAHASFHEA8BNzN+3Vo59uVmGZS/dd0E53n9dWYSGsyQcAb0URBwAfl3e4Qo99ukVfbN6nrslxmj5mgDolxTkdCwDwCyjiAOCjrLX6ZH2unpiTqrLKaj00tIPGnZOikGBGwQHAF1DEAcAH7S0q1yMfbdKitAPq0zJBzwzvrraNYpyOBQD4FSjiAOBDrLV6d3W2/v75VlW5rR67orNuOqMVB/MAgA+iiAOAj8guKNPDH23Ut+kHdUZKop4Z3l0tEqOcjgUAqCGKOAB4Obfb6u0VmXrmyzQFBxn9/bfddG3/5jKGUXAA8GUUcQDwYjsPlOiPH27UmqxDOrdDQ/39t93UtF6k07EAALWAIg4AXqiq2q3Xv9mlF+dvV2RosF64uod+2yuZUXAA8CMUcQDwMtv2FevBDzZqU26RLu7SWE9e2VWNYiOcjgUAqGUUcQDwEpVVbk1anK6XF6UrLiJUL1/XW5d2a8IoOAD4KYo4AHiB1D3F+v3767Vt32Fd2bOpHr2ii+pHhzkdCwBQhyjiAOAgt9tqyrJd+udXaYqPCtUbN/XVkM6NnY4FAPAAijgAOGRvUbn+7/0NWr7zoC7q3FhPD+/OKDgABBCKOAA44PONe/XIx5vkqnbrmeHddHVf9gUHgEBDEQcADzpc4dLjs1P1n7U56tG8niZc01OtG0Q7HQsA4ACKOAB4yJrMAj3w/nrlHirXvRe00+/Ob6vQ4CCnYwEAHEIRB4A65qp2618LdujlRelKTojUB3ecoT4t6zsdCwDgMIo4ANShXfmluv+99dqQXagRfZrpsSs6KzYi1OlYAAAvQBEHgDpgrdW7q7M1fk6qwkKCNOn63rq0W5LTsQAAXoQiDgC1rKC0Ug//Z6O+Tt2vs9om6rmRPZQUH+l0LACAl6GIA0AtKixz6eIJS1VU5tJfLuukW89qraAgtiUEAPwvijgA1IIKV7UyD5ZqX1GFEhJCNe2W/urcNM7pWAAAL0YRB4DTtGVPke5/d73aFFUoKT5CE+45WxGhwU7HAgB4OYo4ANSQ2231xrIM/fOrNNWLClOnpDjFR4ZSwgEAp4STJACgBvYUluv6N1bq73O36fyOjfTV/YMUH8m2hACAU8eIOAD8Sp9t3KNHPtqkKrfVM8O76eq+zWUMCzIBAL8ORRwATtHhCpce+3SLPlqXq57N62nCNT3VqkG007EAAD6KIg4Ap2B1ZoEeeG+99hSW694L2ul357dVaDCz+wAANUcRB4Cf4ap266X5OzRpcbqaJUTpgzvOVJ+WCU7HAgD4AYo4APyEjAMleuC99dqQU6SRfZrpsd90UUw4fzYBALWDbxQAOIG1VrNWZevJz1IVFhKkSdf31qXdkpyOBQDwMx6b4GiMGWqMSTPGpBtjHj7J4+HGmPeOPb7SGNPquMe6G2NWGGO2GGM2GWMiPJUbQGA5WHJEY9/+Xo98vEm9W9bTV/cPooQDAOqER0bEjTHBkl6WdKGkHEmrjTGzrbWpx102RtIha21bY8woSc9IusYYEyLpHUk3Wms3GGMSJbk8kRtAYFmUlqcHP9io4nKX/nJZJ916VmsFBbEtIQCgbnhqakp/SenW2gxJMsa8K2mYpOOL+DBJjx+7/aGkieboxrwXSdpord0gSdbagx7KDCBAVLiq9Y+5WzVtRZY6NI7V9DH91SkpzulYAAA/56kiniwp+7ifcyQN+KlrrLVVxpgiSYmS2kuyxpivJDWU9K619tkT38AYM07SOElKSuI/IwM4NVv2FOm+d9crPa9Et57VWg8N7cAR9QAAj/BUET/Zf9u1p3hNiKSzJfWTVCZpgTHme2vtgv+60NrJkiZLUtOmTU98bQD4L2631evfZOi5r9OUEBWmt2/tr0HtGzodCwAQQDxVxHMkNT/u52aS9vzENTnH5oXHSyo4dv8Sa22+JBlj5krqLWmBAKAG9hSW6/fvr9d3GQW6uEtj/eOq7qofHeZ0LABAgPHUrimrJbUzxrQ2xoRJGiVp9gnXzJZ087HbIyQttNZaSV9J6m6MiTpW0Afrv+eWA8Apm71hj4ZOWKqNOUV6dnh3vXpDH0o4AMARHhkRPzbn+x4dLdXBkt601m4xxoyXtMZaO1vSFEnTjTHpOjoSPurYcw8ZY17Q0TJvJc211n7uidwA/EdxhUuPfbpFH6/LVa8W9TThmp5qmRjtdCwAQADz2IE+1tq5kuaecN+jx92ukDTyJ577jo5uYQgAv9qqXQV64L312ltUrvsuaKffnd9WIcEeO0YBAICT4mRNAH6rqtqtCfN3aNLidDVLiNIHd5ypPi0TnI4FAIAkijgAP5VXXKF7Zq3Tql0FGtGnmR7/TRfFhPMnDwDgPfhWAuB3Vuw8qN/NWqeSIy49P7KHhvdp5nQkAAD+B0UcgN9wu61eWbJTz3+dplYNojXjtgHq0CTW6VgAAJwURRyAXygsq9Tv39+ghdvydHn3JD09vDtTUQAAXo1vKQA+b0N2oe6asVZ5hys0flgX3TiwpYw52WG9AAB4D4o4AJ9lrdU732Xpyc+2qmFsuD6440z1bF7P6VgAAJwSijgAn1R6pEoPf7RJczbs0XkdGuqFq3sqgRMyAQA+hCIOwOfs2H9Yd7zzvXbll+rBizvozsFtFBTEVBQAgG+hiAPwKZ+sy9WfPtqk6PBgvTNmgM5s28DpSAAA1AhFHIBPcFurrINlmvbeevVvVV//vq6XGsdFOB0LAIAao4gD8HrZBWXasqdYpUeqdPvgFD14UQeFBAc5HQsAgNPCNxkArzYvdb8u+9c3OuKqVocmsfrTJZ0o4QAAv8CIOACvVFXt1j+/TtNrSzLUNTlO3SLqKTyUAg4A8B98qwHwOnnFFbrujZV6bUmGrhvQQh/ecSYlHADgdxgRB+BVlu/M172z1qv0SJVevKaHfturmdORAACoExRxAF7B7bZ6ZclOPf91mlo3iNbMsQPUvnGs07EAAKgzFHEAjissq9QD763XorQDuqJHUz19VTdFh/PnCQDg3/imA+CoDdmFumvGWuUdrtCTw7rohoEtZQynZAIA/B9FHIAjrLWa/l2WnvwsVY1iI/ThHWeqR/N6TscCAMBjKOIAPK7kSJX+9NEmzdmwR+d3bKQXru6helFhTscCAMCjKOIAPGr7/sO6853vtSu/VA9e3EF3Dm6joCCmogAAAg9FHIDHfLwuR498tFnR4SF657YBOrNNA6cjAQDgGIo4gDpX4arW+M9SNXPlbvVvXV8Tr+2lRnERTscCAMBRFHEAdWr3wTLdNfN7bc4t1h2D2+gPF7VXSDCnZAIAQBEHUGfmpe7X799fLyPpjZv6akjnxk5HAgDAa1DEAdS6qmq3/vl1ml5bkqFuyfGadH1vNa8f5XQsAAC8CkUcQK3aX1yh381cp1WZBbphYAv95bLOiggNdjoWAABehyIOoNYs35mve2etU+mRar00qqeG9Ux2OhIAAF6LIg7gtLndVq8s2annv05TSsMYzRrbW+0axzodCwAAr0YRB3BaCssq9cB767Uo7YCG9Wyqv/+2m6LD+dMCAMAv4dsSQI2l5x3WmGlrtLewQk9d2VXXD2ghYzglEwCAU0ERB1AjS7Yf0D0z1io8NEizxg1Un5YJTkcCAMCnUMQB/CrWWk1dnqknP0tVhyZxeuPmvkquF+l0LAAAfA5FHMApc1W79einWzRr1W5d1LmxXrymJ/PBAQCoIb5BAZySQ6WVumvGWq3IOKi7zm2jP1zUQUFBzAcHAKCmKOIAflF6XonGTFutvYUVeuHqHrqqdzOnIwEA4PMo4gB+1tLtB3T3zLUKD2FRJgAAtYkiDuCkrLWatjxT4z9LVfvGsXrj5r5qlhDldCwAAPwGRRzA/3BVu/XY7C2auXK3LuzcWBNYlAkAQK3jmxXAf6lyW9385iot33lQd57bRg+yKBMAgDpBEQfwo3JXtdL2HdaaikMsygQAoI5RxAFIOrooc0tukYwxmjVugPq0rO90JAAA/FqQ0wEAOOuHRZm3TF2tsJBgdU2Op4QDAOABjIgDAcxV7dbjs7doxsrdGtKpsbpU71cw88EBAPAIRsSBAFVYVqmb31ylGSt3647BbTT5xj6UcAAAPIgRcSAA7TxQojFTV2tPYYWeH9lDw/uwKBMAAE+jiAMB5psdB3TXjLUKCw5iUSYAAA6iiAMB5O0VmXpiTqraNYrhpEwAABxGEQcCgKvarSfmbNE73x1dlDlhVE/FcFImAACO4psY8HOFZZW6e+ZafZt+ULcPTtFDF3dkUSYAAF6AIg74sZ0HSnTbtDXKPVSu50b20AgWZQIA4DUo4oCfOn5R5syxA9S3FYsyAQDwJhRxwA8dvyjz9Zv6qnl9FmUCAOBtKOKAH/nvRZmNNGFULxZlAgDgpfiGBvwEizIBAPAtFHHAD/ywKDPnUJn+OaK7RvZt7nQkAADwCyjigI/7ZscB3T1jrUKCgzRz7ED1Y1EmAAA+gSIO+LAfFmW2bXj0pEwWZQIA4Dso4oAPclW7NX5OqqZ/l6ULOjbSS9eyKBMAAF/DNzfgY4rKXLpr5vdHF2UOStFDQ1mUCQCAL6KIAz5kT2G5bnpzlbIOlrIoEwAAH0cRB3xEel6JbpqyUocrqjR9zAANTEl0OhIAADgNFHHAB6zPLtQtb61ScFCQ3r19oLo0jXc6EgAAOE0UccDLfbPjgG6f/r0SY8L0zpgBapkY7XQkAABQCyjigBf7bOMePfDeerVpGKO3b+2vRnERTkcCAAC1hCIOeKnpKzL16Owt6teyvl6/ua/iI0OdjgQAAGrRzxZxY0wDSQ0k7bDWVnsmEhDYrLV6acEOTZi/Q0M6NdLE63orIjTY6VgAAKCW/dKI+BJJWyR1MMZskTTJWrtMkowxQyS1t9ZOquOMQMBwu62emLNF01ZkaXjvZnpmeDeFBAc5HQsAANSBX/qGnyuppaQ3JX0h6WpjzKfGmK7W2vmSrqnrgECgqKxy67731mvaiiyNPae1/jmiOyUcAAA/9rMj4tbaB40x7SVdKel6Sd0lGUnrjDF7JE2p+4iA/6t2W42Ztlrf7MjXw5d01B2D2zgdCQAA1LFfXKxprd0u6dlj/8gYEy0pTlI/SSPrNB0QAKqqrbbtK9a3Jfl6dnh3Xd2P0zIBAAgEv3rXFGttqaRSSbONMXNrPxIQOPYUlmvL3iIdcbn1yg19dHGXJk5HAgAAHnJaE1CttVW1FQQINOl5JRrxynJVVrnVsUksJRwAgADDPuKAAzZkF2r0W6sUHGTUOSlO0eH8KgIAEGjYkgHwsGU78nXt698pJiJEH95xJiUcAIAARQMAPOjzjXt1/3vrOLIeAABQxAFPmf5dlh79dLP6tEjQlJv7KT6KI+sBAAhkFHGgjllr9a8F6Xpx/nad37GRXr6utyLDOLIeAIBARxEH6tDxR9Zf1TtZzwzvrlBOywQAAKKIA3WmssqtP3ywQbM37NFtZ7fWI5d2UlCQcToWAADwEhRxoA6UVVbpjnfWaun2A/rj0I66Y3CKjKGEAwCA/48iDtSywrJK3TJ1tTZkF+qZ4d10Tb8WTkcCAABeiCIO1KK9ReW6acoqZRWUadL1fTS0K6dlAgCAk6OIA7Vk54ES3TRllYrKXZp2S3+d0SbR6UgAAMCLUcSBWrAxp1Cj31otI+ndcQPVNTne6UgAAMDLUcSB07RsR75un75GCdFhmj5mgFo3iHY6EgAA8AEUceA0zN20V/e/u16tG0Tr7TH91Zgj6wEAwCmiiAM1NGNllv7yCUfWAwCAmqGIA7+StVYTF6br+XkcWQ8AAGqOIg78Cm631fjPUjV1eaau6pWsZ0ZwZD0AAKgZijhwiqrdVg99uFH/WZujMWe31p85sh4AAJwGijhwCtxuqz99dLSE33dBO90/pB1H1gMAgNNCEQd+gdtt9edPNuv9NTm69/y2euDC9k5HAgAAfoDJrcDPsNbqsdlbNGvVbt15bhtKOAAAqDUUceAnWHt0Yeb077I09pzWeujiDkxHAQAAtYYiDpyEtVb/+GKb3vo2U7ec1UqPXNqJEg4AAGoVRRw4gbVWz36VpslLM3TTGS316OWdKeEAAKDWUcSBE7w4b7teWbxT1/Zvocev6EIJBwAAdYIiDhznXwt26F8L03V132b625Vd2SccAADUGbYvBI7JLSzX69u266reyXr6qu6UcAAAUKcYEQck7S0qV3ZBmYb1bKp/juhBCQcAAHWOIo6AN2XZLmUdLFNiTLieH9lDwZRwAADgARRxBLS3V2Tqyc9SlRgdprYNYxQSzK8EAADwDOaII2DNWJmlRz/doiGdGqut+4DYHAUAAHgSw38ISO+t3q0/f7xZ53dspJev70UJBwAAHkcRR8D58PscPfzRJg1q31CTru+t8JBgpyMBAIAARBFHQPlkXa4e/HCDzmrTQJNv7KOIUEo4AABwBkUcAWPOhj36/fvrNaB1fb1+U19KOAAAcBRFHAHhi017df9769W3ZX29ObqfIsMo4QAAwFkUcfi9r7fs0+9mrVPP5vX05i39FBXGZkEAAMB5FHH4tQVb9+vumWvVNTleU2/pp5hwSjgAAPAOFHH4rcVpebrznbXqlBSnabf2V2xEqNORAAAAfkQRh19atiNf46Z/r7aNYvT2rf0VH0kJBwAA3sVjRdwYM9QYk2aMSTfGPHySx8ONMe8de3ylMabVCY+3MMaUGGP+4KnM8E3Ld+ZrzLTVSmkQrRm3DVC9qDCnIwEAAPwPjxRxY0ywpJclXSKps6RrjTGdT7hsjKRD1tq2kl6U9MwJj78o6Yu6zgrftjLjoMZMXaOWiVGacdsAJURTwgEAgHfy1Ih4f0np1toMa22lpHclDTvhmmGSph27/aGkC4w5evC4MeZKSRmStngoL3zQmswC3TJ1tZrWi9CM2wYqMSbc6UgAAAA/yVNFPFlS9nE/5xy776TXWGurJBVJSjTGREv6o6Qnfu4NjDHjjDFrjDFrai01fMba3Yc0+q3VahIXoVljB6phLCUcAAB4N08VcXOS++wpXvOEpBettSU/9wbW2snW2r7W2r41zAgftTGnUDdPWaXEmDDNHDtQjeIinI4EAADwizy1qXKOpObH/dxM0p6fuCbHGBMiKV5SgaQBkkYYY56VVE+S2xhTYa2dWPex4e025xbphjdWKj4qVDPHDlSTeEo4AADwDZ4q4qsltTPGtJaUK2mUpOtOuGa2pJslrZA0QtJCa62VdM4PFxhjHpdUQgmHJKXuKdYNU1YqNiJUs8YOVHK9SKcjAQAAnDKPFHFrbZUx5h5JX0kKlvSmtXaLMWa8pDXW2tmSpkiaboxJ19GR8FGeyAbflJ5XohumrFRkaLBmjR2o5vWjnI4EAADwq3jsvG9r7VxJc0+479HjbldIGvkLr/F4nYSDTzlw+IhGv7VKQUaaOXagWiRSwgEAgO/xWBEHakNZZZXGTFutgyWVenfcQLVuEO10JAAAgBrhiHv4jGq31b2z1mlzbpH+dW0v9Whez+lIAAAANcaIOHyCtVZPzNmi+VvzNH5YF13YubHTkQAAAE4LI+LwCW98s0tvr8jSuEEpuumMVk7HAQAAOG0UcXi9zzfu1d/mbtVl3ZL08NCOTscBAACoFRRxeLU1mQV64P316tsyQc9f3UNBQSc7gBUAAMD3UMThtTIOlGjs22uUXC9Sr9/UVxGhwU5HAgAAqDUUcXglV7Vbo99aLWOMpt7STwnRYU5HAgAAqFXsmgKv47ZW2/cf1v6yCs0aN1AtE9krHAAA+B9GxOFVqt1W6XklKjlSpZdG9VLvFglORwIAAKgTFHF4lb99vlUFpZVqWT9aQ7s2cToOAABAnaGIw2u8uWyX3vx2l5LiI9QkPsLpOAAAAHWKOeLwCl9u3qcnP0/VxV0aq0X1AafjAAAA1DlGxOG4tbsP6b5316lHs3qacE0vsVM4AAAIBBRxOCrrYKlum7ZGTeIjNOXmvooMY69wAAAQGCjicExBaaVGv7Vabmv11uh+SowJdzoSAACAxzBHHI6ocFVr3NtrlFtYrpm3DVBKwxinIwEAAHgUI+LwOLfb6v/e36A1WYf04tU91bdVfacjAQAAeBxFHB73zJfb9PmmvXrk0o66rHuS03EAAAAcQRGHR01fkanXlmbopjNaauw5KU7HAQAAcAxFHB4zP3W/Hpu9RUM6NdJjV3SRMWxUCAAAAhdFHB6xMadQv5u1Tl2T4/Wva3spOIgSDgAAAhtFHHUuu6BMt05drcSYME25uZ+iwtisBwAAgEaEOlVU5tLot1bJVW317rh+ahjLXuEAAAASI+KoQ0eqqjVu+hplF5Rr8o191LZRrNORAAAAvAYj4qgTbrfVgx9s1MpdBXppVE8NSEl0OhIAAIBXYUQcdeK5r9M0e8MePTS0g4b1THY6DgAAgNehiKPWzVy5W5MW79R1A1rozsFtnI4DAADglSjiqFWL0/L0108367wODTX+N+wVDgAA8FMo4qg12QVluu/d9WrfOFYTr+utkGD+9QIAAPgpNCXUiiNV1bp75lq5rdWrN/RWdDjrgAEAAH4ObQm14qnPtmpjTpEm39hHLROjnY4DAADg9RgRx2n7dH2upn+XpXGDUnRRlyZOxwEAAPAJFHGclvS8w/rTR5vUr1WCHry4g9NxAAAAfAZFHDVWVlmlO99Zq8jQYP372t4KZXEmAADAKWOOOGrEWqtHPtqk9AMlemfMADWJj3A6EgAAgE9hCBM1MnPVbn2yfo9+P6S9zmrbwOk4AADg/7V371F2loW5wJ93JplcyR0JQgIIQS4BCeZgqfdiK6Gr0mVtlVOPttLqWdVebJcW2x57alutdrlsu6TtsUes2lakWiw9DcW20qNSVLBBIGAkAiaBhEBu5DqTmXnPHzN6YhqSIczsb8/ev99as7Iv3977WTPv2vvJu9/v+5h0zIjztO3tH8wf33RfXnr2iXnry89qOg4AwKRkRpynZXC45oGte7Jodl8+9NqL0tPjzJkAAMdDEWfMaq359tY9GRgczod/+uIsmNXXdCQAgElLEWfMPvLFB7Nj30CWLpiZi5fObzoOAMCkpogzJl99cFs+cMu6LJw1zRFSAADGgZ01OabHd/fnFz+1JksXzMxzFjh9PQDAeDAjzlENDdf80qfW5MkDB/Nnr784vXbOBAAYF4o4R/Whf/5Wbn9wW37vxy/IOYvnNB0HAKBjKOI8pVu/uTUfvnV9XrtySV7z/FObjgMA0FEUcY5o0459efsNd+Xck+fkd648v+k4AAAdRxHnP+kfHMpb/2ZNhoZq/uynL870qb1NRwIA6DiOmsJ/8t5/vD/f2Lgzf/76i3P6IkdJAQCYCGbE+T7/8I1H8/Hbv5Ofe9EZuXz5yU3HAQDoWIo437N+655c89m78/zT5ufXV53TdBwAgI6miJMk2TcwmF/4669n2tTefPi/rsjUXkMDAGAiWSNOaq35rRvvzQNb9+QTb7okJ8+d0XQkAICOZ9qTXH/Hxvzdmkfyy5cty4uXndh0HACArqCId7l7H9mV375pbV68bFF+8YeWNR0HAKBrKOJdbNf+g/mFv/6PLJjZlz967UXp7SlNRwIA6BrWiHepWmve8bffyKM79+fTb/mBLJw9relIAABdxYx4l/rfX3oon7/vsVyz6pw8/7QFTccBAOg6ingXuuPh7fmDf/pmLj9/ca5+0RlNxwEA6EqKeJfZ2z+YX73hrpw6f0Y+8JMXphTrwgEAmmCNeJf5w1vWZeP2/bnhLZdmzvSpTccBAOhaZsS7yJ0Pb8/Hb384b7z0tFxyhnXhAABNUsS7xIGDQ3nnZ+7Os+fOyDsvP6fpOAAAXc/SlC7xR//yQB58Ym/+6uoXZNY0f3YAgKaZEe8C39i4Mx/54rfz2pVL8qJli5qOAwBAFPGONzA4nHd+5u6ceMK0/MaPntt0HAAARlmj0OGuvXV91j22Ox9948rMneEoKQAA7cKMeAfbNzCUa29dnx+/6Nm57NyTmo4DAMAhFPEOVWvy4ON7Mm/m1Pz2j53fdBwAAA6jiHeozbv2Z0//YH7nVcszf1Zf03EAADiMIt6B1m/dk0079mfBrL5cccHipuMAAHAEiniHGRqu+fXP3p2enpIzFs1KKaXpSAAAHIEi3mE+cfvD+fp3duT0hTMztdefFwCgXWlqHWTDtn35wD+ty8ufe2IWzZ7WdBwAAI5CEe8QtY4sSZnSU/LeV1/QdBwAAI5BEe8Qn/raxtz+4La864pzc/LcGU3HAQDgGBTxDvDozv157+r784NnLsxVlyxpOg4AAGOgiE9ytdb85o33ZGi45g9efaGjpAAATBKK+CR345pHcuu6x/POy5+bpQtnNh0HAIAxUsQnsa27D+R3/uG+rDxtft546elNxwEA4GlQxCexd39ubfYfHMr7X3NhenosSQEAmEwU8Ulq9T2b809rt+Ttrzg7Z544u+k4AAA8TYr4JLR970De/ff35oJT5ubnX3xG03EAADgOU5oOwNP3nn9Ym137D+avfu4FmeI09gAAk5IWN8n86/2P5XN3PZq3vvysnLN4TtNxAAA4Tor4JLJr/8H8xo335JzFJ+QXXnZW03EAAHgGLE2ZRN63+v48vrs/f/GGlemb4v9QAACTmTY3SXzpgcdz/R0b8+aXnJkLT53XdBwAAJ4hRXwS2Ns/mGs+e0+ec+Ks/MorljUdBwCAcWBpyiTwJ194II/u2p+/fculmT61t+k4AACMAzPibe6RnfvzsdsezqtXnJqVpy9oOg4AAONEEW9zH/z8uiTJr/3I2Q0nAQBgPCnibWztvfrWWAAAFAxJREFUo7ty45pH8rMvPD3Pnjej6TgAAIwjRbyN/cHN38zcGVMdMxwAoAMp4m3qSw88ni898ETe9vKzMnfG1KbjAAAwzhTxNvW+1d/MqfNn5L9delrTUQAAmACKeBt6Yk9/7tv8ZN7xyudm2hSHKwQA6ESKeJsZrjUbt+/PBafMzY9d+Oym4wAAMEEU8Tbz2JMH0j84lHddcU56ekrTcQAAmCCKeBvZuW8gj+zYn3kz+/KDZy5qOg4AABNIEW8j1966PkO1ZumCmU1HAQBgginibWLj9n35+L9/JyfOnpaZfXbQBADodIp4m/jg59ellOTU+WbDAQC6gSLeBu59ZFc+d9ejufpFZ6Rvij8JAEA30PoaVmvN+26+P/NnTs1/f9mZTccBAKBFFPGGffGBJ3Lb+m35pcuWZc50p7IHAOgWiniDhoZr3rf6/ixdMDM//QKnsgcA6CaKeINuXPNIvrlld97xyudaGw4A0GW0v4YcODiUD35+XZ536tz86AUnNx0HAIAWU8Qb8rHbHs7mXQdyzapzncoeAKALKeIN2L53IH966/pcds6zcumZC5uOAwBAAxTxBnz4C+uzd2Awv77qnKajAADQEEW8xTZu35dPfuXh/NTKJTn7pBOajgMAQEMU8Rb7w1vWpben5O0/fHbTUQAAaFDLingp5fJSyrpSyvpSyjVHuH9aKeXTo/d/tZRy+ujtP1xK+Xop5Z7Rf3+oVZnH292bduambzyan3vRc3LSnOlNxwEAoEEtKeKllN4k1yZZleS8JFeVUs47bLOrk+yotZ6V5ENJ3j96+xNJfqzWekGSNyb5ZCsyj7daa967+v4smNWXt7z0OU3HAQCgYa2aEb8kyfpa64O11oEk1ye58rBtrkzy8dHLn0lyWSml1FrX1FofHb19bZLppZRpLUk9jv5t3eP5yoPb88uXLcsJTmUPAND1WlXET0my8ZDrm0ZvO+I2tdbBJLuSHH5sv59IsqbW2j9BOSfE0HDN+26+P6cvnJmrLlnadBwAANpAq4r4kc5YU5/ONqWU8zOyXOUtR3yBUt5cSrmzlHLncaecIJ/9+qZ867E9eefl5ziVPQAASVpXxDclWXLI9VOTPPpU25RSpiSZm2T76PVTk9yY5A211m8f6QVqrR+pta6sta4c5+zPyMDgcD70L9/KRUvmZdXyxU3HAQCgTbSqiN+RZFkp5YxSSl+S1yW56bBtbsrIzphJ8pokX6i11lLKvCT/mORdtdbbWpR33Nx87+Zs3nUgv3zZspTiVPYAAIxoSREfXfP9tiS3JLk/yQ211rWllPeUUl41utlHkywspaxP8qtJvnuIw7clOSvJ/yil3DX686xW5H6maq356JcfynNOnJWXnn1i03EAAGgjU1r1QrXW1UlWH3bbuw+5fCDJTx7hcb+X5PcmPOAE+I8NO3L3pl353R9fnp4es+EAAPx/9hycQNd9+eHMmT4lP3Hx4QeIAQCg2yniE2TTjn25+d7NueoFSzOzr2VfPAAAMEko4hPkk7d/J6WUvOHS05uOAgBAG1LEJ8DQcM2nvrYhly9fnFPmzWg6DgAAbUgRnwBP7OnPkwcG86YXntF0FAAA2pQiPgG27DqQ5y2Zl4uXzms6CgAAbUoRH2c79w1k/8GhvOmFpzuBDwAAT0kRH2ebdx1I35SeXHHByU1HAQCgjSni42jdlt3Ztf9gTpozPVN7/WoBAHhq2uI4+thtD6WnlJx0wvSmowAA0OYU8XGybU9//m7NIznxhGmZ0mttOAAAR6eIj5NPfW1DBgaHs3iO2XAAAI5NER8HA4PD+cTt38mLly3KjL7epuMAADAJKOLjYPU9m7N1d3+ufpET+AAAMDaK+DNUa811tz2UM0+clZcsO7HpOAAATBKK+DN053d25O5Nu/KzLzwjPT120gQAYGwU8Wfoui8/lLkzpubVF5/SdBQAACYRRfwZ2Lh9X25ZuyVXXbI0M/umNB0HAIBJRBF/Bj5x+8MppeQNl57WdBQAACYZRfw47ekfzPV3bMyq5Yvz7Hkzmo4DAMAko4gfp89+fVN2HxjMmxyyEACA46CIH4fh4ZqP3fZQLloyLxcvnd90HAAAJiFF/Djcum5rHt62z2w4AADHTRE/Dtfd9lAWz5meVcsXNx0FAIBJShF/mr655cnctn5b3vCDp2Vqr18fAADHR5N8mj725YczfWpPrvovS5uOAgDAJKaIPw079w3kxrseyasvPjXzZ/U1HQcAgElMEX8aPr/2sQwMDpsNBwDgGVPEn4bV927OkgUzsvyUOU1HAQBgklPEx2jXvoO5bf0TuWL5ySmlNB0HAIBJThEfo3++/7EcHKpZdcHJTUcBAKADKOJjdPM9m3PKvBl53qlzm44CAEAHUMTHYGi45ksPPJFVyxdblgIAwLhQxMdgx76BDAwNW5YCAMC4UcTHYPvegSyeMz0rlsxrOgoAAB1CET+GoeGanfsOZtUFi9PTY1kKAADjQxE/hp37Dma41lxhWQoAAONIET+GbXv709fbk+cvnd90FAAAOogifhR7+wezc9/BLJjVZ1kKAADjShE/in9b93iGa82C2X1NRwEAoMMo4kex+t7NmdrbkxOmT206CgAAHUYRfwr7B4byhfu3ZsGsvliUAgDAeFPEn8L//dbW7D84lIWzLEsBAGD8KeJPYfU9W7JgVp9lKQAATAhF/AgOHBzKv97/WF55/kkp1qUAADABFPEj+OK3Hs/egaGsWu4kPgAATAxF/AhuvndL5s2cmkvPXNh0FAAAOpQifpj+waH8y32P5UfOOylTe/16AACYGJrmYW5b/0R29w9m1QWWpQAAMHEU8cP8491bMmf6lLzwzEVNRwEAoIMp4ocYGBzOP9+3JT983uL0TfGrAQBg4mibh/j3bz+RJw8M5ooLFjcdBQCADqeIH+Lme7bkhGlT8qJllqUAADCxFPFRB4eGc8t9W/KK807KtCm9TccBAKDDKeKjvvLgtuzcdzCrlluWAgDAxFPER62+Z0tm9fXmJWef2HQUAAC6gCKeZGi45vNrt+Syc0/K9KmWpQAAMPEU8SQPbN2dbXsH8rLnmg0HAKA1FPEkazbsTJJcvHR+w0kAAOgWiniSNRt2ZP7MqTlt4cymowAA0CUU8YzMiK9YOj+llKajAADQJbq+iA8O1zywdU9WLJnXdBQAALpI1xfxvf2DSZIV1ocDANBCXV/E9xwYTCnJhUvmNh0FAIAuooj3D2bZs2ZnzvSpTUcBAKCLdH0R390/mBVLLEsBAKC1urqIHzg4lMGh4axYakdNAABaq6uL+B47agIA0JDuLuIHBtPbU3LWs2Y3HQUAgC7T1UV8d/9gZk+bkt4eJ/IBAKC1uraI7x8Yyr6BocyeNqXpKAAAdKGuLeL3PLIrtdbMnq6IAwDQel1bxNds2JEkZsQBAGhEFxfxnZk+tTdTe7v2VwAAQIO6soXWWvMfG3aYDQcAoDFdWcQ37zqQrbv7rQ8HAKAxXVnE12zYmSQ5wYw4AAAN6dIiviPTpvRkZp8iDgBAM7qziG/cmQtOmZviPD4AADSk64r4wOBw7nlkV1Ysndd0FAAAuljXFfH7Nz+ZgcHhrFg6v+koAAB0sa4r4t89kY8ZcQAAmtR9RXzjziyeMz0nz53RdBQAALpY9xXxDTvNhgMA0LiuKuIHh4azYfs+RRwAgMZ1VRHf0z+YJHbUBACgcd1VxA8MZkpPyfJnz206CgAAXa67inj/YM49eU5m9PU2HQUAgC7XNUW8ZqSIWx8OAEA76Joivn9gKEPDVREHAKAtdE0R/96OmkvsqAkAQPO6p4gfOJipvT05beHMpqMAAEAXFfH+wcyeNiWllKajAABAdxTxXfsPZt/AUGZPm9J0FAAASNIlRfzuTTuTJLOnK+IAALSHrijiazbsTCkxIw4AQNvokiK+IzOm9qa3x/pwAADaQ3cU8Y07M3va1KZjAADA93R8ET9wcCg79x20PhwAgLbS8UX8uyfysT4cAIB20vlF/MDI8cNn9PU2HQUAAL6n44v47v7BPG/J3NhNEwCAdtLxRXzfwFBWLJnfdAwAAPg+HV3Eh4Zraq1ZsXRe01EAAOD7dHwRT5KLlijiAAC0l44v4tOn9mbh7GlNRwEAgO/T2UW8VoctBACgLXV0ER+u1Yl8AABoSx1bxOvI8vCcYEYcAIA21LFFPElKSmb2KeIAALSfji3iNTW9PSXFmXwAAGhDnVvEa9KrhQMA0KY6tognSW+PIg4AQHtSxAEAoAGdW8RLrA8HAKBttayIl1IuL6WsK6WsL6Vcc4T7p5VSPj16/1dLKacfct+7Rm9fV0p55Zheb/yiAwDAuGtJES+l9Ca5NsmqJOcluaqUct5hm12dZEet9awkH0ry/tHHnpfkdUnOT3J5kj8dfb6jv6YqDgBAG2vVjPglSdbXWh+stQ4kuT7JlYdtc2WSj49e/kySy0opZfT262ut/bXWh5KsH32+o9PDAQBoY60q4qck2XjI9U2jtx1xm1rrYJJdSRaO8bEppby5lHJnKeXORA8HAKC9taqIH6kX1zFuM5bHptb6kVrrylrryuPIBwAALdWqIr4pyZJDrp+a5NGn2qaUMiXJ3CTbx/hYAACYVFpVxO9IsqyUckYppS8jO1/edNg2NyV54+jl1yT5Qq21jt7+utGjqpyRZFmSr7UoNwAATIgprXiRWutgKeVtSW5J0pvkulrr2lLKe5LcWWu9KclHk3yylLI+IzPhrxt97NpSyg1J7ksymOSttdahVuQGAICJ0pIiniS11tVJVh9227sPuXwgyU8+xWN/P8nvT2hAAABooc49syYAALQxRRwAABqgiAMAQAMUcQAAaIAiDgAADVDEAQCgAYo4AAA0QBEHAIAGKOIAANAARRwAABqgiAMAQAMUcQAAaIAiDgAADVDEAQCgAYo4AAA0QBEHAIAGKOIAANAARRwAABqgiAMAQAMUcQAAaIAiDgAADVDEAQCgAYo4AAA0QBEHAIAGlFpr0xnGXSlld5J1Teeg7SxK8kTTIWg7xgVHYlxwJMYFR/LcWusJx/PAKeOdpE2sq7WubDoE7aWUcqdxweGMC47EuOBIjAuOpJRy5/E+1tIUAABogCIOAAAN6NQi/pGmA9CWjAuOxLjgSIwLjsS44EiOe1x05M6aAADQ7jp1RhwAANqaIg4AAA2Y1EW8lHJ5KWVdKWV9KeWaI9w/rZTy6dH7v1pKOb31KWm1MYyLXy2l3FdKubuU8q+llNOayElrHWtcHLLda0optZTiEGVdYCzjopTyU6PvGWtLKX/T6oy03hg+R5aWUm4tpawZ/Sy5oomctE4p5bpSytZSyr1PcX8ppfzJ6Ji5u5Ry8Vied9IW8VJKb5Jrk6xKcl6Sq0op5x222dVJdtRaz0ryoSTvb21KWm2M42JNkpW11guTfCbJB1qbklYb47hIKeWEJL+U5KutTUgTxjIuSinLkrwryQtrrecn+ZWWB6Wlxvh+8VtJbqi1rkjyuiR/2tqUNOAvk1x+lPtXJVk2+vPmJH82liedtEU8ySVJ1tdaH6y1DiS5PsmVh21zZZKPj17+TJLLSimlhRlpvWOOi1rrrbXWfaNXv5Lk1BZnpPXG8n6RJL+bkf+YHWhlOBozlnHx80murbXuSJJa69YWZ6T1xjIuapI5o5fnJnm0hfloQK31i0m2H2WTK5N8oo74SpJ5pZSTj/W8k7mIn5Jk4yHXN43edsRtaq2DSXYlWdiSdDRlLOPiUFcnuXlCE9EOjjkuSikrkiyptf6fVgajUWN5vzg7ydmllNtKKV8ppRxtRozOMJZx8T+TvL6UsinJ6iS/2JpotLGn2z+STO5T3B9pZvvwYzGOZRs6y5j/5qWU1ydZmeSlE5qIdnDUcVFK6cnI8rWfaVUg2sJY3i+mZOSr5pdl5NuzL5VSltdad05wNpozlnFxVZK/rLV+sJRyaZJPjo6L4YmPR5s6rs45mWfENyVZcsj1U/Ofvxr63jallCkZ+froaF8rMPmNZVyklPKKJL+Z5FW11v4WZaM5xxoXJyRZnuTfSikPJ/mBJDfZYbPjjfVz5O9rrQdrrQ8lWZeRYk7nGsu4uDrJDUlSa709yfQki1qSjnY1pv5xuMlcxO9IsqyUckYppS8jO0vcdNg2NyV54+jl1yT5QnUGo053zHExugThf2WkhFvv2R2OOi5qrbtqrYtqrafXWk/PyL4Dr6q13tlMXFpkLJ8jn0vy8iQppSzKyFKVB1uaklYby7jYkOSyJCmlnJuRIv54S1PSbm5K8obRo6f8QJJdtdbNx3rQpF2aUmsdLKW8LcktSXqTXFdrXVtKeU+SO2utNyX5aEa+LlqfkZnw1zWXmFYY47j4wySzk/zt6L67G2qtr2osNBNujOOCLjPGcXFLkh8ppdyXZCjJO2qt25pLzUQb47j4tSR/UUp5e0aWH/yMib7OVkr5VEaWqC0a3Tfgt5NMTZJa659nZF+BK5KsT7Ivyc+O6XmNGwAAaL3JvDQFAAAmLUUcAAAaoIgDAEADFHEAAGiAIg4AAA1QxAEAoAGKOABJklLKI6WUi5rOAdAtHEccgO+eNXJzktm11v6m8wB0AzPiAF2ulHJWko0Z+UzYVkrZVkqZtGdeBpgszIgDkFLKW5O8pNb62qazAHQLM+IAJMnzktzVdAiAbqKIA5AkFyX5RtMhALqJpSkAXa6U0pNkT5Kzaq2PNp0HoFuYEQdgxuiPzwSAFvKmC9Dlaq17k/x5kvtKKZuazgPQLSxNAQCABpgRBwCABijiAADQAEUcAAAaoIgDAEADFHEAAGiAIg4AAA1QxAEAoAGKOAAANOD/Ae/rO5xPEahPAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x864 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig,ax = plt.subplots(figsize=(12,12))\n",
    "\n",
    "ax.plot(pc.t_grid, l_star, label=r'$\\ell^*$')\n",
    "ax.set_xlim(0.0, 1.0)\n",
    "for s in t:\n",
    "    ax.axvline(x=s, c=\"0.5\")\n",
    "    \n",
    "ax.set_ylabel(r'$\\ell$', fontsize=12)\n",
    "ax.set_xlabel(r'$t$', fontsize=12)\n",
    "ax.legend(fontsize=12)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compare to johns code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "from scipy.optimize import fminbound\n",
    "from interpolation import interp"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "class ProductionChain_j:\n",
    "\n",
    "    def __init__(self,\n",
    "            n=1000,\n",
    "            delta=1.05,\n",
    "            c=lambda t: np.exp(10 * t) - 1):\n",
    "\n",
    "        self.n, self.delta, self.c = n, delta, c\n",
    "        self.grid = np.linspace(1e-04, 1, n)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_prices(pc, tol=1e-5, max_iter=5000):\n",
    "    \"\"\"\n",
    "    Compute prices by iterating with T\n",
    "\n",
    "        * pc is an instance of ProductionChain\n",
    "        * The initial condition is p = c\n",
    "\n",
    "    \"\"\"\n",
    "    delta, c, n, grid = pc.delta, pc.c, pc.n, pc.grid\n",
    "    p = c(grid)  # Initial condition is c(s), as an array\n",
    "    new_p = np.empty_like(p)\n",
    "    error = tol + 1\n",
    "    i = 0\n",
    "\n",
    "    while error > tol and i < max_iter:\n",
    "        for i, s in enumerate(grid):\n",
    "            Tp = lambda t: delta * interp(grid, p, t) + c(s - t)\n",
    "            new_p[i] = Tp(fminbound(Tp, 0, s))\n",
    "        error = np.max(np.abs(p - new_p))\n",
    "        p = new_p\n",
    "        i = i + 1\n",
    "\n",
    "    if i < max_iter:\n",
    "        print(f\"Iteration converged in {i} steps\")\n",
    "    else:\n",
    "        print(f\"Warning: iteration hit upper bound {max_iter}\")\n",
    "\n",
    "    p_func = lambda x: interp(grid, p, x)\n",
    "    return p_func"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "def optimal_choices(pc, p_function, s):\n",
    "    \"\"\"\n",
    "    Takes p_func as the true function, minimizes on [0,s]\n",
    "\n",
    "    Returns optimal upstream boundary t_star and optimal size of\n",
    "    firm ell_star\n",
    "\n",
    "    In fact, the algorithm minimizes on [-1,s] and then takes the\n",
    "    max of the minimizer and zero. This results in better results\n",
    "    close to zero\n",
    "\n",
    "    \"\"\"\n",
    "    delta, c = pc.delta, pc.c\n",
    "    f = lambda t: delta * p_function(t) + c(s - t)\n",
    "    t_star = max(fminbound(f, -1, s), 0)\n",
    "    ell_star = s - t_star\n",
    "    return t_star, ell_star"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def compute_stages(pc, p_function):\n",
    "    s = 1.0\n",
    "    transaction_stages = [s]\n",
    "    while s > 0:\n",
    "        s, ell = optimal_choices(pc, p_function, s)\n",
    "        transaction_stages.append(s)\n",
    "    return np.array(transaction_stages)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Iteration converged in 1000 steps\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAD5CAYAAAAqaDI/AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dd3xV9f3H8dcnYYYNYYQNsocCIopoHSgKYnELWhUcqHW0rlatrVb9WWtb90SliFrELVVUQFFwYNmyIbISAgQICYGQcZPP749ca4oJhNwk9yb3/Xw88uCec77nnA/fR5J3zvoec3dERCR6xYS7ABERCS8FgYhIlFMQiIhEOQWBiEiUUxCIiEQ5BYGISJSrcagGZtYOmAy0AgqACe7+hJk1BaYCHYGNwEXuvruY9a8A7glOPujurxxqn3Fxcd6yZUuaNWvGrl27AGjWrBnAz6Z/VNL80i4vbZuytA1lnVDWK6/1y3s75b2titxmRW63srZf2fuJlv1Fyr5/3P/GjRt3unvzsqxvh3qOwMwSgAR3X2RmDYCFwDnAWCDN3R82szuBJu7++wPWbQosAAYCHlz36OICo6jWrVv7Qw89xNixY5k0aRIAY8eOBfjZ9I9Kml/a5aVtU5a2oawTynrltX55b6e8t1WR26zI7VbW9it7P9Gyv0jZ94/7Hzdu3EJ3H1iW9Q95asjdt7r7ouDnTGAV0AYYBfz41/0rFIbDgc4AZrp7WvCX/0zgzLIUKiIiFeOwrhGYWUegP/Ad0NLdt0JhWAAtilmlDZBUZDo5OE9ERCJEqYPAzOoD7wC/dfc9pV2tmHnFnosys/FmtsDMFpS2JhERCV2pgsDMalIYAq+7+7vB2duD1w9+vI6QWsyqyUC7ItNtgZTi9uHuE9x9YFnPcYmISNkcMgjMzICXgVXu/miRRdOAK4KfrwA+KGb1T4FhZtbEzJoAw4LzREQkQpTmiGAIcBlwqpktCX6NAB4GTjezdcDpwWnMbKCZvQTg7mnAA8D84Nf9wXkiIhIhDvkcgbt/RfHn+gGGFtN+AXB1kemJwMSyFigiIhVLTxaLiEQ5BYGISBWWX+BsTssKaRuHPDUkIiKRaW9OgN9MWUzN9P0hbUdHBCIiVVBSWhbnP/sNX6zdQaf4eiFtS0EgIlLFzN+YxqhnvmZrxn5eGTeIlg3rhLQ9BYGISBXy5oIkLnlxHo3r1uT9G4ZwQtf4kLepawQiIlVAfoHz8MereHHuBk7oEs8zlwygUVzNctm2gkBEJMJlZufxmzeW8PnqVK4Y3IE/juxFjdjyO6GjIBARiWCbd2Vx9eT5/LBjHw+c04fLjutQ7vtQEIiIRKivE3dyw78W4Q6TrxzEkC6hXw8ojoJARCTCuDuvfLORBz5aRef4erx0xUA6NAvtFtGDURCIiESQ3EABf/pgOW/MT+K0ni147OJ+NKhTPheFS6IgEBGJEDsyc7j+tYUs2LSbG0/pwq2ndyMmpqQxP8uPgkBEJAIs35LB+MkLSMvK5akx/Tn7qNaVtm8FgYhImP17aQp3vL2UpnG1ePu64+nTplGl7l9BICISJgUFzqMz1/L07EQGdmjCc786muYNald6HQoCEZEwyMzO45apS5i1KpWLB7bj/nN6U7tGbFhqURCIiFSyjTv3cc3kBazfuY/7zu7FFcd3pPD18OFxyCAws4nASCDV3fsE500FugebNAbS3b1fMetuBDKBfCDg7gPLqW4RkSrpq3WFD4mZVexDYoejNEcEk4Cngck/znD3i3/8bGb/ADIOsv4p7r6zrAWKiFQH7s6kbzby4EerOKJ5PV68vGIfEjscpXl5/Rwz61jcMis8lrkIOLV8yxIRqT5yAvn86f0VTF2QxGk9W/L46H7Urx05Z+ZDreREYLu7rythuQMzzMyBF9x9Qoj7ExGpUrZlZHPdawtZkpReqQ+JHY5Qg2AMMOUgy4e4e4qZtQBmmtlqd59TXEMzGw+MB0hISAixLBGR8FuwMY3rX1/EvpwAz146gBF9I/N3W5kHtDazGsB5wNSS2rh7SvDfVOA9YNBB2k5w94G6oCwiVZ278+q8TYyeMI96tWJ5/4YhERsCENoRwWnAandPLm6hmdUDYtw9M/h5GHB/CPsTEYl42Xn53PtB4fWAk7s354mL+5fbm8QqyiGPCMxsCvAt0N3Mks3squCi0RxwWsjMWpvZ9OBkS+ArM1sK/Af4yN0/Kb/SRUQiy7aMbC6eMI+pC5K48ZQuvHzFMREfAlC6u4bGlDB/bDHzUoARwc/rgaNCrE9EpErIzA4w8qmv2J8b4PlfDeDMPpF7KuhAkXP/kohIFeTubN+TzcZdWTRoUIMp1xxL15YNwl3WYVEQiIiUUXZePn98fzmZO/fRJK4W798whEZ1I/9U0IEUBCIiZbA1Yz/XvbaIpUnpXN2yLm2bxFXJEAAFgYjIYfvPhjR+/fpC9ufm8/yvjmbbgh3hLikkZX6OQEQk2rg7k7/dyCUvzqNhnZp8cOMQzuzTKtxlhUxHBCIipfDj9YC3FiYztEcLHhvdj4YV/FL5yqIgEBE5hKS0LH79+iKWbcng5qFd+e3QrhE3XlAoFAQiIgfx5dod/OaNxeQXOC9dPpDTerUMd0nlTkEgIlKMggLn6dmJPDZrLd1bNuD5Xx1Nx/jIeH9AeVMQiIgcICMrj1veXMLnq1M5t38bHjq3L3Vrhed9wpVBQSAiUsSKlAyuf20RWzP288Co3vzquA5hfZ9wZVAQiIgEvbMwmbvfW0aTuFpMvXYwA9o3CXdJlUJBICJRLyeQzwMfruS1eZsZ3LkZT13Sn/j6tcNdVqVREIhIVEtJ38/1rxcOFXHtSZ25Y1h3asRG17O2CgIRiVpfJ+7kpimLyQ0U8NylAxgewW8Rq0gKAhGJOu7O81+u52+fruaI5vV5/rKjOaJ5/XCXFTYKAhGJKnuy87jjraV8umI7I49M4K/nH0m92tH9qzC6//ciElXWbMvkutcWsjktiz+O7MWVQzpW+1tDS6M07yyeaGapZra8yLz7zGyLmS0Jfo0oYd0zzWyNmSWa2Z3lWbiIyOF4e2Eyo575ir05AaZccxxXndBJIRBUmiOCScDTwOQD5j/m7n8vaSUziwWeAU4HkoH5ZjbN3VeWsVYRkcOWnZfPvR+sYOqCJI7r3JQnx/SnRYM64S4ropTm5fVzzKxjGbY9CEgMvsQeM3sDGAUoCESkUmzYuY/rX1vI6m2Z3HDKEdxyWreouzW0NELpkRvN7PvgqaPiHr9rAyQVmU4OziuWmY03swVmtiCEmkREAJi+bCtnP/UV2/Zk889xx3DHGT0UAiUoa688BxwB9AO2Av8opk1xJ9+8pA26+wR3H+juA8tYk4gIuYEC7pu2gl+/voguLerz0c0nckr3FuEuK6KV6a4hd9/+42czexH4sJhmyUC7ItNtgZSy7E9EpDRyAgVc9MK3LElKZ9yQjtw1vCe1augo4FDKFARmluDuW4OT5wLLi2k2H+hqZp2ALcBo4JIyVSkicgjpWXkk7thLYsFenr10ACOi9CnhsjhkEJjZFOBkIN7MkoF7gZPNrB+Fp3o2AtcG27YGXnL3Ee4eMLMbgU+BWGCiu6+okP+FiEStQH4Bj81ayw/b9hBXqwb/vuEEOlXTF8hUlNLcNTSmmNkvl9A2BRhRZHo6ML3M1YmIHERqZjY3T1nMvPVpXN6sDh3j4xQCZaAni0WkSpq3fhc3TVlMZnYef7/wKPYu2xHukqosXUURkSqloMB59otELnlxHg1q1+D9G4ZwwdFtw11WlaYjAhGpMnbuzeHWN5cyZ+0ORh6ZwMPnH0n9KB8wrjyoB0WkSpi3fhc3T1lM+v48HjynD5ce215jBZUTBYGIRLT8AufpzxN54rO1dGxWj0njBtGrdcNwl1WtKAhEJGKl7snmt1OX8M0Puzi3fxsePKdP1L87oCKoR0UkIs1dt4Nbpi5hb06ARy44kguPbqtTQRVEQSAiESWQX8Djs9bxzBeJdG1RnynXHEfXlg3CXVa1piAQkYixNWM/N09ZzPyNuxl9TDvuPbs3dWvFhrusak9BICIR4fPV27ntzaXkBgp4YnQ/RvUrcdR6KWcKAhEJq9xAAX+fsYYJc9bTK6EhT1/Sn87N64e7rKiiIBCRsElKy+KmKYtZkpTOZcd14A9n9aROTZ0KqmwKAhEJi0+Wb+N3by/FHQ0bHWYKAhGpVNl5+fzfR6t4dd4mjmrbiKfGDKB9s7hwlxXVFAQiUmnWbMvk5imLWbM9k/G/6Mztw7rrDWIRQEEgIhXO3Xntu808+OFKGtSpyStXDuKkbs3DXZYEKQhEpELt3pfL79/5nhkrt3NSt+b8/cKjaN6gdrjLkiJK86rKicBIINXd+wTn/Q04G8gFfgDGuXt6MetuBDKBfCDg7gPLr3QRiXR79ucx/Im57NqXwz1n9eTKIZ2IidEwEZGmNCfnJgFnHjBvJtDH3Y8E1gJ3HWT9U9y9n0JAJHoE8gtI2p3Fqm17iKsVy3u/HsLVJ3ZWCESoQwaBu88B0g6YN8PdA8HJeYBeDyQiQOGzARe98C1bdu+nef3a/PumE+jTplG4y5KDKI9rBFcCU0tY5sAMM3PgBXefUA77E5EI9e+lKdz97jIAhrRoQLP6tTRsdBUQ0n1bZvYHIAC8XkKTIe4+ABgO3GBmvzjItsab2QIzWxBKTSJS+fblBLjjraXcNGUxXVrWZ/pvTqRZ/VrhLktKqcxRbWZXUHgReai7e3Ft3D0l+G+qmb0HDALmlNB2AjABoHXr1sVuT0Qiz/ItGdw8ZTEbdu3jxlO68JvTulIzVs8GVCVlCgIzOxP4PXCSu2eV0KYeEOPumcHPw4D7y1ypiESUggJn4tcb+Osnq2lWrzavX30sxx8RH+6ypAxKc/voFOBkIN7MkoF7KbxLqDYwM/jGoHnufp2ZtQZecvcRQEvgveDyGsC/3P2TCvlfiEil2r4nm9veXMpXiTs5vVdLHjn/SJrU06mgquqQQeDuY4qZ/XIJbVOAEcHP64GjQqpORCLOx8u2ctd7y8jJK+Chc/syZlA7vUKyitPlfBEplb05Af48bQVvLUzmyLaNePzifnpvQDWhIBCRQ1q0eTe3TF1CUlqWLghXQwoCESlRIL+Apz5P5OnZibRqWIc3xg9mUKem4S5LypmCQESKtWnXPn47dQmLN6dzbv82/HlUbxrWqRnusqQCKAhE5H+4O28tTObP01YQG2M8OaY/vzyqdbjLkgqkIBCR/9q9L5e731vGx8u3cVznpjx6UT9aN64b7rKkgikIRASAr9bt5La3lpC2L5c7h/fgmhM7E6vRQqOCgkAkymXn5fO3T9fw8lcb6NKiPi9fcYxGC40yCgKRKLYiJYNbpy5lzfZMLh/cgbuG96RurdhwlyWVTEEgEoUC+QW8MGc9j89aS5O4Wvxz3DGc0r1FuMuSMFEQiESZDTv3cdubS1i0OZ2zjkzgwVF9NE5QlFMQiESR7XuyGfHEXGrG6rZQ+YmCQCQKbMvIZvXWPaTvz+OYTk155PwjadWoTrjLkgihwUJEqjF354MlWxj22Jdk5gToFF+PV8YdoxCQ/6EjApFqave+XO75YDkffb+VAe0b07dmI+rUjNWQ0fIzOiIQqYY+X72dYY/PYcaKbfzuzO68dd3x1Kmp20KleDoiEKlG9uUEePCjVUz5z2Z6tGrAK+MG0at1w3CXJRFOQSBSTczfmMZtby4laXcW157UmVtP70btGjoKkEMr1akhM5toZqlmtrzIvKZmNtPM1gX/bVLCulcE26wzsyvKq3ARKZSdl8+DH67kohe+BeDNawdz1/CeCgEptdJeI5gEnHnAvDuBz9y9K/BZcPp/mFlTCl92fywwCLi3pMAQkcO3cNNuRjwxl5e+2sClx7bn49+cyDEd9eIYOTylOjXk7nPMrOMBs0cBJwc/vwJ8Afz+gDZnADPdPQ3AzGZSGChTylStiACFRwGPzlzLS3PXk9CoLq9ffSxDusSHuyypokK5RtDS3bcCuPtWMytuoJI2QFKR6eTgPBEpo0Wbd3P7W0tZv2MflxzbnrtH9KR+bV3uk7Kr6O+e4m5Y9mIbmo0HxgMkJCRUZE0iVVJ2Xj6PzVrLi3PW06phHV69ahAndm0e7rKkGgglCLabWULwaCABSC2mTTI/nT4CaEvhKaSfcfcJwASA1q1bFxsWItFqSVI6t7+1lMTUvYwZ1I67R/Skgd4fLOUklAfKpgE/3gV0BfBBMW0+BYaZWZPgReJhwXkiUgo5gXz++slqznv2a/blBHjlykH85bwjFQJSrkp1RGBmUyj8yz7ezJIpvBPoYeBNM7sK2AxcGGw7ELjO3a929zQzewCYH9zU/T9eOBaRg1saPApYl7qXiwe24w8je9JQASAVoLR3DY0pYdHQYtouAK4uMj0RmFim6kSiUE4gnyc/W8fzX66nef3aTBp3DCfrpTFSgXSrgUgE+T45nTve+p412zO58Oi23DOyF43q6ihAKpaCQCQCFLjzl+mreHHuepo3qM0/xx7DKT10FCCVQ0EgEmZ7svPYsGMf761Zz5hB7bhzeE8dBUilUhCIhElmdh6PfLKGHSl7qFMzln9dfSzH6+lgCQMFgUgYzF6Tyh/eXcbWPdmMja9Lu6Z1FQISNnoxjUgl2r0vl1unLmHcP+dTr3YN3rn+eDo0iyNGbw2TMNIRgUglcHemL9vGvdOWk56Vx82nduGGU7tQu0Ys34e7OIl6CgKRCpa6J5s/frCcT1dsp2+bRrx61bH0TNBbwyRyKAhEKoi789bCZB78cCU5gQLuHN6Dq0/oRI1YnZGVyKIgEKkASWlZ3P3eMuau28mgjk15+Py+dG5eP9xliRRLQSBSjgL5BUz6ZiP/mLGWGIMHRvXm0mM7EBOji8ESuRQEIuVk+ZYM7nz3e5Zv2cOpPVrwwDl9aNO4brjLEjkkBYFIiLJyAzw2cy0vf7WBpvVq88wlAxjRtxWmW0KlilAQiIRg9ppU7nlvOVvS9zNmUHvuPLMHjeI0PIRULQoCkTLYkZnD/R+u5N9LU+jSoj5vXTeYYzo2DXdZImWiIBA5DO7OmwuSeGj6avbn5nPLad247uTO1K4RG+7SRMpMQSBSSj/s2Mvd7y7juw1pDOrUlIfO7UuXFrolVKo+BYHIIbhDSvp+hj8+lzo1Y3j4vL5cNLCdbgmVaqPMQWBm3YGpRWZ1Bv7k7o8XaXMyhS+13xCc9a6731/WfYpUtgUb01i2JZ2s3HyG9W7Jn87uRYsGdcJdlki5KnMQuPsaoB+AmcUCW4D3imk6191HlnU/IuGwe18uD3+8mqkLkhgVBz1aNeS3lwwId1kiFaK8Tg0NBX5w903ltD2RsCgocN5emMxfPl5FZnaAa3/RmRbbUonVaSCpxspr9KvRwJQSlg02s6Vm9rGZ9S6n/YmUu9Xb9nDRC9/yu3e+p0uL+nx084ncNaKnQkCqvZCPCMysFvBL4K5iFi8COrj7XjMbAbwPdC1hO+OB8QAJCQmhliVSavtyAjz52Tpe+moDDevU4JHzj+SCo9vqYrBEjfI4NTQcWOTu2w9c4O57inyebmbPmlm8u+8spu0EYAJA69atvRzqEjkod2fGyu38edoKUjKyuXhgO+4c3oMm9WqFuzSRSlUeQTCGEk4LmVkrYLu7u5kNovBU1K5y2KdISJLSsrhv2go+W51Kj1YNeHJMfwbqyWCJUiEFgZnFAacD1xaZdx2Auz8PXABcb2YBYD8w2t31176ETW6ggBfnruepz9cRY8YfRvRk7JCO1NTLYiSKhRQE7p4FNDtg3vNFPj8NPB3KPkTKy7z1u7jn/eUkpu7ljN4tuffs3rTWMNEierJYqr8dmTn85eNVvLtoC22b1GXi2IGc2qNluMsSiRgKAqm2AvkFTP52E4/NXEt2IJ9fn3wEN53albq1NECcSFEKAqmW5q3fxb0frGDN9kx+0a05953dS+8MFimBgkCqldxAATdPWcy0pSm0aVyXFy47mmG9WuptYSIHoSCQaiE3UEBK+n62pO9nZmAbNw/tyvUnHaHTQCKloCCQKm/uuh3cO20F3TKyaBJXi1nXnET7ZnHhLkukylAQSJW1JX0/D364ko+Xb6NDszh6tGpA47haCgGRw6QgkConJ5DPi3PW8/TsRABuH9aNq0/szBuva/BbkbJQEEiV8vnq7dz/75Vs3JXF8D6tuGdkL9rooTCRkCgIpEpITM3kgQ9X8eXaHXRuXo9XrxrEiV2bh7sskWpBQSARLSMrj8c/W8ur326ibq1Y7jmrJ5cP7kitGhobSKS8KAgkIgXyC3hjfhL/mLGGjP15jB7UnttO70az+rXDXZpItaMgkIjzTeJO7v9wJau3ZXJc56b8aWRverVuGO6yRKotBYFEjM27svi/6Sv5dMV22japy3OXDuDMPq30VLBIBVMQSNjtzQnw7OxEXpq7gRqx9t/bQevU1FPBIpVBQSBh9fbCZP76yWp2ZOZwXv82/O7MHrRqVCfcZYlEFQWBhEVmdoBNu/bx9uql9GvXmAmXHU3/9k3CXZZIVFIQSKXauHMfD3+8GlIyqFUjhkcvOopz+rUhJkbXAUTCJeQgMLONQCaQDwTcfeAByw14AhgBZAFj3X1RqPuVqiU9K5cnP0vk1XkbqRkbw6VN4khoXIfzBrQNd2kiUa+8jghOcfedJSwbDnQNfh0LPBf8V6JAbqCAyd9u5KnPE8nMzuOige249fRuTH93W7hLE5Ggyjg1NAqY7O4OzDOzxmaW4O5bK2HfEibuzifLt/HwJ6vZtCuLE7vGc/eInvRM0PMAIpGmPILAgRlm5sAL7j7hgOVtgKQi08nBef8TBGY2HhgPkJCQUA5lSbgsSUrn/z5ayfyNu+nWsj6Txh3Dyd1bhLssESlBeQTBEHdPMbMWwEwzW+3uc4osL+4qoP9sRmGATABo3br1z5ZL5EvencUjn6xh2tIU4uvX4qFz+3LRwLbUiNW4QCKRLOQgcPeU4L+pZvYeMAgoGgTJQLsi022BlFD3K5FjT3Yez87+gYlfb8CAG0/pwnUnH0H92ropTaQqCOkn1czqATHunhn8PAy4/4Bm04AbzewNCi8SZ+j6QPWQE8jn9XmbeXp2Imn7cjlvQBvuOKM7CY30fgCRqiTUP9laAu8Fx4KpAfzL3T8xs+sA3P15YDqFt44mUnj76LgQ9ylhVlDgTFuawt9nrCF5936OP6IZdw3vSd+2jcJdmoiUQUhB4O7rgaOKmf98kc8O3BDKfiRypO/PY+RTX7Fy6x56JTRk8pV9ObFrvAaGE6nCdBJXSuX75HRWbd1Dxv48Muvn8cTofpx9ZGs9ESxSDSgI5KA27tzH32as4aPvt3J23Xw6NqvHkzecrDeEiVQjCgIp1o7MHJ78bB1T/rOZWjViuHloVxolpRIbYwoBkWpGQSD/Y29OgAlz1vPS3PXkBgoYM6g9Nw3tQosGdZg06ZtwlyciFUBBIABk5+Xz2rxNPPfFD+zal8tZfRO4/YzudIqvF+7SRKSCKQiiXF5+AW8uSOKpzxLZtiebE7rEc/sZ3enXrnG4SxORSqIgiFL5Bc60pVt4bOY6NqdlMaB9Yx69+CiOPyI+3KWJSCVTEEQZd+fTFdv4x4y1rEvdS6+EhkwcO5BTurfQswAiUUpBEEW+WJPKP2asZdmWDI5oXo9nLhnA8D6t9CyASJRTEESBzOwASWlZ3PfP+bRtUpe/XXAk5/Zvo1FBRQRQEFRr3yen8/cZa4lLyaBWbAwPjOrNxce013MAIvI/FATV0PfJ6Twxax2frU6lcVxNLmoaR6tGdbhscMdwlyYiEUhBUI0sS87gic/WMmtVKo3q1uT2Yd244viOvPOGXv8gIiVTEFQDy7dk8Pisdcxatf1/AqBBnZrhLk1EqgAFQRVWNAAa1qnBbad344ohHWmoABCRw6AgqIKWb8ngic/WMXOlAkBEQqcgqEJWpGTwxKx1zAgGwK2nd2OsAkBEQlTmIDCzdsBkoBVQAExw9ycOaHMy8AGwITjrXXc/8J3Gcgh7cwJs2b2f+578SgEgIuUulCOCAHCbuy8yswbAQjOb6e4rD2g3191HhrCfqOTuzFufxjOzE6m/JYMasTHcclphADSqqwAQkfJT5iBw963A1uDnTDNbBbQBDgwCOQzuzhdrdvD07EQWbtpNfP3anN8sjhYN6nDVaV3DXZ6IVEPlco3AzDoC/YHvilk82MyWAinA7e6+ojz2Wd0UFBQOBvf07ERWpOyhTeO63D+qNxcNbMcbryeHuzwRqcZCDgIzqw+8A/zW3fccsHgR0MHd95rZCOB9oNg/a81sPDAeICEhIdSyqoxAfgHTlqbw7Bc/kJi6l07x9XjkgiM5p18bDQUhIpUipCAws5oUhsDr7v7ugcuLBoO7TzezZ80s3t13FtN2AjABoHXr1h5KXVVBTiCfdxZu4bkvE0lK20+PVg14ckx/zuqbQKxGAxWRShTKXUMGvAyscvdHS2jTCtju7m5mg4AYYFdZ91kd5Bc4L3z5AxO/3sD2PTkc1bYRfxrZm6E9Wmg4aBEJi1COCIYAlwHLzGxJcN7dQHsAd38euAC43swCwH5gtLtX+7/2i7N9Tzab07JI3ZPNh2tXM6RLM/5+4VGc0CVeL4QRkbAK5a6hr4CD/gZz96eBp8u6j+ogMTWTCXPW897iLZxWYz9N69Xm31efQN+2jcJdmogIoCeLK8yCjWk8/+V6Zq3aTp2aMYwZ1J42qduoXTNGISAiEUVBUI4KCpxZq7bzwpz1LNy0m8ZxNfnN0K5cPrgDzerXZtKkBeEuUUTkZxQE5SA7L5/3F29hwtz1rN+xj7ZN6vLnX/bmwoFtiaulLhaRyKbfUiHIyy/g0RlreO27zaTty6VXQkOeHNOfEX1a6X3AIlJlKAjKYEVKBj/s2Muuvbl8nJjI0B4tuPKETgzu3Ex3AIlIlaMgKKX8AuezVduZ+PUG5q1P46w6ubRoUJvPx55Mp/h64S5PRKTMFASHsC8nwFsLkvjnNxvZtCuLNo3rcveIHsQkplIjxhQCIlLlKQhKkJSWxeRvN1QD/qUAAAa3SURBVPLG/CQyswMc3aEJvzujB2f0bkmN2BgmrZ8b7hJFRMqFguAA6Vl5XDlpPrPXpBJjxoi+CVx1Qif6tWsc7tJERCqEggBIz8rlrQXJrEhKJzsvn2W1M7jp1K5cMqg9rRrVCXd5IiIVKqqDYPmWDCZ/u5EPlqSQEyjgooYxtGsSx5PXnaohoEUkakRdEOQE8pm+bCuTv93E4s3p1K0Zy3kD2nL54A589+kOAIWAiESVqAmC/Xn57MjMYfBfPidtXy6d4+vxp5G9OP/otv99B3Bxr1cTEanuqnUQZOfl88nybfzrP5tpuiUdM2NQl6Zccmx7TugSr/H/RUSopkGwZlsmU/6zmfcWbyFjfx4dmsXRr2kczRvU5prLjg53eSIiEaXaBEGBO7v25nLes1+zaHM6tWJjOKNPK8Yc047jOjdj8uRN4S5RRCQiVekgKChwvtuQxruLksnatJv8AiejUR73nNWT8wa0pWm9WuEuUUQk4lXJIMjOy+fRGWt4d/EWknfvp37tGpxXvzbNG9TmuWtP0sBvIiKHIaT7JM3sTDNbY2aJZnZnMctrm9nU4PLvzKxjWfe1JzuPKf/ZzIqUPSxJSuep2Yl0iq/HE6P7Mf8Pp9G5eT0a1KmhEBAROUxlPiIws1jgGeB0IBmYb2bT3H1lkWZXAbvdvYuZjQb+Clxc2n1k5QbYtTeXXftyGPjgLHIDBZxbv4D2TeP4dtxQPfUrIlIOQjkiGAQkuvt6d88F3gBGHdBmFPBK8PPbwFAr5Z/s61L3cvQDs1iXmsne7ACXDGrPBzcM4ai2jWnduK5CQESknIRyjaANkFRkOhk4tqQ27h4wswygGbDzUBvfsz+P8wa0oeW2JBrUrcm4X/YGYHEIBYuIyM+Zu5dtRbMLgTPc/erg9GXAIHe/qUibFcE2ycHpH4JtdhWzvfHA+OBkH2B5mQqrfuIpRXBGAfXDT9QXP1Ff/KS7uzcoy4qhHBEkA+2KTLcFUkpok2xmNYBGQFpxG3P3CcAEADNb4O4DQ6it2lBfFFI//ER98RP1xU/MbEFZ1w3lGsF8oKuZdTKzWsBoYNoBbaYBVwQ/XwB87mU9BBERkQpR5iOC4Dn/G4FPgVhgoruvMLP7gQXuPg14GXjVzBIpPBIYXR5Fi4hI+QnpgTJ3nw5MP2Den4p8zgYuLMOmJ4RSVzWjviikfviJ+uIn6ouflLkvynyxWEREqge9gUVEJMqFLQgqc3iKSFeKvrjVzFaa2fdm9pmZdQhHnZXhUH1RpN0FZuZmVm3vGClNX5jZRcHvjRVm9q/KrrGylOJnpL2ZzTazxcGfkxHhqLMymNlEM0s1s2JvsbdCTwb76nszG3DIjbp7pX9ReHH5B6AzUAtYCvQ6oM2vgeeDn0cDU8NRa4T0xSlAXPDz9dHcF8F2DYA5wDxgYLjrDuP3RVcKn7FsEpxuEe66w9gXE4Drg597ARvDXXcF9scvgAHA8hKWjwA+Bgw4DvjuUNsM1xFBhQ5PUcUcsi/cfba7ZwUn51H4zEZ1VJrvC4AHgEeA7MosrpKVpi+uAZ5x990A7p5ayTVWltL0hQMNg58b8fNnmqoNd59DCc9jBY0CJnuheUBjM0s42DbDFQTFDU/RpqQ27h4AfhyeoropTV8UdRWFaV8dHbIvzKw/0M7dP6zMwsKgNN8X3YBuZva1mc0zszMrrbrKVZq+uA/4lZklU3gn401Er8P9nRK29xEU95f9gbcvlaZNdVDq/6eZ/QoYCJxUoRWFz0H7wsxigMeAsZVVUBiV5vuiBoWnh06m8Chxrpn1cff0Cq6tspWmL8YAk9z9H2Y2mMLnl/q4e0HFlxdxDvt3Z7iOCA5neAoONTxFFVeavsDMTgP+APzS3XMqqbbKdqi+aEDhOFRfmNlGCs9/TqumF4xL+zPygbvnufsGYA2FwVDdlKYvrgLeBHD3b4E6FI5DFI1K9TulqHAFgYan+Mkh+yJ4OuQFCkOgup4HhkP0hbtnuHu8u3d0944UXi/5pbuXeYyVCFaan5H3KbyRADOLp/BU0fpKrbJylKYvNgNDAcysJ4VBsKNSq4wc04DLg3cPHQdkuPvWg60QllNDruEp/quUffE3oD7wVvB6+WZ3/2XYiq4gpeyLqFDKvvgUGGZmK4F84A4vZmTfqq6UfXEb8KKZ3ULhaZCx1fQPR8xsCoWnA+OD10TuBWoCuPvzFF4jGQEkAlnAuENus5r2lYiIlJKeLBYRiXIKAhGRKKcgEBGJcgoCEZEopyAQEYlyCgIRkSinIBARiXIKAhGRKPf/HSDIIy/EqT4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "pc_j = ProductionChain_j()\n",
    "p_star = compute_prices(pc_j)\n",
    "\n",
    "transaction_stages = compute_stages(pc_j, p_star)\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(pc_j.grid, p_star(pc_j.grid))\n",
    "ax.set_xlim(0.0, 1.0)\n",
    "ax.set_ylim(0.0)\n",
    "for s in transaction_stages:\n",
    "    ax.axvline(x=s, c=\"0.5\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAELCAYAAAAoUKpTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXxV1bn/8c9DCJCAjGEMYZ5EEIUAWucJAStoxYpD1VZLa7Xe2tZb7b1tqdXeqx20tvqr1IHigAMOpUKlTnVEIJExjDEEMhAghDEEMj2/P87Ri2kgAZJ9zsn5vl8vXpxz9to5z870zdpr7bXN3RERETmSZpEuQEREop/CQkRE6qSwEBGROiksRESkTgoLERGpk8JCRETqFFhYmNl4M1tnZtlmdlct2882s0/NrNLMptTYdoOZbQj/uyGomkVEJMSCuM7CzBKA9cBFQD6wBLja3Vcf0qYP0Bb4MTDX3eeEX+8IZADpgAOZwCh333m490tJSfE+ffo0xqGIiDRZmZmZxe7eubZtzQOqYQyQ7e45AGb2PDAZ+CIs3D03vK26xr4XA2+6e0l4+5vAeGD24d6sT58+ZGRkNGT9IiJNnpltOty2oE5DpQJ5hzzPD7/W2PuKiEgDCCosrJbX6nv+q177mtk0M8sws4zt27cfVXEiInJkQYVFPpB2yPOeQGFD7uvuM9w93d3TO3eu9ZSbiIgco6DCYgkw0Mz6mlkLYCowt577LgDGmVkHM+sAjAu/JiIiAQkkLNy9EriN0C/5NcCL7p5lZveY2SQAMxttZvnAlcBjZpYV3rcE+BWhwFkC3PP5YLeIiAQjkKmzQUtPT/cjzYbas2cP27Zto6KiIsCqokNiYiJdunShbdu2kS5FRKKMmWW6e3pt24KaOhs19uzZw9atW0lNTSUpKQmz2sbPmyZ3p6ysjIKCAgAFhojUW9wt97Ft2zZSU1NJTk6Oq6AAMDOSk5NJTU1l27ZtkS5HRBrQ9r0HeXphLvNXbmmUjx93PYuKigqSkpIiXUZEJSUlxeUpOJGmZse+g7yRVcS8FVv4JGcH1Q6TRvRg4vDuDf5ecRcWQNz1KGqK9+MXiWU7S8tZkFXE6yu2sDBnB1XVTr+U1tx23gAuObkHg7ud0CjvG5dhISISS3KLS5m3cgvzVmxhTdEe3KF3p2S+e04/LhnegxO7n9DofwQqLEREotDmHftDAbGykFUFewAY1bsDd1w4iPOHdOGkHm0DPUugsBARiRL5O/czP9yDWJ6/G4BT0trz35ecyITh3UltH7nxVoVFE+HumBnTp09n+vTpXzwXkeiWV7Kff6zawryVRSzP2wXAyT3bcfeEIUwc3p20jskRrjBEYdFEPPLII7Rs2ZLS0lLuvvtuxo0bx3nnnRfpskSkFiWl5cxdVsCrywq/CIjhqe34yfghTBzejd6dWke4wn+nsIgxxcXF3HHHHSQnJ/PYY4998fptt93GAw88wMMPP8zbb7/NmWeeGcEqRaSmiqpq/rVuOy9n5vP22q1UVDlDu7flrglDmDisO706RUcP4nAUFjHmZz/7GWlpaTz88MNfCotHH32U9u3bc/vttzNv3jwqKys599xzI1eoiFBV7XyYXcwrn+bz3vrt7NpfQUqbFtxweh+uGNWTE7vHzioKCosYUllZyUsvvcSrr77KnDlzvrTtlltu+bcxCxGJjPVb9/JyZj6vLi1g296DtEtK5MITuzJxeDfOHtSZxITYWzxDYRFDFi5cyKBBg8jMzOSMM8740rbPB7OnT5/+peciEowNW/fy6tIC3l6zjXVb99K8mXHu4C5MGZXKeUO60LJ5QqRLPC4KixjywQcfcNpppzF37lzuuOOOSJcjEve27TnA3OWFvLasgFUFe0hoZozt25Gfjx7KpFN6kNKmZaRLbDAKixiSkZHBwIED2bBhAxMmTIh0OSJxqfRgJf9cXcQrnxbwUXYx1R6a6vqLS4fy1ZN70PmEphMQh1JYAL/8exarC/cE+p5De7TlF5eedFT75ObmsmbNGn74wx/SvLm+dCJBqayq5sPsYl5bWsCCrK2UVVTRs0MSt543gMmnpDKgS5tIl9jo9BsnhmzdupVmzZrxne98J9KliDR5FVXVvLduO/NXbuH9Ddsp3ldOu6RELh+ZyuWnppLeu0NcjQ0qLOCo/8KPlObNm3PfffeRnBzd87FFYtnaoj3MycjntWUFFO8rp0NyImcN7MwlJ3fn3MGdY36g+lgpLGLE66+/zubNm2nZsiWrV69m6dKlXHvttZEuS6RJyN+5n9eWFvCPVUVkFe4hMcG4YEhXpozqyTmDY3Oqa0NTWMSAqqoqZs6cyTPPPMPPf/5zunfvzhNPPBHpskRiWll5FW9kbeGljHw+/mwHAKf2as/0S4cy6ZRUOrZuEeEKo4vCIgYkJCR8cRGeehMix87dydy0kzmZ+by+Ygv7DlaS1jGJOy4cxNdGpkbNon3RSGEhIk3etr0HeDmzgBcz8thYXEpyiwQmDu/OlaN6MrpPR5o1i5+B6mOlsBCRJml/eSVvrCpi/soi3l23japqZ0yfjtx63gAmDOtG65b69Xc09NkSkSbD3VmSu5M5mXnMW7GF0vIqurZtyc1n9eWq9DT6dW7610M0FoWFiMS8gl1lvJKZz5xP89m0Yz+tWyRwycnduTI9Le6uh2gsCgsRiUll5VUsyCpiTmY+H31WjDuc3q8Tt58/kAnDu5HcQr/eGlJcfjbj/ZajWr5cYtXByioW5ZQwb8UW5q/cwt6DlfTskMR/XDCQK0b21GymRhR3YZGYmEhZWVlcXwVdVlZGYmJipMsQqbeNxaU8t2gTL2Xms2t/Ba1bJHDxsG5cOSqNsX01mykIcRcWXbp0oaCggNTUVJKSkuKqh+HulJWVUVBQQNeuXSNdjsgRlR6s5M3VW3kpM4+PsnfQvJkx7qSuXDGyJ1/pn0JSi/hcdiNS4i4s2rYN3cawsLCQioqKCFcTvMTERLp27frF50Ekmrg7Kwt2M3txHn9fXsi+g5Wktk/ix+MG8fX0NLq0bRXpEuNW3IUFhAJDvyxFosfusgr+tqyA2YvzWLNlD60Sm3HJ8B5cNTo0m0mnmSIvLsNCRCLv86U3nlu8mXkrtnCwspqTerTlV5cNY9KIHrRL0rhaNAksLMxsPPAHIAF43N3/t8b2lsAsYBSwA7jK3XPNLBF4HBgZrneWu/9PUHWLSMPatb+clz8tYPbizWRv20ebls25YlRPrh7di+E920W6PDmMQMLCzBKAR4CLgHxgiZnNdffVhzS7Cdjp7gPMbCpwP3AVcCXQ0t2Hm1kysNrMZrt7bhC1i8jxc3cWbyxh9uLNzF9VRHllNSPS2nP/FcP56sk9tPRGDAjqKzQGyHb3HAAzex6YDBwaFpOB6eHHc4A/WWiqkgOtzaw5kASUA8HeA1VEjklJaTkvZ+Yze8lmcraXckLL5kwdncbU0b0Y2kPjhrEkqLBIBfIOeZ4PjD1cG3evNLPdQCdCwTEZ2AIkA3e4e0mjVywix8TdWZizg9mL81iwqojyqmpG9mrPb6aczCUnd9eV1TEqqK9abVMZal5GfLg2Y4AqoAfQAfjAzN76vJfyxc5m04BpAL169TrugkXk6GzZXcbcZYU8vyS0DHjbVs25Zmwvrh7Ti8HdToh0eXKcggqLfCDtkOc9gcLDtMkPn3JqB5QA1wBvuHsFsM3MPgLSgS+FhbvPAGYApKenaz0LkQCUV1bzRlYRsz7OJWPTTgDG9OnI988fwMTh3WmVqAvnmoqgwmIJMNDM+gIFwFRCIXCoucANwEJgCvCOu7uZbQbON7NnCJ2GOg14KKC6RaQWW/cc4NlFm5m9eDPb9x6kd6dkfjxuEOOHdWNAF/UimqJAwiI8BnEbsIDQ1Nkn3T3LzO4BMtx9LvAE8LSZZRPqUUwN7/4I8BSwitCpqqfcfUUQdYvI//l8RtOsTzaxYFURVe6cN7gL3zi9N+cM7KwL55o4a4orkKanp3tGRkakyxBpEvaXV/La0kJmLcxlbdFe2rZqzlWj07jutN707tQ60uVJAzKzTHdPr22bpiWISK02Fpfy9MJNvJSZx94DlQzt3pb7rxjOpBGpWsQvDiksROQLVdXOv9ZtY9bCTby3fjvNmxkTh3fn+tN7M0p3nItrCgsRYdf+cl7MyOPpTzaRV1JG17YtuePCQVw9Riu9SojCQiSOrSrYzayFufxtWSEHK6sZ07cjd40/kXEndSUxoVmky5MoorAQiTPlldX8Y9UWZi3cROamnSQlJnDFqJ5cf3pvhnTTEhxSO4WFSJwo2n2A5xZt4rnFeRTvO0ifTsn87KtDmTKqp5YDlzopLESasKpq583VRfz1400szi2h2p3zB3fh+q/04awBKbo2QupNYSHSBB2oqOLlT/N5/IONbCwupVfHZG4+qy/XjulNr07JkS5PYpDCQqQJ2bW/nGc+2cTMj3Mp3lfOiJ7tePTakVx8UjcS1IuQ46CwEGkCCnaV8cQHG3l+yWb2l1dx7uDOfOfs/pzWr6OujZAGobAQiVHuznvrt39xAZ0Bk0b04Ntn9+PE7prVJA1LYSESYyqqqnl9RSGPvZfD2qK9dG3bkpvP7Mv1X+lDavukSJcnTZTCQiRGlB6s5PkleTzxQQ6Fuw8wsEsbfnvlCCaN6EGL5rqAThqXwkIkym3fe5CZH2/k6YWb2HOgkjF9O3Lv5cM4d1AXTX2VwCgsRKKQu7NoYwmzF2/mH6uKqKiqZvxJ3Zh2dj9O7dUh0uVJHFJYiESRiqpq5q/cwuMfbGRlwW7aJydy5aie3HxWP/qm6N4REjkKC5EoULCrjGc/2cQrnxZQtOcA/Tq35teXD+drI1N1H2uJCgoLkQgq2n2AP76zgRcz8qisdk7v14lff03jERJ9FBYiEbCuaC/PfLKJFzPyqHbn6+lp3HreAHpo6qtEKYWFSIA2bN3Lg2+tZ/7KIhKaGZedksoPLhxIWket1yTRTWEh0sg+n9n0/OLNzF1eSHKL5tx+/gC+dWZf2ie3iHR5IvWisBBpRIs3lvC7f65j0cYSWiU249tn9eM75/SnY2uFhMQWhYVIA6uudt5dt42/LtzE++u30/mElky/dChfH51Gcgv9yEls0neuSANxdxZkFfHgmxtYt3Uv7ZMTuXvCEK4/vQ9JLTT9VWKbwkLkOB2srOKfWVuZ8X4OKwt2069zax666hTGD+umaySkyVBYiByj/eWVzF6cx1/ez6FozwF6dkjiN1NO5vJTU2meoIX9pGlRWIgcpd1lFcz6OJenPs6lpLSc0/p15H+uGM4Z/VO0+qs0WQoLkXoqK69i5se5/Pm9z9hdVsH5Q7rwvXP7k96nY6RLE2l0CguROhysrGL2os088q/P2L73IOcN7syPxg1mWGq7SJcmEhiFhchhVFZV88qnBTz01noKdx9gbN+OPHLNSMb0VU9C4o/CQqSGqmrn9RWF/OGtDeQUlzIirT0PTBnBGQM6YabF/SQ+KSxEwtydN1dv5Xf/XM+6rXsZ0u0EHvvGKMYN7aqQkLgXWFiY2XjgD0AC8Li7/2+N7S2BWcAoYAdwlbvnhredDDwGtAWqgdHufiCo2qVpq6yq5oMNxTz09gaW5+2iX0pr/nj1qVwyvLuWCRcJCyQszCwBeAS4CMgHlpjZXHdffUizm4Cd7j7AzKYC9wNXmVlz4BngG+6+3Mw6ARVB1C1Nm7vzz9VbuW/eGjaX7Ce1fRL/+7XhTBnVU9dJiNQQVM9iDJDt7jkAZvY8MBk4NCwmA9PDj+cAf7JQ338csMLdlwO4+46AapYmyt1ZmreL+/+xlkUbSxjYpQ0PXjWCS4b30HUSIocRVFikAnmHPM8Hxh6ujbtXmtluoBMwCHAzWwB0Bp539wcav2Rpatydf63bzv1vrGVt0V46tW7Bry4bxtWj09STEKlDUGFR24lfr2eb5sCZwGhgP/C2mWW6+9tf2tlsGjANoFevXsddsDQtqwp2c9+8NSzM2UHvTsncd/kwLh3Rg7atEiNdmkhMCCos8oG0Q573BAoP0yY/PE7RDigJv/6euxcDmNl8YCTwpbBw9xnADID09PSaQSRxqmBXGb9dsI5XlxbQITmR6ZcO5ZqxvXW6SeQoBRUWS4CBZtYXKACmAtfUaDMXuAFYCEwB3nH3z08//aeZJQPlwDnAgwHVLTFq9/4KHn0vm6c+ysWAW87tzy3n9ldPQuQYBRIW4TGI24AFhKbOPunuWWZ2D5Dh7nOBJ4CnzSybUI9ianjfnWb2e0KB48B8d58XRN0Se0oPVvKnd7OZ+VEuByqruPzUVH40bjCp7ZMiXZpITDP3pnfGJj093TMyMiJdhgQodOOhrdzz9ywKdx/gslN68J1z+nNi97aRLk0kZoTHg9Nr26YruCXmffxZMffNW0NW4R6GdDuBh68+VSvBijQwhYXErKLdB7h33mpeX7GF1PZJ3H/FcL42sieJmgYr0uAUFhJzdu0v58E31zN7SR4G/ODCgXz3nP66halII1JYSMyornZeyMjjgTfWsrusgqtGp3HLOQPo1Sk50qWJNHkKC4l6VdXO35cX8sd3NvDZ9lLG9OnIPZedxJBuGrwWCYrCQqJa5qYS/vu1LNZs2cOJ3dvyp2tCq8FqyXCRYCksJCplb9vL/8xfy9trt9G9XSsevvpUvqolw0UiRmEhUaWq2nn8gxx+9+Z6klskcPsFA/nO2f1o3VLfqiKRpJ9AiRpvrt7Kg2+uZ/WWPVw0tCv3XT6MLie0inRZIoLCQqLA7v0V/PLvWbyytIB+Ka156KpTmHxKD41LiEQRhYVETEVVNTM/yuUvH+Swo7Sc2y8YyG3nDdCKsCJRSGEhEbGuaC8/emkZqwr2MKZvR564YTTDe7aLdFkichgKCwnUztJyHn5nA89+spkTWjXnz9eNZPyw7pEuS0TqoLCQQLg7b6wqYvrfs9ixr5xJI3rwX5ecSKc2LSNdmojUg8JCGl1Gbgl/fCeb99Zv58TubXn8ep1yEok1CgtpNPvLK3ngjXXM/DiX5BYJ/HTiEL51Rl+aa1VYkZijsJBGsXhjCXfOWc6mHfu58St9+M/xg0luoW83kViln15pUGXlVfxmwTqe+ngjPTskMfvbp3F6/06RLktEjtMRw8LMUoAUYIO7VwVTksSqjNwS7pyzgo3FpVx/em9+Mn6IlukQaSLq+kl+D8gCBptZFvCou38IYGYXAoPc/dFGrlGi3IGKKn67YB1PfLSR1PZJPPftsXylf0qkyxKRBlTXSON8oDfwJPAP4Otm9jczG+bubwFXNXaBEt0yN+1k4h8+4PEPN3Lt2F688YOzFRQiTdARexbufqeZDQIuA64FTgYMWGpmhcATjV+iRKMDFVX8/s31PP5BDt3bJfHszWM5Y4BCQqSpqvOEsruvBx4I/8PMWgNtgdHAlY1anUSlTzfv5McvLSdneylXj+nFTycO4YRWiZEuS0Qa0VGPPrp7KVAKzDWz+Q1fkkSrAxVVPPjWev7yfg7d2rbi6ZvGcNbAzpEuS0QCcFxTVdy9sqEKkejl7ry1Zhv3v7GW7G37mDo6jf+65ET1JkTiiOY1yhHtPVDBj19azoKsrfTulMxfvzWGcwapNyESbxQWcljL8nbxwxeXsWnHfn46cQjfPKMviVqqQyQuKSzk31RWVfPAgnXMeD+HlDYtePbmsZzWT1dhi8QzhYV8yY59B7ntuaUszNnB1NFp/PjiwaRoGXGRuKewkC8sz9vFLc9ksqO0nN9eOYIpo3pGuiQRiRIKC6G62pm9ZDO/nLuazie05OVbvsKwVN1vQkT+j8IizuWV7OfW5z5lRf5uzhqYwsNTT6VD6xaRLktEokxgU1vMbLyZrTOzbDO7q5btLc3shfD2RWbWp8b2Xma2z8x+HFTNTd3ba7by1T9+yMbiUn512TBmfnOMgkJEahVIz8LMEoBHgIuAfGCJmc1199WHNLsJ2OnuA8xsKnA/X16o8EFCixnKcaqsqub3b67n0X99xtDubfl/142kd6fWkS5LRKJYUKehxgDZ7p4DYGbPA5OBQ8NiMjA9/HgO8CczM3d3M7sMyCG0zIgch617DvCD55exMGcHV49J4xeXnkSrxIRIlyUiUS6osEgF8g55ng+MPVwbd680s91AJzMrA35CqFeiU1DHyN2Z+XEuv12wjip3zXYSkaMSVFhYLa95Pdv8EnjQ3feZ1dYkvLPZNGAaQK9evY6xzKbJ3bl33hqe+HAj5w3uzC8uPYk+KTrtJCL1F1RY5ANphzzvCRQepk2+mTUH2gElhHogU8zsAaA9UG1mB9z9T4fu7O4zgBkA6enpNYMobpWUljNtVgYZm3Zy41f68POvDqVZs8OHrohIbYIKiyXAQDPrCxQAU4FrarSZC9wALASmAO+4uwNnfd7AzKYD+2oGhdRu+96DXPf4IjbuKOVnXx3Kt87ow5F6ZyIihxNIWITHIG4DFgAJwJPunmVm9wAZ7j6X0F33njazbEI9iqlB1NZUvbt2Gz96aTll5VU8deNo3cVORI6Lhf54b1rS09M9IyMj0mVEzN+XF3LHC8vok9Ka3145glPS2ke6JBGJAWaW6e7ptW3TFdxNSEVVNffNW8PMj3MZ06cjT9yYrhsUiUiDUFg0EbvLKrj12U/5MLuYr52ayr2XDyO5hb68ItIw9NukCcgr2c83Zy5h045SfjPlZK5MT6t7JxGRo6CwiHGZm0qYNiuTympn1rfGcnp/3aRIRBqewiJGVVc7v39zPTPez6F7+1Y8eeNo+nduE+myRKSJUljEqN/+cx2P/uszJg7vxr2XDaejVosVkUaksIgx5ZXV/GJuFrMXb+aasb2477JhutBORBqdwiKGuDt3v7KSlz/N55qxvbhn0kkKChEJhMIiRpQerOS/Xl3Ja8sK+cGFA/nBhYMiXZKIxBGFRQyoqnZun72Ut9du47vn9Oc/LhgY6ZJEJM4oLKLcztJy7nplBW+v3cavLhvGN07rHemSRCQOKSyi2IGKKm6elUHmpp3cefFgBYWIRIzCIkodqKjitueWkrlpJ49eO5KJw7tHuiQRiWMKiyhUUlrOt2dl8Onmnfxq8kkKChGJOIVFlMnILeHOOSso2FXGI9eoRyEi0UFhEUWyt+3lmscX0S4pkeduHkt6n46RLklEBFBYRI2CXWV8f/YyklskMO/7Z9KlbatIlyQi8gWFRRQoPVjJVY8tZNveg8z4xigFhYhEHYVFhJWVV/GTl0NjFM/dfJqWGBeRqKSwiCB353vPZvLuuu3cefFgBYWIRC2FRYS4O39+L4d3123nrglD+O45/SNdkojIYSksIuShtzbwh7c3cNHQrkw7q1+kyxEROaJmkS4gHi3K2cEj72YzYVg3/nzdKJo10zLjIhLdFBYB+yRnBzc8tZhenZL59eXDSVBQiEgMUFgE6IMN2/nWzCX07JDMC9NOp4NuhSoiMUJhEZCPsou5/snF9OyQxHM3j6XzCS0jXZKISL1pgDsAFVXV/OGtDXRMbsFrt55Bcgt92kUktqhnEYBvPrWExbkl/GTCEAWFiMQkhUUj+zi7mA+zi7n1vP58PT0t0uWIiBwThUUjWpSzgxtnLqFvSmu+d+6ASJcjInLMFBaNpLraeeitDbRq3ozZ3z6N1i11+klEYpfCopFMezqThTk7uHP8ELq10yqyIhLbAgsLMxtvZuvMLNvM7qple0szeyG8fZGZ9Qm/fpGZZZrZyvD/5wdV87FaklvCW2u2cvOZfblubK9IlyMictwCCQszSwAeASYAQ4GrzWxojWY3ATvdfQDwIHB/+PVi4FJ3Hw7cADwdRM3HasPWvdw0cwmp7ZO4/cKBmOkKbRGJfUH1LMYA2e6e4+7lwPPA5BptJgN/DT+eA1xgZubuS929MPx6FtDKzKLyirZ9Byv57jOZtGjejFk3jaFtq8RIlyQi0iCCCotUIO+Q5/nh12pt4+6VwG6g5g0ergCWuvvBRqrzmLk733xqMRuLS/nj1SPp37lNpEsSEWkwQU3Rqe1cjB9NGzM7idCpqXG1voHZNGAaQK9ewY8TLM3bxZLcnXz3nP66iZGINDlB9SzygUOvSOsJFB6ujZk1B9oBJeHnPYFXgevd/bPa3sDdZ7h7urund+7cuYHLP7KS0nK+/9xSUtq04BbdxEhEmqCgwmIJMNDM+ppZC2AqMLdGm7mEBrABpgDvuLubWXtgHnC3u38UUL1HZd6KQgp2lXHvZcNol6xxChFpegIJi/AYxG3AAmAN8KK7Z5nZPWY2KdzsCaCTmWUDPwQ+n157GzAA+JmZLQv/6xJE3fWxtmgPv56/lhE923HxSd0iXY6ISKMI7LJid58PzK/x2s8PeXwAuLKW/e4F7m30Ao/RUx/mUlZRxUNTT9U0WRFpsnQF93FYkFXECxl5XD0mjb4prSNdjohIo1FYHCN3D92jonUL7ppwYqTLERFpVAqLY/T7N9ezesse7rhoEO2SNKgtIk2bwuIYZBXu5o/vZHPOoM5cM0ZrP4lI06ewOErV1c5PX1lJh+REHp56KgnNNKgtIk2fwuIoZRXuYXn+bv5z/BBdUyEicUNhcZQee/8zWiQ046KhXSNdiohIYBQWR2HrngO8vmIL3zqzLyltonLhWxGRRqGwOApPfZQLwNdG1lwwV0SkaVNY1NOBiiqe+mgjk0b0YFDXEyJdjohIoBQW9fTsos0crKzm6+lpdTcWEWliFBb14O7MeP8zzhyQwhkDdK8KEYk/Cot6eHfdNrbuOcikU3posUARiUsKizpUVTvT565mQJc2TBrRI9LliIhEhMKiDm+u3srmkv388KJBtEpMiHQ5IiIRobCow8uf5tO9XSvG6SI8EYljCosjKK+sZvHGEs4YkELzBH2qRCR+6TfgEXz8WTG7yyqYMEy3SxWR+KawOIKFn+0goZlxen9NlxWR+KawOILXV2zh7IEpJLcI7FblIiJRSWFxGJt2lFKwq4yzB3WOdCkiIhGnsDiM+SuLALQUuYgICovDemftVgZ2aUPPDsmRLkVEJOIUFrV4Y1URS3J3cv6QLpEuRUQkKigsarGyYBcAd1w0KMKViKGnBJoAAAbqSURBVIhEB4VFLTYWl9IvpbWW9xARCVNY1CJneyl9U1pHugwRkaihsKihutrJ3aGwEBE5lMKihqI9BzhQUU3fzgoLEZHPKSxq2FhcCqCehYjIIRQWNeSEw6JfSpsIVyIiEj0UFjVs3F5KUmICXdu2jHQpIiJRI7CwMLPxZrbOzLLN7K5atrc0sxfC2xeZWZ9Dtt0dfn2dmV3cmHVuLN5H35TWute2iMghAgkLM0sAHgEmAEOBq81saI1mNwE73X0A8CBwf3jfocBU4CRgPPBo+OM1ipziUg1ui4jUEFTPYgyQ7e457l4OPA9MrtFmMvDX8OM5wAUW+vN+MvC8ux90941AdvjjNbjyymrySvbTT4PbIiJfElRYpAJ5hzzPD79Waxt3rwR2A53quW+D2Fyyn2rXTCgRkZqCCovaBgC8nm3qsy9mNs3MMswsY/v27cdQYsjE4d04sXvbY95fRKQpCios8oG0Q573BAoP18bMmgPtgJJ67ou7z3D3dHdP79z52G5YNKBLGx69dpTCQkSkhqDCYgkw0Mz6mlkLQgPWc2u0mQvcEH48BXjH3T38+tTwbKm+wEBgcUB1i4gIEMjNpd290sxuAxYACcCT7p5lZvcAGe4+F3gCeNrMsgn1KKaG980ysxeB1UAlcKu7VwVRt4iIhFjoj/emJT093TMyMiJdhohITDGzTHdPr22bruAWEZE6KSxERKROCgsREamTwkJEROqksBARkTo1ydlQZrYd2HQcHyIFKG6gcmJBvB0v6JjjhY756PR291qvam6SYXG8zCzjcNPHmqJ4O17QMccLHXPD0WkoERGpk8JCRETqpLCo3YxIFxCweDte0DHHCx1zA9GYhYiI1Ek9CxERqVPchoWZjTezdWaWbWZ31bK9pZm9EN6+yMz6BF9lw6rHMf/QzFab2Qoze9vMekeizoZU1zEf0m6KmbmZxfzMmfocs5l9Pfy1zjKz54KusaHV43u7l5m9a2ZLw9/fEyNRZ0MxsyfNbJuZrTrMdjOzh8OfjxVmNvK439Td4+4foWXSPwP6AS2A5cDQGm2+B/w5/Hgq8EKk6w7gmM8DksOPb4mHYw63OwF4H/gESI903QF8nQcCS4EO4eddIl13AMc8A7gl/HgokBvpuo/zmM8GRgKrDrN9IvAPQncaPQ1YdLzvGa89izFAtrvnuHs58DwwuUabycBfw4/nABeYWW23eI0VdR6zu7/r7vvDTz8hdFfCWFafrzPAr4AHgANBFtdI6nPM3wYecfedAO6+LeAaG1p9jtmBz2+B2Y5a7rYZS9z9fUL3/TmcycAsD/kEaG9m3Y/nPeM1LFKBvEOe54dfq7WNu1cCu4FOgVTXOOpzzIe6idBfJrGszmM2s1OBNHd/PcjCGlF9vs6DgEFm9pGZfWJm4wOrrnHU55inA9eZWT4wH/h+MKVFzNH+vNcpkDvlRaHaegg1p4XVp00sqffxmNl1QDpwTqNW1PiOeMxm1gx4ELgxqIICUJ+vc3NCp6LOJdR7/MDMhrn7rkaurbHU55ivBma6++/M7HRCd+Uc5u7VjV9eRDT476947VnkA2mHPO/Jv3dLv2hjZs0JdV2P1O2LdvU5ZszsQuC/gEnufjCg2hpLXcd8AjAM+JeZ5RI6tzs3xge56/u9/Td3r3D3jcA6QuERq+pzzDcBLwK4+0KgFaE1lJqqev28H414DYslwEAz62tmLQgNYM+t0WYucEP48RTgHQ+PHMWoOo85fErmMUJBEevnsaGOY3b33e6e4u593L0PoXGaSe4ey/fkrc/39muEJjNgZimETkvlBFplw6rPMW8GLgAwsxMJhcX2QKsM1lzg+vCsqNOA3e6+5Xg+YFyehnL3SjO7DVhAaCbFk+6eZWb3ABnuPhd4glBXNZtQj2Jq5Co+fvU85t8AbYCXwmP5m919UsSKPk71POYmpZ7HvAAYZ2argSrgTnffEbmqj089j/lHwF/M7A5Cp2NujOU//sxsNqHTiCnhcZhfAIkA7v5nQuMyE4FsYD/wzeN+zxj+fImISEDi9TSUiIgcBYWFiIjUSWEhIiJ1UliIiEidFBYiIlInhYWIiNRJYSESIDMrMLNTIl2HyNHSdRYiAQlfLb0FaNMEllKROKOehUgAzGwAoVVAmwE7zGxHeM0xkZignoVIQMzsVuBsd78q0rWIHC31LESCMwJYFukiRI6FwkIkOKcQuuWnSMzRaSiRAIRvtLQPGODuMX1LT4lP6lmIBCMp/E8/cxKT9I0rEgB3LwX+DKwO339AJKboNJSIiNRJPQsREamTwkJEROqksBARkTopLEREpE4KCxERqZPCQkRE6qSwEBGROiksRESkTgoLERGp0/8H8SccRaqrW/0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "\n",
    "ell_star = np.empty(pc_j.n)\n",
    "for i, s in enumerate(pc_j.grid):\n",
    "    t, e = optimal_choices(pc_j, p_star, s)\n",
    "    ell_star[i] = e #johns ell*\n",
    "\n",
    "\n",
    "    \n",
    "    \n",
    "fig, ax = plt.subplots()\n",
    "\n",
    "ax.plot(pc_j.grid, ell_star, label=\"$\\ell^*$\")\n",
    "\n",
    "\n",
    "ax.set_ylabel(r'$\\ell$', fontsize=12)\n",
    "ax.set_xlabel(r'$t$', fontsize=12)\n",
    "\n",
    "ax.legend(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAELCAYAAAAoUKpTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3deXhU5dnH8e+dlYRddllkV4KoQNhlE0VQBFFUkLor1sprq13cqqVa69K3taJYpYK4VEWpaEQtLrigIhBQkE2JiKwSw76FZGbu94+MvEiDCZDMZDK/z3VxMWfL3CfbL895nnMec3dERER+SkK0CxARkYpPYSEiIiVSWIiISIkUFiIiUiKFhYiIlEhhISIiJYpYWJjZIDP70sxyzOyWYrb3MbOFZhYwsxEHbbvMzFaG/10WqZpFRKSIReI+CzNLBL4CzgDWAfOBUe6+7IB9mgM1gN8AWe4+Lbz+GCAbyAQcWAB0dveth3q/unXrevPmzcvjVEREKq0FCxbkuXu94rYlRaiGrkCOu68CMLMXgGHA/rBw99XhbaGDjj0TeNvdt4S3vw0MAp4/1Js1b96c7OzssqxfRKTSM7NvD7UtUpehGgNrD1heF15X3seKiEgZiFRYWDHrSnv9q1THmtkYM8s2s+zvv//+sIoTEZGfFqmwWAc0PWC5CbChLI9194nununumfXqFXvJTUREjlCkwmI+0MbMWphZCjASyCrlsTOBgWZW28xqAwPD60REJEIiEhbuHgDGUvRLfjnworsvNbO7zGwogJl1MbN1wAXA42a2NHzsFuBuigJnPnDXD53dIiISGREZOhtpmZmZ/lOjoXbs2EFubi6FhYURrKpiSE5Opn79+tSoUSPapYhIBWNmC9w9s7htkRo6W2Hs2LGDTZs20bhxY9LS0jArrv+8cnJ39u7dy/r16wEUGCJSanH3uI/c3FwaN25Menp6XAUFgJmRnp5O48aNyc3NjXY5IlKWduUS+HQiOxdOK5cPH3cti8LCQtLS0qJdRlSlpaXF5SU4kUpndx7Bpa+yY8FL1Nw0lyRCrKxxOp06jSj52MMUd2EBxF2L4mDxfv4iMW3PFoLLstiR/RI1vptDIkG2hBox1Yazo+UQevTsXS5vG5dhISISUzZ/TWjpK+z5bBrpW5eTiLM91IAXbSjbWpxNpy69ubxtPaokJ5ZbCQoLEZGKaOtqQkums/ezaVTdsoQEYEWoLZ9wAXuaD+CULn257IT65RoQB1JYiIhUFNvWEFryCns/f4mqeYtJAL4KtWYml7Cz5Vn07NSRa06oT1pKZALiQAqLSsLdMTPGjRvHuHHj9i+LSAW39Vt82avkL/o3abmfkwDkhFryHx/N9hZn0b1zJ/7nhPpUTY3ur2uFRSUxYcIEUlNT2b17N7feeisDBw6kf//+0S5LRIqzZwt8MY19C58jddNnGLAy1II3Q6PIazaIHpmZ/CKjAdWrJEe70v0UFjEmLy+PG2+8kfT0dB5//PH968eOHcsDDzzA+PHjeffddzn11FOjWKWI/JdgIeS8w77sZ0jKmUmiB8gJHcdrwZGsaTiQHpmdubpDI+pUS412pcVSWMSYO+64g6ZNmzJ+/PgfhcWjjz5KrVq1uOGGG3j99dcJBAL069cveoWKCISCsOo9Chc+R2jlu6QWbmOn12B6cCALaw/ipM6n8rOTG9Gkdnq0Ky2RwiKGBAIBXnrpJaZPn860aT++S/O66677rz4LEYmS3OUEFv6LwOcvUCX/e3Z7Vd4JdWZ+2qnU73g2QzoexzUNq0e7ysOisIghc+bMoW3btixYsIBevXr9aNsPndnjxo370bKIREjuCkKLprJ3yQyqbv8K90Rmh07hreTLqHHSEM7u1Jzzm9aK2Z9NhUUMmT17Nt27dycrK4sbb7wx2uWIyM7v8C9eYu+C50nfvBQngc+D7ZidcDn5JwxnQGZ77m1Zh6TE2H8Mn8IihmRnZ9OmTRtWrlzJ4MGDo12OSHzatwtWzGB39nOkrZ1NAiFWhlryml/K9pZDOS3zRH4VwZvlIkVhAfzxtaUs27Ajou+ZcWwN/nBO+8M6ZvXq1SxfvpybbrqJpCR96UQiJhiAb95nT/ZzJH/1OsmhfLZ6XZ4MDmVNkyF0yezODSc2pEYFGupa1vQbJ4Zs2rSJhIQErr322v3rZsyYwcyZM3n44YejWJlIJfTDUNdFLxPKeYe0gi0UeFWmB3uxpO4g2mYO4MKTG1O/epVoVxoRCgs47L/woyUpKYl77rmH9PT/H2a3ePFiTj755ChWJVLJbFpKMDySKXXfZnZ7NT4KdWBBtX406HwOZ3dszug6VaNdZcQpLGLEjBkzWLNmDampqSxbtozPPvuM0aNHs2jRIlq2bEnPnj1Zs2YNWVlZdOrUKdrlisSWbWvxRS+Qv3g6aZuXEiKRWcFOvJV8DfU7ns2wzs05p1H1mB3JVBYUFjEgGAwyZcoUnn32We68804aNWrEpEmTgKKWxYABA7j33nt59NFHmT59usJCpDQK9sDy19iX/Qwpaz/CcFaEWjPDL2dX23MZ1LU9f2ldt1KMZCoLCosYkJiYuP8mvNGjR+9fn5+fT2FhIWPGjAGKQqVOnTpRqVEkJrjD2rkEFj6LL3mZ5MBuNnk9pgXO56tGQ+jTpTM3nNSImmmVt6P6SCksYtiSJUvIzMzcv/zFF19w4YUXRrEikQpq5yb88+fInzeFtJ2rKfBUXg924/30M2jV+QzO69yM5nXjrx/icCgsYtiiRYvo0KHD/uXFixdz9913R7EikQqkYDcsf409n79MlW/eIYEgi0MnkMV1kDGMc7q25eHmx5CQEL/9EIdDYRHDfuivgKL5LHJzc2nQoEGUqxKJIndYM4eC7GewZa+QHNzDDq/N08Gz+LrxufTs3p3b2zckPUW/+g6XPmMx7KGHHtr/2sxYtWpVFKsRiaJtawl9/jz7sp8hbdcaCr0KM4LdmVN9IG27DmR4pyY0qpkW7SpjmsJCRGJTwR5Y8Tq75z1F+rqPSMD5PJjB60ljqXLSMM7p0pYLm9SM6+GuZUlhISKxI7APVn9E/qJ/k7D8VVICu9gSqsckP48Nxw2nX/dM7jihPqlJleu5TBWBwkJEKr7NX+PZTxJY+CzJ+7YS8CrMDHUhu/Zg2nYdxOhTmlTYGeYqC4WFiFRMBbthxRvsy36a1DUfEiSRt4KdeSv5NBqdfCbDu7Xh/BibQCiWKSxEpOJwhw2fEVzwFL74JZICu9nsdXgucAGrjzuPwd078kCGLjNFg8JCRKJv7zb44iX2zXuS1LylFJLCjGB33k4dyPFdTueiLsfR9JiKP091ZaawEJHoCD96I5g9BV86naRgPjmh45gauoLtrc5lWI92TGhTT89mqiAiFhZmNgh4CEgEnnD3+w7ango8DXQGNgMXuftqM0sGngA6het92t3vjVTdIlLG9myBxVPJnzuZKlu/Yi9pvBroxfvVBtOxe3/GdmpC/RrxMUdELIlIWJhZIjABOANYB8w3syx3X3bAblcBW929tZmNBO4HLgIuAFLdvYOZpQPLzOx5d18didpFpAyE76wunDeZhOWvkhgqYEWoFS/6GArbncvwbsczqkUdPXqjAotUy6IrkOPuqwDM7AVgGHBgWAwDxoVfTwMesaK7aRyoamZJQBpQAER2DlQROTK7N+OLniN/7pOkbf+afE/j5WBfPq11Dl179OF3HRtTKz0l2lVKKUQqLBoDaw9YXgd0O9Q+7h4ws+1AHYqCYxiwEUgHbnT3LeVesYgcGXdY/RH75k0macVrJHohy0Jt+DfXkdhhOOd1a8ulTWvpzuoYE6mwKO67wku5T1cgCBwL1AZmm9k7P7RS9h9sNgYYA9CsWbOjLlhEDtP29YS+mMa+uZNJ27mafK/Ky8HT+LzeUHr06MNtJx9LtVSNqYlVkfrKrQOaHrDcBNhwiH3WhS851QS2ABcD/3H3QiDXzD4GMoEfhYW7TwQmAmRmZh4cRJWeu2NmjBs3jnHjxu1fFilXwUJYnkXhnMdJXj+XBOCL0PG8mjCWtI7nM6J7a65oWCPaVUoZiFRYzAfamFkLYD0wkqIQOFAWcBkwBxgBzHJ3N7M1wGlm9ixFl6G6A3+PUN0xY8KECaSmprJ7925uvfVWBg4cSP/+/aNdllRWOzbCgikUzp9M8p5cNnp9XgxcwPpGp9Pv1N7c0b4hVZJ141xlEpGwCPdBjAVmUjR0drK7LzWzu4Bsd88CJgHPmFkORS2KkeHDJwBPAksoulT1pLsvjkTdFVFeXh433ngj6enpPP744/vXjx07lgceeIDx48fz7rvvcuqpp0axSqmUwiOagnMnYstfI8EDfBQ8mal2JXVPOYuf9WzBCWpFVFrmXvmu2GRmZnp2dnax25YvX067du0iXFHZue6666hduzbjx49n165d+9c/+uijJCUlsXLlSpKSkjjzzDPp16/fIT9OrH8eJIIKdsMXL1Ew53FS8paxg6pMDfTlw5pDOb1XD4Z3akyNKpqzujIwswXunlncNvU2xZBAIMBLL73E9OnTmTZt2o+2XXfddf/VZyFyVDZ/jc97gsDCZ0gu3ElO6DieCV3D7rbDGdnzeK5uVUf9YnFEYRFD5syZQ9u2bVmwYAG9evX60bYffmjHjRv3o2WRwxIKQs47FM55jORvZhEkkTeDXclKOZuMngO5oVszzTgXpxQWMWT27Nl0796drKwsbrzxxmiXI5XJni3w2bPs+/SfpO5cwxavzXOB81ne+DzO6dWJR9s3JCVJz2iKZwoLgDdvge++iOx7NuwAg+8reb8DZGdn06ZNG1auXMngwYPLqTCJKxsXEfh0IvbFSySG9vF56ASm8iuqdTyXi3u04kZ1WEuYwiKGrF69muXLl3PTTTeRlKQvnRyhQAEsz2LfJ4+RunE+haQyPdCL92oO49Re/fhjp8ZUV4e1HES/ceCw/8KPlk2bNpGQkMC11167f92cOXN4++23ufPOO6NYmcSEHRsIzZ9M4fwnSc3PY0OoAc+GLmFrmwsYcWp7JrZUh7UcmsIihiQlJXHPPfeQnv7/k8D06NGDHj16RLEqqdBCQVjxOoFPHydh7RzwELODp5CV8nOa9zyHa7o1p2FNPQ5cSqawiBEzZsxgzZo1pKamsmzZMj777DNGjx7NiBEjuOWWW8jMLHZotMSrwnxY9DyFH40nedsqvvN6zAiexeIGwzmrbw/+2r4hyZpUSA6DwiIGBINBpkyZwrPPPsudd95Jo0aNmDRpEgBLly4lIyMjyhVKhbFnC2RPovCTx0jOz2NFqAUTQzeQmDGUK3q34edNa0W7QolRCosYkJiYuP8mvNGjR+9fv2fPHoAfXZaSOLVtLaE5EwhlP0VScA8fB0/m6cRf0KbbYG7t1YJja+neCDk6CosYtmTJEtq3bx/tMiRa3OHrdwnMeYyEVbMIOWQFe5BVdQT9evfj4cymVNUjwaWM6Dsphi1evJgOHTpEuwyJtGAhLJ1OwYcPkpK3jM1em+nBwSxqdCHD+nZjUkYDEjU9qZQxhUUM++KLL+jbt2+0y5BI2bcLFj7Nvo8eIXX3elaHGvNE6OcUZpzP5b3bqj9CypXCIobNnj2bW265JdplSHnblUvo08cIzH2ClMLtfB46gacTbqZJ92H8qmdL9UdIRCgsYlB+fj49e/Zk0KBBNGrUKNrlSHkIzx1ROG8StjyLhFAhs4KZZFU9n259BvNA5ybqj5CI0ndbDKpSpQoLFy6MdhlSHoKFsOxVCmaPJyV3EXu8KjOCvZnXaCRn9evDw+3UHyHRobAQqQi2r4PsyRQsfJ6U3RtYG2rEk6GryM+4gEv7tGN0E/VHSHQpLESiacdGmP2/hBY8BaEA84MZ/CvhUpp2H8r1vVpq7gipMBQWItGwaRmePZnQwqfxYICpgX48lzKCs/t15b7ux2maUqlwFBYikZS7An/vz9jyVwmRwCvBXjyfNopzBvZiWmZT0lISo12hSLEUFiKRsPpjgtlTsCXTyCeFJwLnMqvW+Yzu15HnTmmsWeikwovLsHD3uH5uv7tHu4T48e0cgu/eTeKajwmQwpTAYD6odzGXnNaZ69s31MgmiRlxFxbJycns3bs3rh++t3fvXpKTdU283IRCkPM2hXMeJ/mbd9lCLSYUXsqqJudx9YATGdOmblz/sSKxKe7Con79+qxfv57GjRuTlpYWVz+07s7evXtZv349DRo0iHY5lY87rJhB4N17SMpbzm6vxj8Co/i21WiuPq09mc2PiXaFIkcs7sKiRo2iCeg3bNhAYWFhlKuJvOTkZBo0aLD/8yBlILAPVrxOweyHSNn0OWu9EeMDv4B2Q7n6tAzaH1sz2hWKHLW4CwsoCgz9spSjVrAHFj5FYPbfSdr9Hblel0eC15Jwykhu6Hc8LepWjXaFImUmLsNC5Kjkb4d5Ewl88g+S8jeTHWrHk345TTPP4pf92upGOqmUFBYipVW4F+ZNJPjh30jct40Pg6cwJfGXnHLqIO7t1YJjqqZEu0KRcqOwEClJYB8seIrAB/9L0p5NzA6dzCPcTI/eA3j41JbUTNfIMqn8FBYihxIMwOIXCMz6M0k717MwdALj/ee07zmYiX1aqSUhcUVhIXKwUAiWvkxg1r0kbc1haagVD/pttOwyhAf7t6Ze9dRoVygScQoLkR+4w5dvEnjnLpLylvO1N+Wh4E0c03k4957WRh3XEtciFhZmNgh4CEgEnnD3+w7ango8DXQGNgMXufvq8LaTgMeBGkAI6OLu+ZGqXSq5YABWvUfhrHtJ3riAdd6QvwfHUuXkC7h1QFuaHhO/d/uL/CAiYWFmicAE4AxgHTDfzLLcfdkBu10FbHX31mY2ErgfuMjMkoBngUvcfZGZ1QHi7246KXvu8OUbBP9zO4nbviHP6zA+cDWFHUbxq9Pb0Vz3SYjsF6mWRVcgx91XAZjZC8Aw4MCwGAaMC7+eBjxiRc/iGAgsdvdFAO6+OUI1S2XlDusXEHjrDySt+YhvaczDBdcRyjiX/xnYntb1q0e7QpEKJ1Jh0RhYe8DyOqDbofZx94CZbQfqAG0BN7OZQD3gBXd/oPxLlkrHHXLeIfj2OBJzl7CD6vyt8Ao2tx3JDQMzaNdId/WLHEqkwqK4p/Ud/JzsQ+2TBJwKdAH2AO+a2QJ3f/dHB5uNAcYANGvW7KgLlkpm4yJCM+8gYfUHbKQBjxZexbaWQ7juzM50aKJnN4mUJFJhsQ5oesByE2DDIfZZF+6nqAlsCa//wN3zAMzsDaAT8KOwcPeJwESAzMxMTdggRbavw2fdDYumspOqPFh4KauaX8gvB55I5+NqR7s6kZgRqbCYD7QxsxbAemAkcPFB+2QBlwFzgBHALHf/4fLT78wsHSgA+gIPRqhuiVV7t8HHfyc051ECwRCTA0N4+5iL+eWQLoxrWy/a1YnEnIiERbgPYiwwk6Khs5PdfamZ3QVku3sWMAl4xsxyKGpRjAwfu9XM/kZR4Djwhru/Hom6JQYV7IbZfyU051ESAnt5JXgqU1JH87MhvXmxUxPNTCdyhKwyTrGZmZnp2dnZ0S5DIumHYbBv/I7EHevICvXiSYbSr88ArunTgvQU3X8qUpJwf3Bmcdv0EySxb/VHhP5zGwnfLWIVTbm94E5aZZ7B46e3pX6NKtGuTqRSUFhI7NqxEX/r99iSaeRaPf5eeDXftzqfu8/qwPENda+ESFlSWEjs2bsV3r+P0PzJBELOPwLn8W6di/ntkFPo3Uad1yLlQWEhsSMUgs+fJfjWH7D8bbwY6MOLaRcwakhfpqvzWqRcKSyk4gs/Mjz4/v0kbv6KhaHj+TM3c1r/ATzbW53XIpGgnzKp2NbOJ/T6TSR8t5ivacYjhWOp1vkiHj+jLfWrq/NaJFIUFlIx5a2Et+6Ar95ks9Xh7oLr2dV6KLec3Z62DdR5LRJpCgupWEJB+PQf+Lt3s8eTeSIwnDdrXMjNIzPpf3z9aFcnErcUFlJxfDUTn/Un7LvFvO+d+X3wai7s14VX+rakSnJitKsTiWsKC4m+vdvgP7fCoudYl9CYvxVcx/Y25/H80BNpVkez1IlUBAoLiZ5gAOY9TvDjh2FXLo8EhvNy1VHcfv7JnJHRgKK5r0SkIlBYSHTkrsCn/xzb+Bmfezv+FLyenr3P4D/925CWoktOIhWNwkIia88W+PAvhOY9wU6vwq0FN7Cj5RD+d1h7WtWrFu3qROQQFBYSOctnEHr9N/iuXF4N9mBilSsZO7wHZ3dopEtOIhWcwkLK35q5+If/i+W8xUqO4+bCP9Kl5wCmnd6Waqn6FhSJBfpJlfJTsAdm3Y1/+g/2kcKDhaNY3HQ09597ip4KKxJjFBZSPtZ8SuiVX5Cw5WueDZ7BxJRLuXFoJ27p2FiXnERikMJCylbhXpj1J3zOBHKtHjcW3E6TjmcyY0gGNdOSo12diByhnwwLM6sL1AVWunswMiVJzFo7j9D060jYksNzwQFMqnIld4zsqsd0iFQCJbUsPgCWAseb2VLgUXf/CMDMTgfauvuj5VyjVHSFe+G9e/A5E/ieOtxUcCuNOw3mlSEZ1Kii1oRIZVBSWLwB9AEmA1uAC83st8Dt7v6Omd0BKCzi2bpsQtN/TsLmlTwfPI1JaVdyx8hu9FNrQqRS+cmwcPffmllb4FxgNHASYMBnZrYBmFT+JUqFVJgP7/8Z/+Rh8jiGXxfcwrGdzmb6kHZqTYhUQiV2cLv7V8AD4X+YWVWgBtAFuKBcq5OKaf2CotZE3ldMDfbnibQruWNUD/q21fzXIpXVYY+GcvfdwG4gy8zeKPuSpMIK7IP378M//jt51OY3BTdzbOchvHy2WhMild1RDZ1190BZFSIVmDt8NZPQO38g4fsVvBTsxz/TruKOUT3oo9aESFzQfRby0/bthFevh2Wv8p015LaC39Eo8xz+fZZaEyLxRGEhh7Z+IT7953jeSv4SGMmM9PO4Z1QntSZE4pDCQv5bMACz7oaP/87OhJpcW3Ar9U86gzfOPZHqak2IxCWFhfzY7s3w7yth1ftkWX/uL7yIXw4/lQsym+iZTiJxTGEh/2/DZ/jUnxHckctthdewqO5QplzckTYN9IRYkXinsJCi0U4Ln8Lf+B15XoOr8u+gfZf+vDIkQ1OcigigsJBta+ClK2B9NvM4kd/4r/jdyF6cc/Kx0a5MRCqQhEi9kZkNMrMvzSzHzG4pZnuqmU0Nb59rZs0P2t7MzHaZ2W8iVXOlt/Id/PE+5G9czu8Lr+DeOn/m2RvOUlCIyH+JSFiYWSIwARgMZACjzCzjoN2uAra6e2vgQeD+g7Y/CLxZ3rXGhVAQZt2D/2sEqwpqcebeu0ntMYYXr+vNcXWqRrs6EamAInUZqiuQ4+6rAMzsBWAYsOyAfYYB48KvpwGPmJm5u5vZucAqih4zIkdj5yaYPgZWvc9078f9oav586VdGdCuQbQrE5EKLFJh0RhYe8DyOqDbofZx94CZbQfqmNle4GbgDECXoI7GvH/i74wjUFjAbYVj+Lbpebwy6hQa1UyLdmUiUsFFKiyKG6Dvpdznj8CD7r7rp8b5m9kYYAxAs2bNjrDMSsod3r4TPhnPgqRO/G7faM7u15t7B7QhKTFi3VYiEsMiFRbrgKYHLDcBNhxin3VmlgTUpGjCpW7ACDN7AKgFhMws390fOfBgd58ITATIzMw8OIji154tMPVn8O3HvGhnck/gCh65IpPebfTIDhEpvUiFxXygjZm1ANYDI4GLD9onC7gMmAOMAGa5uwO9f9jBzMYBuw4OCjmE3Xnw9LkEv1/BA8HRvF3jAqZf3oWW9apFuzIRiTERCYtwH8RYYCaQCEx296VmdheQ7e5ZFM2694yZ5VDUohgZidoqrZx38ZfHEMjfyZX5v4ZWpzF9VCdqpuvZTiJy+Kzoj/fKJTMz07Ozs6NdRvQsfQX/99V8l9CQX+y+mpN7nMHvz26n/gkR+UlmtsDdM4vbpju4K5NgoKgj+9MJLE08gUv2/IbfnNuV0d2Oi3ZlIhLjFBaVRf4OmHYF5LzD6/Tm7uA1TLiqFz1b1Y12ZSJSCSgsKoNta+C5iwh9/xW/D1zD3NpDeOGyLjSvq7uxRaRsKCxi3boF+PMj2Ze/hyv3/Y7k1v2ZfnFHTXkqImVKYRGr3OH9e/GPHyKX2ly85w/07Xkqt511gjqyRaTMKSxi1Xt/hg8f4KOkHvx6z+XcNLwnI7vqznURKR8Ki1gTLIQ3b4bsScxIHMAthWOYeFWmOrJFpFwpLGKJO7x+Eyx8mlcTBvDH4DX86+runNy0VrQrE5FKTmERKwr2FAXFoud5ws7nsYSLeW5MV05oWCPalYlIHFBYxIJQCF6+BlbMYApDeDJ5FC9e013PeBKRiFFYVHR7t8Jrv4QVM/hz6DLeqj6cqVd3o0nt9GhXJiJxRGFRkQUK4IXR8O3HPBS8gPdrn8eLV3enfo0q0a5MROKMwqKiKsyHl6+Gbz/mV4GxrGo4mKlXdKV21ZRoVyYicUhhURHt2VLUoljzCXcFLmFDkyH86/JMquuubBGJEoVFRbN2Hrx6PcHN3/DLgv9he6tzeOqSTNJSEqNdmYjEMYVFRbL5a3hqKHsS0rks/1aOyejLE6M6kpqkoBCR6NJDhCqK7etg2hUUkET/HX+k3on9mHBxJwWFiFQIallUBAV7YMoQQtvXcf2+X9GqVRsevOgUPRBQRCoMhUW0Fe4tuo9i6zdcHbyd7+r3YOolndWiEJEKRWERTe4w7Sr48nUetYtYWT2Tf1/ZRaOeRKTCUVhEizvMeQS+fJ2JSRczyc/j31d2o3513XAnIhWPwiJaZv8VZt3N3KRMHi44h+fGdNU0qCJSYSksomHNXPyD+8lO7c5lu8byz8u70KFJzWhXJSJySBpuE2lrPsWfGU5uQj2u3X45D1zYmd5t6kW7KhGRn6SwiKRVH+DPjiAv4RiG7LyNsUO6M/TkY6NdlYhIiRQWkbL6Y3jmXLYm1WXI9ps5v28mV57aItpViYiUisIiEoKF8MF9FCZXp8/W39Pt5BO5edDx0a5KRKTUFBaR8MLF8M2HPBAYSYO6dbn3vIFBt2MAAAwzSURBVA6YWbSrEhEpNY2GKm/ffgIr3yKr6gie2dGfV0d3pmqqPu0iElvUsihPa+fBM8PZltqY2zcP4k/nduD4htWjXZWIyGFTWJSXUAjev5cAiZy942YGdW7LiM5Nol2ViMgRUViUl2mXw9ez+FvwQqrVb8Fdw06MdkUiIkcsYmFhZoPM7EszyzGzW4rZnmpmU8Pb55pZ8/D6M8xsgZl9Ef7/tEjVfMTWZcOyV3kjfShTgmcyYXQnzXQnIjEtImFhZonABGAwkAGMMrOMg3a7Ctjq7q2BB4H7w+vzgHPcvQNwGfBMJGo+Ynk58K8RbE9pwC1bhnDveSfRun61aFclInJUItWy6ArkuPsqdy8AXgCGHbTPMOCp8OtpwAAzM3f/zN03hNcvBaqYWWpEqj5cBbvhxUsoCDrn7/w1Z3fNYNgpjaNdlYjIUYtUWDQG1h6wvC68rth93D0AbAfqHLTP+cBn7r6vnOo8cu7w3EV47nJuKBxLcsMM/nDOwY0nEZHYFKkB/8XdgeaHs4+Ztafo0tTAYt/AbAwwBqBZs2ZHVuXR2LgIVs9mVq3zeS/vRN68uCNVktVPISKVQ6RaFuuApgcsNwE2HGofM0sCagJbwstNgOnApe7+dXFv4O4T3T3T3TPr1YvwU1z3boUXL6UwtTa/+e4Mru/fmpb11E8hIpVHpMJiPtDGzFqYWQowEsg6aJ8sijqwAUYAs9zdzawW8Dpwq7t/HKF6D8+yLNj2LfdxJTXrNGRMn5bRrkhEpExFJCzCfRBjgZnAcuBFd19qZneZ2dDwbpOAOmaWA9wE/DC8dizQGrjDzD4P/6sfibpL5fsvYeZtfF+1LZO3d+SPw07U5ScRqXQi9pAid38DeOOgdXce8DofuKCY4/4E/KncCzxScx+Hgl1cuvdaBp14LH3baiIjEal89ES7o/HVTMiexMfVzuTbHU2YNESjn0SkctLjPo6UO3zwAIUpNflF3vncMKANx9ZKi3ZVIiLlQmFxpD78C6zP5lE/n3r1G3JlL816JyKVly5DHYncFfD+vXxbswvjN/Xn2ZEnkpKk3BWRyku/4Q6XO7z2S4Ip1bkw7xrOOaUpPVodfKO5iEjlorA4XLnLYe2nvHbMFWyzGtx2VrtoVyQiUu4UFofr44fwhCT+uqZoMqP6NapEuyIRkXKnsDgcuzfD4hdY0GAE60K1uLq37tQWkfigsDgcc/8BwF82duTMjIa0qFs1ygWJiESGwqK0AgXw6T/4pv7pzM1vypi+alWISPxQWJTWwqegYBcTtnajS/PadGpWO9oViYhEjMKitD55mLw6mUzbmcGYPq2iXY2ISEQpLErj61mw7Vte2NuVVvWqMeCEivPQWxGRSFBYlCQUgjdvZm/15jyypQvX9G5JQkJxk/qJiFReCouS5LwDeV8xJWUU1arV4NyOB08dLiJS+SksSrLoeULp9fjbhnZc3K2ZJjYSkbiksPgpwQCsns3aWl0p9CTObN8g2hWJiESFwuKnrJkDu7/nzWAXjq1ZhYxGNaJdkYhIVCgsfso3H+AYkzc0Y0C7BpipY1tE4pPC4qcsnc7WBj3ILazCgHYaLisi8UthcSjb18PmHOYmnkJ6SiLdW2rOChGJXwqLQ1n2KgCTv8+gT5t6GgUlInFNYXEoX73JvhrHMX/nMboEJSJxT2FRnK/egm8+ZFnVHphBfz3eQ0TiXFK0C6iQNn4OwH37zqdj06rUrZYa5YJERKJLLYvibM4hWL0JczcUcnqGbsQTEVFYFGdzDrkpTQA4vZ3CQkREYXEwd8jLYXlBfZoek0ab+tWiXZGISNQpLA62Ow/2bWfOtlqcrru2RUQAhcV/25wDwMpgQ12CEhEJU1gcLBwWm5Kb0KX5MVEuRkSkYtDQ2YN4Xg4BkmjTNoOUJGWpiAhEsGVhZoPM7EszyzGzW4rZnmpmU8Pb55pZ8wO23Rpe/6WZnVmedW5fv5xvQ/U5LePY8nwbEZGYEpGwMLNEYAIwGMgARplZxkG7XQVsdffWwIPA/eFjM4CRQHtgEPBo+OOVi0DuSr6hEf2Or1debyEiEnMi1bLoCuS4+yp3LwBeAIYdtM8w4Knw62nAACsaijQMeMHd97n7N0BO+OOVvVCQGnvXsqd6C2qlp5TLW4iIxKJIhUVjYO0By+vC64rdx90DwHagTimPLRMb1+aQQiG1m7Yrjw8vIhKzIhUWxd2s4KXcpzTHYmZjzCzbzLK///77IygR9hUUMK9qf1qd1OuIjhcRqawiNRpqHdD0gOUmwIZD7LPOzJKAmsCWUh6Lu08EJgJkZmb+V5iURvM2HWj+21eO5FARkUotUi2L+UAbM2thZikUdVhnHbRPFnBZ+PUIYJa7e3j9yPBoqRZAG2BehOoWEREi1LJw94CZjQVmAonAZHdfamZ3AdnungVMAp4xsxyKWhQjw8cuNbMXgWVAALje3YORqFtERIpY0R/vlUtmZqZnZ2dHuwwRkZhiZgvcPbO4bbpFWURESqSwEBGREiksRESkRAoLEREpkcJCRERKVClHQ5nZ98C3R/Eh6gJ5ZVROLIi38wWdc7zQOR+e49y92KeoVsqwOFpmln2o4WOVUbydL+ic44XOuezoMpSIiJRIYSEiIiVSWBRvYrQLiLB4O1/QOccLnXMZUZ+FiIiUSC0LEREpUdyGhZkNMrMvzSzHzG4pZnuqmU0Nb59rZs0jX2XZKsU532Rmy8xssZm9a2bHRaPOslTSOR+w3wgzczOL+ZEzpTlnM7sw/LVeambPRbrGslaK7+1mZvaemX0W/v4+Kxp1lhUzm2xmuWa25BDbzczGhz8fi82s01G/qbvH3T+KHpP+NdASSAEWARkH7fML4LHw65HA1GjXHYFz7g+kh19fFw/nHN6vOvAh8CmQGe26I/B1bgN8BtQOL9ePdt0ROOeJwHXh1xnA6mjXfZTn3AfoBCw5xPazgDcpmmm0OzD3aN8zXlsWXYEcd1/l7gXAC8Cwg/YZBjwVfj0NGGBmxU3xGitKPGd3f8/d94QXP6VoVsJYVpqvM8DdwANAfiSLKyelOedrgAnuvhXA3XMjXGNZK805O1Aj/Lomxcy2GUvc/UOK5v05lGHA017kU6CWmTU6mveM17BoDKw9YHldeF2x+7h7ANgO1IlIdeWjNOd8oKso+ssklpV4zmbWEWjq7jMiWVg5Ks3XuS3Q1sw+NrNPzWxQxKorH6U553HAz8xsHfAG8D+RKS1qDvfnvUSRmoO7oimuhXDwsLDS7BNLSn0+ZvYzIBPoW64Vlb+fPGczSwAeBC6PVEERUJqvcxJFl6L6UdR6nG1mJ7r7tnKurbyU5pxHAVPc/a9m1oOiWTlPdPdQ+ZcXFWX++yteWxbrgKYHLDfhv5ul+/cxsySKmq4/1eyr6EpzzpjZ6cDtwFB33xeh2spLSedcHTgReN/MVlN0bTcrxju5S/u9/aq7F7r7N8CXFIVHrCrNOV8FvAjg7nOAKhQ9Q6myKtXP++GI17CYD7QxsxZmlkJRB3bWQftkAZeFX48AZnm45yhGlXjO4Usyj1MUFLF+HRtKOGd33+7udd29ubs3p6ifZqi7x/KcvKX53n6FosEMmFldii5LrYpolWWrNOe8BhgAYGbtKAqL7yNaZWRlAZeGR0V1B7a7+8aj+YBxeRnK3QNmNhaYSdFIisnuvtTM7gKy3T0LmERRUzWHohbFyOhVfPRKec5/AaoBL4X78te4+9CoFX2USnnOlUopz3kmMNDMlgFB4Lfuvjl6VR+dUp7zr4F/mtmNFF2OuTyW//gzs+cpuoxYN9wP8wcgGcDdH6OoX+YsIAfYA1xx1O8Zw58vERGJkHi9DCUiIodBYSEiIiVSWIiISIkUFiIiUiKFhYiIlEhhISIiJVJYiESQma03s1OiXYfI4dJ9FiIREr5beiNQrRI8SkXijFoWIhFgZq0pegpoArDZzDaHnzkmEhPUshCJEDO7Hujj7hdFuxaRw6WWhUjknAx8Hu0iRI6EwkIkck6haMpPkZijy1AiERCeaGkX0NrdY3pKT4lPalmIREZa+J9+5iQm6RtXJALcfTfwGLAsPP+ASEzRZSgRESmRWhYiIlIihYWIiJRIYSEiIiVSWIiISIkUFiIiUiKFhYiIlEhhISIiJVJYiIhIiRQWIiJSov8Dag3LQvG3Cg0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "l_star = σ_star(pc.t_grid) # my ell*\n",
    "\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(pc.t_grid, l_star, label=r'$\\ell^*_h$')\n",
    "ax.plot(pc_j.grid, ell_star, label=\"$\\ell^*_j$\")\n",
    "\n",
    "\n",
    "ax.set_ylabel(r'$\\ell$', fontsize=12)\n",
    "ax.set_xlabel(r'$t$', fontsize=12)\n",
    "\n",
    "ax.legend(fontsize=14)\n",
    "plt.show()"
   ]
  },
  {
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}