RandomForest_MNIST_CIFAR10.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification of images in CIFAR-10 database by Random Forest Classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "from collections import OrderedDict\n",
    "from sklearn.ensemble import RandomForestClassifier\n",
    "from sklearn.metrics import precision_recall_fscore_support\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "from sklearn.metrics import accuracy_score\n",
    "import matplotlib.image as mpimg\n",
    "import matplotlib.pyplot as plt\n",
    "from time import time\n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings('ignore')#mute the warning messages"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preprocessing\n",
    "\n",
    "The dataset labeled as \"train\" has 50000 images saved as png files and a csv file with labels to identify each of these files are also available. We will use this dataset and split it into a training set and testing set for this project.\n",
    "\n",
    "The major task in the data preprocessing is to extract RGB from the png files and format the data as formats that are readily readable by the classifer."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Loading the dataset...\n",
      "Raw data loaded!\n",
      "Labels loaded!\n"
     ]
    }
   ],
   "source": [
    "#total number of images\n",
    "img_num = 50000\n",
    "#read the train set\n",
    "print (\"Loading the dataset...\")\n",
    "X=[]\n",
    "for i in range(1, img_num+1):\n",
    "    filename ='train/'+str(i)+'.png'\n",
    "    data = mpimg.imread(filename)\n",
    "    X.append(data)\n",
    "print (\"Raw data loaded!\")\n",
    "\n",
    "#load the labels \n",
    "label_o = pd.read_csv(\"trainLabels.csv\")\n",
    "class_le = LabelEncoder()\n",
    "label_ = class_le.fit_transform(label_o['label'])\n",
    "y = np.array(label_)\n",
    "print (\"Labels loaded!\")\n",
    "\n",
    "#split the training and testing datasets (80% training, 20% testing)\n",
    "train_X=np.array(X[:40000]).reshape(-1,3072)\n",
    "test_X=np.array(X[40000:]).reshape(-1,3072)\n",
    "train_y=np.array(y[:40000])\n",
    "test_y=np.array(y[40000:])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Model optimization\n",
    "\n",
    "We first try to understand how some the of the key parameters affect the performance of a random forest classifier with the training set in this project, and based on the result a optimal set of parameter are selected.\n",
    "\n",
    "Within the evaluated range, it seems that after n_estimators is larger than 400, the reduction of error rate becomes relatively incremental. Evaluating the error rate for even larger n_estimator quickly becomes expensive in terms of time, and may be less \"cost-effective\". For this analysis, we picked n_estimator = 500. We would like to keep the evaluation of even larger n_estimators for future work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RandomForestClassifier, max_feature='sqrt'  fit for  10 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  60 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  110 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  160 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  210 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  260 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  310 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  360 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  410 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  460 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  510 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  10 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  60 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  110 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  160 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  210 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  260 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  310 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  360 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  410 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  460 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  510 estimators \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAELCAYAAAAoUKpTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3XecVOXZ//HPtb1TtlAEWcqySEea\niJFdC4qgmIgaAwlIHo0lsSV5jHliTZ7HJJJoMCa2KLafGk1UYjBWimADBJW2tKXD7lK2wfa9fn+c\nM7OzfYDdnWXner9e5zUzZ87M3HNg97t3OfctqooxxhjTlJBAF8AYY0z7Z2FhjDGmWRYWxhhjmmVh\nYYwxplkWFsYYY5plYWGMMaZZFhbGGGOaZWFhjDGmWRYWxhhjmhUW6AK0lKSkJE1NTQ10MYwx5pSy\nevXqg6qa3NxxHSYsUlNTWbVqVaCLYYwxpxQR2enPcdYMZYwxplkWFsYYY5plYWGMMaZZHabPwjSv\noqKCPXv2UFpaGuiiGGPaWFRUFL169SI8PPyEXm9hEUT27NlDfHw8qampiEigi2OMaSOqyqFDh9iz\nZw99+/Y9ofewZqggUlpaSmJiogWFMUFGREhMTDypVgULiyBjQWFMcDrZn/2OExbFuVBhbfHGGNMa\nOk5YFO6FncsDXQpjjOmQOk5YSAhseT/QpTDNCA0NZeTIkQwdOpRLL72U/Pz8FnnfHTt2MHTo0BZ5\nrzlz5tC3b19GjhzJyJEjmT9/fou8b0OWLFnCJ598Umvf888/z9ChQxkyZAiDBw9m3rx53nK9/vrr\nLfK5+/btY8aMGd7H11xzDcOHD+fhhx/mnnvu4YMPPmiRz2lLmzZtYuTIkYwaNYpt27Yd9+sfeeQR\njh071gola1xLTVG0du1aFi1a5H28YMEC7rvvvhZ5b4+OExaRcbDlvUCXwjQjOjqatWvXsm7dOrp2\n7cpjjz0W6CI16KGHHmLt2rWsXbuWW265xe/XVVVVHdfn1A2Ld955h0ceeYT33nuP9evX8+WXX9Kp\nU6fjek9/9OzZ0xs8Bw4c4JNPPuHrr7/m9ttv54EHHuCCCy7w+70qKytbvHwn4s0332T69OmsWbOG\n/v37H/frTyQs2sN3r6ysrBcWraHjDJ2NTIDD2+HQNkg8/v8oweb+f61nw77CFn3PwT0TuPfSIX4f\nP2HCBL7++msAiouLmT59OkeOHKGiooLf/OY3TJ8+nR07djBlyhTOOeccPvnkE0477TTeeustoqOj\nWb16NXPnziUmJoZzzjnH+76lpaXceOONrFq1irCwMP74xz+SmZnJggULePPNN6mqqmLdunX89Kc/\npby8nBdeeIHIyEgWLVpE165dGy3vyy+/zP/93/+hqkydOpXf/e53AMTFxXHHHXfw7rvv8oc//IHo\n6GjuuOMOiouLSUpKYsGCBfTo0YP58+fz+OOPExYWxuDBg/ntb3/L448/TmhoKC+++CKPPvooDz74\nIPPmzaNnz56AMzb+uuuuq1eWBx54gH/961+UlJRw9tln88QTTyAi9T7jlVdeYenSpdx6662A08m5\nbNkyDh06xLRp01i3bh2TJ08mNzeXkSNH8uijj/K3v/2NadOmMWPGDFavXt3gd8nIyODss89mxYoV\nXHbZZfz0pz9t8Jz5e86feuopnnzyScrLyxkwYAAvvPACMTExTJ8+nSuuuIIf/OAHPPHEEyxbtoyX\nXnqp3ucsWrSIRx55hNDQUJYtW8bixYt58cUXmT9/PuXl5YwfP56//OUvhIaGcuONN7Jy5UpKSkqY\nMWMG999/P/Pnz2ffvn1kZmaSlJTE4sWLiYuLo7i4GIDXX3+dt99+mwULFjBnzhy6du3KmjVrOPPM\nM3nggQf4yU9+wjfffENlZSX33Xcf06dP9+dHgORkZ/6+o0ePctVVV7Fnzx6qqqq4++67ufrqq/nP\nf/7DbbfdRlJSEmeeeSbbt2/n7bff5r777mPfvn3s2LGDpKQkli9fTklJCcuXL+euu+4iOjqauLg4\nv8rgN1XtENvoEUNV701Q/fSvahq2YcMG7/37Fq7Tqx7/pEW3+xaua7YMsbGxqqpaWVmpM2bM0Hfe\neUdVVSsqKrSgoEBVVfPy8rR///5aXV2t2dnZGhoaqmvWrFFV1SuvvFJfeOEFVVUdNmyYLlmyRFVV\nf/azn+mQIUNUVXXevHk6Z84cVVXduHGj9u7dW0tKSvTZZ5/V/v37a2Fhoebm5mpCQoL+9a/O/5fb\nbrtNH374YVVVnT17tqampuqIESN0xIgR+vXXX+vevXu1d+/empubqxUVFZqZmalvvPGGqqoC+uqr\nr6qqanl5uU6YMEFzc3NVVfWVV17Ra6+9VlVVe/TooaWlpaqqeuTIEVVVvffee/Whhx7ynp8uXbpo\nfn5+g+du9uzZ+tprr6mq6qFDh7z7Z82apQsXLmz0M6ZNm6bLly9XVdWioiKtqKjQ7Oxs7/nyve/7\nOU19l0mTJumNN97YYDl9+XvODx486H3N//zP/+j8+fNVVfXAgQPav39/XbZsmaalpdX63nX5nssN\nGzbotGnTtLy8XFVVb7zxRn3uuedqnbvKykqdNGmSfvXVV6qq2qdPH83Ly/O+n+f/qqrqa6+9prNn\nz/aen6lTp2plZaWqqt51113e/5NHjhzRtLQ0LS4u1k2bNnn/D9XdPP82Hq+//rr+13/9l/dxfn6+\nlpSUaK9evXTz5s1aXV2tV155pU6dOtX7Xc8880w9duyY9zzffPPNTf1TeM9LXcAq9eN3bMepWYRF\nQmKa0xR11g2BLk27dzw1gJZUUlLCyJEj2bFjB6NHj+bCCy8EnD9afvnLX7Js2TJCQkLYu3cvOTk5\nAN7+A4DRo0ezY8cOCgoKyM/PZ9KkSQB8//vf55133gFg+fLl/OQnPwFg0KBB9OnTh82bNwOQmZlJ\nfHw88fHxdOrUiUsvvRSAYcOGeWs54DRD+bbpv/XWW2RkZHj/Epw5cybLli3j8ssvJzQ0lCuuuAKA\nrKws1q1b5/1eVVVV9OjRA4Dhw4czc+ZMLr/8ci6//PKTOo+LFy/m97//PceOHePw4cMMGTKESy+9\ntMHPmDhxInfccQczZ87kO9/5Dr169fLrM5r6LgBXX321X+/jzzlft24dv/rVr8jPz6e4uJiLLroI\ngG7duvHAAw+QmZnJG2+80WTNz9eHH37I6tWrGTt2LOD8v0tJSQHg73//O08++SSVlZXs37+fDRs2\nMHz4cL/e1+PKK68kNDQUgPfee4+FCxd6+5ZKS0vZtWsXZ5xxBmvXrvXr/YYNG8bPfvYz7rzzTqZN\nm8a3vvUt1q5dS9++fUlLSwNg1qxZPPnkk97XXHbZZURHRx9XuU9GxwkLgLTJsPJpKD8KEbGBLo1p\ngKfPoqCggGnTpvHYY49xyy238NJLL5GXl8fq1asJDw8nNTXVewFRZGSk9/WhoaGUlJSgqo2OG3f+\nWGqY73uFhIR4H4eEhDTZ/tzUe0ZFRXl/cagqQ4YM4dNPP6133L///W+WLVvGwoUL+fWvf8369evr\nHTNkyBBWr17Neeed1+jnlZaWctNNN7Fq1Sp69+7Nfffd5z1XDX3GL37xC6ZOncqiRYs466yz+OCD\nD4iKimr0/X2/c2PfBSA21r+fMX/O+Zw5c3jzzTcZMWIECxYsYMmSJd7XfPPNNyQmJrJv3z6/Ps9T\n9tmzZ/Pggw/W2p+dnc28efNYuXIlXbp0Yc6cOY1eqOb7/6vuMb7fXVX5xz/+QXp6eq1jsrKyGg3U\nJUuW0LlzZ+/jgQMHsnr1ahYtWsRdd93F5MmTueyyy5q8NsLf899SOk4HN0DaBVBVBtkfB7okphmd\nOnVi/vz5zJs3j4qKCgoKCkhJSSE8PJzFixezc2fTU+x37tyZTp06sXy5M1zatx373HPP9T7evHkz\nu3btqveDfLzGjx/P0qVLOXjwIFVVVbz88sveWo2v9PR08vLyvL9gKyoqWL9+PdXV1ezevZvMzEx+\n//vfe/+Cjo+Pp6ioyPv6u+66i//+7//mwIEDAJSVldUbjeX5xZWUlERxcbG3o7qxz9i2bRvDhg3j\nzjvvZMyYMWzatMmv79zYd2nIG2+8wV133eXX+zakqKiIHj16UFFRUevf8osvvuCdd95hzZo1zJs3\nj+zsbL/e7/zzz+f1118nNzcXgMOHD7Nz504KCwuJjY2lU6dO5OTkeGujQL1/i27durFx40aqq6t5\n4403Gv2siy66iEcffdT7B8WaNWsA5/x5BknU3XyDApzRaTExMcyaNYuf/exnfPnllwwaNIjs7Gzv\nyK6XX3650TLULXtr6DBhkVtURmnPsyA8BrbaENpTwahRoxgxYgSvvPIKM2fOZNWqVYwZM4aXXnqJ\nQYMGNfv6Z599lptvvpkJEybUqo7fdNNNVFVVMWzYMK6++moWLFhQ66/bE9GjRw8efPBBMjMzGTFi\nBGeeeWaDnZgRERG8/vrr3HnnnYwYMYKRI0fyySefUFVVxaxZsxg2bBijRo3i9ttvp3Pnzlx66aW8\n8cYbjBw5ko8//phLLrmEm2++mQsuuIAhQ4YwevToejWezp07c9111zFs2DAuv/xyb1NLY5/xyCOP\nMHToUEaMGEF0dDRTpkzx6zs39l0asm3bNhISEo7zrNb49a9/zfjx47nwwgu9//ZlZWVcd911PPPM\nM/Ts2ZM//OEPzJ07t8lansfgwYP5zW9+w+TJkxk+fDgXXngh+/fvZ8SIEYwaNYohQ4Ywd+5cJk6c\n6H3N9ddfz5QpU8jMzATgt7/9LdOmTeO8886r1fxW1913301FRQXDhw9n6NCh3H333cf9/b/55hvG\njRvHyJEj+d///V9+9atfERUVxZNPPsnUqVM555xz6NOnT6Ovz8zMZMOGDYwcOZJXX331uD/fH+LP\niT8VRPZI0/8sWUHml7dCzjq49WuwqS1q2bhxI2eccUagi2E6oFmzZvHwww97+3RMy1uyZAnz5s3j\n7bffPuH3aOh3gIisVtUxzb22w9QsRGBpVh6kXQj5u+Dg5kAXyZig8eKLL1pQdHAdpoM7LiKMJVm5\ncK4zcoMt70HyybVTG2Pan5tvvpkVK1bU2nfrrbdy7bXXBqhEbSMjI4OMjIyAfX6HCYv4qHB2HDpG\ndmVX+iaf4YTF2T8JdLGMMS2svV7139F1mGao+Cgn95Zk5TpNUTs/hbLWHR1gjDHBosOERURYCP2S\nYlmSledcb1FdAduXBrpYxhjTIXSYsACYlJ7MZ9sPUdpjLETE28SCxhjTQjpUWGSkp1BWWc2nOwuh\nf6YzZXkHGRpsjDGB1KHCYnzfrkSFh7Bkk9tvUbQPchq+4tQEhq1nUZutZ9FyTvX1LE5mltiZM2eS\nnp7O0KFDmTt3LhUVFS1Quto6VFhEhYcyoV8iSzbnwQB3CK1dzd2u2HoWtdl6Fi0nWNezACcsNm3a\nxDfffENJSQlPP/10i39GhwoLgMxBKew8dIzs8gToPsxWz2vMO7+AZ6e27PbOL46rCBMmTGDv3r2A\ns57F+eefz5lnnsmwYcN46623AKfGcMYZZ3DdddcxZMgQJk+eTElJCQCrV69mxIgRTJgwoVbolJaW\ncu2113qnvVi8eDHgrK1w+eWXc+mll9K3b1/+/Oc/88c//pFRo0Zx1llncfjw4SbL+/LLLzNs2DCG\nDh3KnXfe6d0fFxfHPffcw/jx4/n0009ZvXo1kyZNYvTo0Vx00UXs378fgPnz5zN48GCGDx/Od7/7\nXXbs2MHjjz/Oww8/7J3u43jWsxg7dixDhw7l+uuv906BUfczAJYuXeqtJY0aNYqioqJaNTHf9Sw+\n/vjjWjWYxr5LRkYGv/zlL5k0aRJ/+tOfGj1n/p7zp556irFjxzJixAiuuOIK7y/t6dOn8/zzzwPw\nxBNPMHPmzAY/x7OexdNPP+2druPFF1/0TqHxox/9yBvkN954I2PGjGHIkCHce++93vPmWc/C83rf\nv/Rff/115syZAzg1vDvuuIPMzEzuvPNOjh49yty5cxk7diyjRo3y/t/1R0MXMqoqP//5zxk6dCjD\nhg3zTt9RXV3NTTfdxJAhQ5g2bRqXXHKJ99/pkksuQUQQEcaNG8eePXv8LoPf/JnH/FTYRo8eraqq\nOw8e1T53vq3PLN+u+sH9qvd1UT1We+74YFVrLvtFd6o+c0nLbovubLYMtp6FrWdh61nUX8/C9/Ne\nf/11veCCC7SyslIPHDigvXv31n379ulrr72mU6ZM0aqqKt2/f7927tzZ+//Bo7y8XEeNGqXLli1r\n8BzZehY+Tk+M8Q6hvfaCyfDxH2D7Yhjy7UAXrX2Z8tuAfKytZ2HrWdh6Fk1bvnw511xzDaGhoXTr\n1o1JkyaxcuVKli9fzpVXXklISAjdu3f31oB83XTTTZx77rl861vfOu7PbU6HCwtwhtC+9PkuSlLO\nIzqqk9MUZWHRLth6Fraeha1nUX89i7plP579Hvfffz95eXk88cQTTR53ojpcnwVAZnoK5ZXVfLaj\nAPqfD1s/gOrqQBfL+LD1LGw9C1/Bvp6Fr3PPPZdXX32Vqqoq8vLyWLZsGePGjeOcc87hH//4B9XV\n1eTk5NQK1Keffpp3332Xl19+mZCQ1vm13iFrFuP6diU6PJQlWblkpk2G9f+EA19Dz5GBLprxUXc9\ni0svvZQxY8YwcuRIv9ezmDt3LjExMd5mC3Cq4jfccAPDhg0jLCysxdezUFUuueSSJtezuOWWWygo\nKKCyspLbbruNgQMHMmvWLAoKClDVWutZzJgxg7feeotHH32USy65hJycHC644AJv7Wnu3Lm1PsN3\nPYvU1NR661nU/Yy7776bxYsXExoayuDBg5kyZYq3o7opjX2XIUPqL8nbUutZ9OnTh2HDhlFUVORd\nz+LZZ5+ttZ7FRx991OQKclB7PYvq6mrCw8N57LHHOOuss7zrWfTr16/B9Sx69OjB4sWLvetZ9O7d\nm6FDh1JcXNzgZ919993cdtttDB8+HFUlNTX1pKYR//a3v82nn37KiBEjEBF+//vf0717d6644go+\n/PBDhg4dysCBAxk/frx3pNwNN9xAnz59mDBhAgDf+c53uOeee064DA3pMOtZjBkzRletWuV9PHfB\nSrbmFrP0piHIvDTI/BVM+nkASxh4tp6FaS22nkXbKC4uJi4ujkOHDjFu3DhWrFhB9+7d/X79yaxn\n0SFrFgAZ6cl8tCmX7JIY+vU805n6I8jDwpjW8uKLLwa6CEFh2rRp5OfnU15ezt13331cQXGyOm5Y\nDEwB1rMkK49+aRfCsofg2GGI8W80hTGmfQrW9SyAWv0Uba3DhsXpiTH0S45lyeY85k6eDEt/B9s+\ngmEzmn9xB9bUKCJjTgXt9ar/9u5kuxw65Ggoj4yBKXy2/RAlScMhJjHoZ6GNiori0KFDJ/2fxhhz\nalFVDh065NeQ6cZ02JoFOP0Wz6zI5rMd+WQOuMAdQlsFIaGBLlpA9OrViz179pCXlxfoohhj2lhU\nVJTfF2Q2pEOHhWcI7WLPENqvX4V9a6BXsx3/HVJ4eDh9+/YNdDGMMaegDt0MFRUeytn9E1mSlYf2\nywQJsYkFjTHmBLRqWIjIxSKSJSJbRaTelKQi8rCIrHW3zSKS7/PcbBHZ4m6zT7QMGenJ7Dp8jOxj\nkXDamKDvtzDGmBPRamEhIqHAY8AUYDBwjYgM9j1GVW9X1ZGqOhJ4FPin+9quwL3AeGAccK+IdDmR\ncmSkO5OHedfm3vclFOee4Lcyxpjg1Jo1i3HAVlXdrqrlwCtA/fkRalwDvOzevwh4X1UPq+oR4H3g\n4hMpRO+uzhDaxVnu6nkAWz88kbcyxpig1ZphcRqw2+fxHndfPSLSB+gLfHS8r/VHxsAUPs8+TEni\nEIjrZk1RxhhznFozLBq68quxAf7fBV5XVc+alH69VkSuF5FVIrKqqeGgmYOSKa+s5tPswzDgAufi\nvKr2sRyiMcacClozLPYAvX0e9wIam5D+u9Q0Qfn9WlV9UlXHqOqYpiYwq5mFNs9piirNh72rGj3e\nGGNMba0ZFiuBNBHpKyIROIGwsO5BIpIOdAF8V1h5F5gsIl3cju3J7r4TEhnmDKFdnJWL9ssACbWm\nKGOMOQ6tFhaqWgn8GOeX/Ebg76q6XkQeEJHLfA69BnhFfeagUNXDwK9xAmcl8IC774RlpCez+3AJ\n24vD4fSzLCyMMeY4tOoV3Kq6CFhUZ989dR7f18hrnwGeaamyOENonVlo+6ddCB/cB4X7IaFHcy81\nxpig16Gv4PbVu2sM/ZNjWZKVCwM8Q2g/CGyhjDHmFBE0YQFO7eLz7MMc65IO8T2tKcoYY/wUZGHh\nDqHdftgZFbVtMVRVBLpYxhjT7gVVWNQeQjsZyotg12eBLpYxxrR7QRUWkWGhTByQyJLNuWjfcyEk\n3JqijDHGD0EVFgCT0lOcIbRFIdBngnVyG2OMH4IuLDIGOld6L96U6zRF5W6A/N3NvMoYY4Jb0IWF\nZwjt0s1uvwXAVlsQyRhjmhJ0YQHuENrthzmW0A86n26r5xljTDOCMiwy01Mor/IMoZ0M25dCZVmg\ni2WMMe1WUIbF2L5diIlwh9AOuBAqjsLOTwJdLGOMabeCMixqzUKbeg6ERlpTlDHGNCEowwKcIbR7\njpSwrQBIPceutzDGmCYEbVh4htAuyXKH0B7aAoe3B7hUxhjTPgVtWPTuGsOAlDh3CK07C+0Wu0DP\nGGMaErRhAU7t4vPthzka1we69rPrLYwxphHBHRaeIbTbDjlNUdnLoKIk0MUyxph2J6jDwjuEdnOu\n0xRVWQo7lge6WMYY0+4EdVg4Q2iTWJKVh/aZCGHRNirKGGMaENRhAc6CSHuOlLDtSCX0m+SEhWqg\ni2WMMe2KhUW6ZwhtHgy4AI7sgEPbAlsoY4xpZ4I+LHp1cYbQOqvneYbQWlOUMcb4ajYsRKSbiPxN\nRN5xHw8WkR+2ftHaTmZ6Ml9kH+ZoTC9ISrewMMaYOvypWSwA3gV6uo83A7e1VoECofYQ2gth5woo\nKw50sYwxpt3wJyySVPXvQDWAqlYCVa1aqjY2JtUZQrvYM/VHVTns+DjQxTLGmHbDn7A4KiKJgAKI\nyFlAQauWqo3VGkJ7+lkQEWdNUcYY48OfsLgDWAj0F5EVwPPALa1aqgDISE9mb34J2w6XQ78MZ8py\nG0JrjDGAf2GxHpgEnA38CBgCbGrNQgVCrSG0aRdCwW7I63Bf0xhjTog/YfGpqlaq6npVXaeqFcCn\nrV2wttarSwxpKXFOv8UAG0JrjDG+Gg0LEekuIqOBaBEZJSJnulsGENNmJWxDGenJrMw+wtGobtBt\nqK2eZ4wxrqZqFhcB84BewB+BP7jbHcAvW79obc8zhPaTbYecq7l3fQqlhYEuljHGBFyjYaGqz6lq\nJjBHVTN9tstU9Z9tWMY2Mya1C7ERoTWr51VXwvYlgS6WMcYEXFhzB6jqP0RkKk7HdpTP/gdas2CB\nEBkWytkD3CG0l56DRHZy+i0GXxboohljTED5M93H48DVwE8AAa4E+rRyuQLGM4R266Ey6J8JWz+w\nIbTGmKDnz2ios1X1B8ARVb0fmAD0bt1iBU5GegrgGUI7GYr2Q866AJfKGGMCy5+wKHVvj4lIT6AC\n6Nt6RQqs0zpHk5YS56yeN+ACZ6cNoTXGBDl/wuJfItIZeAj4EtgBvNyahQq0zEEpziy0EYnQY4QN\noTXGBL0mw0JEQoAPVTVfVf+B01cxSFXvaZPSBUjGwGQqqpQVWw86TVG7P4eSI4EuljHGBEyTYaGq\n1TjXVngel6mq35MIisjFIpIlIltF5BeNHHOViGwQkfUi8v989leJyFp3W+jvZ7aEMaldnSG0m91+\nC62GbR+1ZRGMMaZd8acZ6j0RuUJE5HjeWERCgceAKcBg4BoRGVznmDTgLmCiqg6h9joZJao60t3a\ndOxqRFgIZw9IYmlWHtrzTIjuAls+aMsiGGNMu+LvrLOvAWUiUigiRSLiz2XN44CtqrpdVcuBV4Dp\ndY65DnhMVY8AqGrucZS9VWWmpzhDaA+WQP/zYev7UF0d6GIZY0xANBsWqhqvqiGqGqGqCe7jBD/e\n+zRgt8/jPe4+XwOBgSKyQkQ+E5GLfZ6LEpFV7v7L/fi8FuWZhda7INLRPNi/tq2LYYwx7YI/NYsT\n1VCzVd2r28KANCADuAZ42h15BXC6qo4Bvgc8IiL9632AyPVuoKzKy8truZIDPTtHM7BbnHO9xYDz\nAbFRUcaYoNWaYbGH2hfv9QL2NXDMW6paoarZQBZOeKCq+9zb7cASYFTdD1DVJ1V1jKqOSU5ObvEv\nkJGewsodhykO6wynjbbrLYwxQas1w2IlkCYifUUkAvguzop7vt4EMgFEJAmnWWq7iHQRkUif/ROB\nDa1Y1gZ5htB+4hlCu3c1HD3U1sUwxpiAa/Y6CxE5obkuVLUS+DHwLrAR+LuqrheRB0TEM7rpXeCQ\niGwAFgM/V9VDwBnAKhH5yt3/W1Vt87CoPYT2AkBh24dtXQxjjAm4JmedVdVqEflKRE5X1V3H++aq\nughYVGffPT73FWe01R11jvkEGHa8n9fSIsJCmDggiSWbctHpGUhMktMUNfyqQBfNGGPaVLNTlAM9\ngPUi8gVw1LOzra99CJSM9BTe25DDlrxjDEy7EDb/B6qrICQ00EUzxpg2409Y3N/qpWjHPENol2Tl\nOmHx1ctO30XvcQEumTHGtB1/rrNYCmwC4t1to7svKPTsHE16t3hnCG3/80BCbAitMSbo+LP40VXA\nFziLHl0FfC4iM1q7YO1JRnqyM4Q2JB56jbMhtMaYoOPP0Nn/Acaq6mx3EaRxwN2tW6z2ZVK67yy0\nFzpXchflBLpYxhjTZvwJi5A6czYd8vN1HcaYPu4QWs/qeeAst2qMMUHCn1/6/xGRd0VkjojMAf5N\nneGwHZ1nCO3SrFy021CI6+5MLGiMMUHCnw7unwNPAMOBEcCTqnpnaxesvckclMK+glI25x51mqK2\nfgRVlYEuljHGtIkmh866a1K8q6oXAP9smyK1T75DaNPTLoQ1L8CeL6DP2QEumTHGtL7mVsqrAo6J\nSKc2Kk+71aOTzxDafhkQEmZZTRA7AAAgAElEQVSjoowxQcOfPotS4BsR+ZuIzPdsrV2w9igjPZlV\nOw9TLLFw+gS73sIYEzT8CYt/4wyVXQas9tmCTkZ6Su0htDnroLDurOvGGNPxNDfrbChwoao+V3dr\no/K1K2NSuxAXGcYSz+p5YLULY0xQ8KfPItldjyLohYeGMHFAIkuy8tCkdEjoZf0Wxpig4M9EgjuA\nFSKykNqzzv6xtQrVnmWkp/Du+hw25x51RkV98xpUlkOY5akxpuPyp89iH/C2e2y8zxaUfIfQkjYZ\nyoth16cBLpUxxrSuZmsWqno/gIjEqurR5o7v6Hp0imZQ93gWZ+Xyo7POhdAI52rufpMCXTRjjGk1\n/sw6O8Fd9nSj+3iEiPyl1UvWjk1KT2bVjiMUaST0mWid3MaYDs+fZqhHgItwJhBEVb8Czm3NQrV3\nGQNTqKxWVmw95AyhzdsER3YGuljGGNNq/Jo9VlV319lV1QplOWV4htAu3ewzhNYmFjTGdGD+hMVu\nETkbUBGJEJGf4TZJBavw0BDOGZDkDKHt2h+6pFpTlDGmQ/MnLG4AbgZOA/YAI93HQS0jPZn9BaVk\n5RY7tYvsZVBRGuhiGWNMq/BnivKDqjpTVbupaoqqzlLVQ21RuPZskncIrbsgUsUx2LkiwKUyxpjW\nEVQr3rUkzxDaJVm5kHoOhEVZU5QxpsOysDgJ3iG0VWGQ+i2b+sMY02FZWJyEzHTPENqDTlPU4W1w\naFugi2WMMS2uuVlnJ4nIcPf+VSLyZxG5XUQi26Z47dvoPl2Ijwxz+y0udHZu/SCwhTLGmFbQaFiI\nyGPAb4CnReRF4HvAOmAU8EzbFK99c2ahdYfQdkmF5EGw4k9wODvQRTPGmBbVVM0iU1W/hXO19hTg\nClV9HJgNDG+Lwp0KMtKTOVBYSlZOEXznKWdU1IJpcHh7oItmjDEtpqmwKAVQ1VJgp7u2BaqqQEUb\nlO2U4BlCu3hTHvQYDj9YCBVHLTCMMR1KU2GRIiJ3iMhPfe57Hie3UfnavVpDaMEJjNn/gooSJzCs\nw9sY0wE0FRZP4axbEedz3/P46dYv2qkjIz2F1TuPUFTqVri6D4PZCy0wjDEdRqPrWXjWsTDNy0hP\n5vGl21ix9SAXD+3h7Ow+zKlhPH+ZExhz3obE/oEtqDHGnKDmhs5OEZFlInJQRPJEZKmIXNJWhTtV\n1BpC66v7UKcPo6rMahjGmFNaU0NnrwN+DdwH9AP6A/cD94nI9W1SulNEeGgI56S5Q2hVaz/ZfahT\nw6gqgwVTLTCMMaekpmoWtwOTVfUjVS10t49whtHe3jbFO3V4htBuOlBU/8luQ9zAKHcC4+DWti+g\nMcachKbCQlT1cN2dNuNswyYNTAFgsWdUVF3dhsDst6GqwgLDGHPKaSosCkVkRN2d7r4G/nwObt07\nRTHq9M48+uFW3t+Q0/BB3QY7NYzqSgsMY8wppamw+CmwUETuE5FLRWSaiNwPvAXc0TbFO7U8MWs0\nA7vFcf0Lq3hq2fb6/RfgBMact30CY0vbF9QYY45To2GhqsuB8e4xc4C57v2z3OeaJSIXi0iWiGwV\nkV80csxVIrJBRNaLyP/z2T9bRLa42+zj+E4Bk5IQxSvXT2DK0O7876KN/PKNb6ioqm7gwDOcwNAq\nZ5SUBYYxpp2TBv/69T1AJAoYACiwzZ3+o/k3FgkFNgMX4izHuhK4RlU3+ByTBvwdOE9Vj4hIiqrm\nikhXYBUwxv3c1cBoVT3S2OeNGTNGV61a5U/RWl11tfKH97N4bPE2Jg5I5C/fG02nmPD6B+Zuguem\ngYQ4/RnJA9u+sMaYoCYiq1V1THPHNTV0NkxEfg/sBp4DXgR2i8jvRaSB33z1jAO2qup2VS0HXgGm\n1znmOuAxTwioqqd3+CLgfVU97D73PnCxH5/ZLoSECD+/aBDzrhzBF9mH+fZfV7Dz0NH6B6YMckJC\n1QmNvM1tX1hjjPFDU30WDwFdgX6qOlpVR+Fca9EZmOfHe5+GEzQee9x9vgYCA0VkhYh8JiIXH8dr\n270Zo3vx4g/Hc/hoOZc/toIvsusNLnMCY44bGAumWmAYY9qlpsJiGnCdqnpHPqlqIXAj4M9V3NLA\nvrptXmFAGpABXIOzdkZnP1+LiFwvIqtEZFVeXl4DLwm88f0SefOmiXSJiWDm05/xj9V76h+UnO4E\nBriBkdW2hTTGmGY0FRaqDXRouFOVN93R4dgD9PZ53AvY18Axb6lqhapmA1k44eHPa1HVJ1V1jKqO\nSU5uvxPhpibF8sZNExmb2pWfvvYV897Norq6zimsFRjTnP4MY4xpJ5oKiw0i8oO6O0VkFuDPb7KV\nQJqI9BWRCOC7wMI6x7wJZLrvm4TTLLUdeBeYLCJdRKQLMNndd8rqFBPOc3PH8d2xvfnz4q385OU1\nlFZU1T4oOR3m/BtEnD4MCwxjTDvR6KyzwM3AP0VkLs5oJAXGAtHAt5t7Y1WtFJEf4/ySDwWeUdX1\nIvIAsEpVF1ITChuAKuDnnivEReTXOIED8EBDV5OfasJDQ3jwO8PolxzLg+9sYk9+CU/9YDQp8VE1\nByUPdDq9n5vmbLP/5Qy1NcaYAPJn6Ox5wBCcfoT1qvphWxTseLWnobP+eHf9AW57ZS1dYyN4evYY\nzuiRUPuAg1uc5qjqSqd5ygLDGNMKTnrorIc7keCjqjq/vQbFqeiiId157YYJVFZXM+Ovn7B4U505\npZLSnJAICXP7MDYGpqDGGIMfYWFaz9DTOvHWzeeQmhTLD59bybMrsmtPEVI3MHI2NP5mxhjTiiws\nAqx7pyheu2ECF5zRjfv/tYF73lpPpe8UIUlpTqd3aDg8d6kFhjEmICws2oGYiDAenzWaH03qxwuf\n7eTaBSsp9KznDZA0wCcwpkHO+sAV1hgTlCws2omQEOGuKWfwuyuG8em2Q1zxl0/YffhYzQGJ/d3A\niHRqGAfWBa6wxpigY2HRzlw99nSe/+E4covKuPyxFaze6TNiOLG/04cRGgnPX2aBYYxpMxYW7dDZ\n/ZN446aziY8K45qnPuettXtrnvQNDKthGGPaiIVFO9UvOY43bprIyN6dufWVtTz8/uaakVKewAiP\ndgPjm8AW1hjT4VlYtGNdYiN48YfjmTG6F3/6cAu3vrK2ZoqQWoFxmQWGMaZVWVi0cxFhITw0Yzj/\nfXE6C7/ax/ee+oyDxWXOk137uYER49Qw9n8d2MIaYzosC4tTgIhwU8YA/jrzTDbsL+Tyx1awOced\nOd4bGLFOp/f+rwJbWGNMh2RhcQqZMqwHr14/gbLKaq74yycs3eyu4dG1r09gTLfAMMa0OAuLU8yI\n3p156+aJ9Ooaw9wFK3nh0x3OE57AiIhzpgb5162w+V2o8GvJdGOMaVKzs86eKk61WWdPVnFZJbe+\nvIYPN+Uy5+xU7p42mNAQgSM74YN7Ycv7UF7s1DYGnAfpl0DaRRCbGOiiG2PaEX9nnbWwOIVVVSv/\nt2gjf1ueTWZ6Mo9+70ziIt0lSirLYMfHsGkRZL0DRftAQqD3WTDoEic8EvsH9gsYYwLOwiKIvPjZ\nTu5duJ60lDj+Nmcsp3WOrn2AKuxf6wbHIshxL+RLSneDYyqcNhpCrFXSmGBjYRFkPt6Sx00vfUlk\nWChPzx7DyN6dGz/4yE6ntpH1b9ixArQKYlMg/WInOPpNcq7fMMZ0eBYWQWhLThFzn1tJTkEZ3x51\nGteek8qg7glNv6jkCGz5wAmOLR9AeZFz3UZ/t59j4EUQm9Q2X8AY0+YsLILUoeIy/vD+Zv755R5K\nK6qZOCCRa8/uy3mDUggJkaZfXFkGO5Y7TVVZ70DhXrefY7wTHOmXONOlG2M6DAuLIJd/rJyXv9jN\n85/uYH9BKamJMcw5O5UZY3rXdII3RdW5XiPL7efwTCeSNLAmOHqNgZDQVv0expjWZWFhAKioquY/\n6w7wzIps1uzKJz4yjKvH9mb22an07hrj/xvl73L7ORY5tY/qSohNhoEXO8HRLwMijuP9jDHtgoWF\nqefLXUd4dsUOFn2zH1Vl8uDuzD2nL2NTuyDSTBOVr5J82PqBExxb3oeyQgiLhv6Zbj/HxRCX3Hpf\nxBjTYiwsTKP2F5Tw/Kc7+X+f76KgpIKhpyUwd2Jfpg3vSUTYcQ6frSyHncudWsemRVC4BxDoPc4J\njkFTnXXEjTHtkoWFaVZJeRX/XLOHZ1fsYGtuMcnxkXz/rD58b/zpJMVFHv8bqjp9G1mLYNO/4YA7\nC27igJrg6DXW+jmMaUcsLIzfVJWPtxzkmRXZLMnKIyIshMtH9uTaiX05o0czQ2+bUrCnpp8j+2Oo\nroCYJLefY4rTbBUR23JfxBhz3CwszAnZmlvMgk+y+cfqvZRUVHF2/0TmTvRz6G1TSgvcfo53YPN7\nUFYAYVFOx3j6JU54xKW01NcwxvjJwsKclPxj5byycjfPfXKCQ2+bUlUBO1fU9HMU7ALEGYrr7ecY\nCMfT6W6MOSEWFqZFVFRV8+76AzyzPJsvT2bobWNUIWd9TT/H/rXO/q79aq7n6D0eQk8yoIwxDbKw\nMC1ujc/Q2+qTGXrblIK9sNmtcWQvc/o5orv69HOcB5FxLfNZxhgLC9N6DhSU8sJnO3jp813kH6sZ\nejt1eA8iw1pwpFNpIWz70O3neBdK8yE00pno0NPPEd+95T7PmCBkYWFaXUl5FW+u3cszy7PZ4g69\nnTW+DzPPOsGht02pqoBdn7r9HP+G/J3O/tNG1/RzJA+yfg5jjpOFhWkzqsryrQd5Znk2i92ht9NH\n9GTq8B4M7pFAcnxkyzVTOR8IuRudmXI3LYJ9Xzr7u6Q6U6ynT4HTJ1g/hzF+sLAwAbEtr5gFK3bw\n+uo9lFRUAdA1NoIzesQzqHsCg7rHc0aPBAakxBEV3kJNVoX7ffo5lkJVOUR1dqZXT5vsrAjYuQ9E\nd7GahzF1WFiYgCouq2T93gI27i9k04EiNu4vJCuniNKKagBCQ4T+ybFOgPRwAuSM7gl0SzjJWkhZ\nsU8/x3+c9To8IuKg8+m1t0693ft9IKarhYkJOhYWpt2pqlZ2HDrKpv1FbDpQyMb9hWzcX8Te/BLv\nMV1iwusFSFq3E6yFVFVC7ganfyN/tzNzrnfb6UyA6Cs8tk6Y9Pa53wdiEi1MTIdjYWFOGQUlFWQd\nqB0gWQeKvM1YoSFC36RYbxPWGW6QdE+IOrlaSEm+ExwFDQRJ/i7nqnNf4TF1aiM+QdL5dGdFQQsT\nc4qxsDCntKpqZdfhY04z1v5CNri1kT1HamohnaLD6wXIwG7xLdcXUlrQcI3Ec780v/bxYdF1aiM+\nYZLQ01nn3DrdTTtjYWE6pMJStxbiEyBZB4o4Vu7UQkIEUpNi3SYsp1M9rVscvbrEEHoyc1s1pLSw\n8VpJ/q7a/SUAiDP/VXx3iO/hc9uj9uOYRAg5zqnijTlBFhYmaFT71EI2up3pmw4UsvtwTS0kMiyE\nfslxDEiJIy2l5rZPYuzxr+Hhr7Iit2ayEwr3QdEBKNrv3rr3jx2s/7qQMIjr7oZH3WDxubXRXaYF\ntIuwEJGLgT8BocDTqvrbOs/PAR4C9rq7/qyqT7vPVQHuws/sUtXLmvosCwtTV1FpBZtzitiaW8zW\n3GK2uLe+TVlhIUKfxBg3POIZ4AZJ/+Q4oiPaYN2NynIozqkTJA3c1m3yAudq9qbCxHMbGW+hYhoV\n8LAQkVBgM3AhsAdYCVyjqht8jpkDjFHVHzfw+mJV9XsSIAsL469j5ZVszzvKllwnSLbkFLM1r5id\nh45RVe38PIhAry7RDEiOI61bPAOS4xjQzQmShKjwti90RUntGkmDtwegvKj+a8NjndBI6AnJ6c6V\n7p7b2GQLkiDnb1i0Zm/bOGCrqm53C/QKMB3Y0OSrjGllMRFhDD2tE0NP61Rrf1llFTsPHXPCI7fY\nGyYrth2ivLLae1y3hEhvTaS/T7NWYmxEy16p7is8Grr2dbamlBVBUU4DQbLPaRL7+u+1hwxHd3FC\nI2lg7RBJ6GkhYmppzbA4Ddjt83gPML6B464QkXNxaiG3q6rnNVEisgqoBH6rqm/WfaGIXA9cD3D6\n6ae3ZNlNEIoMC2Vgt3gGdouvtb+qWtl9+Ji3GWtLbhHbcot5bdVujrod6+BcI+I0Y8XX6hvp0ekk\nh/ge15eId7akAQ0/r+oESN4myNvs3mbBxoXw5XM1x0XE16mFuFun063zPUi1ZjPUlcBFqvpf7uPv\nA+NU9Sc+xyQCxapaJiI3AFep6nnucz1VdZ+I9AM+As5X1W2NfZ41Q5m2pqrsLyj1hoizFbElt5j8\nYxXe4+IiwzitczTJ8ZGkxEeSnBBJclwkKQlR7q2zPy4yrO1Cpf6XgaMH3fBwA+RglnNbnFNzXFg0\nJNephSSlO/Ny2bDgU1J7aIbaA/T2edwL2Od7gKoe8nn4FPA7n+f2ubfbRWQJMApoNCyMaWsiQs/O\n0fTsHM2kgcne/arKoaPl3k71bbnF7MsvIbeojOyDR8krKqO8qrre+0WFh5ASH1UTKu6tZ5/ncWJc\nZMsPAxaBuGRn6/ut2s8dOwwHfWoheZtgx3L4+tWaY0IjIDGtfp9I134QFtGyZTUB0Zo1izCcpqXz\ncUY7rQS+p6rrfY7poar73fvfBu5U1bNEpAtwzK1xJAGfAtN9O8frspqFOVWoKoUlleQWlZJbVEZe\nURm5RaXubRm5hWXkFZeRW1hKYWllvdeHCCTGRXprJbVvo3yCJqp1R3SVFsLBLbVrI3mbaqaPB2cY\ncNd+PiEyyOlsj+rkTPYY1cmZs8uatgIm4DULVa0UkR8D7+IMnX1GVdeLyAPAKlVdCNwiIpfh9Esc\nBua4Lz8DeEJEqoEQnD4L6xg3HYKI0CkmnE4x4aTV6R+pq7SiirwiT3iUkecTKp7bjfsLOVhc7h3J\n5SsuMswbHt0SokhNiqVfUiypSbH0TYqlU/RJjOyKSoBeo53NV/lRJ0R8ayO5G511SLR+jQoJgcgE\nJzii3QDxbp1rQsV38z0uPMY649uAXZRnTAdQXa0cPlZeJ0hKvbWUvMIy9heWsPdICb6ZkhgbQV+f\n8PAESWpibMvXSirL4NA2OJrrTKXiu5Xk199X6u6rONb0+4aE1a6pNBYq3uc7O/vjUpyQCvKgCfh1\nFm3NwsKY5pVVVrH78DG25x1lx6GjZB886r2fU1hW69ienaK8IeK79e4aQ3hoGzYbVZY7w309AdJU\nsNQLoHxnfZPGhEXXTMESlwJx3Zyr5+vu68DzegW8GcoY0/5EhoW6Q3vrN38Vl1Wy46ATIJ7b7QeP\n8vbX+ykoqRndFRoi9O4S7YZHHH2TY+mbGEvf5Fh6JEQR0tKd72EREJbkzOp7IipK64fKscPOKK/i\nHCjOheIDTtNZ9scNXy2POHN21QoVd4vvVvtxB71i3sLCGAM4/RsNXawIcORoOdt9QsQTJJ9tP+yd\nSh6cObhSE91aiE+I9E2Kbd2LFpsSHuVs8d38O76yzA2QnJqpWHwfF+c4wVKc03CtJSzaJ0BS3JpK\nt/r7YpNPqdqKNUMZY06YqpJTWMb2g8XsOHiM7IPF3jDZdfgYFVU1v1/iI8PomxxLry7RdEuIoken\nKLp3iqa7ez8lIZLIsDaYj6ulqDozC/vWUIoO1K6tePY1W1vpVnPrDZbuNQETEdtqX8P6LIwxAVVZ\nVc3e/JJ6NZJ9+SXsLyj1TivvKzE2gu6dnPCoGyjdOzlbXOSp89e4V93aSnGOMy1L8QH31merrj9c\nmoj4mgDx9qd082kW637CMxFbn4UxJqDCQkPokxhLn8RYSK/9nKpSVFZJTkEp+wtKOVBQyoFCz/0S\n9uaXsnrnEY74XAnvER8Z5g0OT62kmxsw3ROi6d4pii4x4YG7Gr4hYZHuwli9mz6uuhpKDvs0f/k2\ng7nBsn8tbM6BiqP1Xx8S3nDNxLfmEu9pAju+YdMWFsaYNiciJESFkxDV9LUmpRVV5LghklPoEywF\npewvLGVzTh55RWXUvcQkIizEWxvp4YZK7RpLNElxEYS15aguf4SEOB35sUnQbUjTx5YV+4TJgdrN\nYEUH4Eg27P4Mjh1q4MU+TWB+srAwxrRbUeGhNbWTRlRWVZNXXOYEiqemUlgTKmt25XOgoLTeFCsh\nAsnxkXRPcAKku2/TV4JTW+meEEVse232ioxztsT+TR9XWe5c2+Jt7jpQuwnMT9ZnYYzp8FSVw0fL\nvSHiqal4mr889xuaXiU+MswbHN3rBIlzP5Kk2MiWHzLcRqzPwhhjXCJCYpwzCeOQnvWHBnscK6+s\nEyBl3iDZX1jK1i0HySsuqze1SliI0C0him4Jkd4aSvc6tZVuCVFEhZ9Co73qsLAwxhhXTEQY/ZLj\n6Jfc+CKdVdXKweKyWqHiaQI7UFjKpgNFLM3Kq7XWiUfnmPCaZq+EKDrHhhMfGUZcZBhxUeHERYYR\nH+V5HOY8FxVGdHhowDvsLSyMMeY4hHprEVGMaOK4otIKb+1kf0GJc9+ntrJhfyEFJRW1VmFsTIjg\nBkm4N0hqBYrPPidswus8dp6PjQg74entLSyMMaYVxEeFEx8V3uDUKr7KKqs4WlZFcWklRWUVFJdW\nUlzmbEWe+6W+jysoLqsk/1g5u48c8z7X0HUrDYmNCPUJG/+Hz1pYGGNMAEWGhRIZFkrX2JNbJKqq\nWr0hU+yGSsNh4/O4rIELABthYWGMMR1AaIjQKTr8uNcoeeGH/h3Xzq5IMcYY0x5ZWBhjjGmWhYUx\nxphmWVgYY4xploWFMcaYZllYGGOMaZaFhTHGmGZZWBhjjGlWh5miXETygJ2BLkcbSwIOBroQ7YSd\nixp2LmrYuajR2Lnoo6rJzb24w4RFMBKRVf7MQx8M7FzUsHNRw85FjZM9F9YMZYwxplkWFsYYY5pl\nYXFqezLQBWhH7FzUsHNRw85FjZM6F9ZnYYwxpllWszDGGNMsC4t2SkSeEZFcEVnns6+riLwvIlvc\n2y7ufhGR+SKyVUS+FpEzA1fylicivUVksYhsFJH1InKruz/ozoeIRInIFyLylXsu7nf39xWRz91z\n8aqIRLj7I93HW93nUwNZ/tYgIqEiskZE3nYfB+W5EJEdIvKNiKwVkVXuvhb7GbGwaL8WABfX2fcL\n4ENVTQM+dB8DTAHS3O164K9tVMa2Ugn8VFXPAM4CbhaRwQTn+SgDzlPVEcBI4GIROQv4HfCwey6O\nAJ4lbX4IHFHVAcDD7nEdza3ARp/HwXwuMlV1pM8Q2Zb7GVFV29rpBqQC63weZwE93Ps9gCz3/hPA\nNQ0d1xE34C3gwmA/H0AM8CUwHudiqzB3/wTgXff+u8AE936Ye5wEuuwteA56ub8EzwPeBiSIz8UO\nIKnOvhb7GbGaxamlm6ruB3BvU9z9pwG7fY7b4+7rcNymg1HA5wTp+XCbXdYCucD7wDYgX1U9Cyr7\nfl/vuXCfLwAS27bEreoR4L+BavdxIsF7LhR4T0RWi8j17r4W+xmxNbg7BmlgX4cb5iYiccA/gNtU\ntVCkoa/tHNrAvg5zPlS1ChgpIp2BN4AzGjrMve2w50JEpgG5qrpaRDI8uxs4tMOfC9dEVd0nIinA\n+yKyqYljj/tcWM3i1JIjIj0A3Ntcd/8eoLfPcb2AfW1ctlYlIuE4QfGSqv7T3R205wNAVfOBJTj9\nOJ1FxPPHn+/39Z4L9/lOwOG2LWmrmQhcJiI7gFdwmqIeITjPBaq6z73NxfkjYhwt+DNiYXFqWQjM\ndu/Pxmm79+z/gTvC4SygwFP17AjEqUL8Ddioqn/0eSrozoeIJLs1CkQkGrgAp3N3MTDDPazuufCc\noxnAR+o2Up/qVPUuVe2lqqnAd3G+20yC8FyISKyIxHvuA5OBdbTkz0igO2Vsa7Sz6mVgP1CB81fA\nD3HaVz8Etri3Xd1jBXgMp+36G2BMoMvfwufiHJwq8tfAWne7JBjPBzAcWOOei3XAPe7+fsAXwFbg\nNSDS3R/lPt7qPt8v0N+hlc5LBvB2sJ4L9zt/5W7rgf9x97fYz4hdwW2MMaZZ1gxljDGmWRYWxhhj\nmmVhYYwxplkWFsYYY5plYWGMMaZZFhbGGGOaZWFhzEkQkZEiconP48tE5BdNveY43vs2EYlpifcy\n5mTZdRbGnAQRmYNzQdOPW+G9d7jvffA4XhOqztxRxrQoq1mYoCAiqe7iSU+5iwa9506X0dCx/UXk\nP+7snR+LyCB3/5Uiss5deGiZu6jOA8DV7oIzV4vIHBH5s3v8AhH5qzgLN20XkUniLGq1UUQW+Hze\nX0VkldRezOgWoCewWEQWu/uucRe3WSciv/N5fbGIPCAinwMTROS3IrLBXdRmXuucURN0An2Zum22\ntcWGszZIJTDSffx3YFYjx34IpLn3x+PMIQTOtAinufc7u7dzgD/7vNb7GGcBq1dwplaYDhQCw3D+\nSFvtUxbPFAyhOBMDDncf78BdnwAnOHYByTizRX8EXO4+p8BVnvfCWZtAfMtpm20nu1nNwgSTbFVd\n695fjRMgtbjToJ8NvOauGfEEzqIxACuABSJyHc4vdn/8S1UVJ2hyVPUbVa3Gmb/H8/lXiciXOHM+\nDQEGN/A+Y4ElqpqnzloMLwHnus9V4czIC04glQJPi8h3gGN+ltOYJtl6FiaYlPncrwIaaoYKwVk8\nZ2TdJ1T1BhEZD0wF1opIvWOa+MzqOp9fDYSJSF/gZ8BYVT3iNk9FNfA+jS7eAZSq20+hqpUiMg44\nH2cm1h/jTN1tzEmxmoUxPlS1EMgWkSvBu7D9CPd+f1X9XFXvwVmSszdQBMSfxEcmAEeBAhHphrM2\nsofve38OTBKRJBEJBa4BltZ9M7dm1ElVFwG34azTbcxJs5qFMfXNBP4qIr8CwnH6Hb4CHhKRNJy/\n8j909+0CfuE2WT14vB+kql+JyBqcZqntOE1dHk8C74jIflXNFJG7cNZqEGCRqr5V/x2JB94SkSj3\nuNuPt0zGNMSGzhpjjFwJdEgAAABNSURBVGmWNUMZY4xpljVDmaAlIo/hrOPs60+q+mwgymNMe2bN\nUMYYY5plzVDGGGOaZWFhjDGmWRYWxhhjmmVhYYwxplkWFsYYY5r1/wF4fu3Or/gEugAAAABJRU5E\nrkJggg==\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1ecad4b4400>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "#define a function to explore the optimal n_estimators setting\n",
    "def mod_optimization (min_estimator, max_estimator, step, train_X, train_y):\n",
    "\n",
    "    RANDOM_STATE = 123\n",
    "\n",
    "    #define a list of classifiers with different splitting options (max_features)\n",
    "    ensemble_classifiers = [(\"RandomForestClassifier, max_feature='sqrt'\", \n",
    "                             RandomForestClassifier(warm_start=True, oob_score=True, max_features='sqrt', \n",
    "                                                    random_state=RANDOM_STATE, n_jobs=-1)),\n",
    "                           (\"RandomForestClassifier, max_feature='log2\", \n",
    "                           RandomForestClassifier(warm_start=True, oob_score=True, max_features='log2',\n",
    "                                                 random_state=RANDOM_STATE, n_jobs=-1))]\n",
    "\n",
    "    error_rate = OrderedDict((label, []) for label, _ in ensemble_classifiers)\n",
    "\n",
    "    #model evaluation\n",
    "    for label, classifier in ensemble_classifiers:\n",
    "        for i in range(min_estimator, max_estimator + 1, step):\n",
    "            classifier.set_params(n_estimators = i )\n",
    "            classifier.fit(train_X, train_y)\n",
    "            print(label,\" fit for \", i, \"estimators \")\n",
    "\n",
    "            oob_error = 1 - classifier.oob_score_\n",
    "            error_rate[label].append((i, oob_error))\n",
    "\n",
    "    #plot the error vs the number of estimators\n",
    "    for label, error in error_rate.items():\n",
    "        x, y = zip(*error)\n",
    "        plt.plot(x, y, label=label)\n",
    "\n",
    "    plt.xlim(min_estimator, max_estimator)\n",
    "    plt.xlabel(\"n_estimators\")\n",
    "    plt.ylabel(\"OOB error rate\")\n",
    "    plt.legend(loc=\"upper right\")\n",
    "    plt.show()\n",
    "    \n",
    "\n",
    "mod_optimization(10, 510, 50, train_X, train_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Fit the model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Based on the model optimization above, n_estimators = 500 was selected for the random forest classifier in thr formal analysis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model trained!\n",
      "Training accuracy: 1.0\n",
      "Elapsed time:  5.4784034530321755\n"
     ]
    }
   ],
   "source": [
    "#train the random forest classifier\n",
    "start = time()\n",
    "rf_cifar = RandomForestClassifier(n_jobs=-1, n_estimators=500)\n",
    "#The first 35000 labels correpond to the images used as training set\n",
    "rf_cifar.fit(train_X, train_y)\n",
    "end = time()\n",
    "print (\"Model trained!\")\n",
    "print(\"Training accuracy:\", rf_cifar.score(train_X, train_y))\n",
    "print(\"Elapsed time: \", (end-start)/60.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Precision:\n",
      "airplane : 0.567398119122\n",
      "automobile : 0.521781219748\n",
      "bird : 0.388607594937\n",
      "cat : 0.404109589041\n",
      "deer : 0.400801603206\n",
      "dog : 0.425349087003\n",
      "frog : 0.450826121164\n",
      "horse : 0.520134228188\n",
      "ship : 0.585454545455\n",
      "truck : 0.47762345679\n",
      "\n",
      "Recall:\n",
      "airplane : 0.53550295858\n",
      "automobile : 0.531558185404\n",
      "bird : 0.322478991597\n",
      "cat : 0.290354330709\n",
      "deer : 0.401203610832\n",
      "dog : 0.386341463415\n",
      "frog : 0.584693877551\n",
      "horse : 0.475946775844\n",
      "ship : 0.642073778664\n",
      "truck : 0.605675146771\n",
      "\n",
      "Overall accuracy:  0.4781\n"
     ]
    }
   ],
   "source": [
    "#evaluate the performance of the classifier\n",
    "\n",
    "pred_y = rf_cifar.predict(test_X)\n",
    "result = precision_recall_fscore_support(test_y, pred_y)\n",
    "classes = [class_le.inverse_transform(i) for i in rf_cifar.classes_]\n",
    "print (\"Precision:\")\n",
    "for i in range(len(classes)):\n",
    "    print(classes[i],':', result[0][i])\n",
    "print (\"\\nRecall:\")\n",
    "for i in range(len(classes)):\n",
    "    print(classes[i],':', result[1][i])\n",
    "\n",
    "print (\"\\nOverall accuracy: \", accuracy_score(test_y, pred_y))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Classification of the handwritten digits in the MNIST database by Random Forest Classifier"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Data preprocessing\n",
    "\n",
    "Images have been preprocessed and readily available in the tensorflow package. For this analysis, we just import data as numpy arrays."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "WARNING:tensorflow:From C:\\Users\\isaac\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\tensorflow\\contrib\\learn\\python\\learn\\datasets\\base.py:198: retry (from tensorflow.contrib.learn.python.learn.datasets.base) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Use the retry module or similar alternatives.\n",
      "WARNING:tensorflow:From <ipython-input-27-ce30900f3751>:3: read_data_sets (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Please use alternatives such as official/mnist/dataset.py from tensorflow/models.\n",
      "WARNING:tensorflow:From C:\\Users\\isaac\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\tensorflow\\contrib\\learn\\python\\learn\\datasets\\mnist.py:260: maybe_download (from tensorflow.contrib.learn.python.learn.datasets.base) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Please write your own downloading logic.\n",
      "WARNING:tensorflow:From C:\\Users\\isaac\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\tensorflow\\contrib\\learn\\python\\learn\\datasets\\mnist.py:262: extract_images (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Please use tf.data to implement this functionality.\n",
      "Extracting MNIST_data\\train-images-idx3-ubyte.gz\n",
      "WARNING:tensorflow:From C:\\Users\\isaac\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\tensorflow\\contrib\\learn\\python\\learn\\datasets\\mnist.py:267: extract_labels (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Please use tf.data to implement this functionality.\n",
      "Extracting MNIST_data\\train-labels-idx1-ubyte.gz\n",
      "Extracting MNIST_data\\t10k-images-idx3-ubyte.gz\n",
      "Extracting MNIST_data\\t10k-labels-idx1-ubyte.gz\n",
      "WARNING:tensorflow:From C:\\Users\\isaac\\AppData\\Local\\Continuum\\anaconda3\\lib\\site-packages\\tensorflow\\contrib\\learn\\python\\learn\\datasets\\mnist.py:290: DataSet.__init__ (from tensorflow.contrib.learn.python.learn.datasets.mnist) is deprecated and will be removed in a future version.\n",
      "Instructions for updating:\n",
      "Please use alternatives such as official/mnist/dataset.py from tensorflow/models.\n"
     ]
    }
   ],
   "source": [
    "import tensorflow as tf\n",
    "from tensorflow.examples.tutorials.mnist import input_data\n",
    "mnist = input_data.read_data_sets('MNIST_data', one_hot=False)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data extracted from the tensorflow package are typically in tuples. Each element, which is an array, in the tuples represents an image. Because the images are gray, the dimension of each image is 28x28. There 55000 images in the training set and 10000 images in the testing set."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Optimize the model\n",
    "\n",
    "Similar optimization practice was conducted by using the mod_optimization function defined in the previous steps. The primary purpose was to identify the optimal n_estimators."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RandomForestClassifier, max_feature='sqrt'  fit for  10 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  60 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  110 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  160 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  210 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  260 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  310 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  360 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  410 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  460 estimators \n",
      "RandomForestClassifier, max_feature='sqrt'  fit for  510 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  10 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  60 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  110 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  160 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  210 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  260 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  310 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  360 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  410 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  460 estimators \n",
      "RandomForestClassifier, max_feature='log2  fit for  510 estimators \n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYsAAAELCAYAAAAoUKpTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMS4wLCBo\ndHRwOi8vbWF0cGxvdGxpYi5vcmcvpW3flQAAIABJREFUeJzt3Xt8FPW5+PHPs5uEJFwCCQGCICQI\nKLcEBRKsFahXEMRT8XbgVGqPnqqttdYepT1apP1Vq1Q9WI6X2kqrPWrFo1KL1VZBihcEFOUmiIDc\nIQQIBBKS3X1+f8zsZhM22U3IJmHzvF+vee3Md74z892B5MnMM/P9iqpijDHG1MfT0g0wxhjT+lmw\nMMYYE5UFC2OMMVFZsDDGGBOVBQtjjDFRWbAwxhgTlQULY4wxUcU1WIjIpSKyQUQ2icjdEdafLyIf\ni4hPRKbUWne6iLwlIutFZJ2I9I1nW40xxtQtbsFCRLzAXGA8MAi4TkQG1aq2DZgO/G+EXfwReEhV\nzwJGAfvi1VZjjDH1S4rjvkcBm1R1M4CIvABMBtYFK6jqVnddIHxDN6gkqerf3Xpl0Q7WtWtX7du3\nb1O13Rhj2oSVK1fuV9XsaPXiGSxOA7aHLe8ACmPcdgBwSET+D8gF/gHcrar+ujbo27cvK1asaGxb\njTGmTRKRr2KpF8+chUQoi7UjqiTg68CdwEggD+d2Vc0DiNwkIitEZEVxcXFj22mMMSaKeAaLHUDv\nsOVewK4GbPuJqm5WVR/wKnB27Uqq+pSqjlDVEdnZUa+ijDHGNFI8g8VyoL+I5IpICnAtsKAB23YR\nkWAE+AZhuQ5jjDHNK245C1X1icj3gDcBL/B7VV0rIrOAFaq6QERGAq8AXYBJInKfqg5WVb+I3Am8\nLSICrAR+G6+2thVVVVXs2LGDioqKlm6KMaaZpaam0qtXL5KTkxu1vSTKeBYjRoxQS3DXb8uWLXTs\n2JGsrCycGGyMaQtUlZKSEo4cOUJubm6NdSKyUlVHRNuHvcHdhlRUVFigMKYNEhGysrJO6q6CBYs2\nxgKFMW3Tyf7sJ06wOFoMFYdbuhXGGJOQEidYlO6APZ+1dCuMMSYhJU6wANhtwaK183q9FBQUMGTI\nECZNmsShQ4eaZL9bt25lyJAhTbKv6dOnk5ubS0FBAQUFBcyZM6dJ9hvJ4sWLef/992uU/fGPf2TI\nkCEMHjyYQYMGMXv27FC75s+f3yTH3bVrF1OmVPfded111zFs2DAeeeQR7r33Xv7xj380yXGa0+ef\nf05BQQHDhw/nyy+/bPD2jz76KMeOHYtDy+rWVF0UrVq1ioULF4aW582bx8yZM5tk30Hx7O6jeXmT\nYfenLd0KE0VaWhqrVq0C4Prrr2fu3Ln89Kc/beFWneihhx6q8cs0Vn6/H6/XG3P9xYsX06FDB849\n91wA3njjDR599FHeeustevbsSUVFBc8++2yD2xFNz549Q4Fnz549vP/++3z1VUy9PpzA5/ORlNTy\nv0peffVVJk+ezH333deo7R999FGmTZtGenp6zNu0hu/u8/lYtWoVK1asYMKECXE7Tsv/CzeV5HQL\nFg1w31/Wsm5X0+Z4BvXsxM8mDY65/ujRo/nsM+dqsKysjMmTJ3Pw4EGqqqr4xS9+weTJk9m6dSvj\nx4/nvPPO4/333+e0007jtddeIy0tjZUrV3LDDTeQnp7OeeedF9pvRUUFN998MytWrCApKYmHH36Y\ncePGMW/ePF599VX8fj9r1qzhRz/6EZWVlTz77LO0a9eOhQsXkpmZWWd7n3/+eX75y1+iqlx22WX8\n6le/AqBDhw7ccccdvPnmm/z6178mLS2NO+64g7KyMrp27cq8efPIyclhzpw5PPHEEyQlJTFo0CAe\neOABnnjiCbxeL8899xyPPfYY999/P7Nnz6Znz56A82z8jTfeeEJbZs2axV/+8hfKy8s599xzefLJ\nJxGRE47xwgsv8O677/KDH/wAcJKcS5YsoaSkhIkTJ7JmzRouvvhi9u3bR0FBAY899hi/+93vmDhx\nIlOmTGHlypURv8vYsWM599xzee+997j88sv50Y9+FPGcxXrOf/vb3/LUU09RWVnJGWecwbPPPkt6\nejqTJ0/myiuv5Fvf+hZPPvkkS5Ys4U9/+tMJx1m4cCGPPvooXq+XJUuWsGjRIp577jnmzJlDZWUl\nhYWF/M///A9er5ebb76Z5cuXU15ezpQpU7jvvvuYM2cOu3btYty4cXTt2pVFixbRoUMHysqcPkzn\nz5/P66+/zrx585g+fTqZmZl88sknnH322cyaNYvvf//7rF69Gp/Px8yZM5k8eXIsPwIEe544evQo\nV199NTt27MDv93PPPfdwzTXX8Le//Y3bb7+drl27cvbZZ7N582Zef/11Zs6cya5du9i6dStdu3Zl\n6dKllJeXs3TpUmbMmEFaWhodOnSIqQ0xU9WEmM4Z0FN1ZmfV40fVRLZu3brQ/MwFa/TqJ95v0mnm\ngjVR29C+fXtVVfX5fDplyhR94403VFW1qqpKS0tLVVW1uLhY+/Xrp4FAQLds2aJer1c/+eQTVVW9\n6qqr9Nlnn1VV1aFDh+rixYtVVfXOO+/UwYMHq6rq7Nmzdfr06aqqun79eu3du7eWl5frM888o/36\n9dPDhw/rvn37tFOnTvr444+rqurtt9+ujzzyiKqqXn/99dq3b1/Nz8/X/Px8/eyzz3Tnzp3au3dv\n3bdvn1ZVVem4ceP0lVdeUVVVQF988UVVVa2srNTRo0frvn37VFX1hRde0G9/+9uqqpqTk6MVFRWq\nqnrw4EFVVf3Zz36mDz30UOj8dOnSRQ8dOhTx3F1//fX60ksvqapqSUlJqHzatGm6YMGCOo8xceJE\nXbp0qaqqHjlyRKuqqnTLli2h8xU+H36c+r7LmDFj9Oabb47YznCxnvP9+/eHtvnpT3+qc+bMUVXV\nPXv2aL9+/XTJkiXav3//Gt+7tvBzuW7dOp04caJWVlaqqurNN9+sf/jDH2qcO5/Pp2PGjNFPP/1U\nVVX79OmjxcXFof0F/6+qqr700kt6/fXXh87PZZddpj6fT1VVZ8yYEfo/efDgQe3fv7+WlZXp559/\nHvo/VHsK/tsEzZ8/X//93/89tHzo0CEtLy/XXr166caNGzUQCOhVV12ll112Wei7nn322Xrs2LHQ\neb711lvr+6cInZfacF6Sjvo7NrGuLLQM9q2DXlHfL2nzGnIF0JTKy8spKChg69atnHPOOVx00UWA\n80fLT37yE5YsWYLH42Hnzp3s3bsXIJQ/ADjnnHPYunUrpaWlHDp0iDFjxgDwb//2b7zxxhsALF26\nlO9///sAnHnmmfTp04eNGzcCMG7cODp27EjHjh3JyMhg0qRJAAwdOjR0lQMn3oZ67bXXGDt2bOgv\nwalTp7JkyRKuuOIKvF4vV155JQAbNmxgzZo1oe/l9/vJyckBYNiwYUydOpUrrriCK6644qTO46JF\ni3jwwQc5duwYBw4cYPDgwUyaNCniMb72ta9xxx13MHXqVL75zW/Sq1evmI5R33cBuOaaa2LaTyzn\nfM2aNfzXf/0Xhw4doqysjEsuuQSA7t27M2vWLMaNG8crr7xS75VfuLfffpuVK1cycuRIwPl/161b\nNwD+/Oc/89RTT+Hz+di9ezfr1q1j2LBhMe036KqrrgrdbnzrrbdYsGBBKLdUUVHBtm3bOOuss0K3\nXKMZOnQod955J3fddRcTJ07k61//OqtWrSI3N5f+/fsDMG3aNJ566qnQNpdffjlpaWkNavfJSKBg\n4Z603assWLRiwZxFaWkpEydOZO7cudx222386U9/ori4mJUrV5KcnEzfvn1DLxC1a9cutL3X66W8\nvBxVrfO5ca2nV4LwfXk8ntCyx+PB5/PVuV19+0xNTQ394lBVBg8ezAcffHBCvb/+9a8sWbKEBQsW\n8POf/5y1a9eeUGfw4MGsXLmSb3zjG3Uer6KigltuuYUVK1bQu3dvZs6cGTpXkY5x9913c9lll7Fw\n4UKKior4xz/+QWpqap37D//OdX0XgPbt20fdB8R2zqdPn86rr75Kfn4+8+bNY/HixaFtVq9eTVZW\nFrt2xdoPqdP266+/nvvvv79G+ZYtW5g9ezbLly+nS5cuTJ8+vc4X1cL/f9WuE/7dVZWXX36ZgQMH\n1qizYcOGOgPq4sWL6dy5c2h5wIABrFy5koULFzJjxgwuvvhiLr/88nrfjYj1/DeVxHkaypsCaV0s\nb3GKyMjIYM6cOcyePZuqqipKS0vp1q0bycnJLFq0KGqytXPnzmRkZLB06VKAGvexzz///NDyxo0b\n2bZt2wk/yA1VWFjIu+++y/79+/H7/Tz//POhq5pwAwcOpLi4OPQLtqqqirVr1xIIBNi+fTvjxo3j\nwQcfDP0F3bFjR44cORLafsaMGfznf/4ne/bsAeD48eMnPI0V/MXVtWtXysrKQonquo7x5ZdfMnTo\nUO666y5GjBjB559/HtN3ruu7RPLKK68wY8aMmPYbyZEjR8jJyaGqqqrGv+VHH33EG2+8wSeffMLs\n2bPZsmVLTPu74IILmD9/Pvv2OQNsHjhwgK+++orDhw/Tvn17MjIy2Lt3b+hqFDjh36J79+6sX7+e\nQCDAK6+8UuexLrnkEh577LHQHxSffPIJ4Jy/VatWRZzCAwU4T6elp6czbdo07rzzTj7++GPOPPNM\ntmzZEnqy6/nnn6+zDbXbHg+Jc2UBkJNvj8+eQoYPH05+fj4vvPACU6dOZdKkSYwYMYKCggLOPPPM\nqNs/88wzoQR38LYFwC233MJ3v/tdhg4dSlJSEvPmzavx121j5OTkcP/99zNu3DhUlQkTJkRMYqak\npDB//nxuu+02SktL8fl83H777QwYMIBp06ZRWlqKqvLDH/6Qzp07M2nSJKZMmcJrr73GY489xoQJ\nE9i7dy8XXnhh6OrphhtuqHGMzp07c+ONNzJ06FD69u0butXi9/sjHuOee+5h0aJFeL1eBg0axPjx\n49m9e3fU71zXdxk8+MRbmF9++SWdOnVq5NmFn//85xQWFtKnTx+GDh3KkSNHOH78ODfeeCPPPPMM\nPXv25Ne//jU33HAD77zzTtS3kQcNGsQvfvELLr74YgKBAMnJycydO5eioiKGDx/O4MGDycvL42tf\n+1pom5tuuonx48eTk5PDokWLeOCBB5g4cSK9e/dmyJAhoWR3bffccw+33347w4YNQ1Xp27cvr7/+\neoO+/+rVq/nxj3+Mx+MhOTmZxx9/nNTUVJ566ikuu+wyunbtynnnnceaNWsibj9u3DgeeOABCgoK\nmDFjRsy3CBsisToSvH8CfPg4zNgJSSkt3aRWZ/369Zx11lkt3QyTgKZNm8YjjzyCjSsTP4sXL2b2\n7NkNDkThIv0OaJsdCebkg78SimO7zDbGNI3nnnvOAkWCS6zbUD3ync/dn0JOw55uMMacGm699Vbe\ne++9GmU/+MEP+Pa3v91CLWoeY8eOZezYsS12/MQKFpl5kNLB+ogyJoHNnTu3pZvQJiXWbSiPB3oM\nsyeijDGmiSVWsAAnb7FnNQT8Ld0SY4xJGAkYLIZB1TEo2dTSLTHGmISRgMEimOS2vIUxxjSVxAsW\nXQdCUqrT7YdpdWw8i5psPIumc6qPZ3EyvcROnTqVgQMHMmTIEG644QaqqqqaoHU1JV6w8CZB98GW\n5G6lgn1DrVmzhszMzFb7ZMtDDz0U6prhtttui3k7v79hubLawSJ8PIu1a9fy8ccfk5GR0aB9xiLS\neBafffYZP/zhD5k1axYXXnhhzPuqr0+t5hQcz+KTTz6hX79+Dd6+McGitXz3qVOn8vnnn7N69WrK\ny8t5+umnm/wYiRcswHkias9nkCBvp8fFG3fDM5c17fTG3Q1qwujRo9m5cyfgjGdxwQUXcPbZZzN0\n6FBee+01wLliOOuss7jxxhsZPHgwF198MeXl5QCsXLmS/Px8Ro8eXSPoVFRU8O1vf5uhQ4cyfPhw\nFi1aBDhjK1xxxRVMmjSJ3NxcfvOb3/Dwww8zfPhwioqKOHDgQL3tff755xk6dChDhgzhrrvuCpV3\n6NCBe++9l8LCQj744ANWrlzJmDFjOOecc7jkkktCXWvMmTOHQYMGMWzYMK699lq2bt3KE088wSOP\nPEJBQQH//Oc/GzSexciRIxkyZAg33XRTqF+i2scAePfdd0NXScOHD+fIkSM1rsTCx7P45z//WeMK\npq7vMnbsWH7yk58wZswY/vu//7vOcxbrOf/tb3/LyJEjyc/P58orrwz90p48eTJ//OMfAXjyySeZ\nOnVqxOMEx7N4+umnGTduHOC8KDhq1CgKCgr4j//4j1Agv/nmmxkxYgSDBw/mZz/7Wei8BcezCG4f\n/pf+/PnzmT59OuBc4d1xxx2MGzeOu+66i6NHj3LDDTcwcuRIhg8fHvq/G4tILzKqKj/+8Y8ZMmQI\nQ4cO5cUXXwScvr9uueUWBg8ezMSJE5kwYULo32nChAmICCLCqFGj2LFjR8xtiFks/Zg3dgIuBTYA\nm4C7I6w/H/gY8AFTIqzvBOwEfhPtWOecc051B+3Lf6/6s06qB7ZE7d+9LanRl/3Cu1R/P6Fpp4V3\nRW2DjWdh41nYeBYnjmcRfrz58+frhRdeqD6fT/fs2aO9e/fWXbt26UsvvaTjx49Xv9+vu3fv1s6d\nO4f+PwRVVlbq8OHDdcmSJRHPUascz0JEvMBc4CJgB7BcRBao6rqwatuA6cCddezm58C7DT54Ttib\n3F36NnjzNmH8Ay1yWBvPwsazsPEs6rd06VKuu+46vF4v3bt3Z8yYMSxfvpylS5dy1VVX4fF46NGj\nR+gKKNwtt9zC+eefz9e//vUGHzeaeL7BPQrYpKqbAUTkBWAyEAoWqrrVXReovbGInAN0B/4GNGyA\nim6DwJPkBItBsQ1vaJqHjWdh41nYeBYnjmdRu+0NKQ+67777KC4u5sknn6y3XmPFM2dxGrA9bHmH\nWxaViHiAXwM/btSRk1Mh+0xLcrdiNp6FjWcRrq2PZxHu/PPP58UXX8Tv91NcXMySJUsYNWoU5513\nHi+//DKBQIC9e/fWCKhPP/00b775Js8//zweT3x+rcfzyiLSn32xZpxvARaq6vb6+q0XkZuAmwBO\nP/30mitz8uGLt5wkd5S+703LsPEsbDyLoLY+nkW4f/mXf+GDDz4gPz8fEeHBBx+kR48eXHnllbz9\n9tsMGTKEAQMGUFhYGHpS7rvf/S59+vRh9OjRAHzzm9/k3nvvbXQbIoolsdGYCRgNvBm2PAOYUUfd\neYQluIE/4eQztgL7gcPAA/Udr0aCW1X1wyecJHfproiJnrYoUnLLmKYwderUUCLcxM+RI0dU1Xkg\nIC8vT3fv3t2g7VtlghtYDvQXkVycJ5quBf41lg1VNfR8nIhMB0aoasOeywxPcnfKqb+uMeakPPfc\ncy3dhDZh4sSJHDp0iMrKSu655x569OjRbMeOW7BQVZ+IfA94E/ACv1fVtSIyCyeSLRCRkcArQBdg\nkojcp6onXuM2RvchgDjBYuClTbJLY0zLa6vjWQA18hTNLa7jWajqQmBhrbJ7w+aXA/U+x6eq83Bu\nUzVMuw6QdYaNbVGL1vMUkTGngtb61n9rpyf5knJivsEdlJNvT0SFSU1NpaSk5KT/0xhjTi2qSklJ\nSUyPTNclsUbKqy0nH9bMh6Ml0D6rpVvT4nr16sWOHTsoLi5u6aYYY5pZampqzC9kRpLgwcJ9K3PP\np9Cv7pec2ork5GRyc3NbuhnGmFNQYt+G6uEGCxvbwhhjTkpiB4v0TOh8uuUtjDHmJCV2sABLchtj\nTBNI/GDRIx8OfAkVh1u6JcYYc8pK/GARfJN775qWbYcxxpzC2k6wsFtRxhjTaIkfLDp2hw49LFgY\nY8xJSPxgAc77FhYsjDGm0dpIsMiH4g1QVd7SLTHGmFNS2wkW6oe966LXNcYYc4K2EywAdjd88HRj\njDFtJVhk9IbUzpa3MMaYRmobwULEubqwsS2MMaZREiZY+ANa/zgNOfmwdy34q5qvUcYYkyASJlis\n232YTfvK6q6Qkw/+Sij+vPkaZYwxCSJhggXAh1sO1L3S3uQ2xphGS5hgkez18OHmkrorZPaDlA42\ntoUxxjRCwgSL9ilelm0+UHfewuOBHkPtysIYYxohcYJFuyT2lx1n8/6jdVfKyYc9qyHgb76GGWNM\nAkioYAHUfyuqxzCoOgolXzZTq4wxJjHENViIyKUiskFENonI3RHWny8iH4uIT0SmhJUXiMgHIrJW\nRD4TkWuiHatdkoduHduxbHMMSW5738IYYxokbsFCRLzAXGA8MAi4TkQG1aq2DZgO/G+t8mPAt1R1\nMHAp8KiIdI52zKK8LD7cXFJ33iJ7IHjbWbcfxhjTQPG8shgFbFLVzapaCbwATA6voKpbVfUzIFCr\nfKOqfuHO7wL2AdnRDliYl8m+I8fZWnIscgVvMnQfbEluY4xpoHgGi9OA7WHLO9yyBhGRUUAKEDXR\nUJSXBcCy+vIWwbEt6nvb2xhjTA3xDBYSoaxBv6FFJAd4Fvi2qgYirL9JRFaIyIri4mLyurana4d2\n9Se5c/KhohQObWtIU4wxpk2LZ7DYAfQOW+4F7Ip1YxHpBPwV+C9V/TBSHVV9SlVHqOqI7OxsRITC\nvEyWbannfQt7k9sYYxosnsFiOdBfRHJFJAW4FlgQy4Zu/VeAP6rqSw05aFFeFrtLK9h2oI68RbfB\nIF4LFsYY0wBxCxaq6gO+B7wJrAf+rKprRWSWiFwOICIjRWQHcBXwpIisdTe/GjgfmC4iq9ypIJbj\nFuVmAtT9CG1yKmSfacHCGGMaICmeO1fVhcDCWmX3hs0vx7k9VXu754DnGnPMM7p1IKt9Ch9uKeHq\nkb0jV8rJhy/fbszujTGmTUqYN7iDQnmLaC/nle2FI3uar2HGGHMKS7hgAVCYm8XOQ+VsrytvYUlu\nY4xpkIQMFsH3Lep8hLbHEECsu3JjjIlR1GAhIt1F5Hci8oa7PEhEvhP/pjVe/24d6JKezLK6BkNq\n1xGy+lm3H8YYE6NYrizm4TzR1NNd3gjcHq8GNQWPRyjMzWLZligv59mVhTHGxCSWYNFVVf+M23+T\n+0hsqx8QojAvk+0Hytl5qDxyhZx8KN0Gx+pJhBtjjAFiCxZHRSQLt6sOESkCSuPaqiZQmBuln6ge\nw5xP667cGGOiiiVY3IHz5nU/EXkP+CNwW1xb1QTO7NGRjLTkupPc9kSUMcbELJaX8tYCY4CBOJ0D\nbuAUeIrK4xFG5WbWneROz4SM0y1YGGNMDGL5pf+BqvpUda2qrlHVKuCDeDesKRTlZfFVyTF2l9aV\ntxhmwcIYY2JQZ7AQkR4icg6QJiLDReRsdxoLpDdbC09CYbR+onLynfG4jx9pxlYZY8ypp77bUJfg\nDHnaC3g4rPwI8JM4tqnJnJXTiY6pSSzbUsIVwyOMu5STDyjsWQN9Rjd7+4wx5lRRZ7BQ1T8AfxCR\nK1X15WZsU5PxeoTC3Ew+rO/KApxbURYsjDGmTlET3Kr6sohcBgwGUsPKZ8WzYU2lMDeLf6zfx97D\nFXTvlFpzZcce0KG75S2MMSaKWLr7eAK4Bvg+ztNQVwF94tyuJhO9n6hh9q6FMcZEEcvTUOeq6reA\ng6p6HzCamsOltmqDenaiY7ukuh+hzcmHfeuhqqJ5G2aMMaeQWIJF8LfoMRHpCVQBufFrUtPyeoQR\nfbvU/SZ3Tj6oH/atjbzeGGNMTMHiLyLSGXgI+BjYCjwfz0Y1taK8LL4sPsq+IxGuHuxNbmOMiare\nYCEiHuBtVT3kPhHVBzgzfGjUU0Ghm7f4KNKtqM6nQ2pn64HWGGPqUW+wUNUA8Ouw5eOq2uo7Eaxt\nSM9OtE/xRk5yi9ib3MYYE0Ust6HeEpErRUTi3po4SfJ6GNG3nnG5c/Jh71rwVzVvw4wx5hQRa6+z\nLwHHReSwiBwRkcNxbleTK8zL5It9ZewvO37iypwC8B+H4g3N3zBjjDkFRA0WqtpRVT2qmqKqndzl\nTs3RuKZUVF/ewsa2MMaYesW1q3ERuVRENojIJhG5O8L680XkYxHxiciUWuuuF5Ev3On6k23L0NMy\nSE/xRn6ENqsfJLe3vIUxxtQhbsFCRLzAXGA8MAi4TkQG1aq2Daezwv+ttW0m8DOgEBgF/ExEupxM\ne5K9Hs7p0yVyP1EeL/QYasHCGGPqEM8ri1HAJlXdrKqVwAvA5PAKqrpVVT/DHd87zCXA31X1gKoe\nBP4OXHqyDSrKy2LD3iMcOFp54sqcfNizGgK1m2KMMSbqexYisqaR+z4N2B62vMMta7JtReQmEVkh\nIiuKi4uj7jQ4vkXEvEXOMKgsgwObY2yiMca0HbG8Z/GpiJzeiH1HetRWm3JbVX1KVUeo6ojs7Oyo\nOx3WqzOpyZ7I71uE3uReFWMTjTGm7YjlNlQOsFZE3haRBcEphu12ULPDwV7ArhjbdTLb1iklyclb\nROxUMPtM8KZY3sIYYyKIOp4FcF8j970c6C8iucBO4FrgX2Pc9k3gl2FJ7YuBGY1sRw1FuVk8/I+N\nHDpWSef0lOoV3mToPtiChTHGRBDLexbvAp8DHd1pvVsWbTsf8D2cX/zrgT+r6loRmSUilwOIyEgR\n2YEzRsaTIrLW3fYA8HOcgLMcmOWWnbTCvCxU63nfYs9noLHeLTPGmLYhlsGPrgY+wvmFfjWwrPY7\nEXVR1YWqOkBV+6nq/3PL7lXVBe78clXtpartVTVLVQeHbft7VT3DnZ5pzJeLJL93Bu2SPJFvReXk\nQ/lBKN1+4jpjjGnDYrkN9VNgpKruAxCRbOAfwPx4Nixe2iV5Ofv0LnUkuQucz92fOr3RGmOMAWJL\ncHuCgcJVEuN2rVZhXibrdh+mtLxWx4HdB4F4LW9hjDG1xPJL/28i8qaITBeR6cBfgYXxbVZ8Fbl5\ni+W1b0Ulp0H2QBvbwhhjaoklwf1j4ElgGJAPPKWqd8W7YfFU0LszKUkelm2p430Lu7Iwxpga6s1Z\nuP07vamqFwL/1zxNir/UZC8FvTvXneT+9Hk4sgc69mj+xhljTCsU7Q1uP3BMRDKaqT3NpigvizU7\nSzlcUStvEXqT225FGWNMUCzRfCG7AAAdw0lEQVQ5iwpgtYj8TkTmBKd4NyzeinIzCSis3Hqw5oru\nQ5zPPXYryhhjgmJ5dPav7pRQhp/ehRSv00/UuDO7Va9I7QSZ/SxvYYwxYWLJWVykqtOaqT3NJi3F\nS37vDD6sK2+xc0XzN8oYY1qpWHIW2SKSUl+9U1VhrpO3KDvuq7kiJx8ObYNjTdLDiDHGnPJiyVls\nBd4TkXtE5I7gFOd2NYuivCz8AWXF1lpBISc4Jvfq5m+UMca0QrEEi13A627djmHTKe/sPp1J8siJ\nj9D2CD4RZXkLY4yBGBLcqnofgIi0V9Wj8W9S80lPSWJYr4wT+4lqnwUZvS1YGGOMK5ZeZ0eLyDqc\nbsYRkXwR+Z+4t6yZFOVlsXpHKccqI+QtLFgYYwwQ222oR4FLcDoQRFU/Bc6PZ6OaU2FeFr6AsvKr\nWu9b9BgGJZvgeFnLNMwYY1qRmHqPVdXaAzz449CWFjGiTxe8HjnxVlROPqCwd02LtMsYY1qTWILF\ndhE5F1ARSRGRO3FvSSWC9u2SGHpaBss2134iypLcxhgTFEuw+C5wK3AasAMocJcTRmFeJp/uOER5\nZdgFU8ce0L6bBQtjjCG2Lsr3q+pUVe2uqt1UdZqqRujb+9RVlJdFlV/5eFtY3kLEed/COhQ0xphT\ne8S7pjKiTxc8Assi5S2K10NVRcs0zBhjWgkLFkDH1GSGnpbBh5HyFgEf7FvXMg0zxphWwoKFqzAv\ni1XbD1FRFZa3sCS3McYAUYKFiIwRkWHu/NUi8hsR+aGItGue5jWfwtxMKv2BmnmLzn0gNQP2WN7C\nGNO21RksRGQu8AvgaRF5DvhXYA0wHPh9LDsXkUtFZIOIbBKRuyOsbyciL7rrl4lIX7c8WUT+ICKr\nRWS9iMxo+FdrmBF9M928RditKBHn5Ty7sjDGtHH19Q01TlUHiUgqsBPopqp+EXkSiPqntjsWxlzg\nIpxHbpeLyAJVDU8AfAc4qKpniMi1wK+Aa4CrgHaqOlRE0oF1IvK8qm5tzJeMRUZaMoN6dmLZlghJ\n7o9+C/4q8CbH6/DGGNOq1XcbqgJAVSuAr9yxLVBVBarq2S5oFLBJVTeraiXwAjC5Vp3JwB/c+fnA\nBSIigALtRSQJSAMqgcOxfaXGK8rN4uNttfMWBeA/Dvs3xvvwxhjTatUXLLq5Y1f8KGw+uJwdw75P\nA8K7CdnhlkWso6o+oBTIwgkcR4HdwDZgtqqeMBKRiNwkIitEZEVxcXEMTapfYV4Wlb4An24/VF0Y\nHNvC3rcwxrRh9QWL3+KMW9EhbD64/HQM+5YIZRpjnVE4/U/1BHKBH4lI3gkVVZ9S1RGqOiI7O5b4\nVb9RfTMRoeYjtFlnQHK65S2MMW1anTmL4DgWJ2EH0DtsuRfOQEqR6uxwbzllAAdwkul/U9UqYJ+I\nvAeMADafZJvqlZGezFk9gnmL/k6hxws9hlqwMMa0adEenR0vIktEZL+IFIvIuyIyIcZ9Lwf6i0iu\nO4b3tcCCWnUWANe781OAd9ycyDbgG+JoDxQBn8f6pU5GYV4mH287yHFfWN6ixzBniNVAoDmaYIwx\nrU59j87eCPwcmAnkAf2A+4CZInJTtB27OYjvAW/i9FL7Z1VdKyKzRORyt9rvgCwR2QTcAQQfr52L\nc7trDU7QeUZVmyVpUJSXRUVVgM92lFYX5uRD5RE4uKU5mmCMMa1OfY/O/hA4r1Zi+R0RGQ8sBZ6K\ntnNVXQgsrFV2b9h8Bc5jsrW3K4tU3hxG9c0EnH6iRrrz1W9yr4Ksfi3RLGOMaVH13YaSSE8gJVqP\ns7V1aZ/CmT061kxyZ58J3hTLWxhj2qz6gsVhEcmvXeiWHYlfk1peUV4WK786SJXfzVEkpUC3s+zx\nWWNMm1VfsPgRsEBEZorIJBGZKCL3Aa/h5BcSVmFuJuVV/hPzFrs/Ba399K8xxiS+OoOFqi4FCt06\n04Eb3Pkid13CGpXr5CpqjMudkw/lB6B0Rwu1yhhjWk59CW5UdY+I/BI4A+dluS/dpHRCy+rQjgHd\nO7BsywFuHecW5hQ4n7s/hc6969zWGGMSUX2PziaJyIM43XH8AXgO2C4iD4pIwveoV5ibxYqtB6rz\nFt0GgXisu3JjTJtUX87iISATyFPVc1R1OM67Fp2B2c3RuJZUlJfFsUo/a3a6eYuUdOg60J6IMsa0\nSfUFi4nAjaoaevJJVQ8DNwOxvsV9ygrmLZZtCXuENpjkNsaYNqa+YKFu1xu1C/2c2CFgwsnu2I4z\nunU4Mcl9ZDeU7Wu5hhljTAuoL1isE5Fv1S4UkWk0Uz9NLa0wN5MVWw/iC+YtrLtyY0wbVV+wuBW4\nVUQWi8ivRWS2iLwL3IZzKyrhFeZlUXbcx9pd7rhLPYY6n7tXtVyjjDGmBdTXRflOoFBEvgEMxhl7\n4g1Vfbu5GtfSikJ5ixLye3eG1AzIzLO8hTGmzan3PQsAVX0HeKcZ2tLqdOuUSl7X9izbfICbznc7\nEMzJh50ft2zDjDGmmdU7noVxbkV9tOUA/oCb0+8xDA59BeUHW7ZhxhjTjCxYRFGUl8mR4z7W73bz\nFsHuyvesbrlGGWNMM7NgEUVhbhYQ1k9UaGwLy1sYY9oOCxZR9MhIpW9WevX4Fu27QqdeFiyMMW2K\nBYsYFOZmsXzrAQLBvEXOMHvXwhjTpliwiEFRv0xKy6tYvycsb7F/I1QebdmGGWNMM7FgEYNg3mJZ\n8FZUTj6gsGdNyzXKGGOakQWLGPTsnEbvzDRLchtj2iwLFjEqys3io2DeomMOpHeFPRYsjDFtQ1yD\nhYhcKiIbRGSTiNwdYX07EXnRXb9MRPqGrRsmIh+IyFoRWS0iqfFsazSFeVkcOlbFxn1HQMS6KzfG\ntClxCxYi4gXmAuOBQcB1IjKoVrXvAAdV9QzgEeBX7rZJOCPzfVdVBwNjgap4tTUWhcFxub8MuxW1\nbz34jrdgq4wxpnnE88piFLBJVTeraiXwAjC5Vp3JOEO2AswHLhARAS4GPlPVTwFUtcQdR6PF9M5M\n57TOadWDIeXkQ8AH+9a1ZLOMMaZZxDNYnIYzfnfQDrcsYh1V9QGlQBYwAFAReVNEPhaR/4xjO2NW\nmJfJsi0HUFUb28IY06bEM1hIhLLaI+zVVScJOA+Y6n7+i4hccMIBRG4SkRUisqK4uPhk2xtVUV4W\nB45W8sW+MuiSC+0yLG9hjGkT4hksdgC9w5Z7AbvqquPmKTKAA275u6q6X1WPAQuBs2sfQFWfUtUR\nqjoiOzs7Dl+hpqLQ+xYlbpJ7mAULY0ybEM9gsRzoLyK5IpICXAssqFVnAXC9Oz8FeMcd9/tNYJiI\npLtBZAzQ4smB3plp5GSkVvcTlZMPe9eA39eyDTPGmDiLW7BwcxDfw/nFvx74s6quFZFZInK5W+13\nQJaIbALuAO52tz0IPIwTcFYBH6vqX+PV1liJCEV5WSzbUuLkLXoMA18FlHzR0k0zxpi4ijpS3slQ\n1YU4t5DCy+4Nm68Arqpj2+dwHp9tVQpzM3nlk518WVzGGeFvcnc7q2UbZowxcWRvcDdQUV5wfIsD\n0LU/JKVZ3sIYk/AsWDRQn6x0undq57xv4fFCj6EWLIwxCc+CRQOJCIW5WXy4uaT6fYs9qyEQaOmm\nGWNM3FiwaISivCyKjxxny/6jzhNRxw/DwS0t3SxjjIkbCxaNUJjn9hO1+YB1V26MaRMsWDRCXtf2\nZHdsx7ItJZB9FniSLVgYYxKaBYtGcPIWmSzbfAD1JjuPze6xPqKMMYnLgkUjFeZlsedwBV+VHKse\n20Jrd31ljDGJwYJFI4128xbLtpQ4weJYCRze2cKtMsaY+LBg0Uj9sjvQtUOKm+QucAotb2GMSVAW\nLBop+L7Fss0laPdBIB4b28IYk7AsWJyEwrxMdpVWsKNMoOsAu7IwxiQsCxYnodAd3+KDzSXVSW5j\njElAFixOQv9uHchsn8Ky4Mt5R3ZBWfxH7DPGmOZmweIkeDzCqL6ZfLi5xBnbAmCPXV0YYxKPBYuT\nVJSXyc5D5exMO8MpsFtRxpgEZMHiJBUGx7fY6Yes/vDPh+Hlf4f1r0NVeQu3zhhjmkZcR8prCwZ2\n70jn9GQ+3FzClVf/EZY97gSK1S9BSgcYcAkMmgxnXAQp6S3dXGOMaRQLFicpmLdYtuUAXDUOLn8M\nLnsYti6Fda/C+r/AmpchuT0MuNgJHP0vhpT2Ld10Y4yJmQWLJlCYl8Vb6/ay61A5PTungTcZ+o1z\npgm/hq/eqw4ca1+B5HTof5EbOC6Bdh1a+isYY0y9LFg0gaKwfqL+ZXivmiu9SZA3xpkmzIav3ncC\nx7oFsO41SEp1A8cVzi2rdh1b4BsYY0z9LFg0gTN7dKJTahLLNh84MViE83gh9+vONP5B2PahEzDW\nveZcdXjbwRkXwuArYMClkNqp+b6EMcbUw4JFE/B6hFG57vsWsfJ4oe/XnOnSB2D7surAseGv4E2B\nfhc4gWPgeEjNiN8XMMaYKOL66KyIXCoiG0Rkk4jcHWF9OxF50V2/TET61lp/uoiUicid8WxnUyjK\ny2JryTH2lFY0fGOPB/qMhvEPwA/Xwg1vwcgbYc9qeOU/4MF+8KerYdX/Qvmhpm+8McZEEbcrCxHx\nAnOBi4AdwHIRWaCq68KqfQc4qKpniMi1wK+Aa8LWPwK8Ea82NqVgP1HLtpQwueC0xu/I44HTC53p\n4l/Aro+dpPi6BfDFm84Qrnlj3SuOCZCe2STtN8aY+sTzNtQoYJOqbgYQkReAyUB4sJgMzHTn5wO/\nERFRVRWRK4DNwNE4trHJDOrZiY7tknh88ZccrvAxdkA2vTNP8r0Kjwd6jXCmUOB41UmQv3YreJIg\nd4wTOM6caIHDGBM38QwWpwHbw5Z3AIV11VFVn4iUAlkiUg7chXNV0upvQYGTt/jpZWcxd/Em7nl1\nDQB52e0ZO6AbYwdmMyo3k9Rkb+MPIAKnneNMF82C3auqA8eC78Nfbofc853Hcc+cCB2ym+ibGWNM\nfIOFRCirPUh1XXXuAx5R1TKRSFXcjUVuAm4COP300xvZzKZz7ajTuWZkbzbvP8riDcW8u7GY55Z9\nxe/f20JqsofReVmMHegEjz5ZJ/FSngj0HO5MF86EPZ9VB47Xb3emtC7Q+XTo3KfWpzvZux3GmAYQ\n1dq/v5toxyKjgZmqeom7PANAVe8Pq/OmW+cDEUkC9gDZwBKgt1utMxAA7lXV39R1vBEjRuiKFSvi\n8l1ORnmlnw83l/DuxmIWb9jH1pJjAPTNSmfswG6MGZhNUW4WaSkncdURpAp718CX78DBrXBoW/Xk\nq5V4T8t0gkaXCMEko7cFE2PaCBFZqaojotaLY7BIAjYCFwA7geXAv6rq2rA6twJDVfW7boL7m6p6\nda39zATKVHV2fcdrrcGitq37j4YCxwebS6ioCtAuyUNhXhZjB2QzdmA2uV3bU98VVYOpwtFiN3B8\nBQe/qhlIDm0D//Ga26RnnXg1Elrubd2VGJMgWjxYuI2YADwKeIHfq+r/E5FZwApVXSAiqcCzwHDg\nAHBtMCEeto+ZJFCwCFdR5eejLQdYvKGYxRv3sbnYyeX3zkwL5TpG98siPSXOr8MEAjWDyaHawWT7\nicGkfXatQBIWTDJ6W6eJxpwiWkWwaE6nYrCobfuBYyzeWMy7G4p5/8v9HKv0k+L1MCo3k7EDnauO\nftkdmvaqIxaBABzdFxY8al2dlG4Hf2XNbdp3g47dIb2rc5WSngXtuzpPbKVn1SxPz3T60zLGNDsL\nFqe44z4/K7YeDN2y2ri3DIDTOqcxZmA2Ywdkc+4ZXenQrhW8hB8IQNnemsHk0FdQtg+OlTjT0RI4\nXlr3PlIzagWWrLBg0vXEYNOuk5PoN8acFAsWCWbnoXLe3VDMuxv38d6mEsqO+0j2CiP6OFcdYwZm\nM7B7x+a/6mgIXyWUHwgLIPvd+QNwbH9YWbDO/hOvWII8yfUHlvRMJ7ikdXH63PImQ1I7pxsVb7L7\nmeK8q9Kaz5kxcWbBIoFV+gKs/Kr6quPzPUcAyMlIZcwAJ8/RrWMqGWnJdE5PJiMtmfQUb+sOJJGo\nQmXZiQGkRrA5ULOs/GADDyLVgSMUUJJrlnlrB5lIgSfCduH7Sk53ppT2zqBYKe2dvE5w3ptiQcu0\nCAsWbcie0gqWbHSS5P/8Yj9HKnwn1EnyCBlpTuDo5H5GmmqscwNN+1Mp0Ph9TsAIBpbyg87Vib/K\n/ax0rnBqlB2vud5fBb7jtepUufXCykL7qbWvxvAkOQNkpdQz1VjfwQ027nxdgchyQSYKCxZtlM8f\n4It9ZRw8VknpsSpKyyNPh2stB+r5b5DkkVAQqRlokpyrl7SUiIGmQ7sk0lO8JHvb0FDvqhDwnRh4\nqsqh6ihUBqcyqDxWPV8VNl95tOa6yqPu+jI4Xgbqj7093pTqQJOc5vR2LF4Qj3Ml4wnOu5+h5bCp\nRp36tpEIZcHlWvsEQMKupmKcD30E56Vx8yJ1fMfa36ue79So8xd+LrzV/x7BOp6ksLKk6jbEUazB\nohVkR01TSvJ6OCunYeNgqCplx30cOhY5kJwwHatkW8nRmAINQLJXSE9xAkdaipf0FG9oOXw+LcVL\nenIS7dvFUC8lifRkLx5PK7viEXFvPSUDcXgXRdUJQpVHa05VR08sqx2Iqo5BwA8aqDmFyvxusPOD\nVrnL4XW0uqzGNoGw7QK16tTexl/9PYIdOoTPmwgkQhDx1AooXqfMk1QdfMIDUX3lMbJgYRAROqYm\n0zE1OfTafKwCAaWs0he6igkPNGXHfRyr9HOs0k95ZfX8MXf+wNFKth/wUV7p51iVn2PH/VT6Aw06\nfmqyp0ZASXODiBNwkkhL9pDk9ZDsEZK8HpK8QrLH/fR6SPIIXo87H7auxjYecco8HpK91WXJXo+7\nrdR5jCSPNO0tPBEnF5LULnE7jtQIQSR0B0SbZj4Y4EJBs3YQrCuYBpyn/xqzXY1t/dXrA77q+gFf\n9bpgcI2pPFCrTj379FfWLI+RBQtzUjweoVNqMp0aEWgi8fkDHKvyU17p56gbbMqrnPnyYLCp8nMs\nyrpdh6oor3ICk8+v+AKKzx+gyv2MdjXUlLxuQPIIeEXwiOAJLrvBJLhOpLquU0fccpz5OtZ5PTX3\n6wkep9a6ZDfweYOBLxQATwyGSW7wC9YLBthQgKwVdIMB19m25rrw/QS/Q51BVKTmp4mvm2M7zxYs\nTKuS5PXQyeuhU2p8E7OBgFIVCNQIJL6AUuUPlgWo8mtovt51NcqCdavnfYEA/gAEVAkEFL+qczcp\noE6ZOu0JaKR1SiCAW+7UjbSuyh9wy53bin53XbBecF31d6hue5XbTn9zRlDcIOoGwPDJIxK64vN4\nIMnjCQU9r8eD10ON7TziBLtgME3yVM+HpvD6btBM9npITvKQ4vWQ4n4me4WUJK/7GSxz1jufQorX\nS3KShNa1c9clh+3jlHkgpAEsWJg2yeMR2nm8tIZ3GluLQCAsiATcYBh2NVZVIzieWFYVCo5uoPRX\nB+TQOr8TOP2q+MOCqM9fHdR8geqg6neDmF8Vvz8YBN06YeuPVwVqrAsG1Br78gf3iRPcfU77G3rr\nMxbVgccNJGFBxSlzgpWIezdMnfOguJ9u0A8uBwLODTl114XXRQn90aE4dQnbT3AdVP8xEfxsyPNN\n9qNijAGcAJriEVLiO9pyq6PqXAlW+QNU+gJU+QMcdz8r/QGqfEql30+lr2adSnfeqeN++jW0bXVZ\ncH9Kpc/vfjrrNEDoNmKSe+tQ3Ft0AqFbijXLqm9LErrlSGhdsCy47PEABG/9Va8Tdz9rYjxPFiyM\nMW2aiDi3l5I8tG/X0q1pfjNjrNe2/oQwxhjTKBYsjDHGRGXBwhhjTFQWLIwxxkRlwcIYY0xUFiyM\nMcZEZcHCGGNMVBYsjDHGRJUw41mISDHwVUu3o5l1Bfa3dCNaCTsX1excVLNzUa2uc9FHVbOjbZww\nwaItEpEVsQxa0hbYuahm56KanYtqJ3su7DaUMcaYqCxYGGOMicqCxantqZZuQCti56KanYtqdi6q\nndS5sJyFMcaYqOzKwhhjTFQWLFopEfm9iOwTkTVhZZki8ncR+cL97OKWi4jMEZFNIvKZiJzdci1v\neiLSW0QWich6EVkrIj9wy9vc+RCRVBH5SEQ+dc/FfW55rogsc8/FiyKS4pa3c5c3uev7tmT740FE\nvCLyiYi87i63yXMhIltFZLWIrBKRFW5Zk/2MWLBoveYBl9Yquxt4W1X7A2+7ywDjgf7udBPweDO1\nsbn4gB+p6llAEXCriAyibZ6P48A3VDUfKAAuFZEi4FfAI+65OAh8x63/HeCgqp4BPOLWSzQ/ANaH\nLbflczFOVQvCHpFtup8RdQeCt6n1TUBfYE3Y8gYgx53PATa4808C10Wql4gT8BpwUVs/H0A68DFQ\niPOyVZJbPhp4051/Exjtzie59aSl296E56CX+0vwG8DrOKOLttVzsRXoWqusyX5G7Mri1NJdVXcD\nuJ/d3PLTgO1h9Xa4ZQnHvXUwHFhGGz0f7m2XVcA+4O/Al8AhVfW5VcK/b+hcuOtLgazmbXFcPQr8\nJxBwl7Nou+dCgbdEZKWI3OSWNdnPiI3BnRgkQlnCPeYmIh2Al4HbVfWwSKSv7VSNUJYw50NV/UCB\niHQGXgHOilTN/UzYcyEiE4F9qrpSRMYGiyNUTfhz4fqaqu4SkW7A30Xk83rqNvhc2JXFqWWviOQA\nuJ/73PIdQO+wer2AXc3ctrgSkWScQPEnVf0/t7jNng8AVT0ELMbJ43QWkeAff+HfN3Qu3PUZwIHm\nbWncfA24XES2Ai/g3Ip6lLZ5LlDVXe7nPpw/IkbRhD8jFixOLQuA693563Hu3QfLv+U+4VAElAYv\nPROBOJcQvwPWq+rDYava3PkQkWz3igIRSQMuxEnuLgKmuNVqn4vgOZoCvKPuTepTnarOUNVeqtoX\nuBbnu02lDZ4LEWkvIh2D88DFwBqa8mekpZMyNtWZrHoe2A1U4fwV8B2c+6tvA1+4n5luXQHm4ty7\nXg2MaOn2N/G5OA/nEvkzYJU7TWiL5wMYBnzinos1wL1ueR7wEbAJeAlo55anusub3PV5Lf0d4nRe\nxgKvt9Vz4X7nT91pLfBTt7zJfkbsDW5jjDFR2W0oY4wxUVmwMMYYE5UFC2OMMVFZsDDGGBOVBQtj\njDFRWbAwxhgTlQULY06CiBSIyISw5ctF5O76tmnAvm8XkfSm2JcxJ8veszDmJIjIdJwXmr4Xh31v\ndfe9vwHbeNXpO8qYJmVXFqZNEJG+7uBJv3UHDXrL7S4jUt1+IvI3t/fOf4rImW75VSKyxh14aIk7\nqM4s4Bp3wJlrRGS6iPzGrT9PRB4XZ+CmzSIyRpxBrdaLyLyw4z0uIiuk5mBGtwE9gUUissgtu84d\n3GaNiPwqbPsyEZklIsuA0SLygIiscwe1mR2fM2ranJZ+Td0mm5pjwhkbxAcUuMt/BqbVUfdtoL87\nX4jThxA43SKc5s53dj+nA78J2za0jDOA1Qs4XStMBg4DQ3H+SFsZ1pZgFwxenI4Bh7nLW3HHJ8AJ\nHNuAbJzeot8BrnDXKXB1cF84YxNIeDttsulkJ7uyMG3JFlVd5c6vxAkgNbjdoJ8LvOSOGfEkzqAx\nAO8B80TkRpxf7LH4i6oqTqDZq6qrVTWA039P8PhXi8jHOH0+DQYGRdjPSGCxqharMxbDn4Dz3XV+\nnB55wQlIFcDTIvJN4FiM7TSmXjaehWlLjofN+4FIt6E8OIPnFNReoarfFZFC4DJglYicUKeeYwZq\nHT8AJIlILnAnMFJVD7q3p1Ij7KfOwTuACnXzFKrqE5FRwAU4PbF+D6frbmNOil1ZGBNGVQ8DW0Tk\nKggNbJ/vzvdT1WWqei/OkJy9gSNAx5M4ZCfgKFAqIt1xxkYOCt/3MmCMiHQVES9wHfBu7Z25V0YZ\nqroQuB1nnG5jTppdWRhzoqnA4yLyX0AyTt7hU+AhEemP81f+227ZNuBu95bV/Q09kKp+KiKf4NyW\n2oxzqyvoKeANEdmtquNEZAbOWA0CLFTV107cIx2B10Qk1a33w4a2yZhI7NFZY4wxUdltKGOMMVHZ\nbSjTZonIXJxxnMP9t6o+0xLtMaY1s9tQxhhjorLbUMYYY6KyYGGMMSYqCxbGGGOismBhjDEmKgsW\nxhhjovr/VAkbeGHmLX0AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x1ecad4a3080>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "train_X = np.array([img.flatten(1) for img in mnist.train.images])\n",
    "train_y = mnist.train.labels\n",
    "\n",
    "mod_optimization (10, 510, 50, train_X, train_y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fit the model\n",
    "\n",
    "After n_estimator is larger than 200 or 300, the classification error becomes relatively stable. Thus, we selected 300 as the n_estimator after balancing the time cost and increased accuracy."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Model trained!\n",
      "Mean accuracy: 1.0\n",
      "Elapse time:  0.7186755140622457 minutes\n"
     ]
    }
   ],
   "source": [
    "start = time()\n",
    "rf_mnist = RandomForestClassifier(n_jobs=-1, n_estimators=300)\n",
    "rf_mnist.fit(train_X, train_y)\n",
    "end = time()\n",
    "\n",
    "et = (end - start)/60.\n",
    "print(\"Model trained!\")\n",
    "print (\"Mean accuracy:\", rf_mnist.score(train_X,train_y))\n",
    "print (\"Elapse time: \", et, \"minutes\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Evaluate the model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Precision:\n",
      "0 : 0.972945891784\n",
      "1 : 0.991166077739\n",
      "2 : 0.963391136802\n",
      "3 : 0.962413452028\n",
      "4 : 0.971544715447\n",
      "5 : 0.975113122172\n",
      "6 : 0.974039460021\n",
      "7 : 0.9736328125\n",
      "8 : 0.963768115942\n",
      "9 : 0.961\n",
      "\n",
      "Recall:\n",
      "0 : 0.990816326531\n",
      "1 : 0.988546255507\n",
      "2 : 0.968992248062\n",
      "3 : 0.963366336634\n",
      "4 : 0.973523421589\n",
      "5 : 0.966367713004\n",
      "6 : 0.979123173278\n",
      "7 : 0.969844357977\n",
      "8 : 0.955852156057\n",
      "9 : 0.95242814668\n",
      "Overal scores: 0.9711\n"
     ]
    }
   ],
   "source": [
    "test_X = np.array([img.flatten(1) for img in mnist.test.images])\n",
    "test_y = mnist.test.labels\n",
    "\n",
    "test_y_pred = rf_mnist.predict(test_X)\n",
    "\n",
    "\n",
    "result = precision_recall_fscore_support(test_y, test_y_pred)\n",
    "classes = [i for i in range(10)]\n",
    "print (\"Precision:\")\n",
    "for i in range(len(classes)):\n",
    "    print(classes[i],':', result[0][i])\n",
    "print (\"\\nRecall:\")\n",
    "for i in range(len(classes)):\n",
    "    print(classes[i],':', result[1][i])\n",
    "    \n",
    "print(\"Overal scores:\", accuracy_score(test_y, test_y_pred))"
   ]
  }
 ],
 "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.6.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}