{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Notebook 6: Phases of the Ising Model with Logistic Regression\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Learning Goal\n", "The goal of this notebook is to show how one can employ Logistic Regression to classify the states of the 2D Ising model according to their phase. We will discuss overfitting, regularization, and learn how to use the scikit-learn library. We will also examine the role of the optimizer in making predictions.\n", "\n", "## Overview\n", "\n", "The energy function of the classical Ising model is given by\n", "\n", "$$H = -J\\sum_{\\langle ij\\rangle}S_{i}S_j,\\qquad \\qquad S_j\\in\\{\\pm 1\\}$$\n", "\n", "where the lattice site indices $i,j$ run over all nearest neighbors of a 2D square lattice, and $J$ is some arbitrary interaction energy scale. We adopt periodic boundary conditions. Onsager proved that this model undergoes a thermal phase transition in the thermodynamic limit from an ordered ferromagnet with all spins aligned to a disordered phase at the critical temperature $T_c/J=2/\\log(1+\\sqrt{2})\\approx 2.26$. For any finite system size, this critical point is expanded to a critical region around $T_c$.\n", "\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "An interesting question to ask is whether one can train a statistical model to distinguish between the two phases of the Ising model. If successful, this can be used to locate the position of the critical point in more complicated models where an exact analytical solution has so far remained elusive. \n", "\n", "In other words, given an Ising state, we would like to classify whether it belongs to the ordered or the disordered phase, without any additional information other than the spin configuration itself. This categorical machine learning problem is well suited for logistic regression. Notice that, for the purposes of applying logistic regression, the 2D spin state of the Ising model will be flattened out to a 1D array, so it will not be easy to learn information about the structure of the contiguous ordered 2D domains [see figure below]. Such information can be incorporated using other methods such as multi-layer deep convolutional neural networks (CNNs), see Secs. IX, X and XI of the review and the corresponding notebooks.\n", "\n", "## The 2D Ising Dataset\n", "\n", "To this end, we consider the 2D Ising model on a $40\\times 40$ square lattice, and use Monte-Carlo (MC) sampling to prepare $10^4$ states at every fixed temperature $T$ out of a pre-defined set $T\\in[0.25,0.5,\\cdots,4.0]$. Using Onsager's criterion, we can assign a label to each state according to its phase: $0$ if the state is disordered, and $1$ if it is ordered. Our goal is to predict the phase of a sample given the spin configuration.\n", "\n", "It is well-known that, near the critical temperature $T_c$, the ferromagnetic correlation length diverges which, among other things, leads to a critical slowing down of the MC algorithm. Therefore, we expect identifying the phases to be harder in the critical region. With this in mind, consider the following three types of states: ordered ($T/J<2.0$), critical ($2.0\\leq T/J\\leq 2.5)$ and disordered ($T/J>2.5$). We use both ordered and disordered states to train the logistic regressor and once the supervised training procedure is complete, we evaluate the performance of our classification model on unseen ordered, disordered and critical states. " ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "import warnings\n", "#Comment this to turn on warnings\n", "#warnings.filterwarnings('ignore')\n", "\n", "np.random.seed() # shuffle random seed generator\n", "\n", "# Ising model parameters\n", "L=40 # linear system size\n", "J=-1.0 # Ising interaction\n", "T=np.linspace(0.25,4.0,16) # set of temperatures\n", "T_c=2.26 # Onsager critical temperature in the TD limit" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Loading in the Ising dataset\n", "We now load in the data which is hosted on Pankaj Mehta's [website](http://physics.bu.edu/~pankajm/MLnotebooks.html)." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "import pickle, os\n", "from urllib.request import urlopen \n", "\n", "# url to data\n", "url_main = 'https://physics.bu.edu/~pankajm/ML-Review-Datasets/isingMC/';\n", "\n", "######### LOAD DATA\n", "# The data consists of 16*10000 samples taken in T=np.arange(0.25,4.0001,0.25):\n", "data_file_name = \"Ising2DFM_reSample_L40_T=All.pkl\" \n", "# The labels are obtained from the following file:\n", "label_file_name = \"Ising2DFM_reSample_L40_T=All_labels.pkl\"\n", "\n", "\n", "#DATA\n", "data = pickle.load(urlopen(url_main + data_file_name)) # pickle reads the file and returns the Python object (1D array, compressed bits)\n", "data = np.unpackbits(data).reshape(-1, 1600) # Decompress array and reshape for convenience\n", "data=data.astype('int')\n", "data[np.where(data==0)]=-1 # map 0 state to -1 (Ising variable can take values +/-1)\n", "\n", "#LABELS (convention is 1 for ordered states and 0 for disordered states)\n", "labels = pickle.load(urlopen(url_main + label_file_name)) # pickle reads the file and returns the Python object (here just a 1D array with the binary labels)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Constructing the training and the test sets" ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "X_train shape: (65000, 1600)\n", "Y_train shape: (65000,)\n", "\n", "65000 train samples\n", "30000 critical samples\n", "65000 test samples\n" ] } ], "source": [ "from sklearn.model_selection import train_test_split\n", "\n", "###### define ML parameters\n", "num_classes=2\n", "train_to_test_ratio=0.5 # training samples\n", "\n", "# divide data into ordered, critical and disordered\n", "X_ordered=data[:70000,:]\n", "Y_ordered=labels[:70000]\n", "\n", "X_critical=data[70000:100000,:]\n", "Y_critical=labels[70000:100000]\n", "\n", "X_disordered=data[100000:,:]\n", "Y_disordered=labels[100000:]\n", "\n", "del data,labels\n", "\n", "# define training and test data sets\n", "X=np.concatenate((X_ordered,X_disordered))\n", "Y=np.concatenate((Y_ordered,Y_disordered))\n", "\n", "# pick random data points from ordered and disordered states \n", "# to create the training and test sets\n", "X_train,X_test,Y_train,Y_test=train_test_split(X,Y,train_size=train_to_test_ratio,test_size=1.0-train_to_test_ratio)\n", "\n", "# full data set\n", "X=np.concatenate((X_critical,X))\n", "Y=np.concatenate((Y_critical,Y))\n", "\n", "print('X_train shape:', X_train.shape)\n", "print('Y_train shape:', Y_train.shape)\n", "print()\n", "print(X_train.shape[0], 'train samples')\n", "print(X_critical.shape[0], 'critical samples')\n", "print(X_test.shape[0], 'test samples')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Visualizing the states" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0IAAAERCAYAAABW2xcEAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XvUZWV92PHvj+EyAUzyTpjgKohAvaSgoxOpIZgGxSYQq4Mp3mJrDFliaqNNtCZRqpZ4WSa1mrRpdInaxhjiPcCoXVWRES8D1gEUGBUkDjcDdbiI4eIo8PSPvV89nDnnffd+9+Xsc57vZ62zZt5znr33s2+/c35n7+d3IqWEJEmSJOVkn1l3QJIkSZL6ZiIkSZIkKTsmQpIkSZKyYyIkSZIkKTsmQpIkSZKyYyIkSZIkKTsmQpIkSZKyYyIkSZIkKTsmQnMmIj4bEZ+ddT+Wtd2fiDgrIlJE7NvWPCXVt3wujvz9zIh4RZW2XfZj1obWH6lP4+/5Qzof/HyktTARkiRN8m7gF0f+fiYwMRGa0HaR5bSu0mo8HzTXzCoHJiIOSCntWbRlSZoPy3EhpXQTcFOVaeq07UtX8W2I6yrNSp/ng5+P1AWvCHUoIk6JiIsj4t6IuDMizouIR4+8vnyZ8zER8cmIuAv40Mjrz4uIb0TEnojYGRG/PmU5j4uIrRFxR7msL0bEvxhrM3VZVaav058pfVxe/mMjYltE3BMRN0fE6yNi0nF4VER8IiLuiojrI+J1o+0i4hER8b6I2FX2+VsR8Y6IWBpb7qMi4tyI+E5EfD8iboiID49eWq66/tI8Ko/vcyPitvL4vjoiXl2+tlJc+NEtLxHxV8ALgcPK9ikirhtZxl63x6y03PL1SudwxXVcLZZWjXG/Uca470fElRGxJSrcCrRarB/r4yNXim3SUFR5zx8/Hyq+59Y5X/x8tLejVoohVWNrxX218J+PDL4diYhTgE8AdwHPBV4CPAb4QkQcNtb8fOAiYAvwZ+X0/xL4W+CbwL8G3gL8N2A8WPw8sB3YAJwBnAbcBlwQEU+Y0LUHLavq9FX7U8F5wAUUt9n8LfBa4HUT2p0LXFi2Ow/4Y4oPYsv+CcW3UL8PnAy8Hngq8L/H5vNx4DCK7X8y8CpgD+Wxv4btJ82NiHgicDHwT4GXA/8KeBtw+FjTvWLQmDdQnFu7KW6D+UVg6ht9xeVWPYfrmBRLq8a4XwHOAb5RtvmvwJ8Dj1ppgTVjPawe26SZa/Cev9p7bt3zxc9He/PzUZtSSj46eAA7KE6KfUeeOwr4IfC28u+zgAT83oTpvwh8Ddhn5LlfKNt/duS5zwBfB/YfeW5d+dx5I89NXFaN6Sv1Z4Xtsbz8V409/y7gH4GfHmt3+li7K4FPrTD/fYFfKqfdXD53SPn3lhWmq7T+PnzM4wP4HHAjcOCU11eKQWcVbxE/+vuvgJtWmk/V5U6Zx17n8KR5r2E9qsa47cBVQIw89/MTYu74uq4a68f6WCu2+fAxi0fV9/zR86Hie27d88XPR34+6vThFaEORMRBFG+gH0wp3bf8fEppF8UJc+LYJOeOTb8O+OfAR1JKD4xM/yXgupF2P1HO68PAAxGxb3lJMyi+VfjlCd07t+70VftT0YfG/v4AcDDFN0KjPjH291XAESN93z8iziwvRd9LEUQ/X768/C3MbcC3gD+JiDMi4pGjM1zj9pPmQkQcCDwJOCeldM8qzc9d5fXWl1vxHK5rPJbWiXHHAR9d/kQHkFK6DNi1wjrUjfWwSmyTZq3Be/5q77lrOV/8fOTno06ZCHVjieJguXnCa7dQXGYcNd7uEGA/4P9NmH70uQ0U2flrKQ700cdLgaUJ95fevIbpq/anivH2y3+PXxK/fezvPcD6kb/fTPHtyN9Q3HbzRIpL0iy3Kz/Q/ArFN1BvBq4p75V9SdluLdtPmhdLFDG+ykDmSbGq6+Wueg6vwfh61I1x35kwz5ViXN1YD6vHNmnW1vSeX+E9dy3ni5+P/HzUKavGdeMOikuOD53w2kMpMvFR4zX4b6U42A6dMP2hwPXl/78LPAD8JfDXkzoy+g3FhGVVmj4iqvanikMpvoUY/Rvg2zXmAfA84K9TSm9cfiIiDh5vlFL6FvCbERHA4yhO4LdHMdD7c9TfftK8uIPi+J503/24Nn8HpOpyK53DNY2vR90Y97MTmhwK3DBleXVjvTQP1vyeX+E9t+754ucjPx91aiGyuaFJKd0NXAo8u7xsCkBEPBw4gWIw3krT3w98GXhWPLgSyC8AR44t5/MUB/BlKaUd448K/Vx1+qr9qeg5Y38/j2LQ5FU153MgRfAZdfq0xqnwFX78OyiPabr9pCErb0v7AvBvy9scmtoDrDqfGsutdQ6vRc0YtwM4rfxQAEA5IPioVea/5lgvDVEb7/krvOc2Ol/8fFSJn49q8IpQd15LcR/nxyPi7RT3ef4xcCfw1grT/2fgU8B5EfFOYGM5/S1j7V5Bkbl/MiLeQ3Fp9xCK+3DXpZRetcpyqk5ftT+rOaMMFl+mqFLyIuCslNJ3a87n/wAvjIgrgWspLvueMNogIjZRVG75YNlmHfBbwH0UFVeg+faThuyVFB8uLo6It1LcrnY08PiU0stqzutrwIby1okdwPdTSlc2WO6q53BL6sa4cyPi7LLNWRQxbqVvPpvGemmIar/nV3zPbeN88fPRyvx8VMdq1RR8rP0BnEJRQvZeipP8fODRI6+fRXEpdt8p0/8GcDXFN7E7KcrVfpaxKiTAP6MYVPedsu1NwFbgaVWWVWX6Ov2Zsi7Ly38MsK3cJrdQlOXdZ7V+UlSsum7k70PKPt9RPs6hGLCYgN8q2/ws8F7gGuAeivtqLwJOXsv6+/Axjw9gM/Axils97qUoD/1H5WsrxYWzeHB1tIOA9/Pj28Gum9Z2teWWr696Dk+b97S+TlqP8vWqMe75E2Lc5cC5q6zrirF+pT6OxzYfPobyoMJ7Pg+uGlf1PXfN50udvpXt/Hzk56MVH1GuqNSpiDiL4luT/dJItRhJGqqIOJzi29I3pZTeMOv+SFo8fj6aLW+NkyRlrxzP9DaK0rC3UtzK94cU35a+e4ZdkyR1xERIkiS4n6Jy1f8AfgZYHiz87JRSm+XFJUkD0eqtcRHxMODPKGqTL//o0u+nlKaVHpWkRow7kvpm3JEWQ2uJUPlr4l+lGEz1GopBWW+kKOO3KRWl+CSpNcYdSX0z7kiLo81b486guKf60SmlawEi4grgm8DvUNx7vaKfOWRdOuLh3q0ndeWG6+/jtlvvj9Vbzo1GcWdSzLni8o0T227avHuv5+q07cqkPvS5/L5N2+aTLPJ2mCdfuewHt6aUqu+44WsUdyIOSvvEUuedbEvV86jOuVlnWV3Nd5o2ltfX8me9XnV1tX8neSB9u1LcafOK0GeA9SmlJ409fxFASunE1eax+QkHpG3bq/wIuqS1eMoJ3+byS/csTCLUNO5MijkP+6kzJra98c537fVcnbZdmdSHPpfft2nbfJJF3g7zZGn9rktTSsfNuh9taRp31u1zeDpw/9/tsIftqnoe1Tk36yyrq/lO08by+lr+rNerrq727yR37TmzUtzZZ7UGNRzL5F+/3Qkc0+JyJGmZcUdS34w70oJoMxHaQPHDTeNuB6ZeA46IF0fEjojYcevu+1vsjqQM1I47xhxJDTWKOw4hkoajzUQIigGD41a8DSeldHZK6biU0nGHbFzXcnckZaBW3DHmSGrBmuNOxEEddktSHW1WJriD4luScUtM/uZEkppqPe7M27iSWY9HmqbPfs3bPtPcm2ncGcJYj0mGvA6zHkvZ1frWGafa5/icvmPypOUtra82bZtXhHZS3Dc77hjgay0uR5KWGXck9c24Iy2INhOhrcDxEXH08hMRcSTwpPI1SWqbcUdS34w70oJoMxF6F3AdcH5EnBoRW4DzgRuBd7a4HElaZtyR1DfjjrQgWkuEyl9SPgm4BngfcA6wCzgppXRXW8uRpGXGHUl9M+5Ii6PNYgmklG4ATmtznpK0EuOOpL4Zd6TF0GoiJEk5mVYZp2klnq7m27QPbVQCGmrFKyl3Tc/vNqpKdhU7u1hWHW2sbxdV3+r0q4116HP/VtX27whJkiRJ0uCZCEmSJEnKjomQJEmSpOyYCEmSJEnKjsUSJKllXQxq7Wq+XQw+1Xyadix5jLRr0+bdbNv+4G3axrbvs9hBU13Fw6EaQpGYrvrQ9H2puwIXZ1Zq5RUhSZIkSdkxEZIkSZKUHRMhSZIkSdkxEZIkSZKUHRMhSZIkSdmxapwk9aCNCkddVP0ZRnWf6uapUtS8cdv244rLN+51fkzb9pPOozptp5k0j6bz7Wod+o6dfVZBq6NOH5rus6baeK/oKx55RUiSJElSdkyEJEmSJGXHREiSJElSdkyEJEmSJGXHYgmSpL00Haja1YBuDcOsBzjPs02bd7Nt+4O3UxcD4aHfAiddrUMbx1qfx2XT7TiEc6irYjeTtLEfm/TXK0KSJEmSsmMiJEmSJCk7JkKSJEmSsmMiJEmSJCk7JkKSJEmSsmPVOEnZuuLyjXtVmxlCxZ46mva37+pNTaef1Id522eLwG0+3+qcR11UEGujml3TWNA0nnW1vdrY3lXXoc/qcHX70FeM8YqQJEmSpOyYCEmSJEnKjomQJEmSpOyYCEmSJEnKjsUSJGVr0+bdbNvez4DMNgaE9jlAvWkRhWl97bOwgjTvhnDOT+pDn4Ps+yziUEed5bdRxKGOqjG5jeOraX/73jbjvCIkSZIkKTsmQpIkSZKyYyIkSZIkKTsmQpIkSZKyYyIkSZIkKTtWjZOkHgyh+lNT81axrU41O2mWujouu6ps1tSQK65V7cOsq9bV1TQetrG+XVUNbdI3rwhJkiRJyo6JkCRJkqTsmAhJkiRJyo6JkCRJkqTsWCxBkuaYhREkTVNncPlQB/9P61dXBQyaFhUYamGFocb0NvbvpLZL66st3ytCkiRJkrJjIiRJkiQpO5USoYg4PCL+IiIujoh7IiJFxJET2q2PiLdExM0RcW/Z/pfb7rSkxWfckdQ3446Ul6pXhB4BPAe4A/j8Cu3eA5wBvA54OnAz8MmIeHyTTkrKknFHUt+MO1JGqhZL+FxK6VCAiHgR8KvjDSLiccDzgd9OKf2v8rmLgJ3A64EtrfRYUi6MO2OaDrht+ivibfwKeZ8DdpsOWh7q4GJ1Ktu4U6ewQp996OqcrTPfrtpO0tV+GELs63fbnFlpnpWuCKWUHqjQbAvwQ+CDI9PdB3wAODkiDqjUI0nCuCOpf8YdKS9tFks4FtiVUrpn7PmdwP4Ul5slqU3GHUl9M+5IC6LNRGgDxT21424feX0vEfHiiNgRETtu3X1/i92RlIHacceYI6kh4460INpMhAJIU56fKqV0dkrpuJTScYdsXNdidyRloHbcMeZIasi4Iy2INhOh25l81Wdp5HVJapNxR1LfjDvSgqhaNa6KncCvR8SBY/fNHgP8ALi2xWVJEixo3OmzOlxX2qh8NOsKc9IUjeLOFZdvrHwedHVc1plv1UpuXVVsG0I1u0XQtCJf3/OdpIv5tnlFaCuwH/Ds5SciYl/gucCnUkp7WlyWJIFxR1L/jDvSgqh8RSginlX+9wnlv78WEbuB3Smli1JKX4mIDwJ/HhH7AbuAlwBHAf+mzU5LyoNxR1LfjDtSPurcGvfhsb/fXv57EfDk8v+nA28C3gj8NPBV4JSU0mUN+igpX8YdSX0z7kiZqJwIpZRWrP5WtrkXeEX5kKRGjDuS+mbckfLRZrEESZordQYtN9V0cHJXpq1/n4O0+xz0nFtRhL73r9ZmCPujzrHSVWGEptO30a+m+6LP9W2qjeOuq8IZfW2HNoslSJIkSdJcMBGSJEmSlB0TIUmSJEnZMRGSJEmSlB0TIUmSJEnZsWqcJGWs72pVfVZEqqNp5aN50+9+OLPHZc2nviuYNa0E15V5qrg2TZ/7rI35Vl3WtOU33eZdzbcqrwhJkiRJyo6JkCRJkqTsmAhJkiRJyo6JkCRJkqTsWCxBkkZ0NXCz6UDVRdbn+rYxALerwcFdmPVAZM1WV8dan8VF6hzD87a+Qz0P66xvne3QdDvWabu0vlo7rwhJkiRJyo6JkCRJkqTsmAhJkiRJyo6JkCRJkqTsmAhJkiRJyo5V4yRpRN9VfHKrEDfJtG3edNv0WelvGvevxm3avJtt26sdF3WqbHVVzbCLmNi00lgbmsadruJW3eXNk64qxDXhFSFJkiRJ2TERkiRJkpQdEyFJkiRJ2TERkiRJkpQdiyVIUgUOep8/7rPpmg5WrzP90vrKTTWmq8Hlk/ZpV4Px+yyMUKeAgfGhnjaKQzQt/tHFPvOKkCRJkqTsmAhJkiRJyo6JkCRJkqTsmAhJkiRJyo6JkCRJkqTsWDVOUrY2bd7Ntu0PrkLTVeWkrvRVWactQ+7bvOvqWHCfzU6deNTVfupz/09aVp1qZW3E764q1zXVRfW9riq21dkGbezfJseoV4QkSZIkZcdESJIkSVJ2TIQkSZIkZcdESJIkSVJ2LJYgSSP6HhjedFCqA9m1zGMhb00LK0ybvmqM6mIg+0rT91nYpqtCAZp93PKKkCRJkqTsmAhJkiRJyo6JkCRJkqTsmAhJkiRJyo6JkCRJkqTsWDVOkmaoi2pEdaoszbpij2ajaYUxrd0Vl28cZGWxOnGjTv+HGne6Woc6FfWaajrfptugDbM+F7wiJEmSJCk7JkKSJEmSsrNqIhQRz4qIj0bE9RFxb0RcHRFvjoiHjLVbioh3R8StEXF3RFwQEY/truuSFpVxR1LfjDtSfqpcEXolcD9wJnAK8A7gJcCnI2IfgIgIYGv5+suA04D9gG0RcXgH/Za02Iw7kvpm3JEyU6VYwjNSSrtH/r4oIm4H3gs8GbgQ2AL8EnBSSmkbQERcDOwC/hD4D212WtLCM+5M0HSA8RAGKKuergaaeyxMNLO400XRlCHrex26WN60edbZl33G9K5iSVf7sk7RicnPn1lpOateERoLCsu+XP57WPnvFuAfloNCOd2dwMeAUyv1RJJKxh1JfTPuSPlZa7GEE8t/v17+eyxw1YR2O4EjIuLgNS5HkpYZdyT1zbgjLbDaiVBEHAa8HrggpbSjfHoDcMeE5reX/y6tML8XR8SOiNhx6+7763ZHUgbajDvGHElVdBV3Urq7/c5KWpNaiVD5Tcf5wH3A6aMvAWnSJKvNM6V0dkrpuJTScYdsXFenO5Iy0HbcMeZIWk2XcSfioPY6KqmRKsUSAIiI9RSVUo4GTkwp3TTy8u0U35KMW/5mZNK3J5K0oq7jzqRfeHcQ+XAM9Rfp+zSEwdB1NBm0PBRD+rxTZyB6ncHldZbV9BgaQnGHLopRTJtnnfNw3goYTNL0uJumr3WodEUoIvYDPgo8EXhaSunKsSY7Ke6bHXcMcENK6a5GvZSUHeOOpL4Zd6S8VPlB1X2Ac4CnAqemlC6Z0GwrcFhEnDgy3U8Czyhfk6TKjDuS+mbckfJT5da4vwSeDbwJuDsijh957abykvFW4GLgbyLiDyguDb+a4p7Z/9JulyVlwLgjqW/GHSkzVW6N+7Xy3/9EcfKPPl4EkFJ6AHg68Gng7cC5FL/O/JSU0o0t91nS4jPuSOqbcUfKzKpXhFJKR1aZUUrpduC3y4ckrZlxR1LfjDtSfipXjZMkDU/TqlJDNm/9nbUhbK9JfVhaP4OOLIguqp210bZOv7qqOtdnRcWuKvLVaTvr/tapfFfHrKvk1f5BVUmSJEmadyZCkiRJkrJjIiRJkiQpOyZCkiRJkrJjsQRJ2dq0eTfbts9+gPm4OoNlhzBAXout6UBzPViduNPnth/qPu2qX00H+dfRRtGLRSha0cby2uYVIUmSJEnZMRGSJEmSlB0TIUmSJEnZMRGSJEmSlB0TIUmSJEnZsWqcJA3MUKs3KU8ej+264vKNe1XEGsI2rlMVbKiVBLuqNFZnvk0rxE2bflLbriq5DWFfNu3D0vpq7bwiJEmSJCk7JkKSJEmSsmMiJEmSJCk7JkKSJEmSsmOxBElSJV0NzFV3g7E134awr7sqQDDrZQ3VvMXZpoUzZh37vCIkSZIkKTsmQpIkSZKyYyIkSZIkKTsmQpIkSZKyYyIkSZIkKTtWjZMkVTLUqkWLwG2bj02bd7Ntez/7u43jqmpVrzrVztqoFNa0WlkdXVWza7pthqDpNq+zf7uoqOcVIUmSJEnZMRGSJEmSlB0TIUmSJEnZMRGSJEmSlB2LJUiSJGmQ2igeUKewwrwVK5h3Xe1fOLPStF4RkiRJkpQdEyFJkiRJ2TERkiRJkpQdEyFJkiRJ2TERkiRJkpQdq8ZJ0pyoU/lI3XE/qIkrLt+41zFU5/ipU2VrWtuhHq9dbYc6853Utmlls+4qo03WdP82jXF1jrtZV+/zipAkSZKk7JgISZIkScqOiZAkSZKk7JgISZIkScqOxRIkaU4MdYDzPGlj8Lj7QU1s2rybbdvXfgx1VVCgzvLaKNjQtG3T6fs8j7taVlcFBfou3tHUpP4ura82rVeEJEmSJGXHREiSJElSdiolQhFxckRcGBG3RMSeiLgpIj4UEceMtXtYRHwkIu6MiO9FxN9FxBHddF3SIjPuSOqTMUfKT9UxQhuAS4G3A7uBI4BXAZdExGNTStdHxIHAhcAe4IVAAt4IbIuITSmlu1vvvaRFZtyR1CdjjpSZSolQSun9wPtHn4uI/wt8A3gW8FbgDOBo4NEppWvLNlcA3wR+B3hbe92WtOiMO+0b6kDiPi3qeqm5vmLOFZdv7GTQ+KRju6vjvelg+r4LPswT43S/mowRuq3894flv1uAS5YDA0BKaRfwReDUBsuRpGXGHUl9MuZIC6xWIhQR6yJi/4h4JPBO4BbgA+XLxwJXTZhsJ3DMhOclaVXGHUl9MuZI+ah7RehLFPfFXgNsAk5KKX2nfG0DcMeEaW4HlqbNMCJeHBE7ImLHrbvvr9kdSRloNe4YcyStotPPOg4jkoajbiL0AuB44PnA94BPR8SRI6+nCdPESjNMKZ2dUjoupXTcIRvX1eyOpAy0GneMOZJW0elnnYiD2uqnpIZqJUIppa+nlL5UDih8KnAwRUUVKL4h2TBhsiUmf3siSasy7kjqkzFHykfV8tl7SSl9NyKuBR5RPrWT4t7ZcccAX1vrciRpWR9xZ1rFnkWozrMI6zBvmlbM0mwN8bPOtOOnq2OtahWzNirB1ZlH03Wr04dJz3VVza7P+FBnG7RRzW6IFQDXXDUuIg4Ffg74+/KprcDxEXH0SJsjgSeVr0lSI8YdSX0y5kiLrdIVoYg4F7gMuILiftlHAS8H7qOoqw/wLuClwPkR8RqKe2jfANxIUXVFkioz7kjqkzFHyk/VK0KXAM8E3gt8AngFcBHw+JTSNQDlrymfRFFl5X3AOcAuimord7Xcb0mLz7gjqU/GHCkzla4IpZT+FPjTCu1uAE5r2ilJMu5I6pMxR8rPmoslSNIimreB7A7GHzb3hYaoTtzoYuD8tOnrzLdpH4ZwbjYtHtB0HdrYBn1uxy6WteZiCZIkSZI0r0yEJEmSJGXHREiSJElSdkyEJEmSJGXHREiSJElSdqwaJ0kta1rNqM58h1D5aJ5M2zddVMbqSp110HxoGgfaaNvVPOrErTb6O2tN43SdbTPrqnNt9KG7ZZ1ZqZVXhCRJkiRlx0RIkiRJUnZMhCRJkiRlx0RIkiRJUnYipTTrPvzI5icckLZtP2zW3ZAW1lNO+DaXX7onZt2Poegz5gxhgHtuxRa6GsQ7T9tsCMfd0vpdl6aUjuttgQPXVdwZanGPSdo4N4da4KRpH9oo6OL+hbv2nFkp7nhFSJIkSVJ2TIQkSZIkZcdESJIkSVJ2TIQkSZIkZcdESJIkSVJ29p11ByQpB7Ou4jOUPnRhyNXh+qzUl1tVQK1u2v7vorJZ38daV+d9F+psm67iTtNltVGRr85x07RtVV4RkiRJkpQdEyFJkiRJ2TERkiRJkpQdEyFJkiRJ2bFYgiRpbszTAGnodwC5hRHyNuv9P+3cnHW/pqnT36bFIbraNk2LB3RVAKHOfGcd070iJEmSJCk7JkKSJEmSsmMiJEmSJCk7JkKSJEmSsmMiJEmSJCk7Vo2TJM21rqpSzVsVLM2vvo+1LqqVzbr6V1u6qLjW5TyGqE6FuTaq0U2ax9L6lXr4Y14RkiRJkpQdEyFJkiRJ2TERkiRJkpQdEyFJkiRJ2bFYgiT1YBEG3g9hHeZpe0lN1RlMX+fcaDpovenyh6Dp9lKhznasU2Sjr4IcXhGSJEmSlB0TIUmSJEnZMRGSJEmSlB0TIUmSJEnZMRGSJEmSlJ1IKc26Dz+y+QkHpG3bD5t1N6SF9ZQTvs3ll+6JWfdjKJrGnCFUUctNV1W0clOnalhTS+t3XZpSOq6Tmc+hdfscng7c/3fXPH3T6m59a3pc9Vk5r+48ujDU/diVro7nu/acWSnueEVIkiRJUnZMhCRJkiRlx0RIkiRJUnZMhCRJkiRlZ1DFEiJiN3A9cAhw64y705VFXTfXaz48PKW0cdadGIqRmAOLt69HLeq6uV7zwbgzIpO4s6jrBYu7bou2XpXizqASoWURsWNRK8ws6rq5Xpp3i7yvF3XdXC/Nu0Xd14u6XrC467ao67Uab42TJEmSlB0TIUmSJEnZGWoidPasO9ChRV0310vzbpH39aKum+ulebeo+3pR1wsWd90Wdb1WNMgxQpIkSZLUpaFeEZIkSZKkzpgISZIkScrOYBKhiHhYRHwkIu6MiO9FxN9FxBGz7lcdEXF4RPxFRFwcEfdERIqIIye0Wx8Rb4mImyPi3rL9L/ff42oi4lkR8dGIuL7s79UR8eaIeMhYu6WIeHdE3BoRd0fEBRHx2Fn1ezURcXJEXBgRt0TEnoi4KSI+FBHHjLWb+2NTky3CvjXuGHc0XxZh3xp3jDuLYhCJUEQcCFwI/BzwQuAFwCOBbRFx0Cz7VtMjgOcAdwCfX6Hde4AzgNcBTwduBj4ZEY/vvIdr80rgfuBM4BTgHcBLgE9HxD4AERHA1vL1lwGnAfvXryDdAAAEa0lEQVRR7MPDZ9HpCjYAlwIvBX4VeDVwLHBJRDwcFurY1JgF2rfGHeOO5sQC7VvjjnFnMaSUZv4Afo/iwHvEyHNHAfcBr5h1/2qsxz4j/38RkIAjx9o8rnz+9JHn9gWuBrbOeh2mrNfGCc/9ZrkeJ5V/n1r+/ZSRNj8F3A7891mvQ411fXS5Hv+x/Hshjk0fE/f1Quxb445xx8f8PBZl3xp3jDuL8hjEFSFgC3BJSuna5SdSSruAL1IccHMhpfRAhWZbgB8CHxyZ7j7gA8DJEXFAR91bs5TS7glPf7n897Dy3y3AP6SUto1MdyfwMeZoHwK3lf/+sPx3IY5NTbQQ+9a4Y9zRXFmIfWvcMe4siqEkQscCV014fidwzITn59mxwK6U0j1jz+8E9qe43DwPTiz//Xr570r78IiIOLiXXq1BRKyLiP0j4pHAO4FbKAI15HVs5ianfWvcGRjjTrZy2rfGnYEx7uxtKInQBor7TMfdDiz13JeurbSuy68PWkQcBrweuCCltKN8erX1GvJ+/BKwB7gG2ERx+fs75Ws5HZu5yWnfGneGx7iTp5z2rXFneIw7Y4aSCEFxn+K46L0X3QvmeF3LbzrOp7hn9PTRl5jf9XoBcDzwfOB7FIMijxx5fV7XS6vLZd/O8/lp3PmxeVgvrS6XfTvP56dx58fmYb3WbCiJ0B1M/mZgicnZ6Ty7nenruvz6IEXEeopKKUcDJ6eUbhp5ebX1Gux+TCl9PaX0pZTS+4GnAgcDrypfzunYzE1O+9a4MzDGnWzltG+NOwNj3NnbUBKhnRT3Jo47Bvhaz33p2k7gqLJM4ahjgB8A1+49yexFxH7AR4EnAk9LKV051mSlfXhDSumujrvYipTSdyn2wfK9yzkdm7nJad8adwbMuJOVnPatcWfAjDuFoSRCW4HjI+Lo5SfKS3VPKl9bJFsp6s0/e/mJiNgXeC7wqZTSnll1bJqydv45FN8enJpSumRCs63AYRFx4sh0Pwk8gznahxFxKEUN/b8vn8rp2MxNTvvWuDNgxp2s5LRvjTsDZtwpRFknfLadKH6o6avAvcBrKO5RfAPwEGDTvGTXUPwqcfnfpwL/Dvj3wG5gd0rporLNB4CTgT8AdlH8WNfTgRNSSpf13ulVRMQ7KNblTcDHx16+KaV0Uxk8vgA8jGK97qD4wa5NwONSSjf22OVKIuJc4DLgCop7ZR8FvBx4KPDElNI1i3Rs6sEWad8ad4w7mg+LtG+NO8adhTDrHzJafgBHUFyK/B7wj8B5jP041zw8KA6cSY/PjrT5CeBtFGULv09RxePJs+77Cut03QrrddZIuw3A/6S4f/Ye4DMUQWHm6zBlvf6I4peWv1v292qKcpJHjrVbiGPTx8RjYCH2rXHHuONjfh6Lsm+NO8adRXgM4oqQJEmSJPVpKGOEJEmSJKk3JkKSJEmSsmMiJEmSJCk7JkKSJEmSsmMiJEmSJCk7JkKSJEmSsmMiJEmSJCk7JkKSJEmSsvP/ARUDwXcrRjQmAAAAAElFTkSuQmCC\n", "text/plain": [ "