quantshah/mlphases.ipynb

## mlphases.ipynb
{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Project II: Machine Learning the Ising Phase Transition\n",
    "## Learning from data [TIF285], Chalmers, Fall 2019\n",
    "\n",
    "#### Christian Forssén and Shahnawaz Ahmed, Chalmers\n",
    "Last revised: 16-Oct-2019 by Christian Forssén [christian.forssen@chalmers.se]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Instructions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- See deadline on the course web page\n",
    "- This project is performed in groups of two students. \n",
    "- The second part of the project is optional and only gives extra points. See examination rules on the course web page.\n",
    "- Hand-in your written report via Canvas."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Written report\n",
    "- Page limit: 6 pages (excluding title page and list of references). 3 extra pages are allowed when doing also the optional extra task.\n",
    "- Give a short description of the nature of the problem and the methods you have used.\n",
    "- Include your results either in figure form or in a table. All tables and figures should have relevant captions and labels on the axes.\n",
    "- Try to give an interpretation of you results.\n",
    "- Upload the source code of your program as a separate file (.ipynb or .py). Comment your program properly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Import modules\n",
    "\n",
    "We need TensorFlow 2 for parts of this project. Make sure to update your `tif285-env` environment using the following command (making sure that you have the most recent version of the `environment.yml` file that includes `tensorflow` on the last line):\n",
    "\n",
    "```\n",
    "conda env update -f environment.yml\n",
    "```\n",
    "\n",
    "Please note that Keras has become a part of Tensorflow in the latest release. Therefore many of the online TensorFlow tutorials which use Keras can now directly use imports such as \n",
    "\n",
    "```\n",
    "from tensorflow.keras.models import Sequential\n",
    "```\n",
    "\n",
    "instead of the earlier:\n",
    "\n",
    "```\n",
    "from keras.models import Sequential\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "%matplotlib inline\n",
    "\n",
    "import os\n",
    "import time\n",
    "import pickle\n",
    "from multiprocessing import Pool as ThreadPool\n",
    "\n",
    "import numpy as np\n",
    "\n",
    "from scipy.optimize import curve_fit\n",
    "import matplotlib.pyplot as plt\n",
    "from sklearn.model_selection import train_test_split\n",
    "\n",
    "from tensorflow import keras\n",
    "from tensorflow.keras.models import Sequential\n",
    "from tensorflow.keras.layers import Dense, Dropout, Flatten\n",
    "from tensorflow.keras.layers import Conv2D, MaxPooling2D\n",
    "\n",
    "from tqdm import tqdm\n",
    "\n",
    "np.random.seed(42)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Background\n",
    "\n",
    "The Ising model is arguably the most famous model in (condensed matter) physics. It is described by the simple Hamiltonian\n",
    "\n",
    "$$\n",
    "H=−J \\sum_{\\langle i,j \\rangle} s_i s_j.\n",
    "$$\n",
    "\n",
    "Here, the $s_i=\\{−1,1\\}$ are classical, binary magnetic moments (spins) sitting on a two-dimensional square lattice and the $\\langle i,j \\rangle$ indicates that only interactions betweens neighboring spins are taken into account. For simplicity, we will set $J=1$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most importantly, the Ising model shows a phase transition between a paramagnetic and a ferromagnetic phase as a function of temperature. The critical temperature $T_c$ at which this change of magnetic character occurs has been calculated exactly by Lars Onsager. He found\n",
    "\n",
    "$$\n",
    "T_c = \\frac{2}{\\log \\left( 1 + \\sqrt{2} \\right) }\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Critical temperature: Tc = 2.2692\n"
     ]
    }
   ],
   "source": [
    "Tc = 2 / np.log(1+np.sqrt(2))\n",
    "print(f\"Critical temperature: Tc = {Tc:.4f}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this project we aim to reproduce this result (roughly) using increasingly sophisticated neural networks. You will use\n",
    "1. a simple binary classifier (single neuron) that you implement yourself.\n",
    "1. a relatively simple neural network (one hidden layer with 100 neurons) implemented using keras / tensorflow.\n",
    "1. a convolutional neural network\n",
    "1. (extra task) a Bayesian neural network\n",
    "\n",
    "At the end, the results that you obtain made it all the way into a Nature Physics publication just a few years ago: [Nature Physics (2017) 13, 431–434](https://www.nature.com/articles/nphys4035)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We will start by quickly simulating the Ising model using the Monte Carlo method to obtain representative sets of spin configurations for a bunch of temperatures. Afterwards, we do not take the traditional approach of inspecting the magnetization, the order parameter of the transition, and its susceptibility. Instead we use supervised learning to train a Neural Network to automagically learn the transition temperature."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Monte Carlo simulation\n",
    "The Monte Carlo method for the Ising model is very straightforward: take a random configuration of spins to begin with and propose individual spin flips until you fall asleep. To decide whether a spin should be flipped we use the Metropolis criterium\n",
    "$$\n",
    "p=\\min \\left( 1, e^{-\\beta\\Delta E} \\right)\n",
    "$$\n",
    "where $\\Delta E = E′−E$ is the energy difference between the new (spin flipped) and the old configuration according to $H$ above and $\\beta = 1/T$ is the inverse of the temperature $T$. Since $\\Delta E$ only depends on the local environment of the spin to be flipped (nearest neighbors), we can evaluate it locally. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Generate spin configurations and study the phase transition\n",
    "\n",
    "In the python file attached with this notebook we have the definition of a `Lattice` class which can be used to generate a 2D lattice for `N` spins at a temperature `T`. Here, we simply import the `Lattice` class and use the `step` method to generate a lattice after a few hundred iterations to simulate a thermalization of the lattice. \n",
    "\n",
    "At every iteration, we select $N^2$ random points to try a flip attempt. A flip attempt consists of checking the change in energy due to a flip. If it is negative or less than $e^{-E/(k_b T)}$, then perform the flip. After a few steps the lattice with thermalize.\n",
    "\n",
    "\n",
    "\n",
    "#### *You need the `lattice.py` file in the same directory to get this to work which contains the definition of `Lattice`*\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from lattice import Lattice"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1  1  1  1 -1 -1 -1 -1 -1  1]\n",
      " [ 1  1  1  1  1  1 -1 -1 -1  1]\n",
      " [ 1  1  1 -1 -1 -1  1 -1 -1  1]\n",
      " [ 1  1  1 -1  1  1  1  1 -1  1]\n",
      " [ 1  1  1  1  1 -1  1  1 -1 -1]\n",
      " [ 1  1  1 -1 -1  1 -1  1  1 -1]\n",
      " [ 1  1 -1 -1 -1 -1 -1  1 -1 -1]\n",
      " [ 1  1  1  1 -1 -1  1 -1  1  1]\n",
      " [ 1 -1  1  1  1  1 -1 -1  1  1]\n",
      " [ 1  1  1  1 -1  1  1  1  1  1]]\n"
     ]
    }
   ],
   "source": [
    "# Initialize a lattice\n",
    "lat = Lattice(N=10, T=4.5)\n",
    "\n",
    "# Make 30 iterations (N**2 spin flip attempts)\n",
    "for i in range(30):\n",
    "    lat.step()\n",
    "\n",
    "print(lat.lattice) # (or even `print lat` to use the convenient repr)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Illustrate some spin configurations, and plot macroscopic quantities as a function of temperature"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "6it [00:05,  1.09it/s]\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAW8AAADtCAYAAABwM/RzAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAO70lEQVR4nO3dX6gkZ5nH8e8Tj0k2GQcJM2ZBYQ5IJpKwk4FzbsTAzkVEokGXBCRkLndyoUS82AUhKAybwHghLF6s5mIDE0zQFUHBrIPLBkbWQNRzNsY/YRwQJhGM4AGTzEzMROPrRfWwtc053dVd3V31dH8/UDDV5+2qt+qp+p3q6vfURCkFSVIu13TdAUnS5AxvSUrI8JakhAxvSUrI8JakhAxvSUrI8JakhAxvSUqocXhHxKXa9JeI+GNt/visOhQRZyPizdqyfzWi7U0R8e2IuBwRL0XEA7Pqx7Sa9ikirouIxwdtLkbE8xFx96L729Si6l9b3y2D4+DJEW3S1n/Q9qGI2IqIKxFxeoHdbGSB5/yTEfFKRLweEecj4sSItr2r+bAJj4Hpc6CUMvEEXADumua9DZZ9FjjRsO3Xgf8A9gF3Aq8Bt8+jXxP0v1GfgBuBk8A61S/Re4CLwHqX/e+6/rV1/BfwP8CTy1j/Qdt7gX8Avgqc7rquXdUcuB24bvDvDwC/Azay1LzlMTB1DvSxkI3Ce7DRbwGHa699Dfhiw/WcAL4/OHH+AJwHbgM+C7wM7AD3Ttj3tn36GXBf1wdfl/UfLP9+4JuDg3rX8F6m+gOPrnJ4D63nVuAV4JMZaj7rPg7aN8qBud3zjoinI+LVPaanx7z9VETsRMSzEXFsjzaHgbdLKedrr71A9Vu8iSPAJvAt4ADwc+DM4GfvBx4BPj/h9kzdp4i4efD+Xzbsf69NW/+I2A/8C/BPY1axVPVfBm3O+Yj4SkS8AZyjCu/v7dKsjzWfaR8nyYG1JgucRinlninf+jngRarfXvcD342Io6WUXw+120f1caTuNeBdDddzB3CqlPIMQES8SPXR7cuD+V9Q2z8Nt2eqPkXEO4GngCdKKeca9r/XWtT/EeDxUspvImJUu6Wp/7JoUXNKKZ+OiM8AHwSOAVd2adbHms+sj5PmQO9Gm5RSflRKuVhKuVJKeQJ4FvjoLk0vAfuHXttPdb+oiSNA/bfnbbvMTxqkE/cpIq6h+lj1FvDQhOtbKhFxFLgL+NcGzZei/vo/pZS3Syk/BN4HfGqXJr2reUQcr32Je2baPk6TA/O8bXJm6NvqS0Mb2VQBdrsEOw+sRcQttdfuoMHHjYg4BFw7WMZVR4Gf1uaP1Ocbbs9EfYrq0vJx4Gaqe1x/Gtf3LKas/zGqL25ejojfAf8M3BcR/7tL2/T1XzYzPOfXqG5jDOtdzUspT5VS9g2mu6fp49Q50KcvL4B3Ax8Brqcq4HHgMnDrHu2/QfXN7o3Ahxj6Vhc4zS5fBAEfB56rze8H/gzcUHvtx8DHptiGkX0aavsY8Bywb9b7cp7THOt/A/C3telLVPcnDy5p/dcGx/opqquu64G1ruu74Jq/h+r26D7gHYPz/zLwiSw1b3MMDNpPlQN9K+RB4CdUHzFeHWzQh2s/PwM8XJu/CfjOoNgvAw8MLe8Z4MFd1vMF4LHa/J3Audr8NcAbwHun2IZxfToDPAwcovpU8SbVR62r0/FZ79c51Gku9d9lPSepjTZZpvrXtq8MTSe7ru8iaz44538wON9fp/oS8cHaz3tf85bHwNQ5EIMFLJ2IuJbqW94jZYluR6gZ6796Vq3mSxvekrTMejfaRJI0nuEtSQkZ3pKUkOEtSQkZ3pKUkOEtSQkZ3pKUkOEtSQm1fiTsgQMHyvr6+tTv397ebtuF7HZKKQe77sS0IqLTv/La2NjocvVjj99R/btw4QI7Ozsjn3vbZ8t+7s/72Nre3m517rf+C8vNzc2ytbU19fvHPLN5FWyXUja77sS0ug7vrv9CeNzxO6p/m5ubbG1tpT0Blv3cn/exFRGtzn1vm0hSQoa3JCVkeEtSQoa3JCVkeEtSQoa3JCXUepx31/o+VGzZbWxsMGq4WJuhdH3Qtr6rfnz02QKGAs51+V55S1JChrckJWR4S1JChrckJWR4S1JChrckJWR4S1JC6cd5j5N9nPGy63t9xq3fcdzqilfekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpTQ3Md5z3scrONsl9u8x4F7/Exve3t75P4bV5uux9DP+9ia9/Z55S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCbUe5912rGfX42y7Xn924+rfd/MeJz5q+Zubm63W3bWNjQ22tra67sbK8spbkhIyvCUpIcNbkhIyvCUpIcNbkhIyvCUpIcNbkhLq/fO8ux4n7vOiRxs31jf7M5nHmffyM2tbm3nv2+znplfekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpTQ3Md5dz3Otu1YznmPVc0+1rTvz/Puc9+WnWPg58srb0lKyPCWpIQMb0lKyPCWpIQMb0lKyPCWpIQMb0lKqPU473k/z7nvY0W7fp50dl2Pg++yPpubm52tO4N5Z8e8n9U/72PLK29JSsjwlqSEDG9JSsjwlqSEDG9JSsjwlqSEDG9JSqjz53lnH8c7rv8+T3q0vtff+u1t3LPcu/4bh7bjsLsexz2OV96SlJDhLUkJGd6SlJDhLUkJGd6SlJDhLUkJGd6SlFDr8L461nPaqa3sy89uY2ODUsrUU9+17f+o925sbCxgC+ZnXO37bty53bb2884Or7wlKSHDW5ISMrwlKSHDW5ISMrwlKSHDW5ISMrwlKaGYwfOOfw+8NJvurKRDpZSDXXdiWta/FWu/2lrVv3V4S5IWz9smkpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCY0N74i4VJv+EhF/rM0fn0UnIuK6iHg8Il6KiIsR8XxE3D2i/ZMR8UpEvB4R5yPixCz6MUsRcVNEfDsiLg+264Ex7c9GxJu1ffurRfV1lJ7Wf6J924WIeCgitiLiSkScbtC+d9u0iNoP1tN4X/VxPw1bVO3XxjUopeyrreQCcKKU8t9NFj6BNeA3wN8DLwMfBb4ZEX9XSrmwS/tTwD+WUq5ExAeAsxHxfClle8b9auPfgLeAm4GjwH9GxAullF+OeM9DpZR/X0jvGupp/afZt4v2W+BR4CPA3zRo37ttWlDtYbJ91bv9tIvF1L6U0ngCLgB3TfKeaSfgZ8B9DdrdCrwCfLLhck8A3we+CvwBOA/cBnyWKjh2gHtb9v3GQTEO1177GvDFEe85S3VyzH3fttiuzus/zb5ddP2H1vcocHrWx8sy1n7cvrL2/3+a+T3viHg6Il7dY3q64TJuBg4De/7miYivRMQbwDmq8P5ewy4eATaBbwEHgJ8DZwY/ez/wCPD5lttzGHi7lHK+9toLwO1j+nYqInYi4tmIONZwe3plAfWfdt9etYj6T6rtNvXCAvaVta8Ze9tkUqWUe9q8PyLeCTwFPFFKOTdiPZ+OiM8AHwSOAVcaruIO4FQp5ZnB+l4EriulfHkw/wtq+2XK7dkHvDb02mvAu0a853PAi1S/he8HvhsRR0spv55i/Z1ZQP2n2bd1i6j/pNpuUy8sYF9Z+5pejTaJiGuoPjK8BTw0rn0p5e1Syg+B9wGfariaI0D9t+Ztu8zv+UtjNxFxvPZFzhngErB/qNl+4OJeyyil/KiUcrGUcqWU8gTwLNW935XRsP4T79shi6j/pNpu06qw9jXzuG1yZuhb6ktNNi4iAnic6qb9faWUP02w2jWqjz3j+nYIuJbqXtdVR4Gf1uaP1OebbE8p5alSyr7BdPdg+WsRcUttuXcw4jbQLgoQE7TvhQXUf+p9u8D6T2oWx0vnpq39BKx9XV++tAAeA54D9o1p9x6q2wr7gHdQfaN7GfhErc1pdvmiAPg48Fxtfj/wZ+CG2ms/Bj42g+35BvB1qi8kPkT1Uej2Pdq+e7Ad11P9Ijo+2KZb57GvW2xT5/Vvsm97Uv+1QT1PUX2auB5Ym8XxsoS1b7yvrH3tfX0oIHCI6krzTaqPEVen44OfnwEeHvz7IPAD4FXgdaovHR4cWt4zw68NXv8C8Fht/k7gXG3+GuAN4L0z2KabgO9QhfDLwANDPx/epp9QfVR6lSrEPjzr/TyDbeq8/g33bR/qf3KwTfXp5G71b7JNXU/zqv24fWXt955i8OalERHXUn1be6RMdutFS8D6r65Vq/3ShbckrYJejTaRJDVjeEtSQoa3JCVkeEtSQoa3JCVkeEtSQoa3JCVkeEtSQq0fCXvgwIGyvr4+g67ktL3d+j/v2SmlHJxFX7qw6vVv48KFC+zs7KR7+NhV1r6d7e3tVud+6/BeX19na2ur7WLSqh6G18pLs+hHV1a9/m1sbm523YVWrH07EdHq3Pe2iSQlZHhLUkKGtyQlZHhLUkKGtyQlNPP/PX7VtH0e+gxGq0haQV55S1JChrckJWR4S1JChrckJWR4S1JChrckJWR4S1JCvR/nPe9x0G3HaY/jOO7VZv01L155S1JChrckJWR4S1JChrckJWR4S1JChrckJWR4S1JCvR/n3da8x3G3Xb/jgJfbqPpn/9/j1S2vvCUpIcNbkhIyvCUpIcNbkhIyvCUpIcNbkhIyvCUpofTjvNuO4247zrrrceSrblz9rI+WlVfekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpRQ78d5O05Xo3h8aFV55S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCfV+nHfb52235ThiSX3klbckJWR4S1JChrckJWR4S1JChrckJWR4S1JChrckJdT7cd7jxll3PQ5ckrrglbckJWR4S1JChrckJWR4S1JChrckJWR4S1JChrckJdR6nPf29vbIsdZdPw+76/VL0jx45S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCbUO742NDUope05tRcTIadS6HeMtaVl55S1JCRnekpSQ4S1JCRnekpSQ4S1JCRnekpSQ4S1JCUXbsdAR8Xvgpdl0ZyUdKqUc7LoT07L+rVj71daq/q3DW5K0eN42kaSEDG9JSsjwlqSEDG9JSsjwlqSEDG9JSsjwlqSEDG9JSsjwlqSE/goum1TEC4BUagAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 6 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# 10 x 10 lattice\n",
    "# six temperatures, 500 thermalization iterations\n",
    "# Plot the spin configurations for varying temperatures.\n",
    "# Display the magnetization too\n",
    "\n",
    "nrows, ncols = 2, 3\n",
    "fig, axs = plt.subplots(nrows, ncols)\n",
    "fig.subplots_adjust(wspace=0.6)\n",
    "\n",
    "for (ip, T) in tqdm(enumerate([5.0, 4.0, 3.0, 2.3, 2.0, 1.0])):\n",
    "    lat = Lattice(N=10,T=T)\n",
    "    for k in range(500):\n",
    "        lat.step()\n",
    "\n",
    "    idx = ip // ncols, ip % ncols\n",
    "\n",
    "    axs[idx].matshow(lat.lattice,cmap=plt.cm.gray_r)\n",
    "    axs[idx].set_title(f\"T = {T:.1f}, $m$={lat.get_avg_magnetization():.1f}\")\n",
    "\n",
    "    axs[idx].get_xaxis().set_visible(False)\n",
    "    axs[idx].get_yaxis().set_visible(False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "100%|██████████| 60/60 [01:19<00:00,  1.32s/it]\n"
     ]
    }
   ],
   "source": [
    "# 10 x 10 lattice\n",
    "# 60 temperatures, 500 thermalization iterations\n",
    "\n",
    "# For a temperature range, thermalize a lattice, then\n",
    "# take a few hundred steps, recording energy and magnetization.\n",
    "# Store the means to plot next.\n",
    "# This takes about 60s with one modern core.\n",
    "\n",
    "# Thermalization and measurement steps\n",
    "ntherm = 500\n",
    "nmeasure = 200\n",
    "\n",
    "# points = array with (T, mean(E), abs(mean(M)), var(E))\n",
    "# with the mean and variance evaluated for a list of many temperatures\n",
    "points = []\n",
    "# Storing nmeasure / nsparse data points\n",
    "nsparse = 10\n",
    "# points_full = array with (T, E, abs(M))\n",
    "# for several different configurations per temperature\n",
    "points_full=[]\n",
    "for T in tqdm(np.arange(4.0,1.0,-0.05)):\n",
    "    lat = Lattice(N=10,T=T)\n",
    "    for _ in range(ntherm):\n",
    "        lat.step()\n",
    "    Es = []\n",
    "    Ms = []\n",
    "\n",
    "    for istep in range(nmeasure): \n",
    "        lat.step()\n",
    "        Es.append(lat.get_energy())\n",
    "        Ms.append(lat.get_avg_magnetization())\n",
    "        if (istep%nsparse==0):\n",
    "            points_full.append((T,Es[-1],np.abs(Ms[-1])))           \n",
    "    Es = np.array(Es)\n",
    "    Ms = np.array(Ms)\n",
    "    points.append((T,Es.mean(),np.abs(Ms.mean()),Es.var()))\n",
    "points = np.array(points)\n",
    "points_full = np.array(points_full)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAf0AAAHwCAYAAACsUrZWAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3gc1fXw8e/ZVe+yZclFtuXeu8GGADZgwIHQIRBqwNSEkPIjIYQ3CSkkpEMCITEdYwcIIXQSQ8CAwcYF27h3y5K7ZKv33fP+MSOxllayJEvalXQ+z7PP7t47d+bOVTlz79yZEVXFGGOMMV2fJ9QVMMYYY0zHsKBvjDHGdBMW9I0xxphuwoK+McYY001Y0DfGGGO6CQv6xhhjTDdhQd8YYxohIvtFpExEHg91XZpDRL4iIiUi4heRU0JdHxN+LOibLklEdolIufsPsPb1cKjr1d5E5EcB+1shIr6A7+tDXb/jISK3ici7Idj02ap6k1uHGBFREclsr42JyNki8oGIFIvIpiD5Q0TkI/dgZL2InFabp6pvqGoCcLC96mc6Nwv6pis7X1UTAl53tPUGRCSirdd5PFT1V7X7C9wGLAnY/zGhrl9jOqIdw+1n1YQSYC5wTyP5LwEfAT2AXwKviEhKB9XNdHIW9E23IyJfF5HFIvJ7ETkiIjtF5MsB+cki8oSI7BORPSLySxHxBpT9WET+JCKHgftExCsifxCRPHddd7i9wQgRuVxEVtbb/v+JyCtB6nWliKyol/ZdEXnN/XyuiGxwe4B7ROSuVu7/WBF5z933jSJyUUDe8yLykIi8IyKlIrJIRNJF5K8iUuD2LMcFLL9fRH4gIptE5LCIzBWR6ID8i0Xkc7fsRyIyul7Zu9wRiCI37SduGxaLyDoROc9NnwQ8CMx0Ry32u+lLReSagHXWjQYE9MpvF5HtwLpj7X84UNVPVHU+sKt+noiMB4YDv1DVClX9B7AdCKt9MOHLgr7prqYBm4E04LfAEyIibt4zQA0wFJgEnA3cVK/sDiAduB+4GfgyMBGYzNH/gF8DBonIqIC0a4B5Qer0GjBCRIYFpF0FLHA/PwHcqqqJwFjgvRbsLwAikgS8464rDbgOeFJEhgYsdgVwl5sfASwFPgB6Am/htFegrwFnACNw2uv77ramA38FbnDLzsPplQb2uK8AznLzwfmZnAwkA78BnheRNFVdBXwHWOSOWvRuwW5/BZgCTGrm/rcJEbnBPdhp7JXeitWOAbaoanlA2ho33ZhjsqBvurJX6v2TvTkgL1tVH1NVH06Q7wNkiEgGTgD/jqqWqupB4E/AlQFl96rqX1S1xv3n+1XgIVXNVdUjwAO1C6pqJfACTqBHRMYAWcAb9SurqmXAqzhBFDf4j8Q5GACoBkaLSJKqHlHVz1rRJhcD61R1vqr6VHU58DpwacAy/1TVNe6+vQoUquoLblu9iBPYAz2kqntV9RDw69r6A7cCD6vqSndbc4FonABc609u2XK3DV5Q1X2q6lfVecCeesu3xv2qWuBuozn73yZU9SlVTWni1Zrz7glAYb20QiDx+GtsugML+qYru6jeP9nHAvL2135wgy04/1AHApHAvtqDBeDvOL36Wjn1ttO3Xlr9/GeAq9yRhGuBF92DgWAW8EXQvAp4JaB+lwLnAtniTPQ6qZF1NGUgcFrgwZC73j4ByxwI+Fwe5HtCvXUG7m82TnvUbutH9bbVC+jXSFlEZE7A6YACnNGWtJbtYgOB22jO/oezEiCpXloSUByCuphOqLNMbDGmo+QAlUCaqtY0skz9R1PuAwJnc/c/amHVpSJSBZyKE8ivamL7C4E0EZmIE/y/G7Ce5cCFIhIJ3IHT6+4fdC2NywEWqur5LSzXlMA6DAD2BmzrTVX9QxNl69pSRIYDf8E5VbBMVf3izF6X+ssGKAXiAr4HG/YPLNce+x+UiMwBHmpikcGt6O2vB4aLSIyqVrhpE3BOoxhzTNbTNyaAqu7DCbx/EJEkEfGIc4nUjCaKvQh8W0T6iTOL+u4gyzwLPAzUqOriJrZfgzM7+3c4s7PfARCRKBG5WkSSVbUaZ+KbrxW7+ArOue0rRCTSXe90N+C21p0i0kdE0oAf4pzOAGcG+rdEZKo4EkTkAhGJa2Q9CYAfOAR4ROQ2nJ5+rQNAf/egp9Zq4DJ30t5I4OvHqGtb7X+0u83aV4P/par6RL2rR+q/ggZ893cuBmfESdz1R7rr/BzYCvxYRKJF5Ks4bfRqC+tvuikL+qYre12Ovk7/380sdx0QBWwAjuAE4aaGfx/DOVD4HFiFM9mthqOD8jycyXfBJvDVtwCYhXNuPXC04Vpgl4gU4VyOd02wwk1x5xycgzO5bh9Or/yXOAGmtZ4H3scJRmtxJ/qp6sfAnTinRwqALTijHMF67LhzFP4GrHDrNsj9XOs/ODPaD4pIrpv2W5wRy0M4BxnPNVXRNtz/bTinOmpfTY3etNTZ7jpfxpmpX44z76DW5cBpOG16H3Cxu1/GHJOoBv37M8a0kjiX//1NVQcGpMXi3DBlsqpuDVnl2ph76dxlTY1edGYisgtIBZ5X1VtDXJ1jEucSx/k4EybPVNVPQlwlE2bsnL4xx8kN6Kfj9PYzgJ8C9UcVbgeWd6WA3x2oalao69ASqvomYDfqMY2yoG/M8RPgZzjnssuBN4Gf1GU6vUXBbqBijAkxG943xhhjugmbyGeMMcZ0EyEJ+uLcj3y9OI9/nFov7x4R2SYim0XknID0KSKy1s37c8AtU40xxhjTDKE6p78OuATnUp464jyM40qc+0j3Bd4VkeHu7T8fBW7BuQ/4W8Bs4O1jbSgtLU2zsrJaVLnS0lLi4+NbVKarszZpKOza5LD7XJ8ex3vX2tYJu/YIA9YmDVmbNNTWbbJy5co8Ve0VLC8kQV9VNwIE6axfiHNpTCWwU0S2ASe6E6GSVHWJW+5ZnElRxwz6WVlZrFix4liLHWXRokXMnDmzRWW6OmuThsKuTRa4f09Xtez3va2EXXuEAWuThqxNGmrrNhGR7Mbywm32fj+cnnytXDet2v1cPz0oEbkFZ1SAjIwMFi1a1KJKlJSUtLhMV2dt0lC4tclM9z1UdQq39ggH1iYNWZs01JFt0m5BX5xnWge7D/a9qtrYLSODnafXJtKDcp/mNRdg6tSp2tIjKDsSbcjapKGwaxP3AbyhqlPYtUcYsDZpyNqkoY5sk3YL+qo6qxXFcjn64R2ZOLfJzOXoB5rUphtjjDGdlt+vVPs77tL5cBvefw1YICJ/xJnINwznaVs+ESkWkenApzj3Rv9LazdSXV1Nbm4uFRUVQfOTk5PZuHFja1ffJTWnTWJiYsjMzCQy8nhu426MMV3HoeJKPtmex668MvJLK8kvreJwSRX5pZUcLq3icGkV52RFctYZHVOfkAR9EbkYJ2j3At4UkdWqeo6qrheRF3EedFIDfNOduQ/ObUyfBmJxJvAdcxJfY3Jzc0lMTCQrKyvYZEKKi4tJTExs7eq7pGO1iaqSn59Pbm4ugwYN6sCaGWNM+CirquHTnYdZvDWPj7flsWl/cV1ecmwkPROi6BkfxaC0eKZm9aBnfBTRRblNrLFthWr2/r9peG/y2rz7gfuDpK/AeUrZcauoqGg04JvWERF69uzJoUOHQl0VY4xpN6pKUXkNB4srOFhc6bwXVXKwuJK1ewpZtfsI1T4lKsLDCVmp3D17JKcMTWNkn0QivcFvjbNo0b4Oq3+4De93GAv4bc/a1BgTLvYVlvPQu1sRgZ7x0fSIj3J72c7njKRoeiZEN2tdB4sqeHvdft5cu481OQVU1vgbLBMb6WVIejw3njKIU4f2YmpWKjGR3rberePWbYN+qHm9XsaNG1f3/corr+SHP/xhCGtkjDFdQ3Z+KVc99il5JZUkxkRwuLSKYHPl+iTHMCEzhQn9U5jQP5lx/ZJJjHHmJAUG+uW7DqMKQ9MTuGb6QPokx9ArMZr0xBjSk6JJT4wmITqiU3R8LOiHSGxsLKtXrw51NYwxpkvZcqCYax7/lGqfn3/edhLjM1Pw+5XC8mpnIl1JFfmlVewtKGftnkLW5BTwn/X7ARCBob0SSIyJYFVOAaowLD2BO88Yxnnj+zA8o/PP9bKgb4wxpktYk1PA9U8tI8rr4YVbT6oL0h6PkBofRWp8FEPTG5Y7UlrFmtwC1uQUsia3gPySSr595jDOG9eHYV0g0Afq9kH/Z6+vZ8PeoqPSfD4fXm/rz8WM7pvET88f0+Qy5eXlTJw4se77PffcwxVXXNHqbRpjTFfk9yvbD5WQnhhDclzjlwN/uiOfOc+sICUukgU3TWdAz7hmbyM1PoqZI9KZOSLIEUEX0+2DfqjY8L4xxjSkqmw5UMKS7Xks2ZHPpzsPU1BWjdcjTBmYypkj0zlzVDpDeiXUnUN/f/NBbpu3kv494nhuzjR6J8eEeC/CV7cP+sF65HadvjHGtD9V5VBJJVv2l7D5QDErsw/z6Y7D5JdWAZCZGstZozI4YVAPdueX8b9NB/n125v49dubGNAjjjNGptMvJZbf/ncTwzMSefbGE5s9I7+76vZB3xhjTPPkHinjwy15fLDlIIdLq7jxS4OYPbZ3s2etbztYzLvZ1bz7ylq27C9hy8FiCsqq6/L7JscwY0Qvpg/uyUmDe9K/x9FD9HedM4K9BeW8t+kg7206yD+W7aayxs8JWak88fUTSIqxu4EeiwX9EKl/Tn/27Nk88MADIayRMcYcraLax7Kdh/lgyyE+2HKIbQdLAOiXEkukV7h9/meM7ZfE/509gpnDewUN/qrKB1sO8cTinXy0NQ+AxJi9DM9I5Mtj+zA8I4ERGYkMy0gkLSHqmAcQfVNiuWb6QK6ZPpDyKh8b9xcxuk9SWF4TH44s6IeIz+c79kLGGNMBaofZN+8vrnttOVDMpv3FVNb4iYrwMG1QD648oT8zR/RiSK8EfH7lldV7efDdLdzw1HJOyErlrrNHMG1wT8A5YHj5sz08+fFOth0sIT0xmrvOHk6fyhwumX16m1zTHhvlZfKA1ONeT3diQd8YY7qwT3fk870X11BcUU1MpJfoSA/REV5i3HcBduSVctg9jw6QlhDFiN6JXDt9IKcMS2PaoJ7ERh3dk47wCpdNyeSCCX15Yflu/vLeNq6Yu5RTh6Uxtl8yzy/bzZGyasb0TeKPX53AV8b3JSrCw6JFezrFTWy6Kgv6xhjTRb2+Zi//9+IaZ0Lc5Ewqqn1U1viprPFRWe2nssZPjd/P2aMzGJ6RyMjeiQzvnUhaCybDRUV4uPakLC6b0p95S3fx6KLtLN6Wx6xRGcw5ZRDTBvWwIB9GLOgbY0wXo6o89tEOfvXWJk7ISuWx66aSEhfVrtuMjfJyy2lDuGb6QEoqa0hPtMvmwpEFfWOM6UJ8fuVnr6/n2SXZnDe+D3+4fEKHTnKLi4ogLspCS7iyn4wxxoS5qho/728+yGtr9hId4eHUYWl8aWhag950eZWPO59fxTsbDnDLaYP54eyReDw2tG6+YEHfGGPCkKqyfm8RL63M5bU1ezlcWkVaQhQ+v/LyZ3sAGNk7kVOGpnHKsDSGpidwx4JVrMkt4GcXjOH6k7NCuwMmLFnQN8aYDqaqFFXUUO3zO68apar2s8/PpzsO89LKXDYfKCbK6+Gs0RlcOqUfpw3rhUeE9XuL+GjbIRZvzePZJdk8vngnANERHv52zRTOGdM7xHtowpUFfWOM6UBLtufzyzc3sL7eg77qm9g/hV9cNJbzx/dpMAlvXGYy4zKT+cbMoZRX+Vi26zCfZR9h1qgMxmUmt2f1TSdnQT9EHnvsMR555BEAPv/8c8aPHw/AGWecwR//+MdGy5WXlzN79mzee+899u7dy3XXXcf+/fvxeDzccsstfPvb3z5q+ZycnCaXKSgo4KabbmLdunWICE8++SQnnXQSDz30EI899hiqys0338ycOXOoqqpi1qxZvPfee0RE2K+OMS2xK6+UX721kYUbDtA3OYbvnzOCxJgIIr0e9yVEuZ8H9YpnSK+EZq03NsrLjOG9mDG8VzvvgekK7D93iNx8883cfPPN7Nmzh5NPPrnZT9x78sknueSSS/B6vURERPCHP/yByZMnU1xczJQpUzjrrLMYPXp03fLHWubb3/42s2fP5qWXXqKqqoqysjLWrVvHY489xrJly4iKimL27NnMmDGDSZMmceaZZ/LCCy9w9dVXt0u7GNPVFJZX8/B7W3n6k11Eej3cdfZwbjp1sN021oSEJ9QV6O7WrVvHuHHjmr38/PnzufDCCwHo06cPkydPBiAxMZFRo0axZ8+eo5ZvapmioiI+/PBD5syZA0BUVBQpKSls3LiR6dOnExcXR0REBDNmzOCNN94A4KKLLmL+/PnHt9PGdAM1Pj//213N6b9fxOOLd3LxpH4sumsmd5wxzAK+CRnr6S9oeDlLmzxU9ypt1mJr165l7NixzVq2qqqKHTt2kJWV1SBv165drFq1imnTpjVavv4yO3bsoFevXtxwww2sWbOGKVOm8NBDDzF27Fjuvfde8vPziY2N5a233qo7/TB27FiWL1/erPoa01moKnsLK+idFIP3OC9xU1X+s24/v1u4mR2Hqpg2qAc//spoxvazc+0m9KynH2It6enn5eWRkpLSIL2kpIRLL72UBx98kKSkpKBlgy1TU1PDZ599xu23386qVauIj4/ngQceYNSoUdx9992cddZZzJ49mwkTJtSdw/d6vURFRVFcXNzKPTYmfPj8yltr93HRIx/zpQfe44q/LyE7v7TV6/t4Wx4XPfIxt8//DI8I35oUzfO3TLeAb8KG9fSD9MiLi4tJTGyT/v4xrV27lu9+97t132tqavjBD36AiDBw4EDuvPPOurzY2FgqKiqOKl9dXc2ll17K1VdfzSWXXBJ0G40tk5mZSWZmZl3P/7LLLqt7vO+cOXPqhv1/9KMfkZaWVleusrKSmBi7xabpvCqqffzrs1we+3AHu/LLyOoZx+0zh/Dc0mxmP/gRPzpvFNdMG9Dse8avzS3kt//dxEdb8+ibHMPvLhvPJZMz+ejDD+y+8yasWNAPIb/fz9atWxk5cmRd2qOPPsqFF17IjBkzGiyfmpqKz+ejoqKCmJgYVJU5c+YwatQovve97wXdRlPL9O7dm/79+7N582ZGjBjB//73v7oJfgcPHiQ9PZ3du3fz8ssvs3DhQgDy8/Pp1asXkZGRbdUMxnSYwrJq5i3dxdOf7CKvpIrxmcn89erJnDOmN16PcO30gdz9r8/58SvrWLh+P7+9bDx9kmODrqv2WfMvLM/hzbX7SI2L5MdfGc3V0wbYOXsTtizoh9C2bdvIzMwkOvqLJ1rVDrc35uyzz2bx4sXMmjWLjz/+mHnz5jFu3DgmTpwIwK9+9SvOPfdczj33XB5//HF27NjR6DIAf/nLX7j66qupqqpi8ODBPPXUUwBceuml5OfnExkZySOPPEJqqvPM6vfff7+urDGdxc68Up7+eCf/XJlLWZWPGcN7cduMIUwffPQT4PqmxPLsjSfy3NJsfvXWJs7504f87MIxXDSxHyLC7vwyFm05yKLNh1iyPZ/yah/xUV7uPHMYN586iMQYOxg24c2CfggNHz6cDRs2HJV20UUXceutt9KjRw/uueceevTocVT+HXfcwR//+EdmzZrFKaecgmrwCYNvvfUWAH379m10GYCJEyeyYsWKBukfffTRUd9rz+EvWLCAX//618feOWNCTFVZsj2fJz/eyf82HSTS4+H8CX2Zc8ogRvcNPvcFQES49qQsTh3Wi//75xq++8IaXlyey4HiCnYccs73D+wZxxUn9GfGiF6cNLin9exNp2FBP8xceOGFdZfkBTNp0iROP/10fD4fXm/H/qOpqqrioosuYsSIER26XWNaoqLax2tr9vLk4p1s2l9Mz/govnXGMK6ZPqBFj3vNSovnxVtPYu6HO3jq452M7JPEtdMHMnNEOoPS4ttxD4xpPxb0O6Ebb7wxJNuNioriuuuuC8m2jWmOLQeKmfPMcnIOlzOydyK/vXQ8F0zs2+qeuNcj3D5zCLfPHNLGNTUmNCzoG2O6hA+3HOKb8z8jJsrLMzeeyGnD0mzmvDH1WNA3xnR685Zmc99r6xmWnsCTXz+BvinBZ9wb092F5OY8IvI7EdkkIp+LyL9FJCUg7x4R2SYim0XknID0KSKy1s37s9ghvDHdns+v/Pz1Dfz4lXXMGN6Ll24/2QK+MU0I1R353gHGqup4YAtwD4CIjAauBMYAs4G/ikjtybhHgVuAYe5r9vFUoKkZ7aZ1rE1NRyqtrOGWZ1fw5Mc7ueFLWTx23VQSom3w0pimhCToq+pCVa1xvy4FMt3PFwLPq2qlqu4EtgEnikgfIElVl6gTWZ4FLmrt9mNiYsjPz7cg1YZUlfz8fLtTn+kQewvKufxvS3h/80F+fuEYfnr+mOO+Z74x3UE4HBbfCLzgfu6HcxBQK9dNq3Y/108PSkRuwRkVICMjg0WLFtXPJz4+npycnKDlVdUmANXTnDbx+XyUlpaSnZ3dQbUKrZKSkga/W6E0030PVZ06oj12F/l4d3cNS/bW4BX4zuRoBlTuYtGiXe263dYKt9+RcGBt0lBHtkm7BX0ReRfoHSTrXlV91V3mXqAGqH1Wa7Cook2kB6Wqc4G5AFOnTtWZM2c2v+I4/zRbWqarszZpKOzaZIHzFqo6tVd71Pj8LNxwgKc/2cWynYeJifRw2dQB3HLa4LC/Xj7sfkfCgLVJQx3ZJu0W9FV1VlP5InI98BXgTP1inD0X6B+wWCaw103PDJJujOmiCsurmf9pNs8tyWZvYQWZqbH86NyRfHVqf1LiokJdPWM6pZAM74vIbOBuYIaqlgVkvQYsEJE/An1xJuwtU1WfiBSLyHTgU+A64C8dXW9jTMc4UlrFV/++hK0HS/jS0J7cd8EYzhyVYeftjTlOoTqn/zAQDbzjnideqqq3qep6EXkR2IAz7P9NVfW5ZW4HngZigbfdlzGmiymrquHGZ5aTfbiM5+ZM45RhaccuZIxplpAEfVUd2kTe/cD9QdJXAGPbs17GmNCqqvFz23OfsSangEevmWIB35g2Fg6z940xBr9fueufa/hwyyF+c+k4zhkTbB6wMeZ4NOs6fRH5l4icJyKhupmPMaaDLN2RT0llzbEXDPB5bgH3Li5j/d7CVm1TVfn5Gxt4bc1efjB7BFecMKBV6zHGNK25QfxR4Cpgq4g8ICIj27FOxpgQyTlcxpVzl/Lrtza2qNwj729jT4ly978+p8bnb/F2H3l/G09/soubThnE7TPsiXbGtJdmBX1VfVdVrwYmA7twJuB9IiI3iEhke1bQGNNx3tlwAICXVuaSX1LZrDK788tYuOEAg5M9rNtTxFMf72rRNud/ms3vF27hkkn9+NG5o+zGWMa0o2YP14tIT+DrwE3AKuAhnIOAd9qlZsaYDrdww37SE6OprPEzb2nz7qz4zJJdeEW4Y1I0s0Zl8Id3NrM7v+yY5QDeXruP//fKOs4Ymc5vLhuPxy7JM6ZdNfec/svAR0AccL6qXqCqL6jqt4CE9qygMaZjFJRVsXzXEb46tT9njkzn2SXZlFf5mixTXFHNC8tzOHdcH3rEePjFRWOI8Hi495W1x3y2xSfb8vj286uZPCCVR66aTKTXpgwZ096a+1f2sKqOVtVfq+q+wAxVndoO9TLGdLD3Nh3E51fOGp3BLacN5nBpFf/6LLfJMi+tzKWksoYbTxkEQJ/kWH4wewQfbc3j36v2NFpubW4hNz+7gkFp8Tx5/QnERnkbXdYY03aaG/RTROSSeq8zRSS9XWtnjOkwC9cfICMpmnH9kjlxUA8mZCbz+Ec78PmD99h9fuXpT3YxeUAKE/un1KVfM20gkwek8Is3NgSdF7Azr5SvP7WMlLgonrnxRJLjbFqQMR2luUF/DvA4cLX7egz4HvCxiFzbTnUzxnSQimofH249xFmjM/B4BBHhltOGsCu/rG5yX33vbTpIdn5ZXS+/lscjPHDpeEoqa/jlm0dfBXCgqIJrn/gUBebNOZHeyfYoZmM6UnODvh8YpaqXquqlwGigEpiGcw99Y0wn9sn2PMqqfJw1+osb4pwzJoP+PWKZ++H2oGWeXLyTvskxzA5yE53hGYncPmMI/161hw+2HAKcB+hc/+QyjpRW8fQNJzC4l00HMqajNTfoZ6lq4OH+QWC4qh7Geda9MaYTW7j+AAnREUwf3KMuLcLr4aZTBvPZ7gJWZh8+avkNe4tYsiOf607OIqKRCXjfOH0og3vFc++/13K4tIqbnlnO9kMl/P3aqYzPTAlaxhjTvpob9D8SkTdE5Hr3kbivAh+KSDxQ0H7VM8a0N79feXfjQWaO6EV0xNET6i6fmklybCR//2DHUelPfbyT2EgvX2viznkxkV4euGQ8uUfKOeuPH7Ai+wh/umKi3U/fmBBqbtD/JvAUMBGYBDyL8wS8UlU9vb0qZ4xpf6tyCsgrqeSs0RkN8uKiIrh2+kDe2XiAHYdKAMgrqeTV1Xu5bErmMSfhnTioB1dNG0B+aRU/v2AMXxnft132wRjTPMcM+iLiBd5R1X+p6ndV9Tuq+pIe6yJcY0yn8M6GA0R4hJkjgl+Mc/3JWUR6PDy+eCcA85fupsrn5+tfymrW+n92wRjeuvNUrj2pecsbY9rPMYO++zz7MhFJ7oD6GGM62MIN+5k+uCfJscF77b0So7lkcj/+tTKXfYXlzFuazekjejGkmRPxIr0eRvdNassqG2NaqbmP1q0A1orIO0BpbaKq3tkutTLGdIjth0rYcaiUr5+c1eRyN506mOeX53DDU8vJK6lscJmeMaZzaG7Qf9N9GWO6kNpr8GeNang+P9DQ9ARmjUrn3Y0HGZaewClDbTKeMZ1Rs4K+qj4jIrHAAFXd3M51MsZ0kIXr9zO2XxJ9U2KPuextM4bw7saD3HzqYHsSnjGdVHMfuHM+sBr4j/t9ooi81p4VM8a0r0PFlazKKeDs0Q1vrhPM1KwefPD9mVw+NbOda2aMaS/NvWTvPuBE3GvyVXU1YCf1jOnE/rfxAKoEvVSvMQN7xlsv32PilG4AACAASURBVJhOrLlBv0ZVC+ul2SV7xnRiCzccIDM1lpG9E0NdFWNMB2lu0F8nIlcBXhEZJiJ/AT5px3oZY9pRaWUNi7flcdboDOu5G9ONNDfofwsYg/OQnX8ARcB32qtSxpj29dHWQ1TV+Jt9Pt8Y0zU0d/Z+GXCv+zLGdHIL1x8gJS6SE7JSQ10VY0wHalbQF5HhwF1AVmAZVT2jfapljGkvpZU1vLvxALNGZTT6hDxjTNfU3Jvz/BP4G/A44Gu/6hhj2ts/lu2mqKKGa04aGOqqGGM6WHODfo2qPtquNTHGtLuqGj9PLN7JtEE9mDzAhvaN6W6aO7b3uoh8Q0T6iEiP2le71swY0+ZeXb2HfYUV3D5zSKirYowJgeb29K93378fkKbA4LatjjGmvfj9yt8+2M6oPknMGN4r1NUxxoRAc2fv2933jOnk3tl4gO2HSnnoyol2bb4x3VSTw/si8oOAz5fXy/tVe1XKGNM86/YUsmr3kWMup6o8umg7/XvEct64Ph1QM2NMODrWOf0rAz7fUy9vdms3KiK/EJHPRWS1iCwUkb4BefeIyDYR2Swi5wSkTxGRtW7en8W6KqYbyy+p5O6XPucrf1nMFX9fyrKdh5tc/tOdh1mdU8Atpw2xy/SM6caO9dcvjXwO9r0lfqeq41V1IvAG8BMAERmNc6AxBueg4q8i4nXLPArcAgxzX60+6DCms/L5lXlLdnH67xfxr89yuemUQWT2iOWWeSuaLPfoou2kJURx+RR7Qp4x3dmxzulrI5+DfW82VS0K+BofsK4LgedVtRLYKSLbgBNFZBeQpKpLAETkWeAi4O3W1sGYzmZl9hF+8uo61u8t4uQhPfnZBWMYlpHIdSdlcfFfP2603Pq9hXyw5RDfP2cEMZHeRpczxnR9otp47BYRH1CK06uPBcpqs4AYVY1s9YZF7geuAwqB01X1kIg8DCxV1efcZZ7ACey7gAdUdZabfipwt6p+pZF134IzKkBGRsaU559/vkV1KykpISEhoVX71VVZmzTUUW1SXKW8sLmKxXtqSI0WvjYqihMyvEdNxttW4OOmslkALMx4jyjvF3l/W1PB6oM+/jAzjvjI9jsrZr8jDVmbNGRt0lBbt8npp5++UlWnBstrsqevqq3uFojIu0Cwp3ncq6qvquq9wL0icg9wB/BTgp8y0CbSg1LVucBcgKlTp+rMmTNbVPdFixbR0jJdnbVJQx3RJuVVPi7/+yds3l/O7TOHcMfpQ4mPbvhnOxNggfP5lf1JPPy1yXg8wu78Mpb9931uOnUw5501ql3rar8jDVmbNGRt0lBHtklzr9NvsdpeeTMsAN7ECfq5QP+AvExgr5ueGSTdmC5LVfn+S2tYv7eIx6+bypmjMppV7q21+/lNj03c8+VRzP1oOxEeD3NOsatujTHNvyNfmxKRYQFfLwA2uZ9fA64UkWgRGYQzYW+Zqu4DikVkujtr/zrg1Q6ttDEd7OH3tvHG5/u4e/bIZgd8gGumD+DvH+zgz//byosrcrlkcj8ykmLasabGmM6i3Xr6x/CAiIwA/EA2cBuAqq4XkReBDUAN8E1VrX3Az+3A0zhzC97GJvGZLuw/6/bxh3e2cMmkftx6WstufHnf+WPYc6ScP76zBRG4pYXljTFdV0iCvqpe2kTe/cD9QdJXAGPbs17GhIMNe4v47gtrmNg/hV9dMq7Fd8+L8Hp4+KrJ3PDUcoakxzO4l02aMsY4QtXTN8YEkVdSyc3PriA5NpK5105p9SV28dERvHjbSTR1dY4xpvuxoG9MmKis8XHbvJXkl1byz1tPJr0NzsPbjSuNMYEs6BsTBlSVH7+yjhXZR3j4qkmMy0wOdZWMMV2Q3YTbmDbWmiH1Jxbv5MUVudx5xlC+Mr7vsQsYY0wrWNA3pg2tzS3khPvf5fllu5td5n8bD3D/Wxv58tjefGfW8HasnTGmu7Ogb0wb8fuVH7+6jrySKu7591peXb3nmGU27ivizn+sYmzfZP741Yl4PHYO3hjTfizoG9NGXvosl9U5Bdx/8VimDerB915cw8L1+xtd/lBxJTc9s4KEmAgeu24qsVH2MBxjTPuyoG9MGygsr+Y3b29iysBUvnbCAB6//gTG9UvmjgWr+GjroQbLV1T7uHXeCvJLK3n8uhPonWx3zDPGtD8L+sa0gQff3cLhsip+dsEYPB4hITqCZ244kcG94rn52RUs33W4bllV5e5/fc5nuwt48IqJNlPfGNNhLOgbc5w27S/i2SXZXHXiAMb2+yKAJ8dFMm/ONPomx3LjU8tZm1sIOPfUf3X1Xr5/zghmj+0TqmobY7ohC/rGHAdV5aevricxJoK7zh7RIL9XYjTP3TSNpNhIrn3yUx55f5tzT/3J/fjGzCEhqLExpjuzoG/McXjj8318uvMw3z9nBKnxUUGX6ZsSy4KbpxHl9fC7/25m6sBUft2Ke+obY8zxsjvyGdNKFTXK/W9uZGy/JK48YUCTyw7sGc+Cm6fx1Me7+O5Zw4mOsJn6xpiOZ0HfmFZ6fXs1+4uqeeTqSXibcX390PRE7r94XAfUzBhjgrPhfWNaYcehEv6zq5pLJ2cyZWCPUFfHGGOaxYK+MS1UXuXjp6+tJ9IDd3+54eQ9Y4wJVza8b0wzlVbW8NzSbB77aAd5JVVcMyqK9ES7qY4xpvOwoG/MMRRVVPPMx7t44uOdFJRVc+qwNO48cxiluz4PddWMMaZFLOibLq3G5+f9zYc4UlbFSYN70r9HXLPLHimt4qmPd/LUJ7sorqjhzJHp3HHGUCYNSAVg0a52qrQxxrQTC/qmSzpUXMnzy3azYNlu9hVW1KX3S4nlpCE9OXlIT04a0pM+ybEAFFdUs35vEev2FDqvvUVsP1SCKswe05s7zhh61N32jDGmM7Kgb7oMVWVF9hHmLcnm7XX7qPYppw5L474LxjAoLZ6lO/L5ZFs+7248wEsrcwHI6hmHiLAzr7RuPb2TYhjbL4nzxvXh3HF9GNE7MVS7ZIwxbcqCvunUCsurWZNTwOqcAt5au49N+4tJjIng2ulZXD19AEN6JdQtOzwjketOysLvVzbuL2LJ9nyW7jiMR+CSSf0Y2y+Zsf2S6ZUYHcI9MsaY9mNB33QKfr9SVFFN7pFyVuUUsHp3AatyjrDjkNNDF4Fx/ZJ54JJxXDCxL3FRjf9qezzCmL7JjOmbzE2nDu6oXTDGmJCzoG/CRkW1j/c2HeSjrYfIK6niSGkVR8qqOFJWTUFZFX79Ytm0hCgm9k/l0smZTOyfwvjMZBJjIkNXeWOM6QQs6JuQ8vmVT7bn8erqvfx33X6KK2tIjo2kT3IMqXFRjOydRGp8JKlxUaTGRZGeFM2EzBQyU2PtgTXGGNNCFvRNSKzJKeCV1Xt4fc0+8koqSYyOYPbY3lw4sR8nDenZrHvZG2OMaRkL+qbDqCrvbz7IX9/fzorsI0R5PZwxMp0LJ/bl9JHpxETak+eMMaY9WdA37a7G5+etdfv56/vb2LS/mH4psdx3/mgunpxJcqydhzfGmI5iQd+0SmF5Na+t3sMrq/ciwICecQzsEU9WWhwDesQxsGc88dFe/rVyD3//cDvZ+WUM6RXP7y+fwIUT+xLptWc9GWNMR7Ogb5rN71eW7sznxeU5vL1uP5U1fkb2TiQpNpJPtuXzctGeo5aP8Ag1fmVCZjL3XDOFs0dn4LFz9cYYEzIhDfoichfwO6CXqua5afcAcwAfcKeq/tdNnwI8DcQCbwHfVlUNtt7urNrnZ19BBUUV1RRX1FBSWUNJpfO59ntxRTUldZ+/SI/wChmJMWQkRZORFFP3SkuI4vXtVfx0+SKy88tIjIngq1P7c8UJ/Y+6NW1FtY+cw2Vk55eRfbiM/YXlzBiezpeG9rSZ9sYYEwZCFvRFpD9wFrA7IG00cCUwBugLvCsiw1XVBzwK3AIsxQn6s4G3O7re4URVyTlczurcgrq70q3bU0hljb/RMhEeITEmgsSYSBKiI0iIiaBvSgwJ0RFU+5QDRRWs3H2EA0WVVNVbz/TBiXxn1jBmj+lDbFTDSXcxkV6GZSQyLMNuW2uMMeEolD39PwE/AF4NSLsQeF5VK4GdIrINOFFEdgFJqroEQESeBS4iDIO+qlJZ46e8ykd5tfuq8lFW5aO0qoayytr3GkqrfJRV1VBYXk1BWXXde0F5FYVl1RRX1hAT4SU+OoL4aC/xUe57dAQ+v7J+bxGHS6sAiI7wMK5fMtdMH8iI3omkxEaSEBNBYrTznhAdQWJMBNERnmb1ulWVgrJqDhRXcKCokv1b13LFeSe1d/MZY4xpRyEJ+iJyAbBHVdfUC0D9cHrytXLdtGr3c/30DnXfa+vZsLeISp+fqho/VTU+qtzPlTV+Kqv9lFf7WrROr0dIjo0kJTaS5LhI0hKiGJqeQHKs0xOvrPFRWuWjtLKG0krn/UhpFT5VZo1KZ0L/FCZkpjCid2KbTo4TEVLjo0iNj2Jkb1i01ybeGWNMZ9duQV9E3gV6B8m6F/gRcHawYkHStIn0xrZ9C86pADIyMli0aNGxqnuUkpKSoGVycispKvET4RFiPJDghYhIiPQIER6I8niI9nqI8kKUV4jyQrTHeY+JEKK9EOMVoiO+eI8QAnre1e6r7IuNRgPxjdX0CJQfIW8r5G1t0S62WGNt0p2FW5vMdN9DVadwa49wYG3SkLVJQx3ZJu0W9FV1VrB0ERkHDAJqe/mZwGciciJOD75/wOKZwF43PTNIemPbngvMBZg6darOnDmzRXVftGgRwcq0cDVdSmNt0p2FXZsscN5CVaewa48wYG3SkLVJQx3ZJh0+Zquqa1U1XVWzVDULJ6BPVtX9wGvAlSISLSKDgGHAMlXdBxSLyHRxjhSu4+i5AMYYY4w5hrC6Tl9V14vIi8AGoAb4pjtzH+B2vrhk723CcBKfMcYYE85CHvTd3n7g9/uB+4MstwIY20HVMsYYY7ocm5JtjDHGdBPS1W9qJyKHgOwWFksD8tqhOp2ZtUlD1iZHs/ZoyNqkIWuThtq6TQaqaq9gGV0+6LeGiKxQ1amhrkc4sTZpyNrkaNYeDVmbNGRt0lBHtokN7xtjjDHdhAV9Y4wxppuwoB/c3FBXIAxZmzRkbXI0a4+GrE0asjZpqMPaxM7pG2OMMd2E9fSNMcaYbsKCvjHGGNNNWNA3xhhjugkL+sYYY0w3YUHfmFYQkbdF5Pp2WO96EZnZ1us14UVE9otImYg8Huq6NIeIfEVESkTELyKnhLo+pvUs6JuQE5FdIlIlImn10leLiIpIVmhqVleP+0TkucA0Vf2yqj5znOt9WkR+WW+9Y1R10fGstz2IyI/cf/olIlIhIr6A7+tDXb/jISK3ici7Idj02ap6k1uHGPd3PbO9NiYiZ4vIByJSLCKbguQPEZGP3IOR9SJyWm2eqr6hqgnAwfaqn+kYFvRNuNgJfK32i4iMw3mMsgkDqvorVU1w//HfBiyp/a6qY0Jdv8aISLs/SbQjttFGSnCuB7+nkfyXgI+AHsAvgVdEJKWD6mY6iAV9Ey7mAdcFfL8eeDZwARE5T0RWiUiRiOSIyH318q8TkWwRyReRH7sjCLPcvPtE5EURedbt6awXkakBZfuKyL9E5JCI7BSRO9302cCPgCvcXu0aN32RiNT20tYE9HpL3B7bTDfvn+5QbqGIfCgiY9z0W4CrgR+4ZV530wPrHC0iD4rIXvf1oIhEu3kzRSRXRP5PRA6KyD4RuSFYw4rIlSKyol7ad0XkNffzuSKywW2XPSJyVzN/ZvW3M1ZE3hORIyKyUUQuCsh7XkQeEpF3RKTUbb90EfmriBS4P49xAcvvF5EfiMgmETksInNr993Nv1hEPnfLfiQio+uVvcsdgShy037i/lyLRWSdiJznpk8CHgRmuj+H/W76UhG5JmCddaMBAb3y20VkO7DuWPsfDlT1E1WdD+yqnyci44HhwC9UtUJV/wFsB8JqH8zxs6BvwsVSIElERomIF7gCeK7eMqU4BwYpwHnA7bX/WN1/+n/FCaR9gGSgX73yFwDPu+VfAx52y3qA14E1bpkzge+IyDmq+h/gV8ALbq92Qv2Kq+qEgF7w94DNwGdu9tvAMCDdTZvvlpnrfv6tW/b8IG1yLzAdmAhMAE4E/l9Afu+A/ZwDPCIiqUHW8xowQkSGBaRdBSxwPz8B3KqqicBY4L0g62iSiCQB77jrSsP5OT0pIkMDFrsCuMvNj8D5mX8A9ATeAn5bb7VfA84ARgCTgO+725qO87O+wS07D6dXGtjjvgI4y80H52dyMk57/QZ4XkTSVHUV8B1gkftz6N2C3f4KMAWY1Mz9bxMicoN7sNPYK70Vqx0DbFHV8oC0NW666UIs6JtwUtvbPwvYBOwJzFTVRaq6VlX9qvo58A9ghpt9GfC6qi5W1SrgJ0D9200uVtW3VNXnbqs2gJ8A9FLVn6tqlaruAB4DrmxJ5cWZ4PRL4AJVLXLr/KSqFqtqJXAfMEFEkpu5yquBn6vqQVU9BPwMuDYgv9rNr1bVt3CGb0fUX4mqlgGv4p4+cYP/SJyDgdr1jBaRJFU9oqqf1V9HM1wMrFPV+arqU9XlOAdSlwYs809VXeMGlleBQlV9wf15vIgT2AM9pKp73X3/NV+c/rkVeFhVV7rbmgtE4wTgWn9yy5a7bfCCqu5zf3fm4fxuBS7fGveraoG7jebsf5tQ1adUNaWJV2vOuycAhfXSCoHE46+xCScW9E04mYfTA/069Yb2AURkmoi87w7BF+KcW66d/NcXyKld1g10+fVWsT/gcxkQ4/YOBwJ9A3tLOEP6Gc2tuIj0xwlc16vqFjfNKyIPiMh2ESnii2HVtEZWU19fIDvge7abVitfVWvq7VNCI+tawBdB8yrgFbeNwAlM5wLZ4kz0OqmZ9Qs0EDitXhteijPqUutAwOfyIN/r1z0n4HPgvg8EflRvW704emQnsCwiMifgdEABMJTm/xwaE7iN5ux/OCsBkuqlJQHFIaiLaUedZQKK6QZUNVtEduIEoDlBFlmAMyT/ZVWtEJEH+eIf9z4CerkiEssXQ7vHkgPsVNVhjeQ3+YAKd1uvAA+q6tsBWVcBFwKzcAJ+MnAEkOasF9iLE0xqZ8cPcNNaYyGQJiITcYL/d2sz3F7phSISCdyBc/DSv4XrzwEWNnKaorUC6xC47znAm6r6hybK1rWtiAwH/oJzqmCZqvrFmb3e1M+hFIgL+B5s2D+wXHvsf1AiMgd4qIlFBreit78eGC4iMapa4aZNwDmNYroQ6+mbcDMHOENVS4PkJQKH3YB/Ik5QrfUScL6InCwiUThD4RJkHcEsA4pE5G4RiXV76GNF5AQ3/wCQ5Z77D+ZJYJOq1j8nnQhU4ow4xOHMDQh0ABjcRL3+Afw/EeklzuWMP6HhPIdmcUcEXgJ+hzM7+x0AEYkSkatFJFlVq3EmvvlasYlXcM5tXyEike56p7sBt7XuFJE+7r7/EHjBTZ8LfEtEpoojQUQuEJG4RtaTAPiBQ4BHRG7D6enXOgD0dw96aq0GLnMn7Y3EGX1qSlvtf7S7zdpXg985VX0i4MqJYK+gAV9EPCISA0Q6XyWmdp/d02VbgR+LM4H0qzht9GoL62/CnAV9E1ZUdbuqrmgk+xvAz0WkGCcAvhhQbj3wLZyJevtwhiUP4gTdY23TB5yPM2FuJ5AHPI7TMwf4p/ueLyLBzndfCVwsR8/gPxXnFEU2zvnjDTgT1wI9gXMuvUBEXgmy3l8CK4DPgbU4EwF/GWS55lqAM+rwz3qnBa4FdrmnIG4DrglWuCmqegQ4B2dy3T6cXvkvcQJMaz0PvI8TjNbiTvRT1Y+BO4G/AwXAFpwDwKAjJ+4chb/htOU+YJD7udZ/cEZiDopIrpv2W5yR0EM4BxlNHmy14f5vwznVUfu6qunFW+Rsd50v48zUL8eZd1DrcuA0nDa9D7jY3S/ThdijdU2XJCIJOP+8hqnqzlDXx7SMOJfOXaaqi0Ndl/YgIruAVOB5Vb01xNU5JnEucZyPM2HyTFX9JMRVMq1k5/RNlyEi5wP/wxnW/z1O73BXKOtkTDCqmhXqOrSEqr6Jc6mr6eRseN90JRfiDKvuxbk2/kq1oSxjjKljw/vGGGNMN2E9fWOMMaab6PLn9NPS0jQrK6tFZUpLS4mPj2+fCnVS1iYNhV2bHF7pvPc43hvNtU7YtUcYsDZpyNqkobZuk5UrV+apaq9geV0+6GdlZbFiRWNXgAW3aNEiZs6c2T4V6qSsTRoKuzZZ4N6W4KqW/b63lbBrjzBgbdKQtUlDbd0mIpLdWF7YDO+LyJPiPC1sXSP5IiJ/FpFt7u00J3d0HY0xxpjOLGyCPvA0MLuJ/C/jzMgeBtwCPNoBdTLGGGO6jLAJ+qr6IXC4iUUuBJ5Vx1IgRUQ6y8MsjDHGmJDrTOf0+3H0U61y3bR9HVWB+Z9mk3ukHK8IHo/gEQI+CxEewRvwCvwe4fXUfY/0Cl6P8z3CI0RGeIjyeoiO8BBV+/J6iIzw4PMp1X4/1T6lxue++/0A9IiLokd8FBHesDl2M8aYBqqrq8nNzaWiooLk5GQ2btwY6iqFlda2SUxMDJmZmURGNv9uz50p6Ad7eErQmwyIyC04pwDIyMhg0aJFLdpQSUlJ0DLPLStnyxE/CvjD5PYGAiREQlK0kBQlJEcLCZHiHlBApAe8Hoj0CF6BCA94pPZVe+DifPcrVPuh2qdU+aHar1T7nLRIrWJ93v/ITPSQHN3c59h0bY39noTKTPc9VHUKt/YIB9YmjoSEBDIyMujXrx9+vx+v1xvqKoUVn8/X4jZRVQoLC1mzZg0lJSXNLteZgn4uRz9qM5NGHjOqqnNxHpLB1KlTtaWzIhubSVk/SVXx+RW/gt/9XOOvffc7774v0nx+pdrnr1uuxuenxq9U+fxU1QS83O/VPr87MuAh0itEeDxEeIUorwcF8kuryCuuJK+k9lXFvpJKjhypotrnrMfXBkcnER6hxi+w03niZlpCFCN7JzGidyIjeifiEaGgrIqCsmqOuO8F5VUUV9SQ1TOeSQNSmDQgldF9koiK6DqjEmE3C3mB8xaqOoVde4QBaxPHxo0byczMREQoLi4mMTEx1FUKK61tk8TEREpKSpg6dWqzy3SmoP8acIeIPA9MAwpVtcOG9oMRESK84d3rrT3QCDyQ8PkVvx9q/H73YMX57BEhJtJLTKSH6Igv3r0e4bWF75M2eBwb9xezaV8Rmw8U89zSbCpr/HXb8gikxEWREhdJalwUybGRLN91mNfWOMdmUREexvZNYtKAVCYNSGF8vxT694hFJLzb0Bhz/OzvvO21pk3DJuiLyD9wRijT3Mdb/hT3sZSq+jfgLeBcnEdPluE8wtIcgzOnwEtM5PENpyVFCScPTePkoWl1aT6/knO4DHGDfWJ0BB5Pw1/CfYXlrN5dwKqcAlbtPsJzS7N5YrHz4LukmAjG9ktmXL/kuveBPePsH4QxxrSDsAn6qvq1Y+Qr8M0Oqo5pBq9HyEo79l2k+iTH0mdcLF8e51xsUe3zs2lfMWv3FLJ2TyHr9hTy1Me7qPI5owYxkR5iIr1EeDxEed1JkF4h0uMhOtJDUkwkSbER7nskSTERJMVGkpEUw8T+KWQkxbTrfhtjTGcVNkHfdB+RXg/jMpMZl5lcl1ZV42fLgWLW7Slk28ESqnxfXLFQO++hxuenssZPcUUNB4oqKKqopqi8hvJq31Hr75scU3cKYdKAFMb0TSYm0kthWTU5R8rIOVxG7pFyco447wDpidGkJ0bTKymm7nNGUgy9k2KCjl4YY0xnZEHfhIWoCA9j3SH+lqqs8VFcUcPuw2Ws3l3AZ7uPsGp3AW+udaZ8RHqduQrFFTVHlUuKiSAzNQ4RWLunkLySSuo/dDIpJoJJA1KZMjCVyQNSmTgghYRo+7MxpjOaOXMmTz/9NC15HktryoQz++9lOr3oCC/RCV7SEqKZPCCVGxkEwMGiCnceQQFlVTX0T42jf49YMlPj6N8jjuTYo69trfH5yS+t4mBRJQeLK9hfVMG6PUV8ln2EP727BVVnsuKI3klMHZjK5Bh/sOoYY0zYsqBvuqz0pBjOGdObc8b0btbyEV4PGUkx7pyAo0ccCsurWZ1TwMrsI6zafYQXVuSwPAkuPFtt+N+YFvjNwu1szStv03WO7pvET88f06Iyl19+ORkZGaxevZqcnBzmz5/P3LlzWbp0KaeeeipPPPFEm9YxXHSdi6aNaUfJsZHMGN6L7501nHlzpvHzC8aw6bCffyzf3aL1tMV9E4wxx2/t2rUMHjyYxYsXc/311zNnzhx+85vfsG7dOl5++WUqKytDXcV2YT19Y1rhihP68+wH6/n1W5s4fUQ6fVNij1nm+WW7uf/Njfz5a5M4fWR6B9TSmPBz99lDQn5znoqKCgoKCvjOd74DQGxsLHPmzKFPH+cKo7i4OKKiokJZxXZjPX1jWkFEuGFMND6/cu+/16L1ZwDWs3hrHve+so6KGh/fmP8Za3IKOqimxpj61q9fz+TJk/F4nBC4Zs0apk2bBkBubi59+/btsvcKsaBvTCv1ivPw/XNG8P7mQ7yyek+jy207WMzt81cyLD2Bhd+dQc+EKG58ejm78ko7sLbGmFpr165lwoQJdd8///xzxo8fDzgHALWfuyIL+sYch+tPzmLygBR+9voGDhU3PAeYX1LJDU8vJzrCyxNfP4FBafE8c+OJ+FW5/qll5JV0zfOGxoSztWvXMnHiRMAZ6i8vLyc1NRU4+gCgK7Kgb8xx8HqE3142nrJKHz99bd1ReRXV/5+9O4+PqjofP/55MslkX0kIIWEJEJYAsgpWQIIr0ioIWK1btajVX61U22r1237b2kW7uKCl+lVrqxbcl7pQNzSyyA5hFxISlhAgC5B9z/n9cSdxyDpJJpkJed6v4S/RfQAAIABJREFUV17J3HvuvWdOJnnuOfcstdz+8hZyiyp5/vuTiXc89x8aE8I/bj6XE0UVLPrXJsqqapo7tVKqizz66KNcc801gLU8bVZWVsO+Bx54gMWLF3sqa11Og75SnTSsbyiLL05ixc7jfLTLmhDIGMN9b+5gy6FTPH7NeMYPiDjjmIkDI3nqexPZebSQHy3bSk2tjvlXSnU9DfpKucHtFwwhOS6MX767m9NlVTzxWTrvbc/h55eNYI5jzYHGLkmO5XfzxvDFvjx++e6uNjsDKqU65+abbyYiIqLthJ08xpvpkD2l3MDP5sOfF57D3KVruf75DezOKWLhpAT+X8rQVo+7fuogjhdW8NTnGfQNC+Cei5PO2l7DSnnazTff3C3HeDOt6SvlJmPiw7lj5hB25xQxNTGKP1411qUAfu8lw1k4KYEnV6Zz7+vbKal0zzP+jNxijp5278xnSqmeTWv6SrnR3RclkRAZxJwxcdh9XbunFhH+tOAcBkQGsWTlfrYdPsXfrpvYocWH6u06WsjCZ74iqW8o7/94eofPo5Q6u3hNTV9EZovIPhHJEJFfNLM/XETeF5HtIrJbRG7xRD6Vao2/r43vTRlIeJBf24md2HyExRcn8cpt51FRXcdVf1/LC2uyOvScP6+4kttf2kx1rWHn0UK+Pl7U7nMopc5OXhH0RcQGLAUuB5KB74lIcqNkPwL2GGPGASnAoyJyds6TqHqtqUP68N/FM5g5PIaHPtjDbS9t5mRpVbvO8cOXN3OyrIp/3XIufjbhrS3ZXZRbpVRP4xVBH5gCZBhjMo0xVcCrwNxGaQwQKtZD0hDgJKADnNVZJzLYznM3TebXVySzan8+c5asZn1mgcvHbz18mkevHs+MpBguHNmXd7blUK1DApVSgHjDMCERWQjMNsbc6nh9IzDVGHOXU5pQ4D1gJBAKXGOM+bCF890O3A4QGxs76dVXX21XfkpKSggJCenIWzlraZk01R1lcrCwlqe3V3KizDAzwZerh9sJsTffOTAlZxYAi0s/5qokqxFsW24NS7ZW8pOJ/ozv27VdePQz0pSWiSU8PJxhw4YBUFtbi81m83COvEtnyiQjI4PCwsIzts2aNWuLMWZyswcYYzz+BVwNPO/0+kbgqUZpFgKPAwIMA7KAsLbOPWnSJNNeX3zxRbuPOdtpmTTVXWVSUlFt/vDhHjPkgQ/NhIc+Ma9vOmzq6urOSPP51yeMWYYxyzC1td/sq6qpNRMf+sTc+e/NXZ5P/Yw0pWVi2bNnT8PPRUVFHsnDs88+a8aNG2fGjRtnRKTh53vuuafV48rKyswFF1xgampqzOHDh01KSooZOXKkSU5ONk888UST9G2lOXXqlFmwYIEZMWKEGTlypPnqq69MUVGReeKJJ8zo0aNNcnKyefzxx40xxlRWVpoZM2aY6urqFvPnXLb1gM2mhZjoLc372cAAp9cJQE6jNLcAbzveUwZW0B/ZTflTymOC/X15cM4oPvjxdBKjg/n5mzu45v/Ws/9EMQAZuSXcvXxbQ3ofn29aAvxsPswdH89ne3I51c6+AUqdTW677TbS0tL48MMPGTBgAGlpaaSlpfHYY4+1etwLL7zA/Pnzsdls+Pr68uijj7J3717Wr1/P0qVL2bNnzxnp20qzePFiZs+ezddff8327dsZNWoUe/bs4bnnnmPjxo1s376dDz74gPT0dOx2OxdddBGvvfaa28rBW4L+JiBJRBIdnfOuxWrKd3YYuAhARGKBEUBmt+ZSKQ8aFRfGGz/8Fn9aMJb9ucXMWbKaP67Yy20vbW51eODCSQlU1dbx/o7G99FK9T67du1i7NixLqdftmwZc+daXczi4uKYOHEiAKGhoYwaNYqjR89cYbO1NEVFRaxatYpFixYBYLfbiYiIYN++fZx33nkEBQXh6+vLzJkzeeeddwCYN28ey5Yt69ybduIVQd8YUwPcBXwM7AVeN8bsFpE7ROQOR7LfAeeLyE5gJXC/MSbfMzlWyjN8fIRrzh3I5z9NYf7EeJ5dlUn2qTL+78ZJLR6T3D+M5Lgw3tRe/Eqxc+dOxowZ41LaqqoqMjMzGTx4cJN9Bw8eZNu2bUydOrXF4xunyczMJCYmhltuuYUJEyZw6623UlpaSnJyMqtWraKgoICysjJWrFjBkSNHABgzZgybNm1q/xttgddMzmOMWQGsaLTtGaefc4BLuztfSnmjqGA7f144ju9NGUhNnWHy4Cj4quX0CyYl8LsP9rD/RDHDY0O7L6NKNRL6fljXnPg61zql79q1i0suucSltPn5+c3Ou19SUsKCBQt44oknCAtr/v00l6ampoatW7fy1FNPMXXqVBYvXswjjzzCfffdx/33388ll1xCSEgI48aNw9fXCs82mw273U5xcTGhoZ3/2/WKmr5SqmMmDIzk3MFRbaabO74/vj46Zl+pxjX9mpoa7r33Xn7605/y5JNPnpE2MDCQioqKM7ZVV1ezYMECrr/+eubPn9/sNVpKk5CQQEJCQkPNf+HChWzduhWARYsWsXXrVlatWkVUVBRJSUkNx1VWVhIQENC5N+7gNTV9pVTXiQ7xZ9bIvry97Sg/v2wEvja931eeUXxFkVtqrB1RV1dHeno6I0d+0wf86aefZu7cucycObNJ+sjISGpra6moqCAgIABjDIsWLWLUqFHce++9zV6jtTT9+vVjwIAB7Nu3jxEjRrBy5UqSk6156HJzc+nbty+HDx/m7bffZt26dQAUFBQQExODn1/7Zvlsif7lK9VLLJyUQF5xJavTtSuM6p0yMjJISEjA39+/YdvWrVuZNm1ai8dceumlrFmzBoC1a9fy8ssv8/nnnzN+/HjGjx/PihXWU+k5c+aQk5PTahqAp556iuuvv55zzjmHtLQ0HnzwQQAWLFhAcnIyV1xxBUuXLiUyMhKAL774gjlz5ritDLSmr1QvMWtEXyKD/HhzSzazRvb1dHaU6nbDhw9vMsRu3rx5/PCHPyQqKooHHniAqKgzH5fdddddPPbYY1x88cVMnz69xfUw6gN7//79W10zY/z48WzevPmMbcXFxaxevbrZ9MuXL+fhhx9u8725SoO+Ur2E3dcas798w2EKy6rbvSiQUmejuXPnNgzJa86ECROYNWuWR2YSrKqqYt68eYwYMcJt59TmfaV6kfox++/pmH2lXPaDH/zAI1MH2+12brrpJreeU4O+Ur3I6P5hjOwXqmP2leqlNOgr1YuICAsnJbD9yGkycos9nR2lVDfToK9ULzNvQjw2H+G9NG3iV6q30aCvVC8THeLPwKggDuSXejorqhdprUe76piOlKkGfaV6odgwf04UVrSdUCk3CAgIoKCgQAO/GxljKCgoaPdMfTpkT6leqF9YAJsPnfJ0NlQvkZCQQHZ2Nnl5eQ2z26lvdLRMAgICSEhIaNcxGvSV6oX6hQeSW3ScujqDj494OjvqLOfn50diYiIAqampTJgwwcM58i7dWSbavK9UL9QvzJ+q2jpOllV5OitKqW7kNUFfRGaLyD4RyRCRX7SQJkVE0kRkt4h82d15VOps0S/cako83k3P9fVZrlLewSuCvojYgKXA5UAy8D0RSW6UJgL4O3ClMWY0cHW3Z1Sps0RsmBX0TxR1fdD/ePdxxv32E4orqrv8Wkqp1nlF0AemABnGmExjTBXwKtB4MuTrgLeNMYcBjDG53ZxHpc4aceGBABzvhqC/MeskRRU1HCoo6/JrKaVa5y1BPx444vQ627HN2XAgUkRSRWSLiLh3QmKlepHoEDs+0j3N+wfySgDIPlXe5ddSSrXOW3rvN9d9uPFDQF9gEnAREAisE5H1xpj9TU4mcjtwO0BsbCypqantykxJSUm7jznbaZk05W1lkuL47mqewuxC2r6DpNqPueX6LZXH7sNWDf/LzTsJyP/aLdfqKbztM+INtEya6s4y8Zagnw0McHqdADSeIzQbyDfGlAKlIrIKGAc0CfrGmGeBZwEmT55sUlJS2pWZ1NRU2nvM2U7LpCmvK5Pl1jdX8zRo1xoI9CMlZapbLt9ceVRU15L/8UcABPTpT0rKaLdcq6fwus+IF9Ayaao7y8Rbmvc3AUkikigiduBa4L1Gaf4DzBARXxEJAqYCe7s5n0qdNfqFB3R5R76s/FLqO+4f1eZ9pTzOK2r6xpgaEbkL+BiwAS8YY3aLyB2O/c8YY/aKyEfADqAOeN4Ys8tzuVaqZ+sXFsC6AwVdeo365/lx4QEcPa1BXylP63TQF5GBLiY9bYwpammnMWYFsKLRtmcavf4L8Jd2Z1Ip1URseABFFTWUVdUQZO+a+/8DuaWIwPRh0Xy690SXXEMp5Tp3/KW/iNXprrW5PA3wL+AlN1xPKeUG/cK+maBnSExIl1wjI6+EhMhAhsSEcHpLNiWVNYT4e0UDo1K9Uqf/+owxs9yREaVU92oI+kVdF/QP5JYwNCaE+EhrXoCjp8oZ0S+0S66llGqbWzvyiYifO8+nlOo69VPxdlVnvro6Q2a+FfQT6oP+aZ2gRylPcls7m4g8D8wXkVKs4XY7gB3GmKfcdQ2llPvUB/1jXTRBT05hORXVdVbQj/impq+U8hx3PlybAcQaY6pFJB5rDP05bjy/UsqNguy+hAb4cqKLgv6BvFIAhsYEEx3ij93mo7PyKeVh7gz664FIINcYcxQ4SqPe+Eop79IvLKDL5t8/kGsN1xvaNwQfH6F/RADZOmxPKY9y5zP9Z4EvReRnIjJDRMLdeG6lVBfoFx7A8aLKLjn3gbwSwgP96BNsByAhMkib95XyMHcG/X8Dr2O1Hvw/4CsROeDG8yul3KxfWADHC7smEB/IK2FY3xBErNG88RGBOkGPUh7mzub9bGPMr503iIi/G8+vlHKzfuEB5BVXUlNbh6/NvbNyH8grZdaImIbX8ZGB5BVXUlFdS4Cfza3XUkq5xp1/5Wkisth5gzGma9oNlVJuERsWQJ2B/JIqt563sLyavOJKhjqN/4939ODP0dq+Uh7jzqAfC9whIjki8oGI/EFErnbj+ZVSbhYX/s0EPe5UP+e+c9D/Zqy+Bn2lPMVtzfvGmO9CQ5P+aGAsMAV4w13XUEq5V2zDVLzlMCDCbed17rlfz3lWPqWUZ3TVgjv5wBfAF077W11wRynV/eon6Dnu5rH6B/JK8bMJAxyBHqxOgzYf0bH6SnlQVy+4U79dF9xRygtFBdnxs4nbh+0dyCthcJ/gMzoH+tp86BemS+wq5Ules+COiMwGlgA24HljzCMtpDsXayKga4wxb7rj2kr1Vj4+Qt/QALfPv38gr4ThfZsurBMfGajN+0p5kHvH6HSQiNiApcDlQDLwPRFJbiHdn4CPuzeHSp294sIDOObGsfrVtXUcLihjaN/gJvsSdKy+Uh7lFUEfq8NfhjEm0xhTBbwKzG0m3Y+Bt4Dc7sycUmez2PAATrixef9QQRk1deaMnvv14iMDOVZYTnVtnduup5RynbcE/XjgiNPrbMe2Bo5FfK4CnunGfCl11rNm5avAGOOW89UP1xvWt5mgHxFInXF/x0GllGvcOSNfZ7TUCdDZE8D9xpja+mk9WzyZyO3A7QCxsbGkpqa2KzMlJSXtPuZsp2XSlLeVSYrje7s/73nVlFfXsuKzVIL9Wv/bavU8jvL49IA10c/Rr7dxMuPM8xXk1wLwYeo6Rkad/bPyedtnxBtomTTVnWXiLUE/Gxjg9DoByGmUZjLwqiPgRwNzRKTGGPNu45MZY57FWgCIyZMnm5SUlHZlJjU1lfYec7bTMmnK68pkufWtvXkqjszhtX3bSDpnMsNjm3a+c1V9ebyXm0a/sAIuv7hpH99B+aX8ZXMqMYNGkDIpocPX6im87jPiBbRMmurOMvGW5v1NQJKIJIqIHbgWeM85gTEm0Rgz2BgzGHgT+H/NBXylVPvUj9U/5qYm9wN5pc124oNvZgDUsfpKeYZXBH1jTA1wF1av/L3A68aY3SJyh4jc4dncKXV26+eYle+EG4K+MYbM3JJmO/EBBPjZiAn15+jpsk5fSynVft7SvI8xZgWwotG2ZjvtGWNu7o48KdUb9A2zFsN0x/z7ecWVFFfWtBj0wZqDvycN2zPGsD7zJJMGRWL39Yp6klIdpp9gpXo5f18bfYLtbgn6Gc0stNNYfETPmqDny/15fO+59Sx85iuy8ks9nR2lOkWDvlKKWMewvc46kGcFxZae6YM1Vj/ndAV1de4ZItjVth8pRMSaf+DbT67mjc1H3Da8UanupkFfKUW/cDcF/dwSgu22hn4CzUmICKSqto68EvfO999VducUktgnmP8unsHY+HB+/uYOfvzKNgrLqz2dNaXaTYO+UorYMPfMv38gr4ShfUNobS6N+iV2e0oP/j3HikjuH0b/iECW33YeP79sBP/ddZw5S1az+eBJT2dPqXbRoK+UIi48gILSKiprajt1ngOt9NyvlxAZBNAjOvMVllWTfaqc5P5hANh8hB/NGsabd3wLm4/w3f9bx1Mr0z2cS6Vcp0FfKdXQHJ/biTn4K2oMOYUVDI1p+Xk+WB35ALJPef+wvd3HCgEY3T/8jO0TBkby4d3TuXxsHI9+up/dOYWeyJ5S7aZBXylFrGPSnM704D9eai2i01ZNP9jfl4ggvx7Rg39PThEAox01fWehAX7873esxUBXp+d3a76U6igN+kqphpp+ZzrzHSu1erQPbWahncbie8gSu3tyiogN8yc6xL/Z/bFhAYyIDWV1el4350ypjtGgr5T6Zla+TtT0j5XW4SMwqE9Qm2kTInvGWP3dOUUkxzWt5TubnhTNpoOnKK/qXH8IpbqDBn2lFGGBvgT62To1//6x0joGRgXh79v26nnxEUFknyrv9vHuu44WUlpZ41LaiupaMvJKmjzPb2xGUjRVNXVs1J78qgfQoK+UQkSssfqdqemX1LX5PL9efGQg5dW1nCrrvrHup0qrmLd0LX9PzXAp/f4TxdTWmWaf5zubmtgHu82HNdrEr3oADfpKKQBiw/w7vOhObZ3heJlx6Xk+fNODvzub+DcePElNnXG5091uRye+5DaCfqDdxuTBkdqZT/UIGvSVUoD1XL+jNf2jp8qpqaPN4Xr1EhwT9HTnansbs6zm951HCzldVtVm+j05RYT6+zIgsu0+CjOSYvj6eDG5bpjgSKmupEFfKQVAv/BAThR1bE78+ufZSbGhLqVP8MCsfBuyCogI8sMYWJ9Z0Gb63TmFjOofho9Py7ML1puRFA3Amgyt7SvvpkFfKQVAvzB/qmsNJ12oBTurqa3j719kMCDUh/EJES4dEx7oR7Dd1m1Bv6iimj05RVw/dSDBdhtrM1oP+rV1hr3HitvsuV8vOS6MPsF2beJXXs9rgr6IzBaRfSKSISK/aGb/9SKyw/H1lYiM80Q+lTpb9Qvv2Fj997bnkJlfytyhfi7VisHqOBgf2X1j9bccPEWdgWlDo5mSGMXaNmrkBwtKKa+ubbMTXz0fH2HasGhWp+frCnzKq3lF0BcRG7AUuBxIBr4nIsmNkmUBM40x5wC/A57t3lwqdXaL7cBY/ZraOp5cmU5yXBgTY9sequcsITKo2zryrc8qwM8mTBgYybRh0WTml3KssOVr726Yia/14XrOZiRFk19SydfHizudX6W6ilcEfWAKkGGMyTTGVAGvAnOdExhjvjLGnHK8XA8kdHMelTqrxYVbz9nb05nvnW1HOVhQxk8uTsKnlZX1mtOds/JtyDzJuIQIAu02pg2znr+31sS/O6cQu82HYS6ORgCrMx+gs/Mpr+br6Qw4xANHnF5nA1NbSb8I+G9LO0XkduB2gNjYWFJTU9uVmZKSknYfc7bTMmnK28okxfG9o3mqrTMIsH7HPuLLs9pMX1Nn+PPqcgaF+eCXu5eS0tJ2XbvyZBWF5dX897MvCPRt3w1De1TUGHZml3F5oh+pqanUGUOoHd5au4vo4ubH7K/dVUFcMHy1ZlW7rtU/RPjPhv0Mr7P+nXnbZ8QbaJk01Z1l4i1Bv7m/+GYfjInILKygP72lkxljnsXR/D958mSTkpLSrsykpqbS3mPOdlomTXldmSy3vnUmT33Xf0ZARAwpKW13mXl90xHyynfw/NWTmJUc2+7yKI7M4fX920gcM4mR/Vx7dt4Rq9PzqDUbuTplAjOHW7Xxmce2sjHrJDNnzkQatVAYY/jp6s+4aFRfl8rB2eziPSzbcIjzps0gwM/mfZ8RL6Bl0lR3lom3NO9nAwOcXicAOY0Ticg5wPPAXGNM22NulFLt4upY/eraOp78PJ1zEsK5aFTfDl2rYax+Fz/X35B5EpuPMGlQZMO26cOiyS2u5EBeSZP0J4oqKSitatfz/HozkqKprKljk07Jq7yUtwT9TUCSiCSKiB24FnjPOYGIDATeBm40xuz3QB6VOuvFhgW41JHvzS3ZZJ8q556LhzepKbsqvmGCnq4N+huzTjKmfxgh/t80bNY/11/TzBC73TmFQNsz8TVn6pAo/GyiQ/eU1/KKoG+MqQHuAj4G9gKvG2N2i8gdInKHI9n/An2Av4tImohs9lB2lTprxYUHtLnoTlVNHX/7PIPxAyJIGRHT4WtFB/tj9/Xp0rH6FdW1pB05zdQhfc7YPiAqiAFRgaw90LTBcE9OESIwysUx+s6C7L5MHhSlQV95LW95po8xZgWwotG2Z5x+vhW4tbvzpVRvEhseQHFFDWVVNQTZm//38MaWIxw9Xc4frhrT4Vo+WGPb4yO6dondbYdPU1Vbx9TEqCb7pg+L5oPtx6iprcPX9k39Z3dOEYP7BJ/RMtAe05Oi+cvH+8grruxwvlvy7raj9AsP4LxGNzFKucoravpKKe/QL6z1CXoqa2pZ+nkGEwZGNHSK64yEyECyu7B5f2PWSURg8uCmQf/8odEUV9aw82jhGdt3Hyt0eSa+5lzgGLrX1gRA7bU6PY+fvJbGDc9vYMXOY249t+o9NOgrpRo0BP0Wnuu/vukIOYUVnXqW7yw+IpAjJ8uoqa3r9LmasyGrgFH9wggP9Guy7/yhVm35K6cm/sLyao6cLO/Q8/x6o/uHERnkxyo3jtc/XVbFz97YztCYYCYMjOCu5Vt5a0u2286veg+vad5XSnle/VS8728/xoHcEipr6qyv6loqa+p4e9tRJg2KbFhgprNmDo/h1U1HeHJlOvdeOsIt56xXVVPH1sOnuPbcgc3u7xPiz6i4MNak5/OjWcMA2Husfia+jgf9+il516Tnc0VM+2YpbMmv/rObgpIqnr/pXIb2Deb2l7bw0ze2U15dyw3nDXLLNVTvoEFfKdWgf0QgwXYbr2w83GSfv68PkUF2Hpwzyi21fIDLx8Zx9aQEnvoig3MToxpmtXOHnUdPU1Fdx3lDmjbt15s2tA8vrT9ERXUtAX62Dk2/25wLkmL4YMcxjpYEduo8AP9JO8r723P46SXDGZtg5ev570/mR8u28st3d1FRXcutM4Z0+jqqd9Cgr5RqEOBnY+0vLqS0qha7zQd/Px/8fX2w23zcFugbe2juGLZnn+Ynr6axYvGMhjUAOmt9pjVW/txmnufXmzYsmufXZLH54CmmJ0WzJ6eImFB/YkL9O3Xt6Y6WkF35tU32GWM4WVqF3deH0ICmjx2cHSss51fv7mLCwAjuTBnasD3Az8bTN0zintfS+P2HeymrquXHFw7rst+Rq46eLuftLdmIwO0XDMXuq0+QvY0GfaXUGSKC7EQEdd/1Au02ll43kSv/tpa7X9nGslunntGbvqM2ZJ0kqW8IfUJaDuBTEqPw9RHWZOQzPSma3TmFnWrar9c/IpChMcFszS3ngx05ZOWVkpVfyoH8UrLySiiqqCHYbuOBOaO4bsrAZlcnrKsz3PfmDqprDY9/d3yTMrH7+rDk2vH4+/nw2Kf7Kauq5c6ZQ8krqSS/pJKCkiryHT8XV9Rw64xEEiLd/4utqK7lo13HeXNLNmsP5FO/yODq9HyeuWESkcF2t19TdZwGfaWUxyXFhvL7eWP46Rvb3fJ8v6a2ji0HT3LVxPhW0wX7+zJhYARfHcinsqaWjNySDs8w2NjM4X15YW0Wdy3fBkD/8AASY4K5cnx/BvcJJnVfHr98dxfvb8/hTwvOYXB08BnHv7TuIKvT8/nDVWOa7Kvna/PhrwvHEehn45kvD/DMlweapKm/nziQV8JLP5jiltYAYwzbjpzmjc3ZfLA9h+LKGhIiA1l8URILJiaw5dAp7ntzB1f9fS3/uPlchsa4vnCR6loa9JVSXmHBpATWZxa45fn+7pwiSqtqmZrY9nj284dG8+Tn6Ww+eIqaOkNyXOee59dbfFESwWU5XD5jCoOjg5rMe7BoeiKvbz7C7z/Yy+wlq/jZpSO4ZVoiNh8hI7eEh//7NbNGxHDdlOY7Itbz8RF+P28MEwdGcqqsiphQf6JD/OkTYic6xJ/IIDsvrzvIb97fw4qdx/n2OXGdel+llTX8/M3trNh5nAA/H+aMiWPh5ATOS+zT0GJRP/nR7S9t4aqla3n6hkkNsyAqz9IHLkopr/HQ3DEk9Q3hJ6+muTQdcEs2ZFnD8JqblKex6UnRGAP/WGOtLOiO5n2A8CA/JsX6ktw/rNmJjkSEa84dyKf3zmTa0Gh+/+FeFjz9FXtyirjntTSC7Db+tOAcl2rmIsKCSQncOmMIc8fHM21YNCP7hREd4o/NR7jhvEEkx4Xxuw/2UFJZ0+H3dKiglPl//4qPdh3nZ5cOZ9P/XMxj14zn/KHRTR5RTBoUxbs/mka/8AC+/8LGZjuHeitjml3v7aygQV8p5TXqn++XVdVy9yvbOjx+f2PWSRKjg+nrQqfAcQkRBNltfP51LiH+vgyM6sYODVjDJJ///mSWXDueQwWlzHlyNTuPFvLHq8a6lH9X+Np8+N28MRwvquDJlekdOseq/Xlc+be1HC+q4MUfTOGuC5Pa7Ig4Be2PAAAgAElEQVQ4ICqIt+48n+lJ0Tzw9k5+/8Ee6twcUI8VllNcUe2WcxljeOyTfZz7h5W8vTX7rAz+GvSVUl6l/vn+hqyT/Pb9PRS18x96bZ1hY9ZJl2r5YHWIm+JIOyoutNlOdV1NRJg7Pp5P753JwkkJ3DFzKJeP7VwzfGOTBkVyzeQBvLAmi33Hi10+zhjDs6sOcPM/NxIXHsD7d01v16OX0AA/nr9pMjefP5jn12Tx180VZOWXduQtNPGvtVl86+HPmfDQpyx8+iue+Gw/mw+epLqDN4tPfJbOk59nYPOBe1/fzqIXN7c4O6WrDuSVcPtLm/nxK9uaXdWxu2nQV0p5nQWTErjxvEG8vP4Q3/rjSn7z3m4OFbgWKL4+XkRRRU1DIHfFdMfz5s6Oz++s6BB//nr1OH5x+cguOf/9l48kJMCXX/1nl0u12PKqWha/msYfV3zN7DH9eOvO8xnYp/0tIb42H35z5WgemT+WrMI6Lnt8FY99up+K6qZDGl1hjGHJZ+n85v09XDyqL7dfMISq2jqWrExn4TPrmPDQp9z64mZeWneQsirXHmf87fN0lqxM57uTE1h7/4X86jvJfHUgn0se/5I3Nh9pd62/orqWxz/dz+VPrGZdZgGf7z3BpY+v4n/e2UluceduJDpDO/IppbzS7+aN4buTB/DC2iyWbTjEi+sOcvGoWH4wLZHzhkS1+Kx7Y5Y1Pr/xynqtuWB4DHy4l/EDItyRda8VFWznvstG8uA7O3ln21HmT0xoMe2BvBJ+vHwbe48Xcd/sEdw5c2ine/5fO2UgAScz+OJ0JE+uTOedbdk8dOUYZo10fcREXZ3hdx/u4Z9rD7JgYgJ/WjAWX5sP92FNV/zVgQLWZOSzJj2fz/ae4NlVmfzhqrGtrhXxzJcH+Osn+5k/IZ6H55+DzUdYND2RC0f25b43t/PzN3fw4c5jPDx/LHHhbU+49NWBfH75zi4y80uZO74/v/x2MiLw5Mp0lm84zDvbjnLrjCHcfsGQDi/s1FEa9JVSXmtsQjiPXzOeX1w+kn+vP8S/1x/i0z0nSI4L48KRfQnytxHkZyPI7kug3UaQ3cane06QEBlIfITrs+ENjw1lxd0zGNEvtAvfjXe49twBvLb5CH9csZeLRsU2WZegorqWp1MP8HTqAQLtNl64+VxmjXDPMEaAiAAfllw7gWvOHcCv3t3FLf/axKXJsfzvFcltziNQU1vHL97eyZtbsrll2mB+9e3kMx7HRATZmTM2jjmORyMbMgt44J2dfP+FjcyfEM8vv5NMVKN5A55fnckj//2aK8b15y9Xj8PmdL7E6GBeu/1bvLjuIH/+aB+XPraKO2cNZWS/UAZEBhEfGXhGJ82TpVX84cO9vLU1m4FRQbz0gynWDaXDQ3PHcMu0RP768T7HDcAh7r4oif513dd3wGuCvojMBpYANuB5Y8wjjfaLY/8coAy42RiztdszqpTqdrFhAfz00hH8aNYw3t12lH99dZClqRm01OJ6zeQB7b5GZxbZ6Ul8fIQ/zBvDlX9bw2Of7OO3c8c07HOuoc4b35//+XZyp2cnbMn5Q6P57+IL+MeaLJ5cmc7Fj33JNZMHMHNEDFMT+xDcqAZcWWN17vx49wnuuXg4d1/U9gyEU4f0YcXdM/j7Fxn8PfUAqfvz+PUVyVw5rj8iwotfHeT3H+5lzth+PP7dMwN+PR8f4ZZpVq3//rd28OeP9p2xv0+w3brJjAxk3YECiitq+NGsofz4wiQC/JquvZAYHczS6ydy25HTPLxiL//7n91ckODLxRd2oBA7wCuCvojYgKXAJUA2sElE3jPG7HFKdjmQ5PiaCjzt+K6U6iUC/GxcO2Ug104ZiDGGypo6yqpqKauqobyqlrKqWiqqa3tNAO+oMfHh3ODoM3H15AH0jwhsqKEO6hPEy4umuHUdhJbYfX24M2UoV47vz8Mr9vLa5iO8uO4QfjZh8qAoZgyP5oKkGAb1CeKOf29hbUYBv74imVumJbp8jQA/G/deOoI558Txi7d2svjVNN7eepQpiVH85eN9XJIcy5JrJ7Q5C+SgPsG8ctt55BVXcuRUOdmnysh2+r4np4jR/cP53yuSGR7bdovR+AERvHr7eaTuy+No+i6X309neUXQB6YAGcaYTAAReRWYCzgH/bnAS8bqTbFeRCJEJM4YowtLK9ULiQgBfjYC/GxNmmxV23566QhW7DzG3a9u42RpFaWVNdw1axh3XTis2RpqV4qPCORv102korqWzQdPsTo9jy/35/Hnj/bx54/24WcT6gw8evU4FkxquR9Ca0b2C+OtO8/npXUH+cvH+/hyfx4XjuzL366bgJ+L0z6LCH3DAugbFsCkQZEdykfj880a2ZfU493Xp95bgn48cMTpdTZNa/HNpYkHNOgrpVQ7hQf68T/fHsU9r21n8qBI/jh/rEs11K4U4GdjelK0Na5/zihyiypYk5HPpoOnuGx0LCmd7FtgczTVXzq6H5/tOcE15w7A37d7b3A8Tbxh8gERuRq4zBhzq+P1jcAUY8yPndJ8CDxsjFnjeL0SuM8Ys6WZ890O3A4QGxs76dVXX21XfkpKSggJ0bminWmZNOVtZZKSMwuA1P5feOT63lYe3qAnlMmxkjpigwWfblqhryeUSXdzd5nMmjVrizFmcnP7vKWmnw0497xJAHI6kAYAY8yzwLMAkydPNikpKe3KTGpqKu095mynZdKU15XJcuubp/LkdeXhBbRMmtIyaao7y8RbJufZBCSJSKKI2IFrgfcapXkPuEks5wGF+jxfKaWUcp1X1PSNMTUichfwMdaQvReMMbtF5A7H/meAFVjD9TKwhuzd4qn8KqWUUj2RVwR9AGPMCqzA7rztGaefDfCj7s6XUkopdbbwio58XUlE8oBD7TwsGsjvguz0ZFomTWmZnEnLoyktk6a0TJpyd5kMMsY0O9HCWR/0O0JENrfU87G30jJpSsvkTFoeTWmZNKVl0lR3lom3dORTSimlVBfToK+UUkr1Ehr0m/espzPghbRMmtIyOZOWR1NaJk1pmTTVbWWiz/SVUkqpXkJr+koppVQvoUFfKaWU6iU06CullFK9hAZ9pZRSqpfQoK+UFxKR/4rI9z2dD+V+IrJeRCpE5BNP58UVIjJWREpEpFZEbvB0flTnaNBXXkFEDorIxY223Swia9x0fiMiw9xxru5gjLncGPMiuLccOkJErnf80y8RkXIRqXN6XeKpfLmDiMwWkQwPXPpWY8ylTvk4LiLTu+piIjJRRD4RkQIRqWhmf4yIvC8ipSKSJSIL6/cZY3YaY0KwVkNVPZwGfaVUq4wxy4wxIY5//JcDOfWvHdu8koj4iEiX/o8TEa9ZtKwNlcArwB0t7H8WOAX0BRYBL4hIUjflTXUjDfqqxxCR/iLylojkOWojdzvtmyIi60TktIgcE5G/iYjdsW+VI9l2R+30mhbOf5uI7BWRYhHZIyITHdt/ISIHnLZf5XTMzSKyVkSeEpFCEflaRC5y2n+L0zkzReSHja45V0TSRKTIcY3Zju2pInKriIwCngG+5cj7aRE5V0ROOAccEVkgImnNvKfzHLVIm9O2q0Rkh1O5bXZc/4SIPNaOX4nzdQaIyH9EJN/xPu9w2veIiCwTkdcc7yFNRBJF5NeO9AdFZJZT+vUi8jsR2eIo07dEJNxp/wwR2eAoi60iMq3RsQ+JyAasJbj7i8gPHb+XYhHJEJEfONL2Ad4Bhji1XPQRkVdF5JdO5zyjNcBRnj8Tkd1AUVvv3xsYY3YbY/4J7G28T0QigSuAXxljSo0xn2Mtc359N2dTdQMN+qpHcNTY3ge2A/HARcBPROQyR5Ja4B6s1aq+5dj//wCMMRc40oxz1E5fa+b8VwO/AW4CwoArgQLH7gPADCAc+C3wbxGJczp8KpDpuPavgbdFJMqxLxf4juOctwCPO91MTAFeAn4ORAAXAAed82WM2YtVO1vnyHuEMWaTI2+XOCW9AXi58fsyxqwHSoELnTZfByx3/LwEWGKMCQOGAq83PkdbHDcUK4CvgP7AbOBBEZnplOwqrJuXCGAf8LkjX/2AR4G/NzrtTVhBJx6wO9IgIoOBd4H/AaKAXwLvOgJXvRscx4cCx4FjWC0UYVhluVRERhtjChz5ynRquSjANddglX8fF9+/W4jIRY6bnZa+OrJoy0ig2BjjvBrpdmC0e3KtvIkGfeVN3nX+B8aZgeBcIMYY85AxpsoYkwk8B1wLYIzZYoxZb4ypMcYcBP4PaM8/3VuBPxtjNhlLRv0/QWPMG8aYHGNMneOGIR2Y4nRsLvCEMabasX8f8G3HsR8aYw44zvkl8AnWDQQ4mlGNMZ86zn3UGPO1i/l9ESu44bjBuIxvAnljrwDfc6QNBeY4tgFUA8NEJNoYU+K4SWiv6UCAMeZPjt/NfuCfOH43DiuNMV8YY2qAN7EC8KOO168CI0Uk0Cn9P40xXxtjSrBupL7n2P594G1jzGeOMlsB7AEudTr2eWPMPsfvo8YY854xJsvxO/gM+NKR58543PGZKHfx/buFMWal48avpa/NHThtCFDYaFsh1k2TOsto0FfeZJ7zPzAcNXWHQVhNtc43BQ8CsQAiMlxEPnA0vRYBf8SqebtqAFaNvgkRucnRJF1/3TGNzn3UnDmf9SGsGh8icrmjyfmk49g5Tse2eE0X/Bu4QkRCgO8Cq40xx1pIuxyYLyL+wHxgq1OtbhEwHPhaRDaJyHc6kJdBwOBGv5t7sWrx9U44/VwO5DmVWbnje7BTmiNOPx8CghxN/IOAGxpdazKO8m7mWETkShHZ6PQ7uJD2fTaa43wNV96/NyvBuglzFgYUeyAvqov1lE4oSh0BsowxLXUuehrYBnzPGFMsIj8BFraQtqXzD228UUQGYbUoXITVxF7reHYuTsniRUScgthA4D1HkH0Lq6n5P8aYahF51+nYZq/ZjCYLZBhjjorIOqzm6Rux3n/zBxuzR0QOYTVxOzftY4xJB77neHwyH3hTRPoYY0pdyFe9I8DXxpix7TimLQOcfh4IlBljCkXkCFZN/setHNtQXiISDLyB9Vn4rzGmRkQ+4pvfQXOLj5QCQU6vmwvezsd1xftvllgjXN5tJcksx+Of9vgaCBORgcaYw45t44DdHcmj8m5a01c9xUagSETuF5FAEbGJyBgROdexPxSrU1WJiIwE7mx0/AlgSCvnfx74mYhMEsswR8APxvoHnwdWxzysmr6zvsDdIuLn6BswCusZrx3wdxxbIyKXc2Yz9D+AWxzPaX1EJN6R98ZOAAni6Jjo5CXgPmAsVoe01iwH7sbqN/BG/UYRuUFEYowxdcBpx+baNs7V2BrHuX4iIgEi4isi59T3Xeigmx2tNyFYfS3q+2G8CFztKDOb47NwkYi0VKsOBPywHsHUiciVQIrT/hNAX8d16qUB3xGRCBGJB1q7wQD3vX+74/j6L1vjBI7HGiGtfDUb8B2f6QCszySO89sd5zwFfAA8JCJBIpKC1S9hWTvzr3oADfqqRzDG1GL1MB4PZAH5WIG6vlf3z7BqscVYNfPGnfV+A7zoaH79bjPnfwP4A1ZwLMaqTUUZY/ZgdSJbhxUgxgJrGx2+AUhy5OkPwEJjTIExphgr0L6ONRzqOuA9p2tuxNG5D+sZ6pdYTcWNfY5V6zouIvlO299xpH/HhZr5K1jB7nNjjPM5ZgO7xRpvvwS41hjTZBx3a4wx1ViPLc7HaorPw2p56MxwvpcdeT4K1AE/dVwrE1iA1aEy33G9xbTwv8zxXn+G1Qm0AJiHdUNWbzvW7+SQ47MRBbwAZACHsYLhK7TCje9/JdajjvqvB9p5fGtGOM65BetGtBzY4bT/NqxHHvnAv4BFjlYgdZbRpXWV6gQRuRlropUum1iljesfAH7o6KB2VhCR9cDfjDH/9nReuoKIfAlMANYaYy73dH7aIiJjgdVYrQSLjDGt3gQp76bP9JXqoURkAdajh889nRflOmOM24fydSVjzE6soZbqLKBBX6keSERSgWTgRsfzeKWUapM27yullFK9hEc78jl6x74p1hSZe0XkWyISJSKfiki643ukU/oHxJpGc598MxObUkoppVzg6d77S4CPjDEjscaF7gV+gTV7VxJWb9ZfAIhIMtYMV6Oxehz/vbkhLUoppZRqnsea90UkDGu4zBDn2cxEZB+QYow5Jtb85qnGmBEi8gCAMeZhR7qPgd8YY9a1dp3o6GgzePDgduWttLSU4ODgthP2IlomTXldmZzcYn2PmuSRy3tdeXgBLZOmtEyacneZbNmyJd8YE9PcPk925BuCNZ71nyIyDmv86GIgtn46UUfg7+tIHw84zwue7djWhIjcDtwOEBsby1//+td2ZaykpISQEK9dMdQjtEya8rYyScmxFqpL7d++z7u7eFt5eAMtk6a0TJpyd5nMmjXrUEv7PBn0fYGJwI+NMRtEZAmOpvwWSDPbmm2mMMY8i7U+NJMnTzYpKSntylhqairtPeZsp2XSlNeViWNyXU/lyevKwwtomTSlZdJUd5aJJ5/pZwPZxpgNjtdvYt0EnHA06+P4nuuU3nk+7gQgp5vyqpRSSvV4Hgv6xpjjwBERGeHYdBHWEpnvYS2fieP7fxw/vwdcKyL+IpKINe3pxm7MslJKKdWjeXpynh8DyxwLP2RizUPuA7wuIouw5r6+GsAYs1tEXse6MagBfuSYj12pXuuh9/eQfaqMZ2+a7OmsKNWi6upqsrOzqaioIDw8nL1793o6S16lo2USEBBAQkICfn5+Lh/j0aBvjEnDWgu7sYtaSP8HrAVNlFLAp3uPc6q0GmNMs51elPIG2dnZhIaGMnjwYEpKSggNDfV0lrxKcXFxu8vEGENBQQHZ2dkkJia6fJynx+krpTqooKSSIyfLKamsIa+k0tPZUapFFRUV9OnTBxG9NXUXEaFPnz5UVLRrUUwN+kr1VNuzTzf8nJXX1sq6Z8rMK6GwvNrdWVKqRRrw3a8jZapBX6keKu3wN0H/YIHrQd8Yw8Jn1vHkSl0uXaneRoO+Uj1UWnYhw2NDsNt8yMx3PegfL6rgZGkVB9txjFLq7KBBX6keyBjD9iOnmTgwkkF9gtrVvJ+RWwLAscL2PQtU6mzwxhtvMHXqVM455xyGDRvGb3/7W09nqVtp0FeqBzpYUEZheTXjB0SQGB1MVjtq7d8E/fKuyp5SXunFF1/kT3/6E2+99RY7duwgLS2NoKAgT2erW2nQV6oH2n7Eep4/zhH0DxWUuXxsfdA/VVZNRbVOdaF6h6KiIu69915ef/11EhISAAgJCeHnP/+5h3PWvTw9OY9SqgPSjpwmyG5jeGwoidHBVNXWuXxsuiPog9XEnxitK56p7vOnTw6Qnu/eVqbk/mH8+orRraZ55513mDp1KkOGDHHrtXsarekr1QOlHTnNmPhwbD7S7qB9ILeEgVFWk+ax09rEr3qH3bt3M378eE9nw+O0pq9UD1NZU8uenCJunjYYgMQY14P+qdIqCkqrmD2mH8s2HNbOfKrb3X/pUI/MyBccHEx5ud7kak1fqR7m62PFVNXWMX5ABAAxIf6E+Lt2/56RZzXtz0iKBrQzn+o95syZwxtvvMGJEycAqKys5LnnnvNwrrqf1vSV6mHqZ+Ib5wj6Iq438dd34hvdP5yoYLvW9FWvce655/Kb3/yGyy67jNraWmpqarjhhhs8na1up0FfqR4m7fBpYkL96R8e0LAtMToYTNvHZuSWEOhnIz4ikH5hARr0Va9y4403cuONN3o6Gx6lzftK9TBp2acZlxBxxrzbrtb003NLGBITjI+P0D9Cg75SvY0GfaV6kMKyajLzShk/IPyM7a4G/QO5JQzrGwJAv/AAfaavVC+jQV+pHmTHUet5/vgBkWdsdyXol1bWcPR0OcNirKAfFx7I6bJqyqt0gh6legsN+kr1IPUz8Y1NOLOmP9iFoJ/pmJ8/KbY+6Ft9ArS2r1TvoUFfqR4k7chphsYEEx7od8b2xq+bk5FXDNDQvB8XHgjowjtK9SYa9JXqIYwxpB0pbBiq114ZuSX4+giD+litAt/U9DXoK9VbaNBXqoc4erqc/JLKhkl52iv9RAmD+gThZ7P+7PvVB32dilepXkODvlI9xPYjhQBtBv3iiupmt2fkfdNzHyDAz0afYDvHirSmr1RvoUFfqR4i7cgp7L4+jOwX1mq6g/lNl9mtqqnjUEHZGUEfHMP2tKaveoHnnnuO8ePHM378eHx8fBp+vvfee1s9rry8nJkzZ1JbW8uRI0eYNWsWo0aNYvTo0SxZsqRJ+rbSnD59moULFzJy5EhGjRrFunXrAFiyZAljxoxh9OjRPPHEEwBUVVVxwQUXUFNT46ZS0KCvlFsZY/hgRw7V7Vjq1lXbjxQyun8Ydt/W/2wz80uabDtUUEptnSGp75kLncSFB7r8TL+4oppb/rmRjNxi1zOtlJe47bbbSEtL48MPP2TAgAGkpaWRlpbGY4891upxL7zwAvPnz8dms+Hr68ujjz7K3r17Wb9+PUuXLmXPnj1npG8rzeLFi5k9ezZff/0127dvZ9SoUezZs4fnnnuOjRs3sn37dj744APS09Ox2+1cdNFFvPbaa24rBw36SrnRtiOnuWv5Nj7efdyt562prWPn0ULGJbT9PD8rv7TJtvo59xvX9OPCXZ+Vb+vh03yxL49nvsx0Kb1S3mjXrl2MHTvW5fTLli1j7ty5AMTFxTFx4kQAQkNDGTVqFEePHj0jfWtpioqKWLVqFYsWLQLAbrcTERHBvn37OO+88wgKCsLX15eZM2fyzjvvADBv3jyWLVvWuTftROfeV8qNDhVYATcrr2ng7Yz9J0oor65lwsDOBf0hjZbhjYsIoLC8mrKqGoLsrf87SD9h1fDf357Dr76dTHhQ28MElWos9P3WH0912HUuLD4B7Ny5kzFjxriUtqqqiszMTAYPHtxk38GDB9m2bRtTp05t8fjGaTIzM4mJieGWW25h+/btTJo0iSVLlpCcnMzvf/97CgoKCAwMZMWKFUyePBmAMWPGsGnTJpfy6wqt6SvlRocLrOfjBwuaPlfvjIaV9TpY00/PLSE+IrBJYG/PsL30EyXYbT5U1tTx9rZsV7KtlNdpT00/Pz+fiIimf3MlJSUsWLCAJ554grCw5m9imktTU1PD1q1bufPOO9m2bRvBwcE88sgjjBgxgvvvv59LLrmE2bNnM27cOHx9rb9Vm82G3W6nuNg9j9W0pq+UGx05ZQX7+hq/u6QdPk1EkB+D+gS1mTYrrxRjzBkL8mTkljRp2genCXpOVzA0pul+Z/tzi5k4KILyqlqWbzjMzecPPuMaSrmi+IoiQkND207YRXbu3Mk999zT8Lqmpob77rsPEWHQoEHcfffdDfsCAwOpqDjzhri6upoFCxZw/fXXM3/+/Gav0VKahIQEEhISGmr+Cxcu5JFHHgFg0aJFDc3+Dz74IAkJCQ3HVVZWEhAQgDto0FfKjQ6ftIL+QTcH/e3NrKzXkuLKGgpKq4gO8Qegrs6QmV/C+UP7NEnr6lS8xhgyTpRw1cR4xvQP5763drDp4CmmJEZ14N0o5Rl1dXWkp6czcuTIhm1PP/00c+fOZebMmU3SR0ZGUltbS0VFBQEBARhjWLRoEaNGjWqx139rafr168eAAQPYt28fI0aMYOXKlSQnJwOQm5tL3759OXz4MG+//XZDr/6CggJiYmLw83PP4zRt3lfKjbIdQT+/pKrF8fLtVVpZw/4Txe2aic+5if/o6XIqquuaren3c7F5/1hhBcWVNSTFhvKdcXGE+vuyfMMhl/OjlDfIyMggISEBf3//hm1bt25l2rRpLR5z6aWXsmbNGgDWrl3Lyy+/zOeff94w5G/FihUAzJkzh5ycnFbTADz11FNcf/31nHPOOaSlpfHggw8CsGDBApKTk7niiitYunQpkZHWolpffPEFc+bMcVsZaE1fKTepqqnjWFEFI2JD2XeimEMFZYyJD2/7wDbsPFpInYEJ7Qn6eaWcO9iqhbfUcx/A39dGdIi9zaC/39GJb3jfEILsvsyfGM8rm47wv6VVRAXbXc6XUp40fPjwJkPs5s2bxw9/+EOioqJ44IEHiIo6s/Xqrrvu4rHHHuPiiy9m+vTpGNN8h8H6wN6/f/8W0wCMHz+ezZs3n7GtuLiY1atXN5t++fLlPPzww22+N1dpTV8pNzl6uhxjYHpSNACH3NSZL82xst45Ca7dQPjZhEynmn5rQR8cE/S00bxff46kWOtZ7HVTB1FVU8dbW7RDn+rZ5s6dyz/+8Q/+8pe/NAn4ABMmTGDWrFnU1nb/EtRVVVXMmzePESNGuO2cHg/6ImITkW0i8oHjdZSIfCoi6Y7vkU5pHxCRDBHZJyKXeS7XSjV1xNG0P32YFfTd9Vz/y315DOsbQp8Q/7YTAwOjgshymqAnI7eE6BA7EUHN18jjwgM5drrtmn50iL2hVj+iXyiTB0WyfOPhVms1Sp0NfvCDH2Cz2br9una7nZtuusmt5/R40AcWA3udXv8CWGmMSQJWOl4jIsnAtcBoYDbwdxHp/t+CUi2o78Q3Mi6UmFB/DjYzdK698ksq2ZBVwJwx/Vw+JjE65Ixn+um5xa32zI9zoaa//0RJk9n8rps6kKz8UtYdKHA5b0opz/Jo0BeRBODbwPNOm+cCLzp+fhGY57T9VWNMpTEmC8gApnRXXpVqy5FTZdhtPsSGBjC4T5Bbmvc/2X2COgOXj41z+ZghMcEcLCijts5Yve5zS0iKbS3oB1JUUUNpZfPze9efY3ijc8wZG0d4oB/LNh52OW9KKc/ydE3/CeA+wHmi8lhjzDEAx/e+ju3xwBGndNmObUp5heyT5SREBuLjIwzuE+yW5v3/7jrG4D5BjOzn+rjmxOhgqmrqyDldTl5JJUUVNQxro22sc8UAACAASURBVKYPLffgzymsoMTRc99ZgJ+NhZMS+GT3cfKKK13On+qd9DGQ+3WkTD3We19EvgPkGmO2iEiKK4c0s63ZdywitwO3A8TGxpKamtquvJWUlLT7mLOdlklTjctk96FyQuxCamoqdUVV5BZX8/FnX+Dv27EJbEqqDGszyrh8sB9ffvllm+lTHN8Ls9MB+M/n67A5Ll16LJPU1OaH2J04aXVQ+mjVBsZEN31itiPPagEozckgNTXrjH1DqaO61vCnN74kpW+VfkYa0b8bS0hICNnZ2YSHh1NXV+e22eXOFrW1te0uE2MMhYWFlJaWtusz5skhe9OAK0VkDhAAhInIv4ETIhJnjDkmInFAriN9NjDA6fgEIKe5ExtjngWeBZg8ebJJSUlpV8ZSU1Np7zFnOy2TphqXyelVnzBtVBwpKWMpicrhrfRtDBw9iVFxHZtr/PXNR6gzO7jj21MZ60rP/eXWt/kXT+ORjSsJix/quFPezYJLpjWMyW9sSEEZD2/8gr6Dh5MyeUCT/emrMoG9fPeyGUQ2MzzvvZx1bMyvYE6in35GGtG/G0t1dTXZ2dkcPXq0YaIb9Y2OlklAQADjxo1r18Q9Hgv6xpgHgAcAHDX9nxljbhCRvwDfBx5xfP+P45D3gOUi8hjQH0gCNnZ3vpVqTlFFNafLqhkYZU2TO7iPtbDNwfzSDgf9/+48RkJkIGPi23d8TKg/wXYbmY7peEP8fYkNa7nnf2y4ta+lHvxWz33/ZgM+WMP37n5lG3sK/LmwXTlVvYWfnx+JiYmAdSM0YcIED+fIu3RnmXj6mX5zHgEuEZF04BLHa4wxu4HXgT3AR8CPjDHdP3BSqWbUD9cb4Aj69XPkd3ThncLyatZk5DNnbFy757cXERJjgsnKLyXdMed+a+ewJujxb7EH//5mOvE5u2x0LH3+f3v3HR91fT9w/PW5Sy57X3ZCEkjYmzBkSdBaQNSfSFVUXHVVrdZqbbW1rbW2VltH1VoRnGhxD5QhIiAge4QNCQkQyB6EhOy7z++PS8JIApd1l3Dv5+ORx5H7fu97n3z4JO/7fsb742NiRVbzEwGFEF1Hlwj6WuuVWuvp9f8u0lpforVOqn8sPu28p7XWvbTWfbTWi51XYiHOlFVsC5gNd/p+nu6YfU1t3njn+3151Fo0U1qxVO90Dcv2Wtpo52y2ZXtN7/RtOffLSDrHNTzcjMxMjmFbvoW8E+ffrU8I4TxdIugL0d013ukHndoFL64dM/gX7cwlMsCToXZspduchBBvskoqyC+rbkXQb3qnn11axckaS5OZ+2e7fmQPrBq+3H6sTeUVQjiGBH0hOkBWSQV+nm4EeJ+aUBPXxrX65dV1rDpQwJSBERgMbZv5nxDqQ8NqnnMt12vQ0p1+Y8798wT9BLMPPQMMfLZVgr4QXZkEfSE6wJHiisau/QbxIT7klFZRVdu6qSff78unps7KtFYk5DlbgvlUoLfrTj/Qi7KqOsrPStCT1hj0z3+NsVFu7MstY2/OiVaWtqmqWgtLduVKtj8hOpjssidEB8gqrmiSprZhMt/hogr6tCK5zuKdOYT6eTCiR9D5T25BQv3qAZOboXFy4bk0JOjJLa0k8bSf40BeOaF+Hi3m7T/d6Eg3PjxQy+fbjrVpxYLFqll3sIgvtx9jya5cyqrrCPByZ+sTP8HYxh4PIcSZ5E5fiHayWjVHSyrpEXJmcE0w1y/ba8W4fkVNHSv3FzBlQNu79gECvN0J8THR0+xjV8CMDPACIPusZXtp55nEdzo/k2JSnzC+2HYMi9W+TGFaa1KzjvOXhXsY8/fl3DRvA4t35fLTgRHcNbEnpZW17DpWate1hBDnJ3f6QrRTQXk11XVWYoO8zng+LtgW9Fszg3/V/gIqay1MHdS2WfunmzYoErOdO/OdutM/FfS11qTll3NtMwl7WjJjWDTL9uSxNr2Qib1Dz3v+S8vTePG7NExGAyl9Q7lqaDST+4bh6W6ksLyaOT9ksCa9kCGxbZvQKIQ4kwR9IdqpYeZ+zFnd6AHe7gR5u7dqrf6iXbmE+JgYFd90X+/Weur/Btp9bri/J0pB9mkz+I8dr6SixnLOzXrONrlfGP6ebny+7dh5g37+iSr+u+ogl/UP57mZQ86YBAlg9vWgf6Q/Pxwo4L6URLvLIIRomXTvC9FOWSW2oH72RD6wLduz906/qtbC93vzuGxABG5Gx/5qmtwMtgQ9p3Xvp+WVA+efuX86Dzcjlw+OYsmu3BZ37Wvwyop0ai2ax6f1axLwG0xIMrP1SMl5ryWEsI8EfSHa6UiR7e44OtCrybH4EG8OFdp3p//DgQJO1liY2saEPO0VGeBJzmnJdRqW69k7pt9gxvBoKmstLN2d2+I5WcUV/G/jEa5NjiW+fu5Dc8Ynmam1aDZmFrd4jhDCfhL0hWinrJIKIvw98XRvukNdXIgP2aWVdi3bW7wrlwAvdy7qFdIZxTyvyABPco6f6t5vzcz90yXHBREb7MXn21pes//S8jSUUjxwybm77UfGB2NyM7A6rbBVZRBCNE+CvhDtdKS4gtjgpnf5YJvBrzUcLTn33X51nYXv9uZxWf9w3B3ctd8gMsDrjIl86flldq3PP5tSiquHRrM2vbDZtLzp+eV8tvUos8fENa4aaImnu5HRCcGsSS9odTmEEE1J0BeinY4WV5yRfvd0jRvvnKeL/8f0Isqq6tqVkKe9IgM8Kauuo6yqFqvVNnP/7NwD9rp6eEyLaXlf+O4AXu5G7p3Uy65rjU80cyCv/IwPJEKItpGgL0Q71NRZyTlR1WICnMYtds8zmW/hjmz8PN0Ym+icrn2AiNOW7TXM3G/NJL7TJZh9GBob2CQt765jpXyzI4fbxycQYudywvFJZgDWpEsXvxDtJUFfiHY4drwSrWkx6Ad6u+Pv6XbOHPxVtRa+3Z3H1IEReLg1nRfgKFH1ExGzS6tIy6+fxNeG7v0GM4ZHN0nL+/yyAwR4uXPHhJ52X6dfhD8hPibWpEkXvxDtJUFfiHZoWKPf3HI9sI1vx5vPvdvein35lFfXccWQqE4po70i/G13+jnHKznQsFyvjd37ANMHR+FmUI0T+rYcLub7ffncfXFPAryaX6LXHINBMT7JzJr0Iqx2ZvoTQjSvXUFfKfWEUurhjiqMEN3NkYYtdVuYyAe2Lv5zBf2FO7Ix+5q4qKfzuvbB1r2vFOSUVpGWV06Yn0eL6+ftEexjYlKfML7cbkvL++yS/Zh9Pbh1bHyrrzU+0UxheTX7csvaXB4hRPvv9GcDr539pFLqDqXUY+28thBdXlZJBSajgXA/zxbPiQ/x5lhJJTV11ibHyqpqWb43n8sHRTo8Ic/Z3I0GQn09yCmtJC2/rM3j+aebMTyavBPVPLt0Hxsyi7k/pRfeptYnAp2QZMvuJ7P4hWif9v6VqdRaNzdY+R5wUzuvLUSXd7S4kpggr3NujhMX4oO1hWV7y/bkUV1ndXrXfoPIAE+yj9vu9Nsznt9gct8w/DzdeH1VBtGBXswa3aNN14kI8CQxzFfW6wvRTu0O+kqpJmuMtNbVgOTNFBe8I8UVTXLuny3efGqL3bMtTM0mOtCL4e3YRrcjRQZ4sT3rOJW1ljYv1zudp7uR6YNtfyIevCSpXRMVJySZ2ZhZbFeiIyFE89ob9P8FfKmUijv9SaVUGNC0L1OIC0xWSQU9zjGeD7Y7fYDMwjPH9UtO1rA6rZDpQyLbtY1uR4oI8KS8Ps99WxLzNOfeSYk8MDmRGcOj23WdCUlmquusbD5U0iHlEsIVtWuXPa31x0opb2CLUmo9sB3bB4mfAX9uf/GE6LoqajXHK2pbTMzTIMTHhJ+HW5ONdxbtyqHOqrlicNfo2geICjw1NyGpA8b0wbac8deX9Wn3dUYnhOBuVKxOL2hcuy+EaJ12zxzSWr8DJAAfAe5AFTBLa/1+e68tRFdWUGnrzGppjX4DpRRxZu8mW+wuTM2mZ6gPA6L8O62MrRVRnxY33N+jVcvqHMHHw41hPYJYI+P6QrRZu+70G2ity4B3O+JaQnQXhZW2NeMtrdE/XVyID7uPlTZ+n1taxYbMYh68JAmlukbXPkBUfVa+jhjP7wwTk8z889sDFJVX253RTwhxiiTnEaKN8itsQf983ftgW7Z3tKSSWoutd+DrHdloTZeZtd+gIRVvR8zc7wzj65furT1Y5OSSCNE9SdAXoo0KK634e7rZlcAmLsSHOqvmWIlt69qFO3IYEOVPr9CuFVwjA7y4YkhU44z7rmZQdAABXu6sPiDr9YVoiw7p3hfCFRVUaGKDfew6N8F8auMdpSA16ziPTe3bmcVrE6NB8fKsYc4uRouMBsW4xBDWpBeite5SQyNCdAdypy9EGxVUWu3q2odTW+weLqpgYWo2ANO7WNd+dzE+MZSc0ioOFpx750IhRFMS9EW38d2ePG54Yz3FJ2ucXRSsVk1hpaZHiH1BP9TXA2+TkUNFJ1mYmkNyXBDRgede3y+aN6Fhq13ZdU+IVpOgL7qFmjorf/pqNz8eLOKRj1OdvttaQXk1tVaIDbIvcCuliAvxYfnefPbnlXHlULnLb6vYYG/iQ7z5KjXb6e1AiO5Ggr7oFj7cnMWx45VMHxzJ9/vymbsmw6nlyWrcXc++O32wzeA/UlyBQcHUgV1zolx38YtJvdh65Djvrjvk7KII0a1I0BddXlWthVe+T2NkfBAvzxrGlAERPLtkP1sOOy8d65G2BP36yXzjEs2E+ska8/a4NjmWSX1CeWbJPjIKyp1dHCG6DQn6osubv/4weSeqefiyPiil+MfMwUQGevLA/7ZxvMI54/tZxbald60Zl4+vH//vamvzuyOlFP+4ZjAebkYe+TgVSyd28+eXVfHmmky+2ZFDen5ZY64FIbojpy3ZU0rFYsviF4Ftc545WuuXlFLBwIdAPHAIuFZrXVL/mseAnwMW4AGt9VInFF040MnqOl5beZDxiWbG9AwBIMDLnVdmDWfmf3/kkY938MbNIxy+dCurpIIgD4Wnu/27xl3SL5xbLjrRZdfAdzfh/p48eeUAfvXhdt5YncE9F/fq0OvX1Fl5a20mL3+f3rgJEYDJaKBnqA+9w/3oE+HHsNhAxibKXgCie3DmOv064GGt9VallB+2TXuWAbcCy7XWzyilfgf8DvitUqo/cD0wAIgCvlNK9dZayz6bF7C3fzxE0ckafn1Z7zOeHxIbyO+m9uOpr/cwb00md0zo6dByHSmuINS7dR80zL4ePHnVwE4qkWu6amgUS3bl8vy3B0jpE0afiI5JH7xiXz5/+XoPmYUnuaRvGI9O6UutxcqBvDIO5JVzIK+MLYdL+Kp++eVvftqH+1ISO+S9hehMTgv6WuscIKf+32VKqb1ANHAVMKn+tHeAlcBv659foLWuBjKVUunAKGCdY0suHKW0spbXVx3kkr5hze43f/u4eNZnFPGPJftIjg9maGygw8p2tLiCBB8ZHXM2pRR/vXogl73wAw9/vJ3P7x2Hu7Ht/y8ZBeU89fUeVuwvoGeoD2/fNpJJfcIajw+MDjjj/LKqWv745W6eW7ofk9HAnRMd++FTiNbqEhn5lFLxwDBgAxBe/4EArXWOUqrhNy4aWH/ay47WPycuUPPWZHKiqo6HftK72eNKKZ6bOZjL/72G+97fyqIHJtiVEre9ik/WkHOiipHmrrULnasy+3rwt6sHcs/8rby6Ip1fXdp8e7FaNbuySykqr6G6zkJVrZXqOgvVdVaqa60cKa5gwaYjeLgZ+f20ftwyNh6T27k/QPh5uvPczMHUWKw8vWgv7kbFreMSOuPHFKJDKK2du85VKeULrAKe1lp/ppQ6rrUOPO14idY6SCn1KrBOaz2//vl5wCKt9afNXPMu4C6A8PDwEQsWLGhVmcrLy/H17Vo50Z3N0XVSVqP5zaoKBpqN3D/M85znHjxu4W8bqkgMNDA1wZ0BZiPuhs4Z4y+v0Ty7qYrsk1YeGqQZENl12smk7BQAVkatcMr7O/v35vXUKjbmWnhijCfxAba5FlprMk9Y2ZhTx8ZcC8VVLf+9U8D4aDdm9jYR4NG69lNn1byWWs2WPAu3DjAxKdb2gdDZddIVSZ001dF1kpKSskVrndzcMafe6Sul3IFPgfe11p/VP52nlIqsv8uPBPLrnz8KxJ728hggu7nraq3nAHMAkpOT9aRJk1pVrpUrV9La11zoHF0nf1+8lxprBs/cOI7E82zzOgnwiTrC04v28uLWavw83Li0fzjTBkUyIcncqsl251J8soYb524gtxLm3ToKnb27a7WTD2wPziqTs39vho2q5ScvrOKDDDf++bMhLNmVy9c7cjhSXIW7UTExKZTLB0fSM9QXDzeD7cvdiGf9o4eboV1DAxMnWrln/hbe3p3PgH59+VlyrNPrpCuSOmnKkXXizNn7CpgH7NVaP3/aoa+AW4Bn6h+/PO35D5RSz2ObyJcEbHRciYWj5JdV8c6Ph/i/odHnDfgNrh/Vg6uHR/NjehGLdubw7Z48Pt92DB+TkUv6hXPD6B6Ns//bovhkDTe8sZ7MwpPMvTmZib1DWdnsR07hLAHe7vzjmsHc9vYmrnxlLUaDYmyvEO5PSeSnAyI6fejH5GbgPzcO5853N/PopzswuRkIOP/LhHAoZ97pjwNmAzuVUtvrn3scW7D/SCn1c+AI8DMArfVupdRHwB5sM//vk5n7F6b/rDhIrUXz4KVJrXqdh5uRlL5hpPQN428WK+sOFrF4Vw5Ld+fxVWo296X04qFLe+PWyru5MwL+LclMqN/TXXQ9KX3DeGbGIOqsmqkDIwjxdWwSJE93I3NmJ3P725v49Uep3D3Y1DgrWYiuwJmz99dgG0ZrziUtvOZp4OlOK5RwuqziCj7YcIRrk2OIC7Fv29rmuBsNTOwdysTeofzpCgt//mo3r644yOZDJfx71jDC/c89T6CBBPzu5/pRPZz6/l4mI3NvSebWtzbyemoJ0yeeoH+Uv1PLJEQDWXMkuozSylrueGczJjcD909u3V3+uXi6G3nmmsE8f+0Qdhwt5fJ/r2ZteuF5X1dUXt0Y8OfdMlICvrCbj4cbr89OxtekeOjD7VTVSqek6Bok6IsuobrOwt3vbSajsJzXZ4/olG1nZwyP4av7xxHkbeKmeRt4YdmBJulbj5ZU8Pm2ozz22U6mv7ymMeCPT5KMa6J1gn1M/Hygif15Zfxz6X5nF0cIoIus0xeuzWrVPPxRKuszinnxuqGM68SUpknhfnx5/zj+8MUuXlqexubDxUwbFMnmQyVszCzm2HFbTn0/TzeS44L4xaRERiUEd1p5xIVtcKgbs8dEMHdNJpP7hnVaul6LVfPK9+mMTzIzIq5pIishGkjQF07398V7+XpHDo9N7cv/Dev8fEveJjf+9bMhjEkI4Ykvd7E2vQizrwejE4K5c0ICIxOC6Rvhj7GT1voL1/L4tH6sTS/kkY9TWfyriQR4dfwqgmeX7uP1VRm8sTqDD+8ew4AoWTcgmidBXzjV3NUZvLE6k1vHxnOXA1OYKqW4dqRte9aTNRbiQ7wdvmmPcA1eJiMvXDeUGa/9yJ++3MWL1w/r0Ot/lZrN66syuGpoFJsyi7n1rU18es9YeoTYv+2zcB0ypi+cZmFqNn/9Zi/TBkXwxPT+Tgm6Yf6eJJh9JOCLTjUkNpAHJifxxfZsFqZ2XIKHXcdKefSTVEbGB/HczCG8+/NR1NRZufnNDRSWV3fY+4gLhwR94RTrDhbx8EepjIoP5vlrh0pXurjg3ZfSi6Gxgfzhi13kllY1e051nYW16YXsyz1x3usVlVdz93tbCPI28Z8bR2ByM5AY5sebt44k90QVt7+9iZOnbQksBEjQF05wpKiCu97bTFyIN2/cnNxhaXKF6MrcjAZeuG4oNXVWfvNJKtb6lSM5pZV8sOEId767mWF/WcaNczcw9aXV/OGLnZRW1jZ7rVqLlfs+2EpheTWvzx5BqN+pJEQj4oJ49Ybh7M4+wT3zt1BTZ3XIzye6BxnTFw735MLdWK2at24b6ZBd8YToKhLMPjwxvT+Pf76Te9/fyqGik+zLLQMgOtCLGcOjubh3GOsOFvH2j5ks2ZXLE9P7c+WQqDOGoJ7+Zi/rM4p54bohDI5puqX0Jf3C+fuMQTz6yQ5+80kqL1w7FEMX703bn1tGkLc7YXYmzhJtI0FfONTyvXks35fP49P6EhMkE42E65k1KpYV+/NZtjeP5LggHpval5S+YSSF+TYG9p/0D2fG8Gh+//lOHlywnY82Z/HUVQPpGerLR5uzePvHQ9wxPoGrh8W0+D7XJsdSUFbNc0v3Y/b14A+X93PI3JWKmjpKK2uJDLAv10ZBWTX/WLKPT7YcJczPgw/uHENimOzC11kk6AuHqaq18OTCPSSG+XKb7DkuXJRSitduHE6NxYq3qeU/wQOjA/js3nF8sOEwzy7Zz5QXV3P9qFgWbMxifKKZ303te973undSLwrKqpm3JpM1aYWk9A1jct8whvcIbPUeFPbYcriEX36wlezSKiYkmbltXDyTeoc128tQa7Hyzo+HeOm7NKrqLNxyURzf7Mzl+jnrmH/HaPpGSOriziBBX7RLVa0FDzeDXXcQc37I4EhxBR/cMbpdW5gK0d25GQ12BV2jQTH7onh+OiCCv36zl3fXHSY22IuXZw2z6/VKKf44vT/xId4s2Z3L3NUZ/HfVQQK83JnYO5SUPqFc3Du03RsTWa2aOaszeG7pfqICPbk/JZGPt2Rx+9ubiQ/x5pax8cwcEYOfp204b01aIX9euJv0/HIu7h3KH6/oT69QX24eG88Nb6zn+jnrmf/z0QyMlnwDHU2Cvmiz4pM1XPbCKob1COKVG4bh4dbyhLys4gpeXZHO5YMjOy0rmRAXqjB/T/49axi3jYsnIsCTIB+T3a81GBS3jkvg1nEJnKiqZU1aId/vy2fl/vzG5YO9w31Jjg8mOS6I5LhgYoO97B4KKCqv5uGPU1m5v4CpAyN45prBBHi58+ClSSzZlctbazN5cuEe/vXtAWaOiGHnwSq2LNlAj2DbRN5L+4U1vlevUF8+uvsibnhjAze8sZ53bh/FsB6SYbAjSdAXbfby92kUnaxh2Z48fjF/K/+5cXiLM/Gf+noPBqX4w+X9HFxKIS4c7Q2A/p7uTBsUybRBkVitml3ZpazaX8CmwyUs3J7NBxuOABDq50FyXBAj4oIYGB1A/yh//D2bTrrdkFHEAwu2UXKylqeuGsBNY+IaA7i70cAVQ6K4YkgUqVnHefvHQ7y/4TAKzSOX9eaOCT2b/XsRF+LDh3eP4YY3NjB73kbeum0kI+MlFXZHkaAv2uRIUQXz1x/muuRYBscE8vjnO7n7vS28PntEk1/klfvz+XZPHo9O6WP35B4hROcyGBSDYwIbZ/9brJoDeWVsPlzClkPFbDpUwuJduY3nx4V4MyDKnwFRtg8Bu46W8sJ3B+gR7M28e0eesyt+SGwgL1w3lCem9+fHtWuZfp5dNGOCvOvv+Ndzy5sbmXtLMmN7SQ9hR5CgL9rkuW/3YzQoHvpJb8L9PTEa4Hef7eTOdzefsfa+us62l31Psw93jHdcml0hROsYDYp+kf70i/Rn9pg4APLLqtidfYI92SfYdayUXcdOsGjnqQ8CVw6J4m8zBuHrYV8oCfYx4Wuyb9ggIsCTBXeP4cY3NnDbW5v4xaRejE80MzgmEJNb++YEVdVa2J19gkHRAe2+VncjQV+02o6jx1mYms39KYmE16+pvW5kD5RS/PbTHdz57mbmzE7Gy2Rk7upMDhVV8O7to1zul0uI7i7Mz5OwPp6k9AlrfO5EVS17sk9gtWou6hXSqcsAw/w8WXDXGO7/YBsvLU/jxe/S8DYZSY4PZmyvEMb2CmFAVIDdGT0tVs2nW4/ywrID5JRWYfb14PqRscwa3aNTtvPuiiToi1bRWvP3RfsI9jFx98Vn3rlfmxyLQSl+80kqP39nE3+5aiAvf5/GlAERTOwd6qQSCyE6kr+nO2N6hjjs/UJ8PfjfXWMoOVnDhswi1h0s4seDRTyzeB9g2wY7pU8Y0wZFMqlPaLPzBLTWLN+bz7NL93Egr5whsYH86tIklu3J49WV6fxnZTqT+4Zx05g4JiaFdvlERu0hQV+0ysoDBazLKOJPV/RvXH5zupkjYjAoeOTjVC7/92qUgieu6O+EkgohLiRBPiamDIxkysBIwDb0sD6jmDVpBSzbk8dXqdn4mIxM7hfO5YMimNQnDE93I1sOF/PM4n1sOlRCT7MPr904nCkDI1BKcd3IHhwtqeB/G4/w4aYsvtubT2ywF9cMj6FvhB/xZh/ign3wMl04qcIl6AsAjh2vxN2oCPNrOQWmxar5x+J99Aj25sbRcS2eN2N4DEaD4tcfpfLIZX1cpttMCOE4YX6eXDkkiiuHRFFrsbI+o4hFO3NYujuPhanZeJuM9A73Y3vWcUL9PHj66oFcmxzbJEdITJA3v/lpXx68pDdLducyf/1hXvwu7YxzIvw9iTd7k2D2ITLAC2+TES+TES93I94mI57uRrxNblTXWcgprSK3tKr+sdL2eKKKUF8PUvqGkdInjOT4IKflKpGg78Kq6yws25PH/zYeYW16EZ7uBn5/eX9uGt2j2XG6z7cdY19uGS/PGnbe8fmrhkYzqXeY5NYXQnQ6d6OBCUmhTEgK5amrrGzILOabnTlsOVTCI5f15vbxCefMfghgcjM0fog4UVXL4cIKDhWd5FDhSTKLTnK4qIJvd+dRdLLGrjKZfU2E+3sSE+TF8LggDhed5K21mcz5IQM/TzcmJoWS0jeMix089ClB3wVlFJSzYFMWn2w5SvHJGqIDvXjo0t5sPlzME1/sYvnePJ69ZvAZG1/UWDTPf7ufwTEBXD4o0q73kYAvhHA0N6OBcYlmxrUjCZi/fUxrkAAACR5JREFUpzuDYgIYFNN0GWKtxUplrYWqGgsVNRYqa22PVbUW3AyKqEAvwvw9mk1WVl5dx5q0Qlbsy2fF/ny+2ZmDUjA13p1Jk9pc3FaRoH+BWXWggLfWZmJQCpPRgMnNgHv9o4ebgb05J9iQWYybQXFpv3Bmje7B+EQzRoNCa8276w7zt0V7+emLP/D3GYMax8++O1xLdmkt/7x2yAU9yUUIIc7F3Wj7m9pcsqLz8fVwY8rACKYMjEBrze7sE6zYl4+l6HAnlLR5EvQvIGvSCrnznc2E+JoI8TVRU2c99WWxUl1nJdTXg0en9GHmiJgm4/dKKW4ZG8+4RDMPfbide+Zv5ZrhMfzq0iS+zqglpU+oJMgQQogOoJRiYHQAA6MDWLnymMPeV4L+BWLL4WLufHczCWZbCstAb/tzc58tMcyXz+4dy8vL03hlRToLU7OptcBv7djVSwghRNcl2VK6KK01OaWVaK3Pe+7u7FJufWsTEQGevHfHqHYF/AbuRgO/vqwPH98zlnizN5fGuclWl0II0c3JnX4XtO5gEc8u3ce2I8cZERfE49P6MiKu+Q0n0vPLuXneRvw83Jh/x+hzLrlrixFxQXz70MWsXLmyQ68rhBDC8eROvwvZdayUm9/cyKw31pNbWsUvJvUiq7iCa15bx93vbeZgQfkZ52cVV3DT3A0opXj/zjGyHl4IIcQ5yZ1+F5BRUM6/lh3gmx05BHq78/tp/Zh9URye7kZ+OTmReaszef2HDC574QeuHxnLg5cmoTXcOHcDlbUWFtw1hgSzj7N/DCGEEF2cBH0Hqa6zUFReQ2F5NUXlNRSUV1NYXk1aXjlfpWbj4Wbgl5MTuXNizzOWgnib3PjlJUnMGt2Dl5en8f6GI3y+7RjBPiZKTtYw/47R9IuUsXYhhBDnJ0G/E1itmrT8cjZkFrE+o4hNh0ooKKtu9lw/Dzdmj4njvpREQv08Wrym2deDJ68ayK3jEvjn0v38cKCAubeMZFiPoM76MYQQQlxgJOh3gDqLlX25ZWzMLGZDZhEbM4spqagFICrAk/GJZhLMPph9PTD7mjD7eRDq64HZ16PVGzkkmH149cbhWK1akuQIIYRolW4X9JVSU4CXACMwV2v9jCPf37aUrortWcdtX0eOs+PYcapqrQDEBHkxuW84Y3oGM6ZnCDFBXp2y37QEfCGEEK3VrYK+UsoIvAr8BDgKbFJKfaW13uOI9//DFztZtiePvBO2rnqT0UD/KH+uH9mDYT0CGREXREyQtyOKIoQQQrRatwr6wCggXWudAaCUWgBcBTgk6LsbDYzpGcLQ2ECG9QiiX6Rfs5sqCCGEEF2RsifjW1ehlJoJTNFa31H//WxgtNb6/rPOuwu4CyA8PHzEggULWvU+5eXl+Pr6dkyhLxBSJ011tTqZlJ0CwMqoFU55/65WH12B1ElTUidNdXSdpKSkbNFaJzd3rLvd6Tc3kN3kU4vWeg4wByA5OVlPauWehStXrqS1r7nQSZ001eXq5APbg7PK1OXqowuQOmlK6qQpR9ZJd8vIdxSIPe37GCDbSWURQgghupXuFvQ3AUlKqQSllAm4HvjKyWUSQgghuoVu1b2vta5TSt0PLMW2ZO9NrfVuJxdLCCGE6Ba61US+tlBKFQCHW/kyM1DYCcXpzqROmpI6OZPUR1NSJ01JnTTV0XUSp7UObe7ABR/020IptbmlmY+uSuqkKamTM0l9NCV10pTUSVOOrJPuNqYvhBBCiDaSoC+EEEK4CAn6zZvj7AJ0QVInTUmdnEnqoympk6akTppyWJ3ImL4QQgjhIuROXwghhHARLhv0lVJvKqXylVK7WjiulFL/VkqlK6V2KKWGO7qMjmZHnUxSSpUqpbbXf/3R0WV0JKVUrFJqhVJqr1Jqt1LqwWbOcal2YmeduFo78VRKbVRKpdbXyZPNnONq7cSeOnGpdgK2nWKVUtuUUl83c8whbaRbJefpYG8DrwDvtnB8KpBU/zUaeK3+8UL2NueuE4DVWuvpjimO09UBD2uttyql/IAtSqllZ23l7GrtxJ46AddqJ9XAZK11uVLKHVijlFqstV5/2jmu1k7sqRNwrXYC8CCwF/Bv5phD2ojL3ulrrX8Ais9xylXAu9pmPRColIp0TOmcw446cSla6xyt9db6f5dh+2WNPus0l2ondtaJS6n/vy+v/9a9/uvsyVKu1k7sqROXopSKAS4H5rZwikPaiMsGfTtEA1mnfX8UF//jVu+i+i67xUqpAc4ujKMopeKBYcCGsw65bDs5R52Ai7WT+m7b7UA+sExr7fLtxI46AddqJy8CjwLWFo47pI1I0G+ZXdv4upit2NI7DgFeBr5wcnkcQinlC3wK/EprfeLsw8285IJvJ+epE5drJ1pri9Z6KLadP0cppQaedYrLtRM76sRl2olSajqQr7Xecq7Tmnmuw9uIBP2WyTa+Z9Fan2jostNaLwLclVJmJxerU9WPR34KvK+1/qyZU1yunZyvTlyxnTTQWh8HVgJTzjrkcu2kQUt14mLtZBxwpVLqELAAmKyUmn/WOQ5pIxL0W/YVcHP9jMoxQKnWOsfZhXImpVSEUkrV/3sUtvZT5NxSdZ76n3UesFdr/XwLp7lUO7GnTlywnYQqpQLr/+0FXArsO+s0V2sn560TV2onWuvHtNYxWut4bFvCf6+1vums0xzSRlx29r5S6n/AJMCslDoK/AnbZBO01v8FFgHTgHSgArjNOSV1HDvqZCbwC6VUHVAJXK8v7OxO44DZwM76sUmAx4Ee4LLtxJ46cbV2Egm8o5QyYgtcH2mtv1ZK3QMu207sqRNXaydNOKONSEY+IYQQwkVI974QQgjhIiToCyGEEC5Cgr4QQgjhIiToCyGEEC5Cgr4QQgjhIlx2yZ4QouMppUKA5fXfRgAWoKD++1Fa6xqnFEwIAciSPSFEJ1FK/Rko11r/09llEULYSPe+EEII4SIk6AshhBAuQoK+EEII4SIk6AshhBAuQoK+EEII4SIk6AshhBAuQpbsCSGEEC5C7vSFEEIIFyFBXwghhHAREvSFEEIIFyFBXwghhHAREvSFEEIIFyFBXwghhHAREvSFEEIIFyFBXwghhHAR/w/jDZVFdnTQQAAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# Plot the energy, magnetization, and heat capacity vs temperature\n",
    "\n",
    "fig, axs = plt.subplots(3,1,sharex=True,figsize=(8,8))\n",
    "\n",
    "axs[0].plot(points[:,0],points[:,1], label=\"E\")\n",
    "axs[0].set_ylabel(\"Energy\")\n",
    "axs[0].set_title(\"Energy vs Temperature [L = 10]\")\n",
    "\n",
    "axs[1].plot(points[:,0],points[:,2], label=\"$|m|$\")\n",
    "axs[1].set_ylabel(\"$|m|$\")\n",
    "axs[1].set_title(\"Magnetization vs Temperature [L = 10]\")\n",
    "\n",
    "# heat capacity\n",
    "# C = var(E) / ( k_B T**2)\n",
    "heat_capacity = points[:,3] / (points[:,0]**2)\n",
    "axs[2].plot(points[:,0],heat_capacity, label=\"$C$\")\n",
    "axs[2].set_xlabel(\"T\")\n",
    "axs[2].set_ylabel(\"$C$\")\n",
    "axs[2].set_title(\"Heat capacity vs Temperature [L = 10]\")\n",
    "\n",
    "for ax in axs:\n",
    "    ax.axvline(x=Tc,linestyle='-', color=\"orange\",linewidth=2.0, label=\"$T_c$ ({:.3f})\".format(Tc))\n",
    "    ax.legend(loc=\"best\", numpoints=1)\n",
    "    ax.grid(True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Energy vs magnetization\n",
    "Blue data is low temperature ($<T_c$) and red data is high temperature ($>T_c$)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYYAAAEICAYAAABbOlNNAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df3CdV5nfv4+u5IyswCpOjIvsXClgonK9/OjGDQTC4F2pS6JtCzsL7WQ0lNIFNTKhnc5kS4ynCx0qZpmdnYZfsquQbKASMJ3d2S3LejdM3A2bQFnWQALEJDtOLMmKKfnhqdtYbmxLT/94r+J7r85zXr1H5773vdL3M3NGes/7nnOe99yr99H7POc8j6gqCCGEkBU6Wi0AIYSQYkHFQAghpA4qBkIIIXVQMRBCCKmDioEQQkgdna0WYL1cc801OjAw0GoxCCGkrfjBD37wvKpud51re8UwMDCAY8eOtVoMQghpK0RkzjpHUxIhhJA6qBgIIYTUQcVACCGkDioGQgghdVAxEEIIqYOKgRBCSB1UDIQQQuqgYiCEEFJHbopBRO4TkWdF5KfGeRGRz4nICRH5sYj8Sj5yrS6EELKZyfON4X4At3jO3wrgddUyBuBQswWylACVAyFkM5ObYlDVvwZwxnPJuwF8RRO+B6BXRF6dj3SEEEJWKJKPYSeAUzXHC9W6VYjImIgcE5Fjzz33XC7CEULIZqFIisFlwHEmpFbVKVXdq6p7t293BgdsKvRLEEI2MkVSDAsArq053gXgdItkMaFfghCy0SmSYvgGgH9RXZ30VgBnVfXnzRxQne8jdj0hhGwGcsvHICJfA7APwDUisgDgEwC6AEBVDwM4AmAEwAkAiwA+mIdcVAKEEFJPbopBVW9LOa8APpKTOE3DZVJSBUSWUe9GUagW6YWNEEIS+GSKiO1/WFEK9SWpJ4SQYkHFkJEwv8SKMkirI4SQ1tP2OZ9bAf0ShJCNDBVDAaBfghBSJPiUyQXF6r16SR39EoSQosE3hogk/+W7zrj0b5qPwfJLEEJIc+EbQ2RUV5dmsEUuQERfLlvkwsvn9u95CJ1yCSKKTrmE/Xseao4QhJANCRVDG7JFLuAiulBrYrqILmyRC9i/5yEcOv5OLKETgGAJnTh0/J1UDoSQNUPF0IZcVgq1JMph6vjNznNJPSGEpEPFsMFYQilTPSGENELFkAMl45lcKtnngsfCklm/c+sLdX6JnVtfiDs4IWRDQMWQA2Njdv3goNXKvcS1q2PJPAcoBruecp5TLOP0+W2o9UucPr+NyoEQsgouV82Bycnk59QUsLSUvCWMjSX1neYn4F7i2ndtJ+bmXNcnD/snL74WLh/DsuGXSJQFIYRcRrTN4zvs3btXjx071moxgsma4EckbQmswk6G565X5f4IQjYbIvIDVd3rOkdTUovJ6mMol1P6M3wMaezcWZ+qdGc12/b+4Sfq90QMPxHUPyGkfaBiaDGW/2FoCNi6tb5u61ZgYgLo63O36esDdnSfheV/cPoscAE7dwKnG5Konj4NbN1yAYeODtbviTg6SOVAyAaHiqHFTE4C4+OX3xxKpeT4wQcTn0R/f/IffH9/cjw6CnR1ufvq6gJOn78adohv196HLauUwgrnL7r9ElNHd2e4Q0JIu0EfQxvS0eH2M/j9D1l9D/429EsQ0t7Qx7DBsPwMaf6HWIT6MQgh7QEVQxsysvsJuPwFI7ufQKXibtPd8f+cbYb6foYO81vg9kt0UDEQsqGhYmhDjjx0JVy2/yMPXYlz59xtXiUvYAjfQq0jegjfwoNdI1g2UzzYfglCyMaFG9zakPkl97Kk+aU+YN5uM4trHSfoKyCE1MM3hjakXHIvIyqXTtv+B6PNehwTM8P3YUDm0CHLGJA5zAzf560HgP37k93eIsnP/fuDhyeENAtVbetyww036GZjevxh3YoX69IBbcWLOj3+sE5Pq27dWp8qaOtW1emhe105hFTHx7VScZ8ClhRYbqhb1i68pNND9zplGO+dccs2dK+Oj5siEEJyBsAxNZ6rLX+wr7dsRsWgmiiH/tIpFSxpf+mUTo8/fPnctGp/v6pI8nN6WpNfXE/l/n5V1VXKoVJRLZXcTUol1X7Mus/honsYzHr7I4Tki08xcB/DZsG3+cHwPvviOAmWoU5LpHvvg319tVV7fw0JaTu4j4EA5TL24/PoxEUIltGJi9iPz7/sY5iZAQYGEv0xMJAc+/JIlHHKfc5YylrGKW9/M/sfwUDnQuKX6FzAzP5HLl/gEo4Q0jS4KmmTsH/3X+LQ3CBW/ptfQicO4SPA7mG8fSaJ2bS4mFw7N5cc79ixOoYSkOSQePXpZzH3v8uofztQ7MD/wmnsXFU/UjmJv0Mnjh7vW3Vu8BXPYOzQP8AiepLxl3Zh7NBVAB7B6Nvn3MIBSXwQQkh0aEraJHR2JrkgGimVgF27YOR4cFMqAVi6VA2u14jblNTfD2BhAXNLu1b3B3df/aUFzO662S1cfz8wO7t2oQkhdRTGlCQit4jIkyJyQkTucpz/JRH5MxF5TEQeF5EP5infRsalFFbq5429D76+suaQnp+3919YfSX7MqyNGRmFJoSsmdwUg4iUAHwRwK0AKgBuE5HGAA4fAXBcVd8EYB+APxARbrONgNdfkHErQ6mUPV5SuWzvpTD9EqXTrQ8MRcgmJM83hhsBnFDVp1X1AoCvA3h3wzUK4BUiIgCuBHAGwKUcZdyw+PJOj4y4z1l5H/btA/ZVfgFXHKW+jp8760cu/ikmxmaxFfUxO7biHMYqjzjrJ8ZmbeGsekLIuslTMewE6payLFTravkCgNcDOA3gJwD+raqakXzI2rHyPkxOAkeOuNv84hfu+hMngBPnGh3MACD4xfKrnPVHTr8Zo5M3Y2r8R+gvLUCwjP7SAqbGf4TJx/c560cnb7aFs+oJIesmT8VgJQOo5V0AHgXQB+DNAL4gIq9c1ZHImIgcE5Fjzz33XHxJNyiTk8ClS8megUuXkmPANtdbfon5eU8by1+AxPQzOnkzZi/twrJ2YPbSruThD2D07XOY3XUzlqUTs7tuTlYjVQebwW0YwEl0YAkDOIkZ3PayADN7Po0Bma2G35jFzJ5PJ/Vc4UpIONbOt9gFwE0AHqg5PgDgQMM1fw7gHTXH/wPAjb5+N+vO55hYm6Ktncr9/apXX3neea4Dl4ydzydtAcw4HtM63fNhd4iNng/rdGXCHZaj70+s7gghVVCEkBhI9kw8DeA6AFsAPAZgT8M1hwB8svr7DgDPALjG1y8Vw/qxnsvj4+bzWq+W550KoAdn3Q/yyoQtgCdcR3/HvPtUx7z242S2sBz9ec0oIcXHpxhyMyWp6iUAdwB4AMDPAPw3VX1cRG4Xkdurl30KwNtE5CcAjgL4mKo+n5eMm5XRUXd+6clJO+/0Gb3K2dcirsRU5bPox2ziL8AspiqfxejjH7cF8CxJnV9udENVTy3vfNk81YhpzuIKV0LWRK47n1X1CIAjDXWHa34/DeDX85SJ+BkddW8wLpdOOzerlUunMfr4x3G5yQAAj1IAkqWnrk1s5TLKC/Y4WLqEOQysOlfCknPD3MsrXGdmgIMHE01RLgMTE9xFTUgNjJVEMFMNiTE3lxhdVqJO+By2I/tehHNZ6r4XswvgWZI6MnjCPc7gCUxUvupe/tr3TWzdWt9i69bk+R90s4RsNiwbU7sU+hjWT0pE7mhtQjrrL51ynyqdUlXV6cqE9uNkEn4cJ1/2ZThDj0cXnJD2BQy7TXwEROQOahMiQIdecobrFixjWQNeeKMKTkj7UphYSaSYhESdiBqpwtOZL40pgOy5Qhlig5BUqBgIJiZg2+QjtjHx+Rh8voz9+4FDhy7vxFtaSo59yiGq4IRsUCwbU7sU+hjiYNrkI7dx4vMx2Kf8uUdzEZyQ9gX0MZBC4/MxYNl2Cagn92ibf68JaTb0MZDi4Api5Ek76nUJ+GKJ+7D8Ej5/RUjwJQZsIu2K9SrRLoWmpDbCiL0xXvkrBZYbLELLOj70M18YJdW+Prcpqa/PlmF83N2mUnHXj497YzllvVearUhRAE1JpBAMDDh3OHfionOncqmURIE1NypLgCnJynFq4ct96ksvatwrU5KSouAzJVExkPwwfAmCZbijsqe4CkIUg6+NbxzT0ZHHRg9C4kMfAykG5bIzt4KV2nPFVTAzfB8GZK6ac2EOM8P3pQ41s/8RDHQuJG06FzCz/5H6TteKL/dpuWz7EbhfgrQzlo2pXQp9DO2DlT9hqPf7tnl/6F53GO+he00fw3T3b7vbjD9s+yW6u7P7GHxxyeljIAUHRcjH0KxCxdA+9GPWvScBszo+fnlbQqmUPHPT2lj7GKw8Df2lU+6Hf60ScAmh6t77kBZ3ifslSIHxKQb6GEhudMhy5rhH3jZG3oUOLGVuAyD73gf6EUgbQx8DyRfD7l7GKeflVn1qm1LJ6bMow52Rx4q7FEyaH4H7GEibQsVA4uLJd7C77xzgiHuU1LuZqMw4cy5MVGYws28KY7gHcxiAogNzGMAY7sFI73fdbfZ9C+jrcw9k1fvwxV1i3gfSzlg2pnYp9DEUDI/dPSi0UX+/TuO2+pwLuM0fR6l0ymwTHF/JwvIjMO8DKTgoQs5nsknw5G+29pUtLcE2u8zPYxRfwyyuwzJKmMV1GMXXknzQ1lBLfWYbrxC+kBhZw3t75iGoP0LyxNIY7VL4xlAwQt4YOpbspZ2+yKtX/1/3qY55+7/1jg73Od9KJSuMxtBQkNxmf7WroAhpMuByVZIbnvX75vOw5377Ierpb/rqj7r3K/R82H5g9/RkUwylkm1+skqK3NHNWYQEQMVA8sWzft+5VUDE/aAU8fcn4vYliHjbZHrIh5Q0uX1tCckJn2Kgj4HkyuRkEhhPNfk5OYnUZZ/7vzOKzoVZiC6jc2EW+78zWnfe2W50NAlWt7yc/BxNaQM4l76iVMoeRmNljO98B1hYSG52YSE5BvzhwrnElRQBS2O0S+EbQ8GIHKLaZ46fHn/YDn2RcazpvjvdfVUm7JDcfX32vfoED/FZEBIZ0JREciN0maZhdvGZ44NXhDrG6i+dcvdVOuUXwjIXpfkRXDY1LnElOeJTDAyJQeISOUyEL0p2SDRsi5DQGwDcAqwIkbUNQ2yQHGFIDNIUnObwyOGmfeZ471A+W73jnBUuo1w67RfC2o+QlnbUSHFq31DKPWWFvgziw3qVaJdCU1JrMN0C4w9HtZN7fQwhMlg+Biu89/jDYelAQwTPK4w3Q4IT9ZuSWv5gX2+hYmgNXnN45HDTWaNhe4XznJsef1j7S6eSpa+lU/VObJcQIX6EtMnLI8QGfRlE/YqBPgYSRKHN4T7hgHiCh/gR0uTLI1VooT88kheF8TGIyC0i8qSInBCRu4xr9onIoyLyuIh8O0/5yNrJM3NlZnO4T7iYghthv1P3PYTIEFNuph0laVivErELgBKApwC8BsAWAI8BqDRc0wvgOIBy9fhVaf3SlNQa8jJTB43jaxQxTpGVqnS6MhH/pmJOOGM1EfWbkvJUDDcBeKDm+ACAAw3X7Afwn7L0S8XQOvLIXBlzr8L6OnTI5tv7kEbI5MWacPoYiPoVQ24+BhF5L4BbVPVD1eP3A3iLqt5Rc83dALoA7AHwCgCfVdWvOPoaAzAGAOVy+Ya5ubkc7oC0gujm8IgdhqQqLQT0MRAUx8fg8tQ1fjs7AdwA4DcAvAvAfxCR61c1Up1S1b2qunf79u3xJSWFIbo5PLRDx34F794HwO8cGR5O+lopw8PmONGhj4GkkKdiWABwbc3xLgCNf1kLAP5SVc+p6vMA/hrAm3KSjxQQX/bMIEZGstUDycP50KHLSX6WloBDhzAx+BV3CtGxWX9qz+Fh4OjR+jGOHgV27nSOE105hMwB2VxYNqbYBcnbwNMArsNl5/OehmteD+Bo9dqtAH4K4Jd9/dLHsPGJ6ssIsa979iuYex9847jqfSV2ngb6GIgWxMcAACIyAuBuJCuU7lPVCRG5vaqgDlev+R0AHwSwDOBLqnq3r0/uYyCZCLGvx457FPI3F/PvlD4GguL4GKCqR1T1elV9rapOVOsOryiF6vHvq2pFVX85TSkQkpkQ+3pa3COjP+ceh8D9El6y+iXoYyApFHjpBCFNIMS+PjaWrR7AzMX3YQz3YA4DUHRgDgMYwz2Yufg+oFJxN+rtzTyO5f/wKgf6GEgalo2pXQp9DCQTofZ1X8Am1zA46R4GJ/0yZBwnKH80fQxEC+RjaAb0MZBM5GRf9+5xkM54MjDvAwmkMD4GQnLD2kMQal+3+jPqy5h3D4P5uDKk+D9m9nwaAzKLDlnGgMxiZs+nw/M+hOZwYO6H9sN6lWiXQlMSWUXsWEkB+RPMHNJ9d8aNldTX576fSsWO5dR3Z/a8D75cEaGfBWkpKEKspGYVKgayisB8DJn7s+z71XGmcZv242SyxwEndRq3XR4n6+aMrPsfSiW/nyNrLCnfvYZ+FqSl+BRDZ6vfWAhJZWYGOHgQmK+aYSYmgNFR+/p5txnHrA89t7ISyLh+FHMYxdcaznl8Aj588rlYWsI83CajeZSB0QH3HAbeq0nIZ0FaTqqPQUTKayyvzENgssnwhZaw6Omx60Ps+9u2ues7jD+fctlus21b2D1l3WMg4vdzWDTGH0kjTS7umWhL1uJ8/jKA+6s/rXI/gPc0RUKyuTl4EFhcrK9bXEzqLc6ds+tjBl/q7g7rK+SeLLktVDFR+ao7llPlq3a78+fd9SJh9xo92BXJBcvG1C6FPoYNjojbRi1it/HZ3lWz2/d9Mlh9+dqE3JMld8q9Tlcm6v0caUmEfP2FBq3KI3EHyQxiOZ8BdGW5Po9CxbDBiRz0zktIch9PG9P5HNJf7Hu1NtKF9kfajiiKAcCXAJwBcArA3wC4B8BH19q+WYWKYYMTstwx5pLUwCWc0+MPu5eKjj8cd0motVy1ry9sfoaG3OeGhtI/K9JWxFIMT668MQDYCWAEwF1rbd+sQsWwCQgxRWQNLZH2ZuKSwdMm9UUnY38mKaYkJ763Ai4v3TT4FMOaQ2KIyJcB/I6qPhvPw7F+GBKDRCEkTISnTQeWM3eXW0hwXxsrLDjDZWw4YoXEmALwbRG5U0TeISK/FEc8QpqALxS1K0RDyLJKT5ugqBPNWNqZNYwGl5cSIJMp6SSA/wjgLgBfA/A4gKfW2r5ZhaYksgqfDT2mfd/jlwhxWUz33emWu1KxZfD5BALCaAT5Z0hbgkg+hocddVestX2zChUDWUWoDT3El+Fpk3mRE066T/j8BaqrlcOKozggjAZ9DJsHn2LI4mP4PIATqvrZ6K8t64A+BrKKgtvQTVcClrEMw8yzxr/TNQ3kowDzQ/Ihlo9hB4DbReS0iHxTRCZE5H1xRCQkIgW3oZsi+EJV+MgaYtwibX6yphAlbcuaFYOq/jNVfT2A6wD8LoC/A3BjswQjJJh9++z6AoRomBh5xB2qou+L7gZDQ3ZnvrhLVqpOK7Xo2JjdpqcnewpR0r5YNqaVAqC8xvLKtL6aUehjIKsI2ZOQs3zmrmjLX+Dpy7zXkBSiIX4J0pZgPT4GEfkrAArAZbhdqVcA96vqV2IoqyzQx0BWUfTUlTHl8/UFZB8nxC8R4v8gLWddPgZV/VVV/bXqz8byazU/c1cKhDgJTV3pI2a6y5ipPX19+c5ddVWiIFbKVVetTYZGLH9OmtyhbZgmNB+sV4l2KTQlkVUEbSII2K8Qmu4yp/Si5jhdXe763l57X0Rvr7s+zdRVhPkmTsDUnmTTERIp1aII6S4DI7xm8hf47skqed5r6HwTJz7FsOZ9DEWFPgaSiZgxkSzSfAWR4zJF9RdkJc97DZWBOIm1j4GQ9idmTCTLvr6edJdZ9ySEtPGxFp9B4/gh58tlYHi43s8xPOxvEzrfJDNUDGRzEbKPwWpj7Zew9gKknd+9296TYMkwMpK9jZWruqvLvqfeXltmH9b5s2eBo0fr644eTZSDNT/79rV8D8qmwbIxtUuhj4FkJlZMpNC4QqE29BAZXG1CYkn5fBI+svosVmTP6k8hmUFRfAwicguAzwIoAfiSqv6ecd0/BPA9AP9cVf/I1yd9DKRlhO5HiGlDj53DwYqV5EM1eUM5eBCYn09MOxMTwOiof6ysMtCXEBWfj6EzRyFKAL4I4B8BWADwtyLyDVU97rjuMwAeyEs2QoIolxPTjavex7ZtwAsvrK63Hoi+/rq6gAsX3PUWvgdvTw/w4ot2W1eblbAci4tJ3Yo5yzeWD2t+tm3L1g8JJk8fw41IorM+raoXAHwdwLsd130UwB8DKFSmOEJWETvuUk9P9v5cSsFXD9gPalXg3Dn3OV9fBw9eVgorLC4m9TlaJEg88lQMOwGcqjleqNa9jIjsBPCbAA77OhKRMRE5JiLHnnvuueiCErImRkeBqSmgvz/5z7i/PzkeHfW3O3PGXX/uXFh/MQl5kM8bUWGt+jSs+bHqSXTyVAxWrKVa7gbwMVVd8nWkqlOquldV927fvj2agISYWEtCR0eB2dnE9j07W/8QD1l6mge+sORZl6umheqO3V/WlK0kDMsrHbsAuAnAAzXHBwAcaLjmJIDZankRiTnpPb5+uSqJNJ2YYR18YTlCQj50d7tX8HR322186UCttJ/WOL50qb6wHNY4vpSkltxpMhAnKEJIDCSO7qeR5HPYAuAxAHs8198P4L1p/VIxkKYTO4SFarzlr75lnyGy+fqzQnVb97OCq11ey2wZLsPEpxjyXq46gsRcVAJwn6pOiMjtAKCqhxuuvR/AN5XLVUmraXUIi9Clp9bftm8c3/Mg5rMir2W2XOJqUpiQGKp6RFWvV9XXqupEte5wo1Ko1v/LNKVASC4hnUNs3iGhv5vhe3CFnWjGOL7PwTVHIfMTmrKVvofsWK8S7VJoStrE5BXSObbNO6Q/C8tW39dn92e1GR9XrVTc5yqVsM/B8jFY4/jmzpK7UmGo7gBQBB9DswoVwyYmz5DOWW3eVpv1yJB1Dlz1vpKWDjREhpAw3qFhOWKGLdkE+BQDw26T9qXVIZ1D7PuxZQj1F1jjAPHmNLYMPqzri57mtYUUxsdASFRaHULbZ/P2ESpDs/0FaelAQ3wwIWG8Q9OLhvgyiBMqBtK+hISkiBnSeSUe0Frr1yPD8LA7TLUvjMbQkHucvj5brpjhvScmgB073GNZsZx277b7s+QeHEyUwKFDwFJ1b+zSUnLc02PfK7GxbEztUuhj2ORkDcMcO6Szb21/TBl8Nnaf3I0O6KGh5uyxCJE7q7/A59OJnZJ0EwD6GAipUgSbc177FWKOHzpvIWG3rfsJ6cuCPgb6GAh5mVCbc8z9Eq1ec9+MPRa+GEZZWPEXZN3HEDslaRFo5f4L61WiXQpNSSQTseMehbQJWXMfEhMppmy+PRbWXgXfvojeXrtNiAzWuUqlPfcx5BD7CdzHQEgNMf0SoW2y2vF99vKYcxCyxyI0hpHln4ndph3Tgeaw/8KnGOhjICSNIsRK8v2dxvwbDtlj4bs+JIZRnj6QopLD/dDHQMh6iLlfohmxkkLs+1n9H749FrFjGDVjvtuNVt+P9SrRLoWmJNJ08vIx+GzoIXkSQmQL8Qn48jtY/fX1ZY+vFHpP7Qh9DFQMpA0IsVPHtOOHxg9ykWd8pdj9xf6MikyT74eKgZAiIZLtYSmSvQ2QjOVy1lp9hYzhk60Z/YVsKPQ9YDeaMsmATzHQ+UxI3gwMJCElGimVLod0qKW/H3jhBeDFF9c+hghw++1JWIhGtmwBLlxYXX/11ck4WbjyyuSxfe7c6nM9Pe56Hz09wEsvAZcurT5nObLHx4HJSXd/MzNJ6I7Fxct1W7cCU1PJ79a52tzdGxSf85mKgZC8sR5WH/gA8OUvux9U73+/+6How1I0FiGKYSWgn2ulTEdH4gx3KSGLjo7sq25KJbciAWwl3N+f/LTOzc5mk6EN4aokQorE6GjysO/vTx6q/f3J8eSku350NLtSALIpBQA4cyb7GKr2g3x5Gbh4MVt/IUsxffc5P2/X+85tcqgYCFkPscMWjI4m/60uLyc/V0waWUM+5BUmwjeOb7mqr78QGSx8yz7zTnHaRlAxEBLKiknIFYo6drvBQXe9Fb56cDAJ4+2iUrHDZFuhrS127LDH2bfPDm9dqbjrx8aA7m73uQ7jceULc26N7wsx7gu17iP0+1BELK90uxSuSiItI3RZZUi7rGGl08JRWKtxso6z0l/I0lNrhZFvrKyrkkLDhYfQZmlEwVVJhDSB0LAFscNuW4SEo4g9DhBXhqzPqzxDZbRZWA46nwlZC1ntw6GhqENs2yE+Bt84lmwhNv5QO36IPd4X/qOV4cJr+1xrfZGxXiXapdCURKIQO+yFL6xDyFhW2ImuLnd9X58tgxWmwhd6o6PD7ssXEsOSYWjIngMr/IZ1r6GhukNCb8T+DrUQcOczISnEDsPgC0Xta5dVPl8JCaMRMkZI2G1f+I+sMoSOk/YZhdBGO6l9ioE+BkKA+PbhmHZyIHs47BUZmv33HRp2O9b1zRqnzZ+La4E+BkLSiG0f9q3tB+L5M2K2CfUvhOxj8IXxzkroOGmf0SaGioEQIP6admtt/cq69qzr3a31+Na+g6Eh+556e91tduzIvodgZMTeYzE4CJw96z5n7Vbu6Ulkd2HJvW8fsHu3LYP1ufo+o82OZWNql0IfA4lGbPtwSBpKC1+bRgfr0JD/nrLa8NfjYwjpU9V9T7458Mng+1xDorVuEEAfAyEFotWpK0P2Kvj6iv0MsforQurTDURhfAwicouIPCkiJ0TkLsf5URH5cbV8V0TelKd8hJgErLk3m2yk1JV52up9cxBbhg0S8ygY61UidgFQAvAUgNcA2ALgMQCVhmveBuCq6u+3AvibtH5pSiJNJ2B9urdJzH0MtWajteLry9rHYBXffonxcXtPglV6e8Pktsbx9WfRZvsRQkER9jEAuAnAAzXHBwAc8Fx/FYBn0vqlYiBNJ8AnkNokqz8j9pp7yy8Rso9ANbs/Ja2/rHMQ0p9Fm8U8CsWnGHLzMYjIewHcoqofqh6/H3/MJbkAAA74SURBVMBbVPUO4/o7Afz9lesbzo0BGAOAcrl8w5wr2QYhsQiw70cPmxN7X4RFyH4J3/Ux+wv1jWQdv81iHoVSFB+D61N1fmIi8qsAfhvAx1znVXVKVfeq6t7t27dHFJFsCmLHRIrTxE+oDd2616xxnNKwxonZX+w5sM6lfXi+/kJiLxXRn2G9SsQuWKMpCcAbkfgirl9LvzQlkUyE2I8DYupEN1OHxPWJGT/IipXU3e2/2RAfQ4jc1jjd3bZsIeOExsfK7YuydlAQH0MngKcBXIfLzuc9DdeUAZwA8La19kvFQDIRew+Bh+hhc2LlIrCKL+aQr/jmJ2tfaf1Zc5B1z4QvV4QvvlLoXoqsn1EO/gyfYsh1H4OIjAC4G8kKpftUdUJEbgcAVT0sIl8C8FsAVpwGl9Swga3AfQwkE63eQ5AnecVX8uVjCHm+5JFHwpcrIqRN6F6KFn63iuJjgKoeUdXrVfW1qjpRrTusqoerv39IVa9S1TdXi1cpEJKZouwhiJkHwCIkVlLMmEyh85NHHglfrghffKXQvRRZfTAt3p/CWElkcxESEyl2HKX9+4FDhy7HC1paSo5jKwcrvpIVc2hwMIlVlIXeXjtO0e7dduwli0rFluHsWXveXvEKdxsrJ/bIiC33jh12G18OaSv39eCgHRvL118rsWxM7VLoYyCZCTH+x3QYNCMPgIs88jH42pVKYbb/WOOnjRPTL9EMn0WTQRGcz80qVAwkGnklWfE9fHwyWJvSLIesSPYHZqtLXjKHjCNit/OdC+2vyfgUQ2dr31cIKQgrobAXF5Pjldd9ABgdjTuWz1FpyfCHfwgcPVp/7dGjwM6dwOnTl+tWzCtAYu46dy6OXKF0dDTfQR8id7mczG8Wtm1Lfr7wgvvclVe6+yyV3GHGV/wIrjb0MRBSAA4evPxAXmFxMamPje8hZsnQqBRWqFUKtUxNAefPZ5fLyoVgsWWL/3xWH0OIYkprY/mHYkaZBWxf1NiYLUNs/1UsrFeJdik0JZEo5PlKX1TzCpDI5zJZ+eanKLJbxTLPhXwOad8TayyfibBFeaLhMSXxjYEQYH1hELISsrQyZIys46xc/+CD9Y+9Bx8MX1aZh0kkNLR2zCWuafc5OgrMziZmtdnZevOk71yLoGIgBPC/0oek4vRhpY4cGrJlsEw8VmrPsTF7+aQvRaaFb0mqtSS0qyv7ssuOjrDUnta57m77s7Put1KxPwff8tLY35NWYr1KtEuhKYlEw3qlb8aSQmslUcxVSSHLJy1CQ16Hht3Omtoz6xhpS0VDvgttFq4bRQmJ0QwYEoM0nc0SEiNmyInadlmfMdb1oak9LbmssULDo4T010IKExKDkEKTZ9iCrOGwfW2yyu0L+WD1FxryOmbY7ZifQ5q/IGtY8vX4HzZz2O1mFZqSSBRih1MOGSsk3PP4ePaw0lb6zkrFblOpuNukpdW0+rPa9PWFtbHmzmozNGTfk9UmLex2SPrVzR52u1mFioFEIU/bccxw2Gn+Apet3DdWyDhpIT5cMoTMt6+ouv0Sob6RrPejyrDbRYI+BhKFPG3HeYTDDvUXhIzju946FzLfPqzrY29i88kVkn6VYbcJKTDNsB1nHcvCFw47zV+QlZBx0vwPw8PJg26lDA+HzXcIob4RXxvX/YSOxbDbhBQY3z6G2GELrLXwlYq7fmzMlsFai+/bP2CNU6nY7fbts+dgcNDdZnAweWi6YjxdvOhuMzJi36sVXsPawwDYITu2bLHnwepvbMy+n+Fhe3+KVQ8wJEazCn0MJBp5hS0ISV1pyRA7VWnI2v6YdnyfbyTEhh/iT/F9Dml+jqzpV617zQHQx0BIgYhpV46dqhTI3l9MO37oOCE+htgpRNvsWUofAyGtIGRfhG9Ne8ja/qxtQmSLacdfGd9lx08bZ8+e+jZ79qxtLFd9EfcW1NJs+axXiXYpNCWRQuJbnx5zT0LI3ofY+yV86/etPRPd3bYMVn9dXe76vj57T4JVfHssfPdqjVOpFOO7lQFwHwMhOZNm+4+1tj8kbWTI3odQv0QeqT1DSjPuNS8iyeBTDPQxENIM8rL9W4S2ie2XyPp8aUYWOWucdo2HFMlHRR8DIXkTsj49xPbv218QuvfBZasP9UvkkXsihNifQ57xkHLY+0DFQEgzCFmf7st5YPW3Y4e7TU+PvSfB2ncwMpIogePH6+uPHweefdaW7exZ97mzZ235rBwOvvwOFkND2VOI9vTY53yfgy8fg0XsPA157H2wbEztUuhjIIUl6/r0kJhDIfZ6n48hq60+ZK/CevqzclKEjBXyOcTeNxJKhL0PYGpPQlqAL2Wjy7SwtOTux6pPY34+W3/W9T5CZQvt77vf9R9nJevnYM3RSr0rXHdamxCanQ7U0hjtUvjGQNoOa7lhR4f9n6rVxvcf8ZVXZvsP+uqrs//XbckcWkSyt7GWvqaVrHPa0WHP0dVX2+HZr7jCbtNCwDcGQgrEwYPA4mJ93eKibScfG7PbWFxxBXDu3PrkXAvd3fZu4JAd0arZ25w/n70NYM+pJXeaH2Nqyl3/0kvZZWsxVAyE5I1lQlhcBMbHL6/kKZWS48nJ7GaHCxeyP2TPnMl2PeBXTkXHmlNr3hYX7Tk6cya7WS1kvnMiV8UgIreIyJMickJE7nKcFxH5XPX8j0XkV/KUj5Bc8C03fPvbgV27kv9ad+1Kjn1tQsJhZ5UrrU3M5ZOlUlgojZBxYi4BbsZ8tzAsR26KQURKAL4I4FYAFQC3iUhj3NtbAbyuWsYAHMpLPkJyw1rauHu3vawxJOy2Fe55aMhe7pg1FPXIiB1C++JFuz/LLLNvn72c1qKry5bPWvo6NhYWYty3VDRkvi1iL3HNiuV8iF0A3ATggZrjAwAONFzzXwDcVnP8JIBX+/ql85m0HXmFqlC1w0D7ljs2xgOqVPzj+By8MVOF+sZQXZ2rubfXPwchIcbT5i5kvrN8RyKG3kARQmKIyHsB3KKqH6oevx/AW1T1jpprvgng91T1kerxUQAfU9VjDX2NIXmjQLlcvmFubi6XeyAkCllTe4aGb4gZosE3ju9eYqcKtdrETLGaZ3gLixxkK0pIDJerv/HO13INVHVKVfeq6t7t27dHEY6Q3IiZpjOv1JCh4xTVN+Jr1+K0ml4ZcpItT8WwAODamuNdAE4HXENIe2PZqcfGstui80oN6RvHlyo05F5j2up9FDWtJtB62SwbU+wCoBPA0wCuA7AFwGMA9jRc8xsA/gLJm8NbAXw/rV/6GEhbYtmcQ0Id5JUaMqtfIq1dHrb69dxTq2mybCiCjwEARGQEwN0ASgDuU9UJEbm9qqAOi4gA+AKAWwAsAvigNvgXGmHYbUIIyY7Px9CZpyCqegTAkYa6wzW/K4CP5CkTIYSQerjzmRBCSB1UDIQQQuqgYiCEEFIHFQMhhJA6qBgIIYTUQcVACCGkDioGQgghdeS6wa0ZiMhzAGJE0bsGwPMR+mkmlDEOlDEOlDEOrZKxX1WdwebaXjHEQkSOWbsAiwJljANljANljEMRZaQpiRBCSB1UDIQQQuqgYrjMVKsFWAOUMQ6UMQ6UMQ6Fk5E+BkIIIXXwjYEQQkgdVAyEEELq2LSKQUTuFBEVkWtq6g6IyAkReVJE3lVTf4OI/KR67nPVhELNlO1TIvJjEXlURL4lIn3V+gEROV+tf1REDte0KYSM1XNFmcffF5EnqnL+iYj0VuuLNI9OGavnijKP7xORx0VkWUT21tQXaR6dMlbPFWIeHTJ/UkSeqZm/kTSZc8NK7baRC5K80g8g2Rh3TbWugiTd6BVI0o8+BaBUPfd9ADchSTn6FwBubbJ8r6z5/d8AOFz9fQDAT402RZGxSPP46wA6q79/BsBnCjiPloxFmsfXAxgE8BCAvTX1RZpHS8bCzKND5k8CuNNRb8qcV9msbwz/GcC/B1DreX83gK+r6kuqehLACQA3isirkTwE/6cmn9pXALynmcKp6v+pOexpkHMVBZOxSPP4LVW9VD38HoBdvusLJmOR5vFnqvrkWq8vmIyFmccMOGXOU4BNpxhE5J8CeEZVH2s4tRPAqZrjhWrdzurvjfVNRUQmROQUgFEAv1tz6joR+ZGIfFtE3lGtK5KMhZrHGv4Vkv8KVyjMPNZQK2NR57GRIs5jLUWfxzuqZsT7ROSqap0lc27kmvM5L0TkQQB/z3HqIICPI3l9X9XMUaee+nXhk1FV/7uqHgRwUEQOALgDwCcA/BxAWVVfEJEbAPypiOwpmIyFmsfqNQcBXAIwUz1XqHk0ZCzcPDoo3Dy6mhmyNEXGVYP7n0WHAHyqOu6nAPwBkn8OcpHNx4ZUDKo67KoXkTcgsdk9VvUz7QLwQxG5EYlWvrbm8l0ATlfrdznqmyKjg68C+HMAn1DVlwC8VG3/AxF5CsD1RZIRBZtHEfkAgH8MYKhqMkDR5tElIwo2j0abQs2jQa7z2MhaZRaRewB8s3poyZwfeTo0ilYAzOKy83kP6h0+T+Oyk+pvAbwVl51UI02W63U1v38UwB9Vf99eI9NrADwDYFvBZCzSPN4C4DiA7Q31RZpHS8bCzGONTA+h3rFbmHn0yFi4eayR7dU1v/87JH4Fr8y5yZbnYEUrqFEM1eODSFYAPImaFQoA9gL4afXcF1DdMd5Euf64Ot6PAfwZgJ3V+t8C8Hj1S/NDAP+kaDIWbB5PILHVPlotKyunijSPThkLNo+/ieS/2JcA/ALAAwWcR6eMRZpHh8z/FcBPqn9D30C9onDKnFdhSAxCCCF1bLpVSYQQQvxQMRBCCKmDioEQQkgdVAyEEELqoGIghBBSBxUDIYSQOqgYCCGE1EHFQEhkRORfi8jPa+LsP1oNx0JIW8ANboRERkS+COCHqnpvq2UhJAS+MRASnzcgCW9BSFvCNwZCIiMiLyAJKLdcrZpU1akWikRIJjZk2G1CWoWIXAvgWVV9Y6tlISQUmpIIicsbATzRaiEIWQ9UDITE5Q2gYiBtDn0MhERERGYAvBPA89UqBfAOVX2xdVIRkg0qBkIIIXXQlEQIIaQOKgZCCCF1UDEQQgipg4qBEEJIHVQMhBBC6qBiIIQQUgcVAyGEkDr+Pyh8L37l3OddAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "high_T = points_full[:,0]>Tc\n",
    "low_T = points_full[:,0]<Tc\n",
    "\n",
    "E_M_high = points_full[high_T][:,1:]\n",
    "E_M_low = points_full[low_T][:,1:]\n",
    "\n",
    "fig, ax = plt.subplots(1,1)\n",
    "ax.scatter(E_M_high[:,0],E_M_high[:,1],c='r')\n",
    "ax.scatter(E_M_low[:,0],E_M_low[:,1],c='b')\n",
    "ax.set_xlabel(\"$E$\")\n",
    "ax.set_ylabel(\"$|m|$\");"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 1: Single neuron binary classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create a binary classifier that can take $(E,|m|)$ as input data and predict a binary label (0=below Tc, 1=above Tc)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Hints:\n",
    "* Build your own binary classifier from a single neuron. Study the lecture notes and the exercise on logistic regression / neural networks.\n",
    "* Normalize the data before training / testing (mean=0, standard deviation=1).\n",
    "* Split into 70 % training data and 30% test data.\n",
    "* Use weight decay alpha=1.0, learning parameter eta=0.01\n",
    "* A rather large number of training iterations will be needed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**How well does it perform? Plot the decision boundary.**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Neural network classifiers based on spin configurations"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should now use a neural network to learn the phase transition based on images of the spin configurations.\n",
    "We will use special tricks to generate new data for free from the configurations that we have generated (e.g. by flipping all spins, and by mirroring)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Generate spin configurations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      ">>> Loaded from ising_config_data_big.pkl\n"
     ]
    }
   ],
   "source": [
    "# Use 6 threads to run some lattice evolutions\n",
    "# With current parameters, this took a minute or two on an i7\n",
    "# Dump out the lattice configurations/temperatures to a file\n",
    "# and just load it if the file already exists, since this\n",
    "# is a CPU intensive cell.\n",
    "\n",
    "fname = \"ising_config_data_big.pkl\"\n",
    "\n",
    "def get_lattices(T, N=10, Nlattices=25, Nthermal=200):\n",
    "    \"\"\"\n",
    "    Generates a set of lattices at a given temperature\n",
    "    \n",
    "    Args:\n",
    "        T (float): temperature of the lattice\n",
    "        N (int): the lattice size\n",
    "        Nlattices (int): the number of lattices to generate\n",
    "        Nthermal (int): the number of steps to simulate thermalization\n",
    "    \"\"\"\n",
    "    lattices = []\n",
    "    for _ in range(Nlattices):\n",
    "        lat = Lattice(N=N,T=T)\n",
    "        for _ in range(Nthermal):\n",
    "            lat.step()\n",
    "        lattices.append(lat.lattice)\n",
    "    return round(T, 4), lattices\n",
    "\n",
    "if not os.path.exists(fname):\n",
    "    pool = ThreadPool(6)\n",
    "    Ts = np.arange(5.0,1.0,-0.1)\n",
    "    d_data = {}\n",
    "    for T, lattices in pool.imap_unordered(get_lattices, Ts):\n",
    "        print(T)\n",
    "        d_data[T] = lattices\n",
    "\n",
    "    with open(fname,\"wb\") as fhout:\n",
    "        pickle.dump(d_data,fhout)\n",
    "    print(f\">>> Dumped to {fname}\")\n",
    "\n",
    "else:\n",
    "    d_data = pickle.load(open(fname,\"rb\"))\n",
    "    print(\">>> Loaded from {}\".format(str(fname)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "# make vector of input matrices, vector of temperatures\n",
    "X_data = []\n",
    "y_data = []\n",
    "T_data = []\n",
    "\n",
    "for T,configs in d_data.items():\n",
    "    for config in configs:\n",
    "        # flip spins to double dataset keeping E same\n",
    "        # this is also needed so ML algorithm doesn't learn\n",
    "        # to prefer one magnetization sign over another\n",
    "        # also make truth labels (0 is low T phase, 1 is high T phase)\n",
    "        # and also mirror lattice horizontally/vertically to get free data\n",
    "        target = 0\n",
    "        if T > Tc:\n",
    "            target = 1\n",
    "            \n",
    "        X_data.append(config)\n",
    "        X_data.append(np.flip(config,0))\n",
    "        X_data.append(np.flip(config,1))\n",
    "        X_data.append(-config)\n",
    "        X_data.append(-np.flip(config,0))\n",
    "        X_data.append(-np.flip(config,1))\n",
    "        \n",
    "        for _ in range(6):\n",
    "            T_data.append(T)\n",
    "            y_data.append(target)\n",
    "\n",
    "\n",
    "X_data = np.array(X_data)\n",
    "y_data = np.array(y_data)\n",
    "T_data = np.array(T_data)\n",
    "\n",
    "# convert spin matrices from -1,1 to 0,1\n",
    "X_data = 0.5*(X_data+1)\n",
    "\n",
    "# Thus, our training and test sets will consist of \n",
    "# the lattice images and the targets will be 0 or 1.\n",
    "# If the lattice is at low (T<Tc) or high (T>Tc) temperature\n",
    "# It's up to the NN / CNN to learn the concept of temperature/magnetization/etc."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We implement two different ways to split our data in training and test sets:\n",
    "\n",
    "1. Randomly\n",
    "1. Using only spin configurations at very high and very low temperatures as training data\n",
    "\n",
    "In the second approach our ultimate goal will be to see if we can correctly predict the phase for all intermediate temperatures. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "# function to split the data\n",
    "def split_data(X_data, y_data, T_data, Thi=4.7, Tlo=1.3, test_size=0,random_state=42):\n",
    "    \"\"\"\n",
    "    Splits the data into train and test sets according to the settings.\n",
    "    \n",
    "    Either using sklearn train_test_split (test_size).  \n",
    "    Or, alternatively, using T>Thi and T<Tlo data as training data and Tlo < T < Thi as predictions/test.\n",
    "    \"\"\"\n",
    "    if test_size:\n",
    "        print(f\"Using sklearn to split data. Test size: {test_size*100:}%\")\n",
    "        X_train, X_test, y_train, y_test, T_train, T_test = \\\n",
    "          train_test_split(X_data, y_data, T_data, test_size=test_size, random_state=random_state)\n",
    "    elif (Thi and Tlo):\n",
    "        print(f\"Using T>{Thi} and T<{Tlo} as training data.\")\n",
    "        train_set = np.logical_or(T_data>=Thi,T_data<=Tlo)\n",
    "        test_set = np.logical_not(train_set)\n",
    "        X_train = X_data[train_set]\n",
    "        X_test = X_data[test_set]\n",
    "        y_train = y_data[train_set]\n",
    "        y_test = y_data[test_set]\n",
    "        T_train = T_data[train_set]\n",
    "        T_test = T_data[test_set]\n",
    "    else:\n",
    "        print(\"No rule to split data.\")\n",
    "        return None\n",
    "        \n",
    "    print(f\"...Training samples: {X_train.shape[0]}\")\n",
    "    print(f\"...Testing samples:  {X_test.shape[0]}\")\n",
    "    return(X_train, X_test, y_train, y_test, T_train, T_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 2: Tensorflow neural network with a single hidden layer\n",
    "\n",
    "1. Implement a neural network with a single hidden layer using tensorflow.\n",
    "1. Flatten the 10x10 input images into arrays of length 100.\n",
    "1. Use random splitting of the data (70% for training)\n",
    "\n",
    "Study in particular:\n",
    "* What is the accuracy of predictions on the training / test data?\n",
    "* Try to predict the critical temperature\n",
    "\n",
    "*Hint*: Plot the average predictions for lattices in the test set at different temperatures. The point at which the prediction passes through 0.5 is the critical temperature. I.e., at a phase transition, the CNN is confused about what phase the system is in. So we fit a sigmoid to the points and use this to estimate the crossing point (`Tc_fit`)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using sklearn to split data. Test size: 30.0%\n",
      "...Training samples: 4200\n",
      "...Testing samples:  1800\n"
     ]
    }
   ],
   "source": [
    "X_train, X_test, y_train, y_test, T_train, T_test = split_data(X_data, y_data, T_data,test_size=0.3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### A note about the output encoding, loss function and number of outputs\n",
    "\n",
    "If you check the target `y_train`, you notice that we have an `array([0, 1, 1, ..., 0, 0, 1])` representing the phases. This represents a binary classification problem where we expect the output of the final layer of the network to give a 1/0 output. In particular, if the final neuron has a sigmoid activation, it will give outputs very close to 0/1 as `array([0.01, 0.99, 0.95, ..., 0.2, 0.1, 0.86])`.\n",
    "\n",
    "While this works for a binary classification problem, if you have more than two classes as output, e.g., {\"red\", \"green\", \"yellow\"}, this simple 1/0 encoding does not work and you need to have a one-hot encoding:\n",
    "\n",
    "red = [1, 0, 0]\n",
    "green = [0, 1, 0]\n",
    "yellow = [0, 0, 1]\n",
    "\n",
    "In this case, think about how many neurons you should have in your output layer. You need to be also careful about the loss function you choose in this case."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### One hot encoding with keras\n",
    "\n",
    "If you do choose a one-hot encoding then you can use the following functionality in `keras`\n",
    "\n",
    "```\n",
    "keras.utils.to_categorical(y_train, 2)\n",
    "```\n",
    "\n",
    "where the 2 represents that you have two classes in your y_train."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "y_train = keras.utils.to_categorical(y_train, 2)\n",
    "y_test = keras.utils.to_categorical(y_test, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 3: Convolutional neural network (CNN) \n",
    "Now use a convolutional neural network (CNN) instead. Those are optimal for working with images.\n",
    "\n",
    "1. Implement a CNN with several hidden layers using tensorflow.\n",
    "1. Now we will use only low- and high-temperature configurations as training data.\n",
    "\n",
    "Study in particular:\n",
    "* What is the accuracy of predictions on the training / test data?\n",
    "* Try to predict the critical temperature\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Using T>4.7 and T<1.3 as training data.\n",
      "...Training samples: 1050\n",
      "...Testing samples:  4950\n"
     ]
    }
   ],
   "source": [
    "X_train, X_test, y_train, y_test, T_train, T_test = split_data(X_data, y_data, T_data, Thi=4.7, Tlo=1.3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "# prepare for CNN input\n",
    "img_rows, img_cols = 10, 10\n",
    "X_train = X_train.reshape(X_train.shape[0], img_rows, img_cols, 1)\n",
    "X_test = X_test.reshape(X_test.shape[0], img_rows, img_cols, 1)\n",
    "\n",
    "input_shape = (img_rows, img_cols, 1)\n",
    "\n",
    "y_train = keras.utils.to_categorical(y_train, 2)\n",
    "y_test = keras.utils.to_categorical(y_test, 2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 4: Understanding the neural network (extra)\n",
    "Description to be added."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Task 5: A simple Bayesian neural network (extra)\n",
    "Description to be added."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}