Markov Chains.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "g74VDfctMYox"
   },
   "source": [
    "# Markov Chains\n",
    "**Author:** Gavin DeBrun"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "code_folding": [
     3
    ],
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 1075,
     "status": "ok",
     "timestamp": 1604601300421,
     "user": {
      "displayName": "Gavin Debrun",
      "photoUrl": "",
      "userId": "16932255542580517286"
     },
     "user_tz": 360
    },
    "id": "HdFTVg6hMYo0",
    "outputId": "cbb6c69c-071b-403a-a4b6-ae8d6799ce26"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "datetime.datetime(2021, 2, 28, 15, 14, 9, 788109)"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import random\n",
    "def resetMe(keepList=[]):\n",
    "    ll=%who_ls\n",
    "    keepList=keepList+['resetMe','np','plt','random']\n",
    "    for iiii in keepList:\n",
    "        if iiii in ll:\n",
    "            ll.remove(iiii)\n",
    "    for iiii in ll:\n",
    "        jjjj=\"^\"+iiii+\"$\"\n",
    "        %reset_selective -f {jjjj}\n",
    "    ll=%who_ls\n",
    "    return\n",
    "import datetime;datetime.datetime.now()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "JdMN0yLQMYo5"
   },
   "source": [
    "## Weather Markov Chains"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "4ok7btSiMYo-"
   },
   "source": [
    "### Weather Probabilities"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2GrzbTxfMYo_"
   },
   "source": [
    "|*state*|sunny|cloudy|rainy|\n",
    "|---------|-------|--------|-------|\n",
    "|*sunny*  |  0.9  |  0.075 | 0.025|\n",
    "|*cloudy* | 0.15  |  0.8   | 0.05 |\n",
    "|*rainy*  | 0.25  | 0.25   | 0.5  |"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 380
    },
    "executionInfo": {
     "elapsed": 1051,
     "status": "ok",
     "timestamp": 1604601300423,
     "user": {
      "displayName": "Gavin Debrun",
      "photoUrl": "",
      "userId": "16932255542580517286"
     },
     "user_tz": 360
    },
    "id": "Xoe-L9BHMYpC",
    "outputId": "42f7d61e-0065-4d76-ac22-f99d9c0d2234"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Sunny Probability is: 0.6318\n",
      "Cloudy Probability is: 0.3107\n",
      "Rainy Probability is: 0.0575\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nOzdd5xcVfn48c8zfXvLpm4apJFeKQIJNfTeq6Gj8AULCooIKjZQfoqKiEpR6SCokNBLaCEkEEJI79nU7b1MOb8/zt3J7Ga2ZiebZZ736zWvmbn13Dsz5znlzrlijEEppVTycvV0ApRSSvUsDQRKKZXkNBAopVSS00CglFJJTgOBUkolOQ0ESimV5DQQ9EIiYkRkRBfX3Sgix7Uy70gRWRVvWRH5oYj8rWspThwRmSciX+/pdHyViMgcEXm/B/f/DRHZKSLVIpLXU+lIJhoI9hEnU61zvtw7ReQREUnv6XTFMsa8Z4wZ3cq8XxhjrgYQkWFOMPJ0ZT9ORhN2zkXsY2A7690lIv9qka6TjDGPdSUd7eyry8coIgERKReRY+LM+38i8lwnt/eok5aDY6aNEJGv3J+ARMQL3AfMNsakG2NK4ixzlYisFJEq57f0sohkOPMeFZG793W6ezsNBPvWacaYdGAqMAP4UcsFupq59kIfOT/02Me2nk5UdzDG1ANPA5fHThcRN3AR0JXAVQr0ugyuC9/nfkAA+LKV7c0CfgFcZIzJAA4CntmrRCoNBD3BGLMVmAeMh2hTzw0isgZY40y7RkTWikipiPw3Tmn5ZBFZLyLFInKviLic9Q4UkbdEpMSZ97iIZLdYd4aILBeRMqdmEnDWPUpECuOluUVpfL7zXO6U5Gc56ZwQs3xfpwaU39nzIyK3ishWp8S3SkSOFZETgR8CFzj7/NxZ9h0RaaqpzBGRD5xSd7lzfr7mTN8iIrtim5FE5BQR+UxEKp35d8Uko+UxHuasc6WIrHDO3asiMrSVw3gMOEdEUmOmnYD9zc1r7TjbOC2PAROdjDDeOWvW5Bf7ecXUbq5wjrNMRK4XkRkistQ5V3/cc5PyBxGpcErfx8bMyBKRv4vIdif9dztBruVnUArcxZ4b9ovI70Rkm/P4nTNtFNDUNFkuIm/FOdQZ2ELEZwDGmFJjzGPGmCoRuRa4BPi+85n9z9nfQBF5XkSKRGSDiNzU4jw9JyJPO5/DpyIyKWZ+Zz6j3ssYo4998AA2Asc5rwdjSzw/c94b4HUgF0gBjgGKsTUHP/AHYH7MtgzwtrP8EGA1cLUzbwRwvLNePjZD+12LdCxz0pALfADc7cw7CihsJc13Af9yXg9z0uCJWfYB4Ncx728G/tfKuZgDvN/KvNHAFmBgzL4ObJmGmOXfiTn2OUAIuAJwY0vQm4E/OedjNlAFpMcc7wRs5jwR2Amc2cYxngmsxZZCPdga3YdtfOargUtj3j/Z9Fm0dZxxtvOocyw3NZ0353M28T6rNj6vB7Gl7dlAPfAi0BcYBOwCZrU4j98GvMAFQAWQ68x/EfgLkOasvxC4rsW6/+eco5Q4x/NTYIGzbj7wIbt/C3uc9xbrHgnUAT8BDgf88c5VzHsXsBj4MeADDgDWAyfEnKcgcK5zrLcAG5zXHf6MevujxxOQLA/nh1oNlAObsBlnijPPAMfELPt34J6Y9+nOl3VYzPInxsz/JvBmK/s9E/isRTquj3l/MrDOeX0UXQ8Ehzg/GpfzfhFwfitpasosymMeTWkYgc2UjgO8LdaLpiFm2js0DwRrYuZNcNLZL2ZaCTC5lXT9Dvh/bRzjPOCqmPcuoBYY2sr2fgS85rzOdJad0t5xxtnOo9hA4McGtpPoWiAY1OI8XBDz/nngWzHncRsgMfMXApdhm24aiMngsc1db8esu7md41kHnBzz/gRgY2vnPc76JwH/c7431dg+BXfsuWrxvdzcYv0fAI/EnKcFLT7T7diA0+HPqLc/tGlo3zrTGJNtjBlqjPmmMaYuZt6WmNcDscECAGNMNfaHO6iV5Tc56zQ1yTzlVGcrgX8BfVqkI+66e8MY8zFQA8wSkTHYH9F/21hlgXMumh4HOttZC3wL+wPd5RxLZ9K3M+Z1nbPNltPSAUTkEBF522kyqACuZ89zFWso8HunKaUc224vNP9cYv0DOFpEBmFLnGvN7iaNTh+nMaYB+JnzkLaWbUXL8xD3vDi2GidndDR9T4ZiS8vbY87DX7Cl+yax3694mn2/6eR30BgzzxhzGrZGewY2+FzdyuJDgYFNaXXS+0NsQNsjvcaYCFCIrQXs7Xex19BAsP+I/dFtw36BARCRNCAP2BqzzOCY10OcdQB+6WxrojEmE7iUPTON1tbtSlpjPebs7zLgOWM7TTvNGPOEMeYI7DkwwK/b2W9XPYENVoONMVnYppOmcxVvX1uwTSCxASzFGPNhvI0bYzYD72HbrS/DBobY+a0dZ1seAbKAs1pMrwFi+yP6d2BbbRkkIrHfm6bvyRZsjaBPzDnINMaMi1m2vc+p2febrn0HMcZEjDFvAm/h9LfF2fcWYEOLzyzDGHNyzDLR34PYvraCpvR08TPqdTQQ7J+eAK4Qkcki4sdeJfGxMWZjzDLfE5EcERmMbY9/2pmegdME5ZREvxdn+zeISIGI5GJLR0/HWaYtRUAE294a65/YDOpSWmR6HSUio0XkGOe467El1bAzeycwzPmxdocMoNQYUy/20syLY+bFO8YHgR+IyDgnrVkicl47+3gMuBHbnv1408R2jrNVxpgQtoR6a4tZS4ALRcQrItOxNZC90Re4ydneedh+kbnGmO3Aa8BvRSRTRFxiL1CI24ndiieBH4lIvoj0wbbf/6uddQAQkTNE5ELnuy/O5zYL2+cA9jsS+5ktBCqdTt8UEXGLyHgRmRGzzDQROVvsFU7fwga6BV39jHojDQT7IaeUcwe23XY7cCBwYYvF/oPtBFsCvIztVwDbiTYV27n3MvDvOLt4AvtjXu88OnVZojGmFvg58IFT3T7UmV4IfIotOb3XzmYOkz3/RzAD2w7+K2xn+Q5shvRDZ51nnecSEfm0M2luxTeBn4pIFTYzil6GGO8YjTEvYEuETznNbsuw7dVteQ7IwfbhbI+Z3tZxtudJ7Pci1h3Y70kZ9jvwRAe31ZqPgZFO+n4OnGt2X9N/Obbjdbmzv+eAAZ3Y9t3YPqSlwBfY70xHv4NlwDXYq+uamj7vNcY0Bdm/A2Odz+xFY0wYOA2YjO0ELgb+hq1VNfkPtkO8DFtzO9sYE2TvPqNeRZo3Ayq1d0TkYWCbMWaP/0gotb8Re8nwCGPMpT2dlp6ULH9eUvuAiAwDzgam9GxKlFKdoU1DqluIyM+wTSX3GmM29HR6lFIdp01DSimV5LRGoJRSSa7X9RH06dPHDBs2rKeToZRSvcrixYuLjTFxx/7qdYFg2LBhLFq0qKeToZRSvYqIbGptnjYNKaVUktNAoJRSSU4DgVJKJble10eglPpqCgaDFBYWUl/fpbEKlSMQCFBQUIDX6+3wOhoIlFL7hcLCQjIyMhg2bBjNBz5VHWWMoaSkhMLCQoYPH97h9RLWNCQiD4u9NeCyVuaLiNwv9naMS0VkaqLSopTa/9XX15OXl6dBYC+ICHl5eZ2uVSWyj+BR4MQ25p+EHd1wJHAt8OcEpkUp1QtoENh7XTmHCQsExpj52Ds4teYM4B/GWgBki0hnhrLtlFU7qnjw3XXUB7+Sw4krpVSX9eRVQ4Nofku7Qlq55Z+IXCsii0RkUVFRUZd29tqXO/jVvJUs2VLepfWVUl99brebyZMnM378eM477zxqa2sBqKurY9asWYTDPVuQPO644ygrK+v27fZkIIhXf4k7Ap4x5iFjzHRjzPT8/Lj/kG7XwcNzAYhEdJA9pVR8KSkpLFmyhGXLluHz+XjwwQcBePjhhzn77LNxu909mr7LLruMBx54oNu325OBoJDm986N3ic0kTQMKKU64sgjj2Tt2rUAPP7445xxxhkAbN++nZkzZ0ZrDu+9Z2/Gl56eHl33ueeeY86cOQDMmTOHm266ia997WsccMABPPfccwC88847HHXUUZx77rmMGTOGSy65BGMMb775JmedtfuW1K+//jpnn302AKeffjpPPvlktx9rT14++l/gRhF5CjgEqGhxK79upZ1QSvUeP/nflyzfVtmt2xw7MJM7TxvXoWVDoRDz5s3jxBNPpLGxkfXr19M02OUTTzzBCSecwO233044HI42H7Vl+/btvP/++6xcuZLTTz+dc8+1t5T+7LPP+PLLLxk4cCCHH344H3zwAccccww33HADRUVF5Ofn88gjj3DFFVcAkJOTQ0NDAyUlJeTl5XXtRMSRyMtHnwQ+AkaLSKGIXCUi14vI9c4ic7H3y10L/BV7/1illOoxdXV1TJ48menTpzNkyBCuuuoqiouLyc7Oji4zY8YMHnnkEe666y6++OILMjIy2t3umWeeicvlYuzYsezcuTM6/eCDD6agoACXy8XkyZPZuHEjIsJll13Gv/71L8rLy/noo4846aTdt8bu27cv27Z1b+NJwmoExpiL2plvgBsStf/W97uv96iU6qyOlty7W1MfQctpsdflz5w5k/nz5/Pyyy9z2WWX8b3vfY/LL7+8WatDy+v4/X5/9HXszcBip7vdbkKhEABXXHEFp512GoFAgPPOOw+PZ3dWXV9fT0pKyl4eaXNJM9aQtgwppboiJyeHcDgczdw3bdpE3759ueaaa7jqqqv49NNPAejXrx8rVqwgEonwwgsv7NU+Bw4cyMCBA7n77rujfQ1gg8iOHTvo7nuyJN0QE0a7i5VSnTR79mzef/99jjvuON555x3uvfdevF4v6enp/OMf/wDgV7/6FaeeeiqDBw9m/PjxVFdX79U+L7nkEoqKihg7dmx02uLFizn00EOb1RC6Q6+7Z/H06dNNV25Ms2hjKec++BH/vOpgjhzZtUtQlVKJs2LFCg466KCeTkZcn332Gffddx///Oc/99k+b7zxRqZMmcJVV10VnXbzzTdz+umnc+yxx7a5brxzKSKLjTHT4y2fdDUCpZTqrClTpnD00UcTDof3yX8Jpk2bRlpaGr/97W+bTR8/fny7QaArki4Q9LIKkFJqP3HllVfus30tXrw47vRrrrkmIfvTzmKllEpySRMImmiFQCmlmkuiQKBVAqWUiieJAoFSSql4ki4Q9LbLZZVS+87eDEM9Z86c6IBye2vYsGEUFxe3Or+oqIgTT2zrvl+dkzSBQDuLlVLt2d+HoW6Sn5/PgAED+OCDD7ple0kTCJpofUAp1RGtDUMNcM899zBhwgQmTZrEbbfdtse6b775JlOmTGHChAlceeWVNDQ0AM1L+osWLeKoo44CoKSkhNmzZzNlyhSuu+66aMvFHXfcwe9///vodm+//Xbuv/9+wA5k9/jjj3fLsSbN/wi0QqBULzLvNtjxRfdus/8EOOlXHVq0rWGo582bx4svvsjHH39MamoqpaXN78hbX1/PnDlzePPNNxk1ahSXX345f/7zn/nWt77V6v5+8pOfcMQRR/DjH/+Yl19+mYceegiAq666irPPPpubb76ZSCTCU089xcKFCwGYPn06P/rRj7pwIvaUdDUCpZRqTUeGoX7jjTe44oorSE1NBSA3N7fZNlatWsXw4cMZNWoUAF//+teZP39+m/udP38+l156KQCnnHIKOTk5gK1B5OXl8dlnn/Haa68xZcqU6H0IunM46qSpEURp25BS+78Olty7W0eGoTbGtHmjq7YuSPF4PEQiEWDPoapb2+bVV1/No48+yo4dO5r9u7k7h6NOmhqB3qFMKdUVLYehnj17Ng8//HD0iqKWTUNjxoxh48aN0f6Ff/7zn8yaNQuwJfym4SOef/756DozZ86MtvfPmzev2Q3qzzrrLF555RU++eQTTjjhhOj01atXM378+G45xqQJBE10GGqlVGc1DUMNcOKJJ3L66aczffp0Jk+ezG9+85tmywYCAR555BHOO+88JkyYgMvl4vrr7Y0Z77zzTm6++WaOPPLIZlcg3XnnncyfP5+pU6fy2muvMWTIkOg8n8/H0Ucfzfnnn99snbfffptTTjmlew7QGNOrHtOmTTNdsWRzmRl660vmzRU7urS+Uiqxli9f3tNJaNWnn35qLr300h7ZdzgcNpMmTTKrV69uNv3II480paWlcdeJdy6BRaaVfDXpagRKKdVZscNQ70vLly9nxIgRHHvssYwcOTI6vaioiO985zvRTuW9lXSdxfrHYqVUV+zLYaibjB07lvXr1+8xPT8/nzPPPLPb9pM0NQLtK1ZKqfiSJhA00RqBUko1l3SBQCmlVHNJEwhEB5lQSqm4kiYQNNGWIaVUa/ZmGOp4Tj75ZMrLy7uUlgsvvJA1a9Z0ad3OSppAoJ3FSqn2dHYY6vYCw9y5c5uNU9QZ3/jGN7jnnnu6tG5nJU0gUEqpzmhtGOp33nmHo48+mosvvpgJEyYAdkjoadOmMW7cuOjIobB72OmNGzdy0EEHcc011zBu3Dhmz55NXV0d69atY+rUqdHl16xZw7Rp06L7f+ONNwiFQgk/1iT8H4E2Dim1v/v1wl+zsnRlt25zTO4Ybj341g4t29Yw1AALFy5k2bJlDB8+HLA1htzcXOrq6pgxYwbnnHNOdJTQJmvWrOHJJ5/kr3/9K+effz7PP/88l156KVlZWSxZsoTJkyfzyCOPMGfOHABcLhcjRozg888/jwaHRNEagVJKOToyDDXAwQcfHA0CAPfffz+TJk3i0EMPZcuWLXHb9ocPH87kyZMBmDZtGhs3bgTs6KKPPPII4XCYp59+mosvvji6TncONd2W5KsR9HQClFLt6mjJvbt1ZBhqgLS0tOjrd955hzfeeIOPPvqI1NRUjjrqqD2WB/D7/dHXbreburo6AM455xx+8pOfcMwxxzBt2rRmNYnuHGq6LUlTI9DOYqVUV7QchrqliooKcnJySE1NZeXKlSxYsKBT2w8EApxwwgl84xvf4Iorrmg2b/Xq1YwbN67Lae+ohAYCETlRRFaJyFoR2ePGniKSJSL/E5HPReRLEbki3naUUqonxQ5D3dKJJ55IKBRi4sSJ3HHHHRx66KGd3v4ll1yCiDB79uzotJ07d5KSksKAAQO6nO6OSljTkIi4gT8BxwOFwCci8l9jzPKYxW4AlhtjThORfGCViDxujGlMVLq0r1gp1Zrq6uq402+88Ubuu+8+jjvuOI466qjoTefBNvnMmzcv7npN/QB9+vRh2bJl0em33HJLs+Xef/99rrzyymaXpz7xxBNcd911XTySzklkH8HBwFpjzHoAEXkKOAOIDQQGyBB7+7B0oBRIyLVS+s9ipVRXxQ5D3fK/BHvrrLPOYt26dbz11lvNpmdnZ3PZZZd1675ak8hAMAjYEvO+EDikxTJ/BP4LbAMygAuMMZGWGxKRa4FrgWZ37ukarRIopTovUcNQv/DCC3Gnt+wvSKRE9hHEK4K3zIVPAJYAA4HJwB9FJHOPlYx5yBgz3RgzPT8/v2uJ0QqBUvs9/Z/P3uvKOUxkICgEBse8L8CW/GNdAfzbuZPaWmADMCaBaVJK7acCgQAlJSUaDPaCMYaSkhICgUCn1ktk09AnwEgRGQ5sBS4ELm6xzGbgWOA9EekHjAb2vB1PN9LvmFL7p4KCAgoLCykqKurppPRqgUCAgoKCTq2TsEBgjAmJyI3Aq4AbeNgY86WIXO/MfxD4GfCoiHyBbUq61RhTnIj0aNOQUvs3r9fb7N+6at9J6D+LjTFzgbktpj0Y83obMLvleglN077cmVJK9QLJ889ivXxUKaXiSppAoJRSKr6kCwTaWayUUs0lTSDQzmKllIovaQJBE6PdxUop1UzSBAKtECilVHxJEwiUUkrFl3SBQDuLlVKquaQJBNpZrJRS8SVNIGiiFQKllGouiQKBVgmUUiqeJAoE1vef+7ynk6CUUvuVpAsE9cE9boCmlFJJLWkCgXYWK6VUfEkTCJRSSsWXNIFAKwRKKRVf0gQCpZRS8WkgUEqpJJc0gUC0t1gppeJKmkCglFIqvqQJBFof6D5GR+5T6islaQKB2juRiGH+6iKeWriZcXe+ymV//7jVgBCJGN5YvpOKuuA+TqVSqis8PZ2AZFBW00hDKEL/rABgM8oXl2xle0U9Nxw9IrrM795YzU3HjuTRDzfSNzPAJQcPweXaXZd5f00xH64r5txpBRyQnx6dvnJHJR+sLeGlpdv408VTGZidEp3XGIrgdgluZzurd1Zx/b8W8+tzJjJjWC4Au6rqCYYNOaleUn2e6Hpet7C+uIZ5X2znd2+sIRTZnfG/t6aYa/+5mHVF1Rw1qi9Th2ZT1xhm4YZSnl1cCED/zADPXn8YJTWNTB6cHffc1AfD7KioZ1iftL0+z0qprpHeVs2fPn26WbRoUafX21RSw6x73wFg469O6eZUxWeMQUQY++NXqG0M7zHf6xZe+ObhfLKxlLdXFTF/ddEey7x80xFk+L3MvPftZtPPmjKIX5w1gfteX8Vf39vQbN7/bjyCt1ft4r7XV0en3XXaWBpCEX45b2V02jPXHUYoEuHiv34MQIbfw1+/Pp0vCiv4+dwVcY/p1IkD+MXZE/j6wwv5bHN5q8c+ODeFLaV10fej+2Vw+uSBnDCuP/2zAry8dBv/WbKNJVvKqQuGOXZMXxrDhpuPHcG0obmtbrczms6/UgpEZLExZnrceckcCFbvrKKwrJZjxvTrziRGnfvnD1m0qSwh2wbIS/NRUtMYfT+8TxobimvaXe+8aQXRUntHnD1lEDcdO7JZqX1XVT33vbaakyYM4K/z17N6ZxXHjOnL0WP6cvxB/XC5hD+9vZa3V+4iFDEs2bJn0OiT7mPswCzW7apma3kdLoGIge+dMJqrjhhOwOvucBrB1mLWFVWzaFMZjy/YxModVVw4YzA/O3M8XnfXW0FD4QhV9SFy0nzRafXBMG6XsGZnNR+tL+GgARmMHZBJdqqvjS0p1XM0EACbS2qjper1vzgZl0sYdtvLwN7VEMIRw9cfXshZUwZxzrQCAF5ZtoONJTX8Kqb0neH3cOqkAdQ0hFm8qYzTJg3kwXfXNdvWgh8cS2V9kFH9MiipbuCsBz5kc2ktADNH5fPw16fjcbvYVl7H1371FgDD8lL548VTGT8oC2MMFzy0gIUbSpk8OJtvHz+KI0b0Yfm2Sk774/tMHZLNk9ceit/j5n+fb+P/nvwMn9vFhz84hj7pfv7y7jp+OW8l50wt4J5zJ/LZ5jImFGTh93QuQ44nFI7w7OJClhZWsKuynssOG8qsUfmICMFwhE82lDI4N5Wbn/qMTzeXk+H38J3Zozh7SgF+r4sV2yv5cF0Jb6/cxeDcVC46eAj9MwM8s2gLq3ZWEQpHWLypjMr6EGBrJHlpfpZsKWdIbiqXHzaU/Aw/dY02Az9nakGzZrfGUITNpbWU1zaycGMpW8vqyEzx4nUJz3+6la3ldQzKTkEEiqoaCEdMs6YysDW8Sw4ZyvhBWWworqZPup+3Vu5iyeZyjj2oL9fOPJAF60v4cF0JKT43o/ul4/O4yE7xMW1YDoOyU5oFv8r6IAvWlVAfipDh91DVEGJLaS0DsgL0zwwwJC+VgpzU6PLltY2ICFkpXsDWiBpCkU4HVLDfa5foZdfdJRiOsKG4hvx0P+kBT7OCSSgcIWLA50lsl60GApoHgieuPoSvjejTLYHgw7XFXPw327Ty729+jbMf+HCPZaYOyebf3zx8j+llNY2sLapmyuBs6oJhMgLeuPsIhiN7lGifXLiZt1fu4oFLpuKJmReJGKQX/4AbQxFeXLKVv7y7jnVFe9ZuctN8lMbUgpoyq4yAh8MOyGPS4GxG989g5sh83C7hpaXbuPfVVWwqqW22ndMnDeQXZ0/g1WU7eGbRFj7ZWEqklZ/CwcNzmTEsh8+3VBDwuhAR8tJ8+DwuRvZN57AD+7B8eyXvrNzFfz7fRjhmQ/0zA4zqn8GC9SU0huzIt3lpPirqgnsEEo9LuPTQoYQiEYqrGnl71S4aQq2Plut2CaP7ZRDw2s//88IKAMYPzCQ/I8DybRUUVTcwKDuFhlAElwipPjej+2cgIvTL8BMMO4FCwCXCjop6Pl5fwvbKegZmpTBhUBbVDSF2VdXj87jomxFga1kdpbWNVNeHGNYnDbcLRvfLZNzATNwuYWTfdMYXZJEZ832ORAyN4Qh+j4vCsjrW7qrG4xY2ldTy5bZK+mX66Z8ZINXvIcPvYWNJDbuqGhiRn05FXZBxAzMZlJNCeW2QoqoGMlM8rNtVQ1F1Ax6XLUzUBcOk+T2ket2sLapmQFYKkYhhZL8MJhRk0S/DTyhiCHjdFFc3sGhjKdUNYQbnpFDn1PAiBpZsLqestpFwxDAkNxWDYWJBNi4RSmsaOWR4Li6XRL+H5bWNbC6tpaiqAb/HRWV9iHVF1Wwtq2NreR1ul7C9oj76+af7PYzpn0GKz82andXsqKzH73FRkGMLL+V1jWwqqSUj4KFPup/cNB9ZKV4aQxFOGNef82cMbvU70RYNBDQPBOdOK+Decycy/Adzgd01hK745dwV/GX+egD6ZfrZWdkQnTdpcDZ3nHIQo/tntJrJq/iMMby3pph3VxexvaKO8YOyOHp0X8b0z6C4upF/fLSRxlCEOYcPo3+m7YRvLfgZY1i9s5ri6gaMgU82lvL7N9dE5w/OTWHmyHyG5aXRN9PPtKE55DqZdcTAoJjO9/bsqqqnsi5In3Q/RVUNDO+ThsftYu2uKt5bU8ysUfkM75NGKGKobbD9RhtKavhscxkL1pfw6pc78bldBLwuTp88kJPHD8DlEqrqQ/TPDNA308+W0lqqGkK8t7qYL7aWU1LTSGbAy6EH5GEwfLC2mO3l9UwenE3fzAA7K+vJCHgIRww1DSGWFlYQMYay2iA+t4vGcASPS4gYQ3aqj8MOzGNIbiqrd1SxckcVGQEPg3NT2VFRT10wzLC8NDJTPKT63KzeUY3XIyzeVNZsiPeA10VBTioRYwh43GworiEUiZDidUdrbU3S/R6qG5pP6wqR3fckT/W54/bLgS15N7YRYJzhekYAACAASURBVJvS7xahppVttKd/ZoCCnBQG5aQQihgKslMY0TedkppGNhbXsKG4hsr6EGP6ZzAkN5WKuiDbyusorwuSleJlcE4q1Q1BymqD7Kysj56fa488gAsPHtKlNGkgALaU1nLkPbs7XBfefiwH//xNAN74zixG9E1vbdVm6hrDVNYHSfW5+dlLy3lmUWGzLyDAxIIslhZW8NEPjmFAVsczEbXvzF9dxDuripg5qk+0iWp/UFLdQI7Tz9DVwklHNDUbeVyCwf7PJhQx+NyuLu3XGENhmb04YPXOKt5dXcTm0tpohjuqXwYRY2gMRRg7MJMD89Nxu8RmltkplNUGKattpD4YprYxzNDcVNIDHnZU1BPwullaWEFZrQ14uWk+ymsbGTMgk74xtRq/x0V1Q4j6YIQ+6bb/LM3nYcmWctYVVbOjoh6v20VVfZA+GX5mDMshze9hZ2UDqT434YghHDFMcGozwXCE2oYwjeEIn24uI+B1k+pzs3BDKT63i9w0H6FIhNw0P4NzU+ifGaAxFCHgczerDXWnvbkAYq8DgYicCsw1xvT4XV26KxC8f+vRHPHr3e/n3nQkYwdmtrudpuakI0f24b01xQBccfgwlmwp57PN5aT53Hz50xP1ihWl1H6lrUDQ0d6JC4E1InKPiBzUiR2fKCKrRGStiNzWyjJHicgSEflSRN7t6Lb3VmwQADj5/vfaXae8dnfbdFMQAHuVyyWHDAXgrVuOAnpvG71SKvl06A9lxphLRSQTuAh4REQM8AjwpDGmKt46IuIG/gQcDxQCn4jIf40xy2OWyQYeAE40xmwWkb57dzjdLxSOcNf/vuSig4dw6/NLm81zCXx6x/Gk+jycM3UQZ0weuFeXKSqlVE/o8D+LjTGVIvI8kAJ8CzgL+J6I3G+M+UOcVQ4G1hpj1gOIyFPAGcDymGUuBv5tjNns7GNX1w4jcS58aAGLNpWxYH1ps6tVAH534ZTodeMigtettQClVO/ToeKriJwuIi8AbwFe4GBjzEnAJOCWVlYbBGyJeV/oTIs1CsgRkXdEZLGIXN7K/q8VkUUisqioaM9/33bsGDq3fFV9kO8+83n0D2E7KupxiXD21N2HcNrEAV1Ki1JK7U86WiM4F/h/xpj5sRONMbUicmUr68TLelv2THuAacCx2JrGRyKywBizutlKxjwEPAS2s7iDae60hlAYY2DMHa+Qn2Ev/2tS3RCiuiHEpIJsBmQFMEb7AZRSXw0dDQTbWwYBEfm1MeZWY8ybraxTCMT+86EA2BZnmWJjTA1QIyLzsbWM1XSzjmTa28rriThXUcUGgViTBme3OoCaUkr1Rh3t2Tw+zrST2lnnE2CkiAwXER/2yqP/tljmP8CRIuIRkVTgECD+aGf7QEMozOV/X7jH9NjmoIMGZOzLJCmlVMK1GQhE5Bsi8gUwRkSWxjw2AEvbWtcYEwJuBF7FZu7PGGO+FJHrReR6Z5kVwCvOthYCfzPGLNv7w+q4SQVZ0dcVtUG2ltftscwZk3cHgu4Yd0cppfYn7TUNPQHMA34JxP4PoMoYU9rexo0xc4G5LaY92OL9vcC9HUrtXmitYegfVx7CS19s4/YXlnHBQwui0686Yjh/f98O7zxzZB/7PCo/0clUSql9rr1AYIwxG0XkhpYzRCS3I8Fgf5eV6uWIEX2aTTtz8kDuOHUsPzrloGin8Gd3HE+qX2sDSqmvno7UCE4FFkN0SJImBjggQenaJ66baZPfNGxvk6aSv4hELzuNHYteKaW+StoMBMaYU53n4fsmOYkT76KhI5wmn3R/89OQ0oXx25VSqrdqMxCIyNS25htjPu3e5Ow7A7ICHDnSlvw9bhenTRrI/z63V7ceNKD9weeUUuqror2mod+2Mc8Ax3RjWhJKEK5xv8TrkWlsNAPIz/A3m/+Hi6ZEA4HeSF0plUzaaxo6el8lJNGkoYrbvU9wtZnLIQ0P4IrTVrSvbmqvlFL7k/aaho4xxrwlImfHm2+M+XdikpUI9h/DKdh/DCfwnh9KKdWrtNc0NAs70NxpceYZoNcEgpYVALdGAqWUAtpvGrrTeb5i3yRn34nXNKSUUsmoo8NQ54nI/SLyqTNc9O9FJC/RietOLbN9DQRKKWV1dNC5p4Ai4BzskNRFwNOJStS+oE1DSilldXQY6lxjzM9i3t8tImcmIkH7ilYIlFLK6miN4G0RuVBEXM7jfODlRCas27XI+AtyUnsmHUoptZ9p7/LRKnaPMfQd4F/OLBdQDdyZ0NQlyJ8unsqxB/Xt6WQopdR+ob2rhr4yd2ER5waXApyi9xpWSqmojvYRICI5wEgg0DSt5e0rewPZ47bJSimV3DoUCETkauBm7H2HlwCHAh/Ri8YaavpnsQYCpZRqrqOdxTcDM4BNzvhDU7CXkPYaTVcJ6cVCSinVXEcDQb0xph5ARPzGmJXA6MQlKwGMrQm4iPRwQpRSav/S0T6CQhHJBl4EXheRMmBb4pLV/aTFs1JKKatDgcAYc5bz8i4ReRvIAl5JWKoSQvsIlFIqns5cNTQVOAKbo35gjGlMWKoSSgOBUkrF6uigcz8GHgPygD7AIyLyo0QmrLtp05BSSsXX0RrBRcCUmA7jXwGfAncnKmGJok1DSinVXEevGtpIzB/JAD+wrttTk0BNAcClgUAppZppb6yhP2Ab1RuAL0Xkdef98cD7iU9e93OJBgKllIrVXtPQIud5MfBCzPR3EpKaRDIaAJRSKp72Bp17rOm1iPiAUc7bVcaYYCIT1t30/gNKKRVfR68aOgpYA/wJeABYLSIzE5iuBIipEexa0XPJUEqp/UxHO4t/C8w2xswyxswETgD+X+KS1f0k9sLRBw6Fso09lhallNqfdDQQeI0xq5reGGNWA97EJClBWvYRfPFsz6RDKaX2Mx0NBItF5O8icpTz+Cu2A7lNInKiiKwSkbUiclsby80QkbCInNvRhO817TxWSimg44HgeuBL4CbskNTLnWmtEhE3tk/hJGAscJGIjG1luV8Dr3Y82d0gEt6nu1NKqf1Vu/8sFhEXsNgYMx64rxPbPhhYa4xZ72znKeAMbBCJ9X/A89j7HSRQixpAJJTY3SmlVC/Rbo3AGBMBPheRIZ3c9iBgS8z7QmdalIgMAs4CHmxrQyJyrYgsEpFFRUXddD+cloHgLzPh39fBrpXds32llOolOjrW0ADsP4sXAjVNE40xp7exTrwr91s2zP8OuNUYE5Y2LvQ3xjwEPAQwffr0rjXut+wTiA0E4RBs/9w+lj4FPy4Fl7tLu1FKqd6mo4HgJ13YdiEwOOZ9AXvezGY68JQTBPoAJ4tIyBjzYhf21zmxgaChsvm8tW/AqBMSngSllNoftDfWUADbKTwC+AL4uzGmo43rnwAjRWQ4sBW4ELg4dgFjzPCYfT0KvJSwINCywhEbCN77bfN5tSUJSYJSSu2P2usjeAxbav8Ce/XPb9tefDcnYNyIvRpoBfCMMeZLEbleRNq84mifWPIkRCLwyd/goz82n/efG3smTUop1QPaaxoaa4yZACAifwcWdmbjxpi5wNwW0+J2DBtj5nRm253Wso8gWAML/wKvxPl7gwlDXRmk5CQ0SUoptT9or0YQHViuE01Cvccrt4HEdArfUbz7dcn6fZ8epZTqAe0FgkkiUuk8qoCJTa9FpLKddfcrGf5WDtWbYp/FBW4vHPEd+75kbfPlVs2zVxcppdRXTJuBwBjjNsZkOo8MY4wn5nXmvkpkd2j18lSP3z6biH0+6jYbFBb+BX7WF+7KgsLF8OSF8Pi+GwFDKaX2lY5ePtr7tTa2UMsrhDx+GxS2xgylNP9e+7x9SWLSppRSPaijYw0lt9Xz7HNDlQ5Wp5T6ykmeGkF7N60fObv9TURCtu9gzWu2k/nQnr8KViml9pbWCJpc8K/dr3+4HcaeGX+5tW/Cqz+EV27VEUyVUl8JyRMI2mrSufad3Z3GAL5UOP8xuGEhpPez01LzIG8ErH5l93KFnyQipUoptU8lTyBoS+6B8afnj4YbPoZDvgGXvQgHHgvr3949/+ET4C+z7GWl5Zvh7V/qJaZKqV5H+wgAfOmtz0vJgZN+ZV9XHGUvK421fQm8cSes+B+Ub4L+E+CgU22zkY5gqvZnkQiEGyBYZ/u/GqvBnwVpeXZ+QxXUV9rfgC8VakqgvhwyB9kadMtLshtr7f9yIiHY+ilU77BX4Pkz7TZEbN9a37Hg7sGsp6l1oI0Rj5NNEgWCNrg6WDEadvju12NOhZUv2dexYxV9+g8I1cPzV8FNSyB3OKoLImEbXHcth+whUHCwPZfubrhVtjGw8CHYvhRGHGsvFPC3URjoLu0VDiq2wvp3IHuwXXbIYeANtL78zuX2MufS9eDywICJUF8BK+dC0Qonw/VBxRZw+wFjv5u7VtqMuKEKEDukSqzsoZCaC7tW2OVdHltYqq8gWqDyZUC/sVC1w86vLbbzA1nOcm3wpYM/A7yptrm13zgoWmUDR6jBLjNgMmT0g2A9pObA1s+gcqsNNP4Me1z5o8GXBg3VznkVqCu1ASlnqN1WqAFK1tgae+V2+/2pLbHBacBE8ASgotD+dygSsoGr31h7HL50+3C5YcN79rMIZMG2zyBjoN2Wyw1p+ZA50D57/HZbWYPt+7R8+z41d/d3t67cfmbVOyFnmD3fjdVQtskWJrMKIG+knVZXas9xY41Nc98x9vfQzcT0ssshp0+fbhYtWtT5FUs3wP2T7esjb4H3frN73l3tfHFjvfAN+wGd9ygUr4bVr9oaQWvOeACGz7Q/btW+UCN8+hi8/zuoLGw+L3sIHHA07FhqS6q+NJhwrs1QsgbDyOPjZ7T1FVC+xf6oAL541gZst9+WiMUNuQfYDMblgb4HQcYAm5lGwna7I2fbH296f7sPE7HbzRlmf+RgM7OP/mgz2JRcm8EOPQIGTYW3fw5fvmgz9+EzoWIzINBnpL0PxpZPnGkxvGkw9Gs2My7fZEvVgSwo22ibIKu2E7emm5YPfUbZY64ttrVUE7HH6XJD3oH2tS/NloqbMj23zx7T9s/tueoz2p6X6p02M0rPtxlUzS77fueX9jhdLghk2360hko7rd84m8GF6uw5rCvbndFuWWiXC9VD0WooXmXPd/YQm9GL2wa4ulKbUYfqIWsI9BlhM/b6SrvdknX2+MXllPKNrdH40uy5aRotIHuIDTgZA2waUnN3H2dDlW0aNmHnP0QGdi6z04N19oGBfhPsuaveaYNUbQmEG+28mhK7v5YBNZbLa4NTYw1U72p72bYcfjMc/9MurSoii40x0+POS8pA0JTxf/pPG91P7cwdOFso2wS/n2hfX/wMPHH+nst4UuDcv8PIE5pXiR8/35Zy+oy070/+Lez4HA48xr4Ph+yXzZe65zYjYQjW2swrVkOVffam7s4Uw0G7rDdt9/4ba21/x+BDIK2PnbZyLrx1NwyeAcfeaYfnPvSbkDWIblFban+82z6Fg06HzAF2en2FLf0PmAwvfxe2LLCZ2TF32Ey4ZC1s+tDeJ6LwE5vh5B5gM4TiVbu3n5JrM8DsITD5YpsZvfETWPI4e2SYh94As38Gmz6ww4fsXLa7JFu81g5KCHZfbQ5LLjYtldts5uQJ2AsMKrfZoBKqcxZzwcQLYNsSG2BcHvsZNVTaTHjo4bYkOPECuz8ThnVv2eM2xhYkTMRm7tmDbeadMwzGn2s/v0A2bF1kv2sDp9jM2Ri7zv7eRBlqsEEotqnGGFvg8qXbc+TP3LMpJxJxvtcpRMeZb6rdN1RjA0OL30dnRcI28w60M5BCOGS/P43V9n3lVqgphpoi+xlUboXiNTZIZQ6EQdMgra+tqVRsscefe4D9bCu32WX96ZDax34HU3LseUjNs4G8CzQQQPxA0F0+f9p+GceeDmvegHd+YTPQ569qfZ0Bk2yJJB5vKnz9JfibExBmft9mbLnDd7dv/mEalK6DMx+EEcfZ0lrJOvjD1LbTetFTMOpE+PPXbLMLwNl/hXVvw+dPxF9n2JF26I1hR+yetuhheOnbey57wFEw+mSYdoUtQZeut6Wyz/7ljNfkVP0Rew7S+tj25KbSOsCxP4av3bxnO7Ix9seVmrs7c6sotCXk2lJY9rztsynbxO6MX2DsGfYceQI2g3V5YPw5rbcRN9bYH3VKjl1n66ew4R3ILLA/xroy+3m7PLaav+kDyOhvayWHXGdfG2MfWxbY4HXgMbZkDrY07U21mVxdmc3Ue7LNXCUFDQRgM6T7p9jX3R0I2trn9s/h2Tnds712S6fd4JT7bKnlheug33ibKZZt2D0/kAVTLtvzHg4tNVXpY+UMh/Fnw6DptvS6+jW77XAjHH27LQlNusCWlvZGsH53/8KYU6Ag7ndfqaSigQB6JhA0WfuGrQIuftSWBofPgi+esfOa7o9cWwof3g+Hf8v2X3z4Bzv/pHth3vfib/fqt+DZr9uqZZOpX4fT77evQ422TTMlx1Yza0vh77Nt5xnAD53mi2fnwKq5cP37u0utsXYsg3m3wqb3m0+/cZFtc80fY0vX4aA9zuX/sX0nfUbAsJlQtQ3GnW3baD2+5tuIhG3pu6mdXSmVEBoIoHmzyb4OBPGUbYTqItse3xHhkB3aongVTL/SlsyblKyzJXcTtu3sibosLhy0GffHf7Z9C0O/lpj9KKW6XVuBQBsme0rOMPvoKLcHxpwMnLznvC52HnWa2+vcsyFO34BSqtfSfxYrpVSS00CglFJJLnkCQS/rC1FKqX0leQKBUkqpuJIoEGiNQCml4kmiQKCUUiqe5AkETX0EJ93Ts+lQSqn9TPIEgiapeT2dAqWU2q8kUSDQPgKllIoniQKBQ+9KpJRSzSRPIND/ESilVFzJEwiitEaglFKxEhoIROREEVklImtF5LY48y8RkaXO40MRmZS41GiNQCml4klYIBARN/An4CRgLHCRiIxtsdgGYJYxZiLwM+ChRKUnJmEJ34VSSvUmiawRHAysNcasN8Y0Ak8BZ8QuYIz50BhT5rxdABQkLDXaR6CUUnElMhAMAmJunUWhM601VwHz4s0QkWtFZJGILCoqKtrLZGmNQCmlYiUyEMTLceMWy0XkaGwguDXefGPMQ8aY6caY6fn5+d2YRKWUUom8Q1khMDjmfQGwreVCIjIR+BtwkjEmgXdm16YhpZSKJ5E1gk+AkSIyXER8wIXAf2MXEJEhwL+By4wxqxOYltid7pPdKKVUb5GwGoExJiQiNwKvAm7gYWPMlyJyvTP/QeDHQB7wgNgMOtTazZW7IUEJ2axSSvV2Cb15vTFmLjC3xbQHY15fDVydyDTsSWsESikVK4n+Waw1AqWUiieJAoFD+wiUUqqZ5AkE2keglFJxJU8giNIagVJKxUqiQKA1AqWUiieJAoFD+wiUUqqZ5AkE2keglFJxJU8giNIagVJKxUqiQKA1AqWUiieJAoFD+wiUUqqZ5AkE2keglFJxJU0gCEXCfKdvH1bW7DEStlJKJbWkCQTra7fyeloqP1j7VE8nRSml9itJEwg84gYgZMI9nBKllNq/JE0gcDuHWhdp5Pvvfp9NlZt6OEVKKbV/SJ5A4NQIdjZWMG/jPH7w3g96OEVKKbV/SJpA4JHmh7qlaksPpUQppfYvSRMIXC3+P1DeUN5DKVFKqf1L0gSCeKoaq+JO31Gzg6VFS4mYCKFICICt1VupDdbuy+QppdQ+kdB7Fu9X4vyh7Lz/nccr57yyx/Tjnzu+2fu3znuL0184ncZII7MKZvHNyd9kbN7YhCVVKaX2peQJBHFsrd7aoeWOefaY6Ot3C9/l3cJ3Abhg9AV8a+q3SPelE4qEeHXjq8wvnM/1k65neNbwbknjqtJVlDWUMaXvFPxuPwAN4Qa+KPoCv9vP6NzR+Ny+btmXUio5JU0gMK0MOldaX0puIJfy+nLuXXQvE/pMINOXyayCWYzNG4vf4+enH/0UgFum30JtqJaHv3iY+nA9T696mrc2v0VRXVGzbc7dMJdZBbM4evDRLC9ZzhGDjmBN+RpeWv8SRw8+mhfWvMCM/jPwu/2MzRvLzIKZDEwfiMdlP4715eu555N7SPWm8vqm16PbHZg2kJE5I6OBCMDn8nFQ3kFMzJ/IyOyRVDRUICKMzh3N6JzR5ARyKK8vZ+GOhTSEG+ib2pcBaQMoqiuipK4EEWFIxhCGZg4l4Al092knHAnjEhd1oTpEhIA7gHRyvCfj1OaqglV4xIPH5eGVja8QMRGyfFn4PX6K64qpbqwmy59Flj+LdG861cFqCqsK2V6znfUV62kMNyIipHpSyfRl0hhuxO/xk+PPIdufTU7APqf70hGENG8aGb4MUjwp1AZrqQnWUNlYSZo3jXRfOhm+DNK96YRNmIqGCopqi/C5fRhjqAnVEI6ECXgCuMVNqieVVG8qbnHTEG7A5/ZRWl/KtuptGAw+ty+6j7AJ43P5SPel23157XOKJwW3uPG6vURMhHAkTKo3Fa/L2+lzmmhNx9H0+XtcHhrCDURMBK/Li9flJWzCNIQbCIaDeN3e6GdbG6olGAlSVFtEbbCWkAkRioRwixuPy4Pb5aa8oRyvePG6vfjdfvxuPz63D6/Li8FgjKEx3AhAyISoC9aB2N+Lz+08XD5EBI/YbXpcHgTBJS7c4sZg2Fm7kx01O6LT6sP10SbjpoeI4BY3IoIxBre4d+/D2Y/P7SNiItSH6psdq9fljf7uXeJCRKgP1RM2YQS73ZAJETZhMJCXkkf/tP7d/nmJ6WVj8EyfPt0sWrSo0+vtXPsqx31wS/T9LdNv4TeLfsP3Z3yfy8Zexp0f3sm/1/w7On/OuDl8d/p3AagL1bG0aCkz+s/A5Vx9FIqEWFGygrs/vpvlJcuj63172repDdbyzKpnKGso63D6PC4PB2QdwJqyNXsEreOHHk9uIJfPiz5nZ81OyhrKGJY5jLNGnkVZfRlLi5byZcmXNIQb9thuv9R+lNWX0RhpbHP/gtA3tS+5gVwGpA1gYPpA0n3pZPuz2Va9Da/LS4onhZxADn1T+5LmTQPA7/ZTF6ojGAlSVl9GxETICeRQG6zliZVPsGTXErwub3T/XpeXgDtAXkoeeSl55Kfk0yelDwFPgNL6Uorrikn3pkczthUlK1hfsR6PeAiZUIfPZ0uDMwaT488hbMLUhmqpbqwmFAnhEhflDeX2h9ZLecRDiieFFG8KqZ5U+7rFI+AJEDERgpEgtcHa6Otmy7VY3+/2Ux2spi5Uh8/lI+AJRDPbhnADO2t3sqt2F4KwtXor5Q3l7KzZSVVjVbvfN9U1V46/km9P+3aX1hWRxcaY6fHmJU2NoGUfwZjcMUzsM5HHVzzOzIKZzYIAEI3SACmeFA4ZcMge8yfkT+DJU55kya4lHJh9IOnedNwu+3+FqydczbwN86gP19M3pS/PrnmWC0ZdwOjc0Wys3Mi4vHEsK15GVWMVy0uWs6FyA8uKl0WDwEnDTuJ7M76HiNAnpU90vxETsaUOZz9NguEg6yvWUx+uZ1D6INaWr2VlyUpWlq1kR80OZhbMZEa/GdSEathatRWXuKK1kKLaItZXrGdj5UYqGyrZVLmJj7Z/RF2obq9P+8nDTyY3kEumPxO3uKkOVlMbrKW8oZyi2iKWlyynqK6IulAdbnEzOGMwoUiI2lAtHvGQn5rP1ROuxhhDtj8bg6E6WM0BWQcwMX8ilY2VVDZUkhvIJTeQS3WwmoqGCiobK/G6vAxIG8CQzCHRAB6PMYaqYBXl9eWUNZTZWhVCTaiGmsYaakO1+N1+0r3ppHpTCUVCVDVW2UewimA4SN/UvuSn5kdLbh6Xh1RPKnXhOiImQl2ojtpgLfWhelK9qTSGG8nyZ1GQUYBb3AQjQfxuP5m+TEQkmmFXB6upaayxz8EaDLakGzZhAu4AdaE6u+1QbXQftSG7n8rGSnbW7owuAzZwNwVxr8sbndf0CEaCHf5sBSEnkANAujedIZm2Ztk/rT9Zvix8bh8ucWGMIWxszcDr8hKKhAhGgggSrSWFTZhgJEgwHMTvsSX8vEAeWf4sW1J3SttN66Z703GJi4ZwAw3hBhrDjfY50hj9z1BTU6pHPNHabmOkkcZwI8FwkMZIIxET2V3CNyGMMba25RQMsv3ZDMkYgsFO97l9tlbmlOTdLnezdZpK8MFwsFmaguEgYRMmw5dBMBKMHkfsuYiYiK0dunzRWl/ERPCIJ1pbGJIxpJO/wA5+lslSI9ixeh7Hf/T96PvfzPoNgvDdd7/bbLn3L3yfez65h4vHXMy4PuP2Or2d0fSFapnJ94SmH2/TZba5gVwawg3RJpDyhnJc4iIYCeJ1eXGJizRvGpm+TJs5N1ZyUO5B5KXkdWp/sQFY7XuhSCgaFOpD9aT70kn1pNIYaaQhtDvTTfGkkOXPItWb2tNJVh2kNQKg5Y1pRmSP4ICsA5pNu/2Q28nyZ/HzI36+LxMW1dTWuD9oajuNrY00NRkkpI3S2Z/qWR6XhwxfBhm+jGbTAwRAr0n4ykrK/xFcPOZiDsw+EBHhw4s+5NqJ1xJwBzh80OE9nTSllNrnkqcIFtMEVla/uxM3w5fB/035P26cfON+d+WFUkrtC0lZI1hRumKPaRoElFLJKokCwe4awdDMoT2YDqWU2r8kTdNQUxi4aMBMbuihzmCllNofJbRGICInisgqEVkrIrfFmS8icr8zf6mITE1YYpw+goPSC8jyZyVsN0op1dskLBCIiBv4E3ASMBa4SERajtR2EjDSeVwL/DlR6YlJWeJ3oZRSvUgiawQHA2uNMeuNMY3AU8AZLZY5A/iHsRYA2SIyIBGJWbB+biI2q5RSvV4iA8EgIPY2YIXOtM4ug4hcKyKLRGRRUVFRy9kdMnzIkZzi7csho8/q0vpKKfVVlcjO4nhtMC3Hs+jIMhhjHgIeAjvERFcSM3n8ht53vgAACJJJREFUxUwef3FXVlVKqa+0RNYICoHBMe8LgG1dWEYppVQCJTIQfAKMFJHhIuIDLgT+22KZ/wKXO1cPHQpUGGO2JzBNSimlWkhY05AxJiQiNwKvAm7gYWPMlyJyvTP/QWAucDKwFqgFrkhUepRSSsWX0D+UGWPmYjP72GkPxrw2wA2JTINSSqm2JdEQE0oppeLRQKCUUklOA4FSSiU5DQRKKZXket09i0WkCNjUxdX7AMXdmJzeQI85OegxJ4e9Oeahxpj8eDN6XSDYGyKyqLWbN39V6TEnBz3m5JCoY9amIaWUSnIaCJRSKsklWyB4qKcT0AP0mJODHnNySMgxJ1UfgVJKqT0lW41AKaVUCxoIlFIqySVNIBCRE0VklYisFZHbejo9XSUig0XkbRFZISJfisjNzvRcEXldRNY4zzkx6/zAOe5VInJCzPRpIvKFM+9+Edmvb+gsIm4R+UxEXnLef6WPWUSyReQ5EVnpfN6HJcExf9v5Xi8TkSdFJPBVO2YReVhEdonIsphp3XaMIuIXkaed6R+LyLB2E2WM+co/sMNgrwMOAHzA58DYnk5XF49lADDVeZ0BrAbGAvcAtznTbwN+7bwe6xyvHxjunAe3M28hcBj2TnHzgJN6+vjaOfbvAE8ALznvv9LHDDwGXO289gHZX+Vjxt6mdsP/b+/8Y62uyzj+eiM/NF0QAu4iLtERNhcLsk2rKXqpP8yhOZggDNvYyqisNtd2g5z+ZzOUP5zKoGkFpWQuzWWwITD9B0h2Bwy9kUF1lZTNrGCLcXef/nieY9/Ozj0XrkfPPd/v89q++36+z+fH93nOufc8nx/f7/MBzovrLcBXy2YzcC0wDzhYkLXMRmAV8FiklwBPDatTuz+UD+mDvwbYWrjuAXrarVeLbHsW+CLQB3SFrAvoa2Qrvj/ENVHmtYJ8KbC+3fY0sXMGsB24oeAISmsz8NH4UVSdvMw21/Ywn4yHyH8e+FIZbQYurXMELbOxVibSY/E3kdVMn6pMDdX+wGr0h6yjiSHfXGA3cJHF7m5xnhbFhrL94kjXy0cr64DvA4MFWZltvgw4Djwe02EbJZ1PiW02szeAHwN/BY7hOxZuo8Q2F2ilje/VMbMB4J/Ahc1uXhVH0Gh+sKOfm5V0AfBr4Ltm9q9mRRvIrIl81CHpJuBtM3vlTKs0kHWUzXhPbh7wqJnNBU7iUwZD0fE2x7z4zfgUyHTgfEnLm1VpIOsom8+Akdh41vZXxRH0A5cUrmcAb7ZJl/eNpHG4E9hsZs+E+C1JXZHfBbwd8qFs7490vXw08nlgoaSjwJPADZI2UW6b+4F+M9sd10/jjqHMNi8AjpjZcTM7DTwDfI5y21yjlTa+V0fSWGAi8E6zm1fFEewFZkmaKWk8voDyXJt1GhHxZMBPgFfN7MFC1nPAHZG+A187qMmXxJMEM4FZwJ4Yfv5b0tXR5opCnVGFmfWY2QwzuxT/7l40s+WU2+a/A3+TNDtE3cAhSmwzPiV0taSPhK7dwKuU2+YarbSx2NYi/P+l+Yio3YsmH+LizI34EzavA6vbrc/7sOML+DBvP9Abx434HOB24HCcJxfqrA67+yg8PQFcBRyMvIcZZkFpNBzAfP63WFxqm4FPA3+I7/o3wMcqYPN9wGuh78/xp2VKZTPwS3wN5DTee1/ZShuBc4FfAX/Cnyy6bDidMsREkiRJxanK1FCSJEkyBOkIkiRJKk46giRJkoqTjiBJkqTipCNIkiSpOOkIko5BkklaW7i+W9K9LWr7CUmLWtHWMPdZHJFEd9TJx0QEyYMRUXJvPDeOpB980Hol1SYdQdJJnAJulTSl3YoUkXTOWRRfCawys+vr5LfhYRXmmNmngK8A70ZeOoLkAyUdQdJJDOB7tn6vPqO+Ry/pRJznS9olaYukP0q6X9IySXui5315oZkFkl6KcjdF/XMkPRA99P2Svl5od4ekXwAHGuizNNo/KOlHIbsHfyHwMUkP1FXpAo6Z2SCAmfWb2T8k3Q+cJ6lX0uZoZ3no3ytpfc0RSTohaa2kfZK2S5oa8rskHQr9nxzB556UnXa/ZZdHHmd6ACfw8MxH8fgpdwP3Rt4TwKJi2TjPx3vWXfhbqm8A90Xed4B1hfq/xztHs/A3Ps8FvgasiTIT8Dd9Z0a7J4GZDfScjodLmIoHj3sRuCXydgJXNagzI+zqBdYCc+ttifQngd8C4+L6EWBFpA1YFul7gIcj/SYwIdKT2v095jH6jhwRJB2FeaTVnwF3nUW1vWZ2zMxO4a/jbwv5ATwufI0tZjZoZoeBPwNX4PHwV0jqxcN9X4g7CvCYL0ca3O+zwE7z4GkDwGZ8M5JmdvUDs/H484PAdkndDYp2A58B9oZO3XjIaqLeU5HehI8+wENUbI5IngPN9Eiqydh2K5AkI2AdsA94vCAbIKY6IwjX+ELeqUJ6sHA9yP//D9THW6mF+/22mW0tZkiaj48IGjGibRHDUb0AvCDpLeAWPO5Mfds/NbOeM2kyzl/GHdFC4IeSrgwHlSRArhEkHYiZvYNvY7iyID6K95TBY9qPG0HTi+PpncvxXnYfvtvTNyL0N5I+Id8gphm7geskTYn5+6XArmYVJM2TND3SY4A5wF8i+3Tt/rhjWCRpWpSdLOnjkTcGjzYJcDvwcrR1iZntwDf2mQRcMOwnkVSKHBEkncpa4FuF6w3As5L24D+WQ/XWm9GH/2BfBNxpZv+RtBGfPtoXI43jeE99SMzsmKQeYAfeg/+dmQ0XBnkasEHShLjeg0eUBF8g3y9pn5ktk7QG2BY/8qeBb+JO4yRwpaRX8F2pbsP3694kaWLo8pCZvUuSFMjoo0lSEiSdMLPs7SdnTU4NJUmSVJwcESRJklScHBEkSZJUnHQESZIkFScdQZIkScVJR5AkSVJx0hEkSZJUnP8CK67sUUNjHz0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 206 ms, sys: 13.9 ms, total: 220 ms\n",
      "Wall time: 222 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "size = 10000\n",
    "sunnyDays = np.zeros(size)\n",
    "cloudyDays = np.zeros(size)\n",
    "rainyDays = np.zeros(size)\n",
    "\n",
    "sunnyDays[0] = 1\n",
    "for i in range(size - 1):\n",
    "\n",
    "  if (sunnyDays[i] == 1):\n",
    "    random = np.random.uniform(0,1)\n",
    "    if (random <= .9):\n",
    "      sunnyDays[i + 1] = 1\n",
    "    elif (random <= .975):\n",
    "      cloudyDays[i + 1] = 1\n",
    "    elif (random <= 1):\n",
    "      rainyDays[i + 1] = 1\n",
    "    \n",
    "  if (cloudyDays[i] == 1):\n",
    "    random = np.random.uniform(0,1)\n",
    "    if (random <= .8):\n",
    "      cloudyDays[i + 1] = 1\n",
    "    elif (random <= .95):\n",
    "      sunnyDays[i + 1] = 1\n",
    "    elif (random <= 1):\n",
    "      rainyDays[i + 1] =1\n",
    "  \n",
    "  if (rainyDays[i] == 1):\n",
    "    random = np.random.uniform(0,1)\n",
    "    if (random <= .5):\n",
    "      rainyDays[i + 1] = 1\n",
    "    elif (random <= .75):\n",
    "      cloudyDays[i + 1] = 1\n",
    "    elif (random <= 1):\n",
    "      sunnyDays[i + 1] = 1\n",
    "\n",
    "\n",
    "print(\"Sunny Probability is:\", np.sum(sunnyDays)/size)\n",
    "print(\"Cloudy Probability is:\", np.sum(cloudyDays)/size)\n",
    "print(\"Rainy Probability is:\", np.sum(rainyDays)/size)\n",
    "plt.plot(np.linspace(1, size, size), np.cumsum(sunnyDays) / np.linspace(1, size, size), label = \"P(sunny)\")\n",
    "plt.plot(np.linspace(1, size, size), np.cumsum(cloudyDays) / np.linspace(1, size, size), label = \"P(cloudy)\")\n",
    "plt.plot(np.linspace(1, size, size), np.cumsum(rainyDays) / np.linspace(1, size, size), label = \"P(rainy)\")\n",
    "plt.xlabel(\"Number of Steps\")\n",
    "plt.ylabel(\"Probability\")\n",
    "plt.title(\"Probability Estimate Vs Number of Steps\")\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "iIm91fHWMYpI"
   },
   "source": [
    "## 2. The thermodynamics of a paramagnetic material"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 295
    },
    "executionInfo": {
     "elapsed": 1441,
     "status": "ok",
     "timestamp": 1604601300837,
     "user": {
      "displayName": "Gavin Debrun",
      "photoUrl": "",
      "userId": "16932255542580517286"
     },
     "user_tz": 360
    },
    "id": "y3aoR98UMYpL",
    "outputId": "93328033-685d-43a3-d66f-f3671aff95cb"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZAAAAEWCAYAAABIVsEJAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3deXxU9fX/8dchLElYEvYtbCLKIggY3MUFse5brYKthVp3canLr7a2iGitVau21dZaW5eqILb6Fbcq4FasCzsiyCIgCSA7YUuAJOf3x72J45iEIclkksz7+XjMI3M/d5lzk8mc+XzuufeauyMiIrK/GiQ6ABERqZuUQEREpFKUQEREpFKUQEREpFKUQEREpFKUQEREpFKUQORbzOxuM9toZl+H0+eZWY6Z7TCzQQmMq8I4zMzN7MAaiOMpM7s73q+TbMzsl2b2RBy2+5iZ/bq6tysBJZAkY2YrzSw//CAueTwSzusC3Az0dfcO4SoPAGPcvZm7z6nC61b1A7664hgZ/g4sqr2hma03szOrEGOtECY5N7Ozo9ofDttHJyi0kjhOMLPcyDZ3v8fdL6vidkeb2fSo7V7l7ndVZbtSPiWQ5HRW+EFc8hgTtncDNrn7+ohluwGf13yI31FdcbwMZALHR7WfCjjwn2p4jdpgCTCqZMLMGgI/AL5MWERS7yiBCABmdjIwBegU9kommNkOIAWYZ2Zfhst1MrN/m9kGM1thZtdHbCMlHIr40sy2m9ksM+tiZh+Ei8wLt31RGa/fwMx+ZWZfhT2BZ8wsw8yalBVHOU42s6VmtsXMHo3uZQC4ewEwCfhx1KwfA8+5e6GZvWhmX5tZnpl9YGb9yvmdfecbb2RPK4z9ATNbZWbrwuGUtHK2Veb+h/O6h9sdFW5ro5ndXsHvAeBV4BgzaxlOnwrMB76OeM2eZvaOmW0Kt/mcmWVGzB9sZnPCv+WLZvZCyfBdSS/CzG4O411rZj+JWLfMfTezpsCbfPM+2xG+p8aZ2bPhuo9E9ZALzWxcOO+2iPfXQjM7L2zvAzwGHBWuszVs/9aQo5ldbmbLzGyzmU02s05Rf7ur9vUekm8ogQgA7j4VOA1YE/ZKRrp7s3D2oe7e08waEHwwzQM6A8OAG83se+FyNwEjgdOBFsClwC53HxqxnWbu/kIZIYwOHycCBwDNgEfcfXd0HBXsxpnAEOBQ4ELge+Us9zRwQcmHefhBfRbwTDj/TaAX0A6YDTxXwWtW5HfAQcBA4ECC39nYcpYdTRn7H7XMscDBBL/3seGHZnkKgMnAiHD6x3yzfyUM+C3QCegDdAHGAZhZY4Le2lNAK2ACcF7U+h2AjHC/fgo8GpGwytx3d9/Jt99nzdx9TeRG3b1kqLJZuM9bgFfC2V8Cx4WveyfwrJl1dPdFwFXAR+G6mUQxs5PC/b0Q6Ah8BUyMWizW95AAuLseSfQAVgI7gK0Rj8vDeScAuVHLO3Bg+PwIYFXU/F8AT4bPFwPnlPO6pdspZ/404JqI6YOBvUDDGNd34NiI6UnAbRUsvxS4OHx+OTCvnOUyw21nhNNPAXeHz0cD08vaT4IP551Az4h5RwEr9nf/ge7hdrMi5n8KjChnW08BdxN8+H5E8GG7DkgDpgOjy1nvXGBO+HwosBqwiPnTI/b9BCC/5O8Ttq0HjtzXvpfzPhsHPBvV1jZ8v5a5n+Eyc0vec+X8PSL/Xn8H7ouY1yz8HXevzHtID6chkozO9aDHsb+6EQw9bI1oSwH+Gz7vQuXH2DsRfCMs8RXBh2d7gg+yWHwd8XwXwQdEeZ4h+Fb+PHAJQa8EM0sBfkNwvKAtUBwu3wbIizEOwnXTgVkRoyBG8PsqS0X7X2J/9g93n25mbYFfAa+5e37kiIyZtQP+SPCNvjnBiMSWiHhWe/hJGsqJeolN7l5YRkz7u+/fYWaNgH8Bz7v7xIj2HxP0dLuHTc0I/jax6ETQowTA3XeY2SaC3tHKsHm/fsfJTkNYsj9yCL5FZkY8mrv76RHzKxpiqsgaggRVoitQSPDNOR6eAYaZ2VEE35qfD9svBs4BTib45t49bC9rLHwnwQdlsIBZh4h5Gwm+ofeL+F1l+DfDcdHitf/PElTWRQ9fQTCc48AAd28B/Ihv9nMt0DnqGECXGF9zX/seyyXA/wRsJ0h+AJhZN+BvwBigtQfDVAsiYt7Xdr/1Ow6Px7Qm9i8oEkUJRPbHp8A2M/t5eEA0xcwOMbMh4fwngLvMrJcFBphZ63DeOoKx/fJMAH5mZj3MrBlwD/BC1DfcauPuXxEMyUwAprh7yTfP5sBuYBNBcrings3MA/qZ2UAzSyU8fhBuv5jgw+6h8Js+ZtY54nhRtHjt/x+B4cAHZcxrTjicaWadgVsj5n0EFAFjLChxPgc4PJYXjGHf1wGtS4oEopnZlQRVcheH2yrRlCBJbAiX+wlwSMT8dUBWePymLM8DPwn/Xk0IfsefuPvKWPZLvksJJDm9GlXl8nIsK7l7EcHB5oHACoJvmk8QfFMHeJBg3PhtYBvBmHNJ1dE44Gkz22pmF5ax+X8A/yT4oFtBcBD4ukrs2/54muAbaeS382cIho9WAwuBj8tb2d2XAOOBqQTHVKZHLfJzYBnwsZltC5c7uJzNxWX/3X2zu0+LGooqcScwmGBo7nXgpYj19gDnExwc30rQO3mNILnGotx9d/cvCBLm8vD90Clq3ZEEXzbWRLxHf+nuC4HfEyS3dUB/4MOI9d4hKPX+2sw2lvG7mAb8Gvg3QQ+rJ98UGUglWNnvKxGRbzOzT4DH3P3JRMcitYN6ICJSJjM73sw6hENYo4AB1J8TLaUaqApLRMpzMMGQZDOC6roL3H1tYkOS2kRDWCIiUikawhIRkUpJqiGsNm3aePfu3RMdhohInTJr1qyN7t42uj2pEkj37t2ZOXNmosMQEalTzOyrsto1hCUiIpWiBCIiIpWiBCIiIpWSVMdAyrJ3715yc3MpKChIdCgJlZqaSlZWFo0aNUp0KCJSRyR9AsnNzaV58+Z0796dZL35mLuzadMmcnNz6dGjR6LDEZE6IumHsAoKCmjdunXSJg8AM6N169ZJ3wsTkf2T9AkESOrkUUK/A5F6av4keOgQGJcZ/Jw/qdo2nfRDWCIicTF/EkwbD3m5kJEFw8bCgLLuZFC93J1tBYWszcunaO4kDv70dhoWh6MLeTnw6vXB82qIRQmkFkhJSaF///6l0yNGjOC2226rlm3PnTuXNWvWcPrpp+97YRGpHvMnBR/Ue/OD6Wr84N5esJe1eQXBY2s+a/IK+Dovn7V5BazZms/XeQXs3FMEwPTGv6Vhg6ih6b35QWJTAqkf0tLSmDt3bly2PXfuXGbOnKkEIlKTpo3/JnmUiOGDe+fuwjA55LN2a0Hp89IksbWA7bu/fZNKM2jbrAkdM9Po1a45Qw9qS6eMNDpkpNL55U1lv1BeblX3EFACqbXy8vI4/PDDmTx5MgcffDAjR47kpJNO4vLLL+fqq69mxowZ5Ofnc8EFF3DnnXcCMGPGDG644QZ27txJkyZNmDJlCmPHjiU/P5/p06fzi1/8gosuuijBeyaSBMr5gPa8XP63bGO5vYdtBd+9g3GbZk3omJFK99ZNObpnGzpkpNIxI5VOmWl0zEilXfNUGjcs53D2O1lB7ydaRlZV9q6UEkiEO1/9nIVrtlXrNvt2asEdZ/WrcJn8/HwGDhxYOl3yQf/II48wevRobrjhBrZs2cLll18OwG9+8xtatWpFUVERw4YNY/78+fTu3ZuLLrqIF154gSFDhrBt2zbS09MZP348M2fO5JFHHqnW/RKRsu3cXUhK046k7lzznXmri1vzwyc+KZ1u1bQxHTNSyWqZzpDureiYmVrae+iUkUb7jCY0aZhS+WCGjf32UBpAo7SgvRoogdQC5Q1hDR8+nBdffJFrr72WefPmlbZPmjSJxx9/nMLCQtauXcvChQsxMzp27MiQIUMAaNGiRY3FL5Ksioud5Rt3MmfVFubkbGXOqq0s/nobZ9q53NvoCdJtT+myexuksnrQrTx/yBGlSSK1URWSQyxKhsvidDBfCSTCvnoKNa24uJhFixaRlpbG5s2bycrKYsWKFTzwwAPMmDGDli1bMnr0aAoKCnB3leKKxFnerr3Mzd3KnFVbmL1qK3NXbSkddmrepCEDu2Yy/MQDGdRtCMXb+sP0e0o/uBsNG8sRNVCF9R0DLoxb9ZcSSC320EMP0adPH+655x4uvfRSPvroI7Zt20bTpk3JyMhg3bp1vPnmm5xwwgn07t2bNWvWMGPGDIYMGcL27dtJS0ujefPmbN++PdG7IlJzqql8tqjYWfz1dubkbGHOqiBpfLlhJxAcuD64fXPOGNCRQV1aMqhrJj3bNqNBg8gvcRfDkIuraadqJyWQWiD6GMipp57KpZdeyhNPPMGnn35K8+bNGTp0KHfffTd33nkngwYNol+/fhxwwAEcc8wxADRu3JgXXniB6667jvz8fNLS0pg6dSonnngi9957LwMHDtRBdKn/qlA+u2H7bubmBIlizqqtzMvdyq6wHLZV08YM6pLJeYM6M7hrS/pnZdA8VdeNS6p7omdnZ3v0DaUWLVpEnz59EhRR7aLfhdR5Dx1STtVRF/jZgtLJPYXFLFy7rTRZzMnZQs7mIOk0bGD07dSCQV0yGdQ16F10bZWe1EPEZjbL3bOj29UDEZH6o4Ly2Tfmry092P3Z6jz2FBYD0L5FEwZ3bcklR3ZjUNeWHNIpg7TGcT64XU8ogYhI/ZFR9nkPq4tbc+3zs2ncsAEDOmcw6qhupb2LjhlpCQi0flACEZF6Y3X2rbR551aa+O7Stj3WhC8H3MQrRxxDn44tyj/pTvabEoiI1HlfbtjBw1OX8tr8NlzY+ApuT32R5rvXYRlZNB42luMTUT6bBBKaQMzsVOAPQArwhLvfGzX/VuCH4WRDoA/Q1t03m9lKYDtQBBSWdYBHROq3nM27+MO0pbw0O5cmDVO4+vieXDF0OC3S7050aEkhYQnEzFKAR4HhQC4ww8wmu/vCkmXc/X7g/nD5s4CfufvmiM2c6O4bazBsEakF1ubl88g7y3hhRg4NGhg/OaYHV5/QkzbNmiQ6tKSSyMHAw4Fl7r7c3fcAE4FzKlh+JDChRiKrYWbGJZdcUjpdWFhI27ZtOfPMM0vb3nzzTbKzs+nTpw+9e/fmlltuAWDcuHGYGcuWLStd9qGHHsLMKClZ3rFjB1deeSU9e/akX79+DB06lE8++eZ6PCJ1xYbtu7nz1c85/v73mDQzh5GHd+WDW0/k12f2VfJIgEQOYXUGIsslcoEjylrQzNKBU4ExEc0OvG1mDvzV3R8vZ90rgCsAunbtWvWo43CTmKZNm7JgwYLSEwCnTJlC586dS+cvWLCAMWPG8Prrr9O7d28KCwt5/PFvdrd///5MnDiRX/3qVwD861//om/fvqXzL7vsMnr06MHSpUtp0KABy5cvZ9GiRVWKWaQmbdm5h79+sJyn/7eSPUXFfH9wZ647qRddWqUnOrSklsgeSFln5ZR3VuNZwIdRw1fHuPtg4DTgWjMbWtaK7v64u2e7e3bbtm2rFnHJWa55OUGoJWe5VsMtIk877TRef/11ACZMmMDIkSNL5913333cfvvt9O7dG4CGDRtyzTXXlM4/99xzeeWVVwBYvnw5GRkZlOzrl19+ySeffMLdd99NgwbBn/uAAw7gjDPOqHLMIvG2rWAvD01ZwnH3vctfP/iS7/Vrz9Sbjue+Cw5V8qgFEplAcoEuEdNZwHevfxwYQdTwlbuvCX+uB14mGBKLr4puElNFI0aMYOLEiRQUFDB//nyOOOKbztiCBQs47LDDyl23RYsWdOnShQULFjBhwoRvXa7k888/Z+DAgaSk6MQoqTt27i7k0XeXcdzv3uUP05ZyXK82vHXjUB4eMYgebZomOjwJJXIIawbQy8x6AKsJksR3rjxmZhnA8cCPItqaAg3cfXv4/BSg6p/i+1LeXbyq4e5eAwYMYOXKlUyYMKFSdw8sSUBvvfUW06ZN48knn6xyTCI1rWBvEc9+/BV/ee9LNu3cw7De7fjZ8IM4pHNGokOTMiQsgbh7oZmNAd4iKOP9h7t/bmZXhfMfCxc9D3jb3XdGrN4eeDm8Nk1D4Hl3/0/cgy7nLNfqurvX2WefzS233MJ7773Hpk3f3IqyX79+zJo1i0MPPbTcdc866yxuvfVWsrOzv3UvkH79+jFv3jyKi4tLh7BEaps9hcW8MGMVj7y7jHXbdnPsgW246ZSDGNy1ZaJDkwok9DwQd38DeCOq7bGo6aeAp6LalgPlf5rGS5zv7nXppZeSkZFB//79ee+990rbb731Vs4//3yOPfZYDjroIIqLi3n44Ye56aabSpdJS0vjd7/7HQcddNC3ttmzZ0+ys7O54447GD9+PGbG0qVLWbhwIeecU1HRm0j8FRYV89Ls1fxh2lJWb81nSPeWPHzRII7q2TrRoUkMdCb6/ojz3b2ysrK44YYbvvuyAwbw8MMPM3LkSHbt2oWZlXkQfMSIEWVu94knnuDmm2/mwAMPJD09ndatW3P//fdXS8wilVFU7Lw2fw0PT13Kio07OTQrg3vO78/QXm2S+qq3dY0u565LmJfS70Lizd156/OveXDKEpas20HvDs25+ZSDOblPOyWOWkyXcxeR+KrgHCl3593F6/n920v4fM02erZtyiMXD+L0QzpG3cVP6hIlEBGpunLuBOjA/9JP4oG3FzNn1Va6tkrnwQsP5ZyBnUlR4qjzlEAIvh0le/c5mYYyJQ7KOUdq4yu388OdD9MpI5Xfnt+fCw7LolGKqgHri6RPIKmpqWzatInWrVsnbRJxdzZt2kRqamqiQ5G6qpxzoVoXbeDOs/sx4vAuNGmok1nrm6RPIFlZWeTm5rJhw4ZEh5JQqampZGVVz/kskoTKO0eqRWdGHd29xsORmpH0CaRRo0b06NEj0WGI1G3DxrL3/66jUXHBN22N0mhw8h2Ji0niToORIlJl7zY5gVt3X8qmhu1wDDK6wFl/rLZzpKR2SvoeiIhUzbL127n++TlktTud1KvuwproYyVZqAciIpW2Zecefvr0TJo0SuGJUdk0VfJIKvpri0il7Cks5qpnZ7E2r4CJVxxJ58y0RIckNUw9EBHZb+7OHZMX8MmKzdz3/QG6am6SUgIRkf325IcrmfBpDtec0JNzB3Xe9wpSLymBiMh+eW/xeu5+fSGn9G3PLaccnOhwJIGUQEQkZsvWb+e65+dwcIcWPHTRQF0IMckpgYhITL6puGqgiisBVIUlIjHYW1TM1c/NYu3WAiao4kpCSiAiUiF3Z+wrn/Px8s08dNGhHNZNFVcSSOgQlpmdamaLzWyZmd1WxvwTzCzPzOaGj7Gxrisi1eOp/61kwqeruPqEnpw3SBfclG8krAdiZinAo8BwIBeYYWaT3X1h1KL/dfczK7muiFTB+0s2cNdrCxnetz23quJKoiSyB3I4sMzdl7v7HmAicE4NrCsiMVi2fgdjnp/NQe2b87AqrqQMiUwgnYHIGwjkhm3RjjKzeWb2ppn12891MbMrzGymmc1M9nt+iMRq6649XPb0DJo0VMWVlC+RCaSsrzPR91WdDXRz90OBPwH/tx/rBo3uj7t7trtnt23bttLBiiSLvUXFXPPcbNZsLeCvlxxGVsv0RIcktVQiE0gu0CViOgtYE7mAu29z9x3h8zeARmbWJpZ1RWT/uTvjJn/O/77cxG/P789h3VolOiSpxRKZQGYAvcysh5k1BkYAkyMXMLMOFt6o3MwOJ4h3Uyzrisj+e+ajr3juk1VcdXxPvn+YKq6kYgkb2HT3QjMbA7wFpAD/cPfPzeyqcP5jwAXA1WZWCOQDI9zdgTLXTciOiNQTHyzZwPjXFnJyn/b8v++p4kr2zYLP4+SQnZ3tM2fOTHQYIrXOsvU7OO/PH9I5M41/XX00zXTQXCKY2Sx3z45u17WwRJJcScVV45Sg4krJQ2Kld4pIEousuHr+8iNUcSX7RQlEJInd+WpQcfXADw4lu7sqrmT/aAhLJEk989FKnv14FVcOPYALVHEllaAEIpKE/rt0A3e+upCT+7Tj/53aO9HhSB2lBCKSZL7csINrn5tNr3bNeHjEIFJ0jSupJCUQkSQSVFzNpFFKA/72Y1VcSdXo3SOSJPYWFXPt87PJ3bKL5y8/ki6tVHElVaMEIpIkxr+6kA+XbeL+CwYwRBVXUg00hCWSBP750Ur++fFXXDH0AH6Q3WWfy4vEQglEpJ6bvnQj415dyLDe7fi5Kq6kGimBiNRjyzfs4JrnZtGzbVMeHjFQFVdSrZRAROqpvF17uezpmTRMacDfRw2heWqjRIck9YwSiEg9VFJxlbNlF3+95DBVXElcqApLpB6667WFTF+2kftUcSVxpB6ISD3zz4+/4pmPvuLy43pwoSquJI7UAxGp6+ZPgmnjIS+XgqYdmZV3Hif1Po/bTuuT6MiknospgZjZMcA4oFu4jgHu7gfELzQR2af5k+DV62FvPgCpO9dwb8Mn8P79SWkwJMHBSX0X6xDW34EHgWOBIUB2+LNKzOxUM1tsZsvM7LYy5v/QzOaHj/+Z2aER81aa2WdmNtfMdJ9aSU7TxpcmjxKp7Cbtg98kKCBJJrEOYeW5+5vV+cJmlgI8CgwHcoEZZjbZ3RdGLLYCON7dt5jZacDjwBER8090943VGZdInZKXu3/tItUo1gTyrpndD7wE7C5pdPfZVXjtw4Fl7r4cwMwmAucApQnE3f8XsfzHgO56IxIpIwvycspuF4mzWBNIybf+7Ig2B06qwmt3BiLf+bl8u3cR7adAZC/IgbfNzIG/uvvjZa1kZlcAVwB07dq1CuGK1D4+bCy7XxpD6jff66BRGgwbm7igJGnElEDc/cQ4vHZZ11TwMhc0O5EggRwb0XyMu68xs3bAFDP7wt0/+M4Gg8TyOEB2dnaZ2xepq97gON7a81N+2+JlmhZ8HfQ8ho2FARcmOjRJArFWYWUAdwBDw6b3gfHunleF184FIovUs4A1Zbz2AOAJ4DR331TS7u5rwp/rzexlgiGx7yQQkfqqYG8R97yxiObtTiP1+ntA17mSGhZrFdY/gO3AheFjG/BkFV97BtDLzHqYWWNgBDA5cgEz60pw3OUSd18S0d7UzJqXPAdOARZUMR6ROuXv01ewems+Y8/qq4skSkLEegykp7t/P2L6TjObW5UXdvdCMxsDvAWkAP9w98/N7Kpw/mPAWKA18GczAyh092ygPfBy2NYQeN7d/1OVeETqkvXbCnj03WWc0rc9R/dsk+hwJEnFmkDyzexYd58OpScW5u9jnX1y9zeAN6LaHot4fhlwWRnrLQcOjW4XSRb3v7WYvUXF3H6GzjaXxIk1gVwNPB0eCzFgMzA6XkGJSPk+y83jX7NzueK4A+jWummiw5EkFmsV1lzgUDNrEU5vi2tUIlImd+eu1xbSKr0x1550YKLDkSRXYQIxsx+5+7NmdlNUOwDu/mAcYxORKG989jWfrtzMPef1p4VuECUJtq8eSEn/uHkZ83ROhUgNKinb7d2hORcN0WXaJfEqTCDu/tfw6VR3/zByXnggXURqSEnZ7vOXH6GyXakVYj0P5E8xtolIHKzfVsCfVbYrtcy+joEcBRwNtI06DtKC4NwNEakBD7y9mD1FxfzydJXtSu2xr2MgjYFm4XKRx0G2ARfEKygR+caC1Xm8OCuXy487gO5tVLYrtce+joG8D7xvZk+5+1c1FJOIhNyd8a8GZbtjVLYrtUysJxLuCu8H0g9ILWl096pczl1E9uHNBSrbldor1oPozwFfAD2AO4GVBBdDFJE4Udmu1HaxJpDW7v53YK+7v+/ulwJHxjEukaT39+kryN2Sz9gzdbVdqZ1iHcLaG/5ca2ZnENy3Q/fMFImTkrLd4X3bc/SBKtuV2inWBHJ3eCHFmwnO/2gB/CxuUYkkuZKy3dtVtiu1WKwJZF5498E84EQAM+sQt6hEkpjKdqWuiPUYyAozm2Bm6RFtb5S7tIhUirsz/jWV7UrdEGsC+Qz4L/BfM+sZtumonkg1e3PB13y6YjM3nXKQynal1ot1CMvd/c9mNg941cx+jq7GK1KtIst2RwzpmuhwRPYp1gRiAO7+oZkNA14AesctKpEk9I8Pg7Ld5y/T1Xalboh1COv0kifuvhY4CTi1qi9uZqea2WIzW2Zmt5Ux38zsj+H8+WY2ONZ1ReqS9dsLePQdle1K3RLTHQmBkSV3IYzyQWVf2MxSgEeB4UAuMMPMJrv7wojFTgN6hY8jgL8AR8S4rkid8cBbutqu1D2JvCPh4cAyd18OYGYTgXOAyCRwDvCMuzvwsZllmllHoHsM64rUCZFluz1Utit1SCLvSNgZyImYziXoZexrmc4xrlsS5xXAFQBdu+rApNQuKtuVuiyRdyQsa0wsuldT3jKxrBs0uj/u7tnunt22bdv9DFEkvv6jsl2pwxJ5R8JcIPISo1kE19iKZZnGMawrUqsV7C3iNyVX283W1Xal7tlXDyT6joQlj+q4I+EMoJeZ9TCzxsAIYHLUMpOBH4fVWEcCeWEVWCzritRqJWW7vz6zLw1TYh0MEKk99uuOhGbW1N13VscLu3uhmY0B3iLozfzD3T83s6vC+Y8RXC7ldGAZsAv4SUXrVkdcIjUhsmz3GJXtSh0V64mEnczsTYLeSFczOxS40t2vqcqLu/sbRF1TK0wcJc8duDbWdUXqit+/tURlu1Lnxdpvfhj4HrAJwN3nAUPjFZRIfbZgdR6TZuUw+ujuKtuVOi3mgVd3z4lqKqrmWETqvZKy3ZbpjRlzUq9EhyNSJbEmkBwzOxpwM2tsZrcAi+IYl0i9VFK2e/MpB5GRprJdqdtiTSBXERyL6ExQWjuQco5NiEjZCvYWcc+bKtuV+iOmg+juvhH4YZxjEanXnvxwJTmb83nusiNUtiv1QkwJxMzaApcTXIOqdB13vzQ+YYnUL+u3F/DIO0s5uY/KdqX+iLWM9xWCOxJORQfPRfZbSdnu7WeobFfqj1gTSLq7/zyukYjUUyVlu5cd20Nlu1KvxDoQ+/FzcNMAABKASURBVJqZnb7vxUQkkrtzl8p2pZ6KNYHcQJBE8s1sm5ltN7Nt8QxMpD546/Ov+WTFZm4arrJdqX9ircIq64ZSIlKBkqvtHty+OSOGqGxX6p9Yq7AGl9GcB3zl7oXVG5JI/aCyXanvYj2I/mdgMPBZON0fmAe0NrOr3P3teAQnUlet317Ao+8uU9mu1Guxfi1aCQxy98Pc/TCCM9EXACcD98UpNpHabf4keOgQGJcZ/Jw/qXTWg28vYXdhkcp2pV6LtQfSO/J+G+6+0MwGuftys7LuLitSz82fBK9eD3vzg+m8nGAaWND6e7wwM4efHqOyXanfYk0gi83sL8DEcPoiYImZNQH2xiUykdps2vhvkkeJvfn4tPHcld6FlumNuW6Yynalfot1CGs0wV0BbwR+BiwP2/YCJ8YjMJFaLS+33HaV7UqyiLWMNx/4ffiItqNaIxKpCzKygmGrKOusjcp2JWnE1AMxs15m9i8zW2hmy0se8Q5OpNYaNhYapX2raW+DVO7Z/QN+dWYfle1KUoj1Xf4k8BegkGDI6hngn5V9UTNrZWZTzGxp+LNlGct0MbN3zWyRmX1uZjdEzBtnZqvNbG740GVWpGYNuBDO+iNkdAGMouZZ3F50ObsOPp/jerVNdHQiNSLWBJLm7tMAc/ev3H0ccFIVXvc2YJq79wKmhdPRCoGb3b0PcCRwrZn1jZj/kLsPDB9vVCEWkcoZcCH8bAGM28rt3SfwcuHR3H5G332vJ1JPxJpACsysAbDUzMaY2XlAuyq87jnA0+Hzp4Fzoxdw97XuPjt8vp3gFrqdq/CaInHx+Zo8XpiZw6ijuqtsV5JKrAnkRiAduB44DLgEGFWF123v7mshSBTsIxmZWXdgEPBJRPMYM5tvZv8oawgsYt0rzGymmc3csGFDFUIW+S53Z/yrC8lMa6SyXUk6MSUQd5/h7jvcPdfdf+Lu57v7xxWtY2ZTzWxBGY9z9idAM2sG/Bu40d1LrgD8F6AnwRnxaym7Oqwk9sfdPdvds9u21di0VK8XZuQEZbunHKyyXUk6FZbxmtnkiua7+9kVzDu5gu2uM7OO7r7WzDoC68tZrhFB8njO3V+K2Pa6iGX+BrxWUZwi8fDJ8k38+pUFHNerDSNVtitJaF/ngRwF5AATCIaPquu6JZMJhsDuDX++Er2ABddI+TuwyN0fjJrXsWQIDDiP4LpcIjVm1aZdXPXsLLq0SueRiwerbFeS0r7e9R2AXwKHAH8AhgMb3f19d3+/Cq97LzDczJaG27wXwMw6mVlJRdUxBMdaTiqjXPc+M/vMzOYTlBX/rAqxiOyX7QV7+enTMyh2+PuoIRq6kqRVYQ/E3YuA/wD/Ca97NRJ4z8zGu/ufKvui7r4JGFZG+xrg9PD5dMrp8bj7JZV9bZGqKCp2rp8whxUbd/LMpYer6kqS2j4vZRImjjMIkkd34I/ASxWtI1Jf3fvmIt5dvIG7zz2Eo3WfD0ly+zqI/jTB8NWbwJ3urmMNkrQmzczhb/9dwaijuvGjI7slOhyRhNtXD+QSYCdwEHB9xL0/DHB3bxHH2ERqjU9XbOb2lz/juF5t+PWZOttcBPZ9DESlJZL0cjZHVFyNVMWVSAn9J4hUoKTiqqjYg4qrdFVciZSI9Y6EIkmnqNi5YeJcvtywk3+q4krkO9QDESnH7/7zBe98sZ5xZ/dTxZVIGZRARMrw4swcHv9gOT8+qhuXqOJKpExKICJRZqzczC9f/oxjD2zDWFVciZRLCUQkQs7mXVz5z1l0aZnOo7rGlUiF9N8hEtqxu5DLnp5JYVExT4zKVsWVyD6oCkuEsOJqwhyWbdjB0z85nAPaNkt0SCK1nnogIsB9//mCaV+sZ9xZfTm2lyquRGKhBCJJ78WZOfz1g+VccmQ3Ljmqe6LDEakzlEAkqZVUXB1zYGvGnqWKK5H9oQQiSauk4iqrZTp/vvgwGqniSmS/6D9GkpIqrkSqTlVYknSKip0bJwYVV0/9ZAg9VXElUikJ6YGYWSszm2JmS8OfLctZbmV47/O5ZjZzf9cXKct9b33B1EXrueOsvhzXq22iwxGpsxI1hHUbMM3dewHTwunynOjuA909u5Lri5T616xc/vr+cn50ZFd+rIorkSpJVAI5B3g6fP40cG4Nry9JaNZXm/nlS0HF1R1n9Ut0OCJ1XqISSHt3XwsQ/mxXznIOvG1ms8zsikqsj5ldYWYzzWzmhg0bqil8qWtyt+ziimdm0SkzlUcvHqyKK5FqELeD6GY2FehQxqzb92Mzx7j7GjNrB0wxsy/c/YP9icPdHwceB8jOzvb9WVfqh5KKqz1FxTwxagiZ6Y0THZJIvRC3BOLuJ5c3z8zWmVlHd19rZh2B9eVsY034c72ZvQwcDnwAxLS+SHGxc+PEuSxdv4MnRw/hwHaquBKpLonqx08GRoXPRwGvRC9gZk3NrHnJc+AUYEGs64sA3PfWYqYuWsevz+jD0INUcSVSnRKVQO4FhpvZUmB4OI2ZdTKzN8Jl2gPTzWwe8Cnwurv/p6L1RSL9e1Yuj73/JT88oiujju6e6HBE6p2EnEjo7puAYWW0rwFOD58vBw7dn/VFSsz6ajO/eOkzju7ZmnFn98PMEh2SSL2jUhSpd3K3BNe46pSZyp9/qIorkXjRpUykXtkZVlztLixm4hWquBKJJ301k3qjuNi58YW5LFm3nUcvHqyKK5E4UwKReuP+txczZeE6fn1mX1VcidQADWFJ3TR/EkwbD3m5kJHFpz3H8Jf/dWHk4V0ZrYorkRqhHojUPfMnwavXQ14O4JCXwyGzfs3NHeYy/hxVXInUFCUQqXumjYe9+d9qSrc9XFP8vCquRGqQ/tuk7snLLbM5ZdvqGg5EJLkpgUjdk5G1f+0iEhdKIFKn5O8pYlrnK8kn6vyORmkwbGxighJJUkogUifsLiziqQ9XMPT+d/np7AN4qtVN7GnaGTDI6AJn/REGXJjoMEWSisp4pVbbW1TMv2bl8qdpS1mTV8ARPVrx5x8OZkj3M4BfJDo8kaSmBCK1UlGx88rc1Tw8dSmrNu9iUNdM7v/BoRzds7XKdEVqCSUQqVWKi503FqzloSlL+HLDTvp1asE/Rmdz4sHtlDhEahklEKkV3J2pi9bz+7cX88XX2zmofTMe+9FgTunbgQYNlDhEaiMlEEkod+eDpRt58O3FzMvNo3vrdP4wYiBnDuhEihKHSK2mBCIJ8/HyTfz+7cXMWLmFzplp3HfBAM4f1JmGOptcpE5QApEaN3vVFh58ewnTl22kfYsm3HXuIVyU3YXGDZU4ROoSJRCpMQtW5/HglCW888V6WjdtzK/O6MOPjuxGaqOURIcmIpWQkARiZq2AF4DuwErgQnffErXMweEyJQ4Axrr7w2Y2Drgc2BDO+6W7vxHnsKWSlqzbzkNTlvDmgq/JSGvErd87mNFHd6dpE31/EanLEvUffBswzd3vNbPbwumfRy7g7ouBgQBmlgKsBl6OWOQhd3+ghuKVSlixcSd/mLqEV+atoWnjhtwwrBc/Pa4HLVIbJTo0EakGiUog5wAnhM+fBt4jKoFEGQZ86e5fxTcsqQ45m3fxp3eW8u/Zq2mc0oArh/bkyqEH0LKp7k8uUp8kKoG0d/e1AO6+1sza7WP5EcCEqLYxZvZjYCZwc/QQWAkzuwK4AqBr165Vi1oq9HVeAY++u4yJM1ZhZow6qjtXn9CTts2bJDo0EYkDc/f4bNhsKtChjFm3A0+7e2bEslvcvWU522kMrAH6ufu6sK09sBFw4C6go7tfuq+YsrOzfebMmfu9LxIh6layDBvLxgPO4S/vfck/P/6K4mLnoiFdGHPSgXTMSEt0tCJSDcxslrtnR7fHrQfi7idXEMw6M+sY9j46Ausr2NRpwOyS5BFuu/S5mf0NeK06YpZ9KLmVbMndAPNy2Pt/1/G7wnn8e+/RnD84ixuG9aJLq/TExikiNSJRQ1iTgVHAveHPVypYdiRRw1clySecPA9YEI8gJUoZt5JtVFzALxq/yFXX/YKebZslKDARSYREnbl1LzDczJYCw8NpzKyTmZWW45pZejj/paj17zOzz8xsPnAi8LOaCTv5uDtrtubz2vw1eDm3km1VuF7JQyQJJaQH4u6bCCqrotvXAKdHTO8CWpex3CVxDTCJ5e8p4rPVecxZtYU5q7YyJ2cL67btBmBQk9Z0to3fXUm3khVJSjqTK4m5O19t2sXsiGSxaO12ioqDwopurdM56oDWDOrakkFdM2m/8R54/YZvD2PpVrIiSUsJpK4oo/ppf2/huq1gL/Nzwt5FzlbmrNrCll17AWjaOIWBXTO5+vieDOqaycAumbRuFlV+m3URNLAqxyEi9YMSSF1QRvUTr14fPC/nw7u42Fm6fse3hqKWrt9BSdV2r3bNGN63fWnvole75rFdPn3AhUoYIgIogdQNZVQ/sTc/aA8/zDfv3MPcnDBZrNrKvJytbN9dCEBmeiMGdcnkzAGdGNQ1kwFZmWSk6XIiIlI1SiB1QTnVT56Xy00vzGXOqi2s3LQLgJQGRu8OzTlnUCcGdWnJ4G4t6d46XbeDFZFqpwRSS+0uLGJd3m7W5uXTP60D6flrv7PM6uLWTF+2kcFdMxlxeFcGdcmkf1YG6Y31ZxWR+NMnzb5Uw8HraHuLilm3rYC1eQWs2ZrP13nfPF8bPt+4Y3fp8mc3OI97Gz1Buu0pbStMSSX91Dv59Ihh6l2ISEIogVSkEgevC4uKWb896DmszStg7dYC1uTls3ZrAWu3FbB2az4bduwm+hJkzZo0pGNGKh0z0+jXqQUdMlLplJFGx8xUOmYMpUHuofD+3aWJrOGwsbTSwWwRSSAlkIqUc/B6z1vjmGrHfbv3ECaJ9dsLKI5KDumNU+iYkUqnzDQOOqgtHTPT6JSRGiSJzDQ6ZqTSfF/3yGg3AgaPqN79ExGpAiWQipRz8LrhjjVc89xsAFIbNaBTRhodMlI55sA2dMpMjeo9pNEitaGGmUSk3lECqUhGVjBsFWV30468fuWxdMpIIzO9kZKDiCSlRF1MsW4YNja4VEekRmmknXon/Tpl0LJpYyUPEUlaSiAVGXAhnPVHyOgCWPDzrD/qTGwRETSEtW+6dIeISJnUAxERkUpRAhERkUpRAhERkUpRAhERkUpRAhERkUpRAhERkUoxj76qXz1mZhuArxIdRyW0ATYmOogalGz7C9rnZFFX97mbu7eNbkyqBFJXmdlMd89OdBw1Jdn2F7TPyaK+7bOGsEREpFKUQEREpFKUQOqGxxMdQA1Ltv0F7XOyqFf7rGMgIiJSKeqBiIhIpSiBiIhIpSiB1CFmdouZuZm1SXQs8WZm95vZF2Y238xeNrPMRMcUL2Z2qpktNrNlZnZbouOJNzPrYmbvmtkiM/vczG5IdEw1wcxSzGyOmb2W6FiqixJIHWFmXYDhwKpEx1JDpgCHuPsAYAnwiwTHExdmlgI8CpwG9AVGmlnfxEYVd4XAze7eBzgSuDYJ9hngBmBRooOoTkogdcdDwP8DkqLqwd3fdvfCcPJjICuR8cTR4cAyd1/u7nuAicA5CY4prtx9rbvPDp9vJ/hQ7ZzYqOLLzLKAM4AnEh1LdVICqQPM7GxgtbvPS3QsCXIp8Gaig4iTzkBOxHQu9fzDNJKZdQcGAZ8kNpK4e5jgC2BxogOpTrqlbS1hZlOBDmXMuh34JXBKzUYUfxXts7u/Ei5zO8GQx3M1GVsNsjLakqKXaWbNgH8DN7r7tkTHEy9mdiaw3t1nmdkJiY6nOimB1BLufnJZ7WbWH+gBzDMzCIZyZpvZ4e7+dQ2GWO3K2+cSZjYKOBMY5vX3hKVcoEvEdBawJkGx1Bgza0SQPJ5z95cSHU+cHQOcbWanA6lACzN71t1/lOC4qkwnEtYxZrYSyHb3unhFz5iZ2anAg8Dx7r4h0fHEi5k1JCgSGAasBmYAF7v75wkNLI4s+Cb0NLDZ3W9MdDw1KeyB3OLuZyY6luqgYyBSWz0CNAemmNlcM3ss0QHFQ1goMAZ4i+Bg8qT6nDxCxwCXACeFf9u54bdzqWPUAxERkUpRD0RERCpFCURERCpFCURERCpFCURERCpFCURERCpFCUQkgcysu5ktSHQcIpWhBCIiIpWiBCKSeClm9rfw3hhvm1laogMSiYUSiEji9QIedfd+wFbg+wmORyQmSiAiibfC3eeGz2cB3RMYi0jMlEBEEm93xPMidJVsqSOUQEREpFKUQEREpFJ0NV4REakU9UBERKRSlEBERKRSlEBERKRSlEBERKRSlEBERKRSlEBERKRSlEBERKRS/j/moH4S7lzGrQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.special as sps \n",
    "\n",
    "def ProbabilityRatio(h, new_theta, old_theta):\n",
    "  return np.exp(h * (np.cos(new_theta) - np.cos(old_theta) ) )\n",
    "\n",
    "N = 2000\n",
    "h_values = np.linspace(-5,5,10)\n",
    "magnetizations = []\n",
    "\n",
    "\n",
    "for h in h_values:\n",
    "\n",
    "  thetas = np.empty(N)\n",
    "  theta_initial = np.random.uniform(0, 2*np.pi)\n",
    "  thetas[0] = theta_initial\n",
    "\n",
    "  for t in range(N - 1):\n",
    "\n",
    "    r = np.random.uniform(0,1)\n",
    "    new_theta = thetas[t] + 0.75*np.random.randn()\n",
    "\n",
    "    if (r < ProbabilityRatio(h, new_theta, thetas[t]) ):\n",
    "      thetas[t + 1] = new_theta\n",
    "    else:\n",
    "      thetas[t + 1] = thetas[t]\n",
    "\n",
    "  magnetizations.append( np.sum( np.cos(thetas) ) / N)\n",
    "\n",
    "\n",
    "exact = sps.iv(1, h_values)/sps.iv(0,h_values)\n",
    "plt.plot(h_values, exact, label = \"Exact\")\n",
    "plt.plot(h_values, magnetizations, 'o', label = \"MCMC\")\n",
    "plt.title(\"Effect of h Value on Magnetization\")\n",
    "plt.xlabel(\"h\")\n",
    "plt.ylabel(\"Magnetization\")\n",
    "plt.legend()\n",
    "plt.show()  \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "65JKv19LMYpN"
   },
   "source": [
    "## 3: Thermodynamics of rubber bands"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "72__XSmAMYpP"
   },
   "source": [
    "### Rubber bands at infinite temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/"
    },
    "executionInfo": {
     "elapsed": 1997,
     "status": "ok",
     "timestamp": 1604601301415,
     "user": {
      "displayName": "Gavin Debrun",
      "photoUrl": "",
      "userId": "16932255542580517286"
     },
     "user_tz": 360
    },
    "id": "c9p4xXeDMYpQ",
    "outputId": "23516ad5-8e0a-4077-e4ea-ed7d09d0666a"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Average length of a polymer configuration is: 49.93230000000056\n"
     ]
    }
   ],
   "source": [
    "average_length = 0\n",
    "sample_size = 10000\n",
    "N = 100\n",
    "\n",
    "for sample in range(sample_size):\n",
    "  configuration = np.random.randint(2, size=N)\n",
    "  average_length += float(np.sum(configuration)) / sample_size\n",
    "\n",
    "print(\"Average length of a polymer configuration is:\", average_length)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "irzNcD8oMYpS"
   },
   "source": [
    "### Rubber bands at finite temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "colab": {
     "base_uri": "https://localhost:8080/",
     "height": 295
    },
    "executionInfo": {
     "elapsed": 5929,
     "status": "ok",
     "timestamp": 1604601320065,
     "user": {
      "displayName": "Gavin Debrun",
      "photoUrl": "",
      "userId": "16932255542580517286"
     },
     "user_tz": 360
    },
    "id": "7hcYOaYPMYpT",
    "outputId": "ca6fa26b-6dd1-4b18-87dc-4611f75857c1"
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYIAAAEWCAYAAABrDZDcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd3xV9f3H8dc7k4SEsDcCIkOWKy7U1i3uWheOOmq1WrW21V+1arHLtlpbra224ta6lSpa66yjiivIVkAEZEMYhrACST6/P84JvYSMA+TmJLmf5+NxHjn3rPs+997czz3re2RmOOecS11pcQdwzjkXLy8EzjmX4rwQOOdcivNC4JxzKc4LgXPOpTgvBM45l+K8ELgdImkXSWslpcedBUDSvyWdH3eOpqipvVf1qZ5X0tuSvhd3rvpIOkfSa3Hn2BFeCKqRNE/SJkkdqw2fJMkk9YknWfJIGixpnKQSSaWS3pI0oq55zGy+meWZWUUDPP92/aNL+oWkf1TLc6yZPbyzWWp4rofCz8PahO7MnVxmK0lfSzq8hnG3S3p2Z5ZfXUO+V9XV8vpM3pllJjNvomrZV0l6XdKgiPP2Cb8PMqqGmdljZnZ08hInjxeCms0Fzqp6IGkYkBNfnIaR+KFNGNYPeB+YCvQFugP/BF6TdGDU5bRwt4ZfTFXdU9szc/XXy8w2Ak8B51WbLp3gc7ddBa0JvB/VX589Ys6zPW41szygB7AIuD/mPLHwQlCzR9n6n/R84JHECSQVSHpEUrGkryTdKCktHHeBpPck3SZptaS5ko5NmLd7+At8laTZki5OGPcLSc9I+kf463yqpAGSfiZpuaQFko6uluN+SUskLZL0m4RN6gskvR/+ylwF/KKGdf0F8IGZ3WBmq8ys1MzuDF+DW8LlVP36uUjSfOA/1X8RRchR4+sh6WbgEOCv4S+zv4bD/xyu6xpJEyQdEg4fCVwPnJn46zNxq0JSWvh+fBW+Zo9IKqi2LudLmi9phaQbIn8ytv4MXBy+f6vC97N7wjiTdLmkL4Avapj9YeBUSbkJw44h+J/8t6QLJX0efgbmSPp+wrIPlbRQ0rWSlgIPSpom6cSEaTLDdduzhvfqbUm/Dj8bpZJeU8IWsKTzwtdupaSfK9hKPnIHXp+q571E0uLws3F1wvj9JBWF7/EySX+qNl9NP1yS8t6a2QbgaWDPhOc6XtLEMN8CSb9ImOXd8O/X4efwwKrPecL8IyR9omBL+xMlbGWH084JX/+5ks6J9KImi5l5l9AB84AjgZnA7kA6sADoDRjQJ5zuEeAFIB/oA8wCLgrHXQBsBi4O578MWAwoHP8OcDfQiuCDVwwcEY77BbCR4EshI3yeucANQGa4zLkJeZ8H7gFaA52Bj4HvJ+QoB64Ml5VTw/ouBS6sYfhhQAWQG66fhVlaE2wdVQ3LiJijrtfjbeB71Z7/XKBDmPvqMGerhNfoH9Wm37IM4LvAbGBXIA8YCzwajqvKfW+4HnsAZcDutXweHgJ+U8Pww4EVwN5ANvAX4N2E8Qa8DrSv6XUPp5kFnJvw+AngjrD/eKAfIOCbwHpg73DcoeH7ekv43DnAT4GnEpZ1MjC12jpnJLxWXwIDwnnfBn4fjhsMrAUOBrKA28L37sjteX2qPe8T4ediGMFn/chw/AfAd8L+POCAOvIm9b0N8z0KTE4Yf2iYOQ0YDiwDvlVTxoTP+Xthf3tgNfAdgs/wWeHjDuFzrQEGhtN2A4bE+r0X55M3xY7/FYIbgd8BIwn+oTPCN74PwZdZGTA4Yb7vA28nfCBmJ4zLDeftCvQi+ILNTxj/O+ChsP8XwOsJ404M/zHTw8f54bLaAl3CHDkJ058FvJWQY34961sOjKxh+KDweXokfOh3TRi/5R8hYo4aX4/w8dtUKwQ15FkN7JHwGtVVCN4EfpAwbiDBl1lGQu6eCeM/BkbV8rwPERTmr8NuRTj8foLdClXT5YXP0Sd8bMDh9azTjcBrYX8bgi/7vWqZ9nngqrD/UGATYWEMh3UHSoE24eNngZ9Wf68SXqsbE+b9AfBK2D8aeKLae7WJugtB4uvzNfBwtecdlDD9rcD9Yf+7wC+BjtWWWVPeZL+3lQQ/uIbX8X7dAdxeU8aEz3lVIfgO8HG1+T8Ip2kdPuep1PIjobE73zVUu0eBswneuEeqjetI8Gvpq4RhXxF8aVZZWtVjZuvD3jyCf9hVZlZax7zLEvo3EHz5VCQ8rlpWb4KthCUKDj5+TfCrvHPC/AtqX0Ug+FXbrYbh3Qj+OVZHWFaUHLW9HjWSdHW4a6QkXF4BweseRXe2fW+qCtY2eQi+gGvNAtxmZm3DrirDVs9hZmuBlWz9Ptb32j8CHCapB3AaQbGcCCDpWEkfhrudvgaOY+v1L7bgWEPV8y8mONZzqqS2wLHAY3U8d23r3z0xd/heraxnPRJfn7Zmdn618Ymvw1fhcwBcRLBVMiPcdXJCPc9Tla/B31uCL/YNBIUFAEn7KzhxolhSCXApO/4ZrMraw8zWAWeGy1si6V+KeJA6WbwQ1MLMviL4hXAcweZnohUEv0J6JwzbheBgU30WA+0l5e/AvNUtIPgl3jHhn7CNmQ1JmMbqWcYbwOk1DD+D4NjB+oRhtS0rSo66bLVcBccDrg0ztAv/UUsIdpPUlaPKYrZ9b8rZusDurK2eQ1Jrgs3+xPexzpxmNh/4L3AOwS/IR8JlZQPPEeyW6RKu/8v8b/1rW/bDBLvUTid473bkM7UE6Fn1QFIOwXrtjF4J/bsQvHaY2RdmdhbBD4ZbgGfD17EuSXlvw/fiKuDP4ToDPA6MA3qZWQHwd3b8M1iVdVH4fK+a2VEEP7hmEOzOio0XgrpdRLB5vy5xYPjr/GngZkn5knoDPwH+UcMytmJmC4DxwO8UnEY4PHyeun691basJcBrwB8ltQkPpPWT9M3tWMwvgRGSbpbUPlyfKwkOll/bSDmWEezzrZJP8M9dDGRIGk2w6yRx+j4KD87X4Angx5L6SsoDfkuw/7w8Yp4oHgcuDA/GZofP8ZGZzdvO5TwMXAEcxP8+A1kE+/6LgXIFB9ajnJb4PMExi6vYdis2qmeBE8MDnVkEnw/VM099fi4pV9IQ4EKCM6aQdK6kTmZWSbCrBILdpnVJ2ntrZq8TfIFfEg7KJ9h63yhpP4I9BFWKCbaYd6VmLwMDJJ0tKUPBKceDgZckdZF0Ulj0ygh2/Sb1VNn6eCGog5l9aWZFtYy+ElgHzAHeI/hieCDios8i2BRdTHCq5k3hh3BHnEfwxfEZwW6cZ6l5V0+NzOwLggODexAcH1lCsO/yGDN7v5Fy/Bk4TcEZRXcCrwL/JjiY+hXBftzE3QvPhH9XSvq0huU9QLBr712CrbqNBO9XgzGzN4GfE/xyX0JwYHfUDizqWaAd8GZYUAl3G/6Q4MfGaoIvoHERMm0I8/Rl263YSMxsOsFr9STBepUCywm+sGrzU219HcGKauPfITjA+ybBrpiqi65GAtMlrSX4DIxK3N1Vi2S/t38gWJ9sgmMnv5JUSnDs5OmqicIt5ZuB98PdoQckLsTMVgInEJzosJLgYP4JZraC4Hv3aoL//1UEJwP8oAHXYbtVnbXhnGsBwq2nAWZ2bgMtL4/g13p/M5u7nfP2IfiyzmzgrTHXwHyLwLkWQlJ7gt2MY3ZyOSeGu3JaExynmEqwtehaKC8EzrUACi5KXAD828zerW/6epxMsNtiMdCfYJeN7zpowXzXkHPOpTjfInDOuRQXd2NV261jx47Wp0+fuGM451yzMmHChBVm1qmmcc2uEPTp04eiotrO6HTOOVcTSdWvdN7Cdw0551yK80LgnHMpzguBc86luKQVAkkPhDeOmFbLeEm6U8GNPaZI2jtZWZxzztUumVsEDxG0JVKbYwkuVulP0MjT35KYxTnnXC2SVgjCqxtX1THJycAjFvgQaCspcmNpzjnnGkacxwh6sHWLkgvZ+qYeWyi452mRpKLi4uJGCeecc6kizusIamrjvMb2LsxsDGFDWoWFhTvUJsasZaW8NHkxudkZ5Galk5tV9Xfr/tbZGeRkpZObmU5Guh9Ld861fHEWgoVsfeeinoR3LkqGWctKufM/s7drnqyMNFqHhaJ1djpdC3Lo1S6HXdrn0qt9Lr3a5bJL+1wKcjOTlNo555IvzkIwDrhC0pPA/kBJ1Y05kuGE4d05bmg3NpZXsK6sgg2bKli3qZz1mypYX/1vWUW14RWUbtzMkpKNTFn4NV+v37zVsvNbZQTFoV0uvdoHhaJn+6BI9GibQ6vM9GStlnPO7bSkFQJJTwCHAh0lLQRuIrjBOWb2d4JbuR1HcOei9QS3sEuqtDSFu4F2brXXbNzMglXrWbBqQ/B39Xrmr1rPF8tLeWvmcsrKK7eavkubbAZ0yWdojwKG9ShgaPcCerXPQdrZOwA659zOa3bNUBcWFlpTbmuostIoXlv2vwKxcgPzV61nxtI1zFpWyuaK4PUuyMlkaI82W4rDsB4F7NI+14uDcy4pJE0ws8KaxjW7RueaurQ00aVNK7q0aUVhn/ZbjSsrr2Dm0lKmLiph2qISpi4q4YH35m4pDm1aZfxvqyH827uDFwfnXHJ5IWhE2RnpDO/ZluE9224ZVlZewayla5kaFoZpi0p48P15bKoIdi/lt8pgSPc2W4rD0B4F9O3QmrQ0Lw7OuYbhhSBm2RnpDOtZwLCeBVuGbSqvZNay0i1bDdMWlfDwB1+xKTz20DorncHdg91KQ7sHxaFfp9Z+uqtzbod4IWiCsjLStvz6HxUO21xRyezlwZbD9EUlTFu8hic/XsCGzfMAaJWZxu7d2jC0e7BLaUiPNvTvnE9WhhcH51zd/GBxM1ZRacwpXhtuNaxh2uISPlu8hrVl5QBkpacxtEcbjh3ajeOHd6N725yYEzvn4lLXwWIvBC1MZaUxb+U6pi1ew7RFJYz/cgXTFq0BYN8+7Thxj+4cO7QbnfKzY07qnGtMXghS3NwV63hp8mJemrKEmctKSRMc2K8DJw7vzsihXWmbmxV3ROdcknkhcFvMXFrKS1MW8+LkxcxbuZ6MNPGNAZ04YXg3jhrchfxW3lyGcy2RFwK3DTNj+uI1vBhuKSz6egNZGWkcPrAzJ+7RncMHdSYny5vGcK6l8ELg6lRZaUxcsJoXJy/hX1OXUFxaRm5WOkcN7sIPDt2NgV3z447onNtJXghcZBWVxkdzV/LSlCW8OHkx68rKOWu/XfjxUQPomOcHmJ1rrrwQuB2yet0m/vzmFzz64VfkZKZz+WG7ceFBfbw1VeeaoboKgV9t5GrVrnUWvzhpCK/+6BscsGt7bnllBkf+6R1emrKY5vYDwjlXOy8Erl67dc7jvvP35bHv7U9edgZXPD6RU/82nonzV8cdzTnXALwQuMgO2q0j//rhIdxy6jAWrN7AKXeP56onJ7Lo6w1xR3PO7QQ/RuB2yNqycu5550vGvDsHgO8d0pfLDt2NvGxvvsq5psiPEbgGl5edwdVHD+Staw7l2KFdueutLzn0D2/z5MfzqahsXj8unEt1XgjcTuneNoc7Ru3F85cfRJ8OuVw3dirH3/lf3vtiRdzRnHMReSFwDWLPXm155tIDufucvVm3qZxz7/+I6/85lfKKyvpnds7FynfougYjieOGdeOI3Tvzp9dncc87c1jy9Qb+evbetPZjB841Wb5F4BpcdkY6Pzt2d24+ZSjvzCpm1JgPWV66Me5YzrlaeCFwSXPO/r257/xCZi9fyyl3jWf28tK4IznnalBvIZDUSdL1ksZIeqCqa4xwrvk7fFAXnvr+AZSVV3Dq3z7g47mr4o7knKsmyhbBC0AB8Abwr4TOuUiG92zLP39wEB3ysjj3vo94cfLiuCM55xJEOYKXa2bXJj2Ja9F6tc9l7GUjuPiRIq58YiJLSjZw8SG7IinuaM6lvChbBC9JOi7pSVyL1zY3i0cv2p/jh3fjty/P4KZx0/3iM+eagChbBFcB10vaBGwOh5mZtUleLNdStcpM5y+j9qJH2xzGvDuHJSUbuXPUXn43NOdiVO8WgZnlm1mambUK+/O9CLidkZYmrj9ud3550hDe+HwZo+79kBVry+KO5VzKinT6qKSTJN0WdickO5RLDeeP6MPfz92HGUvW8O27xzOneG3ckZxLSVFOH/09we6hz8LuqnCYczvtmCFdeeKSA1hbVs6pfxvPhK/89FLnGluULYLjgKPM7AEzewAYGQ5zrkHsvUs7xl42goKcTM6+9yNembYk7kjOpZSoVxa3TegviLpwSSMlzZQ0W9J1NYwvkPSipMmSpku6MOqyXcvSp2NrnrtsBIO7t+Gyxz7l0Q/mxR3JuZQRpRD8Dpgo6SFJDwMTgN/WN5OkdOAu4FhgMHCWpMHVJrsc+MzM9gAOBf4oKWs78rsWpENeNk9cfABHDOrCz1+YzuMfzY87knMpIcpZQ08ABwBjw+5AM3sywrL3A2ab2Rwz2wQ8CZxcffFAvoKrivKAVUD5duR3LUyrzHTuPmdvDhvYiRuen8pzExbGHcm5Fq/WQiBpUPh3b6AbsBBYAHQPh9WnRzh9lYXhsER/BXYHFgNTgavMbJsG7CVdIqlIUlFxcXGEp3bNWVZGGn87dx9G9OvA/z072ZukcC7J6rqg7CfAJcAfaxhnwOH1LLumtgOqX0Z6DDApXFY/4HVJ/zWzNVvNZDYGGAPBPYvreV7XArTKTOfe8wq54IFP+NFTk8jOSOPoIV3jjuVci1TrFoGZXRL2HmtmhyV2RDtraCHQK+FxT4Jf/okuBMZaYDYwFxgUPb5ryXKzMrj/gkKG9Sjgiscn8vbM5XFHcq5FinKweHzEYdV9AvSX1Dc8ADwKGFdtmvnAEQCSugADgTkRlu1SRH6rTB6+cD9265zH9x+dwPjZfi9k5xpaXccIukraB8iRtJekvcPuUCC3vgWbWTlwBfAq8DnwtJlNl3SppEvDyX4NjJA0FXgTuNbM/D/dbaUgN5N/fG9/enfI5aKHiyia5xedOdeQZFbzLndJ5wMXAIVAUcKoUuAhMxub9HQ1KCwstKKiovondC3O8tKNjLrnQ5aXlvHY9/Znj15t65/JOQeApAlmVljjuNoKQcLMp5rZc0lJtgO8EKS2JSUbOOOeD1izoZwnLj6Awd29/UPnoqirEES5juA5ScdL+qmk0VVdw8d0rn7dCnJ4/HsHkJuVzrn3f8QXy/w+yM7trCiNzv0dOBO4kuCU0NOB3knO5VyterXP5fGLDyA9TZx930fMXbEu7kjONWtRzhoaYWbnAavN7JfAgWx9Wqhzja5vx9Y8/r39qag0zrn3QxasWh93JOearSiFYEP4d72k7gR3KeubvEjORdO/Sz6PXrQfa8vKOfu+D1lSsqH+mZxz24h6z+K2wB+AT4F5BO0GORe7Id0LePSi/Vm9bjPn3PsRy0s3xh3JuWYnSiG41cy+Ds8c6k1w5e9vkhvLuej26NWWBy/clyUlGzn3vo9YtW5T3JGca1aiFIIPqnrMrMzMShKHOdcU7NunPfefX8hXK9dzwYMfU1ZeEXck55qNpF1Z7FxjG7FbR/48ai+mLCzhT6/PijuOc81GXa2PHkNwZXFP4E8Jw0uB65OYybkdNnJoV87abxfGvDuHwwZ25oBdO8Qdybkmz68sdi3OurJyjr/zv2yuMF6+6hAKcjLjjuRc7HbqymKCs4bOlnS9X1nsmoPW2RncfuaeLF2zkZtemBZ3HOeavCiF4AWCW0yWA+sSOuearL12aceVh+/G85MWM87vcOZcneo6RlClp5mNTHoS5xrYFYftxtszi7nxn1Mp7N2O7m1z4o7kXJMU6cY0koYlPYlzDSwjPY07ztyT8krjmmcmU1npdzl1riZRCsHBwARJMyVNkTRV0pRkB3OuIfTp2JrRJwxm/JcreeD9uXHHca5JirJr6Nikp3Auic7ctxdvfL6cW1+ZycH9OzKoq9/DwLlEUe5H8BVBa6OHh/3ro8znXFMhiVtOHUabnEx+9OQkNm72q46dSxTlfgQ3AdcCPwsHZQL/SGYo5xpah7xsbj1tGDOWlvLH12bGHce5JiXKL/tTgJMITxk1s8VAfjJDOZcMhw/qwrkH7MK9/53L+Nkr4o7jXJMRpRBssuDyYwOQ1Dq5kZxLnhuOG8yuHVtz9TOTKVm/Oe44zjUJUQrB05LuAdpKuhh4A7g3ubGcS46crHTuGLUnxaVl3PD8VOprYsW5VBDlYPFtwLPAc8BAYLSZ/SXZwZxLluE92/KjI/vz0pQlvDDJrzp2Lsrpo5jZ68DrSc7iXKO59Jv9eGtmMT9/YRqFfdrRs523rO5SV133I3gv/FsqaU1CVyppTeNFdK7hZaSncfsZe1JZaVz99GQq/Kpjl8JqLQRmdnD4N9/M2iR0+WbmV+S4Zm+XDrncdNIQPpq7ivv+OyfuOM7Fps5jBJLSJHk7vq7FOn2fnowc0pXbXpvJ9MUlccdxLhZ1FgIzqwQmS9qlkfI416gk8dtvD6NtbhY/fsqvOnapKcrpo92A6ZLelDSuqkt2MOcaS/vWWdx2+h7MWraWW16ZEXcc5xpdlLOGfpn0FM7F7JsDOnH+gb158P15HD+sG4V92scdyblGE2WLYD5BkxJ5wAIze8fM3kluLOca37XHDqJbQSt+/sJ0P4vIpZS6Th9tI+lpgiuJvwt8D3hD0jOS/Kwh1+LkZmVw4/GD+XzJGh776Ku44zjXaOraIrgT+Azob2bfNrNTgH7AVOCvURYuaWR4Q5vZkq6rZZpDJU2SNF2Sb2m4WB03rCsj+nXgtldnsnJtWdxxnGsUdRWCg8zsF+GZQwBY4FfAgfUtWFI6cBfBjW0GA2dJGlxtmrbA3cBJZjYEOH0H1sG5BiOJX508hPWbKrj1FW+u2qWGugqBdnLZ+wGzzWyOmW0CngROrjbN2cBYM5sPYGbLd/I5ndtpu3XO57sH9+WpogVMnL867jjOJV1dheB9SaMlbVUQJP0c+DDCsnsACxIeLwyHJRoAtJP0tqQJks6raUGSLpFUJKmouLg4wlM7t3N+eER/OudnM9oPHLsUUFchuBIYBsyW9JykZyV9CewBXBFh2TVtUVT/j8oA9gGOB44Bfi5pwDYzmY0xs0IzK+zUqVOEp3Zu5+RlZ3DD8bszdVEJT32yoP4ZnGvGar2OwMzWAKdL6kewj1/AtWb2ZcRlLyS413GVnkD1Nn8XAivMbB2wTtK7BIVmVsTncC5pTtqjO499NJ9bX53BsUO70q51VtyRnEuKKPcj+NLMXjSzcdtRBAA+AfpL6ispCxgFVL8i+QXgEEkZknKB/YHPt+M5nEuaqgPHpRvLuc3vc+xasCgXlO0QMysn2IX0KsGX+9NmNl3SpZIuDaf5HHgFmAJ8DNxnZt7InWsyBnVtw3kH9ubxj+czdaE3SudaJjW3W/UVFhZaUVFR3DFcCinZsJkj/vg2PdvlMvayEaSl7ewJdc41PkkTzKywpnHeDLVz9SjIyeS6Y3dn0oKvefbThXHHca7BeTPUzkXw7b16sE/vdtzy7xmUrN8cdxznGpQ3Q+1cBGlp4pcnDWH1+k3c/oaf1OZaFm+G2rmIhvYo4Jz9e/PIB/M4o7AXg7t724uuZYhy+ug7wDwgM+z/BPg0ybmca5KuPnoAbXOzuGncNJrbiRbO1abeQiDpYuBZ4J5wUA/g+WSGcq6papubxU+PGcgn81bz/KRFccdxrkFEOUZwOXAQsAbAzL4AOiczlHNN2RmFvdijV1t++/IMSjf6gWPX/EUpBGVh66EASMpg2zaDnEsZaWniVycNYcXaMu5444u44zi306IUgnckXQ/kSDoKeAZ4MbmxnGva9ujVllH79uKh8fOYtaw07jjO7ZQoheA6oJjgzmTfB14GbkxmKOeag/87ZhB52RmMfsEPHLvmLcpZQ5Vmdq+ZnW5mp4X9/ql3Ka996yyuOWYgH85ZxYtTlsQdx7kdFuWsoRMkTZS0StIaSaWS1jRGOOeaurP324WhPdpw878+Y11ZedxxnNshUXYN3QGcD3QwszZmlm9mfiWNc0B6mvjlSUNZtqaMO//jB45d8xSlECwApvnuIOdqtk/vdpy2T0/u/+9cZi/3A8eu+YlSCH4KvCzpZ5J+UtUlO5hzzcl1xw6idXYG14+dRqXf49g1M1EKwc3AeqAVkJ/QOedCHfOyueH43fl43ioe/3h+3HGc2y5RGp1rb2ZHJz2Jc83c6fv05PmJi/j9v2dw5O5d6FrQKu5IzkUSZYvgDUleCJyrhyR+9+1hlFdWcuPzfm2Baz6itjX0iqSN4amjfvqoc7Xo3aE1PzlqAG98vox/T1sadxznIolyQVm+maWZWauw308fda4O3z2oL0N7tGH0C9P5ev2m+mdwLmZRtgiQdJKk28LuhGSHcq45y0hP4/ffHs7q9Zv47cufxx3HuXpFubL498BVwGdhd1U4zDlXi6E9Crj4kF15umgh789eEXcc5+oUZYvgOOAoM3vAzB4ARobDnHN1+NGR/enTIZfr/zmVDZsq4o7jXK0i7RoC2ib0FyQjiHMtTavMdH777WF8tXI9d/gN710TFqUQ/A6YKOkhSQ8DE4DfJjeWcy3DiH4dGbVvL+797xymLSqJO45zNYpy1tATwAHA2LA70MyeTHYw51qKnx27Ox3ysvnps1PYXFEZdxznthF111AasAJYDQyQ9I3kRXKuZSnIzeTXJw/hsyVruP+9uXHHcW4b9TYxIekW4ExgOlD1c8aAd5OYy7kWZeTQbhwzpAu3vz6LY4Z0pW/H1nFHcm6LKFsE3wIGmtnxZnZi2J2U7GDOtTS/OnkoWelp/GzsFG9+wjUpUQrBHCAz2UGca+m6tGnFz47bnQ/nrOLpogVxx3Fuiyitj64HJkl6EyirGmhmP0xaKudaqFH79uKFSYu4+V+fc9jAznRu4y2UuvhF2SIYB/waGE9w6mhVVy9JIyXNlDRb0nV1TLevpApJp0VZrnPNVVpa0ELpxvJKbho3Pe44zgERtgjM7GFJWcCAcNBMM9tc33yS0oG7gKOAhcAnksaZ2Wc1THcL8Or2hneuOdq1Ux5XHdGfP863MToAABSESURBVLw6k1emLWXk0K5xR3IpLkpbQ4cCXxB8qd8NzIp4+uh+wGwzm2Nmm4AngZNrmO5K4DlgedTQzjV3l3xjV3bv1obRL0yjZEO9v6ucS6oou4b+CBxtZt80s28AxwC3R5ivB8GN76ssDIdtIakHcArw97oWJOkSSUWSioqLiyM8tXNNW2Z6GrecOowVa8u45ZUZccdxKS5KIcg0s5lVD8xsFtHOIlINw6qfM3cHcK2Z1dkil5mNMbNCMyvs1KlThKd2rukb3rMtFx3cl8c/ms+Hc1bGHcelsCiFoEjS/ZIODbt7iXaweCHQK+FxT2BxtWkKgSclzQNOA+6W9K0Iy3auRfjxUQPo1T6Hn42dysbN3kKpi0eUQnAZwVXFP+R/9yW4NMJ8nwD9JfUNDzaPIjgDaQsz62tmfcysD/As8AMze3478jvXrOVmZfC7U4Yzd8U67nzzi7jjuBQV5ayhMuBPYReZmZVLuoLgbKB04AEzmy7p0nB8nccFnEsVB/fvyGn79GTMu3M4YXh3Bnf3O8G6xqXaLnWXNKWuGc1seFIS1aOwsNCKiorieGrnkmb1uk0cdfs7dG+bw9jLRpCRHrU9SOeikTTBzAprGlfXp60SqAAeBc4ATqzWOecaSLvWWdx04hCmLCzhofHz4o7jUkythcDM9gTOAvKAx4GbgSHAIjP7qnHiOZc6ThjejSN378xtr81k/sr1ccdxKaTO7U8zm2FmN5nZ3sCLwCPAjxslmXMpRhK//tZQMtLSuOH5qd5CqWs0dRYCST0kXS3pPeBcgiLwt0ZJ5lwK6laQw7UjB/LfL1bw3KeL4o7jUkSthUDSOwRbAZnABcD5wL+ALEntGyWdcynonP17U9i7Hb9+6TOKS8vqn8G5nVTXFkFvoB3wfeA1oCjsJoR/nXNJkJYmfn/qMDZsquBXL31W/wzO7aRaryMIL/JyzsVgt875XHH4bvzp9Vl8a8/uHLF7l7gjuRbMT1Z2rom69Jv9GNAljxufn8basvK447gWzAuBc01UVkYavz91OEvXbORWb6HUJZEXAueasL13accFI/rw6IdfUTRvVdxxXAsVqRBIOljShWF/J0l9kxvLOVflmqMH0r0gh+vGTqWs3FsodQ0vyh3KbgKuBX4WDsoE/pHMUM65/2mdncFvThnK7OVrufutL+OO41qgKFsEpwAnAesAzGwxkJ/MUM65rR02sDPf2rM7d789m1nLSuOO41qYKIVgkwXXuhuApNbJjeScq8noE4eQ3yqTa5+bQkWlNz/hGk6UQvC0pHuAtpIuBt4A7k1uLOdcde1bZzH6hMFMnP81j34wL+44rgWptxCY2W0Edw97DhgIjDazvyQ7mHNuWyfv2Z1vDujEra/OZNHXG+KO41qISGcNmdnrZvZ/ZnaNmb2e7FDOuZpJ4uZThgJw4z+9hVLXMKKcNVQqaU21boGkf0ratTFCOuf+p2e7XP7vmIG8NbOYcZMXxx3HtQBRtgj+BPwf0APoCVxDcIzgSeCB5EVzztXmvAP7sGevtvzyxc9YtW5T3HFcMxelEIw0s3vMrNTM1pjZGOA4M3uKoHVS51wjS08Tt5w6nDUbNvMbb6HU7aQohaBS0hmS0sLujIRxvoPSuZgM7JrPDw7tx9iJi3hnVnHccVwzFqUQnAN8B1gOLAv7z5WUA1yRxGzOuXpcfvhu7NY5j2uemczC1X6fY7djopw+OsfMTjSzjmbWKeyfbWYbzOy9xgjpnKtZdkY6fztnbzZuruDCBz+hZMPmuCO5ZijKWUOtJF0u6W5JD1R1jRHOOVe//l3yuec7+zBv5ToufXQCm8or447kmpkou4YeBboCxwDvEJw55I2dONeEjOjXkVtOHc4Hc1Zy3XNT/PoCt12iFILdzOznwDozexg4HhiW3FjOue317b178pOjBjB24iJuf+OLuOO4ZqTWexYnqNrp+LWkocBSoE/SEjnndtiVh+/GwtXrufPNL+jZLoczCnvFHck1A1EKwRhJ7YAbgXFAHvDzpKZyzu2QoAmKYSwp2cj1Y6fSvSCHg/t3jDuWa+Lq3DUkKQ1YY2arzexdM9vVzDqb2T2NlM85t50y09O4+5y92a1zHpf9YwIzlq6JO5Jr4uosBGZWiV8r4Fyzk98qkwcv3Jfc7HQufPATlpZsjDuSa8KiHCx+XdI1knpJal/VJT2Zc26ndCvI4YEL9mXNhs1c+NAnrC0rjzuSa6KiFILvApcD7wITwq4oysIljZQ0U9JsSdfVMP4cSVPCbrykPbYnvHOubkO6F3D3ufswa1kplz/2KZsr/BoDt60oVxb3raGrt/lpSenAXcCxwGDgLEmDq002F/immQ0Hfg2M2f5VcM7V5ZsDOnHzt4byzqxiRr8wza8xcNuIcmVxrqQbJY0JH/eXdEKEZe8HzA6bqNhE0Gz1yYkTmNl4M1sdPvyQ4GI151wDG7XfLlx+WD+e+HgBd7/9ZdxxXBMTZdfQg8AmYET4eCHwmwjz9QAWJDxeGA6rzUXAv2saIekSSUWSioqLvZVF53bENUcP5OQ9u/OHV2fywqRFccdxTUiUQtDPzG4lvLDMzDYAijBfTdPUuE0q6TCCQnBtTePNbIyZFZpZYadOnSI8tXOuOkncetpw9u/bnv97ZgofzlkZdyTXREQpBJvCJqcNQFI/oCzCfAuBxMsaewLb3FdP0nDgPuBkM/NPpnNJlJ2RzpjvFNKrfQ6XPFLE7OXebJiLVgh+AbwC9JL0GPAm8NMI830C9JfUV1IWMIrgyuQtJO0CjAW+Y2aztie4c27HFORm8tCF+5GVkcYFD35CcWmU33WuJYty1tBrwLeBC4AngEIzezvCfOUEF6O9CnwOPG1m0yVdKunScLLRQAfgbkmTJEU6LdU5t3N6tc/l/vP3ZcXaMi586GOWr/ELzlKZ6juVTNI4ggIwzszWNUqqOhQWFlpRkdcL5xrCf2Ys4/LHJpLXKoO7zt6b/fr6taItlaQJZlZY07gou4b+CBwCfCbpGUmnSWrVoAmdc7E4fFAXnr/8IPKyMzjr3g+5/725fp1BCoqya+gdM/sBsCvBBV9nENy/2DnXAgzsms8LVxzEEYM68+uXPuPKJyayzpujSClRtggIzxo6FbgU2Bd4OJmhnHONq02rTO75zj5cO3IQL09dwrfuep8vi9fGHcs1kihXFj9FcLD3cIImI/qZ2ZXJDuaca1ySuOzQfjx60f6sXLeJk//6Pv+euiTuWK4RRL2yuJ+ZXWpm/wEOlHRXknM552Jy0G4deenKg+nXOY/LHvuU3738OeXeWF2LFuUYwSvAMEm3SJpH0LzEjGQHc87Fp3vbHJ7+/gGce8Au3PPuHM69/yO/3qAFq7UQSBogabSkz4G/ElwpLDM7zMz+0mgJnXOxyM5I5zffGsYfT9+DifO/5sS/vMeEr1bXP6NrduraIpgBHAGcaGYHh1/+FY0TyznXVJy6T0/G/mAEWRlpjBrzAY98MM9PMW1h6ioEpwJLgbck3SvpCKI1Nueca2GGdC/gxSsO5pD+nRj9wnR+8vRkNmzy34UtRa2FwMz+aWZnAoOAt4EfA10k/U3S0Y2UzznXRBTkZnLfeYVcfdQAnp+0iFPufp95K2JvbMA1gCgHi9eZ2WNmdgJBC6KTgG1uO+mca/nS0sSVR/TnoQv3Y+majZz4l/e4++3Zfj/kZq7etoaaGm9ryLmmYcGq9Yx+YRpvzSymXW4mF39jV847sA952RlxR3M1qKutIS8EzrmdMmnB1/z5jVleEJo4LwTOuaTzgtC0eSFwzjUaLwhNkxcC51yj84LQtHghcM7FxgtC0+CFwDkXu8SC0DY3k4sP2ZVz9+9NQW5m3NFSghcC51yTkVgQsjPSGDm0K2cU9uLAXTuQluaNFySLFwLnXJMzfXEJT32ygOcnLmLNxnJ6tc/h9H16cdo+PeneNifueC2OFwLnXJO1cXMFr05fytNFC3h/9kokOKR/J84s7MWRgzuTnZEed8QWwQuBc65ZWLBqPc8ULeDZCQtZXLKRdrmZfGuvHpy5by8GdW0Td7xmzQuBc65Zqag03pu9gqeLFvD69GVsqqhkeM8CzijsxUl7dqdNKz/AvL28EDjnmq1V6zbx/MRFPF20gBlLS8nOSOO4Yd04YXg39uvbnnwvCpF4IXDONXtmxtRFwQHmcZMWU1pWTnqaGN6zgBH9OnBQv47s3bsdrTL9mEJNvBA451qUjZsr+HT+asbPXsn4L1cweWEJFZVGVkYa++zSjoN268CB/ToyvGcBmen1trafErwQOOdatLVl5Xw8d2VYGFby2ZI1ALTOSmf/XTswol8HDuzXgd27tknZaxXqKgR+jbdzrtnLy87g8EFdOHxQFyA4rvDhnGBrYfzslfxnxnIA2uVmcmC/DuzVqx1dC1rRraAVXdoEXVZG6m45eCFwzrU47Vtncdywbhw3rBsAS0o28MGXwdbC+NkreHnq0m3m6ZiXRdeCVnRt0yrhb87/Hhe0arHtI/muIedcSjEz1mwoZ+majUFXsoGlJWUsXbOBpSUbWVKykWVrNrJ6/eZt5s3PzqBLWCS6tGlF14LshP5geIe8bNKb4O4n3zXknHMhSRTkZlKQm8nArvm1TrdxcwVLS6qKRcLfsP/LL1ewvLSMisqtf0ynp4nO+dlBcQgLROc2QcHo2ibo75TXijY5GUhNo2B4IXDOuRq0ykynT8fW9OnYutZpKiqNlWvLthSJZVu2MspYtmYjs4vX8v7sFZSWlW8zb1Z6Gh3zsuiUn03HvGw65WfX2t86Kz2pRSOphUDSSODPQDpwn5n9vtp4heOPA9YDF5jZp8nM5JxzDSU9TXRu04rObVoxvGft060rK99SJIpLy1ixdhPFpWVhfxlLSjYyZVEJK9eWUVnD3vqczHQ65mdx3gF9uPgbuzb4eiStEEhKB+4CjgIWAp9IGmdmnyVMdizQP+z2B/4W/nXOuRajdXYGu3bKY9dOeXVOV1FprF6/aUuBSCwWxaVldG6TnZR8ydwi2A+YbWZzACQ9CZwMJBaCk4FHLDhi/aGktpK6mdmSJOZyzrkmKT1NdMwLdgk1pmSeONsDWJDweGE4bHunQdIlkookFRUXFzd4UOecS2XJLAQ1HdmovvcryjSY2RgzKzSzwk6dOjVIOOecc4FkFoKFQK+Exz2BxTswjXPOuSRKZiH4BOgvqa+kLGAUMK7aNOOA8xQ4ACjx4wPOOde4knaw2MzKJV0BvEpw+ugDZjZd0qXh+L8DLxOcOjqb4PTRC5OVxznnXM2Seh2Bmb1M8GWfOOzvCf0GXJ7MDM455+qWus3tOeecA7wQOOdcymt2rY9KKga+ijvHDugIrIg7RCPzdW75Um19ofmuc28zq/H8+2ZXCJorSUW1NQHbUvk6t3yptr7QMtfZdw0551yK80LgnHMpzgtB4xkTd4AY+Dq3fKm2vtAC19mPETjnXIrzLQLnnEtxXgiccy7FeSGIgaRrJJmkjnFnSSZJf5A0Q9IUSf+U1DbuTMkiaaSkmZJmS7ou7jzJJqmXpLckfS5puqSr4s7UWCSlS5oo6aW4szQULwSNTFIvgtt3zo87SyN4HRhqZsOBWcDPYs6TFAm3ZT0WGAycJWlwvKmSrhy42sx2Bw4ALk+Bda5yFfB53CEakheCxnc78FNquAFPS2Nmr5lZefjwQ4L7TbREW27LamabgKrbsrZYZrbEzD4N+0sJvhi3ubtgSyOpJ3A8cF/cWRqSF4JGJOkkYJGZTY47Swy+C/w77hBJEumWqy2VpD7AXsBH8SZpFHcQ/JCrjDtIQ0pqM9SpSNIbQNcaRt0AXA8c3biJkquu9TWzF8JpbiDYlfBYY2ZrRJFuudoSScoDngN+ZGZr4s6TTJJOAJab2QRJh8adpyF5IWhgZnZkTcMlDQP6ApMlQbCb5FNJ+5nZ0kaM2KBqW98qks4HTgCOsJZ70UpK3nJVUiZBEXjMzMbGnacRHAScJOk4oBXQRtI/zOzcmHPtNL+gLCaS5gGFZtYcWzGMRNJI4E/AN82sOO48ySIpg+Bg+BHAIoLbtJ5tZtNjDZZECn7NPAysMrMfxZ2nsYVbBNeY2QlxZ2kIfozAJdNfgXzgdUmTJP29vhmao/CAeNVtWT8Hnm7JRSB0EPAd4PDwvZ0U/lJ2zZBvETjnXIrzLQLnnEtxXgiccy7FeSFwzrkU54XAOedSnBcC55xLcV4IXIsiqSLhdMZJO9ISqKRCSXeG/RdI+msDZ5sm6cX6WmOVtGfiKZmSTkqFlk1d4/PTR12LImmtmeU14PIuILjw74oGWNaWbJIeBmaZ2c2N8dzO1cW3CFxKkDRP0i2SPg673cLhp4e/0CdLejccdmhNbc1L6i3pzfD+Cm9K2iUc/pCkOyWNlzRH0mkRIn1A2DCdpP3CeSeGfwdKygJ+BZwZbkWcmbh1UlsW53aEFwLX0uRU2zV0ZsK4NWa2H8EVz3eEw0YDx5jZHsBJ9Sz7r8Aj4f0VHgPuTBjXDTiYoF2l39e1kPD+BUcA48JBM4BvmNleYZ7fhs1ZjwaeMrM9zeyp7cji3HbxRudcS7PBzPasZdwTCX9vD/vfBx6S9DRQX8NpBwLfDvsfBW5NGPe8mVUCn0nqUsv8OZImAX2ACQQ37gEoAB6W1J+g1dLMenLUl8W57eJbBC6VWPV+M7sUuJGg9dBJkjrs4PLKEvprapYa/lekegNZwOXh8F8Db5nZUOBEgpYtt5cf7HM7zAuBSyVnJvz9AEBSPzP7yMxGAyvYujnp6sYDo8L+c4D3diSEmZUAPwSuCZtyLiBotRTggoRJSwka7UtaFufAC4FreaofI0jcX58t6SOCe87+OBz2B0lTJU0D3gXqunvcD4ELJU0haHlzh2/YbmYTw+caRbBb53eS3gfSEyZ7Cxhcw7GOBs3inJ8+6lJCKtz/wbkd5VsEzjmX4nyLwDnnUpxvETjnXIrzQuCccynOC4FzzqU4LwTOOZfivBA451yK+39SyKaxIbz91AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "def computeAvgLength(epsilonRatio): \n",
    "  N = 10000\n",
    "  polymer_length = 100\n",
    "\n",
    "\n",
    "  configuration_length = 0\n",
    "  initial_configuration = np.random.randint(2, size = polymer_length)\n",
    "\n",
    "  for step in range(N):\n",
    "    \n",
    "    M_initial = np.sum(initial_configuration)\n",
    "    configuration_length += M_initial\n",
    "\n",
    "    index_to_change = np.random.randint(0,polymer_length)\n",
    "    proposed_configuration = initial_configuration.copy()\n",
    "    proposed_configuration[index_to_change] = (initial_configuration[index_to_change] + 1) % 2\n",
    "    M_proposed = np.sum(proposed_configuration)\n",
    "    r = np.random.uniform(0,1)\n",
    "\n",
    "    if r < np.exp(-(epsilonRatio * (M_proposed - M_initial) ) ):\n",
    "      initial_configuration = proposed_configuration\n",
    "\n",
    "  return configuration_length/(N*polymer_length)\n",
    "    \n",
    "\n",
    "epsilon_ratios = np.linspace(-5.0, 5.0, 21)\n",
    "average_monomer_lengths = []\n",
    "for epsilonRatio in epsilon_ratios:\n",
    "  average_monomer_lengths.append(computeAvgLength(epsilonRatio))\n",
    "\n",
    "plt.plot(epsilon_ratios, average_monomer_lengths)\n",
    "plt.xlabel(\"Epsilon Ratio\")\n",
    "plt.ylabel(\"Average Monomer Orientation\")\n",
    "plt.title(\"Monomer Orientation For Varying Epsilon Ratios\")\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "name": "Markov Chains.ipynb",
   "provenance": []
  },
  "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"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}