43_policy_iter.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "# Foundations of Computational Economics #43\n",
    "\n",
    "by Fedor Iskhakov, ANU\n",
    "\n",
    "<img src=\"_static/img/dag3logo.png\" style=\"width:256px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "fragment"
    }
   },
   "source": [
    "## Solving DP problems with policy iterations\n",
    "\n",
    "<img src=\"_static/img/lab.png\" style=\"width:64px;\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "subslide"
    }
   },
   "source": [
    "<img src=\"_static/img/youtube.png\" style=\"width:65px;\">\n",
    "\n",
    "[https://youtu.be/m9R_io8mQoI](https://youtu.be/m9R_io8mQoI)\n",
    "\n",
    "Description: Policy iterations solution method for infinite horizon dynamic models. Solving stochastic inventory management problem with policy iterations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Policy iterations\n",
    "\n",
    "Also known as **Howard policy improvement** algorithm\n",
    "\n",
    "Area of application: dynamic decision problems with\n",
    "\n",
    "- infinite horizon, discrete time  \n",
    "- particularly well suited for problems with finite state space (discrete states)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Main idea of policy iterations\n",
    "\n",
    "Break the search of the fixed point (of Bellman or Coleman-Reffett operators) into two steps:\n",
    "\n",
    "1. Policy evaluation: compute value function for fixed policy function  \n",
    "1. Policy improvement : find best policy for given value function  \n",
    "\n",
    "\n",
    "Iterative (back and forth) approach instead of jointly finding the value and policy functions in the fixed point problem."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### General form of Bellman equation\n",
    "\n",
    "$$\n",
    "V(\\text{state}) = \\max_{\\text{decisions}} \\big[ U(\\text{state},\\text{decision}) + \\beta \\mathbb{E}\\big\\{ V(\\text{next state})  \\big| \\text{state},\\text{decision} \\big\\} \\big]\n",
    "$$\n",
    "\n",
    "- $ V(\\text{state}) $ is **value function** = maximum attainable (discounted) expected reward/utility/payoff  \n",
    "- $ U(\\text{state},\\text{decision}) $ is per-period/flow/instantaneous reward/utility/payoff  \n",
    "- (*next state*) is the *stochastic* next period state resulting from current state and decision  \n",
    "- expectation $ \\mathbb{E}\\{\\cdot\\} $ is taken over the distribution of the next period state conditional on current state and decision  \n",
    "- $ \\beta $ is a discount factor to measure future rewards in terms of current ones  \n",
    "\n",
    "\n",
    "The optimal choices are revealed along the solution of the Bellman equation as decisions which solve the maximization problem in the right hand side!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Policy iterations step 1\n",
    "\n",
    "For fixed policy function (decisions) solve for function $ V(\\text{state}) $:\n",
    "\n",
    "$$\n",
    "V(\\text{state}) = U(\\text{state},\\text{decision}) + \\beta \\mathbb{E}\\big\\{ V(\\text{next state})  \\big| \\text{state},\\text{decision} \\big\\}\n",
    "$$\n",
    "\n",
    "- functional equation, although simpler than Bellman equation  \n",
    "- becomes non-linear system of equations with discrete or discretized states  \n",
    "- becomes *linear* system of equations in many applications where $ \\mathbb{E}\\{\\cdot\\} $ can be expressed as matrix multiplication\n",
    "  - Markov discrete choice problems (all randomness in transition probabilities, no interpolation)\n",
    "  - when quadrature integration and linear interpolation are used in discretized state spaces  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Policy iterations step 2\n",
    "\n",
    "For fixed value function $ V_0(\\text{state}) $ find the optimal policy while computing:\n",
    "\n",
    "$$\n",
    "V(\\text{state}) = \\max_{\\text{decisions}} \\big[ U(\\text{state},\\text{decision}) + \\beta \\mathbb{E}\\big\\{ V_0(\\text{next state})  \\big| \\text{state},\\text{decision} \\big\\} \\big]\n",
    "$$\n",
    "\n",
    "- standard evaluation of the Bellman operator  \n",
    "- any previously considered approaches for implementation of the Bellman operator are applicable  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Why policy iterations are useful?\n",
    "\n",
    "*Rate of convergence!*\n",
    "\n",
    "- known to be considerably faster than VFI and time iterations  \n",
    "- yet, step 1 subproblem is comparable in complexity to VFI  \n",
    "- so, efficient implementation is necessary to justify the use  \n",
    "- *exact theoretical properties on complexity of the algorithm are unknown for the general case*  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Which problems can be solved efficiently by policy iterations?\n",
    "\n",
    "- want to simplify the policy evaluation step as much as possible!  (otherwise, functional equation of the same complexity as solved by VFI)  \n",
    "- best is to present the Bellman equation as a linear system (vectorization, linear algebra solver)  \n",
    "- discrete choice, discrete finite state space problems are best  \n",
    "- with other problems, it’s not completely straightforward, but not impossible  \n",
    "\n",
    "\n",
    "Question: can linear interpolation of the value function on a fixed grid be represented as matrix multiplication?  Gaussian quadrature?\n",
    "3-dim array algebra as in stochastic consumption-savings problem with discretized choices?"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### General policy iterations algorithm\n",
    "\n",
    "1. Initialize policy or value function (in the latter case reverse the order of the next two steps)  \n",
    "1. Compute the value of the current policy by solving the Bellman equation without max operation (assuming that current policy is optimal)  \n",
    "1. Re-compute the policy function by applying the Bellman operator to the found value function  \n",
    "1. Repeat until convergence in policy or value function space  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Inventory dynamics model\n",
    "\n",
    "Recall the inventory management problem (video 37 and exercise 11)\n",
    "\n",
    "- discrete state, discrete choice, infinite horizon in discrete time  \n",
    "- choice: $ q\\ge 0 $ is the order of new inventory  \n",
    "- state: $ x\\ge 0 $ is inventory in the current time period  \n",
    "- state: $ d\\ge 0 $ is stochastic demand in the current time period  \n",
    "- payoff in each period is given by  \n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "\\pi & = & p \\cdot \\text{sales} - r \\cdot \\text{storage} - c [\\text{if order is made}] \\\\\n",
    "& = & p \\min\\{x,d\\} - r (x - \\min\\{x,d\\} + q) - c \\mathbb{1}\\{q>0\\} \\\\\n",
    "& = & p \\min\\{x,d\\} - r (\\max\\{x-d,0\\} + q) - c \\mathbb{1}\\{q>0\\}\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Inventory dynamics model\n",
    "\n",
    "- transition rule is given by  \n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "x' &=& \\max\\{x-d,0\\} + q \\\\\n",
    "d' &-& \\text{idiosyncratic}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "Parameters in the model:\n",
    "\n",
    "- $ p $ is the profit per one unit of (supplied) good  \n",
    "- $ c $ is the fixed cost of ordering any amount of new inventory  \n",
    "- $ r $ is the cost of storing one unit of good  \n",
    "- $ d $ follows truncated geometric distribution with the support $ i \\in \\{0,1,\\dots,N\\} $ and corresponding probabilities $ pr_i = (1-\\lambda)^i \\lambda $ (with last probability corrected for truncation)  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Bellman equation\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "V(x,d) &=& \\max_{q \\ge 0} \\Big\\{ \\pi + \\beta \\mathbb{E}\\Big[ V\\big(x', d' \\big) \\Big| x,d,q \\Big] \\Big\\} \\\\\n",
    "&=& \\max_{q \\ge 0} \\Big\\{ sp - x' r - c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta \\mathbb{E}\\Big[ V\\big( x', d' \\big) \\Big] \\Big\\}\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "s &=& \\min\\{x,d\\} \\\\\n",
    "x' &=& \\max\\{x-d,0\\} + q\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "- the expectation is taken over the distribution of the next period demand $ d' $  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Bellman equation in expected value function space\n",
    "\n",
    "- random demand is *idiosyncratic* = distribution does not depend on the previous period variables  \n",
    "- can reduce the dimensionality of the fixed point problem by rewriting the Bellman operator in expected value function terms  \n",
    "\n",
    "\n",
    "$$\n",
    "EV(x') =  \\mathbb{E}\\Big[ V\\big(x', d' \\big) \\Big| x,d,q \\Big] =  \\mathbb{E}\\Big[ V\\big(x', d' \\big) \\Big]\n",
    "$$\n",
    "\n",
    "$$\n",
    "V(x,d) = \\max_{q \\ge 0} \\Big\\{ sp - x' r - c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta EV(x') \\Big\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Bellman equation in expected value function space\n",
    "\n",
    "Taking the expectation with respect to $ d $ on both sides, we get\n",
    "\n",
    "$$\n",
    "EV(x) = \\mathbb{E}\\Big[ \\max_{q \\ge 0} \\Big\\{ sp - x' r - c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta EV(x') \\Big\\} \\Big]\n",
    "$$\n",
    "\n",
    "$$\n",
    "EV(x) = \\sum_{d} \\Big[ \\max_{q \\ge 0} \\Big\\{ sp - x' r - c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta EV(x') \\Big\\} \\Big] pr(d)\n",
    "$$\n",
    "\n",
    "- demand is a discrete random variable with probabilities $ pr(d) $ for different values $ d $  \n",
    "- functional equation in $ EV(\\cdot) $ which is also a contraction mapping!  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Policy function\n",
    "\n",
    "- fixed point of the Bellman operator $ EV(x) $ in *expected value function* space defines the solution of the model  \n",
    "- but the policy is still given by the function of two variables  \n",
    "\n",
    "\n",
    "$$\n",
    "q^\\star(x,d) = \\arg\\max_{q \\ge 0} \\Big\\{ sp - x' r - c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta EV(x') \\Big] \\Big\\}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "s &=& \\min\\{x,d\\} \\\\\n",
    "x' &=& \\max\\{x-d,0\\} + q\n",
    "\\end{eqnarray}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Further refinements\n",
    "\n",
    "Sequence of events in each time period (timing assumptions in the model):\n",
    "\n",
    "- start period with inventory $ x $  \n",
    "- demand $ d $ is realized  \n",
    "- decision $ q^\\star(x,d) $  \n",
    "- current profit is computed as function of $ x,d,q $  \n",
    "- the next period inventory $ x' =  \\max\\{x-d,0\\} + q $ is stored for the rest of the period  \n",
    "\n",
    "\n",
    "Note that ordering decision is made based on the stock *remaining after trading is over* $ \\max\\{x-d,0\\} $ rather than on $ x $ and $ d $ independently!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Post trade stock (remaining stock)\n",
    "\n",
    "- let $ y = \\max\\{x-d,0\\} = x - \\min\\{x,d\\} $ denote the remaining inventory after the trading phase  \n",
    "- now the optimal policy can be written as  \n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "q^\\star(x,d) &=& \\arg\\max_{q \\ge 0} \\Big\\{ sp - x'r - c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta EV(x') \\Big\\}\\\\\n",
    "q^\\star(x,d) &=& \\arg\\max_{q \\ge 0} \\Big\\{ p \\min\\{x,d\\} - (y + q)r  - c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta EV(y+q) \\Big\\}\\\\\n",
    "q^\\star(x,d) &=& \\arg\\max_{q \\ge 0} \\Big\\{ - q r - c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta EV(y+q) \\Big\\} = q^\\star(y)\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "- thus, policy function is function of *one* variable $ q^\\star(y) $ which is convenient  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Resulting compact solution for the inventory model\n",
    "\n",
    "- Bellman equation in expected value function space  \n",
    "\n",
    "\n",
    "$$\n",
    "EV(x) = \\sum_{d} \\Big[ p \\min\\{x,d\\} - yr + \\max_{q \\ge 0} \\Big\\{ -qr -c \\mathbb{1}\\{q>0\\}\n",
    "+ \\beta EV(y+q) \\Big\\} \\Big] pr(d)\n",
    "$$\n",
    "\n",
    "- Definition of the policy function, assuming $ EV(x) $ satisfies the equation above  \n",
    "\n",
    "\n",
    "$$\n",
    "q^\\star(y) = \\arg\\max_{q \\ge 0} \\Big\\{ -qr - c \\mathbb{1}\\{q>0\\} + \\beta EV(y+q) \\Big\\}\n",
    "$$\n",
    "\n",
    "$$\n",
    "y = \\max\\{x-d,0\\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Preparing for policy iterations\n",
    "\n",
    "- Policy $ q(y) $ evaluation (in terms of expected value function) is $ EV(x) $ that solves  \n",
    "\n",
    "\n",
    "$$\n",
    "EV(x) = \\sum_{d} \\Big[ p \\min\\{x,d\\} - (y+q(y))r -c \\mathbb{1}\\{q(y)>0\\}\n",
    "+ \\beta EV(y+q(y)) \\Big] pr(d)\n",
    "$$\n",
    "\n",
    "- Improved policy $ q'(y) $ given expected value $ EV(x) $ is  \n",
    "\n",
    "\n",
    "$$\n",
    "q'(y) = \\arg\\max_{q \\ge 0} \\Big\\{ -qr - c \\mathbb{1}\\{q>0\\} + \\beta EV(y+q) \\Big\\}\n",
    "$$\n",
    "\n",
    "Policy iterations solver is applicable!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Examine the problem in step 1 (policy evaluation)\n",
    "\n",
    "$$\n",
    "EV(x) = \\sum_{d} \\Big[ p \\min\\{x,d\\} - (y+q(y))r -c \\mathbb{1}\\{q(y)>0\\}\n",
    "+ \\beta EV(y+q(y)) \\Big] pr(d)\n",
    "$$\n",
    "\n",
    "- state space $ S = \\{0,1,\\dots,N\\} $, $ x,d,y,q \\in S $ (all function are defined on finite $ S $)  \n",
    "- $ EV(x) $ is given by a vector $ EV $ with $ N+1 $ elements corresponding to points in $ S $  \n",
    "- $ EV(y+q) $ is a particular indexing of the elements of $ EV $, can be expressed as matrix multiplication  \n",
    "- otherwise RHS is linear function of the re-indexed vector $ EV $  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "EV=\n",
      "[11 21 32 43 54 65 76]\n",
      "y+q(y)=\n",
      "[1, 0, 2, 2, 4, 3, 5]\n",
      "re-indexing matrix =\n",
      "[[0 1 0 0 0 0 0]\n",
      " [1 0 0 0 0 0 0]\n",
      " [0 0 1 0 0 0 0]\n",
      " [0 0 1 0 0 0 0]\n",
      " [0 0 0 0 1 0 0]\n",
      " [0 0 0 1 0 0 0]\n",
      " [0 0 0 0 0 1 0]]\n",
      "EV(y+q(y))\n",
      "[21 11 32 32 54 43 65]\n"
     ]
    }
   ],
   "source": [
    "# example of re-indexing as matrix multiplication\n",
    "import numpy as np\n",
    "ev = np.array([11,21,32,43,54,65,76])  # vector EV\n",
    "j = [1,0,2,2,4,3,5]                    # represent y+q(y)\n",
    "n = ev.size\n",
    "M = np.zeros((n,n),dtype=int)  # start with all zeros\n",
    "M[range(n),j] = 1              # mark\n",
    "ev1 = M@ev\n",
    "print('EV=',ev,'y+q(y)=',j,'re-indexing matrix =',M,'EV(y+q(y))',ev1,sep='\\n')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Plan for the implementation\n",
    "\n",
    "- let $ EV $ be the vector of values of $ EV(x) $ for all $ x $  \n",
    "- let $ M(d) $ denote the matrix such that $ M(d) \\cdot EV = EV(y+q(y)) $  \n",
    "- let $ C(d) = p \\min\\{x,d\\} - (y+q(y))r -c \\mathbb{1}\\{q(y)>0\\} $ denote the constant term  \n",
    "\n",
    "\n",
    "$$\n",
    "\\begin{eqnarray}\n",
    "EV &=& \\sum_{d} \\Big[ C(d) + \\beta M(d) \\cdot EV \\Big] pr(d) \\\\\n",
    "&=& \\sum_{d} pr(d) C(d)  + \\beta \\sum_{d} pr(d) M(d) \\cdot EV  = \\bar{C} + \\beta \\bar{M} \\cdot EV\n",
    "\\end{eqnarray}\n",
    "$$\n",
    "\n",
    "- where $ \\bar{M} = \\sum_{d} pr(d) M(d) $ is a square matrix  \n",
    "- *linear system* $ (I-\\beta\\bar{M}) \\cdot EV = \\bar{C} $!  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "plt.rcParams['figure.figsize'] = [12, 8]\n",
    "\n",
    "# COPY of the code from exercise 11, only infinite horizon version, combined into single class\n",
    "class inventory_model:\n",
    "    '''Small class to hold model fundamentals and its solution'''\n",
    "\n",
    "    def __init__(self,label='noname',\n",
    "                 max_inventory=10,  # upper bound on the state space\n",
    "                 c = 3.2,            # fixed cost of order\n",
    "                 p = 2.5,           # profit per unit of good\n",
    "                 r = 0.5,           # storage cost per unit of good\n",
    "                 β = 0.95,          # discount factor\n",
    "                 dp = 0.5,          # parameter in geometric distribution of demand\n",
    "                 demand = 4         # fixed demand\n",
    "                 ):\n",
    "        '''Create model with default parameters'''\n",
    "        self.label=label # label for the model instance\n",
    "        self.c, self.p, self.r, self.β, self.dp= c, p, r, β, dp\n",
    "        self.demand = demand\n",
    "        # created dependent attributes (it would be better to have them updated when underlying parameters change)\n",
    "        self.n = max_inventory+1    # number of inventory levels\n",
    "        self.upper = int(max_inventory)  # upper boundary on inventory\n",
    "        self.x = np.arange(self.n,dtype=int)  # all possible values of inventory and demand (state space)\n",
    "\n",
    "    def __repr__(self):\n",
    "        '''String representation of the model'''\n",
    "        return 'Inventory model labeled \"{}\"\\nParamters (c,p,r,β) = ({},{},{},{})\\nDemand={}\\nUpper bound on inventory {}' \\\n",
    "               .format (self.label,self.c,self.p,self.r,self.β,self.demand,self.upper)\n",
    "\n",
    "    def sales(self,x,d):\n",
    "        '''Sales in given period'''\n",
    "        return np.minimum(x,d)\n",
    "\n",
    "    def next_x(self,x,d,q):\n",
    "        '''Inventory to be stored, becomes next period state'''\n",
    "        return x - self.sales(x,d) + q\n",
    "\n",
    "    def profit(self,x,d,q):\n",
    "        '''Profit in given period'''\n",
    "        return self.p * self.sales(x,d) - self.r * self.next_x(x,d,q) - self.c * (q>0)\n",
    "\n",
    "    def demand_pr(self,plot=False):\n",
    "        '''Computes stochastic demand probs'''\n",
    "        k = np.arange(self.n)  # all possible values of demand\n",
    "        pr = (1-self.dp)**k *self.dp\n",
    "        pr[-1] = 1 - pr[:-1].sum()  # update last prob to ensure sum=1\n",
    "        if plot:\n",
    "            plt.step(self.x,pr,where='mid')\n",
    "            plt.title('Distribution of demand')\n",
    "            plt.show()\n",
    "        return pr\n",
    "\n",
    "    def bellman(self,ev0):\n",
    "        '''Bellman operator for inventory model\n",
    "           Inputs: model object\n",
    "                   next period EXPECTED value function\n",
    "        '''\n",
    "        pr = self.demand_pr()\n",
    "        ev1 = np.zeros(shape=self.x.shape)\n",
    "        for j,d in enumerate(self.x):  # over all values of demand\n",
    "            # create the grid of choices (same as x), column-vector\n",
    "            q = self.x[:,np.newaxis]\n",
    "            # compute current period profit (relying on numpy broadcasting to get the matrix with choices in rows)\n",
    "            p = self.profit(self.x,d,q)\n",
    "            # indexes for next period value with extrapolation using last value\n",
    "            i = np.minimum(self.next_x(self.x,d,q),self.upper)\n",
    "            # compute the Bellman maximand\n",
    "            vm = p + self.β*ev0[i]\n",
    "            # find max and argmax\n",
    "            v1 = np.amax(vm,axis=0)  # maximum in every column\n",
    "            ev1 = ev1 + pr[j]*v1\n",
    "        q1 = self.optimal_policy(ev1)\n",
    "        return ev1, q1\n",
    "\n",
    "    def optimal_policy(self,ev):\n",
    "        '''Computes the optimal policy function as function of post trade stock for the stochastic\n",
    "        inventory dynamics model for given EV function'''\n",
    "        # idea: 2-dim array with q in axes 0, y = max(x-d,0) in axis 1\n",
    "        q = self.x[:,np.newaxis]  # choices\n",
    "        y = self.x[np.newaxis,:]  # post trading stock\n",
    "        # indexes for next period value with extrapolation using last value\n",
    "        i = np.minimum(y+q,self.upper)\n",
    "        # compute the Bellman maximand\n",
    "        vm = -self.r*q -self.c*(q>0) + self.β*ev[i]\n",
    "        # find argmax and argmax\n",
    "        return np.argmax(vm,axis=0)  # maximum in every column\n",
    "\n",
    "    def solve_vfi(self,tol=1e-6,maxiter=500,callback=None):\n",
    "        '''Solves the model using value function iterations\n",
    "        '''\n",
    "        ev0 = np.zeros(self.n) # initial point for VFI\n",
    "        for i in range(maxiter):  # main loop\n",
    "            ev1, q1 = self.bellman(ev0)  # update approximation\n",
    "            err = np.amax(np.abs(ev0-ev1))\n",
    "            if callback != None: callback(iter=i,err=err,ev1=ev1,ev0=ev0,q1=q1,model=self)\n",
    "            if err<tol:\n",
    "                break  # break out if converged\n",
    "            ev0 = ev1  # get ready to the next iteration\n",
    "        else:\n",
    "            raise RuntimeError('Failed to converge in %d iterations'%maxiter)\n",
    "        return ev1, q1\n",
    "\n",
    "    def policy_evaluation(self,q):\n",
    "        '''The expected value of the given policy\n",
    "           Returns a vector ev for all points in the state space\n",
    "        '''\n",
    "        pr = self.demand_pr()\n",
    "        Cbar = np.zeros(self.n)\n",
    "        Mbar = np.zeros((self.n,self.n))\n",
    "        for j,d in enumerate(self.x):          # over all values of demand\n",
    "            y = self.x - self.sales(self.x,d)  # post trade stock for all values of x | d\n",
    "            p = self.profit(self.x,d,q[y])     # profits for all values of x | d, q()\n",
    "            Cbar += pr[j]*p\n",
    "            i = np.minimum(y+q[y],self.upper)  # next period stock for all x | d, q()\n",
    "            M = np.zeros((self.n,self.n),dtype=int)  # start with all zeros\n",
    "            M[range(self.n),i] = 1             # transformation matrix\n",
    "            Mbar += pr[j]*M\n",
    "        A = np.eye(self.n) - self.β*Mbar\n",
    "        ev = np.linalg.solve(A,Cbar)\n",
    "        return ev\n",
    "\n",
    "    def solve_policyiter(self,tol=1e-6,maxiter=500,callback=None):\n",
    "        '''Solves the model using policy iterations\n",
    "        '''\n",
    "        q0 = np.zeros(self.n,dtype=int) # initial point for policy iterations\n",
    "        for i in range(maxiter):  # main loop\n",
    "            # step 1: policy evaluation\n",
    "            ev = self.policy_evaluation(q0)\n",
    "            # step 2: policy improvement\n",
    "            q1 = self.optimal_policy(ev)\n",
    "            err = np.amax(np.abs(q0-q1))\n",
    "            if callback != None: callback(iter=i,err=err,ev1=ev,q1=q1,model=self)\n",
    "            if err<tol:\n",
    "                break  # break out if converged\n",
    "            q0 = q1  # get ready to the next iteration\n",
    "        else:\n",
    "            raise RuntimeError('Failed to converge in %d iterations'%maxiter)\n",
    "        ev1 = self.policy_evaluation(q1)\n",
    "        return ev1,q1\n",
    "\n",
    "    def solve_show(self,method='vfi',maxiter=1000,tol=1e-6,**kvargs):\n",
    "        '''Illustrate solution'''\n",
    "        fig1, (ax1,ax2) = plt.subplots(1,2,figsize=(14,8))\n",
    "        ax1.grid(b=True, which='both', color='0.65', linestyle='-')\n",
    "        ax2.grid(b=True, which='both', color='0.65', linestyle='-')\n",
    "        ax1.set_xlabel('Inventory')\n",
    "        ax2.set_xlabel('Post trade inventory')\n",
    "        ax1.set_title('Value function')\n",
    "        ax2.set_title('Optimal new orders')\n",
    "        def callback(**argvars):\n",
    "            mod, ev, q = argvars['model'],argvars['ev1'],argvars['q1']\n",
    "            ax1.step(mod.x,ev,color='k',alpha=0.25,where='mid')\n",
    "            ax2.step(mod.x,q,color='k',alpha=0.25,where='mid')\n",
    "            print(argvars['iter'],end=' ')\n",
    "        if method=='vfi':\n",
    "            ev,pk = self.solve_vfi(maxiter=maxiter,tol=tol,callback=callback,**kvargs)\n",
    "        elif method=='policyiter':\n",
    "            ev,pk = self.solve_policyiter(maxiter=maxiter,tol=tol,callback=callback,**kvargs)\n",
    "        else:\n",
    "            raise RuntimeError('Unknown method in solve_show()')\n",
    "        # add solutions\n",
    "        ax1.step(self.x,ev,color='r',linewidth=2.5,where='mid')\n",
    "        ax2.step(self.x,pk,color='r',linewidth=2.5,where='mid')\n",
    "        plt.show()\n",
    "        return ev,pk"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAskAAAHiCAYAAAAatlGFAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAdwElEQVR4nO3df9jdd13f8dfb1MKUHys2Okh70yJlo/4KMzTXLhXcJWBRIXpdQEqGK4rrUNFNZBMBW+x0Iv6YbjJNlU7BYUV0GncVKwhMd1mwKfLDhgGh1DYm0kKQ8Esg5b0/zgmefriTnDT3yX3nzuNxXblyzvfn5863p9cz33zOOdXdAQAA/sEXrPYAAABgrRHJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQycC6VVW/UlU/tkLHWqqqj1XVhunzN1XV96zEsafHe21VXb5SxzuB8/5EVX2wqv52zu27qh6+6HGdjKq6YDrOs1Z7LMDpy/9AgNNSVd2W5MuSHE5yd5I9SV6R5Jru/mySdPezT+BY39Pdrz/aNt19e5L7ndyoP3e+Fyd5eHc/Y+b4T1yJY5/gOM5P8sNJHtrdd57q8wOsZe4kA6ezJ3X3/ZM8NMlLkvxIkpev9EnW8R3Jhyb5kEAG+HwiGTjtdfdHuntXku1JLq+qr0ySqvr1qvqJ6eNzq+p/V9XfVdXBqvqzqvqCqnplkqUkfzidTvEfZ/65/llVdXuSNxzln/C/vKr+oqo+UlV/UFUPmp7rG6tq3+wYq+q2qnpcVV2a5AVJtk/P9/bp+s9N35iO60VV9ddVdWdVvaKqHjhdd2Qcl1fV7dOpEi882p9NVT1wuv9d0+O9aHr8xyV5XZKHTMfx60fZ/z9U1YGq2l9V3z2su09V/ex0HB+YTm/5R7N/BtM/zzunx/j2qvqWqnrP9Bq8YOZYl1TVjdPrc6Cqfqmqzp5Z31X17Kp6b1V9uKpeVlU1XbdhOo4PVtWtSb71aH8eAPMSycC60d1/kWRfkm9YZvUPT9dtzGSaxgsmu/R3Jrk9k7vS9+vul87s89gkj0zyzUc55b9O8t1JHpLJtI//OscY/yjJf07y29Pzfc0ymz1z+utfJnlYJtM8fmnY5uuT/NMk35Tkyqp65FFO+d+SPHB6nMdOx/xd06klT0yyfzqOZ447ToP+eUken+SiJI8bNvnpJI9IsjnJw5NsSnLlzPp/kuS+M8t/NckzknxtJtfoyqp62HTbu5P8UJJzk/yL6c/1fcP5vi3Jo5N8TZKn5R+uy7+ZrntUki1JnnKUPwuAuYlkYL3Zn+RByyz/TJIHZzL/9jPd/Wfd3cc51ou7++Pd/cmjrH9ld/9Vd388yY8ledqRN/adpH+V5Oe7+9bu/liSH01y2XAX+8e7+5Pd/fYkb88kHO9hOpbtSX60uz/a3bcl+bkk3znnOJ6W5H/M/Iwvnjl2ZRKnP9TdB7v7o5nE/2Uz+38myU9292eSXJdJAP/idCy3JLklyVcnSXff3N1v7u7D03HuzCTqZ72ku/9uOj/8jZnE+ZFx/kJ339HdB5P81Jw/H8BRiWRgvdmU5OAyy38myd4kf1xVt1bV8+c41h0nsP6vk3xhJiF4sh4yPd7ssc/K5A74EbOfRvGJLP+mwnOTnL3MsTadwDjGn/GIjUm+KMnN0ykSf5fkj6bLj/hQd989fXzkLxofmFn/ySPjrqpHTKfD/G1VHcokuMc/y6P9zMcaJ8C9IpKBdaOqHp1JAP7fcd307uUPd/fDkjwpyXOr6puOrD7KIY93p/n8mcdLmdw5/WCSj2cSkEfGtSH3jMfjHXd/Jm+qmz324dwzMOfxwemYxmP9zZz7H8jn/4yzx/5kkq/o7n88/fXA7r63nwDyy0n+X5KLuvsBmUyHqRUYJ8C9IpKB015VPaCqvi2Tf9L/ze5+5zLbfFtVPXw6TeBQJnNgj9zl/EAmc3ZP1DOq6uKq+qIkVyd5zfTO6XuS3LeqvrWqvjDJi5LcZ2a/DyS5oKqO9v/g30ryQ1V1YVXdL/8wh/nwiQxuOpZXJ/nJqrp/VT00yXOT/Oach3h1kmfO/IxXzRz7s5nMMf4vVfWlSVJVm6rqaPO3j+f+mVyXj1XVP0vyvSew76uT/GBVnVdV5ySZ518JAI5JJAOnsz+sqo9m8k/tL0zy80m+6yjbXpTk9Uk+luTGJP+9u980XfdTSV40nTbwvBM4/yuT/Hom0wDum+QHk8mnbWTyprNfy+Su7cczedPgEb8z/f1DVfXWZY577fTYf5rk/Un+PskPnMC4Zv3A9Py3ZnKH/VXT4x9Xd782yS8keUMmU1XeMGzyI9Plb55OkXh9Jm8mvDeel2RHko9mEt+/fQL7/mqSGzKZm/3WJL93L8cA8Dl1/PetAADAmcWdZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYnHX8TU6tc889ty+44ILVHgYAAOvczTff/MHu3rjcujUXyRdccEF279692sMAAGCdq6qjfo296RYAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMJgrkqvq0qp6d1XtrarnL7P+uVW1p6reUVV/UlUPnVl3d1W9bfpr10oOHgAAFuGs421QVRuSvCzJ45PsS3JTVe3q7j0zm/1lki3d/Ymq+t4kL02yfbruk929eYXHDQAACzPPneRLkuzt7lu7+9NJrkuybXaD7n5jd39i+vTNSc5b2WECAMCpc9w7yUk2Jblj5vm+JFuPsf2zkrx25vl9q2p3ksNJXtLdv3/CozwFfvwPb8me/YdOaJ9tmzdlx9alBY0IAIDVMk8k1zLLetkNq56RZEuSx84sXuru/VX1sCRvqKp3dvf7hv2uSHJFkiwtnR7RuefAJKhFMgDA+jNPJO9Lcv7M8/OS7B83qqrHJXlhksd296eOLO/u/dPfb62qNyV5VJJ7RHJ3X5PkmiTZsmXLsgG+aFc96StOaPvtO29c0EgAAFht88xJvinJRVV1YVWdneSyJPf4lIqqelSSnUme3N13ziw/p6ruM318bpKvSzL7hj8AAFhzjnsnubsPV9VzktyQZEOSa7v7lqq6Osnu7t6V5GeS3C/J71RVktze3U9O8sgkO6vqs5kE+UuGT8UAAIA1Z57pFunu65NcPyy7cubx446y358n+aqTGSAAAJxqvnEPAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGc0VyVV1aVe+uqr1V9fxl1j+3qvZU1Tuq6k+q6qEz6y6vqvdOf12+koMHAIBFOG4kV9WGJC9L8sQkFyd5elVdPGz2l0m2dPdXJ3lNkpdO931QkquSbE1ySZKrquqclRs+AACsvHnuJF+SZG9339rdn05yXZJtsxt09xu7+xPTp29Oct708TcneV13H+zuDyd5XZJLV2boAACwGPNE8qYkd8w83zdddjTPSvLae7kvAACsurPm2KaWWdbLblj1jCRbkjz2RPatqiuSXJEkS0tLcwwJAAAWZ547yfuSnD/z/Lwk+8eNqupxSV6Y5Mnd/akT2be7r+nuLd29ZePGjfOOHQAAFmKeSL4pyUVVdWFVnZ3ksiS7Zjeoqkcl2ZlJIN85s+qGJE+oqnOmb9h7wnQZAACsWcedbtHdh6vqOZnE7YYk13b3LVV1dZLd3b0ryc8kuV+S36mqJLm9u5/c3Qer6j9lEtpJcnV3H1zITwIAACtknjnJ6e7rk1w/LLty5vHjjrHvtUmuvbcDBACAU8037gEAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDA4KzVHsDpbM+BQ9m+88aFHX/b5k3ZsXVpYccHAGB5Ivle2rZ500KPv+fAoSQRyQAAq0Ak30s7ti4tNGAXeYcaAIBjMycZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAZzRXJVXVpV766qvVX1/GXWP6aq3lpVh6vqKcO6u6vqbdNfu1Zq4AAAsChnHW+DqtqQ5GVJHp9kX5KbqmpXd++Z2ez2JM9M8rxlDvHJ7t68AmMFAIBT4riRnOSSJHu7+9YkqarrkmxL8rlI7u7bpus+u4AxAgDAKTXPdItNSe6Yeb5vumxe962q3VX15qr69hMaHQAArIJ57iTXMsv6BM6x1N37q+phSd5QVe/s7vfd4wRVVyS5IkmWlpZO4NAAALDy5rmTvC/J+TPPz0uyf94TdPf+6e+3JnlTkkcts8013b2lu7ds3Lhx3kMDAMBCzBPJNyW5qKourKqzk1yWZK5Pqaiqc6rqPtPH5yb5uszMZQYAgLXouJHc3YeTPCfJDUneleTV3X1LVV1dVU9Okqp6dFXtS/LUJDur6pbp7o9Msruq3p7kjUleMnwqBgAArDnzzElOd1+f5Pph2ZUzj2/KZBrGuN+fJ/mqkxwjAACcUr5xDwAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAZzfeMeq2PPgUPZvvPGhZ5j2+ZN2bF1aaHnAAA43YjkNWrb5k0LP8eeA4eSRCQDAAxE8hq1Y+vSwuN10XepAQBOV+YkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDA4KzVHgCra8+BQ9m+88aFHX/b5k3ZsXVpYccHAFgEkXwG27Z500KPv+fAoSQRyQDAaUckn8F2bF1aaMAu8g41AMAimZMMAAADkQwAAAORDAAAg7kiuaourap3V9Xeqnr+MusfU1VvrarDVfWUYd3lVfXe6a/LV2rgAACwKMeN5KrakORlSZ6Y5OIkT6+qi4fNbk/yzCSvGvZ9UJKrkmxNckmSq6rqnJMfNgAALM48d5IvSbK3u2/t7k8nuS7JttkNuvu27n5Hks8O+35zktd198Hu/nCS1yW5dAXGDQAACzNPJG9KcsfM833TZfM4mX0BAGBVzBPJtcyynvP4c+1bVVdU1e6q2n3XXXfNeWgAAFiMeSJ5X5LzZ56fl2T/nMefa9/uvqa7t3T3lo0bN855aAAAWIx5IvmmJBdV1YVVdXaSy5LsmvP4NyR5QlWdM33D3hOmywAAYM06biR39+Ekz8kkbt+V5NXdfUtVXV1VT06Sqnp0Ve1L8tQkO6vqlum+B5P8p0xC+6YkV0+XAQDAmnXWPBt19/VJrh+WXTnz+KZMplIst++1Sa49iTECAMAp5Rv3AABgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYHDWag+A9W3PgUPZvvPGhR1/2+ZN2bF1aWHHBwDOTCKZhdm2edNCj7/nwKEkEckAwIoTySzMjq1LCw3YRd6hBgDObOYkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMBAJAMAwEAkAwDAQCQDAMDgrNUeAJyMPQcOZfvOGxd6jm2bN2XH1qWFngMAWFtEMqetbZs3Lfwcew4cShKRDABnGJHMaWvH1qWFx+ui71IDAGuTOckAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwmCuSq+rSqnp3Ve2tqucvs/4+VfXb0/VvqaoLpssvqKpPVtXbpr9+ZWWHDwAAK++4X0tdVRuSvCzJ45PsS3JTVe3q7j0zmz0ryYe7++FVdVmSn06yfbrufd29eYXHDQAACzPPneRLkuzt7lu7+9NJrkuybdhmW5LfmD5+TZJvqqpauWECAMCpM08kb0pyx8zzfdNly27T3YeTfCTJl0zXXVhVf1lV/6eqvmG5E1TVFVW1u6p233XXXSf0AwAAwEqbJ5KXuyPcc25zIMlSdz8qyXOTvKqqHvB5G3Zf091bunvLxo0b5xgSAAAszjyRvC/J+TPPz0uy/2jbVNVZSR6Y5GB3f6q7P5Qk3X1zkvclecTJDhoAABZpnki+KclFVXVhVZ2d5LIku4ZtdiW5fPr4KUne0N1dVRunb/xLVT0syUVJbl2ZoQMAwGIc99MtuvtwVT0nyQ1JNiS5trtvqaqrk+zu7l1JXp7klVW1N8nBTEI6SR6T5OqqOpzk7iTP7u6Di/hBAABgpRw3kpOku69Pcv2w7MqZx3+f5KnL7Pe7SX73JMcIAACnlG/cAwCAgUgGAIDBXNMt4Ey258ChbN9548KOv23zpuzYurSw4wMAJ04kwzFs2zx+b87K2nPgUJKIZABYY0QyHMOOrUsLDdhF3qEGAO49c5IBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBgcNZqDwDOdHsOHMr2nTcu9BzbNm/Kjq1LCz0HAKwnIhlW0bbNmxZ+jj0HDiWJSAaAEyCSYRXt2Lq08Hhd9F1qAFiPzEkGAICBSAYAgIFIBgCAgUgGAICBSAYAgIFIBgCAgUgGAICBSAYAgIFIBgCAgUgGAICBr6WGM8CeA4cW+vXU2zZvWvjXawPAqSSSYZ3btnnTQo+/58ChJBHJAKwrIhnWuR1blxYasIu8Qw0Aq8WcZAAAGIhkAAAYiGQAABiIZAAAGIhkAAAYiGQAABiIZAAAGPicZOCk+UY/ANYbkQycFN/oB8B6JJKBk+Ib/QBYj8xJBgCAgUgGAICBSAYAgIFIBgCAgTfuAWveoj9iLvExcwDck0gG1rRFf8Rc4mPmAPh8IhlY0xb9EXOJj5kD4POZkwwAAAORDAAAA5EMAAADc5IBsvhP0PDpGQCnF5EMnPEW/QkaPj0D4PQjkoEz3qI/QcOnZwCcfsxJBgCAgTvJAKeAOc8ApxeRDLBg5jwDnH5EMsCCmfMMcPoRyQDrwKKncySmdABnFpEMcJpb9HSOxJQO4MxT3b3aY7iHLVu29O7du1d7GADM2L7zxuw5cCgXP/gBCzuHO9XAqVZVN3f3luXWuZMMwHF58yFwpnEnGYBV5041sBrcSQZgTVv0neq3vP9g3vL+g/mDt/3NQs8jxGH9cCcZgHXvVW+5feGB/Jb3H0ySbL3wQQs7hwiHlXWsO8lzRXJVXZrkF5NsSPJr3f2SYf19krwiydcm+VCS7d1923TdjyZ5VpK7k/xgd99wrHOJZABOR4sOcREOK++kIrmqNiR5T5LHJ9mX5KYkT+/uPTPbfF+Sr+7uZ1fVZUm+o7u3V9XFSX4rySVJHpLk9Uke0d13H+18IhkAPt96iPD1wF8kju3e/nd68UMekKue9BULGNGxneyc5EuS7O3uW6cHuy7JtiR7ZrbZluTF08evSfJLVVXT5dd196eSvL+q9k6P5+uhAOAELPqbG0/FlJTT3ama2346W09/2ZonkjcluWPm+b4kW4+2TXcfrqqPJPmS6fI3D/t+3rszquqKJFckydKSv50BwKm26AhfD/xF4vi2XvigdXO3fZ5IrmWWjXM0jrbNPPumu69Jck0ymW4xx5gAAE4pf5E4s3zBHNvsS3L+zPPzkuw/2jZVdVaSByY5OOe+AACwpswTyTcluaiqLqyqs5NclmTXsM2uJJdPHz8lyRt68o7AXUkuq6r7VNWFSS5K8hcrM3QAAFiM4063mM4xfk6SGzL5CLhru/uWqro6ye7u3pXk5UleOX1j3sFMQjrT7V6dyZv8Dif5/mN9sgUAAKwFvkwEAIAz0rE+Am6e6RYAAHBGEckAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADAQyQAAMBDJAAAwEMkAADCo7l7tMdxDVd2V5K9X6fTnJvngKp2bU8d1Xv9c4zOD63xmcJ3PDKt1nR/a3RuXW7HmInk1VdXu7t6y2uNgsVzn9c81PjO4zmcG1/nMsBavs+kWAAAwEMkAADAQyfd0zWoPgFPCdV7/XOMzg+t8ZnCdzwxr7jqbkwwAAAN3kgEAYCCSk1TVpVX17qraW1XPX+3xsBhVdVtVvbOq3lZVu1d7PKyMqrq2qu6sqr+aWfagqnpdVb13+vs5qzlGTt5RrvOLq+pvpq/pt1XVt6zmGDk5VXV+Vb2xqt5VVbdU1b+bLvd6XkeOcZ3X3Ov5jJ9uUVUbkrwnyeOT7EtyU5Knd/eeVR0YK66qbkuypbt93uY6UlWPSfKxJK/o7q+cLntpkoPd/ZLpX3zP6e4fWc1xcnKOcp1fnORj3f2zqzk2VkZVPTjJg7v7rVV1/yQ3J/n2JM+M1/O6cYzr/LSssdezO8nJJUn2dvet3f3pJNcl2bbKYwLm1N1/muTgsHhbkt+YPv6NTP4HzGnsKNeZdaS7D3T3W6ePP5rkXUk2xet5XTnGdV5zRPLkwtwx83xf1ujF4qR1kj+uqpur6orVHgwL9WXdfSCZ/A85yZeu8nhYnOdU1Tum0zH8M/w6UVUXJHlUkrfE63ndGq5zssZezyI5qWWWndlzUNavr+vuf57kiUm+f/rPt8Dp65eTfHmSzUkOJPm51R0OK6Gq7pfkd5P8++4+tNrjYTGWuc5r7vUskid3js+feX5ekv2rNBYWqLv3T3+/M8n/ymSqDevTB6bz3o7Mf7tzlcfDAnT3B7r77u7+bJJfjdf0aa+qvjCTcPqf3f1708Vez+vMctd5Lb6eRfLkjXoXVdWFVXV2ksuS7FrlMbHCquqLp28QSFV9cZInJPmrY+/FaWxXksunjy9P8gerOBYW5Eg4TX1HvKZPa1VVSV6e5F3d/fMzq7ye15GjXee1+Ho+4z/dIkmmHzPyC0k2JLm2u39ylYfECquqh2Vy9zhJzkryKtd5faiq30ryjUnOTfKBJFcl+f0kr06ylOT2JE/tbm/6Oo0d5Tp/Yyb/NNtJbkvyb4/MXeX0U1Vfn+TPkrwzyWeni1+QyXxVr+d14hjX+elZY69nkQwAAAPTLQAAYCCSAQBgIJIBAGAgkgEAYCCSAQBgIJIBAGAgkgEAYCCSAQBg8P8BxxmJnke2NcgAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 864x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 "
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 "
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 "
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy8AAAHwCAYAAACv2lDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde5zc913f+/dnd2Zndna12l3Z2viGRaI6NBHEaYVLE6BpQhGFXCgE6IWc0Kbo0PMoDgcoCZzSxgTa0Bv2OeVAlfLISUpoLtAUyIFsUsCieSg4wbEcO7GVs5sotrXRypF29qLV3r/nj5kxsvY3s7/vd2d+87u8no+HH7J25ze/7/y0l99nPpevOecEAAAAAGk30O8FAAAAAEAcBC8AAAAAMoHgBQAAAEAmELwAAAAAyASCFwAAAACZQPACAAAAIBMIXpB5ZnbEzJyZlXrw3C82s0fMbNnM7u3283c479eZ2YqZDSZ1TgCAn17+rDazd5jZb3b7efup+bv6aL/XgWwjeEHfmdm0mf1CxMffYGYXexGUePgZSQ865w445/7PXp3EzM6b2Xe0/u6ce8o5N+qc2+7VOQGgaMzsR8zsMTNbbf5++TUzG/c4np/VQJ8RvCAN/h9JbzIzu+Hjb5L0fufcVvJLes6dkj7fx/MDALrAzH5K0i9L+meSDkr6FjV+xn/CzIb6uba86fabjn1+ExMpQ/CCNPjvkiYlfVvrA2Y2Iem1kt7X/Pv3NMu3lszsaTN7R7snu/GdsRtT72b2LWZ2xszqZvaomb2qzfP8saS/Kek/NssC7jKzB83sH1/3mB8xs09e93dnZj9mZv+fmS2Y2a9eH5SZ2Y+a2RPNMrQvmNlfMbP/IunrJP1+8zw/c2MpnJndama/Z2ZXzGzGzH70htf3ITN7X/N5P29mx/e+7ABQDGY2Juk+ST/unPuYc27TOXde0g+qEcD8cPNx7zCz3zazDzZ/nn7WzF7W/Fycn9UPmtkvNn/HrJjZ75vZITN7f/P312fM7Mh163qg+TttycweNrNvUwxm9ioze8bMfsrMLpnZV83sH173+YqZ/Tsze8rM5s3s181suPm502b2/c3//9bm+r+7+ffvMLOzbc5ZMbP7zWyu+d/9Zla5YT1vM7OLkt7T/Pg/a65tzsz+UcTztVvjruczs5vM7KPN391XzOx/mhn3sQXEPzr6zjl3TdKHJP0v1334ByU96Zx7tPn3q83Pj0v6Hkn/xMy+1/dcZnabpP9X0i+qETD9tKTfMbObI9b1akn/U9I/bZYFfDHmaV4r6Zslvaz5Ok40z/0Dkt7RfB1jkl4v6bJz7k2SnpL0uuZ5/k3Ec/5XSc9IulXSGyX9KzN7zXWff72kD6hxfX5P0n+MuVYAKIJXSKpK+m/Xf9A5tyLpDyX9res+/AZJH1bjd8RvSfrvZlaO+bNakv6uGpUDt0l6kaRPqXEzPynpCUn/8rrHfkbS3ded68NmVo35ml6gRgbpNklvkfSr1njjT2pkmO5qPvfR5mP+RfNzpyW9qvn/3y7pS5L+xnV/P93mfP+HGtmqu9X4/XaPpH9+w3om1QgGT5rZd6nxO/ZvSfpLkr5Dz9dpjbueT9JPqfF78GZJU5J+TpJrs1bkGMEL0uK9kn6g9a6LGjf472190jn3oHPuMefcjnPuc2rczP+NiOfZyw9L+gPn3B80n+sTkv5c0nfvc/3Xe5dzru6ce0rSn6jxg1mS/rGkf+Oc+4xrmHHOfWWvJzOzOyR9q6S3OefWnHNnJf1nNX45tnyy+Zq2Jf0XNX6xAAAabpL0tTZlyF9tfr7lYefcbzvnNiX9BzWCnm/xONd7nHOzzrlFNQKjWefc/2ie+8OSXt56oHPuN51zl51zW865fy+pIunFMc+zKekXmlmkP5C0IunFzWz/j0r6351zV5xzy5L+lRpBldQITq4PVv71dX//G2ofvPyD5vkuOeeeVSOTdf3voR1J/9I5t958U/IHm9ficefcVTXevJMkxVhj1PNtSrpF0p3N1/w/nXMELwVE8IJUcM59UtKzkt5gZi9UI3PxW63Pm9lfM7M/MbNnzWxR0o/p+b9s4rpTjSCp3vpPjcDglv2/iudcvO7/VyWNNv//DkmzAc93q6TWD/eWr6jxLlW7c1aNGmEAaPmapJva/Fy8pfn5lqdb/+Oc29FfZL3jmr/u/69F/L31O0HNsq8nzGyx+fvooOL/brt8QzDW+n1zs6SapIev+z33sebHpUYm6C4zm1LjzbX3SbrDzG5SI5vyp23Od6sav3tavqLnX5dnnXNrNzz+6Rse37LXGqOe799KmpH0cTP7kpm9vc06kXMEL0iT96mRcXmTpI87567/gf9bapRD3eGcOyjp1yXd2ODfclWNH4otL7ju/5+W9F+cc+PX/TfinHtXzDV2eu69PK1GCUGUTu8ezUmaNLMD133s6yRd8Dg3ABTZpyStS/q+6z9oZiOS/rakP7ruw3dc9/kBSber8XNY6mKZUrO/5W1qZCgmnHPjkhbV/ndbXF9TI0h66XW/5w4650YlyTm3KulhSW+V9LhzbkPSGUk/qUaW6GttnndOjTcAW75Of3FdpN3X5qu67lo2Hx9rjVHP55xbds79lHPuhZJeJ+knbyifRkEQvCBN3qdGTeyP6rqSsaYDamQf1szsHkl/v8PznJX0d82s3Gxcf+N1n/tNSa8zsxNmNmhm1WZj4O0x13hW0veZWc0as+rfEvM4qVHq9dNm9let4aiZtX4RzEt6YdRBzrmn1fjF8q+b6/2m5nnf73FuACisZgnXfZL+LzP7rubvhyNqlHE9o0a5bctfNbPva2ZpfkKNoOfPmp9r+7M6wAFJW2pUHZTM7F+o0Q+5L81s0bsl/YqZHZYa/Z5mduK6h52W9E/1FyViD97w9yj/VdI/N7Obm1maf6HG79R2PiTpR8zsJWZW03W9PjHX+Dxm9trm702TtCRpu/kfCobgBanRnPxyRtKIGlmW6/1vkn7BzJbV+IH5oQ5P9fNqZDgW1Phl9Vz5WTMQeIMajX7PqpEN+WeK/73wK5I21PgF9l55BBDOuQ9L+qXmepb1F1PWpEbN8T9vps9/OuLwvyfpiBrvcn1EjTrgT8Q9NwAUXbPB/uck/Ts1bn4fUuN3wGucc+vXPfR3Jf2QGr9D3iTp+5r9L9LeP6t9TKvRE/NFNUqq1vT8Mqv9eJsaJVZ/ZmZLkv6Hnt9Lc1qN4OlP2/w9yi+q0SP6OUmPSfps82ORnHN/KOl+SX/cXMsfe67xRn+p+ZgVNTJp/7dz7sEOj0dOGb1OAAAAjVHJko46536432sBEI3MCwAAAIBMIHgBAAAAkAmUjQEAAADIBDIvAAAAADKB4AUAAABAJiS6A/dNN93kjhw5Enz88vKyDhw4sPcDc6zo16Dor1/iGhT99Uv7vwYPP/zw15xzN+/9yOLh99T+Ff0aFP31S1yDor9+qbe/pxINXo4cOaI///M/Dz5+enpaJ0603b+oEIp+DYr++iWuQdFfv7T/a2BmX+nicnKF31P7V/RrUPTXL3ENiv76pd7+nqJsDAAAAEAmELwAAAAAyASCFwAAAACZQPACAAAAIBMIXgAAAABkAsELAAAAgEwgeAEAAACQCQQvAAAAADKB4AUAAABAJhC8AAAAAMgEghcAAAAAmUDwAgAAACATCF4AAAAAZALBCwAAAIBMIHgBAOSSmb3YzM5e99+Smf1Ev9cFAAhX6vcCAADoBefcOUl3S5KZDUq6IOkjfV0UAGBfyLwAAIrgNZJmnXNf6fdCAADhyLygp+4x089LOtCl56tIerBLz7Us6XUPPSTdc0+XnhFAiv1dSf+1lye430x3q/3PKH7mAMD+Ebygp35e0uv6vYhO/tpf8wqGuPkAssfMhiS9XtLPtvn8SUknJWlqakrT09NB57lb0qv2eMyle+/VI/fdF/T8WbG4uBh8DfOg6K9f4hoU/fVLvb0GBC+ILSSL8m3NP+uSznqer1za/eW545wGzHZ9fHtnR3JOFvG5yHXt7Dz3/6/yXJfe+U7p93/f9ygA/fO3JX3WOTcf9Unn3ClJpyTp+PHj7sSJE2FneetbdeWP/1iTk5O7PlU/fVrjkg5Xqwp+/oyYnp7O/WvspOivX+IaFP31S729BgQviG0/WZTPjozon9x2W6zHLi0taWhoSLfffvuuzy0sLGhiYmLXx+v1ura2tjQ8PBzrHKNf+ILevrmpsYH4bV/ftLOjcUn1j35UZ2MGSRLZGiAF/p56XDImSbr/fn2mzS/ss2b+b5QAAHYheCmo/WRRFiU9NjgY65jt7W1dHRjQ9Dd9k1548GCsY65cuaKhoSHdcccduz5XLpd166237vr4wYMHtb6+rqGhoVjneKJe18mNDR0+fDjW4yXpnY8+qtdJGhfZGiArzKwm6W9J+l/7vRYAwP4RvBTUfrIoZ8fG9Nav//pYj7148aKq1apeceSIdudLopmZKpWKbr755l2fu3r1auTHWxmX0dHRWOdYXV3V2tqaXvSiF8VclfRvn3xStrGhcY9szbHtbbI1QB8551YlHer3OgAA3UHwUlCtjEtd0qMxb6qdc1qW9NC3fItecijevcDw8LBGR0f1spe9LPbaLl68KEl6wQtesOtzy8vLkZmXq1evSpJGRkZineOpp57S1tZWZCDUzqdvv10/evWqpqamYh9DtgYAAKB7CF5yIKQE7O7mn+eGh3UyojwryjPPPKNKpaI3fcM3KO4t//LyskoRjffdtr29Ledc7McfOHBAq6urXuc4dOiQxsbG9A3f8A2xj3nno482zudxnrvVCHa0vOyzPAAAgNwjeMmBfY0jHh3VC1/4wlgPvXbtmmq1WmTmo52vfe1riQQvgzF7cFoqlYoqlYrXazl8+LDW19e9sjW/NzKiN25uRmaR2nnvU0/R2AsAABCB4CUHri8B8xlHvCzpi69+te6KWQb17LPPasCj30NqNNJXKhUdOXIk9jEDAwMys8iG/Xq9HvlcV65c0dLSUuxzlMtl73V96Utf0urqqlfAc9NNN2l1dTVyclpbTz0lqTFaNapPpt1GnfTJAACAvCN4yZEnKhW9OWYgMjc3p1KppJMe/Rvlcjn2NK+WSqUSu4l+PyYnJyP3Vmjn4sWLWl9f1/nz52Mfs7Ky4p1FGh8fV61W05133hn7mOUzZxrHij4ZAACA68W6EzOzcUn/WdIxSU7SP5J0TtIHJR2RdF7SDzrnFnqySsQyUqvpJS95SazHrqyseJdNTU5OqlqtBmVRkjAxMRF79PHVq1f17LPPer2Wy5cvq16ve63p4MGD2tnZ8RpY8K4Pf1i2taWJNoFS1EadL93aok8GAADkXty3kR+Q9DHn3BvNbEhSTdLPSfoj59y7zOztkt4u6W09Wmdh7NV8H1Uy1Gq+L5XLkRs4RhkfH9fw8LDXzfvjjz8e+7H7MTExofHxcR09enTX5y5fvqwXv/jFuz7uG1SMjIxocnLS6/VfuXJFa2trXgHf6Oiotra2vNb2xIED+nubm23HOC8tLWlsbOx5H7v/0UfpkwEAALm3Z/BiZmOSvl3Sj0iSc25D0oaZvUF/UdXyXjXuqQle9mk/zfc7tVrsZvKxsbHYY4WTduDAgV0353vZ2NjQwsKC5ufnYz3+6aeflnNO6+vrsc9x9epVr4lmklStVrW2tuZ1zOTkpDY2Ntpm0ebm5nYHUM2pZu36ZNqhTwYAAGRJnMzLCyU9K+k9ZvYySQ9LequkKefcVyXJOfdVM4us1zGzk5JOStLU1JSmp6eDF7u4uLiv47Pg+ub7yF3snZNuuDnd3t7WsqRHjh3TcsyyITPT+vq6zp6N3+J/+fJlVatVXb58OfYxV65cUblc9jqmXq+3fR3Xrl3T008/Hfu52nn22WefW19cly9f1tDQkFfQ55yTc84rMzQwMKCBgYHn9q65Ua1W2/W51tUK6ZO5dO+9euS++zyP6p8i/BzYC9cAAFBUcYKXkqS/IunHnXMPmdkDapSIxeKcOyXplCQdP37cnThxImihkjQ9Pa39HJ8FDzb//NLYmN5+7Niuzy8sLOwqDTt79qzK5bJ+9lu/VbsLqqLNz8+rWq3q7rvv3vvBTRsbG6pWq3rFK14R+5jz589rfX09stSrndaN/vj4+K7PnTlzxuv87czMzKher+vAgfg7sDz55JM6cOCA1/kvXLigtbU1r+v86U9/WhsbG23Pc+7cuV3X81999KPSzo4mPEZGH9ve1rikw9Vqpr6vivBzYC9cAwBAUcUJXp6R9Ixz7qHm339bjeBl3sxuaWZdbpF0qVeLLKJqtRo5oapcLu8qGZqdnfUe+xvav7K2tqZz587FfvxTTz0V1LAft3cn1NjYmMbGxmI3+EuNQGR5ednr9fuWjEmNf3vfqWZfGB3VD21u6q677op9DH0yAAAga/a8Q3LOXTSzp83sxc65c5JeI+kLzf/eLOldzT9/t6crzaC9mu+jtN6fL5fLkf0rV69e3fXxAwcOqFKpBK8zrrGxMR04cMAriyLJO7BKgu/IZ0kaHh5OZHJaSJ/M2NiYNjc3Y0+bk/RcnwwAAEBWxH1798clvb85aexLkv6hpAFJHzKzt0h6StIP9GaJ2bWf5ns7cCByqtXy8vKuj4+OjqpcLns9/8jIiEZGRrwCEZ/+kJZKpeLVFJ+UjY0N72NuueUW7xK4VrbGR61WU61WaxvwRW3UefDgQa2vr8ce2PC856PJHwAAZESs4MU5d1bS8YhPvaa7y8mX65vvPxdzZ/qdnR0tS7rymviXNqTMaGRkRNVq1esY38dL0vr6eiJZIV8hmZeQ1zI8PKzh4WGvgOfTn/5022b9dmq1WmRJYSf7afJnM0wAANAPfne8CDIzMqJ7I/YsifL5z39epVJJ77jjjtjPXyqVgm7GfftX6vW6dwCTp8xLyGsJCfhCjqlUKhr0aNaXpH89OCjb3tYhj4DsL6+vsxkmAADoG4KXBNRGRmL3Inz5y19WqVTyKhl6wQteoGq16vXu/vnz51Uul4P6V3wUPfOytrbmHSSura2pUqm0/beJ2qjz8OHDWlpa8lrb50dG9ENbW/r2b//22Me87WMfo8kfAAD0DcFLAto130cZHh72fge9Wq1634yXy2Vtbm56HROSeSh65iWkT+bJJ5/0XZrGx8c79slEGR4e1ubmZs8nuwEAAHQLwUvKDA0NeTffS42bcZ9395eXl732OJHCMg9Fz7yEHLOxsdExW/P000/v+lxIMDYyMqLNzU2a/AEAQGYQvMS0n7HHlaGh2I3UBw8e7NjgHVUydObMGZVKpVRODrvttts0OjrqtUnj3NycLl3q7bZBSWVeQo45dOiQNjc3vb4GPv3pT2tra8vra+DQoUNaWVnxWhtN/gAAoJ8IXmLaz9jjwYmJ2OU84+Pj3pPDSqWStra2vI5JanLY6uqqVldXdfbs2djH1Ot1SY1r0Stpzrxsbm56Z9+2tra8v24mJiZUqVS8JpS9s/nneMzpeZL0jTs7NPkDAICuIHiJ6Xljj2OWy+w41xh7/Nf/ui6ePx/rmFYQ4lMydPXqVY2MjMR6/paQnd9DsgghQcL4+LgOHz4ceVM9Pz/vlcVpJyS7k1TmJaQfqVQq6erVq16lg6VSSWNjY159Mo9Vq/q+rS299KUvjX3M/Y8+SpM/AADoCoIXTz5jjx9//HENDg7qfd/2bToS8/lbjdc+JUOPPvqo987vSWVeOgUi7bQCi6jgYnZ21iuL04nvuiR5l1kllXkJ+fcMUavVtLm5qVtuuSX+QY8+2rsFAQCAQiF48eQz9nhmZsZ7HPFNN92k7e1trzWFNMTXajVNTEwEZTFC+ld8shxJlI3V63XV63WvdbUeOzo6GvuYpDIve/XJRJmcnPTOwI2OjmpjY0N33XVX/IM+9jGvcwAAALRD8OLJZ+xxa+NA3z0+fN9FX19f9868VCoVLSwseGUxzp8/730jHhKIpLVsrFareQ8faPE55ty5c1pdXfV6/pBsTUhJX2sst2/GSmJCGQAA2D+Clx4ql8uJ7NkSknlJqvk8T2VjZ8+e1crKitf5z507533NQibBhWRrxsfHvSeUHTx4UNeuXfM6DxPKAABAtxC8ePIZe1ypVLwnQI2MjHg331erVY2NjXm9u3/p0iUtLCx4nWc/Y48pG4svpH8lqQllIyMj3gF5a0LZQY8JZd/EhDIAABCB4MVTtVqNPZ1pbGzMe/+VkOb7crmspaUlr4zAhQsXvM4hJTf2uOhlYyGBZVITygYGBrwDnkeHhvR3trd17Nix2McwoQwAAEQpZPCynw0nfZRKJe8bvZASsJB33UMa9ltlUz6SLBvzOc+tt97q3beRVNnYwsJCYhPKKpWKV3B96623am1tzWu8crVa1fb2duxBF5KYUAYAACIVMnjZz4aTm8PDsR9bqVS8byhDmu8PHz6sWq0W9O6+z434yspKasvGQsrAfIOqw4cPx35sS5r3hglp2B8aGtLQ0JBXwDM6OqrNzc3Ygy4AAADaKWTwcv2Gkz6t38uSqt///bEfXy6XE8m8rK2taW1trdBlY75lYCHBjhS2N4zkVza21zFR1yBkQtnGxobW1ta8ysYWFxc17BHAS40+mZBNUQEAAG5UyOCl5YlKRW+emor12KeeekqS9IEXvSj284dMDgvJvIQ0eKe5bCyEbxlYSM9LSMATUja21zFRpXMhE8pC9oYZHh72DniGhoZkZl0br1yR9GCbxzNeGQCAfCt08FKpVHT77bfHeuyFCxc04DEtSWrctPmWjSW94aSv0IZ1H50Ci36OSg4JeJIqG0tqQllIqVm5XFa5XPbqk2G8MgAAiFLo4KVarerOO++M9dhHHnnEO3ip1WreN5V7jT2OKhkK6V8J2XCy1fPiIzTDIaVvVLKU3rKxpCaUhewNMzo6qmvXrun8+fOxj/md5p+TEWWXO85pICI7+dKtLcYrAwCQc4UOXkZHRvSyl70s1mOnp6c1ODjoddM2Njbmvaa9xh5HZR5C+ldCNpxMc8+Lr5CgqlUy5xu8dHtCWdTXwMLCgnf2rVX65XPM9PS0d/9KrVZTqVTyyrw8Pjys79/a0jd/8zfv+lzrtd7ol86cYbwyAAA5V+jgxUe5XPbenG9oaEjOOa9jkhp73JKXnhdfSY1KTmpjy5DA8vz58yqXy17HLC0teX9NHzhwQGbmFfiPjIxoY2MjMjNaLpej/+3OnPFaFwAAyB6Cl5jK5bJ3vf/Gxob3efYae5yGsrGi9ryEjEruxcaW/SwbC+l52dra0ubmpneT/+DgYOR45atXrzJ2GQCAgip08FKtVmOXsrRKX3yEZF72Gnvcz7Kx0dFR7xv4ove8hAZanTI83SwbCxmvHJJ58S01q1ar2trairzWy8vLqcvkAQCAZBQ6ePExNDQUFLz4vuseckOZVNlYUhtOFr3nJaTULKRsLGS8cujGlr5KpZL3gIyWqPHKnTBeGQCA7CB4ialUKnkHL6VSyfvGdXV11btsLCmXLl3ynjhGz4t/ILJXqVm3ysaq1WrHr7Uon/rUp7S0tOR1npAJZQcPHtS1a9e8zsN4ZQAA8o/gJaZKpeLdSL+1teUd8ISUjSU19lhKZp+XtEqq5yWk1Ozw4cM6fPhwUJbPt2E/ie+DkZERVSqVyLLOer0e+fF3Nv886JGx+aadHcYrAwCQIZkPXu4x089LOuBxTMitd7lc9r4Bm5ycVLVa7XnZWFJjj5Pa56WTfjbsh2ZrQvj2vOw1XjlKSNlYSJN/qVTS1taW1zG1Ws174MUjpZK+d2dH3/iN3xj7mPsffZTxygAAZEjmg5efl/S6wGMrN98cu5RldHTUuwZ/fX1dGxsbiZWN9XrscUjDfkipWRKS3KQyreOVQ8rGQoLrkMyLb3ZHalyD7e1t3XLLLfEPevRR7/MAAID+yXzw0sq41CU9FjO42N7Z0bKku06ejH2e0dFRlcvlRDb0i9qArxd8y5lCGva7PV45aw37vSg161YAG1I2FpKtCcm8HDhwQNvb25FvLly+fDny46Ojo9rY2NBdd90V/0Qf+5jXugAAQH9lPnhpmRkZ0Y8fPRrrsY823239QKkkxdx7YnV11XtDv7W1Ne+ysb32bOlnz0tIFiUkW5OEpBr2peTGK/sKCZJDsjUhQfz29rb3prAjIyMaGhpK3XAIAADQPbkJXoaGhmKXizz++OMaGBjwnn7kezNVrVY1NjbmdczCwkLqe15o2I+vF+OVu9Xzsp+9YXyb/H33hhkcHNTa2lrkxpZPP/102w0vzWNEMgAAyJ7cBC/DtVrscpE/+ZM/8e5fGRwc1Pb2ttcx5XJZS0tLhe55CdGppCuLDftp7XlJ+94wQ0NDXmVjY2Nj3k3+LT57w7AvDAAA/ZOb4KXiUS5SLpe9syghZSybm5vejcdp7nkJEdJbErKxpa+kGvbT3POS5r1hyuWydxZleHjY6/HSPvaGYV8YAAD6IjfBi4+Q3bsPHjyooaEh7xIbKeyG0qdkKKlpXkkFIp02tsxaw36S45V9pXlvGN9xzFKjh2dkZCRyD5h2WnvDjMf8efCN7AuTOWY2Luk/SzomyUn6R865T/V3VQCAUIUMXkIyLxsbG95jjy9cuKDNzU3vY3wltWdLSMN+p0Ckn9LcsB/Kd23nzp3zzvKFZAVD9oYJybyUy2WVy2WvXrZHh4b0d7a3dezYsViPZ1+YTHpA0secc280syFJtX4vCAAQrpDBS6lU8n4neGhoyLvpeD9lYz4lQyH9KzTs+0tzw35In0zI10Bok7/kl4E8dOiQ99d0qVTyDnhqtZo2Nzf1kpe8JN4B7AuTKWY2JunbJf2IJDnnNiSFNUYBAFIhN8FLtVqNXS5Sq9W8My9DQ0PeN+4hN22hfNeW5ob9tEqqBKwXfTJRQnpefCfhXX+cb6Dsu7Hl1taWNjc3204ii1IulzUwMKCbb77Z61zIjBdKelbSe8zsZZIelvRW59zVXpxsZmZmz6EiDIcAgP3JTfBSq9Vil4uEbDiZtrKnbkhrw34n/Zw2FqrbpWbtsm++kup5SXJjS9/9ZMrlskqlUi6/vyGp8Tvur0j6cefcQ7oS6IAAACAASURBVGb2gKS3S/r56x9kZiclnZSkqakpTU9PB51sbm5OKysrOnPmzK7P3dH803c4xKV779Uj990XtJ5+WVxcDL6GeVD01y9xDYr++qXeXoPcBC8+zMy7F+XSpUva2trqec/LXhtO9rNhP0RIn0wS0jxtLLTULAlp3tiy03jldkZHR7W9ve3V5I9MeUbSM865h5p//201gpfncc6dknRKko4fP+5OnDgRfMLp6WlFHf/FsTEtvvvdOhgz61I/fVrjkg5Xq5HPl2btrkFRFP31S1yDor9+qbfXIDfBy+rqauxyka2traCSFN9jQnpekt5w0kfohK5uZniYNtZeSMN+nja2DGnyL5VK2tnZ0fnz52M9fqr5p0/pj0T5T7845y6a2dNm9mLn3DlJr5H0hX6sZfXYMX35gQdif++cNWM4BABEyE3w4mNkZETDw8NBZWNJ9bz0umE/pH8lJIuS1j6ZtE8bS6JhP28bW4aMVx4ZGfEamx68L4zE3jD98+OS3t+cNPYlSf+wz+sBAOxDboIXn56X1mZ2vS4bW11dVa2WzFRO3+xGK/Pgk31g2lgyJWBJNey3FHljy5GREVUqldhlY619YQ54nONuib1h+sg5d1bS8X6vAwDQHbkJXnyUSiXvkpSQsrHW2OM0ylMWJSlpydZ0q2E/RGiTv6/JyUlNTk5Gnqdd6WDIeOVaraadnZ3Yb3w8PDCg1+/s6K677op9jv/0xS9S/gMAQJdkJni530x3S3rwho+H5ABKpZIqlUrPy8Za71KntWE/rVmUTv0oWZs2lreG/VC+pWbnz59XuVyOPKbd10DIeOWdnR2vsrFKpaLt7W298IUvjH+SL37Ra00AAKC9zAQvd6tzjfn2yEjs5xoYGND6+rrXzdTa2pqq1Wrsx0uN5mbffTHS3LAfIqQxvl6vS5LGx8d7saTnzpHEtLG8NeyHHBPS81Iul717WELGKw8MDGhjYyN2dqi1L0xaM6oAAORdZoKXZw4d0icXFzVwQ037xuamliXddfJk7OfyfbdVagQVGxsbQT0veWnYTyoQGR8fbxskZG3aWFJCMjwhDfshx4T0vHQadtHuayBkvLLvNMBW5oVNLQEA6I/MBC+//5rX6LHHHtv1jmdrM7BHXvrS2M8VskllyI1Rkj0vSZSAhfTJdApE+ikt/Svt+E4bCz2P1N1peFFCRiWHlI2FjFc+cOCAtre3Y/e8HDhwQBsbG6n7egYAoCgyE7wcOXJE8/Pzu24a/uzP/sx7wtDm5qb35pGtsrE09oiESMOeLVmT5mljac4KhQTwIWVjIeOVt7e3NTg4GPvx1WpVpVKJTS0BAOiTzAQvN910kw4ePLirXGNwcNC7BKxcLnu/Q1utVjU2NuZ1TMhY2aR6UZg25i+p/pW9ztPPaWMh9jOhzKdsLGS88uDgoLa3t2M/fmhoyPvnTQsbWwIAsH+ZCV6mpqZ06NChXTd1g4OD3pmXcrnsvUllyLjXkIb9vZrvizBtLG+SKjULkVSTf0j2pd33QaeysbGxMa+v6cnJSW1sbMR+vO80M4mNLQEA6KbMBC/thGReDh486F1icuXKFS0tLfW8YT+p5vs0Y1RyWAmYbyDSusY+gW9Iw37INLyQCWVra2taW1vznrrno1qtepezsbElAADdk/ngZWBgwKtmXZI2Nja8J4f59odI4Q37nTIiUSUzrZt9nzUyKjlfo5JDgqRarRacfQvZ78hHpwll3SwbK5fLXpnbkMzLZyS9XmJjSwAAuiAzwcttt92mqampXY2ytVpNg4ODXjdTo6OjQeNefUtSkpJU/wqjkpNrik9i2lhSfTL76XnxMTk5qcnJSa/zHDp0yCv7UqvVvMrMpEafzM7ODhtbAgDQBZkJXtpplYz53IhtbGyoUqn0/GYqSYxK9pPmUcl5mzYWyrfnpdN45XZWVla8sinlclnlcjn2aGWp8cbH1tYWG1sCANAFsX5rm9l5NfpOtyVtOeeOm9mkpA9KOiLpvKQfdM751YZ0QUjPy9DQUGITvdI6bYxRyf6SCiqSnDaW1ob9pMYrl0olbW1teT3ed0DI8PCwNjc32dgSAIAu8Mm8/E3n3Neu+/vbJf2Rc+5dZvb25t/f1tXVXefOO+/UHXfcsesdz1qtplKp5HVDPTIy4n3TduHCBe+9YS5cuBD7sS1JTRtjVLK/pPpXkhLyb5lUw/7CwoImJia8el46jVdux3fz2a2tLW1ubnplYltvsOQpKwYAQL/sp2zsDfqLyZ/vlfSgehi8tGNm2tra8r4x8t3nZXNzU+Vy2euYVsN+WqeNFTmLkqSkRiUnOZK51w37viPGQ/nu3XTgwAFtb297lY2Njo5qe3ubjS0BAOiCuMGLk/RxM3OS/pNz7pSkKefcVyXJOfdVM+vpW/Jf+cpX9PTTT+96x3Nra8t72lirbKxbm+Z1m++0MaRXUqVmae6TCW3YT6LnxXdvGN99YaSwCWUAACBa3N+qr3TOzTUDlE+Y2ZNxT2BmJyWdlBobTU5PTwcss9GjMTExocuXLz/v4zs7OzIzzc/Px36u5eVlLS8ve62ldUPkMzVrdnZWkrzWttcxi4uLu9Z9+fLl56Z6xXXt2jUNDw97rS1EyNo6WVlZee4a7df4+LgOHTrUlefqppCvgV6cp9/HLC8v69q1a7s+t7KyojNnzkQes7W1pWq1Gvs8jz32mCTFzsQ++eSTcs5Fnr+dtbU1bW9vewVVJ5p/1k+f1tmIHpuKpK987/fqyR/7sdjPCQBAHsQKXpxzc80/L5nZRyTdI2nezG5pZl1ukRTZ/d3M0pySpOPHj7sTJ05EPWxPMzMz+vjHP65XvOIVz/v44cOHNTg4KJ/nPX36tK5evaqpqanYx9x1112q1Wpe52ndrPiWjXU6Znp6etcaWqVmvn0vSUwBC11bO2fOnNn1NRAidC+ZNFyz2dnZyD1DQkcld/Prs1vHtL43o46J+h6QpA9+8IOSpB/6oR+KfR4z09LSUuyvqc997nMyM6+vwfe85z3a2tryev2trSnHpfb7vSws6M7An6cAAGTVnsGLmY1IGnDOLTf//zsl/YKk35P0Zknvav75u71caLuysfX1dZVKpZ73vIRsOJnktLEk+lfSMKGsCPu8pLkErHWuNE4ok9qPV27Ht2zMd18YqTEq2XcK2jubf060KYm1gQG9kvJRAEABxcm8TEn6SHM8aEnSbznnPmZmn5H0ITN7i6SnJP1A75bZedpY6KjkkKbjXk8bS/NErzxNKEvzPi9pHpXcCvh8vgZCJ5T5Cgl21tbWtLa2Fvsa+O4LI4X1vHx2cFDf65y+8dixyM+Xy2V95v77vZ8XAICs2/O3qnPuS5JeFvHxy5Je04tF+QgZlRzy7vTCwoL3DVXItLFWRqBdVqBbo5LTkEXJmrRnRHyFvJ5arRb8NdDrhv1O45Xb+dSnPqWlpaXYj/fdF0ZqXLOdnZ2gjS1f8pKXRH7eN5MDAEBeZH4MztbWlveo5LW1Na+m3pZarZbKXexDsht5yqIkJc3ZGimZUclJ7VsTWjLma3JyUpOTk7G/r1uDEnx+Dhw8eNA72KhUKhocHGy7saVv6RoAAHmRmeDl6NGjmp2d3XXTEFJPXq1Wvfd3SFJSo5KLnEVJCqOSw3Qar9xpk0rfnhff8cpzc3PePXarq6td39jy2Wefjf1cAADkSWaCl5mZmchyka2tLZXL5aBSEd8bkFqtFvvxSdpPCRh6KzRbk9bzpFlItqZcLnu9+RFSNhayseXw8LD3eQAAKILMBC/tlEol70lGS0tLWltb8zomZNpYUvJWAtYpGGvX7yAlM8YY6Z02FroZphS/DGx6etr7Z8f29rb3Rrq+1wsAgKLIfPAyMDCg9fX1oJ6XXk8bO3/+vPeUpTSPSg4RkhUK2YOlXq+rXq97nSvNwU6nAKFdAJfE6wkJeEOnjfk27IeOV/YRUm46ODio7e1tr2NGRkZUqVR05MiRyM/7jnoHACAvMh+87OzseI9KDul5CZk2tr6+7v0OalYzIu2EBCLj4+Ntb8Tb9Tv4ri0k2JHSGyCE9q8kNUxA8mtyD9kjKaRhv3Vc3Gvguy+M1BgKsLGx4bWmkPHKAAAUQWZ+Q7Zr2G+90xpSKuIrdNpYGjMiIULK0zoFIt3k2/OR5k0qk9rnJc1N/kk17F+5csVrXb77wkhhk8FCBpEAAFAEmQle2jXsb25uqlwuez3XlStXUtuwv9dNdbf2eQmV1vI0X2kfe5yEvDX5h5SMVatVrzclfPeFkRpDAZqb/MZWKpVUKpXaNvmHZpkAAMi6zAQv7ZTLZQ0PD3sdk+aG/aSa75lQ5i/NmYrWufISWIXYT8N+XL77wkjSoUOHvLMvg4OD3gEPAABFkJngpV3Z2PDwsK5du9bzhv0kJbHPS94mlCUhzdmavAVWSTTft/iUmp0/f967X25sbMy7T2Z0dFRXr15tG1zNz897rQEAgLzITPDSrmxsaWnJe/JO2jepTEpeSsDSLKmgIskSMN9ApJXh8wmUk5o25vv4kJ6XkI0tpcbEMcrGAAB4vswEL+20ysZ63bAfMv0oqTIrSsDSK299JSHBWK1W8w6UQ6eN+TbsLywsaGJioqc9LyEbW5ZKJcrGAACIkJngpV3Z2M0335zI+UNGJYfcTIUEFZSA5U9a93lJKhgLmTaWxHjlkJ6XkI0tfYMdAACKIjPBS7uysY2NDQ0NDSWyBt9Rya1+Bx+hQUWeSsA6ZZLa3biHSGuzel73eem1JBr2Q4SUqJJ5AQAgWmaCl3Y2Nja0sbHhdQN24cIFbW5uJjIqOU9BRYikNrb0leZNKtnnJVm+DfvlctnreodsbHnw4EFdu3Yt9uMBACiKzAQv7crGRkdHtb6+7vVcIXvDJDUqOc37vCQViHTa2LJb09bSvEllUvLWjxPC93u6XC57bx4Z0uR/7do1DQ4Oep0HAIAiyEzw0kmlUgkqFUljRiTN/Ssha+sUiPRTmsce51HIqORKpRL5uW5OG/MtNQv52RHS5D84OKjt7W2vYwAAKILMBC9p6HlJShL7vOxnQlkag74kpL3MKq2BVch1q1Qq3hnVkAEZSewnE9LkPzk5qY2NjR6uCgCAbMpM8NLO0NAQI38DpDnDk1ZpLrPKa5N/2qaNtY7xef0hG1uOjo56D/sAAKAIMhO8tOt5SaocJ837vIQqchYlSUnc7Cc5wthXUtmn0GljPsHIlStXvNcV0vOysrKiUikzP54BAEhM5n87Li4u6urVqz2fHBayz0vRMxVI981+iDRnn0L4loxVq1XvkelJbWwJAEARZCZ4adfzsrCwIOec13OFvBPaCnh8blpafSU+vSUh2Zr99K+gt9J+s5/WPpmkhDbs+wjd2FJqPxhgfn7eex0AAORBZoKXdlo9LyGlIj5CRiUn1VdC/0r+dAoq2k3b8g0q8pYVSkoSPS+SOr7BMjc35/18AADkQWaCl273vPhmUUIl1VdC/0p+JBVUZDUr1GlUcrvxyu34ThsLmUwWkuldWlryzigDAFAEmQlekE0hJW2dtLtxDZHW8qe9goqoaVtp7iuR/MvTWl8zPtnE0PHKPkKGAoT0vAwNDalarVI2BgDADQheUmavm/2om/c096+ElLQloV6vq16vewdWaQ140iwkk1Sr1dpmE/fa66iX08ZC9oUJ7XmhbAwAgN0yE7y0a9i/dOmStra2ej5tLClp7l9Jw8aW/d6kU0pvv0dam+/TXJ7mG4iE7AsTIqRHBgCAIshM8NLO1taW934IIc33Sep0s9+tm/cQaQ6sfIXeUKc1QMjrJpW9FjJtLKRhv1wue/e8jI2NUTbWBWZ2XtKypG1JW8654/1dEQBgPzITvOzVsO9zQ9/acNL3BsS3nj6N5VItaciiZE2ap3OxSWUyQt70KJfL2tzc9DpmryZ/ysa8/U3n3Nf6vQgAwP5lJnjpppANJ9fX170nGaU16yDlK4uSlLxla0KkuQQsCSEN+629YXrd5A8AQBEUMniRwkclp3GTylBFzqIkpeiZiqT5BoohDfhJ2KvJn7IxL07Sx83MSfpPzrlTvTpRu95MKWyUtyTVT5/WWbPYj1+W9LqHHpLuucf7XACQBZkJXtr9UlhbW1O1Wu3Tqjoju4G0ZyrSnBXy3eclZLzy6upqz6eNtc6TxMaWiPRK59ycmR2W9Akze9I596fXP8DMTko6KUlTU1Oanp4OOtHc3JxWVlZ05syZXZ+7ePGiNjc3Y5f83dH8c1zSqzzXcenee/XIffd5HtU9i4uLwdcwD4r++iWuQdFfv9Tba5CZ4KWdarWa6l/ySWQ39tO/guJKc1ao2+OV22n1v8UVMm2sFxtb0vMSn3NurvnnJTP7iKR7JP3pDY85JemUJB0/ftydOHEi+HzT09OKOr71bxn36/OLY2NafPe7ddAj61I/fVrjkg5Xq5FrSEq7a1AURX/9Eteg6K9f6u01yEzw0q5hv1UbnpdRySHI8ORPp4xIu8yDb0YkyaxQNzM87SbuhQ4g6PW0sYWFBU1MTNDz0gdmNiJpwDm33Pz/75T0C31eViyrx47pyw884PV1c9bMO0sDAFmTmeClnaWlJa2trXkdk/ZRySHoX8mPNGdEQuTp9YSWjPmi56VrpiR9xBrZi5Kk33LOfay/SwIA7EfmgxepUTrGjTvyYq+MSFTmIc0TvdLe9+NjP9PGkDzn3Jckvazf6wAAdE8ugpc82at/JapkKM39KyH9OJ20K5kKkdZxxGmX5ib/tOr2xpb0vAAAiiozwUs3p435NulKyW1Smbf+lZDXk4R6va56ve4dWKX5JjyJoCJPJWBJSWpjSwAAiiAzwUs7IdPG0r5JZaf+lXbNyknYz1Szbq25W68/9LVI6bwJTyqoSHsJWMg+Lz7f1yGjknuxsSU9LwCAospM8NJu2lhoPXkSm1TmTZ6yQqE34WktmUp7UBEiiX1eKpWKV0bVd18YKb0bYQIAkEWZCV7QXWnIomQNJVPJSWqfl5a4x4SUnIZMGwMAANEyH7xcuXIltfu8pHnzyDxlUZKSt2xNmnW61t3c58VXktPGOmV4aNgHABRV5oOXNO/zkvYAochZlKSkPVuTVGBV9ADOt9TsypUrPVwNAADZlZngZa9pY2m9CU8iQEhzhqfo0tyLklRglfYArtdC3iipVqsd+/Jo2AcAFFVmgpd2QqaN5U3aMzxIp6QCqzQHcElgY0sAALonM8FLt6eN5Q0lYPnSqcyq3bStPJVZAQAARMlM8AIURdHLrJLm04+T5NhjGvYBANitkMFLyLhTekSQlL3KrKKmbaV5I0gpuayQ79p894ZJap8X9oUBACBaIYOXhYUF770XQm5aQgKevZrvo0qGCKyQlDRnhZLYGyZ0nxffnx0LCwuamJigYR8AgBsUMniR1HGST5TWfh0+Qpri89Z8HzIJrZN2/R4h8tYjkkRGJMnm+3avp9PXgO/rCblevs33bGwJAED3ZCZ4aTcq+cKFC9rc3Exkk8qkmuI7nafdBn1pFRKMJaFer6ter3sHVmkNeNKcEQmRp9fDtDEAALonM8FLO5ubmyqXy17HJLVJZd7sZz+ZbgVc3QreQl+LlL6bYyl/44g7vZ52XwNpfj0AAKA7MhO87DUqOUvZiKzKU0lb6M1+mpvVkS9MGwMAYLfMBC9ZtJ9MRa+lIYuSNXkqZUoaQZ8fMsMAAEQjeOmhNGcq0ry2tMpbaZaUTFCRt6DP95qFjEreq0+GaWMAgKKKHbyY2aCkP5d0wTn3WjP7ekkfkDQp6bOS3uSc2+jNMrMrzZmKNK8tT9KadUgqqEh70Ofz7+O7L4yU3N4wAAAUgU/m5a2SnpA01vz7L0v6FefcB8zs1yW9RdKvdXl9QKalOeuQ9qAiRK83qfTdF6Z1DkYlAwDQHbGCFzO7XdL3SPolST9pZibp1ZL+fvMh75X0DmUkeAm5mUjjuN+WNPfWFF0eA4S0SmKTypB/G0YlAwDQPXEzL/dL+hlJB5p/PySp7pzbav79GUm3dXltz9Nun5eQPVsWFha839lMc78H/Sv50ymD0G6TxiI3uEv5CxSZNgYAwG57Bi9m9lpJl5xzD5vZq1ofjnioa3P8SUknJWlqakrT09NBC52bm9PKyorOnDnzvI9/4QtfkCRdvnw59nPNzMyoUql43+jMz8/rsccei/342dnZ547rlsXFxV3XsHWeF73oRV7P5ft60iDq9efN5cuXVa/X234+6vvg2rVrGh4e9v4a8DU7O6tr167tOv9exsfHdejQIa/zSNHfO938GvD9Hg35ng49Znl5WdeuXYv8/NbWVu6/DwAAiBIn8/JKSa83s++WVFWj5+V+SeNmVmpmX26XFPlWoHPulKRTknT8+HF34sSJ4MVOT0/rxuM/+MEPanV1VS9/+ctjP8+hQ4c0MTGh7/zO74x9TEhp1sjISNeb4qOuQevd2SI030e9/qLp59dAUiO2W3vqRGUTZ2dnddddd0Ue55t98r1uIdc55JipqamOx/B9AAAoqj2DF+fcz0r6WUlqZl5+2jn3D8zsw5LeqMbEsTdL+t0errOjWq3W85u2pEqz9ro5jCoZon8FSUmqNCvNgw4AAED/7Gefl7dJ+oCZ/aKkRyT9RneWlF5JjBamfwWh0jqSOUSnIGl+fj7y+zDN/SsAAKA7vIIX59yDkh5s/v+XJN3T/SWhU5DU7sYNxUamIjkhm1RWKhWvc7DPCwAA0faTeQH2FNIj0Um7SVsh0pp1CJG3SVtJ6vUmlZVKRevr615rYp8XAACiEbygp0LK4JJQr9dVr9e9A6s8BTxJSmtJm2/GKmSTyhbffV4YlQwAwG6ZCV66uc9LUvK2eWRSk6Y66VbZXOhrkSiz8pXmkjbfjFVS2SpKxgAAiJaZ4KWdWq2W2l/0eWu+z9PrCS2zSmsGIVQSryfJkra8/PscPnxYhw8f7tj7BgBAEWUmeDl69KhmZ2cz16yexISyEGnIomRNmjMIIXg96X49AABgt8wEL9106dIlLSwseB2T5nKuEHnKoiQlb9mavDX55+31AACA3TIfvIQEIhcuXPA+Tx5v3IucRUkK2QAAAIDuyXzwsrCw4D1WtNUnk5cb97wNBsgTsgEAAADdk/ngRWoEI3kJREJQAoYkJVUG1+48nfb6SWPzvRS2sWVaB5EAANBPuQhekpD27AYlYPnS6Wa33c17GvdFkcLK4NJcbucbiIRsbLm6uso+LwAARCB4iYnsBpKS5hv3pMrgOp2n3V4/SZTbhfzbhGxsGdLLBwBAERC8eCC7gSTsFSBE3bzTJxPON5Pim+EK+bdhnxcAAKJlJniZmZmJLJdZXV1VrVbr06q6b6/ytKhrQPM9EMY3k8IkOAAA+iszwUs7rclheUF5GkKldT+ZNPMtgyPDBQBAf2UmeDl69KhmZ2cLUbbVqTytXb0/ii3NfTKtcxFYAQCA/cpM8AKgvTTvJ5NkYEWQBABAvhG8oKdCRkx30mmPD1/ctCYjqcAq7dknAACwfwQvPZT2vWGSENLDk4R6va56ve7970PAk15pzj4BAIDuyHzwErIfQlI303lrvt9PMNatPp1u9fyEvhYpX+/SU2aVDN/rfO7cOVUqlbafZ5NKAEBRZT54WVhY0Orqqtcxe+1eHSU04MnT3jB5CsZC36XP080+ZVZhfL8GWkGyz/dNpVLR+vq699oQzcwGJf25pAvOudf2ez0AgHCZD16kxrhknwChdRPqI6034aHSkEXJmrzd7FNm5S/ka6BWqwV/37BJZde8VdITksb6vRAAwP7kIngJUeSbcClfWZSkcLOfrHYZjk5DG3qd5eJrIHvM7HZJ3yPplyT9ZJ+XAwDYp8IGLyCAS0qeSs2SkrcsF/rqfkk/I+lAvxfi87Ngr76nTuqnT+usWezHL0t63UMPSffcE3Q+AEhSZoKXmZmZyHdcV1dXVavV+rQqoLM83oQnEYx1ynC0G9oQmuHwfT1FDyyzxMxeK+mSc+5hM3tVh8edlHRSkqampjQ9PR18zsXFxcjjL1++rHq9Hvt5Ll68qM3NTa/hDHc0/xyX9KrYRzVcuvdePXLffZ5HRWt3DYqi6K9f4hoU/fVLvb0GmQle2qnVapqYmOj5edI89jjNayu6vJUZ5S0Y8309aX4tiPRKSa83s++WVJU0Zma/6Zz74esf5Jw7JemUJB0/ftydOHEi+ITT09Paz/EtrZ8DPtnxL46NafHd79ZBj6xL/fRpjUs6XK12Zd1S965BVhX99Utcg6K/fqm31yAzwcvRo0c1OzvbtzKnNPeIpHltCNMpG9Cu5yOJjEDegjHf15Pm14LdnHM/K+lnJamZefnpGwOXPFk9dkxffuABvwE2Zt5ZGgDop8wEL2mQ5h6RNK8NfvKW3QAAAOgWghcgZfbKBkT1fJARyJ9O2Tc2qfTnnHtQ0oN9XgYAYJ8IXlJmr/6VqJIh+leAhrROdktiY0sAAIqA4CVl6F9BqLTeuEvJrC2t5Xa92NiSTSoBAEVF8JJCe9200NuCG6X1xl1Kbm1pHSaQ1nUBAJBFmQ9eLl26pIWFBa9jKLNC3qT5BjnNawuR5gwXAAB5l/ngZWFhQaurq17HUGYFIESaM1wAABRB5oMXqVEf7lNK1WqK99nYMSRbw+aRQL7kLYsEAEDW5CJ48ZVUUzzN92E6BX3tNmgMQSkPJVAAACBbMhO8zMzMRN64rq6uqlareT9fUps6Fn3zyJDsU71elySNj4/3YknPnaNer3uvLU837pRAAQCArMlM8NJOrVbTxMREv5dRCEkFIuPj422DhG5NWwst6ZPyc+Oe9hKodlmhTtm3kODSN/uUpwAWAICsyUzwcvToUc3OzhY6i9FvIWVwnQKRfgq9cafMKhlJZYV8z5O3ABYAgKzJTPCC7trPMIGiBpCUWSWnU3DZLvsWkhXyDWLTkHmSywXr0AAAIABJREFUGt+/AAAUEcFLQTFMwF/ay6yQH0X+PgMAoBOCl5jyOPa4yFmUJOWt1CxvrycJ3b5m8/Pz3VgWAACZQ/ASE5kKhMhbqVneXk8SuGYAAHQPwYuHtGYq8pgVyou8lZrl7fUkgWsGAED3ELzkAFmh/OlUZtRuVHDRS7MkStoAAMg7gpecSGtWCP4oMwrDdQMAIP8IXoCU2avMKGpUMGVGyZVnkd0BAKB/Mh+8XLp0SQsLC17HpLnfY6/+laiSoTS/HkDKzw0/2R0AAPor88HLwsKCVldXvY5Jc78H/SvImzzd8NN8DwBAf2UmeJmZmYnMOqyurqpWq+Wq36NT/0q73cWBtGY3uOEHAADdkpngpZ1araaJiYl+LwPoqzxlNwAAANrJTPBy9OhRzc7OdiXrwL4oyBuyG/nTKZM2NzeX8GoAAEiHzAQv3ZRUXwlBEpB+vuV2SZTa0cMGAEC0QgYvUjL7otB8D6Sb7/daUqV2ccZlAwBQRIUNXpLC5pFIs7Q2+Uvt1xY1uKPFd22+5XaU2gEA0F97Bi9mVpX0p5Iqzcf/tnPuX5rZ10v6gKRJSZ+V9Cbn3EYvF4ti6FRu180b16JLc5N/mtcWIs1BIgAAWRIn87Iu6dXOuRUzK0v6pJn9oaSflPQrzrkPmNmvS3qLpF/r4VqRQSF9P/V6XZI0Pj7udUy9Xvc6V9FvDtPc5N9pbe3Ghac1K5K3QAwAgH7aM3hxzjlJK82/lpv/OUmvlvT3mx9/r6R3iOAl15IKRMbHx9sGFu1uXH3XFhLsSAQ88JfmIBEAgKyJ1fNiZoOSHpZ0VNKvSpqVVHfObTUf8oyk23qyQqRGyACCToFIN/neIIZOgmudCwAAAMmLFbw457Yl3W1m45I+IukvRz0s6lgzOynppCRNTU1peno6cKnS4uLiruNnZ2cl+U3fCTnm8uXLz2UR4rp27ZqGh4e7Ohko6hokpXXdXvSiF3kdNz8/r8cee6wra+j367927ZrOnDnjddz4+LgOHTrUtXUU/Rq0e/0h39e+kjhHnPP082sAAIB+8po25pyrm9mDkr5F0riZlZrZl9slRe6a5pw7JemUJB0/ftydOHEieLHT09O68fhWeYXPRK/QY0L2YOl21iHqGoQIyTyMjIz0fXpat15/iP3s29PNa1b0a9Du9Yd8X/tK4hxxztPPrwEAAPopzrSxmyVtNgOXYUnfIemXJf2JpDeqMXHszZJ+t5cLTYN+37h3E3vQ+AvtXcjTpKm0928kca3z9O8JAEDWxMm83CLpvc2+lwFJH3LOfdTMviDpA2b2i5IekfQbPVwneiBPwVhaMWkqOUlca/49AQDorzjTxj4n6eURH/+SpHt6sSgfly5d0sLCgtcxIeVfabafUh70VtozFXmSxLXm3xMAgP7y6nnpp5mZmcgNCi9cuOD9XHkrf6IEDAAAAEWQmeClnVqtpomJicKXP1ECli+d+iqigngp3X0V9In463TN5uYi56MAAJB7mQlejh49qtnZWW7QkXt566vI2+tJAllRAACiZSZ4KYq9+lei3nWnfyVf9uqrmJ+f3xXEp7mvgj4Rf3G+BgAAKKKBfi8Az9fqX/FB/woAAACKoJCZl7RP5+rUvxL1rjsg0VcSyve6cc0AAOifQgYvTOdC3tBXEsb3uoVeMwJLAAC6o5DBi8R0LuQLfSVhfK9byDUjsAQAoHsKG7wkIe3laUCatctWtBsVLaUzW0FgCQBA99Cw30M03wNhDh8+7B3Er6yseL9ZAAAAsoXMS49RnoY0S2svRqdsRbuhFWQrAADIP4IXoKDoxQAAAFlD8AIUFL0YAAAga+h5AQDkkplVzezTZvaomX3ezO7r95oAAPtD5gU9FTJxrZNOk6Z8pXEyFdIvrX1CiLQu6dXOuRUzK0v6pJn9oXPuz/q9MABAGIKXmBh7HCZkQ9Ak1Ot11et1739TbkKLjT6hbHHOOUmtkY/l5n+ufysCAOxXZoKXmZmZyHfdk7oxDrkJz9vY4/0EcN2auNZu0pSv0NcicROa5sxDr9eWZJ9Qp9cyNzcX9JxFZGaDkh6WdFTSrzrnHurzkmLz/Xo+d+6cKpVK0Lnqp0/rrFnsxy9Let1DD0n33BN0PgAIlZngpZ3V1VWtrq56/YAPDXiKPvY4TwEczeph0px5SPPafKXxeyarnHPbku42s3FJHzGzY865x69/jJmdlHRSkqampjQ9PR18vsXFxX0d33L58mXV63WvYy5evKjNzU2v4PaO5p/jkl7ldTbp0r336pH7drcRdesaZFXRX7/ENSj665d6ew0yE7wcPXpUs7Ozu4KHs2fPshFkgDRkUbIozVmHJKQ56Evz2nzt9Vrm5+cTXE0+OOfqZvagpO+S9PgNnzsl6ZQkHT9+3J04cSL4PNPT09rP8fvR+nr2+Rn9xbExLb773TrokXWpnz6tcUmHq9XI19rPa5AGRX/9Eteg6K9f6u01yEzw0knRb6hD5CmLkpQ8vbMPFIGZ3Sxpsxm4DEv6Dkm/3OdlpcrqsWP68gMPeP0OPWvmnaUBgG7JRfCCMAR9fvL0zj5QELdIem+z72VA0oeccx/t85oAAPtA8JIDTEJLt6KXmgH94pz7nKSX93sdAIDuIXjJAUrA0otSswYCOAAA0A0ELzlBCVg6UWpGACcRvAEA0C0EL0AKdbrZjdrvSErvzW7aAzjfwML3OhO8AQDQPYUMXugRQZpxs5sc32sdcp3THrwBAJAlhQxe0twjsldgFfWuO4FVvsTZ4yNqvyP48w0suM4AAPRXIYMXKb09ImkOrAAAAIB+KmzwkmadAquod90BAACAIiB4AXIibxOt2r2edgMLpHS/HgAAsH8EL0AO5K3JP2+vBwAAdAfBSw8x1QxJydtEq06vp13pZJpfDwAA6A6Clx6i+R5AqE5lgHNzcwmvBgCAdCB46bG0TjUDECaJ3iLewAAAIFpmgpeZmRn2OAG6LG9N/r2WVC9OnL1+AAAooswEL+1QZgWEoSneX956iwAAyJrMBC9Hjx7V7OwsJVhAl3AjDgAAsiYzwUs7TPQC0EIZHAAA+Zb54CWpiV4ESUC6UQYHAED+ZT54kZKZ6MXY4zAhQV8n7K6OdiiDAwAg/3IRvCSFscf+QoK+EPV6XfV63StQItgJQ2kWAADoF4IXxLaf0rluBX3tdlf3XVtIsCNxE05pFgAA6CeCF8SW5tI535Kh0ECsda6iojQrDNkqAAC6g+CloNKQReknbsKRFLJVAAB0D8FLQaU5i5JmvIMOXwTKAAB0D8FLgeUli5IU3kHPJ9+AlGAUAID+IXjJAfagSQbvoIdLa8bKNyAlGAUAoL8KGbzk7WafErB0S+uNe1LSnLHyDUgJRgEA6K9CBi95vNmnBCyd0nzjnhQyVgAAoFsKGbxI3OwjGdy4AwAAdE9hg5e02qukbXZ2dtfNbZpL2gAAAIBuIXhJmTyWtMFfpz6ZqABWylefDDp/DczNzSW8GgAA0oHgJYU6lbTNz89T7pZz9MmANyMAAIhG8AKkzF59MlEBbB77ZNplHtplnqRksk9JTI+L8zUAAEARZSZ4mZmZod8DKIi0Zp/Sui4AAIpiz+DFzO6Q9D5JL5C0I+mUc+4BM5uU9EFJRySdl/SDzrmF3i01Gv0eQEOe9pPplHloVzqZRPaJ6XEAAPRXnMzLlqSfcs591swOSHrYzD4h6Uck/ZFz7l1m9nZJb5f0tl4t9OjRo5qdnc1Uv0feNsNEepERAAAARbBn8OKc+6qkrzb/f9nMnpB0m6Q3SHpV82HvlfSgehi8ZBGTw5AUMgIAAKAIvHpezOyIpJdLekjSVDOwkXPuq2bGHXcENsMEAAAAuiN28GJmo5J+R9JPOOeWzCzucSclnZSkqakpTU9Ph6xTkrS4uLiv41tmZ2cl9X5iTy/O061rkFVFf/1S9r4Puq3d60/r6+HnAAAA3RMreDGzshqBy/udc/+t+eF5M7ulmXW5RVJkc4dz7pSkU5J0/Phxd+LEieDFTk9P68bjQ/pKRkZGEsmItMpyunmeqGtQJEV//VL3rkEvvj6T0O71nz17NqhnrNdDC/g5AABA98SZNmaSfkPSE865/3Ddp35P0pslvav55+/2ZIV7oK8ECJenCWUMLQAAIP/iZF5eKelNkh4zs9Zdzs+pEbR8yMzeIukpST/QmyXujb4SwF/ebvbTPLQgT0EiAAD9FGfa2CcltWtweU13l5NejD1G3qT5Zj9P8hYkAgDQT17TxoqM8jQAIQgSAQDoHoIXD5SnAQAAAP1D8ALAC/0bAACgXwheAMRG/wYAAOgnghf0VMigg05mZ2fbvuvPu/u9l8f+Dd9MEl9nAAD0D8ELYgsJROr1uiRpfHy8F0t63nnq9brX+rgJhW8miSxStpjZHZLeJ+kFknYknXLOPdDfVQEA9oPgBbGFTFwbHx/vapAwPz8fOTTBN7DiJhSSfyYp7Vkk7LIl6aecc581swOSHjazTzjnvtDvhQEAwhC8FNR+9q1J48S1kJtQGs+BfHPOfVXSV5v/v2xmT0i6TRLBCwBkVCGDFzacZN8aGs+TRaDor9M1m5ubS3g12WdmRyS9XNJD/V1Jb/l+r507d06VSiXoXPXTp3XWdu9hXZH0YMTjlyW97qGHpHvuCTofAEgFDV6KfuPektYsShJCG8+5CfdHoOgvbz9r+s3MRiX9jqSfcM4tRXz+pKSTkjQ1NaXp6engcy0uLu7r+P24fPnyc32GcV28eFGbm5teAfEdzT/HJb3K62zSpXvv1SP33ed5VLb082sgLYp+DYr++qXeXoNCBi9Svm7cySQlg5vwMHmcUNZre12z+fn5BFeTbWZWViNweb9z7r9FPcY5d0rSKUk6fvy4O3HiRPD5pqentZ/jk9b6XvP5ffjFsTEtvvvdOhiRdZGkK1euaHJy8nkfq58+rXFJh6vVTF2fEFn7GuiFol+Dor9+qbfXoLDBS56QSUoGN+FAtpiZSfoNSU845/5Dv9eTF6vHjunLDzzQNuD5TMRNy1kz7ywNAETJTPAyMzMTuccHGYSGPGWSgDSjdDBTXinpTZIeM7PWP9rPOef+oI9rAgDsQ2aCl3bIICDtuNnND0oHs8U590lJ0bVNAIBMykzwcvToUc3OzuY+u7BX/wrZp2zhZjdfKB0EAKC/MhO8FAX9K/nCzW64dhmrqAC+hYwVAAD5RvCSQp36V9rtMA/kCRkrAAAQheAFSKFOfTLtMg95yjp0yli1C+DJWAEAkH8EL0DKkHUAAACIRvDSQ2weiRBxNii8MfNA1gEAABTBQL8XkGet5nsfNN8DAAAA0TKfeUl7doPNI5EU9pNJL/5tAADojswHL4wWBuiTaUljkMC/DQAA3ZP54EUiuwGwn0x6gwT+bQAA6J5cBC8AQJAAAED+0bAPAAAAIBPIvMSU9sEAQIg09ogAAAC0Q+YlJsYeI28OHz7sHVyvrKx4B/EAAADdQubFA4MBkCf0iAAAgKwheAFQaL6lc0mVzXVa19zcXM/PDwBAGhG8ACgs37LOpPZfodwUAIBoBC8AvOSpyd+3dC6psrm91jU/P5/IOgAASBuCFwCxpXUjSAAAUAwELwBio8kfAAD0E6OSAQAAAGQCwQsAAACATChk2djc3Jz3RnsrKyveG/rlTch1C5XWBm8gTwMLAADImswELzMzM5qdnd110xASVFy6dMn7uNHR0cKPLw25biHq9brq9XpkoBT1NdDCDWJ65eWGn4EFAAD0V2aCl3ZCg4rR0VHdfffdPVhRNuwn+9Tr6xa6NokbxDTK0w0/Awvw/7d390GSXWUdx78PASFkKYdN2K3wuuCCVETlZQmvsUJEQhArpEpeIiUBKYIUSAJlFSm1BEFLFASKAtEFYniHIG8JIgMFBIKLYUMckgAuzFARlknNEtZZM2SVIjz+0XfCZLd7Zrunu+/L+X7+mZ7b093nnL17z/31OfdcSVK9WhNedu7cycLCQtGBY5yaPPq03gni0tJS331gbm6usXdKL50n/JIkaVxaE140WJNHUaalqXdKlyRJ0vgYXjqgyaMo0zLKndK7ch2GJElSKQwvHdGlUZRp6NJ1GJIkSaUwvKhIXofRfINGxlxxTpKkchleGmaj61fGtVy0RuNUs+no2siY+40kSeNheGkYr19prq6dUDfZqCvONZH7jSRJ42N4aaD1rl8ZdOKmyXOqmUbhfiNJ0vgYXqQJc8qQJEnSeLQ+vGzmHifSpDllSJIkaXxaH168RkRN5pSh6XKUS5Kkbmt9eIHm3uPEUSGNar2T8EFLBZd+Eu4olyRJ3deJ8NJUjgppFJ6Ej8ZRLkmSus/wMmFNHRVSc210Et5vxTlPwrtnvdG3xcXFKZdGkqRmMLxIUsM4+ipJUn+GF6kjvFi9O45l9E2SpBJtGF4i4mLgqcCBzHxItW0r8CFgB3AD8IzM/O/JFVPSerxOZnTDhj4DnyRJ9TmWkZdLgLcA716z7SLgc5n52oi4qPr9FeMvnqRj4cXqoxk29Bn4JEmq14bhJTO/FBE7jth8NnB69fhdwBV0PLy47LHUPcOGPgOfJEn1GvWal+2ZeSNAZt4YEQO/voyI84HzAbZv387s7OyIHwmHDh066vULCwvA5OeALywscPjwYY4//vihXjczMzPWsvVrg5KUXn8YXxus7tN79uwZ6nUzMzOceOKJm/78UdW5D0zreLMR/x9Ikko18Qv2M3M3sBtg165deeaZZ478XrOzsxz5+tVvQie9HPG0Pmcj/dqgJKXXH8bXBpsZTazz/0Gd+4DHAUmS6jVqeFmKiJOrUZeTgeHOgCTVzutkJElS24waXi4DzgNeW/38xNhKJEkN5pLUkiTV51iWSv4AvYvzT4qI/cAr6YWWSyPi+cD3gKdPspCSmqPkk3eXpG6ffsv9S5La61hWGzt3wFO/OeaySGq40k/enWrXSpdw9HL/kqSWmvgF+5K6w5N3tc2A5f4lSS3VmvAyPz/PwsLCUSdC3ktFkiRJKkNrwssgW7ZsGXoqizeclKar5Otk1HyTvh9Zk41y76KNXtOvDe5c/Tx48CB7W9Q+o2jbPjAJpbdB6fWHybZBa8LLzp07WVhYGMv9FQ4cODB0GBklJEnyOhk136TvR9Zko9y7aKPX9GuDK6qfW7dubVX7jKJt+8AklN4GpdcfJtsGrQkv41b3jfakUnidjCNPkiSNyx3qLoAkddm2bduGnnK6srIy9NRW9Vct9/8V4JcjYn+1xL8kqaWKHXmRpGlw5Kle6yz3L0lqIcOLpEYaNNWq36qDq5xqJUlStxleJDWOF/lLkqR+DC+SGme9qVZLS0t9F9twqpUkSd1neCnUKPe6GZVTeaThrbdC2eLi4pRLI0lSMxheOmCUILK8vAzAzMzMJIp0u89ZXl4eunwGHo2iK0sSe08pSZL6M7x0wCg33ZyZmZnKSduowWpQ4PFibQ3SpetkNlqhbJi7oUuS1CWGl4bZ6GS/38n7anBp4k03R1kmdpTA09STUE2PSxJLktR9hpeGGWUUZcuWLZ2aZjLqxdrDThlypEaSJKldWh9eRv2Wftg7Xk/TeqMog07eSzdsePNaHEmSpPZpfXhp8khFF4NVUw07ZcipaVrliJ0kSe3R+vAC649U1KnJwap0o14f0ZXVrNQz7P81A6wkSfXqRHhpsqYGKw2vS6tZqWfYEOsF/pIk1cvwIh0jV7OSJEmql+FFmjCnmkmSJI2H4UWaIKeadY9hVJKk+hhejpErh2kUk1gYoN+NSsET5GkwjEqSVC/DyzFy5TBNiyfIzeV1T5Ik1cvwMgRXDtM0bHSC3O9GpS7jLEmSSmB4kTrA0ZpmM1hKkjQehhepA7zpZnMZLCVJGh/Di1QoT6qnw+tkJEkaH8OLVChPqiVJUtsYXiQNxalmkiSpLkWGF+/ZIo3GqWbTs15IXFxcnHJpJElqhiLDi/dskUbjwgDT4bFGkqT+WhNe5ufn+95ZfNQREe/ZIk2HozXDO5Z7/UiSVKLWhJdBHBGRms2FASRJ0ri0Jrzs3LmThYUFR0ukQgyaatZvBHZVyVPNJEkqQWvCi6RyONVMkiT1Y3iR1DjrTTVbWlrqOwLrwgCSJHWf4UVSJzhaI0lS9xleJHWCCwNIktR9hhdJRRt2qpnTzCRJqo/hRVKxhp1q5jQzSZLq1frwsri4yIEDB4Z6zag3tpTULcNONXNRAEmS6nWHuguwWQcOHLjt29Bj5Y0tJY1i27ZtQ3/xsbKyMvQXLJIkqb/Wj7xAL4x480pJk+aiAJIk1asT4UWSmsypZpIkjYfhRZImyPvPSJI0PoYXSZqgUaearTdas7i4OI6iSZLUOoYXSWoYFxSRJKk/w4skNcxGozVLS0tTLI0kSc3R+qWSJUmSJJXB8CJJkiSpFQwvkiRJklrB8CJJkiSpFQwvkiRJklphU+ElIp4cEfsiYj4iLhpXoSRJGgf7KUnqlpHDS0QcB7wVOAs4BTg3Ik4ZV8EkSdoM+ylJ6p7NjLycCsxn5ncz8yfAB4Gzx1MsSZI2zX5KkjpmMzepvBfw/TW/7wcetbniDDY/P8/CwgJzc3O3276yssKWLVsm9bGSpPaaaj/VVisrK0f1revZt28ft9xyC/v27ev7/NzcHMvLy7fbtr36ufzFLzIXMWpRW+HOwBV1F6JmpbdB6fWHXhtwwQXwpjeN/b03E176HX3yqD+KOB84H2D79u3Mzs6O9GGLi4usrKywZ8+eo56bmZkp5o7Thw4dGrkNu6D0+oNtUHr9wTYYwlT7KWjfv82PfvSjo4LGRpaXl7n55psHPn/48OGjwtBjq58zwOnDFVFSSx38/OfZO4Hj4WbCy37gPmt+vzeweOQfZeZuYDfArl278swzzxz5A2dnZ9nM67ug9DYovf5gG5Ref7ANhmA/VYO+bXDOOfCa18A6oacrDh48yNatW+suRq1Kb4PS6w9VG5xxxkSOh5sJL3uBB0bE/YEfAM8Cfm8spZIkafPsp5ri1FPh8svrLsVU7DXAFt8GpdcfJtsGI4eXzPxpRLwEmAWOAy7OzG+MrWSSJG2C/ZQkdc9mRl7IzE8BnxpTWSRJGiv7KUnqlk3dpFKSJEmSpsXwIkmSJKkVDC+SJEmSWsHwIkmSJKkVDC+SJEmSWsHwIkmSJKkVDC+SJEmSWsHwIkmSJKkVDC+SJEmSWsHwIkmSJKkVDC+SJEmSWsHwIkmSJKkVDC+SJEmSWsHwIkmSJKkVDC+SJEmSWiEyc3ofFvFD4L828RYnATeNqThtVXoblF5/sA1Krz9svg3ul5n3GFdhusR+aixKb4PS6w+2Qen1hwn2U1MNL5sVEVdn5q66y1Gn0tug9PqDbVB6/cE2aDL/bWyD0usPtkHp9YfJtoHTxiRJkiS1guFFkiRJUiu0LbzsrrsADVB6G5Ref7ANSq8/2AZN5r+NbVB6/cE2KL3+MME2aNU1L5IkSZLK1baRF0mSJEmFakV4iYgnR8S+iJiPiIvqLk8dIuKGiLguIuYi4uq6yzMNEXFxRByIiOvXbNsaEZ+NiO9UP+9eZxknaUD9XxURP6j2g7mIeEqdZZy0iLhPRHwhIr4VEd+IiAuq7UXsB+vUv6j9oC1K76vsp27bVsTxaVXpfZX91PT7qcZPG4uI44BvA78F7Af2Audm5jdrLdiURcQNwK7MLGbd8Ij4DWAFeHdmPqTa9rfAwcx8bXVycPfMfEWd5ZyUAfV/FbCSma+vs2zTEhEnAydn5jURcTfga8DTgOdSwH6wTv2fQUH7QRvYV9lPldhPgX2V/dT0+6k2jLycCsxn5ncz8yfAB4Gzay6TpiAzvwQcPGLz2cC7qsfvovcfpJMG1L8omXljZl5TPb4Z+BZwLwrZD9apv5rHvqpApfdTYF9lPzX9fqoN4eVewPfX/L6fMjvvBD4TEV+LiPPrLkyNtmfmjdD7DwNsq7k8dXhJRFxbDdV3chi6n4jYATwMuIoC94Mj6g+F7gcNZl9lP7WquOPTAMUdo+ynYgdT6KfaEF6iz7Zmz3WbjMdl5sOBs4AXV8O0Ks/bgF8CHgrcCPxdvcWZjojYAnwEuDAz/6fu8kxbn/oXuR80nH2V/ZR+rrhjlP3U9PqpNoSX/cB91vx+b2CxprLUJjMXq58HgI/Rm6JQoqVqfuXqPMsDNZdnqjJzKTNvzcyfAW+ngP0gIu5E74D4vsz8aLW5mP2gX/1L3A9aoPi+yn7qNsUcnwYp7RhlPzXdfqoN4WUv8MCIuH9E/ALwLOCymss0VRFxQnURFBFxAvAk4Pr1X9VZlwHnVY/PAz5RY1mmbvVAWDmHju8HERHAO4FvZeYb1jxVxH4wqP6l7QctUXRfZT91O0Ucn9ZT0jHKfmr6/VTjVxsDqJZXexNwHHBxZv5VzUWaqoh4AL1vsQDuCLy/hDaIiA8ApwMnAUvAK4GPA5cC9wW+Bzw9Mzt5oeCA+p9Obwg2gRuAF67Oqe2iiHg8cCVwHfCzavOf0JtP2/n9YJ36n0tB+0FblNxX2U+V2U+BfZX91PT7qVaEF0mSJElqw7QxSZIkSTK8SJIkSWoHw4skSZKkVjC8SJIkSWoFw4skSZKkVjC8qJMiYmWKn3V6RDx2Wp8nSZqMiLg1IuYi4vqI+HBE3HWE97hw0OvWe26I998REUPdMyMi9mzmM4f8LPtETZThRdq804GhDtQRccfJFEWStAmHM/OhmfkQ4CfAH47wHhcCgwLKwOci4rgRPuuYZOY0w8Tp2Cdqggwv6rTqG6ArIuKfI+I/I+J90XNWRFx6xN9dXj1+UkR8JSKuqb5521JtvyEi/qLafl1EPDgidtDr3F5WfVt3WkTcLyI+FxHXVj/vW73+koh4Q0R8AXhdRHwnIu5RPXeHiJiPiJOm3ESSpP6uBHYCRMTLq9GY6yPiwmrbCRHxLxHx9Wrkr1oBAAADoElEQVT7MyPipcA9gS9Ux/rb9HsuIlYi4tURcRXwmIj484jYW73f7uru5UTEI6rP+Qrw4jXveVxEvK56zbUR8cJ+FVmdjWCfqC4wvKgED6P3bdcpwAOAxwGfBR4dESdUf/NM4EPVgfLPgCdm5sOBq4GXr3mvm6rtbwP+ODNvAP4BeGP1bd2VwFuAd2fmrwHvA9685vUPqt77ZcB7gWdX258IfD0zbxpv1SVJw6pGAs4CrouIRwDPAx4FPBp4QUQ8DHgysJiZv16N1Hw6M98MLAJPyMwnrH3PAc+dAFyfmY/KzC8Db8nMR1bvdzzw1Orv/gl4aWY+5oiiPh84lJmPBB5Zle3+G1TPPlGtZnhRCb6amfsz82fAHLAjM38KfBr4naqT+m3gE/Q6plOAf4uIOeA84H5r3uuj1c+vATsGfN5jgPdXj98DPH7Ncx/OzFurxxcDz6ke/wG9zkmSVJ/jq2P/1cD3gHfSO4Z/LDN/nJkr9PqB04DrgCdGxN9ExGmZeWiEz7sV+Mia358QEVdFxHXAGcCvRMQvAjOZ+cXqb96z5u+fBDynKvNVwInAAzf4TPtEtZpzDFWC/1vz+FZ+vt9/iN7w+0Fgb2beXA3RfzYzz93gvda+z0ZyzeMf37Yx8/sRsRQRZ9D7Ru/ZR71SkjRNhzPzoWs3rE7dOlJmfrsalXkK8NcR8ZnMfPWQn/e/qyfvEXEX4O+BXVX/8CrgLkBw+37kdsUD/igzZ4f4TPtEtZojLyrZFcDDgRfQO2gD/DvwuIhYned814h40AbvczNwtzW/7wGeVT1+NvDldV77DnpD5Zeu+fZJktQcXwKeVvUHJwDnAFdGxD2BWzLzvcDr6fUncHSfsNZ6z92l+nlTdV3J7wJk5jJwKCJWRyzWntTPAi+KiDsBRMSD1kz9GtYV2CeqBQwvKlZ1YPwkvXnNn6y2/RB4LvCBiLiW3oH7wRu81eXAOasXJwIvBZ5Xvf73gQvWee1lwBYcHpekRsrMa4BLgK/Sm5r1jsz8D+BXga9W06n+FPjL6iW7gX898oL9jZ6rQsrb6U1H+ziwd83TzwPeWl2wf3jN9ncA3wSuid7yyf/IiLNq7BPVFpE5aCRS0qRFxC56FzaeVndZJEmqk32ijoXXvEg1iYiLgBfhvF5JUuHsE3WsHHmRJEmS1Ape8yJJkiSpFQwvkiRJklrB8CJJkiSpFQwvkiRJklrB8CJJkiSpFQwvkiRJklrh/wG0EkL6irjb6QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0 1 2 3 "
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAy8AAAHwCAYAAACv2lDjAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dfZhkZXnv++8tL85AjQwTqA6N4DgDssP2xFY7o9lGg29poyaYRI3uxKBxOzs5O6KJJpI3FUl2jGFH4JidnIkaIeILGo2iiS1bxZeDjiPYiopRxkEgjVM4YzdTMiMB7/NHrcZmqO6u6q7q6lXr+7muubpqVa217vVQ1NO/ftZ6VmQmkiRJkrTWPWDQBUiSJElSJwwvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8KLSi4jNEZERcWQftn1GRHwxIg5ExLm93v4i+z01IpoRccRq7VOS1J1+fldHxGsj4u293u4gFX31aYOuQ+VmeNHARcRkRLyuzfKzI+I7/QglXfgD4OrM3JCZl/RrJxFxU0Q8Ze55Zt6cmbXMvKdf+5SkqomIF0bE9RFxZ9G//G1EbOxifb+rpQEzvGgteBvwgoiIw5a/ALg8M+9e/ZLu9RDgqwPcvySpByLiFcBfAr8PHAc8ltZ3/FURcfQgaxs2vf6j44D/iKk1xvCiteCfgU3A4+cWRMTxwDOBy4rnzyhO37ojIm6JiNcutLHD/zJ2+NB7RDw2Iq6JiJmI+FJEnLXAdj4OPBF4U3FawMMi4uqI+G/z3vPCiPjMvOcZEb8VEd+MiO9FxN/MD2UR8ZKIuKE4De1rEfGoiPhH4FTgymI/f3D4qXARMRoRH4yI/RFxY0S85LDjuyIiLiu2+9WIGF+62SWpGiLiQcD5wEsz8yOZ+R+ZeRPwXFoB5teL9702It4bEe8uvk+vi4hHFK918l19dUT8WdHHNCPiyoj4sYi4vOi/dkXE5nl1XVz0aXdExLUR8Xg6EBFnRcStEfGKiGhExG0R8aJ5rz8wIi6MiJsjYm9E/F1ErC9e+2RE/Erx+GeK+p9ePH9KREwtsM8HRsRFETFd/LsoIh54WD2viojvAP9QLP/9orbpiPjNNttbqMb7bS8iToiIDxV99/6I+HRE+HtsBfkfXQOXmQeBK4DfmLf4ucDXM/NLxfPvF69vBJ4B/HZEPKvbfUXEycCHgT+jFZheCfxTRJzYpq4nAZ8Gfqc4LeAbHe7mmcBPAY8ojmOi2PdzgNcWx/Eg4BeBfZn5AuBm4BeK/byhzTbfCdwKjALPBv5nRDx53uu/CLyLVvt8EHhTh7VKUhX8F2Ad8L75CzOzCfwr8NR5i88G3kOrj3gH8M8RcVSH39UAz6N15sDJwFbgs7R+md8E3AC8Zt57dwFj8/b1nohY1+Ex/TitEaSTgRcDfxOtP/xBa4TpYcW2Tyve8+ritU8CZxWPnwB8C/jZec8/ucD+/pjWaNUYrf5tG/Anh9WziVYY3B4RT6PVxz4VOB14Cve1WI332x7wClr94InACPBHQC5Qq4aY4UVrxaXAc+b+6kLrF/xL517MzKsz8/rM/GFmfpnWL/M/22Y7S/l14F8y81+KbV0FfAF4+grrn+/1mTmTmTcDn6D1xQzw34A3ZOaubLkxM7+91MYi4hTgZ4BXZeahzJwC3kyrc5zzmeKY7gH+kVbHIklqOQH47gKnId9WvD7n2sx8b2b+B/DXtELPY7vY1z9k5u7MnKUVjHZn5v8p9v0e4JFzb8zMt2fmvsy8OzP/F/BA4IwO9/MfwOuKUaR/AZrAGcVo/0uA383M/Zl5APiftEIVtMLJ/LDyF/Oe/ywLh5dfK/bXyMzbaY1kze+Hfgi8JjN/UPxR8rlFW3wlM79P6493AHRQY7vt/QdwEvCQ4pg/nZmGlwoyvGhNyMzPALcDZ0fEFlojF++Yez0iHhMRn4iI2yNiFvgt7tvZdOohtELSzNw/WsHgpJUfxb2+M+/xnUCteHwKsHsZ2xsF5r7c53yb1l+pFtrnuvAcYUma813ghAW+F08qXp9zy9yDzPwhPxr17tTeeY8Ptnk+1ydQnPZ1Q0TMFv3RcXTet+07LIzN9TcnAscA187r5z5SLIfWSNDDImKE1h/XLgNOiYgTaI2mfGqB/Y3S6nvmfJv7tsvtmXnosPffctj75yxVY7vt/RVwI/DRiPhWRJy3QJ0acoYXrSWX0RpxeQHw0cyc/4X/DlqnQ52SmccBfwccfoH/nO/T+lKc8+PzHt8C/GNmbpz379jMfH2HNS627aXcQusUgnYW++vRNLApIjbMW3Yq8O9d7FuSquyzwA+AX56/MCKOBX4e+Ni8xafMe/0BwINpfQ9DD09TKq5veRWtEYrjM3MjMMvCfVunvksrJP3nef3ccZlZA8jMO4FrgZcBX8nMu4BrgN+jNUr03QW2O03rD4BzTuVH7QL3b5vbmNeWxfs7qrHd9jLzQGa+IjO3AL8A/N5hp0+rIgwvWksuo3VO7EuYd8pYYQOt0YdDEbEN+K+LbGcKeF5EHFVcuP7sea+9HfiFiJiIiCMiYl1xYeCDO6xxCvjliDgmWnPVv7jD9aB1qtcrI+LR0XJaRMx1BHuBLe1WysxbaHUsf1HU+5PFfi/vYt+SVFnFKVznA/9PRDyt6B820zqN61Zap9vOeXRE/HIxSvNyWqHnc8VrC35XL8MG4G5aZx0cGRGvpnU95IoUo0V/D7wxIurQut4zIibmve2TwO/wo1PErj7seTvvBP4kIk4sRmleTatPXcgVwAsj4syIOIZ51/p0WON9RMQzi34zgDuAe4p/qhjDi9aMYuaXa4BjaY2yzPd/A6+LiAO0vjCvWGRTf0prhON7tDqre08/K4LA2bQu9Lud1mjI79P5/wtvBO6i1YFdShcBIjPfA/x5Uc8BfjTLGrTOOf6TYvj8lW1Wfz6wmdZfud5P6zzgqzrdtyRVXXGB/R8BF9L65XcnrT7gyZn5g3lv/QDwq7T6kBcAv1xc/wJLf1d3Y5LWNTHfoHVK1SHue5rVSryK1ilWn4uIO4D/w32vpfkkrfD0qQWet/NntK4R/TJwPXBdsaytzPxX4CLg40UtH++yxsOdXrynSWsk7X9n5tWLvF9DKrzWSZIkqTVVMnBaZv76oGuR1J4jL5IkSZJKwfAiSZIkqRQ8bUySJElSKTjyIkmSJKkUDC+SJEmSSmFV78B9wgkn5ObNm5e9/oEDB9iwYcPSbxxiVW+Dqh8/2AZVP35YeRtce+21383ME5d+Z/XYT61c1dug6scPtkHVjx/620+tanjZvHkzX/jCF5a9/uTkJBMTC96/qBKq3gZVP36wDap+/LDyNoiIb/ewnKFiP7VyVW+Dqh8/2AZVP37obz/laWOSJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSuHIQReg4fb1yy5jw8UXc9TBgz3Z3unNJo1arSfbeuAJJ3DchRfCtm092Z6karvxxhvZvXs3U1NTbV+v1+uMjo6uclWSNFwML+qrDRdfzMnXXdez7dV7tqXCYx5D4yd+ouO3G3gkLUez2QQwvEjSChle1LHljKKccNNNrQfHHQdjYyuuYf/+/WzatGnF2+GTn7z3Yf2GG7pb94IL4MorV16DpKFy2mmnsXv3bsbafNctNBojSeqO4UUdW9EoyuMf35Nf+HdNTjIxMbHi7fD5z7dCyIEDHa9y165dHH3nndz18Y8zc+aZHa/naI0kSVJvGF4qalVHUTZsgFe/ursC+23btq7D1MEnPIGjP/1pjr7zTkdrJEmSBsDwUlFrYRSlbI678MJlj9bcc/XVHBwf73i9I48/nnV//ueO1kiSJM1jeKmoe0dchmEUZbUsY7Rm9rGP5cSdOzmi2aR27bXd7W/dukqGREmSpIUYXobA7R/+MEe/4Q0c8f3vd7zOcXOngI2NwdVX96UuAX/6p8x2+d/mqK99jQcePNjVCI8kSVIVGF6GwQUXcNzOnctbd8OG3tai+zjxGc+AZzyjq3UaZ57Z/TU1kiRJFWB4GQJ5xx0A3FOrcfCMMzpe78jjj2ddVU8BK4G7du1qO6vZQjfqdFYzSZI07AwvQ+SIRz+amqeAld4DTzgBYMFZzRa9UaezmkmSpCHWUXiJiI3Am4GHAwn8JvBvwLuBzcBNwHMz83t9qVKqkKVmNWt3o865Wc28TkaSJA2zTkdeLgY+kpnPjoijgWOAPwI+lpmvj4jzgPOAV/WpzspY6v4r7U4Zuvfiew2HJWY1a3ejzhmvk5EkSRWwZHiJiAcBTwBeCJCZdwF3RcTZwFnF2y4FrsbwsmJL3X9l0VOGvPi+8ha6TmYhXicjSZLKpJORly3A7cA/RMQjgGuBlwEjmXkbQGbeFhFtf6+OiO3AdoCRkREmJyeXXezs7OyK1i+D02+/HYD/OPZYDmzZcr/X777nHo484oj7Lz/mGHY/9ancMeTtU4XPwFLatcHDis/EQtfJLKZx7rl88fzze1Zfv/kZsA0kSdXVSXg5EngU8NLM3BkRF9M6RawjmbkD2AEwPj6eh5/u0o3JNqfLDJu5U8KOGh9nU5uL7xdrg0VHZYZEFT4DS2nbBscfv+h1Mu3MXSdTX7euVG3qZ8A2kCRVVyfh5Vbg1sycu5HIe2mFl70RcVIx6nIS0OhXkZKWsMR1Mu14nYwkSSqbJcNLZn4nIm6JiDMy89+AJwNfK/6dA7y++PmBvlZaQktdfN+OF99LkiRJ7XU629hLgcuLmca+BbwIeABwRUS8GLgZeE5/SiyvpS6+X3xlL77X6vAif0mSVBYdhZfMnALG27z05N6WM1zuHXE57jgYG+t8xQ0bwDvfq8+WuhnmorwZpiRJGoBOR160EmNj0Obie2mQlroZZjveDFOSJA2S4UWqKi/ylyRJJfOAQRcgSZIkSZ1w5EVS17zIX2tJRLwVeCbQyMyHH/baK4G/Ak7MzO8Ooj5JUu8YXjrktMeSF/lrzXob8CbgsvkLI+IU4Km0ZsSUJA0Bw0uHnPZY8iJ/rU2Z+amI2NzmpTcCf4D3IZOkoWF46ZDTHkt4kb9KIyJ+Efj3zPxSRAy6HACazSZTU1Mdv79erzM6OtrHiiSpfAwv3XLaY0la0yLiGOCPgZ/r8P3bge0AIyMjTE5OLnvfs7Ozbdfft28fMzMzHW/n4MGDrF+/nq1bty67lkFZqA2qourHD7ZB1Y8f+tsGhhdJ0rDZCjwUmBt1eTBwXURsy8zvHP7mzNwB7AAYHx/PiYmJZe94cnKSlaw/Z26EZqybkf41oldtUFZVP36wDap+/NDfNjC8SFoVzlCm1ZKZ1wP1uecRcRMw7mxjklR+hhdJfeUMZeq3iHgncBZwQkTcCrwmM98y2KokSf1geJHUV85Qpn7LzOcv8frmVSpFktRnhhdJ/eUMZZIkqUcqGV5u//CHOfoNb+CI73+/43W84aQkSZI0WJUML1xwAcft3Lm8db3hpCRJkjQQlQwveccdANxTq3HwjDM6Xu/I449nnTeclCRJkgaikuFlzhGPfjQ1bzgprVntplc+vdmkUau1fb/TK0uSNNwqHV4krU2LTa9cb7fCfE6vLEnS0DK8SFpzFpteef/+/WzatOl+y51eWZKk4Wd4kbT2LDK98q7JSSYmJu633OmVJUkafg8YdAGSJEmS1AnDiyRJkqRSMLxIkiRJKgWveZE0VNpNr7wYp1eWJKk8DC+ShsJi0ysvyemVJUkqBcOLpKGw2PTKC3F6ZUmSyqX04eXrl13Ghosv5qiDBzte57ibbupfQZIGY5HplRfi9MqSJJVL6cPLhosv5uTrrlvmyht6W4wkSZKkvil9eLl3xOW442BsrPMVN2yAV7+6P0VJkiRJ6rnSh5d7jY3B1VcPugpJkiRJfeJ9XiRJkiSVwvCMvEjSMnVzbxjvCyNJ0uAYXiRV1rLvDeN9YbQKms0mU1NTXa1Tr9cZHR3tU0WSNHiGF0mV1e29YbwvjFZLvV7vep1mswlgeJE01Awvkqqry3vDeF8YrZbR0dGuQ0i3ozSSVEZesC9JkiSpFAwvkiRJkkrB8CJJkiSpFAwvkiRJkkrBC/YlqUvd3BcGvDeMJEm9YniRpA4t+74w4L1hJEnqAcOLJHWo2/vCgPeGkSSplwwvktSpLu8LA94bRpKkXipNeLn5V36F03ftolGr3Wf5cTfdNJiCJEmSJK2q0oSXDTfeyPG33LLIGzasXjGSJEmSVl1pwsvxT3wi+zPZtGnT/V/csAFe/erVL0qSJEnSqilNeOGii9g1OcnExMSgK5EkSZI0AN6kUpIkSVIplGfkRZJKzBtbSpK0coYXSeojb2wpSVLvGF4kqY+8saUkSb1jeJGkfvLGlpIk9YwX7EuSJEkqBcOLJEmSpFLo6LSxiLgJOADcA9ydmeMRsQl4N7AZuAl4bmZ+rz9lSpIkSaq6bkZenpiZY5k5Xjw/D/hYZp4OfKx4LkmSJEl9sZLTxs4GLi0eXwo8a+XlSJLUnYh4a0Q0IuIr85b9VUR8PSK+HBHvj4iNg6xRktQbnYaXBD4aEddGxPZi2Uhm3gZQ/Kz3o0BJkpbwNuBphy27Cnh4Zv4k8A3gD1e7KElS73U6VfLjMnM6IurAVRHx9U53UISd7QAjIyNMTk4uo8yW2dnZFa0/DKreBlU/frANqnD8pzeb1IGDn/sct5166v1ef1gE337Uo/j6b/3W6he3BmXmpyJi82HLPjrv6eeAZ69mTZKk/ugovGTmdPGzERHvB7YBeyPipMy8LSJOAhoLrLsD2AEwPj6eExMTyy52cnKSlaw/DKreBlU/frANqnD8s5s3wy23sP4HP2DLLbe0f9NDH8pDhrwdeug3aU0wM/SazSZTU1NdrVOv1xkdHe1TRZLUW0uGl4g4FnhAZh4oHv8c8Drgg8A5wOuLnx/oZ6GSVBXHXXghXHABHDjQ9vX9+/ezaWxslasqp4j4Y+Bu4PJF3jMUZwjs27ePmZmZrtY5ePAg69evZ+vWrT2rowqjo4up+vGDbVD144f+tkEnIy8jwPsjYu7978jMj0TELuCKiHgxcDPwnL5UKElVs20bXHnlgi/vqsDoUy9ExDnAM4EnZ2Yu9L4qnyEwN0oz1sMwXLY26LWqHz/YBlU/fuhvGywZXjLzW8Aj2izfBzy5H0VJkrQSEfE04FXAz2bmnYOuR5LUGyuZKlmSpIGLiHcCnwXOiIhbizMC3gRsoDXJzFRE/N1Ai5Qk9USns41JkrQmZebz2yx+y6oXIknqO0deJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRw56AIkSdLgNJtNpqamulqnXq8zOjrap4okaWGGF0mSKqper3e9TrPZBDC8SBoIw4skSRU1OjradQjpdpRGknrJa14kSZIklYLhRZIkSVIpGF4kSZIklYLhRZIkSVIpGF4kSZIklYLhRZIkSVIpGF4kSZIklYLhRZIkSVIpGF4kSZIklYLhRZIkSVIpGF4kSZIklYLhRZIkSVIpGF4kSZIklYLhRZIkSVIpGF4kSZIklcKRgy5AGmbT09M0Go2ebnP37t1MTU31ZFv1ep3R0dGebEuSJKnfDC9Sh5YTRGZmZgDYuHFjP0pakZmZGWZmZro+JgOPJEkaFMOLKmm1gsjGjRt7/sv+3r17GRsbW/F2ltMGzWYTwPAiSZIGwvCiodDtL+JrJYgM0ujoaNfHMjU1RbPZ7Pq0tWFqN0mSNDiGF605iwWRha736DaMDFsQWS31er3rdRytkSRJvWJ4UV+V+fQs3d9yR2skSZJ6wfCivmo0GjSbTWq1WsfrLBZEenW9h1aXp5qpnyLircAzgUZmPrxYtgl4N7AZuAl4bmZ+b1A1SpJ6w/Ciji33Au9arWbgqDBPNdMqeBvwJuCyecvOAz6Wma+PiPOK568aQG2SpB4yvKhjyxlFqdVqy/rlVcPDU83Ub5n5qYjYfNjis4GziseXAldjeOmZxUZTF7o20dFUSb1geKkoR1G01i30y9FiN+n0lyPNM5KZtwFk5m0RseBfUSJiO7AdYGRkhMnJyWXvdHZ2dkXrl8G+ffvuvTaxnWazyTXXXHOfZQcPHmT9+vVs3bq13+UNXBU+A0upehtU/fihv21geKkoR1G0lnmqmVZTZu4AdgCMj4/nxMTEsrc1OTnJStYfBu3aYO4PDlX445efAdug6scP/W0Dw0uFOYqitWqxU80WmrTBU810mL0RcVIx6nIS0N1QsyRpTTK8DIGVnAImSUPqg8A5wOuLnx8YbDmSpF54wKAL0MrNnQLWDU8BkzQsIuKdwGeBMyLi1oh4Ma3Q8tSI+Cbw1OK5JKnkHHkZEp4CJnk/marKzOcv8NKTV7UQSVLfdTzyEhFHRMQXI+JDxfOHRsTOiPhmRLw7Io7uX5mStLh6vd71qZDNZrPrUy4lSdLgdDPy8jLgBuBBxfO/BN6Yme+KiL8DXgz8bY/rq5ylrl9pN02s169I3k9GkqQq6GjkJSIeDDwDeHPxPIAnAe8t3nIp8Kx+FFg1Xr8iSZIktdfpyMtFwB8AG4rnPwbMZObdxfNbgZN7XFtlLXb9ykLTxEpanm6vk/EaGUmSBmfJ8BIRzwQamXltRJw1t7jNW3OB9b1zcRd2794NtEJKO1Vog8VU/fjBNujl8S91p/DDrZW7hFf9MyBJqq5ORl4eB/xiRDwdWEfrmpeLgI0RcWQx+vJgYLrdyt65uDtL3YW4Cm2wmKofP9gGgzz+tXKX8Kp/BiRJ1bVkeMnMPwT+EKAYeXllZv5aRLwHeDbwLrwBWFvePFKSJEnqnZXcpPJVwO9FxI20roF5S29KGh5efC9JkiT1Tlc3qczMq4Gri8ffArb1vqTh4s0jJUmSpN7oKrxIUtV1OzsZOEOZJEm9YniRpA4t55TOuVNHDS+SJK2c4UWSOjQ6Otp1COl2lEaSJC1sJRfsS5IkSdKqceSlQ057LEmSJA2WIy8dctpjSZIkabAceemC0x5LWg5nKJMkqTcML5LUR85QJklS7xheJKmPnKFMkqTe8ZoXSZIkSaVgeJEkSZJUCp42JkmS+s6JKyT1guFFkiT1lRNXSOqVSoYXbzgpaa1b7K/U09PTq1yNtDJOXCGpVyp5zYs3nJS0ltXrdf9YIklSG5UceQFvOClp7Vrqr9R79+5dxWokSVo7KjnyIkmSJKl8DC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSuHIQRewUtPT0zQaja7WaTab1Gq1PlUkSZIkqR9KP/LSaDRoNptdrVOr1ajX632qSJIkSVI/lH7kBVphZGxsbNBlSJIkSeqj0o+8SJIkSaoGw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skaWhFxO9GxFcj4isR8c6IWDfomiRJy2d4kSQNpYg4GTgXGM/MhwNHAM8bbFWSpJUwvEiShtmRwPqIOBI4BpgecD2SpBUYiptUSpJ0uMz894i4ELgZOAh8NDM/OuCy1IVms8nU1FRX69TrdUZHR/tUkaRBM7xIkoZSRBwPnA08FJgB3hMRv56Zbz/sfduB7QAjIyNMTk4ue5+zs7MrWn8Y9KoN9u3bx8zMTFfrHDx4kPXr17N169YV73+5/AzYBlU/fuhvGxheJEnD6inAnsy8HSAi3gf8F+A+4SUzdwA7AMbHx3NiYmLZO5ycnGQl6w+DQbbB3CjN2NjYQPYPfgbANqj68UN/28BrXiRJw+pm4LERcUxEBPBk4IYB1yRJWgHDiyRpKGXmTuC9wHXA9bT6vB0DLUqStCKeNiZJGlqZ+RrgNYOuQ5LUG6UJLzfeeCO7d+++36wjzWaTWq02oKokSZIkrZbSnzZWq9Wo1+uDLkOSJElSn5Vm5OW0005j9+7dA51BRJIkSdLglH7kRZIkSVI1GF4kSZIklYLhRZIkSVIpGF4kSZIklYLhRZIkSVIpGF4kSZIklYLhRZIkSVIpGF4kSZIklcKS4SUi1kXE5yPiSxHx1Yg4v1j+0IjYGRHfjIh3R8TR/S9XkiRJUlV1MvLyA+BJmfkIYAx4WkQ8FvhL4I2ZeTrwPeDF/StTkiRJUtUtGV6ypVk8Par4l8CTgPcWyy8FntWXCiVJkiSJDq95iYgjImIKaABXAbuBmcy8u3jLrcDJ/SlRkiRJkuDITt6UmfcAYxGxEXg/8BPt3tZu3YjYDmwHGBkZYXJycpmlwuzs7IrWHwZVb4OqHz/YBlU/frANJEnV1VF4mZOZMxFxNfBYYGNEHFmMvjwYmF5gnR3ADoDx8fGcmJhYdrGTk5OsZP1hUPU2qPrxg21Q9eMH20CSVF2dzDZ2YjHiQkSsB54C3AB8Anh28bZzgA/0q0hJkiRJ6mTk5STg0og4glbYuSIzPxQRXwPeFRF/BnwReEsf65QkSZJUcUuGl8z8MvDINsu/BWzrR1GSJEmSdLiOZhuTJEmSpEHr6oJ9SZLU3tcmJnjIl7/MzevXt319Xb1O/ZJLYJsnLfRTs9lkamqqq3Xq9Tqjo6N9qkhSLxleJEnqgU0338yPf+c7C79hzx644AK48srVK6pi6vV61+s0m637cBtepHIwvEiS1AM/PjHB/qOOYtOmTfd77dDOnaw7dAgOHBhAZdUxOjradQjpdpRG0mAZXiRJ6oWLLmLXAvfgaWzZwql79gygKEkaLoYXSR2bnp6m0Wisyr48B12SJB3O8CJV1HKCyMzMDAAbN27sR0n32c/MzEzb+nbv3r3gaZlmlUwAABnHSURBVB4GHkmShpvhRaqoRqNBs9mkVqt1vM7GjRtXJSAsN1gtFHgWY+CRJKk8DC/SEFjOL/tzwWVsbKxPVS3fYhfd7t27t23Ny22Duf1JkqS1z/AiDYHljKLUarVlTSu6Vi13lqFu7wnhSI0kSYNjeJGGxFodRVnLug1vjtRIkjRYhhdpjVnq9Kd2F6x3O+qilm5Ha7wfhCRJg2V4kdYYTwFb27o9zQw81UySpF4xvEhr0GKngC10wbr6bzkB0VPNJEnqHcOL1EcrmQVMa89yJwWQJEm98YBBFyANs7lTwLrhKWCSJEntOfIi9ZmzgMnrZCRJ6g3DiyT1kdfJSJLUO4YXSeojr5ORJKl3DC9Sh7z4XpIkabAML1KHvP+KVtNi18lMT0+vcjWSJK0NhhepC158r9Vg4JUkqT3DiyStMUtdJ7N3795VrEaSpLXD+7xIkiRJKgXDiyRpaEXExoh4b0R8PSJuiIifHnRNkqTl87QxVZIzh0mVcTHwkcx8dkQcDRwz6IIkScvnyIsqaW7msG44c5hULhHxIOAJwFsAMvOuzJwZbFWSpJVw5EWV5cxh0tDbAtwO/ENEPAK4FnhZZn5/UAUd2rmTxpYtHb13Xb1O/ZJLYNu2PlelxaYmX0i9Xu/6BrSSVs7wIkkaVkcCjwJempk7I+Ji4DzgT+e/KSK2A9sBRkZGmJycXPYOZ2dn266/5eijAVh36BCn7tnT2cb27KFx7rl88fzzl13PICzUBmvVvn37mJnpbkDu4MGDrF+/nq1bt97vtbIdfz9UvQ2qfvzQ3zYwvEiShtWtwK2ZubN4/l5a4eU+MnMHsANgfHw8JyYmlr3DyclJ2q5/2WVwwQVw4EBH2zm0cyfrDh2ivm5d++2tYQu2wRCZG6VpN3pfheNfStXboOrHD/1tA8OLJGkoZeZ3IuKWiDgjM/8NeDLwtYEUs20bXHllx29vbNnS+QiNJFWI4UWSNMxeClxezDT2LeBFA65HkrQChhcNhW6nPnbaY6kaMnMKGB90HZKk3nCqZA2Fbqc+dtpjSZKk8nHkRUPDqY8lSZKGmyMvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8CJJkiSpFLxgX2vOYtMe7969+947G8/n1MeSJEnDz5EXrTndTnsMTn0sSZJUBY68aE1aaNrjvXv3Oh2yJElSRTnyIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsHwIkmSJKkUDC+SJEmSSsGpktVXi91wciHecFKSJEntOPKivvKGk5IkSeoVR17UdwvdcFKSJEnqhiMvkiRJkkrB8CJJkiSpFDxtTJKkNerQzp00tmzp+P3r6nXql1wC27b1sSpBa3KZqamp+y3fvXt32+UA9Xqd0dHRfpcmDTXDiyRJa8y6eh327GHdoUOcumdP5yvu2QMXXABXXtm/4rSsSWXmJq8xvEgrs2R4iYhTgMuAHwd+COzIzIsjYhPwbmAzcBPw3Mz8Xv9KlSSpGuqXXNIKIQcOdLzOoZ07WXfoUFfraHlGR0cXDCF79+5tO0nNQqMxkrrTycjL3cArMvO6iNgAXBsRVwEvBD6Wma+PiPOA84BX9a9USZIqYtu2rkdPGlu2dDdKI0kltOQF+5l5W2ZeVzw+ANwAnAycDVxavO1S4Fn9KlKSJEmSurrmJSI2A48EdgIjmXkbtAJORHhXwSE3PT1No9Hoap1ms0mtVutTRZIkSaqSjsNLRNSAfwJenpl3RESn620HtgOMjIwwOTm5nDoBmJ2dXdH6w2CQbbB7924OHjzI+vXru1pv48aN7N27tyc1+BmwDap+/GAbSJKqq6PwEhFH0Qoul2fm+4rFeyPipGLU5SSg7Z/kM3MHsANgfHw8JyYmll3s5OQkK1l/GAyyDeYuNmx3IeJq8TNgG1T9+ME2kCRV15LXvERriOUtwA2Z+dfzXvogcE7x+BzgA70vT5IkSZJaOhl5eRzwAuD6iJib5++PgNcDV0TEi4Gbgef0p0RJkiRJ6iC8ZOZngIUucHlyb8uRJEmSpPaWPG1MkiRJktYCw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skSZKkUujkPi8aQtPT0zQaja7WaTab1Gq1PlUkSZIkLc6Rl4pqNBo0m82u1qnVatTr9T5VJEmSJC3OkZcKq9VqjI2NDboMSZIkqSOOvEiSJEkqBcOLJEmSpFIwvEiSJEkqBcOLJEmSpFIwvEiShlpEHBERX4yIDw26FknSyhheJEnD7mXADYMuQpK0coYXSdLQiogHA88A3jzoWiRJK+d9XiRJw+wi4A+ADYMuZLUc2rmTxpYtHb9/Xb1O/ZJLYNu2PlYlgGazydTUVFfr1Ot1RkdH+1SRVD6GF0nSUIqIZwKNzLw2Is5a5H3bge0AIyMjTE5OLnufs7OzK1p/JbYcfTQA6w4d4tQ9ezpfcc8eGueeyxfPP78ndQyyDdaChY5/3759zMzMdLWtgwcPsn79erZu3dqr8laFn4FqHz/0tw0ML0NgenqaRqPR1TrNZpNardaniiRpTXgc8IsR8XRgHfCgiHh7Zv76/Ddl5g5gB8D4+HhOTEwse4eTk5OsZP0VuewyuOACOHCg41UO7dzJukOHqK9b17O6B9oGa0Avj39ulGZsbKwn21stfgaqffzQ3zYwvAyBRqPRdRip1WrU6/U+ViVJg5WZfwj8IUAx8vLKw4PLUNm2Da68sqtVGlu2dDdKI0kDZngZErVarXR/mZEkSZK6YXiRJA29zLwauHrAZUiSVsipkiVJkiSVguFFkiRJUikYXiRJkiSVguFFkiRJUikYXiRJkiSVguFFkiRJUikYXiRJkiSVguFFkiRJUikYXiRJkiSVguFFkiRJUikcOegCdF/T09M0Go0FX9+9ezdTU1P3WdZsNqnVav0uTZIkSRooR17WmEajQbPZ7GqdWq1GvV7vU0WSJEnS2uDIyxpUq9UYGxtr+9revXsXfE2SJEkaZo68SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoF7/MiSVLFHdq5k8aWLR2/f129Tv2SS2Dbtj5WJYBms8nU1FTH76/X64yOjvaxImmwDC+SJFXUunod9uxh3aFDnLpnT+cr7tkDF1wAV17Zv+JEvV7v6v3NZhPA8KKhZniRJKmi6pdc0gohBw50vM6hnTtZd+hQV+toeUZHR7sKIt2M0EhlZXiRJKmqtm3revSksWVLd6M0ktRDXrAvSZIkqRQceemj6elpGo1GV+s0m01qtVqfKpIkSZLKy5GXPmo0GvdePNepWq3W9QV6kiRJUhU48tJntVqNsbGxQZchSZIklZ4jL5IkSZJKwfAiSZIkqRQML5IkSZJKwfAiSZIkqRQML5IkSZJKYcnwEhFvjYhGRHxl3rJNEXFVRHyz+Hl8f8uUJEmSVHWdjLy8DXjaYcvOAz6WmacDHyueS5IkSVLfLBleMvNTwP7DFp8NXFo8vhR4Vo/rkiRJkqT7WO5NKkcy8zaAzLwtIha8JXxEbAe2A4yMjDA5ObnMXcLs7OyK1l9tu3fvBmDv3r0922bZ2qDXqn78YBtU/fjBNpAkVddyw0vHMnMHsANgfHw8JyYmlr2tyclJVrL+apuamgJgbGysZ9ssWxv0WtWPH2yDqh8/2AaSpOpa7mxjeyPiJIDiZ6N3JUmSJEnS/S03vHwQOKd4fA7wgd6UI0mSJEntdTJV8juBzwJnRMStEfFi4PXAUyPim8BTi+eSJEmS1DdLXvOSmc9f4KUn97iWNW16eppGo7uz45rNJrVarU8VSZIkSdWy3NPGKqfRaNBsNrtap1arUa8vOBGbJEmSpC70fbaxYVKr1Xo6c5gkSZKkzjnyIkkaShFxSkR8IiJuiIivRsTLBl2TJGllHHmRJA2ru4FXZOZ1EbEBuDYirsrMrw26MEnS8hheJElDKTNvA24rHh+IiBuAkwHDSw/ctWsXM2eeeb/lpzebNNpMVvPAE07guAsvhG3bVqO8ymo2m/feJLtT9Xqd0dHRPlUk9ZbhRZI09CJiM/BIYGeb17YD2wFGRkaYnJxc9n5mZ2dXtH4ZnPqA1hnnR995J/Ubbrjf64tNU9M491y+eP75fapsbRjkZ2Dfvn3MzMx0tc7BgwdZv349W7du7VkdVfj/YDFVP37obxsYXiRJQy0iasA/AS/PzDsOfz0zdwA7AMbHx3NiYmLZ+5qcnGQl65fBN177Wmb//u85LqLt6/v372fTpk33WXbXrl2tsLNu3dC3T9k+A3OjNL2ckKhsbdBrVT9+6G8bGF4kSUMrIo6iFVwuz8z3DbqeYXDnwx/OnosvXvCX3V1tfmmZOfPMtqM0ktQtZxuTJA2liAjgLcANmfnXg65HkrRyhhdJ0rB6HPAC4EkRMVX8e/qgi5IkLZ+njUmShlJmfgZof2GGJKmUHHmRJEmSVAqGF0mSJEmlYHiRJEmSVAqGF0mSJEmlYHiRJEmSVAqGF0mSJEmlUMmpkqenp2k0Gl2t02w2qdVqfapIkiRJ0lIqOfLSaDRoNptdrVOr1ajX632qSJIkSdJSKjnyAq0wMjY2NugyJEmSJHWokiMvkiRJksrH8CJJkiSpFAwvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8CJJkiSpFAwvkiRJkkrB8CJJkiSpFI4cdAErNT09TaPR6GqdZrNJrVbrU0WSJEmS+qH0Iy+NRoNms9nVOrVajXq93qeKJEmSJPVD6UdeoBVGxsbGBl2GJEmSpD4q/ciLJEmSpGowvEiSJEkqBcOLJEmSpFIwvEiSJEkqBcOLJEmSpFIwvEiSJEkqBcOLJEmSpFIwvEiSJEkqBcOLJEmSpFIwvEiSJEkqBcOLJEmSpFIwvEiSJEkqBcOLJEmSpFIwvEiSJEkqBcOLJEmSpFJYUXiJiKdFxL9FxI0RcV6vipIkqRfspyRpuCw7vETEEcDfAD8PnAk8PyLO7FVhkiSthP2UJA2flYy8bANuzMxvZeZdwLuAs3tTliRJK2Y/JUlD5sgVrHsycMu857cCj1lZOQu78cYb2b17N1NTU/dZ3mw2qdVq/dqtJKm8VrWfqpJms3m//nhOu75606FDABzauZPGli19r2+QHnLwIDevXz/oMjp27F13kT/8Id98QO8ugx49dIhvrlvXs+2VTdWPH+Ah99wDv/qrcNFFPd/2SsJLtFmW93tTxHZgO8DIyAiTk5PL2tn09DTNZpNrrrnmfq9t3LiRvXv3Lmu7ZTM7O7vsNhwGVT9+sA2qfvxgG3RhVfspqMZ/m3379jEzM7Pg6+366icWvxivO3SIU/fs6Wt9ktaG/R//OLv68H24kvByK3DKvOcPBqYPf1Nm7gB2AIyPj+fExMSydzg5OclK1h8GVW+Dqh8/2AZVP36wDbpgPzUAbdtgfBwuuAAOHBhMUato//79bNq0adBlDFTV26Dqxw9FGzzpSX35PlxJeNkFnB4RDwX+HXge8F97UpUkSStnP7VWbNsGV1456CpWxS4DbOXboOrHD/1tg2WHl8y8OyJ+B5gEjgDemplf7VllkiStgP2UJA2flYy8kJn/AvxLj2qRJKmn7Kckabj0bmoJSZIkSeojw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSoFw4skSZKkUjC8SJIkSSqFyMzV21nE7cC3V7CJE4Dv9qicsqp6G1T9+ME2qPrxw8rb4CGZeWKvihkm9lM9UfU2qPrxg21Q9eOHPvZTqxpeVioivpCZ44OuY5Cq3gZVP36wDap+/GAbrGX+t7ENqn78YBtU/fihv23gaWOSJEmSSsHwIkmSJKkUyhZedgy6gDWg6m1Q9eMH26Dqxw+2wVrmfxvboOrHD7ZB1Y8f+tgGpbrmRZIkSVJ1lW3kRZIkSVJFlSK8RMTTIuLfIuLGiDhv0PUMQkTcFBHXR8RURHxh0PWshoh4a0Q0IuIr85ZtioirIuKbxc/jB1ljPy1w/K+NiH8vPgdTEfH0QdbYbxFxSkR8IiJuiIivRsTLiuWV+BwscvyV+hyURdX7Kvupe5dV4vtpTtX7Kvup1e+n1vxpYxFxBPAN4KnArcAu4PmZ+bWBFrbKIuImYDwzKzNveEQ8AWgCl2Xmw4tlbwD2Z+bri18Ojs/MVw2yzn5Z4PhfCzQz88JB1rZaIuIk4KTMvC4iNgDXAs8CXkgFPgeLHP9zqdDnoAzsq+ynqthPgX2V/dTq91NlGHnZBtyYmd/KzLuAdwFnD7gmrYLM/BSw/7DFZwOXFo8vpfU/yFBa4PgrJTNvy8zriscHgBuAk6nI52CR49faY19VQVXvp8C+yn5q9fupMoSXk4Fb5j2/lWp23gl8NCKujYjtgy5mgEYy8zZo/Q8D1AdczyD8TkR8uRiqH8ph6HYiYjPwSGAnFfwcHHb8UNHPwRpmX2U/Nady308LqNx3lP1UbGYV+qkyhJdos2xtn+vWH4/LzEcBPw/8j2KYVtXzt8BWYAy4Dfhfgy1ndUREDfgn4OWZeceg61ltbY6/kp+DNc6+yn5KP1K57yj7qdXrp8oQXm4FTpn3/MHA9IBqGZjMnC5+NoD30zpFoYr2FudXzp1n2RhwPasqM/dm5j2Z+UPg76nA5yAijqL1hXh5Zr6vWFyZz0G746/i56AEKt9X2U/dqzLfTwup2neU/dTq9lNlCC+7gNMj4qERcTTwPOCDA65pVUXEscVFUETEscDPAV9ZfK2h9UHgnOLxOcAHBljLqpv7Iiz8EkP+OYiIAN4C3JCZfz3vpUp8DhY6/qp9Dkqi0n2V/dR9VOL7aTFV+o6yn1r9fmrNzzYGUEyvdhFwBPDWzPzzAZe0qiJiC62/YgEcCbyjCm0QEe8EzgJOAPYCrwH+GbgCOBW4GXhOZg7lhYILHP9ZtIZgE7gJ+O9z59QOo4j4GeDTwPXAD4vFf0TrfNqh/xwscvzPp0Kfg7Kocl9lP1XNfgrsq+ynVr+fKkV4kSRJkqQynDYmSZIkSYYXSZIkSeVgeJEkSZJUCoYXSZIkSaVgeJEkSZJUCoYXDaWIaK7ivs6KiP+yWvuTJPVHRNwTEVMR8ZWIeE9EHLOMbbx8ofUWe62L7W+OiK7umRER16xkn13uyz5RfWV4kVbuLKCrL+qIOLI/pUiSVuBgZo5l5sOBu4DfWsY2Xg4sFFAWfC0ijljGvjqSmasZJs7CPlF9ZHjRUCv+AnR1RLw3Ir4eEZdHy89HxBWHve/K4vHPRcRnI+K64i9vtWL5TRFxfrH8+oj4TxGxmVbn9rvFX+seHxEPiYiPRcSXi5+nFuu/LSL+OiI+AfxVRHwzIk4sXntARNwYESeschNJktr7NHAaQET8XjEa85WIeHmx7NiI+HBEfKlY/qsRcS4wCnyi+K6/V7vXIqIZEa+LiJ3AT0fEqyNiV7G9HcXdy4mIRxf7+SzwP+Zt84iI+KtinS9HxH9vdyBzZyPYJ2oYGF5UBY+k9deuM4EtwOOAq4DHRsSxxXt+FXh38UX5J8BTMvNRwBeA35u3re8Wy/8WeGVm3gT8HfDG4q91nwbeBFyWmT8JXA5cMm/9hxXb/l3g7cCvFcufAnwpM7/b20OXJHWrGAn4eeD6iHg08CLgMcBjgZdExCOBpwHTmfmIYqTmI5l5CTANPDEznzh/mwu8dizwlcx8TGZ+BnhTZv5Usb31wDOL9/0DcG5m/vRhpb4YmM3MnwJ+qqjtoUscnn2iSs3woir4fGbempk/BKaAzZl5N/AR4BeKTuoZwAdodUxnAv9fREwB5wAPmbet9xU/rwU2L7C/nwbeUTz+R+Bn5r32nsy8p3j8VuA3ise/SatzkiQNzvriu/8LwM3AW2h9h78/M7+fmU1a/cDjgeuBp0TEX0bE4zNzdhn7uwf4p3nPnxgROyPieuBJwH+OiOOAjZn5yeI9/zjv/T8H/EZR807gx4DTl9infaJKzXMMVQU/mPf4Hn70uX83reH3/cCuzDxQDNFflZnPX2Jb87ezlJz3+Pv3Lsy8JSL2RsSTaP1F79fut6YkaTUdzMyx+QvmTt06XGZ+oxiVeTrwFxHx0cx8XZf7OzT3y3tErAP+NzBe9A+vBdYBwX37kfuUB7w0Mye72Kd9okrNkRdV2dXAo4CX0PrSBvgc8LiImDvP+ZiIeNgS2zkAbJj3/BrgecXjXwM+s8i6b6Y1VH7FvL8+SZLWjk8Bzyr6g2OBXwI+HRGjwJ2Z+XbgQlr9Cdy/T5hvsdfWFT+/W1xX8myAzJwBZiNibsRi/i/1k8BvR8RRABHxsHmnfnXrauwTVQKGF1VW8cX4IVrnNX+oWHY78ELgnRHxZVpf3P9piU1dCfzS3MWJwLnAi4r1XwC8bJF1PwjUcHhcktakzLwOeBvweVqnZr05M78I/F/A54vTqf4Y+LNilR3Avx5+wf5SrxUh5e9pnY72z8CueS+/CPib4oL9g/OWvxn4GnBdtKZP/n9Z5lk19okqi8hcaCRSUr9FxDitCxsfP+haJEkaJPtEdcJrXqQBiYjzgN/G83olSRVnn6hOOfIiSZIkqRS85kWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJWC4UWSJElSKRheJEmSJJXC/w8Alaw69xN0jQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1008x576 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Diff in value functions = -8.159e-06\n",
      "Diff in policy functions = 0.00000\n"
     ]
    }
   ],
   "source": [
    "mod = inventory_model(max_inventory=25)\n",
    "mod.dp=.25\n",
    "mod.c = .25\n",
    "mod.p = 3.5\n",
    "mod.r = 0.4\n",
    "mod.β = 0.9\n",
    "mod.demand_pr(plot=True)\n",
    "ev1,q1 = mod.solve_vfi()\n",
    "mod.solve_show()\n",
    "ev2,q2 = mod.solve_policyiter()\n",
    "mod.solve_show(method='policyiter');\n",
    "print('Diff in value functions = %1.3e\\nDiff in policy functions = %1.5f'%(np.amax(ev1-ev2),np.amax(q1-q2)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal orders of new inventory for y = max(x-d,0):\n",
      "[7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n"
     ]
    }
   ],
   "source": [
    "def optimal_policy(self,ev):\n",
    "    '''Computes the optimal policy function as function of post trade stock for the stochastic\n",
    "    inventory dynamics model for given EV function'''\n",
    "    # idea: 2-dim array with q in axes 0, y = max(x-d,0) in axis 1\n",
    "    q = self.x[:,np.newaxis]  # choices\n",
    "    y = self.x[np.newaxis,:]  # post trading stock\n",
    "    # indexes for next period value with extrapolation using last value\n",
    "    i = np.minimum(y+q,self.upper)\n",
    "    # compute the Bellman maximand\n",
    "    vm = -self.r*q -self.c*(q>0) + self.β*ev[i]\n",
    "    # find argmax and argmax\n",
    "    return np.argmax(vm,axis=0)  # maximum in every column\n",
    "\n",
    "ev,_ = mod.solve_vfi()\n",
    "q = optimal_policy(mod,ev)\n",
    "print('Optimal orders of new inventory for y = max(x-d,0):')\n",
    "print(q)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "hide-output": false,
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Optimal orders of new inventory for d,x:\n",
      "(d in rows, x in columns)\n",
      "[[7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2 0]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3 2]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4 3]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5 4]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6 5]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 6]\n",
      " [7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7 7]]\n"
     ]
    }
   ],
   "source": [
    "# COPY from exercise 11\n",
    "def optimal_policy_old(m,ev):\n",
    "    '''Computes the optimal policy function for the stochastic\n",
    "    inventory dynamics model for given EV function'''\n",
    "    # idea: 3-dim array with q in axes 0, d in axis 1 and x in axis 2\n",
    "    q = m.x[:,np.newaxis,np.newaxis]  # choices\n",
    "    d = m.x[np.newaxis,:,np.newaxis]  # demand\n",
    "    x = m.x[np.newaxis,np.newaxis,:]  # inventories\n",
    "    # compute current period profit (relying on numpy broadcasting to get the matrix with choices in rows)\n",
    "    p = m.profit(x,d,q)  # 3-dim array\n",
    "    # indexes for next period value with extrapolation using last value\n",
    "    i = np.minimum(m.next_x(x,d,q),m.upper)\n",
    "    # compute the Bellman maximand\n",
    "    vm = p + m.β*ev[i]\n",
    "    # find argmax and argmax\n",
    "    return np.argmax(vm,axis=0)  # maximum in every column\n",
    "\n",
    "q = optimal_policy_old(mod,ev)\n",
    "print('Optimal orders of new inventory for d,x:\\n(d in rows, x in columns)')\n",
    "print(q)\n",
    "# Note the symmetry in the optimal policy!\n",
    "# This implies that knowing both x and d is not necessary for the\n",
    "# optional new order, it's enough to condition on the inventory\n",
    "# remaining after sales, i.e. x-min(x,d) = max(0,x-d)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Convergence rate comparison between policy iterations and VFI\n",
    "\n",
    "- generally, policy iterations converge much faster  \n",
    "- more work is needed for each iteration (policy evaluation step!)  \n",
    "- if first step done efficiently, significant advantage over VFI  "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "#### Modified policy iteration algorithm\n",
    "\n",
    "In problems where the first step (policy evaluation) is hard to perform efficiently, it can be\n",
    "replaced by a finite forward simulation of Bellman operator under current policy to yield an approximate\n",
    "value."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "slideshow": {
     "slide_type": "slide"
    }
   },
   "source": [
    "### Further learning resources\n",
    "\n",
    "- Discrete state dynamic programming and policy iterations at QuantEcon\n",
    "  [https://python-advanced.quantecon.org/discrete_dp.html](https://python-advanced.quantecon.org/discrete_dp.html)  \n",
    "- On theoretical properties of Howard’s algorithm\n",
    "  [https://link.springer.com/chapter/10.1007/978-3-642-17517-6_37#:~:text=Howard’s%20policy%20iteration%20algorithm%20is,Markov%20Decision%20Processes%20(MDPs).&text=The%20theoretical%20complexity%20of%20Howard’s,known%20on%20its%20running%20time](https://link.springer.com/chapter/10.1007/978-3-642-17517-6_37#:~:text=Howard's%20policy%20iteration%20algorithm%20is,Markov%20Decision%20Processes%20%28MDPs%29.&text=The%20theoretical%20complexity%20of%20Howard's,known%20on%20its%20running%20time).  "
   ]
  }
 ],
 "metadata": {
  "celltoolbar": "Slideshow",
  "date": 1612589586.711941,
  "download_nb": false,
  "filename": "43_policy_iter.rst",
  "filename_with_path": "43_policy_iter",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "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"
  },
  "title": "Foundations of Computational Economics #43"
 },
 "nbformat": 4,
 "nbformat_minor": 4
}