{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "[![Open In Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/worldbank/OpenNightLights/blob/master/onl/tutorials/mod6_6_RF_classifier.ipynb)\n", "\n", "# Random Forest Classifier\n", "\n", "Now that we have processed and explored our data, we will try to classify built-up areas with a Random Forest ensemble of decision trees.\n", "\n", "Decision tree models like Random Forest are among the most powerful, easy to use, and simple to understand models in the machine learning portfolio. The resources noted in {doc}`mod6_2_supervised_learning_img_classification` are a good place to start.\n", "\n", "More in-depth context on methods like Random Forest are out of scope for this tutorial, but worth understanding well if you use them for your analysis, even exploration.\n", "\n", "## Training data\n", "\n", "Let's recreate the training data \"image\" for 2015 that fuses Sentinel-2, VIIRS-DNB and GHSL for the Bagmati province." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import pandas as pd\n", "\n", "# reminder that if you are installing libraries in a Google Colab instance you will be prompted to restart your kernal\n", "\n", "try:\n", " import geemap, ee\n", " import seaborn as sns\n", " import matplotlib.pyplot as plt\n", "except ModuleNotFoundError:\n", " if 'google.colab' in str(get_ipython()):\n", " print(\"package not found, installing w/ pip in Google Colab...\")\n", " !pip install geemap seaborn matplotlib\n", " else:\n", " print(\"package not found, installing w/ conda...\")\n", " !conda install mamba -c conda-forge -y\n", " !mamba install geemap -c conda-forge -y\n", " !conda install seaborn matplotlib -y\n", " import geemap, ee\n", " import seaborn as sns\n", " import matplotlib.pyplot as plt" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "try:\n", " ee.Initialize()\n", "except Exception as e:\n", " ee.Authenticate()\n", " ee.Initialize()\n", "\n", "# define some functions and variables\n", "def se2mask(image):\n", " quality_band = image.select('QA60')\n", " cloudmask = 1 << 10\n", " cirrusmask = 1 << 11\n", " mask = quality_band.bitwiseAnd(cloudmask).eq(0) and (quality_band.bitwiseAnd(cirrusmask).eq(0))\n", " return image.updateMask(mask).divide(10000)\n", "\n", "\n", "se2bands = ['B2', 'B3', 'B4', 'B5', 'B6', 'B7','B8','B8A']\n", "trainingbands = se2bands + ['avg_rad']\n", "label = 'smod_code'\n", "scaleFactor=1000\n", "\n", "# create training data\n", "roi = ee.FeatureCollection(\"FAO/GAUL/2015/level2\").filter(ee.Filter.eq('ADM2_NAME','Bagmati')).geometry()\n", "\n", "se2 = ee.ImageCollection('COPERNICUS/S2').filterDate(\n", " \"2015-07-01\",\"2015-12-31\").filterBounds(roi).filter(\n", " ee.Filter.lt(\"CLOUDY_PIXEL_PERCENTAGE\",20)).map(se2mask).median().select(se2bands).clip(roi)\n", "\n", "viirs = ee.Image(ee.ImageCollection(\"NOAA/VIIRS/DNB/MONTHLY_V1/VCMSLCFG\").filterDate(\n", " \"2015-07-01\",\"2019-12-31\").filterBounds(roi).median().select('avg_rad').clip(roi))\n", "\n", "fused = se2.addBands(viirs)\n", "\n", "# create and overlay labels to training data\n", "ghsl = ee.ImageCollection('JRC/GHSL/P2016/SMOD_POP_GLOBE_V1').filter(ee.Filter.date(\n", " '2015-01-01', '2015-12-31')).select(label).median().gte(2).clip(roi)\n", "\n", "points = ghsl.sample(**{\"region\":roi, \"scale\":scaleFactor,\"seed\":0,'geometries':True})\n", "\n", "data = fused.select(trainingbands).sampleRegions(collection=points,\n", " properties=[label],\n", " scale=scaleFactor)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As a gut-check let's look at the stats:" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'max': 1,\n", " 'mean': 0.18796029458853666,\n", " 'min': 0,\n", " 'sample_sd': 0.39071173874702697,\n", " 'sample_var': 0.15265566279472506,\n", " 'sum': 1174,\n", " 'sum_sq': 1174,\n", " 'total_count': 6246,\n", " 'total_sd': 0.39068046053869543,\n", " 'total_var': 0.15263122224672718,\n", " 'valid_count': 6246,\n", " 'weight_sum': 6246,\n", " 'weighted_sum': 1174}" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "data.aggregate_stats(label).getInfo()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In our total training dataset, we have 6,246 observations and our label clasification shows that about 19% (the mean value above) of our data are classified as built-up (1) with the rest not.\n", "\n", "## Cross-validation\n", "\n", "Cross-validation is one of the most important aspects of machine learning development. Last section we talked about several attributes of the training data that may impact classification:\n", "- varying spatial resolution and choices for re-sample rate\n", "- choices about data cleaning\n", "- decisions about which bands to include\n", "\n", "We also may decide to create new features (known as feature engineering) by transforming our data sources mathmatically (getting derivatives, information about neighboring pixels, etc) or even fusing additional data.\n", "\n", "The classification algorithm itself will have hyperparameters than may be adjusted (known as hyperparameter tuning). \n", "\n", "How do we decide these things? \n", "\n", "Often we will experiment empirically and see what works, this is a big advantage to advances in computing resources and machine learning packages. However, if we just tweak our data until we get the best performance on our training data and leave it at that, we are at risk of over-fitting our model. \n", "\n", "Over-fitting a model means making it too specific to our data on-hand in a way that will fail us on unseen data. This could impact the final analysis and ultimately stakeholder trust and their ability to make informed decisions, so we want to think about strategies to avoid this.\n", "\n", "That is our situation here: since we dont have any \"ground truth\" for our data after 2015 in terms of settlements, we will want to validate our classifier as best we can with the labeled data we have before \"releasing it to the wild.\"\n", "\n", "### train test split\n", "A key way to do that is to split our labeled data into two components: training and testing sets (or even train, validation and test sets). There are many strategies for this, including a K-fold sample technique (or a stratified K-Fold technique, which addresses the class imbalance issue we noted earlier). Unfortunately things can get quite complex with time series or sequential data (since observations in time are \"dependent\" on observations before, we cannot fairly randomly split them). \n", "\n", "For our purposes, a simple 80/20 train/test split randomly among the pixels in our 2015 training image will be fine, but this is another great topic to learn more about.\n", "\n", "
\n", "Python's scikit-learn documentation is a good place to start for understanding concepts like cross-validation and comes with programming packages and examples.
" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# we'll create a column of random numbers\n", "data = data.randomColumn(seed=0)\n", "split_thresh = 0.8\n", "\n", "train = data.filter(ee.Filter.lt('random',split_thresh))\n", "test = data.filter(ee.Filter.gte('random',split_thresh))" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'max': 1,\n", " 'mean': 0.1894335074327039,\n", " 'min': 0,\n", " 'sample_sd': 0.3918919561595822,\n", " 'sample_var': 0.15357930530258393,\n", " 'sum': 943,\n", " 'sum_sq': 943,\n", " 'total_count': 4978,\n", " 'total_sd': 0.3918525917924336,\n", " 'total_var': 0.1535484536944476,\n", " 'valid_count': 4978,\n", " 'weight_sum': 4978,\n", " 'weighted_sum': 943}" ] }, "execution_count": 5, "metadata": {}, "output_type": "execute_result" } ], "source": [ "train.aggregate_stats(label).getInfo()" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "{'max': 1,\n", " 'mean': 0.19444444444444445,\n", " 'min': 0,\n", " 'sample_sd': 0.39592955855613005,\n", " 'sample_var': 0.15676021533845202,\n", " 'sum': 245,\n", " 'sum_sq': 245,\n", " 'total_count': 1260,\n", " 'total_sd': 0.39577241246597245,\n", " 'total_var': 0.1566358024691358,\n", " 'valid_count': 1260,\n", " 'weight_sum': 1260,\n", " 'weighted_sum': 245}" ] }, "execution_count": 6, "metadata": {}, "output_type": "execute_result" } ], "source": [ "test.aggregate_stats(label).getInfo()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have ~5,000 observations (pixels) in our training set and ~1,200 in our test set.\n", "\n", "And notice balance of our label is similar for training and testing (~18% and ~20% classified as built-up, respectively), which we'd expect in a random split." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Train the classifier\n", "\n", "Here we fit the Random Forest estimator to our training data using some initial hyperparameters." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "init_params = {\"numberOfTrees\":100, # the number of individual decision tree models\n", " \"variablesPerSplit\":None, # the number of features to use per split\n", " \"minLeafPopulation\":1, # smallest sample size possible per leaf\n", " \"bagFraction\":0.5, # fraction of data to include for each individual tree model\n", " \"maxNodes\":None, # max number of leafs/nodes per tree\n", " \"seed\":0} # random seed for \"random\" choices like sampling. Setting this allows others to reproduce your exact results even with stocastic parameters\n", "\n", "clf = ee.Classifier.smileRandomForest(**init_params).train(train, label, trainingbands)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Evaluation\n", "\n", "We need metrics to determine how well the classifier performs. We just have a simple binary schema so we can visaulize a confusion matrix, which shows the actual labels (y axis) and the predicted labels (x axis)." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "trainResults = clf.confusionMatrix().getInfo()\n", "trainCM = pd.DataFrame(np.asarray(trainResults), index=['not','built-up'], columns=['not','built-up'])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We will use Python's Seaborn plotting package to make our matrices. The Seaborn plotting library is popular and versatile." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAJrCAYAAAAMBPz+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA5BElEQVR4nO3dd7wcZbnA8d+TiDQVBBRS0ISiFAVUigUxgFIURMULgqDxqlxAru0qYqV6LVexohhbwAKiqCAgRToiEpp0kBJIQgDpSE/Oc/+YOclm2TnnbHLO2c3s75vPfvbszDsz78zuyT7nectEZiJJktSLxnS6ApIkSZ1iICRJknqWgZAkSepZBkKSJKlnGQhJkqSeZSAkSZJ6loGQhl1EZEScNwz7OS8ilpr5HSJiUnnu0ztdl6VZRKwbEX+IiLvL6/nQKBxzanmsqSN9rF4REVPKa3pIp+siDcRAqIbK/3zaeUztdJ21UC9/gUTEWOCPwFuBU4BDga92sk69KiIOKT+HUzpdF2kkPafTFdCIOLTFso8DKwHfAR5qWnfVMB9/feDxYdjP+4AVhmE/WnpMBjYAfpyZ+4zicf8AXALMHcVj1t2lFP8X3NfpikgDMRCqocw8pHlZmfVZCfh2Zs4c4ePfOEz7uXM49qOlyvjy+a7RPGhmPgw8PJrHrLvMfBwYlv8LpJFk01iP6++HExHPjYgvRcRNEfFUfz+XiFgpIj4dEedExOyIeDoi/hURJ0fEayv2+aw+Qo1p9oh4d0RcGhGPR8QDEXF8REyoqlvTsgXNRhGxSUScGhEPlfs6PyJeX1GncRHx84i4NyKeiIirIuL9i9MMFRHPj4gjy+vxZETcGBGfpOL3KSJeFhFfjYjLymv3VETcERHTImJiU9npwLnly4ObmjCnlGXafk8GOZ8VIuIzZf0ejYh/R8QNEfHdiFi9qey4iDgqImY2HPf3EfGaFvtd0O8mIrYu389HI+KR8n1bv6l8Aue3OPdD+q9N+XpSi2O1fB8jYq3yOt9Svu8PRMQ1EXF0RKzaqq4t9v2aiDix/Oz0v3c/iIhxLcouqGNE/Fd5rCcj4p6yHitVvxPP2lfj78weEXF5+Tm/q/z8LVuW26a8to9ExIMR8YvGc2vY39ZlHa4vyz4REddGxMERsVxT2ZnAweXLcxs/hy3Oda2I+O+IuLrc53lV70lEvKtcdklELNN0zFc0nN+Lh3qdpCVlRkj9TgQ2A/5M0Ufj3nL5+sCXgQuAU4EHgZcAbwd2jIidM/P0No6zf7ntyRRfelsAuwMbR8QmmfnUEPezKXAg8DfgJ2WddgXOLvdzU3/B8j/Vi4FJ5XlcDKwB/AA4s426U375nE1xrf4B/ApYGfgi8KaKzd4F7EsR4FwMPA1sCHwI2DkiNs3MOWXZP5bP76e4Puc17Gdm+Txs70lEvLCs18bATcDPyvqtDfwn8HvgnrLsZOAiiqzNOcBxwJrAfwBvi4hdM/OUFofZCdiF4rN1NEXT11uBzSJig8zsbzo5lOI9aj73xmswZGWgMgN4AXAaxWd8OYrmt72B7wP3D7KPncrtAvgdcAfwGmA/YJeIeENFhvXrwPbAnyg+Y1sDHwbWAbZp81T+G9iR4rNxHrAd8AlglYg4CTie4nMwDXg9sBewWrlNo88A61F8Bk+luBZvAA4BpkTEmzNzfln228A7KD7Tx7Dws9fKd4A3lvs8DZhfVTAzfx8RRwEfofgMHwhFMA78BlgW2Csz763ahzTsMtNHDzwo/iNLYFLT8vPK5VcDq7XYbqWK5RMpmi9uaLEugfOalh1SLn8EeGXTul+X63ZrVbemZVPKsglMbVr3X+XyHzQt/2m5/GtNyzcGnirXHTLE6/i5svyJwJiG5ZOBB8p105u2mQAs22Jf21F8afyw4hxb1mlx3pMBzqf/2v+w8XzKdc8HVmp4fUZZ9vNN5V4PzKMIKp7XsHxqWX4esG3TNl8p1x041HMHprf6DFdtRxFAJPCxFuVXBJZvUdepDcueR9G/ZT7wxqbtP1OWP7OijncCL2lY/hyKwDWBzYf43hxSln8YWL9h+bLAdWW97gfe1LBuDHBWud0mTftbC4gWxzm8LL97xfGnVNSv/1znAJOH8p401P8KoA/YoVz287LsoUP97PrwMVwPm8bU74u58C/zBTLz4Yrlsyn+Ql4vIl7SxnG+m5nXNC37cfm8eRv7+WtmTm9a9jOKL90F+4mI5wJ7UHyZHNFYODP/ARzbxjEBPkDxH/iBmdnXsK/bge+22iAz52SLTFdmnknxhbZ9OxUYrvekzJTtTtFB+FON51Pu79Es+s4QRRPedhRf8F9vKncxRXZoFYrsV7PjM/PspmXTyud23vPF9UTzgsx8LDOftbzJLsCqwG8y88Kmdd+k+OPiLRXX+rBs6OOWmfMovuyh/XP+bmbe0LCvpyiyJ2OAUzPz/IZ1fcAvy5cbN+4kM2/LzFbTUXy7fG7rc9jg6+Xnf0jK+u8OPAYcGxGfoghELwAOW8w6SIvNQEj9Lq1aERFviIgTImJW2Ueiv6/Af5dFntW/ZwCXtVg2q3x+4ZLsJzOfoWjGadzPy4Hlgasz89EW+7loqAeMiOdTNG3MycxbWxQ5r2K7iIi9IuIvUfSpmddwDV9Je9evf5/D8Z5sRvF/wAWZ+dggZV9VPl9YXudm5zSVazRc73m7Tgb+DRxV9vHZJyI2jIgY4vavLp/PaV5RBjYXlC9H+pxb7au/M/nlLdb1N7M29z9bMSI+FxEzIuLhiOgrPzP9QXXbn8NS5f8dVTLznxQZ3BcB/1fWYc9c2DQnjRr7CKnf3a0WRsQ7KbIMT1Kk3G+l+EuujyL1/SaKVPdQPdRi2bzyeewS7qd/X4376e+cek9F+arlrQy2r5bXEDiSYvqCuRTNS3NYmKWYCry0jToM53uycvk8Z6BCpf5zrxpe3r985RbrHmpekJnzyniknfe8LZl5R0RsTtHEswMLs1WzIuIbmdkyg9dgWM+ZxfucQ+vRbPOGsG5BZ+SyY/I5FNmoaykySv8C+oPag2nv97hR1ed+MGdRNJW/APhtLuwnJ40qAyEBZUec1g6n6Dy7aWN6HiAifkR1B+Fu8Uj5vHrF+qrlrfR/6VRts0bzgrL56aMUXz6vb85KRcQebRy/33C9Jw+Vz0PJBPSf+7POsTSuqdxI6G+6a/X/1sqtNiivz+4R8RyKpqI3U2TNvhMRj2XmTwc4Xjec83DZhSIIOiYzpzauKDuVH9xqoyFqe/b3Mit3LEUQdB+wT0Qcn5kXDLylNPxsGtNg1gGub/GFOwbYsjNVasuNFNmXjcqmrWZDPocyiLkFmBARa7coMqXFsrUofs/ObBEETSzXN+tvHqjKHAzXe3IpRXCxVUSsOEjZK8vnLcugotnW5fMVbRy/XQ+Wz2u2WLfpQBtm5rzMvDwzv0bRZwyKUVED6T/nKc0rymvQf61H8pyHyzrl84kt1lUFzoN9DpfEpymydL+iGEX3DPDriFhtBI4lDchASIOZCawbEf0T3fX/NXcwxTDorpaZT1M0A6wEfKFxXURsTDF7dTt+TvF787Uy8Ojf12SKzE+zmeXzllHcPqK//PMoOom3Cir6h3RXdXieyTC8J5n5L4qh1+OAbzSeT38do5z3puyIfRbF8PaPN5XbAtiTIlD5w1CPvxj6+6J8uOn4rwQ+1lw4IjaPpnmQSv3LBpv9/I8UIwH3iGfPz/RxiiD2L7l0TPw5s3ye0rgwItYCvlaxzWCfw8VSfl6OoPijYr9y8MQnKDKT09vowyUNC5vGNJhvUcz9cmVEnEjxl9sbKL5w/wTs3MG6DdVBFH91Hlj+J3wxxZf/bhTznryDhc0ug/lmWX5X4IqIOIMiyNqdovPs2xsLZ+bdEXE88B7gqog4syz/Foo+PlcBmzQd4yaKfjvviYinKUZqJfCLzLyD4X1PDgBeQTHP0ZTyfJ6mmA5g+/J8zivL7gv8Ffi/iNiOohNv/zxCfcAHKjqkD5eTgH9SBCYTgb9TfEnvUq7bran8nsBHIuJ8ii/dBynmR9qZYtqEbw90sMz8d0T8J/Bb4PyI+C3Fe/EaihF0d1N0+F0a/IniGnyyDByvpLh2O1HM/9Mq2DmX4n39SkS8gjIjl5lHtCg7JBGxMkXwncB7+j8vmXl0RGwLvBv4JMXvmTQqzAhpQJn5I4oh43MpJrp7L8Xoly1YOpoEyMx7KOa6OZZiIsNPUIz02Z8iNQ8L+xINtq+nKPqZfItixMvHKP7KPqLcbysfBP6XYvTaRygCjFPKOj2rf0k5cuadFCPadqOYaPBwiuBkWN+TzHywrMcXKAKqfSgmC9yQYjqC6xvK3kbRBHU0xWi8T1FM2nc68IbMPKmdY7crM58EtgVOoAjeDqDIyuxJMQ9Ss+MoJtt8EcV1/DjFSLDjKfpX/W0IxzyJIsg8jeJ9+xTFhJZHA68pr0nXK0cFbkMxb9SGFNnLjSg+V3tVbHMDxefrborflcPLx5L4KUVW8aDMbB7x9iHgdorAazSmVZCAcnItqVdFxJcpJkncITPP6HR9JEmjy0BIPSEixmfmXU3LXsnCW15MKDMOkqQeYh8h9YrLIuIWimHsjwHrAm+jaB7e1yBIknqTGSH1hIg4mKKT8ySKe2g9BFwCfCMzz+tUvSRJnWUgJEmSepajxiRJUs9aqvsIPXPfbaazpA5YfvwbO10FqWfNe3rOqE46OZrftcusttaoT6hpRkiSJPUsAyFJktSzluqmMUmSNML65g9eZilmRkiSJPUsM0KSJKlaDvWe1EsnM0KSJKlnmRGSJEnV+swISZIk1ZIZIUmSVCntIyRJklRPZoQkSVI1+whJkiTVkxkhSZJUzT5CkiRJ9WQgJEmSepZNY5IkqZo3XZUkSaonM0KSJKmanaUlSZLqyYyQJEmq5oSKkiRJ9WRGSJIkVfKmq5IkSTVlRkiSJFWzj5AkSVI9mRGSJEnV7CMkSZJUT2aEJElSNe81JkmSVE9mhCRJUjX7CEmSJNWTgZAkSepZNo1JkqRqTqgoSZJUT2aEJElSNTtLS5Ik1ZMZIUmSVM0+QpIkSfVkRkiSJFXK9BYbkiRJtWRGSJIkVXPUmCRJUj2ZEZIkSdUcNSZJklRPZoQkSVI1+whJkiTVkxkhSZJUrc95hCRJkmrJQEiSJPUsm8YkSVI1O0tLkiTVkxkhSZJUzQkVJUmS6smMkCRJqmYfIUmSpHoyIyRJkqrZR0iSJKmezAhJkqRqZoQkSZLqyYyQJEmqlOlNVyVJkmrJjJAkSapmHyFJkqR6MiMkSZKqObO0JElSPRkISZKknmXTmCRJqmZnaUmSpHoyIyRJkqrZWVqSJKmezAhJkqRq9hGSJEmqJzNCkiSpmn2EJEmS6smMkCRJqmYfIUmSpHoyIyRJkqqZEZIkSaonM0KSJKmao8YkSZLqyYyQJEmqZh8hSZKkejIQkiRJPcumMUmSVM3O0pIkSfVkRkiSJFWzs7QkSVI9mRGSJEnV7CMkSZJUT2aEJElSNfsISZIk1ZMZIUmSVM2MkCRJUj2ZEZIkSdUyO12DEWVGSJIk9SwzQpIkqZp9hCRJkurJjJAkSapmRkiSJKmezAhJkqRq3mtMkiSpngyEJElSz7JpTJIkVbOztCRJUj0ZCEmSpGqZo/cYRETsEBE3RcQtEXFQi/UrRcSfIuIfEXFdRHxgsH0aCEmSpK4XEWOBo4AdgQ2APSJig6ZiHwGuz8yNgSnANyPiuQPt1z5CkiSpWvf0EdocuCUzbwOIiOOBXYDrG8ok8PyICOB5wAPAvIF2akZIkiQtDSYAsxpezy6XNfo+sD5wF3AN8LHMgSdCMiMkSZKqjWJGKCL2AfZpWDQtM6f1r26xSXPHou2Bq4BtgLWBsyLiwsx8pOqYBkKSJKkrlEHPtIrVs4E1G15PpMj8NPoA8NXMTOCWiLgdWA+4tOqYNo1JkqRq2Td6j4HNANaNiMllB+j3ACc3lbkT2BYgIlYHXg7cNtBOzQhJkqSul5nzIuIA4AxgLPCzzLwuIvYt1x8NHA5Mj4hrKJrSPpOZ9w20XwMhSZJUKfsGn99ntGTmacBpTcuObvj5LmC7dvZp05gkSepZZoQkSVK17plHaESYEZIkST3LjJAkSao2+GiupZoZIUmS1LMMhCRJUs+yaUySJFXrouHzI8GMkCRJ6llmhCRJUjWHz0uSJNWTGSFJklTNjJAkSVI9mRGSJEnV0lFjkiRJtWRGSJIkVbOPkCRJUj2ZEZIkSdWcWVoaXl/43yPZ6m3v4R177dvpqkhLpe23m8J1117AjddfxIGf/kjLMt868jBuvP4irrj8LF61ySuGvO0nP/FfzHt6Dquu+kIANtt0Ey6bcSaXzTiTyy87i1122WFkTkrqEAMhjbp3vPUtHH3kEZ2uhrRUGjNmDN/9zpfZaee9eOXGW7P77u9g/fXXXaTMjjtsw7rrTGa9DbZkv/0+w1Hf/8qQtp04cTxv3nYr7rhj9oJl1153I1u8dkc23Ww73rbTe/nhUV9j7Nixo3Oy6g7ZN3qPDuh4IBQRyw5lmepj001eyUoveH6nqyEtlTbf7FXceutMbr/9Tp555hlOOOEk3r7z9ouU2Xnn7fnFr34HwN8vvYKVVl6JNdZ48aDbfvMbh3DQ575MNgyXfuKJJ5k/fz4Ayy237CLrpDroeCAE/G2IyySp542fsAazZt+14PXsOXMZP36NRcpMGL8Gs2ctLDNn9lwmjF9jwG132uktzJkzl6uvvv5Zx9x8s1fxj6vO4aorzmb/Aw5aEBipR/Tl6D06oGOdpSNiDWACsHxEvAqIctULgBU6VS9J6mYR8axlzVmaqjJVy5dffjk+d9BH2eGte7Y85qUzrmTjTbZhvfXW4ec//Tann34uTz311GKegdRdOpkR2h74BjAROBL4Zvn4JPC5qo0iYp+IuCwiLvvJsceNSkUlqVvMmT2XNSeOX/B64oRxzJ17zyJlZs+Zy8Q1F5aZMHEcd829p3LbtdeexKRJL+GKy87ilpsvYeLEccz4+xmsvvqLFtnvjTfewmOPPcErNnz5CJ2dNPo6lhHKzGOAYyJi18w8sY3tpgHTAJ657zYbqyX1lBmXXcU660xm0qQ1mTPnbnbbbRf2ft+io79OOeVM9t9vKr/5zUlssfmreeThR7j77nv517/ub7nt9dffzPiJGy/Y/pabL2GL1+3I/fc/yKRJazJr1l3Mnz+fl7xkAi972VrMvGPWaJ+2OihrPqFiN8wjdHZEHAlsVb4+HzgsMx/uYJ00gj598FeZceXVPPTQI2z7jr3Y/4N7s2tTZ09Jrc2fP5+PffwLnHbqrxk7ZgzTj/kN119/M/t8eG8Apv34F5z257PZYYdtuOmGv/L4E0/woQ99csBtB/KGN2zOgZ/+CM88M4++vj4O+OjnuP/+B0f8PKXREp0eARARJwLXAseUi/YGNs7Mdw22rRkhqTOWH//GTldB6lnznp7z7M5eI+ixL79v1L5rV/z8saN6btAdGaG1M3PXhteHRsRVnaqMJEnqHd0QCD0REVtm5kUAEfEG4IkO10mSJEHHJjocLd0QCO1H0Wl6pfL1g8D7O1gfSZLUI7ohELoB+DqwNrAy8DDwDuDqzlVJkiQBtb/pajcEQicBDwFXAHM6WxVJktRLuiEQmpiZ3s5YkqRuVPN5hLrhXmMXR8QrO10JSZLUe7ohI7QlMDUibgeeorjnWGbmRp2tliRJso/QyNux0xWQJEm9qeOBUGbe0ek6SJKkCjWfR6gb+ghJkiR1RMczQpIkqYvVvI+QGSFJktSzDIQkSVLPsmlMkiRVSidUlCRJqiczQpIkqZqdpSVJkurJjJAkSapmRkiSJKmezAhJkqRq3mJDkiSpnswISZKkavYRkiRJqiczQpIkqVKaEZIkSaonM0KSJKmaGSFJkqR6MiMkSZKqefd5SZKkejIQkiRJPcumMUmSVM3O0pIkSfVkRkiSJFUzIyRJklRPZoQkSVKlTDNCkiRJtWRGSJIkVbOPkCRJUj2ZEZIkSdXMCEmSJNWTGSFJklQpzQhJkiTVkxkhSZJUzYyQJElSPZkRkiRJ1fo6XYGRZUZIkiT1LAMhSZLUs2wakyRJlRw+L0mSVFNmhCRJUjUzQpIkSfVkRkiSJFVz+LwkSVI9mRGSJEmVHDUmSZJUU2aEJElSNfsISZIk1ZMZIUmSVMk+QpIkSTVlRkiSJFWzj5AkSVI9mRGSJEmV0oyQJElSPRkISZKknmXTmCRJqmbTmCRJUj2ZEZIkSZXsLC1JklRTZoQkSVI1M0KSJEn1ZEZIkiRVso+QJElSTZkRkiRJlcwISZIkdYGI2CEiboqIWyLioIoyUyLiqoi4LiLOH2yfZoQkSVKlbskIRcRY4CjgLcBsYEZEnJyZ1zeUWRn4AbBDZt4ZES8ebL9mhCRJ0tJgc+CWzLwtM58Gjgd2aSqzJ/D7zLwTIDPvHWynBkKSJKlaxug9BjYBmNXwena5rNHLgBdGxHkRcXlEvG+wndo0JkmSukJE7APs07BoWmZO61/dYpNsev0c4DXAtsDywN8i4pLMvLnqmAZCkiSp0mj2ESqDnmkVq2cDaza8ngjc1aLMfZn5GPBYRFwAbAxUBkI2jUmSpKXBDGDdiJgcEc8F3gOc3FTmJOCNEfGciFgB2AK4YaCdmhGSJEldLzPnRcQBwBnAWOBnmXldROxbrj86M2+IiNOBqynukvaTzLx2oP1GZnPz2tLjmftuW3orLy3Flh//xk5XQepZ856eM2iv4uE0d8utR+27dtxF547quYFNY5IkqYfZNCZJkip1y4SKI8WMkCRJ6llmhCRJUqUcfKLDpZoZIUmS1LPMCEmSpEr2EZIkSaopM0KSJKlS9tlHSJIkqZbMCEmSpEpL8Q0ohsSMkCRJ6llmhCRJUiX7CEmSJNWUGSFJklTJjJAkSVJNGQhJkqSeZdOYJEmq5PB5SZKkmjIjJEmSKtlZWpIkqabMCEmSpEqZZoQkSZJqyYyQJEmqlH2drsHIGjAQiog+YHEGzmVmGmRJkqSuNliwcgGLFwhJkqQa6Kt5H6EBA6HMnDJK9ZAkSRp1Nl9JkqRKjhqTJEmqqcXKCEXEOGBbYAKwbIsimZmHL0nFJElS59V9Zum2A6GIOBQ4qGnbYGGn6v6fDYQkSVJXa6tpLCLeC3wRuBB4N0XQcwywJ/BjoA84HthmeKspSZI6IXP0Hp3QbkZoP2A2sENmzosIgJmZeTxwfET8ATgVOG54qylJkjT82u0s/UrgtMyc17BsbP8PmXkGcAbw6WGomyRJ0ohqNyO0DHB/w+sngJWaylwL7LsklZIkSd2h7p2l280IzQXGNby+E9ioqcwEYB6SJEldrt1A6EqK5rF+5wBvjIi9I2LFiHgbsGtZTpIkLeX6Mkbt0QntBkKnABtGxOTy9VeBh4HpwCPAyRQjyb4wXBWUJEkaKW31EcrM6RRBT//rWRGxGfA/wNrATOAHmXnN8FVRkiR1St1vsbHE9xrLzNuBA4ahLpIkSaPKm65KkqRKnZrocLS0FQhFxEuGWjYz72y/OpIkSaOn3YzQTBbeU2wguRj7liRJXaZTo7lGS7vByrG0DoRWBjYBXgqcB9yxJJWSJEkaDe2OGptatS4ixlDckHVf4P1LVi1JktQN6j5qrN15hCplZl9mHkrRfPbV4dqvJEnSSBmJfjwXA+8bgf1KkqRRVvdRY8OWEWqwCrDiCOxXkiRpWA1rRigi3gzsTnEHekmStJRz1FiDiDhngP2sCfTPM3TYklRKkiRpNLSbEZpSsTyBB4EzgG9kZlXANKxWnLDVaBxGUpPH/vHLTldB0iip+6ixdofPj0SfIkmSpI4wsJEkST2rrUAoIs6JiAGHxkfEXgP0JZIkSUuRvoxRe3RCuxmhKcCkQcq8FHjT4lRGkiRpNI3EhIrLA/NGYL+SJGmU1Xw+xcUKhFpek4gIiuHzbwVmLUmlJEmSRsOggVBE9LFo8HNIRBwy0CbA/y5hvSRJUhdwQkW4gIWB0FbAnRQ3Vm02H7gfOBv4yXBUTpIkaSQNGghl5pT+n8vs0M8z05mjJUnqAU6ouKjJwEMjUA9JkqRR124gdC/wooh4IjOfbl4ZEcsCqwP3ZuaTw1FBSZLUOX2drsAIa3ceoS8BNwHPq1i/InAj8LklqZQkSdJoaDcQ2hH4S2Y+0GplufwvwE5LWjFJktR5SYzaoxPaDYQmATcPUuZmBp99WpIkqePa7SO0DIM3Fyaw3OJVR5IkdZO+mk8t3W5G6DYGv4/YFOCOxaqNJEnSKGo3EDoZeE1EHNhqZUQcBLwa+OMS1kuSJHWBPmLUHp3QbtPYN4D3Al+JiN2AM4E5wARge2ATipmnvz6MdZQkSRoRbQVCmflgREwBfgW8jiL7k7AgjLsY2CszHxzGOkqSJI2Itu8+n5kzgTdExKuB1wIrU8w2fUlmXjGclZMkSZ3VqWHto6XtQKhfGfQY+EiSpKXWYgVCETEO2Jaib9CyLYpkZh6+JBWTJEmdV/dbbLQdCEXEocBBTdsGRV+hxp8NhCRJUldra/h8RLwX+CJwIfBuiqDnGGBP4McUgePxwDbDW01JktQJdb/FRrsZof2A2cAOmTkvIgBmZubxwPER8QfgVOC44a2mJEnS8Gt3QsVXAqdl5ryGZWP7f8jMM4AzgE8PQ90kSVKH9Y3ioxPaDYSWAe5veP0EsFJTmWuBjZekUpIkSaOh3aaxucC4htd3Ahs1lZkAzEOSJC316j5qrN2M0JUUzWP9zgHeGBF7R8SKEfE2YNeynCRJUldrNxA6BdgwIiaXr78KPAxMBx6huClrAF8YrgpKkqTOcdRYg8ycThH09L+eFRGbAf8DrA3MBH6QmdcMXxUlSZJGxmLfYqNfZt4OHDAMdZEkSV2mr963Gmu7aUySJKk2ljgjJEmS6quv5nefNyMkSZJ6loGQJEnqWTaNSZKkStnpCowwM0KSJKlnmRGSJEmVvMWGJElSTZkRkiRJlfrC4fOSJEm1ZEZIkiRVctSYJElSTZkRkiRJlRw1JkmSVFNmhCRJUqW+eg8aMyMkSZJ6lxkhSZJUqY96p4TMCEmSpJ5lRkiSJFVyHiFJkqQuEBE7RMRNEXFLRBw0QLnNImJ+RLx7sH0aCEmSpK4XEWOBo4AdgQ2APSJig4pyXwPOGMp+bRqTJEmVumj4/ObALZl5G0BEHA/sAlzfVO6/gROBzYayUzNCkiRpaTABmNXwena5bIGImAC8Ezh6qDs1IyRJkiqN5i02ImIfYJ+GRdMyc1r/6habNPfl/jbwmcycHzG0VJaBkCRJ6gpl0DOtYvVsYM2G1xOBu5rKbAocXwZBqwFvjYh5mfnHqmMaCEmSpEpdNHx+BrBuREwG5gDvAfZsLJCZk/t/jojpwCkDBUFgICRJkpYCmTkvIg6gGA02FvhZZl4XEfuW64fcL6iRgZAkSarURaPGyMzTgNOalrUMgDJz6lD26agxSZLUs8wISZKkSqM5aqwTzAhJkqSeZUZIkiRVMiMkSZJUU2aEJElSpeyiUWMjwYyQJEnqWWaEJElSJfsISZIk1ZSBkCRJ6lk2jUmSpEo2jUmSJNWUGSFJklQpO12BEWZGSJIk9SwzQpIkqVKfEypKkiTVkxkhSZJUyVFjkiRJNWVGSJIkVTIjJEmSVFNmhCRJUiXnEZIkSaopM0KSJKmS8whJkiTVlBkhSZJUyVFjkiRJNWUgJEmSepZNY5IkqZLD5yVJkmrKjJAkSarUV/OckBkhSZLUs8wISZKkSg6flyRJqikzQpIkqVK9ewiZEZIkST3MjJAkSapkHyFJkqSaMiMkSZIq9UWnazCyzAhJkqSeZUZIkiRVcmZpSZKkmjIjJEmSKtU7H2RGSJIk9TADIUmS1LNsGpMkSZWcUFGSJKmmzAhJkqRKDp+XJEmqKTNCkiSpUr3zQWaEJElSDzMjJEmSKjlqTJIkqabMCEmSpEqOGpMkSaopM0KSJKlSvfNBZoQkSVIPMyMkSZIqOWpMkiSppswISZKkSlnzXkJmhCRJUs8yEJIkST3LpjFJklTJztKSJEk1ZUZIkiRV8hYbkiRJNWVGSJIkVap3PsiMkCRJ6mFmhCRJUiX7CEmSJNWUGSFJklTJeYSkIdhuuylce835XH/9RXz6Ux9pWebIIw/j+usv4vLLzmKTTV6xYPm0H32D2bOu4sor/rJI+UMO/hSXX3YWMy49g1NP/RXjxq0+oucgLe0uuuI6dt7/YN627xf56YmnP2v9I/9+jI9/5Yfs+rHD2fPTX+Gfd8xZZP38+X3s9okvc8ARR41WlaWOMxDSEhszZgzf+c4R7Pz2vdl4463ZffddWH+9dRcps8MO27DOOpPZYIMt2W//z/D9731lwbpjf/Fbdtp5r2ft95tHHs1rNn0Lm22+Paeddjaf//zHR/pUpKXW/Pl9/O+PjuOHXzqAP37vYP584QxunXXXImV+/LvTefnkNTnxO1/kyx/7AF/7yQmLrP/VKecweeIao1ltLQVyFP91goGQlthmm23CrbfO5Pbb7+SZZ57hhBNOYuedt1ukzM47b8evfvk7AC699ApWXvkFrLHGiwG46KK/8+CDDz1rv48++u8FP6+4wvJk1rvDnrQkrv3nTF4y7sVMXONFLLPMc9hhy8049+9XL1Lmtllz2WKj9QCYPHEN7rr3fu5/6BEA7r7vQS647Bre9ZY3jHrdpU7qeCAUEWtFxJ8i4r6IuDciToqItTpdLw3dhPHjmD1r7oLXc+bczfgJ4xYpM378GsyavfCv09lz5jJ+/OB/eR526IHcesul7LHHOzn00G8MX6WlmrnngQdZfbUXLni9+qorc+8DDy5S5mWTJnL2JVcCcM3NtzP3Xw9wz31Fma//9AQ++f53MSZi9CqtpULfKD46oeOBEPBr4ARgDWA88FvguI7WSG1p9f9mc/YmWhQaSobnSwd/nbXX2ZzjjvsD++/3gcWuo1R7LX6dgkV/7z646/Y88u/H+Y+PH8Fxp57HemutydixYzl/xtWsstLz2WCdl45SZaXu0Q2jxiIzf9Hw+pcRcUBl4Yh9gH0Axo5dmTFjVxzp+mkQs+fMZeKaCzNAEyaswdy77l6kzJw5c1lz4vgFrydOGMfcufcM+RjH/+aPnPTHYzjs8G8ueYWlGlp91RcuyO4A3HP/Q7xolZUXKfO8FZbn8I++Hyj+ENlxn88zYfVVOf2iGZw342ouuvxannpmHo89/gSf/dbP+Mon/nM0T0FdqlN9d0ZLN2SEzo2IgyJiUkS8NCIOBE6NiFUiYpXmwpk5LTM3zcxNDYK6w2WX/YN11pnMpElrsswyy7DbbrtwyilnLVLmlFPO5L17vRuAzTd/NQ8//Ch3333vgPtdZ53JC37eaaftuOmmW4e/8lJNbLjuS7lj7r3Mvuc+nnlmHqdfNIMpm2+0SJlH/v04zzwzD4ATz7qIV2+4Ls9bYXk+tvc7+ctPv8rpP/5fvv4/H2TzjdYzCFLP6IaM0O7l8381Lf9PimSv/YW63Pz58/n4x7/Iqaf8ijFjx3DM9N9w/Q038+EPFyPBfvzjX/LnP5/DDjtsww03XMQTjz/Jhz78yQXb/+LY77PVVq9jtdVW4bZbZ3DY4d9k+vTj+fIRn+VlL1uLvr7kzjtn85EDPtupU5S63nPGjuVzH96d/Q79LvPn9/GON7+edV4ynhNOvwCA3XbYittn383nv/NzxowZw9prjuPQA/bucK2lzouleSTOc5eduPRWXlqKPXrVLwYvJGlELLv+1qPao/39k3Ydte/aY2aeOOq99TueEYqI97VanpnHjnZdJElSb+l4IARs1vDzcsC2wBWAgZAkSR3WtxS3HA1FxwOhzPzvxtcRsRJg3l2SJI24jgdCLTwOrDtoKUmSNOLqnQ/qgkAoIv7Ewus8FlifYoJFSZKkEdXxQAhovG/CPOCOzJzdqcpIkqSF+mqeE+p4IJSZ5/f/HBE7ZeZfO1kfSZLUO7phZulGh3W6ApIkaaEcxX+d0G2BkLc9liRJo6bjTWMRsWxmPlW+/K8WyyRJUof0dboCI6wbMkJ/6/8hMy9tXiZJkjRSOpYRiog1gAnA8hHxKhY2i70AWKFT9ZIkSQs5amzkbA9MBSYCRzYsfxT4XCcqJEmSekvHAqHMPAY4JiJ2zcwTO1UPSZJUrVOjuUZLJ5vGPtnq536ZeWTzMkmSpOHUyaax53fw2JIkSR1tGju0U8eWJElDU/fh851sGjswM78eEd+jxc1tM/OjHaiWJEnqIZ1sGruhfL6sg3WQJEkDyLSz9IjIzD+Vz8d0qg6SJKm3dcMtNs6lddPYNh2ojiRJauCEiiPvUw0/LwfsCszrUF0kSVKXiogdgO8AY4GfZOZXm9a/F/hM+fLfwH6Z+Y+B9tnxQCgzL29a9NeIOL8jlZEkSYvollFjETEWOAp4CzAbmBERJ2fm9Q3FbgfelJkPRsSOwDRgi4H22/FAKCJWaXg5BtgUWKND1ZEkSd1pc+CWzLwNICKOB3YBFgRCmXlxQ/lLKG7jNaCOB0LA5SzsIzQPmAl8sGO1kSRJC3TRLTYmALMaXs9m4GzPB4E/D7bTbgiENgD2B7akCIguxCH1kiT1nIjYB9inYdG0zJzWv7rFJi2jtIjYmiIQ2nKwY3ZDIHQM8Ajw3fL1HsAvgP/oWI0kSRIwuqPGyqBnWsXq2cCaDa8nAnc1F4qIjYCfADtm5v2DHbMbAqGXZ+bGDa/PjYgBe3hLkqSeMwNYNyImA3OA9wB7NhaIiJcAvwf2zsybh7LTbgiEroyI12bmJQARsQXw1w7XSZIk0T0zS2fmvIg4ADiDYvj8zzLzuojYt1x/NPAlYFXgBxEBMC8zNx1ov52819g1FG17ywDvi4g7y9cvpaEHuCRJEkBmngac1rTs6IafPwR8qJ19djIjtFMHjy1JkoagW+YRGimdvNfYHZ06tiRJEnRHHyFJktSlumgeoRExptMVkCRJ6hQDIUmS1LNsGpMkSZVGc0LFTjAjJEmSepYZIUmSVKlbJlQcKWaEJElSzzIjJEmSKtlHSJIkqabMCEmSpEpOqChJklRTZoQkSVKlPkeNSZIk1ZMZIUmSVKne+SAzQpIkqYeZEZIkSZWcR0iSJKmmzAhJkqRKZoQkSZJqykBIkiT1LJvGJElSpXRCRUmSpHoyIyRJkirZWVqSJKmmzAhJkqRKaUZIkiSpnswISZKkSo4akyRJqikzQpIkqZKjxiRJkmrKjJAkSapkHyFJkqSaMiMkSZIq2UdIkiSppswISZKkSs4sLUmSVFMGQpIkqWfZNCZJkir1OXxekiSpnswISZKkSnaWliRJqikzQpIkqZJ9hCRJkmrKjJAkSapkHyFJkqSaMiMkSZIq2UdIkiSppswISZKkSvYRkiRJqikzQpIkqZJ9hCRJkmrKjJAkSapkHyFJkqSaMhCSJEk9y6YxSZJUKbOv01UYUWaEJElSzzIjJEmSKvXZWVqSJKmezAhJkqRK6YSKkiRJ9WRGSJIkVbKPkCRJUk2ZEZIkSZXsIyRJklRTZoQkSVKlPjNCkiRJ9WRGSJIkVUpHjUmSJNWTGSFJklTJUWOSJEk1ZSAkSZJ6lk1jkiSpkrfYkCRJqikzQpIkqZKdpSVJkmrKjJAkSarkLTYkSZJqyoyQJEmqZB8hSZKkmjIjJEmSKjmPkCRJUk2ZEZIkSZXsIyRJklRTZoQkSVIl5xGSJEmqKTNCkiSpUjpqTJIkqZ4MhCRJUs+yaUySJFWys7QkSVJNmRGSJEmVnFBRkiSppswISZKkSg6flyRJqikzQpIkqZJ9hCRJkmrKjJAkSapkRkiSJKmmzAhJkqRK9c4HmRGSJEk9LOre9qfuFRH7ZOa0TtdD6jX+7kkLmRFSJ+3T6QpIPcrfPalkICRJknqWgZAkSepZBkLqJPsoSJ3h755UsrO0JEnqWWaEJElSzzIQUleJiKkRMb7T9ZC6UURMiohr2yj/9og4qPz5kIj4VPmzv2dSyUBI3WYq4H/Q0jDIzJMz86stVk3F3zMJMBDSCCv/gr0hIn4cEddFxJkRsXxEbBIRl0TE1RHxh4h4YUS8G9gU+FVEXBURy3e6/lIXek5EHFP+7vwuIlaIiJkRsRpARGwaEeeVP0+NiO83bjzY71lETImIUxpefz8ippY/z4yIr0XEpeVjnZE9VWnkGQhpNKwLHJWZGwIPAbsCxwKfycyNgGuAgzPzd8BlwHszc5PMfKJTFZa62MuBaeXvziPA/u1sPAy/Z49k5ubA94Fvt7mt1HUMhDQabs/Mq8qfLwfWBlbOzPPLZccAW3WiYtJSaFZm/rX8+ZfAlqN8/OManl83yseWhp2BkEbDUw0/zwdW7lA9pDponvMkgXks/P98uXZ2FhFblE1kV0XE25v21Wp/WfGztFQyEFInPAw8GBFvLF/vDfRnhx4Fnt+RWklLh5dERH8mZg/gImAm8Jpy2a5D2MeC37PM/HvZRLZJZp4M3AFsEBHLRsRKwLZN2+7e8Py3xT8NqTs8p9MVUM96P3B0RKwA3AZ8oFw+vVz+BPA6+wlJz3ID8P6I+BHwT+CHwKXATyPic8Dfh7CP6VT8nmXmrIg4Abi63P+VTdsuGxF/p/hDeo8lPRmp05xZWpI0JBExE9g0M+/rdF2k4WLTmCRJ6llmhCRJUs8yIyRJknqWgZAkSepZBkKSJKlnGQhJoywisv9eUEuwj0nlfqYPT606r9V1Ke+YnhExZYSOWbvrKKk9BkKSam04Ak9J9eWEipK62feB44E7R2j/c4D1KWY7l9SDDIQkda1y4r4Rm7wvM58Bbhyp/UvqfjaNqXYa+31ExNoR8buIuD8iHo2IMyPiFWW5F0XEtIiYGxFPRsSMiNi6Yp8rRcRXIuKmsuyDEXFGRLy5ovxzI+KLEXFrRDwVEbdHxBERsewA9X5OROwfEZdExCMR8XhEXBkRB0TEEv2uRsSU8pocEhGvi4i/RMTD5TU5IyI2bbHNgv45EbFnRPw9Iv5dzi7cX2aFiPhsecPOx8r1f4uIlrdeaPe6DNRHKCLWi4ifRcTMcl/3RsSFEbFfuX5qRPRPlPamcj/9j0PKMpV9hCJiXEQcVe7/6Yj4V0T8PiJe06Ls1HI/UyNi64g4r7y2j0TEqRGxfottVo+Ib5Sfqcci4qHy5+kRsVar6yFp+JkRUp1Norjv0g0U91aaBLwTOC+Km1aeDjwC/AZYBXgP8OeIeFlmLmiKiYiVgb8CGwAzgG8DqwG7AWdGxH6Z+aOG8gGcAOwC3ErRvPNc4D+BV7aqaEQsA/wJ2B64Cfg18CSwNfA9YAuKm9MuqS2AzwJ/AY4C1gHeBWwVEdtl5oUttvkf4C1l/c4FVirrvDJwDvAq4ArgZxR/XG0P/DoiNszMLzScY9vXpUpEvA34LbAsxft4HLAysDFwIMX9t64CDgUOpriR6PSGXZw3yP4nU9zMdHx5jscBawL/AbwtInbNzFNabLpTeX5/Bo6m+My8FdgsIjbovzVFFPfY+yuwNnAWxbUN4KXl9r+juAefpJGWmT581OpBEfBk+fh807ovlssfoPiiGtOwbu9y3beatvlRufxHlLOxl8vXpehb8hQwqWH5nmX5vwHLNSxfhSIASOC8pmMcUi7/HjC2YflY4Kflul1anOP0IV6TKQ3X5ICmdbuUy//ZdD366/QY8KoW+5xerj+waflyFMFJH7DJMF2XKQ3LViuv+9PAm1rUa2LT62ftd7DrCJxR8fl5PTAPuB94XsPyqWX5ecC2Tdt8pfk6ATu3+qyV654LPL/Tv0c+fPTKw6Yx1dlM4KtNy44pn5cFPp2ZfQ3rfk3xRbZJ/4IyU7MX8G/gs5m54J40mflP4LsUX1zva9jPB8rnz2Xmkw3lHwAOb65k2ex1AHA38InMnN+wzXyKjEwC7x3shIfgFuAHjQsy8yTgfIrs0BtbbDMtMxe5A3lErEpxXS7LzK837e9J4DMUGY49G1a1dV0G8H7gBcAPM/P85pWZObuNfT1LREwEtqPooN18bhdTZIdWocikNTs+M89uWjatfN68Rfknmhdk5tOZ+Wi79Za0eGwaU51d1RhUlO4qn29u/rLJzPkRcQ8wsWHxesAKwF/LL+xm5wBfoGge6vdqimzIRS3Kn9di2cuAVSkyMl8oWpCe5QmK0U1L6sKm4K+xXm+iOI/m4OLSFuU3o8hWLehv02SZ8rmxzu1elyqvLZ//3MY27eh/Ly/MojN1s3MogsBXAcc2rbusRflZ5fMLG5adTzFi7aCIeDVwGkVTWavPrKQRZCCkOnvWkOjMnFcGGlXDpeex8Escyv4wwNyK8v3LV27a5oGKL9G7WyxbtXxel6I/S5XnDbBuqO6pWN5fr5UGWNeov86blY8qjXVu97pUWbl8ntPGNu1YnPe830PNCxo+c2Mblj0SEa+l6MP0dop+VQD3RcQPgCMqrpOkYWbTmDSw/oBpjYr145rK9f+8Stms1qzVfvq3/UNmxgCPye1X/1lWr1jeX69WAWK2WNZf7luD1Hnrpm3auS5VHiqfJ7SxTTsW5z1vW2bOzswPAi8GXgF8lKLv0ZfKh6RRYCAkDewm4HFgk4h4YYv1/V/0VzQsu4Lid2vLFuWntFh2I8WX+2srgoThtGXFUPwp5fOVLda1cilFM1erPkVV2r0uVS4pn3ccYvk+GrIxQ9B/DbaMiFZZ81bv+WLLwnWZ+T2K0XkA7xiOfUsanIGQNIDMfBr4FUUTz2GN6yJibYq/4p8BftGw6ufl85cjYrmG8qtQ9CdqPsY8itFi44DvRsTyzWXKOW02WLKzAYrmt/2b9r0LRf+gW4BWw+efJTPvpbgum5bzAj0rYIhiDqfGLFZb12UAx1BMe7BfRGzV4rgTmxbdTzH0fUjKztZnUYwo+3jTvreg6AD+IPCHNurcXMdXRMSkFqv6M3aPL+6+JbXHPkLS4A6iyHwcEBGbUcyl0z+P0PMphqPf3lD+OGB3ir4f10bESRT9jt5NMQ/R2i2OcTjFHDj7AjtHxDkUfWBeTBG8vAH4PHD9Ep7L6cA3I2JH4B8snEfoSeCDFR2pqxxQ1u0wYO+IuIiiD9J4ik7SmwF7AP3XZnGuy7Nk5n0RsSfFXDvnRsSfgaspRpJtRBH0NAZgZwPviYg/AZdT9AO7IDMvGOAw+1J0Xv6/iNiOohN0/zxCfcAHlnBk15uBIyPiYoqM4L0UnfR3Kff/f0uwb0ltMBCSBpGZD5QTMH6WImj4JMUorkuB/8vMM5vKZ0T8B0UANZUiYJhLkRE5jCLoaD7GMxHxDorRSFMpJuZ7HvAvikDiixQZmCX197IOh5f1CopRUJ/PzBnt7Kjs8PsmYB+KLMmuFHMI3UMxAu4TFJmV/vJtX5cBjn1qFLNhfwbYlmK4+4MUQcVXmop/jKKf07YUkxuOoeikXBkIZeZt5f6/UG4zhSILdTrw5XavVQtnUEzMuRVF8PMCimtxFnBkOUxf0iiIhmlRJNVUeYuKc4FDM/OQjlZGkrqIfYQkSVLPMhCSJEk9y0BIkiT1LPsISZKknmVGSJIk9SwDIUmS1LMMhCRJUs8yEJIkST3LQEiSJPUsAyFJktSz/h9fGfWwmjWTNQAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "fig, ax = plt.subplots(1, figsize=(10,10))\n", "sns.heatmap(trainCM/trainCM.sum(axis=1), annot=True)\n", "ax.set_xlabel('model predictions', fontsize=20)\n", "ax.set_ylabel('actual', fontsize=20)\n", "plt.title(\"Training data confusion matrix\", fontsize=20);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The matrix shows True Negatives (top left), False Positives (top right), True Positives (bottom right) and False Negatives (bottom left).\n", "\n", "Accuracy measures the ratio of hits for both postiive (built up) and negative (not) classes:\n", "\n", "$$ACC = \\frac{TP + TN}{TP + TN + FP + FN}$$" ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our classifier has an accuracy of 0.98884 on the training data.\n" ] } ], "source": [ "acc = (trainCM.loc['built-up','built-up'] + trainCM.loc['not','not']) / trainCM.sum().sum()\n", "\n", "print(f\"Our classifier has an accuracy of {acc:.5f} on the training data.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With a balanced class (~50% built up and 50% not), this metric might be ok, assuming we didnt care at all about false positives or negatives. \n", "\n", "But accuracly is not a great metric.\n", "\n", "Consider an extreme (but realistic) example where 1 out of 100 observations are in the positive class. A classifier could predict that everything is negative and be right 99 of 100 times...99% accuracy might be considered good, but there is no predictive power here and no information is yielded. \n", "\n", "If you do care about false negatives (say, if the test is for a disease or something with high stakes) you do not want a classifier that just predicts negatives naively: 1 miss in 100 could be terrible for real world consequences.\n", "\n", "For binary classification there are many options. The Matthews correlation coefficient (based on Pearson's phi coefficient) is a better score, so let's use that.\n", "\n", "$$MCC = \\frac{TP * TN - FP * FN}{\\sqrt{(TP + FP)(TP + FN)(TN + FP)(TN + FN)}}$$" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our classifier has an MCC score of 0.96270 on the training data.\n" ] } ], "source": [ "def get_mcc(cm, pos_label='built-up',neg_label='not'):\n", " tp = cm.loc[pos_label, pos_label]\n", " tn = cm.loc[neg_label, neg_label]\n", " fn = cm.loc[pos_label, neg_label]\n", " fp = cm.loc[neg_label, pos_label]\n", " \n", " return (tp * tn - fp * fn)/(np.sqrt((tp + fp)*(tp + fn) * (tn + fp)* (tn +fn)))\n", "\n", "mcc = get_mcc(trainCM)\n", "print(f\"Our classifier has an MCC score of {mcc:.5f} on the training data.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Test results\n", "\n", "Either way, the performance seems good, but remember this is the data we trained the classifier on, so it really just learned about this specific data.\n", "\n", "**Never base your model evaluation on your training data alone!**\n", "\n", "We need to look at the out sample that it was not trained on to get a better indication of how this classifier will do on unseen data, which is what we are relying on to get accurate predictions of 2016 through 2019 land use classification." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our classifier has an accuracy of 0.89106 on the test data.\n", "Our classifier has an MCC score of 0.63075 on the test data.\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAJrCAYAAAAMBPz+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8LElEQVR4nO3deZwcZZ348c83IRyC3MqRRG4EFEHl8IQoKqAgKiqHInHXjbjiuT8VrxVE1/tEEKOLAVZBVmQBCYKC3KBECcghGkIgCeG+Ecwx398fVZN0Ol0z08lMd0/3551Xv7qr6qmnnqqZyXzn+zz1VGQmkiRJvWhMuxsgSZLULgZCkiSpZxkISZKknmUgJEmSepaBkCRJ6lkGQpIkqWcZCKmnRcSciJjT7naMZhGxbkR8v7yWiyMiI2LXET7mluVxpo3kcXpNeU0va3c7pFYyENKgyv8cm3lNHoE2TB6puodDj/8C+TrwIeAvwFeA44B729qiHhQRk8rvw2Pb3RZpNFmt3Q3QqHBcg3UfBdYDvgc8Wrdt5sg2Rx3mAOBvmXlgC485H9gReKyFx+wFOwL/aHcjpFYyENKgMvPY+nVlZmY94LuZOafFTVJn2Ry4opUHzMxFwF9becxekJleU/Ucu8Y07CJiz4j4ZUTcGxELI2JuRPwoIjZvUHbriJgaEbMi4umIeDgi/hIRJ0fERmWZy4Cflrv8tK4bbsshtCci4uiIuCUinomI+RHxg4hYr6L8ehHxiYi4NCLmlefwQEScFxEvqys7OSL6n1Ozd13bjq0rd3ZEzC7P8/GIuDoi3j2Ua9qgjYdExCXl9XqmHJ9zRkTsVldujYg4JiJuioh/lMe9MiLe2aDOpeNuys9nRsSDZf0zIuKAuvKXleceded+We21qerObNSdGBHPjojPR8TNZVufiIg7IuIXEfHSRm1tUO9mEXFieU36v3a/qt2/puzSNkbEa8pzeqI89gURsWP1V2GFupZ2TUXEbhHxm4h4LCIeKb/2E8tyW5fX9oHye+H3EbFLg/q2j4ivltf+gYj4Z0TcVf68TKgrOw34fbn4hbrvw0kNznW/8lwfq/n+XeFrEhFbRcSj5ffZFnXHXDsibouIJRGx91Cvk9RpzAhpWEXEe4EfA/8EzgPmAtsB7wMOjIiXZebdZdnNgOuBdYHpwNnAmsBWwBHAD4CHgGkU3W8HAeeyfNfbo0No1neBDwMLgKnAorKuPYHVgYV15XcEvkyR5bgAeAR4HvBmYP+IODAzf1OWnUnRdfgF4K6yrf0uq/n8Q+DWss4FwEbAG4HTI+L5mfn5IZwHEREUQeGRwIPAr4AHgAnAa4DbgRll2dWBi4C9KbInJwLPAt4O/CIids3MzzQ4zBbAH4HZwOnAhsAhwLkR8brM7P+FO608x/pznzOUc6k4t98ArwCuBX4CLAYmApOAK4E/DVLHVsBVFFmqS4Ezyv3fAbwpIg7OzF832PUAiu+JC4GTgZ0ovj67R8ROmflgE6eyO/Ap4HKKn4WdgbcBO0fEm8v2/RU4jeJavw34bURsnZlP1tTzNuAoigDnGorv0xew7Gdpt8ycX5b9v/L9yPK4l9XUM6eufW8H9qs51y2rTiQz74yI9wH/C5wREXtl5uJy80nADsCxmXn5gFdE6mSZ6ctX0y+K/1wT2LJm3fYU/1nPAsbXlX8tsAQ4p2bdh8o6PtKg/rWBtWqWJ5dlJzfZzleU+80CNqxZvybFL9sE5tTtsx6wcYO6JgD3ALc12JbAZQO0Y5sG61YHLqEIzMYP8XymlMf6I7Be3baxwGY1y58uy04HVqtZ/9yar98ratZvWa5L4At1de/bX9dQz32wr1n9fhQBQ9Z+j9RsGwNs0KCt0+rKXVSu/2yD74PFFIH1Og3auBjYp26fr5TbPjnEr82kmuv3rrpt/12uf7hB2z7f6OcAGA+s0eA4b6D4WfphxfGPrWhf/7n2AfsN5WtSs/6kcttXyuX3lMu/B8YM5fr48tWpL7vGNJw+AIyj+A99fu2GzLyUIkN0YEQ8u26/p+srysynMnOF9SvhveX7lzPz4Zr6n6EIFFaQmY9lgwxAZs4DfgnsEBHPa6YRmXlHg3ULKbI0qwH7DLGqD5Xv78/M5QYKZ+aSzFxQs+pfKH5ZfTyX/RVPZt4PHF8uvq/BMe4CvlRX90XA3cAeQ2znqmj0/dCXmY8MtFPZXfQGinZ+vW7/ayiyQxtSZFrqnZmZl9Stm1q+N3vOV2Xmz+rWnVq+PwZ8tW7baeX7rrUrM3N+Zv6zvvLMvBi4hSI4XRnn5rKM5lB9HLgR+FREHE0RGD1AEfD1rWQ7pI5g15iG08vL970jYvcG259LkbXYnqKL4zzgv4ATI2Jfir/mrwZuzcxssP/KeEn53ih1fyVFJmAFEfFK4CMU5/RciuxNrfEUv3CHpAycPkUR8DwPWKtBfYPVsTbwQuC+zLxhkLLPBrYF5mfjAbCXlu8vbrBtZmYuabB+Lsu+xiPhVoquxsPK8SjnUnQjzSiDxsH0n8uVWQymrncp8O6y3Gl122Y0KD+3fN9gCMcerK57yvdG17b/j4b6cT8BvIsik7NL2Y6xNUWGck0a+WOzO2TmMxFxCMW5nUARYL89M+8ZeE+p8xkIaThtVL5/YpBy6wBk5l0RsQdwLMWYhf6/1OdGxDcz8/vD0Kb+AdH31W/IzCUR8VD9+oh4K0Xm5xngt8AdwFMUXQqTKMbcrDHUBkTE1hS/fDagCL4upsgMLKHo4jlyiPWtX77PH6hQqf+8F1Rs71+/foNtj1bss5gRvMGi/Hq8FvhPinEsXys3PRERpwKfzuXH0NQb1nPOzMVFLLJc8DEUjW7pX1y1reY44+o2fZtimooFFH8kzGdZtmwyxfiilbGyczz9DbiJopvxVorvY2nUMxDScOr/T369zHx8KDtk5m3AIRGxGsVfva+j6P75XkQ8lZn/PUxt2oRi8O9SETGWInirDyyOp/hre7eyfbX7/IgiEGrGx8vjvDczp9XVdxhFIDQUj5bvg2aPWHbem1Zs36yu3Ejo7zJZ4f+ZiFi/0Q5l99fHgI9FxLYU1/r9wNEUAcwRAxyvE855WETEcykG+N9MMY7ribrth61C9SubbT2GIgh6kGLQ9qcpbiqQRjXHCGk4XVe+v7rZHTNzcWb+KTO/BvT/J/+WmiL93QnN/nX+5/K9UfDyahr/MbAtRfdcfRA0BnhVxXH6BmjbtuX72Q22DTmoysynKH4xbhIRjbq0ass+QZHJGh8R2zUo8pry/c8Ntg2X/jE9Exts263BuuVk5qwyEN4beJLirq6B9HcXvqoMrOu14pyHy9YU/z9f3CAImlBur7eyPyODiohXAF+kuCvxheX7cRFR9fMgjRoGQhpOP6C4A+o7EbF9/caIWD0iXl2zvEdEbNKgnv51tTPc9ndhNTVImWW3dH82IjasOfaaFHcFNTIH2C5q5j0qx2t8geK26kYeovEv/P76oOhWW6ocF9VosPJA+rsLfxR18yBFxJhySoJ+p1DM8fONMvvVX25jijuV+suMlBkUAeLhEfGsmuNvSN1g5nL9VhHxggb1bEDRdTjg4PlyMPtvKbobP1pX957A4RTB2TnNnESbzCnfX1X3tVuH4pb8RoHeyv6MDCgiNqAYaL4EODQz76OYTmExxS31Gw20v9Tp7BrTsMnMv0bEv1D8cr0lIn5DMa5gHMV/zq+muNNkh3KXw4EPRsTlFLe3PwJsAxxIMQ/Rd2uqv5YiMPpo+Yu0f8zPCfV3T9W16eqIOIGiu+3miPgly+YReoTG40m+QzG/yg0RcXZZ/pUUQdD5ZfvqXQIcGhHnUwwEXwxckZlXUNxh817gf8v65lP8Vb0fcBbFL5Wh+glFVuo9wN8j4lyKa7o5xRQFp1CMuQL4JrB/ea43RsR0inmE3kExAPzrmXlVE8duSmYuiIifUXRnzYyICyjmjHojxXxK9VmtXYBzIuJPFJmve4DnlO0fx7IxQwM5imLA/Tci4g0UwVj/PEJ9FN2TTwywf0fIzHsj4kzgUIprdzHFGKjXU4xdm0ndXWYUWZr5FN+HCykG8ydwembetQrNOYXi5/fDmTmzbN+NEfEfFH/8/JRiji1pdGr3/fu+RueLBvMI1WzbmSITcxdFQPMwxS+2HwGvrSm3J8VEgzeWZZ6mCIh+CrywQb37UQRET7JsvpYVjt9gv6AYY3Jb2Z57KG5bX688jzkN9plM8cvmKYoxEeeU53VsedxJdeWfC/ycIkBbQt18LhRjKy6lCL6eoLgb6i0MMvfLAOf0Loo74R6j+MV4J/Az4CV15dYEPlNe/6drjn1Ygzq3pMHcPDXbLyv+y1hhfeUcShSZnG8A81g2x9SnKf4Iq59HaALFXYRXUwzo/We534XA/kNtK8UYqh+W338Ly6/f/wG7V3ydhzzX0SBfk8qv5RCu7QrHoQhav1xes2co7mI7kWK8WdXXYneKoPwxisBv6ffqYOfaqB0sm+vr3Iryvyq3f6yZ719fvjrpFZnDdZeyJEnS6OIYIUmS1LMMhCRJUs8yEJIkST3LQEiSJPUsAyFJktSzRvU8QosenO0tb1IbrLV505OHSxomixfOj1Yer5W/a8dtvHVLzw3MCEmSpB5mICRJknrWqO4akyRJI6xvyeBlRjEzQpIkqWeZEZIkSdWyr90tGFFmhCRJUs8yIyRJkqr1mRGSJEnqSmaEJElSpXSMkCRJUncyIyRJkqo5RkiSJKk7mRGSJEnVHCMkSZLUnQyEJElSz7JrTJIkVfOhq5IkSd3JjJAkSarmYGlJkqTuZEZIkiRVc0JFSZKk7mRGSJIkVfKhq5IkSV3KjJAkSarmGCFJkqTuZEZIkiRVc4yQJElSdzIjJEmSqvmsMUmSpO5kRkiSJFVzjJAkSVJ3MhCSJEk9y64xSZJUzQkVJUmSupMZIUmSVM3B0pIkSd3JjJAkSarmGCFJkqTuZEZIkiRVyvQRG5IkSV3JjJAkSarmXWOSJEntFxH7RcTtETErIo5psH2DiDgnIm6KiD9GxAsHq9OMkCRJqtYhd41FxFjgROD1wDzg+og4LzNvrSn2GWBmZr41InYoy+8zUL1mhCRJ0miwBzArM2dn5kLgTOCgujI7AZcAZOZfgS0jYpOBKjUQkiRJ1bKvda+BjQfm1izPK9fVuhF4G0BE7AFsAUwYqFIDIUmS1BEiYkpEzKh5Tand3GCXrFv+KrBBRMwEPgTcACwe6JiOEZIkSdX6WjePUGZOBaZWbJ4HTKxZngDcU7f/48B7ASIigDvLVyUzQpIkaTS4HtguIraKiNWBQ4HzagtExPrlNoD3AVeUwVElM0KSJKnjZebiiDgauAgYC5ySmbdExFHl9pOBHYHTImIJcCvwr4PVayAkSZKqddCEipk5HZhet+7kms/XAts1U6ddY5IkqWeZEZIkSdU6ZELFkWJGSJIk9SwzQpIkqVoHjREaCWaEJElSzzIjJEmSqjlGSJIkqTuZEZIkSdXMCEmSJHUnM0KSJKlSZuseutoOZoQkSVLPMiMkSZKqOUZIkiSpO5kRkiRJ1ZxZWpIkqTsZCEmSpJ5l15gkSarmYGlJkqTuZEZIkiRVc7C0JElSdzIjJEmSqjlGSJIkqTuZEZIkSdUcIyRJktSdzAhJkqRqjhGSJEnqTmaEJElSNTNCkiRJ3cmMkCRJquZdY5IkSd3JjJAkSarmGCFJkqTuZCAkSZJ6ll1jkiSpmoOlJUmSupMZIUmSVM3B0pIkSd3JjJAkSarmGCFJkqTuZEZIkiRVc4yQJElSdzIjJEmSqpkRkiRJ6k5mhCRJUrXMdrdgRJkRkiRJPcuMkCRJquYYIUmSpO5kRkiSJFUzIyRJktSdzAhJkqRqPmtMkiSpOxkISZKknmXXmCRJquZgaUmSpO5kICRJkqpltu41iIjYLyJuj4hZEXFMg+3rRcT5EXFjRNwSEe8drE4DIUmS1PEiYixwIrA/sBNwWETsVFfsg8CtmbkLMAn4VkSsPlC9jhGSJEnVOmeM0B7ArMycDRARZwIHAbfWlEng2RERwDrAw8DigSo1IyRJkkaD8cDcmuV55bpaPwB2BO4B/gJ8JHPgiZDMCEmSpGotzAhFxBRgSs2qqZk5tX9zg13qBxbtC8wEXgtsA/w2Iq7MzMerjmkgJEmSOkIZ9Eyt2DwPmFizPIEi81PrvcBXMzOBWRFxJ7AD8MeqY9o1JkmSqmVf614Dux7YLiK2KgdAHwqcV1fmbmAfgIjYBHg+MHugSs0ISZKkjpeZiyPiaOAiYCxwSmbeEhFHldtPBo4HpkXEXyi60j6VmQ8OVK+BkCRJqpR9g8/v0yqZOR2YXrfu5JrP9wBvaKZOu8YkSVLPMiMkSZKqdc48QiPCjJAkSepZZoQkSVK1we/mGtXMCEmSpJ5lICRJknqWXWOSJKlaB90+PxLMCEmSpJ5lRkiSJFXz9nlJkqTuZEZIkiRVMyMkSZLUncwISZKkauldY5IkSV3JjJAkSarmGCFJkqTuZEZIkiRVc2ZpqTlXXTeDAw59H/u/81/4yelnrbD9scef4MOf/iJvfc8HOPR9H+Hvs+cAcOdd8zj4yA8ufe35+rdx+i/OaXHrpdFr3zdM4pabr+Cvt17FJz/xwRW2P//523DVFefx1BOz+fjH3r/cth9P/Rb3zLuRmTdc0qrmSh3BQEjDasmSJXzpWyfyw28dz3k/+xHTf3cZd9x513JlfnzaL9hhu20457Qf8l+f/3989bsnA7DVFhM4+9QTOfvUEznrlO+z5pprss/er2jHaUijzpgxY/j+977MAQe+m513eQ2HHPIWdtxxu+XKPPzwo3z0Y5/n29/50Qr7n3baWbzpgHe1qrkaTbKvda82aHsgFBFrDGWdRoe/3PY3njdhcyaO34xx48ax/z57c+mV1y1X5o45d/Oyl+4CwNZbTGT+gvt48OFHlitz3YyZTBy/GZtvuknL2i6NZnvs/mLuuGMOd955N4sWLeKss87lzQfuu1yZBx54iBl/upFFixatsP+VV/2Bhx95tEWtlTpH2wMh4NohrtMocP8DD7Lpc5+zdHmT527M/Q88tFyZ52+7Nb+7/BoA/nLr7Sy4737uu//B5cpceMnlvPF1e498g6Uusfn4TZk7756ly/PmL2DzzTdtY4vUNfqyda82aFsgFBGbRsRLgbUi4sUR8ZLyNQl4VrvapVXTaN6tiOWX33fEO3j8iSc5+MgP8rNfnscO223D2LFjl25ftGgRl131B97w2lePcGul7hH1P2hAdvlEeNJwaOddY/sCk4EJwLdr1j8BfKZqp4iYAkwBOOlbX+J97zlsBJuoZm3y3I259/4Hli7fd/+DPGfjjZYrs87aa/Olz34cKP6j3vftk5mw+bIusCuvm8GO22/Dxhtu0JpGS11g/rwFTJyw+dLlCeM3Y8GC+9rYIml0aFsglJmnAqdGxMGZeXYT+00FpgIsenC2f+50mBfusD13z7uHeffcyybP2YgLL7mcr3/hU8uVefyJJ1lrzTUYN24cZ5//G166686ss/baS7dP/+1lvPH1k1rccml0u37GTLbddiu23HIi8+ffyzvfeRBHvGfFO8ekZmWXT6jYCfMIXRIR3wb2KpcvB76YmY+1sU1aSautNpbPfOwDvP/jn2PJkiW89YA3sO3WW/CLcy4A4JC3vonZd83lM8d/k7FjxrD1ls/ji5/+6NL9n37mGa69/ga+8MkPt+kMpNFpyZIlfOSjn2P6BT9n7JgxTDv1F9x669+Y8m9HADD1x6ezySbP4Q/XXsi6665DX18fH/7Qv7HzLpN44okn+Z/TT2TvvV7OxhtvyJzZMzjui9/kp9PObPNZSSMv2t2HHBFnAzcDp5arjgB2ycy3DbavGSGpPdba3PFbUrssXjh/xQFhI+ipL7+nZb9r1/7saS09N+iMjNA2mXlwzfJxETGzXY2RJEm9oxMCoacj4lWZeRVARLwSeLrNbZIkSdC2iQ5bpRMCoQ9QDJper1x+BDiyje2RJEk9ohMCoduArwPbAOsDjwFvAW5qX5MkSRLQ9Q9d7YRA6FzgUeDPwPz2NkWSJPWSTgiEJmTmfu1uhCRJaqDL5xHqhGeNXRMRO7e7EZIkqfd0QkboVcDkiLgT+CcQQGbmi9rbLEmS5Bihkbd/uxsgSZJ6U9sDocy8q91tkCRJFbp8HqFOGCMkSZLUFm3PCEmSpA7W5WOEzAhJkqSeZSAkSZJ6ll1jkiSpUjqhoiRJUncyIyRJkqo5WFqSJKk7mRGSJEnVzAhJkiR1JzNCkiSpmo/YkCRJ6k5mhCRJUjXHCEmSJHUnM0KSJKlSmhGSJEnqTmaEJElSNTNCkiRJ3clASJIkVevra91rEBGxX0TcHhGzIuKYBts/EREzy9fNEbEkIjYcqE4DIUmS1PEiYixwIrA/sBNwWETsVFsmM7+Rmbtm5q7Ap4HLM/Phgeo1EJIkSaPBHsCszJydmQuBM4GDBih/GHDGYJU6WFqSJFXrnMHS44G5NcvzgD0bFYyIZwH7AUcPVqkZIUmS1BEiYkpEzKh5Tand3GCXqijtQODqwbrFwIyQJEkaSAszQpk5FZhasXkeMLFmeQJwT0XZQxlCtxiYEZIkSaPD9cB2EbFVRKxOEeycV18oItYD9gbOHUqlZoQkSVKlzM4YI5SZiyPiaOAiYCxwSmbeEhFHldtPLou+Fbg4M58aSr0GQpIkaVTIzOnA9Lp1J9ctTwOmDbVOAyFJklStc+4aGxGOEZIkST3LjJAkSapmRkiSJKk7mRGSJEmV0oyQJElSdzIjJEmSqpkRkiRJ6k5mhCRJUrW+djdgZJkRkiRJPctASJIk9Sy7xiRJUiVvn5ckSepSZoQkSVI1M0KSJEndyYyQJEmq5u3zkiRJ3cmMkCRJquRdY5IkSV3KjJAkSarmGCFJkqTuZEZIkiRVcoyQJElSlzIjJEmSqjlGSJIkqTuZEZIkSZXSjJAkSVJ3MhCSJEk9y64xSZJUza4xSZKk7mRGSJIkVXKwtCRJUpcyIyRJkqqZEZIkSepOZoQkSVIlxwhJkiR1KTNCkiSpkhkhSZKkLmVGSJIkVTIjJEmS1KXMCEmSpGoZ7W7BiDIjJEmSepYZIUmSVMkxQpIkSV3KQEiSJPUsu8YkSVKl7HOwtCRJUlcyIyRJkio5WFqSJKlLmRGSJEmV0gkVJUmSupMZIUmSVMkxQpIkSV3KjJAkSarkPEKSJEkdICL2i4jbI2JWRBxTUWZSRMyMiFsi4vLB6jQjJEmSKmW2uwWFiBgLnAi8HpgHXB8R52XmrTVl1gdOAvbLzLsj4rmD1WtGSJIkjQZ7ALMyc3ZmLgTOBA6qK3M48KvMvBsgM+8frFIzQpIkqVIHjREaD8ytWZ4H7FlXZntgXERcBjwb+F5mnjZQpQZCkiSpI0TEFGBKzaqpmTm1f3ODXeo77lYDXgrsA6wFXBsR12Xm36qOaSAkSZIqtTIjVAY9Uys2zwMm1ixPAO5pUObBzHwKeCoirgB2ASoDIccISZKk0eB6YLuI2CoiVgcOBc6rK3Mu8OqIWC0inkXRdXbbQJWaEZIkSR0vMxdHxNHARcBY4JTMvCUijiq3n5yZt0XEb4CbgD7gJ5l580D1RnbKfXErYdGDs0dv46VRbK3NX93uJkg9a/HC+S0dvXznLq9v2e/arW78bctHZts1JkmSepZdY5IkqVIH3T4/IswISZKknmVGSJIkVco0IyRJktSVzAhJkqRK2dfuFoysAQOhiOhjxemrhyIz0yBLkiR1tMGClStYuUBIkiR1gb4uHyM0YCCUmZNa1A5JkqSWs/tKkiRV8q4xSZKkLrVSGaGI2AzYBxgPrNGgSGbm8avSMEmS1H7dPrN004FQRBwHHFO3b7BsUHX/ZwMhSZLU0ZrqGouIdwGfB64E3k4R9JwKHA78mOKR92cCrx3eZkqSpHbIbN2rHZrNCH0AmAfsl5mLIwJgTmaeCZwZEecAFwBnDG8zJUmShl+zg6V3BqZn5uKadWP7P2TmRcBFwCeGoW2SJEkjqtmM0DjgoZrlp4H16srcDBy1Ko2SJEmdodsHSzebEVoAbFazfDfworoy44HFSJIkdbhmA6EbKLrH+l0KvDoijoiItSPiTcDBZTlJkjTK9WW07NUOzQZCvwZeEBFblctfBR4DpgGPA+dR3En2ueFqoCRJ0khpaoxQZk6jCHr6l+dGxO7AfwDbAHOAkzLzL8PXREmS1C7d/oiNVX7WWGbeCRw9DG2RJElqKR+6KkmSKrVrosNWaSoQiojnDbVsZt7dfHMkSZJap9mM0ByWPVNsILkSdUuSpA7Trru5WqXZYOU0GgdC6wO7AlsAlwF3rUqjJEmSWqHZu8YmV22LiDEUD2Q9Cjhy1ZolSZI6QbffNdbsPEKVMrMvM4+j6D776nDVK0mSNFJGYhzPNcB7RqBeSZLUYt1+19iwZYRqbAisPQL1SpIkDathzQhFxOuAQyieQC9JkkY57xqrERGXDlDPRKB/nqEvrkqjJEmSWqHZjNCkivUJPAJcBHwzM6sCpmG11uavbsVhJNW5aeKu7W6CpBbp9rvGmr19fiTGFEmSJLWFgY0kSepZTQVCEXFpRAx4a3xEvHuAsUSSJGkU6cto2asdms0ITQK2HKTMFsDeK9MYSZKkVhqJCRXXAhaPQL2SJKnFunw+xZUKhBpek4gIitvn3wjMXZVGSZIktcKggVBE9LF88HNsRBw70C7Af61iuyRJUgdwQkW4gmWB0F7A3RQPVq23BHgIuAT4yXA0TpIkaSQNGghl5qT+z2V26KeZ6czRkiT1ACdUXN5WwKMj0A5JkqSWazYQuh94TkQ8nZkL6zdGxBrAJsD9mfnMcDRQkiS1T1+7GzDCmp1H6D+B24F1KravDfwV+MyqNEqSJKkVmg2E9gd+l5kPN9pYrv8dcMCqNkySJLVfEi17tUOzgdCWwN8GKfM3Bp99WpIkqe2aHSM0jsG7CxNYc+WaI0mSOklfl08t3WxGaDaDP0dsEnDXSrVGkiSphZoNhM4DXhoRn2y0MSKOAV4C/N8qtkuSJHWAPqJlr3Zotmvsm8C7gK9ExDuBi4H5wHhgX2BXipmnvz6MbZQkSRoRTQVCmflIREwCfga8nCL7k7A0jLsGeHdmPjKMbZQkSRoRTT99PjPnAK+MiJcALwPWp5ht+rrM/PNwNk6SJLVXu25rb5WmA6F+ZdBj4CNJkkatlQqEImIzYB+KsUFrNCiSmXn8qjRMkiS1Xyc9YiMi9gO+B4wFfpKZX63bPgk4F7izXPWrwR4U33QgFBHHAcfU7RsUY4VqPxsISZKkYRERY4ETgdcD84DrI+K8zLy1ruiVmTnkJ1w0dft8RLwL+DxwJfB2iqDnVOBw4McUgeOZwGubqVeSJHWmDnrExh7ArMycXT74/UzgoFU9v2bnEfoARRS2X2aeU66bk5lnZuZRFM8Yeyew7qo2TJIkqcZ4YG7N8rxyXb2XR8SNEXFhRLxgsEqbDYR2BqZn5uKadWP7P2TmRcBFwCearFeSJHWgvha+ImJKRMyoeU2paUqjlFH9A0D+DGyRmbsAJzCECZ6bDYTGAQ/VLD8NrFdX5mZglybrlSRJPS4zp2bmbjWvqTWb5wETa5YnAPfU7f94Zj5Zfp4OjIuIjQc6ZrOB0AJgs5rlu4EX1ZUZDyxGkiSNeq3MCA3iemC7iNgqIlYHDqV49NdSEbFpRET5eQ+KOOehFWqq0exdYzdQdI/1uxSYEhFHAL+ieODqwcDVTdYrSZJUKTMXR8TRFENwxgKnZOYtEXFUuf1kihu5PhARiyl6rQ7NzPrus+XEINuXLxwxGTgJeEFm3hkREymCow1qii0CJmXmdUOueCWttvr4oTde0rC5aeKu7W6C1LN2uuOClk71fMEmh7Xsd+2b7juj5dNYN/ussWnAtJrluRGxO/AfwDbAHOCkzPzL8DVRkiRpZKz0Izb6ZeadwNHD0BZJktRh+rr7UWNND5aWJEnqGqucEZIkSd2rr8ufPm9GSJIk9SwDIUmS1LPsGpMkSZW6fZ4aM0KSJKlnmRGSJEmVhvDoi1HNjJAkSepZZoQkSVKlvvD2eUmSpK5kRkiSJFXyrjFJkqQuZUZIkiRV8q4xSZKkLmVGSJIkVerr7pvGzAhJkqTeZUZIkiRV6qO7U0JmhCRJUs8yIyRJkio5j5AkSVKXMhCSJEk9y64xSZJUydvnJUmSupQZIUmSVMlHbEiSJHUpM0KSJKmSt89LkiR1KTNCkiSpkneNSZIkdSkzQpIkqZJ3jUmSJHUpM0KSJKmSGSFJkqQuZUZIkiRVSu8akyRJ6k5mhCRJUiXHCEmSJHUpAyFJktSz7BqTJEmV7BqTJEnqUmaEJElSpWx3A0aYGSFJktSzzAhJkqRKfU6oKEmS1J3MCEmSpEreNSZJktSlzAhJkqRKZoQkSZK6lBkhSZJUyXmEJEmSupQZIUmSVMl5hCRJkrqUgZAkSarU18LXYCJiv4i4PSJmRcQxA5TbPSKWRMTbB6vTQEiSJHW8iBgLnAjsD+wEHBYRO1WU+xpw0VDqNRCSJEmjwR7ArMycnZkLgTOBgxqU+xBwNnD/UCo1EJIkSZWyha9BjAfm1izPK9ctFRHjgbcCJw/1/AyEJElSR4iIKRExo+Y1pXZzg13q46fvAp/KzCVDPaa3z0uSpEp9LZxSMTOnAlMrNs8DJtYsTwDuqSuzG3BmRABsDLwxIhZn5v9VHdNASJIkjQbXA9tFxFbAfOBQ4PDaApm5Vf/niJgG/HqgIAgMhCRJ0gA65aGrmbk4Io6muBtsLHBKZt4SEUeV24c8LqiWgZAkSRoVMnM6ML1uXcMAKDMnD6VOAyFJklTJh65KkiR1KTNCkiSpUqeMERopZoQkSVLPMiMkSZIq9TWaxrCLmBGSJEk9y4yQJEmq1MqZpdvBjJAkSepZZoQkSVKl7s4HmRGSJEk9zEBIkiT1LLvGJElSJSdUlCRJ6lJmhCRJUiVvn5ckSepSZoQkSVKl7s4HmRGSJEk9zIyQJEmq5F1jkiRJXcqMkCRJquRdY5IkSV3KjJAkSarU3fkgM0KSJKmHmRGSJEmVvGtMkiSpS5kRkiRJlbLLRwmZEZIkST3LQEiSJPUsu8YkSVIlB0tLkiR1KTNCkiSpko/YkCRJ6lJmhCRJUqXuzgeZEZIkST3MjJAkSarkGCFJkqQuZUZIkiRVch4hqUn7vmESt9x8BX+99So++YkPrrD9+c/fhquuOI+nnpjNxz/2/uW2/Xjqt7hn3o3MvOGSVjVX6hpr7/VStvntj9j20h+z0fvfscL2Z+25M8+feRZbn38CW59/AhsffRgAq281fum6rc8/gefP/F82nHxQq5svtYUZIQ2rMWPG8P3vfZn93ngY8+Yt4Lprp3P+ry/mttv+vrTMww8/ykc/9nkOOmi/FfY/7bSzOOmkn/LTn36vlc2WRr8xY9js2A9w15GfY9G9D7L1Od/hiUuuY+GsucsV+8f1tzD3345bbt3CO+cz+8APLa1n+2tO44mLr2lVy9XhfOiq1IQ9dn8xd9wxhzvvvJtFixZx1lnn8uYD912uzAMPPMSMP93IokWLVtj/yqv+wMOPPNqi1krdY61dtmfhXfewaO69sGgxj/36Cp79upc1Xc/ar9iFhXcvYNE9D4xAK6XO0/ZAKCK2jojzI+LBiLg/Is6NiK3b3S6tnM3Hb8rcefcsXZ43fwGbb75pG1sk9YbVNtmIRQseXLq8+N4HGbfJRiuUW+vFO7D1r0/geaccxxrbPW+F7esesBePnX/5iLZVo0tfC1/t0PZACPg5cBawKbA58L/AGW1tkVZaRKywLrO706pSR2jws1fvmVtm8fe93svsAz7Ew6edz4STP7d8gXGr8ex99uTx6VeNUCOlztMJgVBk5umZubh8/Q8DTGQZEVMiYkZEzOjre6qFzdRQzJ+3gIkTNl+6PGH8ZixYcF8bWyT1hsX3Psi4zTZeurzaphuz6L6HlivT9+TT5D+eAeDJy2YQq63G2A3WXbp9nb1345lb7mDJQ4+2pM0aHbKF/9qhEwKh30fEMRGxZURsERGfBC6IiA0jYsP6wpk5NTN3y8zdxoxZuw3N1UCunzGTbbfdii23nMi4ceN45zsP4vxfX9zuZkld7+mb/sbqW45n3IRNYNxqrHfAXjx5yR+WKzN24w2Wfl7zRdsTY4Iljzy+dN16B9otpt7TCXeNHVK+v79u/b9QZIYcLzSKLFmyhI989HNMv+DnjB0zhmmn/oJbb/0bU/7tCACm/vh0NtnkOfzh2gtZd9116Ovr48Mf+jd23mUSTzzxJP9z+onsvdfL2XjjDZkzewbHffGb/HTamW0+K2kUWNLHvcf9kOdNO54YM4ZHf/lb/vn3u9ngsP0BeOSMC1l3/1eyweFvhCVL6HtmIfM+8vWlu8eaa7D2K1/Mgs/+oF1nILVFjObxG6utPn70Nl4axW6auGu7myD1rJ3uuGDwAWHD6MgtD27Z79pT55zd0nODDsgIRcR7Gq3PzNNa3RZJktRb2h4IAbvXfF4T2Af4M2AgJElSm/WN4p6joWh7IJSZH6pdjoj1gNPb1BxJktRD2h4INfAPYLt2N0KSJA0wn02XaHsgFBHns+w6jwV2pJhgUZIkaUS1PRACvlnzeTFwV2bOa1djJEnSMn1dnhNqeyCUmUtn74qIAzLz6na2R5Ik9Y5OmFm61hfb3QBJkrRMJz1iIyL2i4jbI2JWRBzTYPtBEXFTRMwsH8f1qsHqbHtGqE7LJ1KSJEmdLyLGAicCrwfmAddHxHmZeWtNsUuA8zIzI+JFFGOOdxio3rZnhCJijZrF9zdYJ0mS2qSvha9B7AHMyszZmbkQOBM4qLZAZj6Zyx6ZsTZDuOmt7YEQcG3/h8z8Y/06SZIkYDwwt2Z5XrluORHx1oj4K3ABxXNLB9S2rrGI2JTiBNaKiBezrFtsXeBZ7WqXJElappV3jUXEFGBKzaqpmTm1f3ODXVZoXGaeA5wTEXsBxwOvG+iY7RwjtC8wGZgAfLtm/RPAZ9rRIEmS1D5l0DO1YvM8YGLN8gTgngHquiIitomIjTPzwapybQuEMvNU4NSIODgzz25XOyRJUrWh3M3VItcD20XEVsB84FDg8NoCEbEtcEc5WPolwOrAQwNV2s6usY83+twvM79dv06SJPWmzFwcEUcDF1E8ieKUzLwlIo4qt58MHAy8JyIWAU8Dh9QMnm6onV1jz27jsSVJ0iiTmdOB6XXrTq75/DXga83U2c6usePadWxJkjQ0Q7itfVRrZ9fYJzPz6xFxAo1HfX+4Dc2SJEk9pJ1dY7eV7zPa2AZJkjSAQYbYjHrt7Bo7v3w/tV1tkCRJva3tzxqLiN/TuGvstW1ojiRJqtHKCRXboe2BEPD/aj6vSXHr2+I2tUWSJPWQtgdCmfmnulVXR8TlbWmMJElajneNjbCI2LBmcQywG7Bpm5ojSZJ6SNsDIeBPLBsjtBiYA/xr21ojSZKW6qBHbIyITgiEdgL+HXgVRUB0Jd5SL0mSWqATAqFTgceB75fLhwGnA+9oW4skSRLgXWOt8PzM3KVm+fcRcWPbWiNJknpGJwRCN0TEyzLzOoCI2BO4us1tkiRJOLP0iImIv1CMCRoHvCci7i6XtwBubVe7JElS72hnRuiANh5bkiQNgfMIjZDMvKtdx5YkSYLOGCMkSZI6VLfPIzSm3Q2QJElqFwMhSZLUs+wakyRJlbp9QkUzQpIkqWeZEZIkSZW6fUJFM0KSJKlnmRGSJEmVHCMkSZLUpcwISZKkSk6oKEmS1KXMCEmSpEp93jUmSZLUncwISZKkSt2dDzIjJEmSepgZIUmSVMl5hCRJkrqUGSFJklTJjJAkSVKXMhCSJEk9y64xSZJUKZ1QUZIkqTuZEZIkSZUcLC1JktSlzAhJkqRKaUZIkiSpO5kRkiRJlbxrTJIkqUuZEZIkSZW8a0ySJKlLmRGSJEmVHCMkSZLUpcwISZKkSo4RkiRJ6lJmhCRJUiVnlpYkSepSBkKSJKln2TUmSZIq9Xn7vCRJUvtFxH4RcXtEzIqIYxpsf1dE3FS+romIXQar04yQJEmq1CmDpSNiLHAi8HpgHnB9RJyXmbfWFLsT2DszH4mI/YGpwJ4D1WtGSJIkjQZ7ALMyc3ZmLgTOBA6qLZCZ12TmI+XidcCEwSo1IyRJkip10Bih8cDcmuV5DJzt+VfgwsEqNRCSJEkdISKmAFNqVk3NzKn9mxvs0jBKi4jXUARCrxrsmAZCkiSpUivHCJVBz9SKzfOAiTXLE4B76gtFxIuAnwD7Z+ZDgx3TMUKSJGk0uB7YLiK2iojVgUOB82oLRMTzgF8BR2Tm34ZSqRkhSZJUqVPGCGXm4og4GrgIGAuckpm3RMRR5faTgf8ENgJOigiAxZm520D1GghJkqRRITOnA9Pr1p1c8/l9wPuaqdNASJIkVeqUeYRGimOEJElSzzIjJEmSKnXKGKGRYkZIkiT1LDNCkiSpkmOEJEmSupSBkCRJ6ll2jUmSpEqZfe1uwogyIyRJknqWGSFJklSpz8HSkiRJ3cmMkCRJqpROqChJktSdzAhJkqRKjhGSJEnqUmaEJElSJccISZIkdSkzQpIkqVKfGSFJkqTuZEZIkiRVSu8akyRJ6k5mhCRJUiXvGpMkSepSBkKSJKln2TUmSZIq+YgNSZKkLmVGSJIkVXKwtCRJUpcyIyRJkir5iA1JkqQuZUZIkiRVcoyQJElSlzIjJEmSKjmPkCRJUpcyIyRJkio5RkiSJKlLmRGSJEmVnEdIkiSpS5kRkiRJldK7xiRJkrqTgZAkSepZdo1JkqRKDpaWJEnqUmaEJElSJSdUlCRJ6lJmhCRJUiVvn5ckSepSZoQkSVIlxwhJkiR1KTNCkiSpkhkhSZKkLmVGSJIkVerufJAZIUmS1MOi2/v+1LkiYkpmTm13O6Re48+etIwZIbXTlHY3QOpR/uxJJQMhSZLUswyEJElSzzIQUjs5RkFqD3/2pJKDpSVJUs8yIyRJknqWgZA6SkRMjojN290OqRNFxJYRcXMT5d8cEceUn4+NiP9XfvbnTCoZCKnTTAb8D1oaBpl5XmZ+tcGmyfhzJgEGQhph5V+wt0XEjyPiloi4OCLWiohdI+K6iLgpIs6JiA0i4u3AbsDPImJmRKzV7vZLHWi1iDi1/Nn5ZUQ8KyLmRMTGABGxW0RcVn6eHBE/qN15sJ+ziJgUEb+uWf5BREwuP8+JiK9FxB/L17Yje6rSyDMQUitsB5yYmS8AHgUOBk4DPpWZLwL+AnwhM38JzADelZm7ZubT7Wqw1MGeD0wtf3YeB/69mZ2H4efs8czcA/gB8N0m95U6joGQWuHOzJxZfv4TsA2wfmZeXq47FdirHQ2TRqG5mXl1+fl/gFe1+Phn1Ly/vMXHloadgZBa4Z81n5cA67epHVI3qJ/zJIHFLPv/fM1mKouIPcsuspkR8ea6uhrVlxWfpVHJQEjt8BjwSES8ulw+AujPDj0BPLstrZJGh+dFRH8m5jDgKmAO8NJy3cFDqGPpz1lm/qHsIts1M88D7gJ2iog1ImI9YJ+6fQ+peb925U9D6gyrtbsB6llHAidHxLOA2cB7y/XTyvVPAy93nJC0gtuAIyPiR8DfgR8CfwT+OyI+A/xhCHVMo+LnLDPnRsRZwE1l/TfU7btGRPyB4g/pw1b1ZKR2c2ZpSdKQRMQcYLfMfLDdbZGGi11jkiSpZ5kRkiRJPcuMkCRJ6lkGQpIkqWcZCEmSpJ5lICS1WERk/7OgVqGOLct6pg1Pq9qv0XUpn5ieETFphI7ZdddRUnMMhCR1teEIPCV1LydUlNTJfgCcCdw9QvXPB3akmO1cUg8yEJLUscqJ+0Zs8r7MXAT8daTql9T57BpT16kd9xER20TELyPioYh4IiIujogXluWeExFTI2JBRDwTEddHxGsq6lwvIr4SEbeXZR+JiIsi4nUV5VePiM9HxB0R8c+IuDMivhQRawzQ7tUi4t8j4rqIeDwi/hERN0TE0RGxSj+rETGpvCbHRsTLI+J3EfFYeU0uiojdGuyzdHxORBweEX+IiCfL2YX7yzwrIj5dPrDzqXL7tRHR8NELzV6XgcYIRcQOEXFKRMwp67o/Iq6MiA+U2ydHRP9EaXuX9fS/ji3LVI4RiojNIuLEsv6FEfFARPwqIl7aoOzksp7JEfGaiLisvLaPR8QFEbFjg302iYhvlt9TT0XEo+XnaRGxdaPrIWn4mRFSN9uS4rlLt1E8W2lL4K3AZVE8tPI3wOPAL4ANgUOBCyNi+8xc2hUTEesDVwM7AdcD3wU2Bt4JXBwRH8jMH9WUD+As4CDgDorundWBfwF2btTQiBgHnA/sC9wO/Bx4BngNcAKwJ8XDaVfVnsCngd8BJwLbAm8D9oqIN2TmlQ32+Q/g9WX7fg+sV7Z5feBS4MXAn4FTKP642hf4eUS8IDM/V3OOTV+XKhHxJuB/gTUovo5nAOsDuwCfpHj+1kzgOOALFA8SnVZTxWWD1L8VxcNMNy/P8QxgIvAO4E0RcXBm/rrBrgeU53chcDLF98wbgd0jYqf+R1NE8Yy9q4FtgN9SXNsAtij3/yXFM/gkjbTM9OWrq14UAU+Wr8/Wbft8uf5hil9UY2q2HVFu+07dPj8q1/+Icjb2cv12FGNL/glsWbP+8LL8tcCaNes3pAgAEris7hjHlutPAMbWrB8L/He57aAG5zhtiNdkUs01Obpu20Hl+r/XXY/+Nj0FvLhBndPK7Z+sW78mRXDSB+w6TNdlUs26jcvrvhDYu0G7JtQtr1DvYNcRuKji++cVwGLgIWCdmvWTy/KLgX3q9vlK/XUCDmz0vVZuWx14drt/jnz56pWXXWPqZnOAr9atO7V8XwP4RGb21Wz7OcUvsl37V5SZmncDTwKfzsylz6TJzL8D36f4xfWemnreW75/JjOfqSn/MHB8fSPLbq+jgXuBj2Xmkpp9llBkZBJ412AnPASzgJNqV2TmucDlFNmhVzfYZ2pmLvcE8ojYiOK6zMjMr9fV9wzwKYoMx+E1m5q6LgM4ElgX+GFmXl6/MTPnNVHXCiJiAvAGigHa9ed2DUV2aEOKTFq9MzPzkrp1U8v3PRqUf7p+RWYuzMwnmm23pJVj15i62czaoKJ0T/n+t/pfNpm5JCLuAybUrN4BeBZwdfkLu96lwOcouof6vYQiG3JVg/KXNVi3PbARRUbmc0UP0gqepri7aVVdWRf81bZrb4rzqA8u/tig/O4U2aql423qjCvfa9vc7HWp8rLy/cIm9mlG/9fyyiwGU9e7lCIIfDFwWt22GQ3Kzy3fN6hZdznFHWvHRMRLgOkUXWWNvmcljSADIXWzFW6JzszFZaBRdbv0Ypb9EodyPAywoKJ8//r16/Z5uOKX6L0N1m1Uvm9HMZ6lyjoDbBuq+yrW97drvQG21epv8+7lq0ptm5u9LlXWL9/nN7FPM1bma97v0foVNd9zY2vWPR4RL6MYw/RminFVAA9GxEnAlyquk6RhZteYNLD+gGnTiu2b1ZXr/7xh2a1Wr1E9/fuek5kxwGur5pu/gk0q1ve3q1GAmA3W9Zf7ziBtfk3dPs1clyqPlu/jm9inGSvzNW9aZs7LzH8Fngu8EPgwxdij/yxfklrAQEga2O3AP4BdI2KDBtv7f9H/uWbdnyl+tl7VoPykBuv+SvHL/WUVQcJwelXFrfiTyvcbGmxr5I8U3VyNxhRVafa6VLmufN9/iOX7qMnGDEH/NXhVRDTKmjf6mq+0LNySmSdQ3J0H8JbhqFvS4AyEpAFk5kLgZxRdPF+s3RYR21D8Fb8IOL1m00/L9y9HxJo15TekGE9Uf4zFFHeLbQZ8PyLWqi9Tzmmz06qdDVB0v/17Xd0HUYwPmgU0un1+BZl5P8V12a2cF2iFgCGKOZxqs1hNXZcBnEox7cEHImKvBsedULfqIYpb34ekHGz9W4o7yj5aV/eeFAPAHwHOaaLN9W18YURs2WBTf8buHytbt6TmOEZIGtwxFJmPoyNid4q5dPrnEXo2xe3od9aUPwM4hGLsx80RcS7FuKO3U8xDtE2DYxxPMQfOUcCBEXEpxRiY51IEL68EPgvcuorn8hvgWxGxP3Ajy+YRegb414qB1FWOLtv2ReCIiLiKYgzS5hSDpHcHDgP6r83KXJcVZOaDEXE4xVw7v4+IC4GbKO4kexFF0FMbgF0CHBoR5wN/ohgHdkVmXjHAYY6iGLz8jYh4A8Ug6P55hPqA967inV2vA74dEddQZATvpxikf1BZ/zdWoW5JTTAQkgaRmQ+XEzB+miJo+DjFXVx/BL6RmRfXlc+IeAdFADWZImBYQJER+SJF0FF/jEUR8RaKu5EmU0zMtw7wAEUg8XmKDMyq+kPZhuPLdgXFXVCfzczrm6moHPC7NzCFIktyMMUcQvdR3AH3MYrMSn/5pq/LAMe+IIrZsD8F7ENxu/sjFEHFV+qKf4RinNM+FJMbjqEYpFwZCGXm7LL+z5X7TKLIQv0G+HKz16qBiygm5tyLIvhZl+Ja/Bb4dnmbvqQWiJppUSR1qfIRFb8HjsvMY9vaGEnqII4RkiRJPctASJIk9SwDIUmS1LMcIyRJknqWGSFJktSzDIQkSVLPMhCSJEk9y0BIkiT1LAMhSZLUswyEJElSz/r/rAmzhOkyHaoAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "testResults = test.classify(clf).errorMatrix(label, 'classification').getInfo()\n", "testCM = pd.DataFrame(np.asarray(testResults), index=['not','built-up'], columns=['not','built-up'])\n", "\n", "fig, ax = plt.subplots(1, figsize=(10,10))\n", "sns.heatmap(testCM/testCM.sum(axis=1), annot=True)\n", "ax.set_xlabel('model predictions', fontsize=20)\n", "ax.set_ylabel('actual', fontsize=20)\n", "plt.title(\"Test data confusion matrix\", fontsize=20);\n", "acc = (testCM.loc['built-up','built-up'] + testCM.loc['not','not']) / testCM.sum().sum()\n", "mcc = get_mcc(testCM)\n", "print(f\"Our classifier has an accuracy of {acc:.5f} on the test data.\")\n", "print(f\"Our classifier has an MCC score of {mcc:.5f} on the test data.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see the difference in the confusion matrix and in the lower accuracy -- and far lower MCC scores.\n", "\n", "On unseen data, our classifier seems to be over-predicting built up areas and producing a lot of false positives.\n", "\n", "This might impact our final analysis...this over-prediction translates into a final assessment of increased economic growth that might be in error.\n", "\n", "Let's see if increasing he number of estimators (trees) in our ensemble improves things on our test data results." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "new_params = {\"numberOfTrees\":500, \n", " \"variablesPerSplit\":None, \n", " \"minLeafPopulation\":1, \n", " \"bagFraction\":0.5, \n", " \"maxNodes\":None, \n", " \"seed\":0}\n", "\n", "clf2 = ee.Classifier.smileRandomForest(**new_params).train(train, label, trainingbands)" ] }, { "cell_type": "code", "execution_count": 14, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our classifier has an accuracy of 0.89286 on the test data.\n", "Our classifier has an MCC score of 0.60866 on the test data.\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAJrCAYAAAAMBPz+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA8JUlEQVR4nO3deZwcVbnw8d+TsArIKluCEDYRRVFZ1IsSRTYVUVHZRKLyRlRU3HFBcbuieN1RiIpBrhIX5AISDQqyCrJo2BcDBJgQ9h1BMpnn/aNqkk6na2Y6mZnu6f598+lPd1WdOnWqZibzzFPnnIrMRJIkqRuNa3UDJEmSWsVASJIkdS0DIUmS1LUMhCRJUtcyEJIkSV3LQEiSJHUtAyF1tYiYGxFzW92OsSwinh0R3y+vZW9EZERsP8LH3Kw8zvSRPE63Ka/p+a1uhzSaDIQ0qPI/x2ZeU0agDVNGqu7h0OW/QL4JfAi4Fvg68CXgnpa2qAtFxOTy+/CYVrdFGktWaHUDNCZ8qcG6I4E1ge8Bj9Rtmz2yzVGbeSNwS2buM4rHnAc8H3h0FI/ZDZ4P/LvVjZBGk4GQBpWZx9SvKzMzawLfzcy5o9wktZeNgQtH84CZuQC4aTSP2Q0y02uqruOtMQ27iNg5In4XEfdExDMRcVdEnBgRGzcou3lETIuIORHxVEQ8FBHXRsQJEbFuWeZ84OflLj+vuw232RDaExFxRERcHxFPR8S8iPhhRKxZUX7NiPhkRJwXET3lOdwfEWdGxMvryk6JiP7n1Oxa17Zj6sqdFhG3lef5WERcEhHvHMo1bdDG/SPi3PJ6PV32zzk1InaoK7dyRBwVEddExL/L414UEe9oUOeifjfl5xkR8UBZ/5UR8ca68ueX5x51535+7bWpup3Z6HZiRKwREUdHxHVlWx+PiFsj4tcR8bJGbW1Q70YRcXx5Tfq/dr+v3b+m7KI2RsRrynN6vDz22RHx/OqvwlJ1Lbo1FRE7RMSfIuLRiHi4/NpvUpbbvLy295ffC3+NiBc3qG/riDi2vPb3R8R/IuKO8udlYl3Z6cBfy8Uv1n0fTm5wrnuV5/pozffvUl+TiJgUEY+U32eb1h1ztYi4MSIWRsSuQ71OUrsxI6RhFRHvBn4C/Ac4E7gL2Ao4DNgnIl6emXeWZTcCrgCeDcwETgNWASYBhwA/BB4EplPcftsXOIMlb709MoRmfRf4MDAfmAYsKOvaGVgJeKau/POBr1FkOc4GHgaeC7wJ2Dsi9snMP5VlZ1PcOvwicEfZ1n7n13z+MXBDWed8YF3g9cApEfG8zDx6COdBRARFUHgo8ADwe+B+YCLwGuBm4Mqy7ErALGBXiuzJ8cCzgLcBv46I7TPzsw0OsylwOXAbcAqwDrA/cEZEvC4z+3/hTi/Psf7c5w7lXCrO7U/AK4FLgZ8CvcAmwGTgIuCqQeqYBFxMkaU6Dzi13P/twBsiYr/M/EODXd9I8T3xR+AEYFuKr8+OEbFtZj7QxKnsCHwauIDiZ2E74K3AdhHxprJ9NwG/oLjWbwX+HBGbZ+YTNfW8FTicIsD5G8X36QtY/LO0Q2bOK8v+X/l+aHnc82vqmVvXvrcBe9Wc62ZVJ5KZt0fEYcBvgVMj4tWZ2Vtu/hGwDXBMZl4w4BWR2llm+vLV9IviP9cENqtZtzXFf9ZzgAl15V8LLAROr1n3obKOjzSofzVg1ZrlKWXZKU2285XlfnOAdWrWr0LxyzaBuXX7rAms16CuicDdwI0NtiVw/gDt2KLBupWAcykCswlDPJ+p5bEuB9as2zYe2Khm+TNl2ZnACjXr16/5+r2yZv1m5boEvlhX9579dQ313Af7mtXvRxEwZO33SM22ccDaDdo6va7crHL95xp8H/RSBNarN2hjL7Bb3T5fL7d9aohfm8k11+/gum0/K9c/1KBtRzf6OQAmACs3OM4eFD9LP644/jEV7es/1z5gr6F8TWrW/6jc9vVy+V3l8l+BcUO5Pr58tevLW2MaTu8HVqT4D31e7YbMPI8iQ7RPRKxRt99T9RVl5pOZudT6ZfDu8v1rmflQTf1PUwQKS8nMR7NBBiAze4DfAdtExHObaURm3tpg3TMUWZoVgN2GWNWHyvf3ZeYSHYUzc2Fmzq9Z9R6KX1Yfy8V/xZOZ9wFfKRcPa3CMO4Cv1tU9C7gT2GmI7Vwejb4f+jLz4YF2Km8X7UHRzm/W7f83iuzQOhSZlnozMvPcunXTyvdmz/nizPxl3bqTy/dHgWPrtv2ifN++dmVmzsvM/9RXnpnnANdTBKfL4oxcnNEcqo8BVwOfjogjKAKj+ykCvr5lbIfUFrw1puH0ivJ914jYscH29SmyFltT3OI4E/hv4PiI2JPir/lLgBsyMxvsvyxeWr43St1fRJEJWEpE/BfwEYpzWp8ie1NrAsUv3CEpA6dPUwQ8zwVWbVDfYHWsBrwQuDcz/zlI2TWALYF52bgD7Hnl+0sabJudmQsbrL+LxV/jkXADxa3GA8v+KGdQ3Ea6sgwaB9N/Lhdl0Zm63nnAO8tyv6jbdmWD8neV72sP4diD1XV3+d7o2vb/0VDf7yeAgykyOS8u2zG+pshQrkkjlze7Q2Y+HRH7U5zbDygC7Ldl5t0D7ym1PwMhDad1y/dPDlJudYDMvCMidgKOoeiz0P+X+l0R8a3M/P4wtKm/Q/S99Rsyc2FEPFi/PiLeQpH5eRr4M3Ar8CTFLYXJFH1uVh5qAyJic4pfPmtTBF/nUGQGFlLc4jl0iPWtVb7PG6hQqf+851ds71+/VoNtj1Ts08sIDrAovx6vBb5A0Y/lG+WmxyPiZOAzuWQfmnrDes6Z2VvEIksEH0PRaEh/b9W2muOsWLfp2xTTVMyn+CNhHouzZVMo+hcti2Wd4+kW4BqK24w3UHwfS2OegZCGU/9/8mtm5mND2SEzbwT2j4gVKP7qfR3F7Z/vRcSTmfmzYWrTBhSdfxeJiPEUwVt9YPEVir+2dyjbV7vPiRSBUDM+Vh7n3Zk5va6+AykCoaF4pHwfNHvE4vPesGL7RnXlRkL/LZOl/p+JiLUa7VDe/voo8NGI2JLiWr8POIIigDlkgOO1wzkPi4hYn6KD/3UU/bger9t+4HJUv6zZ1qMogqAHKDptf4ZiUIE0ptlHSMPpsvL9Vc3umJm9mXlVZn4D6P9P/s01RfpvJzT71/k/yvdGwcuraPzHwJYUt+fqg6BxwC4Vx+kboG1blu+nNdg25KAqM5+k+MW4QUQ0uqVVW/ZxikzWhIjYqkGR15Tv/2iwbbj09+nZpMG2HRqsW0JmzikD4V2BJyhGdQ2k/3bhLmVgXW80znm4bE7x//M5DYKgieX2esv6MzKoiHgl8GWKUYkvLN+/FBFVPw/SmGEgpOH0Q4oRUN+JiK3rN0bEShHxqprlnSJigwb19K+rneG2/xZWU52UWTyk+3MRsU7NsVehGBXUyFxgq6iZ96jsr/FFimHVjTxI41/4/fVBcVttkbJfVKPOygPpv114YtTNgxQR48opCfqdRDHHz3Fl9qu/3HoUI5X6y4yUKykCxIMi4lk1x1+Hus7M5fpJEfGCBvWsTXHrcMDO82Vn9j9T3G48sq7unYGDKIKz05s5iRaZW77vUve1W51iSH6jQG9Zf0YGFBFrU3Q0XwgckJn3Ukyn0EsxpH7dgfaX2p23xjRsMvOmiHgPxS/X6yPiTxT9Clak+M/5VRQjTbYpdzkI+GBEXEAxvP1hYAtgH4p5iL5bU/2lFIHRkeUv0v4+Pz+oHz1V16ZLIuIHFLfbrouI37F4HqGHadyf5DsU86v8MyJOK8v/F0UQdFbZvnrnAgdExFkUHcF7gQsz80KKETbvBn5b1jeP4q/qvYDfUPxSGaqfUmSl3gX8KyLOoLimG1NMUXASRZ8rgG8Be5fnenVEzKSYR+jtFB3Av5mZFzdx7KZk5vyI+CXF7azZEXE2xZxRr6eYT6k+q/Vi4PSIuIoi83U38Jyy/SuyuM/QQA6n6HB/XETsQRGM9c8j1Edxe/LxAfZvC5l5T0TMAA6guHbnUPSB2p2i79ps6kaZUWRp5lF8Hz5D0Zk/gVMy847laM5JFD+/H87M2WX7ro6Ij1P88fNzijm2pLGp1eP3fY3NFw3mEarZth1FJuYOioDmIYpfbCcCr60ptzPFRINXl2WeogiIfg68sEG9e1EERE+weL6WpY7fYL+g6GNyY9meuymGra9ZnsfcBvtMofhl8yRFn4jTy/M6pjzu5Lry6wO/ogjQFlI3nwtF34rzKIKvxylGQ72ZQeZ+GeCcDqYYCfcoxS/G24FfAi+tK7cK8Nny+j9Vc+wDG9S5GQ3m5qnZfn7xX8ZS6yvnUKLI5BwH9LB4jqnPUPwRVj+P0ESKUYSXUHTo/U+53x+BvYfaVoo+VD8uv/+eKb9+/wfsWPF1HvJcR4N8TSq/lkO4tksdhyJo/Vp5zZ6mGMV2PEV/s6qvxY4UQfmjFIHfou/Vwc61UTtYPNfXGRXlf19u/2gz37++fLXTKzKHa5SyJEnS2GIfIUmS1LUMhCRJUtcyEJIkSV3LQEiSJHUtAyFJktS1xvQ8QgseuM0hb1ILrDFxcqubIHWtp5++M0bzeKP5u3bF9TYf1XMDM0KSJKmLGQhJkqSuNaZvjUmSpBHWt3DwMmOYGSFJktS1zAhJkqRq2dfqFowoM0KSJKlrmRGSJEnV+swISZIkdSQzQpIkqVLaR0iSJKkzmRGSJEnV7CMkSZLUmcwISZKkavYRkiRJ6kwGQpIkqWt5a0ySJFXzoauSJEmdyYyQJEmqZmdpSZKkzmRGSJIkVXNCRUmSpM5kRkiSJFXyoauSJEkdyoyQJEmqZh8hSZKkzmRGSJIkVbOPkCRJUmcyIyRJkqr5rDFJkqTOZEZIkiRVs4+QJElSZzIQkiRJXctbY5IkqZoTKkqSJHUmM0KSJKmanaUlSZI6kxkhSZJUzT5CkiRJncmMkCRJqpTpIzYkSZI6khkhSZJUzVFjkiRJrRcRe0XEzRExJyKOarB97Yg4PSKuiYjLI+KFg9VpRkiSJFVrk1FjETEeOB7YHegBroiIMzPzhppinwVmZ+ZbImKbsvxuA9VrRkiSJI0FOwFzMvO2zHwGmAHsW1dmW+BcgMy8CdgsIjYYqFIDIUmSVC37Ru81sAnAXTXLPeW6WlcDbwWIiJ2ATYGJA1VqICRJktpCREyNiCtrXlNrNzfYJeuWjwXWjojZwIeAfwK9Ax3TPkKSJKla3+jNI5SZ04BpFZt7gE1qlicCd9ft/xjwboCICOD28lXJjJAkSRoLrgC2iohJEbEScABwZm2BiFir3AZwGHBhGRxVMiMkSZLaXmb2RsQRwCxgPHBSZl4fEYeX208Ang/8IiIWAjcA7x2sXgMhSZJUrY0mVMzMmcDMunUn1Hy+FNiqmTq9NSZJkrqWGSFJklStTSZUHClmhCRJUtcyIyRJkqq1UR+hkWBGSJIkdS0zQpIkqZp9hCRJkjqTGSFJklTNjJAkSVJnMiMkSZIqZY7eQ1dbwYyQJEnqWmaEJElSNfsISZIkdSYzQpIkqZozS0uSJHUmAyFJktS1vDUmSZKq2VlakiSpM5kRkiRJ1ewsLUmS1JnMCEmSpGr2EZIkSepMZoQkSVI1+whJkiR1JjNCkiSpmn2EJEmSOpMZIUmSVM2MkCRJUmcyIyRJkqo5akySJKkzmRGSJEnV7CMkSZLUmQyEJElS1/LWmCRJqmZnaUmSpM5kRkiSJFWzs7QkSVJnMiMkSZKq2UdIkiSpM5kRkiRJ1ewjJEmS1JnMCEmSpGpmhCRJkjqTGSFJklQts9UtGFFmhCRJUtcyIyRJkqrZR0iSJKkzmRGSJEnVzAhJkiR1JjNCkiSpms8akyRJ6kwGQpIkqWt5a0ySJFWzs7QkSVLrRcReEXFzRMyJiKMabF8zIs6KiKsj4vqIePdgdZoRkiRJ1drkERsRMR44Htgd6AGuiIgzM/OGmmIfBG7IzH0i4jnAzRHxy8x8pqpeM0KSJGks2AmYk5m3lYHNDGDfujIJrBERAawOPAT0DlSpGSFJklStffoITQDuqlnuAXauK/ND4EzgbmANYP/Mgcf/mxGSJEltISKmRsSVNa+ptZsb7FJ/325PYDawMbA98MOIePZAxzQjJEmSqo1iRigzpwHTKjb3AJvULE+kyPzUejdwbGYmMCcibge2AS6vOqYZIUmSNBZcAWwVEZMiYiXgAIrbYLXuBHYDiIgNgOcBtw1UqRkhSZJUrU0esZGZvRFxBDALGA+clJnXR8Th5fYTgK8A0yPiWopbaZ/OzAcGqtdASJIkjQmZOROYWbfuhJrPdwN7NFOngZAkSaqUfe0xj9BIsY+QJEnqWmaEJElStfaZR2hEmBGSJEldy4yQJEmq1iajxkaKGSFJktS1DIQkSVLX8taYJEmq5vB5SZKkzmRGSJIkVXP4vCRJUmcyIyRJkqqZEZIkSepMZoQkSVK1dNSYJElSRzIjJEmSqtlHSJIkqTOZEZIkSdWcWVpqzsWXXckbDziMvd/xHn56ym+W2v7oY4/z4c98mbe86/0ccNhH+NdtcwG4/Y4e9jv0g4teO+/+Vk759emj3Hpp7Np991255pq/cv31F/KJT3xgqe1bb70F559/Oo8++i+OPHLqovUTJ27ErFkzmD37XP7xj7/wwQ++ZzSbLbWUGSENq4ULF/LV/zmen3z3v9lw/fXY/7CP8JpddmaLSZsuKvOTX/yabbbagu9//QvcdsddfO1/judn3z+WSZtO5LSTj19Uz2vffAi77frKVp2KNKaMGzeO733vq7zhDQfT0zOfSy45iz/84c/cdNO/FpV5+OFH+PjHv8ib3rTnEvv29i7k05/+KrNnX8fqq6/GpZeezbnnXrTEvupiaR+hERURKw9lncaGa2+8hedO3JhNJmzEiiuuyN677cp5F122RJlb597Jy1/2YgA233QT5s2/lwceeniJMpddOZtNJmzExhtuMGptl8ayHXfcnltvncvtt9/JggUL+O1vz2KfffZYosz99z/IVVddw4IFvUusv+ee+5g9+zoAnnjiSW66aQ4TJmw4am2XWqnlgRBw6RDXaQy47/4H2HD95yxa3mD99bjv/geXKPO8LTfnLxf8DYBrb7iZ+ffex733PbBEmT+eewGvf92uI99gqUNsvPGG9PTcvWh53rz5bLxx839IbLrpRLbf/gVcfvk/h7N5Gsv6cvReLdCyQCgiNoyIlwGrRsRLIuKl5Wsy8KxWtUvLp9G8WxFLLh92yNt57PEn2O/QD/LL353JNlttwfjx4xdtX7BgAedf/Hf2eO2rRri1UueI+h80IJucCG+11Z7FqaeeyCc+8SUef/yJ4Wqa1NZa2UdoT2AKMBH4ds36x4HPVu0UEVOBqQA/+p+vcti7DhzBJqpZG6y/Hvfcd/+i5Xvve4DnrLfuEmVWX201vvq5jwHFf9R7vm0KE2v+cr3osit5/tZbsN46a49Oo6UOMG/efCZO3HjR8oQJGzF//n1D3n+FFVZgxowTmTHjdM44408j0USpLbUsEMrMk4GTI2K/zDytif2mAdMAFjxwW2eP6RuDXrjN1tzZczc9d9/DBs9Zlz+eewHf/OKnlyjz2ONPsOoqK7Piiity2ll/4mXbb8fqq622aPvMP5/P63efPMotl8a2K6+8mi23nMRmm23CvHn38Pa378Ohh354yPufeOJx3HTTHL7//Z+OYCs1FmWHT6jYDqPGzo2IbwOvLpcvAL6cmY+2sE1aRiusMJ7PfvT9vO9jn2fhwoW85Y17sOXmm/Lr088GYP+3vIHb7riLz37lW4wfN47NN3suX/7MkYv2f+rpp7n0in/yxU8N/T9wScVIyyOPPJqzzjqF8ePHc/LJv+bGG2/hsMPeCcBPf/q/bLDBc7jkkj/w7GevTl9fH0cc8V5e8pLd2G6753Pwwftx7bU38ve//xGAL3zhm8ya9ddWnpI0KqLZe8jD3oCI04DrgJPLVYcAL87Mtw62rxkhqTXWmDi51U2QutbTT9+5dIewEfTk1941ar9rV/vcL0b13KA9MkJbZOZ+NctfiojZrWqMJEnqHu0QCD0VEbtk5sUAEfFfwFMtbpMkSYKOn1CxHQKh91N0ml6zXH4YOLSF7ZEkSV2iHQKhG4FvAlsAawGPAm8GrmldkyRJEtDxD11th0DoDOAR4B/AvNY2RZIkdZN2CIQmZuZerW6EJElqoMPnEWqHZ439LSK2a3UjJElS92mHjNAuwJSIuB34DxBAZuaLWtssSZJkH6GRt3erGyBJkrpTywOhzLyj1W2QJEkVOnweoXboIyRJktQSLc8ISZKkNtbhfYTMCEmSpK5lICRJkrqWt8YkSVKldEJFSZKkzmRGSJIkVbOztCRJUmcyIyRJkqqZEZIkSepMZoQkSVI1H7EhSZLUmcwISZKkavYRkiRJ6kxmhCRJUqU0IyRJktSZzAhJkqRqZoQkSZJaLyL2ioibI2JORBzVYPsnI2J2+bouIhZGxDoD1WlGSJIkVWuTp89HxHjgeGB3oAe4IiLOzMwb+stk5nHAcWX5fYCPZuZDA9VrRkiSJI0FOwFzMvO2zHwGmAHsO0D5A4FTB6vUQEiSJI0FE4C7apZ7ynVLiYhnAXsBpw1WqbfGJElStVHsLB0RU4GpNaumZea0/s0Ndqlq3D7AJYPdFgMDIUmS1CbKoGdaxeYeYJOa5YnA3RVlD2AIt8XAQEiSJA2kfYbPXwFsFRGTgHkUwc5B9YUiYk1gV+CdQ6nUQEiSJLW9zOyNiCOAWcB44KTMvD4iDi+3n1AWfQtwTmY+OZR6DYQkSVKlzLbJCJGZM4GZdetOqFueDkwfap2OGpMkSV3LjJAkSarWPn2ERoQZIUmS1LXMCEmSpGpmhCRJkjqTGSFJklQpzQhJkiR1JjNCkiSpmhkhSZKkzmRGSJIkVetrdQNGlhkhSZLUtQyEJElS1/LWmCRJquTweUmSpA5lRkiSJFUzIyRJktSZzAhJkqRqDp+XJEnqTGaEJElSJUeNSZIkdSgzQpIkqZp9hCRJkjqTGSFJklTJPkKSJEkdyoyQJEmqZh8hSZKkzmRGSJIkVUozQpIkSZ3JQEiSJHUtb41JkqRq3hqTJEnqTGaEJElSJTtLS5IkdSgzQpIkqZoZIUmSpM5kRkiSJFWyj5AkSVKHMiMkSZIqmRGSJEnqUGaEJElSJTNCkiRJHcqMkCRJqpbR6haMKDNCkiSpa5kRkiRJlewjJEmS1KEMhCRJUtfy1pgkSaqUfXaWliRJ6khmhCRJUiU7S0uSJHUoM0KSJKlSOqGiJElSZzIjJEmSKtlHSJIkqQ1ExF4RcXNEzImIoyrKTI6I2RFxfURcMFidZoQkSVKldplHKCLGA8cDuwM9wBURcWZm3lBTZi3gR8BemXlnRKw/WL1mhCRJ0liwEzAnM2/LzGeAGcC+dWUOAn6fmXcCZOZ9g1VqICRJkipljt5rEBOAu2qWe8p1tbYG1o6I8yPiqoh412CVemtMkiS1hYiYCkytWTUtM6f1b26wS334tALwMmA3YFXg0oi4LDNvqTqmgZAkSao0mn2EyqBnWsXmHmCTmuWJwN0NyjyQmU8CT0bEhcCLgcpAyFtjkiRpLLgC2CoiJkXESsABwJl1Zc4AXhURK0TEs4CdgRsHqtSMkCRJqtQuo8YyszcijgBmAeOBkzLz+og4vNx+QmbeGBF/Aq4B+oCfZuZ1A9VrICRJksaEzJwJzKxbd0Ld8nHAcUOt01tjkiSpa5kRkiRJlYYwrH1MMyMkSZK6lhkhSZJUqV06S48UM0KSJKlrmRGSJEmVMs0ISZIkdSQzQpIkqVL2tboFI2vAQCgi+lj6gWZDkZlpkCVJktraYMHKhSxbICRJkjpAX4f3ERowEMrMyaPUDkmSpFHn7StJklTJUWOSJEkdapkyQhGxEbAbMAFYuUGRzMyvLE/DJElS63X6zNJNB0IR8SXgqLp9g8Wdqvs/GwhJkqS21tStsYg4GDgauAh4G0XQczJwEPAToA+YAbx2eJspSZJaIXP0Xq3QbEbo/UAPsFdm9kYEwNzMnAHMiIjTgbOBU4e3mZIkScOv2c7S2wEzM7O3Zt34/g+ZOQuYBXxyGNomSZI0oprNCK0IPFiz/BSwZl2Z64DDl6dRkiSpPXR6Z+lmM0LzgY1qlu8EXlRXZgLQiyRJUptrNhD6J8XtsX7nAa+KiEMiYrWIeAOwX1lOkiSNcX0Zo/ZqhWYDoT8AL4iISeXyscCjwHTgMeBMipFknx+uBkqSJI2UpvoIZeZ0iqCnf/muiNgR+DiwBTAX+FFmXjt8TZQkSa3S6Y/YWO5njWXm7cARw9AWSZKkUeVDVyVJUqVWTXQ4WpoKhCLiuUMtm5l3Nt8cSZKk0dNsRmgui58pNpBchrolSVKbadVortHSbLDyCxoHQmsB2wObAucDdyxPoyRJkkZDs6PGplRti4hxFA9kPRw4dPmaJUmS2kGnjxprdh6hSpnZl5lforh9duxw1StJkjRSRqIfz9+Ad41AvZIkaZR1+qixYcsI1VgHWG0E6pUkSRpWw5oRiojXAftTPIFekiSNcY4aqxER5w1QzyZA/zxDX16eRkmSJI2GZjNCkyvWJ/AwMAv4VmZWBUzDatWNXzUah5FU57L1d2x1EySNkk4fNdbs8PmR6FMkSZLUEgY2kiSpazUVCEXEeREx4ND4iHjnAH2JJEnSGNKXMWqvVmg2IzQZ2GyQMpsCuy5LYyRJkkbTSEyouCrQOwL1SpKkUdbh8ykuUyDU8JpERFAMn389cNfyNEqSJGk0DBoIRUQfSwY/x0TEMQPtAvz3crZLkiS1ASdUhAtZHAi9GriT4sGq9RYCDwLnAj8djsZJkiSNpEEDocyc3P+5zA79PDOdOVqSpC7ghIpLmgQ8MgLtkCRJGnXNBkL3Ac+JiKcy85n6jRGxMrABcF9mPj0cDZQkSa3T1+oGjLBm5xH6AnAzsHrF9tWAm4DPLk+jJEmSRkOzgdDewF8y86FGG8v1fwHeuLwNkyRJrZfEqL1aodlAaDPglkHK3MLgs09LkiS1XLN9hFZk8NuFCayybM2RJEntpK/Dp5ZuNiN0G4M/R2wycMcytUaSJGkUNRsInQm8LCI+1WhjRBwFvBT4v+VslyRJagN9xKi9WqHZW2PfAg4Gvh4R7wDOAeYBE4A9ge0pZp7+5jC2UZIkaUQ0FQhl5sMRMRn4JfAKiuxPwqIw7m/AOzPz4WFsoyRJ0oho+unzmTkX+K+IeCnwcmAtitmmL8vMfwxn4yRJUmu1alh7IxGxF/A9YDzw08w8tm77ZOAM4PZy1e8HeyxY04FQvzLoMfCRJEkjLiLGA8cDuwM9wBURcWZm3lBX9KLMHPJ8hssUCEXERsBuFH2DVm5QJDPzK8tStyRJah9t9IiNnYA5mXkbQETMAPYF6gOhpjQdCEXEl4Cj6vYNir5CtZ8NhCRJ0nCZANxVs9wD7Nyg3Csi4mrgbuATmXn9QJU2NXw+Ig4GjgYuAt5GEfScDBwE/IQicJwBvLaZeiVJUnsazUdsRMTUiLiy5jW1pimNOivVT/f4D2DTzHwx8AOGMJ1Psxmh91NEYHtlZm9EAMzNzBnAjIg4HTgbOLXJeiVJUpfLzGnAtIrNPcAmNcsTKbI+tfs/VvN5ZkT8KCLWy8wHqo7Z7ISK2wEzM7O3Zt34moPOAmYBn2yyXkmS1Ib6RvE1iCuArSJiUkSsBBxAMdHzIhGxYZRZmojYiSLOeXCgSpflWWO1FT4FrFlX5jrg8CbrlSRJqlTeiTqCIuEyHjgpM6+PiMPL7SdQdNt5f0T0UsQoB2TmgE9LazYQmg9sVLN8J/CiujITgF4kSdKY10ajxsjMmcDMunUn1Hz+IfDDZups9tbYPyluj/U7D3hVRBwSEatFxBuA/cpykiRJba3ZQOgPwAsiYlK5fCzwKDAdeIziXl0Anx+uBkqSpNYZzVFjrdDss8amUwQ9/ct3RcSOwMeBLYC5wI8y89rha6IkSdLIWOZHbPTLzNuBI4ahLZIkqc30tc+jxkZEs7fGJEmSOsZyZ4QkSVLn6mujp8+PBDNCkiSpaxkISZKkruWtMUmSVGnAaZk7gBkhSZLUtcwISZKkSu30iI2RYEZIkiR1LTNCkiSpUl84fF6SJKkjmRGSJEmVHDUmSZLUocwISZKkSo4akyRJ6lBmhCRJUqW+zh40ZkZIkiR1LzNCkiSpUh+dnRIyIyRJkrqWGSFJklTJeYQkSZI6lIGQJEnqWt4akyRJlRw+L0mS1KHMCEmSpEo+YkOSJKlDmRGSJEmVHD4vSZLUocwISZKkSo4akyRJ6lBmhCRJUiVHjUmSJHUoM0KSJKmSGSFJkqQOZUZIkiRVSkeNSZIkdSYzQpIkqZJ9hCRJkjqUgZAkSepa3hqTJEmVvDUmSZLUocwISZKkStnqBowwM0KSJKlrmRGSJEmV+pxQUZIkqTOZEZIkSZUcNSZJktShzAhJkqRKZoQkSZI6lBkhSZJUyXmEJEmSOpSBkCRJqtQXo/caTETsFRE3R8SciDhqgHI7RsTCiHjbYHUaCEmSpLYXEeOB44G9gW2BAyNi24py3wBmDaVeAyFJklSpbxRfg9gJmJOZt2XmM8AMYN8G5T4EnAbcN5TzMxCSJEljwQTgrprlnnLdIhExAXgLcMJQKzUQkiRJbSEipkbElTWvqbWbG+xSP6jtu8CnM3PhUI/p8HlJklRpNIfPZ+Y0YFrF5h5gk5rlicDddWV2AGZEBMB6wOsjojcz/6/qmAZCkiRpLLgC2CoiJgHzgAOAg2oLZOak/s8RMR34w0BBEBgISZKkAfS1yZSKmdkbEUdQjAYbD5yUmddHxOHl9iH3C6plICRJksaEzJwJzKxb1zAAyswpQ6nTQEiSJFXyoauSJEkdyoyQJEmq1B49hEaOGSFJktS1zAhJkqRK9hGSJEnqUGaEJElSpb5GD7boIGaEJElS1zIjJEmSKrXLzNIjxYyQJEnqWmaEJElSpc7OB5kRkiRJXcxASJIkdS1vjUmSpEpOqChJktShzAhJkqRKDp+XJEnqUGaEJElSpc7OB5kRkiRJXcyMkCRJquSoMUmSpA5lRkiSJFVy1JgkSVKHMiMkSZIqdXY+yIyQJEnqYmaEJElSJUeNSZIkdSgzQpIkqVJ2eC8hM0KSJKlrGQhJkqSu5a0xSZJUyc7SkiRJHcqMkCRJquQjNiRJkjqUGSFJklSps/NBZoQkSVIXMyMkSZIq2UdIkiSpQ5kRkiRJlZxHSFoOe+4xmeuvu5CbbriYT33yg0ttf97ztuDiC8/kycdv42MffV8LWih1jmdPfgkvvOB4Xnjxj9nwg29davsar3gh29/wS7ad9R22nfUdNjryHUsWGDeObf/0bbac/rlRarHUemaENGLGjRvH97/3NfZ6/YH09MznsktnctYfzuHGG/+1qMxDDz3CkR89mn333auFLZU6wLhxPPer7+OWg77IgvkP8vyzj+ORcy7n6X/1LFHsictvYM6UrzWsYoP3vpGn5vQwfvVVR6PFGiN86Kq0jHba8SXceutcbr/9ThYsWMBvfnMGb9pnzyXK3H//g1x51dUsWLCgRa2UOsNq22/Ff+bO55k77yUX9PLQGRez1h47D3n/FTdalzV324EHfvXnEWyl1H5aHghFxOYRcVZEPBAR90XEGRGxeavbpeW38YQNuavn7kXLPfPms/HGG7awRVLnWmmjdXhm/gOLlp+550FW2midpcqt/rLnse0532GrU45mla03WbR+k2PeS8/XTobs7L/+1by+UXy1QssDIeBXwG+ADYGNgd8Cp7a0RRoWEbHUuvQ/WWmENPp5W3L5yWtv5Zqdp3LDHh/lvp/PZMuffQaANXfbgd4HHuXf1946Gg2V2ko7BEKRmadkZm/5+l8GmMgyIqZGxJURcWVf35Oj2Ew1a17PfDaZuPGi5YkTNmL+/Htb2CKpcz0z/0FW2mi9RcsrbbguC+55aIkyfU88Rd+/nwbg0fOuIlZYgRXWXoPVd9yGtfbYke0uncbmx3+cNf7rRUz6/pGj2Xy1sRzFf63QDoHQXyPiqIjYLCI2jYhPAWdHxDoRsVReNzOnZeYOmbnDuHGrtaC5GqorrpzNlltOYrPNNmHFFVfkHe/Yl7P+cE6rmyV1pCev/herTNqIlTZZn1hxBdbZdxce+fPlS5RZ4TlrLfq82vZbwbig9+HHmXfs/3LNjodx7SumctsH/4fHL7mG2z/83dE9AalF2mHU2P7le/3Y6fdQZIbsLzRGLVy4kI8c+Xlmnv0rxo8bx/STf80NN9zC1P93CADTfnIKG2zwHP5+6R959rNXp6+vjw9/6P+x3Ysn8/jjT7S49dIYs7CPO4/+CVv/8oswbjwP/vovPH3LXTznncUAhfv/dxZrv+GVrH/IXuTChfQ9/Qy3feBbLW601HoxlvtsrLDShLHbeGkMu2z9HVvdBKlr7dDzf0t3CBtBh26236j9rj157mmjem7QBhmhiHhXo/WZ+YvRboskSeouLQ+EgNo/LVcBdgP+ARgISZLUYn1j+M7RULQ8EMrMD9UuR8SawCktao4kSeoiLQ+EGvg3sFWrGyFJkgaYz6ZDtDwQioizWHydxwPPp5hgUZIkaUS1PBACasdv9gJ3ZGZPVWFJkjR6+tooJxQRewHfo0ic/DQzj63bvi/wFYondvQCR2bmxQPV2fJAKDMv6P8cEW/MzEta2R5JktR+ImI8cDywO9ADXBERZ2bmDTXFzgXOzMyMiBdR3GHaZqB622Fm6VpfbnUDJEnSYm30iI2dgDmZeVtmPgPMAPZdoq2ZT+TiCRJXYwhdnNotEBr1iZQkSdKYMAG4q2a5p1y3hIh4S0TcBJxN8ZSKAbU8EIqIlWsW39dgnSRJapG+UXzVPli9fE2taUqjZMlSGZ/MPD0ztwHeTNFfaEAt7yMEXAq8FCAzL69fJ0mSukNmTgOmVWzuATapWZ4I3D1AXRdGxBYRsV5mPlBVrmWBUERsSJHSWjUiXsLiSO/ZwLNa1S5JkrRYG40auwLYKiImAfOAA4CDagtExJbArWVn6ZcCKwEPDlRpKzNCewJTKCK6b9esfxz4bCsaJEmS2lNm9kbEEcAsiuHzJ2Xm9RFxeLn9BGA/4F0RsQB4Ctg/B3m6fMsCocw8GTg5IvbLzNNa1Q5JklRtCKO5Rk1mzgRm1q07oebzN4BvNFNnK2+NfazR536Z+e36dZIkScOplbfG1mjhsSVJklp6a+xLrTq2JEkamr5WN2CEtfLW2Kcy85sR8QMazwPw4RY0S5IkdZFW3hq7sXy/soVtkCRJAxhk0NWY18pbY2eV7ye3qg2SJKm7tXxm6Yj4K41vjb22Bc2RJEk12mhCxRHR8kAI+ETN51UoJkPqbVFbJElSF2l5IJSZV9WtuiQiLmhJYyRJ0hIcNTbCImKdmsVxwA7Ahi1qjiRJ6iItD4SAq1jcR6gXmAu8t2WtkSRJi7TTIzZGQjsEQtsCHwB2oQiILsIh9ZIkaRS0QyB0MvAY8P1y+UDgFODtLWuRJEkCHDU2Gp6XmS+uWf5rRFzdstZIkqSu0Q6B0D8j4uWZeRlAROwMXNLiNkmSJJxZesRExLUUfYJWBN4VEXeWy5sCN7SqXZIkqXu0MiP0xhYeW5IkDYHzCI2QzLyjVceWJEmC9ugjJEmS2lSnzyM0rtUNkCRJahUDIUmS1LW8NSZJkip1+oSKZoQkSVLXMiMkSZIqdfqEimaEJElS1zIjJEmSKtlHSJIkqUOZEZIkSZWcUFGSJKlDmRGSJEmV+hw1JkmS1JnMCEmSpEqdnQ8yIyRJkrqYGSFJklTJeYQkSZI6lBkhSZJUyYyQJElShzIQkiRJXctbY5IkqVI6oaIkSVJnMiMkSZIq2VlakiSpQ5kRkiRJldKMkCRJUmcyIyRJkio5akySJKlDmRGSJEmVHDUmSZLUocwISZKkSvYRkiRJ6lBmhCRJUiX7CEmSJHUoM0KSJKmSM0tLkiR1KAMhSZI0JkTEXhFxc0TMiYijGmw/OCKuKV9/i4gXD1ant8YkSVKlvjYZPh8R44Hjgd2BHuCKiDgzM2+oKXY7sGtmPhwRewPTgJ0HqteMkCRJGgt2AuZk5m2Z+QwwA9i3tkBm/i0zHy4XLwMmDlapGSFJklSpjTpLTwDuqlnuYeBsz3uBPw5WqYGQJElqCxExFZhas2paZk7r39xgl4ZRWkS8hiIQ2mWwYxoISZKkSqPZR6gMeqZVbO4BNqlZngjcXV8oIl4E/BTYOzMfHOyY9hGSJEljwRXAVhExKSJWAg4AzqwtEBHPBX4PHJKZtwylUjNCkiSpUrv0EcrM3og4ApgFjAdOyszrI+LwcvsJwBeAdYEfRQRAb2buMFC9BkKSJGlMyMyZwMy6dSfUfD4MOKyZOg2EJElSpXaZR2ik2EdIkiR1LTNCkiSpUrv0ERopZoQkSVLXMiMkSZIq2UdIkiSpQ5kRkiRJlewjJEmS1KEMhCRJUtfy1pgkSaqU2dfqJowoM0KSJKlrmRGSJEmV+uwsLUmS1JnMCEmSpErphIqSJEmdyYyQJEmqZB8hSZKkDmVGSJIkVbKPkCRJUocyIyRJkir1mRGSJEnqTGaEJElSpXTUmCRJUmcyIyRJkio5akySJKlDGQhJkqSu5a0xSZJUyUdsSJIkdSgzQpIkqZKdpSVJkjqUGSFJklTJR2xIkiR1KDNCkiSpkn2EJEmSOpQZIUmSVMl5hCRJkjqUGSFJklTJPkKSJEkdyoyQJEmq5DxCkiRJHcqMkCRJqpSOGpMkSepMBkKSJKlreWtMkiRVsrO0JElShzIjJEmSKjmhoiRJUocyIyRJkio5fF6SJKlDmRGSJEmV7CMkSZLUocwISZKkSmaEJEmSOpQZIUmSVKmz80FmhCRJUheLTr/3p/YVEVMzc1qr2yF1G3/2pMXMCKmVpra6AVKX8mdPKhkISZKkrmUgJEmSupaBkFrJPgpSa/izJ5XsLC1JkrqWGSFJktS1DITUViJiSkRs3Op2SO0oIjaLiOuaKP+miDiq/HxMRHyi/OzPmVQyEFK7mQL4H7Q0DDLzzMw8tsGmKfhzJgEGQhph5V+wN0bETyLi+og4JyJWjYjtI+KyiLgmIk6PiLUj4m3ADsAvI2J2RKza6vZLbWiFiDi5/Nn5XUQ8KyLmRsR6ABGxQ0ScX36eEhE/rN15sJ+ziJgcEX+oWf5hREwpP8+NiG9ExOXla8uRPVVp5BkIaTRsBRyfmS8AHgH2A34BfDozXwRcC3wxM38HXAkcnJnbZ+ZTrWqw1MaeB0wrf3YeAz7QzM7D8HP2WGbuBPwQ+G6T+0ptx0BIo+H2zJxdfr4K2AJYKzMvKNedDLy6FQ2TxqC7MvOS8vP/AruM8vFPrXl/xSgfWxp2BkIaDf+p+bwQWKtF7ZA6Qf2cJwn0svj/81WaqSwidi5vkc2OiDfV1dWovqz4LI1JBkJqhUeBhyPiVeXyIUB/duhxYI2WtEoaG54bEf2ZmAOBi4G5wMvKdfsNoY5FP2eZ+ffyFtn2mXkmcAewbUSsHBFrArvV7bt/zfuly34aUntYodUNUNc6FDghIp4F3Aa8u1w/vVz/FPAK+wlJS7kRODQiTgT+BfwYuBz4WUR8Fvj7EOqYTsXPWWbeFRG/Aa4p6/9n3b4rR8TfKf6QPnB5T0ZqNWeWliQNSUTMBXbIzAda3RZpuHhrTJIkdS0zQpIkqWuZEZIkSV3LQEiSJHUtAyFJktS1DISkURYR2f8sqOWoY7OynunD06rWa3RdyiemZ0RMHqFjdtx1lNQcAyFJHW04Ak9JncsJFSW1sx8CM4A7R6j+ecDzKWY7l9SFDIQkta1y4r4Rm7wvMxcAN41U/ZLan7fG1HFq+31ExBYR8buIeDAiHo+IcyLihWW550TEtIiYHxFPR8QVEfGaijrXjIivR8TNZdmHI2JWRLyuovxKEXF0RNwaEf+JiNsj4qsRsfIA7V4hIj4QEZdFxGMR8e+I+GdEHBERy/WzGhGTy2tyTES8IiL+EhGPltdkVkTs0GCfRf1zIuKgiPh7RDxRzi7cX+ZZEfGZ8oGdT5bbL42Iho9eaPa6DNRHKCK2iYiTImJuWdd9EXFRRLy/3D4lIvonStu1rKf/dUxZprKPUERsFBHHl/U/ExH3R8TvI+JlDcpOKeuZEhGviYjzy2v7WEScHRHPb7DPBhHxrfJ76smIeKT8PD0iNm90PSQNPzNC6mSbUTx36UaKZyttBrwFOD+Kh1b+CXgM+DWwDnAA8MeI2DozF92KiYi1gEuAbYErgO8C6wHvAM6JiPdn5ok15QP4DbAvcCvF7Z2VgPcA2zVqaESsCJwF7AncDPwKeBp4DfADYGeKh9Mur52BzwB/AY4HtgTeCrw6IvbIzIsa7PNxYPeyfX8F1izbvBZwHvAS4B/ASRR/XO0J/CoiXpCZn685x6avS5WIeAPwW2Bliq/jqcBawIuBT1E8f2s28CXgixQPEp1eU8X5g9Q/ieJhphuX53gqsAnwduANEbFfZv6hwa5vLM/vj8AJFN8zrwd2jIht+x9NEcUz9i4BtgD+THFtA9i03P93FM/gkzTSMtOXr456UQQ8Wb4+V7ft6HL9QxS/qMbVbDuk3Padun1OLNefSDkbe7l+K4q+Jf8BNqtZf1BZ/lJglZr161AEAAmcX3eMY8r1PwDG16wfD/ys3LZvg3OcPsRrMrnmmhxRt23fcv2/6q5Hf5ueBF7SoM7p5fZP1a1fhSI46QO2H6brMrlm3XrldX8G2LVBuybWLS9V72DXEZhV8f3zSqAXeBBYvWb9lLJ8L7Bb3T5fr79OwD6NvtfKbSsBa7T658iXr255eWtMnWwucGzdupPL95WBT2ZmX822X1H8Itu+f0WZqXkn8ATwmcxc9EyazPwX8H2KX1zvqqnn3eX7ZzPz6ZryDwFfqW9kedvrCOAe4KOZubBmn4UUGZkEDh7shIdgDvCj2hWZeQZwAUV26FUN9pmWmUs8gTwi1qW4Lldm5jfr6nsa+DRFhuOgmk1NXZcBHAo8G/hxZl5QvzEze5qoaykRMRHYg6KDdv25/Y0iO7QORSat3ozMPLdu3bTyfacG5Z+qX5GZz2Tm4822W9Ky8daYOtns2qCidHf5fkv9L5vMXBgR9wITa1ZvAzwLuKT8hV3vPODzFLeH+r2UIhtycYPy5zdYtzWwLkVG5vPFHaSlPEUxuml5XVQX/NW2a1eK86gPLi5vUH5HimzVov42dVYs32vb3Ox1qfLy8v2PTezTjP6v5UVZdKaudx5FEPgS4Bd1265sUP6u8n3tmnUXUIxYOyoiXgrMpLhV1uh7VtIIMhBSJ1tqSHRm9paBRtVw6V4W/xKHsj8MML+ifP/6ter2eajil+g9DdatW75vRdGfpcrqA2wbqnsr1ve3a80BttXqb/OO5atKbZubvS5V1irf5zWxTzOW5Wve75H6FTXfc+Nr1j0WES+n6MP0Jop+VQAPRMSPgK9WXCdJw8xbY9LA+gOmDSu2b1RXrv/zOuVttXqN6unf9/TMjAFek5pv/lI2qFjf365GAWI2WNdf7juDtPk1dfs0c12qPFK+T2hin2Ysy9e8aZnZk5nvBdYHXgh8mKLv0RfKl6RRYCAkDexm4N/A9hGxdoPt/b/o/1Gz7h8UP1u7NCg/ucG6myh+ub+8IkgYTrtUDMWfXL7/s8G2Ri6nuM3VqE9RlWavS5XLyve9h1i+j5pszBD0X4NdIqJR1rzR13yZZeH6zPwBxeg8gDcPR92SBmcgJA0gM58Bfklxi+fLtdsiYguKv+IXAKfUbPp5+f61iFilpvw6FP2J6o/RSzFabCPg+xGxan2Zck6bbZfvbIDi9tsH6urel6J/0Byg0fD5pWTmfRTXZYdyXqClAoYo5nCqzWI1dV0GcDLFtAfvj4hXNzjuxLpVD1IMfR+SsrP1nylGlB1ZV/fOFB3AHwZOb6LN9W18YURs1mBTf8bu38tat6Tm2EdIGtxRFJmPIyJiR4q5dPrnEVqDYjj67TXlTwX2p+j7cV1EnEHR7+htFPMQbdHgGF+hmAPncGCfiDiPog/M+hTBy38BnwNuWM5z+RPwPxGxN3A1i+cRehp4b0VH6ipHlG37MnBIRFxM0QdpY4pO0jsCBwL912ZZrstSMvOBiDiIYq6dv0bEH4FrKEaSvYgi6KkNwM4FDoiIs4CrKPqBXZiZFw5wmMMpOi8fFxF7UHSC7p9HqA9493KO7Hod8O2I+BtFRvA+ik76+5b1H7ccdUtqgoGQNIjMfKicgPEzFEHDxyhGcV0OHJeZ59SVz4h4O0UANYUiYJhPkRH5MkXQUX+MBRHxZorRSFMoJuZbHbifIpA4miIDs7z+XrbhK2W7gmIU1Ocy84pmKio7/O4KTKXIkuxHMYfQvRQj4D5KkVnpL9/0dRng2GdHMRv2p4HdKIa7P0wRVHy9rvhHKPo57UYxueE4ik7KlYFQZt5W1v/5cp/JFFmoPwFfa/ZaNTCLYmLOV1MEP8+muBZ/Br5dDtOXNAqiZloUSR2qfETFX4EvZeYxLW2MJLUR+whJkqSuZSAkSZK6loGQJEnqWvYRkiRJXcuMkCRJ6loGQpIkqWsZCEmSpK5lICRJkrqWgZAkSepaBkKSJKlr/X/v1CjIn6yTewAAAABJRU5ErkJggg==\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "testResults = test.classify(clf2).errorMatrix(label, 'classification').getInfo()\n", "testCM = pd.DataFrame(np.asarray(testResults), index=['not','built-up'], columns=['not','built-up'])\n", "\n", "fig, ax = plt.subplots(1, figsize=(10,10))\n", "sns.heatmap(testCM/testCM.sum(axis=1), annot=True)\n", "ax.set_xlabel('model predictions', fontsize=20)\n", "ax.set_ylabel('actual', fontsize=20)\n", "plt.title(\"Test data confusion matrix\", fontsize=20);\n", "acc = (testCM.loc['built-up','built-up'] + testCM.loc['not','not']) / testCM.sum().sum()\n", "mcc = get_mcc(testCM)\n", "print(f\"Our classifier has an accuracy of {acc:.5f} on the test data.\")\n", "print(f\"Our classifier has an MCC score of {mcc:.5f} on the test data.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "...it's about the same.\n", "\n", "If you spend a lot of time with machine learning algorithms, you will see that model tuning at best usually delivers incremental results. But thinking more about the data that is put into the model (feature engineering, additional sources, etc) is usually more likely to imrpove your model.\n", "\n", "What happens when we try cleaning the data like we did before by standardizing the VIIRS-DNB band?" ] }, { "cell_type": "code", "execution_count": 15, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Our classifier has an accuracy of 0.90261 on the test data.\n", "Our classifier has an MCC score of 0.62630 on the test data.\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAkIAAAJrCAYAAAAMBPz+AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAAA+uklEQVR4nO3deZxcZZm4/etOCDuyypKEJUAQEQWRRR0QFIGgICo4LA4KypsBRVxmVHBl0XHBcRwFxagI+EPiggwBgqAgq6IBDMiqAQJJCPuOLOn0/f5xTieVSp1eku6u6qrrm099quqc5zznOdXd6bvvZzmRmUiSJHWiUc1ugCRJUrMYCEmSpI5lICRJkjqWgZAkSepYBkKSJKljGQhJkqSOZSCkjhYRsyNidrPbMZJFxCsi4rvlZ9kVERkR2w/xOTcrz3PWUJ6n05Sf6VXNboc0nAyE1KfyP8eBPI4YgjYcMVR1D4YO/wXyTeBjwN+ArwEnAQ81tUUdKCL2KL8PT2x2W6SRZIVmN0AjwkkNtn0CWBP4X+Cpun0zh7Y5ajH7AX/PzP2H8ZzzgFcDTw/jOTvBq4F/NrsR0nAyEFKfMvPE+m1lZmZN4DuZOXuYm6TWMha4ZjhPmJkLgLuG85ydIDP9TNVx7BrToIuIXSLi1xHxUES8HBFzIuKHETG2QdnNI2JKRMyKiBci4omI+FtEnBER65ZlrgJ+Wh7y07puuM360Z6IiGMj4vaIeDEi5kXEaRGxZkX5NSPi0xFxZUTMLa/h0YiYFhFvrCt7RET03Kdm97q2nVhX7vyIuLe8zmci4vqI+Lf+fKYN2nhwRFxRfl4vluNzzouIHevKrRQRx0fErRHxz/K810bEvzaoc9G4m/L11Ih4rKz/xojYr678VeW1R921X1X72VR1ZzbqToyINSLiixFxW9nWZyPinoj4RUS8oVFbG9S7UUScXn4mPV+739QeX1N2URsj4q3lNT1bnvuSiHh19VdhqboWdU1FxI4R8duIeDoiniy/9huX5TYvP9tHy++FP0TEdg3q2yoivl5+9o9GxEsRcX/58zK+ruxZwB/Kt1+u+z7co8G1Tiqv9ema79+lviYRMSEiniq/zzatO+dqEXFnRCyMiN37+zlJrcaMkAZVRBwJ/Ah4CZgGzAEmAkcB+0fEGzPzgbLsRsAM4BXAdOB8YGVgAnA4cBrwOHAWRffbAcCFLNn19lQ/mvUd4DhgPjAFWFDWtQuwIvByXflXA1+lyHJcAjwJbAK8C9g3IvbPzN+WZWdSdB1+Gbi/bGuPq2pe/wC4o6xzPrAu8A7gZxHxqsz8Yj+ug4gIiqDwg8BjwG+AR4HxwFuBu4Eby7IrApcBu1NkT04HVgUOAn4REdtn5ucanGZT4C/AvcDPgHWAg4ELI+LtmdnzC/es8hrrr312f66l4tp+C7wZ+BPwY6AL2BjYA7gWuKmPOiYA11Fkqa4EziuPfx/wzog4MDMvbnDofhTfE5cCZwDbUHx9doqIbTLzsQFcyk7AZ4GrKX4WXgu8F3htRLyrbN9dwDkUn/V7gd9FxOaZ+VxNPe8FjqYIcP5I8X36Ghb/LO2YmfPKsv9XPn+wPO9VNfXMrmvfQcCkmmvdrOpCMvO+iDgK+BVwXkS8JTO7yt3fB7YGTszMq3v9RKRWlpk+fAz4QfGfawKb1WzbiuI/61nAuLrybwMWAhfUbPtYWcfHG9S/GrBKzfsjyrJHDLCdby6PmwWsU7N9ZYpftgnMrjtmTWC9BnWNBx4E7mywL4GremnHFg22rQhcQRGYjevn9Uwuz/UXYM26faOBjWren1CWnQ6sULN9/Zqv35trtm9Wbkvgy3V179NTV3+vva+vWf1xFAFD1n6P1OwbBazdoK1n1ZW7rNz++QbfB10UgfXqDdrYBexZd8zXyn2f6efXZo+az+/9dft+Um5/okHbvtjo5wAYB6zU4Dx7U/ws/aDi/CdWtK/nWruBSf35mtRs/36572vl+w+U7/8AjOrP5+PDR6s+7BrTYDoGGEPxH/q82h2ZeSVFhmj/iFij7rgX6ivKzOczc6nty+DI8vmrmflETf0vUgQKS8nMp7NBBiAz5wK/BraOiE0G0ojMvKfBtpcpsjQrAHv2s6qPlc//nplLDBTOzIWZOb9m04cofll9Khf/FU9mPgKcUr49qsE57ge+Ulf3ZcADwM79bOfyaPT90J2ZT/Z2UNldtDdFO79Zd/wfKbJD61BkWupNzcwr6rZNKZ8Hes3XZea5ddvOLp+fBr5et++c8nn72o2ZOS8zX6qvPDMvB26nCE6XxYW5OKPZX58CbgE+GxHHUgRGj1IEfN3L2A6pJdg1psH0pvJ594jYqcH+9SmyFltRdHFMA/4LOD0i9qH4a/564I7MzAbHL4sdyudGqftrKTIBS4mIfwE+TnFN61Nkb2qNo/iF2y9l4PRZioBnE2CVBvX1VcdqwLbAw5n51z7KrgFsCczLxgNgryyfX99g38zMXNhg+xwWf42Hwh0UXY2HluNRLqToRrqxDBr70nMt12YxmLrelcC/leXOqdt3Y4Pyc8rntftx7r7qerB8bvTZ9vzRUD/uJ4D3U2RytivbMbqmSH8+k0b+MtADMvPFiDiY4tq+RxFgH5SZD/Z+pNT6DIQ0mNYtnz/dR7nVATLz/ojYGTiRYsxCz1/qcyLiW5n53UFoU8+A6Ifrd2Tmwoh4vH57RLyHIvPzIvA74B7geYouhT0oxtys1N8GRMTmFL981qYIvi6nyAwspOji+WA/61urfJ7XW6FSz3XPr9jfs32tBvueqjimiyGcYFF+Pd4GfIliHMs3yl3PRsTZwAm55BiaeoN6zZnZVcQiSwQf/dFoSn9X1b6a84yp2/VtimUq5lP8kTCPxdmyIyjGFy2LZV3j6e/ArRTdjHdQfB9LI56BkAZTz3/ya2bmM/05IDPvBA6OiBUo/up9O0X3z/9GxPOZ+ZNBatMGFIN/F4mI0RTBW31gcQrFX9s7lu2rPeaHFIHQQHyqPM+RmXlWXX2HUgRC/fFU+dxn9ojF171hxf6N6soNhZ4uk6X+n4mItRodUHZ/fRL4ZERsSfFZ/ztwLEUAc3gv52uFax4UEbE+xQD/2yjGcT1bt//Q5ah+WbOtx1MEQY9RDNo+gWJSgTSiOUZIg+mG8nm3gR6YmV2ZeVNmfgPo+U/+3TVFeroTBvrX+c3lc6PgZTca/zGwJUX3XH0QNArYteI83b20bcvy+fwG+/odVGXm8xS/GDeIiEZdWrVln6XIZI2LiIkNiry1fL65wb7B0jOmZ+MG+3ZssG0JmTmrDIR3B56jmNXVm57uwl3LwLrecFzzYNmc4v/nyxsEQePL/fWW9WekTxHxZuBkilmJ25bPJ0VE1c+DNGIYCGkwnUYxA+p/ImKr+p0RsWJE7FbzfueI2KBBPT3bale47enCGtAgZRZP6f58RKxTc+6VKWYFNTIbmBg16x6V4zW+TDGtupHHafwLv6c+KLrVFinHRTUarNybnu7CH0bdOkgRMapckqDHmRRr/JxaZr96yq1HMVOpp8xQuZEiQDwsIlatOf861A1mLrdPiIjXNKhnbYquw14Hz5eD2X9H0d34ibq6dwEOowjOLhjIRTTJ7PJ517qv3eoUU/IbBXrL+jPSq4hYm2Kg+ULgkMx8mGI5hS6KKfXr9na81OrsGtOgycy7IuJDFL9cb4+I31KMKxhD8Z/zbhQzTbYuDzkM+GhEXE0xvf1JYAtgf4p1iL5TU/2fKAKjT5S/SHvG/HyvfvZUXZuuj4jvUXS33RYRv2bxOkJP0ng8yf9QrK/y14g4vyz/LxRB0EVl++pdARwSERdRDATvAq7JzGsoZtgcCfyqrG8exV/Vk4BfUvxS6a8fU2SlPgD8IyIupPhMx1IsUXAmxZgrgG8B+5bXektETKdYR+h9FAPAv5mZ1w3g3AOSmfMj4lyK7qyZEXEJxZpR76BYT6k+q7UdcEFE3ESR+XoQeGXZ/jEsHjPUm6MpBtyfGhF7UwRjPesIdVN0Tz7by/EtITMfioipwCEUn93lFGOg9qIYuzaTullmFFmaeRTfhy9TDOZP4GeZef9yNOdMip/f4zJzZtm+WyLiPyj++PkpxRpb0sjU7Pn7PkbmgwbrCNXsey1FJuZ+ioDmCYpfbD8E3lZTbheKhQZvKcu8QBEQ/RTYtkG9kygCoudYvF7LUudvcFxQjDG5s2zPgxTT1tcsr2N2g2OOoPhl8zzFmIgLyus6sTzvHnXl1wd+ThGgLaRuPReKsRVXUgRfz1LMhno3faz90ss1vZ9iJtzTFL8Y7wPOBXaoK7cy8Lny83+h5tyHNqhzMxqszVOz/6riv4yltleuoUSRyTkVmMviNaZOoPgjrH4dofEUswivpxjQ+1J53KXAvv1tK8UYqh+U338vl1+//wN2qvg693utoz6+JpVfy358tkudhyJo/Wr5mb1IMYvtdIrxZlVfi50ogvKnKQK/Rd+rfV1ro3aweK2vCyvK/6bc/8mBfP/68NFKj8gcrFnKkiRJI4tjhCRJUscyEJIkSR3LQEiSJHUsAyFJktSxDIQkSVLHGtHrCC147F6nvElNsPYmeza7CVLHeu6f98Vwnm84f9eOWW/zYb02MCMkSZI6mIGQJEnqWCO6a0ySJA2x7oV9lxnBzAhJkqSOZUZIkiRVy+5mt2BImRGSJEkdy4yQJEmq1m1GSJIkqS2ZEZIkSZXSMUKSJEntyYyQJEmq5hghSZKk9mRGSJIkVXOMkCRJUnsyEJIkSR3LrjFJklTNm65KkiS1JzNCkiSpmoOlJUmS2pMZIUmSVM0FFSVJktqTGSFJklTJm65KkiS1KTNCkiSpmmOEJEmS2pMZIUmSVM0xQpIkSe3JjJAkSarmvcYkSZLakxkhSZJUzTFCkiRJ7clASJIkdSy7xiRJUjUXVJQkSWpPZoQkSVI1B0tLkiS1JzNCkiSpmmOEJEmS2pMZIUmSVCnTW2xIkiS1JTNCkiSpmrPGJEmSmi8iJkXE3RExKyKOb7B/7Yi4ICJujYi/RMS2fdVpRkiSJFVrkVljETEaOB3YC5gLzIiIaZl5R02xzwEzM/M9EbF1WX7P3uo1IyRJkkaCnYFZmXlvZr4MTAUOqCuzDXAFQGbeBWwWERv0VqmBkCRJqpbdw/fo3ThgTs37ueW2WrcA7wWIiJ2BTYHxvVVqICRJklpCREyOiBtrHpNrdzc4JOvefx1YOyJmAh8D/gp09XZOxwhJkqRq3cO3jlBmTgGmVOyeC2xc83488GDd8c8ARwJERAD3lY9KZoQkSdJIMAOYGBETImJF4BBgWm2BiFir3AdwFHBNGRxVMiMkSZJaXmZ2RcSxwGXAaODMzLw9Io4u958BvBo4JyIWAncAH+6rXgMhSZJUrYUWVMzM6cD0um1n1Lz+EzBxIHXaNSZJkjqWGSFJklStRRZUHCpmhCRJUscyIyRJkqq10BihoWBGSJIkdSwzQpIkqZpjhCRJktqTGSFJklTNjJAkSVJ7MiMkSZIqZQ7fTVebwYyQJEnqWGaEJElSNccISZIktSczQpIkqZorS0uSJLUnAyFJktSx7BqTJEnVHCwtSZLUnswISZKkag6WliRJak9mhCRJUjXHCEmSJLUnM0KSJKmaY4QkSZLakxkhSZJUzTFCkiRJ7cmMkCRJqmZGSJIkqT2ZEZIkSdWcNSZJktSezAhJkqRqjhGSJElqTwZCkiSpY9k1JkmSqjlYWpIkqT2ZEZIkSdUcLC1JktSezAhJkqRqjhGSJElqT2aEJElSNccISZIktSczQpIkqZoZIUmSpPZkRkiSJFXLbHYLhpQZIUmS1LHMCEmSpGqOEZIkSWpPZoQkSVI1M0KSJEntyYyQJEmq5r3GJEmS2pOBkCRJ6lh2jUmSpGoOlpYkSWq+iJgUEXdHxKyIOL7B/jUj4qKIuCUibo+II/uq04yQJEmq1iK32IiI0cDpwF7AXGBGREzLzDtqin0UuCMz94+IVwJ3R8S5mflyVb1mhCRJ0kiwMzArM+8tA5upwAF1ZRJYIyICWB14AujqrVIzQpIkqVrrjBEaB8ypeT8X2KWuzGnANOBBYA3g4Mze5/+bEZIkSS0hIiZHxI01j8m1uxscUt9vtw8wExgLbA+cFhGv6O2cZoQkSVK1YcwIZeYUYErF7rnAxjXvx1NkfmodCXw9MxOYFRH3AVsDf6k6pxkhSZI0EswAJkbEhIhYETiEohus1gPAngARsQHwKuDe3io1IyRJkqq1yC02MrMrIo4FLgNGA2dm5u0RcXS5/wzgFOCsiPgbRVfaZzPzsd7qNRCSJEkjQmZOB6bXbTuj5vWDwN4DqdNASJIkVcru1lhHaKg4RkiSJHUsM0KSJKla66wjNCTMCEmSpI5lRkiSJFVrkVljQ8WMkCRJ6lgGQpIkqWPZNSZJkqo5fV6SJKk9mRGSJEnVnD4vSZLUnswISZKkamaEJEmS2pMZIUmSVC2dNSZJktSWzAhJkqRqjhGSJElqT2aEJElSNVeWlgbmuhtuZL9DjmLff/0QP/7ZL5fa//Qzz3LcCSfzng8cwyFHfZx/3DsbgPvun8uBH/zooscue72Xn/3igmFuvTRyvX2vt3DzzCu45W9/4FP/cfRS+7faanOu+MP5PP7kXRz38f9v0fZx4zZi+qU/56abf8eMGy/jIx85YhhbLTWXGSENqoULF/KV/z6dH33nv9hw/fU4+KiP89Zdd2GLCZsuKvOjc37B1hO34Ltf+xL33j+Hr/736fzku19nwqbjOf/s0xfV87Z3H86eu7+5WZcijSijRo3i2/9zMu/a73DmzXuIa669kOmX/J677pq1qMyTTz7Np//zJPbff+8lju1a2MUJJ3yVW2bezuqrr8a111/ElVdet8Sx6mDpGKEhFREr9WebRoa/3fl3Nhk/lo3HbcSYMWPYd8/dufLaG5Yoc8/sB3jjG7YDYPNNN2be/Id57Iknlyhzw40z2XjcRozdcINha7s0ku2443bce8/9zJ49hwULFvDrX1/EO/fba4kyjz76ODffdCsLFixYYvvDDz3KLTNvB+C5557n7rtnsdHYDYet7VIzNT0QAv7Uz20aAR559DE2XP+Vi95vsP56PPLo40uUedWWm/P7q/8IwN/uuJv5Dz/Cw488tkSZS6+4mne8ffehb7DUJsaO3ZC58+Yvej9v3kOMXYZgZpNNxrHddttw44yZg9g6jWjdOXyPJmhaIBQRG0bEG4BVIuL1EbFD+dgDWLVZ7dLyabTuVsSS7486/H088+xzHPjBj3Lur6ex9cQtGD169KL9CxYs4Krr/szeb9ttiFsrtY+o/0EDcoAL4a222qqce94P+OxnTuHZZ58brKZJLa2ZY4T2AY4AxgPfrtn+LPC5qoMiYjIwGeD7//0VjvrAoUPYRA3UBuuvx0OPPLro/cOPPMYr11t3iTKrr7YaX/n8p4DiP+p9DjqC8WMXd4Fde8ONvHqrLVhvnbWHp9FSG5g3bz7jx2206P24cRsyf/7D/T5+hRVW4Nyf/4BfTL2QaRdeNhRNlFpS0wKhzDwbODsiDszM8wdw3BRgCsCCx+5t7zl9I9C2W2/FA3MfZO6DD7HBK9fl0iuu5ptf/uwSZZ559jlWWXklxowZw/kX/ZY3bP9aVl9ttUX7p//uKt6x1x7D3HJpZLvpplvZYsvN2HTT8Tz44MMcdND+fOjIj/f7+O//4BvcffcsTvveT4awlRqJss0XVGyFWWNXRMS3gbeU768GTs7Mp5vYJi2jFVYYzec+eQz//qkvsHDhQt6z395sufmm/OKCSwA4+D3v5N775/C5U77F6FGj2HyzTTj5hE8sOv6FF1/kTzP+ypc/c1yTrkAamRYuXMh/fOrL/N+0cxg9ehQ/O+dX3HnnP/jwUYcB8JMf/5z1N1iPa6+bxhprrE53d/LRY49kxx32Ztttt+aw97+X2/52F3+8ofhZPfHLp3L5ZVc18Yqk4RED7UMe9AZEnA/cBpxdbjoc2C4z39vXsWaEpOZYe5M9m90EqWM998/7lh4QNoSe/+oHhu137WqfP2dYrw1aIyO0RWYeWPP+pIiY2azGSJKkztEKgdALEbFrZl4HEBH/ArzQ5DZJkiRo+wUVWyEQOoZi0PSa5fsngQ82sT2SJKlDtEIgdCfwTWALYC3gaeDdwK3Na5IkSQLa/qarrRAIXQg8BdwMzGtuUyRJUidphUBofGZOanYjJElSA22+jlAr3GvsjxHx2mY3QpIkdZ5WyAjtChwREfcBLwEBZGa+rrnNkiRJjhEaevs2uwGSJKkzNT0Qysz7m90GSZJUoc3XEWqFMUKSJElN0fSMkCRJamFtPkbIjJAkSepYBkKSJKlj2TUmSZIqpQsqSpIktSczQpIkqZqDpSVJktqTGSFJklTNjJAkSVJ7MiMkSZKqeYsNSZKk9mRGSJIkVXOMkCRJUnsyIyRJkiqlGSFJkqT2ZCAkSZKqdefwPfoQEZMi4u6ImBURxzfY/+mImFk+bouIhRGxTm91GghJkqSWFxGjgdOBfYFtgEMjYpvaMpl5amZun5nbAycAV2fmE73V6xghSZJUrXXuPr8zMCsz7wWIiKnAAcAdFeUPBc7rq1IzQpIkaSQYB8ypeT+33LaUiFgVmASc31elBkKSJKklRMTkiLix5jG5dneDQ6oGFu0PXN9XtxjYNSZJknozjNPnM3MKMKVi91xg45r344EHK8oeQj+6xcCMkCRJGhlmABMjYkJErEgR7EyrLxQRawK7Axf2p1IzQpIkqVqLLKiYmV0RcSxwGTAaODMzb4+Io8v9Z5RF3wNcnpnP96deAyFJkjQiZOZ0YHrdtjPq3p8FnNXfOg2EJElSpczWyAgNFccISZKkjmVGSJIkVWuRMUJDxYyQJEnqWGaEJElSNTNCkiRJ7cmMkCRJqpRmhCRJktqTGSFJklTNjJAkSVJ7MiMkSZKqdTe7AUPLjJAkSepYBkKSJKlj2TUmSZIqOX1ekiSpTZkRkiRJ1cwISZIktSczQpIkqZrT5yVJktqTGSFJklTJWWOSJEltyoyQJEmq5hghSZKk9mRGSJIkVXKMkCRJUpsyIyRJkqo5RkiSJKk9mRGSJEmV0oyQJElSezIQkiRJHcuuMUmSVM2uMUmSpPZkRkiSJFVysLQkSVKbMiMkSZKqmRGSJElqT2aEJElSJccISZIktSkzQpIkqZIZIUmSpDZlRkiSJFUyIyRJktSmzAhJkqRqGc1uwZAyIyRJkjqWGSFJklTJMUKSJEltykBIkiR1LLvGJElSpex2sLQkSVJbMiMkSZIqOVhakiSpTZkRkiRJldIFFSVJkpovIiZFxN0RMSsijq8os0dEzIyI2yPi6r7qNCMkSZIqtcoYoYgYDZwO7AXMBWZExLTMvKOmzFrA94FJmflARKzfV71mhCRJ0kiwMzArM+/NzJeBqcABdWUOA36TmQ8AZOYjfVVqRkiSJFVqoXWExgFzat7PBXapK7MVMCYirgLWAP43M8/prVIDIUmS1BIiYjIwuWbTlMyc0rO7wSFZ934F4A3AnsAqwJ8i4obM/HvVOQ2EJElSpawPNYb0XDkFmFKxey6wcc378cCDDco8lpnPA89HxDXAdkBlIOQYIUmSNBLMACZGxISIWBE4BJhWV+ZCYLeIWCEiVqXoOruzt0rNCEmSpEqtMkYoM7si4ljgMmA0cGZm3h4RR5f7z8jMOyPit8CtQDfw48y8rbd6DYQkSdKIkJnTgel1286oe38qcGp/6zQQkiRJlVolIzRUHCMkSZI6loGQJEnqWHaNSZKkSsM5fb4ZzAhJkqSOZUZIkiRVcrC0JElSmzIjJEmSKmWaEZIkSWpLZoQkSVKl7G52C4ZWr4FQRHSz9C3u+yMz0yBLkiS1tL6ClWtYtkBIkiS1ge42HyPUayCUmXsMUzskSZKGnd1XkiSpkrPGJEmS2tQyZYQiYiNgT2AcsFKDIpmZpyxPwyRJUvO1+8rSAw6EIuIk4Pi6Y4PFg6p7XhsISZKkljagrrGIeD/wReBa4CCKoOds4DDgR0A3MBV42+A2U5IkNUPm8D2aYaAZoWOAucCkzOyKCIDZmTkVmBoRFwCXAOcNbjMlSZIG30AHS78WmJ6ZXTXbRve8yMzLgMuATw9C2yRJkobUQDNCY4DHa96/AKxZV+Y24OjlaZQkSWoN7T5YeqAZofnARjXvHwBeV1dmHNCFJElSixtoIPRXiu6xHlcCu0XE4RGxWkS8EziwLCdJkka47oxhezTDQAOhi4HXRMSE8v3XgaeBs4BngGkUM8m+MFgNlCRJGioDGiOUmWdRBD097+dExE7AfwBbALOB72fm3waviZIkqVna/RYby32vscy8Dzh2ENoiSZI0rLzpqiRJqtSshQ6Hy4ACoYjYpL9lM/OBgTdHkiRp+Aw0IzSbxfcU600uQ92SJKnFNGs213AZaLByDo0DobWA7YFNgauA+5enUZIkScNhoLPGjqjaFxGjKG7IejTwweVrliRJagXtPmtsoOsIVcrM7sw8iaL77OuDVa8kSdJQGYpxPH8EPjAE9UqSpGHW7rPGBi0jVGMdYLUhqFeSJGlQDWpGKCLeDhxMcQd6SZI0wjlrrEZEXNlLPRsDPesMnbw8jZIkSRoOA80I7VGxPYEngcuAb2VmVcA0qFYZu9twnEZSnZnjX9/sJkgaJu0+a2yg0+eHYkyRJElSUxjYSJKkjjWgQCgiroyIXqfGR8S/9TKWSJIkjSDdGcP2aIaBZoT2ADbro8ymwO7L0hhJkqThNBQLKq4CdA1BvZIkaZi1+XqKyxQINfxMIiIops+/A5izPI2SJEkaDn0GQhHRzZLBz4kRcWJvhwD/tZztkiRJLcAFFeEaFgdCbwEeoLixar2FwOPAFcCPB6NxkiRJQ6nPQCgz9+h5XWaHfpqZrhwtSVIHcEHFJU0AnhqCdkiSJA27gQZCjwCvjIgXMvPl+p0RsRKwAfBIZr44GA2UJEnN093sBgyxga4j9CXgbmD1iv2rAXcBn1ueRkmSJA2HgQZC+wK/z8wnGu0st/8e2G95GyZJkpoviWF7NMNAA6HNgL/3Uebv9L36tCRJUtMNdIzQGPruLkxg5WVrjiRJaiXdbb609EAzQvfS933E9gDuX6bWSJIkDaOBBkLTgDdExGca7YyI44EdgP9bznZJkqQW0E0M26MvETEpIu6OiFllzFG/f4+IeDoiZpaPL/VV50C7xr4FvB/4WkT8K3A5MA8YB+wDbE+x8vQ3B1ivJElSpYgYDZwO7AXMBWZExLTMvKOu6LWZ2e9JWwMKhDLzyYjYAzgXeBNF9idhURj3R+DfMvPJgdQrSZLUh52BWZl5L0BETAUOAOoDoQEZ8N3nM3M28C8RsQPwRmAtitWmb8jMm5enMZIkqbU0a1p7A+OAOTXv5wK7NCj3poi4BXgQ+M/MvL23SgccCPUogx4DH0mSNCgiYjIwuWbTlMyc0rO7wSH1c9puBjbNzOci4h0UY5Yn9nbOZQqEImIjYE+K6GylRg3LzFOWpW5JktQ6hvMWG2XQM6Vi91xg45r34ymyPrXHP1PzenpEfD8i1svMx6rOOeBAKCJOAo6vOzZYHJX1vDYQkiRJg2UGMDEiJlBM1DoEOKy2QERsCDycmRkRO1PMjn+8t0oHNH0+It4PfBG4FjiIIug5u2zIjygCx6nA2wZSryRJak2tcouNzOwCjgUuA+4EfpmZt0fE0RFxdFnsIOC2cozQd4FDMrPXJSEHmhE6hiI1NSkzuyICYHZmTgWmRsQFwCXAeQOsV5IkqVeZOR2YXrftjJrXpwGnDaTOgS6o+FpgehmV9Rhd04DLKCK1Tw+wXkmS1IK6h/HRDAMNhMawZF/bC8CadWVuA7ZbnkZJkiQNh4F2jc0HNqp5/wDwuroy44AuJEnSiNesTM1wGWhG6K8U3WM9rgR2i4jDI2K1iHgncGBZTpIkqaUNNBC6GHhNOXUN4OvA08BZwDMUN2UN4AuD1UBJktQ8rTJrbKgM9F5jZ1EEPT3v50TETsB/AFsAs4HvZ+bfBq+JkiRJQ2OZb7HRIzPvo5jXL0mS2kx3y9xqbGgMtGtMkiSpbSx3RkiSJLWv7ta5+/yQMCMkSZI6loGQJEnqWHaNSZKkSr3esbQNmBGSJEkdy4yQJEmq5C02JEmS2pQZIUmSVKk7nD4vSZLUlswISZKkSs4akyRJalNmhCRJUiVnjUmSJLUpM0KSJKlSd3tPGjMjJEmSOpcZIUmSVKmb9k4JmRGSJEkdy4yQJEmq5DpCkiRJbcpASJIkdSy7xiRJUiWnz0uSJLUpM0KSJKmSt9iQJElqU2aEJElSJafPS5IktSkzQpIkqZKzxiRJktqUGSFJklTJWWOSJEltyoyQJEmqZEZIkiSpTZkRkiRJldJZY5IkSe3JjJAkSarkGCFJkqQ2ZSAkSZI6ll1jkiSpkl1jkiRJbcqMkCRJqpTNbsAQMyMkSZI6lhkhSZJUqdsFFSVJktqTGSFJklTJWWOSJEktICImRcTdETErIo7vpdxOEbEwIg7qq04zQpIkqVKrZIQiYjRwOrAXMBeYERHTMvOOBuW+AVzWn3rNCEmSpJFgZ2BWZt6bmS8DU4EDGpT7GHA+8Eh/KjUQkiRJlXIYH30YB8ypeT+33LZIRIwD3gOc0d/rMxCSJEktISImR8SNNY/JtbsbHFIfP30H+GxmLuzvOR0jJEmSKg3nOkKZOQWYUrF7LrBxzfvxwIN1ZXYEpkYEwHrAOyKiKzP/r+qcBkKSJGkkmAFMjIgJwDzgEOCw2gKZOaHndUScBVzcWxAEBkKSJKkXrTJrLDO7IuJYitlgo4EzM/P2iDi63N/vcUG1DIQkSdKIkJnTgel12xoGQJl5RH/qdLC0JEnqWGaEJElSpX5Max/RzAhJkqSOZUZIkiRV6m7znJAZIUmS1LHMCEmSpEqtMn1+qJgRkiRJHcuMkCRJqtTeI4TMCEmSpA5mRkiSJFVyjJAkSVKbMiMkSZIqdUezWzC0zAhJkqSOZUZIkiRVcmVpSZKkNmVGSJIkVWrvfJAZIUmS1MEMhCRJUseya0ySJFVyQUVJkqQ2ZUZIkiRVcvq8JElSmzIjJEmSKrV3PsiMkCRJ6mBmhCRJUiVnjUmSJLUpM0KSJKmSs8YkSZLalBkhSZJUqb3zQWaEJElSBzMjJEmSKjlrTJIkqU2ZEZIkSZWyzUcJmRGSJEkdy0BIkiR1LLvGJElSJQdLS5IktSkzQpIkqZK32JAkSWpTZoQkSVKl9s4HmRGSJEkdzIyQJEmq5BghSZKkNmVGSJIkVXIdIakf9tl7D26/7RruuuM6PvPpjzYs8z/fPpm77riOm2/6Ha/ffttF2z927IeZ+dcruGXmlRz3saMWbT/pxE9z802/48YZl3PpJT9no402GPLrkEay1d+yAxN/fwYTr5zCekcftNT+1XZ5La++5RdscfF32eLi7/LKjx2yaN+oNVZj49NPYOLvfsCWl/+AVV6/9XA2XWoaAyEtt1GjRvHd//0q++3/b7x2u7dy8MHv5tWvnrhEmX0nvY2JW05g62125ZhjPsvpp30NgNe85lV8+MOH8aY3v5Md3rAX73zH29lyywkAfOu/f8AOb9iLHXfam0um/54vfP6Tw35t0ogxahRjTzqG2Ud+mVn7fIQ199+dlbbceKliz8+4nXv2O4579juOR783ddH2jb40meeuvol/7HUM97zzY7w0a85wtl4tLIfxXzMYCGm57bzT67nnntncd98DLFiwgF/+8kLetf8+S5TZf/99+Nm5vwbgz3+5mTXXWpMNN1yfrbeeyJ//fDMvvPAiCxcu5Jprb+DdB0wC4Nlnn1t0/GqrrUpmew/Yk5bHKtttxUv3z2fBnIfJBV08ffE1rLHXG/t17KjVV2G1nV/Dk7+8HIBc0EX3s88PZXOlltH0QCgiNo+IiyLisYh4JCIujIjNm90u9d/YcRsyZ+6Di97PnTefsWM3XKLMuLEbMnfO4jLz5s5n3NgNuf32u9httzeyzjprs8oqK7PvpLcxfvzYReVOOfmz3HfPDA499D2ceNKpQ38x0gg1ZsN1WTD/0UXvu+Y/xpgN1l2q3Kqv35otLvkem555IitN3ASAFTfekK4nnmHcNz/BFhf9L2O/9jFilZWGre1qbd3D+GiGpgdCwM+BXwIbAmOBXwHnNbVFGpCIWGpbffamqsxdd83i1FNP57eXnsf0i8/lllvvYGHXwkVlvvilbzBhi50477wL+OhHjhz8xkvtrO7n8IXbZ/H33T7EPe/8GI+fczGb/PALxY4VRrPKa7bgiXOnc8/+H6f7ny/xyqPf14QGS8OvFQKhyMyfZWZX+fh/9LKQZURMjogbI+LG7m5Tt61g3tz5bFyTxRk/biPmz394iTJz581n/MaLy4wbvxEPlmV+etZUdt5lEm/d80CefPIp/jHrvqXOcd7UC3jPe94xRFcgjXwLHnqcMRu9ctH7FTZajwWPPLFEme7nXqD7ny8C8NxVNxIrjGb02q+ga/5jLHjoMV645e8APPPb61ll2y2Gr/FqaY4RGnp/iIjjI2KziNg0Ij4DXBIR60TEOvWFM3NKZu6YmTuOGrVaE5qrejNunMmWW05gs802ZsyYMfzrvx7ARRdfvkSZiy++nMPfX8xi2WXnHXjm6Wd46KFHAHjlK4v0/cYbj+Xd796Xqb/4P4BFg6YB9t9vb+6++55huBppZHrh1r+z0mZjGTN+A2LMCqy531t49vd/XqLMCuuttej1Kq/bCkYFC598hq7HnmLB/MdYccI4AFZ/83a8+I8HhrP5UtO0wjpCB5fP/163/UMUmSHHC7W4hQsX8vFPfIHpl/yc0aNGcdbZv+COO/7O5P/vcACm/OhnTL/0CiZNeht333k9/3zhBY466lOLjv/VL37EOuuuzYIFXRx33Od56qmnAfivr57AVlttQXd3Nw88MI+PfPT4plyfNCIs7ObBE89gs7NPJkaN4slf/Y6X/vEAax+2LwBP/vxSXrHvrqzz/n3Jhd3kiy8x57hvLjp8/olnsPF3/pMYswIvP/AQcz/znSZdiDS8YiTPxFlhxXEjt/HSCDZz/Oub3QSpY21778VLD7ocQh/c7MBh+1179uzzh/XaoAUyQhHxgUbbM/Oc4W6LJElqXRExCfhfYDTw48z8et3+A4BTKCahdQGfyMzrequz6YEQsFPN65WBPYGbAQMhSZKarLtFeo4iYjRwOrAXMBeYERHTMvOOmmJXANMyMyPidRSz0ntdJr3pgVBmfqz2fUSsCfysSc2RJEmtaWdgVmbeCxARU4EDgEWBUGY+V1N+NXqZhd6jFWaN1fsnMLHPUpIkacjlMD76MA6ovffL3HLbEiLiPRFxF3AJxcSrXjU9IxQRF7H4+kcDr6ZIZUmSpA4SEZOByTWbpmTmlJ7dDQ5ZKn7KzAuACyLiLRTjhd7e2zmbHggB36p53QXcn5lzm9UYSZK0WPcwLnRYBj1TKnbPBWrvJDweeLCiLJl5TURsERHrZeZjVeWa3jWWmVf3PIC1DYIkSVIDM4CJETEhIlYEDgGm1RaIiC2jvKdTROwArAg83lulrZARqnUycHGzGyFJkgrNuvVFvczsiohjgcsohtKcmZm3R8TR5f4zgAOBD0TEAuAF4ODsY8HEVguEhn0hJUmSNDJk5nRget22M2pefwP4xkDqbHogFBErZeZL5dt/b7BNkiQ1SXezGzDEmj5GCPhTz4vM/Ev9NkmSpKHStIxQRGxIMf9/lYh4PYu7xV4BrNqsdkmSpMWGc9ZYMzSza2wf4AiK6W/frtn+LPC5ZjRIkiR1lqYFQpl5NnB2RByYmec3qx2SJKlaq8waGyrN7Br7VKPXPTLz2/XbJEmSBlMzu8bWaOK5JUmSmto1dlKzzi1Jkvqn3afPN7Nr7DOZ+c2I+B6Nb5p2XBOaJUmSOkgzu8buLJ9vbGIbJElSL/q4Q8WI18yusYvK57Ob1QZJktTZWuEWG3+gcdfY25rQHEmSVMMFFYfef9a8XpnizrFdTWqLJEnqIE0PhDLzprpN10fE1U1pjCRJWoKzxoZYRKxT83YUsCOwYZOaI0mSOkjTAyHgJhaPEeoCZgMfblprJEnSIt5iY+htA3wE2JUiILoWp9RLkqRh0AqB0NnAM8B3y/eHAj8D3te0FkmSJMBZY8PhVZm5Xc37P0TELU1rjSRJ6hitEAj9NSLemJk3AETELsD1TW6TJEnClaWHTET8jWJM0BjgAxHxQPl+U+COZrVLkiR1jmZmhPZr4rklSVI/uI7QEMnM+5t1bkmSJGiNMUKSJKlFtfs6QqOa3QBJkqRmMRCSJEkdy64xSZJUqd0XVDQjJEmSOpYZIUmSVKndF1Q0IyRJkjqWGSFJklTJMUKSJEltyoyQJEmq5IKKkiRJbcqMkCRJqtTtrDFJkqT2ZEZIkiRVau98kBkhSZLUwcwISZKkSq4jJEmS1KbMCEmSpEpmhCRJktqUgZAkSepYdo1JkqRK6YKKkiRJ7cmMkCRJquRgaUmSpDZlRkiSJFVKM0KSJEntyYyQJEmq5KwxSZKkNmVGSJIkVXLWmCRJUpsyEJIkSZUyc9gefYmISRFxd0TMiojjG+x/f0TcWj7+GBHb9VWngZAkSWp5ETEaOB3YF9gGODQitqkrdh+we2a+DjgFmNJXvY4RkiRJlVpojNDOwKzMvBcgIqYCBwB39BTIzD/WlL8BGN9XpWaEJEnSSDAOmFPzfm65rcqHgUv7qtSMkCRJqjScK0tHxGRgcs2mKZnZ070VDQ5p2LiIeCtFILRrX+c0EJIkSS2hDHqqxvXMBTaueT8eeLC+UES8DvgxsG9mPt7XOe0akyRJI8EMYGJETIiIFYFDgGm1BSJiE+A3wOGZ+ff+VGpGSJIkVepukVtsZGZXRBwLXAaMBs7MzNsj4uhy/xnAl4B1ge9HBEBXZu7YW70GQpIkaUTIzOnA9LptZ9S8Pgo4aiB1GghJkqRKwzlYuhkcIyRJkjqWGSFJklSpVcYIDRUzQpIkqWOZEZIkSZUcIyRJktSmzAhJkqRKjhGSJElqU2aEJElSJccISZIktSkzQpIkqZJjhCRJktqUGSFJklTJMUKSJEltykBIkiR1LLvGJElSpczuZjdhSJkRkiRJHcuMkCRJqtTtYGlJkqT2ZEZIkiRVShdUlCRJak9mhCRJUiXHCEmSJLUpM0KSJKmSY4QkSZLalBkhSZJUqduMkCRJUnsyIyRJkiqls8YkSZLakxkhSZJUyVljkiRJbcpASJIkdSy7xiRJUiVvsSFJktSmzAhJkqRKDpaWJElqU2aEJElSJW+xIUmS1KbMCEmSpEqOEZIkSWpTZoQkSVIl1xGSJElqU2aEJElSJccISZIktSkzQpIkqZLrCEmSJLUpM0KSJKlSOmtMkiSpPRkISZKkjmXXmCRJquRgaUmSpDZlRkiSJFVyQUVJkqQWEBGTIuLuiJgVEcc32L91RPwpIl6KiP/sT51mhCRJUqVWmT4fEaOB04G9gLnAjIiYlpl31BR7AjgOeHd/6zUjJEmSRoKdgVmZeW9mvgxMBQ6oLZCZj2TmDGBBfys1IyRJkiq10BihccCcmvdzgV2Wt1IzQpIkqSVExOSIuLHmMbl2d4NDljtKMyMkSZIqDWdGKDOnAFMqds8FNq55Px54cHnPaUZIkiSNBDOAiRExISJWBA4Bpi1vpWaEJElSpVYZIZSZXRFxLHAZMBo4MzNvj4ijy/1nRMSGwI3AK4DuiPgEsE1mPlNVr4GQJEkaETJzOjC9btsZNa8fougy67doodHg6jARMbnsD5Y0jPzZkxZzjJCaaXLfRSQNAX/2pJKBkCRJ6lgGQpIkqWMZCKmZHKMgNYc/e1LJwdKSJKljmRGSJEkdy0BILSUijoiIsc1uh9SKImKziLhtAOXfFRHHl69PjIj/LF/7cyaVDITUao4A/A9aGgSZOS0zv95g1xH4cyYBBkIaYuVfsHdGxI8i4vaIuDwiVomI7SPihoi4NSIuiIi1I+IgYEfg3IiYGRGrNLv9UgtaISLOLn92fh0Rq0bE7IhYDyAidoyIq8rXR0TEabUH9/VzFhF7RMTFNe9Pi4gjytezI+IbEfGX8rHl0F6qNPQMhDQcJgKnZ+ZrgKeAA4FzgM9m5uuAvwFfzsxfU9wj5v2ZuX1mvtCsBkst7FXAlPJn5xngIwM5eBB+zp7JzJ2B04DvDPBYqeUYCGk43JeZM8vXNwFbAGtl5tXltrOBtzSjYdIINCczry9f/z9g12E+/3k1z28a5nNLg85ASMPhpZrXC4G1mtQOqR3Ur3mSQBeL/z9feSCVRcQuZRfZzIh4V11djerLitfSiGQgpGZ4GngyInYr3x8O9GSHngXWaEqrpJFhk4joycQcClwHzAbeUG47sB91LPo5y8w/l11k22fmNOB+YJuIWCki1gT2rDv24JrnPy37ZUitYYVmN0Ad64PAGRGxKnAvcGS5/axy+wvAmxwnJC3lTuCDEfFD4B/AD4C/AD+JiM8Bf+5HHWdR8XOWmXMi4pfArWX9f607dqWI+DPFH9KHLu/FSM3mytKSpH6JiNnAjpn5WLPbIg0Wu8YkSVLHMiMkSZI6lhkhSZLUsQyEJElSxzIQkiRJHctASBpmEZE994Jajjo2K+s5a3Ba1XyNPpfyjukZEXsM0Tnb7nOUNDAGQpLa2mAEnpLalwsqSmplpwFTgQeGqP55wKspVjuX1IEMhCS1rHLhviFbvC8zFwB3DVX9klqfXWNqO7XjPiJii4j4dUQ8HhHPRsTlEbFtWe6VETElIuZHxIsRMSMi3lpR55oR8bWIuLss+2REXBYRb68ov2JEfDEi7omIlyLivoj4SkSs1Eu7V4iIj0TEDRHxTET8MyL+GhHHRsRy/axGxB7lZ3JiRLwpIn4fEU+Xn8llEbFjg2MWjc+JiMMi4s8R8Vy5unBPmVUj4oTyhp3Pl/v/FBENb70w0M+ltzFCEbF1RJwZEbPLuh6JiGsj4phy/xER0bNQ2u5lPT2PE8sylWOEImKjiDi9rP/liHg0In4TEW9oUPaIsp4jIuKtEXFV+dk+ExGXRMSrGxyzQUR8q/yeej4inipfnxURmzf6PCQNPjNCamebUdx36U6KeyttBrwHuCqKm1b+FngG+AWwDnAIcGlEbJWZi7piImIt4HpgG2AG8B1gPeBfgcsj4pjM/GFN+QB+CRwA3EPRvbMi8CHgtY0aGhFjgIuAfYC7gZ8DLwJvBb4H7EJxc9rltQtwAvB74HRgS+C9wFsiYu/MvLbBMf8B7FW27w/AmmWb1wKuBF4P3AycSfHH1T7AzyPiNZn5hZprHPDnUiUi3gn8CliJ4ut4HrAWsB3wGYr7b80ETgK+THEj0bNqqriqj/onUNzMdGx5jecBGwPvA94ZEQdm5sUNDt2vvL5LgTMovmfeAewUEdv03JoiinvsXQ9sAfyO4rMNYNPy+F9T3INP0lDLTB8+2upBEfBk+fh83b4vltufoPhFNapm3+Hlvv+pO+aH5fYfUq7GXm6fSDG25CVgs5rth5Xl/wSsXLN9HYoAIIGr6s5xYrn9e8Domu2jgZ+U+w5ocI1n9fMz2aPmMzm2bt8B5fZ/1H0ePW16Hnh9gzrPKvd/pm77yhTBSTew/SB9LnvUbFuv/NxfBnZv0K7xde+XqrevzxG4rOL7581AF/A4sHrN9iPK8l3AnnXHfK3+cwL2b/S9Vu5bEVij2T9HPnx0ysOuMbWz2cDX67adXT6vBHw6M7tr9v2c4hfZ9j0bykzNvwHPASdk5qJ70mTmP4DvUvzi+kBNPUeWz5/LzBdryj8BnFLfyLLb61jgIeCTmbmw5piFFBmZBN7f1wX3wyzg+7UbMvNC4GqK7NBuDY6ZkplL3IE8Ital+FxuzMxv1tX3IvBZigzHYTW7BvS59OKDwCuAH2Tm1fU7M3PuAOpaSkSMB/amGKBdf21/pMgOrUORSas3NTOvqNs2pXzeuUH5F+o3ZObLmfnsQNstadnYNaZ2NrM2qCg9WD7/vf6XTWYujIiHgfE1m7cGVgWuL39h17sS+AJF91CPHSiyIdc1KH9Vg21bAetSZGS+UPQgLeUFitlNy+vauuCvtl27U1xHfXDxlwbld6LIVi0ab1NnTPlc2+aBfi5V3lg+XzqAYwai52t5bRaDqetdSREEvh44p27fjQ3Kzymf167ZdjXFjLXjI2IHYDpFV1mj71lJQ8hASO1sqSnRmdlVBhpV06W7WPxLHMrxMMD8ivI929eqO+aJil+iDzXYtm75PJFiPEuV1XvZ118PV2zvadeaveyr1dPmncpHldo2D/RzqbJW+TxvAMcMxLJ8zXs8Vb+h5ntudM22ZyLijRRjmN5FMa4K4LGI+D7wlYrPSdIgs2tM6l1PwLRhxf6N6sr1vF6n7Far16ienmMvyMzo5TFh4M1fygYV23va1ShAzAbbesr9Tx9tfmvdMQP5XKo8VT6PG8AxA7EsX/MBy8y5mflhYH1gW+A4irFHXyofkoaBgZDUu7uBfwLbR8TaDfb3/KK/uWbbzRQ/W7s2KL9Hg213Ufxyf2NFkDCYdq2Yir9H+fzXBvsa+QtFN1ejMUVVBvq5VLmhfN63n+W7qcnG9EPPZ7BrRDTKmjf6mi+zLNyemd+jmJ0H8O7BqFtS3wyEpF5k5svAuRRdPCfX7ouILSj+il8A/Kxm10/L569GxMo15dehGE9Uf44uitliGwHfjYhV6suUa9pss3xXAxTdbx+pq/sAivFBs4BG0+eXkpmPUHwuO5brAi0VMESxhlNtFmtAn0svzqZY9uCYiHhLg/OOr9v0OMXU934pB1v/jmJG2Sfq6t6FYgD4k8AFA2hzfRu3jYjNGuzqydj9c1nrljQwjhGS+nY8Rebj2IjYiWItnZ51hNagmI5+X03584CDKcZ+3BYRF1KMOzqIYh2iLRqc4xSKNXCOBvaPiCspxsCsTxG8/AvweeCO5byW3wL/HRH7AreweB2hF4EPVwykrnJs2baTgcMj4jqKMUhjKQZJ7wQcCvR8NsvyuSwlMx+LiMMo1tr5Q0RcCtxKMZPsdRRBT20AdgVwSERcBNxEMQ7smsy8ppfTHE0xePnUiNibYhB0zzpC3cCRyzmz6+3AtyPijxQZwUcoBukfUNZ/6nLULWkADISkPmTmE+UCjCdQBA2fopjF9Rfg1My8vK58RsT7KAKoIygChvkUGZGTKYKO+nMsiIh3U8xGOoJiYb7VgUcpAokvUmRgltefyzacUrYrKGZBfT4zZwykonLA7+7AZIosyYEUawg9TDED7pMUmZWe8gP+XHo59yVRrIb9WWBPiunuT1IEFV+rK/5xinFOe1IsbjiKYpByZSCUmfeW9X+hPGYPiizUb4GvDvSzauAyioU530IR/LyC4rP4HfDtcpq+pGEQNcuiSGpT5S0q/gCclJknNrUxktRCHCMkSZI6loGQJEnqWAZCkiSpYzlGSJIkdSwzQpIkqWMZCEmSpI5lICRJkjqWgZAkSepYBkKSJKljGQhJkqSO9f8DVxE7Qbkvb5YAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "vmu = viirs.reduceRegion(reducer=ee.Reducer.mean(),scale=scaleFactor)\n", "vstd = viirs.reduceRegion(reducer=ee.Reducer.stdDev(),scale=scaleFactor)\n", "\n", "vmu = ee.Number(vmu.get('avg_rad'))\n", "vstd = ee.Number(vstd.get('avg_rad'))\n", "\n", "viirsclean = viirs.subtract(vmu).divide(vstd)\n", "\n", "fusedclean = se2.addBands(viirsclean)\n", "\n", "dataclean = fusedclean.select(trainingbands).sampleRegions(collection=points,\n", " properties=[label],\n", " scale=scaleFactor)\n", "\n", "# we'll create a column of random numbers\n", "dataclean = dataclean.randomColumn(seed=0)\n", "split_thresh = 0.8\n", "\n", "trainclean = dataclean.filter(ee.Filter.lt('random',split_thresh))\n", "testclean = dataclean.filter(ee.Filter.gte('random',split_thresh))\n", "\n", "clf3 = ee.Classifier.smileRandomForest(**new_params).train(trainclean, label, trainingbands)\n", "\n", "testResults = testclean.classify(clf3).errorMatrix(label, 'classification').getInfo()\n", "testCM = pd.DataFrame(np.asarray(testResults), index=['not','built-up'], columns=['not','built-up'])\n", "\n", "fig, ax = plt.subplots(1, figsize=(10,10))\n", "sns.heatmap(testCM/testCM.sum(axis=1), annot=True)\n", "ax.set_xlabel('model predictions', fontsize=20)\n", "ax.set_ylabel('actual', fontsize=20)\n", "plt.title(\"Test data confusion matrix\", fontsize=20);\n", "acc = (testCM.loc['built-up','built-up'] + testCM.loc['not','not']) / testCM.sum().sum()\n", "mcc = get_mcc(testCM)\n", "print(f\"Our classifier has an accuracy of {acc:.5f} on the test data.\")\n", "print(f\"Our classifier has an MCC score of {mcc:.5f} on the test data.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Improvement, albeit barely perceptible. Probably not worth adding this cleaning step to our final model pipeline, but worth exploring further -- there are many ways to improve signal to noise ratios in images.\n", "\n", "In a real scenario, you would spend a lot of time in this phase of research trying to validate your model. You should also look to independent sources of economic growth not just for additional data inputs but as validators.\n", "\n", "We'll move on for now, but think about some ways you might improve this model using the data we have (or exploring other data!).\n", "\n", "## Visualize prediction\n", "\n", "Applying that last trained classifier (but on non cleaned training data) on our entire image, we can visualize the predicted land cover." ] }, { "cell_type": "code", "execution_count": 16, "metadata": {}, "outputs": [], "source": [ "result = fused.select(trainingbands).classify(clf2)" ] }, { "cell_type": "code", "execution_count": 17, "metadata": {}, "outputs": [ { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "5826294b2d87459b954706575bca508a", "version_major": 2, "version_minor": 0 }, "text/plain": [ "Map(center=[27.87388743003947, 85.41973735675019], controls=(WidgetControl(options=['position', 'transparent_b…" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "ghslVis= {\"palette\":['000000', 'ffffff']}\n", "se2Vis = {\"min\":0.0, \"max\":0.3,\"bands\": ['B4','B3','B2']}\n", "\n", "# initialize our map\n", "map1 = geemap.Map()\n", "map1.centerObject(roi, 9)\n", "map1.addLayer(se2, se2Vis, \"S2\")\n", "map1.addLayer(viirs, {}, \"VIIRS-DNB\")\n", "map1.addLayer(ghsl, ghslVis, \"GHSL\")\n", "map1.addLayer(result.randomVisualizer(), {}, 'classified')\n", "map1.addLayerControl()\n", "map1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We have some odd truncated edges in our prediction set likely due to the re-sampling process (and maybe errors in preprocessing). Worth investigating further, but at at glance we can see the built up area (green) that appears quite close to what we would expect given the GHSL data. \n", "\n", "Our MCC score on the validation set of ~ 0.64 is far from perfect, but it'll do for our exploratory and learning purposes we can take a look at what we were trying to assess: change in land cover.\n", "\n", "Note that if you wanted to err on the side of caution for false positives (showing built-up land, ie growth, only when it is quite certain that it is there) there are metrics and weights that can \"penalize\" your classifier in training to decrease the false positives. The trade-off may be in higher false negatives of course, but this making the right choice for the task at hand is why {doc}`mod6_1_framing_the_analysis` is so critical!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "hide_input": false, "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.9" }, "toc": { "base_numbering": 1, "nav_menu": {}, "number_sections": false, "sideBar": true, "skip_h1_title": false, "title_cell": "Table of Contents", "title_sidebar": "Contents", "toc_cell": false, "toc_position": {}, "toc_section_display": true, "toc_window_display": false } }, "nbformat": 4, "nbformat_minor": 4 }