ReCoDE_MCMCFF/code/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](link)\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: LINK \n",
    "\n",
    "```bash\n",
    "#make a new conda environment named recode, with python 3.9 and the packages in requirements.txt\n",
    "conda env create --name recode  python=3.9 --file requirements.txt\n",
    "\n",
    "#activate the environment\n",
    "conda activate recode\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(\"assets/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": "iVBORw0KGgoAAAANSUhEUgAAAOsAAADuCAYAAADYx/BmAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAAaKUlEQVR4nO2dP44sPXLEewQZewHZusTge946e+511nyYS8jWBdYbWQ0UqIzIiCTnW1UpfsAAr6vIZJL1h8xIdr+P7+/vVwjh/z7/9q92IISgkYc1hJuQhzWEm5CHNYSbkIc1hJuQhzWEm5CHNYSbkIc1hJuQhzWEm/DvTuF//OMf33/5y19eX19fr9fr9fr8/IRlT5VBZd+fV6621DoOla+dL9fzSl87u8yntc7qgzM+jq+MzicG8p/VnfR57Wt1zU7cP4zPz8/XP//5z9df//rXj+r8h7Pd8Pfv399//PHH6+OjtGXzbvttb/08KXvtDyq79llpj/l/gso31LbSLuqrU4f5pIzh6u9k/B37Dt34XNtDvkzGtuJq5/fv368//vijLJxlcAg3Yeth/f7+/l9/6znHjmP34+Pj9fHxQdtB5951339VndV+ZWu14/ytNpTxWX1T6qD20GzWjddaBo0T83f9XI0vGi92vVlf2TVXfXLo2kN9ZkvtzKwh3IQ8rCHcBEsNfoOWSi5IHGIBPDqniFMOioCFfHLqsPbQeEwEmuq44kPnP/Nl9feUMNn5ogiWK8q9t5ZlviC7qNy77O/fv6HtzKwh3ATrYf36+mqFmTeKCLLWdYSTtR1WtxMpmMChCBqd/5XAwYQrxClhyRFtVKp+qUIWE7AU35AoyMbLuedUUYr5y6731ZcITCE8ACtm/fz8pLPM69UnkCfxg3puYq9j0g6r242PYn8S01fHlTiKHa98U+qzGNwZD7W9ys/OJ0c7qMoiX1ifOzKzhnATRmrwT89qp9S79zlHLUQqs8NEDUYxUeUn83Hn2nSrJgVnllF8QvanfUbtKPfpSVU7anAID2akBr9R1ExHFV7tOvYUtbNTHJkKOVHA13YVRVoZF0flRApp5VM3fkw1V3DqOmq5qsCytpFKz8bS6TO751R7mVlDuAkjNXii0p6KD1EdJXZVjqtxnOLjRLWt7HT2WZ0T/VGumRJzI1+Z8uq0g3xkPu3En8q121XSr2RmDeEm5GEN4SaMUjdvquWLuhy+lnGWZsrWMOQT86WzpyxlJ6khtoSbLCOdMXU2UKC21eNq+6r/7D5Qlp7oGjl1Hdh1V/ucmTWEm2DNrF9fX69fv35ZQpCS/F9xxCJHQJlsGHCYJNiRjaqMsimi80mt59pzRLSJMMNWSermkYrJBhbkM/PlxEaKzKwh3IQfT904aQaFnXRGFy+q51BZZ7ZZyzjxmxMTO33eiTedNNhuSkttZ7WpoNynymplZwWIyMwawk049rCqW7HY9qpu+9t0axza1lb5gXxj28PWz8g31s7qazWWqJ1qDNF4MV9YGTQObHzUe8L1F9VRxqfyk92Ta1tdf1AZtW6+fB7CA9jKs1bxlaOITmI9ZG8n5nNij9Nx7tqegqP+KnnciVLpaBNOnNbFm9U127nnlHtwEgtPfcpX5EJ4AHlYQ7gJx343eMVZ5nVL26kPO2kex2a3jKxSK8rSsCvjLMmVNM9adzL+TgjjpHCYT91S00kRsZBIPV75u3PPvcnMGsJNOPYbTD/9llZnL0UY6Npl9hlK8lxJjaCyp2eObvUzWUFV9btxmYh2zsrq1GyMfFKur1Kn61Nm1hBuwug3mJSk8IqyAWGymYBRbTSo/hQb6+cukY58RRsHWAJ/TagzG4iq385YrnaU8UFtK+12467Uca7vZEyr/nX+snsimyJCeAhbG/mZcqbEJcrmBFRnbY+V7ZLZTkxTqaldn7vZ91r2VHylxIOorOOvc+06m057O6pqVd+5T1md7l6rjl/HMJsiQngAo4e1WnejNbsSp6ixAWpbLdvFFUpZ5v/a5/W8Gy+7/araVGK9HV+Uvqo2nLJX0Lgr8e0kbkftKe1U1+Xa/8SsITyArV/kv4Le9hP1jr2duxlEgb0ZnTcualtRSJFPSjusz90s373l1RljtaeMpeKTClNaWf/WcUKfHaYzdFX/8/MTlsvMGsJNyMMawk049htMK47Er8rdzH5Vt1uWMOkd2WU2lVQRSxGo9h2fOhtunfU6KimVLk3i9LkqOxmHzgd2zZz7lYU2atk3mVlDuAnHviK3MxMhu0o7jM6+80ZkZZ2kOVsJdL4oM1K3ktkdU+RbhTLrrsedGbtrj9nvWMUotyzywenPSmbWEG7C6GF1ZPq1zvVvkiZZPysyfZeScKT3SQKfjQsatyoFsZNeYP6tOCmciW+TtF7ns2OjssP60d2n1T3djQu6ztkUEcID+LGfdXnjqHmsTBdDsjpOPO0oi7vKqtpu13dHgWWcVHiVdiqc2blr27k+qw2m1k70i9XuZANFZtYQbsLWz7pM8lqKneq4OiOxnODUT2Qfwd7wJ+wzhdfJBXa+MPsKk9kD2WArhUkmorPLVhPKyqArO1GoM7OGcBOO/f+s6A11Yp1f2UExhpNHZLNx5xvrB7KllHH8d970yspmJw6drIbWdpi/Tu50J+50Vg4T+91sny+fh/AA8rCGcBO2HlZlg4OSAF/LKue6dq84CfZJMnsngc+S86hNRrehobpmqJ1qjCegutXxyaaCtW/ONWJ9RePDrg+6nutn9JdNESE8gNFX5BgnkuWs7EQgcMQKVEZta8pE2HDEqDdMAHLSDCfSI4qfTCxayzrXDPk23QyxlkX2OlEzAlMID2D0G0ws1uvW++wc+ly1M0GNhxhKTIlivy6ecmPgKubZiS1Xu5P4jdlBfld1JlrHZHyUe8/RR7p7o7KVmDWEhzHabvimin9QTMDW6l0csbZVocQaa2zDZtVJfPJnU/VDUZ9Vu1MfTth16qL4mcX0SnyO7ChxLopVd8YgM2sIN2FrZp2qqJ09ReVE7TC1c6JKMpC/aCZndSofuz7uKpddO5WvE72gG6fKlxMrhIpOzZ6sVhT7yr33/Z3thiE8gjysIdwE62F9b4pgkj6S8plkzVIdSApH7Sl1ulSCmwboUkBKmqgaixXUHkutMB+7FBFLPSEbLBXVjQHzzanr2EfXuUvRTO13ZZO6CeEBHBOYOpHFqVMxEQY6u7sCTSciKIKZMwY7KKkPdNxJm5z230nr7aSPVhuOrR3RbbUXgSmEB7C13ZDFrEp8xc6hmHgS/6C6LCZx4sOuHywWXuvs0sXarK/dWCsxN6PTKpzr6tRR7jlWB40P02FQHyv7V3uJWUN4AH/apggl1lDiks4XZXZyYjCFLl7fsXW1p5RFdZyYD+kBrt2TsOvSxf1MJ0H2qzqo738WmVlDuAnHv3z+ZuctVL3JutnKUW/XdqqZw1GQd9RAZzWhrFoms2TnfxerIroZyVF2J6uhycqj+uwo3juzbdTgEB7CsV/k796erE6X+6rsKzZ2ZozJW9TJs6o2mX1WR9EDuj5OcppKX52ZulspVHWdMmh1odxPp1YGKplZQ7gJeVhDuAk/nrrZSRmcWrKhMsym4hMSvdSkfmWj6nMnLE0FJuSTcnzivzKmK8r13kkXKjgpRuSTMk4dmVlDuAmj7YaMdbsV2wqGzlVbtVbY9i20LQzVrehssD4z+5M+O2O6+q+U7Xyp2lTsIn+VfnTXkt0b7Doo12j1cfUJHa/OoT4y/xGZWUO4CaOHtXojoDchOl7ZUd9617Loc+cvquO87Vjfqnaqvq022AzqrCZQHTajKT5218jxhTFZDaEVQeXDZMVxArRKev9lI38ID+DY7wa/Wd9WJ1Q3x4cdVbjCKYvarRRkVqY7N1E017qKvxMVeJ2RlXZU/xCKL6isc14Z767Pjq2VzKwh3ITj2w2dN8j6JlRm5e6NuAvqhzP7KiizI5oxmE8nfWHj382win228nBWWZPVzsrOirCq140T8iMb+UN4AHlYQ7gJWw/rNOG9llVSEigFwUDpBpbOcJPb3bh0qQPW965f7JySVlKvleITa1/xG/m2flb6w64v+nPGdD2u3u/Il+u5pG5CeABbqZsKJ62ARBxFWHJEl84nJ52h1Geii7LaUH1xRB1HoFFniM5u194krbfLThpvRREHV7vd9Y3AFMIDGM2s1Vuie5NXb5RupmOSPqpT2UcyusMkwV7V7VIESjqGjTXqmzMTTtJUji+sHaevHROflHKKD869rZKZNYSbsPXrhqfjiJ2kv7OBgm3G6OJP1jZKeCv20efKb9Q+O8Zi+bWOE6dPZrxJXaeMo5es9wDbtICuFZu5lXtbHYfMrCHchGMP6yQv6eTdUB2UV3RygtMyXX02Bugc649TZz3X5fdYrpGN6WR8UPvrCoTlSh3UfGrHOpaVLeS3cm8mzxrCQ7Bi1q+vr9evX7+oGryi5ByRWqus6yf5z/W8Es9VNjq/ldymkzNFMJ+YDTVedvKuVZvqvXE91sXc03hazbU7mY6qnuPvtZ3kWUN4AHlYQ7gJx1M3zlIQ1UWf0THXl510yal0Elr+siXhxBdm30mZIZ9WdlI4zI6TjmG+qOPvpLocnPtnJTNrCDdh63eDmeTupG6QvK2kJBT7q7/os1LHaYf5j8ru9KNqU0kfOXVU31iqSMGpO/Fh0tfVvnJ9nfaSugnhIfx4zKrgxE4oLnFSHopvzO5KF7PupqC6OLRrs6vT1WWw1JMzhmsdB2d8UHtMx+hSN4620vmS1E0ID2D0sFbr7m7t78Qe6/reiUGqNlmM1/3txHFVP1DfUF0WpzuxvRJnTXQA5tOKEss61wZdE+VadXWr66eMzwoaa3RfJmYN4QFsbTe80sUJ7I3E4p6ToNjjp+OsKmZFVL509p3YjNXvYjN2zonPJyh6hhNTOnGoMi6dfeY38/NKZtYQbkIe1hBuwtZ/n3EFTe+nUjk7y9N3nZ2lreOTkk5yfFGEubXsJH2htIuWhCdSXawdBWdpjnxSQo9pWgbVvbaZ1E0ID+D4L0UgudtJfTjyOUt9rMfWz6ydLpWDRDZ2XPGpAsn/Ds74KzjX1UlTTe6Xrp3K724smR2Wrppc52tbSd2E8AC2thtW8Vv3mb0t2bnOnhNfKe3sxKFOmoThrETUdnZj7g4ldTO5B9j1dlJbap/RCqkrU7V5isysIdyErf/rxklMV3XUtyg6NmUS902U3VMKqaPWojLKtXIS+A7ORgp19mIzH2q/srvTZ2UjyGTDCSIzawg34fiXz9dzVVmk7DKVTSmDVLaOqfKHyiiqZ6Uud8oosluBfKr8V5V2dEwdy843pW/KdVYyBGufuyyAcm8wv9EYr3U/Pz/Lc69XZtYQbsOP/f+sCo5aOImfUFzixA1OzDdReJ2YSYn5VrXUGdvTCubJWNtRaZGNquwkZnVUYeWaqWRmDeEm5GEN4SYc+6UI5dzrxbdmMdEFiQaKmNAF/YqfrEwnRKjjsNpC/jLxrusz81/pn2tDab8bZ7U/CqpYpAhMin0kqiIRK9sNQ3gAW9sNryBB4I0i5kwEGtT+qXYmAkc3Fqw95kN3nPnA6Px0xLV1VlHsK6KXIpShug6KuMlEwiqVU5VF1yVfkQvhAYw2RVRr+C52UWIAJ4bs4l4WvzE/uuQ1ixnRWDgxeIUarzMqX9Q+T+2r98QV1QdFm6jiwx29RBnvTlfo6iRmDeEBjGJWJV54o8RXOzEkavc0p+2jmGx9219R4sQufq7soc/dcdafyk9kj/V5536a+o3KKD4697Tjw+uVmTWE23B8u+Hk7TyZYd+wt5KqKDqqoVOW9aebNa//nijgzgw7mQWQD2yW7HxTfJnMmpW/XXuO2sxWNmh1MSEzawg3IQ9rCDdhtAyeJNqZAPFGKdstmZ1lHltSIV92xbCuvrP8UsapO87sOaLgzjL7tFiklJkIlM74IKq612PZFBHCA9j6pYgramKdJaSRLZb4VtrtktlVIhyVUUA+uyJD11fWL5TYV8qgz6zsif5VKNd39Un5W+2vx9l9ytpD46Fcq2yKCOEhHNsU4cQwKkp8y+qsdbt4qzrXtcfsTdJKTmqIgfqqpFaqsm4712PKuKM2Hc1gJ76dpPEmMWsXnydmDeEBHNsUMUncO7NKN1PsbGxwZoWpmq36i1TC6+dTswwqw8bHUZnf57oZ1VFikQ2ViVrerbZ2Vo8OmVlDuAnH1OCJioqUxfW8ogpX7XT2FJj6iRQ+RUVF9itlEsHKqep81Vfkc9Un1tcd0LhUvq7HWL861Zb1Q1GU1b5XdaIGh/AQjv0vcm9QfLWeZ3VOxQCT2KjzAc2kbjvMHsLxX/18tXtKiWZtTevs+FjpDDvtTFRtRW+IGhzCQ8jDGsJN2PoNJiYmoMCdCSgsOO9EF1WUqdpjAoGyTQz5XYlGnQChCBvsr/ONCUiKOIL644w/6x8qy3x0xlRth9lR+orGqfL1Wi4CUwgP4Mc3RbDza7C9I3CwNytLb6x1kVjEhIGJ34pdVahyBBQHpa/sGjoiHaqjiI6OwNeJQ6zPajqtsqf0IwJTCA9h68vnyhty8uZldLNxl1qqjk9npu5tOZntnfac8Wd0ZZ3UhAKr68S9qI5zzzn3oDLrd75U5VQfMrOGcBNGM+vOG6Wqo7zt0Nty8oafvNnZW7rzZbpJQo13lNiSzZ7vc2sZZUbtjis+OddbYccXxZ6iM3T3q/IcrGRmDeEmWA8r+1/kECgP+Hrh/BjLH651WT4L2WU4eVaUa3TqsHbUssp4OeOvgHxQQONV+Y/GwPGN5ZRR7tS59xRfWJ3rueRZQ3gAx74i92Z9A7I3YvdmYrMkKlu9ESdv5e7Ny+pM2q3KdjPq6iNrW7kOXT/YDMTG3xkftR+sz8oshqiuv7JiWuv/BJlZQ7gJeVhDuAnHf92wEkkQXRpASWOw9IaavlCX2Suor07qoFuSKX4qYzwZH5Y2UcZSHY9TqRaWHunsrP1xlrKTe7y6t7+/s90whEdgzaxfX1+vX79+0TLdW7J6oyh1lbc9Oq68cbs6zL7zVu78V1Dq7vh/amWgrgTYbI98VM5Nxkm5B9m9OFn9qOOemTWEm7D1FTn2Rjmx5ldiDkb3ltudQbqVAZvllXZOxMSrb5ProqafVvudJlHVWcsoKwR0bnL/7GgsVdvKeKvjm5k1hJuw9euGV9S1/5osV1Fn393Y4ERC23mbOm9iZ4ZzbHQ+KKq84qOi6O7EnxP/1zJsfHbi9BNaRWbWEG7CSA2eqGDKW3SS63LiEyfntZ5TYsoVZZxOzJZT7WCS/+xmJBZ/rp+dmcrJ5ypj7Ki2OzlYZUaNGhzCw8jDGsJN2BKYlK1lJ5Z5ih1l+YWWbM7SVhFbHEEC+cbqTzYVsDZRe6dSXJN0xqTuCQGxasdJg3X3QASmEP4fsPWL/Fe+v/tv/3981N8TfH92qHyYUvVJ8Wmtg3xi46H4sp5bfavGd+3HZIyra4rsVfbXOj91vVefnHsD+cL6iO6Vbjw6Hz4+8ksRITyCY5si3jiSODrH4oZVnnfiKsXXLt5U0iRK+sTxn80o63lUVonp17IoFcJAK4LKX8Xu5H5i7aix4zTdNtEZrsfyFbkQHsCx32DqYtWqXLfOZ7HWJOZzQDF45dOkPSXe7eJoph2sZaYxFPJJ8b+Lb5U+dz6qSjW6jsr4sHu4G/9JfI7IzBrCTfixr8hVZVR7J/JkTh0l/lFtVWV2c4+T+Gc9p2gGOzMAa6eLuSdU995kNTXRFxS/TsS5K5lZQ7gJx36RH8WqLM/XxSdO/Mni5C7XqMShLGZ1fFvbRDFUdQz5zcaJjQWyw1jHEvWjm+2rsZ7qDF1dxacTsT26juzeXH1InjWEB5CHNYSbsCUwqbI5KosCd0fwYEJNZ8fxTRFxdpiME8MR/JDIUrWL7Cplu3aQf6y9iokQxK43Gg8mpjnCpHp9M7OGcBN+7HeDFQl7MoOudnZkdTYboLfpKswobbOZg72tO/vKmx31WcGZ+Sr73SzzU6sVZQZccVJnbNZE10y5zh2ZWUO4CT8es77ZiUNZm04s4GwQQL44b8TJZgWnz1VZxZ5aVolvGV3ZnfjZme2dmNv1F9V1dACVzKwh3IRj/4ucGp84MxMru/OWc97SjkqIbCgKqRKzKjPsidlMQZkZnNXDatdZQSl2kZ2TSnvVzs5MvpKZNYSbMPqKHNtm9QZt26q2fqEtWNeyaNscY2cLG2rXqbt+ruyhfl3LOr4g+04/lLJXP5XYXS2LfNix4dSryqG669iy54HVVcnMGsJNyMMawk04LjC92QnCFTFKEQqQDzsy+mQjQuUn+uzYr+p06apJKmci9DH7jk8nRColVdddF1Zm8jwgoTK/wRTCA9j6j6kY6G022WbFcN7SSrvdm1zZWvbGSfcw36q0DmqfzbqISbpkUkdZGXRlqzrOjKeuppR+oM/X+jsrkZXMrCHchNHDyiRnJF1XMvpOagX5okjiSGZ3U0RKygmlAZDkX6GmxVhZ1i4al8oGwrm+q0+V/e46MJ9Y+qTysxoXRjdu1TnUj9VufikihAfwp/0i/446fD3ntMtmqut5J65T7EwUV6esM04TxVWJsyYKPvLpXxGzVnZQOeSvc73ZdY8aHMLD2HpYlbU6ikVQTHGFxTZdrINirlOxctXHrj9XXzqbStyv9GMtq8SHyAYaZzX+V687Gxd27Z2yyJcT15CVrWLk67HErCE8gNGXz6u1OlrHr593YqddJnadGViJxSbtIHssZkV1nPiQjZfTRyXuXO2iz53NyjflerM63b2sjIHqQ2LWEB5AHtYQbsKx//JxZSKGKMlmRxzqBAjmAypTiRadcMKEMSaGqX1mopQiyqnjpS55UZ/WMsp4OCKag1PH8UkZQ3bvRWAK4QGMBCZHcHAEjoloobwdJ2ma1e4pn06kjBQfFCZCGLIxoVstXO0zoXJyz72ZjAGro6wC1nLXvkVgCuEBfDhvxr///e///Xq9/uvn3Anh/z3/+be//e0/qhPWwxpC+NeRZXAINyEPawg3IQ9rCDchD2sINyEPawg3IQ9rCDchD2sINyEPawg3IQ9rCDfhfwAiwuC0As5JHgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 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": "iVBORw0KGgoAAAANSUhEUgAABG0AAAEICAYAAADlbwmJAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAApy0lEQVR4nO3dMY4ru9UgYGngwBv4AWcTegMNd+bEO/A23nq8Da/gduLQUOgN/JmBmcyJM00w0P3letThOSTVze7+PqABvyoWyWJRukbp8PB8vV5PAAAAAOzlf310BwAAAAD4NS9tAAAAADbkpQ0AAADAhry0AQAAANiQlzYAAAAAG/LSBgAAAGBDXtoAAAAAbMhLGwAAAIANeWkDAAAAsCEvbQAAAAA29JtK4b/+9a/Xf/7zn8/qC/AOfve7353+/Oc/nz+6HzP+9re/XX/7299+dDeACf/+979Pf/zjH7/Ed9Hlcvl57OXl5XQ6nf7j2PFc6/z9uUirrazbta1+VOvKtHNfb/ZYr75Ipo7oGWTrrz7bkb5lys30IzserT5Fc3fFs4pkP2tR31rHfv/73//fP/3pT/811KkN+P9F8DX861//an4XlV7a/POf/zz98ssv63oFvLu//OUvH92Fab/97W9Pf/jDHz66G8CEv//97x/dhWn/+Mc/Tr/88svper3+PHY+595D3a5plY/OZevK1pHpR0+mn9Uxui//+vqaKpepv1pvdlx6Y3+8Zua5ZK47ndr3dzsf3Xuv3uO/vzPzruV2bWvOtK6L7vPerd+tvv348eO/ww5vzv8vgq/h7e2t+V1keRQAAADAhry0AQAAANhQaXkUAAD/38vLy8NlHr3lQNHykWh5SKuN6nKS+7qqy29a11WXelWXc2WX3GSXRUUyY5pdMpVdutW6v8wcmOlbpFc+M3db9UXzrtdmdD47fq26Ruc/a0RLCIH/IdIGAAAAYEMibQAABlwul9Pr62s6Me2j80dRdEAUTVBNhttrsypbfyYKYiQBcLV8JgKkZ8XYjyacnkmSnDUavfSsRMutNm51ZSOERHR8vEz0oecE/0OkDQAAAMCGvLQBAAAA2JDlUQAAE3rLVDJLQKpLf1rXZJOuRm2NLDEZTSqbTUT8qJ1eP3qqCZyzbWbqzd5La4xG5kpGdVlX1PaKpS29ekf79l2SD48m9c5+pqqfm54VS+bgqxJpAwAAALAhL20AAN7B9Xr9+Rc5n8+n8/ncLH9/7FYuquP+73hdZgvtbJTQ8Vf3zLWte7k/dvxrXRu11RqD7L20yrfaHO1b675a5bNzZuaeo7616qj2MTOfW+23jrXuK9PHYz+jufUZXS6X1HOvfO5Hyo3Oudlr4avz0gYAAABgQ17aAAAAAGxIImIAgEWipRbZxLvZOo7HeomIW8dGE4hmk5xmE5lG9Wf7lqlvxVKY7PhlEz5HcyGaMyPPtjoOmfsaeVajiZmz4xzVkb32M3l5efm5fOwmm/A7mlejx3pzYjQZOXxXIm0AAAAANiTSBgBgwO3X7ZZnRdVUt+N9j1+tq/cXbWddjcypjt9MuVb5meeXOTczZ7IRPJHMs61G17SOteb1TBRGda58dpfL5fT6+vofx7LPvZXEeeS6yrEWETbwmEgbAAAAgA15aQMAMCDaZre6vXdvq+Zoy+PWuaiOqt4W10e9LaOjLZdb9UdbOmfHL7qv6B6y9fbOZ7bEjq57tDV4dZvkFc/vUf9XzbvW/DjOv9acyY7VzPjtKor6uxdtez7zl2kz2zfg17y0AQAAANiQlzYAAAAAG5KIGABgQGub3ZuZ7YVXJKut1p/dinpF8uMo+Wy2v61zrWU5kez2xJFoK+psuRXJlKvJoKv3GSUKnpmTM0mVj8e+89KaViLi97Qi6fVXWKYGzyLSBgAAAGBDXtoAAAy4JSLuJUyNEmxmE+8ez1WTAz+qI+rH8dyj5K+jbUaJYHtJZVeMaSYJbfU+e/Uc2340BlHi3Wx/M4l3RxL0VpMZt45lxqN379FcyCZE/irJb1uJiD/i3rJj/JWfBTyDlzYAAAAAG/LSBgAAAGBDEhEDAAxoLUmIRElio0SvrfI92aTAlT4+KpdJCDtzL63ymfvKJhieSRBdTaBcTRT8Hs8vkrk2s8TskWp/swmRVyTs/kxaiYijpMDvuQwpm5w4mleWTfHdibQBAAAA2JCXNgAAA26JiHuJW6PkxJmEr4+Sy0bls/Ue+9gqV72/XlLZqFwkSqA880t8LxlvpY/Z/vbaXyF6flFS4JlnVH0WUT+qdYxYOd4fKRv1917Jfkeeq0TEn8+zvrv4NS9tAAAAADbkpQ0AAADAhiQiBgAYEC1JiBL13luRpLVlNGlur81qQtgoqXK2/agfLTNJdqO6Wn181nKOaCyriXpnEjjfVOdpto1svaNjn/0cfhe7JCKeKcfHy3z/Wuq2lkgbAAAAgA2JtAEAWKS6BXQku9V2dC4bmZDdlrcarVCN4MnKJjaNZCJRRn4tzo7lsY3edc/eBjy6tjqPevVntolv1Rsdi7ZOf9TWV47uiL4Dnh0tlm0zOpaNMmPOyPzIPKuv5KMiLu+JtAEAAADYkEgbAIAN7fyLZTUXyWiOnV65qmw0S/Xa45bZj+odzVXzqP1Z2VxDrWsyv8A/OpY51+vPaPRFq29vb2+lOj6DasTXR/fjeGzn77+voPpcqt/572llhGtv3o2O0SyRNgAAAAAb8tIGAAAAYENe2gAADLhcLuktlB+Va527Xq+/CrNuHbtde/u7lbler2GbLVGb93/V+3vU12N/j23el2vV0Tp2rKvV5gq98Wh51Mf7+qLyj5YiVfsRzY+oH9k6qnMxOteaf9Vn2Su/cl7soDcnqvMlK/qcZfozMpcZk/nMZX3E84u+11v9yH5nVOfuSLkZXtoAAAAAbEgiYgCAAS8vL82oiaNo2+To17leMt7RBIsjSXArbfa2Xo6umUnKe7x2pP9RHdXxriY1nnkGM/ecMRONkkk23Lv36lhm+jNy7Wc3er8z41SJinpm/S27Rfas/vxWx2bnz0M2qnVlvTuNh0gbAAAAgA15aQMAAACwIcujAAAm9ELao+Uh1eUvrWury5Gyy1Si5VSta2aWWq1c2hTV31tykxm31rPtlXvUx16b2Tpa5TPL8lpLzqr9btXfWwpYXZKQmYurlxN+dtnnc5NdPpk5lp0TUX97z3NFP1aIPjczMt+r1bpa9Y08q/eS7Vt16W/L7uMh0gYAAABgQyJtAAAGXC6X0+vrazkaoqUaRdIrU43eiOroHWtt+3wsn/11ciaiI6rjeN19udGx6smOw2hS4mpy5/v/ne1PptyKX7dnzERQ3K59e3ub7sdnkPm8ZD9T1WPV77hn9aN3LDMnV9aVrW8m2qP6XD46Km30+a1OmC0RMQAAAAAhL20AAAAANmR5FADAgJeXl4dh5CuWO7Xqy15XXaYyk4T0eC5bb9boUoTscoKZpWyj49yqr3qfK5JjVudTq62RJXWZ+Vy9p49e0rGb6lw7lpkpN/IsntG33vwe/V6d6Vu2zMrEyVHfdkiym/ERSYEzy4Hfox+nk0gbAAAAgC2JtAEAWGz1r3ErIlai656VELYaZVH9dXkmkijzy+2KiJHW+V7yz+h5V6OiWm2siKyK+t0rFx3LRmxVn1W2T1/ZigiQlvcaxxVRcb36qtFzo32bqW9FhN+KeqP6V4/fs5ICjya1b5V75udApA0AAADAhkTaAAAMaG35XY14qJR5VC7zi+x9fSvyQlRzO1SjSEaMbmG7OgLqeG6kjUx/R3L9VLdVzm71HalG8MxENI328SsajXzbZcxmvkei+rLjsaJvWdXnMhOJeKxvZjxWRoDe2yV/zU6fIZE2AAAAABvy0gYAAABgQ5ZHAQAMuG35PZLsdGXC1ux21tmlPJFMWyNLbjLh5CuSkPbaqSb/jO55xZbIM0uEqqLtbVcukxqxIml1VNePHz/GOvbJzCRbfW/VxMGr28gu3VvRt8y175HoeDQx88h3/mjf3muZVK8f0bFnLOsSaQMAAACwIS9tAAAGXC6X0/l8Pl2v11/93Tufz6lonNtfq/zt2P3fsXzrryU6P3Ivx/5k6+j1M1N/azwy9/noXo7j0DpW7Xfl1/ZHz/jRr/7ROLfqjcq1ROMQncvWm51Prbaje4o+GyP9/cyy8+QzWH0Pmc9Nth/ZvlU/I6tF9Wa/RyrtvMezepbsZ2d0HlV4aQMAAACwIS9tAAAAADYkETEAwIBbIuL3sCIRcZREM3Mua0Uo+Egdx/HoJZTMJBXtJTPOjlc1wXHrutZ9VfuWbStzrjV+0djPJHbNzPHevUfXfudExJ/Jiv6PzJNMfTNJh491Zc0kvG31e+X8qH5+K/VlrPh3K1vHe3yuRNoAAAAAbEikDQDAhN4vdJkIgxWRBi3ZiIeofCaJ8kyb90YjQXptVKNCWv8dPceW7C+3o7/mtuZYr83jM+21HUWxzPy6XI2KykTrtIxESn1l0edgl7GoRpLNPM+V4xH1rVd/NjHxo/ofHasaHY+ZflTHKDo20o/oWR3LzPRjlkgbAAAAgA15aQMAAACwIcujAAAmrFja0Quhziy1GVnKkwkNr4aczyz1epZen27eq2+98Rhd6pBdVja6RKMn2+9oOUGm/uq5R+Vu7b69vaWu/ewyn4OP/qy2+lE9trKNkfGIloxW63rUr8qx1W1Ur3uve32Pe/+oz5BIGwAAAIANibQBAJgwkkDxaCTB5uivitXrViTM7LV/TEg78ktlNfnns+oYjUKqJhbO9i16HtltxqtWjHOvvuw9R+Vv//u7bPld/S7azUds191r89njNfrZnim34j5nylc/76Nt9uqtfuevSAb9iEgbAAAAgA15aQMAMOByuUz/2nu9Xk/X6/V0Pp9//mXK9/7u6zv+3Ze7idqOyh+jFo5/UR1Rva1zvXGI6mjdX1RHS6ve1nWZsewdi9rKjtGxrt6zysyPVp+y99eq41im96yyYxSNQ/W5fzfZ76Jny37vPFuvzY8cq+x4VMvdy3wPrxZ9Rp/1uR39Lm/17Rm8tAEAAADYkJc2AAAAABuSiBgAYEKU9PR0qid1jM5Vk65mr12Z6LXVZjbRbHZcovIzdbSW9jzq60i5bD9aovqqCTBXJMzMJlW+mRnLzLleva1zlka9TxLVVXqfm9FlStnvpJVtVmX7UT224vPwUXPmWUuzdiTSBgAAAGBDIm0AAAa8vLyUfpWrRldUfwFtlc8ey0bVZKIrMskuj/Wu+LU/U0f2V/Hqr+fZZ5WJ7rk/Fj2/6Fk8cnzePZnImV4/qpEzKyMcstFL32XL75bjOO4YibRiro36iDZH+tGSiXjb5RkTE2kDAAAAsCEvbQAAAAA2ZHkUAMA7iELqq0uQHpXLnJsJh88uu4razJSbSfa7YqlXVfX+WqoJn6vPoNVWby5U6x1ddtVbbrJyCcd7JY79rEaSZL+X9+zHLvccaT2rHZe3MU+kDQAAAMCGRNoAACyyItokunYmwqXleG0vmmU0qexMVEg1cqa6jXQ1iXC2b9k6Mm1ny/eeXyZpdHUL+5koo5lop2N9I3PsVu7t7S3V5ncgEunz8Ky+D5E2AAAAABvy0gYAAABgQ5ZHAQAMuFwup9fX1/84tiJZblRXNklodanQzNKVqB8rEtlmkzBnlhfNLIWKys0sQ6s+q5mxrCbArrYZtV1NkDqzpK7SzneU/Zyxj953v+f3tYm0AQAAANiQlzYAAANeXl6ayVCPf9frNbW1cqZcpNVmr95qm5l67/sR9bOqVW/v2PHvUX0jz+x0isevNw6P6pp5ZvfXZtrOjl+rT9l7j+qIrl3xmegZnYtf0bPH+juJvndW1Bkd42vy0gYAAABgQ3LaAAAs0sq9EeUqmflle8XW2dWtn29mcotkcqccIzoy5TL5gVZszT1SLtqeunVsdOxXbK2elclldF9uxTydyeEhGoGPUs3nlJX9/uPzE2kDAAAAsCEvbQAAAAA25KUNAMCAy+Xyq4SpUSLg1YlVjwkoe0loW337iCSWo8mMe4mCV4736mcWJeWN2ooSjWb71nre2ecelYueQe9YdK46T7NJprPtw4hnJQXOzH1Lo74+L20AAAAANiQRMQDAgNaW39nkr89I+jqSyPZ4bS/RazX5a6vNlmjcqkmHj3WuEiVErpbLjsfqe87Mz16S5OhcNWl09tlm+p2dp60k4T9+/HhYL2RE829FJEz2u5mvSaQNAAAAwIa8tAEAAADYkJc2AAADbomI77WSqGb0Eu+2jh3r7yWyfa+EyK1z1aS5I8cetZ1NRtsat6wo2e+jZNHRPWT623q22cSkz5oL1f7e9Ob6yjn8nvfO9xHN4Va5bH2Z+h8lZ+fr8NIGAAAAYEMSEQMADIgSEfeORWVGk/3eyyb0PbY/k8w4uq6arDabmLalei9RHVm9PkaJdzPlW21VorhGZZMvP0Mmkgp2kU0KPJokvloXX4tIGwAAAIANeWkDAAAAsCEvbQAABrQSEY/KJmJ9dM2jhK3Zdqt9q9bVS5I8WkdUbzY5cFZ2nKOxj8pXEiaP3lM0fpnEyNm5NZIEOjM/ojpGEjOPzgUYFX1vt8rB6eSlDQAAAMCWJCIGAJiQTR75nokqo2S1M4ktVyTozdTVu89MwudMZNCjuqLnd4wMqtTROjeaULjX5ujYZ5NpZ5MqR/1ZkWR1RbJmeIbou6t6zFz+3kTaAAAAAGxIpA0AwITV0Sk3rV9do63BZ/IfrMidkI2IqUZcVPs2syX26Fhmy1f7NhNtks1plCn/rLm1cqv5atv3dby9vXXrgqxMFOTIMb4vkTYAAAAAG/LSBgAAAGBDlkcBACzSWv6yItlq5tpqQt3W+dUJZFvjUU0ivDIB58jSraPes60uscoumYqe1cizz7bTqz97n9nPRLZPVVEi6R8/fkzXD0czSw4lHuaeSBsAAACADXlpAwAw4OXlpRntMvIL6fl8/vmXLXcs3zp368/9X6u/t79eP1r1rvDonu7PVY9FY7Wqr9EYtMY+O2at/lbHu1o+mjPRGI7MhWhsoucWtdUrv2JMYZXoM2BOcuSlDQAAAMCGvLQBAAAA2JBExAAAAy6Xy+n19XV50uEo8e9M2HwmSXKv/moi22pS3lY7rYS3Ub2tc5ny9+dHkxQ/Uk34HI1zNrlp1EaUDLpVR7Z8ZOb+ojZm5hbAZyDSBgAAAGBDIm0AAAbcEhFnoxuyWylnI0WeHTEwGsHQ2lL5PaMbZrbazvRzZOvqTGRJa9yyda2IzqpGGWW2bs/2rVd/NQpodLtzvqcVEZLwTCJtAAAAADbkpQ0AAADAhry0AQAYcLlcTufz+XS9Xn/+3ZzP559/x3P3jmVGQvFvddxr1VWtP1u+Wm+rv4/KPCoX3XOrP1H5qK3sc7xvt1XXfT2Pzh37/+j5rXzevTaja1qO99kbt1Y/WnVF5ar1Anw2XtoAAAAAbEgiYgCAAdlExJHRrZ1b5aKEr4/KHY/NJHqN+j2TVLk6ptHW5r26juWySYd791LdKr113TO2fc8+l2w/somFj3VV+92zSzJv9paZf7vMlxVJz/m8RNoAAAAAbEikDQDAgMvlcnp9ff2PY5kohFb5nujX05ntnlu5UTLlo8iLbMTISH9H+xH1rWpky+9M37Ltjj7b3jXZ59LKnZORjc6Knm3Ut0yepEp/+foy3xUfPV+ykYMrvtvYl0gbAAAAgA15aQMAAACwIcujAAAGtBIR3zwraWQ1+euKkPlqkt3VIfvVOrLLfLKJhTOqyZeziX2ziaFb/cjWG9URiZZf9ebJ6PjOjKnlIzwymjx8tep3Ujb5PJ+fSBsAAACADXlpAwAw4HK5nM7n8+l6vf78a7k/f/y7OZ/PP/8irXK3/27VG7XZOt9q6/78ffvHv+q9tPoRnWv9tdqK7jNbb+S+zWjcqvfZ82i8K32LZMs/Gu/efM72LZpr0TPK1h/NXfgIrXmY+Tz2vvP5Ory0AQAAANiQlzYAAAAAG5KIGABgQi/xYyYhbGZJziMrk2hWkw6P9K11bnSJSlRHr/4ogXO2zZHzxzazz6q1DOt4LpuMN2qzNx5RPx6V6bWZaadX7rgUK1Puduzt7a3UN76+7LxeWX/l/JHlUF+bSBsAAACADYm0AQAYcNvyuyUbfROVr9ZxLxvhkpH9xTcbNZSJcJnZdjprNPpm9Za62XGobuGdmUfZaJ1sf1uqEQvVe8mWjyKPfvz48fgG+NaqkWozbWS3tOf7EWkDAAAAsCEvbQAAAAA2ZHkUAMCEFctlsglkH11zLJNdwjK6VKl6z9XktiuMJIhemTR69dKp6nKg0bbu+70yoXVvqdJRNTn3yLI8iYg5ipauPqP+3jHLojidRNoAAAAAbEmkDQDAgMvlcnp9fR1KFptJCHsvu63xo+setRWVq5aP7n3F9t7Pkn0uK9qY2fZ69HnPJJ6uRvdUI2eqcytq/9lbNPM9rf7uykbFSUTMPZE2AAAAABvy0gYAAABgQ5ZHAQBMGElEfAx5zy4JuVddYtU6d+xHb4lJaynKiiVF1aSyK5aGjS7xyt5fdplRS7T0qDq+UR0jS5ai5Uit8hnZJYbVeQ2zeomDn7VsKUqIbKnU9yPSBgAAAGBDIm0AACZkf/XMJvvNRt9Uo3Wq0RszkQ6P+tjrW2QmIXJ2S91MP7JRIVFb1YibrOw8Go3Eypbvlcsc643R6NbMrf7++PEjdS1fXxSNtjIpcPX7+ln94HMQaQMAAACwIZE2AACLZXOntMpkI0Wqv8RW88Bkt3mOzo3m+MlugR7V0VKN+On1vxo9Va0jKt/b/n3lFtgz+TRWRHiN1n8vauvt7a3UJt/He+a0ib5/Pyq3DnsQaQMAAACwIS9tAAAAADZkeRQAwISVy1Dur80mYs2Wz5Qb2RI7swX0TN+qy8WiMr1EvaNboI8uA3vURqatapu9Yyvuubp0ayaRdGYpW/ZeJCLmkZGE5hkrvq8lJ/4+RNoAAAAAbEikDQDAgJeXl9JW1SsS6o5GLfTabBndBjz7C3K2TCuSZ0VkzmjfqhFQj45l2qwm732PrcSjY9l6VyS0ziTzHukb3Mt+plbXlym3um/sS6QNAAAAwIa8tAEAAADYkOVRAAATVoSojyy5OZ5bkdgyuu5RP1befzUZbsuKJUiR3lKd6hKrlc+lV+5YvjdWK553dtleph9RvdVncF/u7e0tVZ7vYyZJdqbe3nfA8dh7fA+zL5E2AAAAABsSaQMAsEg1aqKaiLWaaLa19fJM1MmKRMQt0a/KmfLZdkciS6I2q9tvR/3Inp9JqLtye+LsHHvWlsjRFujVbZJt+U3FSFRXtr5Hx3plJNr+2kTaAAAAAGzISxsAAACADXlpAwAw4HK5nM7n8+l6vf78u2kduxeVP5/P6aU+j/5udbT6cTt3/9eSOTdSR3Ts2P/j0puorcp4P5IZv6z7/rbuK/ucZ7X60ToXHYv05lh0rDW2rbGP/np9evTZaI0HRLLzb7SumfnN1+alDQAAAMCGJCIGAJgwk3R4RRutRK9R4uJWItbMluKPZJIIv8cW18djvV+fq4mZo35UkzaPzIVqQt/MffXGNEq+HF1XNTPvqv0QlcBnYr5yOom0AQAAANiSSBsAgAEvLy8/c2TcZCMpoiiBajTGiq2oq5E5kWw/RqIrKnX07qVlNPom+/yy21OPbpPdE20jHI1l9lx2zDP9yM6P6vbetkkGPhuRNgAAAAAb8tIGAAAAYEOWRwEADLhcLqfX19duudZymWri3+x21Zn6W6KlLjPLtKpLV6KlRdVlLTNJdmeS2laXFGXqap3P9ruaKLuX0Doa++xcH01e/eh8tkyv/wA7EmkDAAAAsCEvbQAABtwSEd87n8+n8/l8ul6vP/9uWsda5Y9l7stFbmXu67s/1qqj1dZRq2/RvfTqaPU3476O1nXHY636q21m76V1vjX2reta5bP3/OjeHx2L5lhVtt/RfI7qiMo/iqrJzPVeOYDdeGkDAAAAsCEvbQAAAAA25KUNAMCEkeUm1aVQo0umqv3piZbh9JbEHNvKlj+282iZ0QqZe+ktr6mOR7T0J1ry1luCNzonW31rHcssZ3q0NK26LC+qt1U+GoeZMQL4CF7aAAAAAGzISxsAgAGXy6UUERP9teroJaHNRG9ky7f62krcGh1r1ZGVKV+NjOhFXmQjXDJJdrNRJNXnEkWK9Opo3V81mXFkps1MW9l51HoGUR0SET/fyLM4HoOW7/q59dIGAAAAYENe2gAAAABs6Dcf3QEAgM/o5eWlnOj1kZFQ79s1rXqryVVn6jpemy0ftTFTx+1cpp1juWMdxyVTGdFYZuuqXpt9blG9o/NoZIxWtT1SvlXu7e2t1BZtrec/eqz1ueR7ySyhPJ2+x/wQaQMAAACwIS9tAAAGtBIR3/SS92YSEfdUEwtHyRuzbUb9zSRcHo0oqiRrzibqzWxxvSrpZbVvs2OVTTacTVKcqatnZovtaB5l+yEB8cdY+X3G9/Xd54eXNgAAAAAbktMGAGDALadNL/fCijwfM3lajuWi/vbupVVHRuteRvPu9PoR5aOp9qPXZnaMopw9o+OxIu/PiGN/e/lJMm1mc/ysuKfv/Gv9s8w822o5vqfMHPvKn22RNgAAAAAb8tIGAAAAYEOWRwEADLhcLqfX19ehLY+jZUmtUO/RZVfZJUIt0bns1tHVZVrZrcczy52yS5tmlnq1ZLYzXvFsR7YNrz7T7Fystpl53r3PRDSmfIyZz15mTvB9rfwO/4xE2gAAAABsSKQNAMCE7C9/2SiBVtLhauLdR+086lM2CudZiZYf9adXLjrWu5dqBEhUfiTaarT8igTO2bYzz3skqfKxjtXJakcTTzOnmhg8urb3HcDXN/Id/lWJtAEAAADYkJc2AAAAABvy0gYAYMDLy0tzuclIyPbtukfXRueu12upH622bnXcn7sdu/+7P/+ojtX3lx3TVt8y17buL+p3VP7+fNRWrx/R+ep9Rv1ttdmrN5oTrWO9PlWWvrTmZKvfWaOfV/Kqc+L435ZGcTrF353fgZc2AAAAABuSiBgAYMJIstNnbTe9ss3RhLrZxMEz2/1mtqLObjHc6ufMc2wlJ37U/165arLmFc92Jtn0iuTV2WcVyc5Bnm8meTRrvOe/Mx/pKyevFmkDAAAAsCEvbQAAAAA2ZHkUAMCE7DKc7JKU1jKB27Fsva1zrXLPWJIwsvRnNBls1FZviUw0ztGSnpF+t5L5Ztqs1jEjM5+e1Y+Z8cguQ4vafXt7S/WTWPb7b8USOPoy33Gf+blEy1m/2nI7kTYAAAAAG/LSBgBgwOVyGdpm+f58dG12m+xWXdU2W/2O/rLXZM1uO/1onI9GtoV+Rh33ojGtbltebas6fsdrMnMhs81zdq5H99Iq35uLtvxeqzWe0bGR74rPIPrues+2M8da1/aO7SJzn1+FlzYAAAAAG5LTBgBgwMvLy88ogaNqXpdsHpFsrpVIL1dOq1zl3Op8Ftl7bkUbZepfLXoeM9tpR3Vky7d+RT+Wm8mbFLXf+5yM1rGifp7ns+VJeYaPvvds3qnsZ+kjn181v9Yu/Z4l0gYAAABgQ17aAAAAAGzI8igAgAGXy+X0+vraPFfdBjmTIDhb16Ny2eU6mTpmE/n2jIxHKyluVG/Uj5klZ9Xz1SVyWZnnnUnc/KiuZy09Gl3C0FtOaKnU83z2JWrV5aHZOnaRvb/qvy+7LTfa+RnMEmkDAAAAsCGRNgAAi2QS5LbKVROyRm2PlKsmz60mSc5akbi4GmW0OvFu5tqZSKIVv4ZnI1FWPNNMWyP9iMa0GiH19vaWKk9sxbz9aNXP6vG61rXvee8zEWcrvn/fU/Ssdo0GGiXSBgAAAGBDXtoAAAAAbMjyKACAAS8vL6fr9TqUfHM0rLuaeLca0r566c/o0qaRJVnRkoSo3zPh8ytD76uJk7NLilYkpY7qnenHjGN9I/d5u+bHjx9L+/ZdZb8Ld1m6sjJx8i73vjoZdKa+j36ON9l736W/FSJtAAAAADYk0gYAYMJIstZMMshs+WwUyYqttqvlRn9pXp0Ut3rts7b7zSYFjurNJguNImFWROZk+vioH1F/qvcykyhVIuLne69Iq2ybvUjAUdVkv+8p+g6oRvjdn98lufSKhOm7E2kDAAAAsCEvbQAAAAA2ZHkUAMCAy+Vyen19bZ7rhcBnlozM1HEs0zMTsj+67Kq1nCUKu6+Gu2fbjI71lgmsWMaVXcoRLWuYWUaV6WNUbmacs8ssoiUc1aWDI8uoyMkmAc+cm7GyzWoS9Up9mXpXeM8k4Dt/pnbuW49IGwAAAIANibQBAJiwIkplJpFidTvrKMFrNnFntj+jiXezUUatfmaTgK7Yenf0uc0k/5ypN3NtNaJoZkvxbNRQSybiZ6Re6jJJbVdEuVWPVdvs9aNlRT9WGhm3TH0zY/QRduzTDJE2AAAAABsSaQMAsEg16iC75fGK/CSZ6IeR/CQZ2Tqe9St0dfyqUSQzzyrqb1Z1zmTbztz/TDTL6nGO+rHbNsWfXXbcR489q/5qFN9MTpvVfcuY6Ue2/cyz+mqRLh9NpA0AAADAhry0AQAAANiQ5VEAAItUk+tml2xUt4yuJsidSVh5PL8y7L5SRyapcq/eKLntTGLo6hbE1aUFKxK1ztxfVKa6BGpUr37LotbKjm11W/dqIvFH51fJfieuGI9nbZ29+nNgqeH7E2kDAAAAsCGRNgAAE0YiFEZ/qYx+fc1GxIxEV0TlM1svt+59ZhvaaqLZzDbjvXaiiKaZKKBqn1ZvUX48V41YqUYY3Ms+x+w4R221zt3a+PHjR6pe6qrfjyuecVU0h1dGoM2Um7nfldE12fqeFTX0XYm0AQAAANiQlzYAAAAAG7I8CgBgwMvLSyoJb0t1uU6m3IokwlkrQvWzSZtHQ/uzy7Ra16xY6nBff2a8solao3MrliSMjFtVNTFzpLrsr3Uvb29vpTb5T9X5Hc3X6jKp1aJ+VJMkR/dQHY+smWTk1TZGE6B/lF36MUKkDQAAAMCGRNoAACwyk+g2qqNaVzXyYuU23dlteav3kmn7UfnsL6wrntWt3LMTT/dk2prZ0nyX7epnEtdKRPwxMgnNd5Ht60z/M/XNRLs9e7xXjNHq6JcV/77sRKQNAAAAwIa8tAEAAADYkOVRAAATZhJKrgjNjhKxVpNo9kLao+S92fKjy3aeleyyupyql/A2U1826XA1UevOyVCj9rNL6lrXRmWy/ZSI+PmqyyY/egnLzDK96ndc67oViZmf9bk5Xtv7zK5YIrnSZ1kSdU+kDQAAAMCGRNoAAEzo/aJYjVLIRhqMbk1b7Uc1ee9IIuLq1ukro2micjO/ILdko2QqdT2qY3T78mwETzYaKRPpMxPdM5PM+zP+4v5V7ZqIOGvFvFrxvTDSxso6W5/pyLMTIvfKfZbvAJE2AAAAABsSaQMAsEg1Z0lWFMFTzUcT5WQZzUHzqP5qzoPsL82jUTK9/lbHtNV2Zryq0VS9eqt5IZ6VS6klutdsPqHRHCG9crb8/hijOZs+2o592k12jJ7xnJ+V92wHIm0AAAAANuSlDQAAAMCGvLQBAJhwPp9/9Xfver02l+4cl+W0ro3qeFTmmIw2Onb864muvfU/W1/vnlfK3mtm/EbafFT/o2OZcen1LZqLUf3V+Zydd9X7W/GZ6NU7WgfrVb8T+dye/dnLftd+Fl7aAAAAAGzISxsAgAGXyyUVWVL9tTgTCdKKcInarkRVRFEW2WNRv1sROdG9tPqUGe9WvVFdj5LiVtq8Lxcd641fVEd0T9nIoFYdrWcVPYNsv1eIIqVm2hTJsY+V36Hw1XhpAwAAALAhL20AAAAANvSbj+4AAMBn9PLykk7eezqdfrUE5f5YNjlidonAo6U+x3NRu63yt2OrlypkxiGTuPaR1jOoXPfoWPba0b61ZJ9L5p5bdbTmR7XenpX3fDvWG+9WuZu3t7dUP1gr+xmB706kDQAAAMCGRNoAAAy4XC6n19fXbmRC1Wi0QjVqoVXXSPljxMNIP6qRR6P3kK03ilRqyUZKtcZo5p4z/ZiJ+BndHrfXZiuBcVS+GvHTUp0zvK/oO0AUDt+dSBsAAACADXlpAwAAALAhy6MAAAbcEhGPLLs4XjMS/j+axLiapDibcLa6fGhlEtpHx6qO12aXhvX6e6w3m9y5uuyq1Y9WuWpy4l65qP7MsquZz1B1WZelNnta8Z0IX5VIGwAAAIANnSu/Rry9vf2f0+n038/rDvAO/vef/vSn//roTszwXQRfgu8iYBef+vvIdxF8Gc3votJLGwAAAADeh+VRAAAAABvy0gYAAABgQ17aAAAAAGzISxsAAACADXlpAwAAALAhL20AAAAANuSlDQAAAMCGvLQBAAAA2JCXNgAAAAAb+n8fmJJw4DAjRAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1440x360 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=\"assets/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": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAJMCAYAAAAYBLcfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+9UlEQVR4nO3dedB9R13n8c83BBM0EFZZBIFikGWGpeDn74fIQBQELEcFjILCAINLKbigomKpTNCpEaVKqpSgKI6AC/siUAqEIIb9l4uBkMgiIjggoOwEzFhAzx/n3OTmpm/3t5d7n/v0835VPZXnOfec7j7LPfme7u/pn4UQBAAAMJpTDroBAAAA20CQAwAAhkSQAwAAhkSQAwAAhkSQAwAAhkSQAwAAhkSQAwAAhkSQAwAAhkSQAwAAhnRqycovfelLw8c//vFttQXAAbnJTW6ihzzkIXbQ7dimCy64IJx++ukH3QwAnV1++eW6973vHb1/FQU5H//4x/W4xz2uT6sA7I1zzz33oJuwdaeffrqOHz9+0M0A0NnJkyc3fsZwFQAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGNKpJSvf7na3++TrXve6D2+rMQAOzC0PugHb9oUvfOGT559/PvcvYDwb718WQthlQwAAAHaC4SoAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADAkghwAADCk4iDHzL5iZu9c+XniNhq2Vud1zeyxFdudY2ZPWFt2HzN769qyU83sE2Z2sw3lnGVmryqtv5SZvdrM3mVml5rZH5jZNebldzGzt5rZu83slWZ2nZVt7jx/dun8+enz8rvPf3/AzH7XzGxefn0zO8/M/mH+7/US7TndzN5rZndaWfYLZvbM7R2FPDP7kJndcAvlPmo+Lv9gZo/asI77+GH/cP/aHu5fPty/diyEUPQj6bLSbVp/JN1K0iUV250j6Qlry06R9H8l3XJl2QMlvT5RzlmSXrWD/bzO/F+T9BJJD5v/vlDSfebfHyPpN+bfT5V0saS7zH/fQNI15t9PSrrHXNZfS/rOeflvS3ri/PsTJf1Wpk0PlPTGuZxvkPSPkq7XsI+ndjhOH5J0w87H/vqSPjj/93rz71fbz9Ljx89+/XD/2up+cv/ylcH9a4c/XYarzOxMM3ufmd1u/vt5Zvaj8++XmdnT5kj9fDO70bz8NnPk/w4ze6OZ3X5efmMze9n8RPAuM7unpKdIus385PXUeb1fMLMLzexiM3vySlt+xczeb2ZvknS79baGEL4q6YWSHray+GGSnmdmx+eniovM7C3L/Vnb16s8XZnZJWZ2q/n3R5jZybmdz1w+yXiFED4//3qqpK+RFOa/v0nSBfPv50n6vvn3+0u6OITwrnn7T4UQvmJmN9V0w3lbmK7m50p60LzN90p6zvz7c1aWb2rTqyV9TNIjJT1N0433TDN7/Xzszzezb5z3/9lmdvZyWzO7bP7vWfM5foWkv18t38x+fHlO578fbWZPn39/+Xx9XGpmP7beNjO7lZldsvL3E8zsnPn36PWV8ABJ54UQPh1C+Iym4/zAyHpFxw/7j/sX9695/7l/jagiYvyKpHeu/Dx0Xv4dkt6q6Qv36pX1g6SHz78/SdLT59/Pl3Tb+fcTmp9EJL1A0uPn368h6UytPQlp+nL8oabo/BRJr5J0b0l3l/RuSV8r6TqSPqC1J6F5+2OSLpp/P03Sv2qKgK+jOVKXdD9JL5l/P0vzk5DWnq4kXTK37w6SXinpmvPyZ0h65Pz7syQdcx7f10j6jKS/0JVPNW+R9KD595+T9IX598dL+tN5m7+T9Isr+/e6lTL/60r7P7uy3Fb/TrTpZpI+Iulv5r9fKelR8++PkfTy+fdnSzp7ZbvLVo7fFyXdOlL2jSR9YOXvv5Z0r/n368//vdZ8nG8w//0hSTeMXBdPkHRO5vr6Hkm/HmnHEyT96srfv7bh2ik+fvzsz4+4f50j7l/cvwqO32H/OVXl/j2EcNf1hSGE88zs+yWdK+kuKx99VdMXX5L+TNJLzewMSfeU9CKbhlql6csqSd+uKepWCOErkj5nVx83vP/8c9H89xmSbivp2pJeFkL4kiTNkffVhBAWZnbG/KRzB0lvDyF82sxuIek5ZnZbTTe3a+YOxor7arpJXTjv07U03XwUQvgRbyEhhAfYNC7955qOxXmavoi/a2a/JukVkv5jXv1USfeS9M2SviTpfDN7h6TPOesKZhYc6/2Lmb1e081Ykr5F0kPm3/9UUxdozskQwj9Fyv43M/ugmd1D0j9Iur2kN88f/7SZPXj+/RaazvGnchWlrq8Qwis0HcNm3uOHvcL9K477Vxr3r0OqJsiJMrNTNH3hvqRpTPAjG1YNmp5ePhu72Xirk/SbIYSrJJCZ2eMLyniepqe2O8y/S9JvaIr2Hzx34b4hst2XddWE7dNX2vScEMIveyqfu4LfMf/5ihDCk5afhRAuN7O/1NS1eF4I4b2abooys2+S9F3zqh+RdEEI4ZPzZ38l6W6absY3X6nu5pI+Ov/+CTO7aQjhY3O38L962qvpZv/VzDpXHJv5evialc++mNju+ZJ+QNJ7Nd3kg5mdpelp9FtCCF8yszfoymN9tfpmy89rrq+PanpiW7q54ue/9vhhj3H/4v4l7l9D6vkK+c9Keo+kH5L0J2a2fIo4RdJynPOHJL0pTGO3/zQ/Ockmy6en8yX9xLz8GmZ2pqQvaHrKWXqNpMfMEa/M7BvM7Os1jfs+yMyuZWbXlvTdifY+T9IjND1t/OW87Exd+WV69IbtPqTpiygzu5ukW6+0++y5Hcss9ltuqjyE8JUQwl3nnyfNT2Y3nbc9VdON4L3z38syT5H0q5L+YOU43MnMvnbe5j6S/j6E8DFJnzeze9j0KPDIlX18haRHzb8/arncpvH8525q75q36MqcgIdrSuxbHpu7z79/j/xPki/TdEP8QU03DGk6F5+ZbxC315SEuO4Tkr7ezG5gZqdJ+m/SFbkBm66vTV4j6f5mdr35yfv+87J10eOHQ4/7F/evD4n713hKx7d09THtp2hKkHuPpGvP6/yOpCfPv182/32JpNdLutG8/NaSXi3pXZqSuZ40L7+xpgP/7rn8b5mX/8VcxlPnv39mXufdmsbSbzMv/xVJ75f0pnmbq41LruzLOyU9f+Xvb5m3vUjS/5L0oXDlmOxyTPhakl4r6VJJ/2fe71vNnz10LvNiTU8595iXZ8e05/2+cN72Ekm/pyvH139mbtf75+NtK9s9Ym7LJZJ+e2X5sXnZP0p6+nIbTW8wnK+pa/V1unLc+GxJz0y079max6sl3XI+lxfPZX3jyj68bT6nv6Wrjmkn3+7Q1JX8wZW/T9M0vv0eSS/X9FRy1vzZhzS/nSDpp+d9vGBu4zmZ6ys6pj1/9hhNeRAfkPQ/VpZfcf42HT9+DsePuH9x/+L+daTuX8sLZ2vM7LIQwhlbrQTNbHpD4E9DCBcfdFuAfcH963Dg/oVNCHIAYAPuX8DhtvUgBwAA4CDwb1cBAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhEeQAAIAhnVqy8gUXXBBOP/10LRaLK5YdO3ZMkq6ybP0zz+cpy2296+e2bSmvR13e+mPHLCZVjvdcectPnceac1x6LHq1xVNerE0113lMy7XnPRaxujZte5Ob3EQPechDrLpRh8Dy/gVgLJdffrnufe97R+9fFkJwF/SMZzwjPO5xj9PqNma+++Jym9j6vcpL1VHTlpQe+51aP7dNrJ0t++1dL1VvrO3bbGfptjXXWWzbdS3XYE153mPqrUOSzj33XD32sY8dOsg5efJkOH78+EE3A0BnJ0+e1PHjx6P3L4arAADAkAhyAADAkIpyco4dO7axCzw3FJDqvvd29+eGQ7xDJDXDINtuZ6qcXHtLh7Vy6623L3fsYnWUDqnE6qg5Pi1DkJ71a85Fj3bm9sF7HGPl1X4fsDu54UoAcfTkAACAIRX15CwWC504cSL7lF6abOp90vQmZHqfeno8EdU8Yaf2tzTZt6RdqR6SlkTc0l61WJu89dX0hnjbUtq2lh6fluTqWB3L8rzHm96Aw8HbK875BOLoyQEAAEMiyAEAAEMqGq5a6j3XTa7s1Preoa6aIbZUGS2Jo97yvHWVDnn0GOLLdaOXzmPkHZ6sGbrzWm9fy5xJvYYPUmW3tI9k40lLwn/Lei3XR4+hTuCooCcHAAAMqaonx6ump6Ilodj7VFY622xM6Su7vdvU8lp5qk2r28Z6DFpeiU/V1as3pMfTbek16J0+oVdPZI9pCdadPHkyunwkyxcncrzXUMt6LQng9OAAfvTkAACAIRHkAACAITUPV/UYZqmZS6X3EIBHTcJsTI/hsm2W56nLO+9PrC2rn/VIZPa2r/RY1MzkXJoMXDOrs7ddsc+O8hwryxnbW4bHd7FsVWkCP4CroicHAAAMqdu/XbXNXpver3T20DJDcUuibmrZpvp6rLe+fst5jClNJC9pX+0r1N4ekF6vcnuT2VumaDgqvTYxscRj77VxUMti6MEB/OjJAQAAQyLIAQAAQyoKchaLhczsKj9LIYQrflLWt/eUt/7Z6uepMlbXq7Esz1tGrE2x/YmVF1sW25/c8Ytts2m/Nq3nPbapdsb2u8e2pfta0pb19WqOZ49rMHatrJadun68x2z9Z7FYuNt3WKWG21fFjmWPH2+9vdYDQE8OAAAYVFXisTdpzvu04Xmy9tTXUkYsiXT5uTcpOLVezSug23otOLZt6ZNh7tVX73rebVNqjkXp/nqTglPtSJW7SSrxuMd5PEq8Mx5vS49k+JL1ANCTAwAABkWQAwAAhlQ0XLXs7m3pJs8NB5XOl1LTxeud32S9jpqhgNJuae+wWsu2pUMpLd3jNfPf9JgbxjsU2DLsWTqM6R2a2lROav2W700I4Uj8A53eGY93YdN5qF0GII6eHAAAMKRuMx7H5J5MU8mcm7ZJaXnC8dTrfcIqLTenpkfD23vgfWV1U/mt26SSgrd1Hr1le7drSSCvaa83CbpHwv5ISmc83nUSt/fa2qc2A/uOnhwAADAkghwAADCkbjMer1qfmXXTzK0xqW1y67fMnltaXmwfYzPJ5o5Bqg5vO1t4ZsP1tjNW7qpYe1vqSMmdx9SMtLXnadM+etvZcj5btjVjxuNVPb5XXjXnf5ftQ1/but9hM3pyAADAkLolHntfcez9Gm1uvdLXfGOfp15p9iaE1tRful7pa9Y5qf3e5pNkab29X7lteZXbq/S18pZzwOvHZfYp8bhlPewP732bHrr+6MkBAABDIsgBAABDKhquivEOs3ilhiNahipivHNQlHbxlw6DlWgZuou1qWVoMVVualmsjty2u5g7J7Vdy/BSzLav1U3LUsNfR5F3eHNbwwg1Q60H3eajrua49xhGP4wOKuVhFT05AABgSM09Ob3tc1Rb8xTfOxm6Rao3oKUXYbmtt8egJqF4W9dFaU9k7ye2mvPd46l9fduj8G9XxdT0Pm6Lty371OajpPRc9Or13ZbeozClL8XsqteRnhwAADAkghwAADCkqhmPc1IzOcZmekzNOrvenb/+k5vZNiVXb2pmUW9dqTbH2hI7PrllsfJi9fZQOtuqd79j2+TOQWlbctdK6nrzllF6XebOT+x6LD2fuTrMjsaMx0u9ryuv3HfS06Zc29GX9/vpdVDnMXXt5f5f6y03VV7Mro4BPTkAAGBIzTMep5Izc8mpMd4kTa+WZNfSenOv76Y+a0nATe2jty2pRNma416ayNxyLlr226ulJ8z7qmjqGHiTGEvbVLPt6FqOxy62ra2jxzW8b3p+32uOzz5/d7Z5PZUmGR8kenIAAMCQCHIAAMCQqubJyXURprrna4YgYtuWDg15hwpSwzbe9XP1tsylEuOtIzX0UTqUVXN+SvfNW0Zsm9x+lx6flvMYO2al3b0118q2hu5GUvpd9J7r3st6tDk3PO4d9m7RMgSe4h36Ly2v1znbpR7nuyZFo0d520BPDgAAGFJRT85isdCJEyeqnrBjWnpeYuuV9grEyklFoTUJwN4ItkePT0tvVWr9Gt5j0TsJOfZZaSJ8adt7PQm18CY3b9r2qM54HOP9bm1rWU3PrWdZ7/a23Nt2UV5LL0JpeQfdeyP1Pd+1dbaWtw305AAAgCER5AAAgCE1z5OztA/zeJR2T5YmlXmHUbzDSyXtSy07iPlVahLtvENynvVXP+8xrOZdz5twXbMspnTf9qHb/DApvSZ7r9cyB1WP9nnvsy33430tL6U0+XwfHVQCcC5VILbeNtGTAwAAhlT1CnlM70itV29IrLz1clv0ei2zdB9bequ8EX6P10y9x91bR0sPSem5qknQ9j5tppZ5nx5repW87TrKevQsxOz6OPfopU2V1dKj0Ts51XsPaul1Tq3XK7m553HcZgJwywswsfW2+d2gJwcAAAyJIAcAAAypeZ6c0uGYVS3rlXZP9kiwq0n08y7rMRRXM39DbTdhrquxpRu3ZWioR3JzaZJxjnfYL1Zu6fBgrl5cVct3dp+Oa+mLE6X3wF5zQbVcz7Xta0kB6H1d9E7DWOV5Uaa3lnSNXSVG05MDAACGVPUKeU0iY++ETG+07H0qTymtK/Z5TS9LqtyclieM1GfeXqger81u86knZr2O3r1V20ycLy17ff1zzz1Xx48fr2vcYFqSKQ+Kp829r7/c96MlkbmlLaXlefXoCe61rae8Xfc69jjf22gzPTkAAGBIBDkAAGBI3f6BzlXeLrce85x4tQyxtSRKxYY+UnXGtq0ZumuZN6I2YdFbv3f9mqRpbxewp30tidzbTFJvWbapLUf9H+hsGc7eVz2HeL0J/70T6XuvV8r7Qk2voZeew0oHlRhf85JR76G7VfTkAACAIRUFOal/u6q3EMLVfszsip/Y58ufWDm5OlbL3lRHTqyMUjXbxtYv3e9Yeanjniu39JjlznesjvUycu3LXT+p4546t9525pZ5zkVu37z7s77tYrHIrndU1H53902P/Uh9F3dRXm49731mvbwasbp6Xysl39mSY+f9qW1vy/+3tv1doycHAAAMqerfrsqNr3lfJ2wZhysdn63JpfDWkSrD+5SRaqc3vyP3WawOz5hyTW5K7hrxluMpr2bc3HNdxM6j9/rIqT0X69uklsU+i9U7Qq9FTy35TbuUu1est7nlnPe+hrzl1dzL1+vInceW0YmWa6WlLZ56t3mMvddAj7a0oicHAAAMiSAHAAAMqWq4qqaLyrtt6jNv96x3m9JhjpohmtIu1t5yQx9Lu2xf7py1HDPvcF/vbuv1umqGR1uG7jyfbVovBF4hX+X9zhz0d3uVp809hldLlpUOvdQsS6lpc6mW8npfU57zfVD7v83rJ4eeHAAAMKSqIGf11a/YK4GxZbGf2Ct53leAd/G6WqqdMS1tb3lVM7b+ets9vVDr5aTKrTkXqfJWpa4fb3m562y9rvVevdakt9xxiZ2f2DJv21PHMbfMjFfIV/X4vu+Dnu3z3kdi12TNPa1HW0s/672e95jV1Oupo+ac7XK/Y+vVHLMcenIAAMCQCHIAAMCQqv7tqhY178H3SMiM1ZurP9atv/6Zt4xeybapbb3t8ybnxj7z1h/T41j0StT1Xiu1bfIui9UbWy93/ZYei/X1STwu15Ik2lvs2tnV8NCqXJ0HfcxakoK966X2reb/f16eerdZZ+3/Vzct64WeHAAAMKTmV8hjy1qi2xhved4eCm8vRyyRs7TOHr0SsTbVRMPenoJUeT16hrz74y2vZ49GiVQPUe/jWfJ5quzYZ/vQE7HvvL1x+2S9zb16Ebw9nb3rLVXaq9qrR7a0fbu0zTr39ftATw4AABgSQQ4AABhSUZBz7Nix6nf4S+dSib27H3uHPjW/y2o5sfVy7+SnPvPO25JqZ83xSbXTOzdLTmx+i57lxtqe25/YeqlzkGu7Zy6GXP2e68OThLf+k7umPfu66TrbVB/z5KSVft8PSuqekruee/Dej3eh9JzF/n/gXYb9RU8OAAAYUlXisVcuQSuVAJxL+MrVt6muGt4kNO+2pevVJKymPs8d71Kl+xjj3e9eiY2exOOaclsSn0uTllvsWw/EYbTvyci7bMs+7XdK7Jx57gU4vOjJAQAAQyLIAQAAQ2oervIOEZV2j7cMS9TMueKd+6Fku01llA59eIemvPsVa2vLUFLpEFFMTVdx6Vw03uHG0qHVluuyZegwV17JMCIzHpdjyO/w4ZwdPfTkAACAITX/21WlSbE1yaneBL/SHo2WJ2tvkmiPp3JvQnGON8HO+zq4p/7evVAtx7PlXLX0RMbK6JGg3NIW5PU6hzgY3hdGOI9joycHAAAMiSAHAAAMqXnG49hMq97ZIHvMGhmrNzcjZWm93nJTs2l6Z8D1btsyA25uPc/xaTkWufI85y73ufd4e49PrE7P/m+6Pr1t2sXMqmbMeOzFbLf9pe5VPcrNLcPY6MkBAABDan6FPDWDZGy9HjMPb6rD277S14dL68zVnyp7tf6WZORYGT1e9S4tqya5u+Uc9Hhl36s0+bymHanroiWRGdhn3pckSnnvxxgHPTkAAGBIBDkAAGBIVfPk5LrTt8U750rLbMnb0jJskxuW2FaXbqxNpXrM67Npm9K5l0r3o+a495gtueYayF1LKcx4jF3b5j2453xXOPzoyQEAAENqfoU89spuj1cCc9vG6vW+0pt6HTj1SrH3NedN9Za+UhxbL3ect/GKZM3r3aWvV8faHttvb1tz7UspffU0d9xT59t7/XqPRayO3HnhFXLsWu467VlHbhnGRk8OAAAYEkEOAAAY0lYSj1vmY0kt8yZ61swfUirVJu+cPKtKk4y988XkpOY58h670n3MtbN0/pmaRPhddll752ryfh+8x9HTpiUSj7FLpfd+qXz+rNJ7C8ZETw4AABhSUU/OpsTjdS09AC0Rdu9egVT5pb1K3h6fXI9Yqn29Xu0v3SbVTm8PUe5V6lhdNa+il+jRY9Kq9OkVOAy8Pf+elwR6lofx0JMDAACGRJADAACGVBTkLBaLrl3kqflIUvPV5Oadqam3dr3cHCmb5iUpmYfFO6dE7Fj0mEPHO8ePdx6d2DLvcczNGZRSM4fM+jEonWunZr9i5eTEjon32DJPDvZV7p4fWw9YRU8OAAAYUlHi8ZL3FeldJ5d5k3xL29AjEddbXu7Y1rxG76kjtW0uKbh3/bltPPX2Pt+lSeWlr8vXtNNbNsmWOAxqpofwLsPRRU8OAAAYEkEOAAAYUtVw1TZnk411NXpnjK2tq4W3m7RmqKK0fS3DQN6hxVy9pXWl9B7Kym27rTlpvEO13uNTehw9Q2LMeIyDkrqXtywDJHpyAADAoJqDnNirtalX/byv7MZ4XwH2vvK8qexN6+VeZ0y9slvz+nLLsYrJvUZc8hp26rOSdsbOT+pVau81Vfp6ea6uVPu8r7m2tK/XK7Kbji2vkOOgxb5j3ntK73slxkFPDgAAGBJBDgAAGNJW/oFOr9JEy1i9NUmnpXOe9EjOjfEmHnvnhSgpu5a3/l0kn6fUnKvSa6ClrpYkZ++2y/U854rEY+w7hqBQg54cAAAwpKKenMVioRMnTrifXHvM1lpSTqrslp6HXb6OHdvW++p8TY9P6vX81HZe3mPtPd8tr1fn6lg/FrGepNT6OS372Hu6AZ6KARwF9OQAAIAhEeQAAIAhVSUe57rdY8MsqfW9wyy76GLf1vDSLrQMm3jb2ZIA7C2vdAi013Cnd6grVW5p+7zDibnvQ0ubgRI9UhSAXaEnBwAADKkoyFksFhtnD47NDhvjXS8nNvtratZgr9IZNkvb2zqzbWq/c+cltr6nrtwxic3iG6s/VpdnRuFcHbFttjVTam691GzJObH1vNd5S5sBYFT05AAAgCER5AAAgCFtJfE4xZuI6i3Hm3SaW+ZJ5qyZibYlkbr02Hr327utN8nYOxuwp4xN2/ZIPi+9LmLXQulcO7HPNrWptO053mufYSzkbPM63ZZeL0fgcKMnBwAADKlqxuNVpU/xvaJr72u53m09vRzeV91j5W7zNeeW9pUqfYW8Vy9CS8+Ut65YUvT6tr17q2Ll5I5Zqn0tvW7AJjX3m4Oy3paa7wnGQU8OAAAYUnNOztI2xz9Lcx96RemlvVS9nxJKy+n9tO/VkovknYAv1ZMSa0uux6fHhGbeyfti5bYc44PuzcPR1tKb2Vvqevb+v2GfeqHQHz05AABgSAQ5AABgSN1mPF4Vm0XWO5trjHem3FT9ufatl7tadmwW21SbvFqOY8vMv6lyPcm5m477puOyXm+N3LH3tC+nZP9zxz22XuzaSrVj07WX2rZmPTPTYrFI7jewL7zfwZjU/w8wJnpyAADAkIoSj5e8SZ3eZNJcHTG9k996vgpfkxxbKpdMmqqjVwJsyWerWl4hj5Wz3ku1qY4e1573Wuh1fZYen/WeI896IQSdPHmyuG3AUs2UDT3q8HwWQ8/N0UFPDgAAGBJBDgAAGFJRkLOcJ8ebsLvKm6gbSxz1JAB7hrdKE8w8dbUkk7a007vf3m1L6+0ldYy9x8y7rTdBuCRh13MsSrfJre9dlvoerh8fEo/Rg/c726sO7zKSi48uenIAAMCQtpJ47OVNEk3VlUu09NabWi9Wfk2Sba78nkqTv1P7VpOw22Pfaq6LlifH9WPWI5F9tZzYtepNUs/V4Vl/U1tIPEaN2PXXu+fGW0fvew/GQU8OAAAYEkEOAAAYUtFw1WKx0IkTJ6LDAzVz4njnplmuVzoctaktqfVS29TsY485cbbJe/5ay99UR0rNMFSP+Wlix6RkzpmSttUMy8bWS31H6L7HrvQaOi4tz/v/JBw99OQAAIAh7SzxOBZdl/Ye1PSkxD6rfVL39kx59ehd2qS05yy1bzVPZ95jHFOakN7SltJeE2+ie81TpLdnqPQ6B3bloJKCU/XSo3O00ZMDAACGRJADAACGVBXkeGeQzM04mZqltXR22tx6sdlwveWl2pTbx5Ztvfud2h/vstJzmtvHWJs85W8q18s782+q3NTxzLUzNltxqozcdZHa1ntuPW1mxmOUyl27qWWt9ZXU6/1OYEz05AAAgCFVJR7H1CTlliaYepOHeyW79mhTSm6/vUmnLbPnel/jL/mstU3e9VOvTfd+lbolibHlOu9Rx6pN9THjMXrYdeJx6ntEMjIkenIAAMCgCHIAAMCQqoKcWFJlS3JZTaKlN6nM26bYfqzXEdvHXP2pBNTcfsfqTR372La5Zd6k5doE15r9KW177ifVltS+eY9TbL1c2z3X26ZzG9v/1DLPtUfiMXrwfu9qpO4V3v//kIx89NCTAwAAhlSUeHzs2DF3r4hUnjjr+dxTR0syp3dbb72l9ecSa72Jdqn1WtqXKq9XIrW3vFgZvfc3VW7LddYjmb3mWJBwiW3q3TPiLa/3ehgHPTkAAGBIBDkAAGBIzf9AZ4vSoY/Y57m2lA7bpLatGRZJ6T0cUzoctFpezXBabRktw2qxtnivgVgbUsesdztj9Xu39Q495b43m+pgnhy08H7HSr8/uaFbz7Je32McTvTkAACAITUHOaWvjXtfpc69ZuxtS+pVXe8ruLFyva8qp+T2O1VvbB9zP6k2pNaLtTPXJu9+5z73nB8v76vZqW1j7YytlztmLfWuf5b7LuWuBV4hxza1fO9y5XmWxbZtuY/gcKEnBwAADIkgBwAADKko8XixWOjEiRPZBNdUkuQ253JJJaH1SNysKSNVbsuyWBt6zNuyWk5LV6436XAX3cXeuWZqkxK915s3qbwmeTjXrlRbSDzGNvT+bte82LCNduBwoScHAAAMaauvkLe8ppfbNvZ0Xvoatvc19RhvL8u2XqPMtbPkKd6z3vr63t6l3LZePV51z5WXegL09ga1qL0Ga9rC0y1GxbWNVfTkAACAIRHkAACAIVX9A501yZq9h3JSwyYtyba5OlJSQ2ib6k3Vn5IrY5dDWKVJxt7ht97dzt6EcM9nq5/XJFd797smsT1Vf+/kZgDYZ/TkAACAIVW9Qp4T6wEpTdzcxavPuR4pT++K9/Xh2DbeHrGap+7SJOzS3rTS1+831dW7vNJejlR5NdeW99pvSVxf/yxWf816ADAaenIAAMCQCHIAAMCQioKcZeLxqtQ/eBZbFvvH2nLrpcTKW10WK8P7D7PF2lfazlhdnv3a1I7Yti3LSqWOXaydsfMT2yZ2HHP7vb5f60NEqfK810CKt+2pa6BmH1Nt91776+vwD3QCGBE9OQAAYEhVMx63vPIdU5NA2SNxsjRBuWUG3Jb27uJ1X++r1N5X7VuObapN3mWbyklJ7XdpWTXtTLXJmxBfk7AviX+7CsCQ6MkBAABDIsgBAABDKgpyFouFK4k3lvyYSoiMJWnm1otpSTJOJYLG2tQjodibXF2TKJtKmvbuR+p4po5Tr/PjTWb3JrjnjndtQnisjJpk8ZY2pfYxt4zE44NTc75Sy4CYo3yd0JMDAACGVPVvV3nl1i2NLHOJm6VJuS3lxRJmeyRX1xwzb3JqKik3lVTdkmybWy9mW8cxd2xjr7d71CTie+ovaUPp+iQeHyxvUnvpsm29mIDDI3UP2sVLLPuGnhwAADAkghwAADCk5sTjpVyipzfxOMW7XqqulvJi25YmWZfokfhbk6ibqr92HzYlOfdI4PZum7v2ahPHvWoSyNfr2nRNlSYyk7S6f7zXRu01hKOl5X4zEnpyAADAkKoSj3PJSz2SL1uSaHPleRNwU2V4eetPiW2bS0hM1VuaIJxKSl7Ve72clgTllp6i9fJqzkWu7HW5tpdeo0f96e4gtJz/mvVwtHmvt9HvBfTkAACAIRHkAACAIRUNVy0WC504caKq2z81RBTr2m8ZBssNh3jnEUiVm/rcO1TSa3inxzKv0iGi3LFoOc+l56p0fh7vdVlTb+kwqnf4a/Su51GUDie3DDvjaOp97z+s6MkBAABDKurJWap5qvQ+zXojzlTPh/eJveXpuEdytbdN3m1z5S2V9jR5e2NqjknPnqRcHaXnsbRXJldXrpzeiaelPYbYHe89zXudcF6xVHrvHx09OQAAYEhFQU7s367a5cRtq0onkfNOBJeabC9XRu999B7bWPu828b2MVV/av3csfB8lvs8ta+5a8XbZu/5jpWbqmvT/qb2u3TbknO/ug3/CvnBK712cstwtNXcW0ZETw4AABgSQQ4AABhSVeJxTSJjS3JqS3dbab2l7atJFG5JCI2t16OO0rbk1ve+au6tw1tvj/Pc+7X2loTnloRS7/WI/dCSUI6+dvFCxL7w3o8PK3pyAADAkLbyCrm398S7rfeJPbVtrP29o9Zd9Gqt17WpvtLX1FPl1Zzb2Lax8lteue5xHL3leV+nb2mL51ysqllvU70nT570NhMdlH4XR3/a3gct5+KwnB9vT/g+70MpenIAAMCQCHIAAMCQioKcxWLhmvcjNb9KzbZLsflIYp97693Uhk0/3vW9Sudo2VRv6Xw6LXqUsSp1bL1zvniPfWy90uNYc11458vxXvve/V7/rqTWY56c3Yqd69SymvvLYeGd72rb9W9z2b6IHeODOu67Qk8OAAAYUlHi8XLG41Qyr1T+6nHs81x5pQlSsXq9r/t6PtvUppbkLu9+9zg+PdQkQ+/iVc1U0nAqibDmWm2p31NGTVtGfUIbzWFLYt2mg06A7f3dPujzWPrixD61vRU9OQAAYEgEOQAAYEhFw1WLxUInTpyIflYzV0jLsFaMd/ipdE4Rb9dladtivENYvebJiWkZCiz5LNemmmsgVW/pcJl3yLS3lu5h71Apdi93Hg7LuSkd5veWsU+8+1j6/6Z9HPrZ93PRip4cAAAwpKoZj1e1JML26N3wJkjVzFIZez14vbxeT/ipNtXsY2r90qd97znz9sa09FbV9OrVXo+9zm2PROaa81jSE8eMx7vV4xo/aDXfbe+2uzwGLeeix317l3q9FHOY0JMDAACGRJADAACGVBTkLOfJiYnNWrn6sz6r6+qyVS0zxsba4pVa37s/qTbF9itWR27b2HqpbWNtjx2/Upv2o5Z3v2Pb5MpLnb9N23jKTS2LbVNzXabEysud71ibzZjxeNdi58a77CDkrrVtldf7O5NqS25Zj/J2sT+lWu4jhwE9OQAAYEhVicc1CZmpz3uVtyynd3JzyTqb6vQmoXmTTlv2y5tM59WjLd52ehPncomAPRJ/W9oZ27Zlf1qSHUk83g/bSjz23nt731NTvOUdVAJs6p7hfQkgV95B6/Vixb6jJwcAAAyJIAcAAAxpKzMex3i77L3DAqXDBzEtXaG9h2hS3Zk1XYjeej3LvG1qGRrzDvnkhmhS2+a6lD3HuWa/U0NJLd323m1bhk+xHd4hVM9nLVrSDHLr97hHlw4d97KLIcPU/xP3yb63L4eeHAAAMKRuicdeuUi/R4TfkuDpLS/VppreKm9U35LU1pKcuqmsXNu823rPRWvZnm1Le6687a3p8UmpqXfExMLDIve9K00oL11W2oPbuiymRx29bSORv2d5u7SPbWpFTw4AABgSQQ4AABhSt3+gc5W3K7Z312XLsIFnCKmmu7R0aKx3l2zNcSwdhund3VzaZdqrK329vF0MB9UM0/X4PjBstTvec7PtZTXblSYP99q2Rx1eLW0pTWXoVR786MkBAABDau7JKU2ibXntNVZOSyJzS2KYt66WSN9bRo9ekx5P/TWvZfd4mumVqFn7CnlNm9Y/6yVXBz04u9fyYsC21ttWj4G3B7Vm6oRdvobdu8eTHtSDQ08OAAAYEkEOAAAYUlWQY3b1f4Z99Se2XmpZTghhY3fk8rPc57H2e8XWj+23t67Yeqv7sf7jrTd3vEuPxfp2q3XUtD3Wztg2nmurRs21UvJZrss8t6+x9bzXT0ubzUyLxSJbHrbLe42nvjuxazH3vSyVKsP7Xd3mej3uFb3vN6XnFv3QkwMAAIZUlHh87NixaE+ARypBzJsY502q9DxR12h5AvImSPdKEt1WQrh3nWV5PRKvc+XUXAOldfVOcPQmxHuVJlKvX48nT54srhNlar4LpYn0u5T73pUmFOf2IVVHy32z9Bj3+n7u87ldta3k7l2hJwcAAAyJIAcAAAxpq/PkxNTM2VDS7b5pWa+2eLb1Dqvlyu3R9pbyWurq2Y1cw1tXyzCdt/u6R7d9y3Bebj0z07nnnqvjx4+7ysDuxM77Piemetvbsg+lddTcR3Zx3Hscq95DSD3+37Rv6MkBAABDqurJaUnQ6hUBphItW5LfVnkTdT3rt/Qs1Czzqt2Plp4kb69ETS9dS2/e+ra9eqFaeiJj5aXKrb0GSDzeH7lrx3vv26XSe17Nd6v0hY3SY7HN71hs29LRgNx62z73h6n3ZhU9OQAAYEgEOQAAYEhVw1XeYR4vb5dg7/kRcm2pTfxtmc+nV4Jpav0ewyE1x7p0+M9bXqycmuvC0z7vEEHNddRyTZd2qx/WJELsd+KxV6/rr9e9pKT83mV757Za1fMaqKnzMN0z6MkBAABD6vYKee8eiFwPUare3HotPRmp7bxP4qVJtLk29Uj6qzm2nvVjWnruYmqSDns8ibS8Qh7rcYqtX3psa14l5xXy/XVYZsWN2cc27aOaaUu21VvV+8WWfUBPDgAAGBJBDgAAGFLzPDkxu5jDJrVe72S03gmcu0wcbJkJtyXJuKWbvceQU+/k79TwUs2QpXcfe881ElufeXIOl5brCYfXLs7tthO5DwI9OQAAYEhFQc5isZCZKYRwxU+MmRU9WebKW/18+ZOqY/mZ5/PUT2p/vPu43t7145fan1ibPMd9db3cPnrK857v2P5497HlmHn3J1fH+me5c+Ftew+x/e9Vr5lpsVh0aSe2z3sPAEBPDgAAGBRBDgAAGFJR4vGxY8fcQwFSv7kdWuYIiSXklSYSx+aS6cF7LGrmYYmVm0pYzG27vqwm+bxlvaXcNdWS+Lu+rCWhOKf3fqfaHKt3/dySeLzfWpLMgaOMnhwAADCkop6cxWKhEydORHtFWp4uanpqSntXcq8Ue7apebJPlVszq2TpftT0FsWSblNKX1PP9bzk2lfaltpepV7TDnh7fFLt7H3tbSob+8/bUw6AnhwAADAoghwAADCkqsTjHkM/m5Z5yyhNWs7NbOtJ3KzZB++MujHeYa3ewyupdqa6xXPHM1ZXatixdBhsU1tS9ca2Kd3HWB01Sfc9hiJ7JbNjf/W4lwJHBT05AABgSFby9H/++ef/m6QPb685AA7ILe973/ve6KAbsU3cv4Bhbbx/FQU5AAAAhwXDVQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEgEOQAAYEjFQY6ZfcXM3rny88RtNGytzuua2WMrtjvHzJ6wtuw+ZvbWtWWnmtknzOxmG8o5y8xeVVq/o32nmdkLzOwDZvZ2M7tVYl0zszeZ2XeuLPt+M3t173aVMLM3mNmxLZT7QDN733xsotdYyfEDJO5fPZnZk83s0vn796Pzsgeb2cVm9l4z+6N52e3N7C1m9m4z+1szu2GizO8ws7eamc1/X8PMLjKze/Zuv9dB3/8990JsVtOT8+8hhLuu/Dyle6uu7rqSim8SG7xR0s3N7JYry+4n6dIQwr90qsPrhyV9JoTwnyQ9TdJvbVoxhBAk/bik3zGz083sDEn/W9LjaiufA6e9680zs2tIOlfSd0q6o6QfNLM7RlZ1Hz9gxv2rn7dJ+i+STkj6TTM7VdKnJH2rpu/tHc3sXvO6jwgh3EnSWzTdx6JCCOdJ+rCm77Yk/ZSkRQjhLbWNnNu1j7L3r4J7ITbo8j84MztzjjRvN//9vJXI/jIze9oc8Z9vZjeal9/GzF5tZu8wszea2e3n5Tc2s5eZ2bvmn3tKeoqk28xPXk+d1/sFM7twfmp48kpbfsXM3m9mb5J0u/W2hhC+KumFkh62svhhkp5nZsfnp4iL5iePq22//nRlZpcsI3Aze4SZnZzb+cz5Ak35XknPmX9/saT7Lp9gYkIIl0h6paRfkvQkSc+V9Bkze/l8HN5mZndOtXP+eZ+ZPVfSJZJusbLOA83sRSt/X/EEY2a/b2aL+TxecbzXjs1lK7+fbWbPnn+/kZm9ZD5fF5rZt2aOy3FJHwghfDCE8B+Snj8fq3VFxw+I4f5Vd/8KIfz1/PB1iqSvTovCBSGEL0gKkk6XdHkI4b0hhA/Om50m6fJUuZJ+VtIvm9l/lvSTkn7JzH7Qpp6gS8zsimAgcc95tpn9gZm9XdJvrx2Dt81lL/9+g5kd29P7v/deiE1CCEU/kr4i6Z0rPw+dl3+HpLdq+sK9emX9IOnh8+9PkvT0+ffzJd12/v2EpNfPv79A0uPn368h6UxJt5J0yUqZ95f0h5JM0xfsVZLuLenukt4t6WslXUfSByQ9IbIPxyRdNP9+mqR/lXT9eZtT5+X3k/SS+fezJL1q/v2c1TI1BQq3knQHTQHINeflz5D0yPn3Z0k6FmnHJZJuvvL3P0q6Yeb4f52k9837eZqk35P0P+fPvl3SOzPtvJWmG9I9ImWfKumfJX3d/Pfva3oCk6Trr5yTN0i68/z3G5b7JumylbLOlvTs+fe/kHSv+fdvlPSelfPwrEg7zl5dLum/L6+b1uPHz9H+EfevTfeF4vvX/Nk152PxuLXlvyHpRWvLHiDpPZKu6zhPPyXp85IeLelmmu5LN9J0j3q9pAfN62265zx7Pq7XiJT9s5KePP9+U0nvm3/fu/u/nPdCfjb/1HTj/XsI4a7rC0MI55nZ92vqWrvLykdf1fTFl6Q/k/RSm4Za7inpRSuB62nzf79d0iPnMr8i6XNmdr216u4//1w0/32GpNtKurakl4UQviRJZvaK2A6EEBZmdsYcqd9B0ttDCJ82s1tIeo6Z3VbTze2auYOx4r6ablIXzvt0LU03H4UQfqSgnKQQwhfN7AWavtz/z6bu4O+bP3u9md3AzK6TKebDIYS3Rcr+sk05Pt9tZi+W9F2SfnH++AfM7Mc03WRuqqnr9GJns++nqet6+fd1zOyMEMJCUrdjAzhw/4qrvX/9hKb7ybnLBWZ2F0kP1hSMLZedIumPJX1bCOGzjvacK+kpIYRnm9n3SnpDCOHf5rL+XFNQ+PJMGS+az8G6F0p6raT/KekHNPWiSFNAutf3f5TrNlY5X8R3kPQlSdeT9JENqy67Nz8bu9l4q5P0myGEZ6614fEFZTxP01PbHebfpenp429CCA+euyDfENnuy7rqMN/pK216Tgjhlwva8FFNw0UfsWnc+ExNY9o5X51/Uja1U5K+mNju+Zq6iD+taSz8C2Z2a0lPkPTNIYTPzF3Cp0e2DRvqO0VTz1Gum3ppeVyWbj4v27Re6fEDroL7V9X9S5LuLOmv15bdSdLfrn3fbybpcyGEf/AUGkL4qpmF/Job7znShvtcCOGjZvYpm4b2H6orc4T28f7vvRdig55Jpz+rqSvyhyT9iZkto+BTNHW5af7sTSGEz0v6p/nJaZkAu3x6Ol/T08Eys/5MSV/Q9JSz9BpJj5mfqGRm32BmXy/pAkkPMrNrmdm1JX13or3Pk/QITU9efzkvO1NXXkCP3rDdhyTdba73bpJuvdLus+d2yMyub1dNDox5haRHzb+franLO8z7c35m26U3Snr4XOdZkj45H99N7cz523m7H9UU8EhTN+4XNT2V3lhTElzMJ8zsDvP/MB68svy1mrqfNbfnrpk2XCjptmZ2azP7Gk0389hTbfT4ZcoGYrh/ld+/JOmPNA3zrXqzpl6bVZ+R9POrC8zsN83swco7Kek+ZnbDOc/lBzXdp6TN95ycF2jqpT4zhLDskd6L+//aOt57ITaoCXKuZVd9BfMpc7fpj0j6+RDCGzV9WX91Xv+Lko6b2SWavpC/Pi9/uKQfNrN3SbpUVyZT/YykbzOzd0t6h6Q7hhA+JenNNiV5PTWE8FpNeR5vndd7saRrhxD+TtPF+y5NTxcXbtqJEMJ75ra9PoSwjPh/W9NbAhdpcy/XSyRd38wu1dTj8f65vL+f9/m1ZnaxpPM0DevIzJ5l8des/1jSDczsA5J+TtLy9cCbanpi8DhH0t3nOp+iK7800XbmzN27r9IUyLxqXvYuTV3r79V03N+8YfMnztu8RdLHVpb/tKRjNiVZ/r3mJyebkv2eFWnDl+c2v0bT/3heGEK4dN7m183se+ZVNx0/YBPuX/3uX9I0pH37tWV30jQct+pMXX1o+k6SPr5pH5dCCB/T9N3+G03H5h0hhGVgt+mek/NiTQHDC1eW7cX938xuZmZ/NZe78V4IH9v2g6+ZXRZCOGOrlQzGzH5S0j+HEIjYgQPE/Wt7zOw1IYQHHHQ7MDaCHADYgPsXcLhtPcgBAAA4CHs32y0AAEAPBDkAAGBIBDkAAGBIBDkAAGBIBDkAAGBIBDkAAGBI/x9OqXO+iyUIrwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x720 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": "code",
   "execution_count": 6,
   "id": "0ad5b92f-521f-48c2-9f5b-7950f5a3f931",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjkAAAJMCAYAAAAYBLcfAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAABC8ElEQVR4nO3debhtV1nn+9+bBBKkCX0nSCiKovEiPHDqHLC4kAIEvZbSCIjCBR9snpKUCopVeFEEtQqUuvJYRRAsvYKWItI3twRCIAakOdkYEkKPCmWQRmmUgLEuybh/zLGSmZWxxnxHs9bee5zv53n2c9aea84xxmzW3O8a453jWAhBAAAAozlpvxsAAACwDQQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSKeUrPya17wmfO5zn9tWWwDsk1vf+tZ61KMeZfvdjm06//zzw2mnnbbfzQDQ2eWXX64HPOAByftXUZDzuc99TmeddVafVgE4MM4+++z9bsLWnXbaaTp69Oh+NwNAZ8ePH9/4HsNVAABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSAQ5AABgSKeUrHyXu9zl7972trd9eluNAbBv7rDfDdi2r371q3937rnncv8CxrPx/mUhhF02BAAAYCcYrgIAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMiyAEAAEMqDnLM7Aoz+8Ds5xnbaNhanTc2s6dUbPdsM3v62rIHmtl71padYmafN7PbbijnTDN7U2n9pczszWZ2kZl9yMxebGYnx+X3NLP3mNkHzeyNZnaj2TbfFt/7UHz/tLj8PvH3T5rZfzEzi8tvambnmNkn4r83ybTnNDP7qJndY7bsZ83sJds7CsvM7FNmdvMt1/Fv4/H7gJm9y8zuHpdf18x+N753kZmdOdvmumb2W2b28Xjcvi8uP9XMXhHPxfvM7IzZNj8Xl3/MzB62oS13jNt9MpZz3W3u+8i4f20P9y8f7l87vn+FEIp+JF1Wuk3rj6QzJF1Ssd2zJT19bdlJkv5a0h1my75T0tsz5Zwp6U072M8bxX9N0qslPS7+foGkB8bXT5b0y/H1KZIulnTP+PvNJJ0cXx+XdN9Y1p9I+q64/NckPSO+foakX11o03dKemcs55sl/YWkmzTs4ykdjtOnJN18F+civv5eSW+Or8+S9Lvx9S0lvV/SSfH350j6ldl1dvP4+imSXhxfP07SK+Lru0u6SNKpku4Yj+3Jibb88exaeLGkH9/2tTjqD/evre4n9y9fGdy/tnwtzn+6DFeZ2ekxkrtL/P3lZvaj8fVlZvaCGKmfa2a3iMvvFCP/95vZO83srnH5rczstTHKvMjMvl3S8yTdKUalz4/r/ayZXWBmF5vZc2ZteWaMRN8l6S7rbQ0hXLk66LPFj5P0cjM7Gr9VXGhm717tz9q+XuPblZldsopszewJZnY8tvMlq28yXiGEf4gvT5F0XUkh/v4vJJ0fX58j6fvi64dKujiEcFHc/oshhCvM7DaaLvL3hunK+j1Jj4jbPFzSy+Lrl82Wb2rTmyV9VtITJb1A0433dDN7ezz255rZt8T9f6mZPXq1rZldFv89M57jN0j68Lz8+I3j+bPff8jMXhhfvy5eHx8ysx9bb5uZnWFml8x+f7qZPTu+Tl5fXrNzIUnX19Xn4u6S3h7X+YKkr0g6Et97sqTnxveuDCH8XVw+P+avkvTg+M304ZL+KITwTyGEv5L0SUlH1/bRJD0obic5zhnKcP/i/hX3n/vXiPevigjxCkkfmP18f1z+HZLeo+kD9+bZ+kHS4+PrZ0l6YXx9rqQ7x9fHFL+JSHqFpKfG1ydLOl1r34Q0fTh+S1N0fpKkN0l6gKT7SPqgpG+SdCNNB/3piX04IunC+PpUSV+QdNO4zSlx+UMkvTq+PlPxm5DWvl1JuiS2726S3ijpOnH5iyQ9Mb7+bUlHnMf3LZK+LOkPdfW3mndLekR8/dOSvhpfP1XS78dt/lzSv5/t39tmZf7vs/Z/Zbbc5r9n2nRbSZdKekf8/Y2SnhRfP1nS6+Lrl0p69Gy7y2bH72uS7pgo+xaSPjn7/U8k3T++vmn893rxON8s/v4pSTdPXBdPl/TshevreyX9kvNcnKXp28lfz8r6MUmv1HQjv6Omm8T3SbpxXO/X47l4paRbza6R283K/YvY/hdKesJs+e/Mj19cdvO143N7VfQK8MP9K75+trh/cf86ge5fp6jcP4YQ7rW+MIRwjpk9RtLZku45e+tKTR98Sfrvkl5jZjeQ9O2SXjkFepKmD6s0RX1PjGVeIenv7drjrg+NPxfG328g6c6SbijptSGEr0tSjLyvJYSwZ2Y3iN907ibpfSGEL5nZ7SW9zMzurOnmdp2lgzHzYE03qQviPl1P081HIYQf8RYSQniYTePSf6DpWJyj6YP4X8zsFyS9QdL/iqufIun+kv6lpK9LOtfM3i/p7511BTMLjvX+xszerulmLEn3k/So+Pr3NXUhLzkepmh/vey/NbO/NLP7SvqEpLtK+rP49k+a2SPj69trOsdfXKood32FEN6g6RguCiGcLelsM/tBST8v6UmS/h9N18yepE9ruoFfoelc3E7Su0MIP21mPy3pP0v6Pz11YWe4f6Vx/8rj/nVI1QQ5SWZ2kqaD93VJN9EUOacETd9evpK62Xirk/TcEMI1EsjM7KkFZbxc07e2u8XXkvTLmqL9R8Yu3PMS231D10zYPm3WppeFEH7OU3nsCn5//PUNIYRnrd4LIVxuZq/X1B14Tgjho5puijKzfyHpu+Oql0o6P8RuRTP7H5LurelmfLtZdbeT9Jn4+vNmdpsQwmdjt/AXPO3VdLO/cmGdq45NvB7mCWZfy2z3R5IeK+mjmm7ywaaEuIdIul8I4etmdp6uPtbXqi9avV98fZnZf1Q8ront/kjSb8b3viHpabPt3i3p45puXl+X9Jr41isl/XB8/RlNN7lLzewUTd/uvzhbvjI/TytflHRjMzsl1p1aB424f3H/EvevIe9fPR8hf5qkj0j6QUm/a2arbxEnSVqNc/6gpHeFabzwr+I3J9lk9e3pXEk/HpefbGanS/qqpm85K2+R9OQY8crMvtnMbqlp3PcRZnY9M7uhpO/JtPflkp6g6dvG6+Oy03X1CfihDdt9StMHUWZ2b01dfqt2Pzq2Y/UUwB02VR5CuCKEcK/486z4zew2cdtTNF2wH42/r8o8SVNE/uLZcbiHmX1T3OaBkj4cQvispH8ws/vGMdEnzvbxDZoiesV/Xx/LPmpmv7epvWveratzAh6vKbFvdWzuE19/r/zfJF+r6Yb4A5o+kNJ0Lr4cbxB31ZSEuO7zkm5pZjczs1Ml/RvpqvHoTddXUgjhmavzEbe58+zt79b0LU3xWF8/vv4OSd8IIXw4TH2xb9TUtS1N34xX4/fzY/5oTV3PIS5/nE1PL9xR0ze942vtCpLeoas/Q1edM3TF/Yv716fE/evMuM0496/S8S1de0z7eZoS5D4i6YZxnV+X9Jz4+rL4+yWaEp5uEZffUdKbNWVnf1jSs+LyW2k6CB+M5d8vLv/DWMbz4+8/Fdf5oKax9DvF5c/UFJm+K25zrTHt2b58QFPi1Or3+8VtL5T0K5I+Fa4ek12NCV9P0lslfUhT199HJJ0R3/v+WObFmr7l3DcuXxzTjvt9Qdz2Ekn/VVePr/9UbNfH4/G22XZPiG25RNKvzZYficv+QtPYqcXlN9N0Q/uEpLfp6nHjR0t6SaZ9L1Ucb5V0h3guL45lfctsH94bz+mv6ppj2tmnOzR1Jf/l7PdTNY1vf0TS6zR9Kz0zzMa04+ufjPt4fmzjakx70/XlGtOW9BvxuH5A04f0W+PyMyR9LLbrbbrmUy53iO1YPy6nafpm9ElNN4F/NtvmmbH9H1N8giQu/x+Sbhtf/7O43SdjOaeWfm754f4VX3P/4v51Qt2/VhfO1pjZZSGEG2y1EjSz6QmB3w8hXLzfbQEOCu5fhwP3L2xCkAMAG3D/Ag63rQc5AAAA+4H/uwoAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAyJIAcAAAzplJKVzz///HDaaadpb2/vqmVHjhyRpGssW3/P837Oalvv+kvbtpTXoy5v/aljlpIrx3uuvOXnzmPNOS49Fr3a4ikv1aaa6zyl5drzHotUXZu2vfWtb61HPepRVt2oQ2B1/wIwlssvv1wPeMADkvcvCyG4C3rRi14UzjrrLM23MfPdF1fbpNbvVV6ujpq25PTY79z6S9uk2tmy3971cvWm2r7NdpZuW3OdpbZd13IN1pTnPabeOiTp7LPP1lOe8pShg5zjx4+Ho0eP7nczAHR2/PhxHT16NHn/YrgKAAAMiSAHAAAMqSgn58iRIxu7wJeGAnLd997u/qXhEO8QSc0wyLbbmStnqb2lw1pL6623b+nYpeooHVJJ1VFzfFqGID3r15yLHu1c2gfvcUyVV/t5wO4sDVcCSKMnBwAADKmoJ2dvb0/Hjh1b/JZemmzq/abpTcj0fuvp8Y2o5ht2bn9Lk31L2pXrIWlJxC3tVUu1yVtfTW+Ity2lbWvp8WlJrk7VsSrPe7zpDTgcvL3inE8gjZ4cAAAwJIIcAAAwpKLhqpXec90slZ1b3zvUVTPEliujJXHUW563rtIhjx5DfEvd6KXzGHmHJ2uG7rzW29cyZ1Kv4YNc2S3tI9l40pLw37Jey/XRY6gTOFHQkwMAAIZU1ZPjVdNT0ZJQ7P1WVjrbbErpI7u929TyWHmuTfNtUz0GLY/E5+rq1RvS49tt6TXonT6hV09kj2kJ1h0/fjy5fCSrByeWeK+hlvVaEsDpwQH86MkBAABDIsgBAABDah6u6jHMUjOXSu8hAI+ahNmUHsNl2yzPU5d33p9UW+bv9Uhk9rav9FjUzORcmgxcM6uzt12p907kOVZWM7a3DI/vYtlcaQI/gGuiJwcAAAyp2/9dtc1em96PdPbQMkNxS6Jubtmm+nqst75+y3lMKU0kL2lf7SPU3h6QXo9ye5PZW6ZoOFF6bVJSicfea2O/lqXQgwP40ZMDAACGRJADAACGVBTk7O3tycyu8bMSQrjqJ2d9e0956+/N38+VMV+vxqo8bxmpNqX2J1Veallqf5aOX2qbTfu1aT3vsc21M7XfPbYt3deStqyvV3M8e1yDqWtlXnbu+vEes/Wfvb09d/sOq9xw+1zqWPb48dbbaz0A9OQAAIBBVSUee5PmvN82PN+sPfW1lJFKIl29700Kzq1X8wjoth4LTm1b+s1w6dFX73rebXNqjkXp/nqTgnPtyJW7SS7xuMd5PJF4Zzzelh7J8CXrAaAnBwAADIogBwAADKlouGrV3dvSTb40HFQ6X0pNF693fpP1OmqGAkq7pb3Dai3blg6ltHSP18x/02NuGO9QYMuwZ+kwpndoalM5ufVbPjchhBPiP+j0zni8C5vOQ+0yAGn05AAAgCF1m/E4ZembaS6Zc9M2OS3fcDz1er9hlZa7pKZHw9t74H1kdVP5rdvkkoK3dR69ZXu3a0kgr2mvNwm6R8L+SEpnPN51Erf32jpIbQYOOnpyAADAkAhyAADAkLrNeDy3PjPrpplbU3LbLK3fMntuaXmpfUzNJLt0DHJ1eNvZwjMbrredqXLnUu1tqSNn6TzmZqStPU+b9tHbzpbz2bKtGTMez/X4XHnVnP9dtg99bet+h83oyQEAAEPqlnjsfcSx92O0S+uVPuabej/3SLM3IbSm/tL1Sh+zXpLb721+kyytt/cjty2PcnuVPlbecg54/LjMQUo8blkPB4f3vk0PXX/05AAAgCER5AAAgCEVDVeleIdZvHLDES1DFSneOShKu/hLh8FKtAzdpdrUMrSYKze3LFXH0ra7mDsnt13L8FLKtq/VTctyw18nIu/w5raGEWqGWve7zSe6muPeYxj9MNqvlIc5enIAAMCQmntyejvIUW3Nt/jeydAtcr0BLb0Iq229PQY1CcXbui5KeyJ7f2OrOd89vrWvb3si/N9VKTW9j9vibctBavOJpPRc9Or13ZbeozClD8XsqteRnhwAADAkghwAADCkqhmPl+RmckzN9JibdXa9O3/9Z2lm25ylenMzi3rryrU51ZbU8VlaliovVW8PpbOtevc7tc3SOShty9K1krvevGWUXpdL5yd1PZaez6U6zE6MGY9Xel9XXkufSU+bltqOvryfT6/9Oo+5a2/pb6233Fx5Kbs6BvTkAACAITXPeJxLzlxKTk3xJml6tSS7lta79Phu7r2WBNzcPnrbkkuUrTnupYnMLeeiZb+9WnrCvI+K5o6BN4mxtE01246u5XjsYtvaOnpcwwdNz897zfE5yJ+dbV5PpUnG+4meHAAAMCSCHAAAMKSqeXKWughz3fM1QxCpbUuHhrxDBblhG+/6S/W2zKWS4q0jN/RROpRVc35K981bRmqbpf0uPT4t5zF1zEq7e2uulW0N3Y2k9LPoPde9l/Vo89LwuHfYu0XLEHiOd+i/tLxe52yXepzvmhSNHuVtAz05AABgSEU9OXt7ezp27FjVN+yUlp6X1HqlvQKpcnJRaE0CsDeC7dHj09JblVu/hvdY9E5CTr1Xmghf2vZe34RaeJObN217os54nOL9bG1rWU3PrWdZ7/a23Nt2UV5LL0JpefvdeyP1Pd+1dbaWtw305AAAgCER5AAAgCE1z5OzchDm8SjtnixNKvMOo3iHl0ral1u2H/Or1CTaeYfkPOvP3+8xrOZdz5twXbMspXTfDkK3+WFSek32Xq9lDqoe7fPeZ1vuxwe1vJzS5PODaL8SgJdSBVLrbRM9OQAAYEhVj5Cn9I7UevWGpMpbL7dFr8cyS/expbfKG+H3eMzUe9y9dbT0kJSeq5oEbe+3zdwy77fHml4lb7tOZD16FlJ2fZx79NLmymrp0eidnOq9B7X0OufW65Xc3PM4bjMBuOUBmNR62/xs0JMDAACGRJADAACG1DxPTulwzFzLeqXdkz0S7GoS/bzLegzF1czfUNtNuNTV2NKN2zI01CO5uTTJeIl32C9Vbunw4FK9uKaWz+xBOq6lD06U3gN7zQXVcj3Xtq8lBaD3ddE7DWPO86BMby3pGrtKjKYnBwAADKnqEfKaRMbeCZneaNn7rTyntK7U+zW9LLlyl7R8w8i95+2F6vHY7Da/9aSs19G7t2qbifOlZa+vf/bZZ+vo0aN1jRtMSzLlfvG0uff1t/T5aElkbmlLaXlePXqCe23rKW/XvY49zvc22kxPDgAAGBJBDgAAGFK3/6Bzztvl1mOeE6+WIbaWRKnU0EeuztS2NUN3LfNG1CYseuv3rl+TNO3tAva0ryWRe5tJ6i3LNrXlRP8POluGsw+qnkO83oT/3on0vdcr5X2gptfQS89hpf1KjK95yKj30N0cPTkAAGBIRUFO7v+u6i2EcK0fM7vqJ/X+6idVzlId87I31bEkVUapmm1T65fud6q83HFfKrf0mC2d71Qd62UstW/p+skd99y59bZzaZnnXCztm3d/1rfd29tbXO9EUfvZPWh67Efus7iL8pbW895n1surkaqr97VS8pktOXben9r2tvzd2vZnjZ4cAAAwpKr/u2ppfM37OGHLOFzp+GxNLoW3jlwZ3m8ZuXZ68zuW3kvV4RlTrslNWbpGvOV4yqsZN/dcF6nz6L0+ltSei/VtcstS76XqHaHXoqeW/KZdWrpXrLe55Zz3voa85dXcy9frWDqPLaMTLddKS1s89W7zGHuvgR5taUVPDgAAGBJBDgAAGFLVcFVNF5V329x73u5Z7zalwxw1QzSlXay9LQ19rOyyfUvnrOWYeYf7endbr9dVMzzaMnTneW/TeiHwCPmc9zOz35/tOU+bewyvliwrHXqpWZZT0+ZSLeX1vqY853u/9n+b188SenIAAMCQqoKc+aNfqUcCU8tSP6lH8ryPAO/icbVcO1Na2t7yqGZq/fW2e3qh1svJlVtzLnLlzeWuH295S9fZel3rvXqtSW9LxyV1flLLvG3PHcelZWY8Qj7X4/N+EPRsn/c+kroma+5pPdpa+l7v9bzHrKZeTx0152yX+51ar+aYLaEnBwAADIkgBwAADKnq/65qUfMcfI+EzFS9S/WnuvXX3/OW0SvZNrett33e5NzUe976U3oci16Jut5rpbZN3mWpelPrLV2/pcdifX0Sj8u1JIn2lrp2djU8NLdU534fs5akYO96uX2r+fvn5al3m3XW/l3dtKwXenIAAMCQmh8hTy1riW5TvOV5eyi8vRypRM7SOnv0SqTaVBMNe3sKcuX16Bny7o+3vJ49GiVyPUS9j2fJ+7myU+8dhJ6Ig87bG3eQrLe5Vy+Ct6ezd72lSntVe/XIlrZvl7ZZ50H9PNCTAwAAhkSQAwAAhlQU5Bw5cqT6Gf7SuVRSz+6nnqHPze8yLye13tIz+bn3vPO25NpZc3xy7fTOzbIkNb9Fz3JTbV/an9R6uXOw1HbPXAxL9XuuD08S3vrP0jXt2ddN19mm+pgnJ6/0875fcveUpeu5B+/9eBdKz1nq74F3GQ4uenIAAMCQqhKPvZYStHIJwEsJX0v1baqrhjcJzbtt6Xo1Cau595eOd6nSfUzx7nevxEZP4nFNuS2Jz6VJyy0OWg/EYXTQk5F32ZaDtN85qXPmuRfg8KInBwAADIkgBwAADKl5uMo7RFTaPd4yLFEz54p37oeS7TaVUTr04R2a8u5Xqq0tQ0mlQ0QpNV3FpXPReIcbS4dWW67LlqHDpfJKhhGZ8bgcQ36HD+fsxENPDgAAGFLz/11VmhRbk5zqTfAr7dFo+WbtTRLt8a3cm1C8xJtg530c3FN/716oluPZcq5aeiJTZfRIUG5pC5b1OofYH94HRjiPY6MnBwAADIkgBwAADKl5xuPUTKve2SB7zBqZqndpRsrSer3l5mbT9M6A6922ZQbcpfU8x6flWCyV5zl3S+97j7f3+KTq9Oz/puvT26ZdzKxqxozHXsx221/uXtWj3KVlGBs9OQAAYEjNj5DnZpBMrddj5uFNdXjbV/r4cGmdS/Xnyp7X35KMnCqjx6PepWXVJHe3nIMej+x7lSaf17Qjd120JDIDB5n3IYlS3vsxxkFPDgAAGBJBDgAAGFLVPDlL3enb4p1zpWW25G1pGbZZGpbYVpduqk2leszrs2mb0rmXSvej5rj3mC255hpYupZymPEYu7bNe3DP+a5w+NGTAwAAhtT8CHnqkd0ejwQubZuq1/tIb+5x4Nwjxd7HnDfVW/pIcWq9peO8jUckax7vLn28OtX21H5727rUvpzSR0+XjnvufHuvX++xSNWxdF54hBy7tnSd9qxjaRnGRk8OAAAYEkEOAAAY0lYSj1vmY8kt8yZ61swfUirXJu+cPHOlScbe+WKW5OY58h670n1camfp/DM1ifC77LL2ztXk/Tx4j6OnTSskHmOXSu/9Uvn8WaX3FoyJnhwAADCkop6cTYnH61p6AFoi7N69ArnyS3uVvD0+Sz1iufb1erS/dJtcO709REuPUqfqqnkUvUSPHpNWpd9egcPA2/PveUigZ3kYDz05AABgSAQ5AABgSEVBzt7eXtcu8tx8JLn5apbmnampt3a9pTlSNs1LUjIPi3dOidSx6DGHjneOH+88Oqll3uO4NGdQTs0cMuvHoHSunZr9SpWzJHVMvMeWeXJwUC3d81PrAXP05AAAgCEVJR6veB+R3nVymTfJt7QNPRJxveUtHduax+g9deS2XUoK7l3/0jaeenuf79Kk8tLH5Wva6S2bZEscBjXTQ3iX4cRFTw4AABgSQQ4AABhS1XDVNmeTTXU1emeMra2rhbebtGaoorR9LcNA3qHFpXpL68rpPZS1tO225qTxDtV6j0/pcfQMiTHjMfZL7l7esgyQ6MkBAACDag5yUo/W5h718z6ym+J9BNj7yPOmsjett/Q4Y+6R3ZrHl1uOVcrSY8Qlj2Hn3itpZ+r85B6l9l5TpY+XL9WVa5/3MdeW9vV6RHbTseURcuy31GfMe0/pfa/EOOjJAQAAQyLIAQAAQ9rKf9DpVZpomaq3Jum0dM6THsm5Kd7EY++8ECVl1/LWv4vk85yac1V6DbTU1ZLk7N12tZ7nXJF4jIOOISjUoCcHAAAMqagnZ29vT8eOHXN/c+0xW2tJObmyW3oedvk4dmpb76PzNT0+ucfzc9t5eY+193y3PF69VMf6sUj1JOXWX9Kyj72nG+BbMYATAT05AABgSAQ5AABgSFWJx0vd7qlhltz63mGWXXSxb2t4aRdahk287WxJAPaWVzoE2mu40zvUlSu3tH3e4cSlz0NLm4ESPVIUgF2hJwcAAAypKMjZ29vbOHtwanbYFO96S1Kzv+ZmDfYqnWGztL2tM9vm9nvpvKTW99S1dExSs/im6k/V5ZlReKmO1Dbbmil1ab3cbMlLUut5r/OWNgPAqOjJAQAAQyLIAQAAQ9pK4nGONxHVW4436XRpmSeZs2Ym2pZE6tJj691v77beJGPvbMCeMjZt2yP5vPS6SF0LpXPtpN7b1KbSti/xXvsMY2HJNq/Tben1cAQON3pyAADAkKpmPJ4r/RbfK7r2Ppbr3dbTy+F91D1V7jYfc25pX6nSR8h79SK09Ex560olRa9v27u3KlXO0jHLta+l1w3YpOZ+s1/W21LzOcE46MkBAABDas7JWdnm+Gdp7kOvKL20l6r3t4TScnp/2/dqyUXyTsCX60lJtWWpx6fHhGbeyftS5bYc4/3uzcOJraU3s7fc9ez923CQeqHQHz05AABgSAQ5AABgSN1mPJ5LzSLrnc01xTtTbq7+pfatlzsvOzWLba5NXi3HsWXm31y5nuTcTcd903FZr7fG0rH3tG9Jyf4vHffUeqlrK9eOTddebtua9cxMe3t72f0GDgrvZzAl9/cAY6InBwAADKko8XjFm9TpTSZdqiOld/Jbz0fha5JjSy0lk+bq6JUAW/LeXMsj5Kly1nupNtXR49rzXgu9rs/S47Pec+RZL4Sg48ePF7cNWKmZsqFHHZ73Uui5OXHQkwMAAIZEkAMAAIZUFOSs5snxJuzOeRN1U4mjngRgz/BWaYKZp66WZNKWdnr327ttab295I6x95h5t/UmCJck7HqORek2S+t7l+U+h+vHh8Rj9OD9zPaqw7uM5OITFz05AABgSFtJPPbyJonm6lpKtPTWm1svVX5Nku1S+T2VJn/n9q0mYbfHvtVcFy3fHNePWY9E9nk5qWvVm6S+VIdn/U1tIfEYNVLXX++eG28dve89GAc9OQAAYEgEOQAAYEhFw1V7e3s6duxYcnigZk4c79w0q/VKh6M2tSW3Xm6bmn3sMSfONnnPX2v5m+rIqRmG6jE/TeqYlMw5U9K2mmHZ1Hq5zwjd99iVXkPHpeV5/ybhxENPDgAAGNLOEo9T0XVp70FNT0rqvdpv6t6eKa8evUublPac5fat5tuZ9xinlCakt7SltNfEm+he8y3S2zNUep0Du7JfScG5eunRObHRkwMAAIZEkAMAAIZUFeR4Z5BcmnEyN0tr6ey0S+ulZsP1lpdr09I+tmzr3e/c/niXlZ7TpX1MtclT/qZyvbwz/+bKzR3PpXamZivOlbF0XeS29Z5bT5uZ8Rillq7d3LLW+krq9X4mMCZ6cgAAwJCqEo9TapJySxNMvcnDvZJde7QpZ2m/vUmnLbPneh/jL3mvtU3e9XOPTfd+lLolibHlOu9Rx9ym+pjxGD3sOvE49zkiGRkSPTkAAGBQBDkAAGBIVUFOKqmyJbmsJtHSm1TmbVNqP9brSO3jUv25BNSl/U7Vmzv2qW2XlnmTlmsTXGv2p7TtSz+5tuT2zXucUusttd1zvW06t6n9zy3zXHskHqMH7+euRu5e4f37QzLyiYeeHAAAMKSixOMjR464e0Wk8sRZz/ueOlqSOb3beustrX8psdabaJdbr6V9ufJ6JVJ7y0uV0Xt/c+W2XGc9ktlrjgUJl9im3j0j3vJ6r4dx0JMDAACGRJADAACG1PwfdLYoHfpIvb/UltJhm9y2NcMiOb2HY0qHg+bl1Qyn1ZbRMqyWaov3Gki1IXfMerczVb93W+/Q09LnZlMdzJODFt7PWOnnZ2no1rOs1+cYhxM9OQAAYEjNQU7pY+PeR6mXHjP2tiX3qK73EdxUud5HlXOW9jtXb2ofl35ybcitl2rnUpu8+730vuf8eHkfzc5tm2pnar2lY9ZS7/p7S5+lpWuBR8ixTS2fu6XyPMtS27bcR3C40JMDAACGRJADAACGVJR4vLe3p2PHji0muOaSJLc5l0suCa1H4mZNGblyW5al2tBj3pZ5OS1dud6kw110F3vnmqlNSvReb96k8prk4aV25dpC4jG2ofdnu+bBhm20A4cLPTkAAGBIW32EvOUxvaVtU9/OSx/D9j6mnuLtZdnWY5RL7Sz5Fu9Zb319b+/S0rZePR51Xyov9w3Q2xvUovYarGkL324xKq5tzNGTAwAAhkSQAwAAhlT1H3TWJGv2HsrJDZu0JNsu1ZGTG0LbVG+u/pylMnY5hFWaZOwdfuvd7exNCPe8N3+/Jrnau981ie25+nsnNwPAQUZPDgAAGFLVI+RLUj0gpYmbu3j0ealHytO74n18OLWNt0es5lt3aRJ2aW9a6eP3m+rqXV5pL0euvJpry3vttySur7+Xqr9mPQAYDT05AABgSAQ5AABgSEVBzirxeC73H56llqX+s7al9XJS5c2Xpcrw/sdsqfaVtjNVl2e/NrUjtW3LslK5Y5dqZ+r8pLZJHcel/V7fr/Uholx53msgx9v23DVQs4+5tnuv/fV1+A86AYyInhwAADCkqhmPWx75TqlJoOyROFmaoNwyA25Le3fxuK/3UWrvo/YtxzbXJu+yTeXk5Pa7tKyaduba5E2Ir0nYl8T/XQVgSPTkAACAIRHkAACAIRUFOXt7e64k3lTyYy4hMpWkubReSkuScS4RNNWmHgnF3uTqmkTZXNK0dz9yxzN3nHqdH28yuzfBfel41yaEp8qoSRZvaVNuH5eWkXi8f2rOV24ZkHIiXyf05AAAgCFV/d9VXkvrlkaWS4mbpUm5LeWlEmZ7JFfXHDNvcmouKTeXVN2SbLu0Xsq2juPSsU093u5Rk4jvqb+kDaXrk3i8v7xJ7aXLtvVgAg6P3D1oFw+xHDT05AAAgCER5AAAgCE1Jx6vLCV6ehOPc7zr5epqKS+1bWmSdYkeib81ibq5+mv3YVOSc48Ebu+2S9debeK4V00C+Xpdm66p0kRmklYPHu+1UXsN4cTScr8ZCT05AABgSFWJx0vJSz2SL1uSaJfK8ybg5srw8tafk9p2KSExV29pgnAuKXmu93pLWhKUW3qK1surORdLZa9banvpNXqif7vbDy3nv2Y9nNi819vo9wJ6cgAAwJAIcgAAwJCKhqv29vZ07Nixqm7/3BBRqmu/ZRhsaTjEO49Artzc+96hkl7DOz2WeZUOES0di5bzXHquSufn8V6XNfWWDqN6h79G73oeRelwcsuwM05Mve/9hxU9OQAAYEhFPTkrNd8qvd9mvRFnrufD+4295dtxj+Rqb5u82y6Vt1La0+Ttjak5Jj17kpbqKD2Ppb0yS3UtldM78bS0xxC7472nea8TzitWSu/9o6MnBwAADKkoyEn931W7nLhtrnQSOe9EcLnJ9pbK6L2P3mObap9329Q+5urPrb90LDzvLb2f29ela8XbZu/5TpWbq2vT/ub2u3TbknM/34b/hXz/lV47S8twYqu5t4yInhwAADAkghwAADCkqsTjmkTGluTUlu620npL21eTKNySEJpar0cdpW1ZWt/7qLm3Dm+9Pc5z78faWxKeWxJKvdcjDoaWhHL0tYsHIg4K7/34sKInBwAADGkrj5B7e0+823q/see2TbW/d9S6i16t9bo21Vf6mHquvJpzm9o2VX7LI9c9jqO3PO/j9C1t8ZyLuZr1NtV7/PhxbzPRQelncfRv2wdBy7k4LOfH2xN+kPehFD05AABgSAQ5AABgSEVBzt7enmvej9z8KjXbrqTmI0m97613Uxs2/XjX9yqdo2VTvaXz6bToUcZc7th653zxHvvUeqXHsea68M6X4732vfu9/lnJrcc8ObuVOte5ZTX3l8PCO9/Vtuvf5rKDInWM9+u47wo9OQAAYEhFicerGY9zybxS+aPHqfeXyitNkErV633c1/Pepja1JHd597vH8emhJhl6F49q5pKGc0mENddqS/2eMmraMuo3tNEctiTWbdrvBNjen+39Po+lD04cpLa3oicHAAAMiSAHAAAMqWi4am9vT8eOHUu+VzNXSMuwVop3+Kl0ThFv12Vp21K8Q1i95slJaRkKLHlvqU0110Cu3tLhMu+QaW8t3cPeoVLs3tJ5OCznpnSY31vGQeLdx9K/TQdx6Oegn4tW9OQAAIAhVc14PNeSCNujd8ObIFUzS2Xq8eD18np9w8+1qWYfc+uXftv3njNvb0xLb1VNr17t9djr3PZIZK45jyU9ccx4vFs9rvH9VvPZ9m67y2PQci563Ld3qddDMYcJPTkAAGBIBDkAAGBIRUHOap6clNSslfOf9Vld58vmWmaMTbXFK7e+d39ybUrtV6qOpW1T6+W2TbU9dfxKbdqPWt79Tm2zVF7u/G3axlNubllqm5rrMidV3tL5TrXZjBmPdy11brzL9sPStbat8np/ZnJtWVrWo7xd7E+plvvIYUBPDgAAGFJV4nFNQmbu/V7lrcrpndxcss6mOr1JaN6k05b98ibTefVoi7ed3sS5pUTAHom/Le1MbduyPy3JjiQeHwzbSjz23nt731NzvOXtVwJs7p7hfQhgqbz91uvBioOOnhwAADAkghwAADCkrcx4nOLtsvcOC5QOH6S0dIX2HqLJdWfWdCF66/Us87apZWjMO+SzNEST23apS9lznGv2OzeU1NJt7922ZfgU2+EdQvW816IlzWBp/R736NKh4152MWSY+5t4kBz09i2hJwcAAAypW+Kx11Kk3yPCb0nw9JaXa1NNb5U3qm9JamtJTt1U1lLbvNt6z0Vr2Z5tS3uuvO2t6fHJqal3xMTCw2Lpc1eaUF66rLQHt3VZSo86ettGIn/P8nbpILapFT05AABgSAQ5AABgSN3+g845b1ds767LlmEDzxBSTXdp6dBY7y7ZmuNYOgzTu7u5tMu0V1f6enm7GA6qGabr8Xlg2Gp3vOdm28tqtitNHu61bY86vFraUprK0Ks8+NGTAwAAhtTck1OaRNvy2GuqnJZE5pbEMG9dLZG+t4wevSY9vvXXPJbd49tMr0TN2kfIa9q0/l4vS3XQg7N7LQ8GbGu9bfUYeHtQa6ZO2OVj2L17POlB3T/05AAAgCER5AAAgCFVBTlm1/5v2Oc/qfVyy5aEEDZ2R67eW3o/1X6v1Pqp/fbWlVpvvh/rP956l4536bFY325eR03bU+1MbeO5tmrUXCsl7y11mS/ta2o97/XT0mYz097e3mJ52C7vNZ777KSuxaXPZalcGd7P6jbX63Gv6H2/KT236IeeHAAAMKSixOMjR44kewI8cgli3sQ4b1Kl5xt1jZZvQN4E6V5JottKCPeusyqvR+L1Ujk110BpXb0THL0J8V6lidTr1+Px48eL60SZms9CaSL9Li197koTipf2IVdHy32z9Bj3+nwe5HM7t63k7l2hJwcAAAyJIAcAAAxpq/PkpNTM2VDS7b5pWa+2eLb1Dqstlduj7S3ltdTVsxu5hreulmE6b/d1j277luG8pfXMTGeffbaOHj3qKgO7kzrvBzkx1dveln0oraPmPrKL497jWPUeQurxt+mgoScHAAAMqaonpyVBq1cEmEu0bEl+m/Mm6nrWb+lZqFnmVbsfLT1J3l6Jml66lt689W179UK19ESmysuVW3sNkHh8cCxdO9573y6V3vNqPlulD2yUHottfsZS25aOBiytt+1zf5h6b+boyQEAAEMiyAEAAEOqGq7yDvN4ebsEe8+PsNSW2sTflvl8eiWY5tbvMRxSc6xLh/+85aXKqbkuPO3zDhHUXEct13Rpt/phTSLEwU489up1/fW6l5SU37ts79xWcz2vgZo6D9M9g54cAAAwpG6PkPfugVjqIcrVu7ReS09GbjvvN/HSJNqlNvVI+qs5tp71U1p67lJqkg57fBNpeYQ81eOUWr/02NY8Ss4j5AfXYZkVN+Ugtukgqpm2ZFu9Vb0fbDkI6MkBAABDIsgBAABDap4nJ2UXc9jk1uudjNY7gXOXiYMtM+G2JBm3dLP3GHLqnfydG16qGbL07mPvuUZS6zNPzuHScj3h8NrFud12Ivd+oCcHAAAMqSjI2dvbk5kphHDVT4qZFX2zXCpv/v7qJ1fH6j3P+7mf3P5493G9vevHL7c/qTZ5jvt8vaV99JTnPd+p/fHuY8sx8+7PUh3r7y2dC2/be0jtf696zUx7e3td2ont894DANCTAwAABkWQAwAAhlSUeHzkyBH3UIDUb26HljlCUgl5pYnEqblkevAei5p5WFLl5hIWl7ZdX1aTfN6y3srSNdWS+Lu+rCWheEnv/c61OVXv+rkl8fhga0kyB05k9OQAAIAhFfXk7O3t6dixY8lekZZvFzU9NaW9K0uPFHu2qflmnyu3ZlbJ0v2o6S1KJd3mlD6mvtTzstS+0rbU9ir1mnbA2+OTa2fva29T2Tj4vD3lAOjJAQAAgyLIAQAAQ6pKPO4x9LNpmbeM0qTlpZltPYmbNfvgnVE3xTus1Xt4JdfOXLf40vFM1ZUbdiwdBtvUlly9qW1K9zFVR03SfY+hyF7J7Di4etxLgRMFPTkAAGBIVvLt/9xzz/1bSZ/eXnMA7JM7PPjBD77Ffjdim7h/AcPaeP8qCnIAAAAOC4arAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkAhyAADAkIqDHDO7wsw+MPt5xjYatlbnjc3sKRXbPdvMnr627IFm9p61ZaeY2efN7LYbyjnTzN5UWr+jfaea2SvM7JNm9j4zOyOzrpnZu8zsu2bLHmNmb+7drhJmdp6ZHdlyHXc1s/eY2T8lzueNzexVZvZRM/uImd0vLr+Xmb03XqN7ZnZ0m23E4cD9qx8ze46ZfSjev340LnukmV0cP4//LS67q5m928w+aGZ/amY3z5T5HfGzbvH3k83sQjP79t7t99rW8UvU82Yzuyge0xeb2clx+WPisivn91ozu5mZvcPMLjOzF24o8w1mdkli+c+YWdh0LszsSWb2ifjzpF77uB9qenL+MYRwr9nP87q36tpuLKn4JrHBOyXdzszuMFv2EEkfCiH8Tac6vH5Y0pdDCP9c0gsk/eqmFUMIQdK/lfTrZnaamd1A0n+SdFZt5TFwOgy9eV+S9JOS/nPivd+Q9OYQwl0l3VPSR+LyX5P0nBDCvSQ9K/4OcP/q572S/jdJxyQ918xOkfRFSf9K0t0l3d3M7h/XfUII4R6S3q3pPpYUQjhH0qc13Rsl6Sck7YUQ3l3byNiuw+CxIYR7ajqmt5D0mLj8EkmPknT+2vqXS/oFSU9Xgpk9StJlieW3l/RQSf9zw3Y3lfSLms7rUUm/aGY3Kd2Zg6LLHzgzO93MPmZmd4m/v3wW2V9mZi+Ikei5ZnaLuPxOMXJ9v5m908zuGpffysxeGyPai2IE/zxJd4rfvJ4f1/tZM7sgfmt4zqwtzzSzj5vZuyTdZb2tIYQrJf2xpMfNFj9O0svN7Gj8FnFh/OZxre3Xv12Z2SUWe2DM7Almdjy28yWrSDzj4ZJeFl+/StKDV99gUkIIl0h6o6T/oOkP9+9J+rKZvS4eh/ea2bfl2hl/PmZmv6fpw3P72TrfaWavnP1+1TcYM/tNm3pEPjQ/3mvH5rLZ60eb2Uvj61uY2avj+brAzP7VwnFZ3+8vhBAukPT/rdV3uqQHSPqduN7/CiF8ZbWZpBvF16dL2vUfABwS3L/q7l8hhD+JX75OknTltCicH0L4qqbP32mSLg8hfDSE8Jdxs1M1/XHOeZqknzOzb5X07yT9BzP7AZt6gi4xs6u+DGbuOS+1qTfkfVr7ghPvk986+/08MzuyD/f/awgh/EN8eYqk62o6hgohfCSE8LHE+l8LIbxLieNp05fgn5b0K4mqXiDp36/KT3iYpHNCCF8KIXxZ0jmSvrNkXw6UEELRj6QrJH1g9vP9cfl3SHqPpg/cm2frB0mPj6+fJemF8fW5ku4cXx+T9Pb4+hWSnhpfn6zpD9QZki6ZlflQSb8lyTR9wN6k6Y/dfSR9UNI3afoD90lJT0/swxFJF8bXp0r6gqSbxm1OicsfIunV8fWZkt4UXz97XqamQOEMSXfTFIBcJy5/kaQnxte/LelIoh2XSLrd7Pe/kHTzheN/fUkfi/t5qqT/KukX43sPkvSBhXaeoemGdN9E2adoiu6vH3//TU3fwCTpprNzcp6kb4u/n7faN0mXzcp6tKSXxtd/KOn+8fW3SPrI7Dz8dsG1t75P95J0XNJLJV0Yj/Oq7XeL+/LXkj4j6Q6l1zo/4/2I+9em+0Lx/Su+d514LM5aW/7Lkl65tuxhmnpab+w4Tz8h6R8k/ZCk28bP8i003aPeLukRcb1N95yXxuN6cqLsp2nq5ZWk20j6WHy98+OXaNtbJH1Z0z3z5LX3zkuVE4/RC9eWvUDSI3Xta+/hkn4jvv6UEn9vNPUM/fzs919Q4jo8LD813Xj/GKYhgGsIIZxjZo+RdLamYYOVKzV98CXpv0t6TYwyv13SK2cdF6fGfx8k6YmxzCsk/b1du6vsofHnwvj7DSTdWdINJb02hPB1aRqPTO1ACGHPzG4QI/W7SXpfCOFLNnXjvczM7qzp5nadpYMx82BNN6kL4j5dT9PNRyGEHykoJyuE8DUze4WmD/c/2dQd/H3xvbfbNE57o3wp+nQI4b2Jsr9hU47P95jZqyR9t6aIX5Iea2Y/pukmcxtN3dEXO5v9EE1d16vfb2RmNwgh7ElqOTanSLq3pJ8IIbzPzH5D0jM0fSh/XNLTQgivNrPHaurteUhDXRgD96+02vvXj2u6n5y9WmBm99T0B3aeP3KSps/gvw5X97bmnC3peSGEl5rZwyWdF0L421jWH2gKCl+3UMYr4zlY98eS3qppSOaxmnrRpSkg3df7fwjhYWZ2mqQ/0HQtnVPQBklTPqKkO4UQnmazPE8z+yZJ/5ema++E0W2sMl7Ed5P0dUk3kXTphlVX3ZtfSd1svNVJem4I4SVrbXhqQRkv1/St7W7xtTR9+3hHCOGR8eI4L7HdN3TNYb7TZm16WQjh5wra8BlNw0WX2jRufLqmMe0lV8afnE3tlKSvZbb7I01dxF/SNBb+VTO7o6bo/l+GEL4cu4RPS2w77/6cv3+Spp6jpW5qSZKZnSXpR+Ov/0fYnGtwqaRLQwjvi7+/SlOQI0lPkvRT8fUrNX2bApK4f1XdvyTp2yT9ydqye0j607XP+20l/X0I4ROeQkMIV5rZpuGUa6w6e71+T0re50IInzGzL9o0tP/9ujpHaGfHLw5lvT/++oYQwrNm7bvczF6vqdelOMiRdD9JR8zsU5r+xt/SzM7T1Dt2R0kXxUDsdpL+3MyOhhA+N9v+M5p6r1Zup/SxOBR6Jp0+TVNX5A9K+l0zW0XBJ2nqRlR8711hGnv8q/jNaZUAu/r2dK6mbwerzPrTJX1V07eclbdIenL8RiUz+2Yzu6WmxKxHmNn1zOyGkr4n096XS3qCpmj59XHZ6ZpOsDR1AaZ8SlPvgczs3poumlW7Hx3bITO7qV0zOTDlDZr+GEvTMXp7CCHE/Tl3YduVd0p6fKzzTEl/F4/vpnYu+dO43Y9qCnikqRv3a5q+ld5K0ndt2PbzZna3+AfjkbPlb9X0AVNsz71yDQghnB2uTgzdmEsTP5h/PRs7f7CkD8fXfyPpgfH1gyS5bq44YXH/Kr9/SdJ/0zTMN/dninlyM1+W9DPzBWb2XDN7pJYdl/RAM7t5DA5+QNN9Stp8z1nyCk291KeHEFY90js7fiGEK2b3uGfFnrnbxG1P0dSL/tGC/ZmX/ZshhNuGEM6QdH9JHw8hnBlC+GAI4ZYhhDPie5dKuvdagCNN1+dDzewmsRfyoXHZoVQT5FzPrvkI5vPiH5kfkfQzIYR3avqw/nxc/2uSjtr0GNuDJP1SXP54ST9sZhdJ+pCmqFWavn3/azP7oKZI9+4hhC9K+jObkryeH0J4q6Yxy/fE9V4l6YYhhD/XdPFepOnbxQWbdiKE8JHYtreHEFYR/69pekrgQm3u5Xq1pJua2Yc09Xh8PJb34bjPbzWzizVF4KuL9rct/Zj170i6mZl9UlOS2KoX4jaavjF4PFvSfWKdz9PVQVOynUti9+6bNAUyb4rLLtLUtf5RTcf9zzZs/oy4zbslfXa2/Cc1fbO42Mw+rPjNyaZkv8UeFjO7tZldqukY/byZXWpXD8n9hKQ/iPt/L01PnElTkPZ/x+vrP0n6McfuY3zcv/rdv6Tpj/Fd15bdQ9ceEjld1x6avoek9T+w1xJC+Kyme8s7NB2b94cQVoHdpnvOkldp6gn749my/Th+K9eX9Ia47Qc0DXW9OG7/yHj/u5+k/9fMrgo4Ym/Nr0v6oXhfvLtv969pfi8OIXxJU6/WBfHnl+KyQ8lC8PQINlRgdlkI4QZbrWQwZvbvJP3PEEJyTB7AbnD/2h4ze0sI4WH73Q6MjSAHADbg/gUcblsPcgAAAPbDYZjtFgAAoBhBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGBJBDgAAGNL/D9aCgF3R4TtLAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x720 with 4 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "## Solution\n",
    "\n",
    "\n",
    "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": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAD4CAYAAADiry33AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAAsTAAALEwEAmpwYAAATaklEQVR4nO3df6xndX3n8edrGZAEWWHkLiKwgltKgt0F2dtRIzXoKMJAxG5MF7LpYtVMpdjIrkkz1qQ09h/Ube0Pmk6mwAq7FH+jpKAynbWLJgJeWEB+SBkphhmRuToKuG1qpr73j+9n1u9cvnfmzv1+7w8+PB/JN9/P+ZzP95z3nDvzumfOOd9zUlVIkvr1L1a6AEnS0jLoJalzBr0kdc6gl6TOGfSS1Lk1K13AKLfffnsdfvjhK12GJD1vPPvssz9Yv3791Kh5qzLoDz/8cNatW7fSZUjS88a2bdu+O988D91IUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnVuU3Y6XV6qRNt6x0Ccvu8SvPX+kSNCb36CWpcwa9JHXOoJekzhn0ktQ5g16SOnfAoE9yYpKvJnkoyYNJ3t/61ybZmuTR9n70PJ+/pI15NMklk/4DSJL2byF79HuAD1TVacBrgcuSnAZsArZV1SnAtja9jyRrgSuA1wDrgCvm+4UgSVoaBwz6qnqyqu5p7WeBh4HjgQuB69qw64C3j/j4W4GtVbW7qn4EbAXOnUDdkqQFOqhj9ElOAl4N3AkcW1VPtlnfB44d8ZHjgSeGpne0PknSMllw0Cd5MfA54PKqemZ4XlUVUOMUkmRjkpkkM7Ozs+MsSpI0ZEFBn+RQBiF/Q1V9vnU/leS4Nv84YNeIj+4EThyaPqH1PUdVbamq6aqanpoa+SBzSdIiLOSqmwDXAA9X1R8NzboZ2HsVzSXAF0d8/CvAOUmObidhz2l9kqRlspA9+tcDvw68Kcm97bUBuBJ4S5JHgTe3aZJMJ7kaoKp2A38AfLO9Ptz6JEnL5IB3r6yqrwOZZ/b6EeNngPcMTV8LXLvYAiVJ4/GbsZLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzh3wwSNJrgUuAHZV1S+1vk8Bp7YhRwE/rqozRnz2ceBZ4J+BPVU1PZGqJUkLdsCgBz4BXAVcv7ejqv7j3naSPwSe3s/n31hVP1hsgZKk8SzkUYK3Jzlp1Lz24PBfA9404bokSRMy7jH6XwGeqqpH55lfwG1J7k6ycX8LSrIxyUySmdnZ2THLkiTtNW7QXwzcuJ/5Z1XVmcB5wGVJ3jDfwKraUlXTVTU9NTU1ZlmSpL0WHfRJ1gD/AfjUfGOqamd73wXcBKxb7PokSYszzh79m4FvV9WOUTOTHJHkyL1t4BzggTHWJ0lahAMGfZIbgW8ApybZkeTdbdZFzDlsk+TlSW5tk8cCX09yH3AXcEtVfXlypUuSFmIhV91cPE//O0f0fQ/Y0NqPAaePWZ8kaUx+M1aSOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1LmFPGHq2iS7kjww1Pf7SXYmube9Nszz2XOTPJJke5JNkyxckrQwC9mj/wRw7oj+j1fVGe1169yZSQ4B/hw4DzgNuDjJaeMUK0k6eAcM+qq6Hdi9iGWvA7ZX1WNV9VPgk8CFi1iOJGkM4xyjf1+S+9uhnaNHzD8eeGJoekfrGynJxiQzSWZmZ2fHKEuSNGyxQf8XwL8BzgCeBP5w3EKqaktVTVfV9NTU1LiLkyQ1iwr6qnqqqv65qn4G/CWDwzRz7QROHJo+ofVJkpbRooI+yXFDk78KPDBi2DeBU5KcnOQw4CLg5sWsT5K0eGsONCDJjcDZwDFJdgBXAGcnOQMo4HHgN9vYlwNXV9WGqtqT5H3AV4BDgGur6sGl+ENIkuZ3wKCvqotHdF8zz9jvARuGpm8FnnPppSRp+fjNWEnqnEEvSZ0z6CWpcwa9JHXugCdjpdXopE23rHQJ0vOGe/SS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5vxkrab9W6lvIj195/oqst0cH3KNvD//eleSBob6PJfl2ezj4TUmOmuezjyf5VpJ7k8xMsG5J0gIt5NDNJ4Bz5/RtBX6pqv4d8HfAB/fz+TdW1RlVNb24EiVJ4zhg0FfV7cDuOX23VdWeNnkHgwd/S5JWoUmcjH0X8KV55hVwW5K7k2ycwLokSQdprKBP8iFgD3DDPEPOqqozgfOAy5K8YT/L2phkJsnM7OzsOGVJkoYsOuiTvBO4APhPVVWjxlTVzva+C7gJWDff8qpqS1VNV9X01NTUYsuSJM2xqKBPci7wO8Dbquof5hlzRJIj97aBc4AHRo2VJC2dhVxeeSPwDeDUJDuSvBu4CjgS2Noundzcxr48ya3to8cCX09yH3AXcEtVfXlJ/hSSpHkd8AtTVXXxiO5r5hn7PWBDaz8GnD5WdZJesFbycZG9fVnLWyBIUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjq3oKBPcm2SXUkeGOpbm2Rrkkfb+9HzfPaSNubRJJdMqnBJ0sIsdI/+E8C5c/o2Aduq6hRgW5veR5K1wBXAaxg8GPyK+X4hSJKWxoKCvqpuB3bP6b4QuK61rwPePuKjbwW2VtXuqvoRsJXn/sKQJC2hcY7RH1tVT7b29xk8DHyu44EnhqZ3tL7nSLIxyUySmdnZ2THKkiQNm8jJ2KoqoMZcxpaqmq6q6ampqUmUJUlivKB/KslxAO1914gxO4ETh6ZPaH2SpGUyTtDfDOy9iuYS4IsjxnwFOCfJ0e0k7DmtT5K0TBZ6eeWNwDeAU5PsSPJu4ErgLUkeBd7cpkkyneRqgKraDfwB8M32+nDrkyQtkzULGVRVF88za/2IsTPAe4amrwWuXVR1kqSx+c1YSeqcQS9JnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnFh30SU5Ncu/Q65kkl88Zc3aSp4fG/N7YFUuSDsqCnjA1SlU9ApwBkOQQBg/9vmnE0K9V1QWLXY8kaTyTOnSzHvhOVX13QsuTJE3IpIL+IuDGeea9Lsl9Sb6U5FXzLSDJxiQzSWZmZ2cnVJYkaeygT3IY8DbgMyNm3wO8oqpOB/4M+MJ8y6mqLVU1XVXTU1NT45YlSWomsUd/HnBPVT01d0ZVPVNVP2ntW4FDkxwzgXVKkhZoEkF/MfMctknysiRp7XVtfT+cwDolSQu06KtuAJIcAbwF+M2hvvcCVNVm4B3ApUn2AP8IXFRVNc46JUkHZ6ygr6r/C7x0Tt/mofZVwFXjrEOSNB6/GStJnTPoJalzBr0kdc6gl6TOGfSS1DmDXpI6Z9BLUucMeknqnEEvSZ0b65uxWh1O2nTLSpcgaRVzj16SOmfQS1LnDHpJ6pxBL0mdM+glqXOTeGbs40m+leTeJDMj5ifJnybZnuT+JGeOu05J0sJN6vLKN1bVD+aZdx5wSnu9BviL9i5JWgbLcejmQuD6GrgDOCrJccuwXkkSkwn6Am5LcneSjSPmHw88MTS9o/XtI8nGJDNJZmZnZydQliQJJhP0Z1XVmQwO0VyW5A2LWUhVbamq6aqanpqamkBZkiSYQNBX1c72vgu4CVg3Z8hO4MSh6RNanyRpGYwV9EmOSHLk3jZwDvDAnGE3A/+5XX3zWuDpqnpynPVKkhZu3KtujgVuSrJ3WX9VVV9O8l6AqtoM3ApsALYD/wD8xpjrlCQdhLGCvqoeA04f0b95qF3AZeOsR5K0eH4zVpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzhn0ktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1btz70UtSd07adMuKrPfxK89fkuW6Ry9JnVt00Cc5MclXkzyU5MEk7x8x5uwkTye5t71+b7xyJUkHa5xDN3uAD1TVPe25sXcn2VpVD80Z97WqumCM9UiSxrDoPfqqerKq7mntZ4GHgeMnVZgkaTImcow+yUnAq4E7R8x+XZL7knwpyav2s4yNSWaSzMzOzk6iLEkSEwj6JC8GPgdcXlXPzJl9D/CKqjod+DPgC/Mtp6q2VNV0VU1PTU2NW5YkqRkr6JMcyiDkb6iqz8+dX1XPVNVPWvtW4NAkx4yzTknSwRnnqpsA1wAPV9UfzTPmZW0cSda19f1wseuUJB28ca66eT3w68C3ktzb+n4X+NcAVbUZeAdwaZI9wD8CF1VVjbFOSdJBWnTQV9XXgRxgzFXAVYtdx2Ks1DfaJGm18puxktQ5g16SOmfQS1LnDHpJ6pxBL0mdM+glqXMGvSR1zqCXpM4Z9JLUOYNekjpn0EtS5wx6SeqcQS9JnTPoJalzBr0kdc6gl6TOjfvM2HOTPJJke5JNI+a/KMmn2vw7k5w0zvokSQdvnGfGHgL8OXAecBpwcZLT5gx7N/CjqvoF4OPARxa7PknS4oyzR78O2F5Vj1XVT4FPAhfOGXMhcF1rfxZYv/dh4ZKk5THOw8GPB54Ymt4BvGa+MVW1J8nTwEuBH8xdWJKNwEaAW2655Sfbtm17ZDFFXfOWwxfzsX3s3r37mLVr1z6nxtVktddofeOxvvGs9vpgdI3btm0bZ5GvmG/GOEE/UVW1Bdiy0nUAJJmpqumVrmN/VnuN1jce6xvPaq8PlrfGcQ7d7AROHJo+ofWNHJNkDfAS4IdjrFOSdJDGCfpvAqckOTnJYcBFwM1zxtwMXNLa7wD+V1XVGOuUJB2kRR+6acfc3wd8BTgEuLaqHkzyYWCmqm4GrgH+R5LtwG4GvwyeD1bFIaQDWO01Wt94rG88q70+WMYa4w62JPXNb8ZKUucMeknq3As+6JOckeSOJPcmmUmyrvUnyZ+22zfcn+TMoc9ckuTR9rpk/qVPrMbfTvLtJA8m+ehQ/wdbfY8keetQ/35vTbFENX4gSSU5pk2viu2X5GNt292f5KYkRw3NWzXbb07NK7r+VsOJSb6a5KH29+79rX9tkq3tZ7c1ydGtf96f9xLWeEiS/5Pkr9v0ye1WK9vbrVcOa/0rciuWJEcl+Wz7+/dwktet2Parqhf0C7gNOK+1NwB/O9T+EhDgtcCdrX8t8Fh7P7q1j17C+t4I/A3wojb9r9r7acB9wIuAk4HvMDgpfkhrvxI4rI05bYm34YkMTsp/FzhmlW2/c4A1rf0R4COrbfvNqXdF1z9Ux3HAma19JPB3bZt9FNjU+jcNbc+RP+8lrvG/An8F/HWb/jRwUWtvBi5t7d8CNrf2RcCnlmkbXge8p7UPA45aqe33gt+jBwr4l639EuB7rX0hcH0N3AEcleQ44K3A1qraXVU/ArYC5y5hfZcCV1bVPwFU1a6h+j5ZVf9UVX8PbGdwW4qF3Jpi0j4O/A6DbbnXqth+VXVbVe1pk3cw+L7H3vpWy/YbttLrB6Cqnqyqe1r7WeBhBt90H76tyXXA21t7vp/3kkhyAnA+cHWbDvAmBrdaGVXbst6KJclLgDcwuPKQqvppVf2YFdp+Bj1cDnwsyRPAfwM+2PpH3eLh+P30L5VfBH6l/Zfzfyf55dVUX5ILgZ1Vdd+cWauivjnexWCvif3UsZL17a+uFdMOdbwauBM4tqqebLO+Dxzb2std9x8z2Ln4WZt+KfDjoV/qw+vf51YswN5bsSylk4FZ4L+3w0tXJzmCFdp+q+YWCEspyd8ALxsx60PAeuC/VNXnkvwag9/Ab15F9a1hcJjjtcAvA59O8splLO9A9f0ug8MjK2Z/9VXVF9uYDwF7gBuWs7bnuyQvBj4HXF5VzwzvCFdVJVn267OTXADsqqq7k5y93OtfoDXAmcBvV9WdSf6EwaGa/285t98LIuirat7gTnI98P42+RnafwWZ/xYPO4Gz5/T/7RLWdynw+RocyLsryc+AY/ZTH/vpn2h9Sf4tgz2X+1oAnADck8EJ7VWx/Vqd7wQuANa37ch+6mM//cthIbcWWRZJDmUQ8jdU1edb91NJjquqJ9uhhb2HEpez7tcDb0uyATicwaHXP2FwuGNN22sfXv/e2nZk+W7FsgPYUVV3tunPMgj6ldl+y3FSYjW/GBx7PLu11wN3t/b57Hty5K7Wvxb4ewYnEo9u7bVLWN97gQ+39i8y+O9dgFex78nExxicyFvT2ifz85N5r1qmbfk4Pz8Zu1q237nAQ8DUnP5Vt/1aXSu6/qE6AlwP/PGc/o+x78nEj+7v570MdZ7Nz0/GfoZ9T8b+Vmtfxr4nYz+9TLV9DTi1tX+/bbsV2X7L+pdnNb6As4C72z+oO4F/3/rD4MEq3wG+BUwPfeZdDE7ebQd+Y4nrOwz4n8ADwD3Am4bmfajV9wjtyqHWv4HBVRLfYXD4Yrm25XDQr5btt53BL8d722vzat1+q2X9rYazGJxcv39o221gcGx7G/Aog6vB1h7o573EdQ4H/SuBu9rP/DP8/Eq1w9v09jb/lctU2xnATNuGX2CwY7Mi289bIEhS57zqRpI6Z9BLUucMeknqnEEvSZ0z6CWpcwa9JHXOoJekzv0/9EmysSYu9h8AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 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",
      "78.1 ms ± 997 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "\n",
      "Your version\n",
      "75.6 ms ± 2.13 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",
      "77.1 ms ± 1.87 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "\n",
      "Numba version\n",
      "205 µs ± 16.6 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Numba Speedup: 417x !\n",
      "\n",
      "Numpy version\n",
      "167 µs ± 30.7 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)\n",
      "Numpy Speedup: 554x !\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",
      "32.3 ms ± 1.54 ms 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"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [conda env:recode]",
   "language": "python",
   "name": "conda-env-recode-py"
  },
  "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.9.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}