Notebooks/MarkovChainExample.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Markov Chain from scratch\n",
    "\n",
    "Last update: Sept 24th, 2020"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example\n",
    "\n",
    "Consider a worker who, at any given time $ t $, is either unemployed (state 0) or employed (state 1).\n",
    "\n",
    "Suppose that, over a one month period,\n",
    "\n",
    "1. An unemployed worker finds a job with probability $ \\alpha \\in (0, 1) $.  \n",
    "1. An employed worker loses her job and becomes unemployed with probability $ \\beta \\in (0, 1) $.  \n",
    "\n",
    "\n",
    "In terms of a Markov model, we have\n",
    "\n",
    "- $ S = \\{ 0, 1\\} $  \n",
    "- $ P(0, 1) = \\alpha $ and $ P(1, 0) = \\beta $  \n",
    "\n",
    "\n",
    "We can write out the transition probabilities in matrix form as\n",
    "\n",
    "\n",
    "$$\n",
    "P\n",
    "= \\left(\n",
    "\\begin{array}{cc}\n",
    "    1 - \\alpha & \\alpha \\\\\n",
    "    \\beta & 1 - \\beta\n",
    "\\end{array}\n",
    "  \\right)\n",
    "$$\n",
    "\n",
    "We assume $\\alpha = 0.5$ and $\\beta = 0.1$ as an example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# An example of Markov matrix construction\n",
    "alpha = 0.5\n",
    "beta = 0.1\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Construct P\n",
    "P = np.array([[1 - alpha, alpha], [beta, 1 - beta]])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Or, fill in the numbers\n",
    "P = np.zeros((2, 2))\n",
    "P[0, 0] = 1 - alpha\n",
    "P[0, 1] = alpha\n",
    "P[1, 0] = beta\n",
    "P[1, 1] = 1 - beta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.5, 0.5],\n",
       "       [0.1, 0.9]])"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chain Simulation, single agent\n",
    "\n",
    "How do we randomly pick a state that we are going to jump into?\n",
    "1. Generate a uniformly distributed random number U(0, 1)\n",
    "2. Suppose we are at state 0. If this number is greater than $1-\\alpha$ then we are going to jump to state 1. Otherwise, we stay at state 0.\n",
    "\n",
    "Then we iterate forward for T times."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# One step iteration\n",
    "\n",
    "cur_state = 0\n",
    "cum_prob = np.cumsum(P[cur_state, :])\n",
    "tmp = 0.4\n",
    "next_state = np.searchsorted(cum_prob, tmp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.5, 0.5])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P[0, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.5, 1. ])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cum_prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "next_state"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Iterate T times\n",
    "T = 1000\n",
    "\n",
    "ini_state = 0 # Starting from state 0 for example\n",
    "res = np.zeros(T, dtype = int) # Initialize the result record\n",
    "res[0] = ini_state\n",
    "prob_draw = np.random.rand(T) # Draw all the probabilities in a row\n",
    "for i in range(T - 1):\n",
    "    cum_prob = np.cumsum(P[res[i], :])\n",
    "    res[i + 1] = np.searchsorted(cum_prob, prob_draw[i])\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Wrap into a function\n",
    "\n",
    "def Markov_simulation_single(P, ini_state, T):\n",
    "    res = np.zeros(T, dtype = int) # Initialize the result record\n",
    "    res[0] = ini_state\n",
    "    prob_draw = np.random.rand(T) # Draw all the probabilities in a row\n",
    "    for i in range(T - 1):\n",
    "        cum_prob = np.cumsum(P[res[i], :])\n",
    "        res[i + 1] = np.searchsorted(cum_prob, prob_draw[i])\n",
    "        \n",
    "    return res\n",
    "\n",
    "simulation_res = Markov_simulation_single(P, 0, 1000)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.5, 0.5],\n",
       "       [0.1, 0.9]])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "# A bit more optimization\n",
    "def Markov_simulation_single_new(P, ini_state, T):\n",
    "    res = np.zeros(T, dtype = int) # Initialize the result record\n",
    "    res[0] = ini_state\n",
    "    prob_draw = np.random.rand(T) # Draw all the probabilities in a row\n",
    "    cum_prob = np.cumsum(P, axis = 1)\n",
    "    for i in range(T - 1):\n",
    "        res[i + 1] = np.searchsorted(cum_prob[res[i], :], prob_draw[i])\n",
    "        \n",
    "    return res\n",
    "\n",
    "simulation_res = Markov_simulation_single_new(P, 0, 1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "T = 100000\n",
    "simulation_res = Markov_simulation_single_new(P, 0, T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "sum_one = np.zeros(T)\n",
    "average_one = np.zeros(T)\n",
    "for i in range(T):\n",
    "    if (i > 0):\n",
    "        sum_one[i] = sum_one[i - 1] + simulation_res[i]\n",
    "    else:\n",
    "        sum_one[i] = simulation_res[i]\n",
    "    average_one[i] = sum_one[i] / (i + 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x216f9595e08>]"
      ]
     },
     "execution_count": 34,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXgAAAD4CAYAAADmWv3KAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAWJ0lEQVR4nO3dfZAkd33f8fe3Z2Z370mnk3RCj+EESLIvBRg4FEjCsw2ScKxKlasi+QFMICoSKzZ2pYywKwkpqmI7tlMuFdiyCstPMQjbECOTA+Fgx1SKWOhUBqFHdHpAOgtJK4TvTne3u/PwzR/dexotq7ne097tdd/7VbV107/+9cz3tzP3md/0dG9HZiJJap9irQuQJB0bBrwktZQBL0ktZcBLUksZ8JLUUt21euAzzjgjt23btlYPL0mNdPvttz+VmVvr9F2zgN+2bRu7du1aq4eXpEaKiG/V7esuGklqKQNeklrKgJekljLgJamlDHhJaikDXpJayoCXpJZqXMAPR8mf3PYog+ForUuRpBNa4wL+E7d+i1/49B38wf+rfay/JJ2UGhfw3947B8BtDz29xpVI0omtcQF/28NlsH/hrsfXuBJJOrE1LuCLiLUuQZIaoXEBb75LUj2NC/iR1wiXpFoaF/AY8JJUS+MCPscSPtO0l6Tn07iAH47to/nKA99Zw0ok6cTWuIAfjAX8j3/8VrZd+794vDo2XpL0rDW7ZN/RGi7zLevrfvlLXPOWl3HPt/fxpXuf5OPv2sHbvv9MwkNuJJ3EGhfwm9f1lm3/6F/vPnz7fX9YXuv1oV++3JCXdNJqXMCfun75gF/OBR/a+ZzlP3rvJbzhwloXI5ekxmtcwE86cObej1zKA7PP8M7r/u+y63/yd78KwIf/xXZ+44vfZP/84DnrH/ivl9MpTpwZ/6GFIaNMDiwMuOPRvRxYGPDwUwfpdYNn5gbM9DqctXmGg/MDXnrmRk7fMM23vnOA809bz7YzNrBxunFPr6RV1LgEWBrwt3zgjXx77yHedNFWIoJ/fM5mHv6Vd/Izn/w7bv76Y8vex4f/4u5l21/6izuXbV9OEdApgv6wLOhNF23l3C3r+A9vv5jTNkwBsG+uz3/+7F088vRBDswPuPfx/bXvfzX0OsH2czZzykyXmV6H9VMdukXBYDTisX84BEBEMN0tmOoUnH/aes7bso6F4YhMGI2SF22eITOZ6XXYON1l87oeo4SnD8wz1x9RFMHCYMSh/pDZfXOcsq7H/GBEtwjm+uWfdN66aZok6RUFm2a6bKjuZzAa0R8mm2bKl+Fcf8T6qQ5T3YIn9s1xcH7I5vU9BsNkYThi43SXLet7bJzpsmm6R6cIprrF4cNlV7o7bq4/ZDhKRpkcXBjSLYK9h/rsnxswPxjRKaA/LNc/U7X1OkVV94gigu8eWGB+MGKuP+Jgf8C+QwMOzA9YGIwOH9K7WH8RwZb1PdZNdQ/X3O0E66e6RMB0p2D9dJdOBAvD0eH65gcjhqPyPhYGI2Z6BcNR+RpcGIw42B+yabrLppkuU92CXqf8GY6S/nDEdLdgptdhptcBylNJhsMR070OB+YHDEbJKTM9ep1gfjBilEl/mBRRfufVKYJ1vQ7TvYIiglEmnaKgV/3+F3+Hi+MYjpLhKJnrjxhmHn79RPUcFQHdTkG3CAajZL4/ZKbXodsJRiPoj0b0q9fUTK98PfSKgqJ6PnpF0O0UzPTKY0QWfz+H+kP6gxEJTHULOhEMRiMGw2TdVIcgiHg2Q0aZDDMZjZJRVsujpIh4trbBkMEwKQqY6nQObz/dLeh0ApLDz1V/WP6+xyeJmcmg+n1MdQqGWf5eXnHeZl677bQV/59eqcYF/LiXbN3AxWdt4uKzNn3PuuuuehXXXfWq57Ttn+vz8g9/8Tlt55+2jkefPrTixx4ljIbPvtv8zTdnAfjErY8ccdtNM132zw0m9jnrlBke3zfHj7zyHLadvp6Lztp0+Avmc05dx8NPHWDL+ikO9ofM7p9n36E+U92CU9f3+IeDfR566gD3Pr6P/XN97n5sLxumu/Q6BbP75w8/xkUv2sip66bYe6jPwmDEl++fPfyG1RRTVeCOsrzd7QS9TsEok16nDKQiOBwkvSoA5/tD9h3hOTiaWk5ZV76BTXcLgvI/eqcIetUb0QOzz5RvjFUG9IcjDi4MySooxkWUf3tpuluG4VT1Rjw3KLdfHPP6qQ775vocmB8yPxh6tncD/Ns3v9SAX87irOj7ztrEFz7wxhVtu2mmx8O/8k5G1WxocUazaDRKPn/n49z3xH6u+9L9/O+ffyMvO3MTmUlEMLt/nvuf2M/eQ32+fP9T/KvXns/LztzIz33qa/zl3U8s+5gvPn09P/PWC/nB7S/ilJnuqn3peyxeHMNR8sS+ucOz/eEo2fPdQ2yY7nBoYcjeQ332HurT6xSctmGKdVMdBsNkqlsw3S3Yummag/PDclaXSbdKsaeemS8/7QySfXPlDHnfXJ/pbhnABxeGQDLd6/DM3IDhKDlj4zQbZ7rV4wVTnYL9cwNmn5lnrj9k36E+o4QDCwO6RVBE+UliUM1aiwj6wzL4F5+/6W5Bf1jO9qa6BadvmD48C9w406U/GHHKuh6bZnrM9ApGCb0iIGDTdI/pXrl9tyjodcqZ7Knrp8rZbbeg23lhRx1nJgeq3XKLn6qO5vWyOHPvFOVMtJxhlrPMAAjoFgVz/SEbprp0OsH+uf7h57JTBL2ifJMsimBUzY4XZ8qdIhhWn74WBuXvuijgwPzw8JtStwime0X1Ggjm+kOgnCWPEobVzLrbCaa7HQ4uDBmMyt9tt1M+/vrqdbf4uMNR0uvE4RnxwYVyPNO9crY+M9VhqvPsrD4zD39SOFR9Glr8FAHlJ6CiCDoRh8ew+AllcZwzvQ7dIhglzA+Gh9+05wdDBtX99TrlJ6TyNRAMxiZJnSIO/yy+Lme65aeS4yHW6mzQHTt25K5du1a83S/+z2/wiVsf4S0Xb+X33nPJMahMkk5cEXF7Zu6o07dxJzr98MvPBuDfvOEla1yJJJ3YagV8RFwaEfdFxO6IuHaZ9Zsj4i8i4usRcVdEvGf1S118sPKf4gQ62kWSTkRHDPiI6AAfAy4DtgNXRcT2Jd1+Grg7M18JvBn4jYiYWuVaJUkrUGcGfwmwOzMfzMwF4CbgiiV9EtgU5bcXG4GngdU9REGStCJ1Av5c4NGx5T1V27iPAt8PPAZ8A/jZzBwt6UNEXB0RuyJi1+zs7FGWLEmqo07AL7eze+mhN+8AvgacA/wA8NGIOOV7Nsq8ITN3ZOaOrVv9kwGSdCzVCfg9wPljy+dRztTHvQf4TJZ2Aw8B37c6JUqSjkadgL8NuDAiLqi+OL0SuHlJn0eAtwFExIuAi4EHV7NQSdLKHPFM1swcRMQ1wC1AB7gxM++KiPdX668HPgL8fkR8g3KXzgcz86ljUrGnYUtSLbX+VEFm7gR2Lmm7fuz2Y8DbV7e0yTwKXpIma9yZrJKkegx4SWopA16SWsqAl6SWMuAlqaUMeElqqcYFvIfBS1I9jQv4Rat16TtJaqvGBrwkaTIDXpJayoCXpJYy4CWppQx4SWopA16SWsqAl6SWalzAp2c6SVItjQv4RZ7nJEmTNTbgJUmTGfCS1FIGvCS1lAEvSS1lwEtSSxnwktRSBrwktVTjAj69ppMk1dK4gF/keU6SNFljA16SNJkBL0ktZcBLUksZ8JLUUga8JLWUAS9JLdW4gPeCH5JUT+MCfpEX/JCkyRob8JKkyWoFfERcGhH3RcTuiLj2efq8OSK+FhF3RcTfrG6ZkqSV6h6pQ0R0gI8BPwTsAW6LiJsz8+6xPqcCvwVcmpmPRMSZx6pgSVI9dWbwlwC7M/PBzFwAbgKuWNLnx4DPZOYjAJn55OqWKUlaqToBfy7w6Njynqpt3EXAloj4PxFxe0S8a7k7ioirI2JXROyanZ09uoolSbXUCfjljldZerBiF3gN8E7gHcB/jIiLvmejzBsyc0dm7ti6deuKi5Uk1XfEffCUM/bzx5bPAx5bps9TmXkAOBARXwZeCXxzVaqUJK1YnRn8bcCFEXFBREwBVwI3L+nzWeANEdGNiPXAPwHuWd1SS57nJEn1HHEGn5mDiLgGuAXoADdm5l0R8f5q/fWZeU9EfAG4AxgBH8/MO49l4V7yQ5Imq7OLhszcCexc0nb9kuVfA35t9UqTJL0QnskqSS1lwEtSSxnwktRSBrwktZQBL0kt1biAT6/4IUm1NC7gF3nBD0marLEBL0mazICXpJYy4CWppQx4SWopA16SWsqAl6SWMuAlqaUaF/Ce5iRJ9TQu4Bd5npMkTdbYgJckTWbAS1JLGfCS1FIGvCS1lAEvSS1lwEtSSxnwktRSzQt4z3SSpFqaF/CV8JJOkjRRYwNekjSZAS9JLWXAS1JLGfCS1FIGvCS1lAEvSS3VuIBPD4SXpFoaF/CLPApekiZrbMBLkiYz4CWppWoFfERcGhH3RcTuiLh2Qr/XRsQwIn509UqUJB2NIwZ8RHSAjwGXAduBqyJi+/P0+1XgltUuUpK0cnVm8JcAuzPzwcxcAG4Crlim378HPg08uYr1SZKOUp2APxd4dGx5T9V2WEScC/xL4PpJdxQRV0fErojYNTs7u9JaJUkrUCfglzsicenB6L8JfDAzh5PuKDNvyMwdmblj69atdWuUJB2Fbo0+e4Dzx5bPAx5b0mcHcFP1N9rPAC6PiEFm/vmqVDkmPc9JkmqpE/C3ARdGxAXA3wNXAj823iEzL1i8HRG/D3zuWIT7OK/3IUmTHTHgM3MQEddQHh3TAW7MzLsi4v3V+on73SVJa6PODJ7M3AnsXNK2bLBn5k+98LIkSS+UZ7JKUksZ8JLUUga8JLWUAS9JLWXAS1JLNS7gPdFJkuppXMAvCq/pJEkTNTbgJUmTGfCS1FIGvCS1lAEvSS1lwEtSSxnwktRSjQt4D4OXpHoaF/CLvOCHJE3W2ICXJE1mwEtSSxnwktRSBrwktZQBL0ktZcBLUksZ8JLUUo0L+PSKH5JUS+MCXpJUjwEvSS1lwEtSSxnwktRSBrwktZQBL0ktZcBLUks1LuA9Cl6S6mlcwC/ygh+SNFljA16SNJkBL0ktZcBLUkvVCviIuDQi7ouI3RFx7TLrfzwi7qh+vhIRr1z9UiVJK3HEgI+IDvAx4DJgO3BVRGxf0u0h4E2Z+QrgI8ANq12oJGll6szgLwF2Z+aDmbkA3ARcMd4hM7+Smd+tFv8WOG91y5QkrVSdgD8XeHRseU/V9nzeC3x+uRURcXVE7IqIXbOzs/WrlCStWJ2AX+6I82XPN4qIt1AG/AeXW5+ZN2TmjszcsXXr1vpVPuc+jmozSTrpdGv02QOcP7Z8HvDY0k4R8Qrg48Blmfmd1Snv+cWy7zuSpEV1ZvC3ARdGxAURMQVcCdw83iEi/hHwGeAnM/Obq1+mJGmljjiDz8xBRFwD3AJ0gBsz866IeH+1/nrgPwGnA78V5d8QGGTmjmNXtiTpSOrsoiEzdwI7l7RdP3b7fcD7Vrc0SdIL4ZmsktRSBrwktZQBL0ktZcBLUks1MOA900mS6mhgwJe8opMkTdbYgJckTWbAS1JLGfCS1FIGvCS1lAEvSS1lwEtSSzUu4L3ghyTV07iAX+Rx8JI0WWMDXpI0mQEvSS1lwEtSSxnwktRSBrwktZQBL0ktZcBLUks1LuA9z0mS6mlcwC8KPNNJkiZpbMBLkiYz4CWppQx4SWopA16SWsqAl6SWMuAlqaUMeElqqcYFvFd0kqR6Ghfwi7yikyRN1tiAlyRNZsBLUksZ8JLUUrUCPiIujYj7ImJ3RFy7zPqIiOuq9XdExKtXv1RJ0kocMeAjogN8DLgM2A5cFRHbl3S7DLiw+rka+O1VrlOStEJ1ZvCXALsz88HMXABuAq5Y0ucK4A+z9LfAqRFx9irXKklagToBfy7w6NjynqptpX2IiKsjYldE7JqdnV1prQCctXmGy19+Fhunu0e1vSSdLOqk5HJHnC893ahOHzLzBuAGgB07dhzVKUuvefEWXvPi1xzNppJ0Uqkzg98DnD+2fB7w2FH0kSQdR3UC/jbgwoi4ICKmgCuBm5f0uRl4V3U0zeuAvZn57VWuVZK0AkfcRZOZg4i4BrgF6AA3ZuZdEfH+av31wE7gcmA3cBB4z7ErWZJUR61vKjNzJ2WIj7ddP3Y7gZ9e3dIkSS+EZ7JKUksZ8JLUUga8JLWUAS9JLRW5RpdIiohZ4FtHufkZwFOrWE4TOOaTg2M+ObyQMb84M7fW6bhmAf9CRMSuzNyx1nUcT4755OCYTw7Ha8zuopGkljLgJamlmhrwN6x1AWvAMZ8cHPPJ4biMuZH74CVJR9bUGbwk6QgMeElqqcYF/JEuAH4ii4jzI+KvI+KeiLgrIn62aj8tIv4yIu6v/t0yts2HqrHeFxHvGGt/TUR8o1p3XURE1T4dEZ+q2m+NiG3He5zLiYhORPxdRHyuWm71mCPi1Ij4s4i4t3q+X38SjPnnqtf1nRHxyYiYaduYI+LGiHgyIu4cazsuY4yId1ePcX9EvLtWwZnZmB/KP1f8APASYAr4OrB9retaQf1nA6+ubm8Cvkl5IfP/BlxbtV8L/Gp1e3s1xmnggmrsnWrdV4HXU15N6/PAZVX7vwOur25fCXxqrcdd1fLzwCeAz1XLrR4z8AfA+6rbU8CpbR4z5SU6HwLWVct/AvxU28YMvBF4NXDnWNsxHyNwGvBg9e+W6vaWI9a71v8RVvjLfT1wy9jyh4APrXVdL2A8nwV+CLgPOLtqOxu4b7nxUf5N/tdXfe4da78K+J3xPtXtLuXZcrHG4zwP+BLwVp4N+NaOGTiFMuxiSXubx7x4XebTqno+B7y9jWMGtvHcgD/mYxzvU637HeCqI9XatF00tS7u3QTVR69XAbcCL8rqCljVv2dW3Z5vvOdWt5e2P2ebzBwAe4HTj8UYVuA3gV8ARmNtbR7zS4BZ4Peq3VIfj4gNtHjMmfn3wK8DjwDfpryq2xdp8ZjHHI8xHlX2NS3ga13c+0QXERuBTwMfyMx9k7ou05YT2idtsyYi4oeBJzPz9rqbLNPWqDFTzrxeDfx2Zr4KOED50f35NH7M1X7nKyh3RZwDbIiIn5i0yTJtjRpzDas5xqMae9MCvvEX946IHmW4/3FmfqZqfiIizq7Wnw08WbU/33j3VLeXtj9nm4joApuBp1d/JLX9M+BHIuJh4CbgrRHxP2j3mPcAezLz1mr5zygDv81j/kHgocyczcw+8Bngn9LuMS86HmM8quxrWsDXuQD4Cav6pvx3gXsy87+PrboZWPxW/N2U++YX26+svlm/ALgQ+Gr1MXB/RLyuus93Ldlm8b5+FPirrHbarYXM/FBmnpeZ2yifr7/KzJ+g3WN+HHg0Ii6umt4G3E2Lx0y5a+Z1EbG+qvVtwD20e8yLjscYbwHeHhFbqk9Lb6/aJjveX1Cswhccl1MeffIA8EtrXc8Ka//nlB+r7gC+Vv1cTrmP7UvA/dW/p41t80vVWO+j+qa9at8B3Fmt+yjPnpU8A/wp5QXQvwq8ZK3HPVbzm3n2S9ZWjxn4AWBX9Vz/OeWRD20f838B7q3q/SPKo0daNWbgk5TfMfQpZ9XvPV5jBP511b4beE+dev1TBZLUUk3bRSNJqsmAl6SWMuAlqaUMeElqKQNeklrKgJekljLgJaml/j9HIunNwXhGWwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(average_one)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Chain simulation, distribution\n",
    "\n",
    "The distribution is characterized by a vector. Then the simulation is just vector-matrix multiplication"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "def Markov_simulation_dist(P, ini_dist, T):\n",
    "    n = len(ini_dist)\n",
    "    res = np.zeros((T, n)) # initialize the result\n",
    "    res[0, :] = ini_dist # starting point\n",
    "    for i in range(T - 1):\n",
    "        res[i + 1, :] = res[i, :] @ P\n",
    "        \n",
    "    return res"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 116,
   "metadata": {},
   "outputs": [],
   "source": [
    "ini_dist = np.array([0.7, 0.3])\n",
    "T = 1000\n",
    "res = Markov_simulation_dist(P, ini_dist, T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 117,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.16666667, 0.83333333])"
      ]
     },
     "execution_count": 117,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "res[T - 1, :]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 113,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAlMAAAFlCAYAAADPim3FAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3df5TV9X3n8ed7ht8/R2SCjaDSrEbEYo2sgZikarb+ql30bDcFTWuybSg9aqTbs+uvs5pNz5p6uvS0OVHRpJS62+hJGpfaKrU9zUbXxhBIYhIJ6rJEwkA0+AOBtAgD7/3jXsZxMsidmfudSz7zfJyD997v9zP3vu98jvjy8/l8P9/ITCRJkjQ4ba0uQJIk6WeZYUqSJGkIDFOSJElDYJiSJEkaAsOUJEnSEBimJEmShmBUqz54+vTpecopp7Tq4yVJkhr2zW9+8+XM7OzvXMvC1CmnnMKGDRta9fGSJEkNi4itRzrnNJ8kSdIQGKYkSZKGwDAlSZI0BIYpSZKkITBMSZIkDYFhSpIkaQgMU5IkSUNgmJIkSRoCw5QkSdIQGKYkSZKGwDAlSZI0BC27N19L7HkRXt0C3W9AHgISMltdlSRJGoqpM+Edc1r28SMnTP3DbfBPf9rqKiRJUrPN/y24/I9b9vEjI0xt/1YtSM1bDGf9OowaD9EGEUDUHyVJ0s+kidNb+vEjI0w982VoHwOX3gnjO1pdjSRJKsjIWID+o+/Az51lkJIkSU03MsLUqz+AaT/f6iokSVKByg9TB/bB7u2GKUmSVInyw9TrXUBCx8mtrkSSJBWo/DC17/Xa44Rpra1DkiQVqfww9UY9TI2d3No6JElSkUZAmNpTexw7pbV1SJKkIo2gMOXIlCRJar7yw9S+3bVHw5QkSapAQ2EqIi6JiOciYnNE3NTP+akR8TcR8Z2I2BgRH2t+qYPkNJ8kSarQUcNURLQDdwGXAmcASyLijD7NrgW+n5lnAecDKyJiTJNrHZw3dsPoCdA+Mu6cI0mShlcjI1PnApszc0tm7gceBBb1aZPA5IgIYBLwKtDd1EoH643dTvFJkqTKNBKmTgS29XrdVT/W22eBOcAO4HvADZl5qO8bRcTSiNgQERt27tw5yJIH6I09hilJklSZRsJU9HMs+7y+GHgaeCfwi8BnI+KnFill5n2ZOT8z53d2dg642EHp3g+jxg3PZ0mSpBGnkTDVBczq9XomtRGo3j4GPJQ1m4EfAKc3p8QhOnQA2ke3ugpJklSoRsLUeuDUiJhdX1S+GHi4T5sfAh8CiIgZwLuBLc0sdNAO7oc2w5QkSarGUS9xy8zuiLgOeAxoB1Zl5saIWFY/vxL4A2B1RHyP2rTgjZn5coV1N+7gAWg/Ni4slCRJ5Wlov4DMfBR4tM+xlb2e7wAuam5pTXJwP4yZ1OoqJElSocrfAd2RKUmSVKEREqZcMyVJkqoxAsLUfsOUJEmqTPlh6pDTfJIkqTrlh6mDB9waQZIkVWYEhCmn+SRJUnVGQJhymk+SJFVnhIQpR6YkSVI1RkCYcppPkiRVp+wwlenVfJIkqVJlh6lD3bVHR6YkSVJFyg5TB/fXHt0aQZIkVWRkhCmn+SRJUkUKD1NO80mSpGoVHqYcmZIkSdUaIWHKkSlJklSNssNUz9V8jkxJkqRqlB2meq7mG9XaOiRJUrHKDlN5qPbY1t7aOiRJUrFGRpiKsr+mJElqnbJTxqGDtUfDlCRJqkjZKSOz9hhO80mSpGoUHqac5pMkSdUqO2Xk4Wm+aG0dkiSpWIWHKa/mkyRJ1RoZYcppPkmSVJGyU4ZhSpIkVazslOHWCJIkqWJlp4yekSnXTEmSpGoUHqYO7zNV9teUJEmtU3bKSKf5JElStcpOGT3TfO4zJUmSqjEywpT7TEmSpIqMjDDlNJ8kSapI2SnDrREkSVLFGkoZEXFJRDwXEZsj4qZ+zv+niHi6/ueZiDgYEdOaX+4AuTWCJEmq2FHDVES0A3cBlwJnAEsi4ozebTLzjzLzFzPzF4Gbgccz89UqCh4Qp/kkSVLFGkkZ5wKbM3NLZu4HHgQWvU37JcADzShuyAxTkiSpYo2kjBOBbb1ed9WP/ZSImABcAnx56KUNzcYdr/P7X/x27YVbI0iSpIo0Eqb6SyJ5hLa/CvzTkab4ImJpRGyIiA07d+5stMbBO+TWCJIkqVqNhKkuYFav1zOBHUdou5i3meLLzPsyc35mzu/s7Gy8ykEIgrZwmk+SJFWrkZSxHjg1ImZHxBhqgenhvo0iYirwS8BfN7fEwYmANrw3nyRJqtaoozXIzO6IuA54DGgHVmXmxohYVj+/st70SuDvM/MnlVU7ALUw5dYIkiSpWkcNUwCZ+SjwaJ9jK/u8Xg2sblZhQxWEI1OSJKlyxaaMt45MFfs1JUlSixWbMgLXTEmSpOoVmzLeMjLVVuzXlCRJLVZwynDNlCRJql6xKcM1U5IkaTgUnTLeHJlyawRJklSNYsOUC9AlSdJwKDZlRAThNJ8kSapYsSkjgHbDlCRJqlixKeMt9+Zrc82UJEmqRrlhiqAtDq+ZitYWI0mSilVumKpvjZDlfkVJknQMKDpptHGIQ66XkiRJFSo2afSsmTJMSZKkChWbNGpbI6TTfJIkqVLFJo3DWyOki88lSVKFig1TUJvmS28lI0mSKlRsmHrzRseOTEmSpOqUG6aI2tYILkCXJEkVKjZpHL6azzAlSZKqVGzSCOphqtyvKEmSjgHlJo2AcJpPkiRVrNikEQTtJOkCdEmSVKFyw9The/O5NYIkSapQuWEKaAtHpiRJUrXKDVPh1giSJKl6xSaNnqv5DFOSJKlCRSeNNg65NYIkSapUsUkjAsKRKUmSVLFik0ZtawTvzSdJkqpVbJjC28lIkqRhUGzS6Jnma3UhkiSpaOWGKQ5P8DnNJ0mSqlNumIoAbycjSZIqVm6Y4vA0n2FKkiRVp9wwFU7wSZKk6jUUpiLikoh4LiI2R8RNR2hzfkQ8HREbI+Lx5pY5cEE4MiVJkio36mgNIqIduAv4ZaALWB8RD2fm93u16QDuBi7JzB9GxDuqKrhREX2fSJIkNV8jI1PnApszc0tm7gceBBb1aXMV8FBm/hAgM3/c3DIHJ9wYQZIkVayRMHUisK3X6676sd5OA46LiK9GxDcj4jf7e6OIWBoRGyJiw86dOwdX8QDUrudzZEqSJFWnkTDVXxrpO+QzCjgH+BXgYuC/RMRpP/VDmfdl5vzMnN/Z2TngYgfi8KadkiRJVTrqmilqI1Gzer2eCezop83LmfkT4CcR8QRwFvB8U6ocBBegS5Kk4dDIyNR64NSImB0RY4DFwMN92vw18IGIGBURE4D3ApuaW+rAHF53bpiSJElVOurIVGZ2R8R1wGNAO7AqMzdGxLL6+ZWZuSki/g74LnAI+HxmPlNl4UdzeNPOgrfSkiRJx4BGpvnIzEeBR/scW9nn9R8Bf9S80oYmIrzRsSRJqlyxwzaHb3Sc7jMlSZIqVG6Y6rmazzAlSZKqU3CYCveZkiRJlSs2TEmSJA2HosNUhMvPJUlStcoOU27aKUmSKlZ4mHLNlCRJqlbhYcp9piRJUrUKD1Nv/lOSJKkKhYcp10xJkqRqGaYkSZKGoOwwZY6SJEkVKztMOTIlSZIqVn6YcnhKkiRVqOgw5ZV8kiSpakWHqXCXKUmSVLHiw1Smo1OSJKk6hYcpXDMlSZIqNarVBVTJ28lIkkp14MABurq62LdvX6tLKca4ceOYOXMmo0ePHtDPFR6mvNGxJKlMXV1dTJ48mVNOOYVwFmbIMpNXXnmFrq4uZs+ePaCfLXuaLxKv6JMklWjfvn0cf/zxBqkmiQiOP/74QY30lR2mnOaTJBXMINVcg/19Fh6mnOaTJEnVKjxMeTsZSZJUreLDlCRJGl5PPvkk73vf+5g6dSrTpk3jvPPOY/369QCsXr2a97///Q2/1wsvvEBE0N3dPeh6nn76ac455xwmTJjAOeecw9NPPz3o9+pP0WGqFqccmZIkabjs3r2byy+/nOuvv55XX32V7du3c/vttzN27NiW1LN//34WLVrERz7yEV577TWuueYaFi1axP79+5v2GUWHKUemJEkaXs8//zwAS5Ysob29nfHjx3PRRRcxb948Nm3axLJly3jqqaeYNGkSHR0dADzyyCOcffbZTJkyhVmzZvHJT36y5/0++MEPAtDR0cGkSZN46qmnAFi1ahVz5szhuOOO4+KLL2br1q391vPVr36V7u5uli9fztixY/nEJz5BZvKVr3ylad+5+DDlyJQkScPntNNOo729nWuuuYa1a9fy2muv9ZybM2cOK1euZOHChezdu5ddu3YBMHHiRO6//3527drFI488wj333MOaNWsAeOKJJwDYtWsXe/fuZeHChaxZs4Y77riDhx56iJ07d/KBD3yAJUuW9FvPxo0bmTdv3luu1Js3bx4bN25s2nceAZt2SpJUtv/6Nxv5/o7dlX7GGe+cwu2/Oveo7aZMmcKTTz7JnXfeycc//nFefPFFLrvsMj73uc8xY8aMfn/m/PPP73k+b948lixZwuOPP84VV1zRb/t7772Xm2++mTlz5gBwyy23cMcdd7B161ZOPvnkt7Tdu3cvU6dOfcuxqVOnsmfPnqN+l0Y5MiVJkppqzpw5rF69mq6uLp555hl27NjB8uXLj9h+3bp1XHDBBXR2djJ16lRWrlzJyy+/fMT2W7du5YYbbqCjo4OOjg6mTZtGZrJ9+/afajtp0iR2735r0Ny9ezeTJ08e/Bfso/CRKXdAlySVr5ERo1Y5/fTT+ehHP8q9994L9L8x5lVXXcV1113H2rVrGTduHMuXL+8JU/21nzVrFrfeeitXX331UT9/7ty5rFixgszsea/vfve7XHvttUP5Wm9R+MiU03ySJA2nZ599lhUrVtDV1QXAtm3beOCBB1iwYAEAM2bMoKur6y1X0+3Zs4dp06Yxbtw4vvGNb/CFL3yh51xnZydtbW1s2bKl59iyZcv49Kc/3bPu6fXXX+dLX/pSv/Wcf/75tLe385nPfIY33niDz372swBceOGFTfvOhYcpp/kkSRpOkydPZt26dbz3ve9l4sSJLFiwgDPPPJMVK1YAtRAzd+5cTjjhBKZPnw7A3XffzW233cbkyZP51Kc+xYc//OGe95swYQK33nor5513Hh0dHXz961/nyiuv5MYbb2Tx4sVMmTKFM888k7Vr1/Zbz5gxY1izZg33338/HR0drFq1ijVr1jBmzJimfefIbM3Yzfz583PDhg2Vfsa2T57Oy1PncvbvfbnSz5Ekabht2rSpZwG2mudIv9eI+GZmzu/vZ0bAyJQkSVJ1GgpTEXFJRDwXEZsj4qZ+zp8fEa9HxNP1P7c1v9SBc5pPkiRV7ahX80VEO3AX8MtAF7A+Ih7OzO/3afp/MvPyCmocvOj5hyRJUiUaGZk6F9icmVsycz/wILCo2rKax2k+SZJUpUbC1InAtl6vu+rH+loYEd+JiLURcUxseOGaKUmSVLVGNu3sb56sb0b5FnByZu6NiMuANcCpP/VGEUuBpQAnnXTSAEsdODftlCRJVWtkZKoLmNXr9UxgR+8Gmbk7M/fWnz8KjI6I6X3fKDPvy8z5mTm/s7NzCGU3prZpp2FKkiRVp5EwtR44NSJmR8QYYDHwcO8GEXFC1Pdoj4hz6+/7SrOLHSin+SRJUtWOOs2Xmd0RcR3wGNAOrMrMjRGxrH5+JfBrwO9GRDfwL8DibNVuoL24NYIkSapaQ/tMZeajmXlaZr4rM/9b/djKepAiMz+bmXMz86zMXJCZX6uy6EY5zSdJ0vB78skned/73sfUqVOZNm0a5513HuvXrwdg9erVvP/972/4vV544QUigu7u7kHXs3TpUt797nfT1tbG6tWrB/0+R1L0Dug4zSdJ0rDavXs3l19+Oddffz2vvvoq27dv5/bbb2fs2LEtq+mss87i7rvv5j3veU8l7190mAqSTEemJEkaLs8//zwAS5Ysob29nfHjx3PRRRcxb948Nm3axLJly3jqqaeYNGkSHR0dADzyyCOcffbZTJkyhVmzZvHJT36y5/0++MEPAtDR0cGkSZN46qmnAFi1ahVz5szhuOOO4+KLL2br1q1HrOnaa6/lQx/6EOPGjavkOxcepsAsJUnS8DnttNNob2/nmmuuYe3atbz22ms95+bMmcPKlStZuHAhe/fuZdeuXQBMnDiR+++/n127dvHII49wzz33sGbNGgCeeOIJAHbt2sXevXtZuHAha9as4Y477uChhx5i586dfOADH2DJkiXD/2XrGtln6mecaUqSVLi1N8GL36v2M074Bbj0D4/abMqUKTz55JPceeedfPzjH+fFF1/ksssu43Of+xwzZszo92fOP//8nufz5s1jyZIlPP7441xxxRX9tr/33nu5+eabmTNnDgC33HILd9xxB1u3buXkk08e+HcbosJHptL7yUiSNMzmzJnD6tWr6erq4plnnmHHjh0sX778iO3XrVvHBRdcQGdnJ1OnTmXlypW8/PLLR2y/detWbrjhBjo6Oujo6GDatGlkJtu3b6/i6xxV0SNTbo0gSRoRGhgxapXTTz+dj370o9x7770A1LelfIurrrqK6667jrVr1zJu3DiWL1/eE6b6az9r1ixuvfVWrr766mqLb1DRI1PgwJQkScPp2WefZcWKFXR1dQGwbds2HnjgARYsWADAjBkz6OrqYv/+/T0/s2fPHqZNm8a4ceP4xje+wRe+8IWec52dnbS1tbFly5aeY8uWLePTn/40GzduBOD111/nS1/60hFr2r9/P/v27SMzOXDgAPv27ePQoUNN+85FhylHpiRJGl6TJ09m3bp1vPe972XixIksWLCAM888kxUrVgBw4YUXMnfuXE444QSmT6/dee7uu+/mtttuY/LkyXzqU5/iwx/+cM/7TZgwgVtvvZXzzjuPjo4Ovv71r3PllVdy4403snjxYqZMmcKZZ57J2rVrj1jTRRddxPjx4/na177G0qVLGT9+fM/C9maIVm1UPn/+/NywYUOln/H6f53JtydfwPn/8X9U+jmSJA23TZs29SzAVvMc6fcaEd/MzPn9/UzhI1PugC5JkqpVeJhymk+SJFWr6DCFYUqSJFWs6DBVm+aTJEmqTuFhynvzSZLK1aqLyEo12N9n0WEKvDefJKlM7e3tHDhwoNVlFOXAgQOMGjXw/cyLDlPhJJ8kqVAdHR289NJLTd18ciQ7dOgQL730ElOnTh3wzxZ9OxkSp/kkSUWaPn06XV1dPPfcc60upRgTJ07s2Uh0IIoOU7WtESRJKk9bWxsnnXRSq8sQI2Caz60RJElSlQoPU+6ALkmSqlV0mMJpPkmSVLGiw5TTfJIkqWqGKUmSpCEoOkyBt5ORJEnVKjpM1W4n0+oqJElSyYoOU+DVfJIkqVpFhyk37ZQkSVUrPEy5ZkqSJFWr6DCFV/NJkqSKFR2m2gxTkiSpYkWHKYBMw5QkSapOuWGqvieCa6YkSVKVRkCYcmRKkiRVp9wwhSNTkiSpeuWGKaf5JEnSMCg3TNVj1CGn+SRJUoUaClMRcUlEPBcRmyPiprdp968j4mBE/FrzSpQkSTp2HTVMRUQ7cBdwKXAGsCQizjhCuzuBx5pd5KAcnuZzawRJklShRkamzgU2Z+aWzNwPPAgs6qfd9cCXgR83sb4h8Go+SZJUvUbC1InAtl6vu+rHekTEicCVwMq3e6OIWBoRGyJiw86dOwda68C4AF2SJA2DRsJUf0M7fTPKnwA3ZubBt3ujzLwvM+dn5vzOzs5Gaxwkw5QkSareqAbadAGzer2eCezo02Y+8GBEAEwHLouI7sxc05QqB8NNOyVJ0jBoJEytB06NiNnAdmAxcFXvBpk5+/DziFgN/G1LgxTgyJQkSRoORw1TmdkdEddRu0qvHViVmRsjYln9/Nuuk2qZ+sjUIa/mkyRJFWpkZIrMfBR4tM+xfkNUZn506GU1g9N8kiSpeuXugO7VfJIkaRiUG6bq3LRTkiRVqeAw5ZiUJEmqXrlh6vAC9BaXIUmSylZumOLw1XwtLkOSJBWt3DDlpp2SJGkYlBum6hyZkiRJVSo3TLlmSpIkDYNywxTugC5JkqpXbphKF6BLkqTqlRumcJpPkiRVr9ww5ciUJEkaBuWGqTrXTEmSpCoVHKayz6MkSVLzlRum6tN8Bx2ZkiRJFSo3THk7GUmSNAzKDVNu2ilJkoZBuWHq8MiUaUqSJFWo3DDVMzLlmilJklSdcsOUa6YkSdIwKDdMpffmkyRJ1Ss3THk7GUmSNAzKDVM9+0y1uA5JklS0csNUXSZkmqgkSVI1Cg5TbwYoF6FLkqSqlBum6qNRSXDIkSlJklSRcsNUnWFKkiRVqdww1TMy1fNUkiSp6coNU7w5zXfQRVOSJKki5YapnuEop/kkSVJ1yg1TvDnN58CUJEmqSrlhqvfVfKYpSZJUkXLDFG6NIEmSqldumEqn+SRJUvXKDVN1SXg7GUmSVJmGwlREXBIRz0XE5oi4qZ/ziyLiuxHxdERsiIj3N7/UgXozQB00TEmSpIqMOlqDiGgH7gJ+GegC1kfEw5n5/V7N/hF4ODMzIuYBXwROr6Lghr3ldjItrUSSJBWskZGpc4HNmbklM/cDDwKLejfIzL355lzaRHoPC7VMrzVTpilJklSRRsLUicC2Xq+76sfeIiKujIhngUeA/9Cc8obATTslSdIwaCRMRT/HfiqdZOb/yszTgSuAP+j3jSKW1tdUbdi5c+fAKh2oKe/kO79wC9/Pk5zmkyRJlWkkTHUBs3q9ngnsOFLjzHwCeFdETO/n3H2ZOT8z53d2dg642AGZOJ0X3vURtuUMR6YkSVJlGglT64FTI2J2RIwBFgMP924QEf8qIqL+/D3AGOCVZhc7UG21ktwaQZIkVeaoV/NlZndEXAc8BrQDqzJzY0Qsq59fCfw74Dcj4gDwL8Cv5zGQYA6HqYOHWlyIJEkq1lHDFEBmPgo82ufYyl7P7wTubG5pQ9deH3dzmk+SJFWl6B3Q6zOPhilJklSZosPU4Wm+Q07zSZKkihQepmqPjkxJkqSqlB2m2pzmkyRJ1So7TPWsmWpxIZIkqViFh6naoyNTkiSpKkWHqfaeBeiGKUmSVI2iw1Q4zSdJkipWdJg6PM13DGzGLkmSClV2mKqnqYOGKUmSVJGyw5TTfJIkqWKFh6nao1fzSZKkqhQepryaT5IkVavoMNXe5jSfJEmqVtFhKpzmkyRJFSs6TB2e5nNrBEmSVJWiw9So+jTfgYOGKUmSVI2iw9To9trX6z50qMWVSJKkUhUdpka110emuh2ZkiRJ1Sg6TI2pj0ztP+jIlCRJqkbRYapnms8wJUmSKlJ2mBpV+3ouQJckSVUpOkwdvprPaT5JklSVosPUm9N8jkxJkqRqFB2m2tuC9rbggCNTkiSpIkWHKahN9RmmJElSVYoPU2Pa21wzJUmSKlN8mBo9qs01U5IkqTLFhymn+SRJUpWKD1OjneaTJEkVKj5MjXGaT5IkVaj4MDW63Wk+SZJUneLD1Ki2NsOUJEmqTPFhavSoNvY7zSdJkipSfJga0x50OzIlSZIqUnyYcppPkiRVqaEwFRGXRMRzEbE5Im7q5/zVEfHd+p+vRcRZzS91cJzmkyRJVTpqmIqIduAu4FLgDGBJRJzRp9kPgF/KzHnAHwD3NbvQwXKaT5IkVamRkalzgc2ZuSUz9wMPAot6N8jMr2Xma/WXXwdmNrfMwRvd3sb+bsOUJEmqRiNh6kRgW6/XXfVjR/JbwNr+TkTE0ojYEBEbdu7c2XiVQzB+TDv/vP/gsHyWJEkaeRoJU9HPsX4XIUXEBdTC1I39nc/M+zJzfmbO7+zsbLzKIZg8dhR79h0Yls+SJEkjz6gG2nQBs3q9ngns6NsoIuYBnwcuzcxXmlPe0E0eN5q9b3STmUT0lwslSZIGr5GRqfXAqRExOyLGAIuBh3s3iIiTgIeA38jM55tf5uBNGjeKQwn/csCpPkmS1HxHHZnKzO6IuA54DGgHVmXmxohYVj+/ErgNOB64uz76052Z86sru3GTxta+4p593UwY08hAnCRJUuMaSheZ+SjwaJ9jK3s9/23gt5tbWnNMHvdmmJoxpcXFSJKk4hS/A/rhMLX3je4WVyJJkkpUfJiaNHY0gFf0SZKkShQfpnpP80mSJDVb8WHqHZPHAvCj1/e1uBJJklSi4sPUtIljmDx2FD985SetLkWSJBWo+DAVEZx0/AReeOWfW12KJEkqUPFhCmD29Ik8++JuDh3q9y44kiRJgzYiwtSFp7+Dl3a/wT/9v5dbXYokSSrMiNgS/OK5JzBr2vN87M/XM3v6RMaObqMtonYH5/qjt+2TJOln0yVzT+B3fuldLfv8ERGmJo4dxRd/ZyH3P7WVLTv3cuBgkpkcSkgg0+k/SZJ+Vo0d1dqJthERpgB+bup4brzk9FaXIUmSCjMi1kxJkiRVxTAlSZI0BIYpSZKkITBMSZIkDYFhSpIkaQgMU5IkSUNgmJIkSRoCw5QkSdIQGKYkSZKGwDAlSZI0BIYpSZKkITBMSZIkDYFhSpIkaQgiM1vzwRE7ga3D8FHTgZeH4XPUOPvk2GS/HHvsk2OT/XLsGY4+OTkzO/s70bIwNVwiYkNmzm91HXqTfXJssl+OPfbJscl+Ofa0uk+c5pMkSRoCw5QkSdIQjIQwdV+rC9BPsU+OTfbLscc+OTbZL8eelvZJ8WumJEmSqjQSRqYkSZIqU2yYiohLIuK5iNgcETe1up6RIiJmRcT/johNEbExIm6oH58WEf8QEf+3/nhcr5+5ud5Pz0XExa2rvnwR0R4R346Iv62/tl9aKCI6IuKvIuLZ+r8zC+2T1ouI36v//fVMRDwQEePsl+EXEasi4scR8UyvYwPuh4g4JyK+Vz/3mYiIZtdaZJiKiHbgLuBS4AxgSUSc0dqqRoxu4Pczcw6wALi2/ru/CfjHzDwV+Mf6a+rnFgNzgUuAu+v9p2rcAGzq9dp+aa0/Bf4uM08HzqLWN/ZJC0XEicAngPmZeSbQTu33br8Mv9XUfqe9DTbHcuwAAALwSURBVKYf7gGWAqfW//R9zyErMkwB5wKbM3NLZu4HHgQWtbimESEzf5SZ36o/30PtPw4nUvv9/0W92V8AV9SfLwIezMw3MvMHwGZq/acmi4iZwK8An+912H5pkYiYAnwQ+DOAzNyfmbuwT44Fo4DxETEKmADswH4Zdpn5BPBqn8MD6oeI+DlgSmY+lbVF4vf3+pmmKTVMnQhs6/W6q35MwygiTgHOBtYBMzLzR1ALXMA76s3sq+HzJ8B/Bg71Oma/tM7PAzuBP69PvX4+IiZin7RUZm4H/jvwQ+BHwOuZ+ffYL8eKgfbDifXnfY83Valhqr/5UC9bHEYRMQn4MrA8M3e/XdN+jtlXTRYRlwM/zsxvNvoj/RyzX5prFPAe4J7MPBv4CfUpiyOwT4ZBfQ3OImA28E5gYkR85O1+pJ9j9svwO1I/DEv/lBqmuoBZvV7PpDZMq2EQEaOpBam/zMyH6odfqg+3Un/8cf24fTU8zgP+bUS8QG3a+8KI+J/YL63UBXRl5rr667+iFq7sk9b6N8APMnNnZh4AHgLeh/1yrBhoP3TVn/c93lSlhqn1wKkRMTsixlBblPZwi2saEepXSfwZsCkz/7jXqYeBa+rPrwH+utfxxRExNiJmU1sc+I3hqnekyMybM3NmZp5C7d+Hr2TmR7BfWiYzXwS2RcS764c+BHwf+6TVfggsiIgJ9b/PPkRt7af9cmwYUD/UpwL3RMSCen/+Zq+faZpRzX7DY0FmdkfEdcBj1K7EWJWZG1tc1khxHvAbwPci4un6sVuAPwS+GBG/Re0vq38PkJkbI+KL1P4j0g1cm5kHh7/sEct+aa3rgb+s/0/fFuBj1P4n1z5pkcxcFxF/BXyL2u/529R2156E/TKsIuIB4HxgekR0AbczuL+zfpfalYHjgbX1P82t1R3QJUmSBq/UaT5JkqRhYZiSJEkaAsOUJEnSEBimJEmShsAwJUmSNASGKUmSpCEwTEmSJA2BYUqSJGkI/j/t3XT9hH/eMgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize = (10, 6))\n",
    "plt.plot(res[:, 0], label = \"State 0\")\n",
    "plt.plot(res[:, 1], label = \"State 1\")\n",
    "plt.legend(loc = \"best\", fontsize = \"large\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Irreducibility Check\n",
    "\n",
    "To check whether a Markov matrix is irreducible, we check whether $Q = (I + Z)^{n-1}$ has all positive elements where $n$ is the number of the states and $Z$ is a zero-pattern matrix of $P$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "True"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Irreducible or not?\n",
    "n = np.shape(P)\n",
    "Z = P > 0\n",
    "Q = np.identity(n[0])\n",
    "for i in range(n[0] - 1):\n",
    "    Q = Q @ (np.identity(n[0]) + Z)\n",
    "\n",
    "np.all(Q > 0)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Aperiodic Matrix Example"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 118,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[0 1 0]\n",
      " [0 0 1]\n",
      " [1 0 0]]\n",
      "[[0 0 1]\n",
      " [1 0 0]\n",
      " [0 1 0]]\n",
      "[[1 0 0]\n",
      " [0 1 0]\n",
      " [0 0 1]]\n"
     ]
    }
   ],
   "source": [
    "P_periodic = np.array([[0, 1, 0], [0, 0, 1], [1, 0, 0]])\n",
    "Trans = P_periodic\n",
    "for i in range(3):\n",
    "    print(Trans)\n",
    "    Trans = Trans @ P_periodic # becomes identity matrix after 3 periods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 121,
   "metadata": {},
   "outputs": [],
   "source": [
    "ini_dist = np.array([0.2, 0.3, 0.5])\n",
    "T = 10\n",
    "res = Markov_simulation_dist(P_periodic, ini_dist, T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 122,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x1e973476fc8>,\n",
       " <matplotlib.lines.Line2D at 0x1e973489208>,\n",
       " <matplotlib.lines.Line2D at 0x1e9734893c8>]"
      ]
     },
     "execution_count": 122,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOy9WWxb25rf+VucNI/UQFmyLIqyfY7P4FGkpXutmwDdQVV3oytBF9BVeWgEQeP2RVLphzx0VfdDvzSqgbwF3V3BRSFIKnnJxUVVV6GSukklSCepe+ue4+H4+Ay2LIkaLMk2B0mkRGritPphcVMUxZl7b57Y/AEHskluah3vtb79rW/9v+8TUkpatGjRosW7i6XZA2jRokWLFsbSMvQtWrRo8Y7TMvQtWrRo8Y7TMvQtWrRo8Y7TMvQtWrRo8Y5ja/YAijE0NCSnpqaaPYwWLVq0+M+GL774YkdKOVzsve+koZ+amuLJkyfNHkaLFi1a/GeDEOJVqfdaoZsWLVq0eMdpGfoWLVq0eMdpGfoWLVq0eMdpGfoWLVq0eMdpGfoWLVq0eMepytALIX5FCLEkhPALIX6nyPt/RQixL4R4lv3vf6/22hYtWrRoYSwV5ZVCCCvwe8B/CWwDj4UQfyqlfFHw0Z9LKf+bOq9t0aJFixYGUY1H7wX8Uso1KWUC+Anwa1V+fyPX1kQineCffvtP+ezNZ0Z8fdVIKfnp4y3ip6mmjoPEEXzxzyCTaeowAocB/t2rf9fUMQA8f7PP52u7zR4GrP5/EFps9ij42drP2D1u7r9HKp3hXzzaJJFq7hzlaA+++gk0uWT7f9j8D/yTb/8JRpSOr8bQjwNbeX/fzr5WyJwQ4ishxL8WQnxU47UIIX4ohHgihHgSDoerGNZ5bBYbf/D8D/hXa/+q5mv15PmbA/6XP/qaf/Fws6nj4MWfwL/8n2H9PzZ1GP/8xT/n7//Hv8/O8U5Tx/F//myRv/cvvjRkEVWNlPBH/yP8+f/WvDEAoaMQv/3z3+YPnv9BU8fxC/8O/+v/+w1/9s2bpo6Dp/8M/vh/guDzpg7jX6//a3669FOEELp/dzWGvthvLVwtT4ErUsqbwP8N/EkN16oXpfx9KeU9KeW94eGiWbxlsQgL90bv8fDtw6Yu5qVADIDPmu09apN2/edNHcZKZAWAJ4HmZjovBWKEY6eshuPNG0Q8BEe7sPk5pBJNG4Y/4gfgUeBR08YAeWtltdlrJRtJ3mjuWvHv+/H0ewz57moM/TZwOe/vE8C5R7CU8kBKGc/++WeAXQgxVM21euIb8xE8CrIZa543vRxUk/fh2i7JdBO3pFp4YP0vmjcGwB9VRuVh4GHTxrATP2UnrgzrL5tpVEJZg5I8gtdfNG0YK1H18F3cXWT/dL9p41jKrpWm3hP4TqyVZCbJ+v46M/0zhnx/NYb+MXBVCOEWQjiA3wD+NP8DQgiXyO43hBDe7PfuVnOtnnhdXqC5noo2eQ8Tab553bxFlJu8b76Ek4OmDCFyEsmFbB69bd490R6+AL/0N9PQ58Xmm2hUtIevRPIk2LydlnZftiPHbO0dNWcQ6RTsLKk/b/wlZNJNGcbmwSapTKp5hl5KmQJ+C/hzYBH4qZTyuRDiR0KIH2U/9uvAt0KIr4D/C/gNqSh6rRH/IwBXeq8w0jnSXKMSiLFwTYWemrYlPY5C7A3M/Bcg07DZnANqzaB8b/x7bMY2CRwGmjKO5WyI4AfXhvl8fZdMpkmhvdAL6ByCsZtNDRP4I37ujNyh3drO48DjpowhnZGsBOP8oNlrZW8N0gm1Vk734e1XTRmGtsu6OnDVkO+vSkcvpfyZlPKalNIjpfzd7Gs/llL+OPvn/0dK+ZGU8qaU8r6U8pflrjUKIQRel5dHgUdNidPHTpK82T/B5x7kA1dP8yZv+KX6eed/AGtb07zH1egqAH/zg78JNG+ntRyK099p57+9eYnoUZLFQHN2OIRfwsiH4F6ArYeQPDZ9CBmZYXV/lRvOG9weuc3Dt80JqW3tHXGayvBffeJiqLuNX6426bA+nN1leX+ofjZprfgjfizCgrvPbcj3v3OZsV6Xl72TvZw3aSbLQXXQd320h3nPEI839jhNNWErqMWCL92GSR+s/yfzx4Dy6HvsPXx//PsMtA00zagsB2JcG+1hfsYJNMl7lFKFbkY+BPcPlBe5Zf6/x5v4G45Tx8z0z+Ad8+KP+puiiNJCnNddvcx7nPxydbc5IorQIiBg6gEMf9BUp2iyZ5I2a5sh3//OGXrfmA9ojve4nJu8Pcx5nJymMny5GTV9HIQWwdENfZdhagEC3yitsMmsRFaYGZhRiijXvabstKSULAVjXB/tYayvA/dQV3MO//a3IBFXhn7yPghrU4yK5gDNDMzkzrSaoYjSwmlXR7qZ8zgJxU5ZDR+aPg5CL2DQDY5OZew3P2uKIsof9RsWn4d30NBf6r7EePd4U+L0S4EYnQ4r4/0deN2DWESTvEfNcxRChQkANn5h6hCklOcmr8/lI3AYYDu2beo4AgcnxE5SXHP1ADDncfJofY+U2Yoo7SB25Aa09cD43aZIXzVD7+nzcMN5gy57V1OcoqVgjMuDHXS12Zj3ZHdazZAkhxbVPQG1VpJH8OapqUM4SZ2wGdtkZqBl6GvCN+bjcfAxaZNP0JeDMa6O9mCxCPo67Hwy3tckQ/9CGXqA8Ttg7zLdewwfhzlIHOQMvXdMeY9myyw1rfb1UWXo5z1O4qcp8xVRWjht+AP1072gJJansdLXGMBKZIVLXZfodnRjs9i4N3qvabtf7Z5MDnYy3t/BZ2bH6VOnsLt6tlamvg8I09fK+v46GZlpefS14nV5iSViLEWWTP29y8E410a6c3+/73Hy5VaE44SJD5x4WCXlDGcnr9UOV+ZMn7y5EEF28k71TjHcMWz6Tmsle25ybVTdl/vTyns0PXwTegk9l6CjX/3d/UApol6Zq4haja6eS8qZdc3y6uCVqYqoRCrDWviQq1lDL4Tg/rSTz1ZNVkTtrKh7oD18OwfB9XHT1srVfmMUN/AOG3owV7u9Gz9lJ37K9WyIAGDeM0QyLXnyysT4uOY5al4KKO9xZwliQdOGoWVfattRIQTeMfMVUUvBGCM9bfR3OgAY6m7j+miP+XVv8ndZAJd9YHXAhnlGJZVJsba/di5EoJ1pmSmz3Ng9JJWROY8e1E4rcpTMHdKaQn44TcP9A9h6ZKoiyh/1Y7PYuNx7ufKH6+SdNPTDncO4+9ymhgmWc57j2eSdnRrAZhHmeo9FJ68WpzcvJuyP+hlsH2SwfTD3ms/lY/dkl7X9NdPGsRyMnXv4gorTm6qIyqQhvHTe0Ns7lLE30XvcjG2SzCTPeY7XBq7R19ZnqiJKC6flr5U5TxN2WqEXYLGBMy9k4l6A9Kky9ibhj/px97mxW+yG/Y530tCD8uqfBp+SzCRN+X35ihuNToeNW5f7zZ+8HYPQPXL2mutTaO8zVWbpj/ovbEVnXbMAphmVTEayHIydMyigjMpJMsMzsxRRe+vKeOQ/fEGpPN5+bZoiKrfLyosFW4SF2dFZU3day8EYVotgergr99ql/g6mnJ3mxulDi+C8CjbH2WuTc6YrovwRYxU38A4bet+Yj6PUEc93zKlItxSM0ddhZ6TnvA523uPkm+0oByfmPHByKoL8CngWK1z5vmkqj4zMsBpdvaAimOiZYLx73LQwwVbkiJNk5lyIAOC+24kQJqo8tKScfI8esjstCa9+eeESI/BHiyfleMe8vD18y3bcHEXUcjDGlLOTdrv13OtzniEerpmoiAovXrwn7b0q/8Sk3W88EefN4RtD4/PwDhv6e6P3APP09CvBGNdGuy+UGJ3zDJGR8HjdBK9Nymz25QcX33MvQGQdosYXfHt7+Jaj1FHRSnxa5nJGGr+YtXDa1dHuc6/3ddr5+FKfeTstLZw2fP386+N3wd5pmvfoj/q53HOZdlv7udd9rmzuiUlnWsvB+IVdFiinKHaa4vkbEzKXE4cQ2bho6MFURdTqvsoeb3n0dTLQPsD1geumTF4pJUuBiyECgNuT/bTZLOYYlYPXcHpQYvI+UD9N8Oq1EEExL8U75uUgccDSnvGKKC2cdrWEUfly0yRFVOgFDEyBo+v86zaHChWYaOiLGRR3n5uhjiFTzrROkmk2dg+LrhVTFVFamZBShj6TUuWkDUYrE9Iy9A3gHfPyLPyM0/Spob8neHDKwUnqwqEfQLvdyt0rA+ZM3mIHsRrDH6qCWiZsSbUCTaU8ejBnp7UUiDEx0EF328WOmfc9TpJpyRevIoaP41xSTiHuByqEEA8ZOoTT9CmbB5tFDYoQglnXLI8Djw2P0/tDcaSk6FoZ7mnj2mi3OSG1cmvlsg8sdlMewCuRFdqt7Yz3FO3HpBvvtqF3eTlNn/J1+GtDf48mCSvmpYDyHhffHhA5NDi1ujApJx+LRRmV9b8wvGWaP+rH1eWix3Hx32Okc4Sp3ilTDH1+Uk4hs1ODWUWUwYd/qQTs+ot7jmCaImpjf4O0TJfMvvS5fOwc77C+v27oOIopbvKZ9wzxeH3P+PaCoUWwtaudViGOTrjsNcXQ+6Oq2YhFGGuK32lDf3f0LhZhMVzlsVxh8s55hgCM126HFqFnTCV+FMO9oMI7e8bKG1ejq2W3or4xH08CTwxVRCXTGVbD8Vzpg0K622zcNEMRtetXYYBSHr3rJrT1GW5UcmVwSxz6mZW5vByK4bBamHJ2Fn1/zuPkOJnmq22DFVGhRXVmYrEWf9+9oEoWHxu74zO6xo3GO23oexw9fOT8yHCVx3IwxlB3G4NdjqLvfzrRR6fDarxRCS0W9+Y1prLeo4Eyy1QmxVp0rezknXXNcpQ64sXuC8PG8Wr3kGRa5jJiizE37eSb1/vEjFREldtlAVhtcGXecEO/Gl3FZrEx2TtZ9P2J7gnGusYMP9NaDsSYHu7CZi1uejRFlOENYkKLZ9njxZh6AEjVjMQgoidRdo53WoZeD2Zds3wd/pqjpHEdbFRSTmmDYrda8LoHjY095pJySniOAE6PSsM38EB2K7ZFIpOoaOjB2GzMpcDFBLZC5j1O0hnJ4w0DFVGhRaXLHiojn3MvqF3WvnHyRn/Ez1TvVMmkHK2Xw+PgY0MVUcvBeNH4vEZfp52PLvXy2ZqBIbXjiGrMUyqcBjBxD2wdhobU8iuJGs07b+h9Lh8pmeLL0JeGfL9KyikuF8tn3uPEH4oTOjgxZBxENiB1XH7yatUsDYzTVzN5B9sHuTZwzdCQ2lIwhkWAZ7j0A/jOlQEcNoux3mNoUWVe2srUGdfi9AY+gFeiKxW12r4xH/un+yxHlg0ZQ+wkyevocRVrZYinr6KcJA1SRIU0xU0Zp8jWpspJG7jTKqwHZSTvvKG/NXILm8Vm2OHfduSY42S65KGfxty0itMb5tWXUxHk434ARzvn+5fqiD/iRyCY7psu+zmvy8uXoS9JpI05oF4OxJga6rqQlJNPu93KnUmD4/SFNW6KMXJDZTMbZFSOkke8jr+u6Dkanbmc35inHHPTThLpjHGKqGL1oIrhfqA+Gw8bMgytMc9o56gh35/PO2/oO+2dfDr0qWGxx5zipsx2FODGpV56223GeY+lknIKMVjloSXldNg6yn5OU0R9FTamR2c5xU0+854hFgMGKaISR9mknAoPX00RtfFzQ3Za1Wq1XV0urvReMSyktlKkTEgxZt2DWC3CuBLf4Zfg6IG+ifKfc/9A/TRorWiNeQqTLI3gnTf0oBQFL/ZecJDQP+Mul5QzUjpEAGC1ZEuxGuXRhxehfxLayo+D/kklKTPIe9TkYpW461KKKCOMipaUUyxRqpB5jxMp4eG6AfdlZwmQxTOVC3EvqC5UEf3ljbWECLwuL0+CT0hlUrqPYykYo8OuGvOUo7vNxs2JPuOkr6FFdU8qGdixW+qBYMBa0RrzVLNW9OD9MPQuLxmZ4YvAF7p/93Iwxnh/Bz3tlSvPzXucbO4dsbVnwMFwuaScQqay3qPOjVkS6QSvDl5VZVB6Hb3cGLxhSJhgNRwnIyuHCAA+neg3ThFVbTgNzrxHA4yKP+pXSTndlZNyvGNeDpOHhiiilrNlQiyWyh7svGeIr7b3iZ/q/MCREoLPK4dtQCmipr5nyD3ZOd4515jHaN4LQ39z+CZt1jZD4vSq9EEFLzqLpqfX3atPJWBnubrJC8qonOyrXrI6sr6/TlqmuTpQXYGm2bFZvt75muOUvrW/zyqJVr4vDpuFe1ODxoQJQi/A2gYD7sqfdc5At8uQA1l/1M90/zTWUprxPGZHVZzemLVSWbSgMWeUIuowDMd7tTlFe6uw/1rXYVTKa9Cb98LQO6wObo3c0n3yJtOqU06l+LzGtdFunF0OPtfbqOytlk/KKSRX90ZfT6VWFYHP5SOV0V8RtRSI47BauOLsqvxh1E5rJRQnFNNZERVahOFryjOshIGKqFrK4Do7nMz0z+h+prV3mLjQmKccd68M4LBa9H8AV3sQq2HQmVZhYx6jeS8MPSijshxZZu9EPw/h1e4hifTFMrilEEIw53Hyy9VdfWuK1Dp5e1wwdE13Q78aXcUmbEz1TlX1+dsjt7EJm+5GZTmoknLsJZJyCsk1p9bdqLys/uELyqgchlQ+hE7sn+4TOg7VFCLwjfl0V0QtVygTUki73cqdK/36x+mrkVbmM/oxdAwY4hQVNuYxkqpWghDiV4QQS0IIvxDid8p8blYIkRZC/HreaxtCiG+EEM+EEE/0GHQ9aCneTwL6DaFYV6lKzHmcBA5OWN851G0chF6CsKgmCtXiXlB10NP6ZYWuRFe40nsFu7W6Tjmd9k4+Gf5E951WsWYj5fjoUh897TZ9Df3JPhxsl89ULsSAnVY91RFnXbOcpE90rRFVq6EHJUl+/uaA6JGOiiitMU/XcHWft1hU03Cdd1qVyoToTUVDL4SwAr8H/CpwA/hNIcSFx2H2c/8A+PMiX/NXpZS3pJT3Ghxv3dxw3qDT1qmrUVkKxBACZioobvKZNyJOH3oBgx6wt1f+rIZ7AZKH8Ea/sIk/4q95K+p1eXm++5xYQp/a3/HTFNuR46pDBKAUUT63zoqoWj1HUGqo/kld+8jmGk9XeW4CqpeDQOiqiFoKxOhttzHaWyZxrID5GU0RpWOcvlhjnkq4f5BVRG3oMoSMzJhW40ajGo/eC/illGtSygTwE+DXinzu7wF/BBhbb7VO7BY7d0fv6qryUJ1yyiflFDLl7GSsr11flUeoSKecSlz5vvqpU92bo+QR2/Htmievb8xHRmZ4GnyqyzhW6vAcQYVvXu0esR3RSRFVazhNw72gDmQz+pQhWIms0G3vrikpp6+tjw+dH+pa4Ezr3VuLZvzmRD8ddqt+Oy0p61srucxlfR7AWmMes+LzUJ2hHwe28v6+nX0thxBiHPgbwI+LXC+BfyuE+EII8cNSv0QI8UMhxBMhxJNw2JhMNN+Yj42DDUJH+jyLloLVK240hBDMTTv5fHWXTEaHrWDyWNVJqcVzBOhywugnuk1ereF3rSqCT4c/xWFx6GZUcoqbGg39nN5x+tAiOLqh73Jt100twEkUgvooojTPsdakHK/Ly1fhr3RRRJVrzFMOpYga0C9Ov78NiVjthn7oGnSP6rZWyjXmMYpqDH2xGVJoof4h8NtSymLC7O9JKe+gQj9/VwixUOyXSCl/X0p5T0p5b3i4yvhZjejZ9OIkmWZj57BmgwLKqOweJlgO6RCuCGtJOTVOXlCeytYjSDauNqm3QFObtY3bI7d1O5BdCsTpsFuZGCiflFPI9dEeBrsc+oVvwtlKopYa9Q46dgJrJCnH6/KSyqR4FnrW8DhCsdKNeSox7xliORgnHNOheVC4jnAaqDDPlH6Zy+Ua8xhFNbNwG8h3SyaANwWfuQf8RAixAfw68I+EEH8dQEr5JvszBPwxKhTUFK4PXqfX0auLUVkLH5KRxdvUVUJX77FcS7RKuBcgdQLbjcdi/RE/bdY2JrorpJUXwTvmZSmyRPSk8RrkK6EYV6tMysnHYlE7rc/0UkRp2Ze10ntJHarr4D3unuwSPY3WFJ/XuDN6RymidHCKtGYjV0fqMfRqrejSyyEXTqvjvrgXIB5U+SoNshpdZbRztGhjHqOoxtA/Bq4KIdxCCAfwG8Cf5n9ASumWUk5JKaeAPwT+jpTyT4QQXUKIHgAhRBfw14Bvdf0/qAGLsDDrmtVl8i5XWbejGBMDnVxxduoTpw+9AKsDBssXESvKlTml1tHBqPijfqb7qkvKKUTbaT0ONv7AqSdEoDHncfJ2/4SN3Qbj9PGwSsyp1XPUyCmiGssKbaQ6Ype9i4+HPtbFKTpT3NQW5gT46FIvPe02ndZKtjFPx0Dt1+oYp/dHaxctNEpFQy+lTAG/hVLTLAI/lVI+F0L8SAjxowqXjwK/EEJ8BTwC/kxK+W8aHXQjeF1eXsdfsx1rrPb3UjCG3SqYqjIpp5C5aSefr+2SbjROH1pUMcQqJY3naO+DS7d1SQZZia7UrSL4aOgjOmwdDR+URw4ThGKndYXTQMedVlgrfVDHLgtU+CYRg7eNhU1ySTl13pdZ1yzPd58TT8QbGsdSQDXmcXZXr7jRsFkt+NyD+nn09d6TgSl13tKgodca85gZn4cqdfRSyp9JKa9JKT1Syt/NvvZjKeWFw1cp5d+SUv5h9s9rUsqb2f8+0q5tJjnvsUHp2HIgxvRQNw5bfTlncx4nsZMUL940WGitHhVBPu4FFbpJ1K/rP0gcEDoK1e2laIqohu9JlZVESzE91MVob1vjh3+11LgpxpQWp29MEaUl5Tg7nHVd7xvzkZZpnoYaU0RVasxTiTnPEOs7h7yJNnAwXE1jnnJomcsbjSmiqmnMYwTvTWashqffw2D7YMMqj+VQrG6DAmfeY0NG5eRA6XsbMfRTD1T5hM3P6v6KepJyCvG5fKztrxE+ql9x1UiIAJQiat4zxOdrDcbpQ4sqPNBdZ53xriGVkdmg97gSXWnowO/m8E2liGpgp5XJSFZC1de4KYYumcuRDXUe1ahTdByBYP3RZzO7SuXz3hl6IQQ+l4/Hbx/XvZgPT1Ns7R1zrYZEqUJGetqZGeluLPaopcqX631Zicn7YLE3pPJYiSgVQSOGfnas8WJay8E4Pe02XL01JI4VMDftZCeeYCXUQLhC60faSJ3xqQew+RBS9alNpJQNZ1+229q5OXKzoZ3W6+gxR4l0Q4b++mgPA532xhRRuX4NDTpF0FCo0x+trjGP3rx3hh6UUQkdh9g42Kjres0INOLRg/JUHm/skUzXuRWsNyknH0cXTMw25D36o346bZ2MdY3V/R0fDHxAj6OnIaOylG020kgjh9xOy1/nTqvepJxC3AuqNeR2fSU7AocBDpOHDYcIvC4vL/de1q2I0hQ3jRh6i0XViGpIEVVtY55y9I2rDPRG1krEz0TPRMXGPHrzXhp6n8sHULeiYDlQX1JOIfMeJ0eJNF9v1ykrDC2CvRP6rzQ0DtwP1MHfcX3j0FQEjRhYq8XK7Ohs3WECKaWqcdPgw/fyYCeXBzvq32kdvIHT/cYN/ZX5hhRRuTK4dUgr8/GN+ZBIngTre+AsNRhO05jzDPE6esxmvb0cQi/UOqnUmKcS7gXY+Mu6FVFmlz7QeC8N/eWey7i6XHWHCZaCMdrtFi4PdjY0Dp/biRDU314w9KK+pJxC3AsgM0rSVwf+iF8XFYF3zMt2fJvX8dprf4djp0SPkg0/fKFBRVSjB7EaHf0wdrNuQ6/FghtNyvnY+XFDiqhaGvOUY25aO9Oqd63U0JinHA0oomppzKM376WhF0LgdXl5HHhMRtYeNlkOxpgZ6cZaY1JOIQNdDj509TZ/8k7Mgq29rtjj7vEukdOILpM3l7lcx05L8xyvNug5gsrGPDhJsfi2DkVUo9LKfHKKqNq9WH/Ez0jnCL2O3oaGYLfauTNyp+6Q2nIw3rA3D+AZ7mKkp62+A9lUAnZX9LknU/VXGK21MY+evJeGHpRRiZxGcgeJtVBrGdxyzHucfLEZ4SRZY1u/w11Vu7yeLL9CbG1w2VfX5NXLcwR1mDvYPliXUdFKRuvi0TeiiAotKrVNpw51xt0LkEnC1uc1X+qP6rPLArXTWt1fZee4tn+PVDrDaoOKGw2liKqzl0OuMY8Ohr57RDlXdawVTZ1mZukDjffa0EPtevroUYLgQf1JOYXMeZwkUhmebkZqu1BPzxGUUQl+C4e1LeZ6yuCWQgjBrGuWh4GHNS/m5UCMoW5HXUk5hYz2tjM93FXfTquRpJxCLt8Hi61mo5LOpFnbX9MtRFDvTmtj94hEOqObUzTncbITP8VfqyJKD9FCPu4F2Py8ZkWUP+rHJmy4e6toLakz762hH+seY7JnsmY9fa7ZSIOHfhpe9yBWi6h9S6pXLFhDa0698YuaLluJrNDf1o+zvb6knEK8Li+hoxCvDl7VdN2SjrssyCqi1mtURGUytXeVKkdbN4zfq1n6uh3f5jR9qptW+4PBD+ix99R8ptVImZBi1N3LIbQIwlpbY55yTD1QiqjXX9R0Wa2NefTkvTX0oFK8nwSekMpUf4K+VGcZ3FL0tNv5ZLyvDkP/QpUw6Klf0niOS7fB0VOz91hvGdxS+MayiqgajEomI1nR3dAPcZhI8/X2fvUXRTeUAdDLcwTlPb55qjpWVYneZXBtFht3XXdrNvT1NOYpx+XBTiYGOmoXL4QWwVljY55yTH0PELWvlToa8+jFe23ofWM+4sk4L/deVn3NSjBGT5uNsT6dJg3Ke3y2FeXwtAbJVj2dcsphtakiZzVMXj2ScgqZ7JlktHO0JqPyOnrMYYNJOYXcn9ayMWsIZdXTVaoSOUVU9ZnLmrTS3adfiMDn8rEV2+Jt/G314wjV3pinEvMeJ5+v19jLQY+8hnw6BmpWRNXbmEcv3mtDP+uqPRtzKaDK4OrlwYKKPaYykscbVbZM05JyaulHWg3uBaVOOKhuMQePgsSTcV0nb74iqto4/UpICxHo4zkCDHY5+MDVU1uYQIsFN5KUU8jELFjbalJErUZXmeieoNPemPw3n7rXik7evAWIWEMAACAASURBVMacx0n0KMlioEpFlNaYp5GM2GK4H9SkiFrfXwcayx5vhPfa0A91DOHp81R9yKQl5egVc9S4d2UQu1VUb1RiAdWFSE/PEc5KsVZpVHKlD3TejnrHvOyd7OUOeiuxFFDnJvX0BijHvGeIJxs1KKJCi9A3CW06jsPeDpO+mgqcGVEG9+rAVQbaBqo29CfJNBu7R7qvlbnpbJy+2lBnI415yuH+AaQTsFXdGZ+2y2oZ+ibhHfPyNPSUZDpZ8bPh+CmRo6SuIQKADoeV25MD1U9evVUEGqOfQHt/1UalkXrn5ai1E9hyMMalvnZ6G0zKKWTe4+Q0leHLzSozhvUOEWi4FyDwDRxV3vEl00k29jd0L4Or9XJ4+LY6RdRa+JB0Ruq+Vlx9NSqi9BYtaEzWpojyR/w4LA4u99TYWlIn3ntD73P5OE4d8+1u5Yp0ywH9tNqFzE07+fb1PvtHlR84Z5NXZ6NiscDU96tWefijfkY6Ruhr69N1GJe6LzHRPVF1NqYKp+l/T7zTg1hElSqPdFK/pJxCprSdVmVF1MbBBimZMkSr7XV5CR4F2YptVfzsWTjNmLXyaH2PVDWKqPBi/Y15ytHWA5fuVL371Vo61tOYRw/ee0N/z3UPgajKqCznsi/1n7zzHicZCQ/XqzAqoUXoGlHlbPXG/QOIvlJlXStgZKcc35iPJ4EnpDPlwybpjMQfjhtiUHpziqgqDmT31tRWXm/PEWD8Dti7qvIejdplgdr9AlVJkpcCjTXmKce8Z4j4aYpvXlehRAotwtB1JTbQG/cCvH6qyoVXoJHGPHrw3hv6vrY+Phj8oKowwXIwxmCXg6Fuh+7juDXZT7vdUp33GK6zH2k1VNmcOp1JsxZdMyzLz+vyEkvGeBkpr4h6tXtIIqVfUk4hc54hvtyMcpSooIhqpB9pJax2VeSsSkNvFVZdFTcaU71TjHSMVHWmtRyM4R7qqrsxTznuT6us46rCN/X27q0G9wLIdMVeDlpjnmZkxGq894YelFH5KvQVJ6mTsp9TSTn6Km402mxW7l0ZrByn1zspp5DhD6BruKJReR1/zUn6xLCWaDmVRwWjsqxzXkMhZ4qoCpnLoUVVbXLomiHjwP0AdpbUQXwZ/BE/V3qv4LDq74wIIZgdUz2XK8Xp9U5gy8fZ3aYUUZXWih6Necpx2asUURXWilb6oBk1bjRahh61JU1kEnwV/qrkZ6SULAdihhkUUEblZSDGbrxMavX+JiQPjZu8+S3Tyixmo1UEw53DTPdNVwwTLAXiuiblFDI7NYCtmszl0AsVB7YbVGfcXV2c3ugyuD6Xj72TvZzxKobWmMfotfLk1R6nqTKhvbABeQ352DuUsa9g6PVozNMoLUMP3B29i1VYy8bpc0k5BsSCNbSWaZ+vlVFXGKUiyGfqAcTewm5peaMZBZq8Li9Pg09JZkofUC8HY0wOdtLhMOaQq9Nh4/Zkf+U4feilcQ9fANenKhO6jCLqOHXMVmzLUINSTZzer1NjnnLMe4Y4SWZ4Vk4RZZRoIZ8qFFF6NOZplJahB7rsXXw09FHZOP2KVuPGQC/lk/E+utts5asm5pJyDIo7wpn3WMao+CN+xrvHdU3KKcQ3phRRz3eel/yMkSECjTnPEN+83ufgpMQDJ3miKiQa+fC1WNUDuIz3uLa/hkQammY/3j3OePd42ZCa3mVCiuF1K0VU2Th9aFEdYvdNGjYOtVZk2Z3WanS14cY8jdIy9Fl8Lh/Pd55zmDws+n6uU86IcZPXZrXgdVeI04deQu8EtDdWZ7wsg9Pqd5Q5kDVDRXBvtLwi6jSVZmPn0FCDAkrOl5HwqNROa2dZlSkw8uELytBHNiC6WfRtPZq0V4PX5eVJsLQiajkQo83WeGOecvR12Pl4vK+8eCH0QmUpN9qYpxyX7qgub2Vkls3qKpVPy9Bn8Y55SckUT4NPi76/HIjh6m2nr9PYynPzHidrO4cE9kscDBuVlJOPEOrwb+Pn6vC3AC0px+jJ29/ez/XB6yV3Wus7h6Qy0tAQAcDtyX7abJbS3qMZ4TTI22kVNypmJeV4x7wcJA5YiiwVfX8pqMqENNqYpxJzHidfbkY4TpSI0+vVmKccNgdMlq4RtXu8y97JXsvQf1e4NXwLu8Ve0qgs6dCPtBq0phefrRUJ36RTSnlhtKEHZVSOds9CRXm8OnhFSqZMqcTndXl5FnrGafriAfWSTr17K9Fut3JvaqB0SC30Aix2VSHRSEY+hM6hkkZlJbrCdP80NosBmvE8KtWn17MxTznmPUMk05Inr4rstHKNeUxaK+GXEAteeMvIvIZaqMrQCyF+RQixJITwCyF+p8znZoUQaSHEr9d6bbNpt7Vzc/hmUUOfzkj8oTjXdWiJVokPXb30d9qLl2I1MimnEK1lWpEtqX9f3zK45fC6soqo0EVF1HIwhs0icA/pn5RTyNy0UkTtHSYuvhl+qWSVRtcZz99pFVFEadmXRjPSOcJU71TRtbJ/lNS1MU857l0po4jSuzFPOdxl1oqOjXkaoaKhF0JYgd8DfhW4AfymEOKCpcl+7h8Af17rtd8VvGNeFncX2T89n3G3uXfEaSpjSEZsIRaL4L7bWTxMYFSNm2L0X1ax+iLeoz+iknKm+qYMH0ZOEVVE5bEUiBuWlFPIXLbpxefFYsJ6dpWqhHsBDl6rh34esUSMwGHANM/RN+bji+AXFxRRy9nSB2bsfrvabNy63F9irZgUTgNw3YS2vuJrJerXtTFPvVSzQryAX0q5JqVMAD8Bfq3I5/4e8EdAqI5rvxN4XV4kkifBJ+deNytEoDE/4+R19JitvYISqOGXgDAuKaeQqQdKTZA+nxXqj/q53HOZNmvjbfsq0e3o5iPnR0VbPq6EzAmnAXw60UeXw3oxfHMaU4ejRmVfFqJ1AitQROWSckzYZYFaK0epowuKKG2tmBG6AXWm9fV29KIiKteYx2X8IKw21YykhFPk6fc0VXED1Rn6cSC/itF29rUcQohx4G8AP6712rzv+KEQ4okQ4kk4HK5iWPrz6dCntFvbLxiVsxo3xoduQIUJoEhz6tALGHSDwzg1wzncC3B6AIHzYRN/1G/qVnTWNcs34W84Sp49+I4SKTb3jkx7+NqtFmbdgxe9x3D2QNIMzxHULqvn0gWjkosFm9TBqFR9+uVgjO42G5d0bMxTjvvZGlGP1wvi9Ho35qmEewEi6xA9M3dSyu+E4gaqM/TF/qUKA4T/EPhtKWXh8Xc116oXpfx9KeU9KeW94eHhKoalP3arnTujdy7I+ZaySTmdDmMPuTRmRroZ6m67aFTMUBHkU0TlcZI6YfNg09TJm1NEhc4UUf5QHCnN8xwhq4gKHxI8yFNEmRlOg7PM5fXzcXqzk3IG2ge4NnDtgqFfChhXJqQYdyYHcNgs5+P0UpobToOiZ1paYx6zdlnlqMbQbwP5eq0J4E3BZ+4BPxFCbAC/DvwjIcRfr/La7xSzrln8UT+7x2cTZzlgjopAQwjBvMfJZ6u7ZzVFkiewu2ru5O0eUdrwPO9xfX9dJeWYaOhvj9zGZrGdMypnIQJzdlmQ15w636iEXoKtA/qnTBuHUkTtnMWhOQsRWIR5QrpCRZRRjXnK0W63cu/KwHmnKBZQPXbNdIpGbkCn89xaMaoxTz1UMyseA1eFEG4hhAP4DeBP8z8gpXRLKaeklFPAHwJ/R0r5J9Vc+13D51LNqR8HVfgmkcqwvnNoqkEBJbMMxU5ZDWcTuHZXVKU8Mw09KKOy+RmklNrE7BABQIetg0+HPj0n51sOxnDYLFwxoAxuKT4c66W3vSBzOfRCxeeNTMopJFdhNM+oNKEMrtfl5TR9ytfhrwHYiScMacxTiblpJy/eHhDRFFFm77Ig28shm7mcdc7MSmCrhoqzU0qZAn4LpaZZBH4qpXwuhPiREOJH9Vzb+LCN40Pnh3Tbu3NGRUvKMdNLgbO6N7kaK1rjab17X1bCvQDJI3ijwiYr0RXsFjuTPQamlRfBN+Zjce9MEbUcjHN1xPiknHysFsH9aef5bMzQovn3pH8SBqZyYYK9k72mJOXcdd3FIiy5nZZ2lmW2oZ+fUWsl18tB2+mYvlYenFNErURXGO4Y1r0xTz1U5YZIKX8mpbwmpfRIKX83+9qPpZSFh69IKf+WlPIPy137XcZmsXFv9F5u8i41afJODnYy3t9xtiUNvVCty5wmewdXvgeInPfoj/hx97kNT8opxOvykpEZvgh+ASijYtZBbD7zHidbe1lF1NEexAPm77LgrMJoJn3mOZocIuh19HJj8EbOKTJbcaPx6UQ/nQ5r3lrRGvOYLGnMKaKya+U7chALrczYosy6Znl18IrAYYDlQAyrRTA9bF6IAFScfs7j5PO1XTIZqSav86pKuTaTzkFwfdL0yfvp8Jkiav84ydv9E9OklfnMz+TF6c3Uahfi/oGKQwe+zsWCm3Ho5x3z8vXO1xwljwxtzFMOe7ZG1DmnqBkPX+cM9IzB+l/kGvN8F+Lz0DL0RfGNZeP0gccsZTvltNnM7/U4N+0kcpTkZSDWvMkLynvcekT8KMzbw7dNyfJzWB3cGrnFw8BDVoLmH8RqXB3pZqjbocI3ZmZfFjL1ffVz/ef4o356Hb0MdRjQWrICXpeXVCbFs9CzbOkD8xQ3+cxNO/GH4oQOjpTktRkPXyGyuSc/53Vs29DGPLXSMvRFuDpwlf62fh6+VUalGQYFzurePFraVH1cm2no06esLv8ZAJ6+5rRE8435WIms8OXrbcD8EAGondb9aSe/XN1BBhehrRd6L5k+Dnpcqhfq+l/kdlnNMLC3R25jEzYevn3IcjDelHsCZ4qor775JtuYx6QEtkLcC3AYxr+pEtpaoZvvMBZhYdY1y8O3j3i1d9i0yXupvwP3UBfbK1+qF5pl6CfnQFjPJm+TtqNaMa2Hbx/T5bAy3m9QN6cKzHuGCB6ccvLmW3VPmpX16F5Avvql6Qls+XTaO/l0+FN+8fpz4qeppq2VG5eUIurNcjbXohkePeRyT/yb6qC8mX1i82kZ+hJ4XV4CR2/BtteUQz+N+9NOTl9nhUrNmrztvTB+B//uCzpsHYx3F01uNpwbzht02btYOfiSa66epqWVq52WxBI2oWR0OdwPCGVOiCViTfUcZ12zrERfguXEdHWahtUi8E07OX3zrXrB6N4ApRi4Av1X8O8tGt6YpxZahr4EWss0W9dqUw79NOY9TibTr8hY25SkrllMPWDlNIKnd8rUpJx8bBYbd0fvspd50dSH75Szk497j2lLmpyUc2EgD/A7VMXMZhp635gPSQZr55qhjXkqMe9xMnKyRqpn3NjGPJVwP2DldI+ZJoU4i9Ey9CVw97ppt/Rj717jioGdcipxf9rJdbHFbodbtZNrFu4FVu02ZmzNW8gAHw3cAXuYMWeRcsEmIYTgv3YpPX9mqEmeI0DnIH6nymdopqH/dPhTLNjpHXhleGOecsx7hrgmXhPumG7aGACSU99nw2Zhxt7Eh00BLUNfAiEEnenr2LvWTE3KKWS4p40bttcsZSaaNgaAyMh1dmxWZo6Lt1o0i0GL8qBTjpWmjmOuWxVp9QtjuzlVYqVnmKF0mn6rOUXEitFmbcORmsbWtdq0MQBcG25nxvKaxVRzQosam043KSGavlbyaRn6Mhzuu8lYDljfX2/eII4jDMk9PouPkEhdbOtnFv7D1wDM7L5q2hgAjuMjyHQHgcQ3TR3HVbHFjuzlF2+aW3521QoziQRsXyzjbBbpjCQeneJEbBM5iTRtHGJvHQcpfnEwfFYjqgmspKIAzOw00W4U0DL0JTg4SbK7o7bFxZpemEa29MHz1DjPtqJNG0auxk1gGY6bt5iXQ0dYTmb4eueLpo0BoCu6zKb1Suk+siaQkRlWT8LMJNIl2wuawebeEadxFS4p1jfANLI1bh4eutjYParwYePwR/xYAPfmU0gnK37eDFqGvgQrwRgyOciAY7Rkb0xTyE7eFXm5eMs0k/BH/PTYOhlJp2DjL5s2juVgjBHHx7w5fMN2bLs5g5ASwi85GbjGw/Vd0pnmeI+v4685Tp8w0zXWVEO/FIiROZ6g3dpRsueyKYQWkQj8crx0f18T8Ef9TLYP0ZaIw5svmzaOfFqGvgRLgTgA90ZneRx8TEY2KWwSWgRHDwNj7qZP3qsD1xC2jqK9Mc1ASslyIKYOZLnY9MI09rcgEafz8ifETlI8f7Nf+RoD8Eeyu6yxu/D6CSSaExNWmcpWbo/caa6hDy/C4DQDvb3NdYqifq46s2qsJj6A82kZ+hIsB2N0Oqz8lcl59k/3WY4sN2cgoUUY+YD5mWG+3Ixykizs7WI8Wqccz8BVmLzftMn7dv+E2GmKu5c+wNnuvNAgxjSyNW6ufHAXoGnhGy2c5vH8KmRSqpx0E1gKxrg82MH8pfus768TOgpVvsgIQouIkQ8v9nIwkZPUCVuxLWaGPoLRj1uG/rvOcjDG1dEefGNaNmYTjEpep5w5j5NEOsMXr8yPj4ePwxwkDpSEz72gxhQ3v92jVgb3A1cv3jEvjwOPm3Polg2nDVz5lKsj3U019Je6LtHt/itgsTfNqGiVRLXck6Z49XmNeeY8TnYPEywH46YPY31/nYzMqIxY9wJsPYTUqenjKKRl6EugJm83o12jTPVONeeQ6TAMx3swcoPZqUGsFtGU8I0WIrg6cPWsFGsTwjfLecXMvC4v4eMw6wdNUDaEFqF3HDr6mfM4ebKx1xRFlD/qV+UoHJ0wMXuu5aNZJFIZ1sKqTMj1gev0OHqas1byGvPMFfZyMBFtl3W1/6oqcJY6aaoiSqNl6IuwEz9lJ57I1e3wurw8CT4hlUmZO5C8TjndbTZuTvQ1xXtciSrNuqffA2M3wdHTFO9xKRBnpKeN/k5HrhNYUw7K8yqJznucHCXSfL1triIqmUmyvr9+VkvFvQBvn8GxuePIb8xjtViZHZ1tzu43r2T0xEAnk4OdTVsrdoudy72X4co8CMt3InzTMvRF0DxHrW7H7Ngsh8lDXuy+MHcgBfXO5z1DfL29T/zU3AfOanQVZ7uTwfZBsNpg6ntNmbz5/UgneiYY6xozP0yQSUN4OWfofW4nQpgfp9862CKZSZ6VwXUvgMzAq1+aOo7CrlLeMS+v46/NV0SFFlX4alA9+OazvRzMVkStRldx97mxW+zQ0Q9jt1qG/rvKSja2p03e2dFZoAmxx9AL6BiErmFAFdNKZySP1/dMHcaFZiPuBdhbhf3Xpo0hk5GshM6atAshmHXN8jhgsiJqbx3Sp7k2dQNdDj509ZoeUtN2Wbn7MnEPbO2mG5Xl4PnGPFqFUdPDN6FF1fgj25hnzuPk4CTFizcHpg7DHymyVrabp4jSaBn6IiwFY/R12BnpaQPA2eHk6sBV88MEoZfKm89Wabx7ZQCH1WKqUcnIzFksWGMq25zaxDj9VuSIk2TmXDEz35iP6Gk012HJFIo0np73OHlqsiJqNbqKRVhw97nVC7Y2pYgy+exkKRBjytmZa8wz0z/DYPtgc5yivHsyN52N06+Zt1biiThvDt8UGPoHkEnC5uemjaMYLUNfhOWAUhHkl8H1uXx8GfqSRNqkYloy2z4wb/K2263cudJvapjgTfwNx6nj85N39GPoGDDVe8z1I82rJHpWn97EmHBoERAwfD330vyMk0Qqw1MTFVH+qJ/JnknabXk1btwLEPwWDs0zbvnhNFA7La/Ly6O3j8xTRJ3Gs415ziqJjvS2M2OyImp1P9u7N3+tTM6pXs9NDt+0DH0BUkqWgjGuuc53lZp1zXKSPuHr8NfmDGR/GxKxC/XO56aHePH2gOiROQ+cXOPp/MlrsSivfv0v1APJBLRY8NWRs/vi6nIx2TNpbpggvKjKRTvOegifKaLMMyorkZWLTS2mVNMLs7z6k2SaV3tHF5qNzLpmCR2H2DjYMGUc7CypnxfWipNH63sk0+aE9nIJbPm7X0dXVhHVMvTfKQIHJ8ROLnbKuee6h0VYzDMqJRpPz884kRI+XzMnTn8hFqzhXlAZopENU8axFIwzMdBBV5vt3OveMZMVUaHFC/ekp93OJ+N9qo+sCZymT9mMbV68J5duZxVR5hh6fyiOlFzoDZDfc9kUQsV7954poszJXPZH/cUb80w9UIqok+ZkUEPL0F9gueAgVqPX0cuHgx+aV+AsFws+X+/85kQ/HXaraRphf9SPq8tFt6Ogb262ZZpZnooWTivE5/IRT8ZZ3F00fhCpU9j1F+1HOu9x8tVW1BRFlJaUc6Glo9WmJH0m3RMtnHa14L5M9kwy2jlqXkgttKgOogsa89yfNldP74/6me6bvtiYp0mKqHxahr6A5cB5uVg+XpeXr8JfcZw6Nn4g4ZfQM6Zi4Xk4bBZm3YOmeY8XVAQaQ9ege9QUo5JMZ1jbiRft9HXPdQ8wSRG161elBop0lZr3DJHKSB5vGL/TOpeUU4h7QSUPHbw1fBzLwRgOq4Up5/nGPEIIfGM+8xRRoRfqzKSgMc9Al4MbY72mhdQuqNM0JmaboojKpypDL4T4FSHEkhDCL4T4nSLv/5oQ4mshxDMhxBMhxPfz3tsQQnyjvafn4I1gKRhjuKeNwS7Hhfe8Y15SmRTPQs+MH0iBiiCfuWkny8E44ZixqdWpTIq1/bXiBkUIZVQ2fm54nH5j55BkWhb16Ic6hpjpnzHH0JcIEYBSRNmtgs9NMCr+iB+bxcZk7+TFN93mKaKWgjE8I93YrBfNyKxrlshpJPdQMpQi4TSNOY+TL15FDFdERU4i7BzvFG/Sbm+Hy96mZC5rVDT0Qggr8HvArwI3gN8UQhT+q/574KaU8hbwt4F/XPD+X5VS3pJS3tNhzIai1e0oxp2RO9iEzXijkklDeKnk5J3XUrwN9uq3Yiop50KIQGPqAcSDsGNswbcl7SB2tLvo+16Xl6fBpySNrv0dWlQKCufFxdzhsHJ7csAU79Ef9TPVO6WScgoZ/QTa+2H9Pxk+jpVgnOtl7gmYkLl8HIHY25JO0bzHyWkqw5ebxmYM5/o1lGrp6F6A4Ddw2Jy6SNV49F7AL6Vck1ImgJ8Av5b/ASllXJ5pqbqA5rV3aYBMRrISjJc0KJ32Tj4e+tj4yRvZUDUySkzejy710tNuMzz2WNXkBcO3pMuBGBYBnuESRmXMqxRROwYrogqScgqZ9zj59s0++0fGPnD8UX/xXRYoRZT7geH3JHaS5HX0uGg4DeBS9yUu91w2/kwr25inlFPkdStFlNFrpag6LZ8m1oiC6gz9OLCV9/ft7GvnEEL8DSHES+DPUF69hgT+rRDiCyHED0v9EiHED7NhnyfhsPmVEQG2I8ccJ9MlPXpQRuX57nPiCQMr42kHscPFDb3NasHnHjS85rY/4kcgzpJyChmYgr5Jw43KUjDG1FAX7fbizdHvjd5DIIzfaYVewHDpZuBz00oR9XDduPtylDzidfx16V0WKJlldNNQRVROtDBSZq24vHwR+IJ0xsCwSW6tFL8vPe12PjZBEeWP+umx9zDSOVL8A5dug6P7O23oizXFvOCxSyn/WEr5AfDXgf8j763vSSnvoEI/f1cIsVDsl0gpf19KeU9KeW94eLiKYemPFiIo5aWAUnmkZZqnoafGDUTzUvKScgqZ8wyxsXvE66hxB8Mr0RUu91ymw9ZR/ANCKO9x4+eQMe7QTYUISt+TvrY+Phj8wNidVuJQGc4SniPArcl+2u0WQ8M3FT1HyNtpGWdUCutBFcPr8hJLxni599KwcWiNeeibKPmReY+TLzejHCWMU0StRFaYGZg5l2R5DqtdJU816UC2GkO/DeS3up8A3pT6sJTyLwCPEGIo+/c32Z8h4I9RoaDvJMWScgq5OXITh8VhrHQs9AL6r0Bb6XHk4vQGGpWSKoJ83AsqThr81pAxnCTTbOweFlVB5eMb8/FV+CtOUieGjIPwEiBLhtMA2mxWZqeM3WmVVdxoDF+HrhFDjcpSQDXmGe8v4QRArj69oeEbLXu8lIFFrRWliDImc1lrzFPVWtlZNkURVUg1hv4xcFUI4RZCOIDfAP40/wNCiBmRfZQJIe4ADmBXCNElhOjJvt4F/DXAGIugA0uBGOP9HfS0FznkytJmbePWyC1jk0HKqAg0ro/2MNBpN8yoJNIJNg82y4cI4KzujUFGxR+Kk5HF5a75zLpmSWaSPAsbpIgKl48Fa9yfdrIUjLETN0YRtRJdod3aznjPhejpGdpOy8DM5ZWQasxjsZQ2sEMdQ0z3TRu308przFOOe1cGsVuN6+VwrjFPObSd1sYvDBlHOSoaeillCvgt4M+BReCnUsrnQogfCSF+lP3Yfwd8K4R4hlLo/PfZw9lR4BdCiK+AR8CfSSn/jRH/I3qwHIxxrcRBbD5el5eXey+Jnhhwkp9KKB10hclrsQjmPE4+W90xpKbI+v46aZku7zkC9I2r0rAGxR7PQgTl78vd0btYhdU4oxJ6AdY2GCxxXpFF22l9blBM2B/xM91fJCmnEPcCxANK+28AS4HSipt8vC4vT0MGKaLyGvOUo8Nh5fblAcOkr7ldVjFpZT6uT6C9zxRFVCFV6eillD+TUl6TUnqklL+bfe3HUsofZ//8D6SUH2UllHNSyl9kX1+TUt7M/veRdu13kWQ62ymnTMxRwzfmQyJ5EjQgLSCXlFPe0IOK07/ZP+HV7pHuw8j1Iy2sp1IM9wJs/CWk9Y+BLmWTcq44u8p+rsvexcdDHxsXJggtwvC1C0k5hXwy3kd3m82wOH1VIQLIi9Prb1R246fsxE8r7rJArZXj1DHf7hqwkS+RPV6MOY+Tb17vs3+s/wNHq3FTca1YrGc1okymlRmb5dXuIYl0puyhn8ZHQx/RYeswJk4fLp2UU4iRenp/1I9N2Jjqnar8YfeCKsD2Vv+wyUowzvRwF/YiWNdqgwAAIABJREFUSTmFeF1enu885zBpQO3vKsJpYKwiav90n/BxuPIuC2DADX2XDTEqpcqEFENTRBmyVkrUgyrGvMdJRsIjA3o5+KN+BtsHVWOeSrgXVKXNyCvdx1GOlqHPshSofvLaLXbujN4xJk4fWgRhLZqUU8j0UBcjPW2GeI/+iJ+pvins1tLnFTkMjNMvBWJV3RNQh39pmeaL4Bf6DuI4Cgevq3r4gvIe13cOebuvryKqpl2WEOq+bPxCd0VUNYobjf72fq4PXjdorbyATmeuMU85bk3202azGPIALpvXUEgTejlAy9DnWAqqpJyZMoqbfHwuH6v7q+wc63zAE1oEp0elTVdACMG8QXH6qkMEAN3DyqvS2dBrSTnVGBSAW8O3sFvs+sfpw1oZ3MqeI5DXnFpfo3KuSXs1uBfgaPcsxKETywWNeSrhdXl5FnqmvyKqoDFPOTRFlN4HskUb85Rj5EPoHDI9fNMy9FmWAzGuOEsn5RRiWIp3FSqCfOY9Q+zEE6yE9EvgOkoesR3frt7QgzIqm5+rCo86of0/VevRt9vauTVyS//EqSJdpcrxoauX/k677jutlegK3fZuRjtHq7vAbcxOSysTUlIzXoBvzEcik+Cr8Ff6DaJIY55KzHmcvAzE2NVREfX28O3Fxjzl0GpEmdjLAVqGPsdyqDrFjcYHgx/QY+/R16gkjlRP0hIZscUwwntc218DKiTlFDL1AFLH8Fq/sIlWSbSacxONWdcsL/desn+qY+3v0KLKauy7XPmzKEXUfbeTz1Z3dd1pabusag0sfRMwOK1rmEBKyVIgVrJMSDHujNxRiig914rWmKdMpnIh2lp5qGOcPtdspFanKPYWdld1G0clWoaebFLOzmFNBsVqsXLPdU/fybuzTKWknEIuD3YyMdCh65ZU68Fa9XYUYOp7gNDVe1wOxumwW5kYKJ2UU4jPlVVEBXRURGmlD6o1sKgGMa+jx2zt6ROnl1KyGl2t7Z5AVhH1C90UUcGDUw5OUlWH0wC6Hd185PxI391vDQexGp/mFFE6rpVsY56qzk00DFRElaJl6IHVcDYpp4bJC2pLuhXb4m1cp0y3OiYvKEXB52t7ZDL6eI/+qJ82axsT3aXTyi/QMQBjN3U29MpzLJeUU8gnQ58oRZSeMssaQwRwpojSy6jsnuwSPY3W5jmCMiqnBxDQJ2ySKxNSg1ME6qD8251v9VNE1SCt1LBZLXjdg7qG1LTGPD2OGv49Bqehd9zUOH3L0JOnIqhx8s66ZgEdm16EXoDVoSZCDcx5nOwfJ3nx9kCXYaxGV5num8ZaQTN+AfcCbD1SISgdWApWr7jRsFvt3B65rZ/KIx6Go52aH76e4W6GdVRE5XZZtRp6nRVRK3Ua+lnXLCmZ4mlQpxpR4ZfQc+lCY55KzE07WQsfEjzQ52C4ZGOecuR6OeiviCpFy9CjpJV2q2BqqHxSTiEz/TMMtg/qaOgXYei6aglXA3PTQ4B+cfqV6Er1yo583AuQScJW49703mGCcOy05ocvqINyf9SvjyKqhryGfIQQzE07+WxNnzh9xZLRpegeUWc+OhU4WwqUbsxTjtsjt7FZbPo9gGsULWjoeaaVyqRY31+vXlqZz9QD5UCETWiBScvQA8pLmR7qriopJx+LsDDrmuXh24f6HLqFFmvaimq4+tqZHu7SJUywf7pP6ChUW8xRY/K+asyhg/e4XEUl0VLo2py6TFepSsx7nIRjp6yGG1dEaUk5zg5n7Re7F2DzM1Veo0GqLRNSSIetg5vDN/UJqeUa89R+T26M9dLXYddlrWzFtkhkEvWtFYMUUaVoGXqyIYI6DAoo7zF4FGQzttnYIE4O4GC7rskLyqg8Wt8jmW5sK1hVGdxStPXA+F1dJu9KneE00FkRFXqhwgPdVUoa85j3qJ2WHuGbmvIaCnEvQPKoYUVUJiNZDsZrDtto+Fw+FncXG1dEVWjMUw6LRe209LonUKNoQaN/UmUvm9Re8L039PHTFNuR46oKNBUjp6dv1KhUWR2xFHPTQxwm0nzzurFFVFUZ3HJMPYA3X6oHVwMsBWP0tNsY7a0uKScfm8XG3dG7+qg8tNIHNShuNC4PdjDe39FwmEBKiT/ir89zhDNFVIMyy2oa85Rj1jWLRDaeuVxjXkMhcx4n25FjtvYaO0vSGvNM99V2ppbDrWUuG9vPFlqGvu7DJY0rvVcY6Rxp3Kg0OHnvT6s6G40aFX/UT5e9C1eXq74vcC+ATKtQQQMsB+I1JeUU4h3zshnbbEwRJWU2+7K+eyJEtsLo2m5Diqi3h285Sh3V79F3DMDYpw3vtBoJpwF8Ovwp7db2xp2iXGOe2sOcoF8vB3/UX74xTyXcP4DTfXirYyJZCd57Q7/coKEXQuBz+XgUeNRYnD60CPYu1ZqvDpzdbXzg6mk49lhzUk4hl72qnG8DRkVK2VA4DXTaaR28UQuxTkMPyqhEj5IsBurf4VRdBrcc7gV1SJ6sX9e/VEVjnnI4rA5uj9xuvMBZ6IVqY+moTTyhMTPSzVB3m25rpW4M7uWQT8vQB+O02y1cHuys+ztmXbPsnezl4tt1EXqhOgNZ6r8lcx4nTzYinKbq2wpKKVVLtEYmr71DGfsGJm8odsr+cbLuEAEoo9jf1t+YodcOYmvIVC5ED5WHJq2sO3QDqo9sOqHkr3WyHKzcmKcS3jGliNo9bsCbDi02dE/yd1r1OmeJdIJXB68auyc9o2pXYkKBs5ahD8a4OtKDtYaknEI0lUdDigKtQFMDzHuGOE1l+HKzvoYodSflFOJegMA3cFRfqnmjuyw4U0Q1tNNqMJwGMNbXgXuoqyFDvxpdZbRzlF5Hb93fwZU5VRW1gQewqiRanzevoe20HgfrVERV2ZinEvMeJ8GDU9Z26kvgyjXmaWSXBcqrf6WPIqoc772hr6UMbikudV9ionui/jj94Q4chhqevF73IBZRv8qjIRVBPu4FQNbdMm0poBn6xoyKz+UjcBhgK7ZV3xeEFqHbBZ1V1Bkvw5zHycP1PVJ1KqIaDhFAw4qoVA2Necpxw3mDLntX/Wsl15inUadIy1xucK3o4RQlD+GNTolkJXivDX3kMEEodlqxTV01eMe8PAk+IV3PCXoDWu18+jrsfDzeV3fLtIaklflcugP2zrq3pMvBGEPdDpzdtStu8pkdazBzOVx76YNizE07iZ+m+PZN7XH6dCatatw0ek9AqTzePIXTWM2XbuweVd2YpxyaIqruHIc6E9gKmRzs5FJfe0NrperGPOWY+j6qRpSx4Zv32tBrIYKrDU5eUFvSg8QBS5Gl2i+us8ZNMeY8Tr7cinCUqL2I1UpkhYG2AZztdSTl5GNzwORc3d7jUgNa7XzcvW6GO4br8x4zGV3CaaAahkN9dW+0pJyGd1mQzVxOqXLSNaJHOE3D6/KycbBB4DBQ+8VaY56hxkImKk4/VLciaiW6Un1jnnJ0DoLrY8MLnL3fhj5b77xRLwUarE8feqGaBvfUKWnMY94zRDItebIRqflaf1RptetW3OTjXlC5AbFgTZdlMpKVOmrcFEMIgXfMW1+cPrqhyi7XkalcyHBPG9dHe+qK0zec15DPZZ+qpVSHUVkK1NaYpxwNZS5rjXlsje32QIVv9g4TOTVRLTSU11CI+wfqkDypc2OWPN5vQx+I0dNmY6yvcjenSgx3DuPuc9d3IBuuvlNOJe5dGcBmETX3kZVS6hML1tBKsdYYvnkdPeYoka6pDG45vC4vuye7uRr7VaPjLgvUTuvxxl7Niih/VCXluPvcjQ/C3gET9SmiloO1NeYpx7WBa/S19dUns6yzxk0x5uqM09fVmKcc7gVIn8K2zg1z8nivDb2m1dbFg0UZlafBpyQzNXSal1LXydvVZuPW5f6aJ2/gMMBh8rBxFYHG2E1o66vZqJyFCBr3HOFsp1WzUdEUN8PXdRnHnMfJSTLDV1u1ZS77o37Gu8fptNcv/z2HewHefg3Hte34luqscVMMi7AwO1qHIkprzKPTw/dSfwdTzs6ad1qa06DLLgtUmLNBRVQl3ltDL6XMFmjSx3MEtSU9Sh3xfOd59RfF3sLJvm6TF9SW9JvtKAcn1T9wdFMRaFisKvW+xsm7pOO5CcBEzwTj3eO1H8iGXqp6JG36jOO+24kQtcfp/ZEa+pFWQ04R9ZdVX3KSTPNq90iXEKeGd8zL28O3bMe3q7+ojsY8lZjzDPFwbbcmRZRu6jSN9l64dLtl6I0gHDslepTUzUsBmB2tQ+Whg1a7kDnPEBkJj9aq17Frk1e3uCMooxJZh2j18sblQIxLfe30NpCUU4jX5eVx4DEZWYO8UatxoxN9nXY+vtRX005LS8rRzXMEJbG0d9ZkVNbCh6QzsmFpZT4+l4rT13SmpXM4DZRTFDtN8bwGRZQ/Ukdjnkq4F1TRuVP9ej/nU5WhF0L8ihBiSQjhF0L8TpH3f00I8bUQ4pkQ4okQ4vvVXtssloP6HcRq9Lf3c33geo2GvvHsy0JuT/bjsFlqitP7o35GOkboa+vTbRy5FO8a4vRLwbiuBgVU5vJB4oClvSoVUemk8h7rrKVSijmPk2ebUY4T1cXpNw42SMmUfrssyCqi7td0T/RU3Gi4+9wMdQzV7hRZHarqo05oiqha10pdjXnK0YAiqhoqGnohhBX4PeBXgRvAbwohCh+p/x64KaW8Bfxt4B/XcG1TWGqwQFMpvGNenoWecZqustN86CV0jUBXg5LGPNrtVu5dGajJe1yJrOgbIgDleXU6q/YeU+kMq+G4rg9fqKPuze6qaqCio+cIytAn0hm+eFVdfDyX16D3fXEvKKMZD1X18aVgTDXmcdZXW6YYQojaM5frbMxTjuGeNq6Ndte2VqINlgkpxmUfWOyGySyr8ei9gF9KuSalTAA/4f9v77zD4yjPvX2/u+q9W7YkW5JV3C0bN3BvlAChhGbAEFNcRT/JCUmu75/Tcs4JEIIxNjWUQCgBkgAHMGAbFzAY3LBlrVa2LDdp1bXq7f3+GI2ssn1ndom893VxYe2+Iz/enXnmmd/7FLim/wIpZZM8/21FAtLVYzWjrRHeXAkH/+LSclOFlcTIEJK8LMoZzOzU2bR3t3PAcsC1AzTciO3PJWMTKTrXSHWT8xtOd083JxpOaH/yGgxKVH98h0sj007WttDR1aOZPq8yInIEmTGZfHXOxY6aGhXlDGZmZgJBBsEus2s6fUldCUZh9L4oZzB9w6lduwGrg3lCgrRVemenzqa6tdr1HlFVnncSdcQlY5P49kQtbZ3On7QaOxqxtFi0v/mGRHjdI8oRrnxzaUB/kfV072sDEEJcJ4Q4BnyIEtW7fGzv8at7ZZ99VVVVrtg+kNBoZSDB9v9SHr2d4Mk8UleYmTqTmJAYXjv6mvPFPT3nUys15tKJSk7+n3aXOV17pukMbd1t2jt6gAk/BetZKPqb06WmCs+HjThj2Zhl7DmzxzWnYikCYYCkPE1tiAoNYkFeMm9+W05Tu/OCNnO9mTExYwgxuje2zykjC5QuqV9vUrK+nOBtJ1F7zE+fT7AhmNeKXLhW2hqh4ZQujv6yiam0dnbz5rfO95I0qx63Rd7lSi2NDv3pXXH0tnIPh5wdUsr3pJTjgGuBf3Pn2N7jn5VSzpBSzkhOTnbBrMFWCljyW8XZH3jd4VIp1aIc7TZiVSKCI/j5xJ+z/fR2DlUdcry4/qQy+UeHkzdvRDRXThnJi7tPUOMkqi+pV7ojapZa2Z8J1ypa97b/cnoCF1daERoV5Qzmjgl3EB4UztMHnna+2HIUEsZCsPf1FYN5YGkudS2dvLTrhNO1mtY19MdghIW/UDb/TB87XNrc3sWpWs8H8zgiJSKFG/Nu5H3z+5Q3OpnQ5uVgHkfMyU5gdlYCG7eZne6fqJ1ENd0gV5l7P9z6pvL9aIwrjv40kNHv53TgrL3FUsovgbFCiCR3j/Wa3EshbQbs+B/osu/cztS30tzRrUuUAnDb+NuID4137lQ06nFjj4eW5dLW2c2WLx0XC5nrlIwbjyflOMJghEWPQnUxHH7H4VJTpZUxCRGEh2h/oseHxXP7hNvZenIrx2qPOV7s4exeV5iaEcey8SN4budxGlrtP3m2drVy2npae4mgz5AVyqbmtv9wKKuV9FaPay2nqdwz+R6CDEFsPrjZ8cK+7DTtvxchBI9cmk+VtZ3Xvj7pcK3Xg3n8hCuO/lsgVwiRJYQIAW4B/t5/gRAiR/RWHQkhpgMhQI0rx2qKGtU3nobvX7G7zOTFPFJXiAiO4O7Jd7Pn7B7HY9NULVjj7A6VnJRori1I45WvyrA02i+v1rwoZzDjfwqpk2HH7xzKaqbKJt0cCsCdE+8kOiTa8Q24sw1qj+sSOao8vDyPxrYuXthp/wZ8vOE4EqlPRA9gDFZuwBWH4dg/7C7TU04DpaJ8xbgVfHjiQ47XOwhIvBzM44xZWQnMz03imR2lDmU1TduE+BCnjl5K2QUUAp8ARcBbUsojQoi1Qoi1vct+BvwghDiAkmVzs1Sweawe/5A+shfBmLnw5e/tTtMprtA3SgG4Kf8mksKT2Lh/o/2sAksRxGYoBRM6cf/SXDq7JZu229emzfVmfR5FVQwGWPwbxYEefMPmkvaubk5UN+vmUABiQmK4c8KdbD+1ncNVh20vqjaB7NHtKQtgwqgYrpw8khd2naC22XYfcvUpSzdHDzD5BiWLZdt/2pXVTJVWrwfzOGPVpFWEGkPZdHCT/UXqU5YXg3mc8cil+dQ2d/DynjK7a8x1Ol8rOuHSpyal/EhKmSelHCul/I/e1zZLKTf3/vm/pZQTpZQFUsqLpZS7HB2rK0IoTqWpAr59weYSU6WVkbFhxIZrV5QzmPCgcO6dfC/7KvfZ739j0aYNriMykyK58aJ0Xt9bztn6oTe+zu5OyhrK9JMIVPIuV4p17MhqehTl2OL2CbcTFxpnP6rXoSjHFg8tz6W1s5stX9q+AZvrzYQYQsiIzrD5viYYjLD4UUX//uFdm0uKNRjM44yEsARuH387n5R9Yr/WwQfXSkFGHMvGp7BlR6lNWa2mtYa69jp9b746MTwrYzPnQvZi2PWEzUqz4gqrrtG8yg15N5AamcpT+58aGtWrRTk6n7wAhUtykEie+sI85L2TjSe1L8qxhXoDbjhlU1bTW05TiQyO5K5Jd7H77G6+r7Qx7EEtyknQYb+iHzkp0VxTkMbLe8qwWIfKaiX1JWTHZRNk0C5n3Cbjr4ERk3qz1YZKFlq3CbHHnRPvJDrYjqzWN5hH/xKch1RZzcZmueatD3zI8HT0oGj1LdXwzZYBL3f3SMxVTbpkEQwmxBjC6imrOVR1iJ1nBlUi1h5XZnhqWBFrj/T4CG6ZOZq3952ivKZlwHua97hxxNglSgOnnY8NkdWKK6wEGQRZSdoV5djjlnG3kBiWyMYDG4e+aSmCxFxFw9aZB3pltWdsyGrmOp0ybgZjMMDiX0NtKRwaWINS39JBZWO7Ltlpg4kNjWXlxJVsO7VtaK8oi757Wf2ZOCqWKyal8uKuE9QNktV8eq1ozPB19OkzFLlg9x+VpmG9nKxppqOrxydRCsC1OdeSHpU+VKvXOeNmMIVLcjAaBH/8omTA6yX1vUU5sZn6G6FullvPwb6XBrxlqmwiKylS86IcW4QHhXPvlHv5tuLboV0tfSARqGQmRXLD9HT+vLeccw3nb3zWDiuVLZXa9h1yRP5PlKZaO/57wOxStU2I3nKaysrxK4kNjR16A/aRnKby0PI8mju6hmSrldSVEBca5/1gHj8wfB09KJFKWz18/UzfS30SgY9O3mBDMGunrqWotogvyr84/4alCBCatcF1xoiYMG6fM4Z3vz9NadV5Oau0vpTRMaMJNWpbIWyXzHnKoIVdj0PH+cHMJp2KcuxxQ94NjIgYMVBWa7dCQ7nPHD3AfUtzkFKysZ+sphbl+GzTTwhY/FuoL4f9r/a97Cs5TSUqJIpVE1ex68yugZXlVUUQFqfJYB5XyBsRzU+njuLlPWVUWc/vJ6l1Df9sGTcw3B39yKkw/mr46mloUTo5Flc06VaUY48rs68kMyaTjQc2nu+gaDmq6MDB4T6zY92isYQGGXnys/NRvW5FOY5Y8ltoroJvngWgpaOL8lpt2+A6I9QYyuopqzlYdZBdZ3pzB6p6NwJ9FDmCIqvdPDODt/ad4lStIqupBWw+1YJzlir9Vr78fd+kI1OldoN5XGXFuBUkhCWwcX+/qF7tJOpDB/vA0lw6unvYvEO56UoptZvd6weGt6MHWPRrJVLb8xSgnLwZ8RFEhOi8ydWPIEMQ6wvWY64380nZJ8qLPpQIVJKiQvn53Ez+cegsxRVW2rraKG8s9326WMYspbht95PQ1kiJKhH40NEDXJdzHWlRaWw80Cur6dAy2hUKF+cihOCPnysO3lxnJiIogpGRI31nRJ+sdha+U2S14gptB/O4QkRwBPdMvoe9FXuVFsYaD+ZxlezkKK6flsarX5+koqGNypZKmjqb9Kke9wHD39GPmACTfgZ7t0BTlc+yCAZzWeZl5MTlsOnAJrram5TNLx+fvACr52cTGRLEE1tNfUU5PtOC+7P418qUo6+f6esk6is5TSXYGMyaKWs4WnOUL0590VuUEwFxY3xqR2psGLfPHsO7+89worq57ynLIHx8eWYtUJrQ7Xwc2dHce6347slX5ab8m0gJT1FuwI1newfz+P5auX9pLj09kqe3mftaH/jlWtGA4e/oARb9Crpa6dr5uFKUk+r7k9cgDBQWFFLWWMZHP7yse1GOPeIjQ7h7XhYfH6lg+3GlF49f0sVGTYNxV8FXGzl15jShQQZG61iUY4+rx17NmJgxPH3gaXoqjyh7JjoW5dhj3aKxhBgNPPmZqa/60i8s+S00W2jetZm6lk6/BEWqrLbfsp89pveUF/1wrWQkKLLaX74t57tzytNeQLr5MZOUC1NXYNj3Agk9tX45eQGWjF7C+ITxPGN6k07wqRbcn7vnZxEbHsx7R74n2BDM6Gh9ysqdsliR1XJLXyYnJUrXohx7BBmCWDd1HSV1JXzaWOK37yQ5OpQ7L8nk7z+YqG2r9Z9DGT0HcpYRuvcpImn16b5Jf67PvZ5RkaN4qvRdpQuiD9KQbVG4RNl8/dh0UPvBPD7kwnD0AAt/CT3dbAh632+OXghB4bRCTnfU8beYGKVDoh+ICQtm9YJsTjcfZ2TEGP2LcuwxYiJMvI5lje8xPVH71qyucnnm5eTEZLEpTNLtoywoW6xZkE1EpNKr3q9FOYt/TXB7HauMH/s0E6o/wUYlW+1IexXbE0ZqOpjHHUbGhnPrrNGUNx0nLVK7yVa+5sJx9PGZHEq5mluM28gOcW3Cjx7MT5vPFBHOlvh4OvyYpfXzSzIJCrNgbfRvTrB1zi8Ile38rM1xZ0s9MRqMrE9fzomQYD4SLc4P0In4yBDmjFNK73vaR/jNDtIu4kjMPNYEf0iS0X+fx9Vjr2Z0j4GNsZHuzfvVmDWLMjGEWKipS/CbDd5y4Th64PXQmxACQnc/5jcbhBAUNrZSYZC8Y/Kfc5OiDYLqqKiOY68b8zK15lh3Ku/3zGPy2bfAWuE3O5b2hDKuvYNNZ7bR2eN8cI1ejEiqg+4IXvrSf98JwHNBK4imRUlN9hNBGFhb34CJDrae3Oo3OzqpQRg6MZ2OoKQ3ceCfjQvK0e+tCWdP7NWw/zWlBYE/aG9iTnU5M0JH8Nzh52jtst1hU2/Ucu5oQzqPbTW5PrdTY0yVVp7suh6D7Iadj/vFBgBD1TEKmzo53XKOv5v166TtjFNNJxgZkcnnRVUcOFXvFxuklHxWk8wPcUuUYsNmP9106k/yk4Y6skMT2HRgE906TF5yBbWuIaR7FH/4rMTJ6h8nF4yjV4tySvLXKE2rdvyPfwypKkYAhdnXUN1azVvFb/nFDLX68o4Zc/jmRC27zf65mE0VVmqCR0HBbUr+dr3zcW66UHWMBbE5TEmawpZDW+jott0+WE+klJjrzFycMZGEyBAe+9ROJ0edOdvQRlN7Fycm369MQNv9B7/YQdUxjMD63Js53nCcj0585Bcz1Gtl5UWz+PDwOY6ebfSLHd5wwTh6s6UJKSE9IxNm3QOH3oQqk+8N6S3KuWjsT7h45MW8cPgFWjp9r4Oa682EB4Wz+pKLGBUbxu8/LfZLVK/OIxULfqG8sPP3PrdBLcoRKRPYMG0D55rP8deSv/rcjMqWSqydVsYn5rF2YTY7S6r55kStz+1Qh42MHDsVJt8I3zwH1kqf26FeK8sn3Ep+fD6bD272i6xmrlMG86xbMJHosCCe+MwPfsNLLhhHP6BB09wHlcKY7f/le0MsRRAUDvGZFE4rpK69jj8X/dnnZpTUl5ATl0N4cDD3Lc3lwKl6thVbfGqDlJLiCquSwheXARf9vFdWcz5PVVOaKpXirZQJXDzyYqanTOe5Q8/R1mV/Kpce9O+OuHJOJsnRoTzmhxuwWsCWOyIaFv6r0mV1lx9kNUsRxI7GEB7HhoINlFvL+Uep/WlYelFSX0JuXC6xEcGsnp/N1qOVHPSTrOYpF5CjtxISZGBMQgREJsHstXDkXaj4wbeGWI5Cch4YjExJnsLC9IX86cifsHb4dpPHXHe+KOeGi9IZnRDBY5/6VquvbuoYWJQz/xEwBPleVus3j1RNga1qreLN4jd9akb/qVLhIUY2LBrL3hO17Cn1raxmqrCSGtM7mCdxLBTcCvtehIYzPrWj/+zeRRmLmJQ4iS0Ht9DpYByl1qiDedRrZdW8LOIjgnl86z9XVH/BOPriCis5yVEEGXv/yZcUQmis76P6qmMDinI2FGygsaORV4++6uAgbalrq6OmraavKCfYaOCBpbkcOdvIJ0d8l/miZjD0OfroVJh5j9IXvdqHm16D2uDOTJ3JnJFzePGHF30Cg4j5AAAgAElEQVQqq5nrzSSFJxEXFgfALbNGMzI2zOdRffHgTqILf6nIW76U1bq7BgzmUW/AZ5vP8m6J7WlYetA3mKe3riEqNIi1C8eyw1TFvjLfy2qecsE4elOldWAvlfB4xdkf+wDO7veNES21Si/2fuXc4xPHs3zMcl45+gr1bb55HFQlgv7NzK6dlkZ2ciSPbzXR3eMbp6JKBHn9W1LMfVCRtnx5A7YUQWSy8qTXS+G0Qmrbann92Os+M2NwJ9GwYCP3Lcnl+/J6thdX+cSG7h6J2TJoME/caLjoTmUyWF2ZT+zoG8zTLyi6ZNQlTEuZxrOHnvWZrGZuGHqt3HFxJklRoTz26T9PVH9BOPqG1k7ONbSRO7hB0+y1isPf9p++MaTqmPL/QWX266eup6WzhT8d+ZNPzLA1Es1oEDy0LA9TZRMfHDrrEztMlVbiI4JJjurXCz8qGWavUWaYVuo7R74PG51EpyZPZUH6Al764SWfyGo9ssdmG9wbZ6STkRDOY1t9E9WX17bQbmswz/xHQBhhx//qbgNgs5OoEIL7pt2HpdXC26a3fWKGuc48ZDBPeIiR9YvG8tXxGvaYq31ih7dcEI7ebLEzQCEsBuY+ACWfwqlv9DfEThvcnPgcrsi6gtePvU51q/4njrnOTExIDMnhyQNev3LySMalRvPkZyV0detfiVhcoXQSHdIG95L7IDTaNzfgnp4hcpqKKqu9dvQ13c04Yz1DW3fbkDa4iqyWxw9nGvnkiP6ZL8UVdjqJxoxSZLWDr0P10NnDmmMpAmGApLwBL89Mncns1Nk8f/h5n8hq5nqzzcE8t84eTWpMmF9rUNzhgnD0xRUO+p3PWq08tn/x7/obYimC0BiISRvy1rqp62jvbufFH17U3Qx7k3IMBsGDy/I4Xt3M+wf0jeqllJRUNtn+TiIS4OINvbLagaHva0nDKehosjmPdELiBJaOXsorR1+hob3BxsHaoRbl2OpaeW3BKLKTInliq4kenWU1daqUzcE88x6CoDDY8TtdbQCUoCg+y+ZgHlVWe+PYG7qbYW8wT1iwkcIlOXx3so4dJt/Iat5wQTh6U6WVyBAjaXE2pjmFRMK8h+HEDjixc+j7WmI5pjgUG4McMmMzuTr7at4qfovKZv0iNyllX2qlLS6bOIJJaTE8+bmJTh2j+nMNbVjbu+w3zZqzThkfp3dU72Qe6fqC9TR3Nusuq6lFOWNjhzr6IKOBB5blUlxp5cPD53S1o7jSyugEO4N5opKVwOjwO+c/N71wMJinIKWAeWnzeOnISzR1NNlcowXqYB5718pNMzJIjw/n8X+CqP6CcPTFFVZyR0RjsNcGd8ZdED0Ktv2Hkl2gBy5Mylk7dS3dPd08d/g5fWwAqlqrsHZY7XZHFELwyPJ8TtW28va+07rZUexsHmlYbK+s9om+slqV6uiHRvQAefF5XJ55OX8u+jM1rfqlOZbUlzAqchRRIbZnJVw9ZRT5I6J54jOTrrJaibPBPHMfgJAofW/AnW3KZqyDltGFBYU0tDfwapF+2WonGk4gkXYdfUiQgfuX5nLodANbj/qhoMwNXHL0QojLhRDFQgizEOJXNt6/TQhxqPe/PUKIqf3eKxNCHBZCHBBC7NPSeFdxOiknOAwWPALlX0HpF/bXeUOTBVprHZ686dHpXJd7HX8t+Stnm/SRTvrnattjUX4y00fH8dQXJbR16tNfRK2+dPi9zFoNEUnKDVgvLEUQk67cWOywrkB/Wc1cb3bYmthgEDy0PJfjVc38TSdZraOrh+NVTgbzRCTAxeuh6O9w7pAudlBTArLbYVA0MWkiSzKW8OqRV3WT1WwlLQzm+mlpZCUp2Wp6y2re4NTRCyGMwNPAFcAEYIUQYrC3OgEslFJOAf4NeHbQ+4ullAVSyhka2OwW1U3t1DR3OO9BP+0OiB2taPV6RPX9inIcsXrKagwYePbQ4I9QG/oGTztw9EIIHrk0n3MNbfzlm3Jd7CiutDIiJpS4iBD7i0KjYP7DcHw7lO3SxQ7lKcvxd5IVm8VV2VfxZvGbWFq0rx7u7OnkRMMJp8NGLpuYysRRMTz5eYkustqJ6ma6eqTza2XOen1ltT45zfGwkfUF67F2Wnn5yMu6mFFSX+J0ME+Q0cCDy3I5VmHlox/0ldW8wZWIfhZgllIel1J2AH8Brum/QEq5R0qpNnn/GkjX1kzPMbk6jzQoRCkMOfs9mD7W3hAnWrBKamQqN+bfyPvm9ylv1N7JmuvNJIYlEh8W73DdJWMTmZ2VwNPbS2nt0D6qd3l274y7IHokfKGDrNbdpfQ7cmFMnSqrPX/4eW1tAMoby+nq6XLq6JUbcB7ltS288532slrx4AI2e4THKZlRpv+D0zo8pFuOgiHY6WCe/IR8Lsu8jNeKXqO2TfviJXOdmazYLKeDea6aMorclCie8GENiru44ujTgP4tBU/3vmaPu4H/6/ezBD4VQnwnhFht7yAhxGohxD4hxL6qKu12sVWJwKWRaFNXQEK24lR6NI6YqoogIlHJ8HHCPZPvIdgQzOaDm7W1AeXkdWV6kRrVV1nbefXrMk1tUItyXHL0weFKDnf5Hji+TVM7qDsB3e0ujQ/MiM7g2txrecf0DueatI3c+ve4ccbi/BQKMuJ46vMS2ru0vQGbKqwYDYLs5Ejni2evVc5nPWQ1S5Ey/jPIwdNeL+unrqe9u52XfnhJczPsZdwMxmgQPLQ8j9KqZv52wMdtIlzEFUdvawfT5m1LCLEYxdH/a7+X50opp6NIPxuEEAtsHSulfFZKOUNKOSM52bkzdJXiyibiIoJJjg51vtgYBIsehcrDigapJZYixaHYyLgZTFJ4EivGreCD4x9wvF67vvk9sofShtIBVX6OmJWVwPzcJDbvOE5Te5dmdpyqbaGts8f1eaTT74DYDO1lNRclApU1U9YAsOXQFu1sQHEoBmEgK9b5qDohBP9yaT5nG9r4yzfatnQ2VVrJSookNMjofHFolJJuWfoFnNyjqR2OMm4Gkx2XzZVZV/LGsTeoatEuQGzqaOJc87khdQ32uHxiKhNG6iereYsrjv40kNHv53RgyG6QEGIK8DxwjZSyLz1BSnm29/8W4D0UKchnmCqt5KXYKMqxx6SfKSmQ2/8LtBp0IKVbJy/AqkmrCA8KZ9PBTdrYAJxtOktrV6tbg6cfuTSf2uYO/rRbu46S51sfuOjog0JhwS/gzHdg+kQzOxRHLyDJtTmxqZGp3JinyGqnGrVzsuY6M6OjRxMWFObS+rk5iczKSuDpbWZNN8tNlVb3hoHPuBuiRmgrq7U3Qf1Jt66VdVPX0dXTpamsVtqgpLu6eq0YDIKHl+dxsqaFd7/XL1vNU1xx9N8CuUKILCFECHALMCDcFUKMBt4FVkopTf1ejxRCRKt/Bi4FfNYuUkqpOHpHWQSDMRhh0a+UaskfNOpJ7qAoxx7xYfHcNv42Pin7hOJabQZQqBKBraIcexRkxLFsfArPfnmchlZtugaqclquraIcu4bcCvGZilSglaxmOar8zpAIlw+5Z/I9BBmC2HxIO1nNVYlARUmBzcNibee1r09qYkNrRzcna1uGtglxREiEIqud3KXUoWhBVe+5nuy6o8+IyeDanGt52/Q2Fc3aNOVTs9PcuVaWjk9hakYcf/zcrLms5i1OHb2UsgsoBD4BioC3pJRHhBBrhRBre5f9PyAR2DQojXIEsEsIcRD4BvhQSqnDTqdtKhrbsLZ1uRelAIy/BkZMhu2/UzbsvMXFjdjB3DnxTqKDo3n6gDZzO93Rgvvz0PI8Gtu6eGGXNlF9caWV9PhwIkMdb3INwBisyGoVh+CYRj3JVTnNDZIjkrkl/xZFVmvwXlZr726n3FrulkMBmJ2dyPzcJDZtL6VZA1lNHczj9rVy0c+V9FStZDU7bUKcsXrKaiRSM1lNHcyTFuVoO3Ig6g34TH0rb33rp0lpdnApj15K+ZGUMk9KOVZK+R+9r22WUm7u/fM9Usr43hTKvjTK3kydqb3/TVSP9RXFFS5mEQzGYIDFv4baUqVlrrdYHBfl2CM2NJY7Jt7BtlPbOFLtfYMvc72ZkZEj7Rbl2GPiqFh+MjmVF3edoK7Z+xF7JZVN7jsUUKYdJeXBNg1kta525ft106EA3DX5LkKNoTxz4BnvbEApyumRPS5tkA/m4eV5iqy2p8xrO0zuymkqQaGw4F/g9LdQosEA76pjfYN53GFU1ChuyL2B90ve55TVeydbUl/C2NixGIR7NaXzc5OYmRnPRo1lNW8Z1pWxJlfTxWyRfwWMmg7b/xu6vHRuliKl8jbccUqjLW4ffztxoXE8deAp72ygN+PGzWhe5cFleTR3dLHlS++i2I6uHkqrmtx3KNArqz2qZDD94GVP8hoz9HR55OgTwhK4ffztfFz2MaY671rVltQpdQ2ubpD3Z9roeJaOU2S1xjbvZLUBg3ncNuT2XllNg6jechSS85Xv2k3unXIvRoNRk2y10vpSj26+arZaZaN2spoWDGtHX1zRRHJ0KPGRztO0hiAELPkNNJTDfi/LrJ20PnBEVEgUqyatYveZ3RyweN7gq6uni+MNxz06eUG5WV4zdRQv7ymjytrusR1lNUpRjkcRPcCEa2HEJGWz3BtZzUM5TUWV1TYd8G6z3FxvJsgQxOgY+0U5jnhoeR4NrZ28sNM7Wa24ctBgHncwBisjB88dVBrReYMHcppKSkQKN+ffzAfHP+BEg+efR11bHdWt1R4HRXOyE5mbk8jmHdrIalowrB19icXNLILBjF0KGXPgy98r/Tc8oadb2WDy0NED3JJ/CwlhCWzcv9Hj31FuLaezp9PjkxfggWV5dHT38Mz2Uo9/h8dymorBoET1taXKgHdPsRxVxhYmevZ5xIbGsnLiSj4v/5wjNZ7LauZ6pSgn2BDs0fGT0mK5fKL3spqpwkmbEGdMvkn5LLf9p+eb5X2DedyTOPtz16ReWe2g57Kap3tZ/Xl4eT7VTR28/FWZx79DS4ato+/pka5XX9pDCFjyW7Cehe88LMioVYtyPHf0EcER3Dv5XvZW7OWbc541+HKlx40zspIiuX5aGq/tPUlFg2c3PlOlFYPAtaIce4y7EkYWKO1yPZXVLEWKY3KhKMceK8evJDY0lqf3e75Zbq4zkxPr+XcCSlTf1NHFszs9k9Ua2zo529DmmZymotagWI4qs5g9wc5gHndIDE/k1nG38vGJj/tkMXfRwtFfNCaexfnJPPvlcaxeympaMGwd/am63qIcd1IrbZE1H7IWwM7HoaPZ/eOr3CvKsceN+TeSEpHCxgMbPWqJWlpfikCQHZvtlR33L81FSsnGbZ5dRKZKK5lJkYQFu6/B9qHegOvL4YCHQ0HcrGuwRVRIFKsmrmLnmZ0eyWrNnc2cbT7rsZymkp8azdVTRvGn3WVUN7kvq5VUKq1+vXr6BZh4veKkPc1Wc7OAzR6rJq0iMjjSY1mttL6U6JBoUiJSvLLj4eX51Ld08uKuMq9+jxYMW0fvtUTQn8W/hWYLfONB+2D15HUjh94WocZQ1kxZw37Lfnaf3e328SX1JYyOcb0oxx4ZCRHcNCODN789xek69yf8mDzNuBlMzjLImK2MtnNXVutoVmafehE5qqwYt0KR1Q64L6upPei9iRxVHlyWS3tXt0eymldJC/1Rs9VqSuDwW+4f72AwjzvEhsaycsJKPiv/jKM1R90+vqSuhNy4XNeLLO0wOT2WyyaO4Pmdx6lv8T5bzRuGraNXT95cLZzK6NmQsxx2Pwltje4d21eU44VU0ct1OdeRFpXGxv3uR/XuFuU4onCJMp3qqc/dGynX1tlNWU2zNjdfIWDxb3pltT+5d2xVMSC9jhxBkdXumXwPe8/t5duKb9061taQdk/JTo7i+unpvPb1SSob3bvxFVc4GMzjLuOugpFTYcd/Q7ebkoX6lOWlgwVYOWElMSExbtegSCk1vVZUWe05D2U1rRjGjr6JtLhwotwpynHE4l8r/eT3upm6ZSlyq8rPEcHGYNZMWcORmiNsP7Xd5ePau9spb3S/KMceI2PDuW32aN75/jRl1a7LWX1FOd5owf3JXgiZ82HnY9DhxtNF31OWNt/LTfk3kRKe4vYNuKSuhDBjGGnR3kWwKg8szaW7R/L0NvduwKZKKzmOBvO4g3oDriuDA392/Th1MI+XT74q0SHRrJq0ii9Pf8nBqoMuH1fVWkVjR6Nm18q41BiunDySl3aXUeOBrKYVw9jRW7VzKABp05VoZc9GaK1zvh6UjcIasyaRo8rVY69mdPRoNh7YSI90LbuhrKGMbtmtSeSosm7RWIKNgic/d12rL3Zl2Ii7LP6NIqt964asZjkKxlBIcN5EzBVCjaHcO+Vevrd8z1dnv3L5OHO9mey4bLeLcuyRkRDBTTMzeOObcrdkNaXHjYbfSe6lkD5TkdW6XHRuLgzmcZdbx93qdraamrTgajMzV3hwWR5tnd1s3uF5tpq3DEtH39ndW5SjhUTQn0WPQnsDfOXi42BfUY52J2+QIYh1Besw1ZnYetK1SkQtsggGkxIdxp2XZPL+gTOU9MpkzjBZrIQYDYxJ9F7G6mPMxUoa7K4/QLtrdlB1zOOiHHtcn3s9oyJHubVZXlpfqul3AlC4OAeBYOMXrkX1NU3tVDe5MJjHHdSovvE0fOfiUBCNkhb6ExEcwV2T7uLrc1+7LKt50g/KGTkpUVw7LY1XvjqJxU1ZTSuGpaMvq26ms1t6n3EzmNRJMPE6+PoZaHZhfqiHfTuccUXmFYyNHcumA5vodqEVgFqUMyZmjKZ2rFkwlohgI3/4zLWo3lRhJTs5kmBPinIcsfg3SjT4tYuymhdFOfYIMYawZuoaDlcfZsdp5w2+6tvqqWqt0vQpC2BUXDi3zh7N29+5JquZ1IwbLZ9+AbIXwZi5sPP3rslqXhaw2ePm/JtJDk92WVZTB/MkhCVoaoensppWDEtHr7bBzU3R+OQFJarvbIHdf3C+1lIEwqgMUdAQo8HI+oL1HG84zkcnPnK63lxnJjMmk2CjZ0U59kiIDOHueVl8ePgcR88636Q2VTZp71AA0i+C/J/AV09Ba73jta310HhG85svnJfVnj7wtFNZzZV5pJ6yvldW+6MLsprJ2ZB2T1Gj+qZK2PeC8/WWo8p84CjtZlEAhAWFnZfVzjmX1ZzN7vWUMYmR3DgjnTe+OcWZ+lbNf78zhqWjN1U2YRDKI5PmJOcrVYDfPAdWJ5PfLUWQOFZp/KQxy8YsIz8+n80HN9PZ4zi7oaS+RHOJQOXu+dnEhAXx+FbHPV+sbZ2cqW/VXk5TWfxraHNBVusrytHe0Qcbglk7dS3Hao/x2cnPHK7VQ05TSYkJ446LFVnNbHEsZxVXWokNd3Ewj7tkzoXsxbDrCaXPvCM0qGuwx89yf0ZqZCpP73/aYVTfI3s0zbgZTOESJeDb+IVnNSjeMDwdfYWVzEQvi3IcsehfobsDdj3ueJ0XPW6cYRAGCqcVUm4t5x+l9tv2tnS2cKbpjG4nb2x4MPfOz+azokoOnrIfTasSgW6OPnWy0gfn602OZTWd5DSVn2T9hOzYbKeymrneTFRwFCMiRuhix5oF2YQHG3nCiaxmqlDahHibM26XJb+FlhrH2WpSguWYbt9JiDGENVPWcKj6EF+e/tLuOk8G87hDWlw4K2Zl8Pa+05TXuF+D4g3D09F72/rAGQnZMO022PciNNiZJtPRollRjj0Wpi9kctJkNh/cTEe37YIMtWe6Ho+jKqvmZREfEewwqi/RSyLoz6JHlWKoPU/aX2M5BiFRymhCHVBltdKGUj4usz96QY0c9XKwiVGhrJqbxYeHzlF0zras5tFgHndJnwF5l8OeP9qX1RpOQ4dVN0cPcE3ONaRHpTuU1bQsYLPHhsU5GA3uZatpwbBz9H1FOXpowf1Z8AslEvny97bfr9auKMceQgg2FGzgXPM53i2x3V/Emza4rhIVGsTahWPZYapiX1mtzTXFlVbCg42kx2tQlGOPlHFKz/q9z9qX1dSnLL0iWGD5mOXkxefxzMFn6OoZ2gqgryhHx5svwL3zs4l2IKtVNrbT6MlgHndRZbWv7bQk0Gkjtj/BhmDWFayjqLaIz8s/t7mmpF65VvR09IqsNob39p/GbHEiZ2nIsHP0ZksTPVLjXG1bxI1Wpuvsf1WJ3Afjg5MX4JJRlzA9ZTrPHXqOtq6hqVvmerNSlOPGpBxPuOPiTJKiQnnsU9tORXnKitKmKMcRi37VK6s9Yft9HbVgFYMwUFhQyMnGkzZlterWahraG3R1KACxEYqstvVoJYdOD42mi7VqfeCMkVNh/E/hq01Kh8rBqHKaRsVS9rgy60qyYrPsymqeDuZxl7ULxxIWbPRpVD/sHH2JxQcSgcr8R5RWtzv+d+h7lqNgDIF4bYpy7CGEoHBaIZZWC2+b3h7yvtoG16hhzrgtwkOMbFg8lq+O17DHXD3k/eKKJm3aUTgjcSwUrOiV1c4MfK+pClqqNauIdcSijEVMSpzElkNb6BzUCsAXkaPKqrmZxNmR1Uxa9oNyxuJfK3OT9/xx6Ht9g3nidDXBaDCyfup6zPVmPikbOmTeXGfWNH/eHolRofz8kkz+cfAsxyrcbKniIcPO0RdXNBFsFGQmaViUY4+YkTDzHjj4OlQPyo+1FEFSvtK+VWdmps5k9sjZPH/4eVo6B27ymOvMmlb5OWLFrNGMjA3jsa2mAdkNtc0dVDe1++bmC7DglyB7lBzu/ui8EdsfIQQbpm3gTNMZ3jO/N+A9LVpGu0p0WDBrF45le3EV350cGE0XV1o9H8zjLinjYfINsHeLUgXbHx2TFgZzaeal5MbnsungpgGymjqYR0+Jsz+rF2QTHRrEE06y1bRi2Dl6U6WVsclR2hfl2GPug8qMyx2/G/i6jlkEtigsKKS2rZY3jr3R91pDewOWVotPHApAWLCRwiU5fHeyjh2mqr7XPZ5H6inxY2D6HfD9q1DXb5ybBv3O3WHuqLlMS5nGlkNbaO8+3wqgtKGUhLAEEsMTfWLHHRePISkqZIisVlLp5WAed1n4K+hqU6qYVXq6odrks2vFIAxsKNjAycaTfHD8/DSsU9ZTymAenfdNVOIiQrh7fhafHKnk8OkG3f++Yefoiyt0zrgZTFQyzF4Nh9+Byt6Isa1BKf/2oaMvSClgXto8XjryEk0dyiaPL7IIBnPjRRmkx4fzeL+oXreiHEcs+BcQBvjyf86/ZjkK4QkQ5V2fcVcRQlBYUIilxcLbxedlNW9m93pCREgQ6xblsKe0hj2liqymDObRoU2II5JyYOoKpYCq8azyWl2Z4vx9dPMFWJKxhAmJE5QalF5ZTc+6BnvcNS+rV1Yr1v3vGlaOvqm9q7coR+eN2MFccj+ERitzTEGJ5sGnJy9A4bRCGtobeLVImXHb1wbXR9INQEiQgQeW5nLodANbjyqZL8UVVmLCghgRo0NRjj1iRsHMu+HAG1DT20xKbX2gY8bNYGaNnMXsVEVWa+1q1b0oxx63zR5NakwYj3+q3IBP17XS2tmtfZsQZyz8pdL/aedjys8+lNNU1Btwf1nNXGfWZDCPO8SEBbN6QTbbiqv47qSLjRI9ZFg5+hJfZREMJiIB5qyHor8rA5L7Tl59swgGMzFxIksylvDKkVdoaG+gpK5E16Ice1w3LY3spEge32oaMNJRt6Ice8x7SKlK3v673qKcIp9/J6DcgGvaavjLsb9wrvkcLV0tPtn0609YsJENS3LYd7KOL0uqz7cJ8fW1Ep8J01Yqzc7qy/u1jM73qRnz0uYxNXkqzx56lvbudkrqS8iIzvB6MI+73HlxJomRIbpH9cPK0fdJBL7Sgvtz8XoIi1OGI1uKIDgSYkf73Iz1Betp6mzi5SMvY65Xsgh87WCDjAYeWJbLsQorH/1wTpHT/PGdRKXArHvh8NtQ+jm0N/o0clQpSClgbtpcXvzhRQ5alN7ovnzKUrl5RgZpceE89mkxxb3ZHrl6tAlxxoJ/UZ6qdvyPpoN53EHNVqtsqeQd0zt+ecoCiAwNYt2isew21/D1cRcaJXqIS45eCHG5EKJYCGEWQvzKxvu3CSEO9f63Rwgx1dVjtaS4oomwYAMZ8RF6/jW2CYuFufeD6WMo+ocSORp8fx/NT8jn8szLea3oNYpri/1y8gJcNWUUeSOi+PcPinxTlGOPuQ8qlbB/u0/52cdymsp9BfdR317PY98pkoWvI3oYKKu99nU5aXHhRIdp2+jOJWLTYcZdcOB1KNvtt+9kdupsZqbO5LlDz1HeWO6zjdjB3D5nDCNiQvtkNT1w6omEEEbgaeAKYAKwQggx+Js5ASyUUk4B/g141o1jNcNUaSU3RaNJOZ4wa43Sgc961i+Ro8q6gnW0d7dj7bT6JXIEMBoEDy3Lo6K3/7bP5TSViASYs075TkD3ohx7TEyayOKMxVhaLIyIGEFMSIxf7Lh+ehqZiRFUNLb558lXZd7DSp1Js8Vv14qq1de01Wg+mMcdwoKNbFicwzdltewsGVqDogWuhJyzALOU8riUsgP4C3BN/wVSyj1SSnU34Wsg3dVjtaRY7x43zgiNUnRh8FuUApAdm81V2VcBvs0iGMxlE1OZOEpxaD7fIO/PxRuUJ67okYrj9xMbCjYA+vYdckaQ0cCDy/IAP958AaJHKLIa+PVamT5iOnNHzQX8e63cPLNXVtuqT1TvSjVPGnCq38+ngdkO1t8N/J+7xwohVgOrAUaPdl/b7uzuYUFuMvNzk9w+VlNm3g1NFUonRT/ywPQHiA2NpSClwG82GAyC310/hW3FFhKjfJhxM5jwOLhmk5L26kfyE/J5dNajZMXqWy3tjKunjqLEYuW6afq2xXDK/IeVDJycZX4141ezfsV75vf8+r2EBhl5aHkeB0/V097Vo3nnXeHs7m2hh3UAAAPbSURBVCGEuBG4TEp5T+/PK4FZUsr7bKxdDGwC5kkpa9w5tj8zZsyQ+/bt8+gfFCBAgAAXIkKI76SUM2y950pEfxro39M1HThr4y+ZAjwPXCGlrHHn2AABAgQIoB+uaPTfArlCiCwhRAhwC/D3/guEEKOBd4GVUkqTO8cGCBAgQAB9cRrRSym7hBCFwCeAEXhRSnlECLG29/3NwP8DEoFNvTnbXVLKGfaO1enfEiBAgAABbOBUo/cHAY0+QIAAAdzDkUY/rCpjAwQIECDAUAKOPkCAAAGGOQFHHyBAgADDnICjDxAgQIBhzo9yM1YIUQWcdLrQNkmAPg0j/vkIfBYDCXweAwl8HucZDp/FGCllsq03fpSO3huEEPvs7TxfaAQ+i4EEPo+BBD6P8wz3zyIg3QQIECDAMCfg6AMECBBgmDMcHf2z/jbgR0TgsxhI4PMYSODzOM+w/iyGnUYfIECAAAEGMhwj+gABAgQI0I+Aow8QIECAYc6wcfS+HEL+Y0cIkSGE2CaEKBJCHBFCPOBvm/yNEMIohNgvhPjA37b4GyFEnBDiHSHEsd5z5GJ/2+RPhBAP9V4nPwgh3hBChPnbJq0ZFo7e10PI/wnoAh6RUo4H5gAbLvDPA+ABoMjfRvxIeBL4WEo5DpjKBfy5CCHSgPuBGVLKSSjt1G/xr1XaMywcPT4eQv5jR0p5Tkr5fe+frSgXsp8HhPoPIUQ6cCXKBLQLGiFEDLAAeAFAStkhpaz3r1V+JwgIF0IEAREMwyl4w8XR2xpCfsE6tv4IITKBacBe/1riV/4A/BLo8bchPwKygSrgpV4p63khRKS/jfIXUsozwO+BcuAc0CCl/NS/VmnPcHH0wsZrF3zeqBAiCvgr8KCUstHf9vgDIcRVgEVK+Z2/bfmREARMB56RUk4DmoELdk9LCBGP8vSfBYwCIoUQt/vXKu0ZLo4+MIR8EEKIYBQn/2cp5bv+tsePzAV+KoQoQ5H0lgghXvOvSX7lNHBaSqk+4b2D4vgvVJYBJ6SUVVLKTpTZ15f42SbNGS6OPjCEvB9CGdz7AlAkpXzc3/b4Eynlo1LKdCllJsp58YWUcthFbK4ipawATgkh8ntfWgoc9aNJ/qYcmCOEiOi9bpYyDDennQ4H/2cgMIR8CHOBlcBhIcSB3td+LaX8yI82BfjxcB/w596g6Diwys/2+A0p5V4hxDvA9yjZavsZhu0QAi0QAgQIEGCYM1ykmwABAgQIYIeAow8QIECAYU7A0QcIECDAMCfg6AMECBBgmBNw9AECBAgwzAk4+gABAgQY5gQcfYAAAQIMc/4/SBjEfSfAyHkAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.plot(res)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Compute Stationary Distribution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "n = np.shape(P)\n",
    "w, v = np.linalg.eig(P.T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.4, 1. ])"
      ]
     },
     "execution_count": 67,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.70710678, -0.19611614],\n",
       "       [ 0.70710678, -0.98058068]])"
      ]
     },
     "execution_count": 68,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 77,
   "metadata": {},
   "outputs": [],
   "source": [
    "tmp = v[:, np.isclose(w, 1)][:, 0]\n",
    "sta_dist = tmp / np.sum(tmp)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 81,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.16666667, 0.83333333])"
      ]
     },
     "execution_count": 81,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sta_dist"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}