{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Notebook 9: Using Random Forests to classify phases in the Ising Model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "\n", "\n", "## Learning Goal\n", "\n", "The goal of this notebook is to show how one can employ ensemble methods such as Random Forests to classify the states of the 2D Ising model according to their phases. We discuss concepts like decision trees, extreme decision trees, and out-of-bag error. The notebook also introduces the powerful scikit-learn Ensemble class.\n", "\n", "\n", "## Setting up the problem\n", "\n", "The Hamiltonian for 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 of side $L$, and $J$ is some arbitrary interaction energy scale. We adopt periodic boundary conditions. Onsager proved that this model undergoes a 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=1/\\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", "We will use the same basic idea as we did for logistic regression. An interesting question to ask is whether one can train a statistical model to distinguish between the two phases of the Ising model. 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 ensemble methods and in particular Random Forests.\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. 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. \n", "\n", "It is well-known that, near the critical temperature $T_c$, the ferromagnetic correlation length diverges which, among others, 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 random forest and, once the supervised training procedure is complete, we shall evaluate the performance of our classifier on unseen ordered, disordered and critical states. \n", "\n", "A link to the Ising dataset can be found at [https://physics.bu.edu/~pankajm/MLnotebooks.html](https://physics.bu.edu/~pankajm/MLnotebooks.html)." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\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": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "import pickle, os\n", "from urllib.request import urlopen \n", "\n", "# path to data directory (for testing)\n", "#path_to_data=os.path.expanduser('~')+'/Dropbox/MachineLearningReview/Datasets/isingMC/'\n", "\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": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "X_train shape: (104000, 1600)\n", "Y_train shape: (104000,)\n", "\n", "104000 train samples\n", "30000 critical samples\n", "26000 test samples\n" ] } ], "source": [ "###### define ML parameters\n", "from sklearn.model_selection import train_test_split\n", "train_to_test_ratio=0.8 # 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", "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": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAA0IAAAEPCAYAAABvB3TvAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3XnUJGV96PHvj3UEJRnCuFwQR66KB+LIRGKI3sjivYEYBY24XKOJ5qiJid4kJsaIMcHtJCYnJldP9IhZjIr7BsEbF2Qg4gDJADo6KoiZYRPiDIMsAsP23D+qGpqe7vftmlq6up/v55w+73T1U1XPU8tv6umq59eRUkKSJEmScrLbrCsgSZIkSV2zIyRJkiQpO3aEJEmSJGXHjpAkSZKk7NgRkiRJkpQdO0KSJEmSsmNHaM5ExDERkSLimFnXZSAiTo2IxvKwR8QHIuKappYnqboyzpw69P7ZEfHaMeVai0kRsbpc9kubXvauKuPTllnXQ+paRJwbEecOve/N9Uif6jLgtdF8sCMkSRrn54G/H3r/bGCnjhBwSVn2ki4q1QNvBZ4z60pIPZDbua8FtMesK6CdRcTeKaUdHawngD1TSne2vS5J82EQf1JKF05TPqV0MzBV2S61FUdTSt9vepnSPOr63PfaSG3wjlCLIuKEiLggIm6PiJsi4nMRcehImXMj4vyIeFZEXBoRO4DfLj9bFREfiYibI+JHEfFB4CcnrOtXIuLCiLitLPvJiDh4pMyWiPhwRPxGRHwXuBP45fKzfSLiHRGxOSLuLP++MSJ2G1nG2oj4akTcERHXRsSbgJhyewzW/4qIuKJcxiURceyE8oN13RYR34uI3xr5fFVEvC8iLi/LXF1urwNHyj0uIj4bET8s13lVuX32GCpzQES8t2zTjoj4bkS8cpp2SX0VEU8sj/0byjh0WUS8YejzpeLPfY/GRcQHgF8HDiynp8HjYZMeSYmI50TE1yLi1jKG/XtEnDj0+avL+Li9jFkXRsQv72I7Ty3r8NMR8cWIuBX4xNDn08THfcoYcENE3FJut6fEyKN5MebRuIh4RER8MCK2lfFjY0S8eKTMS8tlHRURp5fb5AcR8a6IWLEr7ZbaEhEvLP8f3BERmyJip7ug4879iDg+ItZHcc1zaxlz/nRkPq+Nxq/fa6MZ8I5QSyLiBODzwDnAC4AHA28Bzo+II1JK1w4VfxzwLopHLv4T2F5O/wzwROAU4Hvlct49Zl2/BbwX+KdyHQ8BTgXOi4g1KaVbhoofCxwBvBn4IbClPOi/CBxW1uGbwFHAm4D9gT8o13NA2Z7rKS6KdgCvAx4QVJZxNPAk4I3l/K8H/jUinphSumyo3H7AR4C/Ldv0MuC9EXFZSmldWWZ/4A7gDcBW4L+Vdf1aRDw+pXRHWe4s4EfAq4BtwIHAMyi/CIiI/YCvAQ8qt9tm4PhyfXunlHba5lLfRcSTgXOBK4DfB64BHgusGSk6Kf4MeyuwCvhZYNCZmfjNbES8plzm5yhixa3AzwCrh4qtpnj0bgvF/0XPAs6KiGeklP51mjaOcQbwD8A7gHvLukwbH08Dnld+tgF4OnD6ciuMiH2B84CVFLH6auDFwIciYp+U0mkjs3wI+CjwKxSPFZ0K3Aj82S60V2pcRPxPiv9/P0/xf+oq4P8CewKXLTHfIcCZwKcozrU7KWLOIUNlvDYaz2ujWUkp+WrhRfEf6feAPYamPRq4C3jn0LRzKf7DPmJk/v8FJOCFI9P/tZx+TPn+wcBNwD+OlFtNEYR+b2jaFuA24OEjZV9SLvNpI9PfWC7joeX7t5fvDx4qsy/FCZSm2CZbxsz/EIrg9qGhaR8o63Ps0LS9y/WctsTydwceWc77nHLaAeX7E5eY700UQeOxI9PfX65zj+Xa5stX317Av1FclO+zRJmx8af8LAGnDr3/AHDNmHLHjMSk/YBbgM9UqOtuFJ2hLwFnDE1fXS77pcvMf2pZ7ndHpk8VH4FDy+3wRyPl3jW6/nI7bBl6/+rh9g9NP5vigmr38v1Ly3JvHil3FnD5rI8XX74GL4qL328Duw1N+7ny+D13aNrouX9y+X6/JZbttdHO22TLmPm9Nuro5aNxLSi/IfwZ4OMppbsH01NKmykCzNEjs2xJKX19ZNrPA/cAnx6Z/rEx5fYDTo+IPQYvim9/vws8baT8hSml60emnQBcCawfWcaXKL4BOmpoXRemlK4aatOPgX/ZaSNMNjr/LRTfDv38SLnb0v3fbpCK54K/x8g3LBHxqoj4RhSPwtwNDJY9uM1+A8U3SX9R3nZ+7Jg6nQBcBGweaf8XgZ+i+DZImhsRsQ/wVOD0lNJtyxQfF3/qeArFRcjonZAHiIgnRcRZEfFfFOfuXRQXOYcuNd8yPjvyftr4+HMUj7F8cmT+T02xzqcB16aUzh2Z/mGKb9JH48fnR95/k2rfHEutiYjdKe78fiqldO9gekrpIooL9qV8neI8/lhEnBwRDx1ZttdGk3ltNCN2hNqxkuI/1evGfHY9xW3LYePKPQK4MaV018j0/xp5Pwg0Z1MEoOHXEygO1uXW9VDgUWPm//fy88EyHjFm/ePqtJRJ8x84Mu3GMeV2APc9S18+fvMeirb/CvBk7g9MK6D8Kqa4uNoA/DlweUT8Z0S8ami5D6UIiqPtH1wUjW5Dqe9WUsT3aVKtjosJdQzOl4nrjohHAl+hiIWvoeg8/SzwBYbO8V0w2pZp4+Mjyr8/HJl/mti2/5j1QhHrB58PG330cAfFt7pSHxxAcZFf+f/6lNIVFI9O7UbxCOj1EXFRRAw6OF4bTea10Yw4RqgdN1Lccnz4mM8eTtETHzYuz/x1wMqI2HPkhH/YSLnBsl4KbBqznFtG3o9b1w0Uz34+f8xncP+3QNeNWf+4Oi1l0vzXjpm+nBcCX0kp/cFgQkQ8erRQSuk/gV+LiKB4rvjVwHsiYksqxiLcQHEB9LsT1jPxmWipp26keKxk9D/RcRr7nYvStvLvgcC3JpQ5AfgJ4Pkppfs6TOWdrDpG2zJtfBxcBD2UIhYOTBPbtjP+LtYg/o/Ge6nPtlFc7E76v/rKpWYu71asi4i9Ke5KvwX4fESsxmujpXhtNCPeEWpBeUv0YuB55W1mACLiURTffJ43xWIuoHiu87kj01848n49xQn9mJTShjGvaQ7UL1A8P3rrhGUMLmwuAI4qv80dtGlfikHO0xqd/yEU2VkuqLCMgX0oAvawl00qnApf5/7fQvnp8u8XgMcDV01o/2jAlHqtfBzufODFEfGghha7g2LQ7HLWUyRHWCqz0KDDc9/5GxGPo7hwatK08fEiiguh543MP/p+nPOAgyJitO4voriI+E6N+kudSindA/wHcHIMZUaLiJ/jgclOllvOjpTSOcBfUoyXebTXRkvy2mhGvCPUnjdRPN95VkS8h+KZ+TdTDN776+VmTil9OSLOB95XZiQZZEb56ZFyN0fE64C/i4hVFAMGb6L4NvZoioGNH1lmdadTnCRfiYi/Br4B7AX8d4oMUc8uL6z+hiJ95ZeiSKs7yIxy+3LtGfJfI/O/niJIvrXCMga+ALw+Ik6huFV9HMVgzftExBqKbDcfp8ietTvFN0R3U2R5gaJdLwC+GhF/Q/Etx74UAeAXUkon7ULdpFn7Q4oLiwvK8/oaiuxNR6SUXrMLy/s2sH/56MQG4I6U0jdHC6WUbokiRfe7I+LTFPHlFoqMTHekItPQ2RTn4AfLuj2CIj5eRYNf0E0bH1NKl0XER4C3lhd/F1PEk8GFzL1jFj/wAYpvTD8TEW+k2M6/SvHYyW+WF5bSPPkzinEwn4uI91GMdXsz9z/uOVaZpe1pwP+jSNRyAEXmsh9w/91hr43G89poVqbNquCr+ovi8Y8LKE6GmyhSux46UuZc4PwJ86+iSLN6C0WKww8CJzE+Q9EzgHXAzeX6rgD+EThsqMwW4MMT1rWCIvPSdylOwu0U3wqdygOzu/wM8FWKTCLXUgS1NzN9ZpQPAy8Hvl+u51LguJFyH2B8dqpzeWDGmgdRpMbcWm6jsyiyz9yX7YriUZd/Bi6nyAqzneLi8PiRZa+kOOk3U2Rv+WHZzt9brl2+fPX1BaylGLD7ozIufBd4/dDnS8Wf0axx+5bxaPB4y5Zy+jETYtLJFHdabi/j0kXAM4c+f35ZnzsoHl15ITtnZFtNtaxxY7MYTRkf9ynjyXaKO1pnUnwjm4CThso9oI7ltEdQjInYVsa1jcCLR8q8tFzWY8bVfdbHii9fwy/gf1Nc+O4oz8/njPk/+AHnPsXA/jMoOkE7KB4Z+yQ7X/d4bfTAdWzBa6OZvaJsqNS6KH6E8PyU0ouXKytJs1Z+o/wOYHUayugkSU3x2mi2fDROkpS9iHgmxeM1X6d4FO4XKB4v/ISdIElaTHaEJEkqHiF5NvDHFI8BXkvxg6p/NstKSZLa0+ijcWXGi7+hGCQaFANif89v0yS1xbgjqWvGHWkxNNYRKn//4RsUg7z+hGJQ1tsoBqCuSUXaRElqjHFHUteMO9LiaPLRuFdQpGY9NBW/LkxEbKRIbfibwDuXW8BPHbB7OvhRPq0nteWqK+/mhm33xKzr0aBacWdczNl46aqxZdes3brTtCpl2zKuDl2uv2uTtvk4i7wd5snXL7lzW0pp+h3Xf7XiTsS+abdY2XolmzLteVTl3KyyrraWO0kT6+tq/bNuV1Vt7d9x7k3XThV3mrwj9BVgRUrpqSPTzwNIKR293DLWPmnvtG79ND+ELmlXHPuUa7n04h0L0xGqG3fGxZxH/sQrxpa9+qb37zStStm2jKtDl+vv2qRtPs4ib4d5snLF5otTSkfOuh5NqRt3dt/toLTPXr/TYg2bNe15VOXcrLKutpY7SRPr62r9s25XVW3t33Fu3XHKVHGnsR+uAw7n/h/MGrYJOKzB9UjSgHFHUteMO9KCaLIjtD/FD+2N2k7xg0xjRcQrI2JDRGzYttUf4JZUSeW4Y8yRVFOtuOMQIqk/muwIQTFgcNSSj+GklE5LKR2ZUjrygFW7N1wdSRmoFHeMOZIasMtxJ2LfFqslqYomMxPcSPEtyaiVjP/mRJLqajzuzNu4klmPR5qky3rN2z7T3Jtp3OnDWI9x+tyGWY+lbKu9Vcapdjk+p+uYPG59K1dMN2+Td4Q2UTw3O+ow4NsNrkeSBow7krpm3JEWRJMdoTOBoyLikMGEiFgNPLX8TJKaZtyR1DXjjrQgmuwIvR/YApwRESdFxInAGcDVwPsaXI8kDRh3JHXNuCMtiMY6QuUvKR8HXA58CDgd2Awcl1K6tan1SNKAcUdS14w70uJoMlkCKaWrgOc2uUxJWopxR1LXjDvSYmi0IyRJOZmUGaduJp62llu3Dk1kAuprxispd3XP7yaySrYVO9tYVxVNtLeNrG9V6tVEG7rcv9Nq+neEJEmSJKn37AhJkiRJyo4dIUmSJEnZsSMkSZIkKTsmS5CkhrUxqLWt5bYx+FTzadKx5DHSrDVrt7Ju/QO3aRPbvstkB3W1FQ/7qg9JYtqqQ93/l9pLcHHKVKW8IyRJkiQpO3aEJEmSJGXHjpAkSZKk7NgRkiRJkpQdO0KSJEmSsmPWOEnqQBMZjtrI+tOP7D7Tm6dMUfPGbduNjZeu2un8mLTtx51HVcpOMm4ZdZfbVhu6jp1dZkGrokod6u6zupr4v6KreOQdIUmSJEnZsSMkSZIkKTt2hCRJkiRlx46QJEmSpOyYLEGStJO6A1XbGtCtfpj1AOd5tmbtVtatf+B2amMgPHSb4KStNjRxrHV5XNbdjn04h9pKdjNOE/uxTn29IyRJkiQpO3aEJEmSJGXHjpAkSZKk7NgRkiRJkpQdO0KSJEmSsmPWOEnZ2njpqp2yzfQhY08VdevbdfamuvOPq8O87bNF4Dafb1XOozYyiDWRza5uLKgbz9raXk1s72nb0GV2uKp16CrGeEdIkiRJUnbsCEmSJEnKjh0hSZIkSdmxIyRJkiQpOyZLkJStNWu3sm59NwMymxgQ2uUA9bpJFCbVtcvECtK868M5P64OXQ6y7zKJQxVV1t9EEocqpo3JTRxfdevb9bYZ5R0hSZIkSdmxIyRJkiQpO3aEJEmSJGXHjpAkSZKk7NgRkiRJkpQds8ZJUgf6kP2prnnL2FYlm500S20dl21lNqurzxnXpq3DrLPWVVU3HjbR3rayhtapm3eEJEmSJGXHjpAkSZKk7NgRkiRJkpQdO0KSJEmSsmOyBEmaYyZGkDRJlcHlfR38P6lebSUwqJtUoK+JFfoa05vYv+PKrlwx3fq9IyRJkiQpO3aEJEmSJGVnqo5QRBwUEe+OiAsi4raISBGxeky5FRHxVxFxXUTcXpZ/WtOVlrT4jDuSumbckfIy7R2hxwDPB24EvrpEuX8AXgH8KfBM4DrgixFxRJ1KSsqScUdS14w7UkamTZbwbymlhwFExMuBXxwtEBFPBF4E/EZK6Z/KaecBm4C3ACc2UmNJuTDujKg74Lbur4g38SvkXQ7YrTtoua+Di9WqbONOlcQKXdahrXO2ynLbKjtOW/uhD7Gv221zylTLnOqOUErp3imKnQjcBXx8aL67gY8Bx0fE3lPVSJIw7kjqnnFHykuTyRIOBzanlG4bmb4J2IvidrMkNcm4I6lrxh1pQTTZEdqf4pnaUduHPt9JRLwyIjZExIZtW+9psDqSMlA57hhzJNVk3JEWRJMdoQDShOkTpZROSykdmVI68oBVuzdYHUkZqBx3jDmSajLuSAuiyY7Qdsbf9Vk59LkkNcm4I6lrxh1pQUybNW4am4DnRMQ+I8/NHgbcCVzR4LokCRY07nSZHa4tTWQ+mnWGOWmCWnFn46Wrpj4P2jouqyx32kxubWVs60M2u0VQNyNf18sdp43lNnlH6ExgT+B5gwkRsQfwAuBLKaUdDa5LksC4I6l7xh1pQUx9RygiTi7/+aTy7y9FxFZga0rpvJTS1yPi48DfRsSewGbgVcCjgV9tstKS8mDckdQ1446UjyqPxn1y5P17yr/nAceU/34Z8HbgbcBPAt8ATkgpXVKjjpLyZdyR1DXjjpSJqTtCKaUls7+VZW4HXlu+JKkW446krhl3pHw0mSxBkuZKlUHLddUdnNyWSe3vcpB2l4Oec0uK0PX+1a7pw/6ocqy0lRih7vxN1KvuvuiyvXU1cdy1lTijq+3QZLIESZIkSZoLdoQkSZIkZceOkCRJkqTs2BGSJEmSlB07QpIkSZKyY9Y4ScpY19mqusyIVEXdzEfzptv9cEqH65pPXWcwq5sJri3zlHFtki73WRPLnXZdk9Zfd5u3tdxpeUdIkiRJUnbsCEmSJEnKjh0hSZIkSdmxIyRJkiQpOyZLkKQhbQ3crDtQdZF12d4mBuC2NTi4DbMeiKzZautY6zK5SJVjeN7a29fzsEp7q2yHutuxStmVK6Yr5x0hSZIkSdmxIyRJkiQpO3aEJEmSJGXHjpAkSZKk7NgRkiRJkpQds8ZJ0pCus/jkliFunEnbvO626TLT3yTuX41as3Yr69ZPd1xUybLVVjbDNmJi3UxjTagbd9qKW1XXN0/ayhBXh3eEJEmSJGXHjpAkSZKk7NgRkiRJkpQdO0KSJEmSsmOyBEmagoPe54/7bLK6g9WrzL9yxdRFNaKtweXj9mlbg/G7TIxQJYGB8aGaJpJD1E3+0cY+846QJEmSpOzYEZIkSZKUHTtCkiRJkrJjR0iSJElSduwISZIkScqOWeMkZWvN2q2sW//ALDRtZU5qS1eZdZrS57rNu7aOBffZ7FSJR23tpy73/7h1VclW1kT8bitzXV1tZN9rK2NblW3QxP6tc4x6R0iSJElSduwISZIkScqOHSFJkiRJ2bEjJEmSJCk7JkuQpCFdDwyvOyjVgewa8FjIW93ECpPmnzZGtTGQfan5u0xs01aiAM0+bnlHSJIkSVJ27AhJkiRJyo4dIUmSJEnZsSMkSZIkKTt2hCRJkiRlx6xxkjRDbWQjqpJladYZezQbdTOMaddtvHRVLzOLVYkbVerf17jTVhuqZNSrq+5y626DJsz6XPCOkCRJkqTs2BGSJEmSlJ1lO0IRcXJEfDoiroyI2yPisoj484h4yEi5lRHx9xGxLSJ+HBFnR8QT2qu6pEVl3JHUNeOOlJ9p7gj9IXAPcApwAvBe4FXAlyNiN4CICODM8vPXAM8F9gTWRcRBLdRb0mIz7kjqmnFHysw0yRKelVLaOvT+vIjYDvwzcAxwDnAi8D+A41JK6wAi4gJgM/BHwP9pstKSFp5xZ4y6A4z7MEBZ1bQ10NxjYayZxZ02kqb0WddtaGN9k5ZZZV92GdPbiiVt7csqSSfGTz9lqvUse0doJCgM/Ef598Dy74nADwZBoZzvJuBfgJOmqokklYw7krpm3JHys6vJEo4u/36n/Hs48K0x5TYBB0fEg3dxPZI0YNyR1DXjjrTAKneEIuJA4C3A2SmlDeXk/YEbxxTfXv5ducTyXhkRGyJiw7at91StjqQMNBl3jDmSptFW3Enpx81XVtIuqdQRKr/pOAO4G3jZ8EdAGjfLcstMKZ2WUjoypXTkAat2r1IdSRloOu4YcyQtp824E7FvcxWVVMs0yRIAiIgVFJlSDgGOTildM/TxdopvSUYNvhkZ9+2JJC2p7bgz7hfeHUTeH339Rfou9WEwdBV1Bi33RZ+ud6oMRK8yuLzKuuoeQ31I7tBGMopJy6xyHs5bAoNx6h53k3TVhqnuCEXEnsCngScDz0gpfXOkyCaK52ZHHQZclVK6tVYtJWXHuCOpa8YdKS/T/KDqbsDpwNOBk1JKF44pdiZwYEQcPTTffsCzys8kaWrGHUldM+5I+Znm0bi/A54HvB34cUQcNfTZNeUt4zOBC4APR8TrKG4Nv4Himdm/bLbKkjJg3JHUNeOOlJlpHo37pfLvGylO/uHXywFSSvcCzwS+DLwH+CzFrzMfm1K6uuE6S1p8xh1JXTPuSJlZ9o5QSmn1NAtKKW0HfqN8SdIuM+5I6ppxR8rP1FnjJEn9UzerVJ/NW31nrQ/ba1wdVq6YQUUWRBvZzpooW6VebWWd6zKjYlsZ+aqUnXV9q2S+q2LWWfIq/6CqJEmSJM07O0KSJEmSsmNHSJIkSVJ27AhJkiRJyo7JEiRla83araxbP/sB5qOqDJbtwwB5Lba6A831QFXiTpfbvq/7tK161R3kX0UTSS8WIWlFE+trmneEJEmSJGXHjpAkSZKk7NgRkiRJkpQdO0KSJEmSsmNHSJIkSVJ2zBonST3T1+xNypPHY7M2Xrpqp4xYfdjGVbKC9TWTYFuZxqost26GuEnzjyvbVia3PuzLunVYuWK6ct4RkiRJkpQdO0KSJEmSsmNHSJIkSVJ27AhJkiRJyo7JEiRJU2lrYK7aG4yt+daHfd1WAoJZr6uv5i3O1k2cMevY5x0hSZIkSdmxIyRJkiQpO3aEJEmSJGXHjpAkSZKk7NgRkiRJkpQds8ZJkqbS16xFi8Btm481a7eybn03+7uJ42rarF5Vsp01kSmsbrayKtrKZld32/RB3W1eZf+2kVHPO0KSJEmSsmNHSJIkSVJ27AhJkiRJyo4dIUmSJEnZMVmCJEmSeqmJ5AFVEivMW7KCedfW/oVTpprXO0KSJEmSsmNHSJIkSVJ27AhJkiRJyo4dIUmSJEnZsSMkSZIkKTtmjZOkOVEl85Ha435QHRsvXbXTMVTl+KmSZWtS2b4er21thyrLHVe2bmaz9jKjjVd3/9aNcVWOu1ln7/OOkCRJkqTs2BGSJEmSlB07QpIkSZKyY0dIkiRJUnZMliBJc6KvA5znSRODx90PqmPN2q2sW7/rx1BbCQWqrK+JhA11y9adv8vzuK11tZVQoOvkHXWNq+/KFdPN6x0hSZIkSdmxIyRJkiQpO1N1hCLi+Ig4JyKuj4gdEXFNRHwiIg4bKffIiPhURNwUETdHxGci4uB2qi5pkRl3JHXJmCPlZ9oxQvsDFwPvAbYCBwN/DFwYEU9IKV0ZEfsA5wA7gF8HEvA2YF1ErEkp/bjx2ktaZMYdSV0y5kiZmaojlFL6KPDR4WkR8e/Ad4GTgb8GXgEcAhyaUrqiLLMR+B7wm8A7m6u2pEVn3GleXwcSd2lR26X6uoo5Gy9d1cqg8XHHdlvHe93B9F0nfJgnxulu1RkjdEP5967y74nAhYPAAJBS2gx8DTipxnokacC4I6lLxhxpgVXqCEXE7hGxV0Q8FngfcD3wsfLjw4FvjZltE3DYmOmStCzjjqQuGXOkfFS9I3QRxXOxlwNrgONSSj8sP9sfuHHMPNuBlZMWGBGvjIgNEbFh29Z7KlZHUgYajTvGHEnLaPVax2FEUn9U7Qi9BDgKeBFwM/DliFg99HkaM08stcCU0mkppSNTSkcesGr3itWRlIFG444xR9IyWr3Widi3qXpKqqlSRyil9J2U0kXlgMKnAw+myKgCxTck+4+ZbSXjvz2RpGUZdyR1yZgj5WPa9Nk7SSn9KCKuAB5TTtpE8ezsqMOAb+/qeiRpoIu4MyljzyJk51mENsybuhmzNFt9vNaZdPy0daxNm8WsiUxwVZZRt21V6jBuWlvZ7LqMD1W2QRPZ7PqYAXCXs8ZFxMOAxwPfLyedCRwVEYcMlVkNPLX8TJJqMe5I6pIxR1psU90RiojPApcAGymel30c8PvA3RR59QHeD7waOCMi/oTiGdq3AldTZF2RpKkZdyR1yZgj5WfaO0IXAs8G/hn4PPBa4DzgiJTS5QDlrykfR5Fl5UPA6cBmimwrtzZcb0mLz7gjqUvGHCkzU90RSim9A3jHFOWuAp5bt1KSZNyR1CVjjpSfXU6WIEmLaN4GsjsYv9/cF+qjKnGjjYHzk+avsty6dejDuVk3eUDdNjSxDbrcjm2sa5eTJUiSJEnSvLIjJEmSJCk7doQkSZIkZceOkCRJkqTs2BGSJEmSlB2zxklSw+pmM6qy3D5kPponk/ZNG5mx2lKlDZoPdeNAE2XbWkaVuNVEfWetbpyusm1mnXWuiTq0t65TpirlHSFJkiRJ2bEjJEmSJCk7doQkSZIkZceOkCRJkqTsREpp1nW4z9on7Z3WrT8CtKWGAAAIRklEQVRw1tWQFtaxT7mWSy/eEbOuR190GXP6MMA9t2QLbQ3inadt1ofjbuWKzRenlI7sbIU911bc6Wtyj3GaODf7muCkbh2aSOji/oVbd5wyVdzxjpAkSZKk7NgRkiRJkpQdO0KSJEmSsmNHSJIkSVJ27AhJkiRJys4es66AJOVg1ll8+lKHNvQ5O1yXmfpyywqo5U3a/21kNuv6WGvrvG9DlW3TVtypu64mMvJVOW7qlp2Wd4QkSZIkZceOkCRJkqTs2BGSJEmSlB07QpIkSZKyY7IESdLcmKcB0tDtAHITI+Rt1vt/0rk563pNUqW+dZNDtLVt6iYPaCsBQpXlzjqme0dIkiRJUnbsCEmSJEnKjh0hSZIkSdmxIyRJkiQpO3aEJEmSJGXHrHGSpLnWVlaqecuCpfnV9bHWRrayWWf/akobGdfaXEYfVckw10Q2unHLWLliqRrezztCkiRJkrJjR0iSJElSduwISZIkScqOHSFJkiRJ2TFZgiR1YBEG3vehDfO0vaS6qgymr3Ju1B20Xnf9fVB3e6lQZTtWSbLRVUIO7whJkiRJyo4dIUmSJEnZsSMkSZIkKTt2hCRJkiRlx46QJEmSpOxESmnWdbjP2iftndatP3DW1ZAW1rFPuZZLL94Rs65HX9SNOX3IopabtrJo5aZK1rC6Vq7YfHFK6chWFj6Hdt/toLTPXr+zy/PXze7WtbrHVZeZ86ouow193Y9taet4vnXHKVPFHe8ISZIkScqOHSFJkiRJ2bEjJEmSJCk7doQkSZIkZadXyRIiYitwJXAAsG3G1WnLorbNds2HR6WUVs26En0xFHNg8fb1sEVtm+2aD8adIZnEnUVtFyxu2xatXVPFnV51hAYiYsOiZphZ1LbZLs27Rd7Xi9o226V5t6j7elHbBYvbtkVt13J8NE6SJElSduwISZIkScpOXztCp826Ai1a1LbZLs27Rd7Xi9o226V5t6j7elHbBYvbtkVt15J6OUZIkiRJktrU1ztCkiRJktQaO0KSJEmSstObjlBEPDIiPhURN0XEzRHxmYg4eNb1qiIiDoqId0fEBRFxW0SkiFg9ptyKiPiriLguIm4vyz+t+xpPJyJOjohPR8SVZX0vi4g/j4iHjJRbGRF/HxHbIuLHEXF2RDxhVvVeTkQcHxHnRMT1EbEjIq6JiE9ExGEj5eb+2NR4i7BvjTvGHc2XRdi3xh3jzqLoRUcoIvYBzgEeD/w68BLgscC6iNh3lnWr6DHA84Ebga8uUe4fgFcAfwo8E7gO+GJEHNF6DXfNHwL3AKcAJwDvBV4FfDkidgOIiADOLD9/DfBcYE+KfXjQLCo9hf2Bi4FXA78IvAE4HLgwIh4FC3VsasQC7VvjjnFHc2KB9q1xx7izGFJKM38Bv0tx4D1maNqjgbuB1866fhXasdvQv18OJGD1SJknltNfNjRtD+Ay4MxZt2FCu1aNmfZrZTuOK9+fVL4/dqjMTwDbgXfNug0V2npo2Y4/KN8vxLHpa+y+Xoh9a9wx7vian9ei7FvjjnFnUV69uCMEnAhcmFK6YjAhpbQZ+BrFATcXUkr3TlHsROAu4OND890NfAw4PiL2bql6uyyltHXM5P8o/x5Y/j0R+EFKad3QfDcB/8Ic7UPghvLvXeXfhTg2NdZC7FvjjnFHc2Uh9q1xx7izKPrSEToc+NaY6ZuAw8ZMn2eHA5tTSreNTN8E7EVxu3keHF3+/U75d6l9eHBEPLiTWu2CiNg9IvaKiMcC7wOupwjUkNexmZuc9q1xp2eMO9nKad8ad3rGuLOzvnSE9qd4znTUdmBlx3Vp21JtHXzeaxFxIPAW4OyU0oZy8nLt6vN+vAjYAVwOrKG4/f3D8rOcjs3c5LRvjTv9Y9zJU0771rjTP8adEX3pCEHxnOKo6LwW7QvmuK3lNx1nUDwz+rLhj5jfdr0EOAp4EXAzxaDI1UOfz2u7tLxc9u08n5/GnfvNQ7u0vFz27Tyfn8ad+81Du3ZZXzpCNzL+m4GVjO+dzrPtTG7r4PNeiogVFJlSDgGOTyldM/Txcu3q7X5MKX0npXRRSumjwNOBBwN/XH6c07GZm5z2rXGnZ4w72cpp3xp3esa4s7O+dIQ2UTybOOow4Nsd16Vtm4BHl2kKhx0G3AlcsfMssxcRewKfBp4MPCOl9M2RIkvtw6tSSre2XMVGpJR+RLEPBs8u53Rs5ianfWvc6THjTlZy2rfGnR4z7hT60hE6EzgqIg4ZTChv1T21/GyRnEmRb/55gwkRsQfwAuBLKaUds6rYJGXu/NMpvj04KaV04ZhiZwIHRsTRQ/PtBzyLOdqHEfEwihz63y8n5XRs5ianfWvc6THjTlZy2rfGnR4z7hSizBM+20oUP9T0DeB24E8onlF8K/AQYM289K6h+FXi8p9PB34L+G1gK7A1pXReWeZjwPHA64DNFD/W9UzgKSmlSzqv9DIi4r0UbXk7cNbIx9eklK4pg8f5wCMp2nUjxQ92rQGemFK6usMqTyUiPgtcAmykeFb2ccDvAw8HnpxSunyRjk090CLtW+OOcUfzYZH2rXHHuLMQZv1DRoMXcDDFrcibgVuAzzHy41zz8KI4cMa9zh0q8yDgnRRpC++gyOJxzKzrvkSbtizRrlOHyu0P/CPF87O3AV+hCAozb8OEdr2e4peWf1TW9zKKdJKrR8otxLHpa+wxsBD71rhj3PE1P69F2bfGHePOIrx6cUdIkiRJkrrUlzFCkiRJktQZO0KSJEmSsmNHSJIkSVJ27AhJkiRJyo4dIUmSJEnZsSMkSZIkKTt2hCRJkiRlx46QJEmSpOz8fxTo+00slPksAAAAAElFTkSuQmCC\n", "text/plain": [ "