von_neumann_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A von Neumann growth model\n",
    "##### and a generalization\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This notebook uses the class  `neumann` to calculate key objects of a linear growth model of John von Neumann (1937) that was generalized by Kemeny, Moregenstern and Thompson (1956). Objects of interest are the maximal expansion rate (`alpha`), the interest factor (`beta`), and the optimal intensities (`x`) and prices (`p`). \n",
    "\n",
    "\n",
    "In addition to watching how the towering mind of John von Neumann formulated an equilibrium model of   price and quantity vectors in balanced growth, this notebook  shows how fruitfully to  employ the following important tools:\n",
    "\n",
    "  * a zero-sum two player game\n",
    "  \n",
    "  * linear programming\n",
    "  \n",
    "  * the Perron-Frobenius theorem\n",
    "\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "np.set_printoptions(precision = 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The next few lines of code fetch the file `neumann.py` that does most of the calculations."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "from quantecon import fetch_nb_dependencies\n",
    "fetch_nb_dependencies([\"neumann.py\"])\n",
    "from neumann import neumann"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "-------\n",
    "**Notation:** We use the following notation. $\\mathbf{0}$ denotes a vector of zeros. We call an $n$-vector \n",
    "- positive or $x\\gg \\mathbf{0}$ if $x_i>0$ for all $i=1,2,\\dots,n$ \n",
    "- nonnegative or $x\\geq \\mathbf{0}$ if $x_i\\geq 0$ for all $i=1,2,\\dots,n$\n",
    "- semipositive or $x > \\mathbf{0}$ if $x\\geq \\mathbf{0}$ and $x\\neq \\mathbf{0}$  \n",
    "\n",
    "For two comformable vectors $x$ and $y$, $x\\gg y$,  $x\\geq y$ and $x> y$ mean $x-y\\gg \\mathbf{0}$,  $x-y \\geq \\mathbf{0}$, and $x-y > \\mathbf{0}$. \n",
    "\n",
    "By default, all vectors are column vectors, $x^{T}$ denotes the transpose of $x$ (i.e. a row vector). Let $\\iota_n$ denote a column vector composed of $n$ ones, i.e. $\\iota_n = (1,1,\\dots,1)^T$. Let $e^i$ denote the vector (of arbitrary size) containing zeros except for the $i$th position where it is one.  \n",
    "\n",
    "We denote matrices by capital letters. For an arbitrary matrix $A$, $a_{i,j}$ represents the entry in its $i$th row and $j$th column.  $a_{\\cdot j}$ and $a_{i\\cdot}$ denote the $j$th column and $i$th row of  $A$, respectively.\n",
    "\n",
    "------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Model ingredients and assumptions:\n",
    "\n",
    "\n",
    "A pair $(A,B)$ of $m\\times n$ nonnegative matrices defines an economy.\n",
    "\n",
    "  *  $m$ is the number of *activities* (or sectors)\n",
    "  \n",
    "  *  $n$ is the number of *goods* (produced and/or used in the economy)\n",
    "  \n",
    "  *  $A$ is called the *input matrix*;  $a_{i,j}$ denotes the amount of good $j$ consumed by activity $i$. \n",
    "  \n",
    "  * $B$ is called the *output matrix*; $b_{i,j}$ represents the amount of good $j$ produced by activity $i$. \n",
    "\n",
    "Two key assumptions restrict economy $(A,B)$: \n",
    "- **Assumption I:** (every good which is consumed is also produced)\n",
    "    $$ b_{.,j} > \\mathbf{0}\\hspace{5mm}\\forall j=1,2,\\dots,n$$ \n",
    "- **Assumption II:** (no free lunch)\n",
    "    $$ a_{i,.} > \\mathbf{0}\\hspace{5mm}\\forall i=1,2,\\dots,m$$\n",
    "\n",
    "A semipositive $m$-vector $x$ denotes the levels at which activities are operated (_intensity vector_). \n",
    "\n",
    "Therefore, \n",
    "\n",
    "   *  vector $x^TA$ gives the total amount of *goods used in production*\n",
    "   \n",
    "   *  vector $x^TB$ gives *total outputs*\n",
    "   \n",
    "An  economy $(A,B)$ is said to be _productive_, if there exists a nonnegative intensity vector $x \\geq 0$ such that $x^T B > x^TA$. \n",
    "\n",
    "The semipositive $n$-vector $p$ contains prices assigned to the $n$ goods. \n",
    "\n",
    "The $p$ vector implies *cost* and *revenue* vectors\n",
    "\n",
    "   * the vector $Ap$ tells *costs* of the vector of activities\n",
    "   \n",
    "   * the vector  $Bp$ tells *revenues* from the vector of activities\n",
    "-------------------------\n",
    "\n",
    "A property of an input-output pair $(A,B)$ called *irreducibility* (or indecomposability) determines whether an economy can be decomposed into multiple ''sub-economies'': \n",
    "\n",
    "**Definition:** Given an economy $(A,B)$, the set of goods $S\\subset \\{1,2,\\dots,n\\}$ is called an _independent subset_ if it is possible to produce every good in $S$ without consuming any good outside $S$. Formally, the set $S$ is independent if $\\exists T\\subset \\{1,2,\\dots,m\\}$ (subset of activities) such that $a_{i,j}=0$, $\\forall i\\in T$ and $j\\in S^c$ and for all $j\\in S$, $\\exists i\\in T$, s.t. $b_{i,j}>0$. The economy is **irreducible** if there are no proper independent subsets.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We study two examples, both coming from Chapter 9.6 of Gale (1960)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# (1) Irreducible (A, B) example: alpha_0 = beta_0\n",
    "A1 = np.array([[0, 1, 0, 0], \n",
    "               [1, 0, 0, 1], \n",
    "               [0, 0, 1, 0]])\n",
    "\n",
    "B1 = np.array([[1, 0, 0, 0], \n",
    "               [0, 0, 2, 0], \n",
    "               [0, 1, 0, 1]])\n",
    "\n",
    "# (2) Reducible (A, B) example: beta_0 < alpha_0 \n",
    "A2 = np.array([[0, 1, 0, 0, 0, 0], \n",
    "               [1, 0, 1, 0, 0, 0], \n",
    "               [0, 0, 0, 1, 0, 0], \n",
    "               [0, 0, 1, 0, 0, 1], \n",
    "               [0, 0, 0, 0, 1, 0]])\n",
    "\n",
    "B2 = np.array([[1, 0, 0, 1, 0, 0], \n",
    "               [0, 1, 0, 0, 0, 0], \n",
    "               [0, 0, 1, 0, 0, 0], \n",
    "               [0, 0, 0, 0, 2, 0], \n",
    "               [0, 0, 0, 1, 0, 1]])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    " The following code sets up our first Neumann economy or  `neumann` instance  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Generalized von Neumann expanding model:\n",
       "  - number of goods          : 4\n",
       "  - number of activities     : 3\n",
       "\n",
       "Assumptions:\n",
       "  - AI:  every column of B has a positive entry    : True\n",
       "  - AII: every row of A has a positive entry       : True\n"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N1 = neumann(A1, B1)\n",
    "N1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "\n",
       "Generalized von Neumann expanding model:\n",
       "  - number of goods          : 6\n",
       "  - number of activities     : 5\n",
       "\n",
       "Assumptions:\n",
       "  - AI:  every column of B has a positive entry    : True\n",
       "  - AII: every row of A has a positive entry       : True\n"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N2 = neumann(A2, B2)\n",
    "N2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Dynamic interpretation \n",
    "\n",
    "Attach a time index $t$ to the preceding objects, regard aneconomy as a dynamic system, and study  sequences \n",
    "\n",
    "$$\\{(A_t,B_t)\\}_{t\\geq 0}, \\hspace{1cm}\\{x_t\\}_{t\\geq 0},\\hspace{1cm} \\{p_t\\}_{t\\geq 0}$$\n",
    "\n",
    "An interesting special case holds the technology process constant and investigates the dynamics of quantities and prices only. \n",
    "\n",
    "Accordingly, in the rest of this notebook we assume that $(A_t,B_t)=(A,B)$ for all $t\\geq 0$. \n",
    "\n",
    "-------\n",
    "A crucial element of the dynamic iterpretation involves the timing of production.\n",
    "\n",
    "We assume that production (consumption of inputs) takes place in period $t$, while the associated output materializes in period $t+1$, i.e. consumption of $x_{t}^TA$ in period $t$ results in $x^T_{t}B$ amounts of output in period $t+1$. \n",
    "\n",
    "These timing conventions imply the following feasibiliy condition:\n",
    "\n",
    "\\begin{align}\n",
    "x^T_{t}B \\geq x^T_{t+1} A \\hspace{1cm}\\forall t\\geq 1,\n",
    "\\end{align}\n",
    "\n",
    "which asserts that no more goods can be used today than were produced yesterday. \n",
    "\n",
    "Accordingly,  $Ap_t$ tells the costs of production in period $t$ and $Bp_t$ tells revenues in period $t+1$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Balanced growth\n",
    "\n",
    "We follow John von Neumann in studying \"balanced growth\"\n",
    "\n",
    "Let `./` denote elementwise division of one vector by another and let $\\alpha >0$ be a scalar.  \n",
    "\n",
    "Then *balanced growth* is a situation in which \n",
    "\n",
    "$$ x_{t+1}./x_t = \\alpha , \\quad \\forall t \\geq 0$$\n",
    "\n",
    "With balanced growth, the law of motion of $x$ is evidently $ x_{t+1}=\\alpha x_t $ and so we can rewrite the feasibility constraint as\n",
    "\n",
    "$$x^T_{t}B \\geq \\alpha x^T_t A \\hspace{1cm}\\forall t$$\n",
    "\n",
    "In the same spirit, define $\\beta\\in\\mathbb{R}$ as the **interest factor** per unit of time.\n",
    "\n",
    "We assume that it is always possible to earn a gross return equal to the constant interest factor $\\beta$ by  investing \"outside the model\". \n",
    "\n",
    "Under this assumption about outside investment opportunities, a no-arbitrage condition gives rise to the following (no profit) restriction on the price sequence:\n",
    "\n",
    "$$\\beta Ap_{t} \\geq B p_{t} \\hspace{1cm}\\forall t$$\n",
    "\n",
    "This says that production cannot yield a return greater than that offered by the investment opportunity (note that we compare values in period $t+1$). \n",
    "\n",
    "The balanced growth assumption allows us to drop time subscripts and conduct an analysis purely in terms of a time-invariant growth rate $\\alpha$ and interest factor $\\beta$\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Duality\n",
    "\n",
    "The following two problems are connected by a remarkable dual relationship between the technological and valuation characteristics of the economy:\n",
    "\n",
    "**Definition:** The *technological expansion problem* (TEP) for economy $(A,B)$ is to find a semipositive $m$-vector $x>0$ and a number $\\alpha\\in\\mathbb{R}$, s.t.\n",
    "\n",
    "\\begin{align*}\n",
    "    &\\max_{\\alpha} \\hspace{2mm} \\alpha\\\\\n",
    "    &\\text{s.t. }\\hspace{2mm}x^T B \\geq \\alpha x^T A\n",
    "\\end{align*}\n",
    "\n",
    "Theorem 9.3 of David Gale's book assets that if Assumptions I and II are both satisfied, then a maximum value of $\\alpha$ exists and it is positive. It is called the _technological expansion rate_ and is denoted by $\\alpha_0$. The associated intensity vector $x_0$ is the _optimal intensity vector_.\n",
    "\n",
    "--------------\n",
    "\n",
    "**Definition:** The *economical expansion problem* (EEP) for $(A,B)$ is to find a semipositive $n$-vector $p>0$ and a number $\\beta\\in\\mathbb{R}$, such that\n",
    "\n",
    "\\begin{align*}\n",
    "    &\\min_{\\beta} \\hspace{2mm} \\beta\\\\\n",
    "    &\\text{s.t. }\\hspace{2mm}Bp \\leq \\beta Ap\n",
    "\\end{align*}\n",
    "\n",
    "Assumptions I and II imply existence of a minimum value $\\beta_0>0$ called the _economic expansion rate_. The corresponding price vector $p_0$ is the _optimal price vector_."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "---------------\n",
    "\n",
    "Evidently, the criterion functions in  *technological expansion* problem and the *economical expansion problem*  are both linearly homogeneous, so the optimality of $x_0$ and $p_0$ are defined only up to a positive scale factor. \n",
    "\n",
    "For simplicity (and to emphasize a close connection to zero-sum games), in the following, we normalize both vectors  $x_0$ and $p_0$ to have unit length. \n",
    "\n",
    "A standard duality argument (see Lemma 9.4. in (Gale, 1960)) implies that under Assumptions I and II, $\\beta_0\\leq \\alpha_0$\n",
    "\n",
    "But in the other direction, that is $\\beta_0\\geq \\alpha_0$, Assumptions I and II are not sufficient. \n",
    "\n",
    "Nevertheless, von Neumann (1937) proved the following remarkable \"duality-type\" result connecting TEP and EEP.\n",
    "\n",
    "**Theorem 1 (von Neumann):** If the economy $(A,B)$ satisfies Assumptions I and II, then there exists a set $\\left(\\gamma^{*}, x_0, p_0\\right)$, where $\\gamma^{*}\\in[\\beta_0, \\alpha_0]\\subset\\mathbb{R}$, $x_0>0$ is an $m$-vector, $p_0>0$ is an $n$-vector and the following holds true\n",
    "\n",
    "\\begin{align*}\n",
    "x_0^T B &\\geq \\gamma^{*} x_0^T A \\\\\n",
    "Bp_0 &\\leq \\gamma^{*} Ap_0 \\\\\n",
    "x_0^T\\left(B-\\gamma^{*} A\\right)p_0 &= 0\n",
    "\\end{align*}\n",
    "\n",
    "> *Proof (Sketch):* Assumption I and II imply that there exist $(\\alpha_0, x_0)$ and $(\\beta_0, p_0)$ solving the TEP and EEP, repspectively. If $\\gamma^*>\\alpha_0$, then by defintion of $\\alpha_0$, there cannot exist a semipositive $x$ that satisfies $x^T B \\geq \\gamma^{*} x^T A$. Similarly, if $\\gamma^*<\\beta_0$, there is no semipositive $p$ so that $Bp \\leq \\gamma^{*} Ap$. Let $\\gamma^{*}\\in[\\beta_0, \\alpha_0]$, then $x_0^T B \\geq \\alpha_0 x_0^T A \\geq \\gamma^{*} x_0^T A$. Moreover, $Bp_0\\leq \\beta_0 A p_0\\leq \\gamma^* A p_0$. This two inequalities imply $x_0\\left(B - \\gamma^{*} A\\right)p_0 = 0$.\n",
    "\n",
    "Here the constant $\\gamma^{*}$ is both expansion and interest factor (not neccessarily optimal). We have already encountered and discussed the first two inequalities that represent feasibility and no-profit conditions. Moreover, the equality compactly captures the requirements that if any good grows at a rate larger than $\\gamma^{*}$ (i.e., if it is ''oversupplied''), then its price must be zero;  and that if any activity provides negative profit, it must be unused. Therefore, these expressions encode all equilbrium conditions and Theorem I essentially states that under Assumptions I and II there always exists an equilibrium $\\left(\\gamma^{*}, x_0, p_0\\right)$ with balanced growth.\n",
    "\n",
    "Note that Theorem I is silent about uniqueness of the equilibrium. In fact, it does not rule out (trivial) cases with $x_0^TBp_0 = 0$ so that nothing of value is produced. To exclude  such uninteresting cases, Kemeny, Morgenstern and Thomspson (1956) add an extra requirement \n",
    "\n",
    "$$x^T_0 B p_0 > 0$$\n",
    "\n",
    "and call the resulting equilibria *economic solutions*. They show that this extra condition does not affect the existence result, while it significantly reduces the number of (relevant) solutions.  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "----------------------------\n",
    "\n",
    "## Interpretation as a game theoretic problem (two-player zero-sum game)\n",
    "\n",
    "\n",
    "To compute the equilibrium $(\\gamma^{*}, x_0, p_0)$, we follow the algorithm proposed by Hamburger, Thompson and Weil (1967), building on the key insight that the equilibrium (with balanced growth) can be considered as a solution of a particular two-player zero-sum game. First, we introduce some notations.\n",
    "\n",
    "Consider the $m\\times n$ matrix $C$ as a payoff matrix, with the entries representing payoffs from the **minimizing** column player to the **maximizing** row player and assume that the players can use mixed strategies: \n",
    "- row player chooses the $m$-vector $x > \\mathbf{0}$, s.t. $\\iota_m^T x = 1$\n",
    "- column player chooses the $n$-vector $p > \\mathbf{0}$, s.t. $\\iota_n^T p = 1$\n",
    "\n",
    "**Definition:** The $m\\times n$ matrix game $C$ has the *solution* $(x^*, p^*, V(C))$ in mixed strategies, if \n",
    "\n",
    "\\begin{align}\n",
    "(x^*)^T C e^j \\geq V(C)\\quad \\forall j\\in\\{1, \\dots, n\\}\\quad \\quad \\text{and}\\quad\\quad (e^i)^T C p^* \\leq V(C)\\quad \\forall i\\in\\{1, \\dots, m\\}\n",
    "\\end{align}\n",
    "\n",
    "The number $V(C)$ is called the *value* of the game. \n",
    "\n",
    "From the above definition, it is clear that the value $V(C)$ has two alternative interpretations: \n",
    "* by playing the appropriate mixed stategy, the maximizing player can assure himself at least $V(C)$ (no matter what the column player chooses)\n",
    "* by playing the appropriate mixed stategy, the minimizing player can make sure that the maximizing player will not get more than $V(C)$ (irrespective of what is the maximizing player's choice)\n",
    "\n",
    "From the famous theorem of Nash (1951), it follows that there always exists a mixed strategy Nash equilibrium for any *finite* two-player zero-sum game. Moreover, von Neumann's Minmax Theorem (1928) implies that \n",
    "\n",
    "$$V(C) = \\max_x \\min_p \\hspace{2mm} x^T C p = \\min_p \\max_x \\hspace{2mm} x^T C p = (x^*)^T C p^*$$\n",
    "\n",
    "### Connection with Linear Programming (LP) \n",
    "\n",
    "Finding Nash equilibria of a finite two-player zero-sum game can be formulated as a linear programming problem. To see this, we introduce the following notation\n",
    "- For a fixed $x$, let $v$ be the value of the minimization problem:  $v \\equiv \\min_p x^T C p = \\min_j x^T C e^j$ \n",
    "- For a fixed $p$, let $u$ be the value of the maximization problem:  $u \\equiv \\max_x x^T C p = \\max_i (e^i)^T C p$ \n",
    "\n",
    "Then the *max-min problem* (the game from the maximizing player's point of view) can be written as the *primal* LP\n",
    "\n",
    "\\begin{align}\n",
    "V(C) = & \\max \\hspace{2mm} v \\\\\n",
    "\\text{s.t. } \\hspace{2mm} v \\iota_n^T &\\leq x^T C  \\\\\n",
    "x &\\geq \\mathbf{0} \\\\\n",
    "\\iota_n^T x & = 1\n",
    "\\end{align}\n",
    "\n",
    "while the *min-max problem* (the game from the minimizing player's point of view) is the *dual* LP\n",
    "\n",
    "\\begin{align}\n",
    "V(C) = &\\min \\hspace{2mm} u \\\\\n",
    "\\text{s.t. } \\hspace{2mm}u \\iota_m &\\geq Cp  \\\\\n",
    "p &\\geq \\mathbf{0} \\\\\n",
    "\\iota_m^T p & = 1\n",
    "\\end{align}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "---------------------------------------------------------\n",
    "\n",
    "Hamburger, Thompson and Weil (1967) view the input-output pair of the economy as payoff matrices of two-player zero-sum games. Using this interpretation, they restate Assumption I and II as follows \n",
    "\n",
    "$$V(-A) < 0\\quad\\quad \\text{and}\\quad\\quad V(B)>0$$\n",
    "\n",
    "> *Proof (Sketch)*: \n",
    "* $\\Rightarrow$ $V(B)>0$ implies $x_0^T B \\gg \\mathbf{0}$, where $x_0$ is a maximizing vector. Since $B$ is non-negative, this requires that each column of $B$ has at least one positive entry, which is Assumption I. \n",
    "* $\\Leftarrow$ From Assumption I and the fact that $p>\\mathbf{0}$, it follows that $Bp > \\mathbf{0}$. This implies that the maximizing player can always choose $x$ so that $x^TBp>0$, that is it must be the case that $V(B)>0$.\n",
    "\n",
    "In order to (re)state Theorem I in terms of a particular two-player zero-sum game, we define the matrix for $\\gamma\\in\\mathbb{R}$\n",
    "\n",
    "$$M(\\gamma) \\equiv B - \\gamma A$$\n",
    "\n",
    "For fixed $\\gamma$, treating $M(\\gamma)$ as a matrix game, we can calculate the solution of the game.\n",
    "\n",
    "* If $\\gamma > \\alpha_0$, then for all $x>0$, there $\\exists j\\in\\{1, \\dots, n\\}$, s.t. $[x^T M(\\gamma)]_j < 0$ implying that $V(M(\\gamma)) < 0$\n",
    "* If $\\gamma < \\beta_0$, then for all $p>0$, there $\\exists i\\in\\{1, \\dots, m\\}$, s.t. $[M(\\gamma)p]_i > 0$ implying that $V(M(\\gamma)) > 0$\n",
    "* If $\\gamma \\in \\{\\beta_0, \\alpha_0\\}$, then (by Theorem I) the optimal intensity and price vectors $x_0$ and $p_0$ satisfy\n",
    "\n",
    "\\begin{align}\n",
    "x_0^T M(\\gamma) \\geq \\mathbf{0}^T \\quad \\quad \\text{and}\\quad\\quad M(\\gamma) p_0 \\leq \\mathbf{0} \n",
    "\\end{align}\n",
    "\n",
    "That is, $(x_0, p_0, 0)$ is a solution of the game $M(\\gamma)$ so that $V\\left(M(\\beta_0)\\right) = V\\left(M(\\alpha_0)\\right) = 0$.\n",
    "* If $\\beta_0 < \\alpha_0$ and $\\gamma \\in (\\beta_0, \\alpha_0)$, then $V(M(\\gamma)) = 0$. Moreover, if $x'$ is optimal for the maximizing player in $M(\\gamma')$ for $\\gamma'\\in(\\beta_0, \\alpha_0)$ and $p''$ is optimal for the minimizing player in $M(\\gamma'')$ where $\\gamma''\\in(\\beta_0, \\gamma')$, then $(x', p'', 0)$ is a solution for $M(\\gamma)$, $\\forall \\gamma\\in (\\gamma'', \\gamma')$. \n",
    "> *Proof (Sketch):* If $x'$ is optimal for a maximizing player in game $M(\\gamma')$, then $(x')^T M(\\gamma')\\geq \\mathbf{0}^T$ and so for all $\\gamma<\\gamma'$\n",
    "$$(x')^T M(\\gamma) = (x')^T M(\\gamma') + (x')^T(\\gamma' - \\gamma)A \\geq \\mathbf{0}^T $$\n",
    "hence $V(M(\\gamma))\\geq 0$. If $p''$ is optimal for a minimizing player in game $M(\\gamma'')$, then $ M(\\gamma'')p''\\leq \\mathbf{0}$ and so for all $\\gamma''<\\gamma$\n",
    "$$ M(\\gamma)p'' = M(\\gamma'') + (\\gamma'' - \\gamma)Ap'' \\leq \\mathbf{0}$$\n",
    "hence $V(M(\\gamma))\\leq 0$. \n",
    "\n",
    "It is clear from the above argument that $\\beta_0$, $\\alpha_0$ are the minimal and maximal $\\gamma$ for which $V(M(\\gamma))=0$. Moreover, Hamburger et al. (1967) show that the function $\\gamma \\mapsto V(M(\\gamma))$ is continuous and nonincreasing in $\\gamma$. This suggests an algorithm to compute $(\\alpha_0, x_0)$ and $(\\beta_0, p_0)$ for a given input-output pair $(A, B)$. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Algorithm\n",
    "\n",
    "Hamburger, Thompson and Weil (1967) propose a simple bisection algorithm to find the minimal and maximal roots (i.e. $\\beta_0$ and $\\alpha_0$) of the function $\\gamma \\mapsto V(M(\\gamma))$.\n",
    "\n",
    "### Step 1\n",
    "\n",
    "First, notice that we can easily find trivial upper and lower bounds for $\\alpha_0$ and $\\beta_0$. \n",
    " * TEP requires that $x^T(B-\\alpha A)\\geq \\mathbf{0}^T$ and $x > \\mathbf{0}$, so if $\\alpha$ is so large that $\\max_i\\{[(B-\\alpha A)\\iota_n]_i\\} < 0$, then TEP ceases to have a solution. Accordingly, let **`UB`** be the $\\alpha^{*}$ that solves $\\max_i\\{[(B-\\alpha^{*} A)\\iota_n]_i\\} = 0$. \n",
    " * Similar to the upper bound, if $\\beta$ is so low that $\\min_j\\{[\\iota^T_m(B-\\beta A)]_j\\}>0$, then the EEP has no solution and so we can define **`LB`** as the $\\beta^{*}$ that solves $\\min_j\\{[\\iota^T_m(B-\\beta^{*} A)]_j\\}=0$.\n",
    "\n",
    "The `bounds` method calculates these trivial bounds for us "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(1.0, 2.0)"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N1.bounds()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Step 2\n",
    "\n",
    "Compute $\\alpha_0$ and $\\beta_0$\n",
    "\n",
    "* Finding $\\alpha_0$\n",
    "    1. Fix $\\gamma = \\frac{UB + LB}{2}$ and compute the solution of the two-player zero-sum game associated with $M(\\gamma)$. We can use either the primal or the dual LP problem. \n",
    "    2. If $V(M(\\gamma)) \\geq 0$, then set $LB = \\gamma$, otherwise let $UB = \\gamma$\n",
    "    3. Iterate on 1. and 2. until $|UB - LB| < \\epsilon$\n",
    "\n",
    "* Finding $\\beta_0$\n",
    "    1. Fix $\\gamma = \\frac{UB + LB}{2}$ and compute the solution of the two-player zero-sum game associated with $M(\\gamma)$. We can use either the primal or the dual LP problem. \n",
    "    2. If $V(M(\\gamma)) > 0$, then set $LB = \\gamma$, otherwise let $UB = \\gamma$\n",
    "    3. Iterate on 1. and 2. until $|UB - LB| < \\epsilon$\n",
    "\n",
    "\n",
    ">*Existence*: Since $V(M(LB))>0$ and $V(M(UB))<0$ and $V(M(\\cdot))$ is a continuous, nonincreasing function, there is at least one $\\gamma\\in[LB, UB]$, s.t. $V(M(\\gamma))=0$.\n",
    "\n",
    "The `zerosum` method calculates the value and optimal strategies assocaited with a given $\\gamma$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Value of the game with gamma = 2\n",
      "-0.24\n",
      "Intensity vector (from the primal)\n",
      "[ 0.32  0.28  0.4 ]\n",
      "Price vector (from the dual)\n",
      "[ 0.4   0.32  0.28  0.  ]\n"
     ]
    }
   ],
   "source": [
    "gamma = 2\n",
    "\n",
    "print('Value of the game with gamma = {}'.format(gamma))\n",
    "print(N1.zerosum(gamma = gamma)[0])               \n",
    "print('Intensity vector (from the primal)')\n",
    "print(N1.zerosum(gamma = gamma)[1])               \n",
    "print('Price vector (from the dual)')\n",
    "print(N1.zerosum(gamma = gamma, dual = True)[1])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.legend.Legend at 0x7f94bb78a438>"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzcAAAFjCAYAAAAAUfUSAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzs3XmcjeX/x/HXZyxJq63sKtJOlkTFjLTwpSItKER7KWmh\nTXOESjv17ZfKViK0l6J1pFRUWoW+FCJlLUW2uX5/XAfHmMMwZ+a+z5n38/E4j5n7Pve5z2eMc19z\n3dfn+lzmnENERERERCTZpQUdgIiIiIiISCKocyMiIiIiIilBnRsREREREUkJ6tyIiIiIiEhKUOdG\nRERERERSgjo3IiIiIiKSEtS5ERERERGRlKDOjYiIiIiIpAR1bkREdoOZZefhMd/MMqPfF+h11sy6\nm9lcM1tvZisL8r12Ecc5ZtYrl/2ZZrY5iJh2h5mdYGblgzqvmZU0s3vNrOTuPCciIttT50ZEZPc0\nzvFYCkwCTozZ1y56rCvIQMysEjAU+BjIAE4ryPfbhbbADp0b4GmgSWEGYmYlzGyoma0xs5Vm9qSZ\nVYh5voyZvWtmG8xsnJk1Auo655bnOE93M5sUPS6yk/f7r5ltMrNRZtY1Zn/D3M6bG+fcBuAZ4KHd\neU5ERLZXPOgARESSiXNueuy2ma0HljvnZuTY37YQwqmNv0n1rHPu00J4v93mnFsCLCnk99wIXGlm\nZYFyzrmrcjy/ysxGA48A7wNjnXPn5nKe4dGRt7XAEbm9l5k1BSoAnzvnYjs2ewG353bencQ9L9oZ\n+49z7q28PiciItto5EZEpGAdZmZvRkcRfjGzvjkPMLO6ZvZ69I/XtWb2sZmdsrOTmtkI4MPo5gfR\nFLjh0edGmtnPubwmy8w+iNmORF9XK48xvmJmy6MxzjazPjGxdAWqxKbmxb5HLudraWbToudaHT13\n7RzH5Dm+OH4GDs3lvcsCh0Y7CT2BF3ZyjprAu8BhuZynFFAdqAhk5Xh6V+eNZzBw2x48JyIiqHMj\nIlKQDHgZPzpwDvAK0C9H6lJ94BPgQOAy4FxgBfCemdXbybnvBq6Pfn81Ph2uf3TbkXtKXM59W7Z3\nFWMjYBq+o9AT+A8+Rapq9JD+wFvAMral58Wm5m33vmbWEngT+As4H7gKOBaYGk212634duJnoKqZ\nFcux/xa2pXh1jp4/HgfMx3dycrog+nM0YsfOza7Om/ubObcSKGFmh+zOcyIi4iktTUSk4DjgQefc\ns9HtD8ysBdARGBXd9wDwC9DcObcZwMwmAz8AffGdnR1P7NzPZvZjdPPHnOlyCY7xQWA5cKJzbn10\nX1ZMLPPNbBmwIWd6XhwDgHnAf5xz2QBm9hkwF7gJuHk344vnZ/xNvOrR7zGzM4Fpzrk1ZnYEsMI5\ntym3F5vZofiOzXygjJkd4Jz7M/rcscAcoD6+E/txzOtyPa+ZHQUMB44DXgQudc5tNrPiwLXOucHR\nQz8D0vH/L3La2XMiIkWeRm5ERApWzvkR3+P/2N6S1tQM/4cuZlYsOspQDHgv+lzQMe4NnASMjunY\n7DEzKw3UA8Zt6dgAOOd+wY9gpe9OfLvwM77jcWj0vfcBTnPOvRF9vh6+gxJPOr4T9wu+k1Uzeh4D\nGjnnPscXcvjCObcu5nU7nDfaUfo/fGpZG2AT0Cf69NnAuJjDlwJHxolpZ8+JiBR56tyIiBSsnOWZ\n1wOlot+XxXdk+gIbYx4bgB74VLWgYyyDbysWJ+i9yuA7HL/l8txS/L/J7sS3M7/gOyVb5t30wY9C\nbXEwsHonrz/COfe/aIGCX9k27+Y84KXo9+nAlByvy+28HYCznHMvOOeynHOXAdWiz5Vxzi2NOXYF\n/t8pNzt7TkSkyFNamohIcFYD2cDj+BQrS9B5/wVyWxOlHD69bHeswsdYJb9BxZzP4Sfh51SRHTsy\ne8w5t97MlgKHmlkTYL5z7veYQ/bCdyTjniLm+/lATTOrAvztnPszWhHtROCeHK/b4bzOuXtzOf93\nZnYhO87X2QtfoS03O3tORKTI08iNiEhAnHNrgan4tVBmOue+yvnYw1MvAA42s3JbdphZTeKUM95F\njOvw80kujv4xH896YO88nG8t8CVwfjS9a0t8NfDpbx/Ge+0e+hmfxnWRc25kjuf+IPeRIsysOtvP\na5kH1AJaO+feju5rjG9HP97+1fHPm8Mi4Ezn3Lwc+8tFz5GbnT0nIlLkqXMjScm2rf6e87HZzDoF\nHd/uMrMPzWx8As5T3syGmNnn5lesn5+I+KRA3Qg0MLN3zOxCM2tmZuea2QAzyzkikJvcRnsmRL8+\nb2ZnmNlFwKv4amZ74mb8H9WfmdnFZpZhfoHLITHHzALKmtlVZtYwOuE+nr7A4cBEM2tjZh2Bd/Cj\nOg/vYYzx/AycGee8P5P7CBJASyB27aD50X2TYvZlAF9GO2x5PW+sf4EPctlfE/gmzmt29pyEiNqp\nuOc5zcwmmNlCM/vHzL4zs2vNryklkm9KS5Nkthr/R0vOP+7+F0AsYVEFX1r3c3xKzUHBhlMkxCu7\nTF72O+dmmtkJQCZ+svkB+E7IV8CTeXz/7Xf4BR/b46uSvYKvQtYLuD2X4/MS4xdmdjK+/PQQfGrU\nAmBEzPHP4FO0BuLnCi1g2xyV7d7DOTfZzFrjf+Zx+BSuD4E+Oeae5Cm+XZgDZDrncuvoT4vGvVV0\nrZ2++FLWM8ws4pybAswGHnDOLTSzE4FLgLbAP2b2KHBbTFGBHc4bx9HAmFz2nwhcHuc1O3tOwkft\n1I4ux19DbsOPXp6CL81+CL5Mu0i+mHN5bR9EwsPMMvGlU1Pij3cz+xBY5py7IIHnfABo75zbYfFB\nEfHM7HlgkHPu28I+r5m94JzrkGNfI3yJ6CtzOT7ucxI+aqfinqdsdM2m2H0DgRuAA6MFPET2mIYA\nJWWZ2c1mts7MjozZ18jMNprZpdHt0mb2mPnV1v8xs/lm9riZ7ZfjXNlmdoOZPWh+hfZlZnZT9Lmu\nZjbPzFaZ2TAzKxnzukuir21oZh+ZX419jpm1zUP8x5rZRDP7K/oYb2YHJ+5fSESACNsWQy2080bL\nUue2vk4v4K44L9vZc5KEimI7lbNjEzUTXwExL3PVRHZKnRtJahZdFyT2EfP0Q8AMYJSZpZmfDD0K\nmOScGxY9pjQ+PfNOfD79nUBzILe84huBffAlXZ8HHjCzQUAX4Dr8EPtF+LtPW2wZGn0BP+ehHfAt\nMN7MjtvJz1UTP0m5ZPScXYFjgNd3/a8iInnlnPsJWGBmjQv5vEeS4/NsZl2AiTkquu3yOQk3tVN5\nchI+hU/FMiT/nHN66JF0D3yufnYuj81A9ZjjagJ/4RuDh/FlcA/eyXmL4S+ym4GqMfuzgfditg1Y\ngl9zYp+Y/eOAT2O2u0Zf2yfHa38ExsTs+xAYH7P9XPSYYjH7auHv9LbK47/RA/jSt4H/vvTQI+wP\n/IhIhaDOi/8D8ezdfU6P8D7UTuX53+lofHnzvkH/zvRIjYcKCkgyWw20YMeJmku2fOP8xOo+wKP4\nkcrOLsedTzPrjE/3OBx/xwv8naza+IX7ttha1cg558zsZ+Af59w/Mcf8D2iSIx6HvxsW+9rX8AsB\nxtMCGBmNb8tdvl+ij4bA27m9SET2jHPu7iDP65zbQJw73jt7TkJP7dROmNmB+AVxvwZyWwtKZLcp\nLU2S2SaX+9ogOfPYX45+XQm8GPuEmbXDpwB8gr+In4gfkjd2XAE954rjG+Lsy23l9JxD7X8AlXL/\nsQAoj19NPeeq9YeybVVzEREJN7VTcURT8F4HSgDn5PJvIrJHNHIjRcFQfFnag4B+wB0xz50HfOac\nu27LDjNrVgAxHIRfwyN2+7edHL8S39g9zY53/HZ3hXkREQm3ItVOmV/TZix+7tlJzrk9XYNLZAfq\n3EhKi07C/Q8+P/l44P/M7BXn3BfRQ/bGr6we62LyvoZGnsLA32W7LxqTAecA03fymveBY5xzMxMY\nh4iIhEwRbaf+DzgDaOGcK8pr/kgBUOdGkllx84vp5bTIObfEzKric5jvjzYSX5hf2PBZMzs+msf+\nLvC4md2OX/jyP8CpBRDrZWa2Efgev4BZTeDCnRwfAT43s4nAcPxdsKrAacAI59xH8V4Y/RnB52KX\njtnOcs6tyNdPISIiu0PtVA7Rn+Ny4J7oduy/zyzn3Jp8/hxSxKlzI8nsAPxK4Dn1xV80nwYW4i/A\nW1wKfIdfRf0WfCrAofj1KEoB7wAdgc9ynDPeKvR5uXPm8GU5HwX641dkvsDtuLhf7IrwP0VLyA6I\nxrg3sBh/p2xXd7km5IhrS7nQ5kDcTpGIiCSc2qkdnR49z23RRyy1U5Jv5lz+RzXNrCXbqnwMc84N\nyvH82fgPSzZ+wlkv59wn+X5jkZAzs674O1r7OefWBh2PiIhILLVTkmryPXITnRT2OL4k4BJghpm9\n5pybHXPYe86516PHH4e/k3xUft9bRERERERki0SUgm4E/OScW+Cc24hf4fac2ANy3AnYFz+CIyIi\nIiIikjCJ6NxUwedmbvFrdN92zKytmf0IvAF0T8D7ioSec26Uc66YhvpFRCSM1E5Jqim0RTydc686\n544C2uInn4mIiIiIiCRMIqqlLQaqx2xXje7LlXPuYzM7zMzKOudW5nzezBJZt11ERPaQcy7nwnyC\n2ikRkbDIrZ1KROdmBlDLzGrgV7LtgC9RuJWZ1XTOzYt+Xx8omVvHJibQBIQFkUiESCSSkHNJ/kTM\niCTo9yr5o89FOFiG4bLC+5nwa/hJPAlrp3RtDA1dG8NDv4twSNZ2Kt+dG+fcZjPrga+7vqUU9I9m\ndqV/2j0FtI+uwLsBWAdckN/3FRERERERiZWQRTydc5OAI3LsGxrz/f3A/Yl4L0lS6elBRyASKumH\n6DMh6NooIqGVrO1UoRUUCEJGRkbQIUhUhoaXQ0Ofi3CIXBIJOgQJAV0bw0PXxvDQ7yIckrWdskTl\nDSeKmbmwxSQiUtSYmQoKxKF2SkQkePHaqYSkpYlI6jrkkENYsGBB0GFIAalRowa//PJL0GGIiIgk\nhEZuRGSnondGgg5DCki8369GbuJTOyUiErx47VRKz7kREREREZGiQ50bKRyaNCsisiNdG0UkpCJZ\nkaBD2CNKS5PCYQb6vSYlpaWlNqWl7b6EtlO6NopISFk/w2WG9/qUVGlpus6LSF4ceuihfPDBB0GH\nsdWCBQtIS0sjOzs76FBo3rw5w4cPDzoMERGRQhXKzs0FF8Dq1UFHISKy+8w02CEiIhKUUHZuXnwR\n6teH6dODjkREJHebN28OOgQRERHJIZSdmwYN4Oef4ZRT4JFHlKYmIru2YcMGbrjhBqpUqULVqlXp\n1asXGzduBPxq16+88goAn3zyCWlpabz99tsAfPDBB9SrV2/reYYPH87RRx9NuXLlaNWqFQsXLtz6\nXFpaGk888QS1a9emdu3aucbhnGPYsGFUqVKFKlWq8NBDD+UpxlGjRtG0adPtzpWWlsb8+fMB6Nat\nGz169KBNmzbsv//+NGnShJ9//nnrse+++y5HHXUUZcqU4brrrtM8KRERKZJC2bn55BO4/nrYuBFu\nvBHOPhtWrAg6KsmXzMygI5AUN2DAAKZPn863337LN998w/Tp0xkwYAAA6enpZGVlAfDRRx9Rs2ZN\nPvroIwCmTJlCRkYGAK+99hr33Xcfr776KsuWLaNp06Z07Nhxu/d57bXXmDFjBrNmzYobS1ZWFvPm\nzWPy5MkMGjRo67ygncUIO6a05dweN24c/fr1Y/Xq1dSsWZM77rgDgBUrVtC+fXvuueceli9fTs2a\nNfnkk092819QAqFro4iEVGZ6kl6fnHOheviQvFdece7AA50D56pWde7jj52IFLLYz2SuMjP9hzTn\nIzMz78fHO3YXDjnkEPf+++8755yrWbOmmzRp0tbnJk+e7A499FDnnHPvv/++q1u3rnPOuZYtW7ph\nw4a5Jk2aOOecS09Pd6+88opzzrlWrVq54cOHbz3H5s2bXenSpd3ChQudc86ZmcvKyoobzy+//OLM\nzM2dO3frvt69e7vLLrtslzGOHDnSNW3adLvzmZmbN2+ec865Sy65xF1++eVbn3vrrbfcUUcd5Zxz\n7tlnn93682xRtWpVN2zYsLixbhHv9xvdH3ibEMbHLj8TIiJS4OK1U6EcudmibVv4+mto0gR+/RXS\n0+HeeyEEhYhEZItIJLeuTfz1O3I7PgFrfSxZsoTq1atv3a5RowZLliwBoEmTJsydO5c//viDb775\nhi5durBo0SJWrFjB9OnTSU9PB3y1s549e1K2bFnKli1LuXLlMDMWL1689bxVq1bdaRxmtt0xsXHs\nLMa8qFix4tbvS5cuzd9//731vNWqVdvu2JzbIiIiRUGoOzcANWrAlCnQuzds3gy33w6tWsEffwQd\nmYiESeXKlVmwYMHW7QULFlC5cmUA9t57bxo0aMDgwYM59thjKV68OE2aNOHhhx+mVq1alClTBoDq\n1aszdOhQVq5cycqVK1m1ahV///03jRs33nrevFRDW7Ro0dbvFy5cuDWOncW4zz77sHbt2q3PLV26\nNM8/e6VKlbabG5QzBhERkaIi9J0bgBIlYNAgeOstKF8e3nkH6taFDz8MOjIRCYuOHTsyYMAAli9f\nzvLly+nfvz+dO3fe+nyzZs14/PHHt47SZGRkbLcNcOWVV3LPPfdsnU/z559/8uKLL+5WHM45+vfv\nz7p16/jhhx8YMWIEHTp02GWMdevW5YcffuDbb79l/fr19OvXL89lpVu3bs2sWbN49dVX2bx5M4MH\nD+b333/frbhFRERSQVJ0brZo1cqnqTVrBkuXQosWPptFFVlFiqbYP/7vvPNOGjZsSJ06dahbty4N\nGzbcOuEefFGBv//+m2bNmm23Hdu5adu2LbfeeisdOnTgwAMPpE6dOkyaNCnX99tZTOnp6dSqVYvT\nTz+d3r1706JFi13GePjhh3PXXXfRokULateuvUPltJ0pV64cEyZMoE+fPpQvX5558+Zx8skn5/n1\nIiIiqcL8fJzwMDO3q5g2bYL+/f3DOT8XZ8wYiGZ3SBhFIgmZVyGFz8wI23VCEife7ze6XyuS5iIv\n7VSe6dooIiEVyYoQyYgEHUZc8dqppOzcbPH++3DRRfD77z5d7bnnoGXLAg5Q9oyZFixKUurcpDZ1\nbnZfQjs3ujaKSEhZP8Nlhvf6FK+dSqq0tJxatIBvvoHTToPly33a2m23+fVxRERERESkaEnqzg3A\nwQfD5MkwYACkpcF990FGBuQoHCQiIiIiIiku6Ts34Ds1d9wBWVlQpQpMmwbHHw+vvx50ZCIiIiIi\nUlhSonOzRdOmvppa69awahWccw706gUbNgQdmYiIiIiIFLSU6tyALyzw+uvw4INQvDg8+iicfDLM\nnx90ZEVcZmbQEYiIhI+ujSISUpnpyXl9Supqabvy+edw4YWwYAHsvz8MGwbnnZeQU4sUGaqWltpU\nLW33JbRamoiI7JGUrJa2KyeeCDNnQrt28NdfcP75cO218O+/QUcmIiIiIiKJltKdG4AyZeCll2DI\nEChZEp54Apo0gblzg45MRFJZ8+bNGT58eNBhMGrUKJo2bRp0GCIiIoUi5Ts34NdIu+46X0WtZk1f\ndKBBAxgzJujIREQKnpmyy0REpGgoEp2bLRo0gK++8vNw/v4bLroILrsM1q4NOjIRSWbZ2dlBhyAi\nIiIUsc4N+MICY8fC0KFQqpQvMtCoEcyaFXRkKS4SCToCSUFpaWnMjymF2K1bN+666y4ApkyZQrVq\n1bj33nupUKEChx12GGNihmu7devG1VdfzRlnnMH+++9P8+bNWRiz+u/s2bM544wzKFeuHEcddRQT\nJkzY7rXXXHMNrVu3Zr/99iMrKyvX+P73v/9x4okncsABB9CuXTtWr1699bnXX3+dY489lrJly3Lq\nqacye/bs3fq5Hn74YQ4++GCqVKnCyJEjtx67cuVKzj77bA444AAaN27MvHnzdvNfVQqVro0iElKR\nrEjQIeyRIte5AZ+mdsUVvpraEUfADz9Aw4YwYgSoAE4B6dcv6AgkBe0q3Wrp0qWsXLmSJUuWMHLk\nSK644gp++umnrc+PGTOGzMxMVqxYQd26dbnooosAWLt2LWeccQYXX3wxy5cv54UXXuCaa67ZrgMy\nduxY+vbty5o1azjllFNyff/nnnuOkSNHsnTpUooVK8Z1110HwNy5c+nUqRNDhgxh2bJltGrVirPO\nOotNmzbl+edas2YNS5Ys4ZlnnuHaa6/lzz//BOCaa66hdOnS/P777wwbNiwU835kJ3RtFJGQ6jcl\nOa9PRbJzs0WdOvDFF9ClC6xbB927Q9euPmVNRPImkhXB+tkOj3h3fHI7fk/vDu2qHK+Z0b9/f0qU\nKEGzZs1o3bo148eP3/p869atOfnkkylRogQDBw7ks88+Y/Hixbz55psceuihdOnSBTOjbt26tG/f\nfrvRm3POOYfGjRsDULJkyVzfv3Pnzhx11FHsvffe9O/fnwkTJuCcY/z48bRp04ZTTz2VYsWKcfPN\nN7Nu3TqmTZuWp5+rZMmS9O3bl2LFitGqVSv23Xdf5syZQ3Z2Ni+//DL9+/enVKlSHHPMMXTt2jVP\n/5YiIiKpoHjQAQRt331h1Cho3tyXiX7uOZg+HcaP950fEdm5SEaESEakwI7PjzJlylCqVKmt2zVq\n1GDJkiVbt6tVq7b1+3322YcyZcqwZMkSFixYwGeffUbZsmUB39nYvHkzXbp0yfW18cQeU6NGDTZu\n3Mjy5ctZsmQJNWrU2PqcmVGtWjUWL16cp5+rXLlypKVtuzdVunRp/v77b5YtW8bmzZupWrXqdu87\nderUPJ1XREQk2RXpkZtYl1wCM2bAMcfAnDl+Hs7QoUpTEwmz0qVLszamIsjSpUu3e37VqlWsW7du\n6/bChQupXLny1u1FixZt/f7vv/9m1apVVK5cmWrVqpGRkcHKlStZuXIlq1at4q+//uLxxx/fenxe\nKpDFnn/BggWUKFGC8uXLU7lyZRYsWLDDsVs6Jbv6ueKpUKECxYoV2+59Y+cRiYiIpDp1bmIcfbQf\ntbnsMli/Hq66Cjp29AuAikj41KtXjzFjxpCdnc2kSZOYMmXKds8758jMzGTjxo1MnTqViRMncsEF\nF2x9/q233mLatGls2LCBvn370rhxY6pUqUKbNm2YO3cuo0ePZtOmTWzcuJEvvviCOXPm7FZ8o0eP\nZvbs2axdu5bMzEzOP/98zIwLLriAiRMn8uGHH7Jp0yYefPBBSpUqRZMmTfL0c8WTlpZG+/btiUQi\nrFu3jlmzZjFq1KjdillERCSZqXOTQ+nS8PTT8PzzPmVt3DioX9+XkJZ8yMwMOgJJQY8++iivv/46\nZcqUYezYsbRr12675ytVqkSZMmWoXLkynTt3ZujQoRx++OFbn+/UqRORSIRy5coxc+ZMRo8eDcC+\n++7LO++8wwsvvEDlypWpXLkyt956K+vXr89zbGZG586d6dq1K5UrV2bDhg0MHjwYgNq1azN69Gh6\n9OhBhQoVmDhxIm+88QbFixfP08+V23tt8dhjj7FmzRoqVapE9+7d6d69e55jlgDo2igiIZWZnpzX\nJ9vVxNXCZmYuLDHNnevXxPn6ayhZEh58EHr08NXWRIoKM9vlBPcwmjJlCp07d46bltWtWzeqVavG\n3XffXciRhUu83290v652uQhTOyUiUlTFa6cSMnJjZi3NbLaZzTWzPrk838nMvok+Pjaz4xLxvgWt\ndm349FO45hrYsAGuvx7at4dVq4KOTEREREREcsp358bM0oDHgTOBY4COZnZkjsPmA82cc3WBAcDT\n+X3fwlKqFPz3v7562v77wyuv+DS16dODjkxE8iMvBQFEREQkueQ7Lc3MGgOZzrlW0e1bAeecGxTn\n+AOB75xzudZRDfNw/7x5Pk3tyy+heHEYNAh69VKamqS2ZE1Lk7xRWtruC3M7JSJSVBRkWloVYFHM\n9q/RffFcBrydgPctdDVrwiefQM+esGkT3HQTnH02rFgRdGQiIiIiIlKo1dLMrDnQDdhhXk6y2Gsv\nePRRn5524IHw5ptQrx5EFxaXeCKRoCMQEQkfXRtFJKQiWZGgQ9gjiUpLizjnWka3c01LM7M6wEtA\nS+fcvJ2cz2XGlMbMyMggIyMjXzEWlAULfJra559DsWIwcCDccgukqcD2jsy0ImqSUlpaatvy+83K\nyiIrK2vr/n79+iktLY6EtlO6NopISFk/w2WG5/qU13YqEZ2bYsAcoAXwGzAd6Oic+zHmmOrA+0Bn\n59xnuzhfUuUyb9wIt9/uy0QDtGwJzz4LFSoEG1foqAFPWocccggLFiwIOgwpIDVq1OCXX37ZYb/m\n3MSX0HZK10YRCamwdW5yitdOJWSdGzNrCQzGp7kNc87dZ2ZX4kdwnjKzp4FzgQWAARudc43inCup\nOjdbTJwIXbrAypVQuTKMHQvNmgUdVYioARfZTrI2GqLOjYgUDcnaTmkRzwRatAg6dvRFB9LS4O67\n4bbblKYGqAEXySFZGw1R50ZEioZkbaf0Z3cCVasGWVm+Q5OdDXfe6dPUfv896MhERERERFKfOjcJ\nVrw43HMPTJoE5cvDu+/C8cfDBx8EHVnAYibfighkpuszIejaKCKhlaztlNLSCtDixdCpE3z0kc88\nuOsu6NvXV1YTEQkzpaXFl0rtlIhIslJaWgCqVIH33/cdGoB+/eD00+G334KNS0REREQkFWnkppC8\n9x5cfLGff3PQQTB6tO/oiIiEkUZu4kvVdkpEJJlo5CZgp50GX38Np54Kf/wBZ57pCw5s2hR0ZCIi\nIiIiqUGdm0JUsSK8844vEW0GAwdCixZ+bo6IiIiIiOSPOjeFrFgxPwfn/fehUiVfbOD442Hy5KAj\nK2CRSNARiIRKJCsSdAgSBro2ikhIJWs7pTk3AfrjDz8P5913/fatt0L//r6cdMrRQnUi20nWxdFE\ni3iKSNGQrO2URm4CdNBBfj2cgQMhLQ3uuw8yMmDRoqAjExERERFJPurcBCwtDW6/HT78ECpXhk8+\ngXr14K301TmaAAAgAElEQVS3go5MRERERCS5qHMTEs2a+WpqLVvCihXQujX07g0bNwYdmYiIiIhI\nclDnJkQqVICJE316WrFi8MADkJ4OCxcGHZmIiIiISPipcxMyaWnQpw9kZUHVqvDpp76a2htvBB1Z\nPmVmBh2BSKhkpuszIejaKCKhlaztlKqlhdjy5dC167b5NzfeCPfeCyVLBhuXiKQ+VUuLT+2UiEjw\nVC0tCZUv70ds7r/fp6k9/LCfm7NgQdCRiYiIiIiEj0ZuksS0adChgy8TfeCBMHIknHNO0FGJSKrS\nyE18aqdERIKnkZskd9JJMHMmtGkDq1dD27Y+TW3DhqAjExEREREJB3Vukki5cvD66/Dgg1C8ODzy\nCDRtCr/8EnRkIiIiIiLBU+cmyZjBTTfB1KlQvTpMn+4X/Xz11aAj24VIJOgIREIlkhUJOgQJA10b\nRSSkkrWd0pybJLZyJVxyybYy0T17+uIDoaymZgb6vYpsZf0Mlxnez4Tm3MSX0HZK10YRCalkbac0\ncpPEypaF117zVdSKF4fBg+Hkk2H+/KAjExEREREpfOrcJDkz6NULPv4YatSAL76A+vXh5ZeDjkxE\nREREpHCpc5MiTjzRV1Nr2xb+/BPat/dpauvXBx2ZiIiIiEjhUOcmhZQp40dsHn0USpSAIUPglFOU\npiYiIiIiRYM6NynGzI/YfPIJHHJIiNLUMjMDDkAkXDLT9ZlIZitXJuhEujaKSEglazulamkpbPVq\n6N4dXnnFb19/va+mttdewcYlIuGnamnxmZnbd1/Hddf5xZTLlw86IhGRoideO6XOTYpzDh57DG6+\nGTZuhIYNYdw4OOywoCMTkTBT5yY+M3Pg26l99oEePfz6YxUqBByYiEgRolLQRZSZH7EJXZqaiEgS\nmzYNWraEf/6BQYP89fWWW+D334OOTESkaNPITRGSM03tuuvggQeUpiYiO9LITXyx7dT06XD33TBx\non9u773hqqugd2+oWDHAIEVEUpzS0gTYMU2tQQMYP15paiKyPXVu4sutnfriC9/JeeMNv12qFFxx\nhe/kVKkSQJAiIilOaWkC7Jim9uWXUK8evPRSAb9xJFLAbyCSXCJZkaBDkARq2BBefx2++gratYN/\n//Xl+GvW9HNyFi2K80JdG0UkpJK1ndLITRGWM02tRw948MECSlMz88NGIgKA9TNcZng/Exq5iS8v\n7dS338KAAfDii/7SV6KEv97edhvUqLHdyXRtFJFQStZ2SiM3RdiBB/oRmy2Lfj7+OJx8shb9FBHJ\nrzp1fMrvd99Bhw6waRMMHQq1asHll+s6KyJSUNS5KeK2LPr58ceFnKYmIlIEHHMMjB0Ls2bBxRdD\ndjY88wzUrg3dugUdnYhI6lHnRgBo1MjnirdtC3/9Beed5zs969cHHZmISPI78kh47jmYPRu6dvX7\nRo70Xy++GH78MbDQRERSijo3slWZMn79m0ce8WlqQ4ZA06bw889BRyYikhoOP9x3aubOhcsu8/ue\nf96P8HToAN9/H2h4IiJJT50b2Y4Z3HDDtjS1GTN8mtqWogN7LDMzEeGJpIzMdH0mirLDDoOnn4bV\nPTO5+mp/Q2ncODjuOGjfHr7+OugIRaSoS9Z2KiHV0sysJfAovrM0zDk3KMfzRwAjgPrA7c65h3dy\nLlVLC4lVq3x1n1df9ds9e8L990PJksHGJSIFT9XS4iuIdurXX2HQIN/h2ZIOfNZZ0LcvnHBCQt9K\nRCQlFNginmaWBswFWgBLgBlAB+fc7JhjygM1gLbAKnVukodzPj3tllv8op8nnOArAB1ySNCRiUhB\nUucmvoJsp377zZfk/7//g3Xr/L6WLeGuu6BJkwJ5SxGRpFSQpaAbAT855xY45zYCLwDnxB7gnFvu\nnPsS2JSA95NCFFtNrUaNbWlqW0ZzREQkcSpVgocegl9+gd69YZ99YNIkOOkkOO00+OijoCMUEQm3\nRHRuqgCxay//Gt0nKaRRI5g5E845xy/+2a4d3HgjbNgQdGQiIqnnoIN8mtovv8Add8B++8H770N6\nOmRkwAcfaO1PEZHcqKCA5FmZMr6wwMMPQ/Hivqpa06a+8RURkcQrXx4GDIAFCyAS8YsvT5kCLVrA\nKafA5Mnq5IiIxCqegHMsBqrHbFeN7ttjkUhk6/cZGRlkZGTk53SSQGbQq5dPkbjgApg+3aepjRoF\nZ5+9kxdGIv4hIgBEsiJEMiJBh7FVVlYWWVlZQYeRNBLWTuXx2limjC86ecMN8N//+tS1adP8fJwT\nTvBzclq39tdoEZFESNZ2KhEFBYoBc/AFBX4DpgMdnXM7LElmZpnA3865h3ZyPhUUSBIrV8Ill8Ab\nb/jtG2+Ee++NU03NTLcXRWJYP8NlhvczoYIC8SW0ndrDa+OaNb7owIMPwrJlfl+9er662jnnQJry\nMkQkn5K1ncr35c85txnoAbwD/AC84Jz70cyuNLMrom9+sJktAnoBd5jZQjPbN7/vLcEqWxZee803\nrsWL+3S19HRYuDDoyEREUtt++/mCAz//7K+9FSv6eZHnngvHHw8TJkB2dtBRiogUvoSsc5NIGrlJ\nTp9+ChdeCIsW+fSJZ5+FNm1iDtDIjch2kvWOmIRj5Candetg2DC47z5YHE0MP+oouPNOf20uVizf\nbyEiRUyytlMauJaEaNLE3zX8z3/84p9nneXvKm7cGHRkIiKpb++9oUcPmDfPp6tVrw4//ggXXQRH\nH+1vOG3SYgwiUgSocyMJU66cn38zaJC/S/jAA75k6aJFu3ypiIgkwF57wVVXwU8/wTPPwKGHwty5\n0LUrHHkkDB+um04iktrUuZGESkvzIzZZWVCliq/mU68e/NQpM+jQREIlM12fCcGXQCsAJUvCpZfC\nnDkwYgTUquVHdS69FGrXhqef1jplIrJzydpOac6NFJhly6BLF7+6NsCtt0L//r74gIiEm+bcxJeM\n7dSmTfDCC37NnDlz/L7q1f11uXt3P+IjIpJM4rVT6txIgcrO9hNc+/b13zdtCmPH+lEdEQkvdW7i\nS+Z2avNmX0mtf3+YNcvvq1LFd3IuuwxKlQo2PhGRvFLnRgI1ZQp07Ai//eZX3B49Gs48M+ioRCQe\ndW7iS4V2KjsbXnrJd3K++87vq1wZ+vSByy/3BQpERMJM1dIkUOnp8PXXcPrpsHw5tGrlS5Sqeo+I\nSOFLS4Pzz/fX5Zdegrp1YckS6NkTataEwYN9eWkRkWSjkRspVJs3w733+jm02dm+0zN2LFSqFHRk\nIhJLIzfxpWI7lZ0Nr78Od9/ty/qDXxi0Tx+48kqN5IhI+GjkRoIViQC+RPSdd8J77/mGc8oUv5r2\ne+8FG55IYYtkRYIOQcIgem0MWloatG0LX37pOzn168PSpdCrly8n/cgjsHZt0FGKSGFK1nZKIzdS\nOHJZhfv33/0Cc++/75++6y5feEAraUtRkKwrP0uC26lcro1h4BxMnOj7Xl9+6fcdfPC2kZzSpQMN\nT0QKQbK2Uxq5kcAcfDBMnrztxmW/fnDGGf5uoYiIBMcM2rSBGTP84swNGvgbUjfeCIcdBo8+qjk5\nIhJO6txIoIoV8/Nv3nkHDjoIPvjAL/qZlRV0ZCIiEq+T06uXLzwwZAj8+2/QUYqIbKPOjYTCaaf5\nqj3p6X7kpkULv9hcdnbQkYmISGwn57XX/E2o337z1dVq1YInnoD164OOUkREnRsJkUqVfGGBO+7w\nnZq+fX3J6GXLgo5MRETAd3LOPtvPw3nlFV9CevFiuPZaOPxwGDoUNmwIOkoRKcrUuZHCkZmZp8OK\nF/cjNpMm+cU+33nHV1ObOrWA4xMpZJnpeftMSIrL47UxbMx8dbWvvoIXX4RjjoFFi+Cqq+CII2DE\nCK1jJpLskrWdUrU0Ca1ff4UOHeCTT/zcnIED4ZZbfMlSESlYqpYWn9qpHWVnw4QJvkDM7Nl+3+GH\n+0IxF16o67aIJJ6qpUnSqVoVPvwQevf2i3/eeiucdRasWBF0ZCIiEistzXdivv8ennvOz8P56Sfo\n1Mmnrr36aigrXotICtLIjSSFN9+ELl1g1SqoVg3Gj4fGjYOOSiR1aeQmPrVTu7ZpE4waBXffDQsX\n+n0NG8I99/gCMqb/WSKSTxq5kaTWpg3MnAknnujzups29Stm6+8LEZHwKV4cLr0U5s6Fxx+HihXh\niy/8WmanngrTpgUdoYikKo3cSFLZsMGnpz3yiN9u2xaGD4cyZYKNSyTVaOQmPrVTu2/tWt/Jue8+\nPwIPPs144EA47rhgYxOR5KSRGwlWJJKQ05QsCQ8/7EuQHnCAz+OuX9/fERRJJpGsSNAhSBgk6NoY\ndqVL+/mT8+fDnXfCPvv4RUHr1oXOneHnn4OOUERyStZ2SiM3UjjMEp5DNn8+XHCBX2+hZEl46CG/\n1oJyuSUZWD/DZYb3WqeRm/gS2k4VwLUxGfz+u59/83//Bxs3QokScM01fp2zChWCjk5EIHnbKY3c\nSNI67DBfJrpHD5+udt11vlrPn38GHZmIiOzMwQfD4MEwZ44fudm0yW/XrOnXOlu7NugIRSRZqXMj\nSW2vveCxx3z1tP328+ssNGwIX38ddGQiIrIrhx4Kzz7rC8a0agVr1kDfvn6NnGHD/DIAIiK7Q50b\nSQnnn+/T0+rWhf/9z5eJfuqpIpntISKSdOrWhbfegg8+gAYNYMkSuOwyOP54mDw56OhEJJmocyMp\n4/DD4dNP4YorYP16uPJKuPhi+PvvoCMTEZG8aN4cpk+HMWPgkEP8oqAtW/pRnR9+CDo6EUkG6txI\n4cjMLJS32XtvGDoURo/21XjGjIETTvANpEiYZKYXzmdCQq6Qro3JJC0NOnaE2bPhgQdg//1h0iSo\nUweuvhqWLQs6QpGiIVnbKVVLk5Q1ezacd56/27f33vDEE3DJJUFHJZIcVC0tPrVThWvZMujXD558\n0s/BOeAA3ye89lpfKVNEiqZ47ZQ6N5LS/vnHV1MbOdJvd+vmF5IrXTrQsERCT52b+NROBWPWLLjx\nxm1zcI44AoYMgTPOCDYuEQmGOjdSpI0Y4e/yrVsHxx4LL77oG0YRyZ06N/GpnQqOc77wwI03wty5\nft+55/rFnWvUCDY2ESlcWudGirRu3eDzz32H5vvvfTWeMWOCjkpERHaHGbRuDd99B4MG+bmVL78M\nRx4JAwf6YjIiUrRp5EaKlDVrfDW1F17w21ddBY88AqVKBRuXSNho5CY+tVPhsXgx3HILjB3rt488\n0s+vbN482LhEpOBp5EaCFYkEHQHgF/ocM8Y3fiVL+gmqJ50E8+YFHZkUNZGsSNAhSBiE5NqYrKpU\n8df099/3I/OzZ8Opp0KXLrB8edDRiSS3ZG2nNHIjhcMsdCtqfvWVX/xz/nxfanTECJ+7LVIYrJ/h\nMsP1mYilkZv4EtpOhfDamKzWr4cHH4QBA+Dff6F8eXj0UejUyf8zi8juSdZ2SiM3UmTVrw9ffgnt\n2sFff0H79tCrF2zYEHRkIiKyu/baC+64w8/Had7cj9xcfDH85z+waFHQ0YlIYVHnRoq0Aw+El17y\n826KF/d3+Zo1gwULgo5MRET2RK1aPk1t2DB/jZ80yVfJHD5cg2QiRYE6N1LkmcENN8DUqVCtmq+q\nVq8evPlm0JGJiMieMIPu3eHHH+Hss/3o/KWX+kprS5YEHZ2IFKSEdG7MrKWZzTazuWbWJ84xQ8zs\nJzP72syOT8T7iiRS48Ywc6ZPYVi1Cs46C269FTZtCjoyERHZExUrwquvwnPP+VGct9+GOnX8PhFJ\nTfnu3JhZGvA4cCZwDNDRzI7McUwroKZz7nDgSuDJ/L6vJJnMzKAjyJNy5eCNN+C++6BYMb+Owqmn\n+nKjIomUmZ4cnwkpYElybUxmZn7uzQ8/wBlnwIoVfq7lFVfAP/8EHZ1IeCVrO5Xvamlm1hjIdM61\nim7fCjjn3KCYY54EPnTOjYtu/whkOOd+z+V8qpYmofDRR9ChA/z2G1So4MuNnnZa0FGJFA5VS4tP\n7VTyys6Gxx6DPn18dbXatf21vUGDoCMTkd1VkNXSqgCxdUh+je7b2TGLczlGJFSaNYOvv/YdmmXL\n/B2/fv1g8+agIxMRkT2RlgY9e8KMGb7IwNy50KQJ3H+/7/iISPIrHnQAImF20EG+0s6AAb5jE4nA\nxx/D88/750REJPkcd5zv4PTpA0OG+K9jx0LVqkFHlhi1akHv3lCpUtCRiBS+RHRuFgPVY7arRvfl\nPKbaLo7ZyrTaloTYe+/BwQcHHYWIBEntVOr5+mv/SBWPPhp0BCLBSETnZgZQy8xqAL8BHYCOOY55\nHbgWGBedo7M6t/k2WyiXWcJqyRI/D2fqVF9wYOBAuOUWn+ogkkr0x/vOqZ1KLatWwbRpqZGatnkz\njBwJr73mtytV8ml3F13kiyuIpIp47VS+CwpET94SGIyfwzPMOXefmV2JLyzwVPSYx4GWwD9AN+fc\nV3HOpYmaqSgS8Y8UsGkT3Hmnr6QG0KYNjBoFZcsGG5ckl0hWhEhGJOgw4lJBgfgS2k6l0LVRwmXK\nFLjpJvjyS7/drBn8979+rpFIXiRrO5WQzk0iqXOTosxSbmnoN9+ELl38Hb/q1WHCBGjUKOioJFlY\nP8Nlhvczoc5NfAltp1Lw2ijhkZ3tb7717g3Ll/uMg+uv9/3p/fcPOjoJu2Rtp5RMI7KH2rTxi342\nagQLF8Ipp/iJqfo7RUREwiAtDbp181XhrrnGt0+PPAJHHukLKKi9klSkzo1IPtSo4eff9OwJGzf6\nrxdcAH/+GXRkIiIiXpkyPiVtxgxo3Niv39apk1+k+ocfgo5OJLHUuRHJp5IlfVWaCRNgv/3gxReh\nYcPUqrojIiLJr359+OQTGDYMypeHrCw4/nhfGGfNmqCjE0kMdW5EEuS88/zEzbp14X//83fHnn5a\nw/4iIhIeaWnQvTvMmQNXX+2rqz34oE9Ve+EFtVmS/NS5kcKRmRl0BIXi8MPh00/hsstg/Xq44gro\n2hX++SfoyCRsMtOLxmdCdqGIXBslfMqWhSeegOnT/dzRJUugY0do0QJmzQo6OgmDZG2nVC1NpIA8\n+6y/K7Z2LRx9tE9XO+qooKMSyRtVS4tP7ZSkmuxsn6p2222wYgUULw69ekHfvj7dWiSMVApaJAA/\n/ADnnw8//gj77ANPPeUncYqEnTo38amdklS1YoVfx23oUJ+eVrkyPPywL5SjBUAlbFQKWiQAxxzj\nh/wvusinpl10kR/N+fffoCMTERHZXrly8H//59utE07wqWodOsBppylVTZKHOjciBWzffeG55/yd\nsL32giefhJNPhvnzg45MRERkRw0bwmef+aI45crBBx/4Yjm9e6uqmoSf0tJECtFXX/k0tfnz4YAD\nYMQIaNcu6KhEdqS0tPjUTklRkjNVrUoVeOghpapJ8JSWJsGKRIKOIBTq1/flotu18wt9nnsu3Hgj\nbNgQdGRS2CJZkaBDkDDQtVFCLmeq2uLF21LVfvwx6OikICVrO6WRGykcZiqeH8M5GDzYL5y2aRM0\naQLjxkG1akFHJoXF+hkuM7yfCY3cxJfQdkrXRkki2dnwzDO+qtrKlb6q2g03wF13qapaKkrWdkoj\nNyIBMPMNwkcfQdWqfm2cevVg0qSgIxMREcldWppfv23uXLjyyu0XAB03Tv10CQd1bkQC1KQJzJwJ\nLVv6vOZWrXxu8+bNQUcmIiKSu3LlfHGczz9XVTUJH3VuRAJWvjxMnAgDBvi7YgMHwumnw9KlQUcm\nIiIS3wkn+KpqTz0FZctuq6p2yy2qqibBUedGJATS0uCOO+Ddd+Hgg+HDD32a2pQpQUcmIiISX1oa\nXH557qlqY8cqVU0Knzo3UjgyM4OOICmceqpPU0tP9yM3p54K993nJ3FKaslM12dC0LVRUsaWVLXp\n06FRI5+q1qmTb8d++CHo6GRPJGs7pWppIiG0aZOvPnPvvX67dWt49lk/7C9SGFQtLT61UyI7l50N\nw4fDrbf6+aTFi0PPnr4vr6pqkijx2il1bkRCbOJE6NwZVq2C6tVhwgR/R0ykoKlzE5/aKZG8WbnS\np1xvWQC0UiW/AGiHDloAVPJPpaBFklDr1j5NrVEjWLgQTjkFHn9cOcwiIhJ+Zcv6BUBnzIATT4Tf\nfvOpas2bw/ffBx2dpCp1bkRCrkYNmDoVrr8eNm6E667zd73++ivoyERERHatQQOYNs0vAFq+vC+W\nc/zxcOONassk8ZSWJpJEJkyA7t3h77/h8MPhxRehTp2go5JUpLS0+NROiey5lSuhb19ffCA7GypW\n9NXVOnVSqprsHqWlSbAikaAjSAnnnw9ffgnHHQc//eSH+UeMCDoq2RORrEjQIUgY6NooRUzZsvDf\n//pUtcaNfWXQiy+GjAz47rugo5NYydpOaeRGCoeZJook0Nq10KPHto5Nt25+Lk7p0sHGJXln/QyX\nGd7PhEZu4ktoO6VroxRh2dkwciT06QPLl0OxYr5t69cPDjgg6OgkWdspjdyIJKHSpX2ZzeHDoVQp\n38lp3NgvoiYiIpIM0tJ8qvXcuXDttb6fP3gwHHGEX/5A/X7ZE+rciCSxbt3g88/9/JvvvvOTNseP\nDzoqERGRvCtTxmcffPEFnHQS/P47dO0KTZvCN98EHZ0kG3VuRJJcnTq+QTj/fF9o4MILfUW19euD\njkxERCTv6tXz1UFHjoSDDoJPPoH69X210NWrg45OkoU6NyIpYP/9Ydw4GDIESpTwd8CaNYMFC4KO\nTEREJO/S0vyozZw5vlMD8NhjPlVt5Eg/T0dkZ9S5kcKRmRl0BCnPzI/YTJ0K1avD9On+LthbbwUd\nmeQmM12fCUHXRpE4DjzQz7+ZOdMvYP3HHz4Vu2lTv08KXrK2U6qWJpKCVqyALl22dWxuv91Xnyle\nPNi4JHmoWlp8aqdECpdzMHo03HKLn4+TlgZXXQUDBvj5OlI0qVqaSBFSrhy88Qbcc49vBO65B04/\n3a8nICIikkzMoHNnn6p2ww1++4knoHZtGDZMqWqyPY3ciKS4rCzo0MHf7apYEV54AdLTg45Kwk4j\nN/GpnRIJ1vff+/Vwpkzx2yee6BcGbdAg2LikcGnkRqSIysiAr7/2X5cuhVNPhXvv1Z0uERFJTsce\nCx9+CM8/D5Uq+SURTjjBp6qtWBF0dBI0dW5EioCKFeHdd+G223yn5vbb4eyzYeXKoCMTERHZfWbQ\nqRPMng033QTFisHQoT5V7amnYPPmoCOUoKhzI4UjEgk6giKveHE/9+bNN/0EzIkT/foBM2YEHVnR\nFMmKBB2ChIGujSL5sv/+8OCDfrHP5s39Tbsrr4TGjX3VUNlzydpOac6NFA4zX+5EQmHBAr/o54wZ\nfl2cRx6Ba67xvyYpHNbPcJnh/Uxozk18CW2ndG0USRjnYPx4P5KzeLH/eF16qb+xV6FC0NEln2Rt\npzRyI1IE1ajh18Pp0QM2bvRfO3WCNWuCjkxERGTPmMGFF/pUtT59fMbCM8/4BUCfeEKpakWFOjci\nRdRee/lVn194Afbd13894QRfhUZERCRZ7bsv3HcffPutXwZh1Sq49lrfxn36adDRSUHLV+fGzMqY\n2TtmNsfMJpvZAXGOG2Zmv5vZt/l5PxFJvAsvhC++gGOO8WsINGoEzz0XdFQiIiL5c+SRMHkyvPgi\nVKsGM2fCSSdBt27wxx9BRycFJb8jN7cC7znnjgA+AG6Lc9wI4Mx8vpeIFJAjjvClNLt0gXXr/Ncr\nroB//w06MhERkT1nBu3bw48/+kqhJUvCyJG+qtqQIbBpU9ARSqLlt3NzDjAq+v0ooG1uBznnPgZW\n5fO9JJllZgYdgezCPvv4C/7TT/uUtaefhiZNYN68oCNLTZnp+kwIujaKFJJ99oGBA33qdatW8Oef\n0LOnX/hz6tSgowunZG2n8lUtzcxWOufKxtvOcWwN4A3nXJ1dnFPV0kQCNnOmr6Y2bx4ccIDv9LTN\n9daFpCpVS4tP7ZRIcnMO3njDd25++cXvu/hiuP9+vyioJIc9rpZmZu+a2bcxj++iX8/O5XBd7UVS\nQL168OWX0K6dv7vVrh3cfLOvrCYiIpLMzPxC1rNm+cHTvfaC0aN9ivbDD6utS3b5Hbn5Echwzv1u\nZhWBD51zR8U5Ns8jN5kxw/QZGRlkZGTscYwisuecg0cfhd69fV7yySfDuHFQpUrQkUmiZWVlkZWV\ntXW7X79+GrmJQ+2USGqZPx9uuMGP5oAvsPPYY35RUAmPvLZT+e3cDAJWOucGmVkfoIxz7tY4xx6C\n79wct4tzarhfJGSmTYMLLvCLolWoAGPGwGmnBR2VFCSlpcWndkokNU2c6FPVtsw17dABHngAqlYN\nNi7JXUEt4jkION3M5gAtgPuib1bJzN6MefMxwDSgtpktNLNu+XxfESlEJ53k5+GcfjosWwZnnAF3\n3w3Z2UFHJiIikhitW/uCA/37w957+/XfjjzSz8XZsCHo6CSv8tW5cc6tdM6d5pw7wjl3hnNudXT/\nb865NjHHdXLOVXbO7eWcq+6cG5HfwCXJRCJBRyD5VKECvP32tl9lZqavOLNsWaBhJa1IViToECQM\ndG0UCZVSpeDOO33p6Hbt4J9/oE8fqFMH3n036OgKV7K2U/lKSysIGu5PUWZ+AoekhHfegYsuguXL\n/XD9+PG+bLTknfUzXGZ4PxNKS4svoe2Uro0ioTZ5Mlx/Pcyd67fPPRceeQSqVw82rsKQrO1UftPS\nRKQIOuMMn6bWpAn8+is0awaDB+tvNBERSS1nngnffQf33efXynn5ZZ+qNmCAFroOK3VuRGSPVK0K\nWVm+wsymTf7rBRfAX38FHZmIiEjilCzpU9Nmz/bt3Lp10LcvHHssvPVW0NFJTurciMgeK1nSD8+/\n+CLst5//2rAhfPtt0JGJiIgkVtWqfjmE99+Ho4/2VdVat/Zr5syfH3R0soU6NyKSb+3b+0U/69SB\nn7djEi4AAB80SURBVH6CE0+EkSODjkpERCTxTj0Vvv4aHnrI39h74w3f2YlE/KiOBEudGykcMQve\nSWo6/HD47DPo3t3nIXfrBpdeqgt9PJnp+kwIujaKJKkSJeDGG2HOHLj4Yli/Hvr1852c115LjTmo\nydpOqVqaiCTc8OFw7bW+k1O3rk9Xq1Ur6Khkd6haWnxqp0Qkp6lToUePbWnZLVvCkCH+xp8UDFVL\nE5FC0727H8WpVQu++QYaNPAVZkRERFJR06Y+PXvIEDjgAJg0CY45Bm6/3a+VI4VHnRsRKRB168IX\nX/j5OH/95b/edBNs3Bh0ZCIiIolXvDhcd51fE6dbN9/e3XuvLx09YUJqpKolA6WliUiBcg4efRR6\n9/Ylo086yVebqVo16MhkZ5SWFp/aKRHJi88+8ynaX33lt1u0gMceg6OOCjauVKG0NBEJhBn06gVT\npkCVKjBtGtSrB+++G3RkIiIiBadxY5g+HZ58EsqW9SWk69SBm2/WmnAFSZ0bKRyRSNARSMBOOglm\nzoTTT4fly/2qz3ffDdnZQUcWjEhWJOgQJAx0bRRJacWKwZVX+lS1K6+EzZt9Cekjj4Tnnw93qlqy\ntlNKS5PCYRbuT7AUms2boX9/37FxzndyRo+G8uWDjqxwWT/DZYb3M6G0tPgS2k7p2ihSpHz5pU9V\n+/xzv920KTz+uB/RCZtkbac0ciMihapYMX+z+u23oVw5mDwZ6tf3uckiIiKprEEDn549fDhUqOBL\nSNevDz17wurVQUeXGtS5EZFAnHmmT1Nr3BgWLYJmzfxES93EFhGRVJaW5qupzZ3rq6s550tI164N\nI0YU3XTtRFHnRkQCU62aLzRwww2+ZOb110OHDrBmTdCRiYiIFKwDD/Sdmpkz4ZRTYNkyv07cySf7\n9DXZM+rciEigSpaERx7xawDstx+MHw8nnADffx90ZCIiIgWvTh346CM//7RiRZ+mfcIJcPXVsGJF\n0NElH3VupHBkZgYdgYTceef5RT+POw7mzIFGjeC554KOquBkpuszIejaKCKAry1y0UW+/bvpJj8/\n9cknfara0KG+GE9hS9Z2StXSRCRU1q71d6uefdZvX3EFDB4MpUoFG1dRo2pp8amdEpGCNmuWn4/z\nwQd+u359X1WtSZNg4woTVUsTkaRQujSMHAlPPQV77eW/nnQSzJ8fdGQiIiKF4+ij4b33fKp21arw\n1Ve+LezWDf74I+jowk2dGxEJHTO4/HL49FM47DA/2bJ+fXjttaAjExERKRxmcP75MHs23H67n6M6\ncqRPVRsyBDZtCjrCcFLnRkRCq149XzHmnHPgzz+hbVvo00cXdBERKTr22QcGDoTvvoOWLX172LOn\nv+n30UdBRxc+6tyISKgdeCC88grcf7+fYHn//dCiBfz2W9CRiYiIFJ7ateGtt3wWw6GH+s5Oerov\nRLDk/9u78/Aq6rP/4+877IhgCIgJhIAsiiIU5SfwsAWotJRFKYICAmKfn61bXeoCCuakUIWqVNR6\nPfwUCg8UUBQRBRdaTFAp4gLiBQKKGpAIZQkIBtny/f0x8TSEBAJJzpwz+byu61w5w0zO3OHkO3fu\nMzP3N9vv6KKHihuJjFDI7wgkhpnBfffBO+9AYqL3SVW7dpCR4XdkZy+UEfI7BIkGOjaKyBkwgwED\nYP167/BRvTrMnQsXXQSPPQZHjpTdvmI1T6lbmkSGmaaelzKxcycMHeoVOnFx3qn6++/3nscSSzdc\nWvSOCXVLK16Z5ikdG0WkFL7+Gu6++z/3pF58sXc/zlVXlf61YzVPxdifAyJS0TVoAMuWwUMPQV4e\njB3r3YuTk+N3ZCIiIpHVtCksWgRvvAEtWnjNB3r3hkGDICvL7+j8oeJGRGJOpUowcSK8/jrEx8Nr\nr3k3Vn78sd+RiYiIRN4vf+ndg/Poo96UCgsXQqtWXq788Ue/o4ssFTciErP69vV6/7dvD998480B\nMG2arvIREZGKp1o1GDPGO3tz3XVw6BCMHw+tW8OSJX5HFzkqbkQkpjVpAu+9B7fc4t1I+bvfwahR\n8MMPfkcmIiISecnJMH8+LF/uTQa6ZQv06wf9+3vPg07FjURGWprfEUiAVasGzz4Lc+Z4p+Nnz4YO\nHWDTJr8jK15ad40JQcdGESk3PXrA2rXwl79A7drepdyXXuqdzcnNPf33x2qeUrc0EQmUDRu8Gyk3\nboRatWD6dBgyxO+oYo+6pRVPeUpEYs2OHd4la7NmecuNG3tFz8CBXtPGWKRuaSJSIVxyCaxe7V1v\nfPCg9/X3vy/b3v8iIiKx5IILYOZMeP99b564rVu9DwJ/8Qvvw8AgUXEjIoFz7rkwbx48/TRUqeJ9\n7dbNO5iLiIhUVP/1X/Dhh96l3PHx3tQKl13mzRd34IDf0ZUNXZYmIoH2wQcweDBs2wYJCfD3v3uf\nVMmp6bK04ilPiUgQ7N7tzRn33HNel9HERHj8cW+i7Fi4VE2XpYlIhdShg9cu+he/gD17oE8fCIXg\n+HG/IxMREfFPvXre9AmrV3u58rvvYPhw6N4d1q3zO7qzp+JGIiMU8jsCqcDq1YOlSyE93VtOT4df\n/cr71MovoYyQfzuX6KFjo4j4rH17WLnSa8BTvz68+653X06HB0Ls2+d3dGdOl6VJZJhpZkWJCsuW\nwbBhXmHTqBEsWAAdO0Y+Dks3XFr0jgldlla8Ms1TOjaKSBTJyfE61P/1r5D3sFH/r45Jk+DGGyEu\nyk6JlMtlaWYWb2Zvm9kmM3vLzOoUsU0jM1tuZuvN7DMz+31p9ikiUhpXXQVr1kCnTvDtt16jgaef\n1t+XIiIi8fHw1FNengTYtQt+8xsvZ370kb+xlVRpa7AxwD+ccxcBy4GxRWxzDLjHOXcp0Am4zcwu\nLuV+RUTOWqNGkJEBd94JR496raKHDg1OpxgREZHSaNPG+/r3v3uNBlavhiuvhJtv9veS7pIobXFz\nNZA/HRCzgGsKb+Cc2+GcW5v//CDwOdCwlPsVESmVqlXhySfhxRe9yT5feME7cK9f73dkIiIi0WHY\nMNi0Ce67DypV8jqrtWzptZKO1sY8pS1uznfO7QSviAHOP9XGZtYE+BnwQSn3KyJSJgYP9k61t27t\nTWR25ZXeJ1UiIiLizR335z/DZ595l3bn5MBtt3mNCN5/3+/oTnba4sbMlpnZugKPz/K/Dihi82Kv\nWjezWsBLwJ35Z3CkIklL8zsCkWJddBGsWgUjRkBuLtxwA9xyCxw+XH77TOuuMSHo2BglmjRpgpnp\nEdBHkyZN/P4Vi0mF89TFF8Nbb8HChdC4MaxdC126wMiRXhvpaFGqbmlm9jmQ6pzbaWYXAO8451oV\nsV1l4HXgDefc1NO8pksrcLBPTU0lNTX1rGMUESkp57xT7nfcAUeOeJ9KLVgAFSEvZmRkkJGREV5O\nT09Xt7RiKE8Fj3ldl/wOQ8qJ3t+yl5sLkyZ5Z3QOH/bO7oRCXv6sUqV89lnSPFXa4mYysNc5N9nM\nHgDinXNjitjuf4Hdzrl7SvCaagUtIr76+GO49lr45huvc8zs2dC3r99RRZYV02JTlKeCSH/8Bpve\n3/KzZQvcfTe89pq3fMklXre1Xr3Kf9/F5anS3nMzGbjKzDYBvYBJ+TtLNLPX8593BoYDPc1sjZl9\nYma/LOV+RUTKzRVXeAVO377etcX9+sG4cdF786SIiIgfmjWDxYvh9de95xs2wM9/DkOGwNat/sSk\nSTxFRIqRlweTJ3uFTV4e9OwJ8+bB+adsnRIMOnNTPOWp4NEn+8Gm9zcyfvwRpkyBiRPh0CGoUQMe\negjuvReqVSv7/ZXXmRsRkcCKi4OxY+Ef//AKmuXLoV276OwOIyLB1bRpU5YvX+53GGFZWVnExcWR\nl5fndyj06NGDGTNm+B2GANWrw4MPep1HBw/2Cpxx47xupEuWRC4OFTcSGaGQ3xGInLUePbzZmrt0\ngexsSE2Fv/zFa0BwtkIZobIKT2KZjo0So8x0YjfozjZPNW7szSH3z3969+B8+aV3eXf//vDtt2Ub\nY1FU3EhkpKf7HYFIqSQleWdu/vAHOHYM7rnH+2Tq++/P7vXSMzUmBB0bJeod182GFVZp81TPnl67\n6ClTvG5qK1aUXye1glTciIiUUJUq8Pjj8PLLULu297V9e1i3zu/IRKSiOHLkCHfddRcNGzakUaNG\n3H333Rw9ehTw2pK/8sorALz//vvExcXxxhtvALB8+XLatWsXfp0ZM2ZwySWXkJCQQJ8+fdha4O7v\nuLg4nn32WVq2bEnLli2LjMM5x/Tp02nYsCENGzbkiSeeKFGMs2bNomvXrie8VlxcHF999RUAo0eP\n5vbbb6dfv37Url2bTp068fXXX4e3XbZsGa1atSI+Pp477rhD99JEuSpVvG5qmzd7Uys0aFD++1Rx\nIyJyhn79a6+bWps28MUX0LEjzJrld1QiUhFMnDiR1atXs27dOj799FNWr17NxIkTAejevXt4HpAV\nK1bQrFkzVqxYAUBmZmZ4PqZXX32VSZMmsWjRInbt2kXXrl0ZOnToCft59dVX+fDDD9mwYUOxsWRk\nZLBlyxbeeustJk+eHL4v6FQxwsmXtBVefuGFF0hPT2ffvn00a9aMhx56CIA9e/YwaNAgHnnkEXbv\n3k2zZs14XzdBxoQLLoDevSOzLxU3IiJnoXlzWLUKRo/2bpq88Ua4+WavW4yISHmZO3cuaWlpJCQk\nkJCQQFpaGrNnzwa84iYzMxPwipuxY8eGlzMzM+nevTsA06ZNY+zYsbRs2ZK4uDjGjBnD2rVr2bZt\nW3g/Dz74IHXq1KHaKdpchUIhqlevTuvWrRk9ejTz5s07bYxFKXz2ZeDAgVxxxRXExcUxfPhw1q5d\nC8DSpUtp3bo1AwcOpFKlStx1111ccMEFZ/pfKAGn4kZE5CzVqAEzZsDzz3ttLp97Djp3hgJXUIhI\nUIRCYHbyo7imEIW3L6PmEdnZ2TRu3Di8nJKSQnZ2NgCdOnVi8+bN/Pvf/+bTTz9l5MiRbNu2jT17\n9rB69epwcZOVlcWdd95J3bp1qVu3LgkJCZgZ27dvD79uo0aNThmHmZ2wTcE4ThVjSRQsWGrWrMnB\ngwfDr5ucnHzCtoWXRVTcSGSkpfkdgUi5+c1v4F//ggsvhE8+gcsv/89szcVJ664xIejYGEtCIa9F\nYuHHqYqbkmx3hpKSksjKygovZ2VlkZSUBECNGjW44oormDp1Kq1bt6Zy5cp06tSJKVOm0Lx5c+Lj\n4wFo3Lgx06ZNY+/evezdu5ecnBwOHjxIx44dw69bkm5oBc/0bN26NRzHqWI855xzyM3NDa/bsWNH\niX/2xMTEE+4NKhyDlK1YzVMqbiQy1O5UAq5dO/joIxgwAPbt876OHet1VitKKDUU0fgkSunYKGdo\n6NChTJw4kd27d7N7924mTJjAiBEjwuu7devGM888Ez5Lk5qaesIywG9/+1seeeSR8P00+/fv56WX\nXjqjOJxzTJgwgUOHDrF+/Xr+9re/cf311582xrZt27J+/XrWrVvH4cOHSU9PL3Fb6b59+7JhwwYW\nLVrE8ePHmTp1Kjt37jyjuKXkYjVPqbgRESkj8fGwaBFMnuxNADppkncDpXKviJRGwT/+x40bR/v2\n7WnTpg1t27alffv24Rvuwbvv5uDBg3Tr1u2E5YLFzTXXXMOYMWO4/vrrOe+882jTpg1vvvlmkfs7\nVUzdu3enefPmXHXVVdx///306tXrtDG2aNGChx9+mF69etGyZcuTOqedSkJCAgsWLOCBBx6gXr16\nbNmyhc6dO5f4+6VisGhroWdmLtpiEhE5U5mZcN11XmGTmOhNaNali99RlZyZ4ZzTLH1FUJ4Knvzf\nd7/DkHKi9zeYistTOnMjIlIOuneHNWugWzf47jtITYUnnvAuvRcREZHyoeJGRKScJCbCP/8J990H\nx4/DvffCtdfC/v1+RyYiIhJMKm4kMnTTrFRQlSvDn/8Mr7wCtWvDwoXQvj3c8kLI79AkGujYKCJR\nKpQR8juEs6J7biQyzHQ9jlR4X37pnbn59FMgZMxs4hg1yu+oiqZ7bopXpnlKx8aooHsygk3v79mx\ndMOlRe//m+65ERHxWfPm3nw4o0d7yzfeCDffDD/+6GtYIiIigaHiRkQkgmrUgBkzvOfVqsFzz0Hn\nzvD11/7GJSIiEgQqbkREfPKvf8GFF8Inn8Dll8Nrr/kdkYiISGxTcSMi4pN27eDjj2HAANi3z/s6\ndiwcO+Z3ZCIiIrFJxY1ERlqa3xGIRJW07t6YOO88r5PapEkQF+d97d3bm/xTKgAdGyXG9ejRgxk/\nXWvro1mzZtG1a1e/wwiUn/JUrFFxI5GhdqciJwilhsLP4+LggQe8OXEaNIB33vHO6rz3nn/xSYTo\n2ChSZszU4LEsFcxTsUTFjYhIlEhNhTVroGtX+O47b3nKFHUKFpHokJeX53cIIqel4kZEJIokJnpn\ncO69F44fhz/8AQYPhu+/9zsyEfFLXFwcX331VXh59OjRPPzwwwBkZmaSnJzMo48+Sv369bnwwguZ\nO3fuCdvecsst9O7dm9q1a9OjRw+2bt0aXr9x40Z69+5NQkICrVq1YsGCBSd876233krfvn0599xz\nycjIKDK+L7/8kg4dOlCnTh0GDhzIvn37wusWL15M69atqVu3Lj179mTjxo1n9HNNmTKFBg0a0LBh\nQ2bOnBnedu/evQwYMIA6derQsWNHtmzZcob/qxJUKm5ERKJMlSrw2GOwcCHUrg0vvwzt28Nnn/kd\nmYj44XSXW+3YsYO9e/eSnZ3NzJkzufnmm/niiy/C6+fOnUtaWhp79uyhbdu2DB8+HIDc3Fx69+7N\nDTfcwO7du5k/fz633nrrCQXIvHnzGD9+PAcOHKBLly5F7n/27NnMnDmTHTt2UKlSJe644w4ANm/e\nzLBhw3jqqafYtWsXffr0oX///hzL75pSkp/rwIEDZGdn8/zzz3Pbbbexf/9+AG699VZq1qzJzp07\nmT59elTc9yPRQcWNiEiUGjgQPvoI2rSBL76ADh1g9my/oxKRSHOnuTbVzJgwYQJVqlShW7du9O3b\nlxdffDG8vm/fvnTu3JkqVarwpz/9iVWrVrF9+3Zef/11mjZtysiRIzEz2rZty6BBg044e3P11VfT\nsWNHAKpWrVrk/keMGEGrVq2oUaMGEyZMYMGCBTjnePHFF+nXrx89e/akUqVK3HvvvRw6dIiVK1eW\n6OeqWrUq48ePp1KlSvTp04datWqxadMm8vLyWLhwIRMmTKB69epceumljBo1qkT/lxJ8Km4kMnTT\nrMgJQhmhEm3XogWsWgU33giHDsHIkfC738GPP5ZreBIpOjbGjFBGCEu3kx7FjeXC25d0zJ+N+Ph4\nqlevHl5OSUkhOzs7vJycnBx+fs455xAfH092djZZWVmsWrWKunXrUrduXeLj45k7dy47C7RrLPi9\nxSm4TUpKCkePHmX37t1kZ2eTkpISXmdmJCcns3379hL9XAkJCcTF/edP1Zo1a3Lw4EF27drF8ePH\nadSo0Qn7lbJVnr+z5UnFjURGerrfEYhElfTMko+JGjVgxgx47jmoVg2mTYMuXeCbb8ovPokQHRtj\nRig1hEtzJz2K6yhVePvSdJ6qWbMmubm54eUdO3acsD4nJ4dDhw6Fl7du3UpSUlJ4edu2beHnBw8e\nJCcnh6SkJJKTk0lNTWXv3r3s3buXnJwcvv/+e5555pnw9iXpQFbw9bOysqhSpQr16tUjKSmJrKys\nk7b9qSg53c9VnPr161OpUqUT9lvwPiIpG2eSp6KJihsRkRhgBv/937ByJTRt6k3+efnlsHSp35GJ\nSHlr164dc+fOJS8vjzfffJPMzMwT1jvnSEtL4+jRo7z77rssWbKEIUOGhNcvXbqUlStXcuTIEcaP\nH0/Hjh1p2LAh/fr1Y/PmzcyZM4djx45x9OhRPvroIzZt2nRG8c2ZM4eNGzeSm5tLWloagwcPxswY\nMmQIS5Ys4Z133uHYsWM8/vjjVK9enU6dOpXo5ypOXFwcgwYNIhQKcejQITZs2MCsWbPOKGYJLhU3\nIiIx5PLLvcKmf3/IyYG+fWH8eK+zmogE05NPPsnixYuJj49n3rx5DBw48IT1iYmJxMfHk5SUxIgR\nI5g2bRotWrQIrx82bBihUIiEhATWrFnDnDlzAKhVqxZvv/028+fPJykpiaSkJMaMGcPhw4dLHJuZ\nMWLECEaNGkVSUhJHjhxh6tSpALRs2ZI5c+Zw++23U79+fZYsWcJrr71G5cqVS/RzFbWvnzz99NMc\nOHCAxMREbrrpJm666aYSxyzBZqe7mSvSzMxFW0xSBsw0WYdIAZZuuLSzHxN5eTB5Mowb5z3/+c9h\n/nxISCij+MxwzmlGvCKUaZ7SsTEq5P+++x3GWcnMzGTEiBHFXpY1evRokpOT+eMf/xjhyKJHLL+/\nfiptnipvxeUpnbkREYlBcXEwdiwsWwb168P27d79OCIiIhVZZb8DkAoiLc3vCESiSlr3shkTPXvC\nmjXwww9Qq1aZvKREko6NUs5K0hBApChllaciTZeliYjISXRZWvGUp4JHly0Fm97fYNJlaSIiIiIi\nEmgqbkREREREJBBU3IiIiIiISCCouBERERERkUBQcSOREQr5HYFIVAllhPwOQaKBjo1RISUlBTPT\nI6CPlJQUv3/FYlKs5qlSdUszs3jgBSAF+AYY4pzbX2ibasAKoCpe6+mXnHPpp3hNdaEJItNEdSIF\nxerkaKJJPEWkYojVPFXaMzdjgH845y4ClgNjC2/gnDsM9HDOtQN+BvQxsytLud8SycjIiMRupAQy\n/A5AwjQuosTXfgcg0SDD7wAkTMfG6KH3IkrEaJ4qbXFzNTAr//ks4JqiNnLO5eY/rYZ39iYiZaAG\nR/TI8DsACdO4iBLf+B2ARIMMvwOQMB0bo4feiyjxjd8BnJ3SFjfnO+d2AjjndgDnF7WRmcWZ2Rpg\nB7DMOfdhKfcrIiIiIiJygsqn28DMlgENCv4T3pmXcUVsXuQZGedcHtDOzGoDi8zsEufchrOIV0RE\nREREpEilbSjwOZDqnNtpZhcA7zjnWp3me8YDPzjnphSzPnrvXBIRqUDUUKBoylMiItGhqDx12jM3\np7EYuBGYDIwCXi28gZnVA4465/abWQ3gKmDSmQQpIiISLZSnRESiV2nP3NQFXgSSgSy8VtD7zCwR\neM4518/MLsNrNhCX/3jBOfen0ocuIiIiIiLyH6UqbkRERERERKJFabulRQUz+6WZbTSzzWb2QBHr\nu5vZPjP7JP9RVDMEKSUzm25mO81s3Sm2ecrMvjCztWb2s0jGV5Gc7r3QmIgMM2tkZsvNbL2ZfWZm\nvy9mO42LgFOeig7KU9FDeSo6BDFPlfaeG9+ZWRzwDNALyAY+NLNXnXMbC226wjk3IOIBVix/A54G\n/reolWbWB2jmnGthZh2A/wE6RjC+iuSU70U+jYnydwy4xzm31sxqAR+b2dsFj08aF8GnPBVVlKei\nh/JUdAhcngrCmZsrgS+cc1nOuaPAfLzJRQvTDaDlzDn3HpBzik2uJv8g5pz7AKhjZg1Osb2cpRK8\nF6AxUe6cczucc2vznx8EPgcaFtpM4yL4lKeihPJU9FCeig5BzFNBKG4aAtsKLH/LyW8KQKf8U2lL\nzOySyIQmhRR+r7ZT9HslkaExEUFm1gT4GfBBoVUaF8GnPBU7NB6ji8ZEBAUlT8X8ZWkl9DHQ2DmX\nm39qbRHQ0ueYRPykMRFB+af6XwLuzP9kTKQwjUmRE2lMRFCQ8lQQztxsBxoXWG6U/29hzrmDzrnc\n/OdvAFXy21hLZG3Haxv+k5PeK4kMjYnIMbPKeAljtnPupLnA0LioCJSnYofGY5TQmIicoOWpIBQ3\nHwLNzSzFzKoC1+NNLhpW8LpAM7sSrwX23siGWWEYxV8juxgYCWBmHYF9zrmdkQqsAir2vdCYiKgZ\nwAbn3NRi1mtcBJ/yVHRRnooeylPRIVB5KuYvS3POHTez24G38Yq16c65z83st95q9/+Aa83sFuAo\ncAi4zr+Ig8vM5gKpQIKZbQXSgKrkvw/OuaVm9isz+xL4ARjtX7TBdrr3Ao2JiDCzzsBw4DMzWwM4\n4EEgBY2LCkN5KnooT0UP5anoEMQ8pUk8RUREREQkEIJwWZqIiIiIiIiKGxERERERCQYVNyIiIiIi\nEggqbkREREREJBBU3IiIiIiISCCouBERERERkUBQcSMiIiIiIoGg4kZERERERAJBxY2IiIiIiASC\nihuRMmBm/9fMdprZJjO7LP/fEszs137HJiIiojwlFYWKG5FSMrMBwHnAZcD1wJj8Vf2BxX7FJSIi\nAspTUrFU9jsAkQA45Jx7LP/5v81svpm1AI445475GZiIiAjKU1KBmHPO7xhEAsXMGgITgXucczl+\nxyMiIlKQ8pQEmS5LEyl7O4DaShgiIhKllKcksFTciJS9usDLfgchIiJSDOUpCSwVNyJl71fAGr+D\nEBERKYbylASWihuRsvd/nHOf+x2EiIhIMZSnJLBU3IiUve1+ByAiInIKylMSWOqWJiIiIiIigaAz\nNyIiIiIiEggqbkREREREJBBU3IiIiIiISCCouBERERERkUBQcSMiIiIiIoGg4kZERERERAJBxY2I\niIiIiASCihsREREREQmE/w983OpD6O5QRwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7f94e4430b00>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "numb_grid = 100\n",
    "gamma_grid = np.linspace(0.4, 2.1, numb_grid)\n",
    "\n",
    "value_ex1_grid = np.asarray([N1.zerosum(gamma = gamma_grid[i])[0] for i in range(numb_grid)])\n",
    "value_ex2_grid = np.asarray([N2.zerosum(gamma = gamma_grid[i])[0] for i in range(numb_grid)])\n",
    "\n",
    "fig, ax = plt.subplots(1, 2, figsize = (14, 5), sharey = True)\n",
    "fig.suptitle(r'The function $V(M(\\gamma))$', fontsize = 16)\n",
    "\n",
    "ax[0].plot(gamma_grid, value_ex1_grid, lw = 2)\n",
    "ax[0].set_title(r'Example 1', fontsize = 15)\n",
    "ax[0].axhline(0, color = 'k', lw =1)\n",
    "ax[0].set_xlabel(r'$\\gamma$', fontsize = 14)\n",
    "ax[0].axvline(N1.bounds()[0], color = 'r', lw = 1, linestyle = '--', label = 'lower bound')\n",
    "ax[0].axvline(N1.bounds()[1], color = 'g', lw = 1, linestyle = '--', label = 'upper bound')\n",
    "ax[0].legend(loc = 'best')\n",
    "\n",
    "ax[1].plot(gamma_grid, value_ex2_grid, lw = 2)\n",
    "ax[1].set_title(r'Example 2', fontsize = 15)\n",
    "ax[1].axhline(0, color = 'k', lw =1)\n",
    "ax[1].set_xlabel(r'$\\gamma$', fontsize = 14)\n",
    "ax[1].axvline(N2.bounds()[0], color = 'r', lw = 1, linestyle = '--', label = 'lower bound')\n",
    "ax[1].axvline(N2.bounds()[1], color = 'g', lw = 1, linestyle = '--', label = 'upper bound')\n",
    "ax[1].legend(loc = 'best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `expansion` method implements the bisection algorithm for $\\alpha_0$ (and uses the primal LP problem for $x_0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha_0 = 1.2599210478365421\n",
      "x_0 = [ 0.33  0.26  0.41]\n",
      "The corresponding p from the dual = [ 0.41  0.33  0.26  0.  ]\n"
     ]
    }
   ],
   "source": [
    "alpha0, x, p = N1.expansion()\n",
    "print('alpha_0 = {}'.format(alpha0)) \n",
    "print('x_0 = {}'.format(x))\n",
    "print('The corresponding p from the dual = {}'.format(p))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `interest` method implements the bisection algorithm  for $\\beta_0$ (and uses the dual LP problem for $p_0$)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta_0 = 1.2599210478365421\n",
      "p_0 = [ 0.41  0.33  0.26  0.  ]\n",
      "The corresponding x from the primal = [ 0.33  0.26  0.41]\n"
     ]
    }
   ],
   "source": [
    "beta0, x, p = N1.interest()\n",
    "print('beta_0 = {}'.format(beta0)) \n",
    "print('p_0 = {}'.format(p))\n",
    "print('The corresponding x from the primal = {}'.format(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Of course, when $\\gamma^*$ is unique, it is irrelevant which one of the two methods we use. In particular, as will be shown below, in case of an irreducible $(A,B)$ (like in Example 1), the maximal and minimal roots of $V(M(\\gamma))$ necessarily coincide implying a ''full duality'' result, i.e. $\\alpha_0 = \\beta_0 = \\gamma^*$, and that the expansion (and interest) rate $\\gamma^*$ is unique. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uniqueness and irreducibility\n",
    "\n",
    "As an illustration, compute first the maximal and minimal roots of $V(M(\\cdot))$ for Example 2, which displays a reducible input-output pair $(A, B)$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "alpha_0 = 1.2528658034279943\n",
      "x_0 = [ 0.    0.    0.33  0.26  0.41]\n",
      "The corresponding p from the dual = [ 0.56  0.44  0.    0.    0.    0.  ]\n"
     ]
    }
   ],
   "source": [
    "alpha0, x, p = N2.expansion()\n",
    "print('alpha_0 = {}'.format(alpha0)) \n",
    "print('x_0 = {}'.format(x))\n",
    "print('The corresponding p from the dual = {}'.format(p))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "beta_0 = 1.0000000009313226\n",
      "p_0 = [ 0.5  0.5  0.   0.   0.   0. ]\n",
      "The corresponding x from the primal = [ 0.33  0.33  0.33  0.    0.  ]\n"
     ]
    }
   ],
   "source": [
    "beta0, x, p = N2.interest()\n",
    "print('beta_0 = {}'.format(beta0)) \n",
    "print('p_0 = {}'.format(p))\n",
    "print('The corresponding x from the primal = {}'.format(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As we can see, with a reducible $(A,B)$, the roots found by the bisection alhorithms might differ, so there might be multiple $\\gamma^*$ that make the value of the game with $M(\\gamma^*)$ zero. (see the figure above)\n",
    "\n",
    "Indeed, although the von Neumann theorem assures existence of the equilibrium, Assumptions I and II are not sufficient for uniqueness. Nonetheless, Kemeny et al. (1967) show that there are at most finitely many economic solutions, meaning that there are  only finitely many $\\gamma^*$ that satisfy $V(M(\\gamma^*)) = 0$ and $x_0^TBp_0 > 0$ and that for each such $\\gamma^*_i$, there is a self-sufficient part of the economy (a sub-economy) that in equilibrium can expand independently  with the expansion coefficient $\\gamma^*_i$. \n",
    "\n",
    "The following theorem (see Theorem 9.10. in Gale, 1960) asserts that imposing irreducibility is sufficient for uniqueness of $(\\gamma^*, x_0, p_0)$.\n",
    "\n",
    "**Theorem II:** Consider the conditions of Theorem 1. If the economy $(A,B)$ is irreducible, then $\\gamma^*=\\alpha_0=\\beta_0$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "\n",
    "## A special case \n",
    "\n",
    "There is a special $(A,B)$ that allows us to simplify the solution method significantly by invoking the powerful Perron-Frobenius theorem for nonnegative matrices.\n",
    "\n",
    "**Definition:** We call an economy _simple_ if it satisfies \n",
    " 1. $n=m$\n",
    " 2. Each activity produces exactly one good\n",
    " 3. Each good is produced by one and only one activity\n",
    "\n",
    "These assumptions imply that $B=I_n$, i.e., that $B$ can be written as an identity matrix (possibly after reshuffling its rows and columns).\n",
    "\n",
    "The simple model has the following special property (Theorem 9.11. in Gale): if $x_0$ and $\\alpha_0>0$ solve the TEP with $(A,I_n)$, then\n",
    "\n",
    "$$x_0^T = \\alpha_0 x_0^T A\\hspace{1cm}\\Leftrightarrow\\hspace{1cm}x_0^T A=\\left(\\frac{1}{\\alpha_0}\\right)x_0^T$$\n",
    "\n",
    "The latter shows that $1/\\alpha_0$ is a positive eigenvalue of $A$ and $x_0$ is the correponding nonnegative left eigenvector. The classical result of **Perron and Frobenius** implies that a nonnegative matrix always has a nonnegative eigenvalue-eigenvector pair. Moreover, if $A$ is irreducible, then the optimal intensity vector $x_0$ is positive and _unique_ up to multiplication by a positive scalar.\n",
    "\n",
    "Suppose that $A$ is reducible with $k$ irreducible subsets $S_1,\\dots,S_k$. Let $A_i$ be the submatrix corresponding to $S_i$ and let $\\alpha_i$ and $\\beta_i$ be the associated expansion and interest factors, respectively. Then we have \n",
    "\n",
    "$$\\alpha_0 = \\max_i \\{\\alpha_i\\}\\hspace{1cm}\\text{and}\\hspace{1cm}\\beta_0 = \\min_i \\{\\beta_i\\}$$\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### References\n",
    "\n",
    "Gale, David, **The Theory of Linear Economic Models**, New York: McGraw-Hill Book Company, 1960.\n",
    "\n",
    "Hamburger, Michael J.,  Gerald L. Thompson and Roman L. Weil, Jr., 1967, \"Computation of Expansion Rates for the Generalized von Neumann Model of an Expanding Economy\", Econometrica Vol. 35, No. 3/4,  pp. 542-547\n",
    "\n",
    "Kemeny, John G.; Morgenstern, Oskar; Thompson, Gerald L., 1956, \"A generalization of the von Neumann model of an expanding economy\". Econometrica. 24. pp. 115–135.\n",
    "\n",
    "Nash, John. 1951. \"Non-Cooperative Games.\" Annals of Mathematics, Second Series, 54, no. 2, pp 286-95.\n",
    "\n",
    "von Neumann, John, 1928, \"Zur theories der gesellschaftsspiele\". Mathematische Annalen 100, pp 295–320. English translation by S. Bergmann in Contributions to the Theory of Games IV, ed. R.D. Luce and A.W. Tucker. Princeton: Princeton University Press, 1959\n",
    "\n",
    "von Neumann, John, 1937, \"Über ein ökonomisches Gleichuns-system und eine Verallgeneinerung des Brouwerschen Fixpunksatzes,\" Ergebnisse eines mathematischen Kolloquium, 8, , pp. 73-83, translated as \"A Model of General Economic Equilibrium,\" Review of Economic Studies, 13. 1945, pp. 1-9.\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}