2d_ising_model_mcmc.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# <a id='the top'>Monte carlo simulation of 2D Ising model</a>\n",
    "\n",
    "In this notebook we will try to understand the physics of <a href='https://en.wikipedia.org/wiki/Square_lattice_Ising_model'>two dimensional Ising model</a> numerically. The Hamiltonian is given by a claasical variable called spin $s_i\\in\\{\\pm1\\}$, these degrees of freedom sits on the vertex of the square lattice.\n",
    "\n",
    "$$ \\mathcal{H}=-\\sum_{i,j}s_is_j $$\n",
    "\n",
    "We will try to find the ground state configuration of this system using the Markov chain monte carlo (<a href='https://en.wikipedia.org/wiki/Markov_chain_Monte_Carlo'>MCMC</a>) method. Notice that this problem has the complexity of $2^{N^2}$, since in a $(N\\times N)$ lattice there are $N^2$ spins, and every spins has a two dimentional Hilbert space. \n",
    "\n",
    "Here we will show that to find the ground state we don't need to look at the full Hilbert space. For a finite temperature, we will let the system evolve it's energy landscape using the monte carlo dynamics, and at the end we will see which state the system is settling into, that will be the best candidate for the ground state (GS).\n",
    "\n",
    "Of course this problem has a trivial ground state, actually a two fold degenerate GS with all the $s_{i}=\\pm 1$. Still this problem is iteresting on its own, since there is a prediction of the critical temperature ($T_c\\sim 2.27$), using <a href='https://en.wikipedia.org/wiki/Kramers%E2%80%93Wannier_duality'>Krammers-wannier duality</a>.\n",
    "\n",
    "Also the problem of Ising model is usefull to study numerous complex problems such as in <a href='https://en.wikipedia.org/wiki/Ising_model#Applications'>physical systems</a>, <a href='https://aapt.scitation.org/doi/abs/10.1119/1.2779882?journalCode=ajp'>socio economic problems</a>, and <a href='https://en.wikipedia.org/wiki/Monte_Carlo_methods_in_finance'>finance</a>."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "def initial_state(N):\n",
    "    return 2*np.random.randint(2,size=(N,N))-1\n",
    "\n",
    "def mc_move(state,T):\n",
    "    '''given a current state, take decission to go to the next state'''\n",
    "    N=state.shape[0]\n",
    "    a=np.random.randint(N)\n",
    "    b=np.random.randint(N)\n",
    "    s=state[a,b]\n",
    "    #print(state,a,b)\n",
    "    # adding all the spins in the neighbour state\n",
    "    field = state[(a+1)%N,b] + state[a,(b+1)%N] + state[(a-1)%N,b] + state[a,(b-1)%N]\n",
    "    # change in enrgy if we flip the spin at (a,b)\n",
    "    del_energy=2*s*field\n",
    "    if del_energy<0:\n",
    "        state[a,b]=-s\n",
    "    elif np.random.rand()<np.exp(-del_energy*(1./T)):\n",
    "        state[a,b]=-s\n",
    "    return state\n",
    "\n",
    "def calculate_energy(state):\n",
    "    '''calculate energy of a state'''\n",
    "    energy=0\n",
    "    N=state.shape[0]\n",
    "    for i in range(N):\n",
    "        for j in range(N):\n",
    "            increament = state[i,j]*(state[(i+1)%N,j] + state[i,(j+1)%N])\n",
    "            energy = energy + increament\n",
    "    return -energy\n",
    "\n",
    "def calculate_magnet(state):\n",
    "    '''Magnetization of a given configuration'''\n",
    "    N=state.shape[0]\n",
    "    magnetization = np.sum(state)\n",
    "    return magnetization"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "data=32\n",
    "temp_arr=np.linspace(1,4,32)\n",
    "N=8\n",
    "mc_equ=2**14                     # no of iteration to reach the groud state, i.e. equilibration time\n",
    "mc_iter=2*16                    # no of iteration to calcualte the expectation values one gs is reached\n",
    "\n",
    "energy, magnet, spec_heat, suscept = np.zeros(data), np.zeros(data), np.zeros(data), np.zeros(data)\n",
    "loop=16"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.random.seed(8003)\n",
    "\n",
    "for i,T in enumerate(temp_arr):\n",
    "    e,m=0.,0.\n",
    "    state = initial_state(N)\n",
    "    #print(T)\n",
    "    for k in range(loop):\n",
    "        for j in range(mc_equ):\n",
    "            mc_move(state,T)\n",
    "        #\n",
    "        for j in range(mc_iter):\n",
    "            mc_move(state,T)\n",
    "            en = calculate_energy(state)\n",
    "            mag = calculate_magnet(state)\n",
    "            #\n",
    "            e += en\n",
    "            m += mag\n",
    "    #    \n",
    "    energy[i] = e*(1./(2*N*N*loop*mc_iter))\n",
    "    magnet[i] = m*(1./(N*N*loop*mc_iter))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAgMAAAJWCAYAAADFmIQTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nOzde5wddX3w8c83gNpuVkBBiIE1atJWtBZ1i7eoqKAYrahV10sfoY1Pan1sa7Et+KDWy2NLsbZiay8xVcF6ibdKqlgFvMaKkoBVkdKkFCEkclNxs9QL5Pv8MbPu5mT2cs7OOWd2z+f9eu1r9sz8Zua7h0Pme37XyEwkSdLgWtbvACRJUn+ZDEiSNOBMBiRJGnAmA5IkDTiTAUmSBpzJgCRJA+7gfgfQD0cccUSuWrWq32FIktQz27dvvzUzj6w61thkICLuBWwGVgHXAc/PzO+3lDke+DvgnsBdwJszc/Nc1161ahXbtm2rO2RJkhorIr4z07EmNxOcBVyamWuAS8vXre4AXpKZDwZOAd4WEYf1MEZJkha9JicDpwLnl7+fDzyrtUBm/mdm7ih/3w3cDFRWgUiSpGpNTgaOysw9AOX2PrMVjogTgLsB/9WD2CRJWjL62mcgIi4Bjq44dHab11kBvBc4LTP3zVBmA7ABYGRkpM1IJUlauvqaDGTmSTMdi4ibImJFZu4pH/Y3z1DunsAngddk5mWz3GsjsBFgdHTU1ZkkSSo1uZlgC3Ba+ftpwIWtBSLibsA/Axdk5od7GJskSUtGk5OBc4CTI2IHcHL5mogYjYhNZZnnA48HTo+Ir5c/x/cnXEmSFqfIHLwa89HR0XSeAUnSIImI7Zk5WnWsyTUDkiSpBxo7A6EkSYNkfBw2b4YdO2DNGhgbg+Hh3tzbZECSpD7buhXWrYN9+2BiAoaG4Iwz4KKLYO3a7t/fZgJJkvpofLxIBMbHi0QAiu3k/r17ux+DyYAkSX20eXNRI1Bl377ieLeZDEiS1Ec7dkzVCLSamICdO7sfg8mAJEl9tGZN0UegytAQrF7d/RhMBiRJ6qOxMVg2w9N42bLieLeZDEiS1EfDw8WogeHhqRqCoaGp/cuXdz8GhxZKktRna9fC7t1FZ8GdO4umgbGx3iQCYDIgSVIjLF8O69f35942E0iSNOBMBiRJGnA2E0iSBkY/5/9vMpMBSdJA6Pf8/01mM4Ekaclrwvz/TWYyIEla8pow/3+TmQxIkpa8Jsz/32QmA5KkJa8J8/83mcmAJGnJa8L8/03W2GQgIu4VERdHxI5ye3hFmftFxPaI+HpEXBURL+tHrJKkZlvI/P/j47BpE5x5ZrEdH+9NzL0UmdnvGCpFxLnA9zLznIg4Czg8M89sKXM3ir/hxxGxHPgW8JjM3D3btUdHR3Pbtm1di12S1Ex797Y3/3/VcMRlyxbncMSI2J6Zo5XHGpwMXAOcmJl7ImIF8PnM/MVZyt8buBJ4lMmAJGmhxsdh5crqmoDh4WJhoV4tJFSH2ZKBxjYTAEdl5h6AcnufqkIRcWxEfAO4AfjzuRIBSZLmY5CGI/Z1BsKIuAQ4uuLQ2fO9RmbeADw0Iu4LfDwiPpKZN1XcawOwAWBkZKTDiCVJg2KQhiP2NRnIzJNmOhYRN0XEimnNBDfPca3dEXEV8DjgIxXHNwIboWgmWFjkkqSlbnI4YlVCsNSGIza5mWALcFr5+2nAha0FIuKYiPi58vfDgccC1/QsQklS33S7l/8gDUds8kJF5wAfioj1wPXA8wAiYhR4WWa+FHgQ8NaISCCAv8jMb/YrYElSb/Ri0aHJYYczjSZYTJ0H59LY0QTd5GgCSVq8et3Lv93hiJMxNm2p5NlGEzS5ZkCSpAPMp5f/+vX13W/58vautxiXSm5ynwFJkg7Q5F7+i3WpZJMBSdKi0uRFhxbr3AQmA5KkRaXJvfybXGsxG/sMSJL6rp0Od03u5b9Y5yZwNIEkqa86XQyok17+3dbk9QwcTSBJaqTpHe4mTX6rXrdu9odnu738e6HJtRazMRmQJPVNr4cJ9sLatUUS07Rai9mYDEiS+maxdribSxNrLWbjaAJJUt80eZjgIDEZkCT1TZOHCQ4SmwkkSTPqZI79pTJMcJA4tFCSVKmTIX9LaZjgUjPb0EKTAUnSAToZL9/kMfaaPRmwz4Ak6QCdzLG/WOfll8mAJKlCJ0P+luowwUFgMiBJOkAnQ/4cJrh4mQxIkg7QyZA/hwkuXiYDkrQIjY/Dpk1w5pnFtqrT3kJMDvkbHp76tj80NLW/qiNgJ+eoGRo7miAi7gVsBlYB1wHPz8zvz1D2nsDVwD9n5ivmurajCSQtZp0O3+tkzoBOhvw5TLCZFuXQwog4F/heZp4TEWcBh2fmmTOUPQ84sixvMiBpyep0+F6nCYSWjsU6tPBU4Pzy9/OBZ1UViohHAEcBn+lRXJIGRLer4jvRyfC96csET/b2n5iY2r93b/fi1eLQ5GTgqMzcA1Bu79NaICKWAW8F/qjHsUla4rZuLb6Bv/KVcO65xXblymJ/P3UyfM/x/5pLX9cmiIhLgKMrDp09z0u8HLgoM2+IiLnutQHYADAyMtJOmJIGzPRv0pMmH8Dr1vV3Jr3J4XtVCcFMw/cc/6+59DUZyMyTZjoWETdFxIrM3BMRK4CbK4o9GnhcRLwcWA7cLSL2ZuZZFffaCGyEos9APX+BpKVoPt+k61yrvp2OfWNjcMYZ1cdmGr7XSQKhwdLkZoItwGnl76cBF7YWyMwXZ+ZIZq4C/hC4oCoRkKR29PKbdLvNEZ0M33P8v+bS5CWMzwE+FBHrgeuB5wFExCjwssx8aT+Dk7R09eqbdKfNEWvXFsfmO3zPZYI1l8YOLewmhxZKmk2vVt/btKmoCZgp6TjvvHqbIxz/P9hmG1rY5JoBSeqLXn2T7nXHvuXL600utHSYDEhShXar4jthxz41hc0EktQnvWqOkMBmAknqmXaGCdqxT01hMiBJNama//+MM2af/78XzRHSXGwmkKQaWOWvplusCxVJ0qLh/P9azEwGJKkGzv+vxcxkQJJqMDlMsIrDBNV0JgOSVAPn/9diZjIgSTXoZAEhqSkcWihJNXGYoBYrkwFJqpHz/2sxsplAkqQBZzIgSdKAMxmQJGnAmQxIkjTg7EAoaSC0s5qgNGhMBiQteZ2sJigNksY2E0TEvSLi4ojYUW4Pn6HcXRHx9fJnS6/jlNRs4+NFIjA+PrV2wMTE1P69e/sbn9QEjU0GgLOASzNzDXBp+brK/2Tm8eXPM3sXnqTFwNUEpbk1ORk4FTi//P184Fl9jEXSIuVqgtLcmpwMHJWZewDK7X1mKHePiNgWEZdFhAmDpP24mqA0t752IIyIS4CjKw6d3cZlRjJzd0Q8APhsRHwzM/+r4l4bgA0AIyMjHcUrafEZGys6C1ZxNUGp0NdkIDNPmulYRNwUESsyc09ErABunuEau8vttRHxeeBhwAHJQGZuBDYCjI6OZg3hS1oEJlcNbB1NsGyZqwlKk+adDETEIZn5024G02ILcBpwTrm9sCKmw4E7MvPHEXEE8Fjg3B7GKGkRcDVBaXbt1AzcGBHvBt6Zmb3ocnMO8KGIWA9cDzwPICJGgZdl5kuBBwH/EBH7KPo/nJOZ3+5BbJL6qJMJhFxNUJpZZM6vxjwibgXuBSTwWeDvgQsz887uhdcdo6OjuW3btn6HIakDVRMITVb5O4GQNLOI2J6Zo1XH2hlNcF/gN4AvAU8GPgTcEBFvjoj7LzxMSZqdEwhJ3THvZCAzf5KZ78/ME4FfAt5G0czwamBHRFwUEadGRJOHK0paxJxASOqOjh7cmfmfmfkqYCVTtQWnAB8Dro+I10fEfesLU5KcQEjqlgV9i8/MnwCfBP4Z2A0ERXPC64D/joi3RcTdFxylJOEEQlK3dJwMRMSjytEFu4G/AoaAtwPHA78FXAP8LkVzgiQt2NhY0VmwihMISZ1rKxmIiOGIeHlE/DvwZYrx/1dTzOx338x8ZWZ+IzPfQzH5z2eB59Ycs6QBNTmB0PDwVA3B0NDUfucNkDrTzqRDm4Ax4OeBHwPvBf42M79WVT4z7ypnBHxSDXFKWqLanTPACYSk+rUzz8A+iml+/x54d2Z+bx7nPBY4KTPfsKAoa+Y8A1IzOGeA1DuzzTPQzgyET8vMT7dz48z8MkVzgiTtZ/qcAZMmRwqsW1d8+/fbvtQb7cwz0FYiIEmzcc4AqTmcIEhSXzhngNQc804GImJfRNw1x8+dEfG9iPhyRLzKOQYkzcQ5A6TmaKdm4IvANygmFtoHfAf4WrndV+7/JnAjxbDCc4EvR8QM/7tLGmTOGSA1RzvJwAuBQ4EPAg/MzAdk5qMz8wHAA8v99wROBo4C3gU8HPjjekOWtBQ4Z4DUHO0MLbwAOG6mYQllmW3AVZl5WkQcBFwF3JmZD6kl2po4tFBqjr17nTNA6oW6hhY+leLb/mwuppiKeHLSoS8CL27jHpIGzPLlsH59v6OQBls7zQTDFM0Aszm0LDdpzomJJElSf7WTDPwHMDbT0sQRcQzFdMVXT9t9LHBb5+FJkqRua6eZ4K0U6xFcERF/TTGz4E0UnQXXUqxQeBjwlwARcTBwEvClOgOWJEn1mncykJnvK2sF3gy8seVwAHcCZ2fm+8p9hwGvA75aR6CSJKk72qkZIDPfEhEfoegUeDxFH4EfAlcC78/Ma6eVvRX4h04Di4h7AZuBVcB1wPMz8/sV5UaATRRNEgmsy8zrOr2vJEmDpq1kACAz/xv4f12IpdVZwKWZeU5EnFW+PrOi3AXAmzPz4ohYTjEBkiRJmqd2piO+KyLe381gWpwKnF/+fj7wrIqYjgMOzsyLATJzb2be0bsQJUla/NoZTTBOMfVwrxyVmXsAyu19Ksr8AvCDiPhYRFwZEW8pJzuSJEnz1E4zwZXAcXXePCIuAY6uOHT2PC9xMPA4irUQrqfoY3A68I8V99oAbAAYGRnpIFpJkpamdpKBPwf+JSJOnqyWX6jMPGmmYxFxU0SsyMw9EbECuLmi2C7gysmOixHxceBRVCQDmbkR2AjFdMR1xC9J0lLQTjJwH+BfgU+VD93Lge9S9ODfT2ZeUENsW4DTgHPK7YUVZS4HDo+IIzPzFuBJgIsOSH0yPl6sM7BjR7FE8dhYsfCQpGZrZ6GifRQP/mg5NP0CAWRmLrjdPiLuDXwIGKFoAnheZn4vIkaBl2XmS8tyJ1NMiBTAdmBDZv5ktmu7UJFUv61bYd062LcPJiaKFQiXLStWIFy7tt/RSZptoaJ2koHT5nvDzDx/7lL9YzIg1Wt8HFauLLathodh925XIpT6rZZVC5v+gJfUP5s3FzUCVfbtK467MqHUXO0MLZSkSjt2FE0DVSYmYOfO3sYjqT1tz0AYEUcCvw48CBia1nZ/JHB/4JuZ+T+1Rimp0dasKfoIVCUEQ0OwenXvY5I0f23VDETEeop1At5BsUrhb047fBTwFeBFdQUnaXEYGys6C1ZZtqw4Lqm52pmO+GSKcfr/CTwb+LvpxzPzW8BVVEwbLGlpGx4uRg0MDxc1AVBsJ/fbeVBqtnaaCc4E9gBPyMwfRsTDKsp8A3h0LZFJWlTWri1GDWzeXPQRWL26qBEwEZCar51kYBT4YGb+cJYyu6ieXljSAFi+3FED0mLUTp+BuwEz9Bf+mcOAuzoPR5Ik9Vo7ycB1wCPmKPNI4JqOo5EkST3XTjJwIfC4iHhe1cGI+E3gocBH6whMkiT1Rjt9Bs4FXgB8ICKeCxwKEBGvoFhG+DnADuCv6w5SkiR1TzvTEX8/Ip4AXABMrx14e7n9EvCizJyrX4EkSWqQtmYgzMzrgRMj4qEUQwjvDdwOXJaZ27sQnyRJ6rK2pyMGyMxvUMwpIEmSFjkXKpIkacC1VTMQEYcApwInAIcDB1UUy8x02hFJkhaJeScDEXFf4GLgl4CYpWgCJgOSJC0S7dQMvJVi2eIPAO8EbgDu7EZQkiSpd9pJBp4CfDEzX9ytYCRJUu+104HwHsBXuxVIq4i4V0RcHBE7yu3hFWWeGBFfn/bzo4hwCWVJktrQTjLwLeB+3QqkwlnApZm5Bri0fL2fzPxcZh6fmccDTwLuAD7TwxglSVr02kkG3gI8MyKO61YwLU4Fzi9/Px+Y6xv/c4FPZeYdXY1KkqQlpp0+AzcD/wL8W0ScB2wHflBVMDO/WENsR2XmnvJ6eyLiPnOUfwHwlzXcV5KkgdJOMvB5imGDAby2/H0mVfMPHCAiLgGOrjh0dhtxERErgF8GPj1LmQ3ABoCRkZF2Li9J0pLWTjLwRmZPANqWmSfNdCwiboqIFWWtwAqKmomZPB/458z86Sz32ghsBBgdHa3175CabnwcNm+GHTtgzRoYG4Ph4X5HJakp2lm18PVdjKPKFuA04Jxye+EsZV8IvLoXQUmLzdatsG4d7NsHExMwNARnnAEXXQRr1/Y7OklN0OS1Cc4BTo6IHcDJ5WsiYjQiNk0WiohVwLHAF/oQo9Ro4+NFIjA+XiQCUGwn9+/d29/4JDXDrMlARDw+IubdwB4RvxIRL1l4WJCZt2XmkzNzTbn9Xrl/W2a+dFq56zJzZWbuq+O+0lKyeXNRI1Bl377iuCTNVTPwOeD06Tsi4syIuG2G8s8C3l1DXJJqsGPHVI1Aq4kJ2Lmzt/FIaqa5koGqBYnuARzWhVgk1WzNmqKPQJWhIVi9urfxSGqmJvcZkLRAY2OwbIb/y5ctK45LksmAtIQNDxejBoaHp2oIhoam9i9f3t/4JDVDO/MMSFqE1q6F3buLzoI7dxZNA2NjJgKSppgMSANg+XJYv77fUUhqqvk0EzhbnyRJS9h8agZeHxGvb90ZEXfVH46kpnAKY2lwzCcZqBpeOBtrEqRFzimMpcEyazKQmY42kAbM9CmMJ01OXLRuXdEZ0c6H0tLiw17SfpzCWBo8JgOS9uMUxtLgMRmQtB+nMJYGj8mApP04hbE0eEwGJO3HKYylweMMhJIO4BTG0mAxGZBUySmMpcFhM4EkSQPOZECSpAHX2GQgIu4VERdHxI5ye/gM5c6NiKsi4uqIeHtEtDt9siRJA62xyQBwFnBpZq4BLi1f7yciHgM8Fngo8BDgV4En9DJISZIWuyYnA6cC55e/nw88q6JMAvcA7gbcHTgEuKkn0UmStEQ0ORk4KjP3AJTb+7QWyMyvAJ8D9pQ/n87Mq3sapSRJi1xfhxZGxCXA0RWHzp7n+auBBwHHlLsujojHZ+YXK8puADYAjIyMdBawJElLUF+Tgcw8aaZjEXFTRKzIzD0RsQK4uaLYs4HLMnNvec6ngEcBByQDmbkR2AgwOjqadcQvSdJS0ORmgi3AaeXvpwEXVpS5HnhCRBwcEYdQdB60mUCSpDY0eQbCc4APRcR6iof+8wAiYhR4WWa+FPgI8CTgmxSdCf81M/+lT/FKPTE+XkwTvGNHscLg2FixboAkdSoyB6/GfHR0NLdt29bvMKS2bd0K69bBvn0wMVEsILRsWbGA0Nq1/Y5OUpNFxPbMHK061uRmAknTjI8XicD4eJEIQLGd3L93b3/jk7R4mQxIi8TmzUWNQJV9+4rjktQJkwFpkdixY6pGoNXERLHUsCR1wmRAWiTWrCn6CFQZGoLVq3sbj6Slw2RAWiTGxorOglWWLSuOS1InTAakRWJ4uBg1MDw8VUMwNDS1f/ny/sYnafFq8jwDklqsXQu7dxedBXfuLJoGxsZMBCQtjMmAtMgsXw7r1/c7CklLic0EkiQNOJMBSZIGnM0EUk1cM0DSYmUyINWgas2AM86Ye80AEwhJTeBCRdICjY/DypXFttXwcNH7v6q3v4sOSeolFyqSuqiTNQNcdEhSk5gMSAvUyZoBLjokqUlMBqQF6mTNABcdktQkJgPSAnWyZoCLDklqEpMBaYE6WTPARYckNYlDC6UatLtmwGSiMNNoAtcakNRLJgNSTdpdM8BFhyQ1RWOTgYi4F7AZWAVcBzw/M79fUe7PgaeXL9+UmfbD1qLhokOSmqDJfQbOAi7NzDXApeXr/UTE04GHA8cDjwT+KCLu2dMoJUla5JqcDJwKnF/+fj7wrIoyxwFfyMw7M3MC+HfglB7FJ0nSktDkZOCozNwDUG7vU1Hm34GnRcTPR8QRwBOBY3sYoyRJi15f+wxExCXA0RWHzp7P+Zn5mYj4VeDfgFuArwB3znCvDcAGgJGRkY7i1eLkYkCSNLvGLlQUEdcAJ2bmnohYAXw+M39xjnPeD/xTZl40WzkXKhocLgYkSYXZFipq7GgCYAtwGnBOub2wtUBEHAQclpm3RcRDgYcCn+lplGqs6YsBTZqcAnjduplXE5w819oESYOiycnAOcCHImI9cD3wPICIGAVelpkvBQ4BvhQRAD8EfiMzK5sJNHjmsxhQ1bC+qtqEM86wNkHS0tXYZCAzbwOeXLF/G/DS8vcfUYwokA7QyWJAC6lNkKTFqsmjCaQF6WQxIJcWljSITAa0ZHWyGJBLC0saRI1tJpCqtNOxr5PFgCZrE6oSApcWlrRUNXZoYTc5tHBx6nSY4N69818MaHwcVq7cv8/ApOFh+wxIWrwW69BC6WcW0rGvncWAXFpY0iAyGVDftFPl3+kwwU64tLCkQWMyoL5odyx/rzv2ubSwpEHiaAL13PQq/8kH/MTE1P69ew88p5NhgpKk+TEZUM91Mpa/k2GCkqT5MRlQz3VS5T/ZsW94eKqGYGhoar/t+ZLUOfsMqOc6Hctvxz5J6g7nGVDPOZZfknpvtnkGbCZQz1nlL0nNYjOB+sIqf0lqDpMB9Y1j+SWpGUwGdIB2ZgZcyDmSpGawA+ES1+5DupPFgDpdQEiS1DuzdSA0GVjC2n1Id9LL35EBkrQ4OJqggcbHYdMmOPPMYlv1MF3IOZ1M+dvJzICdnCNJapbG9hmIiOcBrwceBJyQmZVf5SPiFOA84CBgU2ae07MgS3VUxc+2SE8n53Syyl8nMwP2egEhSVL9GpsMAN8CngP8w0wFIuIg4B3AycAu4PKI2JKZ3+5NiO0/pKd/Y580+TBdt27mqvh2z+nkId3JzICdziYoSWqOxjYTZObVmXnNHMVOAHZm5rWZ+RPgg8Cp3Y+u0OSq+E5W+etkMSAXEJKkxa+xycA8rQRumPZ6V7mvJzp5SPeqKr6Th3QnMwM6m6AkLX59bSaIiEuAoysOnZ2ZF87nEhX7KodHRMQGYAPAyMjIvGOcTZOr4icfxjONJpjpId3JzIDOJihJi1vjhxZGxOeBP6zqQBgRjwZen5lPLV+/GiAz/2y2a9Y1tHDTJnjlK2d+SJ933oGd9Ho9fG/vXh/SkqTZhxY2uQPhfFwOrImI+wM3Ai8AXtSrm4+NFZ0Fq8xVFd/ON/ZOv+WDU/5KkubW2GQgIp4N/DVwJPDJiPh6Zj41Iu5LMYRwXWbeGRGvAD5NMbTwXZl5Va9itCpekrQUNL6ZoBvqnoHQqnhJUtMt5WaCRrAqXpK0mC32oYWSJGmBTAYkSRpwJgOSJA04kwFJkgacyYAkSQPOZECSpAFnMiBJ0oAzGZAkacCZDEiSNOAGcjriiLgF+E4XLn0EcGsXrrsY+V5M8b3Yn+/HFN+LKb4XU7r1XtwvM4+sOjCQyUC3RMS2meZ9HjS+F1N8L/bn+zHF92KK78WUfrwXNhNIkjTgTAYkSRpwJgP12tjvABrE92KK78X+fD+m+F5M8b2Y0vP3wj4DkiQNOGsGJEkacCYDbYqId0XEzRHxrRmOR0S8PSJ2RsQ3IuLhvY6xV+bxXpwYEbdHxNfLn9f1OsZeiYhjI+JzEXF1RFwVEb9fUWYgPhvzfC8G6bNxj4j4WkT8e/l+vKGizN0jYnP52fhqRKzqfaTdN8/34vSIuGXaZ+Ol/Yi1VyLioIi4MiI+UXGsZ5+Lg7t14SXsPcDfABfMcPxpwJry55HA35Xbpeg9zP5eAHwpM5/Rm3D66k7gVZl5RUQMA9sj4uLM/Pa0MoPy2ZjPewGD89n4MfCkzNwbEYcAWyPiU5l52bQy64HvZ+bqiHgB8OfAWD+C7bL5vBcAmzPzFX2Irx9+H7gauGfFsZ59LqwZaFNmfhH43ixFTgUuyMJlwGERsaI30fXWPN6LgZGZezLzivL3cYr/uVe2FBuIz8Y834uBUf733lu+PKT8ae2sdSpwfvn7R4AnR0T0KMSemed7MTAi4hjg6cCmGYr07HNhMlC/lcAN017vYoD/IQQeXVYJfioiHtzvYHqhrMp7GPDVlkMD99mY5b2AAfpslFXBXwduBi7OzBk/G5l5J3A7cO/eRtkb83gvAH69bEr7SEQc2+MQe+ltwB8D+2Y43rPPhclA/aqytkHNfK+gmP7yV4C/Bj7e53i6LiKWAx8FXpmZP2w9XHHKkv1szPFeDNRnIzPvyszjgWOAEyLiIS1FBuazMY/34l+AVZn5UOASpr4ZLykR8Qzg5szcPluxin1d+VyYDNRvFzA9kz0G2N2nWPoqM384WSWYmRcBh0TEEX0Oq2vKNtCPAu/LzI9VFBmYz8Zc78WgfTYmZeYPgM8Dp7Qc+tlnIyIOBg5liTfBzfReZOZtmfnj8uU7gUf0OLReeSzwzIi4Dvgg8KSI+KeWMj37XJgM1G8L8JKy5/ijgNszc0+/g+qHiDh6sn0rIk6g+Lzd1t+ouqP8O/8RuDoz/3KGYgPx2ZjPezFgn40jI+Kw8vefA04C/qOl2BbgtPL35wKfzSU4Ccx83ouWfjTPpOhzsuRk5qsz85jMXAW8gOK/+W+0FOvZ58LRBG2KiA8AJwJHRMQu4E8oOsGQmX8PXASsA3YCdwC/2Z9Iu28e78Vzgd+JiDuB/wFesBT/gSs9FvhfwDfL9lCA/wuMwMB9NubzXgzSZ2MFcH5EHESR9HwoMz8REW8EtmXmFork6b0RsZPim98L+hduV83nvfi9iHgmxaiU7wGn9y3aPujX58IZCCVJGnA2E0iSNOBMBiRJGnAmA5IkDTiTAUmSBpzJgCRJA85kQJK6JCI+FhHXR8TdF3CNQyPiexHxd3XGJk1nMiDNQ0Rkmz+n9zvmxSwibo0ZlsZeLCLiCcCzgf+XmT+OiG7do5QAACAASURBVId08DkazczbgbcC/7ti6l6pFs4zIM1DRLy+YvcrKaYHPQ/4Qcuxj2fm1w88RfMREbcC383MRfvwi4h/Ax4IHJOZP42I+wAvbyl2N+DVwE+BP624zMbM3B3FUtDfBT6dmc/pZtwaTCYDUofKOcXvB9w/M6/rbzRLy2JPBiLiYRSLMf1VZp4xS7nlwDgwkZnL57jme4EXUizis6vOeCWbCaQuK+dj/4uIuCYifhQR34+IT0fEiRVlX1FWDz83In4tIr4SERMRcVNE/H358CAiHlle4wcR8cOI+GhEHLAcckRsi4i9EfHzEXFuRHwnIn4cETsi4qxy8ZOqmH85It4XETdGxE8iYk9EnB8RD6go+5Ey5qMi4lURcVX5d36iPP5zEfH7ZbzXl/e/LSL+NSKe3HKtZ0REUizT+uCWKvO/Kcs8ZPrrmf7mqutGxB9GxNoylu+X+46YVm5VRPxDRFxXxnlr2e5/fNW9ZrG+3G5u87zZfBA4iKm56qXauDaB1EUR8QvAZynWJf8c8EngnhQLsFwaEf8rM99fceoLgWdQLFTyZeAJwG8Dx0TEWynWObgU2AQ8HHgOxepmJ1SFAVwI/CLwMYolUJ8N/BlwPC3znUfEs4EPlOdtAf6bogbkhcAzIuJxmfntivtsAh5TxvYvFOsvUP7tby3/jk8Dt5b7nglcHBEvzswPlGX/E3gDxRrve4G/nXb9r1Xcs11Povi7P1vGezTFHPhExGTsw8CngA8DR1G8V0+LiKdl5ufneZ+TKNZcuKKGmCf9G8V/u5OBN9d4XQky0x9//OngB7iO4h/nVbOUuZziYfPMlv33plit7YfAYdP2v6K85o+BE6btP4jiYZoUC5Y8q+V6m8tjT27Zv63c/w1geNr+5cCV5bFnT9t/NMVDeA+wuuVajwB+BHypZf9HyutcS9E+3voe/DywomL/vSkWbdoNHNxy7FbgWzO8pw8p7/c3MxzfBuxt2feM8pwEXlxxzj0olovdO/19L4+tAm6hSIoOrrpnS/kjyvt8ZR5ll5dl985Vtiy/o/xvcEi/P//+LK0fmwmkLomIxwKjwHuzWIHsZzLzNuBNFN9Cn1lx+rsy82vTyt8FvK98eVlmfryl/AXldqbq7D/JzPFp19sLvLZ8+VvTyq0HhoCzMnNnS8zbgfcCayNipOIeb86KtuzMvCMrlmou34MLKFaye+gMcddta2a+r2L/r1PUVrxl+vsOkEV/kLdRJAWPmcc9Jt+bbixP/V3g7hQ1FlJtbCaQuufR5fbIGUYjTLbxP6ji2LaKfbvL7faKYzeW22NmiOULFfs+X24fNm3fZMy/GhH3rzhnVbl9EHB9y7EZq/HLDnWvoljeeAXFA226ldRbpT6TmWKc/LvXzPDf6sHl9kHAF+e4x73L7ffbC21evlduj6CoyZBqYTIgdc/kQ+Hp5c9MqnqR316x7855HDuk6lhmfq91Z2bujYgJiuGRkyZj/j8zxDqpKubvVhWMiCcC/0pRHX4J8HGKHvT7KPo4PI0Dk4NuqYyRqb/7xXOcP2uP/9L/lNt7zCui9vxcyz2kWpgMSN0z+dBen5nv6mMcB0fEvVoTgnJkwhBTtQowFfMDM/PaNu8z0zjlP6FIUh6ZmZe3xPBmimSgHfvK7Uz/fh02y7kzxTj5dz85Mz/bZjytbi639561VGcmr3nzrKWkNtlnQOqey8rt4/oaReEJFftOLLdXTtvXjZhXAze0JgKlx89wzl0UnSarTFa/H9t6oBwmuKrdAKn3776WoiPiL9VwrZ+JiGXAGuD6zOxGE4QGmMmA1D1foGgH/42IeGFVgYh4WEQc3oNY3lDOYjd53yGKDowA755WbiMwAfxpRPxK60Ui4uCq+RHmcB1w33KY5fRr/T6wdoZzbgNWRMQBzR5lZ8RdwJOnz3tQln07MycRs9lMUUPyqoh4UuvBKDyuKp6K+O6kGPlxv4g4uoNYZvJgig6nn6vxmhJgM4HUNZmZEfE8ivkA3h8Rr6IYavhDim+1D6P49vjLdKez2aQ7KKqVr4qI6fMM3A/4UGZ+bFrMuyPiBRQPxysi4mLgaoo5B46l6AB4MEUHtvn6K4rhh1+NiA9TfGt+JEV/gY9RzJHQ6lKKYZafiogvU0zXe3lmfro8/haKaaC/GhEfoahJeDLwE4ohmwfUGswmM/8nIp5DOX9DRHyRYjjmjylGB5xA8X4Nl7HM5aPAU8uf89uJZRZPmXZtqVbWDEhdVLa7Pwx4PcVD9CXA71I8DHdQDOXbOdP5dYVBMXxxM8WD9+UUHQ7/LxUd5jLzExRDFN9JUS39O8DpwHEUk/G0NQNeZn4UeC7F3/ki4Dcp2ujXMvO33NcA/1je82yKWoxfm3bNt1MkC7dRvIfPoUggHk9Rs9G2ckjhL1NMkHQf4KUUEz0dT9GM8EKmJlKay/sp/saXdBLLDE6jGMFxUY3XlADXJpCWtIjYBvxSzjHvvepXdo58NfCgzLxmgdd6NMUMhH+QmW+rIz5pOpMBaQkzGeifso/GDuBzmVnZZ6SNa11C0UxxXGbOp5lCaovNBJLUBeWMj78BXB0RHc+jEBGHAl8CTjMRULdYMyAtYdYMSJoPkwFJkgaczQSSJA04kwFJkgacyYAkSQPOZECSpAFnMiBJ0oAzGZAkacCZDEiSNOBMBiRJGnAmA5IkDbiD+x1APxxxxBG5atWqfochSVLPbN++/dbMPLLq2EAmA6tWrWLbtm39DkOSpJ6JiO/MdMxmAkmSBpzJgCRJA85kQJKkAWcyIEnSgDMZkCRpwJkMSJI04EwGJEkacCYDkiQNuMZPOhQR7wKeAdycmQ+pOB7AecA64A7g9My8oqdBjo/D5s2wYwesWQNjYzA83IxzJEmaQ2Rmv2OYVUQ8HtgLXDBDMrAO+F2KZOCRwHmZ+cjZrjk6Opq1zUC4dSusWwf79sHEBAwNwbJlcNFFsHZtf88BEwhJEgARsT0zRyuPNT0ZAIiIVcAnZkgG/gH4fGZ+oHx9DXBiZu6Z6Xq1JQPj47ByZbFtNTwMu3fD8uX9OQd6m0CYdEhSo82WDCyFPgMrgRumvd5V7uu+zZuLB22VffuK4/06Z3y8SATGx4tEAIrt5P69e6uvt3VrkXi88pVw7rnFduXKYv9MOjlHktQYSyEZiIp9B1R3RMSGiNgWEdtuueWWeu68Y8fUg7bVxATs3Nm/c3qVQHSadEiSGmMpJAO7gGOnvT4G2N1aKDM3ZuZoZo4eeWTlCo7tW7OmqHqvMjQEq1f375xeJRCdnCNJapSlkAxsAV4ShUcBt8/WX6BWY2NFG3yVZcuK4/06p1cJRCfnTBofh02b4Mwzi21VnwhJUtc1PhmIiA8AXwF+MSJ2RcT6iHhZRLysLHIRcC2wE3gn8PKeBTc8XHTGGx6eevAODU3tr+rU16tzepVAdHIO2M9AkhpkUYwmqFutQwuhaBffvLn4Frx6dfGgrXpA9/qcdkcTNH10xOS5jlqQpLYt+qGFdas9GWiybicQnZyzaVNRE1DVvDA0BOedB+vX1xNbp0w6JC0xsyUDjZ+BUAu0fHn1g3Uma9cW38zbSSDaPaeTfgbTRy1MLwvF/tlqE9pVlXSccUZ3kg5JagCTAR2o3QSi3XMm+xnMVDNQ1c9gPqMW2o25Si+TDklqiMZ3INQS1EnnxoWMWmiHQyUlDSBrBtR7k6MgZmr/r/rm3UltwqR22v97lXRIUoOYDKg/2u1nMDZWtNtXmak2Adpv/19I0iFJi5SjCbR4NHWopCQtAo4m0NLQbm1CJ50OO2nCkKRFzmRAi0s7oxY6bf/vZHilJC1iJgNauhbS/t/J8EpJWqQcWqilq5MhjJI0gEwGtHR1ssCTJA0gmwm0tNn+L0lzMhnQ0mf7vyTNymRA6idXR5TUACYDUr+4OqKkhrADodQP01dHnBz6ODExtX/v3v7GJ2mgmAxI/eDqiJIaxGRA6gdXR5TUICYDUj9Mzo5YxdURJfWYyYDUD86OKKlBTAakfnB2REkN4tBCqV+cHVFSQ5gMSP3k7IiSGsBmAkmSBpzJgCRJA85kQJKkAWcyIEnSgDMZkCRpwJkMSJI04EwGJEkacCYDkiQNOJMBSZIGnMmAJEkDzmRAkqQB1/hkICJOiYhrImJnRJxVcXwkIj4XEVdGxDciYl0/4pQkabFqdDIQEQcB7wCeBhwHvDAijmsp9hrgQ5n5MOAFwN/2NkpJkha3RicDwAnAzsy8NjN/AnwQOLWlTAL3LH8/FNjdw/gkSVr0mr6E8UrghmmvdwGPbCnzeuAzEfG7wBBwUm9CkyRpaehKMhARoxTf6g8HDqookpn5pvlcqurcltcvBN6TmW+NiEcD742Ih2TmvpaYNgAbAEZGRuZxa0mSBkOtyUBE3BP4GPBEqh/kkxKYTzKwCzh22utjOLAZYD1wCkBmfiUi7gEcAdy83w0zNwIbAUZHR1sTCkmSBlbdNQNvAZ4EfAl4N0UV/50LuN7lwJqIuD9wI0UHwRe1lLkeeDLwnoh4EHAP4JYF3FOSpIFSdzJwKnAF8MTWavpOZOadEfEK4NMUzQ3vysyrIuKNwLbM3AK8CnhnRPwBRY3D6ZnpN39Jkuap7mTgUOC9dSQCkzLzIuCiln2vm/b7t4HH1nU/SZIGTd1DC3cAR9V8TUmS1EV1JwPvAH4tIlbWfF1JktQldTcTfIqiA+GXI+INwHbgB1UFM/P6mu8tSZI6UHcycB1FJ74ANs1SLrtwb0mS1IG6H8gXcOCkQJIkqcFqTQYy8/Q6rydJkrqv6QsVSZKkLutau31EHAM8DDgMuB24IjN3det+kiSpM7UnAxExQrEGwMkVxy4GXpaZ19V9X0mS1Jm6Fyo6GvgyxdLD1wFfBPYAK4C1wFOArRExmpnfrfPekiSpM3XXDLyWIhE4E/jLzLxr8kBEHAT8AXAu8BrgFTXfW5IkdaDuDoRPBz6TmW+ZnggAZOZdmfkXwGeAZ9R8X0mS1KG6k4GjKWYdnM32spwkSWqAupOB24H7zVFmpCwnSZIaoO5kYCvw3Ih4TNXBiHgk8LyynCRJaoC6OxC+maLfwBci4oPA5yhGExwNnAi8ENgH/GnN95UkSR2qezriKyLiucB7gBcDL5p2OIDvAb+VmXP1K5AkST1S+6RDmfmJiLgfcCrwcOBQij4CVwIfz8yJuu8pSZI615XpiMsH/vvLH0mS1GAuVCRJ0oBbUM1ARLyk/PWfM3N82us5ZeYFC7m3JEmqx0KbCd4DJHAZMD7t9WyiLGMyIElSAyw0Gfgtigf7nvL1by7wepIkqccWlAxk5ntaXp+/oGgkSVLP1dqBMCIeHxEjc5Q5NiIeX+d9JUlS5+oeTfA54PQ5yrykLCdJkhqg7mQg5llmrk6GkiSpR/oxz8AIxcgDSZLUAAuegTAiXtey68SIygqCgygSgRfgqoWSJDVGHdMRv37a70mxOuGJs5S/ETirhvtKkqQa1JEMPLHcBvBZiomHqoYY3gXcBlyTmftquK8kSarBgpOBzPzC5O8RcT7FyoRfmOUUSZLUILWuWpiZzkAoSdIi46qFkiQNuNqTgYhYERHviIidEfE/EXFXxc+ddd9XkiR1pu7piFcC24DfBiaAuwPXAzsoOhAG8O/Al9q45ikRcU2ZXFSOQoiI50fEtyPiqoh4/0L/DkmSBkndNQOvA44GTsnMXyn3vTszfwl4APBp4OeA58znYhFxEPAO4GnAccALI+K4ljJrgFcDj83MBwOvrOMPkSRpUNSdDDwV+NfMvKT1QGbuAp5HkQy8YZ7XOwHYmZnXZuZPgA8Cp7aU+d/AOzLz++V9bu40eEmSBlHdycDRwFXTXt9F8fAHIDP3Ahdz4AN9JiuBG6a93lXum+4XgF+IiC9HxGURcUrbUUuSNMBqHVoI/BC427TX3+fAh/ftwJHzvF7VvMatixwdDKyhmPXwGOBLEfGQzPzBfheK2ABsABgZmXWVZUmSBkrdNQPfAY6d9vrfgSdFxM8DRMQy4CkU3/DnY1fL9Y4BdleUuTAzf5qZ/w1cQ5Ec7CczN2bmaGaOHnnkfHMRSZKWvrqTgUuBJ0bEIeXr84H7Av8WEW8Bvgw8GNg8z+tdDqyJiPtHxN0oFjna0lLm45RTIkfEERTNBtcu6K+QJGmA1N1M8I8UTQNHAHsy858i4hHA7wIPLct8EHjzfC6WmXdGxCsoRiEcBLwrM6+KiDcC2zJzS3nsKRHxbYo+Cn+UmbfV+ldJkrSERWZrE3wXbhJxJMXQwusy86au33AOo6OjuW3btn6HIUlSz0TE9swcrTpWd81Apcy8BbilF/eSJEntqXsGwndFxF9HxL1mKXNqRLyrzvtKkqTO1d2B8HTg5RQdBh8wQ5njgdNqvq8kSepQN1YtvJKif8BXIuIxXbi+JEmqUTeSgS3AOuAewCUR8fwu3EOSJNWkG8kA5doEj6XoNPj+iDizG/eRJEkL15VkACAzvwU8kmIWwj+NiI3lKoSSJKlBupYMAGTmd4HHAZ8EXgpcBBzazXtKkqT2dDUZAMjMOyhWKfwb4GTg97p9T0mSNH/dWKjoB607s/B7wBlUr0QoSZL6pNYZCDPz/nMcf1tEfIBipIEkSWqAnkxHPF0T1iaQNA/j47B5M+zYAWvWwNgYDA/3OypJXbCgZCAiRspfb8zMu6a9nlNmXr+Qe0vqoq1bYd062LcPJiZgaAjOOAMuugjWru13dJJqttCageuABB4E/Oe013PJGu4tqRvGx4tEYHx8at/ERLFdtw5274bly/sTm6SuWOgD+QKKB/vtLa8lLVabNxc1AlX27SuOr1/f25gkddWCkoHMPH2215Iaop32/x07pmoCWk1MwM6d3YtTUl9YVS8tde22/69ZU5SpSgiGhmD16u7HLKmnap1nICLuiojXzlHm7Ii4s877SprB9Pb/yYf7xMTU/r17DzxnbAyWzfBPw7JlxXFJS0rdkw4F85tUyImHpF6YT/t/q+HhotZgeLioCYBiO7nfzoPSktOPZoLDgR/14b7S4Om0/X/t2mLUwObNRZnVq4saARMBaUlacDIQEY9v2bWqYh/AQcAI8GLgmoXeV9I8LKT9f/lyRw1IAyIyFzYSMCL2Mf/hhAHsA16Sme9f0I0XYHR0NLdt29av20u9Mz4OK1fuP2fApOFh5wyQBkhEbM/M0apjdTQTvJEiGQjgdcDngS9UlLsLuA34XGb+Rw33lTSXyXb+1tEEy5bZ/i/pZxacDGTm6yd/j4jTgI9n5tsXel1JNbH9X9IcerpqoaQ+sf1f0iy6MpogIg4BnkyxZsHyzHxTuf8ewD2BWzNzhvFOkiSpl+qeZ4CIOIViwaJPAm8FXj/t8PHAHsBZSyRJaoi6ZyAcBT5O0aHwD4D9Rgxk5mXAfwPPrvO+kiSpc3XXDLwWuAMYLTsR7qgocznwKzXfV5IkdajuZOCxFKMJvjtLmRuAFTXfV5IkdajuZGA5cOscZX6+C/eVJEkdqvuhfCPw4DnKHA9cW/N9JUlSh+pOBj4FPDUiKhZJh4h4GvAY4BM131eSJHWo7mTgz4AfAJ+JiD8HjgOIiKeXrz9MMbTwL2u+ryRJ6lCtyUBm3gg8BdgN/BHwPIo1C7aUr/cAp2TmXP0KfiYiTomIayJiZ0ScNUu550ZElsMbJUnSPNU+A2FmXhERvwg8HXg0cG/gduAy4MLMvHO+14qIg4B3ACcDu4DLI2JLZn67pdww8HvAV+v5KyRJGhxdmY44M++iqA3YssBLnQDszMxrASLig8CpwLdbyr0JOBf4wwXeT5KkgdP0IX4rKeYlmLSr3PczEfEw4NjMtFOiJEkd6NZCRaMU3+oPBw6qKJKTixfNdamqc6fdZxnwV8Dp84hpA7ABYGRkZB63liRpMNSaDETEPYGPAU+k+kE+KSmq9ueyCzh22utjKDonThoGHgJ8PiIAjga2RMQzM3PbfjfM3AhsBBgdHU0kSRJQf83AW4AnAV8C3k1RxT/vDoMVLgfWRMT9KSY0egHwosmDmXk7cMTk64j4PPCHrYmAJEmaWd3JwKnAFcATM3PfQi+WmXdGxCuAT1M0N7wrM6+KiDcC2zJzoR0UJUkaeHUnA4cC760jEZiUmRcBF7Xse90MZU+s676SJA2KukcT7ACOqvmakiSpi+pOBt4B/FpErJyzpCRJaoS6mwk+RdGB8MsR8QZgO8VaBQfIzOtrvrckSepA3cnAdRTDBgPYNEu57MK9JUlSB+p+IF/AtEmBJElS89WaDGTm6XVeT5IkdV/T1yaQJEldZjIgSdKAq3ttgnfNo9g+4IfA1cAnMnNPnTFIkqT21N2B8HSmOhDOtOLg9P0/jYjXZOZbao5DkiTNU93NBA8ELgRuA14DnAg8qNy+ttz/z8Ajgd8GbgLOiYhTa45DkiTNU901A88CHgccn5k3Ttt/DfDFiLgAuBL4Uma+LSL+laK54BUUSYQkSeqxumsGNgAfbkkEfiYzbwA+XJabfP0J4OE1xyFJkuap7mRgFXD7HGV+ANx/2uvrgOU1xyFJkuap7mTgVuDkOco8haLvwKTDmDuBkCRJXVJ3MvBR4OER8U8RMTL9QESMRMT7gOOBj0w79AiKpY+lxW18HDZtgjPPLLbj4/2OSJLmpe4OhK+j6ED4ImAsIm6kGDFwFLASOAj4elmOiFgB/BR4b81xSL21dSusWwf79sHEBAwNwRlnwEUXwdq1/Y5OkmZV99oEP4yIxwB/DJwGPACYrCG4lmIho3Mz80dl+T3AY+qMQeq58fEiEZheEzAxUWzXrYPdu2G53WIkNVft0xFn5o8z802ZuRo4FDgWODQzV2fmGycTAWnJ2Ly5qBGosm9fcVySGqzuZoL9ZOY4YMOplrYdO6ZqAlpNTMDOnb2NR5La5EJF0kKtWVP0EagyNASrV/c2HklqU+3JQEQMRcQfRcQlEXF1RFxb8fNfdd9X6puxMVg2w/9Ky5YVxyWpwepetfAwYCtwHMXKhPekmEPgbsDPlcV2U4wgkJaG4eFi1EDraIJly4r9dh6U1HB19xl4DUUisB54D3AX8FfAmygWJ/obYAJ4as33lfpr7dpi1MDmzUUfgdWrixoBEwFJi0DdycAzgS9m5rsBIorVijMzgcsiYh3wTeBsilUMpaVj+XJYv77fUUhS2+ruM3AscMW01/uAu0++yMybgU8BL6j5vpIkqUN1JwN3UDQNTLodOLqlzE0UsxFKkqQGqDsZuIGidmDSt4HHR8RB0/atBb5b830lSVKH6k4GvgA8ISY7C8Bm4IHAJyPi/0TEh4FHARfVfF9JktShujsQnk8xjPAYilqCvweeBDyLYuligC9TjDqQJEkNUPdCRVcAvzPt9Z3AcyLiEcBq4Drg8sycYSJ3SZLUa11dm2BSZm4HtvfiXpIkqT09SQYk1Wh8vJjcaMeOYl2EsbFiFkRJ6tCCk4GIeEkn52XmBQu9tzRwtm49cNrjM84opj1eu7bf0UlapOqoGXgPkG2Uj7K8yYDUjvHxIhEYn7Yq+OTSyevWFdMhO/2xpA7U1UxwJ/AJinkFahURpwDnAQcBmzLznJbjZwAvLWO4BfitzPxO3XFIfbd5c1EjUGXfvuK40yH3n804WoTqSAa+ADyeYvjgfYB3Ah/KzB8t9MLlZEXvAE4GdgGXR8SWzJyedFwJjGbmHRHxO8C5gGvGaunZsWOqJqDVxESxQJL6y2YcLVILnnQoM58I/CLwFxTDB98N7ImIv46Ihy7w8icAOzPz2sz8CfBB4NSW+38uM+8oX15GMceBtPSsWVM8XKoMDRUrJap/pjfjTCZtExNT+/fu7W980ixqmYEwM3dm5pkUUxE/H/gqxXwDV0bE1yJifUTM8K/YrFZSTF40aRezr2uwnmIhJGnpGRuDZTP8L7tsWXFc/TOfZhypoWqdjjgz78zMj2bmKRTTEP8psALYCOyOiEe3ecmo2FfZWTEifgMYBd4yw/ENEbEtIrbdcsstbYYhNcDwcFHdPDw8VUMwNDS1386D/WUzjhaxrs0zUHbie21E/BvwDxTf6I9s8zK72H/ho2OA3a2FIuIk4GzgCZn54xni2UiRlDA6OtrO6AepOdauLUYNbN5cPFxWry5qBEwE+m+yGacqIbAZRw0XmfU/FyPivsBvlT/3A34EfAQ4OzN3tXGdg4H/BJ4M3AhcDrwoM6+aVuZh5bVPycwd87nu6Ohobtu2bb5hSNLcxsdh5cr9h35OGh526Kf6LiK2Z+Zo1bHamgkiYllEPDMitlCsQfBGYBz4feC+mXlaO4kA/Gxtg1cAnwauphilcFVEvDEinlkWewuwHPhwRHy9vL8k9ZbNOFrEFlwzEBH3p+i495sU/QMmKHr9vzMzv7bgCLvAmgFJXbN3r804aqTZagbq6DMw2StmG/AnwAcyc4ZeNJK0xC1f7uRPWnTqSAYC+ClFrcDrgNdFVA0C2E9m5v1quLckSVqgukYTHIKT/UiStCgtOBnIzFrnKpA0YJzLX+q7rs0zIElzci5/qRH8Vi+pP5zLX2oMkwFJ/eFc/lJjmAxI6g/n8pcaw2RAUn+4JLPUGCYDkvrDJZmlxnA0gaT+mJyzv3U0wbJlzuXfLQ7j1AxMBiT1j0sy947DODWLrixh/P/bu/foOcr6juPvTxMKSCK0CQpNgFChlZsSsaEICCW2CmIoSA9BELCx9mhpoaWnBHsBtUjFUi6l1WKgAuUSDRwaKcjhKqDcwqVWpJYcGiEECiTcfgGDgW//eJ7lt2xmf9nd7M7ub+fzOmfP5DfzzOyzz052vvPcZtD5QUVmVil+vLJR0iOMzcxsQG3IMM6XX4YFC+Dkk9OyKKCwcc/NBGZmw67TYZxuWqgM1wyYmQ27ToZxeobISnEwYGY27DoZxukZIivFzQRmVg1VHlbXyTBOzxBZKQ4GzGz4SjHJgAAAD65JREFUue27/WGctaaFooDAM0QOHQ8tNLPh5mF1nXG5DR0PLTSz6nLbd2dqTQuTJ492Ptxss9H1DgSGipsJzGy4ue27c54hMqlAfxMHA2Y23Nz2vWEmTYJ58/qdi/6pSH8T9xkws+E2rG3fFbhb7bshO3fcZ8DMqmsY277vvDNdpE48Ec48My2nTUvrrXsq1N/EzQRmNvyGqe27fmbAmloTyEEHjbu71YFWof4mDgbMrBqGpe27lbvVYficreikqaSdfSrU38TBgJnZeFKhu9UxddKxr919jjgibS/SbBrnccrBgJlZv/lutT2dNJV0sk8n0ziPU+5AaGbWT+12BuzkoUPDppOOfZ12Bqz1Nzn3XJg/Py1XrBiqYYXgmgEzs/7x3WpnOmkq2ZDmlbL6m/RxuKiDATOzZnr949xpZ8BBHx3R63LrpKlk0JtX+jy5kScdMrPxp4w7qKIf59rdd7d+nE8+OTUNNDN/PpxxRnfeqyydlls732knkwEN8gRCJeVtrEmHXDNgZt3Tr4t0t++gyhrLP+h3q+3qtNza/U47aSoZ5OaVARgu6mDAzLpjmC7SZf04D9vQtU7KrdPvtJOmkjKbV9oJjAdguOjABwOSPgKcC0wAFkTE3zVs3xi4BNgDWAkcERHLys6nWaUN20W6rB/nQb5b7UQn5bYh32knHfvK6AzYbmA8ADVEAz20UNIE4J+AA4GdgSMl7dyQbB7wfETsAJwNfKXcXJpZaXO4l3WRrv04F+n2j/MwDV3rpNwG4K64q+oD49rnWr16dP3IyLr7DMBw0YEOBoBZwNKIeCwiXgOuBA5pSHMIcHH+9yJgtiSVmEczG7aLdNk/zrW71TPOSMvxViNQ00m5lRl4laGTwHgAHqY16MHANOCJur+X53WFaSJiLfAiMKWU3JlZMmwX6QH4cR6XOim3Abgr7qpOA+M+1xANep+Bojv8xrGQraRB0meAzwBsu+22G54zMxtVVke4MtvYB30s/6Bqt9yGrd/EhrT/9/FhWgM9z4CkvYDTIuLD+e9TACLijLo0N+Q0d0maCDwNbBljfDDPM2DWA2WMy68ZGfFFetgMy3c6wPMZjOd5Bu4DdpS0PfAkMBf4REOaxcCxwF3A4cAtYwUCZtYjZd5JD8vjiG3UsHyn47SmY6CDgYhYK+l44AbS0MKLIuJhSV8ElkTEYuBC4FJJS4FVpIDBzPphWH7QzTbEOGxiGuhmgl5xM4GZmVXNWM0Egz6awMzMzHrMwYCZmVnFORgwMzOrOAcDZmZmFedgwMzMrOIcDJiZmVWcgwEzM7OKczBgZmZWcQ4GzMzMKq6SMxBKehb4aQ8OPRV4rgfHHY9cFqNcFm/l8hjlshjlshjVq7LYLiK2LNpQyWCgVyQtaTbVY9W4LEa5LN7K5THKZTHKZTGqH2XhZgIzM7OKczBgZmZWcQ4GuuuCfmdggLgsRrks3srlMcplMcplMar0snCfATMzs4pzzYCZmVnFORhok6SLJD0j6UdNtkvSeZKWSvqhpPeVnceytFAW+0t6UdJD+fU3ZeexLJK2kXSrpEckPSzphII0lTg3WiyLKp0bm0i6V9J/5vL4QkGajSUtzOfGPZJmlJ/T3muxLI6T9GzdufHpfuS1LJImSHpQ0rUF20o7Lyb26sBD7JvA+cAlTbYfCOyYX3sCX8vLYfRNxi4LgDsi4uBystNXa4GTIuIBSZOB+yXdGBE/rktTlXOjlbKA6pwba4ADImJE0kbAnZKuj4i769LMA56PiB0kzQW+AhzRj8z2WCtlAbAwIo7vQ/764QTgEeDtBdtKOy9cM9CmiLgdWDVGkkOASyK5G9hC0tbl5K5cLZRFZUTEUxHxQP73y6T/3NMaklXi3GixLCojf98j+c+N8quxs9YhwMX534uA2ZJUUhZL02JZVIak6cBHgQVNkpR2XjgY6L5pwBN1fy+nwj+EwF65SvB6Sbv0OzNlyFV5M4F7GjZV7twYoyygQudGrgp+CHgGuDEimp4bEbEWeBGYUm4uy9FCWQB8PDelLZK0TclZLNM5wF8AbzTZXtp54WCg+4qitqpGvg+Qpr98L/CPwDV9zk/PSZoEXAWcGBEvNW4u2GVoz431lEWlzo2IeD0idgemA7Mk7dqQpDLnRgtl8R1gRkS8B7iJ0TvjoSLpYOCZiLh/rGQF63pyXjgY6L7lQH0kOx1Y0ae89FVEvFSrEoyI64CNJE3tc7Z6JreBXgVcFhFXFySpzLmxvrKo2rlRExEvALcBH2nY9Oa5IWkisDlD3gTXrCwiYmVErMl/fgPYo+SslWVvYI6kZcCVwAGS/q0hTWnnhYOB7lsMHJN7jv8m8GJEPNXvTPWDpK1q7VuSZpHOt5X9zVVv5M95IfBIRPxDk2SVODdaKYuKnRtbStoi/3tT4EPAfzckWwwcm/99OHBLDOEkMK2URUM/mjmkPidDJyJOiYjpETEDmEv6zo9uSFbaeeHRBG2SdAWwPzBV0nLgVFInGCLi68B1wEHAUuAV4FP9yWnvtVAWhwOflbQWeBWYO4w/cNnewCeB/8rtoQCfB7aFyp0brZRFlc6NrYGLJU0gBT3fiohrJX0RWBIRi0nB06WSlpLu/Ob2L7s91UpZ/ImkOaRRKauA4/qW2z7o13nhGQjNzMwqzs0EZmZmFedgwMzMrOIcDJiZmVWcgwEzM7OKczBgZmZWcQ4GzMx6RNLVkh6XtPEGHGNzSaskfa2beTOr52DArAWSos3Xcf3O83gm6Tk1eTT2eCFpP+BQ4G8jYo2kXTs4j94fES8CZwF/UDB1r1lXeJ4BsxZIOq1g9Ymk6UHPBV5o2HZNRDy07i7WCknPAU9HxLi9+En6AfAuYHpE/FzSO4DPNST7ReAU4OfAlwsOc0FErFB6FPTTwA0RcVgv823V5GDArEN5TvHtgO0jYll/czNcxnswIGkm6WFMZ0fEn42RbhLwMrA6Iiat55iXAkeSHuKzvJv5NXMzgVmP5fnY/17STyT9TNLzkm6QtH9B2uNz9fDhkj4m6S5JqyX9n6Sv54sHkvbMx3hB0kuSrpK0zuOQJS2RNCLpbZLOlPRTSWskPSppfn74SVGed5N0maQnJb0m6SlJF0v61YK0i3Ke3ynpJEkP5895bd6+qaQTcn4fz++/UtJ3Jc1uONbBkoL0mNZdGqrMz89pdq3/u9lnLjqupD+XtE/Oy/N53dS6dDMk/YukZTmfz+V2/92L3msM8/JyYZv7jeVKYAKjc9WbdY2fTWDWQ5J+DbiF9FzyW4H/AN5OegDLzZI+GRGXF+x6JHAw6UEl3wf2A/4QmC7pLNJzDm4GFgDvAw4jPd1sVlE2gH8Hfh24mvQI1EOBM4DdaZjvXNKhwBV5v8XA/5JqQI4EDpa0b0T8uOB9FgAfyHn7Dun5C+TPflb+HDcAz+V1c4AbJR0VEVfktP8DfIH0jPcR4J/rjn9vwXu26wDS574l53cr0hz4SKrlfTJwPfBt4J2ksjpQ0oERcVuL7/Mh0jMXHuhCnmt+QPrufhs4vYvHNYOI8Msvvzp4ActIP84zxkhzH+liM6dh/RTS09peAraoW398PuYaYFbd+gmki2mQHljyuw3HW5i3zW5YvySv/yEwuW79JODBvO3QuvVbkS7CTwE7NBxrD+BnwB0N6xfl4zxGah9vLIO3AVsXrJ9CemjTCmBiw7bngB81KdNd8/ud32T7EmCkYd3BeZ8AjirYZxPS42JH6ss9b5sBPEsKiiYWvWdD+qn5fe5qIe2knHZkfWlz+kfzd7BRv89/v4br5WYCsx6RtDfwfuDSSE8ge1NErAS+RLoLnVOw+0URcW9d+teBy/Kfd0fENQ3pL8nLZtXZp0bEy3XHGwH+Ov/5+3Xp5gGbAfMjYmlDnu8HLgX2kbRtwXucHgVt2RHxShQ8qjmXwSWkJ9m9p0m+u+3OiLisYP3HSbUVX60vd4BI/UHOIQUFH2jhPWpl04vHUz8NbEyqsTDrGjcTmPXOXnm5ZZPRCLU2/p0Kti0pWLciL+8v2PZkXk5vkpfvFay7LS9n1q2r5fk3JG1fsM+MvNwJeLxhW9Nq/Nyh7iTS4423Jl3Q6k2ju1XqzTTLY+1z79jku9olL3cCbl/Pe0zJy+fby1pLVuXlVFJNhllXOBgw653aReGj+dVMUS/yFwvWrW1h20ZF2yJiVePKiBiRtJo0PLKmluc/apLXmqI8P12UUNJvAd8lVYffBFxD6kH/BqmPw4GsGxz0SmEeGf3cR61n/zF7/Gev5uUmLeWoPZs2vIdZVzgYMOud2kV7XkRc1Md8TJT0y40BQR6ZsBmjtQowmud3RcRjbb5Ps3HKp5KClD0j4r6GPJxOCgba8UZeNvv92mKMfZvlsfa5Z0fELW3mp9EzeTllzFSdqR3zmTFTmbXJfQbMeufuvNy3r7lI9itYt39ePli3rhd53gF4ojEQyD7YZJ/XSZ0mi9Sq37dp3JCHCc5oN4N093M/RuqI+O4uHOtNkn4B2BF4PCJ60QRhFeZgwKx3vkdqBz9a0pFFCSTNlPRLJeTlC3kWu9r7bkbqwAjwr3XpLgBWA1+W9N7Gg0iaWDQ/wnosA34lD7OsP9YJwD5N9lkJbC1pnWaP3BlxOTC7ft6DnPY8mgcRY1lIqiE5SdIBjRuV7FuUn4L8rSWN/NhO0lYd5KWZXUgdTm/t4jHNADcTmPVMRISk3yPNB3C5pJNIQw1fIt3VziTdPe5Gbzqb1bxCqlZ+WFL9PAPbAd+KiKvr8rxC0lzSxfEBSTcCj5DmHNiG1AFwIqkDW6vOJg0/vEfSt0l3zXuS+gtcTZojodHNpGGW10v6Pmm63vsi4oa8/aukaaDvkbSIVJMwG3iNNGRznVqDsUTEq5IOI8/fIOl20nDMNaTRAbNI5TU552V9rgI+nF8Xt5OXMfxO3bHNuso1A2Y9lNvdZwKnkS6ixwB/TLoYPkoayre02f7dygZp+OJC0oX3c6QOh5+noMNcRFxLGqL4DVK19GeB44CdSZPxtDUDXkRcBRxO+pyfAD5FaqPfh+Z3uX8FXJjf8y9JtRgfqzvmeaRgYSWpDA8jBRAfJNVstC0PKdyNNEHSO4BPkyZ62p3UjHAkoxMprc/lpM94TCd5aeJY0giO67p4TDPAzyYwG2qSlgDvjvXMe2/dlztHngLsFBE/2cBj7UWagfBPI+KcbuTPrJ6DAbMh5mCgf3IfjUeBWyOisM9IG8e6idRMsXNEtNJMYdYWNxOYmfVAnvHxaOARSR3PoyBpc+AO4FgHAtYrrhkwG2KuGTCzVjgYMDMzqzg3E5iZmVWcgwEzM7OKczBgZmZWcQ4GzMzMKs7BgJmZWcU5GDAzM6u4/wdXnU4X/7M/UwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x720 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "f = plt.figure(figsize=(8, 10)); # plot the calculated values    \n",
    "\n",
    "sp =  f.add_subplot(2, 1, 1 );\n",
    "plt.scatter(temp_arr, energy, s=50, marker='o', c='b')\n",
    "plt.xlabel(\"Temperature (T)\", fontsize=20);\n",
    "plt.ylabel(\"Energy \", fontsize=20);         plt.axis('tight');\n",
    "\n",
    "sp =  f.add_subplot(2, 1, 2 );\n",
    "plt.scatter(temp_arr, abs(magnet), s=50, marker='o', c='r')\n",
    "plt.xlabel(\"Temperature (T)\", fontsize=20); \n",
    "plt.ylabel(\"Magnetization \", fontsize=20);   plt.axis('tight');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the low temperature it is quite evident the system is settling to low energy state, with uniform magnetization density. While in the high temperature, the energy increases and magnetization decreases, signaling that the system is becomeing disorder. The analytical value of the critical temperature is $T\\sim 2.27 K$.\n",
    "\n",
    "# Now we will visualize the simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "from __future__ import division\n",
    "%matplotlib inline\n",
    "# Simulating the Ising model\n",
    "import numpy as np\n",
    "from numpy.random import rand\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "class Ising():\n",
    "    ''' Simulating the Ising model '''    \n",
    "    ## monte carlo moves\n",
    "    def mcmove(self, config, N, beta):\n",
    "        ''' This is to execute the monte carlo moves using \n",
    "        Metropolis algorithm such that detailed\n",
    "        balance condition is satisified'''          \n",
    "        a = np.random.randint(0, N)\n",
    "        b = np.random.randint(0, N)\n",
    "        s =  config[a, b]\n",
    "        nb = config[(a+1)%N,b] + config[a,(b+1)%N] + config[(a-1)%N,b] + config[a,(b-1)%N]\n",
    "        cost = 2*s*nb\n",
    "        if cost < 0:\n",
    "            s *= -1\n",
    "        elif rand() < np.exp(-cost*beta):\n",
    "            s *= -1\n",
    "        config[a, b] = s\n",
    "        return config\n",
    "    \n",
    "    def simulate(self):   \n",
    "        ''' This module simulates the Ising model'''\n",
    "        N, temp     = 64, .4        # choosing a ordering temperature\n",
    "        config = 2*np.random.randint(2, size=(N,N))-1\n",
    "        f = plt.figure(figsize=(15, 15), dpi=80);    \n",
    "        self.configPlot(f, config, 0, N, 1);\n",
    "        \n",
    "        msrmnt = 2**20 + 1\n",
    "        for i in range(msrmnt):\n",
    "            self.mcmove(config, N, 1.0/temp)\n",
    "            if i == 2*0:       self.configPlot(f, config, i, N, 1);\n",
    "            if i == 2*12:       self.configPlot(f, config, i, N, 2);\n",
    "            if i == 2**14:      self.configPlot(f, config, i, N, 3);\n",
    "            if i == 2**16:     self.configPlot(f, config, i, N, 4);\n",
    "            if i == 2**18:     self.configPlot(f, config, i, N, 5);\n",
    "            if i == 2**20:     self.configPlot(f, config, i, N, 6);\n",
    "                 \n",
    "                    \n",
    "    def configPlot(self, f, config, i, N, n_):\n",
    "        ''' This modules plts the configuration once passed to it along with time etc '''\n",
    "        X, Y = np.meshgrid(range(N), range(N))\n",
    "        sp =  f.add_subplot(3, 3, n_ )  \n",
    "        plt.setp(sp.get_yticklabels(), visible=False)\n",
    "        plt.setp(sp.get_xticklabels(), visible=False)      \n",
    "        plt.pcolormesh(X, Y, config, cmap=plt.cm.RdBu);\n",
    "        plt.title('Time=%d'%i); plt.axis('tight')    \n",
    "    plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/Users/animeshnanda/anaconda3/lib/python3.7/site-packages/ipykernel_launcher.py:48: MatplotlibDeprecationWarning: Adding an axes using the same arguments as a previous axes currently reuses the earlier instance.  In a future version, a new instance will always be created and returned.  Meanwhile, this warning can be suppressed, and the future behavior ensured, by passing a unique label to each axes instance.\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA7UAAAJvCAYAAACki/bIAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAAMTQAADE0B0s6tTgAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3de7Rt2V0X+O+surcqocpuRcKjLUnUDBUKMICIYqOI+OI1FB/xUSBEWlpaY1uCtt3agAn4LggojitaVx6iMBABBaEhgAICAnJ5CUEhBCo8Eh4BqhIrl2T2H3ufyqld+5y99rpzPeY+n88YZ1Tdc9aac+611v6d/Ttr/eYstdYAAABAj+5aegAAAAAwlqQWAACAbklqAQAA6JakFgAAgG5JagEAAOiWpBYAAIBuSWoBAADolqSWp5RSPrGU8o1LjwNgaeIhcFWId5wCSe0VUUqpB77eL8nfT/KhCw/1aUopv76U8nWllDeUUn6klPKipccE9K3HeFhKeUEp5QtLKT9eSnmilPKdpZQ/csn2X7J9LR8w5ziBdek03j2rlPI5pZQfKKW8uZTy0gu2e49SystLKa8vpfxcKeULz/3sN5dSvrGU8jPbn393KeVP7ez/QCnli0opP11K+flSyleXUt71gr7+0vZ47R0Ly7u29ACYzTuc+/+/muS9k3zYue/9bK31jfMO6XKllOtJvjzJrSTvlc2Yb5RSXlVrffmigwN61l08TPLuSR5L8sIkr07ywUn+VSnlA2qtX39+w1LKRyV59uwjBNaox3h3d5LHk/ydJC/et0Ep5Z2SfG2Sl223eXOSdzq3yRuS/ONsPkO+PskfSPLPSyk/Wmv9hu02n5OkJPndSZ5M8olJvryU8txaa93p6y8k+Z5Gr48JuFN7RdRaf/LsK8kTSd54/nu11jfuPn5SSvn6UsrfK6V8VinlF7d3Sj9w+5etr9neLfimUspzz/dVSnlxKeWHt38Z+7btXwHH+ANJfnWSF9Vav7fW+s+S/MtsAgvAKD3Gw1rrzVrrw7XWb6i1/nCt9dOTvDw7d1e2/X9Skj8zph/gtHQa756otX5srfVmkp+/YLOXJvk3tdZPqLV+X631+2utX3yuje+rtf6L7X9fWWv9zCTfneR9zrXx3kleVmv9rlrrDyT55Gw+d77dudd0LcnnJvmLSX52zOthHpJaDvmzSb43yXtkc9f0c5P80ySfluQ9s/kL1yNnG28fD/6LST42ybtk81ewryilPG/78/ctpTx+4Osdt839liTfVmv9xXPjeXk2QQhgbkvGw33eJuc+ZJVS7kry2Uk+odb6WKsXDVxJa4t3Tyml3J3k9yd55TYB/6nto8PvdsH2ZZtg/4Yk/+ncj745yQtLKfeXUu5J8hFJvivJT53b5v9N8n211n87ZGwsx+PHHPJNtdaXJUkp5SXZBKuvqbX+u+33Pj3JPzq3/V9P8uJa61du//0ZpZQPSfJQNn9V+/YkLzjQ549v//u2SV6z87PXJnnOyNcCcCeWjIdPU0r5w9k8ane+rvYvJXl8e3cD4E6sJt7t8Zwkb5Xk45N8XJJvS/Lnk7y8lPL8WutTd3dLKY9l83myJvmYc48eJ5tyji9O8gvbn/9Qkg84e/S4lPJbskl0D42bFZDUcsj5+oGzv1x938733nr7V7NnJ/k1Sb6glFLPbXNvNrVgqbW+Icl/H9h3GTVigGksGQ+fUkp5nyQ3k3x0rfWV2++9U5K/nOQ3H9sewB6riHcXOHvS9ItqrTeSpJTyMdnMNfCh2dxVPvO+SX5ZkvdP8qmllB+stZ7drf3kJG/abvNkkoeTfGkp5b2z+Qz6OUn+XK31dY3GzYQktRxy++x/aq21lPK072Xzl61k8+a/b/v/fzJPD3xJ8ovJ5vGTJP/+QJ/vXGv90WwC5m/c+dlzsrlbCzC3JeNhtvu8V5KvSPLxtdbPP7fdeyd5+yQ/uh3Xma8qpfyrWuvTZv0EOGDxeHeJn84mGX3FuTHeLqX8cDY1sTn3/Vdu//e7SykPZnNn98NKKc9P8jFJnnvWZynlTyd5XZLfl01S/xuS/NtzMfXuJL+jlPKRtdYHBoyTGUlqaek1SX4yyTvWWr/0gm2OefzkPyf5y6WU+2utj2+/9/5JvvWORwowrdbxMKWUd0/yVUleenZ34pwv2bZ33vdk86HtKwMwnebx7jLbya2+M8nzz763ndDpeUkuS4rvSvJL2/9/q+1/33S+6e3XXdnMMr+7vM/NJP8l52qJWQ9JLc1s/5L3KUleUkp5PMl/TPIrknxAkv9ca/3aIx8/+cpsgsqjpZRPyuZOxJ/IZlZkgNVqHQ9LKe+S5KuzmQH+80opb7/90RtqrT+/fTzudTv7JMmPmDQKmNIEn/9SSnnnJPckuT/J25VSXpDNnAFnbXxqkn9WSvm6bGpqz5b++bLt/h+V5OeyuXN8VzZ3Xz88yR/fbvcDSV6Z5J+UUv7vbB4//rgkb8ymnvh2NhNlnR/TE0leW2t9RVgdSS1N1Vo/o5TyZJK/kuRGkp/JZna5fzOirTeWUj5o2853ZPM48p+zRi3Qg5bxMJsJoX5lNpO1fOy57392ko+8s5EC3JnG8S7ZlFmcLRn0nkk+Osl/SPJ+2/4+v5TynCR/K5sE+tuzmeTpF7b7vDnJS5L82mwem35Fkg8/W/bn3GfMv5fNerd3Z7Om7QfWWn965JhZUDm3tjAAAAB0xTq1AAAAdEtSCwAAQLcktQAAAHRLUgsAAEC3JLUAAAB06+glfe4upT47d1/48yeu33v0IO67/eTBbca0O6bfFuNvNdYhx6WFKY7tPi2O05A2pjofh/oZYncsLdpo2c6xRp2PvOmNtdZ5LrqZnVJ83Ne3+DidKY71vnbEx/HtHGtMfExON0aWu+6uufasp/4913t4iCHnu+eYs9R74BTMFTOPHUcr3Z3n26+/MD4endQ+O3fnoTxw4c8fffD9j20yD9362oPbjGl3TL8txt9qrEOOSwtTHNt9WhynIW1MdT4O9TPE7lhatNGynWONOR838qrX3nHHK3VK8XFf3+LjdKY41vvaER/Ht3OsMfExOeEYee1Zuf7gC5/651zv4SGGnO+eY85S74FTMFfMPHYcrfR2nm/funlhfPT4MQAAAN2S1AIAANAtSS0AAADdKrXW43a45756vibiRQ2e8X70BU9/nntfm0O26cnu69ln9zUO2WdMG4eOZauxjjlnh877VK+nhSGvt8U5HdJOi2Pfyu1bN19da7248LRjc8THfe2Kj9O9l1q8j6f6nTbV7861xMirGB+T042R95dr9bI5B5a0phg6xbXV6vWsOZ600NvrWVvsmsNl8dGdWgAAALolqQUAAKBbkloAAAC6dXRN7bE1Ea2e3Z6i7nPKdo7tZ0wt1JqO7Zh+1lzDMqZ2t8W1s+R1PFfN9o286iTrxZJxNWMt3l9T1X0e28bQdo7tZ83xsVVfU81L0MJc52PN1/Gcc1qcaoxUU3u43yWd2mue6vWs6Zz1pNX5uCw+ulMLAABAtyS1AAAAdEtSCwAAQLcO1tSWUh5O8vBT37jr+q+6/m4PPfXPMTWCU60vOKTvFtZSgzrVPkMste7fUnVnp3gdH9LqHJ5yTe2dxsdkmtpQ8XG6fYaYq90W11crp3YdHzJVfExOJ0buxsd7Un7VR+UdZx+HGshhel/ffCqunzZWUVNba32k1vrA2Vfuvt5kUAC9Ex8B9tuNj9c9HAhMSIQBAACgW5JaAAAAunXt2B3uu/1kHrrkueie1voas9beVH3va6On9Sun6qfFMWixxmKrsU1VBzjVGp2HjDm2uXVzyiEt6tj4mPQVI8XHce3MtXbqVHWeQ/qZa13XMft0FR+Tk46RwLpNNY/EHNypBQAAoFuSWgAAALolqQUAAKBbkloAAAC6dfREUbumKA6ea1KeVpO2jJmEYq5JNMZMJDJmm6n6maKNMcd+qonIxvTd4vwMbafFWK6yqSZPmGtSnhbnu6f4uK+dY38+dCwt+hmipxgpPgJzWssER71bS6xzpxYAAIBuSWoBAADolqQWAACAbh1dU/vE9Xvz6IMXP4Peolal1eLph7ZpVfM4lSlqepasZxtThzvmWhljzNjWVD/V4lpXQ3bnjo2PybhjONf7vEXN41Smul4Pvcap6uvHxNAx18oYU8X3uYiP63AoPjKvuertl7ru1/T7ailXLea4UwsAAEC3JLUAAAB0S1ILAABAt46uqb3v9pN56Nwz2nPV6xzbxhBT1WMO6Wuq8U5RA9fqmfwWtYMttKqfmqtWocXaoFO1O2Tt56tkjvg4tt2lYuSpx8d97Y6xlvi4r2/xcVy74uPVsuR7FscWd2oBAADomKQWAACAbklqAQAA6NbRNbVzGFPvMqaeqlW7Uz3HP8X6lUP2aWFIHeCa6h9arFM7pp8hxqzn26K+cKpatRtHt3o65qrhHHO9io/D9mmlp/rLFuvUjulniFOLj8nVjpGHtPpdzDhrWffUOWcfd2oBAADolqQWAACAbklqAQAA6FaptR63wz331esPvvCpfy+5PuydGlsb0KI2rUV9zpB2d41Zo7XV+o+HXuNUNVhj6mOXqiEb+94Yc2ynqks/5EZe9epa6wN33NAKnVJ8TKYZ/5rj476+WtTXt1obe4oYOSY2TPV75NA4WrZzqM2l4mNyujFyNz6OsWRN51zrac9lLfWxY6zpOPam5/OeXB4f3akFAACgW5JaAAAAuiWpBQAAoFuSWgAAALp17dgd7rv9ZB66pMh4rsl+hrTRYiKOIfu0GH+r49Zq4qFD/ey2O9XC9GPGP+ZYHxpvq+tryFgOWdOkYkP2uUrWGh/3tSM+tpts5FC7Y2LOmH6GGLJPi/MsPg7bh370PjEU85nrve76ejp3agEAAOiWpBYAAIBuSWoBAADo1tE1tU9cvzePPnjxM9xDan6mWJS91QLyY/reNdVi9i20qsGaooZsqvqwIddKq2vjkCVrxI7V6rw/w62bdzSuNTs2PiaHr71W18yYflrVxh/bz1LxcWy7Y45ti5rmMWPbNdecCUOcenwc2s6pxsjdOQdOrR6wt9fT+/ttKT3Vxk81x8Baj4E7tQAAAHRLUgsAAEC3JLUAAAB0647XqT21Gp+p6prG9DPHOFppscblmuoNx9TqjtGq9m6KdZynOranbI74mMy3LuqYNsTHZ2q1BvBctcctauenut7GfBZYS3xkPXqpEWR5a7425hrbmo/Bee7UAgAA0C1JLQAAAN2S1AIAANCtgzW1pZSHkzz81Dfuup5H3+3OaknG1KZMUeMzVV1Qq/qpMZZa77HF8V+yhuyQqWriWvQzVotjO2S8u9vcGDrADvQSH4e006oW8ZA1x8eh2xyr1Tk8NA/BXK9nnxbrE4+JJ0P6GWOp+JicTozcjY/3pDzt53PV6c1V/zzV55E11XO3+J2w5nr0XmpHpzTXnAlTOHinttb6SK31gbOv3H19jnEBrJ74CLDfbny87uFAYEIiDAAAAN2S1AIAANAtSS0AAADdKrXWo3a4v1yrD+WBC3/eYqH3uSw5qcZS5pqkZcixbbHPkHbGTKYyZGwttJq8bIoJfaaa3OZGXvXqWuvFQaRjpxQfk3Hvpd6tJU612mdMGy3G1kKrWNdTfExON0Yeio9TWTLuztX3WmKxiaL609vEUJfFR3dqAQAA6JakFgAAgG5JagEAAOjWtak7WHJx7RaLp++zljrPVsbUN8xVc7Xm4zam3ynqwfbtN2SfFmNZc21ML5Z6n495L4mP4/ZZKj6OHcsUWh2DFtek+Hi1jPm92irmrKkWcQ5r+kw2xJrHtiZr+T1yiDu1AAAAdEtSCwAAQLcktQAAAHTr6JraJ67fm0cfHF4n0qoGq0V9S6s6rTWvOzbXeqtrWa9yyTXRWtQUtBrbmPqgXUNeT5PatFs3D+7Tq2PjY9JmDc0WdYVLrle6lvg4tt1D/SxZW7mWNbjHuJLxMTnpGNmzKdY7XpOp4tZc8a/34z+Guv2nc6cWAACAbklqAQAA6JakFgAAgG7d8Tq1Lep1htRttVgfrNVYxvSz5jWe5jpOU6z71+o4zvV6xmwzxlw1fUvWNPegxfEZUvsvPk5nLfFxyFimWiNyrtdz7M/HEh9pYU1xagpj1veFO3GnnwXcqQUAAKBbkloAAAC6JakFAACgW5JaAAAAunXHE0WNsZZF2aeaxGHMhEdD2zm2jSFaTFQzRIvzPmSCj0OT6LQ6p7tanI8W18W+duaa2GXI+bhxsJerbS3xMZkmBp1afNy3zVriY3J4bK0mg1nL+eg9PiZXJ0ZO9bt4jJ4mfZprsqa5+jGhGhfZe95v3bxwe3dqAQAA6JakFgAAgG5JagEAAOhWqbUetcP95Vp9KA889e8W9S1DLFXjM6beaE31AWNqUFvUh42pw52qPm+uGqtWx2lXi3ZbHadj+93X7o286tW11geeseEJODY+JsvV0LaIU6cWH/f1PSaGDjFVXeeuFvFdfJwvPianGyPLPffV6w++8I7a6Kn2dSq915euKRb0bi3XQqu5GYa4fevmhfHRnVoAAAC6JakFAACgW5JaAAAAujX5OrVT1XAOqXOaqt2l1lvd1VttV4vz0WIcQ/qZop5qzDj2jWXIfms6ts9o55I1xq6iKWo4x6wF2+KaGbrNsW2Ij/vNVZ881+/FVmM5tM+q42MiRl5iqrXcx7TBMD0dy6nqQOdaX3mpGtu11Pa6UwsAAEC3JLUAAAB0S1ILAABAt45ep/bYdcaWXP9tyFhatDumnzHr5LXoZ4i56s7GtDuVueqnxmhRQ9aqtr3FNXmqazAm49ZhnKuee4oaMvFxmKlqUNcSH5PlYuSpxcfkdGNki3Vq57LmOtC56j6nmrNkiLlic4t+W+l9/HOxTi0AAAAnSVILAABAtyS1AAAAdKv5OrVLPRM+dn3PMfuMqTuYa93TMTU9LdavHDKOMe1MsRblkDWNh5ii9njsesstzvuxPx+6z5prk+a2ZK15i3qpU4+P+/o+FB/37TNmLFPEx7H9DHnNh0x1TMb8vhIfmUOrzyNLmWpumrm0iFOtcgmW4U4tAAAA3ZLUAgAA0C1JLQAAAN06mNSWUh4upTx29pU33Z5jXACrJz4C7Cc+AnM6OFFUrfWRJI+c/fv+cq0+dK6QeqkJS8b0s2vspA1TjHeqCViGvJ65zseYSZKm0GriqzGTjyz1mqea3KZFmz075fi4r2/xcbn4OGXfu1rEhjHHX3w8LbvxsdxzXz20z1KTjQ5p49DY1nR+l4rVvRkSd3ua7Ouqn3ePHwMAANAtSS0AAADdktQCAADQrYM1tbueuH5vHn1w+PPYQ+oSTq1maUgba6qfGrPI/JAFq4/td4oapqH7HLoGh5yfqWoUW5iqzmJUu7duHt1PL46Nj8nha2Ku2DDmfdGqn0PbrDk+7ttmivjYsp1j9xlyHYyJocf2O5VVxcfkZGPkfbefzEMrqju9zFTX3pJzJBy7T6uxjnlfL1VffxXMFVeXOGfu1AIAANAtSS0AAADdktQCAADQraNraneNWberRd3MmBqfqZ7jb1FvNKTdIcbUNc1VD7uWNS+nug6mWu9sqXrlq77eWQstrrVW565F3eepx8d921y1+Dik7zE1zq2O9RRzMYiPrNma1r8do/fxT8ExGWbfcbpxyfbu1AIAANAtSS0AAADdktQCAADQraNrag+tM9Zq3dAW+xzaZs61R1us+zfEUuv8HRrH2G12DTluU5zDVutkruX8JPOtX3mVnFJ8TOZbe/SqxcdkXA3qIUOO7VTnsEUs6On8DNlHfOyHGsdn6u39N8QUr2nJddRPzZ2eH3dqAQAA6JakFgAAgG5JagEAAOiWpBYAAIBuHT1R1K6lJnaYapKhMX0vVUzfalKQKY7lVIXzYya3GeLQPlP1O0UbY9sZ8/4Y8164bOHsU9NzfBzbTs/xcd9+U/2umWISoaXi45C+l4xtLdqYKz4mVytGHqv3CXhMSPVMvZ9T1sOdWgAAALolqQUAAKBbkloAAAC6dcc1tYdMVds1VV3NXLWTS9W6tnp9h9qZ67zv62dMbVeLfcZco62OE/2aq6Z9jDF1hD3HxyHttnpfT9HP2DriQ3232GdIXfEu8ZGpuE6WM1UN7VLn9CrWBK/1NbtTCwAAQLcktQAAAHRLUgsAAEC37rimdoramyX32bXkun+H9hlSDzZVXecYY+qnDu3Tqg73kLnqwfaZa33EqermnuHWzcPbnIg1xzrxsa/4mIwb/1xrii91bE8uPiZXKkYestbavUTN91W15mvyqnOnFgAAgG5JagEAAOiWpBYAAIBulVrrUTvcX67Vh/LA4O2H1CJOtcZiC2taW2+uNSKHmKKmoNX6xC3Wf5xKi/Vwp6q9a1FrN6SfG3nVq2utw4NIR46Nj8m4GvC1xEjxcb+51mGcqzZ0Li3qVnuPj8npxsgW8fHUrCWWL2mq+UemsqZ5FpayZD35ZfHRnVoAAAC6JakFAACgW5JaAAAAunXHNbVz1cO2qEWcylyveal6gbnW8JuyRumQKeoD5lzbdsxxmKuWc7ed27dunmS9WLLe+Livb/GxTb9D+m51flqsRzzGXHXQU61t21N8TE43Ro6pqd211PmfyprrJufS+3q/c82h0KrvuX4nTHWO1NQCAABwkiS1AAAAdEtSCwAAQLcktQAAAHTrjieKamEtEz+0GsuSExlMdVzutM0h7Y6ZGGVN57CFIa9n7DYtjJk8ZdepToKSnFZ83NeO+Djd5HQtJlZa8zlspcXvhDXHx+R0Y2RPE0W16nuMNU2CtJSlJgIcY8l4OMZcE6eZKAoAAACOIKkFAACgWweT2lLKw6WUx86+bufNc4wLYPXER4D9xEdgTndcU9uqhuSyNsY61PfYftZcVzFFvdSYGqWp6mPnWjR6iDE1wofaGDu2uermWrzmU60XS46Pj/u2GWLNi6f3Hh+Xqv2fKoaO0eIcjqkRHtLGWuq6p4qPyenGyCnmHEim+6zXwhTxcM31pWsy1e+4Mf2siZpaAAAAWCFJLQAAAN2S1AIAANCta3fawKFnplvU1QzZZl8/a647G9PPXOt2rWktxzH9HNpmTJ1WizUkp+qnVTtz1dHduONe+jHkeLSowxtzTcxVl9VbfBwTP45tc2w7LfoZcpzG1KBOUYvcqp9j29jXzlzxMblaMXIOc73/hmhR3z3VGtwt2jxFa5knorfa3SW4UwsAAEC3JLUAAAB0S1ILAABAt45ep7bcc1+9/uALn/r3XM/xz7X+agut1iedq55hqef0p1rH8FAbU9XEjelnrmtyjKmui1NdgzGZJz7ua3epdTfHEB+Hmau+fqr6vKnWdV1LjJzyujjVGDnVOrU9a3Ud9VRTu6a5alrF2THtHmuqmLOWmDqUdWoBAAA4SZJaAAAAuiWpBQAAoFuSWgAAALp17U4baFFoPqZIufdJeHaNmRyr1T5jthkzKcixbQxpd8hEImMmLGlRkD/XBCwX9XVsvy0muhrS942De5yOFhNViI/zxceL9jvm5/v6GTOWVrHt0FiGjG2K+Dikn1Zt9BQfk6sVI2ljzfF8TLyYa4K+pSYCHGOqiSdPiTu1AAAAdEtSCwAAQLcktQAAAHSr1FqP2uHYxbPnWvx4qgWUh5ir77lqQ6eqnzrUzz5zHbcx+0xR/9XqnLbYZyqXLZzdu2PjYzJfbDvUzlV4n08RH8e2M6afKfqdql55rvrYU4uPyenGyDHxcU16qk+c63P1ED0dt1bWUtO8ttjWwmXx0Z1aAAAAuiWpBQAAoFuSWgAAALo1eU3tPkutpzVVrc2hNva1s9Qz7EvWbB7b5hCtal0PtTvV+Vqy1mSp+uvbt26eZL1YIj4Oseb4mCxXsznEVHWqd9pmq3aH9DOXJeenONUYqaZ2PlPF6jF6Om6trKVOdcnfPVNRUwsAAMBJktQCAADQLUktAAAA3Tq6prbcc1+9/uALB28/Zq2sNdeQtWpnqXUll1z7sEU/La6nJY/bUms3Lml3bKdaL5YcHx+Tce/ztcTIU4uPU7Y7Rz+tft9O0bf4uN++43KqMbKnmtqp3uc91Ze2mudkLnN97lyyPvlYS64/PxU1tQAAAJwkSS0AAADdktQCAADQLUktAAAA3bp27A733X4yD50rEJ5iUqExE6O0mkylRfHzmIkqhhzHuSbaGNJmi9cz1aQgPRnyelq8P8b03WqCn91tbgwdYIfmiI/72p1rQr6rEB+nmEhkzGRMrWJ5zzFyyOs5tfiYnHaMvGpafL4ChnGnFgAAgG5JagEAAOiWpBYAAIBuHV1Te6w1L7g+tpZhSDuHtKgRm6o+bIwW9UWtar+m6GdMbdcQU9WD7RrTbqs6zWf0fevmwXavkjXXVraoIxQfp62/PNZSNaitarin+L2xqviYXOkYOcV7ck115UPGMibuTjWWnoyZo2SMNf1u4encqQUAAKBbkloAAAC6JakFAACgW6XWetwO99xXrz/4wqf+3aJ+sUXdbYu1Bce2M6TdFv1MVTt0qJ8l138cY65rZUy7U9XKLPWax7iRV7261vrAJI0v7Nj4uG+bXWtef3WMNcfHIX1PtebsWuJjMs3xFx+HO9UYuRsf5zJVXf8YvY9lV2+1o6dWRzxG73Xrl8VHd2oBAADolqQWAACAbklqAQAA6NbBdWpLKQ8nefjs3/ek5KMueXa6RT3mvn2mqHOaag3dMZZci3JIu1P006LNuWqu5mpjzn7GvFenWFO1Z3caH5M2sW2qesW1nG/xcXy7c9VBz9HGnP30/H5Zi934mLuuH9xnivrnqWLBEHOtLT3HOMa2u6a62ynWuV6zNR37IQaN95J1vA/eqa21PlJrfeDs67qbuwBJxEeAi+zGx9x9OKkFGMsnMAAAALolqQUAAKBbzdep3bXkOnNT1fe2GMuYfeaq5RpSc9CiLmGutdfGvJ5DbQ6x5HUwxlz1hrdv3TzJNRiT4+Njstzar+LjOGPG0lN8HLrNoX2GaLGu/anFx+R0Y+SQdWrn+kw5xlrqL5esEV6zNdUaL2Wp98dc68Qnl8dHd2oBAADolqQWAACAbklqAQAA6FRfa0MAACAASURBVJakFgAAgG4dPVHU/eVafSh3Nn9Bi8khxljzBBmt2h1jLcX1S45jiuL6MRNsjZ0AosWxazFB2FWfKOqU4uOQsYiP+001tqXGsuZJ/U4tPianGyOHTBQ1hd4mA1pTzOnJmieJm0uL351r/yxuoigAAABOkqQWAACAbklqAQAA6NbRNbW7NRFrWSh7rrqaVu2Mqc8ZMrYpjKnzHPNcf6sa50NtDNHi/IyxZO3GXNfkjbzqJOvFkuPjY7KeGCk+jjdFneeUNZuH+jlkqlg9xFIxcs5r8lRj5Fw1tT3VQO6jpnacuX6X9m7N14qaWgAAAK4sSS0AAADdktQCAADQrWvH7nDf7Sfz0LnnsddU49hinzHmrOc9tp8Wx6nVWoGHjGljyNimqPEbss2YmquxNSFjxnLImNc8aPy3bh49ll6Ij8906vFx3/emOrYtYuQU8XFfO1PNmTBXzfaYcTQ71iccI7lza66JZP1a1P63qF+e4jp2pxYAAIBuSWoBAADolqQWAACAbh1dU7urxdp0c61nN1Wd01xrY81VNzdXPdiYfqaq02pVUzVFv2PqzMYcpzHWtA7rGvUcH/dtM0ZP8XHIfnO9l8b0M9XatlPExyHbiI+MYZ1XLjLFefYeHmeK4+ZOLQAAAN2S1AIAANAtSS0AAADdktQCAADQrTueKGqKQt8hk11MNcnQkLHMNXnKFIvZTzUxypLnY1erCWMOtdGin1YTWiz1PjS5xuWWmuRpqkmGhuzTc3wc0u5UYxny8ykmVmo1mVGL3wniI2vW4nPCqZtq4rkhWkzM2KLfIa7adTEXd2oBAADolqQWAACAbklqAQAA6NYd19S2eJa8xXPwY+unxrTbQov6r6nqPIeYqp7t2H1a1TLMda1MVfMxRU3HmDpAnm6K+Lhvmxbtio/7TXGsh2zTKp6s+Xf0LvHxapnq/beUucbaam6UQ1zP01F/PQ13agEAAOiWpBYAAIBuSWoBAADoVqm1HrXD/eVafSgPPPXvqWocd7WonzrU5lhLrYU1V41Sq2M717VyyJJrOR4yVZ3FVOdwzD63b918da31gWdseAJ242MLU9VWzrWeYE/xceh+U/S9pvkCWqxdPmSfMaaou11TfExON0aWe+6r1x984dLD6J51UffrqQZ4zedj7dfXZfHRnVoAAAC6JakFAACgW5JaAAAAunV0Te1uTUSLOpq51qZrUbPUqp2l1ndcU51Wq/Oxq0U97FJrRraqc5xqTboW7d7Iq06yXiw5Pj4m09TuDbHUGqdj2mhlzeu6TjW2XVPEx33b7Jpr3eDe42NyujHy1Gpq1zwvyBT9rp2a2jamOo6tXrOaWgAAAE6SpBYAAIBuSWoBAADo1sGktpTycCnlsbOvvOn2HOMCWD3xEWA/8RGY07VDG9RaH0nyyNm/yz33PW1mqRaFv1MVTE81thYTW00x4cdU+7QyxQQyU018tZRWk5e1aHeySatu3Tx6LGslPp5WfLyT/e7UmmLbVJP4teh3ili2qviYnEyMPBQfGWfM+3yuidrmspbPbAzX6pzduORnHj8GAACgW5JaAAAAuiWpBQAAoFul1uNKHO4v1+pDabsm+FQLu7cYy1x6qw9ocT5a1CSNqemba+HyMTUsS9YozvU+vJFXXbhwdu/miI/JuGtiyTrVO3Vq8TEZF9sOtdFqn6Vi5JgYtOb4OPa6PdUYOSY+tog5a44frWLqWmrnp/odMVdMatX3IWupTd5nze+X5PL46E4tAAAA3ZLUAgAA0C1JLQAAAN06uE7trieu35tHH7z4WfAW9S2tasjmqq1p0cZcdU1TrUPWok5rrrqmuepwh5ir/mSq9+GQfZ7hRNZg3OfY+Ji0OTfi4zBz1boOaXNMjBnzmueqW+2p7mzV8TE56Rg5hzXXBK553fEWpvqMuZbXN1ZP41/TZ+RjuVMLAABAtyS1AAAAdEtSCwAAQLear1M71VqXU9RpLfkM+FR1q0uZqo7iUD9zrcs41/qxY2vV1nI9DXnv3r518yTXYEyOj4/JdNf0oX3GjG0ua7meW1kqPg5pt1W971zrx7YY25i+5zw/pxojp1qndor1Vadah3hNcWnXmtepXZNTO+9rtu9YW6cWAACAkySpBQAAoFuSWgAAALrVvKZ2KVPVSY7pe8n6qTGmqJ8as89UtYRjTLW2bYva8CGWqu8dss9l9RC9W2t8TOZ7n0/R7xinFh/3bbOr1TrIh8y1PvjQvo+15viYnG6MXHN8nMuaayvV1I4zVdxd8zwRS7pszgF3agEAAOiWpBYAAIBuSWoBAADolqQWAACAbl07docnrt+bRx+8uFi5xWQKU2kx6cmYfvb1NdUkSWNe05jJLcZYquh9rvM81yQLrc7xYpNW3brZpN01OjY+JvPFpUPmitWnFh8v6utO+xEfh7W7q/v4mJxsjNyNj0vFujn1NLnPXBO1nZoWE8/1bq5JCg9xpxYAAIBuSWoBAADolqQWAACAbpVa63E73HNfvf7gCy/8eYtF5ls91z9FndO+dueq/xpSo9SiDqinfob0PaafqepjpxrbFLULU73HLls4u3dzxMd927ToZ4xTi48X7XenY1mynzFtjDm2Y4y5VqY678ea8j12qjFyNz6uua6wp1rYJS01F8BV0NOxHfNeniI+ulMLAABAtyS1AAAAdEtSCwAAQLcmr6ntbf2mqWqhejJVrdAUx22u5/hb1eEuuVbyIXPVs51qvVgyrqa2pxgpPi4XH6dsd4p+xMfx9b6nGiMPxcdkPXW2vcepuax5ThzmM9c65ImaWgAAAE6UpBYAAIBuSWoBAADo1rVjd7jv9pN56JJnp1vU0fSu1dqNc623urvPVHVzY9bWbHHchjjU7pBjMlc92FRaXF+n9l4+1rHxMTn9GtRd4mO7dudag3vI+Tg0/laxeini4zymWj97rdY01qt27E9R7znXkPHfuORn7tQCAADQLUktAAAA3ZLUAgAA0C1JLQAAAN06eqKoXWMmh9jVYlH2FpNsDB3bXAvTt5jgY0wbU/Wz1GQwh9ocss/Qdu7UVJOpzDWpSatjeyrGTDC2q0V8HNvPUu+lpeLjkHZave+neM1jYuiYsQ0x1YQxU0zIJz72rcU1MNV5MHESzMedWgAAALolqQUAAKBbpdZ6+QalPJzk4XPfeockPzHloBq6P8njSw/iCD2Nt6exJn2N99TH+pxa671TDGZuncfH5PSvtaX0NNakr/H2NNbkCsdI8XFWPY016Wu8xjqdpvHxYFL7jB1KeazW+sCRA1hET2NN+hpvT2NN+hqvsfart+PR03iNdTo9jbensSb9jXdKvR2Lnsbb01iTvsZrrNNpPV6PHwMAANAtSS0AAADdGpPUPtJ8FNPpaaxJX+PtaaxJX+M11n71djx6Gq+xTqen8fY01qS/8U6pt2PR03h7GmvS13iNdTpNx3t0TS0AAACshcePAQAA6JakFgAAgG5JagEAAOiWpBYAAIBuSWoBAADolqQWAACAbklqAQAA6JakFgAAgG5JagEAAOiWpBYAAIBuSWoBAADolqQWAACAbklqAQAA6JakFgAAgG5JagEAAOiWpBYAAIBuSWoBAADolqQWAACAbklqAQAA6JakFgAAgG5JagEAAOiWpBYAAIBuSWoBAADolqQWAACAbklqAQAA6JakFgAAgG5JagEAAOiWpBYAAIBuSWoBAADolqQWAACAbklqAQAA6Jak9kSUUj6xlPKNS48DYGniIcCGeMhVIantQCmlHvh6vyR/P8mHLjzUpymlXCulfFIp5UdLKU+WUn6wlPJ7zv38n+95LZ927ufPu+D1/vJz23xaKeWHSilvKKX8VCnl80opb78zjvtLKZ9eSvnJ7XbfU0p5t3mOAtBSj/GwlPKCUsoXllJ+vJTyRCnlO0spf2TPdu9RSnl5KeX1pZSfK6V84bmfvV8p5d+VUl5TSvnFUso3lVJ+187+H1NK+cbt/o8dGNMf3B6vz7vg5/9TKeVV222ujX3twHQ6jYfPKqV8TinlB0opby6lvPSC7f6vbcx8fSnly3Y/253bbm8sK6U8WEr5qlLK60opP1NK+eJSyjue+/lH7jlet879/P0uOKav3+nn15VSvqSU8gvb2PwfxcxlOOh9eIdz//9Xk7x3kg87972frbW+cd4hDXIjyXsl+egkP5jkHZP87M42/zrJnz/37yf2tPPbkvzIuX///Ln//64kX5zkR5O8XTbB+/OSfECSlFJKkn+T5HqSP5rksSTPT/ILI14PsLwe4+G7ZxN7Xpjk1Uk+OMm/KqV8QK3165OklPJOSb42ycuSvDjJm5O807k2fluSb0/y0iQ/neSjknxFKeUFtdZXbLd5VpIvTfKt2772KqW8bZJPTXLZ3ZvPSPL92cRtYJ16jId3J3k8yd/JJtY9Qynlo5L89SQfkeSHk3xaki9I8jt3trssln1ZNjHztya5d7vd5yf5X89t8xNJ3uPcv2+f+///lKcf32TzmfWV5/p/zrbvL96O7fEkL0hS970uJlZr9dXRVzYfaL5+z/c/Mck3nvv31yf5e0k+K8kvZpMUfmCSB5J8TTbJ4zclee5OOy/OJoC8Psm3JXm/keN812yCw6+7ZJt/nuTzLvn587IJDM8/ot8PSfKGnX//XJL/eelz58uXr7ZfvcTDC8b+VUkeOffvf53k5pFtvCLJi/d8/yOTPHbJfl+a5C9cFIOT/KFsEuPfvY3B15Y+1758+br8q8d4uB3LS/d8/78k+eRz//6121j0gp3t9sayJM/Z3X7P58NL4+SeMf3qJG9K8nvOfe8fJPm6pc+9r82Xx49P259N8r3Z/BXqy5N8bpJ/ms1fvN4zSUnyyNnGpZQXJfmLST42ybsk+Zxs7gQ8b/vz9y2lPH7g6+yv+h+U5IeS/LFSyo+VUl5RSvmEUsrdO2P8wFLKa0sp31dK+eRSyrP3vI6vLaX8RCnlq0spv/WiF1s2jyX/yTz9L3YflM1f6j5x+/jx95RSPvbwoQNOzJLxcJ+3yfbJlW1c/P1JXllK+fqyKaX46nJJmUQp5a4kb51nPv1yqe3r+p+S/MMLfv522dwt/shsPsABp2dt8fAppZR7k/ymbJ5cSZLUWn84m+T7vXfGdFEs+5kk/y3Jh5dS7i2l3J/kTyT56p3tnlNKeayU8sqyKV/7VZcM7SOyedLm5ee+90FJbm0fj35N2ZSF/M79uzM1jx+ftm+qtb4sSUopL8kmGH1NrfXfbb/36Un+0bnt/3o2f/X/yu2/P6OU8iFJHsrmL4Dfns1jFZf58e1/n5fk1yT5vUn+SJL/JZvHkW8n+ZTtNv8+m8dJfiybIPl3s3nU7cO3P388myD6n7J5XOXPJPkPpZR3r7X+17MOt0nq301yX5JvySbInHlekt+RzQe/D9r285mllNfVWj//wGsBTseS8fBpSil/OJtHi8/qap+T5K2SfHySj8vmLsifT/LyUsrza60/v6eZv5TNI8r/9sAYzvf7vO3Yf3uttW6qM57hs5J8eq31+7cJLnB6VhMP9/iV2cz585qd7782ydtux/e8XBLLaq1vLqX83mzi4/+ZTZL+bdl8Jj3zA9mUcXzvtt1PSPL1pZR3rbX+jz3j+ogkn1trffO57z0vyZ9L8je3+//RJF9VSnmnWusrn9kEU5LUnrbvOff/P7X97/ftfO+tt3cJnp1NEvoFpZTztQD3ZlMLllrrG5L894F935XkniQfWWt9VZJs/0r34myT2lrrF5zb/ntLKa/J5kPcw7XW19ZafzrJp5/b5ltLKb8xm+B7vg73X2Tz17d3zCao/NO8pabkrmwS6RfVWp9I8h2llBdkkyBLauHqWDIePqWU8j5Jbib56HMfes6emvqiWuuN7XYfk03t7YdmcxflfBsfls0jhR9yQcJ7kc9O8jcv+rC1rWN7m5y7QwOcpFXEwwvs/WvbjkOx7K4kn5nkvyZ5UTafRz8pm899H5QktdZvyeZGyNk+357N/CwfnOSLdtp7nyS/ftvveXdl8weCs5s131lK+QPZJPsvGfA6aEhSe9qeKng/95es80XwZ8GpZHOXM9k8vns+sCWbmouUUt43m7url3nnWuuPZhMQnzxLaLdekU3NxkW+Y/vf52XzF7mLtvn157+x/VD380n+WynlB5I8Vkp5t1rrd2/H8dg2oT0/jvN3c4HTt2Q8zHaf90ryFUk+fudJkZ/O5lHfswmfUmu9XUr54WzquHKujQ/O5tG/F9btJFNH+B1Jfnsp5exxvbu2bf7xbF7z78zm8b43bo/P2YfL/1FK+dha6z85sj9gnRaPh5f46WyeQnnbne8/J2+5e3solr1vkt+V5JfXWp/c/uwjkrx6eyf2e7Kj1vq6UsoPZZPA7/rIJN9ca/3Bne//VM7F7a1XZCduMw9JLWdek+Qnk7xjrfVLL9jmmMdLviXJvaWUB2qtZ0tLPD+bR40v8pu2//2RS7Z5t2weGbnI2R2PXzo3jj9YSnn29i+JZ+MYEliBq6l1PEwp5d2zmRzqpWd3Y8/UWt9YSvnObGLT2fbXsvkD3/mk+PdlU7Lxolrrlw9+NW/xrjv/fmk2pR1/Lckbk/w/2cwgf+a9kjyaTY2dmAlXU/N4eJla65OllO/KJil9eZKUUn5NNvHwW7ebHYplb5VNYn7+UeGz/987n1Ap5Zdlk9D+yM73n5XkjyX5K3t2+5aci9tbz0/yJfv6YFqSWpI89Ze6T0nyklLK40n+Y5Jfkc3SOP+51vq1Rz5e8lXZLAfxWaWUj8tmWvS/ls0EJNkW7f+NbGb8fE02ta6fluQLa62v3W7zEUn+R5Jb2SzJ82ey+evbw9ufv002jyJ/RTZ3dp+bzeMe35m3/OXsX2TzSPI/LqX8rW0/fzabZYYAnqF1PCylvEs2JRL/Msn5tbTfcO7x4U9N8s9KKV+XTe3X2VIXX7Zt43dlszzZS7KZW+CsjcdrrY9vt3n7JG+fTSnG9W2pRZL811rrG2ut37szrtdlM7Px2fdfvf06+/nbbP/3+2qtvxTgypng82FKKe+czSPB9yd5u22serzWetbGP0zyslLKd2Qz4/KnJvmGWuut7ZgujWWllG9O8mSSf1JK+Tvbvj45mwlMv3+7zd/IZs6WH87mrvBLspl/5St2hvuHsnnU+gvyTC/LJh6/eLvfH0vyzkn+8NBjQTuSWp5Sa/2MUsqT2fw16kY2s8d9czYfpI5t65dKKR+U5B9n8wHtNds2/8F2kzdls3bj2ex1P5ZNgnu+BqFkUwPx3GyS2+9O8gG11u/a/vx2kt+c5P9I8suzWW/s/0vyCbXWN23H8brt3Y1/mE1y/GNJ/mqt9QuPfU3A1dEyHmYzIdSvzOaPcOdnX//sbB5rS63188tmzcO/lc0Hxm/PJt6dran9p7OpbfuUvGWyvWQTIz9x+///ezZ/xDvzndv/PuPuA8BQjeNhskkAn7v9//fM5kbDf0jyftv+Ht1OVPeZ2Xy++5ok/9sR431tKeUDk/ztbO7u/tJ2vB9c37Ju71tnsxTQ22bzyPM3JHn/nXK1ZBN7v2Tf/AW11m8spfypbCaK+tvZTDr1+3dK75hJqdX6wAAAAPTJOrUAAAB0S1ILAABAtyS1AAAAdEtSCwAAQLcktQAAAHTr6CV9yl1311x71oU/v+/2k3c0oKGeuH7vLP3M9XqWNNexhCTJ7de/sdZ6khfd3aXUZ+fupYcxiPc9tNXq88ITedNJxshDnx8BDrrkM+Tx69Ree1auP/jCC3/80K2vPbrJMR598P1n6Weu17OkuY4lJMntWzdfu/QYpvLs3J2H8sDSwxjE+x7aavV54UZedZox8sDnR4BDLvsM6fFjAAAAuiWpBQAAoFtHP3583+0nr8QjuQAAAKyfO7UAAAB0S1ILAABAtyS1AAAAdEtSCwAAQLeOX6eWPPqC9us7vmjBybfG9D3FMWD9WlynNxqMA2BJS/7OBuCZ3KkFAACgW5JaAAAAunUwqS2lPFxKeezs63bePMe4AFZPfATYbzc+5k23lx4ScMIO1tTWWh9J8sjZv+8v1+qkI7rAUjWcakf3W6qeyPmYj5qxw9YSHwHWZl98fKiT3ys+a8D8hnzuvGxeFo8fAwAA0C1JLQAAAN2S1AIAANCtVa5Tq5YBAIAl7Kvt89mUq6qXOVbcqQUAAKBbkloAAAC6JakFAACgW6uoqVWnsP8Y9PIMO+vnWkKcfSZ1c8BQu/FCrIB1cacWAACAbklqAQAA6JakFgAAgG5JagEAAOjWIhNFKa4f5tBxugqT/7hWxrkK1wbcKfHlahAPAU6fO7UAAAB0S1ILAABAtyS1AAAAdGuRmlqgLTVjAGIhwFXlTi0AAADdktQCAADQLUktAAAA3Zq8ptY6gADAFNTQspQh157PwKzdKcVQd2oBAADolqQWAACAbklqAQAA6NbkNbVjntVWgzDMkOO0pmflr9p5XdOxB+iNGErvDl3DV+1zEUzJnVoAAAC6JakFAACgW5JaAAAAujV5Te0Yrepo1Cosx7EHAADm4E4tAAAA3ZLUAgAA0C1JLQAAAN2S1AIAANCtVU4U1cruhFMmLxrnKhy3VpOTATCOOAzAWO7UAgAA0C1JLQAAAN2S1AIAANCtk66p5WrUwx5L3RbAtMRZgGVdtTjsTi0AAADdktQCAADQLUktAAAA3TpYU1tKeTjJw2f/vidl0gFNaciz5WpQ+3LV6gVYl1OKjzCUuMsQ4iMwp4N3amutj9RaHzj7uu7mLkAS8RHgIuIjMCcRBgAAgG5JagEAAOiWdWp3tKgVUpcLrM2h2CZuXU3qYwE4Be7UAgAA0C1JLQAAAN2S1AIAANAtSS0AAADdMlHUBPZNvGESlsNMWAIwLXEW1sPnRWjHnVoAAAC6JakFAACgW5JaAAAAuqWmFlZuzfU16vOgL7vxxHsY4DRc9XjuTi0AAADdktQCAADQLUktAAAA3VJTC8CgWpw113cDnILdWCzuwjDu1AIAANAtSS0AAADdktQCAADQraNrap+4fm8effDi5/uv+hpJcJWMqvW5dbP9QAAAuLLcqQUAAKBbkloAAAC6JakFAACgW5JaAAAAunX0RFGHjJk45ipMLmUx7cP2HZOrcG0c4tphLVyLAMAauVMLAABAtyS1AAAAdEtSCwAAQLea19QCXFVPXL83jz54Z3WmPdWRr2ms6nvHMZfBsvZet7duzj8QgM65UwsAAEC3JLUAAAB0S1ILAABAt1ZRU2ttWy6yVJ2c6wuAltR9A0zHnVoAAAC6JakFAACgW5JaAAAAurWKmlpYm93apzXV2A4Zi9qtfrU4d2u6XgEYbzee+/3ORdb82XUO7tQCAADQLUktAAAA3ZLUAgAA0K2DNbWllIeTPPzUN+66PuV4rowWNRJzPSuvfmP/MbhqtQo8k/i4HmPej2Lbfle9LquVq3597cbHe1IWHA1w6g7eqa21PlJrfeDsK3f70AaQiI8AF9mNj9c9HAhMSIQBAACgW5JaAAAAuiWpBQAAoFsHJ4paqyETMPQ0ucWQsS416cSaxwYwlsmlAOA0uFMLAABAtyS1AAAAdEtSCwAAQLe6ram9itZcI7w7NnVny3I+rrZTm3NgTa7iHANjXo/r62peKwBLcacWAACAbklqAQAA6JakFgAAgG4dXVN73+0n89ARtTLqRWB5rWq7WtTJ3bjjFmhB3e101LQDwLzcqQUAAKBbkloAAAC6JakFAACgW5OvU7uvJkt90elbcn2+Q327/vZTPwnT8HuQi+y7Nsw7wGXEE4a6anNnuFMLAABAtyS1AAAAdEtSCwAAQLcktQAAAHRr8omi9plrYfrddk+pGPoUOB+cmvtuP5mHzl3XJu/gInP9HlzKVZugBOZyarECWnGnFgAAgG5JagEAAOiWpBYAAIBuLVJTC1Oaq95EzTaH9F43OWa83gcAw/T2OwHWzJ1aAAAAuiWpBQAAoFuSWgAAALp1pWpqW9UuqBljH+syAkB/1LZyVZ3S/DDu1AIAANAtSS0AAADdktQCAADQrVXU1Pa+liPz6f3aOKXaBY6373z3fk0DDPHE9Xvz6IPiHazZ6j+T3Lp54Y/cqQUAAKBbkloAAAC6JakFAACgW6uoqYWrytq2nNqcAurGGWrMte56AmAfd2oBAADolqQWAACAbklqAQAA6NbBpLaU8nAp5bGzr9t58xzjAlg98RFgv934mDfdXnpIwAk7OFFUrfWRJI+c/fv+cq1OOqKFDZm0ZamJXEyQAety1eIjLM1EZP3YjY/lnvvER2AyHj8GAACgW5JaAAAAuiWpBQAAoFsHa2qXoEbmalqqVnnt1JBdLUPq+oEN8RGAxJ1aAAAAOiapBQAAoFuSWgAAALq1yppaTlPPtYFz1WkNOUZDtlFXdjr2ncue30sAAK25UwsAAEC3JLUAAAB0S1ILAABAt9TUwh6916Rau5G1UAMOAEzNnVoAAAC6JakFAACgW5JaAAAAuiWpBQAAoFsmiuKgIRO90Ma+CXMc/6vtKpz/Q6/xFCeSugrndQljjuspXl8AV407tQAAAHRLUgsAAEC3JLUAAAB0S01tR3ZrhdQBAcCdaVXf7HcywHLcqQUAAKBbkloAAAC6JakFAACgW1e6ptY6gVwV+6519V/rJC4BABzHnVoAAAC6JakFAACgW5JaAAAAunWlamrVql1NPdWOukavHuf8MDXhAMBl3KkFAACgW5JaAAAAuiWpBQAAoFuSWgAAALp1pSaK4mowgQxLeeL6vXn0QRM/wVXUbNK3WzfbtANwhbhTCwAAQLcktQAAAHRLUgsAAEC3jq6pXapmrEWdZKtay2Z1M7Cg3eu4xftj0HtDvRgNTHH9AgB9cqcWAACAbklqAQAA6JakFgAAgG4drKktpTyc5OGnvnHX9SnHc6E11U/t9n3qNbZDO/dXqQAABLpJREFUjvWajsGarpVdazpOu9Y8trVaS3wEWBvxEZjTwTu1tdZHaq0PnH3lbkEJIBEfAS4iPgJz8vgxAAAA3ZLUAgAA0K2j16nlmeaqsV1TbSjAmqy5lh4AmJY7tQAAAHRLUgsAAEC3JLUAAAB0S03tivVUE7bmtWytvwoAAKfLnVoAAAC6JakFAACgW5JaAAAAuiWpBQAAoFsmimrARETDHJpMynEEAACO5U4tAAAA3ZLUAgAA0C1JLQAAAN1SU3vAXHWeh+pNr4K5joHaXWBp5hgAgHbcqQUAAKBbkloAAAC6JakFAACgW93W1A6pNxpTo6mOCQAAoB/u1AIAANAtSS0AAADdktQCAADQrW5raodQH8s++2qtT+1aUU/OVbfveu5pPfCrEKcAoBV3agEAAOiWpBYAAIBuSWoBAADolqQWAACAbp30RFE96X1Sk94dOtYmaIH+7b6PxVgAOA3u1AIAANAtSS0AAADdktQCAADQLTW1QJJn1heqI+bUucYB4DS4UwsAAEC3JLUAAAB0S1ILAABAt9TUrtiYNRXXUiNm/cf+jTmHa7n+AAC4OtypBQAAoFuSWgAAALolqQUAAKBbB2tqSykPJ3n4qW/cdX3K8TR1ajWBax4bJPvfc6d83R4bH6/a8QGurp4/PwL9OXinttb6SK31gbOv3C0oASTiI8BFxEdgTh4/BgAAoFuSWgAAALolqQUAAKBbByeK2nXf7Sfz0IgJmHqxO5GLSVzGWfK4jZkgbFfv533M+FsctyHt3pikl3Wa6pgCAPAW7tQCAADQLUktAAAA3ZLUAgAA0K2ja2ph7Xqvh13KkOOmRvRyY+YcUMcPAHBn3KkFAACgW5JaAAAAuiWpBQAAoFtXuqZW7Rocp8l75tbNO2/jhE1VtyzeAQCnyp1aAAAAuiWpBQAAoFuSWgAAALp1pWtqAQDmMLRe/sbE4wA4Re7UAgAA0C1JLQAAAN2S1AIAANCto2tqn7h+bx598OL1DqdaY7EF6zQCAHNY8+chgFPjTi0AAADdktQCAADQLUktAAAA3ZLUAgAA0K2jJ4o6ZMxkTCZTAJjWbpw1cR4AcCrcqQUAAKBbkloAAAC6JakFAACgW5JaAAAAuiWpBQAAoFuSWgAAALolqQUAAKBbzdepHcK6tAAAALTgTi0AAADdktQCAADQLUktAAAA3VqkpnYKj77g/Sdpd831v1O9ZgDgOGv+vABw6typBQAAoFuSWgAAALolqQUAAKBbpdZ6+QalPJzk4XPfeockPzHloBq6P8njSw/iCD2Nt6exJn2N99TH+pxa671TDGZuncfH5PSvtaX0NNakr/H2NNbkCsdI8XFWPY016Wu8xjqdpvHxYFL7jB1KeazW+sCRA1hET2NN+hpvT2NN+hqvsfart+PR03iNdTo9jbensSb9jXdKvR2Lnsbb01iTvsZrrNNpPV6PHwMAANAtSS0AAADdGpPUPtJ8FNPpaaxJX+PtaaxJX+M11n71djx6Gq+xTqen8fY01qS/8U6pt2PR03h7GmvS13iNdTpNx3t0TS0AAACshcePAQAA6JakFgAAgG5JagEAAOiWpBYAAIBuSWoBAADolqQWAACAbv3/yAnG2R0NaWwAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1200x1200 with 6 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "np.random.seed(2389)\n",
    "Ising().simulate()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Interestingly, this is a $2^{64\\times64}$ dimesional Hilbert space, but only with the $2^{20}$ steps we are pretty colse to one of the ground state, this is the power of MCMC simulations. <b>Run time: 5 mins</b>"
   ]
  }
 ],
 "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.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}