courses/IMTx-Queue-Theory/Week3_Lab_Discrete_Time_Markov_Chains (1).ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<p><font size=\"4\">  MOOC: Understanding queues</font></p>\n",
    "<p><font size=\"4\">  Python simulations</p>\n",
    "<p><font size=\"4\">  Week III: Discrete time Markov chains </p>\n",
    "\n",
    "----------------------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this lab, we consider the Markov chain of the weather forecast example of the course. We check convergence of the probability $\\pi(t)$ of the chain at time $t$ to a steady state distribution $\\pi^*$, independently from the initial distribution $\\pi(0)$ of the chain. We solve the load balance equations to get $\\pi^*$.\n",
    "\n",
    "-------------"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us consider the Markov chain of the weather forecast example of the course. Recall that its states 1, 2 and 3 represent clear, cloudy and rainy states, and the transition matrix is\n",
    "\n",
    "$$\n",
    "P=\n",
    "\\begin{pmatrix}\n",
    "0.7 & 0.3 & 0\\\\\n",
    "0.3 & 0.5 & 0.2\\\\\n",
    "0.1 & 0.4 & 0.5\n",
    "\\end{pmatrix}.\n",
    "$$\n",
    "\n",
    "**1)** Complete below the code of the function that generates trajectories of the Markov chain. The function inputs are the chain initial state $x0$, the transition matrix $P$ and final time index $T$. Its output will be a  trajectory $x$ of the chain observed between instants $0$ and $T$. Draw a trajectory of the evolution of the weather between time 0 and time $T=100$. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline              \n",
    "from pylab import *          "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "P    = array([[.7, .3, 0], [.3, .5, .2], [.1, .4, .5]])\n",
    "\n",
    "def X(x0,P=P,T=100):\n",
    "    # Function X supplies a trajectory of the discrete Markov chain \n",
    "    # with initial state x0 and transition matrix P, till time T\n",
    "    x        = [x0]\n",
    "    for _ in range(T):\n",
    "        #####################\n",
    "        # supply the vector p of probabilities to transit to states\n",
    "        # 1,2,3 from the last calculated state\n",
    "        p = P[ x[len(x)-1]-1 ]\n",
    "        #####################\n",
    "        u = rand()\n",
    "        if u<p[0]:\n",
    "            x.append(1)\n",
    "        elif u<p[0]+p[1]:\n",
    "            x.append(2)\n",
    "        else:\n",
    "            x.append(3)\n",
    "    return array(x)\n",
    "\n",
    "V1 = mean(X(x0=1,T=10**4)) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA4AAAAEWCAYAAAA3lPkgAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAFpJJREFUeJzt3X20pVddH/DvjySAEiBAYqC8DViiYqwDDO8UB5EuI7TIqthQEJOlxioIvi2XdHVV7KpdamlRQbARjNEVeTGAhSxaQcgQSAGZhGkCAUReYnAiGQvhnZiQX/84zzR3Te/cO5lzzpw7d38+a91177PP3vvZ59x99rnf+zzPOdXdAQAAYPu7w6oHAAAAwLEhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAWJKqOqeq3rPqcQDAQQIgANtaVb2oqt56SNnHD1N29hz72VFVXVUnHm0fALBsAiAA291lSR5fVSckSVXdO8lJSR5+SNk/nupuSYIlAIsgAAKw3X0gs8C3c9p+YpJLk3zskLJPdPf+qvr2qnp7VX2uqj5WVT98sKOqempVfbCqvlhV11XVi9fs52B4vLGqvlxVj13T7iVV9fmq+lRVnbWm/O5V9eqqur6q/raq/uOaUHpOVV1eVS+tqs8lWbsvADgqAiAA21p3/0OS92cW8jJ9f3eS9xxSdllV3SXJ25P8SZJvSfKsJK+oqu+c6n0lyXOTnJLkqUl+qqp+cE0fSXJKd5/c3e+dth+dWdg8NclvJnl1VdV024VJbsns6OPDkvyzJD++ZviPTvLJaSy/NsfDAABJBEAAxvCu3BbQ/mlmAfDdh5S9K8nTkny6uy/o7lu6+8okb0jyQ0nS3Xu6++ruvrW7r0rymiTfs8m+r+3u3+/ub2QW+O6T5PSqOj3JWUl+tru/0t03JHlpkrXXIe7v7pdNY/nafA8BACSuJwBgBJcleV5V3SPJad398ar6bJILp7IzpzrPTPLoqrpxTdsTk/xxklTVo5P8+lT/jknulORPN9n33x38obu/Oh38OznJPTM7NfX62w4I5g5JrlvTdu3PADA3ARCAEbw3yd2TnJfk8iTp7i9W1f6pbH93f6qqrkvyru5+ymH6+ZMkL09yVnd/vap+K7NTO5Okb+eYrktyU5JTu/uWw9S5vX0CwIacAgrAtjedPrk3yc9ndurnQe+Zyg6+gcslSc6oqh+pqpOmr0dW1XdMt981yeem8PeoJP96TV8Hktya5MFHOKbrk7wtyX+pqrtV1R2q6lurarNTSgHgqAmAAIziXZm9mcraD2Z/91R2WZJ095cyeyOWs5Psz+z0zd/I7FTPJPnpJP+hqr6U5N8nef3Bjrr7q5m9UcvlVXVjVT3mCMb03MxOJb0myeeTXJzZNYIAsBTV7ewSAACAETgCCAAAMAgBEAAAYBACIAAAwCAEQAAAgEFsi88BPPXUU3vHjh2rHgYAAMBKXHHFFX/f3adtVm9bBMAdO3Zk7969qx4GAADASlTVtUdSzymgAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCCOOgBW1Teqal9Vfaiq3lJVpxxBm/91tPsDAABgPvMcAfxad+/s7jOTfC7J8zZr0N2Pm2N/AAAAzGFRp4C+N8l9k6SqTq6qd1TVlVV1dVU9/WClqvry9H13Ve2pqour6qNVdVHNPLmq3rSm/lOq6o0LGiMAAMDQTpy3g6o6IcmTk7x6Kvp6kmd09xer6tQk76uqN3d3H9L0YUm+M8n+JJcneXySdyb53ao6rbsPJDk3yQXzjpGtY/fuo2+7Z898fW/Ufp62mzle+16lee/XRu2P18dkM8t8bsEijDpHl7kejbjWbVer/FthWfud1zzj3szoz495jgB+U1XtS/J/ktwzydun8kryn6rqqiR/kdmRwdPXaf+X3f2Z7r41yb4kO6aQ+MdJnjNdU/jYJP9jvZ1X1XlVtbeq9h44cGCOu8Gxsnt3sm/f0bXdt2/jhWCzvjdqP0/bzRyvfa/SvPdro/bH62OymWU+t2ARRp2jy1yPRlzrtqtV/q2wkVXOo3nGvRnPj/mOAH6tu3dW1d2TXJLZNYC/k+TZSU5L8ojuvrmqPp3kzuu0v2nNz99YM5YLkrwlsyOJf9rdt6y38+4+P8n5SbJr165Djy6yRe3ceXT/dTmSJ+pGfW/Wfp62mzle+16lee/X4dofz4/JZpb53IJFGHWOLnM9GnGt265W9bfCRlY9j4523JtZ9f3aCua+BrC7v5DkBUl+sapOSnL3JDdM4e9JSR54O/vbn9lpof8uyR/OOz4AAABm5r4GMEm6+4NV9b+TnJ3koiRvqaq9mZ3a+dGj6PKiJKd19zWLGB8AAABzBMDuPvmQ7X++ZvOxG7Xp7j1J9qwpf/4hVZ+Q5PePdmwAAAD8/xZyBHCRquqKJF9J8gurHgsAAMB2suUCYHc/YtVjAAAA2I4W9UHwAAAAbHECIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAM4ogCYFXdu6peW1WfqKprquqtVXVGVX1oEYOoqnOq6uWL6AsAAID1nbhZhaqqJG9KcmF3nz2V7Uxy+pLHBgAAwAJtGgCTPCnJzd39ewcLuntfVe04uF1Vd07yyiS7ktyS5Oe7+9KqOifJru5+/lTvkiQv6e49VXVukhcluT7JXyW5qarumuSqJGd0981Vdbdp+yHdffPc93YL2b17eX3v2bO8fc/T9759yc6dR7/vffsO3/+R9H249vO03czx2ncy/zzaqP0y58m8lvncXJZlPreSzefCRlb5eK5qDi9i31vVPHNhHvOsZcly5/A8fS/zuTfvY7aRUef3su7XvK/ny3xernIeLdMq15St4EhOAT0zyRWb1HleknT3dyV5VpILp1C4rqq6T5JfTfL4JE9J8tCp/ZeS7Eny1Knq2UnesF74q6rzqmpvVe09cODAEdyNrWP37tnEW4bNJvQ8+5637507j/4Js2fPxovjZn1v1H6etps5Xvue93e9UftlzpN5LfO5uUzLfG7N8yK5ysdzlXN4mevwKi3zD8GNzLOWJcudw/P+k21Zz715H7ONjDq/l3m/5nk9X3ZAW9U8WqZVrilbxZEcATwST0jysiTp7o9W1bVJztig/qOT7OnuA0lSVa9bU/9VSX4pyZ8lOTfJT6zXQXefn+T8JNm1a1cv4D4cU8v6I/dIJuTR7nuZfR+Jefudp/0yA8lW7Hve3/Vm7VcZ8jazlce2LPMc6drMqh7PVc7hVa+Vy7KVj+ZuZJlzeN6+l/ncW9b8Gnl+r/J+LWsOHu1+57XqAHW8ntmyKEdyBPDDSR6xSZ06TPkth+xj7VHBdUNbd1+eZEdVfU+SE7p7IW80AwAAMLojCYDvTHKnqvp/R+Kq6pFJHrimzmVJnj3ddkaSByT5WJJPJ9lZVXeoqvsnedRU//1JdlfVvarqpCTPPGSff5TkNUkuuN33CAAAgHVtGgC7u5M8I8lTpo+B+HCSFyfZv6baK5KcUFVXJ3ldknO6+6Yklyf5VJKrk7wkyZVTn9dPfbw3yV8cLF/joiT3yCwEAgAAsABHdA1gd+9P8sPr3HTmdPvXk5yzTrvOdGRwndsuyOGP8D0hycXdfeORjA8AAIDNLepNYBamql6W5KwkP7DqsQAAAGwnWy4AdvfPrHoMAAAA29GRvAkMAAAA24AACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBLCQAVtW9q+q1VfWJqrqmqt5aVWdU1YcW0T8AAADzO3HeDqqqkrwpyYXdffZUtjPJ6fP2vab/6u5bF9EfAADAqBZxBPBJSW7u7t87WNDd+5Jcd3C7qk6oqv9cVR+oqquq6ien8pOr6h1VdWVVXV1VT5/Kd1TVR6rqFUmuTHL/BYwTAABgaHMfAUxyZpIrNqnzY0m+0N2PrKo7Jbm8qt6WWUh8Rnd/sapOTfK+qnrz1Obbkpzb3T+9gDGuxO7d65fv25fs3Lm8/e7bt7x9L7NvtpZ5f9eHa7+IeTLP2Mzh22+jx2yzdqt8PFc5h7frPNuOc2Gzdqv8XW/VebRVxzWv4/V+GTe3V3X3fB1UvSDJg7r75w4p35Hkku4+s6ouTvJPknx1uvnuSX4yyaVJXprkiUluzSz0PSjJnZNc2t0P2mC/5yU5L0ke8IAHPOLaa6+d634sw0YvNnv2rGa/8+57mX2ztcz7u17m/J9nbObw7Xc0fzgftMrHc5VzeLvOs+06Fzayyt/1Vp1HW3Vc8zpe75dxc1BVXdHduzatt4AA+OQkv9LdTzykfEduC4BvSHJ+d//5IXXOSXJWkud0981V9ekku6ebL+nuM49kDLt27eq9e/fOczcAAACOW0caABdxDeA7k9ypqn5izc4fmeSBa+r8eZKfqqqTptvPqKq7ZHYk8IYp/D3pkDYAAAAs0NwBsGeHEJ+R5CnTx0B8OMmLk+xfU+1VSa5JcuX00RD/LbPrDy9Ksquq9iZ5dpKPzjseAAAA1jf3KaBbgVNAAQCAkR3LU0ABAAA4DgiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBACIAAAwCAEQAAAgEEIgAAAAIMQAAEAAAYhAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGIQACAAAMQgAEAAAYhAAIAAAwCAEQAABgEAIgAADAIARAAACAQQiAAAAAgxAAAQAABiEAAgAADEIABAAAGIQACAAAMAgBEAAAYBDV3asew9yq6kCSa1c9jnWcmuTvVz0ItjVzjGPBPONYMM9YNnOMY2GV8+yB3X3aZpW2RQDcqqpqb3fvWvU42L7MMY4F84xjwTxj2cwxjoXjYZ45BRQAAGAQAiAAAMAgBMDlOn/VA2DbM8c4FswzjgXzjGUzxzgWtvw8cw0gAADAIBwBBAAAGIQACAAAMAgBcAmq6vur6mNV9ddV9curHg/bQ1Xdv6ouraqPVNWHq+qFU/k9q+rtVfXx6fs9Vj1Wjm9VdUJVfbCqLpm2H1RV75/m2Ouq6o6rHiPHt6o6paourqqPTmvaY61lLFpV/dz0evmhqnpNVd3Zesa8quoPquqGqvrQmrJ116+a+Z0pE1xVVQ9f3chvIwAuWFWdkOR3k5yV5KFJnlVVD13tqNgmbknyC939HUkek+R509z65STv6O6HJHnHtA3zeGGSj6zZ/o0kL53m2OeT/NhKRsV28ttJ/md3f3uS785svlnLWJiqum+SFyTZ1d1nJjkhydmxnjG/P0zy/YeUHW79OivJQ6av85K88hiNcUMC4OI9Kslfd/cnu/sfkrw2ydNXPCa2ge6+vruvnH7+UmZ/MN03s/l14VTtwiQ/uJoRsh1U1f2SPDXJq6btSvK9SS6eqphjzKWq7pbkiUlenSTd/Q/dfWOsZSzeiUm+qapOTPLNSa6P9Yw5dfdlST53SPHh1q+nJ/mjnnlfklOq6j7HZqSHJwAu3n2TXLdm+zNTGSxMVe1I8rAk709yendfn8xCYpJvWd3I2AZ+K8kvJbl12r5Xkhu7+5Zp25rGvB6c5ECSC6ZTjV9VVXeJtYwF6u6/TfKSJH+TWfD7QpIrYj1jOQ63fm3JXCAALl6tU+azNliYqjo5yRuS/Gx3f3HV42H7qKqnJbmhu69YW7xOVWsa8zgxycOTvLK7H5bkK3G6Jws2XYP19CQPSvKPktwls9PxDmU9Y5m25GuoALh4n0ly/zXb90uyf0VjYZupqpMyC38Xdfcbp+LPHjydYPp+w6rGx3Hv8Un+RVV9OrPT1783syOCp0ynUCXWNOb3mSSf6e73T9sXZxYIrWUs0vcl+VR3H+jum5O8McnjYj1jOQ63fm3JXCAALt4HkjxkepepO2Z2wfGbVzwmtoHpWqxXJ/lId//XNTe9OcmPTj//aJL/fqzHxvbQ3S/q7vt1947M1q53dvezk1ya5IemauYYc+nuv0tyXVV921T05CTXxFrGYv1NksdU1TdPr58H55n1jGU43Pr15iTPnd4N9DFJvnDwVNFVqu6VH4XcdqrqBzL7r/kJSf6gu39txUNiG6iqJyR5d5Krc9v1Wf82s+sAX5/kAZm94D2zuw+9OBlul6raneQXu/tpVfXgzI4I3jPJB5M8p7tvWuX4OL5V1c7M3mjojkk+meTczP4pbS1jYarqV5P8q8zeRfuDSX48s+uvrGcctap6TZLdSU5N8tkkv5Lkz7LO+jX98+Hlmb1r6FeTnNvde1cx7rUEQAAAgEE4BRQAAGAQAiAAAMAgBEAAAIBBCIAAAACDEAABAAAGceLmVQBgLFV1ryTvmDbvneQbSQ5M21/t7setZGAAMCcfAwEAG6iqFyf5cne/ZNVjAYB5OQUUAG6Hqvry9H13Vb2rql5fVX9VVb9eVc+uqr+sqqur6luneqdV1Ruq6gPT1+NXew8AGJkACABH77uTvDDJdyX5kSRndPejkrwqyc9MdX47yUu7+5FJ/uV0GwCshGsAAeDofaC7r0+SqvpEkrdN5VcnedL08/cleWhVHWxzt6q6a3d/6ZiOFAAiAALAPG5a8/Ota7ZvzW2vsXdI8tju/tqxHBgArMcpoACwXG9L8vyDG1W1c4VjAWBwAiAALNcLkuyqqquq6pok/2bVAwJgXD4GAgAAYBCOAAIAAAxCAAQAABiEAAgAADAIARAAAGAQAiAAAMAgBEAAAIBBCIAAAACD+L/yqVdmiaAgZQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa836de15f8>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def step(x,y,Tmax=0,color='b'):\n",
    "    # step function\n",
    "    # plots a step function representing the number\n",
    "    # of clients in the system at each instant\n",
    "    if Tmax==0:\n",
    "        Tmax = max(x)\n",
    "    x = append(x,[Tmax])  # number of clients\n",
    "    y = append(y,[y[-1]]) # instants of events\n",
    "    for k in range(len(x)-1):\n",
    "        vlines(x[k+1],y[k],y[k+1],color=color)\n",
    "        hlines(y[k],x[k],x[k+1],color=color)\n",
    "        \n",
    "T = 100\n",
    "x = X(x0=1)\n",
    "figure(num=None, figsize=(15, 4))\n",
    "step(range(T),x)\n",
    "axis(ymin=0.5,ymax=3.5)\n",
    "xlabel(\"Time\")\n",
    "title(\"Weather\")\n",
    "yticks([1.0,2.0,3.0], [\"Clear\",\"Cloudy\",\"Rainy\"]);"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**2)** Run the following code that computes recursively the state probability vectors $\\pi(t)$ at times $t=0,\\ldots,100$. The state probability vectors can be computed recursively : $\\pi(t+1)=\\pi(t) P$. Check that, when changing the initial state $x0$, $\\pi(t)$ still converges rapidly to the same asymptotic vector $\\pi^*$ as $t$ increases."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmcAAAEYCAYAAAAZLlGYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAIABJREFUeJzs3Xd4XOWZ9/HvPSPNaFRmJFuyLfeCDS4Ug+lhQ0mhkwIJLQRSyJvdkM2mbTabtiS7kJBs2qbAJgQ2pFASQg8htCQQbAzYgI1t3ItcZMnqXfO8f5yRPJIlWbKkOcea3+fKXDNz5syZG8eP5577aeacQ0RERESCIeR3ACIiIiKyn5IzERERkQBRciYiIiISIErORERERAJEyZmIiIhIgCg5ExEREQkQJWciIiIiAaLkTEQOiZltNrO3HeJ7V5nZmSMc0mA+90gze8XM6s3sk4fw/pvM7FOH+NnLzGzhobxXRLKLkjORMS6VRDWbWUPa7X8y/Pk9kjjn3ELn3DOZiiHN54FnnHNFzrkf9H7RzErMzKX+jJrMrKIrGTOzMuAa4Na08+eYWaOZlacduyr1vmm9Lv9t4MaR+g8xs0+Y2XIzazWzOwb5nqiZ/dzMtqQS1FfM7LxBvO8uM9tpZnVmts7MPjLs/wAR6ZeSM5HscJFzrjDt9gm/A/LJDGDVAK8fB1Sm/ozygY8D3zWzqcC1wKPOueauk51zG4CHga4E7lTgf4B3Oee29br2g8BZ6YncMFUA3wBuH8J7coBtwFuBBPBl4B4zm3mQ990EzHTOxYGLgW+Y2QlDDVhEBkfJmUiWMrMvmNl9vY5938x+kHo838yeMbOaVDfkxQNcy5nZEWnP7zCzb5jZL4HpwEOpatTnU6/3qKYN9Fmpcz9rZq+aWa2Z3W1mef3EMdB1ngLOAv4nFcu8Pi5xHPBi2vOlqfsIcB7wbB/v+SbwMTNbBPwe+H/OuWW9T3LOtQAvAe/oK/ahcs793jn3B6BqCO9pdM59zTm32TmXdM49DGwCBky0nHOrnHOtXU9TtzmHGruIDEzJmUj2+g1wvpnFAcwsDLwP+LWZ5QIPAX8CJgA3AL8ysyOH8gHOuQ8AW9lfuftW73MG+VnvA84FZgHH4FWxhnQd59zZwF+BT6RiWddHyIuBZanrFQP/iZdQbQKOBtb28d/4cuo9S4GfOOfuHuCP5A3g2D5ifziVUPZ1e3iA6w2LmU0E5jFwNbHr3B+bWROwBtgJPDpacYlkOyVnItnhD72+8D/qnNsCvAy8K3XO2UCTc+4F4BSgELjZOdfmnHsKr/vuilGIbTCf9QPnXIVzrhovATvuEK9zMMcBnzOzarykzOEllg4oBup7v8HMQkAnkMSroqW/dpaZTU87VJ+6Tg/OuQudc8X93C4cQvyDlkpmfwXc6Zxbc7DznXP/CBQBZ+BVCFsHfoeIHColZyLZ4V29vvD/N3X81+xPXq5MPQeYDGxzziXTrrEFmDIKsQ3ms3alPW7CS8IO5Tr9MrMoMB84yjk3zjk3xzn3EefcztQp+/CSk96+g5dwvQlc1eu1D+EleF2KgJrBxDOaUgnlL4E2YNDjD51znc65vwFT8cbjicgoUHImkt3uBc5MDXh/N/uTswpgWupLvMt0YEc/12kC8tOeT0p77BjYUD9rtK6zCGh0zm3v5/VX8boAu5nZx/D+3N6FVzX7nJlZ6rWLgQuBX5rZ1am3zAdW9r6wmT3WazZt+u2xQcY/KKn4fg5MBN7rnGs/hMvkoDFnIqNGyZlIFnPOVQLPAL8ANjnn3ki9tBRoBD5vZrnmrUl2EfDbfi61ArjSzMJmdi7ebMAuu4HZA4Qx1M8aressZuCxV4+S9t+VmtDwX3jdnruB+/AmDlySOuVh4CXn3JnOubtSlbkTgCd6X9g5d16v2bTptz6XujCznNTEiDAQNrM8M8tJvXbHAMtr/AQvSbwofeZpf+8zswlmdrmZFab+/30nXrX1qQH+rERkGJSciWSHh3pVY+5Pe+3XwNvYXzXDOdeGt2TCecBe4MfANQOMTfpnvESoBq9r7w9pr90EfCk11u2zvd94CJ/VpxG4znHA6wO8/n94EyhiZnYUXtL3Aefca6nP7wT+G/jX1PlH4HV1drkYb421ikHGczBfApqBLwBXpx5/KfXaNOC53m8wsxnAx/D+W3el/X24aoD3ObwuzO14XbvfBj7lnHsg7bqPmdkXR+i/SyTrmTfOVUREDsbM/gvY45z73iDOfTcwo+tcM1sKfNg5N1ACOBIxRvC6To8ZSpflob5PREaekjMRkVFgZguAe4A/O+cOacsnEclOSs5EREREAiQjY87M7HYz22NmfZbzzfMDM1ufWgX8+EzEJeIXtQmRntQmRPbL1ISAO/BW9+7PecDc1O16vNlEImPZHahNiKS7A7UJESBDyZlz7i9A9QCnXAL8n/O8ABTbyG0OLBI4ahMiPalNiOyX43cAKVOAbWnPt6eO7ex9opldj/eriYKCghOOOuqoAy7W3NbJ+soGZozLJx7LHZ2IRXp56aWX9jrnykbociPaJmith6r1UDoXIn0tri8y8oLcJura6thWv405iTnk5eSNUIgiAxtsmwhKcmZ9HOtzpoJz7jbgNoAlS5a45cuXH3DOnvoWTvrPJ/n3ixfywdNmjmScIv0ysy0jebk+jh1ym6DiFbjtTLj8v+Go80cwTJH+BblNLN25lI/86SPc/s7bOXHSiSMYpkj/BtsmgrII7Xa8xQ+7TMXbiuWQlBZEyQ0bFbXNBz9ZJJhGtE2Ql9pru3nfcGIS8dOItol4JA5AXWvd8KISGQVBSc4eBK5JzcY5BahN22x4yEIhozwRY2dNy8hFKJJZI9omiJV49y2+77ktcqhGtE0kognA694UCZqMdGua2W+AM4FSM9sOfBXIBXDO/RRvz7rzgfV4GyhfN9zPLE/kUVGjypkEU8bbRDQOmCpnEliZbhNdlbPa1trhXEZkVGQkOXPOXXGQ1x3wTyP5mZOLYyzbNNDEHxH/ZLxNhEIQK4ZmVc4kmDLdJgpyCwhbWJUzCaSgdGuOuPJEHrvqWuhMagcEEcAbd6bKmQgAZkY8ElflTAJpzCZnk4tjdCYdlfWtfociEgyxEo05E0mTiCZUOZNAGsPJmbduzQ6NOxPxxEpUORNJE4+qcibBNGaTs/JEDICdWk5DxBNTt6ZIungkTm2bkjMJnjGbnE0uTiVnWk5DxBMr0YQAkTSJaELrnEkgjdnkLJ6XQ0EkrG5NkS55xd6Ys2TS70hEAkGVMwmqMZucmRnlxTF1a4p0iZWAS0Jbvd+RiARCIpqgoa2BzmSn36GI9DBmkzPwltPYWatuTRHAG3MGGncmkhKPxHE4Gtob/A5FpIcxnZxNKY5RoTFnIp6uLZw07kwESNvCSePOJGDGdHJWnoixt6GV1g6VrEX2J2eqnIkAJCJecqZxZxI0Yzs5S611tktdmyLehABQciaSEo9qf00JpjGdnE1JLaehrk0R9lfOtEuACLC/cqZdAiRoxnRyVp7wKmeasSmCJgSI9KLKmQTVGE/OuipnSs5EyI1BTp4mBIikxCNecqbKmQTNmE7OYpEwJfm5VGjMmYgnT1s4iXSJhCPEcmKqnEngjOnkDLxtnHaqcibiiZVozJlImngkrsqZBM6YT87KE1rrTKSb9tcU6SERTahyJoEz5pOzycV5VGhCgIgnpm5NkXSqnEkQZUFyFqO+pYOG1g6/QxHxnypnIj2ociZBNOaTs+7lNDTuTEQTAkR6iUfi2r5JAmfMJ2eTUwvR7lByJuJVztoboaPN70hEAiERTahbUwJnzCdn+xei1aQAke6FaDVjUwTwKmctnS20drb6HYpItzGfnE2M5xEydWuKAGmbnys5EwGvcgaoa1MCZcwnZ7nhEBOK8tih5TREtPm5SC9duwRoUoAEyZhPzgDKi/O0v6YIaPNzkV669tfUuDMJkqxIziYnYhpzJgLa/Fykl65uTVXOJEiyIzkrzqOiphnnnN+hiPire8yZkjMRSOvWbFNyJsGRFclZeSJGa0eS6kYtHyBZLs+rEmhCgIhHEwIkiLIiOZtcrOU0RAAIhSGaUOVMJKUwtxDDVDmTQMmS5MxbiLZCy2mIeOPONCFABICQhYhHtUuABEtWJGflCS85U+VMBG1+LtJLPBJX5UwCJSuSs/EFESLhkCpnIqDNz0V6SUS0hZMES1YkZ6GQMSmRR4UqZyKp5EyVM5EuiWhC3ZoSKFmRnIE3KUBbOIng7RKg5EykWzwSV+VMAiV7krNETN2aIuBVzlpqQOv+iQDeLgFahFaCJGPJmZmda2ZrzWy9mX2hj9enm9nTZvaKmb1qZueP5OeXF+exu76VzqS+kCQYfGsTsWJIdkBbw4hcTmSk+NUmuipnSZccicuJDFtGkjMzCwM/As4DFgBXmNmCXqd9CbjHObcYuBz48UjGMLk4RmfSsade487Ef762ie5dAjQpQILDzzaRiCZIuiSN7Y0jcTmRYctU5ewkYL1zbqNzrg34LXBJr3McEE89TgAVIxnA5ETXWmdKziQQ/GsTedpfUwLJtzbRvYWTujYlIDKVnE0BtqU93546lu5rwNVmth14FLihrwuZ2fVmttzMlldWVg46gPLULgEadyYB4V+b6KqcaSFaCRbf2kT3Fk6aFCABkankzPo41nvw1xXAHc65qcD5wC/N7ID4nHO3OeeWOOeWlJWVDTqA/QvRKjmTQPCvTWjzcwkm39pEV3KmypkERaaSs+3AtLTnUzmwHP1h4B4A59zfgTygdKQCiOflUBjNUbemBIV/bSKmbk0JJN/aRFe3pipnEhSZSs5eBOaa2Swzi+AN5Hyw1zlbgXMAzGw+XqMbfL/lQZgZ5Yk8Vc4kKPxrE5oQIMHkW5tQ5UyCJiPJmXOuA/gE8DjwBt5sm1VmdqOZXZw67TPAR81sJfAb4FrnRnYhpvLimCpnEgi+toncfAjlqnImgeJnm1DlTIImJ1Mf5Jx7FG8AZ/qxr6Q9Xg2cPpoxTCnOY3WFfhlJMPjWJsz2L0QrEiB+tYm8nDyi4ai2cJLAyJodAsCbFLC3oY2W9k6/QxHxV0xbOImki0fi1Lbpx7sEQ5YlZ95yGru0Abpku1iJxpyJpNHm5xIkWZWcTS5OLUSrSQGS7bT5uUgPqpxJkGRlcrZTkwIk26lyJtJDPBpX5UwCI6uSs65uTS2nIVlPEwJEekhEEqqcSWBkVXKWlxtmXEGEHaqcSbaLFUNrHXR2+B2JSCCociZBklXJGaCFaEUgbX9NVQpEwKucNXU00d7Z7ncoItmXnE0ujmnMmUietnASSRePegvRqmtTgiD7krNEHhU1qpxJluuunGncmQh4lTPQLgESDFmXnJUXx6hv7aC+RaVryWLa/Fykh67KmcadSRBkXXLWvZyGFqKVbKbNz0V6UOVMgiT7krPUchrq2pSs1p2cqXImAt4OAQC1rRpzJv7LuuSsvGuXAE0KkGyW530RKTkT8cQjqW5NVc4kALIuOZtYFCVkWohWslw4FyJFmhAgklIUKQJUOZNgyLrkLCccYmI8T5UzkZj21xTpEg6FKcotUnImgZB1yRnAlOIYW6oa/Q5DxF+xYk0IEEkTj8bVrSmBkJXJ2YmzxrFiW42W05DslqfKmUi6eCSuypkEQlYmZ2fOK6Mj6Xhu/V6/QxHxjzY/F+khEU2ociaBkJXJ2fEzSiiK5vDM2kq/QxHxT6xElTORNIloQpUzCYSsTM5ywyHeMreUZ9ZW4pzzOxwRf3RNCFAbEAG8bk1VziQIsjI5AzjzyDJ21bWwZle936GI+CNWAp1t0K5lZUQg1a3ZWqcf7eK7LE7OJgCoa1OyV5721xRJF4/E6XAdNHfoB4v4K2uTs4nxPOaXx3lm7R6/QxHxR9cWTpoUIAJoCycJjqxNzsDr2nxpyz7qtKSGZKOYKmci6bq2cKptU3Im/sru5KxrSY03taSGZKHuzc9VOROB/ZWzulZNChB/ZXVypiU1JKtpzJlID6qcSVBkdXKWGw5xxrxSnl2nJTUkC3VXzpSciYAqZxIcWZ2cAZw5b4KW1JDsFC0CC2tCgEiKKmcSFFmfnL31yDJAS2pIFjLbvxCtiBDLiZETylHlTHyX9clZ15IaT2tJDclGsRJNCBBJMTMSkYQqZ+K7rE/OQEtqSBbLU+VMJF08Gtc6Z+I7JWfAWUdOoFNLakg2ipVozJlImkQkof01xXdKzoDjpxdTlKclNSQLacyZSA/xaFxjzsR3Ss6AnHCIM+aW8sy6PVpSQ7KLxpyJ9KDKmQRBxpIzMzvXzNaa2Xoz+0I/57zPzFab2Soz+3WmYgNvSY3dda28sVNLakhmBKJNxEqgpRaSnSN+aZGhCkKbSEQTGnMmvsvJxIeYWRj4EfB2YDvwopk96JxbnXbOXODfgNOdc/vMbEImYuvSvaTGuj0smBzP5EdLFgpMm8grBpyXoOWPG/HLiwxWUNpEPBKnob2BjmQHOaGMfEWKHCBTlbOTgPXOuY3OuTbgt8Alvc75KPAj59w+AOdcRte2mBjPY0F5XOPOJFOC0Sa6dgnQpADxXyDaRDzq/Tivb1MvivgnU8nZFGBb2vPtqWPp5gHzzOw5M3vBzM7t60Jmdr2ZLTez5ZWVI5tIaUkNyaBgtImY9teUwAhEm+jaJUDjzsRPmUrOrI9jvUfe5wBzgTOBK4CfmVnxAW9y7jbn3BLn3JKysrIRDfJMLakhmROMNtG9v6YqZ+K7QLSJrv01Ne5M/JSp5Gw7MC3t+VSgoo9zHnDOtTvnNgFr8RphxnQtqaHdAiQDgtEm8lQ5k8AIRJvo3l9TyZn4KFPJ2YvAXDObZWYR4HLgwV7n/AE4C8DMSvHK1xszFB+wf0mNZ9dVakkNGW3BaBMacybBEYg20VU5U7em+CkjyZlzrgP4BPA48AZwj3NulZndaGYXp057HKgys9XA08DnnHNVmYgv3ZlHektqvLhZlQQZPYFpExpzJgERlDahbk0JgozNE3bOPQo82uvYV9IeO+DTqZtvLji6nG8/vpb/evQNfv/x0wiF+hoGITJ8gWgTOVHIzdeYMwmEILSJokgRoMqZ+Es7BPRSEM3hc+88khXbanhwZe/hDiJjkHYJEOmWG8qlILdAlTPxlZKzPrz3+KkcPSXBzY+toamtw+9wREZXnvbXFEkXj8RVORNfKTnrQyhkfOWiBeyqa+HWZzM6J0Ek82IlmhAgkiYRTWjzc/GVkrN+nDhzHBccU86tf9lARU2z3+GIjJ6YKmci6eKROLVt6tYU/yg5G8C/nXcUSQff/OMav0MRGT2xYo05E0mjypn4TcnZAKaW5HP9GbN5YEUFL21RZUHGqFiJKmciaVQ5E78pOTuIj585hwlFUW58eDXJpBamlTEorxg6mqG9xe9IRAIhHo1T11qnxcjFN0rODqIgmsPnzz2KldtqeGDlDr/DERl52iVApIdEJEFbso2WTv1gEX8oORuE9yyewjFTE3zzsbVaWkPGHu0SINJDPOrtr6lxZ+KXISdnZlZgZuHRCCaoQiHjKxd6S2v8VEtryFjTVTnTpAARwKucARp3Jr45aHJmZiEzu9LMHjGzPcAaYKeZrTKzW8xs7uiH6b8lM8dx4THl3PrsBnZoaQ0ZS/JUORNJp8qZ+G0wlbOngTnAvwGTnHPTnHMTgDOAF4CbzezqUYwxML5w3lEAfPaelerelLFDY85EelDlTPw2mOTsbc65rwO1zrlk10HnXLVz7nfOufcCd49ahAEytSSf/3z30SzdVMU1P19GXUu73yGJDJ/GnIn0oMqZ+O2gyZlzrisDub/3a2Z2Sq9zxrxLT5jKD684npXba7jyf1+gqqHV75BEhieaAEzJmUhKV+VM+2uKXwYz5ux9ZnYzUGRm83tNBrht9EILrguOKee2a5bw5u4G3n/bC+yq1XRrOYyFQtolQCRNQW4BYQtT26puTfHHYLo1nwPeAEqA/wbeNLOXzexhIGtHxp915ATu/NBJ7Kpt4bJbn2drVZPfIYkcurxiaK72OwqRQDAzb5cAJWfik8EkZxXOuTuBS5xz5znnZgNvA74KnA1gZjaKMQbWKbPH86uPnEx9SweX3fo8b+6u9zskkUMzYQFs+it0qJteBGBSwSS21W/zOwzJUoOarWlmNwDdf0udc9XAa8DJZnYn8MFRii/wjp1WzN3Xn0rSwftve4HXd+iXlhyGTvwwNO6BVX/wOxKRQFgwfgGrq1drCyfxxWCSs3OBTuA3ZlZhZqvNbBPwJnAF8F3n3B2jGGPgHTmpiHs/diqx3DDvv/Xv/PDJN2lo1VIbchiZczaUzoOlPwF9GYmwYPwCaltr2dGgbfsk8wYzW7PFOfdjvCRtBnAOsNg5N8M591Hn3IrRDvJwMLO0gPs+fiqnHVHKd55Yx1u/9TQ/++tGWto7/Q5N5ODM4KTroeIV2L7c72hEfLdw/EIAVlet9jkSyUZD2b7pZWAuUOmcqwEws3mjEtVhqjwR43+vWcL9/3ga88vjfOORNzjzlmf49dKttHcmD34BET8dewVE47DsVr8jEfHd3JK55IRyWFW1yu9QJAsNJTkbB/wI2GZmr5nZPcCjoxPW4W3x9BLu+sjJ/PqjJzO5OI8v3v8ab/vvZ7n/le10JtVlJAEVLYTFV8Oq+6Fup9/RiPgqEo4wr2SekjPxxVCSs63OubOcc+XAecD/AteOSlRjxGlzSvndx0/j9muXUBDJ4V/uXsmJ//ln/uXuFTywYgc1TW1+hyjS04kfgWQnvPQLvyMR8d2C8QtYXaVJAZJ5OUM4N2FmpwGrnHPbge2jFNOYYmacfdREzpw3gSfe2M0fX9/Fs+squf+VHYQMjp9ewllHTeDMI8tYUB4nS1clkaAYPwfmvROW3w5nfAZyon5HJOKbheMXct+6+9hev51p8Wl+hyNZZCjJWSHwWWChmUXxFqZ93Tn3uVGJbIwJhYx3LpzEOxdOojPpeHV7DU+vreSZtXu45fG13PL4WsYXRJg7sZAjJhRyRFkhR0woYu7EQiYURZW0Seac/DH45bu9ZTWOfb/f0Yj4pmtSwKqqVUrOJKOGkpyd6pzbAGBmecACYOGoRDXGhUPG4uklLJ5ewqffPo/K+laeXVfJsk1VrN/TwAMrKqhv2b8UR1E0hzkTCplcnEdpYTTtFqG0KEppQZRxhRHyc8OEQkriZJhmn7V/WY1j3ufN5BTJQkcUH0EkFGF11WrOnXWu3+FIFhl0ctaVmKUet+DN3nx5NILKNmVFUS49YSqXnjAVAOcclfWtrN/TwPrKBtbvaWBDZQNrd9XzXEMVtc397zOfHwlTEM2hMJrT43E0J0QkJ0QknLpP3aKp5+FQiJyQEQoZOSEjnHYfMu/ejO7HIfO6bENmGN73txkYlnq8/zikHe8KNHVu99Pu80g71n9ScLB8YaCXh5NrlCdiTC6OHfoFDgddy2o8+llvWY1pJ/od0YhzztGRdHR23ZwjmfSOJZOOpIOkcySdw3U/JvXcO+YgdZ96nvZ4/+ek7nG9nveMpfexvmMe1n/xcN48oGkl+UyI543a9f2UG87VpIBB6m4/vdtW6lh3G0qmtynvGBzYppLJA9sNaa+nHw9ae5oxvoDSwuENCRlK5UwyxMyYEM9jQjyP044oPeD1to4kVY2t7K1vY29DK5UNrdQ0tdHQ2kljawdNbR3djxtbO9hT30Jre5K2ziRtHftvrannMnifefs8bjhnrt9hjL5jr4Anb4SlP/U9OWtp72RfUxvVjW3UNLVT3djGvqY29jW209jm/X1vauukua1z/317By3tSdo7k7R3JGnrdLSn/r63dybp0KzpEfO1ixZw7emz/A5j1CwsXcgjGx8h6ZKEbChz6ILJOUddSwfVjW1pt1b2NbWnvjM6aW737rvaVmNbJ63tnd3fIfvbkvO+T7RUVA///b5jec/xU4d1DSVnh6FITojyRIzyxPArOM452jvTfu10OjqSSTqd94uno9P1qBr0rip0/brZ/6vH7b9P+1XjVRZcj+fdMXDgz5+BvjoP9ovHDfDu4U66mjm+YHgXOFxEC2HxB7w1z+q+AfHyUfuo9s4k2/c1s6WqkS1VTalbI1uqm6ioaaaprf+FnPNyQ+RHcojlhsmPeLdYJMyEojzycr1KcW44RG5O12PznocPrBSH0yrFXoV4f3U4lKoad1Vdux6nV4S917xzYX/lt6tQ210dTqsmd+ujctyX4Yw9Ha3O6TkTCkfpysGwcPxC7l57N9vqtzEjPsPvcAalvqWdzXub2FTVyKbKRjZXNbJxbyM7a5rZ19RGe2ff/xCGDPIjOd1tqetxIpZLXlG0R+9LVzvyHhs5oRDhEITT741Uuwr1aEMhM0Ih7x569rQYXT0zANajVyW9TQW5PR05qeiQr9tFyVmWMzMiORpTJH046SPwwo+9ZTXO+uKIXXZPXQt/31jFCxurWbqpii1VTT3W/4vlhpkxPp85ZQW8dV4Z4woijCuIUJKfS0l+6nFBhOJYLjnhw7+SIcG2YPwCAFbtXRXI5KytI8mKbTU8t34vSzdVsX5PI3sbWnucMzmRx6yyAs48sozxhVHGpdrRuMII41Pta1xBhFhuWJPPAkLJmYj0bdzsEVlWY099C0s3VqcSsio2VjYC3kSXk2aN44Kjy5kxvoAZ4/OZMT6fskLNTpbgmF08m2g4yqqqVZw/+3zvYEMlbHwapp8CxdMzGk8y6XhjVx3Pr6/iuQ17Wbapmqa2TkIGi6YkOOeoCcwsLWBW6jZjfD55ueGMxijDp+RMRPrXvazG/XDs5YN+m3OOFzZW89NnN/DsukoAClPJ2OUnTuOU2eNZODlBWLOLJeByQ7kcOe7I/Xts7lkDv7oMard6z8uPgwUXw/yLoXT0xqNu39fEz/66iQdXVlDd6C1gPqesgEtPmMppc0o5dfZ4Evm5o/b5kllKzkSkf7PPgtIj4YWfwNHvg9DA3YjJpONPq3fz02c3sGJbDaWFET71trmceeQEFk2OqxtSDksLxi3gwQ0Pktz0F0J3Xw3hKFxxN1SugTce8ibPPHkjlB0F8y/ybpOOGZFlaN5A6IsNAAAgAElEQVTcXc9Pnt3AgysqADjv6HLOnFfG6UeUMikxNmfJipIzERmIGZzycXj4U/CTU+H0f4ZFl0JOpMdprR2dPPBKBT/9ywY2VjYyfVw+33jXIi49Yaq6VOSwt7B0Ib9d+1s2/+YyZsdnwFX3QskMOPJceMunoHYHrHkE3ngQ/vod+MstEJ8K894Bc98Js/4BIvlD+swV22r48dPr+dPq3cRyw1xz6kw+csassb+UjwBKzkTkYE64FiKF8Lfvwh8+Dk99A079Jzj+g7hIAb9etpUfPPkmu+taWVAe54dXLOa8RZNUJZOxwTkWblkOwOry+cy+/H6IlfQ8JzEFTr7euzXuhbWPwbo/wsq7vTGbOXlegjb3Hd44zgHGqf19QxU/fOpNnt9QRSKWyyfPmcu1p81kXEGk3/fI2JOx5MzMzgW+D4SBnznnbu7nvEuBe4ETnXPLMxWfSKYdNm3CDI65DI6+FNb/Gf72PXj8i7hnv8XfSt7NdzadzBGzZnLLpcdyxtxSDeaXQxa4NtHZAY9+hlkv3UHerJmsWvAOLuydmPVWUArHf8C7dbTCludg3ePe7c0/eYs7j58L5cfAxIUw8WiYuBBXVM73n1rP9/78JhPjUb50wXwuP2k6hVHVULJRRv5fN7Mw8CPg7Xgbpr9oZg8651b3Oq8I+CSwNBNxifjlsGwTZjD37TD37XRsWcqa+77OGTt/wdLYr8nJPwNbPxuqZnmzPMfN9rp9tHG6DFLg2kRrA9x3Hbz5J3Le8mmOalvP6uq1Q7tGThTmnO3dzr0ZqtZ7FbXNz8H2F+H133Wf2hQq4pSOqfxy8gJOOX4xucV1sGs3FE3ybrnqzswmmUrJTwLWO+c2ApjZb4FLgNW9zvs68C28DdZFxrLDtk20dnTyz3/J4Y+VH+erp36Sa0OPYhWvwPZl0FqXdqZBYprXhZOXSN3i3n00vv95boH3JZaT18d9BEK5EMpJ3cLa63PsClabuO86r1J84XdhyYdYuOxmfv/m7+lMdhIOHcI4SjNvNmfpXDjtBu9YSy01m1dwzyOPU7BvDWePq2RSw+PYn3934PtjJVBUDgVlXruJpm55cYgWpZ4XQW6+135yYwe2p3AUwl1tqatdafhBEGUqOZsCbEt7vh04Of0EM1sMTHPOPWxm/TY6M7seuB5g+vTMri8zKjrboaUOWmq8L7b2Fuho9u7bm6CjBdqbvfvONq/Mnmz33pfs2H+fbMdbrj8JrhOSnWn3ydQS/UnA9Xzefcyl3ad0Het+TP/PexzjwNeGarhL+Y+WJdfB4qtH4kqHZZtobO3g/931En99cy9fvnAB171lFvAO70XnoKkaqjf2vNVug32boaXW+zveI4E7BN2JWg5Y2Ptysa5b2LvvTuK6Nn0N9f24eylx6+NY1wemrzxuBx7rcbyP13obVnIZwMT01H+CRe8ZiSsFq01sXQrHfxCWfAjwFqP91Ru/YnPdZuYUzzm0a/byZm2I6x7opLL+DL77/hsoP7rca0ctNVC/C+oqvPv6rvtd0LDHu7XWed8dbfXDC8JCaclaqt30aEdp7cpIaz+hvttTehvpPtbjA1N3A7SZw709nfkFr5dhGDKVnPX1J9D97WtmIeC7wLUHu5Bz7jbgNoAlS5YE8xs8mYTGSqjd7n0x1W7f/7ipyvuS6rq1NRzCBxiEc/c3qHB6wwqlNapwr8ZFr4YVOvBLqb8vrJ47ku+Po8fztGMHnHsoAvhFFB6xbrrDrk3UNrVz7R3LWLmthlsuPYbLlkzreYIZFIz3bgPtx5nshNb6/cla14+PjtZe9y3Q0Zb68dHr1tm+/wdI14+M7h8ind4XXLKTnj88kn08hh4/Svr9EdLfsbTjfb7W2xj8wTJyXdfBaRMdrdBaC/Ep3YcWjl8IwKqqVSOSnP1lXSX/9KuXyYuEuedjp3LstGLvBTOvShYrgQnzD36hZNJL0FrqvHbV0dyzLbWnPe9s6/WjvqPn8wPaUddzt79d9f6B37ttQa/21Ff7GOSP+sO1PYWHv95cppKz7UD6v+RTgYq050XAIuCZ1GDiScCDZnZx4CcFdLbDrldh2zLYthR2rvQSsc62nudFCr2GXjjBG4+TV+x168RS93nFqZJ0LFWOzut1H02VpFNJmBzuDqs2saeuhWtuX8bGykZ+fNUJnLto0qFfLBT2/t7HikcuQBkLgtMmGvd69wWl3YdmxmcSy4mxumo1F8+5eFiX/+ULW/jag6uYO6GQn197IlOGszxGKLR/2ICMGZlKzl4E5prZLGAHcDlwZdeLzrlaoLsVmNkzwGcDmZi11MKWv3uJ2LalsONl75cKeONrphzvrRSdmOo9T0z1bnkJjZWRdIdNm9hT18Jlt/6dyvpWbr/2RN4yt/TgbxIZuuC0iUZvVwsKyroPhUNh5o+bz6q9q4Z16Z//bRNff3g1Zx81gR9csVizMaVPGflb4ZzrMLNPAI/jTZG+3Tm3ysxuBJY75x7MRBzD0lIHf/+Rd2ur97oUy4/1xiBNOwmmnuStdSMyCIdTm/jJsxuoqGnm7o+dyvHTD7KMgMghClSb6K6clfU4vGD8Au5bdx8dyQ5yQkP/+ly/p4Fv/nENb5s/kVs/cIK2L5N+ZSxld849Cjza69hX+jn3zEzENCjtzfDiz+Cv/w3N1bDgEjjpeph8/JBXfBZJdzi0iYbWDu5dvp0Lji5XYiajLjBtounAbk3wdgq464272FS7ibklQ9tHszPp+Px9K8mPhPmv9yxSYiYDUj21P53t8Mpd8Oy3vJkyc86Bc74Mkxf7HZlIxvz+5e00tHZw7emz/A5FJHP66NYEr3IG3qSAoSZnv3huEy9vreF77z+OCUXaE1MGpuSsN+e8hQGf+gbs2wTTTob3/i/MfIvfkYlkVDLpuOP5zRw7rZjjpmnwvmSRxkpvAla0qMfhmfGZ5Ofks2rvKt51xLsGfbmNlQ3c8vha3jZ/ApccN3mko5UxSMlZby/fCQ/9s7elxpX3eHuhaSC/ZKG/rt/LxspGvvf+4/wORSSzGvd6VbNe//aHLMSC8QtYXd17Xdz+ed2ZrxLNCfGf7z5a25vJoGhp4HQttfDk12H6afCxv3gb1KohSZa647lNlBVFOf/ocr9DEcmsxsoDxpt1WTB+AWur19KebB/Upe58fjPLt+zjKxctZGJc3ZkyOErO0v3lFm+R2HNv0pYWktU27W3k6bWVXHXydCI5aguSZRorDxhv1mXh+IW0draysWbjQS+zeW8j33p8DWcdWcZ7j9dsfhk8/avbpWoDvPBTWHwVTFY3jmS3//v7ZnLDxpUnj4Et0kSGqqtbsw8LS72dAlZXDdy1mUw6Pv+7V8kNh7jpPceoO1OGRMlZlz992VuF/+w+Z22LZI305TM0q0yyjnMDdmtOK5pGYW4hq6oGXoz2ly9sYdmmar58wQImJdSOZGiUnAFsfAbWPgJnfAaKJvodjYivfveSls+QLNbW4O1D2U/lrHtSwACVs61VTXzzj2v4h3llXLZk6mhFKmOYkrPODvjjF6F4Opzyj35HI+KrZNJx5/ObOU7LZ0i26meNs3QLxy9kTfUa2nrvoZzyxftfI2TGze/R7Ew5NErOXvk/2LMK3v51yFXpWbLbX96sZOPeRq47fabfoYj4o5+tm9KdXH4y7cl2/rzlzwe8trqijr+t38s/nzOXycPZ0FyyWnYnZy213mKzM073tmUSyXJ3Pr+ZsqIo5y3S8hmSpborZ32POQM4dfKpzIzP5K437jrgtXtf2kYkHOLSE9SdKYcuu5OzZ78FTdXe0hkqPUuW0/IZIgyqWzNkIa6cfyWv7X2NlZUru4+3dST5wys7ePuCiZQUREY7UhnDsvdf4KoNsPRWWHw1lB/rdzQivrvzeS2fITKYyhnAJXMuoSi3iF+t/lX3sSff2M2+pnZNApBhy97k7E9fSi2d8WW/IxHxXX1LO/e9tJ0Lj5ms5TMkuzXuhWjc+34YQH5uPu+e+26e2PIEuxt3A3DP8m1Miudxxtz+q24ig5GdydmGp2Hto1o6QySle/mM02b6HYqIvwZY46y3K466giRJ7l57N7vrWnh2XSXvPWEK4ZCGycjwZOfG5899HxLTtHSGSMr9r+zg2KkJjtXyGZLtBti6qbepRVM5a9pZ3LvuXqg5h6SDS0+YNsoBSjbIzspZ9QaYfqqWzhBJ2VzVxDFTlZiJDLR1U1+umn8VNa01/Gb1A5w0cxyzSgtGMTjJFtmXnCWTULcTEtqEVgSgqa2D2uZ2yov1Y0VkKN2aAEsmLmFawRHU5j7FpSfoe0VGRvYlZ417INkOcTUiEYCKmhYAJie0YKZkuWQnNFUNqXJmZhS3n0U4bzcTJ24fxeAkm2Rfcla7w7tPaKqzCMDO2mYArWYu0rwPXHJIyVljawcr18wiQpz71v9mFIOTbJKFydk2717JmQgAFTVeclaeULemZLnurZsG36356Gs7aWwN8c7p7+LZbc+yrW7bKAUn2ST7krO6VOVM3ZoigNetaQaTlJxJthvE7gC93fvSdmaVFvCpkz5IOBTm12t+PUrBSTbJvuSsdgfk5kOsxO9IRAJhZ20zZYVRcsPZ98+BSA9DTM427W1k2aZqLlsylQkFE3jnzHdy//r7aWhrGMUgJRtk37/Gtdu8qpn20hQBvMqZxpuJkNatObjk7L6XthEyeO/x3jCZq+dfTWN7Iw9seGC0IpQskX3JWd0OjTcTSVNR28xkLaMh4lXOLDSonpXOpON3L+3grfPKmBj32s+i0kUcV3Ycv37j1yRdcrSjlTEs+5Kz2h1a40wkxTnHzpoWyrWMhoiXnOWPh1D4oKf+9c1KdtW1cNmSnjsCXLXgKrbWb+XZbc+OVpSSBbIrOetog4bdEFflTASgtrmd5vZOdWuKwJC2brr3pe2U5OdyzvwJPY6/bfrbmF40nW+88A32NO0ZjSglC2RXclZfATh1a4qk7EgtozFZMzVFUls3HXwZjX2NbTyxajfvWjyFaE7PKltOKIfvnvVd6tvr+dTTn6Klo2W0opUxLLuSs+4FaNWtKQKwM7U7QLkqZyKDrpw99vou2jqTXHpC3z/055XM46YzbuK1va/xtb9/DefcSEcqY1x2JWfda5ypciYCabsDqHImMuhNz5/fsJdJ8TwWlMf7Peec6edww+IbeGTjI9z++u0jGaVkgexKzmpT+56pciYCQEVtC7lho7Qw6ncoIv7qaIXW2oN2azrnWLapmpNmjcMOsiTTR4/+KOfOPJfvv/x9TRCQIcm+5CxWApECvyMRCYSKmmYmJfIIhbTun2S5Qa5xtrmqiT31rZw8e9xBL2lm3Hj6jcwfP59//eu/sn7f+pGIVLJAdiVndTvUpSmSRstoiKQMcneApRurADh51vhBXTaWE+P7Z32fvHAeNzx1AzUtNcMKU7JDdiVnWuNMpIeK2maNNxOBQVfOlm6qprQwwpyywffATCqYxPfP/j67m3bzmWc/Q3uyfTiRShbIWHJmZuea2VozW29mX+jj9U+b2Woze9XMnjSzGSMeRO02LaMhgeF3m+hMOnbVausmCQ5f20R35WzgMWeDHW/W27Flx/LVU7/Ksl3LuGnpTXQmOw81UskCGUnOzCwM/Ag4D1gAXGFmC3qd9gqwxDl3DHAf8K0RDaKtEVpqvH01RXwWhDaxt6GVjqTTMhoSCL63ia7kLL//5GxbdRM7apoH3aXZ2yVHXMJ1i67j3nX38sE/fpBNtZsO6Toy9uVk6HNOAtY75zYCmNlvgUuA1V0nOOeeTjv/BeDqEY2ge40zVc4kEHxvExVpC9CurFzJq5Wv9ntufk4+75j5DooiRSMZgkg6f9tEYyWEoxDt/+/40k3VAIOaDNCffzn+X7x10JbexGUPXcYNi2/g6vlXEx7EllGSPTKVnE0BtqU93w6cPMD5HwYe6+sFM7seuB5g+vTpg4+gLrWMhipnEgy+t4mdtd4CtJOLY3z2r19ka/3WAc//9vJvc9mRl/GB+R+gLH9wW9yIDIG/baJrjbMBuiuXbqyiOD+XeRMO/UeKmXHh7As5pfwUbvz7jXx7+bd5YssTfP30rzMrMeuQrytjS6aSs77+tve5ZLKZXQ0sAd7a1+vOuduA2wCWLFky+GWXu9c4U+VMAsH3NtFVOSspgG312/jo0R/l2kXX9nnutrpt3LnqTu5cdSd3rb6Li+dczHWLrmNGfOSHhkrW8rdNNFYefLzZ5mpOnDluRJaeKY2V8v2zvs8jmx5RFU0OkKkJAduBaWnPpwIVvU8ys7cB/w5c7JxrHdEIancABvHJI3pZkUPke5uoqGkhPxJmb+tWHI6F4xcSj8T7vC0sXci33votHn7Xw7xn7nt4aMNDXHT/RXz6mU+zau+qkQxLspe/beIgWzftqm1hS1UTJ8869C7N3rqqaA+86wFOm3wa317+bT74xw/y4q4XSbrkiH2OHH4ylZy9CMw1s1lmFgEuBx5MP8HMFgO34jW4PSMeQd12KJwI4dwRv7TIIfC9TeysbaY8kcf6Wm9hzDnFcw76nmnxaXzplC/x+KWP85GjP8ILFS9w+SOXc8XDV3D767ezrW7bQa8h0g9/20RT1YDJ2dJN3vpmp8w+tMkAA+mqot18xs1srtvMhx7/EOf97jx++MoP2VK3ZcQ/T4IvI8mZc64D+ATwOPAGcI9zbpWZ3WhmF6dOuwUoBO41sxVm9mA/lzs0tTvUpSmBEYQ2UZFaRmP9vvVEQhGmFU07+JtSSmOlfPL4T/KnS//E55Z8DoDvvvRdzr//fC576DJue/U2NtZuHMlwZYzztU04d9BuzRc2VlMUzWH+APtpDoeZccHsC3ji0ie4+YybmZWYxc9e+xkX3n8hH3j0A9yz9h7q2upG5bMleMy5wQ/bCpolS5a45cuXD+7kHy6BiQvgff83ukFJ1jKzl5xzS/yMYSht4sT//DNnHzmB2sSPqW6p5p6L7hnWZ1c0VPDnLX/miS1PsKJyBQBzEnM4e/rZLCxdyPxx8ykvKB/y+lBy+Dps2kRrPdw0Fd7+dTj9k32ecs53nmH6uHx+cd1JoxBl33Y37uaRTY/w4PoH2VC7gUgowmmTT+OYsmNYOH4hC0sXkogmMhaPDN9g20SmJgT4yzlv66a57/A7EpFAaOtIsrehlfLiPJbve5OTJw00KW5wJhdO5pqF13DNwmvY3bibJ7c+yZ+3/pmfv/7z7vEz8Uico8Yd1X2bP24+0+LTiIa18br46CBbN1XWt7KhspHLlgy+ujwSJhZM5EOLPsR1C69jdfVqHlz/IM9VPMcz25/pPmda0TQvUUsla9OKplEWK9OkgsNcdiRnzfugvUlbN4mk7K5rwTkoKepgz449HFFyxIhef2LBRK6cfyVXzr+S5o5m1u1bx5qqNbxR/QZrq9fy2zW/pS3Z1n1+SbSECfkTmFgw0bvPn8jE/ImUxkopihR13wpzC4nlxFR9k5F1kK2blnWtbzaCkwGGwsy6EzCAurY6VletZtXeVayqWsXKypX8cfMfu88PW5iy/DIm5U9iUsEkygvKmVgwkbJYGUWRou6JPkWRIgojheSEsiMVOJxkx/8jWkajh6RL0pHsoD3Z3n2fdEmSLkmn6ySZTNLhOrqPJV0Sh8M5R5IkzvV8DHS/nn7fJf2c9HO7X+86N62H3fU9g/6gr422GUUzmBbP7K/n0bAjtYyGy9kNwBHFI5ucpYvlxDi27FiOLTu2+1hHsoNNtZtYU72GioYK9jTtYXfTbvY07eH1va9T3VLd7/VyLIfCSKGXqOXGiIaiRHOi5IXziIQj3ffRcJRwKExOKMe7Wc7+x6EcwhbGMMKhMCELESJEKBTqPm5mhCzU/djY/9z7n5cgpr/e9VqX9HO6j/VKLAd6bSisz5UoRsZAcc1OzGZy4WE+C/4gWzct3VRFfiTMoinB6EKMR+KcUn4Kp5Sf0n2sqrmKNdVr2NGwg12Nu9jdtJudjTtZVbWKp7Y+1ePHUG8FuQXdP3yi4f3tKRqOkpfj3UfCEXJDuX22pa72FLJU+zHr8Txkof1tJO0+xP7jXW1qLLSneSXzmFgw8ZDfD9mSnNWldgeIH57JmXOOfa372Neyj9rWWu/WVtv9uK6tjvq2elo6WmjpbOl5n3rc3tlOe9K7dTrt6XaoPnHcJ/jYsR/zO4xh21nrJWeNzptdObd4bkY/PyeUw9ySucwt6ftz2zrb2NO0h73Ne2lob6C+rZ76tvoDHnf9/W7rbKO+rZ7WztYet07XSUeyo/umv/sj7wsnfYGr5l/ldxjDc5BuzWWbqjlhRgm54YxtRz1k42PjOX3K6X2+5pyjuqWavc17u9tPfbt3X9dat/87pLOF1o5W776ztUebauts29+W3P4f93Kg/3rLf3HRnIuGdY3sSM66K2fB7dZs7mhmQ80GdjTsoKKhgh0NO7ofVzRU0NLZ0uf7QhYiEUlQGPF+9eTl5BELxyiKFJGXk0deOI+8nDxyQ7ndv3pyQ7nkhnPJsRxyw7kH/MIJh8KELdzjF1B/lYTev3B6H+vvV033r6Wu43bgL6KB+NWtNTF/eL+GgqKixvv7VNW2lYLcAiYVTPI5op4i4QhTi6YytWhkf1AlXZLOZCftyXYcrrtSnCTZ/VrX8d6V4PQKco+KMe6A6jAcWDHu/Th1oP/Xehlo8tZoVpMPdu0phcH9d3XQBqic7WtsY82uei48pjzDQY0cM2N8bDzjYyO7DIhzrvsHUHfPS3ovTK/el6RLgmN/D0wfPS5joT0NZeZ7f7IjOavbAaFcKJjgdyTddjXuYsWeFayoXMGKPStYW72WDrf/V0g8EmdK4RRmJ2bzlilvYXLhZMbnjScejZOIJkhEEiSiCQpyCwhZcH/NSTDtrG2mJD+XzfUbOaL4iKwZwxWyEKFwiFytdyjpGvdCNAE5B05MWba5az/NkV/f7HBnZt1dnDKysuNPtHa7tzNAyL8kpr6tnsc2PcbSnUtZWbmS3U3eWJ9YToxFpYu4dtG1LBq/iKlFU5lcOFkbTMuoqqhpYVIijzf3vck508/xOxwRfw2wxtnSjdVEc0IcMzUY480kO2RJcubPArTOOVZVreLedffy2KbHaO5oprygnOMnHs+xZcdy3ITjmFcyj9yQfsVLZlXUNDOxpI3trTWjOhlA5LAwwNZNyzZXsXh6MdEcLU0hmZMdyVnddph2ysHPGyFN7U08sukR7l17L29Uv0EsJ8b5s87nsnmXsbB0YcbiEOnPztoWZk6thkZGfBkNkcNO414YN/uAw3Ut7ayuqOOGszM7YUZk7CdnyU6oq8hI5WxP0x5ue/U2Ht74MI3tjcwrmce/n/zvXDD7AnVTSmA0tnZQ29wOkV1ecqbKmWS7xkqYduBCzMs3V5N0cPJsf9Y3k+w19pOzhj2Q7Bj1mZpPbn2Srz7/VZrbmzl31rm878j3cUzpMVkz0FoOH13LaLTYDkqiJYzP00BnyWLJzn43PV+6sZrcsLF4WokPgUk2G/vJ2Sivcdbc0cwtL97CvevuZf64+XzzH77JrMSsUfkskZHQtYzGvvatHFGSPTM1RfrUvA9csu/kbFM1x04tJhbReDPJrLG/BsMornG2pnoN73/4/dy77l6uW3Qdvzr/V0rMJPC8ypljZ9NmdWmK9LPGWWNrB6/tqFWXpvhi7FfORmHrpqRLctfqu/jey9+jOFrMbW+/jVMnnzpi1xcZTTtqWgjl1tDc2aTkTKSf3QFe2rKPzqTj5Fnq9pfMG/vJWd0OyC2AvOIRudze5r38+9/+necrnuesaWfxH6f9ByV5Go8gh4+dNc2MK66iFfrdPkkka/STnC3dVEU4ZBw/Q/++S+aN/eSsdrvXpTkC42qqmqu46pGrqG6p5sunfJnL5l2m8Tpy2NlZ20JBkZeczSme43c4Iv5q3Ovd9+rWXLqxmkVTEhRGx/7XpATP2P9bVzcyC9C2d7bz6Wc+TVVLFb945y84uuzoEQhOJPMqapvJKdvNxNhE4pG43+GI+KuxEiwEsf0Vsr0Nrby8dR//eKa6/cUf2TEhID78yQA3LbuJl/e8zI2n3ajETA5bzjkqapppC1VovJkIeMlZ/ngI7Z+R+dhrO0k6uOjYyT4GJtlsbCdnHW3eOmfDrJzdveZu7l13Lx9a9CHOn33+CAUnknk1Te20tHdQn9yh5EwEvG7NXuPNHlq5k3kTCzlykhYPF3+M7eSsvgJww6qcvbjrRW5edjP/MPUf+OTiT45cbCI+qKhtxiJVdLp2bdskAqnkbP94s521zSzbXM2Fx6hqJv4Z28nZMJfR2NGwg8888xmmxadx8xk3Ew5pIUI5vFXUtBCO7gZgbrFmaor03vT8kVd3AnDhMeV+RSQy1pOz1O4Ah5CcNbU38cmnPklHsoMfnPUD7Y0pY8LO2mZC0V0YpgWTReCAbs2HXt3JoilxZpcV+hiUZLuxnZzVpSpnQ+zWTLokX3ruS6yvWc8tb72FmYmZIx+biA8qalrIydvN1KKp5Ofm+x2OiL86WqG1trtbc2tVEyu31XCRujTFZ2M7OavdAbFxEBnal9Ctr97KE1ue4NMnfJrTp5w+SsGJZN7O2mZyY3s0GUAE0tY48ypnD71aAcAF6tIUn43x5Gz7kPfUXLpzKT9e8WMunnMx1yy4ZpQCE/HH9po6kuFKJWcicMDuAA+trOCEGSVMLVFVWfw1tpOzuh0QH/x4s45kBzcvu5mphVP5yqlf0er/MubsaNgKltS2TSLQo3L25u561uyq10QACYSxnZwNsXL2u3W/Y33Nej675LNEw9FRDEwk8zqTjn3tWwFUOROBtMpZKQ+9uhMzuOBoJWfiv7GbnLU2QEvNoGdq1rXV8T8r/ocTJ53I2dPPHuXgRDJvb0MrLrKbEGFmxmf6HY6I/1LJmcsv5eFXKzhl1ngmxPN8DuQdn5MAAAeOSURBVEpkLCdndallNAbZrfnTlT+ltrWWz5/4eXVnypi0o8ZbRmNCbCq54Vy/wxHxX2Ml5OSxuirJxspGbdckgTF2k7PuBWgP3q25uXYzv3njN7xn7ns4atxRoxyYiD92phagnRWf43coIsGQWuPsoVd3kRMyzl00ye+IRICxnJzVDX4B2u8s/w7RnCifWPyJUQ5KxD9b9+0jFKlmYemRfociEgyNlbiCUh5aWcFb5pYyriDid0QiwFhOzmq3AwZFAw/ufL7ieZ7Z/gzXH3M9pbHSAc8VOZytqV4PwMKyeT5HIhIQjZXUhYrZUdOsvTQlUMZwcrYDiibBAGNrOpId3PLiLUwtnMrV86/OYHAimbe1fiOAltEQ6dK4l43NMSLhEO9YONHvaES65fgdwKip237QbZvuW3cf62vW870zv0ckrHK2jG2VrZux3FymFg59r1mRMcc5XGMlr7kIZx5ZRjxPk2QkOMZ25WyA8Wa1rbX8aMWPtHSGZI365HaKwlMIh8J+hyLiv9Z6rLOVra0FmqUpgZOx5MzMzjWztWa23sy+0MfrUTO7O/X6UjObecgf5lxqAdr+k7OupTP+9cR/1dIZ4otMtonWjk46wjuZGD3kS4iMuox+T6TWOKsPF3PO/AmHfBmR0ZCR5MzMwsCPgPOABcAVZrag12kfBvY5544Avgt885A/sHkfdDT32625qXYTv13zW94z9z0cOU4z1yTzMt0m1lfuIZRbx/T47EO9hMioynSb6KjfA8D06TPIj4zdET5yeMrU38iTgPXOuY0AZvZb4BJgddo5lwBfSz2+D/gfMzPnnBvqh72+9Vk+PGMqrP85bLzzgNfbk+1aOkP8ltE28dKuNQAcNU6TASSwMtom1mzYxCJg8VHaykyCJ1PJ2RRgW9rz7cDJ/Z3jnOsws1pgPLA3/SQzux64PvW0wczWDvC5pb3fn66MskEFP8IGjMknimnwBoprxhCu40ebKP04F+39+BCCzJAg/n+tmAbvsG4T/MfbDrc/U78opsEbdpvIVHLW16Cu3r90BnMOzrnbgNsG9aFmy51zSwZzbqYopsEJYkwwonFlvE1kwZ/piFFMg6c2MfKCGJdiGryRiCtTEwK2A9PSnk8FKvo7x8xygARQnZHoRDJPbUKkJ7UJkZRMJWcvAnPNbJaZRYDLgQd7nfMg8MHU40uBpw5lHIHIYUJtQqQntQmRlIx0a6bGBnwCeBwIA7c751aZ2Y3Acufcg8DPgV+a2Xq8X0KXj8BHD6r7M8MU0+AEMSYYobh8ahNj+s90hCmmwVObGHlBjEsxDd6w4zL96BAREREJjrG7Q4CIiIjIYUjJmYiIiEiAjMnk7GBbgPjBzDab2WtmtsLMlvsYx+1mtsfMXk87Ns7MnjCzN1P3JQGI6WtmtiP157XCzM7PcEzTzOxpM3vDzFaZ2T+njvv6Z3Wo1CYGjENtYnAxqU2MMrWJIcc0ZtvEmEvObHBbgPjlLOfccT6vy3IHcG6vY18AnnTOzQWeTD33OyaA76b+vI5zzj2a4Zg6gM845+YDpwD/lPp75Pef1ZCpTRzUHahNDIbaRGaoTQw+JhijbWLMJWekbQHinGsDurYAEcA59xcOXBfoEqBrn6s7gXcFICZfOed2OudeTj2uB97AW53c1z+rQ6Q2MQC1icFRm8geahODM5ptYiwmZ31tAdL3DuiZ5YA/mdlL5m0tEiQT3f9v745Vo4iiOIx/p7DSQFAQ0gXTWqitYmmRziewsFMLe18gINpaWIpYqQi+goVdjBFBsTUkYKOl6LGYiWYFww6ue4+z3w+GCQkLJ/fOH87e2dmbuQPdxQacbFzPvhsRsdUvZze7VRIRq8BZ4CV1x+owZmK4qvNsJmbDTAxXdZ5HmYkxNmdTbe/RwPnMPEe3jH49Ii62Lqi4e8AacAbYAe60KCIijgGPgZuZ+blFDTNgJsbBTMyOmRiH0WZijM3ZNFuAzF1mfuzPe8BTumX1KnYjYgWgP+81rofM3M3Mb5n5HbhPg/GKiCN0gXuYmU/6X5cbqymYieHKzbOZmCkzMVy5eR5zJsbYnE2zBchcRcTRiFja/xm4BGwf/qq5OrglyhXgWcNagJ8X9L7LzHm8IiLovo38bWbePfCncmM1BTMxXLl5NhMzZSaGKzfPo85EZo7uANaBd8AH4FaBek4Br/rjTcuagEd0y79f6d49XgVO0D1R8r4/Hy9Q0wPgNbDVX+grc67pAt1tji1gsz/WW4/VX/w/ZuLPtZiJ6WoyE/+2HjMxvKbRZsLtmyRJkgoZ421NSZKk/5bNmSRJUiE2Z5IkSYXYnEmSJBVicyZJklSIzdkCi4jliLjWug6pCjMhTTITbdicLbZlwNBJv5gJaZKZaMDmbLFtAGsRsRkRt1sXIxVgJqRJZqIBv4R2gUXEKvA8M083LkUqwUxIk8xEG66cSZIkFWJzJkmSVIjN2WL7Aiy1LkIqxExIk8xEAzZnCywzPwEvImLbD3pKZkL6nZlowwcCJEmSCnHlTJIkqRCbM0mSpEJsziRJkgqxOZMkSSrE5kySJKkQmzNJkqRCbM4kSZIK+QFCBSxqSNVYnQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa834a73f98>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "T = 20\n",
    "\n",
    "def PI(pi0,P=P,T=T):\n",
    "    # Function PI computes the state probability vectors\n",
    "    # of the Markov chain until time T\n",
    "    pi_ = array([pi0])\n",
    "    for i in range(T):\n",
    "        pi_ = vstack((pi_,pi_[-1] @ P))\n",
    "    return pi_\n",
    "\n",
    "def plot_PI(x0):\n",
    "    # subplot(1,3,n+1) of successive states probabilities \n",
    "    # with initial state x0\n",
    "    pi_0       = zeros(3)\n",
    "    pi_0[x0-1] = 1\n",
    "    pi_  = PI(pi_0)\n",
    "    subplot(1,3,x0)\n",
    "    plot(pi_)\n",
    "    xlabel('t');axis(ymin=0,ymax=1)\n",
    "    if x0==1: ylabel(r\"$\\pi(t)$\")\n",
    "    if x0==2: title(\"Evolution of $P(X_t)=1,2,3$.\")\n",
    "\n",
    "rcParams[\"figure.figsize\"] = (10., 4.)\n",
    "for x0 in range(1,4):\n",
    "    plot_PI(x0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**3)** To compute the steady state distribution $\\pi^*=[\\pi^*_1,\\pi^*_2,\\pi^*_3]$, we must solve the system of load balance equations $\\pi^*=\\pi^* P$ with the normalization condition $\\pi^*_1+\\pi^*_2+\\pi^*_3=1$. The system of equations $\\pi^*=\\pi^* P$ is redundant : the third equation is a straightforward linear combination of the first two ones. Taking into account the normalization condition $\\pi^*_1+\\pi^*_2+\\pi^*_3=1$ and discarding the third redundant equation in $\\pi^*(P-I_3)=0$ yields a full rank system of equations. Complete the code below to solve this system and obtain the steady state ditribution. We will use the [solve](https://docs.scipy.org/doc/scipy/reference/generated/scipy.linalg.solve.html#scipy.linalg.solve)\n",
    "function from the **scipy.linalg** library."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "steady state distribution: pi* = [ 0.44736842  0.39473684  0.15789474]\n"
     ]
    }
   ],
   "source": [
    "from scipy.linalg import solve\n",
    "####################\n",
    "# complete the code to get the steady state distribution\n",
    "# of the discrete time Markov chain\n",
    "pi_ = solve([[-.3, .3, 0.1], [.3, -.5, .4], [1, 1, 1]],[0,0,1])\n",
    "print(\"steady state distribution: pi* =\",pi_)\n",
    "####################\n",
    "V2,V3 = pi_[0],pi_[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Your answers for the exercise"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "---------------------------\n",
      "RESULTS SUPPLIED FOR LAB 3:\n",
      "---------------------------\n",
      "V1 = 1.72\n",
      "V2 = 0.45\n",
      "V3 = 0.39\n"
     ]
    }
   ],
   "source": [
    "print(\"---------------------------\\n\"\n",
    "      +\"RESULTS SUPPLIED FOR LAB 3:\\n\"\n",
    "      +\"---------------------------\")\n",
    "results = (\"V\"+str(k) for k in range(1,4))\n",
    "for x in results:\n",
    "    try:\n",
    "        print(x+\" = {0:.2f}\".format(eval(x)))\n",
    "    except:\n",
    "        print(x+\": variable is undefined\")   "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}