diff --git a/examples/ensemble/ExtraTrees_vs_RF_Gridsearch.ipynb b/examples/ensemble/ExtraTrees_vs_RF_Gridsearch.ipynb new file mode 100644 index 0000000000000..e417dcda0d6d8 --- /dev/null +++ b/examples/ensemble/ExtraTrees_vs_RF_Gridsearch.ipynb @@ -0,0 +1,417 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Demonstration of warmstart grid search to compare classifier performance\n", + "\n", + "An important step in classifier performance comparison is hyperparameter \n", + "optimization. Here, we specify the classifer models we want to tune and a \n", + "dictionary of hyperparameter ranges (preferably similar for fairness in \n", + "comparision) for each classifier. Then, we find the optimal hyperparameters \n", + "through a function that implements warmstart grid search and refit the optimized \n", + "models to obtain accuracies. The performance of each hyperparameter value pairing is visualized in heatmaps.\n", + "\n", + "In this example, we tune hyperparameters for two classifiers, Random Forest and Extra Trees, and compare their performance on an OpenML-CC18 benchmarking suite dataset (https://www.openml.org/d/15). We can see clearly in the resulting plot that the optimized models perform better than or atleast similar to the default parameter models. On the dataset we use in this example, RF performs marginally better than ExtraTrees overall.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Automatically created module for IPython interactive environment\n" + ] + } + ], + "source": [ + "print(__doc__)\n", + "\n", + "import pandas as pd\n", + "import numpy as np\n", + "\n", + "import sklearn\n", + "from sklearn import metrics\n", + "from sklearn.metrics import cohen_kappa_score, make_scorer\n", + "from sklearn.model_selection import cross_val_score\n", + "from sklearn.model_selection import train_test_split\n", + "from sklearn.ensemble import ExtraTreesClassifier\n", + "from sklearn.ensemble import RandomForestClassifier\n", + "from sklearn.datasets import fetch_openml\n", + "\n", + "import matplotlib\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from warnings import simplefilter\n", + "\n", + "simplefilter(action=\"ignore\", category=FutureWarning)\n", + "from warnings import simplefilter\n", + "\n", + "simplefilter(action=\"ignore\", category=FutureWarning)" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "metadata": {}, + "outputs": [], + "source": [ + "def hyperparameter_optimization(X, y, *argv):\n", + " \"\"\"\n", + " Given any number of classifier types and \n", + " a dictionary of two hyperparameters to tune for each classifier, \n", + " find optimal pairs of hyperparameters.\n", + "\n", + " Parameters\n", + " ----------\n", + " X : numpy.ndarray\n", + " Input data, shape (n_samples, n_features)\n", + " y : numpy.ndarray\n", + " Output data, shape (n_samples, n_outputs)\n", + " *argv : list of tuples (classifier, hyperparameters)\n", + " List of (classifier, hyperparameters) tuples:\n", + "\n", + " classifier : sklearn-compliant classifier\n", + " For example sklearn.ensemble.RandomForestRegressor\n", + " hyperparameters : dictionary of hyperparameter values or range\n", + "\n", + " Returns\n", + " -------\n", + " clf_best_params : dictionary\n", + " Dictionary of classifiers and their respective optimal hyperparameters\n", + " \"\"\"\n", + "\n", + " best_models = {}\n", + "\n", + " # Iterate over all (classifier, hyperparameter dict) pairs\n", + " for clf, params in argv:\n", + " best_params = grid_search(X, y, clf, params)\n", + " best_models[clf.__class__.__name__] = best_params\n", + "\n", + " return best_models\n", + "\n", + "\n", + "def grid_search(X, y, clf, params):\n", + " \"\"\"\n", + " Given a classifier and two hyperparameters and the \n", + " range/values to search for each, find optimal hyperparameter \n", + " values using warmstart grid search parameter sweeps.\n", + "\n", + " Parameters\n", + " ----------\n", + " X : numpy.ndarray\n", + " Input data, shape (n_samples, n_features)\n", + " y : numpy.ndarray\n", + " Output data, shape (n_samples, n_outputs)\n", + " clf : sklearn-compliant classifier\n", + " For example sklearn.ensemble.RandomForestRegressor\n", + " params : dictionary of hyperparameter values or range\n", + "\n", + " Returns\n", + " -------\n", + " best_params : dictionary\n", + " Dictionary of best hyperparameters\n", + " \"\"\"\n", + " param1_name = list(params.keys())[0]\n", + " param2_name = list(params.keys())[1]\n", + " param1 = params[param1_name]\n", + " param2 = params[param2_name]\n", + "\n", + " # sweep over all pairs of parameter combinations and collect mean scores\n", + " kappa_scorer = make_scorer(cohen_kappa_score)\n", + " mean_scores = np.zeros((np.shape(param1)[0], np.shape(param2)[0]))\n", + " for idx1, val1 in enumerate(param1):\n", + " clf.max_features = val1 #change .max_features to .name of 1st parameter\n", + " for idx2, val2 in enumerate(param2):\n", + " clf.n_estimators = val2 #change .n_estimators to .name of 2nd parameter \n", + " score = cross_val_score(clf, X, y, scoring=kappa_scorer, cv=5)\n", + " mean_scores[idx1][idx2] = np.mean(score)\n", + "\n", + " # select parameter pair with highest kappa score\n", + " best_idx1, best_idx2 = np.unravel_index(\n", + " np.argmax(mean_scores, axis=None), np.shape(mean_scores)\n", + " )\n", + " best_params = {param1_name: param1[best_idx1], param2_name: param2[best_idx2]}\n", + "\n", + " # generate heatmap\n", + " param_heatmap(params, mean_scores, clf.__class__.__name__)\n", + "\n", + " return best_params\n", + "\n", + "\n", + "def param_heatmap(params, scores, clf_name):\n", + " \"\"\"\n", + " Given a dictionary of two parameter ranges, scores \n", + " for each pair of parameter values, and classifier name, \n", + " generate heatmap showing model performance scores for each \n", + " pair of parameter values.\n", + "\n", + " Parameters\n", + " ----------\n", + " params : dictionary of hyperparameter ranges\n", + " scores : ndarray \n", + " Scores for each parameter value pair\n", + " clf_name : string\n", + " Name of sklearn-compliant classifier\n", + " \"\"\"\n", + " param1_name = list(params.keys())[0]\n", + " param2_name = list(params.keys())[1]\n", + " param1 = params[param1_name]\n", + " param2 = params[param2_name]\n", + "\n", + " scores = -np.array(scores)\n", + " scores = scores.ravel().argsort().argsort().reshape(scores.shape)\n", + " plt.figure(figsize=(8, 6))\n", + " plt.subplots_adjust(left=0.2, right=0.95, bottom=0.15, top=0.95)\n", + " plt.imshow(scores, interpolation=\"nearest\", cmap=plt.cm.Blues)\n", + " plt.xlabel(param2_name)\n", + " plt.ylabel(param1_name)\n", + " plt.colorbar()\n", + " plt.xticks(np.arange(len(param2)), param2)\n", + " plt.yticks(np.arange(len(param1)), param1)\n", + " plt.title(\"Grid Search Kappa Rank \" + clf_name)\n", + " plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Building classifiers and specifying parameter ranges or values to search\n" + ] + }, + { + "cell_type": "code", + "execution_count": 13, + "metadata": {}, + "outputs": [], + "source": [ + "# get some data\n", + "X, y = fetch_openml(data_id=40979, return_X_y=True, as_frame=True)\n", + "y = pd.factorize(y)[0]\n", + "X = X.apply(lambda x: pd.factorize(x)[0])\n", + "n_samples, n_features = np.shape(X)\n", + "\n", + "# build a classifier with warm_start=True\n", + "extraTrees = ExtraTreesClassifier(warm_start=True)\n", + "\n", + "# specify parameters and ranges or values to search\n", + "extraTrees_param_dict = {\n", + " \"max_features\": [\"sqrt\", \"log2\", None],\n", + " \"n_estimators\": [10, 30, 50, 70],\n", + "}\n", + "\n", + "# build another classifier with warm_start=True\n", + "rf = RandomForestClassifier(warm_start=True)\n", + "\n", + "# specify parameters and ranges or values to search\n", + "rf_param_dict = {\n", + " \"max_features\": [\"sqrt\", \"log2\", None],\n", + " \"n_estimators\": [10, 30, 50, 70],\n", + "}" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Obtaining best parameters dictionary and refitting" + ] + }, + { + "cell_type": "code", + "execution_count": 14, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc4AAAFuCAYAAAAI+vheAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3de5wkdXnv8c93EeQO6iIRMYDIwQNRUVBU1OAtIUIw3lAS4zVyPIlR491zNJJookaTaBI1rqLo0RjwEiXGiIjiLYoCQkDxgiJyk6vcVhCXfc4fVaPNsDvbNT09NTXzee+rXttT3V31dE1PP/38fr+qX6oKSZI0nlV9ByBJ0pCYOCVJ6sDEKUlSByZOSZI6MHFKktSBiVOSpA5u13cAkqSVbbPtd6tad+NE26gbrzixqg5ZoJDmZOKUJPWq1t3I7fc+YqJt3HTm21YvUDibZOKUJPUskOH0HJo4JUn9CpD0HcXYhpPiJUlaAqw4JUn9s6lWkqQOBtRUa+KUJPVsWIODhhOpJElLgBWnJKl/NtVKkjSmMKimWhOnJKlnGVTFOZwUL0nSEmDFKUnqn021kiR1YFOtJEnLkxWnJKlnw7oAgolTktSvgc2OYuKUJPVvQBXncCKVJGkJsOKUJPXMPk5JkrpZZR+nJEnjGdi1aocTqSRJS4AVpySpf56OIknSuBwcJElSNwOqOIeT4legJP+c5NVz3F9J7rGYMXWV5EdJHtV3HEORZPf29+qX2gWS5HVJjp3i9r+b5KHt7SR5f5JrkvxXkoOTfGta+1Y/TJyLJMlTkpyaZG2Sy9vbf5xs/GtWVT23ql47z/3tm+QzSa5u/4hPT/KY+b+C6UvyjCRfHvl5+yRfSfLRJFv0GduoJEcn+UWSG0Y+IB/UQxwHJ1nfxjG6bDKWJMcmed0E+/4/I/u7KcktIz8veqJI8ofte3xtkkuT/EeSBy/Gvqtq76r6UvvjwcBvArtU1YOr6pSq2ncx4hi8rJps2dTmk/e0n73njKy7Y5KTkny//f8O44Rq4lwESV4MvBV4E/BrwM7Ac4GDgA0mhCSbTbjbfwdOavd3Z+D5wHUTbvM2plUZtW/gk4ELgCdX1c3T2M8EjquqbYHVwOeBD/cUxyVVte2s5auTbnRTv9eq+uuZ/dG8l786sv/bJIppVtBJXga8GXgtsBOwG7AGeOy09jmH3YDzq+pnk25oRbU6JJMvm3YscMisda8ATq6qvWg+b14xzoZMnFOWZAfgL4E/rqqPVNX11fhmVf1BVf28fdyxSd6R5FNJ1gIPn10VJHlp+236kiTPmmOfq4E9gHdV1c3t8pWqGq3mDkty5kjFdO+R+16R5AdJrk/y7SSPG7nvGW0V+PdJrgKObtc/J8m5I8+530hI+yX57yTXJjkuyZabOGY70SSjc4CnVtW6dv2hSb6Z5LokFyY5euQ5M02cR7XH59IkLxm5/+gkH2n3f32SM5LcZ5zXPJc2tg8Cd23jJskdknwyyRVJftre3nVkX6ckeW17HK9vWwZWb+RYPCFNc/dvjBPPyPPumOSiJL/b/rxtkvOSPC3JUcAfAC9rK8R/bx/zoyQvT/LfwNokt5vvcWmfW2laVc4DvtOu3yfJZ9O0hHwnyRNGnrNlkr9rf7eXJXn7zHslyZ3bv41r2ud+ceZY07wHn1tVH6+qn7Xv909U1cs3ENeq9n3wk3ZbpyT5nyP3HzbyPr4oyZ/Ntf/2vovSVP9HAf8MPLQ9rq9O8qgkPxp57K5J/q19b5yf5E9G7ntd+/78UJLrgaeOc6yXjSlXnFX1ReDqWasfC7yvvf0+4PfGCdXEOX0PAm4PfGKMx/4+8FfAdsCXR+9IcgjwEuDRwF7AXP2GVwHnAR9I8ntJdp61rfsC7wH+F3An4J3ACUlu3z7kB8BDgR2Av2i3c5eRTRwI/JCmcv6rJE+i+fB6GrA9cHgbw4wjaL7p7QHcG3jGHLHfETgF+CrwrKpaP3Lf2nYfOwKHAv87yew3+sNpjs9vAS/PrftXH0tTGd4R+Bfg40k2H/M1b1CaJuSnta/3p+3qVcB7aaqPXwduBP5p1lN/H3gmTWvAFjS/29nbfibwRuBRVXXO7PvnUlVXA88C3pXkzsDfA2dW1furag1Nsv+btkL83ZGnHklzbHdsvxTM67iMOBy4P3CvJNvStIK8v33dfwCsSbJ3+9g38av3yF7A7sD/be97Kc17bieaVpRXtesPohnkOM7f14xPttv/NZovZ/9v5L73As+uqu3aOL6wif3/Untcnwd8qT2ut+pmSbKq3fc3gLvS/C2/NMkjRx72OJr35g7AcR1ek+Zn56q6tL39E5rPtE0ycU7fauDKmaoJIE2Fd02SG5M8bOSxn2grw/VVddOs7RwBvLeqzqmqtbSV3oZUVdEkkB8BfwtcmuSLSfZqH3IU8M6qOrWqbqmq9wE/Bx7YPv/DVXVJG8dxwPeBB4zs4pKq+seqWldVNwJ/RPMh/I22mj6vqi4Yefw/tNu7mqYJeb85jtfdgP8BHNu+jtHXdUpVnd3G9d/Ah2j6k0b9RVWtraqzaT4Ejxy57/S26v8F8HfAlh1e82xHJLmGJik+B3jizO+4qq6qqo+21c/1NF+GZsf53qr6Xnv8jt/AMXkhzYf1wVV13hxx7NK+l0aXbdo4PkPzReFk4DE0X5Q25R+q6sI2rvkcl9n+uqp+2m7vscD32uS9rqpOBz4OPLFNKs8BXtg+/jrg9cBT2u38AtgF+PW2opyp+O4EXD7rC9ZGta/j2Lbl5yaav6P9Z45Zu599kmxXVVdX1Rmb2H8XDwK2b5u5b25/r8eMvEaAL1fVv7dx3jiPfQzX5E21q5OcNrIc1WX37edNbfKBmDgXw1U0v9Bf9ldUM2hgx/a+0d/BhXNsZ5dZ91+wsQe2+7ioqp5XVXvSVD5rab7p0/784tEPW5qEtQtA25x35sh9v0HzBWBjcd6NpjLZmJ+M3P4ZsO0cjz2Lpvr6z7Yy/qUkByb5fNvMdS1N39rsJs7Zx2iXDd3XftBexPivebbj29/hzjRVy/4jcW6d5J1JLkhyHfBFYMfcut96U8fkpcDbquqiOWKA5kvMjrOWtSP3r2lfy7FVddVGtjHqVr/beRyXuba3G3DQrPfdk4G70FRxtwfOGrnvkzSVKcAbaH6fJ7dNxy9t118F3LlNvJuUZLMkf5Pkh+3vZuZLycxrehxNlfzjthn3wE3sv4vdgF+f9fpf1r72GXN9BixjWYim2iur6oCRZc0YO75spgWl/f/ycaI1cU7fV2mquXEGKsz1bedSmgQ149fHDaCqLgTeRvOhB80f51/N+rDduqo+lGQ34F00TU53apPDOcBo7/vsOC8E9hw3njHifSvNB9VJuXXf3r8AJwB3q6odaPqTZo8KmH2MLtnQfe0H7a7AJWO+5o3FeiVNBX/0SBPmi4G9gQOrantgplWhy4lqvwW8KiN9gF21iXoNzRemP86tT13a2Hvtl+snOS4b2c+FNAMxRt9321bV84DLgJuBvUfu26H9PVNV11XVn1XV7jT9UC9P8pvAV4B1NMluHE+jqb4fQdMcOnNM0u7n1Ko6nCZhfxL4103sv4sLge/Pev3b1a2byseqeJal6Q8O2pATgKe3t5/OmE3+Js4pq6praPqG3p7kiUm2SzNAYT9gm008fdTxwDPSDK7YGnjNxh6YZnDKXyS5R7uv1TT9XV9rH/Iu4LltBZck26QZeLNdG1MBV7Tbeia/Srgb827gJUn2b7d3j/ZDd96q6m9oRiJ/dqQPbDvg6qq6KckDaPoJZ3t1W/HtS9OHONpPtH+Sx7fV/wtpvtB8jfm95tFYvwucSFM9zMR5I3BNkjsyx+9qDt+i6Rd+W5Jxk8Js/4fmdT2Lpv/w/SNV72XA3Tfx/ImOywacAOyb5PeTbN4uD0iyd1XdQvM+ekuSndr30a5Jfqvd9+8m2TNJgGuBW4D1VfVTmr+vdyQ5PMlW7XYPTfKGDcSwHc3v/Spga5pmdNp9bNXGtn3bnH89sH6u/Xd8/V8Fbk7y4jQDoTZLcq8k+2/ymZpYkg/R/A72TjOg69k0X9AfneT7NONGNvSeuQ0T5yJok8CLaD5YL2uXdwIvB/5rzG38J/AW4HM0zUufm+PhN9MMrPgszSko59B8WDyj3dZpNP1J/0QzoOW8kfu+TdMv+tU2znvRfKufK7YP03wA/QvNh83HaQbgTKSawRXvpmke2xP4Y+Av04w4/HOaLxOzfaF9PScDb277+WZ8gqZp8KfAHwKPr6pfzOc1b8CbgKPSDMR5C7AVcCVNYv50x20BUFVnAYfRDPD5nY08bJfc9jzOJ7Qfxi8CntYmpTfSJMGZ4fbH0PTlXZPk4xvZ/0Icl9HtXQv8Ns1o0UtpmqtfT9NEC02lfgHwdZrk9BmaQTzQVPCfA25oY3hrtedOVtUbaf6WjqZJiBcC/5vmfTjbe2laIS6h+XIy++/v6cBME/uz+dXI1o3uv8PrX0dT7T6AZvzBlTSfA9t32c6yNDM7ynRH1R5ZVXepqs2rateqOqYdj/DIqtqrqh5VzTiMTYdbtXJbBrR8JNkdOB/YvEYGYo3cfzRwj6paWUP8pQFYteNudfuH3ubsoU5u+uSfnF5VByxQSHNaOSfYSpKWLq9VK0nS8mTFqWWhqn7EHKM9q+roRQtGUndOKyZJUgc21UqStDwtq4pz1Zbb1aptd+o7jGVnz19ztPy03Ly+66mAGscdt1oys9AtKxdc8COuvPLKhS8NE5tq+7Jq253Y8bC/7juMZeedr3x03yEsWz++fuLZp7QBR+x3t00/SJ0ddOAUz/YYUFPtskqckqRhyoAS53BqY0mSlgArTklSr8KwKk4TpySpX6HbnDs9M3FKknqWQVWc9nFKktSBFackqXdDqjhNnJKk3pk4JUnqYEiJ0z5OSZI6sOKUJPXL01EkSRpfBnY6iolTktS7ISVO+zglSerAilOS1LshVZwmTklS70yckiSNa2Cjau3jlCSpAytOSVLvbKqVJGlMnscpSVJHQ0qc9nFKktSBFackqX/DKThNnJKknsWmWkmSli0rTklS74ZUcZo4JUm9M3FKkjSmoZ3HaR+nJEkdWHFKkvo3nILTxClJ6tnATkcxcUqSejekxGkfpyRJHVhxSpJ6N6SK08QpSerfcPKmiVOS1L8hVZz2cUqS1IEVpySpV8mwrhy05BNnkv2AXarqU33HIkmajiElziXdVJvkdsB+wGP6jkWSND0zVed8l8W0KBVnkm2A44Fdgc2A1wLXAm8BfgZ8Gbh7VR2W5GhgT+DuwI+Bg4CtkjwEeH1VHbcYMUuStCGL1VR7CHBJVR0KkGQH4BzgEcB5wOxkuA/wkKq6MckzgAOq6nmLFKskabENp6V20ZpqzwYeneSNSR4K7AGcX1Xfr6oCPjDr8SdU1Y3jbDjJUUlOS3Ja3XT9AoctSVoMQ2qqXZTEWVXfA+5Hk0BfBxy+iaes7bDtNVV1QFUdkC23myBKSVIvMqzEuVh9nLsAV1fVB5JcAzwP2D3JnlX1A+DIOZ5+PWBGlCQtCYvVVHsv4OtJzgReA7wKOAr4jyRnAJfP8dzPA/skOTPJk6cfqiRpMQVIJlvG2k/yZ0m+leScJB9KsuV84l2UirOqTgRO3MBd9wRIcjDwkvaxR8967tXA/acboSSpP9Nvbk1yV+D5wD7twNPjgacAx3bd1pK/AIIkaflbpG7K29Gc3vgLYGvgkvlsZElcAKGqTqmqw/qOQ5K0PFXVxcCbaa4PcClwbVV9Zj7bWhKJU5K0si3AqNrVM6cmtstRs7Z/B+CxNKdD7gJsk+Sp84nVplpJUr86DPCZw5VVdcAc9z+K5voBVwAk+RjwYG57HYFNsuKUJK0EPwYemGTrNCXqI4Fz57MhK05JUq8CrFo13dFBVXVqko8AZwDrgG8Ca+azLROnJKl3izGqtqpeQ3MtgYmYOCVJvXM+TkmSlikrTklSvxZmVO2iMXFKknrVXKt2OJnTxClJ6tniTw02Cfs4JUnqwIpTktS7ARWcJk5JUv+G1FRr4pQk9Wtgo2rt45QkqQMrTklSrzwdRZKkjgaUN02ckqT+DanitI9TkqQOrDglSb0bUMFp4pQk9SzDaqo1cUqSetWMqu07ivHZxylJUgdWnJKkng1rdhQTpySpdwPKmyZOSVL/hlRx2scpSVIHVpySpH4NbHYUE6ckqVdDu8i7TbWSJHVgxSlJ6t2QKk4TpySpdwPKmyZOSVL/hlRx2scpSVIHVpySpH55Okp/1v/sBtae9eW+w1h2Pvej+/UdwrJ1xL679B3CsnTVDTf3HcKytG59TWW78Vq1kiR1M6C8aR+nJEldWHFKknq3akAlp4lTktS7AeVNE6ckqV+J53FKkrRsWXFKknq3ajgFp4lTktS/ITXVmjglSb0bUN60j1OSpC6sOCVJvQrNZfeGYuyKM8mTkmzX3n5Vko8l8SKmkqSJrcpky6LG2uGxr66q65M8BHgUcAzwjumEJUlaMdJc5H2SZTF1SZy3tP8fCqypqv8Atlj4kCRJWrq69HFenOSdwKOBNya5PQ4ukiQtgCGNqu2SOI8ADgHeXFXXJLkL8NLphCVJWinCsC7yPnbFWFU/Ay4HHtKuWgd8fxpBSZJWluZ6tfNfFlOXUbWvAV4OvLJdtTnwgWkEJUnSUtWlqfZxwH2BMwCq6pKZ01MkSZrEcr3k3s1VVUkKIMk2U4pJkrSC9NHcOokuo2KPb0fV7pjkOcBngXdNJyxJkpamsSvOqnpzkkcD1wF7A39eVSdNLTJJ0ooxpFG1YyXOJJsBn62qhwMmS0nSghpO2hyzqbaqbgHWJ9lhyvFIklagxbjkXpIdk3wkyXeSnJvkQfOJtcvgoBuAs5OcBKydWVlVz5/PjiVJWmRvBT5dVU9MsgWw9Xw20iVxfqxdJElaMM2Vg6a8j6bF9GHAMwCq6mbg5vlsq8vgoPfNZweSJM1pcWY42QO4AnhvkvsApwMvqKq1cz/ttrpcOej8JD+cvXTdoSRJsy3AJfdWJzltZDlq1i5uB9wPeEdV3Zemy/EV84m1S1PtASO3twSeBNxxPjuVJGmBXVlVB8xx/0XARVV1avvzR5hn4uxykferRpaLq+otNHNzSpI0kWmPqq2qnwAXJtm7XfVI4NvziXXsijPJ/UZ+XEVTgXapWCVJuo3FGBzU+lPgg+2I2h8Cz5zPRrokvr8dub0OOJ9mjk5JkiayGBd5r6ozuXW347x0SZzPrqpbDQZKssekAUiSNCRdLvL+kTHXSZLUSSZcFtMmK84k9wT2BXZI8viRu7anGV0rSdK8JcvvIu97A4cBOwK/O7L+euA50whKkrSyDChvbjpxVtUngE8keVBVfXURYpIkacnqMjjom0n+hKbZ9pdNtFX1rAWPSpK0oizGqNqF0mVw0P8Dfg34beALwK40zbWSJE1kAS65t2i6JM57VNWrgbXtBd8PBQ6c6wlJbphvYEk+mOS7Sc5J8p4km893W5KkpSuEVZlsWUxdEucv2v+vSfIbwA7AnRc+pF/6IHBP4F7AVsAfTXFfkiSNpUviXJPkDsCrgRNorvH3N+M8MY03tdXj2Ume3K5fleTt7WzcJyX5VJInAlTVp6oFfJ2maViStNxM2Ey72E21XebjfHd78wvA3Tvu5/HAfsB9gNXAN5J8ETgI2B3Yh6Z6PRd4z+gT2ybaPwRe0HGfkqSBWJaDg5LsnOSYJP/Z/rxPkmeP+fSHAB+qqluq6jKa5Hv/dv2Hq2p9e+X6z2/guW8HvlhVX9pIXEfNzL9W624c9+VIkpaQVRMuix3ruI4FTgR2aX/+HvDChQ5oVJLXADsBL9rYY6pqTVUdUFUH5HZbTTMcSZI6Jc7VVXU8sB6gqtYBt4z53C8BT06yWZKdgIfR9Ft+BXhC29e5M3DwzBOS/BHNqS9HVtX6DnFKkgYkTH8+zoXU5QIIa5PcCSiAJA8Erh3zuf8GPAg4q33+y6rqJ0k+yq8mE70QOGNkm/8MXAB8tT0oH6uqv+wQryRpIBZpPs4F0SVxvohmNO2eSb5C04T6xLmeUFXbtv8X8NJ2Gb1/fZKXVNUNbVL+OnB2e5+TZEuSlpxxZkd5UlV9GPgp8Js0F30P8N2q+sWcTx7PJ5PsCGwBvLYdJCRJWkGWW8X5SuDDwEer6n7AtxYygKo6eCG3J0kaluZczOFkznES51VJPgPskeSE2XdW1eELH5YkaSVZbhXnocD9aC7y/rfTDUeSpKVtnPk4bwa+luTBVXXFxh6X5B+r6k8XNDpJ0oowoJbaTpfc22jSbB00YSySpBUosOgznEzCUz4kSb1b7MvmTWJIsUqS1LuFrDiHU2dLkpaUAbXUjp84k2xZVTfNWre6qq5sf3zrgkYmSVoRkgyqj7NLU+032uvTApDkCcB/zfxcVccuYFySpBVkWU5kDfw+8J4kp9BMLXYn4BHTCEqSpKWqy+koZyf5K5oLIVwPPKyqLppaZJKkFWO5XTkIgCTHAHsC9wb+B83F2f+xqt42reAkScvfcj6P82zgj9opws5PciDwd9MJS5K0kgwob3Zqqn3LrJ+vBZ694BFJkrSEdWmq3Qt4PbAPsOXM+qq6+xTikiStFBlWH2eX01HeC7wDWAc8HHg/8IFpBCVJWlky4b/F1CVxblVVJwOpqguq6miaKcckSZq3ZnDQZMti6jI46OdJVgHfT/I84GJg2+mEJUnS0tSl4nwBsDXwfGB/4KnA06YRlCRpZVmuFWfRXPxgN2Dzdt27aM7rlCRp3jKg81G6JM4PAi+lOZ9z/XTCkSStNDN9nEPRJXFeUVUnTC0SSZIGoEvifE2SdwMnAz+fWVlVH1vwqCRJK0cPM5xMokvifCZwT5r+zZmm2gJMnJKkiSzXa9Xev6r2nlokkiQNQJfTUf4ryT5Ti0SStCIt5wsgPBA4M8n5NH2cAaqqPB1FkjSRAbXUdkqch0wtCknSChZWLfL1ZifRZVqxC6YZiCRJQ9Cl4pQkacGF5dtUK0nSwhvYfJwmTklS74Z0HmeX01EkSVrxrDglSb2yj1OSpI6G1FRr4pQk9W5AeXOZJc7Nbw8779l3FMvOcZ/7Yd8hLFvnX3Fj3yEsSxdeeUPfISxLP7xybd8hTCzJZsBpwMVVddh8trG8EqckaXDCoo5UfQFwLrD9fDfgqFpJUr8CSSZaxtpNsitwKPDuScK14pQk9W6RujjfArwM2G6SjVhxSpKWg9VJThtZjhq9M8lhwOVVdfqkO7LilCT1qpmPc+Ka88qqOmCO+w8CDk/yGGBLYPskH6iqp3bdkRWnJKl3mXDZlKp6ZVXtWlW7A08BPjefpAlWnJKkJcDzOCVJWqKq6hTglPk+38QpSerZ+KeULAUmTklSrxb5AggTM3FKkno3pIpzSElekqTeWXFKkno3nHrTxClJ6ltsqpUkadmy4pQk9cpRtZIkdTSkploTpySpd8NJm8OqjiVJ6p0VpySpdwNqqTVxSpL61QwOGk7mNHFKkno3pIrTPk5Jkjqw4pQk9SzEplpJksY3pKZaE6ckqVdDGxxkH6ckSR1YcUqS+hWbaiVJ6sTEKUlSB0MaVWsfpyRJHVhxSpJ6FWDVcApOE6ckqX9Daqo1cUqSejekwUH2cUqS1IEVpySpdzbVSpI0JgcHSZLUybBmR7GPU5KkDqw4JUn98lq1kiR1M6C8aVOtJEldDL7iTHIUcBQAW96h32AkSZ01o2qHU3MOvuKsqjVVdUBVHZAttuk7HEnSPGTCZTENvuKUJC0Dwyk4h19xSpK0mKw4JUm9G9IFEEyckqTeDWhskIlTktS/AeVN+zglSerCilOS1L8BlZwmTklSr5pzMYeTOU2ckqR+Dewi7/ZxSpLUgRWnJKl3Ayo4TZySpCVgQJnTxClJ6lkGNTjIPk5Jkjqw4pQk9W5Io2pNnJKkXvUxp+YkTJySpP4NKHPaxylJUgcmTklS7zLhv01uP7lbks8n+XaSbyV5wXxjtalWktS7RRgctA54cVWdkWQ74PQkJ1XVt7tuyIpTktS7TLhsSlVdWlVntLevB84F7jqfWE2ckqQVJcnuwH2BU+fzfJtqJUn9WpjzUVYnOW3k5zVVteY2u0q2BT4KvLCqrpvPjkyckqTeLcAl966sqgPm3EeyOU3S/GBVfWy+O7KpVpK07CUJcAxwblX93STbMnFKknoVmlG1kyxjOAj4Q+ARSc5sl8fMJ16baiVJvZv22ShV9eWF2o2JU5LUPy+5J0nS8mTFKUnq3ZAmsjZxSpJ653yckiR1MKC8aR+nJEldWHFKkvo3oJLTxClJ6lVzqdrhZE4TpySpX+Nf/WdJsI9TkqQOrDglSb0bUMFp4pQkLQEDypwmTklSzzKowUH2cUqS1IEVpySpd0MaVbusEmddd9GVN534ogv6jmNMq4Er+w5iHN87se8IOhnMcQX4Xt8BdDOoYzsgQzquu01jo2FQXZzLLHFW7dR3DONKclpVHdB3HMuNx3V6PLbT4XFtDShz2scpSVIHy6rilCQN05BG1Zo4+7Om7wCWKY/r9Hhsp8PjioODNIaq8o9lCjyu0+OxnQ6Pa2NAedM+TkmSujBxLoIk70lyeZJzRtbdMclJSb7f/n+HPmMcoiRbJvl6krOSfCvJX7Tr90hyapLzkhyXZIu+Yx2aJD9KcnaSM5Oc1q7zPTuhJHu3x3RmuS7JC1f8sW1nR5lkWUwmzsVxLHDIrHWvAE6uqr2Ak9uf1c3PgUdU1X2A/YBDkjwQeCPw91V1D+CnwLN7jHHIHl5V+42cKuF7dkJV9d32mO4H7A/8DPg3PLb86mzO+S6Lx8S5CKrqi8DVs1Y/Fnhfe/t9wO8talDLQDVuaH/cvF0KeATwkXa9x3bh+J5dWI8EflBVF+CxHRQTZ392rqpL29s/AXbuM5ihSrJZkjOBy4GTgB8A11TVuvYhFwF37Su+ASvgM0lOT3JUu8737MJ6CvCh9vaKPrZhWE21jqpdAqqqklTfcQxRVd0C7JdkR5omr3v2HNJy8ZCqujjJnYGTknxn9E7fs5Np+90PB145+76VemwdVatxXJbkLgDt/5f3HKvQBH4AAAUQSURBVM+gVdU1wOeBBwE7Jpn5UrgrcHFvgQ1UVV3c/n85zReSB+B7diH9DnBGVV3W/rzij+2QKk4TZ39OAJ7e3n468IkeYxmkJDu1lSZJtgIeDZxLk0Cf2D7MY9tRkm2SbDdzG/gt4Bx8zy6kI/lVMy14bAfFptpFkORDwMHA6iQXAa8B3gAcn+TZwAXAEf1FOFh3Ad6XZDOaL4HHV9Unk3wb+NckrwO+CRzTZ5ADtDPwb2m+xt8O+Jeq+nSSb+B7dmLtl5FHA/9rZPWK/zwY0iX3UrXimtIlSUvIfe67f534ha9NtI277LDF6Ys1y4wVpySpd8OpN+3jlCSpEytOSVKv+hgZOwkTpySpd0MaHGTilCT1bzh50z5OSZK6MHFKU5BkvySPGfn58CQLMuNFOw3V1guxLWmpGM7cKCZOaVr2A36ZOKvqhKp6wwJt+4VAp8TZXiRCWrK85J40EEl2T3Jukne1k2F/pr1834Yeu2eST7czhnwpyT3b9U9Kck47ofYX2wt4/yXw5Hay4icneUaSf2off2ySdyT5WpIfJjm4nez83CTHjuzvHUlOmzVJ9/OBXYDPJ/l8u+7IdtLpc5K8ceT5NyT52yRnAQ9K8oYk307y30nePJ0jKs1HJv63mEycEuwFvK2q9gWuAZ6wkcetAf60qvYHXgK8vV3/58BvtxNqH15VN7frjmsnLT5uA9u6A80F6f+M5jqlfw/sC9wryX7tY/5veyWUewO/meTeVfUPwCU0k0w/PMkuNBN3P4Kmyr1/kpm5HLcBTm3jOhd4HLBvVd0beF3XgySpYeKU4PyqOrO9fTqw++wHJNkWeDDw4Xb+z3fSXCsX4CvAsUmeA4zbJPrv1Vzv8mzgsqo6u6rWA98a2f8RSc6gud7uvsA+G9jO/YFTquqKdg7SDwIPa++7Bfhoe/ta4CbgmCSPB342ZpzS1DkfpzQ8Px+5fQuwoabaVTQTZO83+46qem6SA4FDgdOT7N9hn+tn7X89cLske9BUtfevqp+2TbhbjrHdUTe185VSVeuSPAB4JM3MMc+jqVIldWTFKY2hqq4Dzk/yJIA07tPe3rOqTq2qPweuAO4GXA9sN8EutwfWAtcm2Zlm/sYZo9v+Ok0z7up2ANCRwBdmb6ytmHeoqk/RNA/fZ4LYpAVnxSktT38AvCPJq4DNgX8FzgLelGQvmhank9t1PwZe0Tbrvr7rjqrqrCTfBL4DXEjTHDxjDfDpJJe0/ZyvoJmDNMB/VNWG5nLcDvhEki3bx72oa0ySGk4rJknq1X3vd0Cd8pWvT7SNHbfezGnFJEkrhBd5l4YtyduAg2atfmtVvbePeKTlro+r/0zCxCnNUlV/0ncMkpYuE6ckqX8DKjlNnJKk3g1pPk7P45QkqQMrTklS7xxVK0lSBwPKmzbVSpKWgEWYyTrJIUm+m+S8SSaWN3FKkpa99lrOb6O57vM+wJFJNjTj0CaZOCVJvVuEiawfAJxXVT9s58z9V+Cx84nVPk5JUq9m5uOcsrvSTJgw4yLgwPlsyMQpSerVGWecfuJWm2f1hJvZMslpIz+vqao1E25zg0yckqReVdUhi7Cbi2nmyp2xa7uuM/s4JUkrwTeAvZLskWQL4CnACfPZkBWnJGnZq6p1SZ4HnAhsBrynqr41n205kbUkSR3YVCtJUgcmTkmSOjBxSpLUgYlTkqQOTJySJHVg4pQkqQMTpyRJHZg4JUnq4P8DbJ24lo1EFF8AAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAc4AAAFuCAYAAAAI+vheAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZgsZXn38e/vIIogmx5EERVEgi9EOQqKexQ1IUIwcUFJjMGN+CZGjbt5VXCLu9EkxICKGFFk0USSqIhENG4oIAiICi7IJqussgjnfv+oGm2Gc+Z0Tc9MTc18P+eq63RXdVfdXd3Td9/P81RVqgpJkjSeFX0HIEnSkJg4JUnqwMQpSVIHJk5JkjowcUqS1IGJU5KkDu7QdwCSpOVtvU3uW3XLDROto2647Liq2mOOQpqRiVOS1Ku65QbutMM+E63jxtMOWjlH4ayTiVOS1LNAhtNzaOKUJPUrQNJ3FGMbToqXJGkRsOKUJPXPplpJkjoYUFOtiVOS1LNhDQ4aTqSSJC0CVpySpP7ZVCtJ0pjCoJpqTZySpJ5lUBXncFK8JEmLgBWnJKl/NtVKktSBTbWSJC1NVpySpJ4N6wQIJk5JUr8GdnUUE6ckqX8DqjiHE6kkSYuAFackqWf2cUqS1M0K+zglSRrPwM5VO5xIJUlaBKw4JUn983AUSZLG5eAgSZK6GVDFOZwUv0Qk+dckb5hheSW5/0LG1FWSnyV5Yt9xDEWSbdr3ddH9UE1yWJK39h3HYpFkvyRfm8f1fz7JX4zcf2uSy5P8Isl9klyXZL352r7mholzAkmeleSkJNcnubS9/VfJ2n86VdWLquots9zeTkm+mOTKJFclOSXJk2f/Cubf9C+iJJsk+XqSTye5Y5+xjUpyYJJft19cVyX5RpJH9BDH45KsbuO4NskPkzx3oeOY1LTXMTX95wLHcGCSw9cw/w+SfLXdv5cl+UqSvRcipqr6w6r6WBvHfYBXADtW1T2q6udVdZequnUhYll0smKyaV2rTw5tv6fPHJl31yTHJzmn/X/zcUI1cc5SklcAHwDeDdwD2BJ4EfAoYI0JYQ5+Sf4ncHy7vbsDLwGumXCdtzNflVH7oTwBOA94ZlXdPB/bmcCRVXUXYCXwZeDonuK4qI1jE+BvgQ8l2aGnWCZxUZsIpqY/6rqCuf4sJnk6zfv6b8DWNH+3bwQ6xzYH7gNcUVWXTrqixdia0Uky+bRuhwF7TJv3WuCEqtqe5rvpteOsyMQ5C0k2Bd4M/FVVHVNV11bju1X1Z1V1U/u4w5J8MMnnklwPPH5601iSVyW5OMlFSZ43wzZXAtsCH6qqm9vp61U1Ws3tleS0kYrpQSPLXpvkx+2v7O8n+ZORZfu1VeA/JLkCOLCd/8IkZ4885yEjIa1K8r0kVyc5MskG69hnW9AkozOBZ1fVLe38PZN8N8k1Sc5PcuDIc6aaOPdv98/FSV45svzAJMe02782yalJdh7nNc+kje0TwL3auEmyeZL/aiuUX7a3tx7Z1olJ3tLux2vbloGVa9kXT0vT3P2764ijqupzwJXA6Hv5gXZfXZOm1eEx0/bJUUn+rY3jrCS7jix/cLufrk1yJHCb9619z89N06pxbJKtRpZVmhaVc9rnvyXJdu1n7Zp2u+tsRUhypyTvb9/Ti9rbd2qXPS7JBUlek+QXwEfb+TN9tl+T5ML8tkJ/QpI9gL8Dnpmm2j09SYD3AW+pqg9X1dVVtbqqvlJVL1xLrDPt64clOblddkmS97XzN0hyeJIr2ni/k2TLdtmJSV6QpqvjeGCrNr7DMq1JP8mmST7Sfu4vTNOsu167bI1/s4M2zxVnVX2V5m9p1FOAj7W3Pwb88Tihmjhn5xHAnYDPjvHYPwXeBmwM3KbvpP3jfiXwJGB7YKZ+wyuAc4HDk/zx1B/iyLoeDBwK/CVwN+Bg4NipLyTgx8BjgE2BN7XruefIKnYDfkLzC/xtSZ5B88f4HJrKZ+82hin70Px625bmS32/GWK/K3Ai8E3geVW1emTZ9e02NgP2BP5vkukf3sfT7J/fB16T2/avPoWmgrgr8EngP5KsP+ZrXqP2y/857ev9ZTt7Bc2X+H1pKoUbgH+e9tQ/BZ5L0xpwR5r3dvq6nwu8E3hiVZ05ffm0x65I04S4kua9n/IdYNXIaz46t/3hsjfwKZp9euxUnO3r+g/g4+1zjwaeNrK93YG307y396RpGfjUtLD+ANgFeDjwauAQ4NnAvYHfBfad6TW1/l/7/FXAzsDDgNePLL9HG999gf1n+mynqcRfDDy0qjZu4/tZVX0B+HvaVoSq2hnYoY3zmDFinDLTvv4A8IGq2gTYDjiqnf8XNJ+5e7fxvojm8/IbVfUl4A/5bVW+3xq2fRhwC3B/4ME0n/8XjCy/zd9sh9ek39qyqi5ub/+CZl+uk4lzdlYCl09VTQDtr+CrktyQ5LEjj/1sWxmurqobp61nH+CjVXVmVV3PDL8aq6poEsjPgPcCF6fpp9m+fcj+wMFVdVJV3dr2o9xE8wVFVR1dVRe1cRwJnEPzhTXloqr6p6q6papuoPkDfVdVfaetfM6tqvNGHv+P7fqupGlCXjXD/ro38DvAYe3rGH1dJ1bVGW1c3wOOAH5v2vPfVFXXV9UZNMlr9Mv5lLbq/zVNNbFBh9c83T5JrqL5knsh8PSp97iqrqiqT1fVr6rqWpovqulxfrSqftTuv6PWsE9eBrwKeFxVncvabTUSx78DL6+q704trKrD23huqar30vyIG23K/VpVfa7tK/s4TXKi3S/rA++vql9X1TE0iWHKnwGHVtWpbavJ64BHJNlm5DHvqqprquosmtaDL1bVT6rqauDzNF/wt3kdI9M+I9t5c1VdWlWX0fyo+fOR560GDqiqm9p9OdNn+9b29e+YZP2q+llV/Xgt+/Vu7f8Xr2X57axjX/8auH+SlVV1XVV9a2T+3YD7t/GeUlWdulTaH8ZPBl7WfvYvBf4BeNbIw6b/zQ7b5E21K9sWgKlp/y6bb7+bap0PxMQ5W1fQvEm/6VeoqkdW1WbtstH9ev4M69lq2vLz1vbAdhsXVNWLq2o7ml/j19P01dDef8XoFxVNwtoKIMlzRpq6rqKpDkabEqfHeW+aim1tfjFy+1fAXWZ47Ok01dfn2+rhN5LsluTLaZpAr6b5dT69iXP6PtpqTcvaSvYCxn/N0x3Vvodb0iSFXUbi3DDJwUnOS3IN8FVgs9y233pd++RVwEFVdcEMMUDzhbgZTaX/j8DuowuTvDJNE/rV7evadNrrmh7HBu1ndSvgwmk/XkY/c1uN3q+q62g+z/caecwlI7dvWMP90dd8UVVtNjJNVWS32Q63f08vm/Yjc62f7fYHyMtofnRemuRTGWlenmaqxWSdrQ5T1rGvn0/zg/AHbXPsXu38jwPHAZ9K0xT9rpFWkHHdl+ZHzsUjr/lgmtaMKTN9twxM5qKp9vKq2nVkOmSMDV8y1QrV/j9Wf7OJc3a+SfOL9yljPHamXzAX03wBTLnPuAFU1fnAQTTJAJo/ordN+6LasKqOSHJf4EM0TVp3a7+UzwRGe9Snx3k+TfPTnKiqDwDvAI7Pbfv2PknTnHjvqtoU+NdpccHt99FFa1qWZAXNgI+LxnzNa4v1cpoq58CRpt1X0FQau1XTNDfVqtDl4LPfB16f5GnrfGQTx03Aa4AHTjVft31sr6Zprdi8fV1XjxnHxTT9tqOPHf3MXUTzhU27rY1oKqcLx4m3g9tsh9u/p2v6LK7xsw1QVZ+sqke36yyapvA1reeH7brG2v/r2tdVdU5V7UuTzN4JHJNko7aaf1NV7Qg8EtiLpum/i/NpvmNWjrzmTapqp5HHjFUdDcb8Dw5ak2NpmtZp/x+n+83EORtVdRVN89K/JHl6ko3b/qhVwEYdVnUUsF+SHZNsCBywtgemGZzypiT3b7e1EngeMNU89CHgRW0FlyQbpRl4s3EbUwGXtet6Lr9NuGvzYeCVSXZp13f/NhnNWlW9i6Zf6Ev57SjRjYErq+rGJA+j6Sec7g1txbcTTR/ikSPLdkny1LaiehnNl823mN1rHo31hzRVw6tH4rwBuCrJXZnhvZrBWTT9wgdlzMMfqhl5/F6akZ9TcdxC87rukOSNNJXpOL7ZPvclSdZP8lRu23R9BPDcJKvS9I3/PXBSVf1szPWP6wiaHxBbtJ/jNwK3O2xkxFo/20l2SLJ7G++NNO/RVB/6JcA27Q+qqaa4l9N8np6b5tCoFUkenWRN1cmM+zrJs5Ns0bZ0XNXOXp3k8Uke2LZGXEPTdLuaDqrpd/si8N6ROLdLMr17QGNKcgTN38AOaQagPZ/mx/yTkpxDM8bkHeOsy8Q5S20SeDnNF+sl7XQwTYXwjTHX8Xng/cD/0Az++J8ZHn4zsA3wJZo/xjNpksR+7bpOpumX+2eaAS3njiz7Ps2X7zfbOB8IfH0dsR1N04/3SeBamkEldx3nda1jvW+hSconJNkO+CvgzUmupfkCPWoNT/tK+3pOAN5TVV8cWfZZ4Jk0r/nPgae2v/g7v+Y1eDfN4JS707xPdwYup0nMX+i4LgCq6nSaCuRDSf5wzKcdCtwnyR/RJPMvAD+iaeK8kTGb7Nok/FSaz8WVNPvtMyPLvwS8Afg0TXW6HbftU5srbwVOBr4HnAGc2s5bW9xr/WzT9Dm+g+Z9+QVN9fe6dtnU4URXJDm1XdcxNK/7eTRV7iXtttdUaaxrX+8BnJXkOpofhM9q+xrvQTMA6RrgbJrP78dn3CNr9hyaQWbfb1/3MXRoZh6UqaujzO+o2n2r6p5VtX5VbV1VH2n7r59QVdtX1ROrGbOx7nBv290hLR5pBqX8FFi/RgZijSw/kGYAxrMXNjJJc2nFZvetOz3mNROt48b/+utTqmrXdT9ycsM+aFaStDR4rlpJkpYmK04tWu2glJnO+3vgggUjaX55WTFJkjqwqVaSpKVpSVWcG2y8eW20xdpOGqLZ2nrTGc/frglcdv1NfYewJN1zYz+z8+G8837G5ZdfPvelYWJTbV822mIrnvzmT/YdxpLzjr3+T98hLFkHn/SzvkNYkv7uCb/TdwhL0qN2m8ejPQbUVLukEqckaZgyoMQ5nNpYkqRFwIpTktSrMKyK08QpSepX6HadoZ6ZOCVJPcugKk77OCVJ6sCKU5LUuyFVnCZOSVLvTJySJHUwpMRpH6ckSR1YcUqS+uXhKJIkjS8DOxzFxClJ6t2QEqd9nJIkdWDFKUnq3ZAqThOnJKl3Jk5JksY1sFG19nFKktSBFackqXc21UqSNCaP45QkqaMhJU77OCVJ6sCKU5LUv+EUnCZOSVLPYlOtJElLlhWnJKl3Q6o4TZySpN6ZOCVJGtPQjuO0j1OSpA6sOCVJ/RtOwWnilCT1bGCHo5g4JUm9G1LitI9TkqQOrDglSb0bUsVp4pQk9W84edPEKUnq35AqTvs4JUnqwIpTktSrZFhnDlr0iTPJKmCrqvpc37FIkubHkBLnom6qTXIHYBXw5L5jkSTNn6mqc7bTQlqQijPJRsBRwNbAesBbgKuB9wO/Ar4G3K+q9kpyILAdcD/g58CjgDsneTTw9qo6ciFiliRpTRaqqXYP4KKq2hMgyabAmcDuwLnA9GS4I/DoqrohyX7ArlX14gWKVZK00IbTUrtgTbVnAE9K8s4kjwG2BX5aVedUVQGHT3v8sVV1wzgrTrJ/kpOTnHzTNb+c47AlSQthSE21C5I4q+pHwENoEuhbgb3X8ZTrO6z7kKratap2vdMmm08QpSSpFxlW4lyoPs6tgCur6vAkVwEvBrZJsl1V/RjYd4anXwtsvBBxSpK0LgvVVPtA4NtJTgMOAF4P7A/8d5JTgUtneO6XgR2TnJbkmfMfqiRpIQVIJpvG2k7yt0nOSnJmkiOSbDCbeBek4qyq44Dj1rDoAQBJHge8sn3sgdOeeyXw0PmNUJLUn/lvbk1yL+AlwI7twNOjgGcBh3Vd16I/AYIkaelboG7KO9Ac3vhrYEPgotmsZFGcAKGqTqyqvfqOQ5K0NFXVhcB7aM4PcDFwdVV9cTbrWhSJU5K0vM3BqNqVU4cmttP+09a/OfAUmsMhtwI2SvLs2cRqU60kqV8dBvjM4PKq2nWG5U+kOX/AZQBJPgM8ktufR2CdrDglScvBz4GHJ9kwTYn6BODs2azIilOS1KsAK1bM7+igqjopyTHAqcAtwHeBQ2azLhOnJKl3CzGqtqoOoDmXwERMnJKk3nk9TkmSligrTklSv+ZmVO2CMXFKknrVnKt2OJnTxClJ6tnCXxpsEvZxSpLUgRWnJKl3Ayo4TZySpP4NqanWxClJ6tfARtXaxylJUgdWnJKkXnk4iiRJHQ0ob5o4JUn9G1LFaR+nJEkdWHFKkno3oILTxClJ6lmG1VRr4pQk9aoZVdt3FOOzj1OSpA6sOCVJPRvW1VFMnJKk3g0ob5o4JUn9G1LFaR+nJEkdWHFKkvo1sKujmDglSb0a2knebaqVJKkDK05JUu+GVHGaOCVJvRtQ3jRxSpL6N6SK0z5OSZI6sOKUJPXLw1H6c/U1N/D5L53ddxhLzvMfunXfISxZ++y0Vd8hLElHnXZ+3yEsSVfecPO8rDeeq1aSpG4GlDft45QkqQsrTklS71YMqOQ0cUqSejegvGnilCT1K/E4TkmSliwrTklS71YMp+A0cUqS+jekploTpySpdwPKm/ZxSpLUhRWnJKlXoTnt3lCMXXEmeUaSjdvbr0/ymSQPmb/QJEnLxYpMNi1orB0e+4aqujbJo4EnAh8BPjg/YUmSlo00J3mfZFpIXRLnre3/ewKHVNV/A3ec+5AkSVq8uvRxXpjkYOBJwDuT3AkHF0mS5sCQRtV2SZz7AHsA76mqq5LcE3jV/IQlSVouwrBO8j52xVhVvwIuBR7dzroFOGc+gpIkLS/N+WpnPy2kLqNqDwBeA7yunbU+cPh8BCVJ0mLVpan2T4AHA6cCVNVFU4enSJI0iaV6yr2bq6qSFECSjeYpJknSMtJHc+skuoyKPaodVbtZkhcCXwI+ND9hSZK0OI1dcVbVe5I8CbgG2AF4Y1UdP2+RSZKWjSGNqh0rcSZZD/hSVT0eMFlKkubUcNLmmE21VXUrsDrJpvMcjyRpGVqIU+4l2SzJMUl+kOTsJI+YTaxdBgddB5yR5Hjg+qmZVfWS2WxYkqQF9gHgC1X19CR3BDaczUq6JM7PtJMkSXOmOXPQPG+jaTF9LLAfQFXdDNw8m3V1GRz0sdlsQJKkGS3MFU62BS4DPppkZ+AU4KVVdf3MT7u9LmcO+mmSn0yfum5QkqTp5uCUeyuTnDwy7T9tE3cAHgJ8sKoeTNPl+NrZxNqlqXbXkdsbAM8A7jqbjUqSNMcur6pdZ1h+AXBBVZ3U3j+GWSbOLid5v2JkurCq3k9zbU5JkiYy36Nqq+oXwPlJdmhnPQH4/mxiHbviTPKQkbsraCrQLhWrJEm3sxCDg1p/A3yiHVH7E+C5s1lJl8T33pHbtwA/pblGpyRJE1mIk7xX1WnctttxVrokzudX1W0GAyXZdtIAJEkaki4neT9mzHmSJHWSCaeFtM6KM8kDgJ2ATZM8dWTRJjSjayVJmrVk6Z3kfQdgL2Az4I9G5l8LvHA+gpIkLS8DypvrTpxV9Vngs0keUVXfXICYJElatLoMDvpukr+mabb9TRNtVT1vzqOSJC0rCzGqdq50GRz0ceAewB8AXwG2pmmulSRpInNwyr0F0yVx3r+q3gBc357wfU9gt5mekOS62QaW5BNJfpjkzCSHJll/tuuSJC1eIazIZNNC6pI4f93+f1WS3wU2Be4+9yH9xieABwAPBO4MvGAetyVJ0li6JM5DkmwOvAE4luYcf+8a54lpvLutHs9I8sx2/ook/9Jejfv4JJ9L8nSAqvpctYBv0zQNS5KWmgmbaRe6qbbL9Tg/3N78CnC/jtt5KrAK2BlYCXwnyVeBRwHbADvSVK9nA4eOPrFtov1z4KUdtylJGoglOTgoyZZJPpLk8+39HZM8f8ynPxo4oqpurapLaJLvQ9v5R1fV6vbM9V9ew3P/BfhqVf3vWuLaf+r6a3WjY5UkaYhWTDgtdKzjOgw4Dtiqvf8j4GVzHdCoJAcAWwAvX9tjquqQqtq1qnbNBhvPZziSJHVKnCur6ihgNUBV3QLcOuZz/xd4ZpL1kmwBPJam3/LrwNPavs4tgcdNPSHJC2gOfdm3qlZ3iFOSNCBh/q/HOZe6nADh+iR3AwogycOBq8d87r8DjwBOb5//6qr6RZJP89uLiZ4PnDqyzn8FzgO+2e6Uz1TVmzvEK0kaiAW6Huec6JI4X04zmna7JF+naUJ9+kxPqKq7tP8X8Kp2Gl2+Oskrq+q6Nil/GzijXeZFsiVJi844V0d5RlUdDfwS+D2ak74H+GFV/XrGJ4/nv5JsBtwReEs7SEiStIwstYrzdcDRwKer6iHAWXMZQFU9bi7XJ0kaluZYzOFkznES5xVJvghsm+TY6Qurau+5D0uStJwstYpzT+AhNCd5f+/8hiNJ0uI2zvU4bwa+leSRVXXZ2h6X5J+q6m/mNDpJ0rIwoJbaTqfcW2vSbD1qwlgkSctQYMGvcDIJD/mQJPVuoU+bN4khxSpJUu/msuIcTp0tSVpUBtRSO37iTLJBVd04bd7Kqrq8vfuBOY1MkrQsJBlUH2eXptrvtOenBSDJ04BvTN2vqsPmMC5J0jKyJC9kDfwpcGiSE2kuLXY3YPf5CEqSpMWqy+EoZyR5G82JEK4FHltVF8xbZJKkZWOpnTkIgCQfAbYDHgT8Ds3J2f+pqg6ar+AkSUvfUj6O8wzgBe0lwn6aZDfgffMTliRpORlQ3uzUVPv+afevBp4/5xFJkrSIdWmq3R54O7AjsMHU/Kq63zzEJUlaLjKsPs4uh6N8FPggcAvweODfgMPnIyhJ0vKSCf8tpC6J885VdQKQqjqvqg6kueSYJEmz1gwOmmxaSF0GB92UZAVwTpIXAxcCd5mfsCRJWpy6VJwvBTYEXgLsAjwbeM58BCVJWl6WasVZNCc/uC+wfjvvQzTHdUqSNGsZ0PEoXRLnJ4BX0RzPuXp+wpEkLTdTfZxD0SVxXlZVx85bJJIkDUCXxHlAkg8DJwA3Tc2sqs/MeVSSpOWjhyucTKJL4nwu8ACa/s2pptoCTJySpIks1XPVPrSqdpi3SCRJGoAuh6N8I8mO8xaJJGlZWsonQHg4cFqSn9L0cQaoqvJwFEnSRAbUUtspce4xb1FIkpaxsGKBzzc7iS6XFTtvPgORJGkIulSckiTNubB0m2olSZp7A7sep4lTktS7IR3H2eVwFEmSlj0rTklSr+zjlCSpoyE11Zo4JUm9G1DeXFqJc/WvruP607/WdxhLzl8eulnfISxZPz/97L5DWJJ22+NhfYewJF1+7c19hzCxJOsBJwMXVtVes1nHkkqckqThCQs6UvWlwNnAJrNdgaNqJUn9CiSZaBprM8nWwJ7AhycJ14pTktS7BerifD/wamDjSVZixSlJWgpWJjl5ZNp/dGGSvYBLq+qUSTdkxSlJ6lVzPc6Ja87Lq2rXGZY/Ctg7yZOBDYBNkhxeVc/uuiErTklS7zLhtC5V9bqq2rqqtgGeBfzPbJImWHFKkhYBj+OUJGmRqqoTgRNn+3wTpySpZ+MfUrIYmDglSb1a4BMgTMzEKUnq3ZAqziEleUmSemfFKUnq3XDqTROnJKlvsalWkqQly4pTktQrR9VKktTRkJpqTZySpN4NJ20OqzqWJKl3VpySpN4NqKXWxClJ6lczOGg4mdPEKUnq3ZAqTvs4JUnqwIpTktSzEJtqJUka35Caak2ckqReDW1wkH2ckiR1YMUpSepXbKqVJKkTE6ckSR0MaVStfZySJHVgxSlJ6lWAFcMpOE2ckqT+Damp1sQpSerdkAYH2ccpSVIHVpySpN7ZVCtJ0pgcHCRJUifDujqKfZySJHVgxSlJ6pfnqpUkqZsB5U2baiVJ6mLwFWeS/YH9AVj/Lv0GI0nqrBlVO5yac/CJs6oOAQ4BWLHh3avncCRJszCctLkEEqckaQkYUOa0j1OSpA6sOCVJvRvSCRBMnJKk3g1obJCJU5LUvwHlTfs4JUnqwopTktS/AZWcJk5JUq+Cg4MkSRrfwE7ybh+nJEkdWHFKkno3oILTxClJWgQGlDlNnJKknmVQg4Ps45QkqQMrTklS74Y0qtbEKUnqVRhUF6eJU5K0CAwoc9rHKUlSByZOSVLvMuG/da4/uXeSLyf5fpKzkrx0trHaVCtJ6t0CDA66BXhFVZ2aZGPglCTHV9X3u67IilOS1LtMOK1LVV1cVae2t68FzgbuNZtYTZySpGUlyTbAg4GTZvN8m2olSf2am+NRViY5eeT+IVV1yO02ldwF+DTwsqq6ZjYbMnFKkno3B6fcu7yqdp1xG8n6NEnzE1X1mdluyKZaSdKSlyTAR4Czq+p9k6zLxClJ6lVoRtVOMo3hUcCfA7snOa2dnjybeG2qlST1br6PRqmqr83VZkyckqT+eco9SZKWJitOSVLvhnQhaxOnJKl3Xo9TkqQOBpQ37eOUJKkLK05JUv8GVHKaOCVJvWpOVTuczGnilCT1a/yz/ywK9nFKktSBFackqXcDKjhNnJKkRWBAmdPEKUnqWQY1OMg+TkmSOrDilCT1bkijapdU4qwbLrv8xtMOOq/vOMa0Eri87yDG8aPTDuo7hC4Gs18HaDD79ivH9R1BJ4PZr8B952OlYVBdnEsscVZt0XcM40pyclXt2nccS437df64b+eH+7U1oMxpH6ckSR0sqYpTkjRMQxpVa+LszyF9B7BEuV/nj/t2frhfcXCQxlBV/rHMA/fr/HHfzg/3a2NAedM+TkmSujBxLoAkhya5NMmZI/PumuT4JOe0/2/eZ4xDlGSDJN9OcnqSs5K8qZ2/bZKTkpyb5Mgkd+w71qFJ8rMkZyQ5LcnJ7Tw/sxNKskO7T6ema5K8bNnv2/bqKJNMC8nEuTAOA/aYNu+1wAlVtT1wQntf3dwE7F5VOwOrgD2SPBx4J/APVXV/4JfA83uMccgeX1WrRg6V8DM7oar6YbtPV4d9vgUAAAXkSURBVAG7AL8C/h33Lb89mnO208IxcS6AqvoqcOW02U8BPtbe/hjwxwsa1BJQjevau+u3UwG7A8e08923c8fP7Nx6AvDjqjoP9+2gmDj7s2VVXdze/gWwZZ/BDFWS9ZKcBlwKHA/8GLiqqm5pH3IBcK++4huwAr6Y5JQk+7fz/MzOrWcBR7S3l/W+DcNqqnVU7SJQVZWk+o5jiKrqVmBVks1omrwe0HNIS8Wjq+rCJHcHjk/yg9GFfmYn0/a77w28bvqy5bpvHVWrcVyS5J4A7f+X9hzPoFXVVcCXgUcAmyWZ+lG4NXBhb4ENVFVd2P5/Kc0PkofhZ3Yu/SFwalVd0t5f9vt2SBWnibM/xwJ/0d7+C+CzPcYySEm2aCtNktwZeBJwNk0CfXr7MPdtR0k2SrLx1G3g94Ez8TM7l/blt8204L4dFJtqF0CSI4DHASuTXAAcALwDOCrJ84HzgH36i3Cw7gl8LMl6ND8Cj6qq/0ryfeBTSd4KfBf4SJ9BDtCWwL+n+Rl/B+CTVfWFJN/Bz+zE2h8jTwL+cmT2sv8+GNIp91K17JrSJUmLyM4P3qWO+8q3JlrHPTe94ykLdZUZK05JUu+GU2/axylJUidWnJKkXvUxMnYSJk5JUu+GNDjIxClJ6t9w8qZ9nJIkdWHilOZBklVJnjxyf+8kc3LFi/YyVBvOxbqkxWI410YxcUrzZRXwm8RZVcdW1TvmaN0vAzolzvYkEdKi5Sn3pIFIsk2Ss5N8qL0Y9hfb0/et6bHbJflCe8WQ/03ygHb+M5Kc2V5Q+6vtCbzfDDyzvVjxM5Psl+Sf28cfluSDSb6V5CdJHtde7PzsJIeNbO+DSU6edpHulwBbAV9O8uV23r7tRafPTPLOkedfl+S9SU4HHpHkHUm+n+R7Sd4zP3tUmo1M/G8hmTgl2B44qKp2Aq4CnraWxx0C/E1V7QK8EviXdv4bgT9oL6i9d1Xd3M47sr1o8ZFrWNfmNCek/1ua85T+A7AT8MAkq9rH/L/2TCgPAn4vyYOq6h+Bi2guMv34JFvRXLh7d5oq96FJpq7luBFwUhvX2cCfADtV1YOAt3bdSZIaJk4JflpVp7W3TwG2mf6AJHcBHgkc3V7/82Cac+UCfB04LMkLgXGbRP+zmvNdngFcUlVnVNVq4KyR7e+T5FSa8+3uBOy4hvU8FDixqi5rr0H6CeCx7bJbgU+3t68GbgQ+kuSpwK/GjFOad16PUxqem0Zu3wqsqal2Bc0FsldNX1BVL0qyG7AncEqSXTpsc/W07a8G7pBkW5qq9qFV9cu2CXeDMdY76sb2eqVU1S1JHgY8gebKMS+mqVIldWTFKY2hqq4BfprkGQBp7Nze3q6qTqqqNwKXAfcGrgU2nmCTmwDXA1cn2ZLm+o1TRtf9bZpm3JXtAKB9ga9MX1lbMW9aVZ+jaR7eeYLYpDlnxSktTX8GfDDJ64H1gU8BpwPvTrI9TYvTCe28nwOvbZt13951Q1V1epLvAj8AzqdpDp5yCPCFJBe1/ZyvpbkGaYD/rqo1XctxY+CzSTZoH/fyrjFJanhZMUlSrx78kF3rxK9/e6J1bLbhel5WTJK0THiSd2nYkhwEPGra7A9U1Uf7iEda6vo4+88kTJzSNFX1133HIGnxMnFKkvo3oJLTxClJ6t2QrsfpcZySJHVgxSlJ6p2jaiVJ6mBAedOmWknSIrAAV7JOskeSHyY5d5ILy5s4JUlLXnsu54Nozvu8I7BvkjVdcWidTJySpN4twIWsHwacW1U/aa+Z+yngKbOJ1T5OSVKvpq7HOc/uRXPBhCkXALvNZkUmTklSr0499ZTj7rx+Vk64mg2SnDxy/5CqOmTCda6RiVOS1Kuq2mMBNnMhzbVyp2zdzuvMPk5J0nLwHWD7JNsmuSPwLODY2azIilOStORV1S1JXgwcB6wHHFpVZ81mXV7IWpKkDmyqlSSpAxOnJEkdmDglSerAxClJUgcmTkmSOjBxSpLUgYlTkqQOTJySJHXw/wER7+mPzT635wAAAABJRU5ErkJggg==\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + }, + { + "name": "stdout", + "output_type": "stream", + "text": [ + "{'ExtraTreesClassifier': {'max_features': 'sqrt', 'n_estimators': 70}, 'RandomForestClassifier': {'max_features': 'sqrt', 'n_estimators': 70}}\n" + ] + } + ], + "source": [ + "tuned_params = hyperparameter_optimization(\n", + " X, y, (extraTrees, extraTrees_param_dict), (rf, rf_param_dict)\n", + ")\n", + "\n", + "print(tuned_params)\n", + "\n", + "# extract values from dict - seperate each classifier's param dict\n", + "keys, values = zip(*tuned_params.items())\n", + "\n", + "# train test split\n", + "X_train, X_test, y_train, y_test = train_test_split(\n", + " X, y, test_size=0.33, random_state=42\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 21, + "metadata": {}, + "outputs": [], + "source": [ + "def get_accuracy(model, X_train, y_train, X_test, y_test):\n", + " \"\"\"\n", + " Given a model, train, and test data, \n", + " fit model and calculate accuracy of predictions.\n", + "\n", + " Parameters\n", + " ----------\n", + " model : sklearn-compliant classifier\n", + " X_train : numpy.ndarray\n", + " Train input data, shape (n_samples, n_features)\n", + " y_train numpy.ndarray\n", + " Train output data, shape (n_samples, n_outputs)\n", + " X_test: numpy.ndarray\n", + " Test input data, shape (n_samples, n_features)\n", + " y_test:numpy.ndarray\n", + " Test output data, shape (n_samples, n_outputs)\n", + "\n", + " Returns\n", + " -------\n", + " accuracy : float\n", + " An sklearn metric for model performance.\n", + " \"\"\"\n", + "\n", + " model.fit(X_train, y_train)\n", + " predictions = model.predict(X_test)\n", + " accuracy = metrics.accuracy_score(y_test, predictions)\n", + " return accuracy\n", + "\n", + "\n", + "# get accuracies of optimized and default models\n", + "extraTrees_models = [ExtraTreesClassifier(**values[0]), ExtraTreesClassifier()]\n", + "extraTrees_acc = []\n", + "for model in extraTrees_models:\n", + " extraTrees_acc.append(get_accuracy(model, X_train, y_train, X_test, y_test))\n", + "\n", + "rf_models = [RandomForestClassifier(**values[1]), RandomForestClassifier()]\n", + "rf_acc = []\n", + "for model in rf_models:\n", + " rf_acc.append(get_accuracy(model, X_train, y_train, X_test, y_test))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Plotting the result" + ] + }, + { + "cell_type": "code", + "execution_count": 22, + "metadata": {}, + "outputs": [ + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeIAAAEYCAYAAACA4PXcAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8li6FKAAAgAElEQVR4nO3deZgV1bX38e+iEZEAoohGaAQUB1AQI2gM1ziAqARBcUKuIhr1TYLGqJiQa4JIHEgkyVXJNdFEESdQVByiUSLORgEVEUGZRGkkCCJgO4RpvX/sfdrqw+kJuru6Ob/P8/TTp2rXsKpOVa3au4Zj7o6IiIiko0HaAYiIiOQzJWIREZEUKRGLiIikSIlYREQkRUrEIiIiKVIiFhERSVEqidjM9jKzYjMr2Mrxi81s72qO6Xkzu6A6p5k1/T3M7EUz+9zMfr+N0xpqZi9XV2yy7czsZTMbmnYcVWVmp5lZUdynuqQdTz4ws3vMbFT8fLSZvbuV02loZm5m7asxtiIzO7q6pldX5lXXVSoRxwP/O2b2pZn928xuNbMWlZ2JmS0xs96Zbnf/yN2buvumrQk6jrt4a8atCjP7i5ldFJd/UzxYFZvZB2Z2p5ntV4XJXQSsApq7+xXVHKebWcdyyrPjz/y1rsS0t+kExcz+nJjfejPbkOh+amunW5+Y2bWJ5V5jZq+Y2WGJ8t5mtjnru3mkjGndE9djsZmtNrNnqrgdZvs98P/iPvXONkxHtoK7P+/uB6YdR32RY18pMrNJZnZoFaZxrZmNr8EwqzyfChOxmV0B/Ba4EtgZ+C7QDphqZo22Ic764ETgyfj5X+7elLAOegNfAW+Y2UGVnFY7YK6n9waVf8WDbfLv422dqJk1LK/c3X+UmR9wPTApMf8Tqzq9euzeuA5aAS8BD2aVf5T13ZxSzrSuj9NqC6wG7qhqMLE21SBOY2trZFvVoiXpq+f72Udx+28GHAEsAF6p17Vrdy/zD2gOFANnZPVvCqwEzo/do4DJwCTgc+BN4OBYdjewmZC4ioGfA+0BBxrGYZ4HrgVejcM8DrQE7gXWATOA9on5O9ARaB2Hz/x9GRapZLjzgXnAZ8DTQLtE2XHAe8BaYBzwAnBBorwrMDt+Hgq8nGP9PAFMTnR/Ny7DGuBt4OjYfzywAVgf4+wNHAb8Kw67PMbQKA5fav0k1tEF2fEAL8Zhv4jTPjNHnDnjj2X7EA7m34ndreN3ezRwHbAJ+DpOe1xi/Q8j7AAfxH43AUvj9/UGcGSOeY0C7snq1zFO7zzgI2Ba7N8TeC2un1nA9xPjtADujOutCBgNNIhl+8V1spbQAnFfGcvdgLDN/jvO43mgU6L8HuBm4CnCNv0voEOi/ATg/Tifm4BXgKFlzOtaYHzWtuXALrG7N7CkvH0xK65Rie4BwJrEMv0PsCgu+8TEPLLX83PxO81sO+/H4Q4k7AtrgHeAH2TN+0/AP+I4R8d+txD2r+K47veI/dYQ9r+DE9P4FbA4rtN3gf6JsgvivP8Yx10M9EmUtyTsS8sJ+/RDibL+hH1uDfAycFA56/C/gJnxu5sOHJ4oexm4hrAffx6XdddypjWQsH2uAxYCfYCzgNezhvt5Jt7kd5j93RO258vjul8L3A/smCgfQdhmlwE/jN9f+1jWGPgDYT9cAfwf0Dg5n7h9/Bu4s4zlKQJ+wTfHzb9l5h/X/5OE48NnhON0m8quO8Jx6EPCtjkizitzjPwuIW+si7HfWEZ8OfcV4M/Aa4nucXH6mfzxvdi/H+E4vIGwvb6R2PbmxbgXUToX7B6Xew3hWPlioqwQeCSukw+AYeXNp8ztqIKd/gRgI4mEkCi7C7g/cYDdAJwG7AAMj0HtEMuXAL0T47Zny0S8kJAUdgbmAvPjSm8ITEhuOHHcjjliujcR04A4zU5xGr8CXo1lu8UVnon3sricyZU/ArghsQHlSsTnAyvi5zbAp0BfwgHxuNjdKpaPB65NjHsoYeNrGNfHPOBnudZPYh1tkYjLWx9ZO0DORBzLL4zrvAnhgDo213yz5jcV2BXYKfY7m7CjNgSuIOzsjbPGG0XZifjOOP+dCLW0T4Hj47o8gbDztozjPE44yDQhHPTfAH4Yyx4kHEgaEA5MPctY5gZxvTSLw40DZibK74nz7B63kUmZ2Ak7ZjFwSiy7Mm4/FSZiYEdgLPAJUFDewaWMaSUP4s1iXM/F7isIJwRt4jL9Fbi7nPXckNIH8kaE/fbncbl6x+XsmJj3Z4RaSIO4LPfEZTkkzvOFOI3BQAEwBpiaiP8MYM84/uA4/T0SB8MNhP2qALgEWJoY92ngPmCXGN/3Y/8ehIN3jzje+YSDaaMc6283QoI7Ky7/OYRtLXPC8jLhBHPfuJ5eIrHfZk3re4SDc6+4PG2B/eO6XQPsmxj2HWBAju8wVyJ+Dfg2YX+azzf7fT/CSUhn4FvAA1nf3y2EpLALoRL1JPCbxHw2ElqlGhH32xzLVATMJiSY3WIsmVhbEbb5neL0H6Z0RaTMdQd0id91z7jd3BzjySTiGcBZie368DLiKysR9yFUGjInHucQjk8NCceDZXxzQlHqxDj2OwnYGzDgWELFsWssu5FwfNghrrvMdteAcBL2P7F/R0Ku61XWfMrcryvY6c8G/l1GWckORjjAJs9GGsQN5sjYvYSKE/FVifLfA09lraRZWYmgY1Y8vyAckDOJ4SniwTkR05eEJuIhWfFa3ACTifilRPxDyZ2ITwA2JOZ/d1b508C58fN4ytihY/nPgEdyrZ/EOtqWRLyRcHDI/C3KGuYxwsFiNqXPwEvmmzW/YyvYdj4jURNKbCdlJeK9Ev2uIuuMHXgW+G9CkvkqK8Zz+GZbvA+4lcSZeqV2hHDQceBbsfse4M+J8v7AnPj5/Kz1n9ney0vE6+N630Q4e07W8HsTWo2S38/AMqZ1D6GFItOSMoVYUyccBI9KDNs2DtugjPWcnYiPIRywLDHMg8CvEvO+I0c8tya6LwPeSXQfAqwqZ73PIda6CYn4vURZ8xjfbnFZNgI755jG7cDVWf0WkeMkjNAi8GpWvxnA2fHzy8CIRNlPgSfKiP1vlF1zux24Jn7uRjipy1RMKkrEgxLdf+CblqgJlD6Z75z5/uJ3/DWlW/2OBBYk5vM1OU5OsuLOPg72J7aW5Bi2O7Ay0V3muiO0Wt2TKGtK2BeOjt2vAiOJJ9vlxFdWIj4oros9cpQZoeJ1YGJ/HF/BfJ7gm9rt9YSTjn2yhukJLM7q92vg9srOJ/NX0TXiVcBuZVxP2DOWZyzNfHD3zfELrfBmoIQVic9f5ehuWtaIZnYicClwsrt/FXu3A26KN8dkmhSMcCBvnRWvJ7vjjWgHEDaO8rSJ083M7/TM/OI8/4uwnnLFvJ+ZPRFvfltH+LJ3q2B+2+I1d2+R+Nsnq/x2wsZ8i7v/pxLTW5rsMLPhZjbPzNbGZd+Zqi1PcnrtgLOy1uV3Cd9bO8IZ9YpE2Z8INWMItcIdgJnxBsNzc83MzArM7Hdmtjiu/4WxKBnzvxOfv+SbbTB7+8ls7+W5z91bEGo67xMSVNJHWd/Pw+VMa0wcZk93P9ndP4j99wIeT6yXzM1XuyfGLfW9ZWkd4/BEvw8J23l541d63403Dr6diPEAyl/nxPHbEhL62hzzbwf8Imt72TMr7uQyfpjVL3sZy/res7UlJPxc7iKcOEKo0Exy9w1lDJutUtsdpZfj24T9Irlun6D0d7/C3ddnOuKNfpmbns5MDJc9j9Zx+KZm9lcz+yjuM9PYch+v7D5TzDfHTggnSJ2B981supn1pWraEE5m18ZYf25m75nZWkKl4Fs5Yi1hZv3M7PV4A+QaQg07M/wYwnp41swWmdmVsX87YK+s7e7nhO+iSiq6YP8v4D+E6yAPJIJuSriR6X8Sw7ZNlDcgNG1kbgZK7tjVysz2J2z0A909uQEtBa5z93tzjLNvVryW7CY0iU7ziu/qPoVQc87M7253v7CSod8KvEVojvnczH5GaCqHcP0NQvPOuvi5yl9uZcXv838JZ/ijzOwhd8/sJGV9dyX9zexIwgbYC3jX3Teb2WeEE59KyTr4LyXUiH+cI9a2hB1815gAs6eznFCzwsy+T7ip8MVEssoYQriMcCxhJ2tJqKlWJublhNaQTEyZ7b1C7r7SzC4CXjezie6+osKRKq8IGOzur2cXWLyrPms9Z/sYaGtmlhhuL0IrScZW78sWHjm8lbCdvO7um8xsDpVb50sJlYLm7r4uR9k17v7bSkznY+AHWf32IrQsVNVSwuW0Lbj7y2aGmfUkNMEP3IrpZ1tO6ePUXonPKwitLvuXs02V+u7cvU8Zw2XPI3McvxLoABzm7v82s+6E1oTKxt4h0xGPObsmYnkfGBT3pdOBh8xsF3f/upLTPwWY4e5fm9kxhOvsvQiX3CAk6Mx2Vmo9mNlOhPtFBgF/d/cNZvZEZvi4vV0GXGbhEb/nzGw64ftf4O6dyoip0vtKuTXiePZ5DXCLmZ1gZjtYeGbtAcJOf3di8EPNbGCsPf+MkMBfi2UrCO3v1crMmgOPEpq1s5+r/TPwSzM7MA67s5mdHsv+DhyYiPenlE50feMwueZZYGYdzOwWws0q18Sie4CTzOz4OExjC88IlnWAbkZIssVmdgBQknTcfSWhifDsOK3zKWOHj7Z1/d5EuD56AWG5/1zFaTcjNBuuBBqa2UhCs+LWuhs4xcyOS6zLY8ysdTzZegEYa2bNzayBmXWMSRczO8PMMrWbNYSdIdcJVTPCNvop4YTnuirE9wTQzcwGmFnmHoNWlR3Z3ecSmtqHV2GelfFn4Hoz2wvAzHY3s/5VGP9Vwvd4RdzXjyXsC5OqKb6mhO9jZQjPLiTUiCsUv/d/An8ysxYxvu/H4tuBYWbWw4KmZnaSmX0rx6SeIOz7Z1q4c3wwodk+5/5egb8BF8Rts4GZFcaKQcbdhBOPYnd/LfckquQB4HwzOyAu29WZglhp+Cvwv2bWKq6HQjMrK9mW52Iza2NmLYFf8s3334xwEvxZLBtZhWk+CAwwsyPMbEdCs23yZP4cM9stnlyvjWVbnGgnJZbxGsLlt0zFMHM8WkVoHRtFqBFnrADaxwoYhJaERoTtcpOZ9SMk8cx8TjKzfeLwawnHk82Eiup6M7siHqMKzKyLffMoVfZ8ylTh40vu/ru4gGMJieN1wplAr6wmzEeBMwnNAOcQaqiZppgbgF9ZqL5X58HnO4SbI/5oiWcwY9yPEB67mmihGWUOoRaPu68inHWNIRyI9yXc5JKpHR9PuOMv6Yg47XWE66bNgR4en72MB4oBhHW1krCOrqTsdTyccKb8OeFAkn2wuzCO/ynhTtbymslHAXfF9XtGGcMcYVs+R9zDzAYQaneZE4HLge+YWaZZ7SbgNDP7zMxuLmPaTxPW13xC7fJrym8CLZe7LyGc4f6asC4/IjQ5Z9bl2YQday5he3uQb06kDgdmmNkXhOs6w9z9oxyzuZNwpv8x4e7dii5DJONbQdjWbyTs7HsR9ouquBH4sZlV5+WIPxC+h2fN7HPCMvWo7Mhxfz6JsB2vItxQM9jdF1RHcO4+m3BD0XRCDWl/qrbezo7/5xMOcpfE6b5G2H5vJWwP8xPDZsewknDd8xeEfesyoJ+7f1bFxcHdXyXspzcTDtDPUbo2OYFwuefuLceuOnd/nHAZ5gXCMk7NGuQKwv43PcbzDOHYVlX3E056FhEuo1wf+/+BcMnpU8K2Ven3AMTv/lLCycQyQhN2shm7LzAvbrdjCU9/rN9iQsFe8VhcTNh+OhPuuZgWy5+M8S8g3J+0jrC9ZUwiJN7VZjbd3dcQtoNHCM3lpxFO2DL2JzTDFxPyxE3u/pK7b4xxHxbnswr4C99UQkrNp7z1Y+W3VFWOhbfEdHT3nBt/fWLhRQvj3P2wCgcWESlDrLV+QniUKvvSiEgJvWs6t6srHkREpFzDgFeUhKUiSsRZ3H26u+dscjGzO8zsEws3mOQqNzO72cwWmtlsM/tOouxcM1sQ/3LeySsi2wczKyI0l1f3fQCyHaqWpul8EW8OKQYmuPsWr7a0cMv9JYTrBocTriUcbma7Et7k051wE8IbwKFbc11KRES2L6oRV4G7v0jpZ9+yDSAkaY83kLQwsz0JN39NdffVMflOJfH4i4iI5K/6/OLvuqgNpe8WLor9yuq/BQvPmF4E8K1vfevQAw6o1NMdIiIVeuONN1a5e6UftZPaoURcx7j7bcBtAN27d/eZM2emHJGIbC/MLPutYlIHqGm6ei2j9HOEhbFfWf1FRCTPKRFXr8eAIfHu6e8Ca+MrF58G+pjZLma2C+E9pk+nGaiIiNQNapquAjO7n/Bay93i4wlXE16hhrv/mfBGl76EHxD4kvAic9x9tZn9hm/eyzo68S5nERHJY0rEVeDuZ1VQ7oSH+HOV3QHcURNxidRXGzZsoKioiK+/ruy7/aUyGjduTGFhITvssEPaoUglKBGLSGqKiopo1qwZ7du3pxLvxpdKcHc+/fRTioqK6NChQ8UjSOp0jVhEUvP111/TsmVLJeFqZGa0bNlSrQz1iBKxiKRKSbj6aZ3WL0rEIiIiKdI1YhGpM9qP+Hu1Tm/JmB9UOExRURHDhg1j7ty5bN68mX79+nHjjTfSqFGjnMOvWbOG++67j5/85CcAfPzxx/z0pz9l8uTJlY5r5MiRfP/736d3796VHieXpk2bUlxcvE3TkPSpRiwiecvdGThwICeffDILFixg/vz5FBcXc9VVV5U5zpo1a/i///u/ku7WrVtXKQkDjB49epuTsGw/VCMWkbw1bdo0GjduzHnnnQdAQUEBf/zjH+nQoQMdOnTg6aefZu3atSxbtoyzzz6bq6++mhEjRrBo0SK6devGcccdx7Bhw+jXrx9z5sxh/PjxTJkyhS+++IIFCxYwfPhw1q9fz913382OO+7Ik08+ya677srQoUPp168f7du354ILLgBg06ZNzJkzB3dn0aJFDBs2jJUrV9KkSRNuv/12DjjgAD744AMGDx5McXExAwYMSHPVSTVSjVhE8ta7777LoYceWqpf8+bN2Wuvvdi4cSPTp0/noYceYvbs2Tz44IPMnDmTMWPGsM8++zBr1ixuvPHGLaY5Z84cHn74YWbMmMFVV11FkyZNeOuttzjiiCOYMGFCqWG7d+/OrFmzmDVrFieccALDh4efL77ooou45ZZbeOONNxg7dmxJM/ill17Kj3/8Y9555x323HPPGlorUttUIxYRKcNxxx1Hy5YtARg4cCAvv/wyJ598crnjHHPMMTRr1oxmzZqx8847c9JJJwHQpUsXZs+enXOcSZMm8eabb/LMM89QXFzMq6++yumnn15S/p///AeAV155hYceegiAc845h1/84hfbvIySPiViEclbnTt33uL67rp16/joo49o2LDhFo8BVeaxoB133LHkc4MGDUq6GzRowMaNG7cYfs6cOYwaNYoXX3yRgoICNm/eTIsWLZg1a1bO6evRpO2PmqZFJG/16tWLL7/8sqTJeNOmTVxxxRUMHTqUJk2aMHXqVFavXs1XX33FlClT6NmzJ82aNePzzz+vlvmvWbOGs846iwkTJtCqVfiZ4ObNm9OhQwcefPBBINxQ9vbbbwPQs2dPJk6cCMC9995bLTFI+lQjFpE6ozKPG1UnM+ORRx7hJz/5Cb/5zW/YvHkzffv25frrr+f+++/nsMMO49RTT6WoqIizzz6b7t27AyEhHnTQQZx44okMG5bz9fKV8uijj/Lhhx9y4YUXlvSbNWsW9957Lz/+8Y+59tpr2bBhA4MGDeLggw/mpptuYvDgwfz2t7/VzVrbEQu/UyB1Uffu3X3mzJlphyFSY+bNm0enTp3SDiOn8ePHM3PmTMaNG5d2KFsl17o1szfcvXtKIUkZ1DQtIiKSIjVNi4jkMHToUIYOHZp2GJIHVCMWERFJkRKxiIhIipSIRUREUqRELCIikiLdrCUidceonat5emsrHKSgoIAuXbqwYcMGGjZsyJAhQ7jsssto0KD8esqVV17Jk08+Sd++fXO+c7oimZ8wXLJkCa+++iqDBw+u8jRk+6BELCJ5baeddip5neQnn3zC4MGDWbduHddcc0254912222sXr2agoKCbZr/kiVLuO+++5SI85iapkVEot13353bbruNcePG4e5s2rSJK6+8kh49etC1a1f+8pe/ANC/f3+Ki4s59NBDmTRpEo8//jiHH344hxxyCL1792bFihUAjBo1irFjx5ZM/6CDDmLJkiWl5jlixAheeuklunXrxh//+MdaW1apO1QjFhFJ2Hvvvdm0aROffPIJjz76KDvvvDMzZszgP//5Dz179qRPnz489thjNG3atKQm/dlnn/Haa69hZvz1r3/ld7/7Hb///e8rNb8xY8YwduxYnnjiiZpcLKnDlIhFRMrwzDPPMHv27JJfaFq7di0LFiygQ4cOpYYrKirizDPPZPny5axfv36LcpHyKBGLiCQsXryYgoICdt99d9ydW265heOPP77ccS655BIuv/xy+vfvz/PPP8+oUaMAaNiwIZs3by4Z7uuvv67J0KWe0jViEZFo5cqV/OhHP+Liiy/GzDj++OO59dZb2bBhAwDz58/niy++2GK8tWvX0qZNGwDuuuuukv7t27fnzTffBODNN9/kgw8+2GLc6vxZRamfVCMWkbqjEo8bVbevvvqKbt26lTy+dM4553D55ZcDcMEFF7BkyRK+853v4O60atWKKVOmbDGNUaNGcfrpp7PLLrtw7LHHliTcU089lQkTJnDggQdy+OGHs99++20xbteuXSkoKODggw9m6NChXHbZZTW7wFLn6GcQ6zD9DKJs7+ryzyDWd/oZxPpDTdMiIiIpUiIWERFJkRKxiKRKl8eqn9Zp/aJELCKpady4MZ9++qkSRzVydz799FMaN26cdihSSbprWkRSU1hYSFFREStXrkw7lO1K48aNKSwsTDsMqSQlYhFJzQ477KC3UEneU9O0iIhIipSIRUREUqRELCLV6h//+Af7778/HTt2ZMyYMVuUf/jhh/Tq1YuuXbty9NFHU1RUBMBzzz1Ht27dSv4aN25c8harI488sqR/69atOfnkk2t1mURqkt6sVYfpzVpS32zatIn99tuPqVOnUlhYSI8ePbj//vvp3LlzyTCnn346/fr149xzz2XatGnceeed3H333aWms3r1ajp27EhRURFNmjQpVXbqqacyYMAAhgwZUivLtD3Rm7XqJtWIq8DMTjCz981soZmNyFHezsyeNbPZZva8mRUmyn5nZu+a2Twzu9nMrHajF6l506dPp2PHjuy99940atSIQYMG8eijj5YaZu7cuRx77LEAHHPMMVuUA0yePJkTTzxxiyS8bt06pk2bphqxbFeUiCvJzAqAPwEnAp2Bs8ysc9ZgY4EJ7t4VGA3cEMf9HtAT6AocBPQAjqql0EVqzbJly2jbtm1Jd2FhIcuWLSs1zMEHH8zDDz8MwCOPPMLnn3/Op59+WmqYiRMnctZZZ20x/SlTptCrVy+aN29eA9GLpEOJuPIOAxa6+2J3Xw9MBAZkDdMZmBY/P5cod6Ax0AjYEdgBWFHjEYvUQWPHjuWFF17gkEMO4YUXXqBNmzYUFBSUlC9fvpx33nkn528A33///TkTtEh9pueIK68NsDTRXQQcnjXM28BA4CbgFKCZmbV093+Z2XPAcsCAce4+L9dMzOwi4CKAvfbaq3qXQKSGtWnThqVLv9lNioqKSn6nN6N169YlNeLi4mIeeughWrRoUVL+wAMPcMopp7DDDjuUGm/VqlVMnz6dRx55pAaXQKT2qUZcvYYDR5nZW4Sm52XAJjPrCHQCCgkJ/VgzOzLXBNz9Nnfv7u7dW7VqVVtxi1SLHj16sGDBAj744APWr1/PxIkT6d+/f6lhVq1axebNmwG44YYbOP/880uVl1XrnTx5Mv369dOrG2W7o0RcecuAtonuwtivhLt/7O4D3f0Q4KrYbw2hdvyauxe7ezHwFHBE7YQtUnsaNmzIuHHjOP744+nUqRNnnHEGBx54ICNHjuSxxx4D4Pnnn2f//fdnv/32Y8WKFVx11VUl4y9ZsoSlS5dy1FFb3kJR1nVjkfpOjy9Vkpk1BOYDvQgJeAYw2N3fTQyzG7Da3Teb2XXAJncfaWZnAhcCJxCapv8B/K+7P17ePPX4kohUJz2+VDepRlxJ7r4RuBh4GpgHPODu75rZaDPLtL0dDbxvZvOBPYDrYv/JwCLgHcJ15LcrSsK1qSZewPDDH/6Qgw8+mK5du3LaaadRXFxcq8skIlJfqEZch9VGjbimXsCwbt26kkdMLr/8cnbffXdGjNji0WsRqUWqEddNums6zyVfwACUvIAhmYjnzp3LH/7wByC8gCHXyxSyX8CQScLuzldffYXeX1L/tB/x97RDqFZLxvwg7RBEclLTdJ6ryRcwnHfeeXz729/mvffe45JLLqmhJRARqd+UiKVCW/sChjvvvJOPP/6YTp06MWnSpNoOW0SkXlAiznNVeQHDW2+9xXXXhfvPKvMCBoCCggIGDRrEQw89VENLIFKzauJmxnHjxtGxY0fMjFWrVtXq8kjdo0Sc52riBQzuzsKFC0s+P/bYYxxwwAE1vCQi1W/Tpk0MGzaMp556irlz53L//fczd+7cUsMMHz6cIUOGMHv2bEaOHMkvf/lLINxPMWvWLGbNmsW0adNo0qQJffr0AaBnz57885//pF27drW+TFL3KBHnuZp4AYO7c+6559KlSxe6dOnC8uXLGTlyZK0vm8i2qqlfkzrkkENo3759jccv9YPumhb69u1L3759S/UbPXp0yefTTjuN0047Lee47du33+LmrgYNGvDKK69Uf6AitSzXzYyvv/56qWEyNzNeeumlpW5mbNmyZckwEydO5PLLL6+1uKV+UY1YRGQbbMuvSYmAasQiImWqyV+TEslQjVhEpAw1+WtSIhlKxCIiZaipX5O6+eabKSwspKioiK5du3LBBRfU6nJJ3aJ3TddhVX3XtF5JKNVJ29P2R++arptUI5btyta+fAHgo48+ok+fPnTq1InOnTuzZBL1dvYAABCYSURBVMkSAI488siSlzK0bt0657u2RUS2lhKxbDe25eULAEOGDOHKK69k3rx5TJ8+nd133x2Al156qeTFDEcccQQDBw6s1eUSke2bErFsN7bl5Qtz585l48aNHHfccQA0bdq05OULGevWrWPatGmqEYtItdLjS7Ld2JaXL8yfP58WLVowcOBAPvjgA3r37s2YMWNKPQ86ZcoUevXqVfITj1LPjNo57Qiq16i1aUcg1UQ1YskrZb18YePGjbz00kuMHTuWGTNmsHjxYsaPH19qXD2GIiI1QYlYthvb8ktShYWFdOvWjb333puGDRty8skn8+abb5aMt2rVKqZPn84PfqA7b0WkeikRy3ZjW16+0KNHD9asWcPKlSsBmDZtGp07dy4Zb/LkyfTr14/GjRvX0tKISL5QIpbtxra8fKGgoICxY8fSq1cvunTpgrtz4YUXlkx74sSJapYWkRqhF3rUYXqhh5qB07TdbU+NB6cdQvXaipu19EKPukk1YhERkRQpEYuIiKRIiVhERCRFSsQiIiIp0pu1pO7Sm5BEJA+oRiwiIpIiJWIREZEUKRGLiIikSIlYREQkRUrEIiIiKVIiFhERSZESsYiISIqUiEVERFKkRCwiIpIiJWIREZEUKRGLiIikSIlYREQkRUrEVWBmJ5jZ+2a20MxG5ChvZ2bPmtlsM3vezAoTZXuZ2TNmNs/M5ppZ+9qMXURE6iYl4koyswLgT8CJQGfgLDPrnDXYWGCCu3cFRgM3JMomADe6eyfgMOCTmo9aRETqOiXiyjsMWOjui919PTARGJA1TGdgWvz8XKY8JuyG7j4VwN2L3f3L2glbRETqMiXiymsDLE10F8V+SW8DA+PnU4BmZtYS2A9YY2YPm9lbZnZjrGFvwcwuMrOZZjZz5cqV1bwIIiJS1ygRV6/hwFFm9hZwFLAM2AQ0BI6M5T2AvYGhuSbg7re5e3d3796qVataCVpERNKjRFx5y4C2ie7C2K+Eu3/s7gPd/RDgqthvDaH2PCs2a28EpgDfqZ2wRUSkLsvLRGxml5jZLlUcbQawr5l1MLNGwCDgsazp7mZmmXX6S+COxLgtzCxTxT0WmLt10YuIyPYkLxMxsAcww8weiI8kWUUjxJrsxcDTwDzgAXd/18xGm1n/ONjRwPtmNj/O47o47iZCs/SzZvYOYMDt1b1QIiJS/zRMO4A0uPuvzOzXQB/gPGCcmT0A/M3dF5Uz3pPAk1n9RiY+TwYmlzHuVKBrNYQvIiLbkXytEePuDvw7/m0EdgEmm9nvUg1MRETySl7WiM3sUmAIsAr4K3Clu2+I13cXAD9PMz4REckfeZmIgV2Bge7+YbKnu282s34pxSQiInkoX5umnwJWZzrMrLmZHQ7g7vNSi0pERPJOvibiW4HiRHdx7CciIlKr8jURW7xZCwhN0uRvM72IiKQoXxPxYjP7qZntEP8uBRanHZSIiOSffE3EPwK+R3hFZRFwOHBRqhGJiEheysvmWHf/hPCKShERkVTlZSI2s8bAD4EDgcaZ/u5+fmpBiYhIXsrXpum7gW8DxwMvEH5J6fNUIxIRkbyUr4m4o7v/GvjC3e8CfkC4TiwiIlKr8jURb4j/15jZQcDOwO4pxiMiInkqL68RA7fF3yP+FeE3hZsCv043JBERyUd5l4jjDzusc/fPgBeBvVMOSURE8ljeNU3Ht2jp15VERKROyLtEHP3TzIabWVsz2zXzl3ZQIiKSf/KuaTo6M/4flujnqJlaRERqWV4mYnfvkHYMIiIikKeJ2MyG5Orv7hNqOxYREclveZmIgR6Jz42BXsCbgBKxiIjUqrxMxO5+SbLbzFoAE1MKR0RE8li+3jWd7QtA141FRKTW5WWN2MweJ9wlDeFkpDPwQHoRiYhIvsrLRAyMTXzeCHzo7kVpBSMiIvkrXxPxR8Byd/8awMx2MrP27r4k3bBERCTf5Os14geBzYnuTbGfiIhIrcrXRNzQ3ddnOuLnRinGIyIieSpfE/FKM+uf6TCzAcCqFOMREZE8la/XiH8E3Gtm42J3EZDzbVsiIiI1KS8TsbsvAr5rZk1jd3HKIYmISJ7Ky6ZpM7vezFq4e7G7F5vZLmZ2bdpxiYhI/snLRAyc6O5rMh3u/hnQN8V4REQkT+VrIi4wsx0zHWa2E7BjOcOLiIjUiLy8RgzcCzxrZncCBgwF7ko1IhERyUt5mYjd/bdm9jbQm/DO6aeBdulGJSIi+Shfm6YBVhCS8OnAscC8dMMREZF8lFeJ2Mz2M7Orzew94BbCO6fN3Y9x93EVjI6ZnWBm75vZQjMbkaO8nZk9a2azzex5MyvMKm9uZkWJ55dFRCTP5VUiBt4j1H77uft/ufsthPdMV8jMCoA/AScSfjbxLDPrnDXYWGCCu3cFRgM3ZJX/BnhxG+IXEZHtTL4l4oHAcuA5M7vdzHoRbtaqjMOAhe6+OL6beiIwIGuYzsC0+Pm5ZLmZHQrsATyzDfGLiMh2Jq8SsbtPcfdBwAGERPkzYHczu9XM+lQwehtgaaK7KPZLepuQ7AFOAZqZWUszawD8HhheUYxmdpGZzTSzmStXrqx4oUREpF7Lq0Sc4e5fuPt97n4SUAi8BfyiGiY9HDjKzN4CjgKWEZq+fwI86e5FlYjtNnfv7u7dW7VqVQ0hiYhIXZaXjy8lxbdq3Rb/yrMMaJvoLoz9ktP6mFgjju+xPtXd15jZEcCRZvYToCnQyMyK3X2LG75ERCS/5H0iroIZwL5m1oGQgAcBg5MDmNluwGp33wz8ErgDwN3/OzHMUKC7krCIiECeNk1vDXffCFxMePnHPOABd3/XzEYnftv4aOB9M5tPuDHrulSCFRGRekM14ipw9yeBJ7P6jUx8ngxMrmAa44HxNRCeiIjUQ6oRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpERcBWZ2gpm9b2YLzWxEjvJ2Zvasmc02s+fNrDD272Zm/zKzd2PZmbUfvYiI1EVKxJVkZgXAn4ATgc7AWWbWOWuwscAEd+8KjAZuiP2/BIa4+4HACcD/mlmL2olcRETqMiXiyjsMWOjui919PTARGJA1TGdgWvz8XKbc3ee7+4L4+WPgE6BVrUQtIiJ1mhJx5bUBlia6i2K/pLeBgfHzKUAzM2uZHMDMDgMaAYtyzcTMLjKzmWY2c+XKldUSuIiI1F1KxNVrOHCUmb0FHAUsAzZlCs1sT+Bu4Dx335xrAu5+m7t3d/furVqp0iwisr1rmHYA9cgyoG2iuzD2KxGbnQcCmFlT4FR3XxO7mwN/B65y99dqJWIREanzVCOuvBnAvmbWwcwaAYOAx5IDmNluZpZZp78E7oj9GwGPEG7kmlyLMYuISB2nRFxJ7r4RuBh4GpgHPODu75rZaDPrHwc7GnjfzOYDewDXxf5nAN8HhprZrPjXrXaXQERE6iI1TVeBuz8JPJnVb2Ti82Rgixqvu98D3FPjAYqISL2jGrGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxCIiIilSIhYREUmRErGIiEiKlIhFRERSpEQsIiKSIiViERGRFCkRi4iIpEiJWEREJEVKxFVgZieY2ftmttDMRuQob2dmz5rZbDN73swKE2XnmtmC+Hdu7UYuIiJ1lRJxJZlZAfAn4ESgM3CWmXXOGmwsMMHduwKjgRviuLsCVwOHA4cBV5vZLrUVu4iI1F1KxJV3GLDQ3Re7+3pgIjAga5jOwLT4+blE+fHAVHdf7e6fAVOBE2ohZhERqeMaph1APdIGWJroLiLUcJPeBgYCNwGnAM3MrGUZ47bJNRMzuwi4KHYWm9n72x56/WSwG7Aq7TiqzTWWdgR5TdsTAO2qOwzZdkrE1Ws4MM7MhgIvAsuATVWZgLvfBtxW/aHVP2Y20927px2HbB+0PUldpURcecuAtonuwtivhLt/TKgRY2ZNgVPdfY2ZLQOOzhr3+ZoMVkRE6gddI668GcC+ZtbBzBoBg4DHkgOY2W5mllmnvwTuiJ+fBvqY2S7xJq0+sZ+IiOQ5JeJKcveNwMWEBDoPeMDd3zWz0WbWPw52NPC+mc0H9gCui+OuBn5DSOYzgNGxn5RPTfRSnbQ9SZ1k7p52DCIiInlLNWIREZEUKRGLiIikSIlYqszMNpnZrMTfFq/7zBr+f7ZiHo/EaS80s7WJeX1v6yOXfJHYRueY2eNm1iL2b29mX2Vtv43Sjlfym64RS5WZWbG7N93W4c3MCNvg5nLGPRoY7u79yihvGG+kEymR3ObM7C5gvrtfZ2btgSfc/aA04xNJUo1YqoWZ7Rx/EGP/2H2/mV1oZmOAnWLN495YI3nfzCYAc4C2Znarmc00s3fN7JpKzKvIzMaY2VvAKWa2r5k9bWZvmNmLZrZfHG4PM3s4Tnu6mX039j/WzN6OMb1pZt+quTUjdcC/KONNdiJ1gV7oIVtjJzOblei+wd0nmdnFwHgzuwnYxd1vBzCzi929W/zcHtgXONfdX4v9rnL31fGHNZ41s67uPruCGD5x90Pi+M8BF7j7IjPrCYwjPKt9M/A7d38tUxMCDgKuBC5y99fji1e+3uY1InVS3KZ6AX9L9N4nsf2+4u7Daj8ykW8oEcvW+CqTWJPcfaqZnU74laqDyxn/w0wSjs6I79huCOxJ+PGMihLxJIB47e+7wEOhpRv4ZrvuDeyf6L+Lme0EvALcZGb3Ag+5e3EF85L6J3Oy2Ibw3P/URNmiXNuvSFqUiKXaxLeKdQK+BHYh/LhFLl8kxulAeEd3D3f/zMzGA40rMbvMNAxYVcaB1YDD4q9lJV1rZo8BPwBeM7Ne7r6gEvOU+uMrd+9mZk0IL+EZRmghEalzdI1YqtNlhNrHYOBOM9sh9t+Q+JytOSGprjWzPQi/91xp8Wcll5vZKRBOBswsUxv/J+EATCzLNI/v4+6z3f0G4E1g/6rMU+oPd/8S+ClwhZmp4iF1khKxbI3MzVeZvzHxJq0LgCvc/SXCr0/9Kg5/GzA7NgWX4u5vA28B7wH3EZqNq2oQ8CMzext4F8jcYT0M6Glms81sLnBh7D88PtYyGygGntmKeUo94e5vES51nJV2LCK56PElERGRFKlGLCIikiIlYhERkRQpEYuIiKRIiVhERCRFSsQiIiIpUiIWERFJkRKxiIhIiv4/2YTjNqHgNaEAAAAASUVORK5CYII=\n", + "text/plain": [ + "
" + ] + }, + "metadata": { + "needs_background": "light" + }, + "output_type": "display_data" + } + ], + "source": [ + "labels = [\"ExtraTrees\", \"RF\"]\n", + "x = np.arange(len(labels))\n", + "width = 0.35\n", + "fig, ax = plt.subplots()\n", + "rects1 = ax.bar(x - width / 2, extraTrees_acc, width, label=\"Optimized\")\n", + "rects2 = ax.bar(x + width / 2, rf_acc, width, label=\"Default\")\n", + "\n", + "# Add some text for labels, title and custom x-axis tick labels, etc.\n", + "ax.set_ylabel(\"Accuracy\")\n", + "ax.set_title(\n", + " \"Optimized/Default ExtraTrees and RF Performance on cylinder-bands Dataset\"\n", + ")\n", + "ax.set_xticks(x)\n", + "ax.set_xticklabels(labels)\n", + "ax.legend()\n", + "\n", + "def autolabel(rects):\n", + " \"\"\"Attach a text label above each bar in *rects*, displaying its height.\"\"\"\n", + " for rect in rects:\n", + " height = float(\"%.3f\" % (rect.get_height()))\n", + " ax.annotate(\n", + " \"{}\".format(height),\n", + " xy=(rect.get_x() + rect.get_width() / 2, height),\n", + " xytext=(0, 3), # 3 points vertical offset\n", + " textcoords=\"offset points\",\n", + " ha=\"center\",\n", + " va=\"bottom\",\n", + " )\n", + "\n", + "autolabel(rects1)\n", + "autolabel(rects2)\n", + "fig.tight_layout()\n", + "plt.ylim((0.9, 1))\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.7.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/examples/ensemble/ExtraTrees_vs_RF_Gridsearch.py b/examples/ensemble/ExtraTrees_vs_RF_Gridsearch.py new file mode 100644 index 0000000000000..4bc5b3bc9f538 --- /dev/null +++ b/examples/ensemble/ExtraTrees_vs_RF_Gridsearch.py @@ -0,0 +1,291 @@ +""" +============================================================================ +Demonstration of warmstart grid search to compare classifier performance +============================================================================ +An important step in classifier performance comparison is hyperparameter +optimization. Here, we specify the classifer models we want to tune and a +dictionary of hyperparameter ranges (preferably similar for fairness in +comparision) for each classifier. Then, we find the optimal hyperparameters +through a function that implements warmstart grid search and refit the +optimized models to obtain accuracies. The performance of each hyperparameter +value pairing is visualized in heatmaps. + +In this example, we tune hyperparameters for two classifiers, Random Forest +and Extra Trees, and compare their performance on an OpenML-CC18 benchmarking +suite dataset (https://www.openml.org/d/15). We can see clearly in the resulting +plot that the optimized models perform better than or atleast similar to the +default parameter models. On the dataset we use in this example, RF performs +marginally better than ExtraTrees overall. +""" +print(__doc__) + +import pandas as pd +import numpy as np + +import sklearn +from sklearn import metrics +from sklearn.metrics import cohen_kappa_score, make_scorer +from sklearn.model_selection import cross_val_score +from sklearn.model_selection import train_test_split +from sklearn.ensemble import ExtraTreesClassifier +from sklearn.ensemble import RandomForestClassifier +from sklearn.datasets import fetch_openml + +import matplotlib +import matplotlib.pyplot as plt + +from warnings import simplefilter + +simplefilter(action="ignore", category=FutureWarning) +from warnings import simplefilter + +simplefilter(action="ignore", category=FutureWarning) + + +def hyperparameter_optimization(X, y, *argv): + """ + Given any number of classifier types and + a dictionary of two hyperparameters to tune for each classifier, + find optimal pairs of hyperparameters. + + Parameters + ---------- + X : numpy.ndarray + Input data, shape (n_samples, n_features) + y : numpy.ndarray + Output data, shape (n_samples, n_outputs) + *argv : list of tuples (classifier, hyperparameters) + List of (classifier, hyperparameters) tuples: + + classifier : sklearn-compliant classifier + For example sklearn.ensemble.RandomForestRegressor + hyperparameters : dictionary of hyperparameter values or range + + Returns + ------- + clf_best_params : dictionary + Dictionary of classifiers and their respective optimal hyperparameters + """ + + best_models = {} + + # Iterate over all (classifier, hyperparameter dict) pairs + for clf, params in argv: + best_params = grid_search(X, y, clf, params) + best_models[clf.__class__.__name__] = best_params + + return best_models + + +def grid_search(X, y, clf, params): + """ + Given a classifier and two hyperparameters and the + range/values to search for each, find optimal hyperparameter + values using warmstart grid search parameter sweeps. + + Parameters + ---------- + X : numpy.ndarray + Input data, shape (n_samples, n_features) + y : numpy.ndarray + Output data, shape (n_samples, n_outputs) + clf : sklearn-compliant classifier + For example sklearn.ensemble.RandomForestRegressor + params : dictionary of hyperparameter values or range + + Returns + ------- + best_params : dictionary + Dictionary of best hyperparameters + """ + param1_name = list(params.keys())[0] + param2_name = list(params.keys())[1] + param1 = params[param1_name] + param2 = params[param2_name] + + # sweep over all pairs of parameter combinations and collect mean scores + kappa_scorer = make_scorer(cohen_kappa_score) + mean_scores = np.zeros((np.shape(param1)[0], np.shape(param2)[0])) + for idx1, val1 in enumerate(param1): + clf.max_features = val1 # change .max_features to .name of 1st parameter + for idx2, val2 in enumerate(param2): + clf.n_estimators = val2 # change .n_estimators to .name of 2nd parameter + score = cross_val_score(clf, X, y, scoring=kappa_scorer, cv=5) + mean_scores[idx1][idx2] = np.mean(score) + + # select parameter pair with highest kappa score + best_idx1, best_idx2 = np.unravel_index( + np.argmax(mean_scores, axis=None), np.shape(mean_scores) + ) + best_params = {param1_name: param1[best_idx1], param2_name: param2[best_idx2]} + + # generate heatmap + param_heatmap(params, mean_scores, clf.__class__.__name__) + + return best_params + + +def param_heatmap(params, scores, clf_name): + """ + Given a dictionary of two parameter ranges, scores + for each pair of parameter values, and classifier name, + generate heatmap showing model performance scores for each + pair of parameter values. + + Parameters + ---------- + params : dictionary of hyperparameter ranges + scores : ndarray + Scores for each parameter value pair + clf_name : string + Name of sklearn-compliant classifier + """ + param1_name = list(params.keys())[0] + param2_name = list(params.keys())[1] + param1 = params[param1_name] + param2 = params[param2_name] + + scores = -np.array(scores) + scores = scores.ravel().argsort().argsort().reshape(scores.shape) + plt.figure(figsize=(8, 6)) + plt.subplots_adjust(left=0.2, right=0.95, bottom=0.15, top=0.95) + plt.imshow(scores, interpolation="nearest", cmap=plt.cm.Blues) + plt.xlabel(param2_name) + plt.ylabel(param1_name) + plt.colorbar() + plt.xticks(np.arange(len(param2)), param2) + plt.yticks(np.arange(len(param1)), param1) + plt.title("Grid Search Kappa Rank " + clf_name) + plt.show() + + +############################################################################### +# Building classifiers and specifying parameter ranges or values to search +# ---------------------------------------------------------- +# + +# get some data +X, y = fetch_openml(data_id=40979, return_X_y=True, as_frame=True) +y = pd.factorize(y)[0] +X = X.apply(lambda x: pd.factorize(x)[0]) +n_samples, n_features = np.shape(X) + +# build a classifier with warm_start=True +extraTrees = ExtraTreesClassifier(warm_start=True) + +# specify parameters and ranges or values to search +extraTrees_param_dict = { + "max_features": ["sqrt", "log2", None], + "n_estimators": [10, 30, 50, 70], +} + +# build another classifier with warm_start=True +rf = RandomForestClassifier(warm_start=True) + +# specify parameters and ranges or values to search +rf_param_dict = { + "max_features": ["sqrt", "log2", None], + "n_estimators": [10, 30, 50, 70], +} + +############################################################################### +# Obtaining best parameters dictionary and refitting +# ---------------------------------------------------------- +# + +tuned_params = hyperparameter_optimization( + X, y, (extraTrees, extraTrees_param_dict), (rf, rf_param_dict) +) + +print(tuned_params) + +# extract values from dict - seperate each classifier's param dict +keys, values = zip(*tuned_params.items()) + +# train test split +X_train, X_test, y_train, y_test = train_test_split( + X, y, test_size=0.33, random_state=42 +) + + +def get_accuracy(model, X_train, y_train, X_test, y_test): + """ + Given a model, train, and test data, + fit model and calculate accuracy of predictions. + + Parameters + ---------- + model : sklearn-compliant classifier + X_train : numpy.ndarray + Train input data, shape (n_samples, n_features) + y_train numpy.ndarray + Train output data, shape (n_samples, n_outputs) + X_test: numpy.ndarray + Test input data, shape (n_samples, n_features) + y_test:numpy.ndarray + Test output data, shape (n_samples, n_outputs) + + Returns + ------- + accuracy : float + An sklearn metric for model performance. + """ + + model.fit(X_train, y_train) + predictions = model.predict(X_test) + accuracy = metrics.accuracy_score(y_test, predictions) + return accuracy + + +# get accuracies of optimized and default models +extraTrees_models = [ExtraTreesClassifier(**values[0]), ExtraTreesClassifier()] +extraTrees_acc = [] +for model in extraTrees_models: + extraTrees_acc.append(get_accuracy(model, X_train, y_train, X_test, y_test)) + +rf_models = [RandomForestClassifier(**values[1]), RandomForestClassifier()] +rf_acc = [] +for model in rf_models: + rf_acc.append(get_accuracy(model, X_train, y_train, X_test, y_test)) + +############################################################################### +# Plotting the result +# ---------------------------------------------------------- +# + +labels = ["ExtraTrees", "RF"] +x = np.arange(len(labels)) +width = 0.35 +fig, ax = plt.subplots() +rects1 = ax.bar(x - width / 2, extraTrees_acc, width, label="Optimized") +rects2 = ax.bar(x + width / 2, rf_acc, width, label="Default") + +# Add some text for labels, title and custom x-axis tick labels, etc. +ax.set_ylabel("Accuracy") +ax.set_title( + "Optimized/Default ExtraTrees and RF Performance on cylinder-bands Dataset" +) +ax.set_xticks(x) +ax.set_xticklabels(labels) +ax.legend() + + +def autolabel(rects): + """Attach a text label above each bar in *rects*, displaying its height.""" + for rect in rects: + height = float("%.3f" % (rect.get_height())) + ax.annotate( + "{}".format(height), + xy=(rect.get_x() + rect.get_width() / 2, height), + xytext=(0, 3), # 3 points vertical offset + textcoords="offset points", + ha="center", + va="bottom", + ) + + +autolabel(rects1) +autolabel(rects2) +fig.tight_layout() +plt.ylim((0.9, 1)) +plt.show()