discrete_markov_chain.ipynb

{
 "cells": [
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "This code is based on code from IPython Cookbook 2nd edition, chapter 13.1\n",
    "ISBN:9781785888632"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f3fa7a56610>]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeYAAAD4CAYAAAA5OEWQAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/d3fzzAAAACXBIWXMAAAsTAAALEwEAmpwYAABCQElEQVR4nO3dd3wUZf4H8M+T3hMSEgIhEELvLRSpIoggKup5yp2e/bDg2e7O09PzbJyc/rw7vTsLZ8WGFewUUSlKC0ivIQkhCSU9kF6e3x87szuzO9uSze4m+bxfL17MPDO7OxmWfOdp30dIKUFERET+IcDXF0BEREQWDMxERER+hIGZiIjIjzAwExER+REGZiIiIj8S5OsLAICuXbvKtLQ0X18GERGR1+zYsaNYSploXe4XgTktLQ2ZmZm+vgwiIiKvEUIcNypnUzYREZEfYWAmIiLyIwzMREREfoSBmYiIyI8wMBMREfkRBmYiIiI/wsBMRETkR/xiHjMREXVu72w5jkdW7sOg5Gj07xaNayf0wsT0BF9flk+wxkxERD73yMp9AIBDp87ii92FWLB0i4+vyHcYmImIyCe2ZpfgtU05aGxq9vWl+BU2ZRMRkU9co9SKe8dHGB4/WVGD7rHh3rwkv8AaMxER+dT+wkrD8tOVdV6+Ev/AwExERF53trbBvP3Pb48YnvPCuqPeuhy/wsBMREReV1pV7/Sc7w6d8cKV+B8GZiIi8rraBg74soeBmYiIvK6mocnpOVP7d/XClfgfBmYiIvK6/LJqm7JuMaG6/Y1Hi71yLW9vzkXag1/hkZV7vfJ5zjAwExGR19XU29aYb5iUZlNW39j2Td5/+Ww/AOCdLXlt/lmuYGAmIqIWWbXvJP6x5jCklOayn44V4+lvDqLBSdKQP368x6asT0Ikdj86G89eNcJctjm7xHMX3E4wwQgREbXI7e/sBABMH5iIsb3jAQALl+3AubpGDOwWjSvH9HTr/SSA2IhgDO0Ray5b9lMupg9I9Ng1O1NaVY/4yBCvfZ4R1piJiKhVquoszdLn6hoBAIdPnbV7vraGrS83/T2kRwyuzjAF9XovpOvM6N3FvF1a5fukJgzMRETktopqS4IQIWyP55XaDu5S1Wn6jXt2CTe/PjHaMvhLrW3XeaGPWfuYsPN4eZt/njNsyiYiIrdlFZ0zb2cXVWFq/0TUaqZABQUa1/uklMgprjLvr753GvbkV+DgyUqMS7PUXEODTK93Nvirpr4Jh05VYlRqHITRE4ILCstrzNtNdmrz3uRSjVkIkSuE2CuE2CWEyFTKRgkhtqhlQojxSrkQQrwghMgSQuwRQoxpyx+AiIi8q66xCb946Sfz/l8/34/Dp84i46lvzWVf7C40fO2TXx7E3Oc3mvcjQ4NwXt8E3Dyljy6whgYFKp/lODAPfnQVrnjxJ7ywLqtFP8vZ2gacrKg173tjFLgz7jRlz5BSjpJSZij7zwB4XEo5CsCjyj4AzAXQX/mzEMBLHrpWIiLyA0bpNLPOnDP3LwNA99gww9e+/mOOS58RotSY6xqdJyIB7OfbdsZ6oYz2FpitSQAxynYsAPXxaD6AZdJkC4A4IUT3VnwOERH5kUMGA7t+/9Eu3b6z6VLOqE3Z2UVVSHvwK6Q9+JXNoLFaF7KHOWMd+L0x2MwZVwOzBLBGCLFDCLFQKbsXwLNCiBMA/g/AQ0p5CoATmtfmK2U6QoiFShN4ZlFRUYsunoiIvK/4rO3IZevc167UPM8faH8alHYgmPlzz+lr6qv3nzJv9+ka6fTzjFg3lXtjsJkzrgbmKVLKMTA1Uy8SQkwDcAeA+6SUqQDuA/CaOx8spVwqpcyQUmYkJnpvjhoREbXOy+uPAQAuGdEdgQHGA64amyWe/uYg/rHmsN33uWxkD7vH1BqzVrNVjblO8zCQU1yFF3+w7WfOK6lG2oNf4dHP9hl+zu4T5br9dtOULaUsUP4+A2AFgPEAbgDwqXLKR0oZABQASNW8vKdSRkREHcCxItOo6g1HivDE/KGG51TXN+GV9dl44bssw/SbgGWAlxGjEdZ1VrXy7nH6fuxnVh22ae6e9uz3AIBlm48bfs7jXxzQ7beLwCyEiBRCRKvbAGYD2AdTn/J05bQLAKgrWn8O4HpldPZEABVSypMev3IiIvK6I6ct/cvNEvj1+F664+PT4m1e882+k1i2ORd78yt05Ua1YkesV6QyCuz//PYoqutNg9DOnK3VHcvWTPGyx9XBaW3JlXnM3QCsUJ5eggC8J6VcJYQ4B+B5IUQQgFqYRmADwNcALgaQBaAawE0ev2oiIvKJ+z/cZd6+cVKaTc12SI8YbMsttXrNbsP3iosIduuzT1XWYmBytHnfqHb7wrqjqK5rxCOXDLGpDd//4W6sXDTZ6efkl1WjZ5cIt67Nk5wGZillNoCRBuWbAIw1KJcAFnnk6oiIyG9U1DRgX0Glef+O8/vanONOLXhUapxbn2/dx7y/sMLwvFc35eCGSWn4ao++sXaXVX9y8Tnj9Jv5ZTU4WVGL0alxdhOltCWm5CQiIpeMfHyNbj8y1LZuF+JGYHY36NVpmrIraxvw9DeH7J479Znvnb7frH+sNyxfsHQLfvnyZizffsLweFtjYCYiIrc9qRn09e9fjTZvu1pj/tsVw52ec+0Eff+1dkrWqYpa69PdVq7J9/3NPVNtjr+ywTT6vLK2weZYW2JgJiIit/16Qm/z9iUjLDmkggIDsORKS9BdMC4V1mLCgvBrq6BrZPEVw5G7ZB6uyTC9x+78cgBAU7PE7H9uaOmlAwCes5rGNbh7jM05J0prcPOb2zHisTW45c3trfo8dzAwExGR27Tzl7UDwCpqGvCLsZZ1mKf2t81T8fYtE9z6rIoaU401Smk6L6+2TQnqyMJp6ebt5mZTP/W/v7PMeb57Zn8AwPXn9Ya17w6dAQCsU/72BgZmIiI/8vXek5j/3x9xwsGyib5QZJDty0jPLuEIDgxA7pJ5yF0yz2auMQCMdHPQ15jepvNrG5qwv7ACYzWLZUSFBiF3yTzD180clITcJfPw54sHIzLENLWqpKoen+zI1513nVJ7/+ulxnOyVRuOeCdLJQMzEZEfufPdndh9ohz/+a5lqyW1Fe0iERPTbecqqwKspk/1SYhEtGaQ2G8m2tZKnQkLNgXVmoYmvLs1T3dMXTjjkXmDbV6nreVWKUlOcoqr8PuP9NO3gpVBaPaymKkOnap0eNxTuB4zEZEf2n681PlJXlJYXoP3NAHxX9eMtnuudWjrEhmCrQ/PRLM0NUn3sLPqlCNhSiKRXSfKMShZ3xc8Z2gyAODWqem4bFQP3PzmdvOUrpgwS4ib2r8rNh4txtEztgtwaEeSZy2ei8LyWnPGMK1Jfbu6fe0twRozEZEfCnOQrtLbnlmln5YUFmwbOoYog6eG94y1ORYREoSo0CCkxIUbptp0plHpFw4PDjQ3Sau0tdyk6DAc1qx8tWhGP/N2lVKzzi+rsXl/bWAOCgxArwTj5CL2yj2NgZmIyA8Fu5mu0hVNzRJf7inEqYpaLN+Wh/s+2OVSbuhv9p3S7Rulwnzz5nFYceckDO1hG5hbS1056sjpc8gp0fe9B1g1Pzc0WZKQ3DrVMuhLre2uP2zbTxxk0IQ9uV+CeTsuIhir7p2KmDD3MpW1FJuyiYj8kPViDJ7w5Z5C3LN8FxKjQ82DuUICA/D3q0Y4fJ31GsVGSUSSosOQFO1+M7UrTlea5ixX1DTYDMAa1sN2mpNKW5tOiAoBABw4adtPbFSLT44JN2/fNaOfTRN6W2KNmYjIDwW3QSrIH7OKAehHWH+QecLpQ8B56Zba48vXjXU6SMrTyhxMj7ppch/DcutYGx5sW8v/1zWj8MZN45x+/kVKP7a3MDATEfmhHcfLPP6e9lJgvrctz7BcVac0d390+3mYM8y7QQoAkmOMa+K/Gt/LbgpQ60eHUIN+8ctHp2DGwCTD13fXDFKLdXOxjdZiYCYi8kMpceHOT/LQe2adcbwcYq2So9qo1ukNsw1qrJP6JuCPFw20+xrr5ulubjazh2sGmXmrb1nFwExE5CfqGi2LNBSU16Cg3HYEcWscsxOAl29zvFjD/kJTv6zRaGxvMGo6f++3ExEfGWJTri4lmRgVqitPjdePqI4IcfyQ0Sved8s+MjATEfmJsir9Ygm78so9+v5drAJZV2VAlFGA04oNNwW7mHDv1hxb4sn5wxAbHozHLhuiK+9uNX/61RsyHL7PrMHd0DcxErdMMe7DbksMzEREfmLdodO6/Td/ymnV+7279TjSHvwKn+40paDU1sgBy3rIBeU1eHvLcVz10k82axj/nFdmzlUdF+44gPuDS0f2wO6/zsacYd115db9686ShYSHBGLd78/HXy4Z4vC8tsDATETkJx5esU+3b5QMoyXvd/+HphSU1nOWtTXlv6zch8zjZTarLl3x4k/mbXfWWvY07WpUNxgsNuGOcWldWns5bYrzmImI/JQ2WYaRH7OKcbyk2qUlFEvO1ZlHV6uMpmRlF1cBAGb/cz2OnHY8KMybFl8+DJeO6IGuUSFIT4xq0XscfGIOtmSX4PyBtite+RMGZiIiP9XQ5Dgr17WvbgUAjO3dBQOTo3XHKqr1/dUvrDuKugb9+208Wmz4vpuOFvtVUAZMo6zP65vg/EQHwkMCMWOQ8fQof8KmbCIiP6EOUHpgjmkakNq3a0Q7z1nNjKVVdE6/TOM7W/Owar8+tabR6wDgq72FNmXPOMkORp7DwExE5AeamyVOVpgCpXbO7f7CCptzaxua8IuXLH2/S745ZHPOwyv26vabmm2bxSfZqYG+bzB9apBVjZzaDgMzEZEfaGi2NDN3jbbMwd1XYBuYi61qw0b5n7fm2F82MqO3afDT4iuGO53PqxrWBotTkDH2MRMR+YGfNXOWExzMKy6rqseUv9uuFeyOxy4bimEppkAbERKI6nrTNKr+SVE4aicJifUqTtR2WGMmIvID2oCoHcgVHqKvP23Ksh2wZdTMfNnIHgCAwd1tV0UaoikrPmdZIGKo1UpNahatRTP6Orx28izWmImI/MC+fFOT9W3T0hEcGIBfjOmJT3bmY1tOCR76ZA+qlFrtTZPTbF576NRZzP/PJrx43VhzPuxTSn/1laNTsFjT1J27ZJ7da+iXpJ+GtOvR2a36mahlWGMmIvID0WGmelK5Ms2ppqERANDYJM1BGQBOlBonHdmdX4Envzhg3k9U+qmr65sMz1c9ftlQ83ZokG8WqSA9BmYiIj9Qr8xZHtzd1Cw9a3A3AMDy7foR0t8e1Kft1DpbZ5lepSYTGdTd8WjqGyalYdejF+LIU3PRrFmXecmVw924evIkBmYiIj+gpssMUWqt0nHSLwC2/cfaKVFqoHcljWZcRAhCggJ0CTwGGfRNk3ewj5mIyA9YArMpkKZ0cbwe84vXjsGInrG6EdpbsktRXl2PuIgQ1CsLVoQapN20Z0TPOPzv+gwEBlgWuCDvY42ZiMgP1FnVcK3nKlu7aGgykqLDbMr/810WAEugDw4KwEwlDeWUfo5XVAKAC4d0wwWDurl+4eRxrDETEfkBc41ZqeGWVtU7Oh2BAQKBBnOLtyupOs1N2YEBeO7qkfgoMx9Xjknx5CVTG2FgJiLyA2pgDlVqzF0iWrb28e4T5cgvq9Y1jcdFhOC309I9c6HU5tiUTUTkRdLOqC7rPuaZg/WrIH18+3kuf0ZeabXN+1H7wX8xIiIv2V9YgbFPfYs3fsyxObY5uwSAJZCGB+vnFGekxSMs2PZX9uWjetiU1TU0I7ek2vR+bgz+Iv/AfzEiIi9ZuiEbpVX1eFyTCESlNmH3jo8AYFp/OC3BtP3YpUMAAG/eNB7JMWF4+5bx5tfdOjXdJvhW1Teat7vF2A4QI//mUh+zECIXwFkATQAapZQZSvnvACxSyr+SUj6glD8E4Bal/G4p5WrPXzoRUfvy+W7LOsc19U0IV1Z2amhqNicE6RplWVnqhz/O0L1+YnoCtvx5pq5sWEosjiyeCwC4d/nPWLmrEHe99zMAIDo0iE3Z7ZA7g79mSCnN2dOFEDMAzAcwUkpZJ4RIUsqHAFgAYCiAHgC+FUIMkFI6zgtHRNTBabuXn/7mIJ6YPwwA8PXek+by1qziVFbdoNs/W9do50zyZ60ZlX0HgCVSyjoAkFKeUcrnA1iulOcIIbIAjAewuVVXSkTUgby/LQ+jUuPw1ubj2H2i3CPvecGgJKw/UmTe75sY6ZH3Je9ytY1DAlgjhNghhFiolA0AMFUIsVUIsV4IMU4pTwGgTe6ar5TpCCEWCiEyhRCZRUVF1oeJiDqU8mr9vOSGJon7P9ztsaAMAMKqsl1QbrzgBfk3V2vMU6SUBUpz9VohxCHltfEAJgIYB+BDIYTLE+WklEsBLAWAjIwMF7LCEhG1X1tzSr3+mTOZwatdcqnGLKUsUP4+A2AFTE3T+QA+lSbbADQD6AqgAECq5uU9lTIiok6ptqEJt729o80/xzopybCU2Db/TPI8p4FZCBEphIhWtwHMBrAPwEoAM5TyAQBCABQD+BzAAiFEqBCiD4D+ALa1ydUTEbUDZdWO02uq/nLJkFZ9zuyh3TCypykY3zy5D26ektaq9yPfcKUpuxuAFcLUeREE4D0p5SohRAiA14UQ+wDUA7hBmlLa7BdCfAjgAIBGAIs4IpuIOrMTpZa+3rSECHPyD2vTBzhfZMKR0KBAfHbXlFa9B/me08AspcwGMNKgvB7AdXZesxjA4lZfHRFRB1BZY5nG9NTlw3Hda1sNz+ubGOWtSyI/xpnnRERtbGeeacWnS0f2wJT+XXHZSNs0mv+7PgPCelg1dUoMzEREbaxZmXdyrtZUc144LR0hQQH4w+wBmDkoCaFBAZiQHu/DKyR/wmUfiYja2I7jpqlSMwaZVowalhKLPX+djbDgQEgp0dAkmTqTzBiYiYjamBp0m5stKRvClNWjhBAICWITNlkwMBNRu7avoAL5ZdWorG3EzEFJSNAsAuEvGppMAXlAcrSPr4TaAwZmImrXLvn3JvP2laNT8I9rRvnuYuyobTDNGI0I4a9cco6dGkTUYXz6s38mGcw6cw4AEBbMX7nkHL8lRNShlJyr8/Ul6OSVVKO6XqkxB7PGTM4xMBNRuyWl7fo3h0+d1Q2y8rXs4nPm7dT4cB9eCbUXfHwjonarvqnZpuzXr5qyal03sReeuny4ty/Jxg+HTcvaXjw8mQlEyCWsMRNRu1XXaBuYVe9syfPildjXrNTq6xrsXyuRFgMzEbVLO/PKMOKxNQ7PUUdD+5IakC8cwrWRyTUMzETULl354k9Oz1FzVPtSXaPp4YCZvchV/KYQUbujTj9yxlfNx2cqa/HCuqM4c7YWK3cVAjAtyUjkCg7+IqJ259HP9tmUzRmajFX7T+nK1Nqqtz39zSGs+LkAP2YVm8tCWWMmFzEwE1G7cq6uET8dKzHvD+0Rg9dvHIfAAIG5w5MBAI9/cQClVfV47PMDmDOsu9evceUuU6KTrTml5jI2ZZOr+E0honZl8VcHdPvv3joB3WLC0DUqFPNHpWD+qBSkJUQAAE5V1mKln2QDiwhhUza5hoGZiNqV7w6d0e2rqzRp7S+sNG+/8VNuW1+SjYHdbBerGN2ri9evg9onBmYialesc4oY9d0mRIaYt3efKMeW7BKbc9rSsSL94LRrMlIRGMDkIuQaBmYialf6Jkaat2+clGaYTev5X43W7Wtr0N7QNzFKt8/+ZXIHB391In/6eA8+yDwBANjwxxnopfTDkffUNzbjute2YltOKab274qdx8tQVd+Euy/oh/tnD/T15XlVQ1Mzrn9tGyakx+PeWQNcfp06oOqTOyZhbG/j5uFxafG4cEg3rD1wGgDw5JcHEBUaiGvG9Wr9hSs+21WAe5bvwpR+XfHOrRN0xw6dOqvbP1vb4LHPpY6Pj3GdRHOzNAdlAPhiT6EPr6bz2ni0CNuUwLLxaDGqlFWHXvguy5eX5RObjhZjc3YJ/vXtUbdeF6kMouoeG+bwvGsyUnX7f/pkr3sX6MQ9y3cBADZppkQBxlO01LnMRK5gjbmDySupxkvrj+HO8/uioqZBt4i81rOrD2PRjH5evjrKL6vx9SXo1DY04ckvD2BCegIuG9nD7dc/9vl+vPlTLtb9frpN860972/LQ2VNAxKiQt3+vOZmaX6YSY5xHJhnDemGN24ah5ve2O7257gr7cGvAACv3ZDht2tCU/vBwNzB3PbODhw8WYmf88pQVl3v68shK47+TSpqGhAbHuzFqwHWHjiNd7fm4d2teS0KzG8qI55vfSsT3//hfJde89CnpprrXS14MCyusqy1HODCYCqj0dFt6Za3Mg3Lr5vouSZ06vjYlN3BHDxpGuRy6NRZnK50vGB8o8GSedS2vj142u6xp748YPeYK/LLqnHNK5txvKTK5dd4Kpd0TrFrn1lvZzWoj3fko8nBGsprD5zGTW9sw8GTpr5bVx9gesTp1z+e9PQ68/8Rd52ra8TbW44jV/lZtSO/jaTGh+O1GzLwh9kD8Phlw1r0mdQ5MTB3MnfP7G/eXrox24dX0jntK7AfFD7akY/DVoOG3DHl799ja04ppj/7g8uveePH3BZ/nparyTNKqywtBjWalZ/+8NFurD1wyuglAIDfLsvE94eLcMPr2wAAcREta1korKjF3Oc3tui1S9cfw19W7sOi93YCMF4LWquxSWLm4G6464L+nCpFbmFg7mTu1QRmdQF38h/awNUaB09WQkr7NVAjRufvza/AidJqp+ePUZJnnK6sxZnKWptzt2SXILvoHKrqG81lb1ol/jhw0vahpKa+yfDz3fnRPrljkusnO7At1zRob39hJWobmnC2ttHh+c1u3n8iFQNzB2JvSoa22U/bLxccyKd4fyPRsl/mT399ULc/9/mN+Cgz3633OFenDzRHT5/Fpf/ZhKnPfG+4rnGdpll6U1YxahuaMP3Z7zHtWf35m44WY8HSLbjgufU4pwlm1k3X+QYBWP18a3l2HhaM2JtS5a4DmrnQg/6yyun5nLtMLcVvTgdy2qCmAuibDAGga5Spb2zWYC7c7i0NTc26YDUoORqjUuPwyLzBePm6seZye32wzryywbZb4oPME5BS6mq2zc0SzXb6cq0Dszb1ZVl1vfl1Tcp7aANzaFAASqvqUdvQjNqGZpRXWx4Sf9b0YxedtT/uITXedl69q8s7epL1PVONspNSc2r/roblJ8uN/z8SOcNR2R3I/zbkGJYnRoWioLwG8cpglctGpuD1H3OwN7/Cm5fXaX25pxB3vfezeT8kKACr7p2mO2fW4CR8e/BMiwOzkR3Hy9Dnoa8BAKvvnYZe8RGY8/wGhAUF4su7pyA4MAAxYUGoVGqxNfX6B7hczSCyTUeL8fQ3h5AaH4HdJ8ptPquusRlvaZqmD58+i2RlnnHmcUtgvnWZ7ajlB+YMxDOrDqPWav6v9fV4SmK0/WlaUkpc88oWhAYH4O1b9ElDjFoNAODtWyagtKoeY55cqyuPdzI4jMge1pg7kFNKjTkqNAhqi/Vt09LxyLzB6BoViifmDwVgqUFHh/G5zBus8zQbBV+12dPZgCJ7LhnheGnDxV8fxPHSKhwvqcbh02fNrSvaWm9tg/6zB3ePMW+/tikHpVX1hkFZpa21a5vWe8TZn2/8x4sGIlxZhKLWwYOBtZmDkuwes/c5qjo7ARYAquubsC23FBuPFtvUmh29rovBYLSH5w126xqJVPzN3IGsP2IazPXGTeMwLi1ed2zucMsv7jG94vD+tjy8tfk43tp8HIBvUnSerW3APct3YXK/rrhlSh+vfrY3WQfiB+cOsjnneImpz3RXXjkuGWF/PvEH2/OwKasE/7x6JIICLc/V+woct35IKbE127I28KmKWmTmlukCc/E5fTNzY5MlMFmnmHRGfcCorm/E+9tO2D1vXFq8ecGHtzYfx8asYnxzz1SEBgXaraECwC1T3fu+LJrRDzdP7oPBj65CZW0jvtxTaHiftYPT6hqbdStX7XbQwmSdr/vgE3MQzmUeqYVYY+4gtGkAezsJsNUGTYR/+mSPx6/JmTX7T+O7Q2fwZCvn7/o768BslPRCXWTh1U3G3RGqP32yF1/sLsTq/fr50M4GGjVLib9+vt+8X1JVj3s/2KU7x3osQktr7wBwwUBTjXbdwTN2z4kMCcTAbtG61aGyi6qwRwmA1tejFRrkftDT3iNt14KWtl+4zurfzWCtDNw4Kc28Palvgnk7LJi/WqnlWGNux6rqGnH/h7sQFhyom5+cFO04VWFXg1SI1rUlbzhZYUlP+ezqQ6iub0LXqFDcMqWP4Rq77ZUa4G6anIZrJ/RCvyTn2ahyi6vwxe5C3DylDyJDbf+bbjhShHma5mvrZmhrP2bpm9ONkpDc9vYOXU3PUbO1VteoEBSf00/zenVTDgIDBIrsfK8WTkvHHdP7IjYiGEN7xOqOZZ05h4zeXfCY5kHCWkumIlnPJc4vq0bPLvqHWO2DkfaBqqGp2XCK1qOXDDFvv3nTeOSXVaNbTJjhildEruJjXTv21uZcrN5/Gp/tKsTcf7meNKFfkm1OY22zqLc0akYH//f7Y3jjx1w8u/owlhqMMG7P1F/w56UnuBSUAdM0oefWHsFza46Yy7TN1doFSQDLwCSjf1sjf/v6kGG5mjwDAL7ZZz/hh9Y141INy1/ZkI1Pdxrnje7ZJRxdlMFR1rX9A4WV2F9YiSOn9SOy52m6Y5wtYOGKC55bb1P2xW7LYhPaVii19h4VGoQgTYDXTj8MCQpAemKU4YMUkTsYmNupqrpGPK9ZlcedZseBybbBwdHAlrZyxs7UmY1HO1biE7Vv0lFz85OXm1I2ju4Vh6ZmaU5esT3X0i+s9sUaUe/lGzeOwyu/GYsebgSu53450ry9PafUwZnG7ps1wO05u1drVn6yrltuzSkx97lrLb5iGP40ZxCW3TzepqbrqtWa0fDWXQzWg71+0rQyqA8+YcEB2PzQTFw2sgfW3KcfWU/kKS79bxJC5Aoh9gohdgkhMq2O/V4IIYUQXZV9IYR4QQiRJYTYI4QY0xYX3tk9u/qwTR9Ya2S7mOvYkzYdLTYs357rmfzN/qCqrtE8d9dR8FJH9R45dVaXmnKvppb8llWmLJV2kFRSTCguGpqMO91YIEI7+rqbEtDdyRoWFBiA89ITnJ+ooe2qCLXqjz1y+hye/kafMCUyJBBxESG44/y+mDYg0a3P0rJ+KNXO6bb+//SAZtxFWZVlXnZidChe+NVoDPDyAhnUebjzmDtDSjlKSpmhFgghUgHMBpCnOW8ugP7Kn4UAXvLEhZKedTpDT6iobrvF3KvqGnHmrGVgTUV1A8o7wepXFTWWexrqIDAPSjYFx7iIEBwr0j8kqYFX22+pDfLa1JDqoKgF41Jx9wX9cNXYnk6vcUgPS2BWV5iyN/DqNxN74/kFo2zKtYOdnC3HqE2oAgDdY8PxwJyBujLr5TFfv3Gcw/d0x+3T+5q3tS1NdQb99Gp2MrVZu7wN/48QqVrblP1PAA8AujyC8wEskyZbAMQJIRxPsiS3uJsD2Uh610ibspFPrHGYmaml6hqbMPSvqzF+8TrkFlcht7gKI59YY05sYaQtHxK8SdsU7WgksbooQ1V9I55dfVh3bNBfVqGytgE7NIk66hubzQk41MCdollJKSgwAPfPHoj/++VIjOypH1xl5M7zTcHqH2uPoKlZYsijqw3Pe/LyYZg/KsW83y3GNJBQWwNeev1Ym9edP9BSy50zLNng8/thyZXDDT9z0Yy+mOBmjdyRB+cOMqep1bY2zHl+g82564+YRpWrzd4jU+M8dh1E9rgamCWANUKIHUKIhQAghJgPoEBKudvq3BQA2pEp+UqZjhBioRAiUwiRWVTUsfoU25qj5fG0iRQc+ftVIzAurQtWLpqsK3dnyUBXaRdm+HRnvnkxAEd25Zd7/Dp84eMdlnzVjpo+1cBmr0aWa9DVoCYJ0fZ/GhmlCSbWg8PUYFiuqdlX1ztenEErQKnFd9PUktMMHvqeunwYJvSJx9u3jLf7XkYBOy0hAr+ZmOby9bhKvVfaloGTFbYpNNW+brWZ21GrB5GnuPotmyKlHANTM/UiIcQ0AH8G8GhLP1hKuVRKmSGlzEhMbHmfUWdU7WCg1rAU57UjwJTY4aPbJ2FUapwuyb+n+q2bmyWuf30bnvjigK6J8IXvsvA/O6Ou/zB7gHm7sqb915h3nSjHRk0/uqM+5jAnv/Dvft8073aEpvZ7VMkjrU5VsjfF7Bea5uxv75+uO6YmntE2wgx/bI15OzVev56xNTUwzxhoycQVFWI7Krlnlwh8cNt5mNrf/v/1uAjbFJbrfn++ObWnJ6n36nhJNXadKEfag1/pjl87oRcA4PEvDuDDzBPmpmwuTEHe4NK3TEpZoPx9BsAKANMB9AGwWwiRC6AngJ1CiGQABQC08yd6KmXkIRuPGA+aAoCQFkx70uYk9lSu5uziKmw4UoTXf8yxGTF+1GphggXjUhEcKHC1ZtpNeAeYx6xdBMIZZ9PVcpWam7bG9vCKvQBgzjFt3S+rGpgcjYTIEEzpZ1psYdZgSxCNVqb2ZNsZ8X3PzAHmpvDfXWAZUPbMVSMAAH//helv7RThgAChS7zhDrWPW9VW6xgXKPfqzNk6XP7fH22OazOdPfDxHvPDJWvM5A1Ov2VCiEghRLS6DdNgr+1SyiQpZZqUMg2m5uoxUspTAD4HcL0yOnsigAop5cm2+xE6j892FeBPH+/BkdOWXxpHnpqrO6clT/TnaTIW1TV6ZtqUdq6nsxWClvxiBA4+MQdJ0WGYroy49ef8DO9tzcNTXx7ArH+sR9qDX+kScRw9fRYPfrIHheU1eGGdZTrbIx7Km6ztp1anSC1db2qBmNzPuB82NCgQW/8809yM/MAcS0pQdR6uvc6RoT1i8NldU3D4qTn4/WxLN8nVGanIWjwXU+ysrPTYZUMtn+HGv+ULvxpt3r53Vn8HZ7aO2k+utkRYs05gok5Va0nGMSJ3ufJbvBuATUKI3QC2AfhKSuloMdKvAWQDyALwPwB3tvoqCQBwz/Jd+CDzBH44YumTtw7ELXmi1y7/6Kmm7CbNL7Y7391p97xxaaZmdLXGqPb9eXKVJU/784q9eHVTjvmBY76mxnX1K5uxfPsJzHtBn/Blcj/jAOaI0UOWtkxdbnCzskjG13vtJwQJCgwwj+o2ahq2N3pbHSRlFJC0tfz0xCgEBghdH7Zaa3Y3D7pas2/NtChnwkPs/z+5cVIarpvQW1d2VlkSs6wTzCQg33OaokZKmQ1gpJNz0jTbEsCiVl8Z2ZWl1JjHKwtVbHloJiY+vQ4AEBToflXzvL4JGNEzFnvyK1DogTVkv957Es+tOez8RADv3KpfWk/N413up33M9pYiPFFajR+zilGmDN4q0wziWnPfNJfmvF41tqd5sNj6P56PACEw9ZnvdeeEBgXg1xN64b2teejZRd//q66z7UxMWDA2/HEGIkItwfaXY3vigY9t86V3czL1SZUYHYq1900zB3IAeOjiQZg7LBmjesW59B6qF68bgxOl1TapOj0pzEHN95F5gxEgBH7/kWVc6ytKq8Skvu4/YBG5ix0m7VCVEhz6dTPVTrQ1oKgWpgNMSzCNpD10qrJV15ZfVo07391pMxfXyLUTetnUxNTBUg99urdV19FWnvzKeMGNC577AQ/auWZXE1Fol+FMiApFanwEJqbrVwlraJL4TlkYwnrVJm3LhzO9EiJ0OdPt5XZ2p483PTEKCZr3DA0KxIT0BLebf2PCbPNne1qjg5kNQYEBulSbWlGhbMqmtsfA3I5pn/qXL5yI5xeManGqwr6JUeb3PHr6LGobmgyn6DhjlEpR629XDEfXqBCM7BmLR+YNcXju2gOnUeKDxTUcWbPfuLm4oan1c8vjNaOS1UF8z16lb6walhJjXnfb2j0e7JP905xBeP+3Ez32fv4mPtK2dWFYSgw+uv088/7GB2bYnDMqtYtNGZGnMdt6O2E0dzk4yPJUP7GVCRjUaTEfZJ7QLZDwyR3nYWzveHsvs3Htq1vtHpvSryt+PaEXfq1MRXHmt8tM2V9zl8xz+fPbmqOaVmulxlseqoKVLonU+AhcNrIHPlcWV2hskoiLCDbPdz5XZ5lznBBpu2qYO5JjwnCqshbDUmJwx/l9nb+gHesRZzsN7MvfTdXtp8ZHICUuHAXlltHuidGtu8dErmCNuZ0wGgyVYvDLpaXsjeZec8Cy7q+jxCau+M+vRzs/yUBrP9eTmtysGf/7V67/zDMGJuHSkT3w8MWDdU3Lah5twJRJ7M8XW0Z4a1s1WjvH9pmrRmDmoCQ8rhlR3VFp5+4HBwq7rQNv3axPBdoWc6qJrDEwtxPatYtVriYTcUWgnT7G15X1aZ/68gDGPrUW+WWOm6rtGdEz1jCBhNF51vx5hLYz0we6PrI4NiIY//7VaPx2WrquXNsaEhUahERNP+6BQtOYAFeXe3Rk2oBEvHbjOLdaSNqrJE3N99CTc3VTBrX6JUUjQlmfOqiN5lQTWWNgbie0zWmq/h74ZawqtTMNRO07fXVTDsqrG+yur+vMnvwK5ycBurmyKn8KzDM1yTkcuXRkD8wanGRO4NEa2trd4iuGY5JmvrL6vSj2s754fxcREohZg5Nw2cgeTge43TG9L6JDg7DIjRW7iFqDfcztwPeHzuCmN7fblLd0BLYRR8FPm67wH2uP4K4Z/QxHraq5mwMEcOxvF2P59hPm0dVDNEsLOmK0uMb+wgpMcmEe8MajRVi6IRtPXzm8xYPg7KltaMJd7/2Mbw+amvbvndUf984agHUHT+OWt3QroWLBuFQsUTJieUJSTJhNP/ug5GgcOnUWzytJTOYa5Jkm+4QQePUG11as+t3M/vjdzLZLdkJkjTXmduB5TQYpLXtTXFrCnak2uXYWutiiJLpolqZr0yaI+O0015JMGA2uuXVZpsGZtu54Zyc2Hi3Gv9dluXS+O/YXVpiDMgDzSk9GfY7a1KJtxfpBqrLG9YUniMi/scbsxz7MPIG48GDdaj9XjE7BLVP6mKc3eYp2RLAz9rKDqascqfmOU+LCsevRC3GurtHlGqzRQgzVdpJ6aG3LKTWPUM4rrcYnO/JR39SMypoGZKTF65qDTT9DE1764RimD0jE6F62U2Cq6xvx9NeHsDWnBKcr6+wOtFPnfwPAy9eNwfg+CYZTcTztjxcNxB2ajGotySxGRP6JgdlPlZyrM2di0g6Iqm1o8uigr5awF5jVrFhxmlHEcREhLg36cqamvgnhIfaTO1z9ymbLuQ1NuqxNgO2Uq092FOBf3x7Fe1vzsO3hWTbvt2zzcby95bh5v8IqE9kFg0x9zdoHidjwEK8EZdNnBev2HaWYJKL2hf+b/ZS2lqgdOGWdhtEX6qyWnaxtaMLq/afwj7VHANhfftBVGx+YgfsvHKDL+20dGB3ZpVlUwkhtQxNe3WRKsaguBGFt87ESh+9x/XlpAPSZsVxNiekJA5L12cQcpZgkovaFgdlP2es+dqfJua1Y15gfXrEPt729A5W1pqZkZ2sLO5MaH4G7Z/bHbyZaFhL4Oa/M7vn7CpyP+C7SBOCnvjqAbCcpQ3OcZD0zGskbY1WLbUva9J2A82Ujiaj9YFN2O/LbqX1w+eiUNnv/1fdOwxe7C5HSJRy5JVX4Zu8p5JXazls+VnQOU/t3NQ8++2Rnvu64p67x3gsH4FVlHnWVg37mn53UkAGgpKrOPLDs3a15Ts93Zz3ol68bg+Jz9S4v+OAJ1vmn+ybajmYnovaJgdlPqasMaT3sJLd0aw1MjsbAZMs84ofmmjJMaadLAcDjXxzA+9vysOa+6eYBX1p9DKY8tURUaBB+M7E33t5yHKv3n7K7NOFfVu5z+l7aUczSSfKuxqZmHNaseW1tqtUaxHOGdXf6+W2hb2KkebGQ6DDv1daJqG2x/ctP/etb/RSp5xeM8s2F2HHktGkt4tIq28QknpzGpa7PHOlg4Jcr3ElScrbW8dQjZ4tveIt2pLt6n4io/WON2U/sOF6G+z/chYze8dhp0J86f1TbNWE7ExsebHfwVVtn5RrfJwH/25iD7w8XGR53NY+2ep3His45PfcPViO6tfxpQY2bJqdh/RHTfWntgDsi8h98zPYTf191CMdLqvHJznybgUfaRB2+8MT8oQgNCsBgg+xdDU36wDymV5xHPztSWf822U7/7Z78ct3+nKGWDFiXKvOpAaBOuc4l3xxy+pnrDp0xLH9k3mDDcl8Z0j0GseHBGNu7C4I5+Iuow2CN2ccamppxxzs7sC2n1OZYgAAOPDHH57Wh+aNSMG94d7z4wzEcPFlpLr//w124QhnolZ4YiVX3TDMvV+gp6oCqes0DwH++O4rwkCDMHZaMK178yVy++6+zERsejJr6JgQEmAZInattwPeHi/DMqsOYMTAJazWrZamufnkzPlTW4f3hsD4oZ//tYjRLCSGE05zK3pYUE4btD8/i4gpEHQwDs499vqsQ3x40rqE1S/9pogwKDMDKn/ULWHy6swBnKk3TkKLDglu97KARdS6zOne6srYB/7fGNF8626pZWl0wQpuIJCLEVHbwZCXySoxXxtqWW4rTlbXoFhOGG9/Q5yQPCBAIgP8Gvra450TkWwzMPuZO4gxfqzQYFKU2u4d4uKasUqcFFVbU4oGPd+PHLEviD+20pxsnpRkurKENXD8c0T8APX7ZUPz18/0ATH381iPM9zw2u/U/ABGRm/i47WNRYe3n2eiiobYLXag12raquYVqRht/mJlvuPwlAAzuHm1Yrm0Cf/Sz/bpjszU/z53v7sT9H1oGfUWEBCKGU5CIyAcYmH3M032ybenRS4dgVGqcruxEmal5uK0GH4W6GPDHGCxEAQAJdnJXv/Kbsegeaz+96Ye3nefS5xIReRoDs4/d94H9qTn+JjQoECsXTdaVNTSZpiuFtFFgdvV97S2UYS/ZifVqU9Z8vVAIEXVeDMx+7I2bXFvI3R+0VVO2K8lK7r6gn+E6zgAwb4RxVq6uUcbnA8Bjl/pHAhEi6pzaTwdnB9TsIDlG1uK57WphgpJzthnAvOEPswfgrgv62z2eFG07//m89AS7598+vS9unNzHI9dGRNQS7ec3fwdkPSJ71uAk87Y/B+WF09Jtapybsx0vk9gat01Lt3usJQtHZB63zBlfcuVw3bH21OdPRB2T//727+Byi6sw+sm1lv0l83Be364OXuE//nzxYGQ+Mstrn/fQxYORu2Se+c/dMy015DFO+ooBIN2qn1k7UG3B+F6Y3M9Sg2bfMhH5GgOzjzy8cq9NmZrOcmTP9hEcbp/e17z93C9Heu1z6xot841dScDy+9kDdft3nt9Xt3/LFEvTdYYLgZ6IqC0J6WwNPC/IyMiQmZmZvr4Mr8p46lsUn6sz76uLI5w5W4vY8GCb9Xb9kZQS+WU1aJYSvRO8tx7wd4dO4+Y3Td+XnX+5EPF2pkRpFZ+rQ3VdE2oamjCgW5TNoLJTFbUIDw5EbATnLhORdwghdkgpM6zLOfjLRwYmR6E4yxSYfzOxt7ncaLCSvxJCIDU+wvmJHtYjzjL/WF3kwpmuUaFAlP3jybHt574TUcfGwOwjgQGmXoS0hAg8MGegk7NJa1ByDJbdPB7xkSHtomWBiMgdDMw+IKXEBmUd3dum90U0Uz+6zddLYRIRtRUO/vIBbf7mi4cZJ8AgIqLOiYHZB17dmGPe5mAjIiLSYmD2gWdXH/b1JRARkZ9yKTALIXKFEHuFELuEEJlK2bNCiENCiD1CiBVCiDjN+Q8JIbKEEIeFEBe10bUTERF1OO7UmGdIKUdp5lytBTBMSjkCwBEADwGAEGIIgAUAhgKYA+BFIQSHzmqoKxtFhPC2EBGRXoubsqWUa6SUjcruFgA9le35AJZLKeuklDkAsgCMb91ldiynKmoBAM8vGO3jKyEiIn/jamCWANYIIXYIIRYaHL8ZwDfKdgqAE5pj+UqZjhBioRAiUwiRWVRU5M41t3tqpqoArpdARERWXA3MU6SUYwDMBbBICDFNPSCEeBhAI4B33flgKeVSKWWGlDIjMbFzzUndW1ABAOjZxftZs4iIyL+5FJillAXK32cArIDSNC2EuBHAJQCulZak2wUAUjUv76mUEYC8kmrztis5nomIqHNxGpiFEJFCiGh1G8BsAPuEEHMAPADgMillteYlnwNYIIQIFUL0AdAfwDbPX7pvVdc3YuXPBWhqlvgpqxj/XHsERWfrbM5rbGrG/sIKNDebnlsKymvMxxKjQ23OJyKizs2VlJzdAKxQVuMJAvCelHKVECILQCiAtcqxLVLK26WU+4UQHwI4AFMT9yIpZZOd9263hjy6GgDwxJcHUFpVDwB4ft1R8ypRqie+PIBlm4/jwbmDcPv0vnhtkym5yMT0eO9eMBERtQtOA7OUMhuAzWK7Usp+Dl6zGMDi1l2af6lvbIYQQIAQCNSM2lKDsj3LNh8HALyy/hhun94X6mqDg5Jj2uxaiYio/eIiFi744fAZ3PjGdvP+GzeOc/s9yqob0NQsUdtgajw4f2DnGvBGRESuYUpOF7z4/THd/he7C1v0PqVV9ahrMC1gERbM5CJERGSLNWYHvthdiNc25WBPfrmu/NOfXRtkrg74UpVU1WFbbikAIJyBmYiIDDAwO/C79392+zVVdY2IDDXd1s3ZJbpjZ2sbERYcgNqGZo7IJiIiQ2zK9rCaBssA9O8PndEd++XLm1GrNGUnMTATEZEBBmYPq9UE5lc35dg9LyiQt56IiGwxOnhIckwYAKCipgEA0NDU7MvLISKidoqB2Y6jp886PWd4SiwAYEyvOJyqNK0YtXDZDgDAV3tOtt3FERFRh8XAbMeR0+d0+9FhQVj/x/MhBDBveHf87YrheOOmcbhnZn+8dN1Y83kF5TUoKK/BwZOV5rJJfRO8dt1ERNS+cVS2HUs3Zuv2/3zxYPROiETO0/qUm/ddOMDmtZOXfGfevv/CAcguOmdzDhERkRHWmO2ICbM8s0xMj8eMgUktep/AAIG80mrnJxIREaGTB+Zlm3Mx/LHVePqbgzbH1AxdyxdOxPKF5yE5Nszhe/1iTE/D8iHdYzAyNU5XFhwoDM8lIiLq1IH5tU05OFvbiFfWZ9scq200TXtyNXXmbdPTDctjI4Jx3cTeukxfiy8f3oKrJSKizqDT9jHX1DfheIlxE3Nzs8Se/AoArqfOHNAtGleMTsEKq3SdiVGhSI2PwL7HL4IAIAHd6lRERERanTYwb84utnvs0CnLVKluMa5n6Pr24GmbMjX1JoMxERG5otMG5oMnLcE3MsRUK1617xRe25SN313QH4ApbWZcRIjL71lTb8n6tXLRZHSPDeMqUkRE5JZO28f87OrD5u16JUvX7e/swPbcMlz/+jYAQL+kKLfeUzsAbGTPWHSLcTxgjIiIyFqnqDFLKXHgZCX6JkYZ1mAbmiTqGptsyt2t7T58yWCcOVuLy0enQAg2XRMRkfs6RWD+Ys9J3P3+zzh/YCLevGm84Tn3fbDLpqy6vtGtz4kJC8Ybdt6fiIjIFR2+KVtKiWU/5QIAfjhcZC4PDdL/6F/vPWXzWnWJRiIiIm/p8IH5zyv2IvN4mXm/vtEUbLtGOR9tzRzXRETkbR0+MG/JLtXt5xRXAQBqlHWTp/bvave1103s3XYXRkREZKDDB2Y1EKvU6cSlVfUAgJeuG4s/zLZdiAIAesSFt+m1ERERWevQgbm5WdqUCQHdCOyI4ECEBHXo20BERO1Ih45IG7Ms2b2SlTnFDU0SZVUN5vKAAIEJfdiXTERE/qFDB+YblEQhAJCkpNasb2zG6z/m6M7rrqwcpV31KYXN2ERE5AOdYh4zAAQHmp5Bvtp7Et8e0Oe0TooJw6d3TkJ0aBDCQwLx8Y58DvwiIiKf6LCBWUp9/3KIEpiXbrBd4hEAxvTqYt6+d5bxYDAiIqK21mGbsgsras3b//31GAQbDPAa27uLTRkREZEvdcjA/NWek5i85DsAQHRYEOaN6I4NR4pszvvl2J42ZURERL7UIQPzv787at52tBAFl2QkIiJ/0yEDszYPdjgDMxERtSMdMjD3jI8wb+eVVgMA7r6gn815ybFcL5mIiPxLhwzMJefqzNs9u5jmI98/eyCyFs/FwxcPNh/rmxjp9WsjIiJypEMGZu3CFfllNebtoMAA9OsWBQCIjwxBZEiHnS1GRETtVKeLTOcPSMTXd09FUkwoAgKE8xcQERF5kUs1ZiFErhBirxBilxAiUymLF0KsFUIcVf7uopQLIcQLQogsIcQeIcSYtvwBrKmrRqkuGtpNty+EwJAeMS6tx0xERORt7jRlz5BSjpJSZij7DwJYJ6XsD2Cdsg8AcwH0V/4sBPCSpy7WFbtOlJm3H7t0CJ6cP8ybH09ERNQqreljng/gLWX7LQCXa8qXSZMtAOKEEN1b8TlueWdLnnn7xsl9kBTDkddERNR+uBqYJYA1QogdQoiFSlk3KeVJZfsUALXNOAXACc1r85UyHSHEQiFEphAis6jINitXSwWy35iIiNoxVwd/TZFSFgghkgCsFUIc0h6UUkohhLTzWkNSyqUAlgJARkaGW691pGtUiKfeioiIyOtcCsxSygLl7zNCiBUAxgM4LYToLqU8qTRVn1FOLwCQqnl5T6XMK347NR07jpfhhklp3vpIIiIij3HalC2EiBRCRKvbAGYD2AfgcwA3KKfdAOAzZftzANcro7MnAqjQNHm3ufTEKKy5bzquncD1lImIqP1xpcbcDcAKIYR6/ntSylVCiO0APhRC3ALgOICrlfO/BnAxgCwA1QBu8vhVExERdVBOA7OUMhvASIPyEgAzDcolgEUeuToiIqJOpkOm5CQiImqvGJiJiIj8CAMzERGRH2FgJiIi8iMMzERERH6EgZmIiMiPCNPsJh9fhBBFMM2F9pSuAIo9+H6dFe9j6/Eeth7vYevxHrZeW9zD3lLKROtCvwjMniaEyNQsT0ktxPvYeryHrcd72Hq8h63nzXvIpmwiIiI/wsBMRETkRzpqYF7q6wvoIHgfW4/3sPV4D1uP97D1vHYPO2QfMxERUXvVUWvMRERE7RIDMxERkR/pcIFZCDFHCHFYCJElhHjQ19fjT4QQqUKI74UQB4QQ+4UQ9yjl8UKItUKIo8rfXZRyIYR4QbmXe4QQYzTvdYNy/lEhxA2++pl8RQgRKIT4WQjxpbLfRwixVblXHwghQpTyUGU/SzmepnmPh5Tyw0KIi3z0o/iEECJOCPGxEOKQEOKgEOI8fg/dI4S4T/l/vE8I8b4QIozfQ+eEEK8LIc4IIfZpyjz23RNCjBVC7FVe84IQQrh9kVLKDvMHQCCAYwDSAYQA2A1giK+vy1/+AOgOYIyyHQ3gCIAhAJ4B8KBS/iCAvyvbFwP4BoAAMBHAVqU8HkC28ncXZbuLr38+L9/L+wG8B+BLZf9DAAuU7ZcB3KFs3wngZWV7AYAPlO0hyvczFEAf5Xsb6Oufy4v37y0AtyrbIQDi+D106/6lAMgBEK75/t3I76FL924agDEA9mnKPPbdA7BNOVcor53r7jV2tBrzeABZUspsKWU9gOUA5vv4mvyGlPKklHKnsn0WwEGY/oPPh+kXJZS/L1e25wNYJk22AIgTQnQHcBGAtVLKUillGYC1AOZ47yfxLSFETwDzALyq7AsAFwD4WDnF+h6q9/ZjADOV8+cDWC6lrJNS5gDIgun72+EJIWJh+uX4GgBIKeullOXg99BdQQDChRBBACIAnAS/h05JKTcAKLUq9sh3TzkWI6XcIk1RepnmvVzW0QJzCoATmv18pYysKE1ZowFsBdBNSnlSOXQKQDdl29797Oz3+V8AHgDQrOwnACiXUjYq+9r7Yb5XyvEK5fzOfA/7ACgC8IbSHfCqECIS/B66TEpZAOD/AOTBFJArAOwAv4ct5anvXoqybV3ulo4WmMkFQogoAJ8AuFdKWak9pjzlcQ6dHUKISwCckVLu8PW1tGNBMDUlviSlHA2gCqbmQzN+Dx1T+kDnw/SQ0wNAJDpXa0Gb8YfvXkcLzAUAUjX7PZUyUgghgmEKyu9KKT9Vik8rTTBQ/j6jlNu7n535Pk8GcJkQIhemrpILADwPUxNXkHKO9n6Y75VyPBZACTr3PcwHkC+l3KrsfwxToOb30HWzAORIKYuklA0APoXpu8nvYct46rtXoGxbl7ulowXm7QD6KyMTQ2Aa5PC5j6/Jbyh9Sq8BOCil/Ifm0OcA1FGFNwD4TFN+vTIycSKACqW5ZzWA2UKILsqT+2ylrMOTUj4kpewppUyD6fv1nZTyWgDfA7hKOc36Hqr39irlfKmUL1BGy/YB0B+mQSMdnpTyFIATQoiBStFMAAfA76E78gBMFEJEKP+v1XvI72HLeOS7pxyrFEJMVP5drte8l+t8PULO039gGkV3BKbRhQ/7+nr86Q+AKTA10ewBsEv5czFMfU3rABwF8C2AeOV8AeC/yr3cCyBD8143wzRQJAvATb7+2Xx0P8+HZVR2Oky/0LIAfAQgVCkPU/azlOPpmtc/rNzbw2jByM32/AfAKACZyndxJUwjW/k9dO8ePg7gEIB9AN6GaWQ1v4fO79v7MPXLN8DUenOLJ797ADKUf5NjAP4DJcOmO3+YkpOIiMiPdLSmbCIionaNgZmIiMiPMDATERH5EQZmIiIiP8LATERE5EcYmImIiPwIAzMREZEf+X/qEsn3wAi6HgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# Initial values\n",
    "N = 2000  # maximum population size\n",
    "a = .5 / N  # birth rate\n",
    "b = .5 / N  # death rate\n",
    "nsteps = 10000 # number of iterations\n",
    "\n",
    "# initialize memory with zeros \n",
    "# np.zeros(3) -> array([0., 0., 0.])\n",
    "x = np.zeros(nsteps) \n",
    "# set \n",
    "x[0] = N//4\n",
    "\n",
    "# initialize  memory  with  zeros\n",
    "for t in range(nsteps - 1):\n",
    "    if 0 < x[t] < N:\n",
    "        # Is there a birth?\n",
    "        birth = np.random.rand() <= a * x[t]\n",
    "        # Is there a death?\n",
    "        death = np.random.rand() <= b * x[t]\n",
    "        # We update the population size.\n",
    "        x[t+1] = x[t] + 1 * birth - 1 * death\n",
    "    # The evolution stops if we reach $0$ or $N$.\n",
    "    else:\n",
    "        x[t+1] = x[t]\n",
    "fig, ax = plt.subplots(1, 1, figsize=(8, 4))\n",
    "ax.plot(x, lw=2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interactive Version\n",
    "======"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "fff1a07a6086489c82c7930819726048",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "interactive(children=(FloatSlider(value=0.5, description='a', max=0.95, min=0.05, step=0.01), FloatSlider(valu…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Interactive variant\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "from numpy.random import MT19937\n",
    "from numpy.random import RandomState, SeedSequence\n",
    "import ipywidgets as widgets\n",
    "from ipywidgets import interact\n",
    "\n",
    "# Initial values\n",
    "N = 2000  # maximum population size\n",
    "a = .5  # birth probability\n",
    "b = .5  # death probability\n",
    "nsteps = 10000\n",
    "seed = 100\n",
    "\n",
    "w_a=widgets.FloatSlider(value=a, min=0.05, max=0.95, step=0.01)\n",
    "w_b=widgets.FloatSlider(value=b, min=0.05, max=0.95, step=0.01)\n",
    "w_N=widgets.IntSlider(value=N, min=10, max=5000, step=1)\n",
    "w_nsteps=widgets.IntSlider(value=nsteps, min=1000, max=10000, step=1)\n",
    "w_seed=widgets.IntSlider(value=seed, min=1, max=200)\n",
    "\n",
    "@interact(a=w_a, b=w_b, N=w_N, nsteps=w_nsteps, seed=w_seed)\n",
    "def discrete_markov(a=a, b=b, N=N, nsteps=nsteps, seed=seed):\n",
    "    # b=1-a\n",
    "    x = np.zeros(nsteps)\n",
    "    x[0] = N//4\n",
    "    rs = RandomState(MT19937(SeedSequence(seed)))\n",
    "    #prev (x[t-1]) -> transition -> current (x[t])\n",
    "\n",
    "    for t in range(1, nsteps):\n",
    "        prev = x[t-1]\n",
    "        if 0 < prev < N:\n",
    "            # Is there a birth?\n",
    "            birth = rs.rand() <= a * prev/N\n",
    "            # Is there a death?\n",
    "            death = rs.rand() <= b * prev/N\n",
    "            # We update the population size.\n",
    "            x[t] = prev + 1 * birth - 1 * death\n",
    "        # The evolution stops if we reach $0$ or $N$.\n",
    "        else:\n",
    "            x[t] = prev\n",
    "    display(f\"b: {b}\")\n",
    "    fig, ax = plt.subplots(1, 1, figsize=(8, 4))\n",
    "    ax.plot(x, lw=2)"
   ]
  },
  {
   "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.8.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}