PY580_Notebook06_Logistics_Regression_2.ipynb

{
 "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": [
       "<Figure size 432x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "##### plot a few Ising states\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "\n",
    "# set colourbar map\n",
    "cmap_args=dict(cmap='plasma_r')\n",
    "\n",
    "# plot states\n",
    "fig, axarr = plt.subplots(nrows=1, ncols=3)\n",
    "\n",
    "axarr[0].imshow(X_ordered[20001].reshape(L,L),**cmap_args)\n",
    "axarr[0].set_title('$\\\\mathrm{ordered\\\\ phase}$',fontsize=16)\n",
    "axarr[0].tick_params(labelsize=16)\n",
    "\n",
    "axarr[1].imshow(X_critical[10001].reshape(L,L),**cmap_args)\n",
    "axarr[1].set_title('$\\\\mathrm{critical\\\\ region}$',fontsize=16)\n",
    "axarr[1].tick_params(labelsize=16)\n",
    "\n",
    "im=axarr[2].imshow(X_disordered[50001].reshape(L,L),**cmap_args)\n",
    "axarr[2].set_title('$\\\\mathrm{disordered\\\\ phase}$',fontsize=16)\n",
    "axarr[2].tick_params(labelsize=16)\n",
    "\n",
    "fig.subplots_adjust(right=2.0)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Cost function, optimizers, regularizers, and performance metrics\n",
    "\n",
    "In Sec. VII of the review, we have shown that the cross-entropy is a natural cost function used for training a logistic regressor. As we already mentioned, minimizing it requires the use of numerical toolboxes. Here, we compare the performance of two different optimization routines: a `liblinear` [the default one for scikit's logistic regression], and stochastic gradient descent (SGD) [see Sec. IV of the review for more details].\n",
    "\n",
    "It is important to note that all these methods have built-in regularizers. Indeed, we did not discuss the role of the regularisor explicitly in the context of Logistic Regression extensively, yet this concept is crucial in order to prevent overfitting, and we encourage the interested reader to play with the different regularization types and regularization strengths and compare model performances. \n",
    "\n",
    "Below, we define the accuracy of a classification model on a given data set as the percentage of correctly classified data points. Comparing the accuracy on the training and test data, we obtain a good estimate of the degree of overfitting. Well-trained models do not overfit the data, which is reflected in an almost equal performance on the training and test data sets [recall that the test set consists of samples which the model has not been trained on]. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Run the cell below (this may take several minutes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "accuracy: train, test, critical\n",
      "liblin: 0.7273, 0.6924, 0.6228\n",
      "SGD: 0.4849, 0.4691, 0.5056\n",
      "finished computing 1/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7273, 0.6924, 0.6228\n",
      "SGD: 0.4992, 0.4781, 0.4980\n",
      "finished computing 2/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7273, 0.6924, 0.6228\n",
      "SGD: 0.4784, 0.4430, 0.5044\n",
      "finished computing 3/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7273, 0.6924, 0.6228\n",
      "SGD: 0.7216, 0.6872, 0.6318\n",
      "finished computing 4/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7273, 0.6924, 0.6228\n",
      "SGD: 0.6758, 0.6474, 0.6767\n",
      "finished computing 5/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7273, 0.6924, 0.6228\n",
      "SGD: 0.5464, 0.5437, 0.6716\n",
      "finished computing 6/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7272, 0.6924, 0.6232\n",
      "SGD: 0.4616, 0.4614, 0.3333\n",
      "finished computing 7/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7266, 0.6917, 0.6245\n",
      "SGD: 0.4616, 0.4614, 0.3333\n",
      "finished computing 8/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7228, 0.6879, 0.6331\n",
      "SGD: 0.4616, 0.4614, 0.3333\n",
      "finished computing 9/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.7031, 0.6711, 0.6610\n",
      "SGD: 0.4616, 0.4614, 0.3333\n",
      "finished computing 10/11 iterations\n",
      "accuracy: train, test, critical\n",
      "liblin: 0.6949, 0.6669, 0.6611\n",
      "SGD: 0.4616, 0.4614, 0.3333\n",
      "finished computing 11/11 iterations\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEQCAYAAACnaJNPAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzsnXl8TFf7wL8nk0wSERGRRGy174ISaqdaFNVWV910US/aH9VF0c3S7S2lC1WtKm83WtW+XV+lxFqVILWrnRAREllknzm/P85MMomEmWSWhPP9fO5n5p577n2ekzu5zz3nPOd5hJQSjUaj0WjsxcvTCmg0Go2mcqENh0aj0WgcQhsOjUaj0TiENhwajUajcQhtODQajUbjENpwaDQajcYhtOHQaDQajUNow6HRaDQah9CGQ6PRaDQOoQ2HRqPRaBzC29MKuIKaNWvKBg0aeFoNh7l48SIBAQGeVsOt6DZfG+g2Vw62bdt2TkoZeqV6V6XhaNCgAbGxsZ5Ww2Gio6Pp06ePp9VwK7rN1wa6zZUDIcRxe+rpoSqNRqPROIQ2HBqNRqNxCG04NBqNRuMQ2nBoNBqNxiG04ShGQgL07g1nzrhf7vjx7d0u1yr7Wmqzp9prlX2ttVlz9aENRzFmzICNG2H6dPfL3bUryO1yrbKvpTZ7qr1W2ddamz35UqRxDeJqTB3bqVMn6ag7rr8/ZGdfWu7jA8uXO0mxErjrLsjLc79cT8quLHKFKP1alztW0vE77ihd9n//W/J5xa9RlmODBkFu7qVyjUZYtw4MBvDyUltJ3690/HJ1rXqMHQsffSQZPVrw4YeX6nK1UkndcbdJKTtdsZ42HIqEBHjqKVixwkVKaTQaQBmW22+H4GCoXr1wK75vLfP3v7KhLomEBLjvPli2DGrVcn47rsTVbDiuygWAZSEiAmrWVD9QHx/1dnjnnTB5sutlv/GGMlje3mby873cJtdWtifa/N13ABIQ3HWXe//WV2rv5d6nrvSuVdrxt96C778vvM933AETJ5Z8XvFrlOWY9fs778BPPxW2ecgQ9ZJkMoHZrDZXfL9wAX7/HQ4cgPx81QupVQvq1FFlFy6o7eLFy/89fXyubFxKKnvzzcLhuWupp+MOtOGwISkJxoyBUaPg44/VG8v117tertms5HbosI0dO6LcJtdWtjvbXHRYUL1KLl8OP/8MWVmule2J9lqR8tL73KWL6+W+++6lbe7f3/VyQcndtw+MRhP5+QaGDr30IZ6bC6mphYYkJaXwe2n7x48Xfi9pKM6W+fPV5utb8nC0xnH0UFUFojJ2bctCQgI895x6+87KAj8/GDZMvRl7YkjB3Vwr9xnUfY2IgA4dYgqMpbOHg7OyihqWw4eVoYiJKTqv5OcHt9yieplDhkBQkHP1KE5lvM96qEpTYYmIUMMWWVlgMJjJzvZi3To19q25urAaiejoi4wc6RoZ/v5qi4hQ+127wqZNsGWLMha5ucpQ1Kun9Pn+ezX8dfPNyrDddpsaptbYj/5X1XiE/fvV5+TJ++nTR/VC2rdXnj4aTXlJTITRo5XxGD1avajMnQvx8cqojBsHe/fCyJGql9uvnxpCS0jwtOaVA204NB6hQweoVg369DnL2rWwYwdUrQo33qgmsM1mT2uoqcysWAHz5kG7durT2vPx8oJu3WDWLDhyBLZtgxdegFOn4Mkn1cR9jx4wZ46aR9GUjDYcGo8QHQ09e6o3QYDISPVPfM89yhNHGw6NqxFCOUa8/rqawN+9G6ZOhfR0eOYZaNAAoqKUN9zBg57WtmKhDYfG7Zw+Df/8A8XnDQMD4auv4JdfwNtbDTds3uwRFTXXGEJA69bwyivw99/KULz1liqfPBmaNVMvN9Omwa5dV3bJvtrRhkPjdry91Vve4MGXHhMCrEnTpkyBXr1g5kzdA9G4lyZN1BDW1q1qyOrdd9X6kGnTlAFp3lwZlNjYa9OIuM1wCCEGCiEOCCEOCSEmlXB8jhAizrL9I4S4YHNshBDioGUb4S6dNa4hLEwZhZYtL19v9mwKFsnddhskJ7tHP43Glvr1Yfx4WL9e9Zbnz4frrlMvNFFR0LChGtratKnwBedqj8/lFsMhhDAA84BbgFbAcCFEK9s6UsoJUsr2Usr2wAfACsu5NYBXgS5AZ+BVIUSwO/TWuIbVq5W//ZUICoJvvoEPPoCVK9WE+u7drtdPoymNWrWUl9aqVWooddEiaNNGTcD36AF166pJ9n/9SwWznDbN0xq7Bnf1ODoDh6SUR6SUucBS4LbL1B8OfG35PgBYJaVMllKmAKuAgS7VVuMyzpxR/vMLF9pXXwgVHmPzZmjUCGrXdq1+Go29hITAo4+qiAdJSfDll8qYfPihCvEipeCjj9Rv2GBQQ1sLF8KaNWr4y2TydAvKjrsWANYBTtrsx6N6EJcghLgOaAisucy5dVygo8YNWNdpOLqgtlMnWLtWfc/NVd4vzz+vYhJpNJ6mWjW4/37o2xeeflpFPM7JUfN5oaFqIeI77xRdyW40Ks+txo0v3Ro2VIsaKyruMhwlxbYsbUrpPmC5lNJqj+06VwgxChgFEB4eTnR0dBnU9CwZGRmVUm9H+OqrpgQEhJOauonoaFmmNsfFBfH22+347LMcXn11Ly1apLtGWRdxLdzn4lxLbc7Obkpubm18fFQwy6io00yYcBCTCZKS/Dh1yo/Tp/0LtoMH/Vi3zp/MzKKP45o1c6hdO4vatbOoUyeL2rWzC/arVcsvUfb580amT2/Fq6/upUaNKwTxKg9SSpdvQFdgpc3+ZGByKXV3AN1s9ocDC2z2FwDDLyevY8eOssycPi1lr15SJiSU/RpllJsSGel+uRbZ7mpzixZSDh5cKLesbd6yRcr69aX08ZHy/felNJsdONlT99gi2yP3+Vpss4e44w4px46V8pNPtsqxY9X+lTCbpTx7Vso//5Tyiy+knDZNyocflrJ7dylr1ZJS+W4VbtWrS3n99VLefbeUkyZJ+cknUq5ZI+WDD0rp5SXlmDFl0x2IlfY80+2pVN4N1bM5ghqCMgJ/A61LqNccOIYl+KKlrAZwFAi2bEeBGpeTVy7DMWZM+f7y5ZBrFsL9ci2y3dHmhAT1i5s5s1Buedp8/ryUQ4aoaz7/vAMneuoeW2R75D5fi232MGvXrnXatTIypNy1S8offpDynXeUYRowQMomTaT09r7UsFg3Pz/H5NhrONwWHVcIMQh4FzAAi6SUrwshplsU/dFSZyrgJ6WcVOzcx4Aplt3XpZSfXU5WmaLjlpYCUAgVRMlVxMWV7AjuarkekG2Wgr3ZjQjZv4kISggKVAa5ZimYffZBBlbbTBv/w5evbG97S8saVJ7y0hz+hVA+nY7gSFajrVtLluvlpQbkren6rJvtflmPWfdnzFCJOIpjNMLOnYXJM3x9HWu/vXg4k5O7ouPm56uoC1OmwIYNah6lShXlyj5rlmNN1xkAHTUcCQkq8tn33yt3B2vWmVat1MyWq8jOVtHWzpxxr1xPyrbKTUhQju9eXiq0qRPkjt05mrbVjjP6ut8ufb7a097S/h/KW56drSI7JiYWtjk8HFq0cOzB6ej/a06Oknv2bKHcmjXV7Ku3t/o7WDdrJqbi3x3dd1RHP78rpwMsbT8oSLWjJMaOhQULlG+sBzI5uTus+pgxKt+K0agcSMrSbHsNh1uGqty9lXmoavRo1Z3383Nvt94iN99odP9wghvbPGWKGoctkAvSbO1TO0FudraUAweqy917r5RpaSVU8tQ9tpHt9vvs7jabzVLm50uZkyPlyJFSenlJk4+PlEKoAf9ff5Xyyy+lnDdPytdfV+OMI0dKedddUt50k5SdOknZuLGUISFSGgylj8NYt6pVpaxXT8q2baXs2VO10RnjNuXEmUNV9mCdW4mLk3bPrRQHO4eqdD4OW6yxmG1TpblR7vYOHYjascO9sZ3d1OYzZ1TU26Ag6DuocFiwoFMwfz589lm5UgD6+qo4V//+N7z0koq4++23KkREAZ66xzay3X6f3d1m68IFgwHOn4fRo9lm2+ZbbrH/WlKq3LL2pAW0fm/ZEo4ehczMwuv4+amkHF9+qRYShYU5v90exjZB1rx5LhZmj3WpbFu5Jsc9iLvfUNzJ0qXqpe+vv6Ty8LnvPvUGan0brFNHyu3bnSYvOlrKiAgpw8OlzMx02mWdwtV8n0vD7W229rJ8fdXvrEkT1YOx/t6uv17KyZOlXLtW9YxcQGW8z9jZ49BBDjVuITpa5du4/nrUfEb16iAEJqNRvaGeOaNiiziJ3r1Vj2PpUuX3IKXr85lrKhDWXtZff6nB/7ZtVdnWrWrSvkoVePtt5SAQEgJDh6rX9IMHr82ohQ6ih6o0bmHdOpV/o2Aes/iwzcGDaowJ4NAhVbFBg3LJDA9XGyibNH8+LF8OeX4J9PngPtaPW0Zko2sgyfm1SGnjNlFRanvpJUhNVeEIVq5U208/qToNG8KAAWq78Ua1LFxTBN3j0LicixeVc08RB5P8fAgN5WKTJuof+/ffVWwGUN4w1shxToqn3rq1iq4bFQWD/z2D1KCN3L9gulOuramkBAXB7berN4rDh1WSmLlz1W/v88+VP2uNGuqN57XXICamcgeYciK6x6FxOQEB6v+y4H8uPR1+/VUNH5TEJ5+oidynnlLhcT/9VCVIKAf9+sHZUf7gnY11xGpPlfmIafMh3w85I4vfDv7G8dTjhAWEEVollNCAUMICwqjhX6Ncsm1JSE9gfNx4VnZaSa2qurdTYRACmjZV25NPKn/WzZtVT+T33+Hll9UWEgI33aR6I/37q1yz1yDacGjchjVNLJs3KytSmo/7ddfB//4HixfDhAnKLeq339TERTlYN3w7A7/qT5ZPvCowG/AzhfHXY9sBWBS3iOV7lxc5p1bVWiQ8q7yQxv82nl1ndymDUiWM0IBQGgc35oHIBwA4mXoSfx9/avjXwEuU3JmfsX4Gu1J3MX3ddD4c7P61BRo7MRrV77NPH3jzTbUOZtWqQkOybJmq16ZNoRHp2bMwMmFCAu3Hj1f1PbD40NVow6FxOX37qv+tSdZ4ANHRag6jWzfV/S8JIVTM6gED1DCBdXV1Xh74+JRJj9q1fMjxPgtSQL4vGHJpnHd7wTzHl8O+5L2B73H24lmSLiaRlJmEWRYOlVXxqUKuKZcdCTtIykziQvYFOkZ0LDAcdyy7g20J2/ASXtSsUpPQKqH0adCHuYPm4v+6P9n5hZEJ5sfOZ37sfPy8/ch6Uc/aV3jCwuCBB9QmpVr5bp0b+eADFfrWz0+93AwYAFu2ELRrl0oZOH++p7V3OtpwaFxKYqKyE0Vc99etU4bAmiP2ctSuXbj89eJFdd6DD6qY6nYYELM089/9/+X2FrfTpEYTwtL6E+JTn1cGj2L6Lx+TnFu4psFoMFI7sDa1A0tO+vHmTW8W2c815ZKRm1Gw/2rvVzl64ShJF5OU8clMItAYCMCRcUdo9H6jAuNRxbsKd7S8g1n9Z135b6CpWAgB7dqpbeJE9btct67QiKxcqaoBfPQRBUk5unRRcybBwSV/Fi8zGh3XzU1hVrTh0LiU9evVZ8GolJTQvXvZPKZyctTQwIsvKveoRYsuG9sq6WISI34YwW+HfuOX+39hUNNBJMz+qeB49NJ57NnjuBpWjAZjkfmPW5vfWmrdiMAIHo58mI+3f4xAkG3KpppvNT3PcTUQEACDBqlt0iTlBvy//6l5Eqt3YMuWyh88MRH27StcsHil617JuBT/fP112LgRpk93aZgVbTg0LqXI+g1Qb14zZ5btYjVqqMny775TE5hRUSqt2ssvX9L7iD4WzQMrHuB85nnmDZrHLU0uXa3cooUaRdi0SdkyV5OUmUSHWh3YeWYnI68fyZmMqzQh9bVMRITqJefnYzIaMeTnq5XqJT3ETSblEpycrLaUlMt/HjxYuH+lRUnz56vNz88lC5i04dC4lOholYu5YP3G6dPK7baM8xQA3Hmn6sJMmKD88KdOLXJ4zp9zeG7VczSp0YRf7v+F9rVK7pU8/rh6Mfv3v+HHH8uujr2suHcFX+36igdWPMCTUU/SNrwUrzJN5SYxkYSxDzG49mZ+Pd2dWqdKeUEwGAp7EY6SnV1oUJKTldviggUkHIjlvtvyWfaLP7X6D1PhcV2ANhwal2E2K2eTDh1sCu+5RxkNax7YshISAv/5j3qb8vJSK8/nzYMpU2gZ2pKHIh9i7qC5VDVWLfUSAQEqIPKrr8Lu3WoUzNUMbjqYLzt/Seuw1q4XpvEMK1bw6k+jiNt+iBdv6sms/p9CVgrSkrhUWlaml7Zvd50ACQHVkfWCoF1DZMwfTInYwobrYHrnLD6sVs1l8xw6rHoFwt1hmN3OxYtqHPaZZ+CttwDntfm3D8Zx8MsPGJfcVM199Ohh13nnzyvv37vvVjEW3cFVf59L4Fpps+9rvuSaXJiy1UEc9dqzN6y67nFoXMaJE2rIt2BU6s8/lTutEx8geaY8XlzzIjOTP6DDXU0Y82EuPr16qcWDb7yhJlguQ0iImjLpdOUMBE5jzdk1xG2J4+kbnnafUI3LkFISfSyamZtnkmvKxSAMCCHIN+fj4+VDZHgkt7e4nWq+1RCWeNDCJlmMPWUCm2MllKXmpPLt3m+JOxNHrinX5V57bjMcQoiBwHuoDIALpZRvlVDnHmAqIIG/pZT3W8pNwC5LtRNSyqFuUVpTLgYNgkaNbOYP1q1T47pOmok+duEY9y2/j79O/cXojqOZPWA2PqNNyuvKGjDx/feveJ0BA5yijt1sOb+FvfF7teGo5JjMJr7f/z1vb3qbmNMxhAWE8fqNr3Mo+RBL/l6C0ctIvsync53OvNTrJZfrczj5MLGnY/Hz9nO5155bDIcQwgDMA24G4oEYIcSPUsq9NnWaApOB7lLKFCGEbcD8LCmli/OoapzJ2bOwZ49aclFAdDR07AiBgeW+fnpOOlGfRJFryuWbu77h7tZ3qwM+wHvvqbGn5s1V2fHjhdniSiEmBp57TuXvcHWqhkZVG7Hq7CqSs5KdGs5E4x6y8rJY8vcSZm2exeGUwzSp0YQFQxbwcLuH8fP2Y9iyYYzuOJoOsgM7xA4SMtyTdyXxYiKjO45mVMdRfLztY5fKdVePozNwSEp5BEAIsRS4DdhrU+cJYJ6UMgVASnnWTbppXMAl6zdArQDPySnXdU1mEwYvA4G+gcwZMIdu9brRKLjRpRWtcxxSqtW+x4+rNKKDBpV43cBAla/5gw9U1G1X0ihA6bsrcRe9G5QvjIrGfSRnJfNhzIe8/9f7JGUm0blOZ96++W1ua34bBi9DQb0V96rIvNHR0YzsM9Jt+lnlAswb7NpMTu6KjlsHOGmzH28ps6UZ0EwIsUkIscUytGXFTwgRaym/3dXKaspPdLTyWurY0aawd2/lZlVGDpw7QKdPOvHzPz8D8GDkgyUbDVuEgNmzVW9j8GB4+GGVd7x3b+WJZaFFCxUode5cFYPRlTQOaAzAzsSdrhWkcQrHLxzn6f89Tf059Xl57ct0rtOZdY+sY8vjWxjWclgRo3Gt4K4ehyihrLg7lzfQFOgD1AU2CCHaSCkvAPWllKeFEI2ANUKIXVLKw0UECDEKGAUQHh5OdHS0k5vgejIyMiql3iXxyy9RtGqVw6ZN6uEYvG0bAClFLIn9bf498Xfm/DMHo5eRvbv3UvX05Se9iyNmz+a6L76g/ldfIZctwysvj9P/+hcHJ0woqHPTTYF8/31HXnjhEPfcE+/Q9R3BmGsk2CeYLXu20Dbr2ljLURl/24cyDrHs5DLWnF2DEIJ+Yf24r959NAxoiPmomXVH1132/MrYZruxJ01geTegK7DSZn8yMLlYnY+AR2z2/wCiSrjWYuCuy8nTqWM9z5o1Uq5fb1PQq5eUnTpdUu9KbU7PSZcjvh8hmYrs9VkveTL1ZNmV8vMrTB1qu/n5FVTp21fK2rWlzM4uu5grsXbtWpmT75p0pRWVyvLbNpvN8o8jf8gBnw+QTEVWfaOqfOZ/z8gTF044fK3K0mZbqGCpY2OApkKIhkIII3AfUHyt7g9AXwAhRE3U0NURIUSwEMLXprw7RedGNBWQvn1VlGlALdLbsqVMbrg/HviR//z9H17p9Qp/PPwHdavVLbtSR47A/fertKGgQmAPHw5HjxZUmTYNXnlFjXC5EqOhDAHsNC4j35zPN3u+IeqTKPr9px9xZ+J448Y3OPH0Cd4Z8A71gup5WsUKhVuGqqSU+UKIp4CVKHfcRVLKPUKI6SgL96PlWH8hxF7ABDwvpTwvhOgGLBBCmFFzMm9JG28sTcXjp5/UZHOBnfjrLxXwzc58GlJKDiUfomlIU4a3GU7bsLbOCc8REaHSgGZng6+vMmgnThRZXduzp43BcyGxp2N5Ze0rzB0098rzNBqXkZmXyeK4xbzz5zscSTlCs5BmfDzkYx5q9xB+3n6eVq/C4rZ1HFLKX4Ffi5W9YvNdAs9YNts6m4FrYyD4KmHyZJUYrcBwREersCB2rOZOzU7liZ+e4NeDv7J77G4aVG/g3JhOllznPPEEDB2qjNrZs0V8cPPyVBLCxo1dt8ZDIPjt0G/sSNihDYcHOJ95nnkx8/hg6wecyzxHlzpdmHXzLIY2H3pNTnY7il45rnEqSUlq/cYDD9gUbt2qAlZVr37Zc2NOxXDv8ns5kXqC1298nfpB9Z2v4IpCl0V+/12lr50yBRYuLCj28oI5c5S6/fu7ZtiqVWgrvIQXOxN3cmerO50vQFMixy4cY86fc1i4YyGZeZkMaTaEid0m0qN+jyIrtzWXx11zHJprhBLXb/z0U4nhZ635txPSE5j952y6L+qOSZrY8OgGXujxQqnpV51GixYwfryKbWUT28xgUPl5YmNhzRrXiPb38adZSDN2nd115cqaMpGQnkDvxb05k3GGuDNxPLDiAZq834T5sfO5p/U97B6zm5+G/0TP63pqo+Eg2nBonEp0tJp7LhL7yWBQOQqKYc2/PWP9DA4lH2JIsyHE/SuOrvW6uk1fXnlFDVN99FGR4oceUlMfb10SGMd5RIZH6rUcLmT6uulsOL6BqE+i6LCgAz8d+IkJN0zgyPgjfHbbZzpCcTnQQ1Uap7Jtm5rKKAhsuHAh7NoF775bMOZTUv5tUJE8g/2D3atwtWqqm9S4cZFiPz+V7uOFF1TPwxVBELvX607SxSTyzfl4e+l/xbIgpeRMxhkOJR8q2N7a9FaRXPHxaWpNTq4pl5n9y5hETFME/WvVOJX161VemQK+/lolnLEZCjgy7ggTVk5g2Z5lQAXIv92smfpMTlZ5ni0RdUePhtWrVaI2VzCuyzjGdRnnmotfRZilmdPpp4sYB9vtYt7FgrreXt5cF3Qd2fnZJGUqo+zv7c+wlsN0fncnog2Hxql4e9s4KOXkwObNMGZMkToRgREkpKsAbN7Cu2Lk305JUUERH3+8YHyqWjU1f64pH9a5rJWdVpZ6j83STHxaPAfPHyw0Cinq83DyYbLyC3NKGA1GGgU3okmNJvRt0JcmNZoUbPWD6uNj8GHMz2P4ePvH+Hn7kWPK8fzv6ypDG45iJKQncN9397HsrmVu/aHZ88/lStnOaPPMmSoxUsG8wNatas1EsfUb+eZ8YhNiCfEPIdgrmO5Nurs1/3aJ7Q0OhiFDVFyrxx+Hpk0L6iclqbzktzshSlrx+9zrs150qt2J2QNml//iV5Drid81FM5lTY2eygvdX+BQ8iEOJh8s0ms4knKEHFNhAEw/bz8aBzemSY0mDGg8oIhxqFet3hVdZt0ZKfZaRBuOYsxYP4ONJzYyfd10PhxcQoJ5F8rdlbrL7XKtsp3R5s8/L5apMjpaDVEVW1H3zZ5vyMzLpFf9Xqw8vJKbfW5m8e2LyyzXUUpt75tvqqxOEybAzz8XFE+frgLrHj2q1qeUV7btfTZLM7Gny5atUkqJRGKWZszSjJQ2323KzdLMi2teZOPxjUxcNZFpfaaRZ84j15RLnsnyeZl9R+rmmgvLv9v3XZG5hgXbFrBg24KC/So+VWhSowktQ1tya7NbixiHOtXqlMurzp2RYq9FdOpYC8UnbK0IhGvWE1g4kXqiSB5hd8l1tmxTRjDxL++g+qC3CbpZPYzHrE6ly+FsHvlXeJG6x1OPl3iNivC3fmJtKi/+eIFHnggjupU/AHnn63L6jWiq9f6U4KFvOlW2VX6If0jBw774Q780g+BpvIQXRoMRHy8f9WnwKdj3MfggEJzJOENqdipmzHh7edMpohMTu0+kS90uRFSNuKrdYCtjulydOtZBjow7wvj/jWfFvhWYpAmDMFA/qD5RtaPw9/F3mdzMvExiT8dyIvWEW+U6W/bxPzsSD3TtmUtYgz4A7BsJ+1DhjovLjTkVw8m0kxXub33kQRMJsT8y/EwIDOqsTmwA63vEcPLPhxn02N/4Vs10iuwG1RsQ4h/C1tNbGdhkIEF+QQgEXsKrYBOicN/2mG355Y6l56bz4/4f2XV2F3nmPIxeRq6vfT0PRT5EaJXQSx749u77ePnYtcLaOtdgFCobXoeIDtzR8g6H/36aioU2HBYiAiMK3vr8vP3INeUysMlAtwwbFfxzWVJNukuurezytnncGthaBX545iWMRlTcDm/vIt5UUkp2nNnB9RHXe6zNdrX3phQigoOxjTaysxG0awdtjn/Iiy+WT7a1zf0b9+fByAfpvqg797a5lyHNhpS5XZfjTPoZ4hLjCtrcoVYHxkaNdYms4ljnGtydDU/jWrThsMFTE2qe/OdyVpsDA2HYMOXNCsDbb6uAT/v2qSi0wNpja+n3n34sv3u5x9psV3uDLWtJ/vlHJYAKDycyUuWBOl7yKJtDsm3b3CasDfe3vZ8Q/5CyX9hOuZ6YKPZUNjyNi7En9npl23Q+jgrATTdJ2bZtkaJ+S/rJWrNqyay8rIKyCtvmlBQpq1aVcsSIgqLcXOdcusK22YXoNlcOqGD5ODRXMTk5KiNSAbm5av2GzcTgX/F/8cfRP3i267OVI1x19erw5JOwZImKoEvhaviDByE/33mipJQkXUxy3gU1GhejDYem3EycCC1b2hiP2FjIzCzznojvAAAgAElEQVRiON7c+CbBfsH8q+O/PKJjmXjxRZXD4//+D8zKiykmRq0T/OYb54l5fcPrRLwTUaJXn0ZTEdGGQ1NuoqOhXj2beXBrnuVevQC4kH2BTSc3Ma7LOAJ9Az2hYtkIDFRzNTExqucBdOwIrVqpRY7O8mRvHtIckzSxL2mfcy6o0bgYtxkOIcRAIcQBIcQhIcSkUurcI4TYK4TYI4T4yqZ8hBDioGUb4S6dNVfm/HnYubNYGPXu3WHqVKhZE4DqftU5Nv4Yz3Z91hMqlo8HHlAG8NQpQOXqeOEFFbfxt9+cIyIyPBJAR8rVVBrc4lUlhDAA84CbgXggRgjxo7RJASuEaApMBrpLKVOEEGGW8hrAq0AnQALbLOemuEN3zeUpMf9G794FYUbSc9Kp4lOFAGOA23VzCkKopByGwjUL990HL72keh2DBpVfRJMaTfDz9tOGQ1NpcFePozNwSEp5REqZCywFbitW5wlgntUgSCnPWsoHAKuklMmWY6uAgW7SW3MFoqOVt21UlKXg1CmIiyuYE5i4aiLXf3w9+WYnzia7G6vRWLsW/vkHHx947jkViuvoUSdc3stAm7A27DyrDYemcuCudRx1gJM2+/FAl2J1mgEIITYBBmCqlPJ/pZxbzohBGmcxeDA0bGizfuPzz1XS8cREEvxNLIpbxCPtHqn8+SbS0uCOO6BLF/jf/3j8ccE990B4+JVPtYfnuz1f+f9GmmsGd/1SSwpIU3xq0RtoiopQURfYIIRoY+e5CCFGAaMAwsPDibZO0FYiMjIyKp3eRiO0b184Hx65YgW+111HzN69fHT4I/JN+fTy7lVquypTm+s8+CBN581j1xtvcL57d0Ctb8zN9cJotD92VEltDkPFoo9OjL70hKuAynSfncVV3WZ7FnuUdwO6Aitt9icDk4vV+Qh4xGb/DyAKGA4ssClfAAy/nDy9ANA9HDwo5Y4dUppMloLcXLVobuxYeT7zvKz6RlV5/3f3X/YalarNublStmolZaNGUmZlSbNZyltukfKhhxy7TEltzs3PlbGnYmV8arxzdK1gVKr77CQqY5upYAsAY4CmQoiGQggjcB/wY7E6PwB9AYQQNVFDV0eAlUB/IUSwECIY6G8p03iYuXOhWzcVlgqA7dshIwN692ZJ3BIycjOY1L1EB7rKiY8PvPceHDkCc+YghFq/8tVX5QtFAspludMnnQqyImo0FRm3GA4pZT7wFOqBvw/4Rkq5RwgxXQgx1FJtJXBeCLEXWAs8L6U8L6VMBmagjE8MMN1SpvEw0dHKcPj6WgrWrVOfvXszrss4okdE0za8rafUcw033QSPPVbgajxhgnLRfeed8l02NCCUiKoR2rNKUylw22yclPJX4NdiZa/YfJfAM5at+LmLgEWu1lFjP8nJav3G9Ok2hWPGQOfOyLAwDELQu0HvUs+v1Hz6acHXunXhoYdg4UJ4+WUIDS37ZSPDI7Xh0FQK9MpxTZnYsEGtnC6yfiMwkJweXemwoANLdy/1lGruwWyGRYvgzz95/nmVIXf+/PJdMjI8kj1Jeyq367LmmkD7/2nKxCXrN/btg2+/ZUmvKvyd+DehVcrx6l0ZyMqCV1+F0FBaxMTw3/8auPHG8l0yMjySXFMu/5z/h1ahrZyjp0bjAnSPQ1Mmpk1T6+EK5jd++YX8aa/y7x1ziaodxY0Ny/kUregEBMDMmbBjB3z6KbfeqorKw82NbuaPh/+gQfUGTlFRo3EV2nBoykS1amotXAHR0Xx7UwRH0o4zpeeUqzqXdAH33gs9e8KUKZCSwm+/Qd++Ksx8WQivGs6NDW+kik8V5+qp0TgZbTg0DrNpk5oUT0uzFJhMyA3reTMqh1ahrRjafOhlz79qEALefx9SUmDqVAwGNYT3xRdlv+S6Y+uu/vkhTaVHGw6NwyxfDm++aTNMFRcHaenMbfgk7w98Hy9xDf2s2rdXVrR/f26+GTp0UJHYTaayXe6T7Z8wcdVE5+qo0TiZa+g/XOMsLlm/cfgwokoVet0ymn6N+nlSNc/w4osweDBCwKRJKlX5Dz+U7VKR4ZGcTDtJSpYO/qypuGjDoXGI5GT4+++CqOkArOsczlPLRpAS7O85xTxNTg68/DJ3+vyXxo3LnuipbZhaMLnr7C4nK6jROA9tODQOUdL6jdc2vMZ3B77H3+caNhwGA/z3vxgmjGf261k8/3zZDIc1qdOuRG04NBUXuwyHEOIpS5wozTXOqVMQEgKdO6v9rWu/YPWR1Txb/178vP08q5wn8fZWE+XHjzP0n1ncc48KReIotQNrU8O/hu5xaCo09v60a6Gy9n1jSQF7Dfhaakpi7FhITAQ/i414c/O/Cc6Cf3Uc7VnFKgJ9+sDdd8Obb5Kx9wSvvabiPjqCEIId/9rB3EFzXaKiRuMM7DIcUsqXULkyPgUeAQ4KId4QQjR2oW6aCoo1Id6es3v4IX834w4GE9iohWeVqijMnAmA78sTmTVLeZ85Sv2g+jqpk6ZCY3dn2hKE8IxlyweCgeVCiLddpJumgvHLL9CxY2G61KreVXhil5H/C3FC4u2rheuug48/xufVKYwdC999p7ysHGH/uf2M/WUsJ1JPuEZHjaac2DvHMU4IsQ14G9gEtJVSjgE6Ane6UD9NBWLNGti7FyIi1P51J1L5+LtcQnrpFPBFePBBiIxk/Hgw+khrJ8Ru0nPSmR87n22nt7lGP42mnNjb46gJDJNSDpBSfiulzAOQUpqBIS7TTlOhiI6Grl3V/MbC7QuJvbAX7ryzWIhcDQBZWYQ/+yCLu8xnyRLlVGAvrcNaIxA6xLqmwmLvHMcrUsoSc5xJKfc5VyVNReTCBRXPr08fSEhP4Klfn+KT9HVqGXndup5Wr+Lh5wcJCdy98yWG9T7PxYv2n1rFpwpNajRh51ltODQVE3uHqpYIIarb7AcLIXRipWsI2/Ubc7bMIc+cx8SmIzytVsVFCHjvPQwZaSxt+jLNmjl2emR4pF7Loamw2DtUFSmlvGDdkVKmAB0cEWRx4z0ghDgkhLgkEbUQ4hEhRJIQIs6yjbQ5ZrIpL56rXOMGgoOVp2mTtsnMj53PvXUG0Lhld9Xj0JRMmzbKf3nBAk79+jcrV9p/amR4JPnmfPJMeVeurNG4GXsNh5ftAkAhRA0cSAIlhDAA84BbgFbAcCFESZlqlkkp21u2hTblWTbl10jo1YpFjx7wzTewcOdcMnIzmJSmVjgXZnLSlMi0aRAczOlHp/DQQ5CZad9pL/V6iSPjj+Bj8HGtfhpNGbD34f8OsFkIYX29vBt43QE5nYFDUsojAEKIpcBtwF4HrqHxEFlZKkZVnTrg7+3Pg5EPEvnFQWjQQLmfakonOBiWL8d0oQVJd8C778KyZe1ZuRJq1Sr9tGsqwrCm0mGX4ZBS/kcIEQvcCAiUh5UjD/06wEmb/XigSwn17hRC9AL+ASZIKa3n+Fnk5wNvSSkviT0qhBgFjAIIDw8nOjraAfUqBhkZGRVS782bQ3jxxbZ88MF2otpEERXUkbzVd3CuWzcOlFPfitpmZyOD9tO6lR//nmEkPSeIf/3rFBMmHLzsOW/se4N6Verx0HUPuUlL13Gt3Gdbruo2Synt2lAL/joDvaybA+feDSy02X8I+KBYnRDA1/J9NLDG5lhty2cj4BjQ+HLyOnbsKCsja9eu9bQKJfLMM1Iaq2TLFbt/liazScqdO6UEKT/7rNzXrqhtdjZBvllyA93lbMbLaHrJcBIkSOnnV/o5XRd2lb0/6+02HV3JtXKfbamMbQZipR3PdHu9qkYC64GVwDTL51QH7FM8UM9mvy5wupgBOy+ltCbd/AS1uNB67LTl8wgQjYMT85rysW4dNBj6H4YtH8KmE5uU++2SJTBggKdVqzTsO+qHuO46/o8P6MkGphum88ADhavwSyIyPJKdiTutL08aTYXB3oHU8UAUcFxK2Rf14E5yQE4M0FQI0VAIYQTuA4p4RwkhImx2hwL7LOXBQghfy/eaQHf03IjbuHABtsflk9T830TVjqJH/R5q3P7hhwuXkGuuSEQjf7of/wpvzHghGWWazxdfCmo1LD0UfWR4JCnZKZxKd2D1oEbjBuw1HNlSymwAIYSvlHI/0NxeIVLKfOApVE9lH/CNlHKPEGK6EMLqJTVOCLFHCPE3MA4VTBGgJRBrKV+LmuPQhsNNbNwIsuW3pIjDTOk5BQHwySdwQsdRcogjR1hf937yvYwA5Asffg2+fJfDmptDryDXVDTs9aqKtywA/AFYJYRIodhQ05WQUv4K/Fqs7BWb75OBySWctxlo64gsjfPoFGWm7vA3CKzWiqHNh6pgVaNGwaJF8Oijnlav8hARQa8h1eDjfKSXFwZzHqdSA9iVVIu2pXhXtQ1ryw11b9AeVpoKxxUNhyX3xjipFgBOFUKsBYKA/7laOY3nyfWLx6tKGpN7vqYeYFYvEdvcsRr7SEyE0aPZ0aIFLf/aSe3lSdx/P2zdCv4ljFgF+QXx5+N/ul9PjeYKXNFwSCmlEOIHLJPVUsp1LtdKUyFITYVfv67PxrsPERFhyd21bh3UqwcNG3pWucrIihUApEVH4/9//4fx7kx63D6fic+P5oO5pedGM0uz7nVoKhT2/hq3CCH0EuFrjO//iGfMUzkcPuijEgtJqXocffqoWEyacnHzua+Zz1jM8z5k06aS68zbOo/gfweTk59TcgWNxgPYazj6oozHYSHETiHELiGEnrG7ynl1x6OIkd3p0sXiDnrkCCQl6WEqZ/Hoo5gGDuJ972fo6h9XYpWaVWqSlpPG/nP73aycRlM69k6O3+JSLTQVjphTMZzwXk3Di2/j72/pXTRuDOfOgdHoWeWuFry8MPxnMbRvD8Pv5exv26jZoCpeNq9ztp5V7Wq184yeGk0x7DUcpcXPnu4sRTQVi+lr34Ss6tzTaHTRAyEhnlHoaiU0FL78EtmvH7+0eo7Utz7i6acLDzcNaYqvwVe75GoqFPYOVV202UyoHkgDF+mk8TB7k/by8+HvYes4BvQJVIVSwvDh8PPPnlXuaqRPH1jwMVu6P8cLL8Dffxce8vbypnVYa53USVOhsDfI4Tu2+0KIWRRb+a25evh619cE+AQQ9/k46odaCg8cgKVL4cYbParb1YoY+Tiv3QY/tpWMvvcCa3YEF7joPtb+MfLN+Z5VUKOxwe6cGsWoggo4qLkKmd53Og+3e5gmtsNS6yxe2Dq/uMsIDYWYjqM5/+sWJj39F+8t8APgyc5PelgzjaYo9gY53GXxptophNgDHADec61qGk+Qk59Derrg6YebsnGjzYHoaBWbqkkTT6l2TVD3ydtox04Gr30Wk6mw/HzmeS5kXyj9RI3Gjdg7xzEEuNWy9UeFOZ/rMq00HiEhPYE6s+vw5o/f8uuvkGfNWiql6nHo9RuuZ9AgzBOepf/BDzH8Vy0YTLqYRM2ZNVkct9izumk0FuwyHFLK4zbbKUvQQs1Vxpwtc0jJTuHcrg4YjXDDDZYDqamqp3HzzR7V71rB6603ICoK06OPM/2xY4T4hxIeEK49qzQVBnuHqpZYghxa94OFEItcp5bG3SRnJTM/dj73tr6Xv9c2oUsXm/hJ1avD+vU6qKG7MBph6VIyfGvy+2fxvPdeYW4OjaYiYO9QVaQlyCEAUsoUdDKlq4q5W+eSkZvBU+0nsW1bsTnwfN3BdDuNGlHt9H5q3taDSZOglohkT9Ie7V2lqRDYazi8hBDB1h0hRA3K7pFVoUlIT6D34t6cyTjjaVXcxrGUY8xYN4ObG91MqDmS7t2hXz/LQSmhaVN45ZXLXkPjfIS3gYULTEwzvs6xRYLs/GwOJR/ytFoajd0P/3eAzUKI5YAE7gHecJlWHmTaumlsPLGR6eum8+HgDz2tjlt4e/PbmKSJGv41aNpUjUoVcPgwHDsGtWu7RZe8vDzi4+PJzs52izxPEBQUxL59++yrLCW3fn8Dt5hMpAbeRX5iPvuS7DzXQ/j5+VG3bl18fHw8rYrGRdi7APA/QohY4EZAAMMczcInhBiIcuE1AAullG8VO/4IMBOw5smcK6VcaDk2AnjJUv6alHKJI7Ltwf91f7LzCx9W82PnMz92PgZh4LUbX2NEuxFEBEaQZ8rD28sbUYm9izJyMziVdorIjyLJNeUWlC/bs4xle5bh5+1H1otZqtCaf8NN6zfi4+MJDAykQYMGlfpvfDnS09MJDAy0/4SGDZF790HVAESzZhXas01Kyfnz54mPj6ehDr1/1WKX4RBCLAHGW11wrZPjUsrH7DzfAMwDbgbigRghxI8lGJ9lUsqnip1bA3gV6ITq7WyznJtij2x7OTLuCM/9/hwr9q0g25SNQRgI8gsCYPIfk7mt+W1EBEbwUexHTFw9kQbVG6gtSH2OiRpDNd9q5OTnYDQYHX7oJSTA+PHtWbkSapWSEe5KSClJvJjIqbRTnEo/VfDZpU4Xbm1+K2cvnqXpB01Jy0krcp6Plw955jz8vauQE3cH0/rOKjy4bh2Eh0NzuzMFl4vs7Oyr2miUCX9/RP16ZJ86Ts7p41QJa0BFfZkXQhASEkJSUpKnVdG4EHuHqi6ZHBdCODI53hk4JKU8AiCEWArcBtjTaxkArJJSJlvOXQUMBL52QP4ViQiMoJpvNXLNufh5+5FryuXe1vfy4eAPycjNwN9buRh1iOjAk1FPcuzCMY5dOEbMqRjOZ50vWN37ytpXmBsz9xLD8my3Z/ESXmTnZ+Nr8L3kwTjp9QR2Xn8fk15bxuK5pVuOTSc2cSL1BPFp8co4pJ+ic+3OPN/9eczSTN3ZdTHJwpVjBmHg2a7PcmvzWwnxD2FEuxHUCaxD3Wp1qVOtDh9v+7igl5Gdn43MqkbH5hb51vwbvXu79S1XG40SqFmT09kJZJjOEXC4No2aGytsx0Pfv6sfew2HlxAi2PqWX4bJ8TrASZv9eKBLCfXuFEL0Av4BJkgpT5Zybp3iJwohRgGjAMLDw4m2DrE4wJ7je7g14laGRAzh54Sf2X1sd4nXGWIcAmGoDcjMzyR2cywAIekhDAobRGJOIgcTDrLh6AYAovJUHqwZe2ew+fxmavnVItQngtjVHZBnW0PYbui4kSVHHmZJ3654VY+nU9/dnMs9R4OABrzc8mUA7t1yL2dzzgLg5+VHqG8ovhm+ROcpPZ9t9iyB3oHUNNakpm9Ngo3BGIShoB3D/IdBPpCstvgz8QVt/vfq9RwOTCAvbz3R0WZEfj51hwzh4nXXkVyGv6c9ZGRkFPkbBwUFkZ6e7hJZ9hIREUFCQgIJCQlMnDiRzz//nC+//JLt27fzzjvv2FX/cphMpjK10eAXSG7eeXIvenHiRDY1auSxYcMGjEYjXbqU9O9UOtu3b+frr79m5syZDuthD9nZ2UXua/H7fC1wNbe5LJPjAHcDrzsgp6RXEFls/yfgaylljhBiNLCEwjmVK52LlPJj4GOATp06yT5lGJNf36dwVngkIx0+H6APl8q9mHuRAGMA6enQ/shZkjMjOXbqGMfEMej0YdEWNlkFTVZhlnAmoyOt6rahV8Mb6NNTXfenpj9R1ViVOoF1qOZb7ZK3u5LkXw7bNn8yYyS1jDBwoE2Fm25y6HqOEh0dje292rdvn2Pj/y4iMDCQwMBAfvjhB0BN+BqNxlJ1K17/cjg8x2HB5GMi4Wwi1YOzyUoCQ82qbN26lapVq3JTCfcpPz8fb++S/8V79+5Nbxcm5PLz86NDh8JBieL3+Vrgam6zvSvH/wPcCSQCZ1GT45d/rSpKPFDPZr8ucLqYjPNSSmt+zE+w5Di359yKiJTwzz+waBHE/hkAQGIizH78Hva/O4f2B77npZo7+KrTaRpevAdyq6gT8/wJjh9Oo+8TOPFiLFsn/Jfnu04GwGSCznU60yq0FUF+QU4dEkhL49L1GzExKuNfBSchQY2mnXGyB/WxY8do06ZNwf7JkycZOHAgzZs3Z9q0aZetv3jxYoYNG8bAgQNp2rQpEydOLKj3xx9/0LVrV66//nruvvtuMjIyAJg+fTpRUVG0adOGUaNGIaV6P+rTpw9TpkxhaP+hLF24lOrVUmjBfnavi+Wjjz5izpw5tG/fng0bNvDII4/wzDPP0LdvX1544QW2bt1Kt27d6NChA926dePAgQOAeqgNGTIEgKlTp/LYY4/Rp08fGjVqxPvvv+/cP6TmqsPedRwACcBW4G+gpmVIyV5igKZCiIZCCCNwH8XCsgshImx2hwJWn8OVQH/LhHwwKlbWSgdkuw0p4d134c471QR38+bw+OPw1VfqeOPG8NdfKoLH2rUwYwYMHxKBzKoBPtn4YATvHKr7VedQXC127IAPPwTrS2PnznDHHbBsGVy86FzdzWZ46y0YNsymMXfeCWPHOleQC5gxAzZuhOkuTiu2detWvvzyS+Li4vj222+JjY29bP24uDiWLVvGrl27WLZsGSdPnuTcuXPMnDmT1atXs337djp16sTs2bMBeOqpp4iJiWH37t1kZWXxs03ukwsXLrBu3TpGjBlBjreZvIAgOoV5MerRx5gwYQJxcXH07NkTgH/++YfVq1fzzjvv0KJFC9avX8+OHTuYPn06U6ZMKVHX/fv3s3LlSrZu3cq0adPIKwhUptFcir1eVSOB8ai3/TjgBuBP1FDSFZFS5gshnkI98A3AIinlHiHEdCBWSvkjME4IMZTCEfhHLOcmCyFmoIwPwHTrRLknychQRmDjRtUTmD5dzR8vWAC5uWq4p0cPtVkdkoRQD//idOiRyKCqo+kgO7BD7CAhIwEhVEbR9u1Vndxc6NVLGY0ffoAqVeDWW2H8eOjatfztqV4dnnvOpuDYMTh5El54ofwXLyNPPw1xJafiBmDDBmXwrMyfrzYvL7A8Qy+hfXtl3MvCzTffTIgl1PywYcPYuHEjnTp1KrV+v379CApSnnmtWrXi+PHjXLhwgf3799O9e3cAcnNz6Wq5gWvXruXtt98mMzOT5ORkWrduza233grAvffeixCCpiFNMRqM+FT1Qu7di9eFFAiuXkTu3XffjcFgACA1NZURI0Zw8OBBhBClGoTBgwfj6+uLr68vYWFhJCYmUrdu3bL9oTRXPfbOcYwHooAtUsq+QogWwKV99csgpfwV+LVY2Ss23ycDk0s5dxHglthYCQlw333qAV2SW+zcubBkCezYoQyGENC3b+HxmBioWtUxmSvuVVFQo6OjGdmn5LkVoxHmzIFZs5SxWrYMvv1WGY+uXZXecXFqSqIsrpp//AHt2kHNmpYCN6/fKAudO8ORIyoNutmsDEbNmqpn5wqKDw1eaajQ19e34LvBYCA/Px8pJX379mX58uVF6mZnZzN27FhiY2OpV68eU6dOLbIIMiBADXdWNRb+uETDhkiTiZzEFHJzJUajKFIX4OWXX6Zv3758//33HDt2rNQx95J01WhKw17DkS2lzBZCIITwlVLuF0K4x7HfzViHPZ55RgWD3bhR9SxiY8HPTz2kAgJg0iT1VnvDDWB5qQQcNxqOYjCo8fzeveH99wvfuJcuVTrXqKFGmO67T9WxvHhelvR0GDAAJk9W7QfU+o2aNaFVK5e15UrY0zMYMwY+/ljdm9xc1fYPXbTgf9WqVSQnJ+Pv788PP/zAokWOv8vccMMNjB07lkOHDtGkSRMyMzOJj48nLEy56NWsWZOMjAyWL1/OXXfddcn5ufm5JGcnU8O/BsbAQALCa3MqMYOjR6FZs0vlpaamUqeOckJcvHixw/pqNCVhr+GIt0TH/QFYJYRIoRJMUDuCvz/YRrn4+mu1AQwdCsnJKurG1KkeUa9EbB1mxo5Vkc+XLlVzKp98ovQ9dMgmym0pbNqkek9FXkY9sH6jLCQmwujRMGqUMiAJCa6T1aNHDx566CEOHTrE/ffff9lhqtIIDQ1l/vz5DB8+nJwc5Qvy2muv0axZM5544gnatm1LgwYNiIqKKvH8XHMu8Wnx+Bn8MPobuf2h+7n99rv4LXo1b7996aT2xIkTGTFiBLNnz+ZGnfZX4yyklA5tQG/U5LXR0XPdtXXs2FE6yunTUt5/v5S+vlKClEajlIMHS3nqlMOXKhunT8uUyEgpExLKfamLF6X89lspp00rLBsxQsrnnpMyNlZKs7lo/SeflFIIKQ8ftik8cEDK3bvLrcuVWLt2bZH9vXv3ulymp0lLSyvzufmmfBlzKkaeSiv8YZrNUibsPS/TY/bKi+kmZ6hYborfx+L3+VqgMrYZNed8xWesI15VVkOzTkr5o5Qy98q1Kw8REVCtmsp65+enIonXr++22H4wYwZBu3Y5xTWoShW4667CgLYmE1y4AO+9B506qWC3L70EFs9Mvv1WOVHNsok0QrNm0Lp1uXXROBeDlwFfgy9ZeVkFZUJAaLiBqlwk9/DJy5yt0TgHhw3H1Yx12GPLFvXp7HUBJeLvr/7z589HSKncgoS48viSAxgMyhMrMRE+/RQaNVKut23aKFFn1UL0AtGjfT6Fb75xmnyNc/H38ScrP6tImaFGEHkhtaiel6TGVTUaF3JV5tQoKytWFH6fN89FQlatUu5PBw6oFYJVqypf2LQ0yMxULlFDhyr3LScTHAyPPaa2s2fh/Hk1Gf7992p+p0oVtU5k7obXYWl7uOcep+ugKT9VfKqQmp2KWZrxEoXvfj7X1YbsdOTx42QbAvAP8r3MVTSasqMNh7M5fRr27FFGwboB/Pab+nzzTbX6LzxcDQfddhscPAgbN2L28cErLw9++kktAHnyybL51tpBWJjagoKUN5KfnzIe14kTeJ84Cs8+7RK5mvITHhBOraq1ihgNQPkjN2qEefdekg9fILRNOEajZ3TUXN1ow1GcKy3kkFK9qtsahuPH4Ysv1DjPpElgDXJXtaoyDjZhK1i8WD2tbX14hw2D0aPZ1qEDUX/8AWvWwIQJajXhu+8qX4YJVFEAACAASURBVFkXUdwrKeyvdeqAC+MYacqHwesyPta+vuQ3b03iP0YyLC66FdwxTlMJ0YajONaFHC+/DE89pQzDgQNqiXZgoJq8tvXJ9fZWK87S0pQxePppNRbUrJmacS/+X1u//qUyLWNkF6OjYeRIZZx++UUZj0cfVVn4nDjnUYJowDI893g0HAmGtm1dIk/jHE6lncLH4ENYQNglx3yrGqlfH84eu0jycRMhDap5QEPN1YyeHLdiM0mN2QwLF6r4FPfco4zIIUuu5/79YfZs9WA/eBCysmD//sIexPXXqwURtWuX/VVPCBgyBHbvht9/V7rl5akZ7bS0K59fHk6fVr0Nr2v3p1HVsorz9OnTBYvwFi9ezFNPPWV3fVeTlpNGSlYK0dHRbN68+ZLjISGSRoYTVDt3hMzUkh0gjx07xlfWQGoajQNcu0+H4hw5AvffD9bQCz4+ap5h9WoVmMoaIrprV9UTGDRIrbgrJWy1U/D1LRzmWrsWpkxRvrSLFhUN0uRMfvtNDdNpqF279iWhQZxZvyxYQ4FYPavWrl1bouEQQuDdtAEGzPiePqp6scXQhkNTVrThsFJ8IYfJpIZr+vVTMUY8Tf/+sHWrMlaPP64CNW3a5BpZlW1G1UVx1StaWPXevXvz3nvvAcqz6sTxEyxYsKBIWPWkpCTuvPNOoqKi6Nq3F38mJmC4mM6ab5YT2bwl7SLb0aFDB9LT05k0aRIbNmygffv2zJkzx6l/O1sS4xKoOuRFzu50h397Ublx1Xu7Xa5V9lXdZntWCVa2rSwrx6WUUt5xh5Rjx0oZF6c+77ijbNcpI3atNDWbpfzqKynr1JGyfftLl4GXh+eeU0vM3YhTVo6PGSOll5f6dAIBAQFSSimPHj0qW7duLaWU8rPPPpO1atWS586dk5mZmbJ169YyJibmsvUbNmwoL1y4ILOysmT9+vXliRMnZFJSkuzWrZvMyMiQUkr51ltvyWmWJf7nz58v0OHBBx+UP/74o5RSyt69e8sxxdqWlp0mY07FyBdefEHOnDmzoHz48OFyw4YNUkopjx8/Llu0aCHNhw/LIT16yA0LF8q0XcdkWlq6zM3Nk2vWrJWDBw+WZrMsdZOy9GOXO75nz96C49Gtxsh8vGR0qzEyL08WbNbj+fmySLm9x62Udjy6dcly7T3fSknHrnS8oM2tS5adn3/58+09bjaXLrcsYOfKcT05botbFnKUEyFg+HC11iMhQe2npMBHH6mJ+fJMov/wA7Rs6Txdy0sFi6vu6bDqtvj7+OPj5VPQK7GyevVq9u7dW7Cfdv48GSdP0r1dO56dM4cHBu5jWN++VA0PRx46BGlpyLi/OWUKJ5FaeJNHa9T5Bm/185JmOGmKIIkwfMmhBfuLHDeb4YSpDuepiT+ZNOMg5nOJ0FoFyLT65/XeOx985iOBRMIJDQWDF2RmwF0Xl/A7A+jHar7gQQDCwwABF9Phlszv2Ex3buMHPmI0ALXC1XUz0qBP1m/E0YEH+ILPeQhRgtws/JjAHKYyFS8vCAtVx1NToHXuDs4QwbPM4jlmYTBAqCVSdGoyNMg7SAaBvMx0xvIhPt5g+Slw4TyE5icAgjy88cZUKHtPoewqZDGPsQxjBb6+hdHwD54LpZVpFwBLeJj+/I6fH1S3TJtuO9uQG+SfAHzHMLqxGX9/CLL4PKxNbEsPNuJPdlG5Qsn1l0UXizoDbTgqKwEBatgK1AN/yhRlPGbNUvFGHJ2YP3VKOQCMGeN8XV2Fm+OqV4Sw6la8vbxpV6sd3/t8X6TcbDbz559/4m95gcjLzCXnUDzPP/IYg3v04OdNm7nhscf48ZMvkQFVwccHc7XqBAg/avuBlxTkpqknmp8vYABzPlT18sXHF7zMXuSmW477AV5gzoNAb198jWAwGcjNqI4pLZ3UWx8ke+0WQjKO4oOJi1QhJrAfgfWqQ0AVQtqCwQfST8Od9WrRMwJCz4RzauvtAIS2V4Yl9SQ80KQmt4RCRHwdTm2/HQHUul61OeU4PNYqmNRgqHe8Ads3P0jooS3Uzj2KNybyMBBd7Ta2PDiPnpm7ObXndow+EGZxHDx/CJ7t4ke2PzQ50IJTB27HzxdCLRF3zh2AF3t5k+8DLfb8f3t3Hh51dTZ8/HtnnQmkgZRA2VQULJAwGSCoWGWrID4utWK1wkWFojQirY/1qUUFUXjwtYCixYKGqtHnArVKQeujVXlLEF4lTYBBEEU2ZRWSQEJWssx5/5jMkGWyTJKZZIb7c11zJb/1nDO/mdw5v+U+iRw7cCudYuCHA13Lc76ChdcBAv/IeI5hm5+nV7mrzZWEczT6Mv5nxiYW9oS+21I4dsxJlzjoWv0xLf7mByyszj+ZkHkVx07GEB8PXS5xzav8JsGzPPb/XcOxvO4kdIO4iwADEQf68syPX2fomv9ifP7bRFFBCTHsuOTnDHh3KX65H7M53ZJge7X4VFU7a1VStIwMY5KTjQFjRo0yZscO37Zfvdq17bZtLa9DC7T6VFVqqus0lcXSZqerGjr11LNnT5OXl2dKSkrMkCFDmjxVdf/993v2eeONN5qNGzeaU6dOmT59+ph9+/YZY4wpLi42e/fuNWfOnDHdu3c3JSUlprCw0CQmJpr58+cbY1ynqtxl1bV06VLz+OOPe6bvuusus3jxYs/0jh07zNld35p969aZqqxs48zKMjeOG2/WrVtnsrOzzahRo1r9fnnjPo6bBqeaSsJMCZZWnULxVXuV255lt0W5+CvJYUuJyEQR2Ssi+0VkTiPr3S4iRkRSqqcvEZFSEXFUv14MVJ2DyujRrkHDX3oJ9uwBLxdvG5WR4bqlODnZL9XzmwAmGHOnVbfb7UyaNKnVadVtNhtXXXUVX3/9NV26dPGkVb/11lsbTKteU15JHgOvHsi6des8F8f//Oc/k52djc1mY/Dgwbz44otIVQVL3l7PkKlTsU2ZijUqihtuuAGbzUZERATJycl+uzgeeeYkWxJT+fDxv7IlMZWo04G5WOwu9/BbWwNabs2yQ7rNzYkurX3hGi72AHApEIVr3PLBXtaLBT4FtgIp1fMuAXb7Ut4F2eOo6cyZ8/ngv/nGmGXLjCkvb3ybZ591XRwPME2r3nK5xbkm61iWKS4vbpP9tSVNqx6cbaaD9TiuAPYbYw4aVzr2N4GfeVlvIbAYKPOyTDVXly7n88G/8YbruRObDf75z4a3efBBWLIkMPVTbSImMgagVop1pQIhUIGjN1BzoICj1fM8RGQo0NcY876X7fuJyA4R2SQiDdwuo7yaN8+VNLGyEm64wfVEujvxolteHlSPRqeCR3RENILUS7GulL8F6q4qb7efeO4jFJEwYBkwzct6J4CLjDF5IjIcWC8iicaYWrk3RGQmMBOgR48eZGRktFHVA6eoqMg/9e7cGVmxgj5//zsXv/46Jx57jAP33+9ZPPCpp+j+r3+x9a23KHffYxggddscFxdHYWFhQOsQaFVVVW3WxqiwKApLCymUjvWelZWV1Tqufvtsd2Ah3ebmnM9q7QsYCXxUY/oR4JEa03FALvBt9asM15jmKV72leFtfs3XBX+NozHff29Mfr7r982bjVm1ypi4ONcdVW30AJ0v9BpH6xw7e8wczj/cZvtrK3qNIzjbTAd7ADALGCAi/YBjwC+Bye6FxpgCoJt7WkQygP8yxmSLSAJw2hhTJSKXAgOAgwGqd+jp0eP872PGuFKruLkfoLNYXMkbVYfXKzZQYxsrdV5ArnEYYyqB2cBHwFfA34wxX4rIAhG5pYnNRwFfiMhO4B0g1RijY2O2hcOH4eqrz2fCjYmBKVPg0KH2rZfyiTEGp/FT0kulvAjYcxzGmA+MMZcbYy4zxiyqnve4MeY9L+uOMcZkV/++1hiTaIxJNsYMM8b8I1B1Dnm9ernutoLzQwD+4AfeB7C6gARDWnW3Smclju8d5BTnNLrePffc40lF8tRTT9VadvXVV7eo7CeeeIKlS5e2aFsV3DQ77oUugA/QBZuOnFbdLVzCEWn8zqqqqir++te/MniwK3dU3cDhLS27Uo3RwHGh+/vfXQkdk5NdP2smegwSJwpPMDp9NN8XXThp1d2Ki4tZ8OACJlw9AZvNxtq1awFXL+jxxx/nyiuv5PPPP2fMmDFkZ2czZ84cSktLsdvtTJkyxbOu2+LFixkyZAjJycnMmeNK8LBq1SpGjBhBcnIykyZNoqSkpNXvrQpuGjhU0Fv46UK2HN7Cgk0L/FrOv//9b1avXo3D4eDtt98mOzu70fUdDgdvvfUWu3bt4q233uLIkSPk5uayZMkSNmzYwPbt20lJSeHZZ58FYPbs2WRlZbF7925KS0t5//3zjzTl5+ezadMmHnrooVplLFy4kK5duvLm/32TnTt3Mm6cKxtecXExSUlJZGZmcs0113jWf/rpp7FarTgcDlavXl1rXx9++CHr168nMzOTnTt3eoLdbbfdRlZWFjt37mTQoEG8/PLLLX8TVUjQ7Liqw/rPf/4nju8bTqu++fDmWheFV2avZGX2SsIkjGsv8v6cqP1Hdp6bGPxp1d02bNjAildX4DROzlWdo2vXroArG++kSZN8at+GDRuYPn06MTGuJ9Lj4+MB2L17N3PnziU/P5+ioiKuv/56n/arQo8GDhW0ruh1BQfPHCS3NBencRImYXSL6cZlXUM/rbqbMYbOUZ3p0rkLYTVOIFgsFsLDw5vdNve+vLVp2rRprF+/nuTkZNLT00P3oTbVbBo4VIfVnJ7Bfe/fR9r2NCwRFsqrypk0aBIrblzhl/p88sknnD59GqvVyvr163nllVd83sdVV13FrFmz2L9/P/3796ekpISjR4/SvXt3ALp160ZRURHvvPNOs+7QmjBhAqteXMVz1YNTnTlzxtPraEhkZCQVFRVERkbW29eCBQuYPHkyMTExnD59mvj4eAoLC+nZsycVFRWsXr2a3r17N7BndaHQaxwqqJ0sPknq8FS2zthK6vDUNr9AXlNHS6sOMHfuXM6cOUNSUhI2m42NGzc2uc3MmTOx2Wyei+NuEydO5JZbbiElJQW73e651XbhwoVceeWVjB8/noEDB/rcZhV6xH3nRihJSUkxTV247IgyMjIYM2ZMe1cjoOq2+auvvmJQRxq+1g8KCwuJjY1t030eOH2AkooShvQY0qb7bam6x1E/28FBRLYZY5r8j0h7HEqFgJjIGM5VnaPKWdX0ykq1kgYOpUKANdI1srSmWFeBoIFDqRBgjagOHDqokwoADRxKhYCo8CjCJZySCn2qW/mf3o6rVAgQES7pcgnREdFNr6xUK2ngUCpEdLU2/vyGUm1FT1UpVceiRYtITEzEZrNht9vJzMwEXJlpH330UQYMGIDdbsdut7No0SLPduHh4djtdhITE0lOTubZZ5/F6aw/TsZ3333HmjVrWlS3xlKgVzorOVN6hoqqihbtW6nm0h6HUjV8/vnnvP/++2zfvp3o6Ghyc3MpLy8HXA/bff/99+zatQuLxUJhYSHPPPOMZ1t38kCAU6dOMXnyZAoKCupl0j18+DBr1qxh8uTJ1FVZWUlERMNfy8ZSoJdVlnHgzAH6x/enS3gXn9qtlC8C1uMQkYkisldE9ovInEbWu11EjIik1Jj3SPV2e0VEM6ypWk6cgNGj22YokRMnTtCtWzdPnqlu3brRq1cvSkpKWLVqFcuXL8disQAQGxvLE0884XU/3bt3Jy0tjRdeeIG6D9nOnz+fzZs3Y7fbWbZsGenp6fziF7/g5ptvZsKECRQVFfHTn/6UYcOGMWTIEN59913Ptu4U6O6Hy26//XYGDhzIlClTsIS76qV3Vil/C0jgEJFw4C/ADcBg4C4RGexlvVjgd0BmjXmDcY1RnghMBFZU708pABYuhC1bYEEbZFWfMGECR44c4fLLL2fWrFls2rQJgP3793PRRRf59MT3pZdeitPp5NSpU7XmP/nkk1x77bU4HA4efPBBwNXTee211/jXv/6FxWJh3bp1bN++nY0bN/LQQw/VCz4AO3bs4LnnnmPPnj0cPHiQrZ9vJTo8Wu+sUn4XqB7HFcB+Y8xBY0w58CbwMy/rLQQWA2U15v0MeNMYc84YcwjYX70/dQEYM6b+a0V1DkOrFURg5UpwOl0/RSAqyrU8N7f+tk3p3Lkz27ZtIy0tjYSEBO68807S09Prrffqq69it9vp27cvR44caXB/zU3pM378eE8ac2MMjz76KDabjeuuu45jx45x8uTJettcccUV9OnTh7CwMOx2O99++y3WSKs+BKj8LlCBozdQ89t1tHqeh4gMBfoaY96ntia3VRem3buhe3cIq/4Uh4W5ppcsad1+w8PDGTNmDE8++SQvvPACa9eupX///hw+fJjCwkIApk+fjsPhIC4ujqoq72k+Dh48SHh4uCfzbWNqpk1fvXo1OTk5bNu2DYfDQY8ePWqlWHfzlrbdGmGlrLLM60V5pdpKoC6Oexu4wPOvmIiEAcuAab5uW2MfM4GZAD169AjKMQOKioqCst6tUbfNcXFxnj/OAP/4h/ftCgtdQeKmm6JJT4/EYoHycrj55gp+/etzFBZCdHT97Wvs2qt9+/YhIvTv3x+AzMxMevbsSVVVFVOnTuU3v/kNzz//PBaLhaqqKsrKyigqKvLU2f0zNzeXe+65h3vvvdczNKxbTEwM+fn5nnXLysooLy/3TJ88eZIuXbpQVlbGxx9/zHfffVevjJKSEiorKz3zysvLKSsrw2qs9OvUj6KioibHC/GnsrKyWsdVP9uhJVCB4yjQt8Z0H+B4jelYIAnIqP6w/wh4T0Ruaca2ABhj0oA0cGXHDbaslBCc2TRby1t2XF+uI5w5A6mpMHMmpKXBiRNRxMZGtbg+xhhmz55Nfn4+ERER9O/fn7S0NGJjY1myZAnz5s1j5MiRxMbGYrVamT59OpdffjlRUVGUlpZy7bXXUlFRQUREBFOnTuX3v/89YWG1O/Y2m43o6GiuueYapk2bRteuXYmKivK0e8aMGdx8882MHTsWu93OwIED6dy5s2d5bGwsMTExREREeOZFRUVhsViIj4tvcdvbksViYejQoZ5p/WyHGGOM31+4AtRBoB8QBewEEhtZPwNIqf49sXr96OrtDwLhjZU3fPhwE4w2btzY3lUIuLpt3rNnT/tUJIDOnj3r1/2fKjpl8kry/FpGU+oeR/1sBwcg2zTjb3pAehzGmEoRmQ18BIQDrxhjvhSRBdUVfa+Rbb8Ukb8Be4BK4H5jjOaOVqoBuSW5hEkY8daO0ftQoSdgDwAaYz4APqgz7/EG1h1TZ3oRsMjbukqp2qyRVvLL8hscQ1yp1tKUI0qFGGuElUpnJRVOTT2i/EMDh1IhJiYyBtAnyJX/aOBQKsS4B3Uqrypv55qoUKVJDpUKMRHhEQz90VDCwzQzj/IP7XEoVUdHTqsO8NRTTzW5jgYN5U8aOJSqoWZa9S+++IINGzbQt6/r+dO5c+dy/Phxdu3ahcPhYPPmzVRUnL8A7U6r/uWXX/LJJ5/wwQcf1EupDufTqrdUcwJHQVkB+/L24TSaekS1PQ0cKvi1YV719kirXlVVxR/+8AdGjBiBzWbjpZde8tRl1KhR2O12kpKS2Lx5M3PmzKG0tBS73c6UKVMabEels5KCcwWcqzzX6vdEqbo0cKjg14Z51dsjrfrLL79MXFwcWVlZZGVlsWrVKg4dOsSaNWu4/vrrcTgc7Ny5E7vdztNPP+3p2axevbrBsq2RrgvkmmJd+YNeHFcdm7dcP3fcAbNmufKq18wau3Kl6xUZ6cp4mJsLt99ee9smks6506pv3ryZjRs3cuedd/L0008zbNiwWuu9+uqrPP/88+Tl5fHZZ595TmfVVbe34c3HH3/MF198wTvvvANAQUEB+/btY8SIEfz617+moqKCW2+9Fbvd3uS+3CwRFgTRFOvKL7THoYKXn/KqBzqtujGG5cuX43A4cDgcHDp0iAkTJjBq1Cg+/fRTevfuzdSpU3n99deb3YYwCcMSYdFnOZRfaI9DdWyN9RAuuwxuu82VFtedV33SJHjgAdfybt2a7GHUtXfvXsLCwhgwYAAADoeDiy++mJiYGGbMmMHs2bN56aWXPGnV3eOR15WTk0NqaiqzZ8+ul/ajc+fOtVLHX3/99axcuZJx48YRGRnJN998Q+/evcnNzaV3797ce++9FBcXs337dn71q18RGRlJRUUFkZGRjbYlNjqWiip9ely1PQ0cKridPFk3r3qrdldUVMRvf/vbemnVwXWb7rx580hKSvKkVb/77rvp1asXgOeidd206nUlJSURERFBcnIy06ZN44EHHuDbb79l2LBhGGNISEhg/fr1ZGRksGTJEiIjI+ncubOnxzFz5kxsNhvDhg1r9DrHRXEXteq9UKoh0pxzsMEmJSXFZGdnt3c1fBbS+fsb4G08jkGDBrVfhQKgsLDQp4vswajucdTPdnAQkW3GmJSm1tNrHEqFqCpnFV+e+pJTxaeaXlkpH2jgUCpEhUkYFc4KvSVXtTkNHEqFKBHBGmHVO6tUmwtY4BCRiSKyV0T2i8gcL8tTRWSXiDhEZIuIDK6ef4mIlFbPd4jIi4Gqs1LBzhpppbSytFnPkyjVXAG5q0pEwoG/AOOBo0CWiLxnjNlTY7U1xpgXq9e/BXgWmFi97IAxpvlPPymlANfYHE7j5FzVOSwRlvaujgoRgepxXAHsN8YcNMaUA28CP6u5gjHmbI3JToD+i6RUK3WK7OQae1y/TaoNBSpw9AaO1Jg+Wj2vFhG5X0QOAIuB39VY1E9EdojIJhG51r9VVRc6f6dVb64XX3zR8+xGeno6x48f9yy755572LNnT0ObelgjrVza9VIska7eRkZGBjfddFOL66QUBO4BQPEyr97/QMaYvwB/EZHJwFzgbuAEcJExJk9EhgPrRSSxTg8FEZkJzATo0aMHGT4+MdwRFBUVBWW9W6Num+Pi4mo9VR1omZmZvPvuu2zatIno6Gjy8vIoLy+nsLCQ+fPnc/LkST777DMsFguFhYUsX77cU1+r1crmzZsB15PjM2bM4NSpUzz22GO1yqiqqmqyjZWVlZ7st4WFhbz88sv069fP8/zHsmXLPMuaYozBiZNwCaekpITKykq/v8dlZWW1jqt+tkOMMcbvL2Ak8FGN6UeARxpZPwwoaGBZBpDSWHnDhw83wWjjxo3tXYWAq9vmPXv2+LyP42ePm1GvjjInCk+0uj5r1641N910U735xcXFJj4+3pw9e7bBbTt16lRr+sCBAyY+Pt44nc5a88+ePWv+9Kc/maSkJGOz2cwf//hHY4wxo0ePNo888ogZNWqUWbp0qZk/f75ZsmSJefvtt02nTp3M5ZdfbpKTk01JSYkZPXq0ycrKMsYY8+GHH5qhQ4cam81mxo0bZ4wxJjMz04wcOdLY7XYzbMQw8+6Wd40xrvf7xhtvbPkb1Ex1j6N+toMDkG2a8Tc9UKeqsoABItJPRKKAXwLv1VxBRAbUmLwR2Fc9P6H64joicikwADgYkFqroLDw04VsObyFBZuCI636xx9/zPr168nMzGTnzp08/PDDnmX5+fls2rSJhx56yDPv9ttvJyUlhdWrV+NwOLBarZ5lOTk53Hvvvaxdu5adO3fy9ttvAzBw4EA+/fRTduzYwR/m/oHnnnqOKqf3ZIxK+Sogp6qMMZUiMhv4CAgHXjHGfCkiC3BFuPeA2SJyHVABnMF1mgpgFLBARCqBKiDVGHM6EPVW7W9M+ph68+5IvINZI2Zh/W8rZVXn06qvzF7JyuyVRIZFUj6vnNySXG7/W+206hnTMhotLxBp1TMyMpg+fToxMTEAxMfHe5bdeeedjdavrq1btzJq1Cj69etXa18FBQXcfffd7Nu3DydOSspKNMW6ajMBS3JojPkA+KDOvMdr/P5AA9utBdb6t3YqGO2+fzdX//VqcktzcRonYRJGt5huPHrNo63arzut+pgxYxgyZAivvfYad9xxhyetemxsLNOnT2f69OkkJSX5nFbdGFMvY65bp06dfKprQ/uaN28eY8eOZd26dezdv5exY8fqg4CqzeiT46pDy5iWUe81a8QsAC7rehm3DboNwPOMwqRBk3jgKtf/IN1iutXbtil79+5l3759nmlvadXLqgePamla9XHjxvHKK69QUuJKBXL6dNMd6NjYWK8XtEeOHMmmTZs4dOhQrX0VFBTQu7frxsU1/+Ma3/x44XEqqyqbLKutnSg8wQOOB/i+qPVD+/pa7uj00QEv1112KLdZA4cKaieLT5I6PJWtM7aSOjy11V+YoqIi7r77bgYPHozNZmPPnj2eccUXLVpEz549SUpKYujQoVx77bVe06onJiZy3XXXMWHCBObPn1+vjPHjx3PLLbeQkpKC3W5n6dKlTdZr2rRppKamYrfbKS0933NISEggLS2N2267jeTkZM+procffphHHnmEn/zkJzidTkSECmcFuaW5rXp/WmLhpwvZVbCrTa5B+VpuW137aknZodxmTavegQRjGubW0rTq/rXtxLYG04380PpDAH4Q/QN+GPNDnMbJd/nf1Vuvi6ULXa1dqXRWcqTgSL3l8dZ44ixxlFeVc+zsMQCOHjzK2P8dS5WpfxovXMKZPGSyZ/q+lPsY2XckX+d+zVObn6q3/oNXPcjQnkPZcWIHy7Yuq7d8zjVzGJwwmM+PfM7K7JWs2bXGa7lR4VHcmVj/GtKfrvsTPWN78uG+D3lj9xv1lj8/8Xm6Wruy7qt1rPt6Xb3lL930EtZIK2/seoOp66Z6LdsSYWHFf6xg47cba82PDo9m1S2rAFiRtYKtR7fWWh4XHcfy/1gOwLOfP4vje0et5T/q/CMWj1+MdZGVssoy6rJEWCh9rPmnKJubVl0HclIqhA3pPoSjZ4+SX5aP07geRgyTMCLCIigqLwLOn+YzR0SvOQAABXZJREFUxnjm1RQTGdPo8s5RnestP1d5jl6xvThdepoKZwXlVeVYwi2Eh4UTb41ny+Etnu0nDZoEQEFZQa35bncnu+6TOVN2xuvygrICAHJKcthyeIun3JKKEgwGS7iFSYMn8ZO+P2HJZ/WHFXbfNHCs8JjX/ZdXuU5HHi447HW5O1Acyj9Ur2xBuGvIXTwz4RmWZy6vt737vQX4Ju+bessTOiV4ft+Ts6fe8ku6XALAwd8dZOTLIzlccBiDISYihp8P+jlLJzTdm20J7XF0INrj0B6HP3yX/x05JTmIiGuEwZgELu5ysV/LdB/H+96/j7TtaURIBJWmkt8M/w0rblzh17IBT7lR4VGUV5UHrNyaZQdjm7XHoZQCoMJZQUJMAgmdEsgpzqHCGbhxyN3XoIaaoeyQHZwoat3Qvr6WO3P4TNK2pQWs3Jplh3KbtcfRgWiPw/Wf6sCBAxu8XTUUhPrQscYYvv76ax06NgjbrEPHqqBksVjIy8vT8SOClDGGvLw8LBZN4R7K9FSV6lD69OnD0aNHycnJae+q+E1ZWVlI/2G1WCz06dOnvauh/EgDh+pQIiMjPekzQlVGRgZDhw5t72oo1WJ6qkoppZRPNHAopZTyiQYOpZRSPgnJ23FFJAeonzuh4+sGBD6ZUPvSNl8YtM3B4WJjTEJTK4Vk4AhWIpLdnHuoQ4m2+cKgbQ4teqpKKaWUTzRwKKWU8okGjo4lrb0r0A60zRcGbXMI0WscSimlfKI9DqWUUj7RwKGUUsonGjiUUkr5RANHEBGRTiKyTURuau+6BIKI3Coiq0TkXRGZ0N718YfqY/padTuntHd9AuFCOK7ehNL3VwNHAIjIKyJySkR215k/UUT2ish+EZnTjF39Efibf2rZttqizcaY9caYe4FpwJ1+rG6b8rHttwHvVLfzloBXto340uZgPa51teAzHjTf36Zo4AiMdGBizRkiEg78BbgBGAzcJSKDRWSIiLxf59VdRK4D9gAnA135FkqnlW2usenc6u2CRTrNbDvQBzhSvVpVAOvY1tJpfpvdgu241pVO8z/jwfb9bZSOxxEAxphPReSSOrOvAPYbYw4CiMibwM+MMf8HqNeVFZGxQCdcH8ZSEfnAGOP0a8VboY3aLMDTwIfGmO3+rXHb8aXtwFFcwcNBEP8j50ubReQrgvC41uXjce5MEH1/m6KBo/305vx/muD6A3JlQysbYx4DEJFpQG6Qfuh8ajPwW+A6IE5E+htjXvRn5fysobb/GXhBRG4E/tEeFfOjhtocSse1Lq9tNsbMhqD//npo4Gg/4mVek09jGmPS274qAeNTm40xf8b1hzUUeG27MaYYmB7oygRIQ20OpeNaV6Of8SD//noEbdc4BBwF+taY7gMcb6e6BMqF2Ga3C7Ht2uYQbbMGjvaTBQwQkX4iEgX8Enivnevkbxdim90uxLZrm0O0zRo4AkBE3gA+B34sIkdFZIYxphKYDXwEfAX8zRjzZXvWsy1diG12uxDbrm2+MNrspkkOlVJK+UR7HEoppXyigUMppZRPNHAopZTyiQYOpZRSPtHAoZRSyicaOJRSSvlEA4dSSimfaOBQSinlEw0cSgVI9bgj34nIfe1dF6VaQwOHUgFijNmFK3fRr9q7Lkq1hgYOpQLrFJDY3pVQqjU0cCgVWE8D0SJycXtXRKmW0sChVICIyERcw4f+L9rrUEFMA4dSASAiFmAxMAvYBSS1b42UajkNHEoFxlzgdWPMt2jgUEFOA4dSfiYiPwbGA89Vz9LAoYKaDuSklFLKJ9rjUEop5RMNHEoppXyigUMppZRPNHAopZTyiQYOpZRSPtHAoZRSyicaOJRSSvlEA4dSSimf/H/vhm2SzOqhgQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "###### apply logistic regression\n",
    "from sklearn import linear_model\n",
    "from sklearn.neural_network import MLPClassifier\n",
    "\n",
    "\n",
    "# define regularisation parameter\n",
    "lmbdas=np.logspace(-5,5,11)\n",
    "\n",
    "# preallocate data\n",
    "train_accuracy=np.zeros(lmbdas.shape,np.float64)\n",
    "test_accuracy=np.zeros(lmbdas.shape,np.float64)\n",
    "critical_accuracy=np.zeros(lmbdas.shape,np.float64)\n",
    "\n",
    "train_accuracy_SGD=np.zeros(lmbdas.shape,np.float64)\n",
    "test_accuracy_SGD=np.zeros(lmbdas.shape,np.float64)\n",
    "critical_accuracy_SGD=np.zeros(lmbdas.shape,np.float64)\n",
    "\n",
    "# loop over regularisation strength\n",
    "for i,lmbda in enumerate(lmbdas):\n",
    "\n",
    "    # define logistic regressor\n",
    "    logreg=linear_model.LogisticRegression(C=1.0/lmbda,random_state=1,verbose=0,max_iter=1E3,tol=1E-5,\n",
    "                                           solver='liblinear')\n",
    "\n",
    "    # fit training data\n",
    "    logreg.fit(X_train, Y_train)\n",
    "\n",
    "    # check accuracy\n",
    "    train_accuracy[i]=logreg.score(X_train,Y_train)\n",
    "    test_accuracy[i]=logreg.score(X_test,Y_test)\n",
    "    critical_accuracy[i]=logreg.score(X_critical,Y_critical)\n",
    "    \n",
    "    print('accuracy: train, test, critical')\n",
    "    print('liblin: %0.4f, %0.4f, %0.4f' %(train_accuracy[i],test_accuracy[i],critical_accuracy[i]) )\n",
    "\n",
    "    # define SGD-based logistic regression\n",
    "    logreg_SGD = linear_model.SGDClassifier(loss='log', penalty='l2', alpha=lmbda, max_iter=100, \n",
    "                                           shuffle=True, random_state=1, learning_rate='optimal')\n",
    "\n",
    "    # fit training data\n",
    "    logreg_SGD.fit(X_train,Y_train)\n",
    "\n",
    "    # check accuracy\n",
    "    train_accuracy_SGD[i]=logreg_SGD.score(X_train,Y_train)\n",
    "    test_accuracy_SGD[i]=logreg_SGD.score(X_test,Y_test)\n",
    "    critical_accuracy_SGD[i]=logreg_SGD.score(X_critical,Y_critical)\n",
    "    \n",
    "    print('SGD: %0.4f, %0.4f, %0.4f' %(train_accuracy_SGD[i],test_accuracy_SGD[i],critical_accuracy_SGD[i]) )\n",
    "\n",
    "    print('finished computing %i/11 iterations' %(i+1))\n",
    "\n",
    "# plot accuracy against regularisation strength\n",
    "plt.semilogx(lmbdas,train_accuracy,'*-b',label='liblinear train')\n",
    "plt.semilogx(lmbdas,test_accuracy,'*-r',label='liblinear test')\n",
    "plt.semilogx(lmbdas,critical_accuracy,'*-g',label='liblinear critical')\n",
    "\n",
    "plt.semilogx(lmbdas,train_accuracy_SGD,'*--b',label='SGD train')\n",
    "plt.semilogx(lmbdas,test_accuracy_SGD,'*--r',label='SGD test')\n",
    "plt.semilogx(lmbdas,critical_accuracy_SGD,'*--g',label='SGD critical')\n",
    "\n",
    "plt.xlabel('$\\\\lambda$')\n",
    "plt.ylabel('$\\\\mathrm{accuracy}$')\n",
    "\n",
    "plt.grid()\n",
    "plt.legend()\n",
    "\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "## Interpreting the results\n",
    "\n",
    "The first thing we can read off the figure above is the relative degree of overfitting. This information is contained in the difference in accuracy of our model on the training (blue) and test (red) datasets. Notice that the accuracy difference between test and training sets is significant but not unreasonable, within $10\\%$. Interestingly, which optimizer performs better depends on the value of the regularization strength. Moreover, similar to the Linear Regression examples, we find that there exists a sweet spot for the regularization strength $\\lambda$ that results in optimal performance of the logistic regressor, at about $\\lambda\\sim 10^{-1}$.\n",
    "\n",
    "Due to the physics of the Ising model close to criticality, we expect that predicting the phase of a sample will become much more difficult close to the critical point. We can visually see this by looking at the states in the critical region, (see Fig. above and plot other examples). Notice that it is no longer easy even for a trained human eye to distinguish between the ferromagnetic and the disordered phases close to $T_c$. \n",
    "\n",
    "It is an interesting exercise to compare the training and test accuracies in the ordered and disordered phases to the accuracy of the model near the critical point (i.e. critical states). Recall that the model is not trained on critical states. Notice that the accuracy is about $10\\%$ smaller for the critical states (green curves). \n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises: ###  \n",
    "<ul>\n",
    "<li> Change the regularization to $L^1$, or $L^1+L^2$ (i.e. elastic net, see Notebook 2 for more details), and study the performance of the model.\n",
    "\n",
    "<li> Try out different solvers supported by `LogisticRegression()`, [see online scikit documentation](http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html). How do they compare to `liblinear` and `SGD`?\n",
    "\n",
    "<li> One can expect the regressor to be maximally confused exactly at the critical point. Take the data states for the closest temperature, i.e. $T/J=2.5$, and check the performance of the model. Do we get an accuracy of about $50\\%$? Does this depend on the proportion of ordered to disordered states in the training set?\n",
    "\n",
    "<li> Take the regressor trained on the square lattice Ising model, and try to predict the phases on a different lattice geometry (e.g. triangular or honeycomb). What accuracy do you obtain? Can you roughly locate the critical point in this new model? Note that this requires you to generate your own dataset using Monte-Carlo sampling.\n",
    "</ul>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python mlreview",
   "language": "python",
   "name": "mlreview"
  },
  "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}