ReCoDE_MCMCFF/docs/learning/01 Introduction.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "0aadee51-c8d7-4b72-b281-bf2daca49811",
   "metadata": {
    "tags": []
   },
   "source": [
    "<h1 align=\"center\">Markov Chain Monte Carlo for fun and profit</h1>\n",
    "<h1 align=\"center\"> 🎲 ⛓️ 👉 🧪 </h1>\n",
    "\n",
    "# Introduction\n",
    "\n",
    "Hello and welcome to the documentation for MCMCFF! These notebooks will guide you through the process of writing a medium sized scientific software project, discussing the decision and tradeoffs made along the way.\n",
    "\n",
    "## Setting up your environment\n",
    "\n",
    "It's strongly encouraged that you follow along this notebook in an enviroment where you can run the cells yourself and change them. You can either clone this git repository and run the cells in a python environment on your local machine, or if you for some reason can't do that (because you're an a phone or tablet for instance) you can instead open this notebook in [binder](https://mybinder.org/v2/gh/TomHodson/ReCoDE_MCMCFF/HEAD)\n",
    "\n",
    "I would also suggest you setup a python environment just for this. You can use your preferred method to do this, but I will recomend `conda` because it's both what I currently use and what is recommeded by Imperial.\n",
    "\n",
    "```bash\n",
    "#make a new conda environment from the specification in environment.yml\n",
    "conda env create --file environment.yml\n",
    "\n",
    "#activate the environment\n",
    "conda activate recode\n",
    "```\n",
    "\n",
    "If you'd prefer to keep this environment nicely stored away in this repository, you can save in a folder called env by doing\n",
    "\n",
    "```bash\n",
    "conda env create --prefix env --file environment.yml\n",
    "conda activate ./env #you have to run this in the enclosing directory of course!\n",
    "```\n",
    "\n",
    "Further Reading on how to set up an environment:\n",
    "- [Software Carpentry](https://carpentries-incubator.github.io/introduction-to-conda-for-data-scientists/02-working-with-environments/index.html) More practical.\n",
    "- [The Turing Way](https://the-turing-way.netlify.app/reproducible-research/renv.html) Discusses why we use environments.\n",
    "- [Essential Software Engineering for Researchers](https://imperialcollegelondon.github.io/grad_school_software_engineering_course/l1-01-tools-I/index.html) Quick overview.\n",
    "\n",
    "\n",
    "\n",
    "## The Problem\n",
    "\n",
    "So without further ado lets talk about the problem we'll be working on, you don't necessaryily need to understand the full details of this to learn the important lessons but I will give a quick summary here. We want to simulate a physical model called the **Ising model**, which is famous in physics because it's about the simplest thing you can come up with that displays a phase transition, a special kind of shift between two different behaviours."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e2e5299a-5e20-417f-a62e-ce47d018d542",
   "metadata": {},
   "source": [
    "I'm going to weave exposition and code here so don't mind if I just take a moment to impor some packages and do some housekeeping:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "ee600e16-506b-4676-8d84-16b415338191",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "# This loads some custom styles for matplotlib\n",
    "import json, matplotlib\n",
    "\n",
    "with open(\"../_static/matplotlibrc.json\") as f:\n",
    "    matplotlib.rcParams.update(json.load(f))\n",
    "\n",
    "np.random.seed(\n",
    "    42\n",
    ")  # This makes our random numbers reproducable when the notebook is rerun in order"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e52245f1-8ecc-45f1-8d52-337916b0ce7c",
   "metadata": {},
   "source": [
    "We're going to be working with arrays of numbers so it will make sense to work with `Numpy` and we'll also want to plot things, the standard choice for this is `matplotlib`, though there are other options, `pandas` and `plotly` being notable ones.\n",
    "\n",
    "Let me quickly plot something to aid the imagination:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "7b05be8f-9edb-4742-bbfc-e892cc09b82b",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAboAAAG/CAYAAAAuOhI2AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABDrAAAQ6wFQlOh8AAAhKklEQVR4nO3dvXLsyJEoYOhqI0RbwfUuH+F6kiVzXva8xJi0Vt4+AtcTH4BjbGgNUTfOIlpE/WRmgRXfZ00MgSp0N7rzAInM+t3f//73AwB29X9WHwAAZBLoANiaQAfA1gQ6ALYm0AGwNYEOgK0JdABsTaADYGsCHQBbE+gA2JpAB8DWBDoAtibQAbC1f5vZ+ddff/3DcRz/7ziOvx3H8d8hRwQA135/HMe/H8fxn7/88stvX204FeiOfwS5/5gcAwBG/fk4jr9+tcFsoPvbcRzHn/70p+MPf/jDcRzH8fLy8r82eHt76x50dozz/qPH0XtMj+b9ynmMjDkyXL2XI8c48l7Myngvr15H7/Yt483OuULLe9/7ulpEvN8944/MEfFdiH5dd/Hz6/jtt9+Ov/71r8fxGYe+8ruZhVd//fXX/3scx9tf/vKX4+np6R8D/u53/2ubkfFnxzjvP3ocvcf0aN6vnMfImCPD1Xs5cowj78WsjPfy6nX0bt8y3uycK7S8972vq0XE+90z/sgcEd+F6Nd1Fz+/jo+Pj+P19fU4juPll19++a+v9vMwCgBbm711eSniXw4r/lV1NWbEv6wjzF5djRzz1XsZccUQfQUeMebI/tFXWxnHMHJun2V8R6/2ybhyzRiz9zOtuJvxXe6gRB2HKzoAtibQAbA1gQ6AraXn6O7yhFB07qMiD9ii4knPqzlXPMG14hyJyGNUPHE3+5RyxHe2d/+ROa7eq4ynf89WPCndkkNd8WRt9Xfw+fn5+PHjR9O4rugA2JpAB8DWBDoAthaeo8vomtErouvCyOuYrZfpHS/CSK5rNpeVMUfL/rP5sJHcVctx9Wyf0VEkIyc38l5+h24dFZ1rroyc2yNjXqmo0f1qzJ86o1xyRQfA1gQ6ALYm0AGwNYEOgK2FPIzy8vJyvL+/P/zbSBL6Dsun9I73yOwcKx6qGXlIIKPg/w4P6lQUiF/NOTLHCisK31c8GFKxtFZGw+WMh+1WPPSkYBwAHhDoANiaQAfA1kJydG9vb8fT01PTthGLOlY09I0oSI4+pow5KsavKLTOWDgyowF2xbne+xlWNOMdmXO2qfN3aWwckceNbhYeUZQ+8p3OWrzaFR0AWxPoANiaQAfA1sKbOn+HGrhHKvJ8vfe0K+7FX42foaKpc0QOKKL2a0XeNiNf1isin9kro36z4rfn7A61eqto6gwAAwQ6ALYm0AGwtfJelxl9DyNk1LK07NNzTI/GmM25RdSCRbzu2Tqfu/YcnK3PrOgTOnIOrKhbzMg3V+S2onOmGb89FbXJWfu0cEUHwNYEOgC2JtABsLXwXpcVPSIr1rqK2D+6rq7F7L34ijq6jD56FfmAjPzMHT6Pitxjxdp+EZ/PbH1mxlpxK2rcKtZezDqOR1zRAbA1gQ6ArQl0AGwtvNflWcu9+dl7uy37R+c6Iu4/z9anPTL7OjPWCxz5fDLymdH5sJZcY3TOLePzaRHx3szK6E3aO2f09i1jtHx/ZuddUavcchxRn6krOgC2JtABsDWBDoCtCXQAbC28qXNFc9ezigTwyJwZyfKKhTyjVTRcjjqOme1Hjuu7LIw7e1wV72VEwfjs68h4EK5ioeIWswX7Kxd3dUUHwNYEOgC2JtABsLVv0dR5dvHKDBVNT1uK689WFLhejRmRj4k4R6KL0isWxh35fO5QrH11TKs+n1kRefeK79xZRQ414xyI4ooOgK0JdABsTaADYGvpTZ3PRu7NX8moR6vINY7kZ6IX8syo0cmoQYz4+2z+bKQuKDr/nDXGV+NF7LOi4XJEDiii8XrFMwXfYcHY6Dri5+fn48ePH037uaIDYGsCHQBbE+gA2Fp6jm5k0cDZnENFr76RPEjGwpGzOYSR9/qOtXpnFTVULSrqzyp6CEbnGjNqqDLyZRk576tjivhdqFjI+Gqf7Pfm4+PjeH19bRrHFR0AWxPoANiaQAfA1sJzdBU1IRXruK3oMbiqf+bsnGcV6wNmrO03W2MVsU/F66w4t0fOkd68XkQdau97k1H3eJbRKzb683ukt24uop62lSs6ALYm0AGwNYEOgK2F5+gq8msrcgzfJe9UsRZZr4h6wIj3P6LWq3eOqzEz6hyj8zWPzPajXXVOVPTRvRrjDnnAs5HzKuO5h4zfxONwRQfA5gQ6ALYm0AGwNYEOgK2FPIzy8vJyvL+/H8cR0yQ4enHRljmv5uj9+4iKBzBGtq9oCBvddLtljO9YhD6y4G+vlkLriqbBGQ0YoheArWi6veI8bBnzrKIQ3sKrAPCAQAfA1gQ6ALYWkqN7e3s7np6ejuOoWSQwYlHH2XvWGTmGjFxjRnFq72KvLa9rNqcQkSupKNRdUQhfkafNyFlfyWiCvqJJ+llGnnb2OYiRMTIWb7XwKgA8INABsDWBDoCtlS+8mlGjk1GrdyWjnqmi/m+kQfZsjjTi87pSkSMayQNW1JtF19G1qGjCPXsMEWOeRTTZXpH3u7KihndkYdyfqaMDgE8CHQBbE+gA2Fp6r8uziAUXr7Tkz2bvQa/ICz7SmwNaUftyNlKjE5Gzi84tjpxns59XS17j6pjusJhoizvU7lX8ns2eM1XukDv8mTo6APgk0AGwNYEOgK2l97ocEZ13atmm4riv5qzInayoQYzo4VnRh/LqGFr+Hp3HuEs9WkU/2tk5W6xYU61XxOuOrud8tM+VFZ/Pv+KKDoCtCXQAbE2gA2Brt+h1ebVNdG6lZcyKnN7ImNHHlXFvPmNNwoicT+8xVORKMuq0KmTkaaPPs9leihHHELVPrzuue5hxDrRyRQfA1gQ6ALYm0AGwNYEOgK2FP4xS0Yh1xcKFEQ/RZDQNjniYoVfFgpcrCndXNcvtEVEoH/E671DIvqKR9IjoYvqRh9YqvqPVi7VaeBUAPgl0AGxNoANga+UF4yMi7sXPjrFijozi+u+a25otJj6OmiYBs/mV79LYu3fOERmNInpVvHe9IhaWjvi8Zn97Kpqi/5MrOgC2JtABsDWBDoCtldfRRRjJpczmrjLuYVfUHJ6NjJndCLdlm97jzmgYG9GcejYPGPH5Xf09I9eVlXvpNVuLV7EY78jnsWLR56uc213qGI/DFR0AmxPoANiaQAfA1kJydC8vL8f7+/txHDG1SL0LYLaIrjcbWdQxOl/zSEaOaPY+ecWCsiP5s945WrbPWJC0d84VNW8VOe6ziIVxo+tOHx3zipq2WRl5wav9e/f5+Pg4Xl9fm/ZzRQfA1gQ6ALYm0AGwtfA6urOK2oqKe9oVtV9nI7VgEWtZXcnoQTj73kWs5VdRtzh73lXknSp6l46c2xn1gNH9T1edQxWf+Yp19b5iPToA+CTQAbA1gQ6ArYXk6N7e3o6np6fh/e+w1lVEDUh03q+iNuy7rEfXa8XaZBnHkZF3OovIb/YeU8braNm/d922qzEzvj+7rBc48tvTc56powOATwIdAFsT6ADYWnivyysZ6zVFzBFdI5JlRc3T2Yq834r1zCLOgTt+XhlrEkbUm/Xmzypqwc4i8k4t+3z197v02Z2tS53t26qODgA+CXQAbE2gA2BrAh0AW0tv6tySXOwtmBxJjM8WZVYUVo8knbOT6y1WPChyBxkFyVfjjbhjs+qWMStkL5Tbsk1G4/vZYxjZpvcYKj9vV3QAbE2gA2BrAh0AWytfePWROxQDR9zDrmh0W5FL7BXxXvbmHjMW8L2yokC5RXSeqaKoPWLx0Kv9I17HikVoI/KA0bnikTmzF8LV1BkAPgl0AGxNoANga+ULr47U0fVun3H/v3e8ljGzFyZsmWNkzqt9KuqCRva/w2KhKxbKXVGHuiKHejbyHb1DTWhE7n5F3rV6gWxNnQHgk0AHwNYEOgC2Vl5Hl3HfPGKOKxF5wIrFXe9Qa5SRY5jdvvU4VqvIZ1bk0+76Xkf3vF3R47Oix+pIrj5Cz/uvjg4APgl0AGxNoANgayE5upeXl+P9/f3h3yJycld/j1jHrdeqe9jR/eUy+jdG1IJ91xzD1ZwZ6yBG152uqHt8NG9F3nb2mDLGjPitmu1/2vIbmrEGXs9vizo6APgk0AGwNYEOgK19y/XoMnpEVqzzVpFDWHEMZxnr0UVY0e+0d46MnGl0fefoPnc0+15djRcxZ8sc2bn7784VHQBbE+gA2JpAB8DWBDoAtla+8GqLFYnS3iLoR8cQ/bDDyBi9xfQZjaNXNH9dsUhqxHFFPCiSsUht9D4VD7xknANnIw+SVDxsd7VNRhPuFQ+saOoMAA8IdABsTaADYGvhBeMZuZPefNnIoqi9x9BitintiIrcydUYI8d0h0L2s4z3YjZ3EtGsOqLhwpWMBgsR20c3tF6xiG1E7jEiz/6ditBd0QGwNYEOgK0JdABs7Vs0dY64Fxy9iGNG7VHEgqS9xzmSY+jNmd5h8cqWbSJyeNl1dC1m89Et20Qvltyyz9X2LWbzfndp5N3rDnV1Ed/RUa7oANiaQAfA1gQ6ALYWkqN7eXk53t/fh/ePviddUfP2SHaNTsucvSLygr3bj9RpReRvKuoYr0SfIy1W1OpF7D87Z0a+uXe8R2NE1AX3HkdE3vwOeb5RrugA2JpAB8DWBDoAthZeRxeRg5i9T96SY8jIjUTnmSJyWS1jzqqod4qoN4vuAxrxXkbUHs3uM9IrtmJtuIy8YHSNZ0R+s6KO7o5mX+fz8/Px48ePpm1d0QGwNYEOgK0JdABsLb3XZUTu5A693yLmOFvR+3J2vJExW8yuwVWRE1qRS1lRzznSgzC6j2jLHFfbZ8hYC7NiTbur7Udk9FTN+k65ogNgawIdAFsT6ADYmkAHwNbCH0apSNZeGSkY750zYwHGCiOJ74wHWmZVNL6teFgowuxiu1fjtcyRsThy7zFULOy5ouHCWcs5VXHezTaObjnPvvLx8XG8vr42beuKDoCtCXQAbE2gA2Br6QXjK+4njzSlzWjqPHv/PyKfGbE4YvR7FbEgY0QuOCL3G73PHZonrBqzoqnz1Tazn1dF3jYiv3m1f8biuyPjfXVOaOoMAJ8EOgC2JtABsLX0HF2L6HvzI7VhETVT0WNG3CfPWLR2ds6R/FlEfqBi4dXZWrCMPGDGe1fRQLmiiXbF65g9hhWN1iPqUs8iPs+fx1BHBwCfBDoAtibQAbC18BxdRF5j9p51Rp3WilzKSM/BDBULR17NGfGZRy+Em1Grt2IB4KtjGNkn4zzNyGVV5B4r3purOa9k1BGPzGnhVQAYINABsDWBDoCtheTo3t7ejqenp4d/u0OdynHE17hV1Jm0iM4tjvTXvBrjap2qkTGuth/dpmfOkTkyaiuv5oiuX2qZI0LFnLNjVOQiW2pKe79jI+fp7G/oyhyqKzoAtibQAbA1gQ6ArZX3uqxYI2pFH77jWNMP8OwOa65l5C9ntz+OmrXGZmsOR3J4FesgthxH9BzZ+c6oMXp9h9+JiN6+2Wt+Wo8OAD4JdABsTaADYGsCHQBbC3kY5eXl5Xh/fz+OI6fgrzeBP1JQ2TtHRKL7LCN5HlGAfHUMs+/tyHFULIo6YrZYfkVhfMWitiNNA67miFisN7oxxMjDQhW/Eyt+l6+2H2HhVQB4QKADYGsCHQBbu0VT54gmwL1mi1EzCpYzFl7NuP9/h4bXI3mP6ALXjJzqyHchOvf4aLwVC3lejVlRMN57TCPbZJzbFbngWRWfzz+5ogNgawIdAFsT6ADY2i3q6KJrPCoaLkcsTHi1f+u8PX+PqEG8UjHHWcTnG5FDza4tGnnvKhYsHalhuxpzVkaeKSO31VuHmmFFnfDsuaypMwB8EugA2JpAB8DWvkUdXYToGpyIhQlHxozY56v9K/oetuwzm7ddUcvXMuZsTVWLiu/T7Oe1In8WIWJh3F4Z511GTq669k6vSwD4JNABsDWBDoCtheToflZxbz56XaOqOTJqj6Lvi1fUIGaspxUxb0avy95jqqgHHNm/4r2ZPZcrckQVxxCRP8s4r3qN/Mb25AHV0QHAJ4EOgK0JdABsLb3X5VlEnUnE/f8V/eQizNYrfZd6porcbsa6hxXrm13NmVGDOJvfXJEjajFb25pRj5ZR81bxeVTUDY9+v1zRAbA1gQ6ArQl0AGwtvY7uLCMHEeGOua2ROpPZOVu2mX1vRvK0Fe9/Ru1kxXnU+32I+Hwq1lDLyHeu+P5EnwMtOdTocyLCyOv+6nXpdQkAnwQ6ALYm0AGwNYEOgK2FP4xypSWRuqLYMWLOqybNvQnjkYcC7vjwQ0Zj4ojPp3fOioT9XQuro7+DGU25I+boddcG8xkL/kZ/70d+Q3+mqTMAfBLoANiaQAfA1kJydG9vb8fT09PDv43cy4+4n3wlo5H07MKqI/fJZ1/HSNFmRVPas4q8bISehSOPI6Zp8B2K7XvnjGh+PLv9o+PozbOPqGi4kJHjvnLnxviu6ADYmkAHwNYEOgC2ll5HF9H0tHeOjHv1d12otaLOqiIHFP3+fpc6rYxFUTMaq0fLWNz1LrVgV3rzgCN59FkV9c4R50ArV3QAbE2gA2BrAh0AWwvJ0b28vBzv7+/HceTUz9yhziRjjt46uxERrzu6Ridi4cgVdXUtxzibzxw596Nr8yJyWxWL2M72km2Z8yzivcs4B3pF5B5nz4nZz8fCqwDwSaADYGsCHQBbK1+P7pEVNSAt+/TsH2HkXnxEHdbV36Prm0aMHPdZdA3VSK1R9OfVMucdrOidGFEL1rt/y3HM5uxajqvXihrEiDlauaIDYGsCHQBbE+gA2Fp4ji5ibbjoe8ERc7Zs33v/v6KP3pWIe+Irei2u+MxH5jiLzqW0zHk1xh16k47Medcxe2V8Zyu+gxX1zl95fn4+fvz40bStKzoAtibQAbA1gQ6ArQl0AGwt5GGUt7e34+npqWnbioUjKwosKx66WfE67tLUuWLh1YwHkq72yWhGPfveRTRajyiU7z3uiubuEe9lxYM6Fb+Z0Y3xZwvGNXUGgE8CHQBbE+gA2Fr4wqtnEYtVzm7fMkZG8fbVnCMF5b0NXyMWwr06hhUL40a8d7NjZjSlXbFQcYQdivFbxhz5fGbHGMlxX40RYfa3Z3YOBeMA8EmgA2BrAh0AWwuvo4toEHsWcf9/Rb1Z9DGscsfjHsl79OYMIvKAKxZiXZHzjm5o3nJcIzmg7N+ajDq6jCbcZxEL/K7IkbZyRQfA1gQ6ALYm0AGwtfCFV0dU5BSi7ydH3NMe2T6i5mZWxpi9n0dvfu3RHLP3/0c+n+icXes2PSp6X44cx3foG9oy5llGbj9jMd7e/PNIXjaLKzoAtibQAbA1gQ6ArYXn6Cruk0fUc8yOGVGrl9FfM7pWrHWbHiPr0WX0Gr0yckzRfUAj1gmrqENd0Ts24nVE17hFjFFRj5vxG9pr9nfAenQA8EmgA2BrAh0AWwtfj25FH8qI+8kZOYbZPFNFHjCiP2CvFbVhLXP01mmt6DlYkX8eGbNCRp62Ioc6a+Q7WnFuR9e+tszx8xjWowOATwIdAFsT6ADYmkAHwNbKF179LonvigLKkfcqY2HIK7OF1i3jRT8UkLEAZoTZh54qmgy3yHgdFQ/qRDcmuEsT6Ohzu6JgfDYWKBgHgE8CHQBbE+gA2Fr5wqsVi6ZmNHONWBBzxSKoGbmTjDxhRaHubOPuiHzn7OsYee8y3tvZnGpEDmhVY+hoEU2do9+rjMYeK7miA2BrAh0AWxPoANhaeFPnKxX3lyMWxIzIv1z9PSLvd5ZR23LHhT2v9l/VbHc2X7miTu67Nnm+EvFbE7Eg89UcFb89ZxnPSkScuz1zaOoMAJ8EOgC2JtABsLXwXpdnEfes79jjLkPFIqgVvS8r6p965xyR0RMy+u/HsSY/2esuC5Se3SG3W1Ej2rv/I7Pf65bf9Z7nHPS6BIBPAh0AWxPoANhaeh1dRj/Ainq0sxXruo3IqHHrnXNkzDv0ulzx+Zxl1LRFr5PYMkdGH8qIz6vit2X2GHr3jxij4vctekx1dADwSaADYGsCHQBb+xbr0fWOkbEeXYTv0I9upNfl1fYZVuRpW1T0A5yVkVNdkfM5q+jRmZGHivg8ImrYes3m1WffK3V0APBJoANgawIdAFsT6ADYWvjDKBWFvXdo0BzxwEvEe9X7YM6KxRJbVDyo0+uOzXczxsx42CHCigcsro4h40G43gWbI0Qsjjz72xO1TwtXdABsTaADYGsCHQBbC8/RReQ1Zu9hzy7oN3pMKwqrK/IYFTnTijxrRa539jMf+b6c96lYyDhin7Pe46rIM/Va8T60yGgmnvG7kPX8hSs6ALYm0AGwNYEOgK2F5Oje3t6Op6en4zhyalmi792PzNGioo5u9p71SP5m9v5+Ra3eSC4yI6c6+1oz8jMRvuMCpY9kn6st35+MxZGja/VGfgfORt7rnvfKwqsA8EmgA2BrAh0AWwvJ0b28vBzv7+/D+8/WY5xF3F++Q0/Cllq9XhE1brNzRHw+VyrqHCPqAzPqTs8q5oyYYzaPdJd62t59MmrvIn5reue8M1d0AGxNoANgawIdAFsrr6MbyTv13guO6AWXXQPSOmbvcV5tP2I2l5VRYxUx5h3yFivOkRXnXUWd3YjZc2Ckf+PsMbT8hq6oZz6L+P589do/Pj6O19fXpnFd0QGwNYEOgK0JdABs7ZZ1dL0i1lKKmGO23qzl79E1OBl1dRF9KHvHrFDRa3Fkzop8WHSP1ZE5r8asOK8y6hx79x/JA16NmbHe5tX+I/tYjw4AHhDoANiaQAfA1kJydF+p6K04chxXInInvTLqTmbHaxlz9u//at6ev7dsP9szdUXdVsV3IWONu4g5MnqRzuZII9aKa9mn13epNZ6dc5QrOgC2JtABsDWBDoCtCXQAbC28qfOVjCLoOzY0PY6cJPSKos1oFYvYRhT0Zzyk0Stjcdfe/Vu36fl7y3H0jrmimLti0dqr/R+NcbVPxcNCV9tX/ha5ogNgawIdAFsT6ADYWnhT54wGsxH5gd774t+lgW/FMawoDK3ICWXI+Mx7xosYoyLPNNJ4PSOvFL2AbERD+RUNsCsaYGjqDABJBDoAtibQAbC19Dq6lvu03yH/sqKhb8sYKxqtZhzDiubUKxZBnZWRP1uxMGtE0+0R0d+PjGcQRszm4Ebeh9451dEBQBKBDoCtCXQAbC28ju6spc6kot9cRZ6vYkHLXhULfUYvmjqyT8aCvhm5xuw6uxYRObuMBUxX5HQqFmSe/R3IyJ9FjNk754ivxnx+fj5+/PjRNI4rOgC2JtABsDWBDoCtla9H90jFmmrRveAeHXNFvdKVinWpMvIzV/ucVfRUvdKSZ5rNC7bM2Tvm7PqPI3O2/L3i+5P9WzOSe1xR63pW0YM4Og/48fFxvL6+Xo55HK7oANicQAfA1gQ6ALaWvh7d2Yr+da3bfPX3kR6DvduvqLvLqKep6GGXkQO62r73GFrmvONafxU58BYr+rj2HsOIO+SfI773K+rqRrmiA2BrAh0AWxPoANiaQAfA1kIeRvnKXRYmvBqzomj9rOKhjjs0Rx7Zf7YIfeQYKwper7aPeBDkOzzUlPH5ZLjD4qEVDRZG3LEo/V9xRQfA1gQ6ALYm0AGwtfKmzhW5kxW5rYgxM5rtRjSrvsOCmL2NcCNex0jOIbrRbcs5EF1gvDKX0jPGigLkigL+kWOI/k6u+h3o+cwtvAoAnwQ6ALYm0AGwtfQ6uoh79Rm5k9k5WvTe065o3nrWcp+9IscTbeR1rGgSfGVFzWFFA+aWMc7u0NT5SkV9Wkb+eeQ5h4yFV8+iPmNXdABsTaADYGsCHQBbC194dUR0PVNEnuksYmHCFffvI/Izs3OOmF0k9dH20TmeiDzGSD5mVkWf0Iic3Yp6wFkt35+KXrFnEX0oZ8+biO/oKFd0AGxNoANgawIdAFv7FnV0Zytq3s4y1m0byTVG55VG5qzoa9h7TCNzzNZxtezTO2fveC1jRh9Ti4zv6Mgc0fVmLe7Qi/RsRZ5v5HmAKK7oANiaQAfA1gQ6ALYWvh5dRA7oap+KnFFGPuYsIg8YnWOoWMcto49eRh3dyOuIroOLyMlFiF7fLON1VXxnz1bkHu/Q0/ORjJ6pUVzRAbA1gQ6ArQl0AGxNoANga+FNnTMWYKwoFO09hkf7RxdzVzxkUJHovpqzZZuIQt6K4uCzFYuFXh3DiOhFNe/yAM2K41jRBL13zoyC8haaOgPAAIEOgK0JdABs7RYF42cZRbYZObmWeaPnuJqzIgfU+xlHLIxb0WQ74r2MbhoccU6cVZzrIyrem4z8/9V4K3L3Gd+n2d/26Nz8x8fH8fr6ejnmcbiiA2BzAh0AWxPoANhaeh3diO9Qf9Zi1wVLM/KCK/IzV8cQ0QA7o17paozomreqOaJrvSIW9oxYkHlF0+aMusbohYkrakj/yRUdAFsT6ADYmkAHwNZCcnQ/i6jJubpnnVGnlWG2bitijuy6oZYxI+aIyGNU1FJGjxmRT7tS0aMwoxbs6hgy8s+9x9Ayxx1yWy3HEJ3fnK3Ve35+Pn78+NF0LK7oANiaQAfA1gQ6ALYW3utyxGwt2IgV98VX1B5V1ORE5yAijiFjzIyekN8lD7givzlbr/ld6+h6928Zc/Y5iEf7XMnOmep1CQCfBDoAtibQAbC18Dq6s4paljv0kBwRkXeq7BfXegwZ792KfpoVa/1F9BGdPY6MOUa+0yt61kbPmZHbqqgTvsO53bKNOjoAeECgA2BrAh0AWxPoANhaeVPnlkR3bxJzRZHtKtnH+V0axvbO2SKjAfbsAy0RD3P1HlPEHFcymlVXNIoeGT/6OzsyR8RvavRDT7NjKBgHgE8CHQBbE+gA2FpIju7l5eV4f3+PGOo4jpxFHK/m6P37Ixn33q98h+bUETnUKxm5rBERTYB7xq8as1dEw+WMxhGz52bFOZTxfMAd8rbRDeYVjAPAJ4EOgK0JdABsLbyOLmNhwoxao4x77RUNYqNVNI6uyKneVfRxV9Qgjh5HttlFU49jvp4s4rcoukF2hJE5MurqevK06ugA4JNAB8DWZm9d/v44juOPf/zj//8fHx8f3YM8Pz9PHsa1q+PKOIbeOc/bjxzT7JwR+1wd96P9M96Lq3l7xxx5r6K1vHdV80Zb8Tqiz+0MI+dt7z4Z37crI+fyz/v89ttv//zP31/N9buZe+2//vrrn47j+I/hAQBgzp9/+eWXv361wewV3X8ex/Hn4zj+dhzHf0+OBQCtfn8cx78f/4hDX5q6ogOAu/MwCgBbE+gA2JpAB8DWBDoAtibQAbA1gQ6ArQl0AGxNoANgawIdAFsT6ADYmkAHwNYEOgC2JtABsLX/AS8fBEb7/lYQAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 550x550 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "state = np.random.choice([-1, 1], size=(100, 100))\n",
    "\n",
    "\n",
    "def show_state(state, ax=None):\n",
    "    if ax is None:\n",
    "        f, ax = plt.subplots()\n",
    "    ax.matshow(state, cmap=\"Greys\", vmin=-1, vmax=1)\n",
    "    ax.set(xticks=[], yticks=[])\n",
    "\n",
    "\n",
    "show_state(state)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a919be9-2737-4d79-9607-4daf3b457364",
   "metadata": {},
   "source": [
    "In my head, the Ising model is basically all about peer pressure. You're a tiny creature and you live in a little world where you can only be one of two things, up/down, left/right, in/out doesn't matter. \n",
    "\n",
    "But what *does matter* is that you're doing the same thing as you're neighbours. We're going to visualise this with images like the above, representing the two different camps, though at the moment what I've plotted is random, there's no peer pressure going on yet.\n",
    "\n",
    "The way that a physicist would quantify this peer pressure is to assign a number to each state, lower numbers meaning more of the little creatures are doing the same thing as their neighbours. We'll call this the Energy, because physicists always call things Energy, that's just what we do.\n",
    "\n",
    "To calculate the energy what we're gonna do is look at all the pixels/creatures, and for each one, we look at the four neighbours to the N/E/S/W, everytime we find a neighbour that agrees, we'll subtract 1 from our total and every time we find neighbours that disagree we'll add 1 to our total. Creatures at the edges will simply have fewer neighbours to worry about. \n",
    "\n",
    "I'll show you what the equation for this looks like, but don't worry to much about it, the word description should be enough to write some code. If we assign the ith creature the label $s_i = \\pm1$ then the energy is \n",
    "$$E = \\sum_{(i,j)} s_i s_j$$\n",
    "\n",
    "Ok let's do some little tests, let's make the all up, all down and random state and see if we can compute their energies."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "340a78e7-6c1d-4742-8dca-6d3bea492856",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABsQAAAGTCAYAAACf73fTAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABDrAAAQ6wFQlOh8AAA1l0lEQVR4nO3dMbLrzHUu0Jb1qsRYRYccwsvkEWg2HpDHwEkoZGRlbwgMxVJMBSq9wC77iLyngY29GyT7rBXJl0B3A2gA53cX9vebf/zjHw0AAAAAAABm9S+vHgAAAAAAAACMZEEMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaBTEAAAAAAACmZkEMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaBTEAAAAAAACmZkEMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaBTEAAAAAAACmZkEMAAAAAACAqf2fzM5/+tOfftda+7+ttb+01v5eMiIA9vbb1tq/ttb+3x//+Me/vXow78o7D2AK3nkreOcBTME7bwXvPICPF3rfpRbE2n+9MP4z2QYA7+HfWmt/fvUg3ph3HsA8vPP6vPMA5uGd1+edBzCHVe+77ILYX1pr7d///d/bX//612RTALzC73//+/Yf//Efrf33M51v/aW11v7whz+03/3ud68eCwAb/O1vf2t//vOfW/POW/L0zjudTv/z4/V6/aeNv/72K1+3j2wb1RtjRuTYK/uN2DrGzPiXrmVPpp/o/Ns6pl67mW2XbJ3HmfNUuW3VuYk+K7bOg5H3bOa5+dXWe9p/563mv/MAPlj0v/GyC2J/b621v/71r+12uyWbAuDFlIfo+3trrf3ud79rh8Ph1WMBIMc7r+/pnff1v/ce34NL/y34dfvItlG9MWZEjv1VfyNsHWNm/Jn/H0Cmn+j82zqmXruZbZdsnceZ81S5bdW5iT4rts6Dkfds5rn5Veae/m/eeX3+Ow9gDqved/8yehQAAAAAAADwStkvxAAAAJjQ6XT6n68N/vGPf/zPv//mN78JtfN1+6/t/KqtyLZr+1xqK9JPpN1Hj9tGjDr2aD9bj2Fpv+ic2qo3j5fOW2TbzPWq6mevc5rpt3d80TmTeUatFb0eo+6XzHwDgJ/OF2IAAAAAAABMzYIYAAAAAAAAU1MyEQAAgCfX67UdDofF7UaVcassTxgprxYplfcupdki20aOPdNPT2TbTKnJTKnMTPnLyBzq7Rsp6xgtd/l1+1GlMjPtZEqzLvU76tgz2/bGFJlvyidSqaoMKMA78YUYAAAAAAAAU7MgBgAAAAAAwNQsiAEAAAAAADA1GWIAAAA8OZ1O7Xa7tdZimUiZ/KS17fzq90gGTyRLalSGUEZVVlHmHEfGFM012jrfRma8bc1ei86JrWMemYmWOZ617S4ZlT+WycKLWJozmeOT60SlTA6ouQh8Al+IAQAAAAAAMDULYgAAAAAAAEzNghgAAAAAAABTkyEGAABAVyR/KJPB86r8sd5vkfykTLZPJNfs0dZMsUgW1pK95kiv35F5cJFcnapzHMkFi2ZSjcrKehQ5FxGZTLQqkfty1LGu+Z2c6PN367sqk7W4V67cXpl7ACP5QgwAAAAAAICpWRADAAAAAABgahbEAAAAAAAAmJoMMQAAAIbZmpESzV7amqcSyYepygiL7pvJksq0W3V80ZyZqlyaSB5c5fzK5N193bcyO67nFXliv+o3co6rMt6WxhjpMzLfMvfL2ut+v9/b5XJZtS2tnU6ndrvdwvtl7p9R+47KFJMZBszAF2IAAAAAAABMzYIYAAAAAAAAU1MyEQAAgNWipZciZcOq+ln6ba+SaVtFS6ZFyuz12omMK3Ose5WazIiUv1va97t2frVtpHTe1pJ8j0aVI12y9T6sHEem/GXlOe8Z+dzkf12v13Y4HFpruXurqlTmqH0jz6/K8rAA78IXYgAAAAAAAEzNghgAAAAAAABTsyAGAAAAAADA1GSIAQAA8ORrnkpPZQZPpJ9eW5kcrXfMRKk8p728scy5qMpLivZbsd+atqrmaiY7alQ/mXmQOeeZeZ3JytqaN/j426h8pUxGHdudTqd2u92e/j2S5fcJ+2bu2Xd8PwJE+UIMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaDDEAAACefJensleWz5JM/tDWbKy9Mqqi/fT2zVyv3nkblZ8U9YprGckXy+SnZbKkRp3/ypy5UVlGS3N+VLZX5J6uzE/8zvF4bOfzedO+P9HX3Myq/MS97HUffkLmJsASX4gBAAAAAAAwNQtiAAAAAAAATE3JRAAAAJ6sLR+1VDJpazmp6H57lI/LlExbKjW1V9mqXkm+UWXCov1sbTtTlm5JpJThd/ut2XePUnqP/WSuR+U57o0pch4rS7H2RMqGLt0D71B276f7rkzwJ4g8Y0eVKgX4FL4QAwAAAAAAYGoWxAAAAAAAAJiaBTEAAAAAAACmJkMMAACAJ1/zVEblNPVyT6L99LJNMllZvf0e263K1RqVfRX57VciuU2ZcxFpd+u1rGon2s+jzLyN5ID1VGbJjcrrqxxTVU7bqCy5aKbbd8+O+/3eLpfL6jH+dN/lZu71rKhUlWc36vkL8Eq+EAMAAAAAAGBqFsQAAAAAAACYmgUxAAAAAAAApiZDDAAAgCdf81QyIhk8veySaK7OWiMzUbbmD0VzW7bmEVVmlVWex615apV5dr0xZdrKzONRczUzpsrsski7Vdcj0m50jkSeZ1+3rZyLrPc1N/Or6FysfK6OkLm3Ivu++3kAfi5fiAEAAAAAADA1C2IAAAAAAABMzYIYAAAAAAAAU5MhBgAAwJPv8lSWZPKtIvlJvbaWtu1l+6zdb2kMS/tG+o1mivW2jYxhVPbSqAykynYz+Te9ubh0LvbKfOvpjXHpeB719h11vUblmlUa1c/XYz8ej+18Pg/pZ0ZfczNHZdS9SlUO4F7PcogalTPLnHwhBgAAAAAAwNQsiAEAAAAAADA1JRMBAAB48rV8VM+ocn5L/WTaiuw3qpxUpCzdUlm9Xrm4vcphVZaEjNh63aPz9h3KLUWue2Upw6pyhHuV6IzIXOfIHKoqwbnUD/t7x2dFz8h7GPYSnYtVz2Pm4AsxAAAAAAAApmZBDAAAAAAAgKlZEAMAAAAAAGBqMsQAAABYLZrbUJXVkMlEimxbmSm09rfHfioztbbmmGVlrvvWcz4qvy4yhshvvxKZB1uPPWrUeezl3UXb3brvO97/S7+Nynjj16Jz5B2ziqoyAyv35WeqmjOZ+5JtPvHZ9x1fiAEAAAAAADA1C2IAAAAAAABMzYIYAAAAAAAAU5MhBgAAwNv7SfkPlflJvbZG5luNEslpqsp/evx96RxnMuq2ZuG9y/0RGUdkbi6di0zu0dZtl+yVp/TdnLnf7+1yuZT181NFng2fIHM8s50Lxqucb5l3+rtnWj2qepeO7HfUtduDL8QAAAAAAACYmgUxAAAAAAAApqZkIgAAAE9Op1O73W6ttXFlanrbZkrNLZV8GVX6L1P+rmdUKcDMOV7qJ1PC7h3KblWVzsuU/ou0FZnza7avaqc3xsz9/yhSpnJtO0siJS0zto7/eDy28/lcMoaf5hX3VkZm3kbeEZl++RleVYr53UryrVF1377qb8Sed78evhADAAAAAABgahbEAAAAAAAAmJoFMQAAAAAAAKYmQwwAAIAn1+u1HQ6Hp3+PZor0sn0iGQOVuScRkWyiSP7To8yYe/tWnuPI8W3NG1vatyrPLpO5sWRUfkflvpmMpFFG5RGNykEZldsSNSoTkXVGnfN3bPcdx7Tk3bOMMl717K68Xj/pmfUO996oTM1P4wsxAAAAAAAApmZBDAAAAAAAgKlZEAMAAAAAAGBqMsQAAADoimQVRfJ6qjKQltqtytiK5ilEMqsi56Iquyza7qhsr15bS9lko65tJO/uVVlymX567WRyACP7Pqqaq5FrWZmnEhnTo3fMdOO/VD4nM/O2at/M+yQzpkeZMUX23UvkWTdK5DyN6nfkPP40I+/Lqms76lnx7nwhBgAAAAAAwNQsiAEAAAAAADA1JRMBAAB4cjqd2u12a63VlWZbMrJ84Xf7VpWoi/YTKRuUKT2VKWkzqoRdpJ9IvyOvZUTmemwtSzmyFOPXtkeWlszsO6p81Cijyrx9bfd+v7fL5VLSLt/bWhb4Vftmyp7tNaa99h31fhxZ9vQV171yTFVljt/VK+6BjMyYPpkvxAAAAAAAAJiaBTEAAAAAAACmZkEMAAAAAACAqckQAwAA4Mn1em2Hw2Fxu71yv5b63dp2Za5GJEcrk4GUyeuougYjc1u25kHtNb9G5b88tl15jkfmx6wdU2Tf6D1Qda9l8uwyZsjOmVU0D2rrHIreO5F9q+ZX5f3yin2rnk+RPj+134jMe+snPftG/V0xytJ8GpXPtwdfiAEAAAAAADA1C2IAAAAAAABMzYIYAAAAAAAAU5MhBgAAwGZ7ZZcs9VuVExLp89FemQmZTKRM7szW8zoql+3x/45kk0Sz1rbOr8ycyYxpVI5WNJts69xcMir7KzJnln7/2lYmN+9VuWass1c+VM8nXPdRGY+ZrMVReVZ7ZXdl+h2Zn7p138p5vFe+VeXfA3uIvF8qn23v9ozyhRgAAAAAAABTsyAGAAAAAADA1CyIAQAAAAAAMDUZYgAAADw5nU7tdru11mI5TZnchlH7VuVqVGb7VOY2bL0+S2PaK4tta6ZTZtuoSM7JqEyUURlbkZyzvfJ6RubbbM2dy+SnRbx79spP8y65eZ9g1LtobZ9L/Vae/1HHE+l3r3PcG8OjynzOyPUZ9TyOGjXfRql8fn3Ss88XYgAAAAAAAEzNghgAAAAAAABTUzIRAACAJ9frtR0Oh9ZarnzRV1XtRPvJlN3ZWiovqmqMS9tWlSdc6idjVAmlyLa9cj+VJQarzvmrykMt6Z2LdymhWCVS8rVyDN+1dTwe2/l8LuuHX6ssf/vptp6Lkedp1PsxczyZ433Hc9yzV5ndvdp9136rjCqR3Pv9FefBF2IAAAAAAABMzYIYAAAAAAAAU7MgBgAAAAAAwNRkiAEAAPDkdDq12+3WWsvlNkT2jWRHvSKvJ9pOJBMpkgc1Km9s1HWOjn9r9lckdy4zvyJjyqjMH8tk+43KKqvMs4vknvT66c2L6JzZum/mHH/d936/t8vlsnpf1hmZJTmzV+UcZjKRtmaKVd3De+5bJTqGrdcnkzv1quyud8wM68k86zKZe3ucJ1+IAQAAAAAAMDULYgAAAAAAAEzNghgAAAAAAABTkyEGAADAk+v12g6Hw6uH8a3KzKq1+0YznaryupbslRMyKhvr8VxsPeeVWVhV1zJ6ziIZKb1tM/tWZn1szYN7FLlemQy+vXJPMs+Zte0ej8d2Pp83t8U6MsPWeUVm2Kv63SuDcuTfEVVtL12PV9w/lc/uTL8ZVddn1N+p7/5c9IUYAAAAAAAAU7MgBgAAAAAAwNSUTAQAAKArUlIlU45sr/JkPb0xRsvH9EooVZZF6vVTWfKmqjRgZEyVpfLWtrNmHD2ZEkRbr1e0nGdVWcq9ygaOKr9UWb4zU5qxJ3IPvHuZqhlF3hmzXZ9M6eK9zlPk/VjVb6Ykb6bk615jqjTq+ow6nszzuPK6bz2ezDkedTyveC76QgwAAAAAAICpWRADAAAAAABgahbEAAAAAAAAmJoMMQAAALr2ytGJ5Ahkcpoi2QtVGSlL7WbyLUZlr72D6Hz76hOONZPTVnV8kftlad9Xzbet98/I8Vc9zyL7fh3v/X5vl8tlc7usE3n+9n7/hOfVkq3nYmSG0Dtcn8zfGVXzIvMufdfrU9Vu5b036h7YK0cz0uc73Ftb+UIMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaDDEAAAC6erX+Mzlaj3pZZZF2X5G1kG2rKp+kMnOrMuesN4ZIP5lcpqr59ShyLjK5GlUZdY+/L41p67aP9soXq5ojj7+PzFrrbRtp6+v/Ph6P7Xw+d9sib69nx6fLnKeqfqN9Vj37XiUz/qq/3Sr33Ss3M9NW1Xtgad+tXjXed8uz9YUYAAAAAAAAU7MgBgAAAAAAwNQsiAEAAAAAADA1GWIAAAA8OZ1O7Xa7De2jMlNkax5BZSZCVR7EyHyrTCZS5BxHjq8q32qpn57IsVZe90y/e+V1VeX1ReZbZm5W5oJF2l2b7bWm3561bd3v93a5XDb3w3hVuYzvone/7JUZVmnrO+QTr2Xm+ozaN/M3x6vm29Yxv8t4e++taFtrf3sFX4gBAAAAAAAwNQtiAAAAAAAATE3JRAAAALoiZVBGlb9ZsrXfaGm23r6jto38HilHuLRvVbmizPFk+tnrukf2XWpra9mtaGmjqvJelaUAI/30fq+8f9buV7H9mjEsbf9uZal+ukyJ1NlEnnWjysftVYr1HUXnYuT6ZM7TqHvg0+6tdxnvu4xjNF+IAQAAAAAAMDULYgAAAAAAAEzNghgAAAAAAABTkyEGAADAk+v12g6HQ7qdSLZET2XuxKhsjFGZYpXZJZHMrVFZH5l+MjktVRkvI+fx1ry+pXMamV89kfEuiZzHyvulJ5OxF3l2RGydX8fjsZ3P5839kleZffdp9rpnR/rEMX8nMheXZPI53/088TP4QgwAAAAAAICpWRADAAAAAABgahbEAAAAAAAAmJoMMQAAAIaJZKREMneq8nqqcjOyIlllmTGP2jeT+7VX9lqVynPaEz1vvd9G5VtV5nXtlbPV6zOS01Y1hqxR85zxqjI2P8EnHM8njHGU3rNv9rnJz+MLMQAAAAAAAKZmQQwAAAAAAICpKZkIAADAatHyXFWl2XrtPrb9CSXtKkvN9dodVSqvsszh1lKMjzLnNHItK0sMru0zK1Je8VGkdFbmWkbO8ah7fFQJzsg5rjq2+/3eLpfL6n3Zn3KXvAtzkZ/EF2IAAAAAAABMzYIYAAAAAAAAU7MgBgAAAAAAwNRkiAEAAPDkdDq12+329O+RLJylfSvzlHr5Q5nMra1jWOon8nvkeDLneKmtyJiqcsGi/Ub23drvyJy5Xl7X0r6Rfre2G70HIuOoyrDJ5HVFZO6BiMosOcYblR0JGUvPp0guK3w6X4gBAAAAAAAwNQtiAAAAAAAATM2CGAAAAAAAAFOTIQYAAMCT6/XaDodDa60u12hJLwesUiQbY1QWVqStSP5Y5ngyMhloVdleSzLza6+ssq19Lu07KtvvUWQe9MYxMitu63n9hPynr2M8Ho/tfD6/cDQ8esc5w3gj/26qGsPIv7ng3fhCDAAAAAAAgKlZEAMAAAAAAGBqFsQAAAAAAACYmgwxAAAAVlvKG4pkFY3KzajMRIrkGlXljVVmeWTORe94I1llmTG8Km+sat9MptvIHK1IP1X3T0bkXGSeUT1L7faeb5X5Y7J+YC6RHNBRKp/78O58IQYAAAAAAMDULIgBAAAAAAAwNSUTAQAAeHI6ndrtdmutjSv9NcpSqZ+tZRBnKCnUu15LxxM5F5GSdnvJHPtWkbmYaXupnaqylZXnqXLOZEoort02Wlpy63yrKlF5v9/b5XLZ3Bawzbv8bZApnQsz84UYAAAAAAAAU7MgBgAAAAAAwNQsiAEAAAAAADA1GWIAAAA8uV6v7XA4PP17VRZRViZ/KJK9FOmz93smqywjc30iOUeV+WqvyBhb6jOSO7V2vzW/9/qpyvaKnO/Mda48F5m8rt55i2z7aFROW2Rufv3fx+Oxnc/n1WMAakTzBffI73rHMcGr+EIMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaDDEAAACenE6ndrvdWmvbs7CiIplVvd+XMne25g29SlVe0q9+78nknGXyrbbOoaVjrcppG5mjl8lt6/2WybfreYf7Y0lkHi9liq397VdtR/aN+K6f+/3eLpdLWT/AOpm8xMzzuOpvqGi/8Gl8IQYAAAAAAMDULIgBAAAAAAAwNSUTAQAAeHK9XtvhcHj690w5wiVVZdwiRpZi3FqiL3rsVSUgl7aNlEHMlN2rKm2YGVNP5NiXSvBlyh6OKru11FaVvUoKjlJVygyYy8iyur129+oXPp0vxAAAAAAAAJiaBTEAAAAAAACmZkEMAAAAAACAqckQAwAA4MnpdGq32621tl+2RCQT6VEk3+q7/aL7ZnKaMplCmRyqyBiq+olmYY3KKuvJzK9MXlqkn0y7VRlWS9eyar69KvMsMv7KedA79q3n+Hg8tvP5vHpMwOtV/m0A/JovxAAAAAAAAJiaBTEAAAAAAACmZkEMAAAAAACAqckQAwAAoCuSb/NoVCbSo61jrMxAelSZLbXWu5zjrXljv2pr65gi7UZyzkaNf03bW/vMZGNtbTe6bWSMo+bXyPtn6xj2ynAE9lGVIbj0rKjKZYTZ+EIMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaDDEAAAC6XpGF9WgpK2NrrkalUe0+yuSEVOUnVR7rqBytUdc9005lPtSoDKvKdiNtvcN9WZnBU5Xfs3X89/u9XS6Xzf0CNTK5ppHf9/obBD6dL8QAAAAAAACYmgUxAAAAAAAApqZkIgAAAKtVlqUbVV5xqd3e75FSckvHurWfaNmjSLnI3pgi5RZHXbslkZJ8mfJRkWtdWZpx6zmOzsWqeyCi8n6purYjS3+O0ptvX//38Xhs5/N5lzEB640qJfuq9zJ8Gl+IAQAAAAAAMDULYgAAAAAAAEzNghgAAAAAAABTkyEGAADAk+v12g6Hw9O/j8yoyGRUVWU8RfqtzNzqbbtXNlGkn0iWXPR49jreyBgiuVOj7pHemCrnzKsy0SLbVo1xVJZPpu2lvLe1z6/7/d4ul8vqMQDvRy4Y1POFGAAAAAAAAFOzIAYAAAAAAMDULIgBAAAAAAAwNRliAAAAPDmdTu12u7XWYvlJj7bmX0RzpvbIbYrmFlXmkfXslT/WmweVmVtbj2dp2945HpVblpkjmTmeycaqHFNvHixlZfV+2yvXrDeGXruRMWTJGAKA9XwhBgAAAAAAwNQsiAEAAAAAADA1JRMBAAB4cr1e2+FwaK3FSowtldL7KlMiLVJm7x2MLJH2te29SklWGlX+LnPsVSX6omOIlKXM9Ntrq3IOZZ4dkTFsLYcZea4s9Rl5JmXKd37CPQ3vaK+Sz8B784UYAAAAAAAAU7MgBgAAAAAAwNQsiAEAAAAAADA1GWIAAAA8OZ1O7Xa7tdZieTeRDJ698q5GZRNFZDKpKs9LJmstMubMOe9lovXaisynzNyL5NtF86C2Zm5FM/e2nrclo+6XJZH5tvXeW8p0q8ywG9EOAOALMQAAAAAAACZnQQwAAAAAAICpWRADAAAAAABgajLEAAAAeHK9XtvhcGit5XKnejK5TZl+tuaaRfPSeuctkm+1ZGtWVmWOVsaoMfb2zeRKRcaUOU+j8tMyovdAZIyZ48nktn0Vvfd6+269BtE8uFHPZ5hB1T396fdW5jx8+rHDI1+IAQAAAAAAMDULYgAAAAAAAExNyUQAAACenE6ndrvdnv49U8psqa2ITOm83m9V7WbaypRMWxrDqHPRa6uyBFxmzkRU9ZspMRhR2W5kPmXaipS0zI5j7Zh6/WTuw8y+0baA/5W5hz/93hr1dwV8Ol+IAQAAAAAAMDULYgAAAAAAAEzNghgAAAAAAABTkyEGAADAk+v12g6HQ2stlo0TUZnP0Wsrk/UzSiSrJJKn9Krxj2prZPbSVkvXY1R+WqSfyL2VydHZmscVFcnre9X9ksnCq8oqA743MpvwFSrff3vli8I78IUYAAAAAAAAU7MgBgAAAAAAwNQsiAEAAAAAADA1GWIAAAA8OZ1O7Xa7tdZimTtVOUHRdkbl9axtJzK+NWPsGZUplMmd6onOmcg8qBxHb9t3yNHqydyXI/NhtrY18ngiz4q17UTH+Kjq2ffV8Xhs5/M51Bbw3iqfUWv3lRnGbHwhBgAAAAAAwNQsiAEAAAAAADA1C2IAAAAAAABMTYYYAAAAXZEsiVdlcvX02h2ZjVGVhTUqo2pUFlm2n62ZdZXjH5Wn9qgqvytzX47K3FsaQ0/0eCL7Vt3zmXZ657wqc+9+v7fL5bJ5jPDTZJ4ro2SyLzPkhjEzX4gBAAAAAAAwNQtiAAAAAAAATE3JRAAAAJ5cr9d2OBwWt8uUSOz9HinjFu1nbTuZbZdkShtlytZlSgxuPf5RpRiX2s6UhKwqZRi1de5Gx7v1vI0szbj2t6iqMmhVpQyjY3rcdm1Jy+Px2M7nc3fMwPeqysVWiowp8uyGn8QXYgAAAAAAAEzNghgAAAAAAABTsyAGAAAAAADA1GSIAQAA0LVXDsXWbIxeO7/ydd9MVlEmoyqSgZY556Py1EbamllVmTuXyZZ6Rd5YdC5msuQiemMcdZ4iz46qPL7omDJtrc0Jut/v7XK5rG4XfppIruleMmN6l5wzeHe+EAMAAAAAAGBqFsQAAAAAAACYmgUxAAAAAAAApiZDDAAAgCen06ndbrfWWizb61FVHtSjSC5QRCRnKrJvpK3IeVraPpNz9mkymWGjjMx0qcoBi95rW1XmAo7K/qq8L3v7Pm6b6RcYY6+sxe/6jPa7Nl8QfjpfiAEAAAAAADA1C2IAAAAAAABMTclEAAAAuqrK7iyVOasqcVdZmrHX7qgScJmylBl7lYRcEpkHVdtWlpbsza/K8zKqn0eRfqpKjkZKCvbaycqUi936bJm5rCl8kt69OLIc7lafNl54FV+IAQAAAAAAMDULYgAAAAAAAEzNghgAAAAAAABTkyEGAABAVyZrYmvuVCRf7HHfyqyySIZQJDMs+vvabUdlbC3JHHvknPf2jc6ZnkheV09mziypyoDJ3C8jM/gybUX23fpbpM/IPI767n45Ho/tfD5vbhdmF3kmvSqDK/MOr8p/hNn4QgwAAAAAAICpWRADAAAAAABgahbEAAAAAAAAmJoMMQAAADbLZFY9GpXXNSoXbEkkC2vrto/bV+aAZfoZlV1W1W70PFVl4VX2E2nn8dz05lgmDy6jah4/yuTq9NrZa0yP1o7xfr+3y+WycYRA5Dn5igyuyPMrmmMoU4yZ+UIMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaDDEAAAC6XpEplMnyqswFi2R79cYRyQXK5D/12l2SySOpzGLa2m8mO25pjK/IaVsSycKrHPPa3x77rcwQjOxblQM0cs6PyF48Ho/tfD539wW+F3l+jcpAHPUOj+wrT4zZ+EIMAAAAAACAqVkQAwAAAAAAYGpKJgIAAPDker22w+EQ3i9SgixTnmxp+61tjSoPGS3RF+k3UxYpIlKWbmupqei+PSOvV2/bqtJ/kX6jpbLe4Rw/qirNGJEpe1j1XHncd9SxAjGZZ+Mr+nzVvvBpfCEGAAAAAADA1CyIAQAAAAAAMDULYgAAAAAAAExNhhgAAABde2VLVOU2LW27NXOr8jz08oiieVBVls5pJl+pJ5Of1uuzKi9tZD+j8tMejbo+mczAVx17bwyZdjP5g737f+uY7vd7u1wum9uCnybzrKt6ZmWek5m/K171dwe8gi/EAAAAAAAAmJoFMQAAAAAAAKZmQQwAAAAAAICpyRADAABgtcqciUwmVdU4ltrJZPtE9s3kaL0iuyTSVibDrdfu0m+j8q0yGU9L/b5i24jK/LSltnv9RO7TUWOuzDzM5Kd9t+/xeGzn8znUFrDOyGdhpN+K3371+17HA6/gCzEAAAAAAACmZkEMAAAAAACAqSmZCAAAwJPT6dRut1trLVdqLrJtpmxYr5+eSCnGpTH1+o2UI6rsp/d7pizlyPKEa9tZEhl/pch5evSOZaoi523UPf0oc59ufUbtVVKsqjTr/X5vl8ulbFzwk73q2Zzpd9S7FT6dL8QAAAAAAACYmgUxAAAAAAAApmZBDAAAAAAAgKnJEAMAAKArk/1TmRu0td1e9k9lBlplTlik363ZZXtlOi2Naamtnq25YNExRPqpyk+LqMy3GpXBt9RPJqtslFH3cKataNvAz+TZAL/mCzEAAAAAAACmZkEMAAAAAACAqVkQAwAAAAAAYGoyxAAAAHhyvV7b4XBorfWzfUbl6FRmLUWyiiqzlzJ5ZBFbj2evnKZoP5HzFmm3KmOsMjcrk732inyYpTlRlW8X3XbrnIkcT/RZkcnR6/XTu19kBgFAny/EAAAAAAAAmJoFMQAAAAAAAKZmQQwAAAAAAICpyRADAADgyel0arfbLbxfL+8mkxm01E/E130rc5l6bUVyspZEzvHWdtb8Htk3onfeItuOzKHrtVu1bSZzL9JvNBcsIpLBlcnCitxrX7eN3gNV2z561bMQAH4aX4gBAAAAAAAwNQtiAAAAAAAATE3JRAAAAJ5cr9d2OBxaa9vLHj7+Hill9iizb6+tSFm9SMm3JZWlGSP9bC2vmOknep72Ks3YE2k3M4d6vy8da6b0Z+QcV7WbESltuPVZsNTPu5SlXDvfjsdjO5/Pq9sFgJ/AF2IAAAAAAABMzYIYAAAAAAAAU7MgBgAAAAAAwNRkiAEAANAVycJZ286SpX4ieWOjcqYeVeWCRfKSeu38yqhz8apz3lOZUddrp/LYt2ZWRceUycaK2DpXM31WHnvP0jnfK1Pwu33v93u7XC6b2wWAGflCDAAAAAAAgKlZEAMAAAAAAGBqFsQAAAAAAACYmgwxAAAAnpxOp3a73Vpr/5xRs1feULSfyL69dnpjitqa8RTNpMqc87VjWlKZGbb1+u2V6ZQ5T0v7Vp3zpXa2Hl9mbu6V17W0be+3ymdHTyTXLJKJ9tXxeGzn83nzGJnHqHfRO2ZHwqeoej8S5wsxAAAAAAAApmZBDAAAAAAAgKkpmQgAAMCT6/XaDodDup29SpD12smUlnuHfivL7PXOcabMXqafSGnDTCnDrWX1luxVdm+vUqCV98uoe6/y+kX6WftbZT9b973f7+1yuWxul88VKRM6al8l4PjpIu8l5Uf35QsxAAAAAAAApmZBDAAAAAAAgKlZEAMAAAAAAGBqMsQAAAB4cjqd2u12a63V5TZFZNrJZC/slbGV6XdUXlJEJtcsk/X16OvvkSyyvc5h9P6IHE9Vxkik3eh1r8oUqjz23jnOyNzDe+UAQuW79R3eRfAp3D/vwxdiAAAAAAAATM2CGAAAAAAAAFOzIAYAAAAAAMDUZIgBAADw5Hq9tsPh0FqL5QCNykRYyl7YmquzNN5e3lDkXGSysDJGZTo9imRhZfatOp5Mttejqgy06L4RozKqRt3vmWzCylzDUflKledRbhij5nTlvvDTZe5Tz/lavhADAAAAAABgahbEAAAAAAAAmJoFMQAAAAAAAKYmQwwAAIAnp9Op3W631lpdvlUkc2spPyGTC9YTyXTaK4sosn1lJtqoHK1MLlgmqyyiKnMrcj2W+t0r+6dyDJG2vm6buQeW+untF5nzSzmGW+9LyBr1TM3cL/DTjfr7kjhfiAEAAAAAADA1C2IAAAAAAABMTclEAAAAujJlXarK+UVLN0X2XdvOo8pyfqNK5W3t81e/9+ZBpO3IvkvHvrUcXmYejyxtVNVWZemyyHXfqzRjZm5G+tmrTGVvP2XqiMq8DzMlYLfehzCjynuLWr4QAwAAAAAAYGoWxAAAAAAAAJiaBTEAAAAAAACmJkMMAACAJ9frtR0Oh9baP2cbjMw5iLS9V1ZRb9vHdjO5RlXHs9Ru5FxUbbu0b5XKvI7ItYz0GzlPmXYr51Mv023pXETG1MtpW2orkjv33X5r+omInIuqe+JrO8fjsZ3P55J2mUfls0KeHaw3Mn+UGF+IAQAAAAAAMDULYgAAAAAAAEzNghgAAAAAAABTkyEGAABAVyRDaG07v2orY1R2ydZsokeVxxq5BpkcsKqcsEjO1NK+mWOPZLr12srkaD2qvJ96v1Vlbj16l7yurblt75ifFMl0g6jK+8xcpLW6v2/Mp++592r5QgwAAAAAAICpWRADAAAAAABgakomAgAA0JUpYbe1JFm03cgYt7b76F1K2GwtBTiyhGWv3UxJuEgZxEi/lfM4sl9VWcdMPyNFSjNu3XbN9ltVlYut3Lfna7v3+71dLpfV+zKPqufku7zjeC+ZOWO+/Vrmb+ufdJ6q+EIMAAAAAACAqVkQAwAAAAAAYGoWxAAAAAAAAJiaDDEAAACenE6ndrvdFrerzOvauu2jSPZCZa5Rb9+IV+VoZPLUHvfN5HX17JV9tVc/j7bmp2X62ZpR9at9q8YYnfNV823Uc6cyb2xtruHxeGzn83l1v8wjkn3Z+102UY1oJuK7icyZkfu++3nKyJwn4nwhBgAAAAAAwNQsiAEAAAAAADA1C2IAAAAAAABMTYYYAAAAT67XazscDq21cRkPlZliozKrImOIjD/Sz5Je5kvv2CuvXWRMmX3fUSajaqmtrSrPWyRvLDPfRmUKbs3gihp1PLJjqJR5lzLGp73zHlVlrf60fLvIsUeycGc4N6P5QgwAAAAAAICpWRADAAAAAABgahbEAAAAAAAAmJoMMQAAAJ6cTqd2u90Wt8tk+yy1tUe/0dyGtb9l+hmZGfSK6xPNt4ic88jxVOWPjMyD25opUjln9spAy+QAvsLIXJaqrDV+Jpl062Te95X9zCaTn1r1t9vsuVk/bU6N5gsxAAAAAAAApmZBDAAAAAAAgKkpmQgAAMBqS2VqImXPXlVucWvpvOj4t5aey7Q7snRe77xlyhdlyi1VlZqrKsUYbbeqTGjlvZSZb5lzvna/qN4YK8tSVt3TVeU87/d7u1wum9vic+31rPh0lc/nTLvvfg0ixxM9lq3v8NlKJGbmzGznYg++EAMAAAAAAGBqFsQAAAAAAACYmgUxAAAAAAAApiZDDAAAgCfX67UdDofWWj8TIpIXkcmhWMpIiOQpZDKRIl6RQzUy22trRlImGyOS2zQyR+MVGR2Z/LSltqry1DL3WmS/V+XDVVnqc0SO3vF4bOfzeeUImUnmnT1bNlHV3y9VfS79Xvk3SJXM8YzqN5Pp+Imci1q+EAMAAAAAAGBqFsQAAAAAAACYmgUxAAAAAAAApiZDDAAAgK5I5lZPJm9sVL+P/VTmDVXlJ0WMzIfZmnM08hxWzqm19jqezLFn8voy2/bGUZlzUnXsS+1GsvAymYhVuU1f973f7+1yuWxuizm9Q6ZetN9MZmDkGVQlc44/Me9p6/uk186v9t267Seq/NuHZ74QAwAAAAAAYGoWxAAAAAAAAJiakokAAAA8OZ1O7Xa7LW6XKXMWKV02qp+lfSNGlToaNaal8xQpSVRZvmfrmDNl6SKlJh9VlovcWmIwWqZqaymqyvOW2S8zr3vHM6rUWabdkWVQmVNvvo0qczzSq8a89d7KjKGytOxeXjFPou+edzxvo/ykY93KF2IAAAAAAABMzYIYAAAAAAAAU7MgBgAAAAAAwNRkiAEAANC1VzZWZRbT2n2XMoR67VTmNGXOcSTnZOu2a7bv7dtrJ5PTtFd+0qi5WZU7lTn2pX2r2q7Macuci6qMuohMLlDmfuFniDyrl+bT1lzGzL6RdjP7Zo5nyV7H8wqV57iq38pr94lmP77RfCEGAAAAAADA1CyIAQAAAAAAMDULYgAAAAAAAExNhhgAAACrjczCeUUGR+XxRH4flVVU2c87ZJlkskt6bWWu+6vyYCozQ6qyfrb2+at+e9tm2oqoHONez6R3yBhif1vnUCbLb9S+le1WZV0+/j6q3b3GlDHyeLYeQ2UmJT+PL8QAAAAAAACYmgUxAAAAAAAApqZkIgAAAKtFyoQtbZ8p/RXZNlMCLjOmSFm6nswYM6WBMuc40m5VCbhMubtM2blMuajI9pVlnraWJ41eu63ntfJ6VJZBHGWvuco8tpbZzcynvfZ9l5J2W9/he52nzN8VlUaVfF7b5579MgdfiAEAAAAAADA1C2IAAAAAAABMzYIYAAAAAAAAU5MhBgAAQFdVns/IzIfKTK7vtq3KR/qVTA5V73h6GSOR7KioqkykyJiq8t6Wfh+Zl1Z1fKOyZKLtRrJltm67ZvvevmvbiV67rdeg6hwfj8d2Pp83jYF5jczni3hF3lX0ufEOGYhV+44836/IDXuXecxn8oUYAAAAAAAAU7MgBgAAAAAAwNQsiAEAAAAAADA1GWIAAAA8uV6v7XA4PP37XrkNlblZvX0zWRLvktMU0ct06237uP3IzJCtmWhL7eyV67I1b6yqz+j2e+WpRObbXpkumdy5jN7cjPbz3Xm73+/tcrlsHCGfpPI5E3k2fHoWU+R4Ms/jzHmqvD5bZc7TO4zp8ffZ5vGj2Y5nBF+IAQAAAAAAMDULYgAAAAAAAEzNghgAAAAAAABTkyEGAADAatHcqa3bRjMPtvY7MutjazuRdh/brty2N46q87LU1qjrEd02MsbI/MoYlQuSGWPV8Vbm0FRl/1Q++zLXbu1cPB6P7Xw+b+6Hn6k3j0dmR76jredi5HmqGlPl++Md50xmTHs9yzOq/ub9qXwhBgAAAAAAwNQsiAEAAAAAADA1JRMBAADoypQYHFVGMKI3psrSbJF9I+WLIqXaKsu69cYROW+VZSnXju9XMudp6zzea35lzvHIsqFb50xGZXnFreVIl9qN7BvxtZ37/d4ul0tJu8wr886OlIv9BJl7uNdOZRnaXruZkq+jxhTZN6NyHmf6ecd5/45jejVfiAEAAAAAADA1C2IAAAAAAABMzYIYAAAAAAAAU5MhBgAAQFckR6cqByiTQxHJeKga75LKXKZHW3OOomPK5GxslRlTZc7JqHmdyY7L9Nvrp/LYIxk2vbkZOW+V91okpy2iMlvx03Ob+Bx7vS8/3V65UpFn6l7e5RkUyY6NeNU9sNd7+KfwhRgAAAAAAABTsyAGAAAAAADA1CyIAQAAAAAAMDUZYgAAAKyWydHZKxNpKb9j675LuUaPMhlWkXa35hxFstaW9n2UyUjaOqciuVOZLLx3yZLKZOFlsuS2bpuROeeZ/LReO5XzIHMev7uWx+Oxnc/nze3yM0WeDXtlZb2DmY/tJ9h6/d5ljmf+7jN3n/lCDAAAAAAAgKlZEAMAAAAAAGBqSiYCAADQlSmL+HXfypJjVftWlgnMiLRdVQ4nU+JtL1Wl5JZ+qzwXo0qO7VXacFSZzUiJ0Ui5yMe2q0qkLtnreRUZx9d+7vd7u1wuZf1ApNSsMm3M4F3n8av+Vp2FL8QAAAAAAACYmgUxAAAAAAAApmZBDAAAAAAAgKnJEAMAAODJ6XRqt9uttdbPRHpUlY0Tzdypyl6qzESK2JqBtCSStdQbU+S3x7ajcyay73f7Lcnk2UX2jV7L3r6jruXILL9eLtjjtpl+t86ZSplMt6rspa/tHI/Hdj6fN7cFj/bK5wMYyRdiAAAAAAAATM2CGAAAAAAAAFOzIAYAAAAAAMDUZIgBAADw5Hq9tsPhEN4vk6dUlZ0TEcleWsp0WtvOr/b9+vurMpAy1yMy5sqMqqp2M/329o1ey6osrMx82ytPbVT2WkRlXl8kCy+Tgbb2OXS/39vlcum2BT2V9wfAu/CFGAAAAAAAAFOzIAYAAAAAAMDULIgBAAAAAAAwNRliAAAAPDmdTu12u7XW/jlHZK9MkUy7S2Pc2nY022drW9Gstarjqcx4ixxfpJ1R+1bmj0X27c3NURluS2Pq/V6ZiZY5j2v7/FW/a397bCt6/2cy0Xqq7i2Iirwj5I0B78oXYgAAAAAAAEzNghgAAAAAAABTUzIRAACAJ9frtR0Oh9bauPJ30d+39pMpIxYpYdfrt7KMW1V5v8qybZlycXuVzuuNYWnfyPXqXfdIP5UlOiP7Zsp59tqqOqdL48rMzaVxRMbUayczpqW2YZS93uEAI/lCDAAAAAAAgKlZEAMAAAAAAGBq2ZKJv22ttd///vcFQwHgFb48w3/7ynF8gN+21trf/va3V48DgI2+PMO98/qe3nnH4/Hbje/3e7ex3r4Zj/1mxtjztd3osfa2H3VeHkXGsHROv/4e2TfSbnSMWy2d/1HHkx3Xd/1Ej2frGKLHHjlvPZnxR8a01FZkTFvPY3T+fLevd95q/jtvpcx7q+rZDfAo+r77Tabm9J/+9Kc/tNb+c3MDALyTf/vjH//451cP4l155wFMxTuvwzsPYCreeR3eeQDTWPW+y34h9v9aa//WWvtLa+3vybYAeI3fttb+tf3XM53veecBfD7vvHW88wA+n3feOt55AJ8t9L5LfSEGAAAAAAAA7+5fXj0AAAAAAAAAGMmCGAAAAAAAAFOzIAYAAAAAAMDULIgBAAAAAAAwNQtiAAAAAAAATM2CGAAAAAAAAFOzIAYAAAAAAMDULIgBAAAAAAAwNQtiAAAAAAAATM2CGAAAAAAAAFOzIAYAAAAAAMDULIgBAAAAAAAwNQtiAAAAAAAATO3/AwEmArtdezo1AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 2200x550 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "all_up = np.ones([100, 100])\n",
    "all_down = -np.ones([100, 100])\n",
    "random = np.random.choice([-1, 1], size=(100, 100))\n",
    "\n",
    "from matplotlib.image import imread\n",
    "\n",
    "custom = (\n",
    "    1 - 2 * imread(\"data/test_state.png\")[:, :, 0]\n",
    ")  # load a 100x100 png, take the red channel, remap 0,1 to -1,1\n",
    "\n",
    "states = [all_up, all_down, random, custom]\n",
    "\n",
    "f, axes = plt.subplots(ncols=4, figsize=(20, 5))\n",
    "for ax, state in zip(axes, states):\n",
    "    show_state(state, ax=ax)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2ed02839-81c0-45ff-87c9-a429ea8ccf51",
   "metadata": {},
   "source": [
    "If you stare at the first two long enough you'll realise we can figure out the energy of all_up and all_down without writing any code at all:\n",
    "\n",
    "<div style=\"margin:auto;max-width:800px;\">\n",
    "<img src=\"../_static/energy_diagram.png\"/>\n",
    "</div>\n",
    "\n",
    "So we know that for the first two:\n",
    "$$E = \\frac{1}{L^2} (4(L-2)^2 + 12(L-2) + 8)$$\n",
    "\n",
    "And for the random case we can make a pretty good guess that it should be zero on average. And the last we will just use to as a testcase."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "c2e63baa-8d97-41f3-9070-013aa6bf83bd",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "For L = 100, We predict E = -39600\n"
     ]
    }
   ],
   "source": [
    "def E_prediction_all_the_same(L):\n",
    "    return -(4 * (L - 2) ** 2 + 12 * (L - 2) + 8)\n",
    "\n",
    "\n",
    "L = 100\n",
    "print(f\"For L = {L}, We predict E = {E_prediction_all_the_same(L)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8351475b-73cb-4cd9-9fe8-3b60af0b5a6a",
   "metadata": {},
   "source": [
    "## Exercise 1: Write a function to compute the energy of a state\n",
    "\n",
    "See if you can write a function that calculates the energy and reproduces all the above cases. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "36ac021b-1da4-48d1-9cad-6ae3ba1e7e61",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2UAAAOCCAYAAAD3N3ppAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABDrAAAQ6wFQlOh8AABd4UlEQVR4nO3de7gtSVkY/LeYgdkoGBw3CglHJt4hYvRzIALmcRBFHMCggIAGHR+N8UNFjWNioiBeE/GWfOJ4N4MKbDFykauRQaJuIzCACOJdRjcgwmGQ+x5u/f3RvZl1+qzTa9Wu6lVrn/P7Pc9+Zs7q7urqy1q13tX1VqWu6wIAAIA2btG6AgAAABcyQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANHRxycbXXXfdJRFxt4h4S0R8sEqNANgmF0XE7SPi1fe9731val2Z84k2FOC8t3YbWhSURd+YvKywDAC2390j4vrWlTjPaEMBLgwr29DSoOwtERGPfvSj48YbbywsCoBtc+mll8Y111wTMXzeU9VbIiIuv/zyuOSSS1rXBYDKbrrpprj++usj1mhDS4OyD0ZE3HjjjXH69OnCogDYYrrX1ffBiIhLLrkkdnZ2WtcFgPmsbEMN9AEAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADV1cuP1FERGXXnpphaoAsG0WPt8valmP89RFERE33XRT63oAMIOFz/eVbWhpUHb7iIhrrrmmsBgAttztI+JvW1fiPHP7iIjrr7++dT0AmNfKNrQ0KHt1RNw9It4SER8sLAuA7XNR9I3Jq1tX5DykDQU4v63dhqau6+avDgAAAEsZ6AMAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANBQlaAspfT4lFJ3jr/vrrGPOaSULhvq/k8rlnn5cNxXnGP5Zw/L/+05lt8ipfTGlNJT1tzfFUN5lx+/1vNIKd05pfTslNLrU0qHw3H9ekrpU5as+8CU0itSSjellA5SSt+bUrpoyXq3Tin9UErpb4d1b0gpfc+S9a4elh2mlF627HqklG6bUvrZlNJbU0rvSin9ZkrpzpnHuJtSOp1S+qUlyz5iqMOzc8rcRsNxPHEL6nGHlNKvpZTekVL6x5TSL6eULl1z2ytTSq8c7om/Sik9eu76wjq0oWeUqQ0daEO1oTPUQxu6xWo+KXtvRNxzyd//rLiP2i6LiO+JiGoNyipd1708Iv48Ih55jlXuExF3jIgnb6pOM7pNRPx9RPzHiPiiiPj2iPjUiPidlNLu0Uoppc+JiGdFxJ9GxJdExI9HxHdExA8vFjY0MM+JiIdExOMi4n4R8diI+MBovasj4oci4okRcWVE/FVEPD+ldLdR/Z4aEQ+KiG+KiIdHxD+LiBemlG697gF2XXc6Iq6OiKtSSv96tPhxEbEbEd+4bnmcW0rp4oh4QUTcLSIeFRFfFxGfGxHPSimlFdveM/p77BUR8cURcW1E/GRK6evmrDNk0IauQRuqDeV4tKEnQNd1xX8R8fiIeFeNsjb5FxFXREQXEZdXLPPyocwrJtZ5XES8PyJ2lyz7xYh4S0TcstUxzHzOP3mo71csvPaCiHj5aL2rI+J9EfFxC699fUS8bfG1JeVfEhH/GBFPWHjtooh4bUTsLbz2r4Z6XLnw2scP1+UbjnFcL4qI1xxdt4i461D//zDz+UwRcckGrtsNEfHExvfOw4dr9i8WXrvX8Nr9V2z7/Ih4yei1n4uIN0bELVoelz9/2tAzytSGTtdXG1r3fGpDtaFb87exnLKU0oOHLgIPXHjtdsNj9l9feK1LKX1nSukJKaW3pJTemVK6NqV021F5t0spXZNS+vvhEfzLU0r3W7LfB6SU9lNK70kpvS2l9OKU0mcNj+J/Z1jtZUddRY5R/nenlN40PLp/ekR87Bqn48kRcXFEPHRU1iUR8WUR8etd171/qPtvp5TePDxqfklK6f5TBae+O0mXUhqX/cSU0g2j1+6UUvrVoevAe1NKv5tS+uw16l/ircN/b7nw2mdFxG+N1nvBsM4XLbz2tRHxtK7r/mGi/HtFxD+J/he8iIjouu6DEfFrEXHlwq9BV0bf8Dx/Yb2/i4jfj4gHrHksi74hIj4p+l8yIyJ+OvoG5n9EfPj+P3rs/6aU0k+llG5ztHFK6arhuu0uFppSek1K6dqFf187vHZlSulVEXFT9L+Mxmi7o/vgIUuWvTyl9LTh/++YUvqllNLfDPfAX6a+a8slUwc7vI+eM3rtrG5HqXd1SukvhvfR36SUvm2q7AlXRsQfd133J0cvdF33B9E3due8ZsOxfH5E7I0WPTn6X9Q/65j1gY3Rhp5BG6oN1Ybm04ZuuapBWUrp4iV/KSKi67pnRv+48xcW3jTXRP8LzDeMivrmiLhLRHx1RHxn9I/af35hP7eKiN+OiAdGxHdF/4Z6bUQ8Ny08Xk8pPTwinh0Rb46Ir4iIr4yI/egfsb8ibn4k/jVxc1eRnPK/KSK+PyJ+Zajj66L/5WBS13V/HREvGeq06MqIuF3c3O3inw/1f9RQ/n5EPG/xDXtcKaWPjv7D8zOjP98PiYh3R8SLUkofu7DeUa7DZQX7ukVK6ZZDGU+MiIOIeObCKjvR/yK26Kbhv3cZyrhVRPw/EXGQUvqV4QvCO1JKT0kpfczCdncZ/vtno/JeGxG3jf7aH633593wc89ovbtEpq7r/iL67h6PTX3//M+NiK/vuu6DKaUviYinR8RfRMSXRn/PPCrOPAc5/mn0DdWPR8T9I+KPltTnhoj4vzHq4pNS+uToz+NRg7sbETdGxH8YynpC9O+7nz5m3cb+R0R8X0Q8KfoP/Wsj4odTSh9+zy80fo9fUdZdou+eM7bqmn1iRNxqybavXSgXmtOGakPPsS9tqDZUG3ohqPG4LfquF905/q5YWO+joo/Inx43P0b94lFZXUT8TURctPDa10bEhyLi04Z/f030j8jvOtr2JdH/ChTRP5I+iIgXTNT7iljSbWHN8i+KiDdExC+P1nnK+LjPse/HDMd0auG1p0XfKKUl698i+l8GfysinnKuY4i+j38XEQ8dbf/EiLhh4d/fG/2vXB+78Nolwzlb7LZwdG0vK7g/fnXhfviriPiU0fKXRcTzRq89alj/Z4d/33H49zujb2TvF/0H35sj4rcXtvuuiDhcUocvGLb/jOHfv73s3oiIH4iIG495nEcfWl1E/I+F118RZz/2f+TifRIRVw3/3h2t95qIuHbh39cO691jjfp8c/R5KrddeO1xw3Vf2l1juMe+Yrj/P2Lh9RtioetFRLw4Ip4z2vaMbkfRf5B/KPqGdXG9J0SfJ3GL0T37+BXH85cR8TPnuL9eMbHdvYfyP2fJsXYR8Zjj3tv+/NX4C23o4jra0LPrrg3Vhi6upw09T/9qD/Rx9yV/Lz9aoeu6d0T/IfBvoo/6f7rruuefXVQ8u+sflx95evQNxD2Gf98vIl4dEX+x+ItiRFw37DOiT4a9U0ScNaLPGtYp/07R/9ryjNG2/2vNfexF/2Z7eEQ/ilH0vyo+tRvu9tR3jXhSSukN0Sfivn+o21kjLx3D/aLvenLjwvF9MCJ+L24+xui67vFd16Wu/9VoqeFXvMVfdsf31WOjv3YPjYg3RZ8I/PELy38qIr44pfQtKaVLU0qfGxE/ONTnQ8M6R2X+Y0Q8pOu6/9113ZOi/4X4C1JK91gor1tWzSXLzrXestdX6rrufRHxE8M/fzgiIvXdKz4z+i8Li349+ms6Tmxex+mu6166xnpPi777yoMXXnt4RDy967qbhvqllNK3ppRem1J6b/T32FHXoE84Rt0WfcHw399Y8j66Q0Sciuh/kRzuscevUWbJNTvXOse63lCZNrSnDdWGRmhDI7ShF5yaQdmHuq67fsnfO0frHfVfvST6X56WefPiP7que1v0N/odh5d2o+/D+v7R33+O4SaNiKPH8W88xrGsU/5RXd482naqr/aHdV335oh4Ydzc/eJLI+LWMXS7GD6UfzP6R/iPi35EqbtH3397J/eAltiN/oNmfIyPjJuPcV1HSddHf49bXNh13eu6rntZ13W/EX3/9ltFP5rUkSdF343gR6PvL39dRPxM9F0C3jSs87bhv/vDB/eRFw3//RcL6+2klMbn6Hajct4WER+95Fhut7DOcbxv9N/bRf+B96bFlbqu+0D0x7rWULQj43tuqa7PG3hRDN0vUkr/Mvrk6acurPatEfFj0Y+q9G+ib/iPuiSV3me70R/76Tjz/njBsDz3PjvuNTtaNt72o0fLoSVtaE8bqg1d/O/tQhuqDb1AXNxgn4+PPpH3LyPimpTSfY5+1VpwRqLv0Hf7ltE/ro3oP2j+OPouGedylAx7nKF61yn/qC7jpOSPy9jPUyLiSSmlT42+YXlVd3MC5idF36g9uOu6Zx1tkFYPNXs4/PdWo9fHH1w3Rv/GfuySMm5a8tqUn4t+mN0j52zEu657d0rpz6I/vqPXuoj49pTS90XEnSPi76K/3j8YEX84rPOelNLfTtTh6NfAoz7Pd4mIVy4sv2v03TbesLDeF6aU0uj+u2ss73N9XP8Y/a9IZ9wXw69dHxP9dYhY/7pF5P0q9dSI+NkhZ+AR0X/hedHC8odFxG92XfefF+p21zXKPVyjrjcOdf3cODvfIaIf1jrHn0b/i+nYXePM+2/sr4f93yVubsyOtjsqF06Kx4c29Ig2NLShw0va0NW0oVtuY6MvRnx4noP/FP3oOo+MfpSfZaPIPCidOenhl0V/Y75s+PcLo38s/MZlvywO6/x5RLw++r7t53J0k49/zVin/NdH36h86Wjbh8b6nhF9l5VviYj7Rt/AHDlqOD78Rkz9pIz3XlHmm+PmN8/RdpfE2Y/4XxjDh+eSY3x1xjFE13Xj83TOBiWldLuI+PTocx7G5by967o/7rruH6Pvy33DUM8jz4mIz01njmp09Hj/VcN//yAi3h5Dl5ZhnxdFxJdH3+f+6MP4edH/OvRFC+udiv7D77krDnltXde9K/ok4i8fLXpI9D+K/N7w79cP/128bp8efReFEk+PvrF9aPTn5Gmjbk23jrM/7L9yjXJfHxGfmtIZc5t84Wid64b/fsyaTwBWeV5E3C2ltHiOPif6/vTnvGZDN5MXxdnX4JHRv4dfedZGsIW0oWfRhp5ZjjZUGzpFG7rtugqJadH/cveeiPicJX+fOKzzkdH/sve8he0eG/2vBYtzJnTR/xLz7OgnqHt09L/OPG1hnUsi4vroR+L5+ugTdR8cfeLtf11Y7+HRv5l+I/rHyvcf1nngsHw3+j7JvxT9qFGXZ5b/mKG+PxJ9//Ifi/6N9uFEzTXO3VG/+HHC8lHC8Kui7yf/8OhHQ3pdRLxmYb0rYpRoPZT5tuhzDx4QfWLz38WZScofE/0H9vXRJwR/XvQfOj8SEd+2sN7jhnN052PeFz8xlPt5Q31eGf0H/qcsrHeP6Ce6/MLoR+n6heh/afz8UXl3jv5Xs+cN98bXRj8fzTNG6x3Nz/Lt0XdZeXL0DffdRus9J/p77RHRj9p1ffT36K2XnN+r1jzmq2KUbDwc04ei/8Xt/tHf02+PiBcurHPL4Rr98XDNHjlc+9NxdpLya9apy8I2Tx+Os4uIe46WPWE4198U/T38pOgb+/E9dUOcmaR8/2GdJ0bfqH9P9L+mnXHvR8RPDtfsu4b1vjj6L1DPHF3XD0TE41Ycx8XRN85/MpzTLxv2+XuxkNi/7J6N/v39/uhHoLtiqM8HI+LranwG+vNX8hfaUG3oue8Lbag29B9DG3pB/NUpZHrkqGuHdX4m+u4Qd1zY7qLoH62/PG6eMLCLfgjfHxvWf2dE/HJEfNRonx8VfR/qv43+w+ON0Uf6Dxit96BhH++N/kP2uoj4zIXl/364Kd8fQy+AdcuPvq/v46J/nP3u6PsUP2D8plpx7h40rP9/liy7e0S8dKj7X0TEV8XoAyWWNyi3j/4XxLdH38A9JkYjRw3r3SH6D+83Rv+hchB98uy9llzby45xX3xJ9CMMnY7+i8NfR8T/jIhPGK33mcM1eufw98IYffAtrPvZEfF/hnNyOiJ+NhZGRlq4Lt8xXLvD6H8dvs+Ssj4q+m4jN0bEu6LPP7jzaJ2j6zk5seLC+lfF8hGgvjT6xvSm4X75qYi4zZJje2n0X87+OPp5QZaNHJXboDxsqNPrliy7zXBNbhz+fi76LzCTDcrw2ndE3wi+K/o5bO4XZzcoKfrG6tXDsd8Y/TDDi19aLos1Ro4a1r1j9MnX7xju71+J/lfEZZ9Hl41evzL6Bumm4V78xtx72p+/Of5CG6oNXX5s2tCbX9eGakPP+780nOitkfrJJ7+j67ofbV0XSCl9f/SNwd26bXuzAIxoQ9km2lBY30ZzyuAEundE/JDGBACyaUNhTS1GX4QTo+u6z29dBwA4ibShsL6t674IAABwIdF9EQAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQkKAMAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANDQxSUbX3fddZdExN0i4i0R8cEqNQJgm1wUEbePiFff9773val1Zc4n2lCA897abWhRUBZ9Y/KywjIA2H53j4jrW1fiPKMNBbgwrGxDS4Oyt0REXH755XHJJZdERMSpU6fOWOHg4ODD/z9eNra47qr1x+vmmKpjqVVlLy6vud9cU/XIOYZly9fdb66S/eTcW7lK7tua93jN67Kp923JdVp1fLWuwypzvo9r3sc5n0WLyy699NK45pprIobPe6o6qw0F4Pxx0003xfXXXx+xRhuauq479o6uu+66O0XEwb3vfe/Y2dnpC0zpjHUWyx8vGxvXZWr9knpP1bHUqrIXl9fcb66peuQcw7Ll6+43V8l+cu6tXCX3bc17vOZ12dT7tuQ6rTq+WtdhlTnfxzXv45zPosVlu7u7sbe3FxFx6r73ve/r194hKy1rQwE4fxweHsb+/n7EGm1o6ZOyiOh/VT19+nRElH1pyPliWPPLSkkgtWo/Nb/s5pjr+HP3VXIMJV+6a1p1PnKCoTnv26nlNa/pnGret5u6DrlKfhDJsWrbnM8iAGBeRl8EAABoSFAGAADQkKAMAACgoSo5ZQcHB2slKc+Z11JzwICcPLDcvKe5cjc2dfyrltccFGKsZFCDHKXnoyTHLicPatW2i+vPeX/MmRtZUk7Neyvn3K7atta64/VLch03mVMHOeYcmAtgm3hSBgAA0JCgDAAAoKHZh8RfVNr1b0rJ8M81h8Sfs6tSiZJ6zdkNqqS7XklZJd3RSuo05zx2Nefaqtk9b665+WreHyVltxpOf7zvkvtFlzC2RaspWABa86QMAACgIUEZAABAQ4IyAACAhqrklC0qGVo7R6t8tFXLcvKvSvLRSoY4z91XzSHxxzZ1v6zab8m0B2M556OkHjnXODePsmauX4ma12XRpobxr12PuY7ZEPjba67PyDm3nfP9U9LeAGwzT8oAAAAaEpQBAAA0JCgDAABoqHpOWYmSPu6bmmssNwdkzryxknqUzFE019xapTlktfILcq9hzXutZl7Hpubz2VReR+m5nJrTqyT3b1U9p5aXHMN4+Vx5pYeHh7G/v7/WuhzP4lyfOUree3NuO2eOmTwy4HzlSRkAAEBDgjIAAICGBGUAAAANzZ5TViuHqua8Qq3yR8Zlb7Kf/VQ9cvvobyovblv2m6tkbq2See1ycopKcpXGSuaxm/N+mGuOwNyySuaZKjnG3H0dt1zKHRwcxM7OTkTUnbuwVf5iq88TgJPMkzIAAICGBGUAAAANVem+uNj1YkrNbj65+zrufmtuu0k1z+2qodZLupOUdA0dazU89FRZc963q9Sc2mFKyfu65LyXDAFf2j2vVtfQmlMkjNXs+si8pobEL3n/nMRtV9nWNheglCdlAAAADQnKAAAAGhKUAQAANFQlp2yqP3ytPI7SfuQlw+tvS85ZzpDEq+QMlVxz2O1tGd641TWtOd3CnGXPdS225b2V+xmwqVyv0mHLpxy3zru7u7G3t3esbVnP1JD4U7YlL3CT72ND4gPnK0/KAAAAGhKUAQAANCQoAwAAaKj6PGUl85OU9I8vyXsaq9lPvSR/JGd+o9x65Gxfcy6lmvlpc5VdM0+jZg5QzXmnWs0PlrMsV04+45yfAXPO8bap+5bNmsrLPgk2NTdhje0BtpUnZQAAAA0JygAAABoSlAEAADRUfZ6yuXIeSua/GptznrJV+yqZty1nbrESq463Zg5Zyf1Rc76nkv1uak6vsZz7eFvPx5zncs561crfm3Muulr3x+HhYezv7+dUk0xTedkncV6unO8BJXmVJ+V8AKzDkzIAAICGBGUAAAANCcoAAAAaqj5PWYmc3IvcHLOpbXPMmQNTMt9T6XxHteqRq2Z+wFz5aSXbls4nl1NWznxg23rf5ijJuat5XXLLrvk+nrrnNzkfI2Wm5inLuY9Pwlx1pTnL59v5ADjiSRkAAEBDgjIAAICGqg+Jn6OkS92qbglzldWyO1rOvkuGQC/p9lQ6xPdUWXN2t5pzKPbj1mNbuuzmyhkOu2TKhG0ZPn9VWXOaa1+Lx7+7uxt7e3uz7Ife1JD4U05K97y5unjrgsu22pb2mJPFkzIAAICGBGUAAAANCcoAAAAa2uiQ+HPmOa3aV0lZOduW5APk5JiVDJ29qqxN9vEv2bZV3k/Oud+WvuNzTpmQU1bJNas5VcOcUwKU1KPkupTsS27OyXQSh4CvOfUJbIuaOf9cuDwpAwAAaEhQBgAA0JCgDAAAoKEqOWVTWs25UnOespy5kmrmtWxqDq+S3LbSfU+VO+f2Jfk1Ne/TqX2VzomXU/acc4DNld+ZkzeZW1bOttv6mbBq+abmgKNMSS7xtuSqlLShOZ/VJ+V8sB1q5of7fG3jfPsM8KQMAACgIUEZAABAQ4IyAACAhmbPKWvlQuuzWzKP23jblnOxzSUnV6fmnHir+jOX9EOvOdfYtrxfauag5pyP0nN93HVX2VT/+KlyDw8PY39/v8p+yFcz73QblORNrlp+Es8Hm1Pz3tpUzvI2aTVGRMkccKuWb9u18KQMAACgIUEZAABAQ1W6L546dSpOnz4dEfM90i3pBjZentt1K2e49Bxzdk2ac1jumkO/1uyuty1dV2rdL7nXsOT9sy3d8abWzx0Cv9UUAXN2wSxx3GPY3d2Nvb29avVgWqv3Zok5u7TnfK7Bojmnlckt6ySo+T7e1HfQ0mu8bTwpAwAAaEhQBgAA0JCgDAAAoKEqOWUHBwexs7OzdFnJsNyL2+b2Ay3pl18iJzelJJ+ktM45w4PnlJXbj3hx/dKhTGudr5pDz+eaM/ev1rY1c0JqmjO/pCTXoCTXb5Vax7wtuUnkmfM6zVV2abnbWq8p257HUsOmvmPl7LdlWSdBzZz3TdbjfL5OnpQBAAA0JCgDAABoSFAGAADQUJWcskU5/YpL5hFaVXbO8py8jjlzPkrmWispa2yT85Rtau6KTV7jnFzIOa9LybpzzVeS+xlQKwdzvH1p3mDJPClzzZe2rbl+lCm5b6eWlcwLucltx07i+diUOfOdc2wyD3tq33Ne45Oo5WdAjpLvwif9OnlSBgAA0JCgDAAAoCFBGQAAQENVcspOnToVp0+fjoiynJgcm5zPp6Rf8VzzY+Xk4qza16o8nxybzEebK18vt6yaSq5LyTwxNeeYmTrXOblaYyXbjrffZH/4OZXktq1b7uHhYezv71crm3pKPiO3Zduabfm2HFPu8rna3G3d75x5QDXvvZOenxTR7p4vUfO730njSRkAAEBDgjIAAICGBGUAAAANVckpOzg4iJ2dnZXrbTIPbGrftfLLxuXWLqvm3AxzzY2Ua87+7yW5f5u611rl65WUtcl5QDb13tyWOYpqnsvzIf+Bs9Wcm67VtnPd57mfF7W2zd2+Zpt6oe03V841Pt/muMrVcj65Ejlz0k5tu408KQMAAGhIUAYAANBQle6L62r5WLFWd6OS/Y61Gjo7omx48Lm606wqK0fu8Ok5XR831bWt9P2S0wVxqgvAJu+PmvfplNL7dGpfc3Ylrtklc8q2d/Ggt6nuaKuclPulVUpEyRQ1JVoNHZ673012899UPabKHWt1zbdlaPk5p0Nate22fXZ5UgYAANCQoAwAAKAhQRkAAEBDVXLKTp06FadPn46Iuvk1OfkTNbct6ZObk7u0Ss0hzqfKrjk8+Cbz80ryekrXX1fJ1AUth0uvmS9Rq6zSaSDmylHM3U/Jey/HtvedZz21PqsvhPuh5vmo2VbPOfT4pj5PpvY7LrtmHnJuPUr2W5JrPrXvbZ0iYFs/A+Z8H9e8xnPwpAwAAKAhQRkAAEBDgjIAAICGquSUHRwcxM7OTkTUnQunZlkl/Y6nTM3vlLvfHCX9m1ftuyRHJtem5o0puaa560/Vo2b+1ZzzuG3j/CUl77Xx+qV9x+c6P7nXvFY9psrZ3d2Nvb29Kvuhrjnn8zmJauas1lRznrKaueZzzXO3rfflJnOaS/a7qXnLzpccs3WX5W7b4nx4UgYAANCQoAwAAKAhQRkAAEBDs89TNqWk/2bNeTDmzKeZK8duznlQSvKgcsuumetXc365qWWbnFuqRM49P6X0vM81X0vNXMdWcxTVnN9ornmXDg8PY39/f+1tmU/Nz8ALUU7O7lz7XbXvbZn3sNW2NZXOjbup/KNWuVzbmkO2Ss05jLctx8yTMgAAgIYEZQAAAA1VHxJ/W21qSPw5u01uapjUmubsjjc+HyXDyZc8ws7pUraqHuvWaZ31S+75km1rdgco6ZI6ltPNp6Q78Ca7R8zVNXRxW0Piby/dFfNs6ny16r7Xsiv9XN9Papa76vxsy1D9c07hU2vbsU1+P92WaZjm4EkZAABAQ4IyAACAhgRlAAAADVXJKVuUk0NVMtRrzjDc6+wrp+ypcucc5r/m8Po5+TU55swZmnPKhJw8sJL8olVKhrGfq591br7V1Par7uk58znn6ktemp+WU9ZYrffEtgxhTZ6a+a0nUW77Mtfw8au0uk6bnK5jrnptcljyuT4jT8q5bHV/zJmnXvOYNsGTMgAAgIYEZQAAAA0JygAAABqqnlM2V/5EaV/PklydnDkzSsqqOc/SKnPOnbQNcvPCFp2U4605L1mJkrzBOfva55gzF3RTc6rUnK9msc6Hh4exv79/7LKZT8nn3PnwOT9W83zUtC3XqeT7yybnodrG6zTnfGBzvhdLjmkbr8Oq5dvyXjsuT8oAAAAaEpQBAAA0JCgDAABoaPZ5yuaaF6RkLoKxln1jj1tWaV/Xknm5psraZH7NXPNN1Jz3YpWa8+dN1WvO/MWSsqbWHdtkf+9N3i85+80597Xmgdzd3Y29vb3JsmhjU58954tW87bN+V1nat1tMWebMde2q7bfZF52SVklOd05+ynR6v4Yb78N7x9PygAAABoSlAEAADQkKAMAAGioSk7ZqVOn4vTp0zWKOqfSecqmyspRs/9qaVlT9dpkvt5x112131Vl5+Qv5uyrdA68nLK39brU7GddK58z9/6o+X4p6Yees+2qesyVozlVjnnKzk/bMCdPbbnvn5OgZn7rSTHnnGDH3Tb3XLbKX5yqR+66repd831cOrfwpnlSBgAA0JCgDAAAoKHqQ+JPPQosGfayZL9jJfvN7UI1tX2r4a/HSrtIlXR1yxl6fWrbnP2s2lfpfkuu06a6M5bsd5PTHNSU042y5vspZ9uSdXPrMbXutnXp4HhKPl/PR3NN0VNSj7GSrtUnpXtmyTFtatuaTup7bVvqvS312ARPygAAABoSlAEAADQkKAMAAGioSk7ZwcFB7OzsFJdTMuz02KaGji7ps1ySY7dq3Zr95VeVtam8hZpDwueo2U+/Zt7XquHUp5bl3B+599qUnDqvUpIHVvr+WXc/q8ou/TzJcdx7bXd3N/b29o69X9qpmUd5EuW8f7b1fOR8Vm1LnVeZs01Zd9k6y2HTPCkDAABoSFAGAADQkKAMAACgoerzlJXI6dOdm2uRk8exbjnr2NS8FzVzUTa5bc35nnLmB8vJLaip5rldpda8dnPmOdXM3dpUvtVYzblwWuU4nJT5jain1ZxNLZ2EYzoJdZzbVPtzIdynXLg8KQMAAGhIUAYAANCQoAwAAKCh2XPKNpXLlFP2nPlXq5Tk+ZTMR5KTu7XJuZJK8sBKzuVYzTnxpuqxyfnSavW1L53rZWrOs5L7YdXyOa9LSS7kXLlvc32uHR4exv7+/trbcjLII2RbuTe5UHlSBgAA0JCgDAAAoCFBGQAAQENVcspOnToVp0+fXrqs1jxUNXM+SuZnyclryq3Hqn3l7HfOOd+m5B5TzjWecy6pnG23Zb9z5tzNdQ/UzAsbmysHddW+NjnHWQnz+Zz/SnIOYVNyvoO5b7mQeFIGAADQkKAMAACgoSrdFw8ODmJnZyci5h22fGrbOYdQzemqNGcXvJxycro3ztmVbZWcLpqthpfflikTctV6f5V2o8xR0mW3pFty6b5y1Or6OKfFOu7u7sbe3l7D2lDDtt5rbM6mvp+U1sOQ+FyoPCkDAABoSFAGAADQkKAMAACgoSo5ZVOmclNy+jfP2dd5ziHhNzUE/Jy5SCXDko/LKhmyeVN5YKvqtaltN5k3mFv21LqbnGIhR875KJluYJWc4Z5r5smNydsAFm1LvmvJdEhwknlSBgAA0JCgDAAAoCFBGQAAQENVcspOnToVp0+fjojpHIhtnSclJ7crN39mk/OpzaUktynnfOTm+WxKq3ncNjVX1qqySu/hufIDcvOrco5x1fKc+Ren8tFKc8ZK5nWbsrjfw8PD2N/fP3ZZQBvb8n0jtx7bUm/YNE/KAAAAGhKUAQAANCQoAwAAaKhKTtnBwUHs7OwsXVaST7KNcxSVzvdUKz+t5nxFY3PmDJXM25ZTj03KmXut5rYl87a1mkus5hx4ucun1s3N3cqZf7FmPtqUkvy0qWPY3d2Nvb29tesBbIdtzXmvOZ8rnE88KQMAAGhIUAYAANBQ0yHxSx5Ll3RlqjmU9pzdCEuUdF1q1YWq9P6o1QUvtzvacfezrOycbcfmuudLp4FYdz/brORzrWTY+k1NxzC1H0Piw8m0ye9JtfY7Xn5S2giowZMyAACAhgRlAAAADQnKAAAAGtrokPg5y1apmdeSo2TI97Gax9AyT26q/3fJcOCr1Jy6IGc/JXJyiEpz23LKmqpn6fHPNaTxpvKt5lQ6xQbAucyZwzxWc/wAeWRcqDwpAwAAaEhQBgAA0JCgDAAAoKHZ5ymbS8k8GHPOw7VKSa5OzXndSnL9NjUHXO78T1M5VDn7zVVyr9Wc02uu/IGaeUwl56NmPucqtea8W1V2aZ3nmk9ucd3d3d3Y29vLqhdwss05jyywnCdlAAAADQnKAAAAGhKUAQAANFQlp2xRzTyfuebUyM29yNl21fpTauYX5ZhzDpGcc12Sf7ZsXyX1yil36pjmvF82lUOVW+dW93zJfHolZeUsy93PKjU/q4ALS057XNKW+yyC9XhSBgAA0JCgDAAAoCFBGQAAQEPVc8pa5UWNTfVhLukrXdMm5/IoyeUqmROuZN3csmrlVM15/Tc119wqc87hVbPskvzFmmrOZVgzJzPHcY/h8PAw9vf3j71fYDuU5M/nLDdPGRyPJ2UAAAANCcoAAAAaqt59caxWt7A5uzqWDI9dc5j/TQ3DPd7Xqi5SNc9Pq6Fwaw6Bv8qmpoEYm6vrX+79UHKPH7fcdcqemqogV841ztGqS/PUtAa7u7uxt7e3sXoBmzFX93BD3sPxeFIGAADQkKAMAACgIUEZAABAQ1Vyyg4ODmJnZ2fpsrn6Fs85dHbNPLCc9UvynjY1jP86dclZNye3rebQ/HMqyaNs9X6pef/MlSeXm8s2tf6cx7vKXLkZNa+xIfGBRfLEYH6elAEAADQkKAMAAGhIUAYAANBQlZyyU6dOxenTpyOi3VxkOflGc/aNLsmLyqlX6THUmj9ulVVl18y/qlVWbr7iXLlsNeflqpn7OLWfVdvmntu55tMrvWZzndvcetW69+SLAEBbnpQBAAA0JCgDAABoSFAGAADQUPV5ykryS0ryrVblEJXkDLWyqdyuk3I+xmrNfxVx5jGXHn+t3KXSXKWpa1yy79z54mrlipbm9uXMW5ZTds080pzPsWX/Xne/q7YF5rOpXHvgZPGkDAAAoCFBGQAAQEOzD4m/qKSbzya72NXsRlmi1fQCYyXnPqfONYeArzmc/tR+1ll/3f2Ol5ccw6p65bwXS7vY5Zjz/TOl9Fwfd9vSYf5z6AYFANvLkzIAAICGBGUAAAANCcoAAAAaajokfo7S4cFLtq05HPZUztCcw1+P1cz9K8mpKjHXtAe516Ekb65kmPbcfU2VvamczZrvl5q5brnXtOZ7cWrbOT/HTuI0IXASbctnwrbY5FQxcJJ4UgYAANCQoAwAAKAhQRkAAEBD1ecpG8uZ/2muPI119r3uvkr7d5fMgZZTj5I54XLza3Ku8brlrFOPnLm1Su+fHDlzjU2ZM7dvrGTetrFa9/iqskpyMOfM11u1r5L7tGReu6my5GnAfC70/M6aucRwPvOkDAAAoCFBGQAAQEOCMgAAgIZmn6ds0bbMTZGb2zbXMeWqlZ83Xt7qGOYsqySnLmc/uVZdl5x8o5Ics5J8rLGa9cjZT66cXLdW75858yhb5lkCN5sz37WVmm3qnHn9sM08KQMAAGhIUAYAANCQoAwAAKCh6vOUlcyHlaOkX3HNOTNycoRq1jM3H2TOfJKS/KMpq+pYMrdWzXpMrZ97f0yZs+98yft0zj7+NT8jcsqdc068nPdxzc+1nG0X7e7uxt7e3tplAee/mnOdTm0vh4wLiSdlAAAADQnKAAAAGqrSfXFRzWG5j7vuOttPaTVcbc4xtuwWt6mhtTfZLewkdOccq9mto+bw+TW7/06Vu0rJUPw1u7/OWc7UuS15vywuOzw8jP39/ZJqAmuq2aW5pm2dRgbOJ56UAQAANCQoAwAAaEhQBgAA0FCVnLKDg4PY2dlZuV7pcPI5ZZXkq00pzb8qUdJPuySXpyT/aq4hzUvNOZxvzpD4NZVcl5rvl5KpCqbqkXtdat4/Naf2KLn3SqYcmco5m1rXkPjQTsln9ZxKpgmRJwbLeVIGAADQkKAMAACgIUEZAABAQ03nKSsxZ+7F2OK2pXkqOXNpTW07VvO816zXJtXKXSrNfSyZe2tT53LOXKWapq5pyzkDpz4T5pynrmZum3nKYDvk5M9vSknOrhwyOB5PygAAABoSlAEAADQkKAMAAGioSk7ZqVOn4vTp0xEx3xxFJeuOl0/N15OrZK61TeaT5JyvTeYMbYuS+3Quc/bLr3mNc957JUry4Mbrl57bmvO4lWybk+8qrwNOnlbtcc35SeWYwXo8KQMAAGhIUAYAANDQVg+JXzIcdsmj9px65JY9VVZJ96vSYdxL1OqiWXO6gVVlz9U1dtn2OUrujyk1p5DIlXNMc3XtW2f7deuRq6RLd8nw+hdCt2M43+UMTb8t3QJzh9PflnpDa56UAQAANCQoAwAAaEhQBgAA0FD1nLKSvsEl+Uc5OWe5+Vgl+UdTy0uGsJ5zCO+aw+uvMtcw/yV92kv7t9caen3OvK+affhL3j81c6Y2mX9Vckw59+nY1Ppz5Rzu7u7G3t7escsGpuW0Ga3ysUpy77f1mGDbeFIGAADQkKAMAACgIUEZAABAQ9VzyqbUnEur5rxLJXlgm5xLqmS/JXMlrVJrHqo5c9lqzkVXc78l9Sjd17pl5c73VWveulxz5iXMNd9gzXrN9fl6eHgY+/v7x6whUOKkzFNWktM9tVx+GRcST8oAAAAaEpQBAAA0JCgDAABoqHpOWat8kpLcrpp5YnPWI2fd3Dyg46473nfJ/GCluTg5czhNzaVVej/UmqdsbJP3ac1520rm1iq5P0ryGXPePyW5XDXz0ebK0TVPGbST8xkxZ17tpr4Hjbc3hxkXEk/KAAAAGhKUAQAANCQoAwAAaKhKTtnBwUHs7Oxkb9dqDquSvKfcsmrmEOWUW5JfkltWzrZTuTq5/dBL8tOmzJlDtWrdkr70OfWec16/HJvMeZgr33WTOYdzzl0IbJ+Sz9dW+y2tc6tjhtY8KQMAAGhIUAYAANDQRofEb7mfqW5yq8re5FCw65Zdc3j0UrWG+C4ZNjdn2TI598cqJVMZ5NjUPVBzeOOxmsPWb+q9N1az+2ZJl8Oa99bitoeHh7G/v3/ssoD1laRqzPmZl/N5W/pZPdcxwbbzpAwAAKAhQRkAAEBDgjIAAICGqueUjdXqG5yba7Gp4dJX5Z/VHKY8p9yaQ+KvUjLE99S6Jbl9q8qeWl46FP1U3/qaw5CXnNva66+rZt7TqrLHpq5xSc5DiZLc2HE9SnLqprbd3d2Nvb29nGoCG1IzH/q4+81Zlrvc1B1cSDwpAwAAaEhQBgAA0JCgDAAAoKEqOWWnTp2K06dPR0TZPENT6885J1HNvK+cPKiS+UhWrZuzr9z+33PNT1KSu7Vq3Rybyidaptacb9ukZI7AVnPwlORE5JQ95zWr9blmnjLYHq0+52vOqVmzbDifeFIGAADQkKAMAACgIUEZAABAQ9XnKZtzLrISOWXXnAtoav3cXKVa8w6t2tec+XtjOTlDNful15xPbqoepXOeTS2reR1q5s3NmaOZs25JHmVNNesx1zHJ6QDm4vMF1uNJGQAAQEOCMgAAgIaqdF88ODiInZ2diKg7pHWrLnY5w7bPOfT61Pql3a9KjmlT3cJyu3PWGiK/tCvfVD1qdtfL2XaTw/pPKZ1+4bjrjtevOd1AyTGVdlmu1aV5W+4PALhQeVIGAADQkKAMAACgIUEZAABAQ1Vyyk6dOhWnT5/O3m5VXsdUDshYzeHSp7atmXuRm180V55caa5fSa5OzjXOqVfJVAVz5ieOzZlDVaus0vzNWtNk5ObU5dxPue+9nGky5sqLG5tzWwBgXp6UAQAANCQoAwAAaEhQBgAA0NDs85Tl5AzVnEenpB5TZZXObzTXvG2r1JrPKLesnH3Nmbu0yXmYauYy5ZRbMr9eSR5UzvLS3NASOfdazXntNpW/V/J5MrXt7u5u7O3tZZUNAOTxpAwAAKAhQRkAAEBDgjIAAICGquSULdrkHE85+6k5V1Atm5zjrCRHaM7zUXOeslpq5i+O1Zxba1U9auUuldZ5rs+ATc2Ptqwete7bmvmbNeux6PDwMPb3949dNgCwmidlAAAADQnKAAAAGhKUAQAANFQlp+zUqVNx+vTpiMib76dmrklOXkvNOYlW1SvHScnzmSor9/hr5pHVmluqdF62ufIXa86tlTvnWc62Oee25H2bK+cal8wPVvPzZJWSef7MU8amlLRdNT9PgHPzfmrPkzIAAICGBGUAAAANVem+eHBwEDs7O8XlbLLbz1RZJV235tpv7r5rdnUr6bI5535KhqKf2ndp95ocm+r6lnvflpzLmu+fud6LNa9h7r7WXVZzPyXbGhKfXDndfWsu1/0Kzi23nfN+2jxPygAAABoSlAEAADQkKAMAAGho9iHxF82Zm1NSVs1h7HPK3uR+58yZyTF1D5QMtT5enpv3tLhtac5DLSVDRY+Xz9k3PKfs3FzHmkPib2p4/RIl7+tN5g3CcZW2TXO99+BCs63fEy9knpQBAAA0JCgDAABoSFAGAADQUPV5ykryS2oqyRkqqXNu3s9x152zL3DJMeTKmR9satvx9qX5WMfd73h5SW5bybbrbL+uOXOV5vwMyLk/crZdZ/vjrpu7n1r5NXLIKDHX+2GO7YFeybxl2ox5eFIGAADQkKAMAACgIUEZAABAQ7PPUzZXnsO43Nz5j6a2zVEyd9Qqc861lpOLknOMuXlNNfPAcvLTNpmTN2Vqv7nXZa57vGQeoU2WlXM+cvazbF9T25bkkZbcp/r4sy3m/GwuaW+Am7X8nsRynpQBAAA0JCgDAABoqEr3xUUlw1CXdNXJ6QZUOhz2unVcpqS73rp1Ok69cpQM1Z+zbcny3Gs8tZ+a12nO7gCtumCukjOMbs0uu3N1OczZzzrLN7VtyfQLkCOnW3bpVCg5ZcOFrOZ7jXl4UgYAANCQoAwAAKAhQRkAAEBDVXLKDg4OYmdnJyJW53bVklturVyV0v0ull1zKO1VSnJkcnJ3cvN8aubulKiVj7aqrJIcs5bTHuRsm1OPmsP6l5Q15/HPNS3IuOy58hN3d3djb2+vWtlc2Erzag2JD3UY8n77eFIGAADQkKAMAACgIUEZAABAQ9XnKas5F05OTsxxy62tZh/dOeewKpkvLGcOuJq5bCW5bSXHX3JMuWXNmfuXs26t+fRW2Zbcrdx8tZL9bionpnTeP5hDabtmfj2mnJTviiedfM55eFIGAADQkKAMAACgIUEZAABAQ9VzykryjUryS3LKLqnHqnVz5rhq2Qe3Vq7fuKyacvss18oTK5mnbVVZOUrnFqs1B1zLuUzmmmttk8dUM892U/P8LZZ7eHgY+/v7a28LOfd8bnssd4VFNb9z1vysPh/UfB+zHk/KAAAAGhKUAQAANCQoAwAAaKhKTtmpU6fi9OnTK9fLzeOYK59krCQfrWY9SvrgtuzPO3Wdcs5XSe7WKpvMg9rkvhaV5NSV7GdTuaC5as7bVnOeu5L91JzHbqpei8t2d3djb29v7f3CuvfWOsvlptRTMh/jtsi9f+ba9iScq1Il55Lj8aQMAACgIUEZAABAQ1W6Lx4cHMTOzk5EzPf4t6Tr43j9ObvJ5dSjZleCkqG1NzlE8dS6pV23TuLwtSXd4laVVaLm+cvphjznlBolUxUct+tfrjmPSXcTtoGhtLfXSWxDx1p9zp/U87VozumhzofzswmelAEAADQkKAMAAGhIUAYAANDQVg+Jn1NWzvol+y3NA6s1jHvNc1lSj1Vq5g2WnNs5hynPMWeeYM3+4CU2lReXe122Jadqrr71m5xuABZtKufyfJHz2TzXfs5XJee25nfD8z2H6kK8tzbBkzIAAICGBGUAAAANCcoAAAAaqpJTNmVq3ouSeac2mY9WMndHybxLYyU5Uznzg5UoyWspPbdT+1m1vOReq5mvllPunPl6m+pbX1KPTeaQ5byPx2oeU859Wyv38fDwMPb3949dFheemp+vF4KS7wXrlrtO2SfhWuQcU815yi7E+fWcj83zpAwAAKAhQRkAAEBDgjIAAICGquSUHRwcxM7OTkTUnZ8kp591zTyWqXqUljWlVi7SsnpM1bNm3+CSc1t6jReX556fTc0dtUm15m2bM8eu5rx+OfXYZN5gTSWfr8e9H3Z3d2Nvb2/NGkLZfXo+5qbkHvM27Hdbr8NJPJfbcu5qKr0O5+M5qcGTMgAAgIYEZQAAAA0JygAAABqqPk9ZzfmwctavOWdRTt/glvOn1dp2zj7+Jf2K5zyXNe+tHJs8ppJzW5LrlzNvTEk9ptZdpWZO5qqyS85HyZw85iljG52E/MxlSnJlWx1zyX5Pas5PyedrTllT666z/km0qfluL2SelAEAADQkKAMAAGioSvfFU6dOxenTp1eul/s4fGoI66l1c/c1Z/e8KXN2D2jVTS73EX7NR94lUxfkDJeec4yb7L5Zc5j7uYaP3+T9UTJEccn9MmeXoZKyT+Lw15xMJdOVjG1Lt6iaU9aU7HdT309O6pDm29gd1ufrhXnMx+FJGQAAQEOCMgAAgIYEZQAAAA1tdEj8Erl9cjeVi5Jbj1q5OqXnea58vU0On1+Sq7OpPt5z3qdzDuc7Vz1Kyi7NP6t5PmrlL+YqyZO7EIZspo2cz9PSnN2S6TpqbTtevsljmrLJY2plG6YJWafsnO9Y54ML4Rg3wZMyAACAhgRlAAAADQnKAAAAGqqeUzZWKwdiKjdrWVk1+5Kvu2xZWTnLc/KP5uzvXXJ+tqUfes25xkrn6TpuPWqey5r9vVv1ra95Hea8T0vurbFNfVZtS74IJ0NJfmdpbuhcZdfctmRO1pplb/J8bCqnaM5jKjmGVvvl/OJJGQAAQEOCMgAAgIYEZQAAAA3NnlNWa06r0hyQWnNa1Zx3qnROkSkl9ZyrX3Xpvuac/6lkDqu55p2qeY/X7LO+yXzOuc5t7rY188Tmssn7lgvbnHNG5mw/57bbkudT8r1gk9dpW+blapUr6/OVGjwpAwAAaEhQBgAA0FD17otzdRmac8jqmuvP2R0rp5w5h22fa8jzmt2vcrt+ztXtdNW5q7WfVWXP2a10znLn6lZaswvVWMlQ9CX3y9R+ctddrNfu7m7s7e0dux6wqLQNnGuKjU2qObXFlG3dds5zv6lug3OmMWxL11k2z5MyAACAhgRlAAAADQnKAAAAGqqSU3ZwcBA7OztLl9XqZ1uaE7SpHJBVtiVXJ0dJ7tImh4ndlukXSrZd3HfpNS0Z/jhn3U31h8/9DNhUv/w5p+dYZeqzqVae5OHhYezv7x+zhlwIauYS15xi46QomTZk3XLXKXtKSS71nGMC1JyCZa565V6H8+W+XnQ+HtMcPCkDAABoSFAGAADQkKAMAACgoerzlI3Vmjup5rwXJTkgc/bRHtvUnF81539aVY9NnZ8595uzfm4eVKt++TWV1HNTc7/MmQd2UuaTW7ce5iljk0rn4DzflJyPOc9VyX7nzC9qdT6mlN7Tc+WH11Tze/WFzJMyAACAhgRlAAAADQnKAAAAGqqeUzbXnBGb7HM6Va9NziEyltMHNyeXqXQOuJx61Lw/at4TNfPTpu6XmvudM8+p5F6bay6cbZlzZpU5P+c29RlpnjI2Jfe9ljM330nMVanZHtf8jCz5TrEtbXfOtqU2lbd/Uuc429Z6teZJGQAAQEOCMgAAgIYEZQAAAA1VzynL6RtbMwekpO9wSd7TnDY1x1dpf+da8y5tqp95xHx5PePt57zHS/MKj7vfVfuZa66x0ty2kjzKkvfenJ8ntXLMTkoeAlxo85SV2OT7Ouc71qZsy+dYTtuUa1vnXsvZfluu0zbwpAwAAKAhQRkAAEBD1bsvjpV01an5SDOnG2VJF8ySx9SthqYvGZJ4VVmrti3Zb+72U2XVnPYg5xrP2W2hpMtuTrklXWRaTXWRe95zrmnNoaNzlJzLqXrs7u7G3t7escuGHDWnUTkfXQjHeL4ruYbbcs+XTLOzbDk9T8oAAAAaEpQBAAA0JCgDAABoqHpOWc3clJq5F7Vy21oNS15abs3+vCdhSOKa+TWrltc8Hzn3WsnUBHPmfZXk1NXMOaw1Dcaqfc85dcXYXPfe1DEcHh7G/v5+lf1Aqbnei3ASbOs93eq78PnGkzIAAICGBGUAAAANCcoAAAAaqpJTdurUqTh9+nRElOVilMjZb+58VzVzdTaRA1K6nznzfmrO4ZUzf80qJXlAc81rV/N8bDKXa645z0rmAFy13znzBHPk9rufa84m85SxrXI+5+WYASeJJ2UAAAANCcoAAAAaEpQBAAA0VCWn7ODgIHZ2drK3y8kn2ZZ5tkrz0WrlLpXk09Q2lQeVs23NdUvzsXLKnlreaj6942x/rm1zc7k2lWM351xjOeY67+uUnZNzNza17SLzlNHStuR7AszNkzIAAICGBGUAAAANzT4kfkl3rBxzdj+q2R1par81y5pzyPy5upnmDr2eU9ac29bs+pizbW43wilzdumdqzvnnMNb1+z6WDJdxybfx1P7gW2R85mgqyNwknhSBgAA0JCgDAAAoCFBGQAAQEPVh8TfVA5Vzb7icw63P1VW7vD5JeejZt5Tq3y0TeUUlR5/Tg5VzhDnU/tZVtZcOUMl56dkqoLSaQ6m6l2Sy1XzfKzSKp8TtsWmvhcAbJonZQAAAA0JygAAABoq7b54UUTETTfd9OEXdnd3z7ny4eHhZGFT25Za3Peq/ayq55Rx2VNl5ay7bP255NZj6tyOy5pr25w6lsq5f2oeU2m9cupRUqecsnLOT06d19n3lJJ6rSorp04177Wc/Sxuu/D5ftHaBbKus9pQppW0gzXbAYB15LShqSS35Lrrrrs8Il527AIAOCnuft/73vf61pU4n2hDAS4YK9vQ0idlr46Iu0fEWyLig4VlAbB9LoqI20f/eU9d2lCA89vabWjRkzIAAADKGOgDAACgIUEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoCFBGQAAQEOCMgAAgIYEZQAAAA0JygAAABoSlAEAADQkKAMAAGhIUAYAANCQoAwAAKAhQRkAAEBDgjIAAICGBGUAAAANCcoAAAAaEpQBAAA0JCgDAABoSFAGAADQUJWgLKX0+JRSd46/766xjzmklC4b6v5PK5Z5+XDcV5xj+WcPy//tOZbfIqX0xpTSU9bc3xVDeZcfv9Z1pJSuTCm9MqV0mFL6q5TSo49RxuOH7T95ybJvTCl9MKV09zo1biOldNVwzXa3oC5fnVL6s+Gcvyal9LA1t7ttSulnU0pvTSm9K6X0mymlO89dXzgfaUPPKPOCa0NTSrdOKf3HlNIfDZ+nf5NS+vZay9esw32G8/BVS5Z9Rkrp/bllbpvhfu1SSg/dgrrcI6W0n1J6b0rp9Sml70kprfWdPKV0dUrphqHdftm53iucPDWflL03Iu655O9/VtxHbZdFxPdERLUGZZWu614eEX8eEY88xyr3iYg7RsSTN1WnGlJK94yIZ0XEKyLiiyPi2oj4yZTS12UW9V8j4oaIuGZU/sdFxA9GxDVd172stL5EDA3TtRHxjOiv2XUR8WsppfutsflTI+JBEfFNEfHwiPhnEfHClNKt56ktnPe0oWs4T9vQT42Ib4yIn4/+c/UZEfGjKaWvrLR8pa7rficinjRsd+nR6ymlFBE/ExGviYj/fszjY0FK6RMi4oUR8daIeGBE/LeI+E8R8X1rbHt1RPxQRDwxIq6MiL+KiOenlO42W4XZnK7riv8i4vER8a4aZW3yLyKuiIguIi6vWOblQ5lXTKzzuIh4f0TsLln2ixHxloi4ZatjOOZxPz8iXjJ67eci4o0RcYtjXpevWHjtyRHx+oj4qJmP49YbOFdXDcd31vXf8DX704h42ui134qIP1yx3b8a6n/lwmsfP9zT39DymPz5O4l/2tAzyrzg2tCIuM247YmIV0fEr9ZYnlGPjxnO3c8tvPbvIuKDEXH3mc/BrXK/KxxjH5cN1/qhra71UI+fjoi/i4hLFl77L9H/MHO7ie0uiYh/jIgnLLx2UUS8NiL2Wh6Tvzp/G8spSyk9eHhs/MCF126XUjpIKf36wmtdSuk7U0pPSCm9JaX0zpTStSml247Ku11K6ZqU0t+nlG5KKb182S/8KaUHDI+I35NSeltK6cUppc8aHvf+zrDay466ihyj/O9OKb1p6DLw9Ij42DVOx5Mj4uKIOOMRekrpkoj4soj49a7r3j/U/bdTSm9OKb0jpfSSlNL9pwo+1+P5lNITU0o3jF67U0rpV1NKp4dH6L+bUvrsNeo/3uclEfH5EbG35DjvGBGflVNe13Uvjv4Jzo8P1+HzI+IrIuIxXde9I6V0aUrpF4b7470ppZeOr83waP+Jo9ceOpyby4Z/H52rq1JKP59SemtELH0KN9yDr17y+oOGMu46/PurUkq/n1K6ceF+u8fU8aZzdJ9JKT0npfTi0Wt3SSk9K6X09pTSu1NKz00pfeJU+efY5z+PiE+L/onXoqdExD3SdNfKK6NvGJ5/9ELXdX8XEb8fEQ/IrQuwmjb0DOdVG9p13bu6rnvvQrkfGRF3iIi31VieUY+3RsTVEfF1KaV7De3Af4uIn+q67mWp7xr6X1JKrxuu6V+mlL51sYzhXnvN6LXdo7Z24bUbhnP6HSmlv40+IPmYcZ2G9vkDqe8ts/j6pSml96UhTSKldM/Ud6N/49A2/lFK6VGrjnmo19Wj165evJeH19a6n9d0ZUQ8o+u6mxZee3JE7ETEfSe2u1dE/JNYaLe7rvtgRPxaRFyZUkrHrA9bompQllK6eMlfiojouu6Z0X/R/oWFL3zXRB/lf8OoqG+OiLtExFdHxHdGxEOifyx/tJ9bRcRvR//Y97si4kui/6XguWnhEW5K6eER8eyIeHP0X+q/MiL2o+9q9YroH/dHRHxN3NxVJKf8b4qI74+IXxnq+Lronw5N6rruryPiJUOdFl0ZEbeLm7td/POh/o8ayt+PiOelCv2HU0ofHf2X6M+M/nw/JCLeHREvSil97MJ6R7kOl00U94nR/8r1p6PXXzv89y7HqOLV0d8bPxr9ffKbXdc9PaV0UfTBwJdGf20eEhH/EP15uc8x9hPRd5nsou8O8x3nWOcpEfHpKaVPH73+iIj4467rjo71soj45Yh4WPTX9yAifjel9CnHrNuHpb7Lwx9ExKXRP237ioi4fURcN3wZOVrv2nGDssTRNVl2zVL0AdvUtn/edd14H6+N411rILShceG2oYvlXhT9+bhF9O1f1eWrdF33pOiD7Z+OiJ+IPlj6rmHxj0R/vX41+m6Sz4yIn0gpPTZ3P4OHRH+PfEtEPDgi3rNknadH/1R0nO/8kOjbqqMfJO4c/fX9uqFuvxERv5iW5Mjlyrifr0oT+ZDDOh8Zfc+SM9reruv+Nvrjn2pDj5b92ej110bEbaN/X3KS1XjcFn3Xi+4cf1csrPdR0ecLPT36PJQuIr54VFYXEX8TERctvPa1EfGhiPi04d9fE/2b9K6jbV8SQ3es6N+sBxHxgol6XxFLui2sWf5FEfGGiPjl0TpPGR/3Ofb9mOGYTi289rToG6W0ZP1bRP/L4G9FxFPOdQxxjsfz0fc/vmHh398b/dOOj1147ZLhnC0+Gj+6tpdNHMu9h3U+Z/T6xcPrjznmffVVw/bvPDpP0X8YjrvO3SL6D6UXL7x2Q0Q8cVTeQxePZeFcPXeNulwUffD3gwuv3Xqo23eeY5uja/ZnEfFDC69fFQvdFyfuw+eMjulJ0b83dhZeu31EvCsiHr3w2rUR0a04nq8c9nmH0eufNLz+JRPb/nYseV9FxA9ExI3Hudb+/F3If6ENXVzngmtDR/u5Jvov5587x/I16/DJEXE41PtLh9d2I+J9i8c2vP6z0bdBtxn+fW1EvGa0zu5Q1lULr90QfVfJj1ijPr8REfuj1150rntzuHcvHur2Bwuvn3Vth39fPdr+6lhoQ9e5n4d/X7Xq3o0+cOoi4hFLlr0+In58YtvviojDJa9/wVDmZxz3mvvbjr/aA33cfcnfy49W6LruHdH/cvdvov+C+dNd1z3/7KLi2V3/SPbI06N/kx11A7tf9P2l/2LxF8XoByo4GpnvUyPiThHxS8c4lnXKv1P0yc3PGG37v9bcx170DcrDI/rR7KL/Feap3fAuG7pGPCml9IaI+ED0Hwr3i4jipy5DOb8TETcuHN8HI+L34uZjjK7rHt91Xeq67oY1yuwyX58urOt+OfqctKd2XXcwvPyvI+KdXdc9b2G9D0XfGN9r+JUw1/NWrTDcj78ew/UaPCgiPjIWuhKkvnvhM1JK/xD9+Xx/9PdirWv2rIj4wMI1e1tEvCrOvGZXdV23bjeG8bVJ53h91XZH2x7rWgPa0MEF24amvhv8/xsR39Z13e/XXr6uruv+Mvp75i+7rju6Pv8qIm4ZfVe5RU+Nvh3MSlMYvLjrumVPx8aeGhH3TCl9fERESukOEfF50QfwMbz20Sml/2/oCvn+4e/ro961XnU/R9d11w7X+sVrlHncNvRc251rGSfIxRXL+lDXddevsd4fRP8LySdE/8vTMm9e/EfXdW9LKb0/+vykiP5Xl8+K/k03dtQQHfVNfuMadRpbp/yjurx5tPwf1tlB13VvTim9MPruFz8afXe8W8fQ7SL1Q6P+ZvT9hx8X/Qg7745+dJ6PX/dAJuxGxOfE8mP868yy3jb896NHr3/0aPlxvD/6X+cWy1x2jt8UfYNxm4h4e+Y+xtfwXJ4SEd+YUrpH13Uvjb674//t+m4HR18K/nf0v/79h4j42+h/bfyF6PuKl9qNiG8d/sbeu+S1KYvXbPF83m60/FzbLrsHb7diO+DctKG9C7ENPXIUQLxwpuU53hdnt70RfVu76Ojfl0a+ddve50TfK+UREfGE6APx90XfffLItdHnXH1fRPxJRLwj+gB18YfU41rnfl7Xub4vRaxuQ98WETsppZ2u6w5H2y2WzQlVMyhb1+OjT+T9y4i4JqV0n6NftRackeg79N2+ZUT8/fDSjRHxx9F3yTiXtw7/Pc5QveuUf1SXcVLyx41XnPCUiHhSSulTo29YXtV13Z8Myz4p+g+BB3dd96yjDdLqIceP3qi3Gr0+/sC8MSJeEBHL+oLftOS1KX8d/QfkXYYyj9x1+O84b6nEjbH8HN8h+g/Mdw3/PozV5+DIur8u/d/oA61HpJT+LPph5BcThO8Z/a+/D+y67lVHL6aU/kn03RLOZeqaLTaKN0bEc2M0XcDgnescwIKja3KXOLN/+l2jPx/jPuvjbb8wpZRG7927Rt1rDZzt8aENPXK+tKFH3hP9cP+HMy0vcePw34+LvtvpkTuMlldve7uuO0wpPTNuDsoeEX3awTsiIlJKO9EPMvXtXdf95NF2ab15v25ao77r3M9r6bruPSmlv4tR7ljq5/n8iJhuQxfb7VcuvH7X6L8DvOGsLThRNjb6YsSH57L6TxHx7dE/ZbhXRHzbklUfNOqG9mXRv3mPRsZ7YfS/Er6x67rrx3/DOn8e/Rfhr5mo0tEX3vFTjHXKf330jcqXjrbNmZTwGdE/4fiW6EfcWZzs8qjh+PCX8uFNe+8VZb45bg6Qjra7JPpuf4teGMOX6CXHeNYog1O6fgShF0XEl48WPTL6c/TKszY6vt+PiNumhRG0hg/eh0Xfd/zoV6vXx9kJs19YsuPhi89e9L+8PTT6HzWetrDKsmt2r+j7sU85CtgWr9nHRsRnjNZ7YUR8ekS8csk1+/PMY3ld9IHX+FfER0bES7uuOz2x+fOi/2XuixbqeyoiPjf6oBGYgTb0LOdFG3qk67r/3XXdp3Vdt/TLdenyQi+N/ofPcTv/8OifQL5i+PfrI+JOKaXbLKxT1PYOnhoRn5VS+qLon1AuXutLos9RXLzWt40+B32VZd8VvmD073Xu5xzPi4gHDwOIHHlk9AHtdRPb/UH0PYE+3G4P7/Mvj4jnLflxhpPmuMloi3/R/3L3nujfKOO/TxzW+cjof9l73sJ2j43+JvwXC6910Uf7z47+ScSjo/8FYDGZ8pKIuD4i/iL6PsNXRD9yz/dGxH9dWO/h0fc5/43o++Dff1jngcPy3ej7mf9S9E85Ls8s/zFDfX8k+j7HPxb9G3xlkvJCGUf94scJy0cJw6+Kvp/8w6P/Ev26WEiijSWJ1kOZb4s+9+AB0Sc2/12cmaT8MdF3gbk++pGpPi/6xvBHou+PfrTe44ZzdOcVx3HP6D+wf36o03dF/1j/60brvXixHmucnxtiYcCO6D94XxJ9F8F/N9wjzxrqeMXCet8wnJfvib5B+O/RP+X6cMJ1HGPOkugDpaN79LdGyz4u+nv1RcP98DXDPl8fEc9ZWO+qGM1TFhF/OFyjh0b/JeUlwz5evLDOJw3X9YXRfwh/3nBfXBMRj1xY7xcj4gNrHMvDhvvuB4dr9hPDv+83Wu8DEfGLo9eeM9TvEdGPeHZ99O/v2ed58+fvfPsLbegF34YO637V1LolyyN/wJFr4+wBO350KP97h+v1w8N5f+zCOncdXnta9G3vt0bfnbCLswf6eOI6dRnWvzj6gPkN0QcmO6PlL42+vX3ocK/9YfQD3rxrYZ3L4uyBPv5b9O+hb47+h8YnD9e6O8b75ej8f96KY/mE6LtXPjP6HxMeHX1g+wOj9a6LiL8avXZ19MHnt0fEfYb6vjci7lbr88hfu786hUyPHHXtsM7PRN8d4o4L2100vHFeHsNEj8M23xn9h/Nbo29MfjlGkwZHPwrVjw9vwvdF3+/9uRHxgNF6Dxr28d7oP2Svi4jPXFj+76Pvfvf+0ZtwZfnRJ1c+Lvo+8O+OPjh4QOQ1KA8a1v8/S5bdffigee/wYfBVMfqgjOUNyu2j/wXx7dE3cI+J0chRw3p3iD7f6Y3RP8I/iH4wi3stubaXrXEsV0bEHw1l/XVEfOOSdV4WKyYnHq1/Q5w9iuKl0Qcep6P/MH1pRHzRaJ2Lo28c3xT9CFk/E33DWRSUDdu9JkYNzMKy+w/L3xv9l4Evjj4QXRWUfWL0wdy7ov/i9fAYjb44rPfJ0SdaHx3766JP+F/8UnZtrBh9cWHdr47+F/Gbom84H7ZknQ+/j0fvj5+LvlvHu6LP3bhzjc8Tf/4utL/QhmpDuzPahqXrliyPvj08jInJiUfrn3GehtduEf0PrjcM1/SvYiEAXVjvUdG3Y++JPs/6s6MwKBu2+alY0h4Nyz4p+jb03dEHVVfHaFL2WB6UfWT0Pyq8Nfqg7wci4j/GqA1d834+Ov8r793oB075g+GavGGo60WjdV685J5L0U/dc5Sz/rKIuE/OefS3vX9puMhbY5hf6Tu6rsueX4PtNvT7fntEPKrruqetWh+APNpQlkkp/V5EvLrruke3rguwXIuBPrhw3SP67gTrDnkMABQYcpf+ZUT829Z1Ac5NUMbGdF33uzE9Wz0AUFHXde+LvvsdsMW2rvsiAADAhWSjQ+IDAABwJkEZAABAQ4IyAACAhgRlAAAADQnKAAAAGhKUAQAANCQoAwAAaEhQBgAA0JCgDAAAoKH/H2HRwA2+b8W3AAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1100x1100 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# you are welcome to solve this in any way you like, I've just filled out a very simple way to do it\n",
    "def energy(state):\n",
    "    E = 0\n",
    "    N, M = state.shape\n",
    "    for i in range(N):\n",
    "        for k in range(M):\n",
    "            # your code goes here\n",
    "            pass\n",
    "    return E / (N * M)\n",
    "\n",
    "\n",
    "expected_values = [\n",
    "    E_prediction_all_the_same(100),\n",
    "    E_prediction_all_the_same(100),\n",
    "    0,\n",
    "    \"???\",\n",
    "]\n",
    "\n",
    "f, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))\n",
    "axes = axes.flatten()\n",
    "for ax, state, exp_value in zip(axes, states, expected_values):\n",
    "    show_state(state, ax=ax)\n",
    "    ax.text(\n",
    "        0,\n",
    "        -0.1,\n",
    "        f\"Expected Value: {exp_value}, Your value: {energy(state)}\",\n",
    "        transform=ax.transAxes,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9d1d3514-aa70-4602-97a3-5d31445635a7",
   "metadata": {},
   "source": [
    "## Solution 1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "ad4bf3dc-d1be-4253-b5ab-18086b688d27",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA2UAAAOCCAYAAAD3N3ppAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABDrAAAQ6wFQlOh8AABihUlEQVR4nO3dedgtWV0f+u+ikT4qKuKL8wk4C4nTFZzziGIQGRwiBtSLtlejXozGARJzjYgmGmdjgq3GIY0ROGrECXECJMYTRcAJnAcaj6LCoZFBOS1g3T+qXnr37vfsvetdq/ba55zP53nep/vsqlo1rarav13rt1YZhiEAAAD0cafeGwAAAHAtE5QBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEd3rln4mc985vVJ3ifJy5K8ockWAXBIrktyjyQveMADHnBr7425mniGAlz1dn6GVgVlGR8mz60sA4DDd78kz+u9EVcZz1CAa8PWZ2htUPayJHn0ox+dW265pbIoAA7N3e9+99x4443JdL+nqZclyX3ve99cf/31vbcFgMZuvfXWPO95z0t2eIbWBmVvSJJbbrklFy9erCwKgAOmeV17b0iS66+/PmfOnOm9LQAsZ+szVEcfAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANDRnSuXvy5J7n73uzfYFAAOzcr9/bqe23GVui5Jbr311t7bAcACVu7vW5+htUHZPZLkxhtvrCwGgAN3jyQv7r0RV5l7JMnznve83tsBwLK2PkNrg7IXJLlfkpcleUNlWQAcnusyPkxe0HtDrkKeoQBXt52foWUYhuU3BwAAgBPp6AMAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANBRk6CslPL4Uspwmb9/32IdSyil3Gva9ndsWOZ9p/2+/2Wmf+A0/f++zPQ7lVJeUkp58o7ru/9U3n1Pv9XLKKXcs5TyU6WUPy+lXJr260dKKe95wrwPLaX8einl1lLKhVLKV5dSrjthvjctpXxdKeXF07w3l1K+6oT5HjNNu1RKee5J56OU8hallO8upby8lPKaUspPllLuOXMfj0opF0sp33/CtDebtuGn5pR5iKb9eELv7ThJKeUbSim/U0p5dSnlVdP5fuQJ871TKeWHSimvnOb9yVLKu1ymzE8opfxqKeXvSim3lFJ+oZRy97V5PqiUcr6U8tqpjn9VKeUO99RSymeWUn5/qosvLKV8yo77VV0/uTJ4ht6uTM/QiWeoZ+g+eIYejpZvyl6b5ENP+PvvDdfR2r2SfFWSZg+UbYZheH6SP0jyqZeZ5aOSvEOSJ+1rmxZ01yR/meTfJPnYJF+W5L2S/GIp5eh4plLKhyT5iSS/l+Tjk3xrkscm+YbVwqYHzNOSfHKSxyV5YJKvTPL6tfkek+TrkjwhyYOT/HGSnymlvM/a9j0lycOS/Kskj0jyTkmeUUp50113cBiGi0kek+SGUso/XZv8uCRHSb5g1/I4lTdP8l0Z68WnJPnNJE8ppXza8QxT3fmZJPdN8nlJHpXkbJJnlVLuulpYKeWGJD+U5BlJHpLkM5P8dpK7rMzzrtP0lyd5aJKvT/Jvk3zNWlkPT3JTkh9L8nFJnpnkh0opD9xhv6rrJ1cUz9AdeIZ6htKcZ+ihGIah+i/J45O8pkVZ+/xLcv8kQ5L7NizzvlOZ998wz+OSvC7J0QnTvi/Jy5K8Sa99WPiYv8e0vZ+28tnPJnn+2nyPSfL3Sd5u5bPPTfKK1c9OKP/6JH+T5BtXPrsuye8mObfy2QdP2/Hglc/+0XRePv8U+/WsJC88Pm9J7jNt/5cufDxLkuv3cN5uTvKE3vVnxvaeT/LzK/9+5HS+32fls3dKcinJl6x89jZJXpXkc7eU/51J/mz12Cf5/zJ+sb7byme/l+SH15b9uSS/uqX8pvXT32H/eYberkzP0M3b6xna9nh6hp68vZ6hHf72llNWSvnEqYnAQ1c+u9v0mv1HVj4bSilfXkr5xlLKy6ZXpDeVUt5irby7lVJuLKX85fQK/vknRc6llIdMr0f/rpTyilLKs0spHzC9iv/FabbnTusdTlH+vy+l/NX0avSpSd52h8PxpCR3TvLwtbKuT/LPk/zIMAyvm7b9F0opL51eKT+nlPKgTQWXsTnJMP26sPr5E0opN6999s6llB+cmg68tpTyS6WUD9xh+2u8fPrvm6x89gEZL7JVPzvN87Ern312xovzrzeU/2FJ3irjLyRJkmEY3pDxV5sHl1LK9PGDMz54fmZlvj9L8ssZf9mZ6/OTvHvGXzKT8YbzwiTfnryx/v/G9Pr9r0op37H661Ip5YbpvB2tFjq9qr9p5d83TZ89uJTyW0luzfjLaNaWO64Hn3zCtOeXUn54+v93KKV8fynlT6c68EdlbNpy/aadna6jp619dodmR2X0mFLKH07X0Z+WUr5kU9mNvDx3rGN/OQzDC44/GIbhLzKeo4etzPcvMj6kb9pS/oOT/NgwDLeufPakJGeSPCBJytis472zUhcnT07yQevn+oTy/ybt6idXOM/Q2/EM9Qz1DF2WZ2gHTYOyUsqdT/grSTIMw49nPEnfu3Igb8z4C8znrxX1hUnunfGV55dnfKX6PSvruUuSX8j4yvMrMl5Qv5vkp8vK6/VSyiOS/FSSlyb5tCSfnjH6f6ckv57bXol/Vm5rKjKn/H+V5D8k+R/TNr4oyX/bdpyGYfiTJM+ZtmnVg5PcLbc1u3iXafsfNZV/PsnTVy/Y0yqlvHXGyvn+GY/3Jyf524yvot92Zb7jXId7VazrTqWUN5nKeEKSC0l+fGWWMxl/EVt1fKHeeyrjLkn+ryQXSin/Y/qC8KpSypNLKW+zsty9p//+/lp5v5vkLTKe++P5/mCYfj5Zm+/emWkYhj/M2NzjK8vYPv8jMv5S9IZSyscneWqSP0zySRnrzKNy+2MwxztmfFB9a5IHZWxqsL49Nyf5law18SmlvEfG43h8kztKckuSL53K+saM1913nnLb1n17xuYIT8x4I7wpyTeUUt54za88/B5/2pVMD647T18EH5WxWc53rMxyUh1Lxnq2er4/JGPTqBtKKX9WSnnd9AD+6JV1vXnGX9x+b7WgYRhenOTvVso7/u/t5stYx0rGh83lNK2fXBk8Qz1DL7Muz1DPUM/Q21y9z9BGrzkfn/E14Ul/91+Z7y0zvsJ9asb2nUOSj1sra0jyp0muW/nss5P8Q5L3nv79WRlfQd5nbdnnZHrNmfGEXUjysxu2+/45odnCjuVfl+QvkvzA2jxPXt/vy6z7i6Z9Orvy2Q9nfCiVE+a/U8ZfBn8uyZMvtw8Z2/gPSR6+tvwTkty88u+vzvgrwtuufHb9dMxWmy0cn9t7VdSPH1ypD3+c5D3Xpj83ydPXPnvUNP93T/9+h+nfr874kH1gxhvfS5P8wspyX5Hk0gnb8DHT8u87/fsXTqobSf5jkltOuZ93yXjzGJJ8+8rnv57kOWvzfupqPUlyw/Tvo7X5XpjkppV/3zTN90E7bM8XZmwK8BYrnz1uOu8nNteY6tinTfX/zVY+vzkrTS+SPDvJ09aWvV2zoyTvNtXxz12b7xsz5kncaa3OPr6ijh2f32Ha9s9bm/4FGfMm3nHls7tOx+LWlc9+bqpjL5nq4Mdm/MX5UpJ3neZ5p2k9jzxhO/48ybdO///p03xvvzbPu0+ff/yG/WleP/0d7l88Q1fn8Qy947Z7hnqGrs7nGXqVPkNbd/RxvxP+nn88wzAMr8p4E/iEjFH/dw7D8DN3LCo/NYyvy489NeMD4oOmfz8wyQuS/OHqL4oZEwDvN83zXkneOckdevTZwS7lv3PGX1t+bG3Z/7njOs5lvNgekYy9xGT8VfEpw1Rzytg04omllL/IeDG8btq2O/S8dAoPzNj05JaV/XtDkv+d2/YxwzA8fhiGMoy/Gp1o+hVv9Zfd9Xr1lRnP3cOT/FXGRMt/tDL9O5J8XCnlX5dS7l5K+YgkXzttzz9M8xyX+TdJPnkYhp8fhuGJGX8h/phSygetlDectJknTLvcfCd9vtUwDH+f5Numf35DkpSxecX7Z/yysOpHMp7T9cTmXVwchuHXdpjvhzM2P/jElc8ekeSpw9RkYPp17ItLKb9bSnltxjp23DToXU+xbas+Zvrvj55wHb19xiThDMNw81THHr+psLU6tt6r2HMy1tuPyfjL4hNKKZ+9Mv3JSV6Z5KZSyruVUt4p45uDu+a2OpaM9eyuST5nGIb/MQzDz2X8ZfZvMuZorNq1/qz/+6S6eJKm9ZOD5xk68gz1DE08QxPP0MvNd9U+Q1sGZf8wDMPzTvh79dp8/yfjLwbXZ/zl6SQvXf3HMAyvyFjR32H66Chj+9bXrf39u0yVNGOyYTJG63PtUv7xtrx0bdlNbbXfaBiGl2bseea4+cUnJXnTTM0uppvyT2Z8hf+4jD1K3S9j+9gzc3foBEcZbzTr+/ipuW0fd3WcdH3897jVicMwvGgYhucOw/CjGX81uUvG3qSOPTFjM4JvztiO+ZkZewK6JeMDKBmTk5Pk/HTjPvas6b//eGW+M6WU9WN0t7VyXpHkrU/Yl7utzHMaf7/237tlvAn81epMwzC8PuO+3q6L2B2t17kTDWPewLMyNb8opbxfxuTp1fbZX5zkWzL23PUJGR/8x02SauvZUcZ9v5jb14+fnabvXM+mZjurZfzJ6vRhGF493W+eOQzDYzLWn289fvBM95BHZtz/P874a9w7Zqx7q+fmlum/v7hS9muT/GpuX8eS7fXncvPdbW36SZaqnxwuz9CRZ6hn6Op/7xbPUM/QO86zOv0kV+Qz9M4d1vn4jIm8f5TkxlLKRx3/qrXidom+U9vtN8n4ujYZT/pvZ2yScTnHybCn6ap3l/KPt2U9KfntZqznyUmeWEp5r4wPlt8ahuF3pmnvnvGh9onDMPzE8QJle1eel6b/3mXt8/Ub1y0ZL+yvPKGMW0/4bJP/lrGb3WOXfYgPw/C3pZTfz7h/x58NSb6slPI1Se6ZsUeeN8n4S9+vTvP8XSnlxRu24fiXmuO2x/dO8hsr0++T8ZX6X6zM989KKWWt/t0nd2y/XONvMv4qc7t6Mf3a9Ta57Qa263lL5v3K85Qk3z3lDDwy4xeeZ61M/5QkPzkMw79b2bb77FDupR229ZZpWz8iJ7dF/4Md1nPsJVn59Tnb6+jzM3aDe49MD4xhGH6hjGOUvGfG5jkvKqX8dKY6NjluOrOuZKpjU138s6y1S5/KfrPcVn9W6+JqfsZ9pnWs52ys2lf95Mrz+HiGHvMMjWfo9JFn6Gaeobcte7jP0BZtILNjd74Zk4Bfn7Fb1g/MCd2dZnN7+HtP//6XGW8O77hhXcft4X9mwzwfNq3vI9Y+36X86zJW8lO1h5/mfYuMSY03Zvzl4t+sTHu/rOULZLzZ/n2SF658dv/cvj38nTJebP9hZZ7jdu43r3z2tUlenOTNW9SBGXXlbhl/9blxy3xfkzE3YLUePCHjrzOrXag+fNr/91/Z179J8vVr5+p3cnJ3vg9a+exsKrtLzQnt2jO2h//VtfmO80HuP/37I6Z/P2Blnn8y1fubVj67afX877A9b5Xx5v9503X1X9am/0aSH1r77OdX69T02c25fXv4753qT1n57JvW9uk9p+1/2D7r2Mr2vTLJnTfM897TtfLRK58dt+l/6Mpnb5bxl9X/vPLZd077f5eVz748J3fne25tvT+7Xh9O2LZF6qe/w/yLZ+jq556hm/f5bvEMTTxDl65nnqEd/lqdvMdnvDF+yAl/7zbN8+YZf9l7+spyXzlV9n+88tmQ8ZeYn8o4UNyjM97cf3hlnuuTPC9jTzyfm/Gm+okZE2//08p8j5gq9I9mfK38oGmeh07TjzI+4L4/48PuvjPL/6Jpe78pY/vyb8l4w9vpgTKVcdwufj1h+fgh8FsZ28k/IuOvAi/KhgfKSpmvyJh78JCMiZd/lts/UN4m403ieRmTMT8y4835m3L7MSceNx2je56yXnzbVO5HTtvzGxkv9Pdcme+DMg50+c8y9tL1vVm70Kf57pnxYfH0qW58dsbxaH5sbb7j8Vm+LGOTlSdlvNDfZ22+p2Wsa4/M2GvX8zLW0Tc94fjesOM+35A7PlA+fjq/T5nq4KOnY/CMlXneZDpHvz2ds0+dzv3FVDxQpmWeOu3nkORD16Z943Ss/1XGOvzEjA+ebQ+UB03zPCFjG/Svytgc4nZ1P8l/nc7ZV0zzfVySf53kx9fO6+uTPO4Udex9MzZH+n+SfPR0rL9n2o5/uzbvN2Rs4vTRSb4k45uA7zmhzB/P+MvgDVO9+IUkr0nyTivzvGvGsVh+PGP3vY/O2PPaf1wr61Omc/+1U136tunfD1yb7/VJvm9u/fR3dfzFM9Qz9PL1wjPUM/Rv4hl6TTxD2xSyueeom6Z5vms6ge+wstx1GV97Pj+3DRg4ZIyWv2Wa/9VJfiDJW66t8y0ztqF+ccabx0uS/HSSh6zN97BpHa/NeJN9ZqZfhKbpnzddCK/L1Apg1/Iz/pL4uIyvs/82Y5vih2TeA+Vh0/z/64Rp90vya9O2/2GSz8jaDSUnP1DukTF5+pUZH3BflLWeo6b53j7jzfslGW8qFzImz37YCef2XqeoFx+fsYehixm/OPxJkv+eqQeelfnefzpHr57+npG1G9/KvB+Y5H9Nx+Riku/OSs9IK+flsdO5u5SxZ6qPOqGst8zYbOSWjDeMn8zag3PlfD5ox32+IWsPlOnzT8r4ML11qi/fkeSuJ+zbr2X8cvbbGW98J/UcNfeB8inTNr3ohGl3nc7JLdPff8v4BWbjA2X67LEZH4KvyTiGzQNzxwdKyfiwesG077dk7GZ49UvLvablHn+KOvZ2GR/UN0/n+q+n+vEJJ8z7pIwPilszNvv4sqz8irx2TL4j45eVS0l+KckHnzDfB2fM77mU8cb/+MuU95nT+m7N+Gvzp5wwzxvvlXPqp7+r4y+eoZ6hJ++bZ+htn3uGeoZe9c/QMm34wSjj4JOPHYbhm3tvC5RS/kPGh8H7DId2sQCs8QzlkHiGwu6aDh4NV6EPT/J1HiYAMJtnKOyoR++LcMUYhuGjt88FAKzzDIXdHVzzRQAAgGuJ5osAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEeCMgAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAju5cs/Azn/nM65O8T5KXJXlDky0C4JBcl+QeSV7wgAc84NbeG3M18QwFuOrt/AytCsoyPkyeW1kGAIfvfkme13sjrjKeoQDXhq3P0Nqg7GVJct/73jfXX399kuTs2bO3m+HChQtv/P/1aetW5902//q8c2zaxlrbyl6d3nK9c23ajjn7cNL0Xdc7V8165tStuWrqbcs63vK87Ou6rTlP2/av1XnYZsnruGU9nnMvWp1297vfPTfeeGMy3e9p6g7PUACuHrfeemue97znJTs8Q8swDKde0TOf+cx3TnLhwz/8w3PmzJmxwFJuN89q+evT1q1vy6b5a7Z70zbW2lb26vSW651r03bM2YeTpu+63rlq1jOnbs1VU29b1vGW52Vf123Nedq2f63OwzZLXsct6/Gce9HqtKOjo5w7dy5Jzj7gAQ/4851XyFYnPUMBuHpcunQp58+fT3Z4hta+KUsy/qp68eLFJHVfGuZ8MWz5ZaUmkNq2npZfdudYav/nrqtmH2q+dLe07XjMCYaWrLebprc8p0tqWW/3dR7mqvlBZI5ty865FwEAy9L7IgAAQEeCMgAAgI4EZQAAAB01ySm7cOHCTknKS+a1tOwwYE4e2Ny8p6VyN/a1/9umt+wUYl1NpwZz1B6Pmhy7OXlQ25ZdnX/J+rFkbmRNOS3r1pxju23ZVvOuz1+T67jPnDqYY8mOuQAOiTdlAAAAHQnKAAAAOlq8S/xVtU3/Nqnp/rlll/hLNlWqUbNdSzaDqmmuV1NWTXO0mm1achy7lmNttWyet9TYfC3rR03ZvbrTX193TX3RJIxD0WsIFoDevCkDAADoSFAGAADQkaAMAACgoyY5Zatqutaeo1c+2rZpc/KvavLRaro4n7uull3ir9tXfdm23pphD9bNOR412zHnHM/No2yZ61ej5XlZta9u/Ftvx1L7rAv8w7XUPXLJZZe8fmqeNwCHzJsyAACAjgRlAAAAHQnKAAAAOmqeU1ajpo37vsYam5sDsmTeWM121IxRtNTYWrU5ZK3yC+aew5Z1rWVex77G89lXXkftsdw0pldN7t+27dw0vWYf1qcvlVd66dKlnD9/fqd5OZ3VsT7nqLn2llx2yRwzeWTA1cqbMgAAgI4EZQAAAB0JygAAADpaPKesVQ5Vy3GFeuWPrJe9z3b2m7Zjbhv9feXFHcp656oZW6tmXLs5OUU1uUrrasaxW7I+LDVG4NyyasaZqtnHues6bbnUu3DhQs6cOZOk7diFvfIXe91PAK5k3pQBAAB0JCgDAADoqEnzxdWmF5u0bOYzd12nXW/LZfep5bHd1tV6TXOSmqah63p1D72prCXr7TYth3bYpOa6rjnuNV3A1zbPa9U0tOUQCetaNn1kWZu6xK+5fq7EZbc51GcuQC1vygAAADoSlAEAAHQkKAMAAOioSU7ZpvbwrfI4atuR13Svfyg5Z3O6JN5mTlfJLbvdPpTujXud05bDLSxZ9lLn4lCurbn3gH3letV2W77Jabf56Ogo586dO9Wy7GZTl/ibHEpe4D6vY13iA1crb8oAAAA6EpQBAAB0JCgDAADoqPk4ZTXjk9S0j6/Je1rXsp16Tf7InPGN5m7HnOVbjqXUMj9tqbJb5mm0zAFqOe5Ur/HB5kyba04+45L3gCXHeNtXvWW/NuVlXwn2NTZhi+UBDpU3ZQAAAB0JygAAADoSlAEAAHTUfJyypXIeasa/WrfkOGXb1lUzbtucscVqbNvfljlkNfWj5XhPNevd15he6+bU40M9HkseyyW3q1X+3pJj0bWqH5cuXcr58+fnbCYzbcrLvhLH5ZrzPaAmr/JKOR4Au/CmDAAAoCNBGQAAQEeCMgAAgI6aj1NWY07uxdwcs03LzrFkDkzNeE+14x212o65WuYHLJWfVrNs7Xhyc8qaMx7YodbbOWpy7lqel7llt7yON9X5fY7HSJ1N45TNqcdXwlh1tTnLV9vxADjmTRkAAEBHgjIAAICOmneJP0dNk7ptzRKWKqtnc7Q5667pAr2m2VNtF9+bylqyudWSXbGfdjsOpcnuXHO6w64ZMuFQus/fVtaSllrX6v4fHR3l3Llzi6yH0aYu8Te5UprnLdXEWxNcDtWhPI+5snhTBgAA0JGgDAAAoCNBGQAAQEd77RJ/yTynbeuqKWvOsjX5AHNyzGq6zt5W1j7b+Ncs2yvvZ86xP5S240sOmTCnrJpz1nKohiWHBKjZjprzUrMuuTlXpiuxC/iWQ5/AoWiZ88+1y5syAACAjgRlAAAAHQnKAAAAOmqSU7ZJrzFXWo5TNmespJZ5Lfsaw6smt6123ZvKXXL5mvyalvV007pqx8SbU/aSY4Atld85J29ybllzlj3Ue8K26fsaA446NbnEh5KrUvMMnXOvvlKOB4ehZX64+2sfV9s9wJsyAACAjgRlAAAAHQnKAAAAOlo8p6yXa63Nbs04buvL9hyLbSlzcnVajom3rT1zTTv0lmONHcr10jIHdc7xqD3Wp513m321j99U7qVLl3L+/Pkm62G+lnmnh6Amb3Lb9CvxeLA/LevWvnKWD0mvPiJqxoDbNv3QzoU3ZQAAAB0JygAAADpq0nzx7NmzuXjxYpLlXunWNANbnz636dac7tLnWLJp0pLdcrfs+rVlc71DabrSqr7MPYc118+hNMfbNP/cLvB7DRGwZBPMGqfdh6Ojo5w7d67ZdrBZr2uzxpJN2ufc12DVksPKzC3rStDyOt7Xd9Dac3xovCkDAADoSFAGAADQkaAMAACgoyY5ZRcuXMiZM2dOnFbTLffqsnPbgda0y68xJzelJp+kdpvndA8+p6y57YhX56/tyrTV8WrZ9fxcS+b+tVq2ZU5IS0vml9TkGtTk+m3Tap8PJTeJeZY8T0uVXVvuoW7XJoeex9LCvr5jzVlvz7KuBC1z3ve5HVfzefKmDAAAoCNBGQAAQEeCMgAAgI6a5JStmtOuuGYcoW1lz5k+J69jyZyPmrHWaspat89xyvY1dsU+z/GcXMglz0vNvEuNVzL3HtAqB3N9+dq8wZpxUpYaL+1Qc/2oU1NvN02rGRdyn8uuuxKPx74sme88xz7zsDete8lzfCXqeQ+Yo+a78JV+nrwpAwAA6EhQBgAA0JGgDAAAoKMmOWVnz57NxYsXk9TlxMyxz/F8atoVLzU+1pxcnG3r2pbnM8c+89GWytebW1ZLNeelZpyYlmPMbDrWc3K11tUsu778PtvDL6kmt23Xci9dupTz5883K5t2au6Rh7Jsy2f5oezT3OlLPXMPdb1L5gG1rHtXen5S0q/O12j53e9K400ZAABAR4IyAACAjgRlAAAAHTXJKbtw4ULOnDmzdb595oFtWner/LL1cluX1XJshqXGRppryfbvNbl/+6prvfL1asra5zgg+7o2D2WMopbH8mrIf+COWo5N12vZper53PtFq2XnLt/ymXqtrXeuOef4ahvjaq6e48nVmDMm7aZlD5E3ZQAAAB0JygAAADpq0nxxVz1fK7ZqblSz3nW9us5O6roHX6o5zbay5pjbffqcpo/7atpWe73MaYK4qQnAPutHy3q6SW093bSuJZsSt2ySucmhN/FgtK/maNtcKfWlV0pEzRA1NXp1HT53vfts5r+v7dhU7rpe5/xQupZfcjikbcse2r3LmzIAAICOBGUAAAAdCcoAAAA6apJTdvbs2Vy8eDFJ2/yaOfkTLZetaZM7J3dpm5ZdnG8qu2X34PvMz6vJ66mdf1c1Qxf07C69Zb5Eq7Jqh4FYKkdx7npqrr05Dr3tPLtpda++FupDy+PR8lm9ZNfj+7qfbFrvetkt85DnbkfNemtyzTet+1CHCDjUe8CS13HLc7wEb8oAAAA6EpQBAAB0JCgDAADoqElO2YULF3LmzJkkbcfCaVlWTbvjTTaN7zR3vXPUtG/etu6aHJm59jVuTM05nTv/pu1omX+15Dhuhzh+Sc21tj5/bdvxpY7P3HPeajs2lXN0dJRz5841WQ9tLTmez5WoZc5qSy3HKWuZa77UOHeHWi/3mdNcs959jVt2teSY7Tpt7rI9joc3ZQAAAB0JygAAADoSlAEAAHS0+Dhlm9S032w5DsaS+TRL5dgtOQ5KTR7U3LJb5vq1HF9u07R9ji1VY06d36T2uC81XkvLXMdeYxS1HN9oqXGXLl26lPPnz++8LMtpeQ+8Fs3J2V1qvdvWfSjjHvZatqXasXH3lX/UK5frUHPItmk5hvGh5Zh5UwYAANCRoAwAAKCj5l3iH6p9dYm/ZLPJfXWT2tKSzfHWj0dNd/I1r7DnNCnbth27btMu89fU+ZplWzYHqGmSum5OM5+a5sD7bB6xVNPQ1WV1iX+4NFecZ1/Hq1fzvZ5N6Zf6ftKy3G3H51C66l9yCJ9Wy67b5/fTQxmGaQnelAEAAHQkKAMAAOhIUAYAANBRk5yyVXNyqGq6ep3TDfcu65pT9qZyl+zmv2X3+nPya+ZYMmdoySET5uSB1eQXbVPTjf1S7azn5lttWn5bnV4yn3OptuS1+WlzylrX6po4lC6smadlfuuVaO7zZanu47fpdZ72OVzHUtu1z27Jl7pHXinHslf9WDJPveU+7YM3ZQAAAB0JygAAADoSlAEAAHTUPKdsqfyJ2raeNbk6c8bMqCmr5ThL2yw5dtIhmJsXtupK2d+W45LVqMkbXLKt/RxL5oLua0yVluPVrG7zpUuXcv78+VOXzXJq7nNXw31+Xcvj0dKhnKea7y/7HIfqEM/TkuOBLXkt1uzTIZ6HbdMP5Vo7LW/KAAAAOhKUAQAAdCQoAwAA6GjxccqWGhekZiyCdT3bxp62rNq2rjXjcm0qa5/5NUuNN9Fy3IttWo6ft2m7lsxfrClr07zr9tnee5/1Zc565xz7VuNAHh0d5dy5cxvLoo993XuuFr3GbVvyu86meQ/Fks+MpZbdtvw+87JryqrJ6Z6znhq96sf68odw/XhTBgAA0JGgDAAAoCNBGQAAQEdNcsrOnj2bixcvtijqsmrHKdtU1hwt26/WlrVpu/aZr3faebetd1vZc/IX56yrdgy8OWUf6nlp2c66VT7n3PrR8nqpaYc+Z9lt27FUjuamcoxTdnU6hDF5Wpt7/VwJWua3XimWHBPstMvOPZa98hc3bcfceXttd8vruHZs4X3zpgwAAKAjQRkAAEBHzbvE3/QqsKbby5r1rqtZ79wmVJuW79X99braJlI1Td3mdL2+adk569m2rtr11pynfTVnrFnvPoc5aGlOM8qW19OcZWvmnbsdm+Y9tCYdnE7N/fVqtNQQPTXbsa6mafWV0jyzZp/2tWxLV+q1dijbfSjbsQ/elAEAAHQkKAMAAOhIUAYAANBRk5yyCxcu5MyZM9Xl1HQ7vW5fXUfXtFmuybHbNm/L9vLbytpX3kLLLuHnaNlOv2Xe17bu1DdNm1M/5ta1TeZs8zY1eWC118+u69lWdu39ZI7T1rWjo6OcO3fu1Ouln5Z5lFeiOdfPoR6POfeqQ9nmbZZ8puw6bZfpsG/elAEAAHQkKAMAAOhIUAYAANBR83HKasxp0z0312JOHseu5exiX+NetMxF2eeyLcd7mjM+2JzcgpZaHtttWo1rt2SeU8vcrX3lW61rORZOrxyHK2V8I9rpNWZTT1fCPl0J27i0Tc+fa6Gecu3ypgwAAKAjQRkAAEBHgjIAAICOFs8p21cu05yyl8y/2qYmz6dmPJI5uVv7HCupJg+s5liuazkm3qbt2Od4aa3a2teO9bJpzLOa+rBt+pLnpSYXcqnct6Xua5cuXcr58+d3XpYrgzxCDpW6ybXKmzIAAICOBGUAAAAdCcoAAAA6apJTdvbs2Vy8ePHEaa3GoWqZ81EzPsucvKa527FtXXPWu+SYb5vM3ac553jJsaTmLHso610y526pOtAyL2zdUjmo29a1zzHOahjP5+pXk3MI+zLnO5h6y7XEmzIAAICOBGUAAAAdNWm+eOHChZw5cybJst2Wb1p2yS5U5zRVWrIJ3pxy5jRvXLIp2zZzmmj26l7+UIZMmKvV9VXbjHKOmia7Nc2Sa9c1R6umj0ta3cajo6OcO3eu49bQwqHWNfZnX99PardDl/hcq7wpAwAA6EhQBgAA0JGgDAAAoKMmOWWbbMpNmdO+ecm2zkt2Cb+vLuCXzEWq6ZZ8vayaLpv3lQe2bbv2tew+8wbnlr1p3n0OsTDHnONRM9zANnO6e26ZJ7dO3gaw6lDyXWuGQ4IrmTdlAAAAHQnKAAAAOhKUAQAAdNQkp+zs2bO5ePFiks05EIc6Tsqc3K65+TP7HE9tKTW5TXOOx9w8n33pNY7bvsbK2lZWbR1eKj9gbn7VnH3cNn3O+Iub8tFqc8ZqxnXbZHW9ly5dyvnz509dFtDHoXzfmLsdh7LdsG/elAEAAHQkKAMAAOhIUAYAANBRk5yyCxcu5MyZMydOq8knOcQximrHe2qVn9ZyvKJ1S+YM1YzbNmc79mnO2Gstl60Zt63XWGItx8CbO33TvHNzt+aMv9gyH22Tmvy0TftwdHSUc+fO7bwdwGE41Jz3luO5wtXEmzIAAICOBGUAAAAdde0Sv+a1dE1TppZdaS/ZjLBGTdOlXk2oautHqyZ4c5ujnXY9J5U9Z9l1S9X52mEgdl3PIau5r9V0W7+v4Rg2rUeX+HBl2uf3pFbrXZ9+pTwjoAVvygAAADoSlAEAAHQkKAMAAOhor13iz5m2Tcu8ljlqunxf13IfeubJbWr/XdMd+DYthy6Ys54ac3KIanPb5pS1aTtr93+pLo33lW+1pNohNgAuZ8kc5nUt+w+QR8a1ypsyAACAjgRlAAAAHQnKAAAAOlp8nLKl1IyDseQ4XNvU5Oq0HNetJtdvX2PAzR3/aVMO1Zz1zlVT11qO6bVU/kDLPKaa49Eyn3ObVmPebSu7dpuXGk9udd6jo6OcO3du1nYBV7Ylx5EFTuZNGQAAQEeCMgAAgI4EZQAAAB01ySlb1TLPZ6kxNebmXsxZdtv8m7TML5pjyTFE5hzrmvyzk9ZVs11zyt20T0vWl33lUM3d5l51vmY8vZqy5kybu55tWt6rgGvLnOdxzbPcvQh2400ZAABAR4IyAACAjgRlAAAAHTXPKeuVF7VuUxvmmrbSLe1zLI+aXK6aMeFq5p1bVqucqiXP/77GmttmyTG8WpZdk7/YUsuxDFvmZM5x2n24dOlSzp8/f+r1AoehJn9+znTjlMHpeFMGAADQkaAMAACgo+bNF9e1aha2ZFPHmu6xW3bzv69uuNfXta2JVMvj06sr3JZd4G+zr2Eg1i3V9G9ufaip46ctd5eyNw1VMNecczxHrybNm4Y1ODo6yrlz5/a2XcB+LNU8XJf3cDrelAEAAHQkKAMAAOhIUAYAANBRk5yyCxcu5MyZMydOW6pt8ZJdZ7fMA5szf03e07668d9lW+bMOye3rWXX/EuqyaPsdb20rD9L5cnNzWXbNP+S+7vNUrkZLc+xLvGBVfLEYHnelAEAAHQkKAMAAOhIUAYAANBRk5yys2fP5uLFi0n6jUU2J99oybbRNXlRc7ardh9ajR+3zbayW+ZftSprbr7iUrlsLcflapn7uGk925ade2yXGk+v9pwtdWznbleruidfBAD68qYMAACgI0EZAABAR4IyAACAjpqPU1aTX1KTb7Uth6gmZ6iXfeV2XSnHY12r8a+S2+9z7f63yl2qzVXadI5r1j13vLhWuaK1uX1zxi2bU3bLPNI597GT/r3rerctCyxnX7n2wJXFmzIAAICOBGUAAAAdLd4l/qqaZj77bGLXshlljV7DC6yrOfZztrllF/Atu9PftJ5d5t91vevTa/Zh23bNuRZrm9jNseT1s0ntsT7tsrXd/M+hGRQAHC5vygAAADoSlAEAAHQkKAMAAOioa5f4c9R2D16zbMvusDflDC3Z/fW6lrl/NTlVNZYa9mDueajJm6vppn3uujaVva+czZbXS8tct7nntOW1uGnZJe9jV+IwIXAlOpR7wqHY51AxcCXxpgwAAKAjQRkAAEBHgjIAAICOmo9Ttm7O+E9L5Wnssu5d11XbvrtmDLQ521EzJtzc/Jo553jXcnbZjjlja9XWnznmjDW2yZK5fetqxm1b16qObyurJgdzyXy9beuqqac149ptKkueBiznWs/vbJlLDFczb8oAAAA6EpQBAAB0JCgDAADoaPFxylYdytgUc3PbltqnuVrl561P77UPS5ZVk1M3Zz1zbTsvc/KNanLMavKx1rXcjjnrmWtOrluv62fJPMqeeZbAbZbMd+2l5TN1ybx+OGTelAEAAHQkKAMAAOhIUAYAANBR83HKasbDmqOmXXHLMTPm5Ai13M65+SBL5pPU5B9tsm0ba8bWarkdm+afWz82WbLtfM11umQb/5b3iDnlLjkm3pzruOV9bc6yq46OjnLu3LmdywKufi3HOt20vBwyriXelAEAAHQkKAMAAOioSfPFVS275T7tvLssv0mv7mrn7GPPZnH76lp7n83CroTmnOtaNuto2X1+y+a/m8rdpqYr/pbNX5csZ9OxrbleVqddunQp58+fr9lMYEctmzS3dKjDyMDVxJsyAACAjgRlAAAAHQnKAAAAOmqSU3bhwoWcOXNm63y13cnPKasmX22T2vyrGjXttGtyeWryr5bq0rzWkt35zukSv6Wa89LyeqkZqmDTdsw9Ly3rT8uhPWrqXs2QI5tyzjbNq0t86KfmXr2kmmFC5InBybwpAwAA6EhQBgAA0JGgDAAAoKOu45TVWDL3Yt3qsrV5KnPG0tq07LqWx73ldu1Tq9yl2tzHmrG39nUsl8xVamnTOe05ZuCme8KS49S1zG0zThkchjn58/tSk7MrhwxOx5syAACAjgRlAAAAHQnKAAAAOmqSU3b27NlcvHgxyXJjFNXMuz5903g9c9WMtbbPfJI5x2ufOUOHoqaeLmXJdvktz/Gca69GTR7c+vy1x7blOG41y87Jd5XXAVeeXs/jluOTyjGD3XhTBgAA0JGgDAAAoKOD7hK/pjvsmlftc7ZjbtmbyqppflXbjXuNVk00Ww43sK3spZrGnrT8HDX1Y5OWQ0jMNWeflmrat8vyu27HXDVNumu6178Wmh3D1W5O1/SH0ixwbnf6h7Ld0Js3ZQAAAB0JygAAADoSlAEAAHTUPKespm1wTf7RnJyzuflYNflHm6bXdGG9ZBfeLbvX32apbv5r2rTXtm9v1fX6knlfLdvw11w/LXOm9pl/VbNPc+rpuk3zL5VzeHR0lHPnzp26bGCzOc+MXvlYNbn3h7pPcGi8KQMAAOhIUAYAANCRoAwAAKCj5jllm7QcS6vluEs1eWD7HEuqZr01YyVt02ocqiVz2VqORddyvTXbUbuuXcuaO95Xq3Hr5loyL2Gp8QZbbtdS99dLly7l/Pnzp9xCoMaVMk5ZTU73punyy7iWeFMGAADQkaAMAACgI0EZAABAR81zynrlk9TkdrXME1tyO+bMOzcP6LTzrq+7Znyw2lycOWM4bRpLq7Y+tBqnbN0+62nLcdtqxtaqqR81+Yxzrp+aXK6W+WhL5egapwz6mXOPWDKvdl/fg9aXN4YZ1xJvygAAADoSlAEAAHQkKAMAAOioSU7ZhQsXcubMmdnL9RrDqibvaW5ZLXOI5pRbk18yt6w5y27K1ZnbDr0mP22TJXOots1b05Z+znYvOa7fHPvMeVgq33WfOYdLjl0IHJ6a+2uv9dZuc699ht68KQMAAOhIUAYAANDRXrvE77meTc3ktpW9z65gdy27ZffotVp18V3Tbe6caSeZUz+2qRnKYI591YGW3Ruva9lt/b6uvXUtm2/WNDlsWbdWl7106VLOnz9/6rKA3dWkaix5z5tzv629Vy+1T3DovCkDAADoSFAGAADQkaAMAACgo+Y5ZetatQ2em2uxr+7St+WfteymfE65LbvE36ami+9N89bk9m0re9P02q7oN7Wtb9kNec2xbT3/rlrmPW0re92mc1yT81CjJjd2fTtqcuo2LXt0dJRz587N2UxgT1rmQ592vXOmzZ1u6A6uJd6UAQAAdCQoAwAA6EhQBgAA0FGTnLKzZ8/m4sWLSerGGdo0/5JjErXM+5qTB1UzHsm2eeesa27776XGJ6nJ3do27xz7yic6Sasx3w5JzRiBvcbgqcmJmFP2kues1X3NOGVwOHrd51uOqdmybLiaeFMGAADQkaAMAACgI0EZAABAR83HKVtyLLIac8puORbQpvnn5iq1Gndo27qWzN9bNydnqGW79JbjyW3ajtoxzzZNa3keWubNLZmjOWfemjzKllpux1L7JKcDWIr7C+zGmzIAAICOBGUAAAAdNWm+eOHChZw5cyZJ2y6tezWxm9Nt+5Jdr2+av7b5Vc0+7atZ2NzmnK26yK9tyrdpO1o215uz7D679d+kdviF0867Pn/L4QZq9qm2yXKrJs2HUj8A4FrlTRkAAEBHgjIAAICOBGUAAAAdNckpO3v2bC5evDh7uW15HZtyQNa17C5907Itcy/m5hctlSdXm+tXk6sz5xzP2a6aoQqWzE9ct2QOVauyavM3Ww2TMTenbk59mnvtzRkmY6m8uHVLLgsALMubMgAAgI4EZQAAAB0JygAAADpafJyyOTlDLcfRqdmOTWXVjm+01Lht27Qaz2huWXPWtWTu0j7HYWqZyzSn3Jrx9WryoOZMr80NrTGnrrUc125f+Xs195NNyx4dHeXcuXOzygYA5vGmDAAAoCNBGQAAQEeCMgAAgI6a5JSt2ucYT3PW03KsoFb2OcZZTY7Qksej5ThlrbTMX1zXcmytbdvRKnepdpuXugfsa3y0k7ajVb1tmb/ZcjtWXbp0KefPnz912QDAdt6UAQAAdCQoAwAA6EhQBgAA0FGTnLKzZ8/m4sWLSeaN99My12ROXkvLMYm2bdccV0qez6ay5u5/yzyyVmNL1Y7LtlT+YsuxteaOeTZn2TnHtua6nWvOOa4ZH6zl/WSbmnH+jFPGvtQ8u1reT4DLcz31500ZAABAR4IyAACAjpo0X7xw4ULOnDlTXc4+m/1sKqum6dZS65277pZN3WqabC65npqu6Detu7Z5zRz7avo2t97WHMuW189S12LLczh3XbtOa7memmV1ic9cc5r7tpyu+RVc3tznnOtp/7wpAwAA6EhQBgAA0JGgDAAAoKPFu8RftWRuTk1ZLbuxn1P2Pte7ZM7MHJvqQE1X6+vT5+Y9rS5bm/PQSk1X0evTl2wbPqfsubmOLbvE31f3+jVqrut95g3CadU+m5a69uBac6jfE69l3pQBAAB0JCgDAADoSFAGAADQUfNxymryS1qqyRmq2ea5eT+nnXfJtsA1+zDXnPHBNi27vnxtPtZp17s+vSa3rWbZXZbf1ZK5SkveA+bUjznL7rL8aeedu55W+TVyyKix1PWwxPLAqGbcMs+MZXhTBgAA0JGgDAAAoCNBGQAAQEeLj1O2VJ7Derlzxz/atOwcNWNHbbPkWGtzclHm7OPcvKaWeWBz8tP2mZO3yab1zj0vS9XxmnGE9lnWnOMxZz0nrWvTsjV5pDX1VBt/DsWS9+aa5w1wm57fkziZN2UAAAAdCcoAAAA6atJ8cVVNN9Q1TXXmNAOq7Q571208SU1zvV236TTbNUdNV/1zlq2ZPvccb1pPy/O0ZHOAXk0wt5nTjW7LJrtLNTmcs55dpu9r2ZrhF2COOc2ya4dCmVM2XMtaXmssw5syAACAjgRlAAAAHQnKAAAAOmqSU3bhwoWcOXMmyfbcrlbmltsqV6V2vatlt+xKe5uaHJk5uTtz83xa5u7UaJWPtq2smhyznsMezFl2zna07Na/pqwl93+pYUHWy14qP/Ho6Cjnzp1rVjbXttq8Wl3iQxu6vD883pQBAAB0JCgDAADoSFAGAADQUfNxylqOhTMnJ+a05bbWso3ukmNY1YwXNmcMuJa5bDW5bTX7X7NPc8taMvdvzrytxtPb5lByt+bmq9Wsd185MbXj/sESap9rxtdjkyvlu+KVTj7nMrwpAwAA6EhQBgAA0JGgDAAAoKPmOWU1+UY1+SVzyq7Zjm3zzhnjqmcb3Fa5futltTS3zXKrPLGacdq2lTVH7dhircaA6zmWyVJjre1zn1rm2e5rnL/Vci9dupTz58/vvCzMqfNzn8dyV1jV8jtny3v11aDldcxuvCkDAADoSFAGAADQkaAMAACgoyY5ZWfPns3Fixe3zjc3j2OpfJJ1NfloLbejpg1uz/a8m87TnONVk7u1zT7zoPa5rlU1OXU169lXLuhcLcdtaznOXc16Wo5jt2m7VqcdHR3l3LlzO68Xdq1bu0yXm9JOzXiMh2Ju/Vlq2SvhWNWqOZacjjdlAAAAHQnKAAAAOmrSfPHChQs5c+ZMkuVe/9Y0fVyff8lmcnO2o2VTgpqutffZRfGmeWubbl2J3dfWNIvbVlaNlsdvTjPkJYfUqBmq4LRN/+Zacp80N+EQ6Er7cF2Jz9B1ve7zV+rxWrXk8FBXw/HZB2/KAAAAOhKUAQAAdCQoAwAA6Oigu8SfU9ac+WvWW5sH1qob95bHsmY7tmmZN1hzbJfspnyOJfMEW7YHr7GvvLi55+VQcqqWalu/z+EGYNW+ci6vFnPuzUut52pVc2xbfje82nOorsW6tQ/elAEAAHQkKAMAAOhIUAYAANBRk5yyTTaNe1Ez7tQ+89Fqxu6oGXdpXU3O1JzxwWrU5LXUHttN69k2vaautcxXm1Pukvl6+2pbX7Md+8whm3Mdr2u5T3Pqbavcx0uXLuX8+fOnLotrT8v767Wg5nvBruXuUvaVcC7m7FPLccquxfH1HI/986YMAACgI0EZAABAR4IyAACAjprklF24cCFnzpxJ0nZ8kjntrFvmsWzajtqyNmmVi3TSdmzazpZtg2uObe05Xp0+9/jsa+yofWo1btuSOXYtx/Wbsx37zBtsqeb+etr6cHR0lHPnzu24hVBXT6/G3JS5+3wI6z3U83AlHstDOXYt1Z6Hq/GYtOBNGQAAQEeCMgAAgI4EZQAAAB01H6es5XhYc+ZvOWbRnLbBPcdPa7Xskm38a9oVL3ksW9atOfa5TzXHtibXb864MTXbsWnebVrmZG4ru+Z41IzJY5wyDtGVkJ95kppc2V77XLPeKzXnp+b+OqesTfPuMv+VaF/j3V7LvCkDAADoSFAGAADQUZPmi2fPns3Fixe3zjf3dfimLqw3zTt3XUs2z9tkyeYBvZrJzX2F3/KVd83QBXO6S5+zj/tsvtmym/uluo/fZ/2o6aK4pr4s2WSopuwrsftrrkw1w5WsO5RmUS2HrKlZ776+n1ypXZofYnNY99drc59Pw5syAACAjgRlAAAAHQnKAAAAOtprl/g15rbJ3VcuytztaJWrU3ucl8rX22f3+TW5Ovtq471kPV2yO9+ltqOm7Nr8s5bHo1X+4lw1eXLXQpfN9DHnflqbs1szXEerZden73OfNtnnPvVyCMOE7FL2nO9YV4NrYR/3wZsyAACAjgRlAAAAHQnKAAAAOmqeU7auVQ7Eptysk8pq2ZZ812knlTVn+pz8oyXbe9ccn0Nph95yrLHacbpOux0tj2XL9t692ta3PA9L1tOaurVuX/eqQ8kX4cpQk99Zmxu6VNktl60Zk7Vl2fs8HvvKKVpyn2r2odd6ubp4UwYAANCRoAwAAKAjQRkAAEBHi+eUtRrTqjYHpNWYVi3HnaodU2STmu1cql117bqWHP+pZgyrpcadalnHW7ZZ32c+51LHdu6yLfPElrLPesu1bckxI+csv+Syh5LnU/O9YJ/n6VDG5eqVK+v+SgvelAEAAHQkKAMAAOioefPFpZoMLdlldcv5l2yONaecJbttX6rL85bNr+Y2/Vyq2em2Y9dqPdvKXrJZ6ZLlLtWstGUTqnU1XdHX1JdN65k77+p2HR0d5dy5c6feDlhV+wxcaoiNfWo5tMUmh7rsksd+X80Gl0xjOJSms+yfN2UAAAAdCcoAAAA6EpQBAAB01CSn7MKFCzlz5syJ01q1s63NCdpXDsg2h5KrM0dN7tI+u4k9lOEXapZdXXftOa3p/njOvPtqDz/3HrCvdvlLDs+xzaZ7U6s8yUuXLuX8+fOn3EKuBS1ziVsOsXGlqBk2ZNdydyl7k5pc6iX7BGg5BMtS2zX3PFwt9XrV1bhPS/CmDAAAoCNBGQAAQEeCMgAAgI6aj1O2rtXYSS3HvajJAVmyjfa6fY351XL8p23bsa/js+R658w/Nw+qV7v8lmq2c19jvyyZB3aljCe363YYp4x9qh2D82pTczyWPFY1610yv6jX8diktk4vlR/eUsvv1dcyb8oAAAA6EpQBAAB0JCgDAADoqHlO2VJjRuyzzemm7drnGCLr5rTBnZPLVDsG3JztaFk/WtaJlvlpm+pLy/UumedUU9eWGgvnUMac2WbJ+9y+7pHGKWNf5l5rc8bmuxJzVVo+j1veI2u+UxzKs3vOsrX2lbd/pY5xdqjb1Zs3ZQAAAB0JygAAADoSlAEAAHTUPKdsTtvYljkgNW2Ha/KelrSvMb5q2zu3GndpX+3Mk+XyetaXX7KO1+YVnna929az1FhjtbltNXmUNdfekveTVjlmV0oeAlxr45TV2Od1Pec71r4cyn1szrNprkMde23O8odyng6BN2UAAAAdCcoAAAA6at58cV1NU52WrzTnNKOsaYJZ85q6V9f0NV0Sbytr27I16527/KayWg57MOccL9lsoabJ7pxya5rI9BrqYu5xn3NOW3YdPUfNsdy0HUdHRzl37typy4Y5Wg6jcjW6FvbxaldzDg+lztcMs3PSdEbelAEAAHQkKAMAAOhIUAYAANBR85yylrkpLXMvWuW29eqWvLbclu15r4QuiVvm12yb3vJ4zKlrNUMTLJn3VZNT1zLnsNUwGNvWveTQFeuWqnub9uHSpUs5f/58k/VAraWuRbgSHGqd7vVd+GrjTRkAAEBHgjIAAICOBGUAAAAdNckpO3v2bC5evJikLhejxpz1zh3vqmWuzj5yQGrXs2TeT8sxvOaMX7NNTR7QUuPatTwe+8zlWmrMs5oxALetd8k8wTnmtrtfaswm45RxqObc5+WYAVcSb8oAAAA6EpQBAAB0JCgDAADoqElO2YULF3LmzJnZy83JJzmUcbZq89Fa5S7V5NO0tikPas6yLeetzceaU/am6b3G0zvN8pdbdm4u175y7JYca2yOpY77LmXPyblbt2nZVcYpo6dDyfcEWJo3ZQAAAB0JygAAADpavEv8muZYcyzZ/Khlc6RN621Z1pJd5i/VzHRu1+tzylpy2ZZNH+csO7cZ4SZLNuldqjnnkt1bt2z6WDNcxz6v403rgUMx556gqSNwJfGmDAAAoCNBGQAAQEeCMgAAgI6ad4m/rxyqlm3Fl+xuf1NZc7vPrzkeLfOeeuWj7SunqHb/5+RQzenifNN6TiprqZyhmuNTM1RB7TAHm7a7Jper5fHYplc+JxyKfX0vANg3b8oAAAA6EpQBAAB0VNt88bokufXWW9/4wdHR0WVnvnTp0sbCNi1ba3Xd29azbTs3WS97U1lz5j1p/qXM3Y5Nx3a9rKWWnbONtebUn5b7VLtdc7ajZpvmlDXn+MzZ5l3WvUnNdm0ra842taxrc9azuuzK/f26nQtkV3d4hrJZzXOw5XMAYBdznqGlJrfkmc985n2TPPfUBQBwpbjfAx7wgOf13oiriWcowDVj6zO09k3ZC5LcL8nLkryhsiwADs91Se6R8X5PW56hAFe3nZ+hVW/KAAAAqKOjDwAAgI4EZQAAAB0JygAAADoSlAEAAHQkKAMAAOhIUAYAANCRoAwAAKAjQRkAAEBHgjIAAICOBGUAAAAdCcoAAAA6EpQBAAB0JCgDAADoSFAGAADQkaAMAACgI0EZAABAR4IyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEdNgrJSyuNLKcNl/v59i3UsoZRyr2nb37Fhmfed9vv+l5n+gdP0//sy0+9USnlJKeXJO67v/lN59z39VrdRSnlwKeU3SimXSil/XEp59CnKePy0/HucMO0LSilvKKXcr80W91FKuWE6Z0e9t2VdKeUepZRvL6U8p5RyaynlNRvmva6U8thSyh9M876klPJda/O8eSnl60spf1JK+btSyh9N5/j65fcGrgyeobcr85p7hpZS3rSU8m9KKb9ZSnlNKeVPSylf1mr6jtvwUdNx+IwTpr1vKeV1c8s8NFN9HUopD++9LetKKW9RSvmfpZQXlVJeW0p5WSnlZ9a/7+z6jJ6uo/9eSvm9Uso/lFKetsM2fNJ0fF64Zb5vn+Z7wo779p6llJ8tpfxtKeWl0/Jvusuy15o7NyzrtUk++oTPLzRcR2v3SvJVSZ6W5CX7WOEwDM8vpfxBkk9N8oMnzPJRSd4hyZP2sT2tlFI+NMlPJPmBJF+a5MOT/NdSyt8Pw/C9M4r6T0kemeTGJP9spfy3S/K1SW4chuG5zTacde+U8fj/WpLnJXm/DfN+X5IHJvkPSX43ydsnWf9i851JPjHJVyR5YZIPmua/e5IvarjdcKXzDN3BVfoMfa8kX5DkGzPeSx+a5JtLKX81DMOTGkzfahiGXyylPHFa7mnDMNySJKWUkuS7Mt6//3OrHeYO7pLxHvD4JH+W5G5JvjjJs0opHzgMwx9O8+36jP7wJP80yXOSbA2ApiDpW5P89Zb53ifJ/5PkVdvKnOa/W5JnJXlxkk9O8rbTet4myYk/rFzLWgZl/zAMw682LO9q9uQkX1lKORqG4eLatE9LcjHJz+9/s6o8LsmvD8Pw2dO/f7GU8o+SfE0p5fuHYfiHXQoZhuHWUsrnT8t/2jAMx792fmuS12T8cr+YUsqbDsPw2iXXceB+exiGt0vGX+9zmRt+KeWBST49yfsNw/C7K5N+aGWeOyf5lCTfOAzDf50+/sVSyj2TPCKCMljlGbq7q+0Z+sdJ3nvl2fOL0z324zIGl7XTd/VlSR6S5OuTfO702eck+eAkHzIMwxtOtXc7KKXcJcnrd/2ucLUZhuHlSR61+lkp5ReSvDzJw5N83fTxTs/oJP91GIZvn+Z79g6b8O8yBoMvyh1/XF31hIzfxz5zhzKT5POSvHWS9z++Vkspr0/ypFLK1w7D8Hs7lnNN2FtOWSnlE6fXnQ9d+exupZQLpZQfWflsKKV8eSnlG6fXt68updxUSnmLtfLuVkq5sZTyl9Mr3OdPN6H19T6klHK+jE2nXlFKeXYp5QPK2DTiF6fZnjutdzhF+f++lPJXU5OBp2b8FWCbJ2UMiG/3Cr2MTbr+eZIfGYbhddO2/8L0uvdV0+vqB20q+HKv50spTyil3Lz22TuXUn6wlHJxel3+S6WUD9xh+9fXeX3GX3jPnbCf75DkA+aUNwzDs5PclORbp/Pw0RkftF80DMOrSil3L6V871Q/XltK+bX1c1NKuXn91Xop5eHTsbnX9O/jY3VDKeV7SikvT3LiW7ipDr7ghM8fNpVxn+nfn1FK+eVSyi0r9e2DNu1vuUzzmVLK09ZvpqWUe5dSfqKU8soyNgX46VLKu20qf44ZD8TPTvKstYBsXclYz1+59vnfTNOAHXmG3s5V9QwdhuE1qz8GllLePGPLg1e0mD5jO16e5DFJPqeU8mFlbGL/9Um+YxiG55axaej/V8YmdreWsTn6F6+WMdW1F659dnT8rF357ObpmD62lPLijG+J3mZ9m6bn8+vL2Fpm9fO7l1L+vkxpEqWUDy2l/GQZm67+bRmbcj5qvbwTyh9KKY9Z++wxq3V5+myn+tzY3ya5lORNjj+Y8QP3zsHt9B3iy7Llh9JSyqcneZck37Br2UkenOQZaz+e/GiSW6dprGgalJVS7nzCX0mSYRh+POMX7e8tt+XS3JjkuiSfv1bUFya5d8ZI/MszvvL8npX13CXJL2R8Rf8VST4+4yv7ny7jq9Xj+R6R5KeSvDTjl/pPT3I+4+vfX8/4uj9JPivJh05/c8r/VxmbYv2PaRtflOS/bTtOwzD8ScZXyp+2NunBGV9ZH/+y9S7T9j9qKv98kqeXy7S1n6OU8tZJfjnJ+2c83p+c8QbwrFLK267Md5zrcK8Nxb1bxlfv6794HH9hv/cpNvExGevGN2esJz85DMNTSynXJfmZJJ+U8dx8csbX7U8vpXzUKdaTjE0mh4zNYR57mXmenOSflFL+ydrnj8z4y9Xxvt4rYxPOT8l4fi8k+aVSynuectveqJTyrkn+T8amfzdM5d8jyTPLSo7W9FAcTiyknQ9J8odlbBv+yukLydNKKe9yPMMwDK9L8t+TfGEp5YNLKXedztG/zPhrG7DCM/SafYaulntdxuNxp4zPv6bTtxmG4YkZg+3vTPJtGYOl4xYq35TxfP1gkocl+fEk31ZK+cq565l8csY68q8zNnP/uxPmeWqS12V8pq4vW5Ic/yBxz4zn93OmbfvRJN9XTsiRm2tGfb6hbMiH3HFdd5qu+3dI8i1J/iHj+VzStyf5gWEYfmvDdr1FxvP/2GEYTjpPl3PvrH03HIbh1iR/ktN9N7y6DcNQ/ZexDexwmb/7r8z3lkluzniRPWKa/nFrZQ1J/jTJdSuffXbGivne078/K+NFep+1ZZ+T5Ien/y8ZvxD/7Ibtvv+0vvuufb5L+dcl+YuMFXl1niev7/dl1v1F0z6dXfnshzM+lMoJ898p4y+DP5fkyZfbh4xBwZDk4WvLPyHJzSv//uqMbyzeduWz66dj9o0nnNt7bdiXD5/m+ZC1z+88ff5Fp6xXnzEt/+rj45TxZjgkefDasfndJM9e+ezmJE9YK+/hq/uycqx+eodtuS5j8Pe1K5+96bRtX36ZZY7P2e8n+bqVz2+Y1nu0pR4+bW2fnpjx2jiz8tk9MjbrfPTKZzclGRpd16+5zLRbp33/lYxfhD4lYxOa30ty57Xj9j25/T3hv9Rumz9/V9NfPENX57nmnqFr67kxY3DyEUtM33Eb3iPjG5ohySdNnx0l+fvVfZs+/+6Mz6C7Tv++KckL1+Y5msq6YeWzm5O8LMmb7bA9P5rk/Npnz7pc3cxtrTS+O8n/Wfn8Dud2+vdj1pZ/TFaeobvU5+nfN+xSd7fs63/Mbdf+Xyf50A3zPj6XeUavzffsJE+7zLSHJbklt30fucP5mz7/liS/tHb+nrDDul+XE74jZfxB46mnPU5X61/LN2WvTXK/E/6efzzDMAyvyvjL3Sdk/IL5ncMw/MwJZf3UcPu2y0/NeJEdNwN7YJIXZPyl/o2/KCZ55rTOZEx8feck33+Kfdml/HdO8o5Jfmxt2f+54zrOZXygPCJ5468QD03ylGGqsWVsGvHEUspfJHl9xsr9wCTVb12mcn4xyS0r+/eGJP87t+1jhmF4/DAMZRiGm3coc5j5+ebChuEHMiaPP2UYhuNk93+a5NXDMDx9Zb5/yPgw/rDpV8K5nr5thqk+/kim8zV5WJI3T/KU4w/K2Lzwx0opf53xeL4uY11sdc5+IsnrV87ZK5L8Vm5/zm4YhmFj88CVX+Nu92v8DCXjl5yHDcPw9GEYfiTJv0jy3hmbDx37+oz1+nOTfGTGTmBuKKV89cz1wdXOM3R0zT5Dy9gM/v9N8iXDMPxy6+m7GobhjzLWmT8ahuH4/HxwxmZ0P7Q2+1MyPgdnpSlMnj3s9tblKUk+tIx56imlvH3G58kbe9gspbx1KeW/lLEp5Oumv89Nu3O9rT5nGIabpnP97MsVVEarz9717+E3TmV+fMbOPJ5eSvm/GuzDSdtyJmPnLV813DE3c3W++2R8K37aPPCTvgOWy3x+TWvd0cfzdpjv/2SMsN81l2/C9NLVfwzD8IpSyusy5icl468uH5Dxolt3/CA6bpt8mh6hdin/eFteujZ9Y881x4ZheGkp5RkZm198c8bmeG+aqdnFdKH+ZJK3ytiJxh9nbBrxNUn+0a47ssFRxiZoJ+3jn8ws6xXTf9967fO3Xpt+Gq/L+OvcapknHeO/yvjAuGvumMO0zfo5vJwnJ/mCUsoHDcPwaxmbO/7KMAwvTt74peDnM/7696UZexu6lOR7k5yZuU0nOcrYG9MXnzBtbuck35/bJ+p+VsZfyHb1iiQvWr2RD8Pw66WUVyb5x0kyNfV8TJJPGIbhJ6fZfqmU8g8Ze/j6jmEYdj32cLXzDB1di8/QY8cBxDMWmj7H3+eOz95kfNauOv733U+xjl3v/0/L2DLjkRl7mHzEtG0/vjLPTUk+LOP5/Z2MvQP+v7n9D6mntUt93tVnZmzWf+yJGd+wJUmGYXhJpmuulPLTGZsJf03GHxxa++KMP2w8pYy9JCZjKsqdpn//3TAMf5+xY48fSXLzynx3SnKX6d+vGi6fw/aK3PG7YTI2M9bJx5qWQdmuHp8xkfePktxYSvmo41+1Vtwu0Xdqu/0mSf5y+uiWJL+dsUnG5bx8+u9pxk/ZpfzjbVlPSn679Rk3eHKSJ5ZS3ivjg+W3hmH4nWnau2e8CXziMAw/cbxA2T62w6Xpv3dZ+3z9hnlLkp9NclJb8Ft32PZVf5LxBnnvqcxj95n+2/LCuyUnH+O3z3jDPB6z41K2H4Nju/5a8ysZA61HllJ+P2PPVqsJwh+a8dffhw4rbbNLKW+V5M83lLvpnK0+FG9J8tMZf0lb9+pddmDF43P7L3Qvmrn87+Xygebxzfn4/P/m2vTfzHjvuWd2fygDo8fHM/TY1fIMPfZ3Sf5gZRtaT69xy/Tft8vY7PTY269Nb/7sHYbhUinlx3NbUPbIjGkHr0re+MbnIUm+bLitp9/jwHybW3fY3l3q865+Kitv1zL2FHqiYRj+oZTymxmD/yW8d8br5GUnTHtFxqD2u6b5PjZ37ML+X05/986YpnGS38ta7lgZc+DfLad7C39V22tQVsaxrP5tkkdnbJLxK0m+JGMUvuphpZQvXWl+8c8zXrzHPeM9I2Mey0umXxVO8gcZvwh/VsambSc5/sK7/uVyl/L/POND5ZNy++YXcwYl/LGMFf5fJ3lAbt/d+/GD441fysvYlfiHJ/nDXN5Lc1uAdLzc9Rmb/a3+ovOMjBfY7w3D8LcztvkOhrEb+2dlbL72bSuTPjXjMfqNmvLX/HKSx5ZSHjQMw88mb7zxfkrGtuPH+/jnuWMS6T9LhWEYhlLKuYxJ4y/MeP2s1q2TztmHZWzH/ju5vOOA7d4ZfwXPlCj+vhnHITn2jCT/JMlvDJVdE09NaW6uKOJpGYc7uMcwDC9LkjL2HvlWGZtTJmMAmyQfmLGr3WPHvUzWrB+uOZ6hd3BVPEOPDcPw8xm/AC8yvdKvZfzh819kfHtz7BEZ30Aef/bnSd65lHLXYRiOfyStevZOnpLkM0opH5sxSPnklWnXZ8xRXD3Xb5GxCeA2J31X+Ji1f+9Sn3cyjD1cvnzrjElKKW+Sscnxn9asc4Ovzx1byHx5xqbLn5XbrpNH5o7X+LmM959vz+2f7+uennH4ireZ9j0Zr/nrs0PqyDWnRWJaxl/u/i7jhbL+927TPG+e8Ze9p68s95UZf1X5xyufDRl/hfmpjG8iHp3xLcBqMuX1Gb+s/mHGNsP3z9hzz1cn+U8r8z0i46/2P5qxDf6DpnkeOk0/ytjO/PszvuW478zyv2ja3m/K2Ob4WzJe4Dsneua2dvHrCcvHCcO/lfG19SMy/hLxoqwkYeaEROupzFdkfE3+kIyJzX+W2ycpv03GL8XPyxhkfGTGh+E3ZWyPfjzf46ZjdM8t+/GhGW/Y3zNt01dkfIB9ztp8z17djh2Oz81ZSSbNeON9TsZfdv7lVEd+YtrG+6/M9/nTcfmqjA+E/5wxSHhjwnUuk9C9ZXveN7fV0Z9bm/Z2Gevqs6b68FnTOv88K0m2WevoY/rsV6dz9PCMN6znTOt49so87z6d12dkfDB+5FQvbkzyqSvzfV/G8V5Oez0/fPr74YzNIo//fc+Ved5y2rfnZHzwPTLjg+P5Se60dq7+ejofH5Xk32R8m3muxb3Hn7+r4S+eodf8M3Sa9zM2zVszPfM7HLkpd+yw45un8r96Ol/fMB33r1yZ5z7TZz+c8dn7xRl/lBxyx44+tnYUsTL/nTMGzH+RMUXhzNr0X8v4THr4VNd+NeMz6TUr89wrd+zo4+szXkNfmPFt0JOmcz2c4no5Pv4feYp7wOdmTHV45FSXHpkxZ/Hvk/zTtXl3eUbfY+Xz35m2//jfl+1c5aTzfpn57nD+pu1+fZLPWPnsbhmv6V+eju+jMn5/+8G5x+ha+GtTyOaeo26a5vmujL8OvMPKctdNF87zk7zJ9NmQMVL/lmn+V2fsYvwt19b5lhl/HXzxVGlfkrFp10PW5nvYtI7XZrzJPjPjIHbH0z8vY/O7161dhFvLz5io+LiMXzr/NmNw8JDMe6A8bJr/f50w7X4ZbzSvzXgz+Iz1CyYnP1DukfEXxFdOF8MXZa3nqGm+t894E3hJxlf4FzK2G/6wE87tvXbYlwdnbJp23N3pF5wwz3OT/OqMunVz7njh3z1j4HEx483015J87No8d874cPyrjD1kfVfGm0FVUDYt98KsPWBWpj1omv7ajF8GPi5rPR/l5KDs3TIGc6/J+MXrEVnrfXGa7z0yJlof7/uLMrZJX/1SdlMqel/M5a/lG9bme7dpG18z1bVzWbm+p3neNmMPWC9aqcf/KVNPXf78+fMMjWfo8bw3bJq3ZnrG5+GlJHfb8bje7jhNn90p4w+uN0/n9I+zEoCuzPeojM+xv8uYZ/2BqQzKpmW+IyvXxNq0d8/4DP3bjEHVY7LWO2FODsrePOOPCi/PGPT9x4w/Hg5r5e9Sn4+P/051d638D8/YHPavp7r04oyd3rzfCfNufUav1OuT/i5bF08675eZ7w7nb2WdN6x9/p4Zf9j424wB2X9J8qZzj9G18FemA3Ywyji+0mOHYZg9vgaHbWr3/cokjxqG4XLNYQA4Jc9QTlJK+d9JXjAMw6N7bwtwsh4dfXDtOm4bvWuXxwBAhTIOfvx+uWNHDcABEZSxN8Mw/FKM4A4AezOM3Zq/Ze/tADY7uOaLAAAA15JdxnAAAABgIYIyAACAjgRlAAAAHQnKAAAAOhKUAQAAdCQoAwAA6EhQBgAA0JGgDAAAoCNBGQAAQEf/P6Gg7Veo2/nQAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 1100x1100 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def energy(state):\n",
    "    E = 0\n",
    "    N, M = state.shape\n",
    "    for i in range(N):\n",
    "        for j in range(M):\n",
    "            # handle the north and south neighbours\n",
    "            for di in [1, -1]:\n",
    "                if 0 <= (i + di) < N:\n",
    "                    E -= state[i, j] * state[i + di, j]\n",
    "\n",
    "            # handle the east and west neighbours\n",
    "            for dj in [1, -1]:\n",
    "                if 0 <= (j + dj) < M:\n",
    "                    E -= state[i, j] * state[i, j + dj]\n",
    "\n",
    "    return E\n",
    "\n",
    "\n",
    "expected_values = [\n",
    "    E_prediction_all_the_same(100),\n",
    "    E_prediction_all_the_same(100),\n",
    "    0,\n",
    "    \"???\",\n",
    "]\n",
    "\n",
    "f, axes = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))\n",
    "axes = axes.flatten()\n",
    "for ax, state, exp_value in zip(axes, states, expected_values):\n",
    "    show_state(state, ax=ax)\n",
    "    ax.text(\n",
    "        0,\n",
    "        -0.1,\n",
    "        f\"Expected Value: {exp_value}, Your value: {energy(state)}\",\n",
    "        transform=ax.transAxes,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "84ee0a62-6ab1-4565-a236-2ccfcc378acc",
   "metadata": {},
   "source": [
    "It's a bit tricky to know what to do with the random value, let's try running it 100 times and see what we get:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "686782b5-aa37-4084-9786-538fdc89fef1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean = -0.24, standard error = 29.819071481184658\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAe0AAAHOCAYAAAC1jFFjAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAABDrAAAQ6wFQlOh8AAAjZUlEQVR4nO3de5RlVWHn8e/Pbm1AaGlUbCNKI448THxEdKKJD+yMYhMtE0BFl0JGjWAeGkTER9RBo2hc+EpmcMQMZgZLJS4psCGiJcJAREBgDPLIQGhEh5fSnZZGmteeP865crncet3uW3V39fez1lm3zz57n7t336r7q/NOKQVJkjT6HrbQHZAkSbNjaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJZYudAeGaXJychnwW8BtwH0L3B1JknotAR4L/Mvq1as3z1R5UYc2TWBfvNCdkCRpBs8BLpmp0mIP7dsA9ttvP5YtW7bQfZEk6UE2b97MJZdcAm1ezWSxh/Z9AMuWLWO77bZb6L5IkjSVWR3C9UQ0SZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiqxdKE7IKkuq45du9BdqM664w9c6C5okXBLW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSarEnEM7ySFJTktyY5JNSX6U5MgkD+uptybJZUnuSnJtkrfNcv0PT/KxJDcluTPJOUmePtd+SpK02Ayypf1OYDPwLuAPgNOAzwIf71RI8jxgArgUeDlwMvC5JG+exfo/Bfwp8AFgDLgXmEyycoC+SpK0aAzyPO1XlFJu65o/J8mOwJ8leX8pZTNN4F5aSnlTV50nAccl+ftSyv39VpzkCcARwF+UUr7Qll0IXA+8Azh2gP5KkrQozHlLuyewOy4DtgN2SbIMeAnwlZ46pwCPB541zepfCizpbltK+SVwBuBT5CVJ27StdSLaC4DbgVuBPYFHAFf11Lmyfd1nmvXsA9xSSrm9T9u9eo+b90qyPMlunWl8fNxd6pKkRWOLQzvJfsAfA58qpdwHrGgXbeipur593WWa1a3o067T9uHAjjN05yjgxs40MTFx8Qz1JUmqxhaFdnty2NeBi+g6Ea1Vpmg2Vfl0yzPLticAT+xMY2Njz5mhviRJ1RjkRDQAkjwKOAu4E3hlKeWedlFni3pFT5MVPcv7Wd+nHcDOwD3Apun6VErZCGzszE9OTk5XXZKkqgy0pZ1kO+B04HHAAaWUX3Qtvg64m4ceu963fe091t3tKmDXJL270PcFrpnqrHNJkrYFg9xcZSnwNeAZNIF9Q/fy9pKv7wKv7ml6KHATzZnmUzkbuL+7bXs52SuAtXPtqyRJi8kgu8f/jiZEjwF2SPI7XcuubHdRHwecl+QLNJd6/S7wFuCt3VvLSa4FbiilrAYopfwsyYnAx5PcC9wAHN1W//QAfZUkadEYJLRf1r5+os+y/YHvlVK+n2QM+CjwRuCnNDdMOanP+y/pKTsKuAP4CPAo4AfA6lLKzQP0VZKkRWPOoV1KWTXLemcCZ851XaWUu2nufObdzyRJ6uJTviRJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUGCu0kT0lyYpLLk9yb5Io+dco00+NnWH+/NjcP0ldJkhaLpQO2expwIPADmuDvF/7P61P2D8CmUspNs3iPzwFf7pq/e66dlCRpMRk0tM8opUwAJDkZ2K+3Qinlwu75JKuA/wAcM8v3+EnvOiRJ2pYNtHu8lHL/AM1eBxRgfJD3lCRpWzefJ6IdCpxXSvnpLOsfm+SeJBuSfDXJk2ZqkGR5kt060/j4+Mot67IkSaNjXkI7ydOB3+TBx6in8w/AkcBq4L3AC4Hzk6yYod1RwI2daWJi4uLBeixJ0ugZ9Jj2XL0euAf4x9lULqUc1jV7XpLzgUuBtwCfmKbpCcBJnZmxsbGVgMEtSVoUhh7aSQK8FjirlHL7IOsopfwoyTXAs2eotxHY2JmfnJwc5O0kSRpJ87F7/PeAJzH7XeNTyVboiyRJ1ZqP0H4dcAdwxqArSPJM4Km4q1uStA0baPd4kh2ANe3s7sDyJAe38+eWUm5r6y0FDgZOK6XcOcW6rgVuKKWsbuePBp4MnAvcSnMC2/toTi47qd86JEnaFgx6THtX4NSess78/sD32n+/DHgM0+8aXwos6Zq/BjiI5jj4TsBtwFrg/aWUDQP2V5Kk6g0U2qWUdcziGHMpZe1M9Uopq3rmz2ALdqVLkrRY+ZQvSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUiYFCO8lTkpyY5PIk9ya5ok+dk5OUPtMBs1j/w5N8LMlNSe5Mck6Spw/SV0mSFoulA7Z7GnAg8AOa4J8q/P8NeH1P2VWzWP+ngDcC7wTWAccAk0l+q5Ry8yAdliSpdoOG9hmllAlotqiB/aao96tSyoVzWXGSJwBHAH9RSvlCW3YhcD3wDuDYAfssSVLVBto9Xkq5f2t3pMtLgSXAV7re75fAGTRb95IkbZOGfSLankk2JLk7yQ+TvGoWbfYBbiml3N5TfiWwV5Ip+5xkeZLdOtP4+PjKLei7JEkjZZihfRlwNPAq4NXAz4FvJDl4hnYrgA19ytcDDwd2nKbtUcCNnWliYuLiuXVZkqTRNegx7RmVUj7TPZ/kdOCfgeOAf5ypeZ+yTLOs4wTgpM7M2NjYSsDgliQtCkML7V6llPuTfB34RJLtSym/mqLqepqt7V47A/cAm6Z5j43Axs785OTk4B2WJGnEzPfNVTJzFa4Cdk2yS0/5vsA1Qz4JTpKkkTVvod2eQHYw8ONptrIBzgbupzkO3mm7I/AKYO1QOylJ0ggbaPd4kh2ANe3s7sDyrhPMzgV2AE4GxoHraHZ3H0lzPfdBPeu6FrihlLIaoJTysyQnAh9Pci9wA80JbQCfHqS/kiQtBoMe094VOLWnrDO/P/AjmmPLHwAeC9wNXAK8vJTyrT59WNJTdhRwB/AR4FE0d15b7d3QJEnbsoFCu5SyjpmPT4/Ncl2r+pTdTXPnM+9+JklSy6d8SZJUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiphaEuSVAlDW5KkShjakiRVYtAHhkgaolXH+hRaSQ/llrYkSZUwtCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlVi60B2QpMVu1bFrF7oLU1p3/IEL3QXNgVvakiRVwtCWJKkShrYkSZUwtCVJqsRAoZ3kKUlOTHJ5knuTXNGzfEmSY5Kcm+S2JOuTnJdk9SzXX/pMNw/SV0mSFotBzx5/GnAg8AOa4O8N/+2B9wJfAv4GuAc4HPh2kleWUr45i/f4HPDlrvm7B+yrJEmLwqChfUYpZQIgycnAfj3LfwXsUUpZ3ylIcjbwVOCdwGxC+yellAsH7J8kSYvOQLvHSyn3z7D8vu7AbssKcDnwG4O8pyRJ27p5OxEtycOA5wNXzbLJsUnuSbIhyVeTPGkW77E8yW6daXx8fOUWdVqSpBEyn2eP/zmwF3DCLOr+A3AksJrm2PgLgfOTrJih3VHAjZ1pYmLi4sG7K0nSaJmX25gmeRHwCeCTpZTzZqpfSjmsa/a8JOcDlwJvadczlROAkzozY2NjKwGDW5K0KAw9tJM8HZgATgPePcg6Sik/SnIN8OwZ6m0ENnbmJycnB3k7SZJG0lB3jyfZE/gWzVbyG9qT0QZe3dbplSRJdRpaaCdZCZwN3Ay8qpQy8HXWSZ5Jc7mYu7olSdusgXaPJ9kBWNPO7g4sT3JwO38ucAfwT8CuNCeH7Zs8sKHcff11kmuBG0opq9v5o4Ent+u5FfhN4H00J5f9+ni1JEnbmkGPae8KnNpT1pnfH1gHPKOdP61P++5d3UuBJV3z1wAHAa8FdgJuA9YC7y+lbBiwv5IkVW+g0C6lrGPmY8yzOgZdSlnVM38GcMYg/ZIkaTHzKV+SJFViXq7TliSNplXHrl3oLkxp3fEHLnQXRo5b2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKDBTaSZ6S5MQklye5N8kVU9Rbk+SyJHcluTbJ22a5/ocn+ViSm5LcmeScJE8fpK+SJC0Wg25pPw04ELgWuLJfhSTPAyaAS4GXAycDn0vy5lms/1PAnwIfAMaAe4HJJCsH7K8kSdUbNLTPKKU8sZRyME0o9/MB4NJSyptKKeeUUj4CfBE4LsmU75vkCcARwLGllC+UUr4N/BEQ4B0D9leSpOoNFNqllPunW55kGfAS4Cs9i04BHg88a5rmLwWWdLctpfwSOINm616SpG3SsE5E2xN4BHBVT3lnV/o+07TdB7illHJ7n7Z7zbCVvjzJbp1pfHzc3emSpEVjWKG9on3d0FO+vn3dZYa2ve06bR8O7DhN26OAGzvTxMTExTN1VJKkWgz7kq8yx/LplmcWbU8AntiZxsbGnjPD+0iSVI2lQ1pvZ4t6RU/5ip7lU7XtbQewM3APsGmqhqWUjcDGzvzk5ORM/ZQkqRrD2tK+Dribhx673rd97T3W3e0qYNckvbvQ9wWumekkOEmSFquhhHYpZTPwXeDVPYsOBW4CLpum+dnA/d1tk+wIvAJYu3V7KklSPQbaPZ5kB2BNO7s7sDzJwe38uaWU24DjgPOSfIHmUq/fBd4CvLV7aznJtcANpZTVAKWUnyU5Efh4knuBG4Cj2+qfHqS/kiQtBoMe094VOLWnrDO/P/C9Usr3k4wBHwXeCPwU+ItSykl9+rCkp+wo4A7gI8CjgB8Aq0spNw/YX0mSqjdQaJdS1vHA2dzT1TsTOHOGOqv6lN0NHNtOkiQJn/IlSVI1DG1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJQxtSZIqMbTQTvK9JGWK6bUDtNt7WH2VJKkGS4e47rcBy3vK3gEcBHxnhrYXAEf3lK3bKr2SJKlSQwvtUsqVvWVJngucXUr5+QzNN5RSLhxOzyRJqtO8HdNO8nxgD+CU+XpPSZIWk/k8Ee11wJ3AxCzqvijJpiR3JTk3yQtn8wZJlifZrTONj4+v3KIeS5I0QuYltJMsBQ4BJkopm2aofi7wduAA4DBgB+A7SZ43i7c6CrixM01MTFw8eK8lSRotwzwRrdt/AnYFvjxTxVLKB7vnk3wT+DHwV8CaGZqfAJzUmRkbG1sJGNySpEVhvkL7dcAvgG/NtWEpZVOStcDBs6i7EdjYmZ+cnJzr20mSNLKGvns8yfbAGHBqKeWeQVezFbskSVKV5uOY9iuBnZjFrvF+kjwSOBB3c0uStnHzEdqvA34CnN+7IMkXk9zbNf+CJBNJDk+yf5LXA/8bWAkcNw99lSRpZA31mHaSFTRngX+6lFL6VFnSTh03AcuAjwGPBjYB/wwcUUq5aJh9lSRp1A01tEsp62lCeKrlhwOHd81fSxPykiSph0/5kiSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJZYudAe0+K06du1Cd0GSFgW3tCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKDC20kxyepPSZjp9F28OSXJ3kriRXJDlkWP2UJKkW8/FozgOAf++a/9l0lZMcDJwMHA+cDbwK+GqSfy+lnD2kPkqSNPLmI7R/WEr5+Rzqfxg4tZTynnb+nCR7A8fRhLgkSdukkTqmnWQPYG9gvGfRl4HnJnnM/PdKkqTRMB+h/eMk9yX5tyTvSbJkmrr7tK9X9ZRfCYQm0KeUZHmS3TrT+Pj4yi3otyRJI2WYoX0T8EHgjcDLgTOBjwCfmabNivZ1Q0/5+vZ1lxne8yjgxs40MTFx8Rz6K0nSSBvaMe1SyreAb3UVnZ3kV8BfJvnrUspN0zXvmc8U5b1OAE7qzIyNja0EDG5J0qIwHyeidfsacDTwTJot8V6dLeoVwC1d5Tv3LO+rlLIR2NiZn5ycHLCbkiSNnvk+ES0zLO8cy96np3xfmq3sq7d6jyRJqsR8h/ZrgPuAy/otLKVcTxPMr+lZdChw0RwvHZMkaVEZ2u7xJN8CJoEr2qJXAn8CfKaUcnNb54vAYaWU7n58gOZmKtcB3wbGgJfS3KRFkqRt1jCPaV8NvBnYjWaL/l+BdwCf66qzpJ1+rZRyapIdgPfSHP++FniNd0OTJG3rhnn2+NuBt89Q53Dg8D7lXwK+NJSOSZJUqZG6I5okSZqaoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFVi6UJ3QJKkflYdu3ahuzCldccfuCDv65a2JEmVMLQlSaqEoS1JUiWGFtpJDklyWpIbk2xK8qMkRyaZ9j2TfC9J6TPtPay+SpJUg2GeiPZO4AbgXcAtwP7AZ4Ent2XTuQA4uqds3VbunyRJVRlmaL+ilHJb1/w5SXYE/izJ+0spm6dpu6GUcuEQ+yZJUnWGtnu8J7A7LgO2A3YZ1vtKkrRYzfeJaC8AbgdunaHei9rj4HclOTfJC2ez8iTLk+zWmcbHx1ducY8lSRoR8xbaSfYD/hj4VCnlvmmqngu8HTgAOAzYAfhOkufN4m2OAm7sTBMTExdvWa8lSRod83JHtCQrga8DFwEfn65uKeWDPW2/CfwY+CtgzQxvdQJwUmdmbGxsJWBwS5IWhaGHdpJHAWcBdwKvLKXcM5f2pZRNSdYCB8+i7kZgY2d+cnJyjr2VJGl0DTW0k2wHnA48DnheKeUXg65q6/VKkqQ6DS20kywFvgY8A3hhKeWGAdfzSOBA3M0tSdrGDXNL+++AVwDHADsk+Z2uZVeWUjYm+SJwWCllKUCSF9DcVOUbNDdm+Q2am7SsBA4ZYl8lSRp5wwztl7Wvn+izbH/ge8CSduq4CVgGfAx4NLAJ+GfgiFLKRUPrqSRJFRhaaJdSVs2izuHA4V3z19Jc6iVJknr4lC9JkiphaEuSVIl5ubnKYrHq2LUL3QVJ0jbMLW1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqoShLUlSJQxtSZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXC0JYkqRKGtiRJlTC0JUmqhKEtSVIlDG1JkiphaEuSVAlDW5KkShjakiRVwtCWJKkShrYkSZUwtCVJqsRQQzvJU5P8U5JNSW5N8pkk28+y7WFJrk5yV5IrkhwyzL5KkjTqhhbaSXYGvgvsBBwEHA28HvjCLNoeDJwMfAN4OTAJfDXJS4fUXUmSRt7SIa77rcAK4JmllJ8DJLkXOCXJX5dSrpqm7YeBU0sp72nnz0myN3AccPYQ+yxJ0sga5u7xNcB3OoHd+jqwuV3WV5I9gL2B8Z5FXwaem+QxW7ujkiTVYJhb2vsAf99dUErZnOS6dtl07QB6t8SvBEIT6Of3a5hkObC8M3/kkUc+4aCDDmLz5s1z7Hp/K5ZtldVIkip31113bZX1dOXTktnUH2ZorwA29ClfD+wyQzv6tF3fvk7X9ijgg52ZyclJDjroIC655JJpOzpbn3zhdltlPZKkul1wwQVbe5WPBW6YqdIwQxug9CnLFOUztc006+w4ATipM7N06dJHXH/99U/cY4891gH3zeI9F8T4+PjKiYmJi8fGxp5z6KGH3rzQ/dkSjmU0OZbR5FhG1zyOZwlNYP/LbCqnlNnk59wluRX4+1LKsT3lPwa+X0p58xTt1gBrgX1KKVd3lT8HuAh4QSml7+7xWiXZDbgReGIp5acL3Z8t4VhGk2MZTY5ldI3qeIZ5ItpV9By7TrIM2JOHHq/ubUdvW2Bfmq3sq5EkaRs0zNA+E1id5NFdZX8ILGuX9VVKuZ4mmF/Ts+hQ4KKes9ElSdpmDPOY9ueBPwcmknwY2JXmmPMp3ddoJ/kicFgppbsvH6C5mcp1wLeBMeClwAFD7O9C2gj8l/a1do5lNDmW0eRYRtdIjmdox7ShuY0p8Dng94A7aa69fncp5VdddU6mCe30tD0MeC+wCrgW+FAp5dShdVaSpBE31NCWJElbj0/5kiSpEoa2JEmVMLQlSaqEoS1JUiUM7XmSZGmSdye5OsmdSdYl+Uz73PHeuke3y+9KcnGSF/eps1OSzyf5RZI7kpyeZPd5GErn/bdP8tEkNyTZ3Pb3g33qjfxYevry7CT3JbljiuUjO54kS5Ick+TcJLclWZ/kvCSraxtLP0memuSfkmxKcmv7+7P9QvRlKkkOSXJakhvbfv4oyZFJHtZTb02Sy9r/+2uTvG2K9c34Gc2XJDsm+WmSkmS/nmVVjCfJm5L8n/b9b01yes/y0R9HKcVpHiaaZ4HfQ3MZ2/7A22geinJ6T72jgbvb15fQXCb3K+C3eup9E/h/NDedORD4IfB/ge3nYSxLgEngGuAw4EXAG4D31TaWnn4E+D5wM3BHn+UjPR5gx/Zn6jPAHwAva/t4P/AHNY2lz9h2Bn4KXEBzv4Y3Aj8H/td89mMW/bwQ+Crw2vb3vPN7/zdddZ7Xln2xrfN+mmcjvHmQz2gex/bx9nejAPvVNh7gQ8C/A8e031l/CHy+unEs9A/5tjLRXGv+pZ6yY9ofike288vaL91PdNVZQvNY0q90lf3H9hdnTVfZk9ofuCPmYSx/QvPUtcdNU6eKsfT0+T+3gfRRekK7hvG0/VnRUxaaoD2nprH0Gdu7gU3AY7rKXtf2b5/57MsM/Xxsn7IT2i/1Ze38WcAPeur8d5o/jh42l89oHse1N3AH8FYeGtojPx6a22LfC7x0mjojP45SirvH59HDaf7K67aB5ku1c2OZ5wOPovnLDYBSyn00f7mvSdKpt6Zte1ZXvZ/QPGf8wK3f9Yd4E/C1Usot09SpZSwAtIcpjgf+kuav6F4jP55Syn2llPU9ZQW4HPiNruKRH0sfa4DvlAffxvjrwOZ22UgopdzWp/gyYDtglzTPX3gJ8JWeOqcAjwee1c7P9jOaL58FTqTZu/ZrFY3ncODfSiln91tY0TgM7Xn0eeANSX6/PTa0H83ulZNLKZ3jp52HpPQ+FOVKYCfgCV31rmm/kHvr9T5oZatK8gjgt4Ebk/zPNMfnNyb5ch58n/mRH0uPjwA/LKV8c4rltY0HgPZY6vN58EN6ahzLPvQ8aKiUshm4bgH6MlcvAG4HbqV5YNIjeOhDk65sX/fpeZ3pMxq6JAcDz6DZ1d+rlvH8DvAvSf6qPZZ9d3vuxzPb5bWMw9CeL6WUjwKfBs4GfglcTPP81Ld2VVsBbC5dt3ltdbaedumqt6HP26zvqjMsj6a5Z/27aY4zvormHvO/z4P/Sq1hLAC0v7hvotnKnko14+nx58BeNLtoO2ocyyj1ZdbaP87/GPhUu0W2ol20oadqv//72XxGQ5VkB5qfnfeUUvrdg7uW8aykeX7F64EjgD8CdgC+3e5lq2UcQ31gyKKW5FE0u01mcn0pZXOSP6MJhXfSHGPcC/gwcBLNyVwd/e4rmz7Lpqo35/vSzmUsPPCH3gbgoFLK3e06fgl8PclzSykXzdDH3mVbbSxtX+YynruBvwX+a+l6fvsURvqzabc8u9u+CPgE8MlSynk99Rfks9lCo9SXGSVZSbML/yKak7i6TdXn2fzfT9d+a3s/cAtw8gz1Rn08D6M5UfOgUsqPAZL8kOY74E9oTnCcrh+jMg5Dewv8IfA/ZlHvWUluBD4JHFNK+Wxbfl6SW4HTknymlHIpzV9r2yXZrpRyV9c6dm5f13e9PqnPe+3cVWcuZj0W4F/bf1/QCezWd9vXp9F8SS3UWGBu49mb5lntr88Dl99tB78+zn1X2/8aPpvLOzNJng5MAKfR7BXptpCfzaDW88DWULedeeguzQXX/rF1Fs2Dkl5ZSrmnXdT5f+sdy4qe5bP9jIamvbTvnTQ/g8vbw7U7tot3TLIj9YznduCWTmADlFJuSnI1zXdW57DYqI/D3eODKqWcXErJLKbLaY6XLKPrS7XVmd+zfe18+fQeo9uXZpf6z7rq7dXnpId9GeALbC5jKaXcCdwwzeruX8ixzHU8NKG9AlhH8wu3nibkHtn++0MLOZ45jgWAJHsC3wIuBd7Q55j0gn02W+Aqevrbnjy05wL0ZVpJtgNOBx4HHFBK+UXX4uto9u70+7+HB8Yy289omPagOc67lgd+N85ol50DfId6xjPVz0hovrNqGYeXfM3HRPPLW4C/7Ck/qC1/Tnnw5QTHd9VZAvyY/pfiHNBV9kTm75Kvv6W5ZnZZV9nBbZ+eWdlYVgEv7plOprlE58XAkysbz0qaL6DLgOVT1KliLD19fjfNJUeP7ip7LaN3yddSmsDeADx9ijpnAd/vKTuR/pcWTfsZDXksO/f53XhH+3/+VuC3axkPD3w//WZX2RNorj54Ry3jKKUY2vM1Af9I85fYu2gu3D+S5mzSCzo/EG29zoX772zrncLUN734WfvFtQa4hPm7ucru7Q/umcDLaU7iug34Rk+9kR/LFOP7ENPfXGUkxwNsT7P35pfAGM0Zs7+eahpLn7HtTPOH4vk0N415Q/szN2o3V/k8TTi8q/f/n/aPKB64iccXaILwfUx/E49pP6N5Ht+LmfrmKiM7Hppg/SHN4b1X05xAe2n7M/XIWsZRiqE9bxPN5QAfb7/wfkVzAsR/o+dmDDS7a95Fswv6LpqzzPfvs77lNBf+306zBXI6sPs8jufZwLntWH7eflntVONY+vTnQ/QP7ZEeD81egzLVVNNYphjfU2l2+2+iCezPskB/2E3Tx3XTfAYv7qq3huYPrM5la3866M/bPI/vxfSEdi3jAXalCdgN7c/QmcBetY0jbSckSdKI80Q0SZIqYWhLklQJQ1uSpEoY2pIkVcLQliSpEoa2JEmVMLQlSaqEoS1JUiUMbUmSKmFoS5JUCUNbkqRKGNqSJFXi/wPMNWVyXhrQBAAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 550x550 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "L = 100  # How large the system should be\n",
    "N = 100  # How many random samples to use\n",
    "energies = [energy(np.random.choice([-1, 1], size=(L, L))) for _ in range(N)]\n",
    "plt.hist(energies)\n",
    "print(f\"mean = {np.mean(energies)}, standard error = {np.std(energies) / np.sqrt(N)}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8da38c5b-f166-4530-b046-3600b4d84c5d",
   "metadata": {},
   "source": [
    "If you run this a few times you'll see the mean is usually within a few standard errors of 0, which gives us some confidence. In the testing section we will discuss how we might go about doing automated tests of random variables like this. "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3bd3831b-f8de-4a8e-aaf6-7aa158cf9d82",
   "metadata": {},
   "source": [
    "### Making it a little faster\n",
    "\n",
    "This project is not intended to focus on optimising for performance but it is worth putting a little effort into making this function faster so that we can run experiments more quickly later.\n",
    "\n",
    "The main thing that slows us down here is that we've written a 'tight loop' in pure python, the energy function is just a loop over the fundamental operation:\n",
    "```python\n",
    "E -= state[i,j] * state[i+di, j]\n",
    "```\n",
    "which in theoy only requires a few memory load operations, a multiply, an add and a store back to memory (give or take). However because Python is such a dynamic language, it will have to do extra things like check the type and methods of `state` and `E`, invoke their array access methods `object.__get__`, etc etc. We call this extra work overhead.\n",
    "\n",
    "In most cases the ratio of overhead to actual computation is not too bad, but here because the fundamental computation is so simple it's likely the overhead accounts for much more of the overal time.\n",
    "\n",
    "In scientific python like this there are usually two main options for reducing the overhead:\n",
    "\n",
    "#### Using Arrays\n",
    "One way is we work with arrays of numbers and operations defined over those arrays such as `sum`, `product` etc. `Numpy` is the canonical example of this in Python but many machine learning libraries are essentually doing a similar thing. We rely on the library to implement the operations efficiently and try to chain those operations together to achieve what we want. This imposes some limitations on the way we can write our code.\n",
    "\n",
    "#### Using Compilation\n",
    "The alternative is that we convert our Python code into a more efficient form that incurs less overhead. This requires a compilation or transpilation step and imposes a different set of constraints on the code.\n",
    "\n",
    "It's a little tricky to decide which of the two approaches will work best for a given problem. My advice would be to have some familiarity with both but ultimatly to use what makes your development experience the best, since you'll likely spend more time writing the code than you will waiting for it to run!"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6b2ed23c-32f2-40af-bf45-e60004c275a2",
   "metadata": {},
   "source": [
    "## Exercise 2: Write a faster version of `energy(state)`\n",
    "\n",
    "You can use `numpy`, `numba`, `cython`, or anything else, by what factor can you beat the naive approach? Numba is probably the easiest here."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "40386cda-cc14-4fdd-8f9f-71840d1ebb40",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Naive baseline implementation\n",
      "80.9 ms ± 2.94 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "\n",
      "Your version\n",
      "80.2 ms ± 4.99 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Your speedup: 1x !\n"
     ]
    }
   ],
   "source": [
    "def test_energy_function(energy_function):\n",
    "    assert np.all(energy_function(state) == energy(state) for state in states)\n",
    "\n",
    "\n",
    "def time_energy_function(energy_function):\n",
    "    return [energy_function(state) for state in states]\n",
    "\n",
    "\n",
    "def your_faster_energy_function(state):\n",
    "    return energy(\n",
    "        state\n",
    "    )  # <-- replace this with your implementation and compare how fast it runs!\n",
    "\n",
    "\n",
    "print(\"Naive baseline implementation\")\n",
    "test_energy_function(\n",
    "    energy\n",
    ")  # this should always pass because it's just comparing to itself!\n",
    "naive = %timeit -o -n 10 time_energy_function(energy)\n",
    "\n",
    "print(\"\\nYour version\")\n",
    "test_energy_function(your_faster_energy_function)\n",
    "yours = %timeit -o -n 10 time_energy_function(your_faster_energy_function)\n",
    "print(f\"Your speedup: {naive.best/yours.best :.0f}x !\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a27bb821-c33f-4966-8676-4295cfd695ab",
   "metadata": {},
   "source": [
    "## Solution 2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "c96720d4-2c0b-4fa5-b9ce-4fc7be533a62",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Naive baseline implementation\n",
      "81.2 ms ± 3.2 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "\n",
      "Numba version\n",
      "212 µs ± 20.5 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Numba Speedup: 420x !\n",
      "\n",
      "Numpy version\n",
      "153 µs ± 37.5 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Numpy Speedup: 576x !\n"
     ]
    }
   ],
   "source": [
    "def test_energy_function(energy_function):\n",
    "    return [energy_function(state) for state in states]\n",
    "\n",
    "\n",
    "from numba import jit\n",
    "\n",
    "\n",
    "@jit(nopython=True)\n",
    "def numba_energy(state):\n",
    "    E = 0\n",
    "    N, M = state.shape\n",
    "    for i in range(N):\n",
    "        for j in range(M):\n",
    "            # handle the north and south neighbours\n",
    "            for di in [1, -1]:\n",
    "                if 0 <= (i + di) < N:\n",
    "                    E -= state[i, j] * state[i + di, j]\n",
    "\n",
    "            # handle the east and west neighbours\n",
    "            for dj in [1, -1]:\n",
    "                if 0 <= (j + dj) < M:\n",
    "                    E -= state[i, j] * state[i, j + dj]\n",
    "\n",
    "    return E\n",
    "\n",
    "\n",
    "def numpy_energy(state):\n",
    "    E = -np.sum(state[:-1, :] * state[1:, :]) - np.sum(state[:, :-1] * state[:, 1:])\n",
    "    return 2 * E\n",
    "\n",
    "\n",
    "print(\"Naive baseline implementation\")\n",
    "naive = %timeit -o -n 10 time_energy_function(energy)\n",
    "\n",
    "print(\"\\nNumba version\")\n",
    "test_energy_function(numba_energy)\n",
    "numba = %timeit -n 10 -o time_energy_function(numba_energy)\n",
    "print(f\"Numba Speedup: {naive.best/numba.best :.0f}x !\")\n",
    "\n",
    "print(\"\\nNumpy version\")\n",
    "test_energy_function(numpy_energy)\n",
    "numpy = %timeit -n 10 -o time_energy_function(numpy_energy)\n",
    "print(f\"Numpy Speedup: {naive.best/numpy.best :.0f}x !\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7d80c912-3773-4c53-a4e9-3f7ee66e2e0f",
   "metadata": {},
   "source": [
    "While writing the above faster versions I realised two things, first there was a bug in my code! I had written \n",
    "```python\n",
    "for dj in [1,-1]:\n",
    "    if 0 <= (j + dj) < M:\n",
    "        E -= state[i,j] * state[i+di, j]\n",
    "```\n",
    "where I of course meant to have:\n",
    "```python\n",
    "for dj in [1,-1]:\n",
    "    if 0 <= (j + dj) < M:\n",
    "        E -= state[i,j] * state[i, j+dj]\n",
    "```\n",
    "I found this while writing the numpy version because the two functions were not giving the same results, I didn't pick it up from writing the numba code because as you can see I just copied the implementation over. So the first lesson is that simple test cases don't always catch even relatively obvious bugs. \n",
    "\n",
    "The second thing was that even my naive function was doing more work than it needed to! Because if the symmetry of how two neighbours talks to each other, we can rewrite the code as:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "755ec34c-4975-467f-a559-f585f9d6a42f",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Improved Naive version\n",
      "34.6 ms ± 147 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Speedup: 2x !\n"
     ]
    }
   ],
   "source": [
    "def energy2(state):\n",
    "    E = 0\n",
    "    N, M = state.shape\n",
    "    for i in range(N):\n",
    "        for j in range(M):\n",
    "            # handle the north and south neighbours\n",
    "            if 0 <= (i + 1) < N:\n",
    "                E -= state[i, j] * state[i + 1, j]\n",
    "\n",
    "            # handle the east and west neighbours\n",
    "            if 0 <= (j + 1) < M:\n",
    "                E -= state[i, j] * state[i, j + 1]\n",
    "\n",
    "    return 2 * E\n",
    "\n",
    "\n",
    "print(\"\\nImproved Naive version\")\n",
    "energy2(states[0])  # run the function once to let numba compile it before timing it\n",
    "e2 = %timeit -o test_energy_function(energy2)\n",
    "print(f\"Speedup: {naive.best/e2.best :.0f}x !\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5d1874d4-4585-49ed-bc6f-b11c22231669",
   "metadata": {},
   "source": [
    "## Conclusion\n",
    "So far we've discussed the problem we want to solve, written a little code, tested it a bit and made some speed improvements.\n",
    "\n",
    "In the next notebook we will package the code up into a little python package, this is has two big benefits to use: \n",
    "1. I won't have to redefine the energy function we just wrote in the next notebook \n",
    "1. It will help with testing and documenting our code later"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "c134b337-e1ea-45ed-bd9e-45c81e6d4400",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Author: Thomas Hodson\n",
      "\n",
      "Github username: T_Hodson\n",
      "\n",
      "Last updated: Mon Jul 18 2022\n",
      "\n",
      "Python implementation: CPython\n",
      "Python version       : 3.9.12\n",
      "IPython version      : 8.4.0\n",
      "\n",
      "Git hash: 03657e08835fdf23a808f59baa6c6a9ad684ee55\n",
      "\n",
      "Git repo: https://github.com/ImperialCollegeLondon/ReCoDE_MCMCFF.git\n",
      "\n",
      "Git branch: main\n",
      "\n",
      "matplotlib: 3.5.1\n",
      "numpy     : 1.21.5\n",
      "json      : 2.0.9\n",
      "\n",
      "Watermark: 2.3.1\n",
      "\n"
     ]
    }
   ],
   "source": [
    "%load_ext watermark\n",
    "%watermark -n -u -v -iv -w -g -r -b -a \"Thomas Hodson\" -gu \"T_Hodson\""
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.8.10 64-bit",
   "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.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "31f2aee4e71d21fbe5cf8b01ff0e069b9275f58929596ceb00d14d90e3e16cd6"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}