Lab 4 Supervised regression_2021.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Lab 4: Supervised regression"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Purpose:**\n",
    "The purpose of this lab is introduce you to concepts of supervised classification and regression: prediction of nominal or numeric values of a geographic variable from other geographic variables.  You will explore processes of training data collection, classifier selection, classifier training, image classification and accuracy assessment.  At the completion of the lab, you will be able to perform supervised classification and regression in Earth Engine.\n",
    "\n",
    "**Prerequisites:** Lab 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 1. Introduction to classification and regression\n",
    "For present purposes, define prediction as guessing the value of some geographic variable of interest *g*, using a function *G* that takes as input a pixel vector **p**:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\\begin{equation}\n",
    "G_{T}(p_i) = g_i \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *i* in this equation refers to a particular instance from a set of pixels.  Think of *G* as a guessing function and *gi* as the guess for pixel *i*.   The **T** in the subscript of *G* refers to a *training set* (a set of known values for p and the correct g), used to infer the structure of G.  You have to choose a suitable *G* to train with **T**.  When *g* is nominal (e.g. {'water', 'vegetation', 'bare'}), call this setup classification.  When g is numeric, call this setup regression.  This is an incredibly simplistic description of a problem addressed in a broad range of fields including mathematics, statistics, data mining, machine learning, etc.  Interested readers may see [Witten et al. (2011)](http://www.cs.waikato.ac.nz/ml/weka/book.html), [Hastie et al. (2009)](http://statweb.stanford.edu/~tibs/ElemStatLearn/) or [Goodfellow et al (2016)](http://www.deeplearningbook.org/)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###  2. Regression\n",
    "In the present context, regression means predicting a numeric variable instead of a class label.  No lab on regression would be complete without the requisite introduction to least squares regression.\n",
    "\n",
    "#### a. Ordinary Least Squares (OLS)\n",
    "\n",
    "A very [ordinary regression](https://en.wikipedia.org/wiki/Ordinary_least_squares) is when G is a linear function of the form G(**p**) = **βp** where **β** is a vector of coefficients.  Once G is trained by some training set **T**, guess the value for some unknown **p** by multiplying it with **β**.  Suppose the goal is to estimate percent tree cover in each Landsat pixel.\n",
    "\n",
    "- Import data to use as known values for g.  Search 'vegetation continuous fields' and import 'MOD44B'.  Get the most recent image out of this collection:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[20, 0], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(Togg…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import ee\n",
    "ee.Initialize()\n",
    "\n",
    "# importing ipygee for dynamic mapping\n",
    "import geemap\n",
    "\n",
    "Map = geemap.Map() # from ipygee\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "MOD44B EarthEngine dataset link: [https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD44B](https://developers.google.com/earth-engine/datasets/catalog/MODIS_006_MOD44B)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Julian date:  65\n",
      "Year:  2020\n"
     ]
    }
   ],
   "source": [
    "# The MOD44B Version 6 Vegetation Continuous Fields (VCF) yearly product is a global representation of surface vegetation \n",
    "# cover as gradations of three ground cover components: percent tree cover, percent non-tree cover, \n",
    "# and percent non-vegetated (bare).\n",
    "# https://lpdaac.usgs.gov/products/mod44bv006/\n",
    "    \n",
    "mod44b =ee.ImageCollection(\"MODIS/006/MOD44B\")\n",
    "\n",
    "# this data is the referential information, or \"true data\", the \"output data\", the \"y\" that we use to calibrate a model.\n",
    "\n",
    "tree = ee.Image(mod44b.sort('system:time_start', False).first()) #this command is asking for the first image (in 2020)\n",
    "\n",
    "# Get the timestamp and convert it to a date.\n",
    "julian_date = ee.Date(tree.get('system:time_start')).format('DD').getInfo();\n",
    "year = ee.Date(tree.get('system:time_start')).format('Y').getInfo();\n",
    "\n",
    "print('Julian date: ', julian_date); # ee.Date\n",
    "print('Year: ', year)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Check which calendar date it is here: [https://www.scp.byu.edu/docs/doychart.html](https://www.scp.byu.edu/docs/doychart.html)"
   ]
  },
  {
   "cell_type": "raw",
   "metadata": {},
   "source": [
    "I found that The Julian date of 65  is May 5th, 2020"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since water is coded as 200 in MOD44B, replace the 200's with 0's and display the result:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[20, 0], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(Togg…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "percentTree = tree.select('Percent_Tree_Cover').where(tree.select('Percent_Tree_Cover').eq(200), 0)# replacing 200 with 0\n",
    "\n",
    "Map.addLayer(percentTree, {'min': 0,'max': 100, 'palette':['red','yellow', 'green', 'blue', 'purple']}, \n",
    "             name='percent tree cover')\n",
    "\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We could benefit from a stronger or more detailed inland water occurrence such as *JRC/GSW1_0/GlobalSurfaceWater* (which is Landsat resolution). Let's create a water mask for it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[20, 0], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(Togg…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# we are adding this filter to minimize water pixels occurrence at Landast scale, which may not be vissible at MODIS scale\n",
    "waterOcc = ee.Image(\"JRC/GSW1_0/GlobalSurfaceWater\").select('occurrence')\n",
    "jrc_data0 = ee.Image(\"JRC/GSW1_0/Metadata\").select('total_obs').lte(0)\n",
    "waterOccFilled = waterOcc.unmask(0).max(jrc_data0)\n",
    "waterMask = waterOccFilled.lt(50)\n",
    "\n",
    "Map.addLayer(waterMask,[],'WaterMask')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let's update the *percentTree* layer with the *waterMask* and plot it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[20, 0], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(Togg…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "percentTreeWM = percentTree.updateMask(waterMask)\n",
    "\n",
    "Map.addLayer(percentTreeWM, {'min': 0,'max': 100, 'palette':['red','yellow', 'green', 'blue', 'purple']},\n",
    "             name='percent tree cover with WM')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that this image represents percent tree cover at 250 meter resolution in 2000.\n",
    "\n",
    "- Import data to use as predictor variables (**p**).  Search 'landsat 5 raw' and import 'USGS Landsat 5 TM Collection 1 Tier 1 Raw Scenes'.  Name the import l5raw.  Filter by time and the [WRS-2](https://drive.google.com/file/d/1Uv5H3cbJXbrLN0J0RSPLa18kz8ut0H_z/view?usp=sharing) (Open this file in Google Earth Pro) Path and Row to get only scenes over the San Francisco bay area in 2000:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "l8 = ee.ImageCollection(\"LANDSAT/LC08/C02/T1_L2\")\n",
    "# looking for Landsat images around the month of the MODIS dataset (forestry does not change much in few months)\n",
    "\n",
    "# For San Francisco: from Google Earth Pro and WRS2 file\n",
    "Path =44\n",
    "Row = 34\n",
    "\n",
    "l8filtered = l8.filterDate('2020-04-01', '2020-08-30').\\\n",
    "filterMetadata('WRS_PATH', 'equals', Path).\\\n",
    "filterMetadata('WRS_ROW', 'equals', Row)\n",
    "\n",
    "# while we are using a large large of months to merge them and find non cloudy pixels, the correct technique would be to \n",
    "# remove the clouds in a single image. The technique is presented in Lab 6."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use an median to get a cloud-free composite of Landsat imagery in 2020. Also apply the *waterMask* mask to the *landsat* composite."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[20, 0], controls=(WidgetControl(options=['position', 'transparent_bg'], widget=HBox(children=(Togg…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "landsat = l8filtered.median()\n",
    "\n",
    "landsat = landsat.select(['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7'])\n",
    "\n",
    "\n",
    "landsat = landsat.multiply(0.0000275).add(-0.2)\n",
    " \n",
    "landsat = landsat.updateMask(waterMask)\n",
    "# print(landsat.getInfo())\n",
    "Map.addLayer(landsat, {'bands': ['SR_B5', 'SR_B4', 'SR_B3'],'min': 0.0,\n",
    "  'max': 0.3}, 'False color composite')\n",
    "\n",
    "Map.centerObject(l8filtered.first().geometry()) #centering the map to landsat scene\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Specify the bands of the Landsat composite to be used as predictors (i.e. the elements of p):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "predictionBands = ['SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now that all the input data is ready, build T.  It's customary to include a constant term in linear regression to make it the best linear unbiased estimator.  Stack a constant, the predictor variables and the image representing known g:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "trainingImage = landsat.select(predictionBands).addBands(percentTreeWM)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Sample this stack at 1000 locations to get T as a table:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'type': 'Feature', 'geometry': {'type': 'Point', 'coordinates': [-121.20608001850448, 37.08640098326397]}, 'id': '0', 'properties': {'Percent_Tree_Cover': 4, 'SR_B1': 0.03611500000000001, 'SR_B2': 0.057166249999999974, 'SR_B3': 0.0921325, 'SR_B4': 0.13462000000000002, 'SR_B5': 0.22901375000000002, 'SR_B6': 0.32127625, 'SR_B7': 0.21655625}}\n"
     ]
    }
   ],
   "source": [
    "region = l8filtered.first().geometry()\n",
    "\n",
    "# // Create 5000 random points in the region.\n",
    "# randomPoints = ee.FeatureCollection.randomPoints(region);\n",
    "\n",
    "trainingTesting = trainingImage.sample(\n",
    "  region = region, \n",
    "  scale = 30, \n",
    "  numPixels= 1000,\n",
    "  seed=0,\n",
    "  geometries= True,  # Set this to False to ignore geometries\n",
    ")\n",
    "\n",
    "\n",
    "print(trainingTesting.first().getInfo())\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "# // Randomly split the data into 60% for training, and 40% for testing\n",
    "trainingTesting = trainingTesting.randomColumn('random',0);\n",
    "\n",
    "training = trainingTesting.filter(ee.Filter.lt('random', 0.6));\n",
    "\n",
    "testing = trainingTesting.filter(ee.Filter.gte('random', 0.6));"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[37.46862098849245, -122.11682996798393], controls=(WidgetControl(options=['position', 'transparent…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "image1 = ee.Image().paint(training, 0, 2)\n",
    "image2 = ee.Image().paint(testing, 0, 2)\n",
    "\n",
    "Map.addLayer(image1, {}, 'Training Points')\n",
    "Map.addLayer(image2, {}, 'Testing Points')\n",
    "\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Inspect the first element of T to make sure it's got all the data you expect.  \n",
    "\n",
    "The next step is to train G.  Make a list of the variable names, predictors followed by g:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['constant', 'SR_B1', 'SR_B2', 'SR_B3', 'SR_B4', 'SR_B5', 'SR_B6', 'SR_B7', 'Percent_Tree_Cover']\n"
     ]
    }
   ],
   "source": [
    "trainingList = ee.List(predictionBands).insert(0, 'constant').add('Percent_Tree_Cover')\n",
    "print(trainingList.getInfo())"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{'type': 'Feature', 'geometry': {'type': 'Point', 'coordinates': [-121.20608001850448, 37.08640098326397]}, 'id': '0', 'properties': {'Percent_Tree_Cover': 4, 'SR_B1': 0.03611500000000001, 'SR_B2': 0.057166249999999974, 'SR_B3': 0.0921325, 'SR_B4': 0.13462000000000002, 'SR_B5': 0.22901375000000002, 'SR_B6': 0.32127625, 'SR_B7': 0.21655625, 'constant': 1, 'random': 0.5108621334513251}}\n"
     ]
    }
   ],
   "source": [
    "# // For the linear regression model, we need to add a constant term to each\n",
    "# // training and testing data point\n",
    "def add_constant(feature):\n",
    "    return feature.set('constant',1.0)\n",
    "\n",
    "training_const = training.map(add_constant)\n",
    "\n",
    "testing_const = training.map(add_constant)\n",
    "                             \n",
    "print(training_const.first().getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In Earth Engine, linear regression is implemented as a Reducer.  This means that training G is a reduction of the T table, performed using the list of variables as an input.  The argument  tells the reducer how many of the input variables are predictors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[31.740127529197363, 108.03651912682838, -67.19241994600918, -176.0148952645832, 9.247451693023855, 62.498372582659364, -38.06823535317738, -65.1105474036357]\n"
     ]
    }
   ],
   "source": [
    "# // Define a linear regression classifier\n",
    "\n",
    "linRegression = training_const.reduceColumns(\n",
    "    reducer= ee.Reducer.linearRegression(8,1),\n",
    "  selectors= trainingList\n",
    ");\n",
    "\n",
    "coefficients = ee.Array(linRegression.get('coefficients'))\\\n",
    "    .project([0])\\\n",
    "    .toList();\n",
    "\n",
    "print(coefficients.getInfo())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Observe that the output is an object with a list of coefficients (in the order specified by the inputs list).  To use the coefficients to make a prediction in every pixel, first turn the output coefficients into an image, then perform the multiplication and addition that implements βp:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "# // Perform the linear regression on the landsat image\n",
    "# // (remember to add the constant band)\n",
    "\n",
    "linRegressionImage = ee.Image(1).addBands(landsat.select(predictionBands))\\\n",
    "    \t.multiply(ee.Image.constant(coefficients))\\\n",
    "    \t.reduce(ee.Reducer.sum())\n",
    "\n",
    "linRegressionImage = linRegressionImage.clip(region).updateMask(waterMask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[37.46862098849245, -122.11682996798393], controls=(WidgetControl(options=['position', 'transparent…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "Map.addLayer(linRegressionImage, {'min': 0, 'max': 100,\n",
    "                                  'palette':['red','yellow', 'green', 'blue', 'purple']}, 'OLS prediction')\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### - Non-linear regression functions"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the garden variety linear regression isn't satisfactory, Earth Engine contains other functions that can make predictions of a continuous variable.  Unlike linear regression, other regression functions are implemented by the classifier library.  \n",
    "\n",
    "- For example, a Classification and Regression Tree (CART, see Brieman et al. 1984) is a machine learning algorithm that can learn non-linear patterns in your data.  Reusing the T table (without the constant term), train a CART as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "cartRegression = ee.Classifier.smileCart(minLeafPopulation=10).setOutputMode('REGRESSION').train(\n",
    "      features= training, \n",
    "      classProperty= 'Percent_Tree_Cover', \n",
    "      inputProperties= predictionBands\n",
    "    );"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Make predictions over the input imagery (classify in this context is a misnomer):\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[37.46862098849245, -122.11682996798393], controls=(WidgetControl(options=['position', 'transparent…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "cartRegressionImage = landsat.select(predictionBands).classify(cartRegression, 'y_hat');\n",
    "\n",
    "Map.addLayer(cartRegressionImage, {'min': 0, 'max': 100,\n",
    "                                   'palette':['red','yellow', 'green', 'blue', 'purple']},'CART regression');\n",
    "\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "example https://code.earthengine.google.com/ade42831c23ae03ab9b26289d89b31f2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Use a Random Forest regression classifier"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "trees =110 #this value needs to be optimized\n",
    "\n",
    "rfRegression = ee.Classifier.smileRandomForest(numberOfTrees=trees, seed=0).setOutputMode('REGRESSION').train(\n",
    "      features= training, \n",
    "      classProperty= 'Percent_Tree_Cover', \n",
    "      inputProperties= predictionBands\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[37.46862098849245, -122.11682996798393], controls=(WidgetControl(options=['position', 'transparent…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# // Perform the RF regression on the landsat image\n",
    "rfRegressionImage = landsat.select(predictionBands).classify(rfRegression, 'y_hat');\n",
    "    \n",
    "# // Visualize the RF regression\n",
    "Map.addLayer(rfRegressionImage, {'min': 0, 'max': 100,'palette':['red','yellow', 'green', 'blue', 'purple']}, 'RF regression');\n",
    "\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using Support Vector Machines"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ee.Classifier.libsvm(decisionProcedure, svmType, kernelType, shrinking, degree, gamma, coef0, cost, nu, terminationEpsilon, lossEpsilon, oneClass)\n",
    "\n",
    "gamma_v =21 #this value needs to be optimized\n",
    "\n",
    "# // Create an SVM classifier with custom parameters.\n",
    "svrRegression = ee.Classifier.libsvm(svmType ='EPSILON_SVR',kernelType='RBF',gamma=gamma_v,cost=100).setOutputMode('REGRESSION').train(\n",
    "      features= training, \n",
    "      classProperty= 'Percent_Tree_Cover', \n",
    "      inputProperties= predictionBands\n",
    "    )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "5333a3eb84734b00802e84e750f3ddf4",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[37.46862098849245, -122.11682996798393], controls=(WidgetControl(options=['position', 'transparent…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# // Perform the RF regression on the landsat image\n",
    "svrRegressionImage = landsat.select(predictionBands).classify(svrRegression, 'y_hat');\n",
    "    \n",
    "# // Visualize the RF regression\n",
    "Map.addLayer(svrRegressionImage, {'min': 0, 'max': 100,'palette':['red','yellow', 'green', 'blue', 'purple']}, 'SVR regression');\n",
    "\n",
    "Map"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We need indicators of the performance of these models. We need a goodness of fit statistics. The most used is Root Mean Square Error (RMSE)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# ////////////////////////////////////////////////////////////////////////////////\n",
    "# // Function to compute the sum of squared differences for the results from \n",
    "# // any Regression. This function uses the actual values stored in the \n",
    "# // 'Percent_Tree_Cover' property of the 'regression' feature collection, and the predicted \n",
    "# // values stored in the 'prediction' property.\n",
    "def set_sqError(feature):\n",
    "    feature = ee.Feature(feature)\n",
    "    actualValue = ee.Number(feature.get('Percent_Tree_Cover'))\n",
    "    predict = ee.Number(feature.get('y_hat'))\n",
    "    sqError = (actualValue.subtract(predict)).multiply(actualValue.subtract(predict))\n",
    "    feature = feature.set('sqError',sqError)\n",
    "    return feature\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accuracy Assessment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Linear Regression RMSE: 11.686813744188472\n"
     ]
    }
   ],
   "source": [
    "# // Perform the linear regression on the test set (with constant term)\n",
    "\n",
    "def linearReg_applied(feature):\n",
    "    feature =ee.Feature(feature)\n",
    "    const =ee.Number(feature.get('constant')).multiply(coefficients.getInfo()[0])\n",
    "    b1 =ee.Number(feature.get('SR_B1')).multiply(coefficients.getInfo()[1])\n",
    "    b2 =ee.Number(feature.get('SR_B2')).multiply(coefficients.getInfo()[2])\n",
    "    b3 =ee.Number(feature.get('SR_B3')).multiply(coefficients.getInfo()[3])\n",
    "    b4 =ee.Number(feature.get('SR_B4')).multiply(coefficients.getInfo()[4])\n",
    "    b5 =ee.Number(feature.get('SR_B5')).multiply(coefficients.getInfo()[5])\n",
    "    b6 =ee.Number(feature.get('SR_B6')).multiply(coefficients.getInfo()[6])\n",
    "    b7 =ee.Number(feature.get('SR_B7')).multiply(coefficients.getInfo()[7])\n",
    "    y_hat=const.add(b1).add(b2).add(b3).add(b4).add(b5).add(b6).add(b7)\n",
    "    feature=feature.set('y_hat',y_hat)\n",
    "    return feature\n",
    "\n",
    "linRegressionTesting = testing_const.map(linearReg_applied)\n",
    "\n",
    "# // Calculate the RMSE for linear regression and print the result\n",
    "regression=linRegressionTesting.map(set_sqError)\n",
    "regression = ee.FeatureCollection(regression)\n",
    "MSE = regression.reduceColumns(ee.Reducer.mean(),['sqError'])\n",
    "RMSE =  ee.Number(MSE.get('mean')).sqrt()\n",
    "\n",
    "print('Linear Regression RMSE:',RMSE.getInfo())\n",
    "\n",
    "# print(coefficients.getInfo()[6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CART Regression RMSE: 11.022595887004154\n"
     ]
    }
   ],
   "source": [
    "# // Perform the CART regression on the test set\n",
    "cartRegressionTesting = testing.classify(cartRegression, 'y_hat');\n",
    "\n",
    "# // Calculate the RMSE for CART regression and print the result\n",
    "regression=cartRegressionTesting.map(set_sqError)\n",
    "regression = ee.FeatureCollection(regression)\n",
    "MSE = regression.reduceColumns(ee.Reducer.mean(),['sqError'])\n",
    "RMSE =  ee.Number(MSE.get('mean')).sqrt()\n",
    "\n",
    "print('CART Regression RMSE:', RMSE.getInfo());"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "RF Regression RMSE: 10.09936436245924\n"
     ]
    }
   ],
   "source": [
    "# // Perform the RF regression on the test set\n",
    "rfRegressionTesting = testing.classify(rfRegression, 'y_hat');\n",
    "\n",
    "# // Calculate the RMSE for RF regression and print the result\n",
    "regression=rfRegressionTesting.map(set_sqError)\n",
    "regression = ee.FeatureCollection(regression)\n",
    "MSE = regression.reduceColumns(ee.Reducer.mean(),['sqError'])\n",
    "RMSE =  ee.Number(MSE.get('mean')).sqrt()\n",
    "\n",
    "print('RF Regression RMSE:', RMSE.getInfo());"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "SVR Regression RMSE: 10.021654299969228\n"
     ]
    }
   ],
   "source": [
    "# // Perform the SVR regression on the test set\n",
    "svrRegressionTesting = testing.classify(svrRegression, 'y_hat');\n",
    "\n",
    "# // Calculate the RMSE for RF regression and print the result\n",
    "regression=svrRegressionTesting.map(set_sqError)\n",
    "regression = ee.FeatureCollection(regression)\n",
    "MSE = regression.reduceColumns(ee.Reducer.mean(),['sqError'])\n",
    "RMSE =  ee.Number(MSE.get('mean')).sqrt()\n",
    "\n",
    "print('SVR Regression RMSE:', RMSE.getInfo());"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### BONUS: Hyperparameters tuning\n",
    "\n",
    "A random forest is a collection of random trees the predictions of which are used to compute an average (regression) or vote on a label (classification).  Note that the only parameter to the classifier is the number of trees (10).  How many trees should you use?  Making that choice is best done by hyperparameter tuning.  For example, \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min RMSE = 10.045687889039481\n",
      "Min Trees = 135\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA4DklEQVR4nO3deXhU5fXA8e/JvicQQjaWsMkugVBxwwKKInVrq1W7oVLpYmtrbV1au2htf7ZVa1ttq1aFagtWrRWpC4hB3BEQNMi+CCGEQEiAJGQ/vz/uDUzCZM8kk5nzeZ555q5zz525c+ad9773vaKqGGOMCVwhPR2AMcYY37JEb4wxAc4SvTHGBDhL9MYYE+As0RtjTICzRG+MMQHOEn0QEZHPi8geESkTkYl+EM8KEflGD207WkReFJHDIvJMT8TQViIyyP3MQns6ltaIyC4ROa+n4zCNWaJvhj8esCKiIjK8Ey9xL/BdVY1T1Q+bef2PRSTEY9rdIjK/E9v0V5cDqUCyql7hOUNE/uYm1jIRqRaRGo/xl7s7UFXd7X5mde1dV0SmiUi9G/tREdksItf6Is7uJCLz3c+mzONxZTdu/5ci8lR3ba+zLNEHl8HAhlaWyQCu6oZYuow42nssDwa2qGpt0xmq+i03scYBvwGebhhX1Qs9thvWuci7TYG7LwnATcCjIjKyh2PqCr/z+FziVPXp9qzciz6/TrNE3wYico2IvC0ifxCRUhHZISJnutP3iEiRiMzxWH6+Wypc5pai3hCRwR7z/+iud0RE1ojIVI95oSLyExHZ7q67RkQGishKd5H1zZVeRCRERO4QkU/dmP4hIokiEikiZUCou/72Fnb3d8Cd3r4Ebukwv8m04/983FLOMyLylBv7xyJyiojc7sazR0TOb/Kyw0RklVuF8oKI9PV47dNF5B33PV8vItM85q0QkV+LyNtABTDUS7yj3eVKRWSDiFziTr8T+Dlwpftezm3h/Wj6mrtE5FYR+QgoF5GwVuJMFJHHRGSfiOwV5x9SqDtvuHtsHBaRgyLiNVGJSJY4/7bCPPb9V+4xeVRElopIv9ZiV8dLwCHgVPe1+ojIEhE5ICIl7vAAj223uC0R+Zp7vBWLyE+bxB0pIg+ISIH7eEBEIt1500QkX0RucY+NfSJymYjMFpEtInJIRH7S1s+lyXavF5Ft7mssFpEMj3kqIjeIyFZgqzvtIhFZ535+74jIqR7L3+p+bg3/hs4VkVnATzhx/KzvSJzdSlXt4eUB7ALOc4evAWqBa3GS5d3AbuAhIBI4HzgKxLnLz3fHz3Hn/xF4y+O1vwokA2HAzUAhEOXO+zHwMTASEGACTvUCgALDW4j5OmAbTtKLA/4DPOkxv7X1FRgBrAG+4U67G5jvDk8D8lt4n34JVAIXuPv2D2An8FMgHLge2Omx7gpgLzAOiAWeA55y52UCxcBsnALJTHc8xWPd3cBYd1vhTeIKd9+LnwARwAz3MxnpEetTbTgOGi3n7u86YCAQ3YY4/ws87O5ff2AV8E133kL3vQkBooCzm4khy/1swjz2fTtwihvDCuCeZtY9/pm527kEqAcmutOSgS8CMUA88Azw3yafkddtAWOAMk4c5/fjfE8ajoe7gPfc/U4B3gF+5RFXLc4PbsOxcQD4lxvHWJxjaWgz+zUfuNvL9BnAQWCSG9OfgZVNjvFlQF93fyYBRcAUnO/2HPczjsT5Du4BMjw+h2HtOX785dHjAfjrg5MT/VaPeePdAybVY1oxkO1xEC7ymBcH1AEDm9lWCTDBHd4MXNrMcq0l6uXAdzzGRwI1nEgQbUn0w3GS1m73YG9vol/mMe9iNxGEuuPx7jaS3PEVeCQonMRR7X7hbsXjR8qd/yowx2Pdu1rYl6k4P6AhHtMWAr/0iLWjif46j/Fm48Q5B1AFRHvMuxrIdYf/ATwCDGglhixOTvR3eMz/DvBKM+tOw0nspW4sdcAPWthWNlDiMd7stnCStOdxHut+fg3Hw3Zgtsf8C4BdHnEd83JsTPFYfg1wWTNxzsf5ISh1Hwfd6Y/hVOl4fvdqgCyPY3yGx/y/4v74eEzbDHwW57tQBJzHyQWJNh0//vKwqpu22+8xfAxAVZtOi/MY39MwoKplOH+XMwBE5GYR2ej+ZS8FEoGGv8MDcb4gHZEBfOox/ilOaTe1PS+izt/73cC8DsTQ9D05qCdOIh5zn72+TzjxhuO8F4OBK9y/06Xu+3Q2kN7Muk1lAHtUtb7J62e2dUda4LndluIcjLM/+zzmPYxTwgW4Bedf2yq3aum6dsRQ6DFcQeP3tKkCVU3CqaP/E06pFwARiRGRh93qlyPASiBJGrfwaW5bGTQ+zstxCjx4zG96PGZ4jBd7OTZa+k41da+qJrmPhu9Po226371iGn/uTT+/m5t8fgNxSvHbgB/gJPUiEVnkWQ3Um1ii952BDQMiEofzV7FAnPr4W4EvAX3cL+BhnC88OAfhsA5uswDnwG0wCOfv8X7vi7foDpxqhRiPaeWe424ySOnAa3sa6DE8CKf0dRDnfXjS44ucpKqxqnqPx/Itdb1aAAyUxidpB+FUFXWW53ZbinMPTim6n8e8BFUdC6Cqhap6vapmAN8E/iKda1XVctCqVTjH3ngRucydfDPOP78pqpqAUw0DJ47Hluyj8XEeg1MV1MDb8VjQoeDbrtE2RSTWjcnzc2/6+f26yecXo6oLAVT1X6p6tvuaCvzWy2v4PUv0vjNbRM4WkQjgV8D7qroH5y9qLU59ZJiI/BynpNXg78CvRGSEOE4VkYYvz368nHT0sBC4SUSGuD8uDS1GTmpZ0hpVXYFzrmCOx+QtQJSIfE5EwnF+DCLb+9pNfFVExrhJ4i7gWbeU9xRwsYhcIM4J6ij3BN6All/uuPdxfphuEZFwcU6QXgws6mS8TTUbp6ruA5YC94lIgjgny4eJyGcBROQKj/0pwUke7W5C2R6qWg3ch1PtAs7xeAwoFedE+C/a8XLPAhd5HOd30TinLATuEJEU9wTuz3HeL1/6F3CtiGS7J35/g/Pd29XM8o8C3xKRKe73LdY9vuNFZKSIzHBfpxLnfWr4fPYDWdL+1l49olcE2Uv9C+dLcwjIAb7iTn8VeBknaX6KcwB5/pW8H/g3ToI4glPnGO3O+yWwwP2L+SUv23wceBLn7/dO97W/14l9uAPnnwgAqnoYp4727zglpHIg3/uqbfYkTn1rIc4JyRvdbe0BLsU5mXoA5z36MW08Zt2EdglwIc4/hL8AX1fVTZ2Mt+l2Wovz6zgngz/BSebPcqL66TPA++K0iFoMfF9Vd3ZlfM14HBgkIhcDD+AcXwdxTpy+0tYXUdUNwA04x/o+nP3zPB7uBlYDH+EUGta603xGVZcDP8M5sb8P599xs82FVXU1zongB3Hi34ZzTg6cQsw9OO9NIU6VW0NLoIaL7IpFZG2X7oQPiHtiwXQhcS4wylfVO3o6FmOMsRK9McYEOEv0xhgT4FqtuhGRx4GLgCJVHedO6ws8jdO+dxfwJVUt8bLuTcA3cE4yfQxcq6qVXRi/McaYVrSlRD8fmNVk2m3AclUdgXORzm1NVxKRTJwTa5PdH4hQelkfKsYYEwha7dRHVVeKSFaTyZfiXNkGsADn6rlbm3n9aBGpwWl/3aY2tP369dOsrKabbFl5eTmxsbHtWqen9KZYweL1NYvXd3pTrNC5eNesWXNQVb1e19LR3ttS3TbCqOo+EenfdAFV3Ssi9+JcYXkMWKqqS5t7QRGZh3slZmpqKvfee2+7AiorKyMurqWL6PxHb4oVLF5fs3h9pzfFCp2Ld/r06Z82O7Mt/STg1MXneYyXNplf4mWdPsDrOFdOhuN07vTVtmwvJydH2ys3N7fd6/SU3hSrqsXraxav7/SmWFU7Fy+wWru4r5v9IpIO4D4XeVnmPJyeCg+oag1OT4pndnB7xhhjOqijiX4xJy6NnwO84GWZ3cDpbqdJApwLbOzg9owxxnRQq3X0IrIQ58RrP3FuOvELnMuC/y3ODRt2A1e4y2YAf1fV2ar6vog8i3PZcy3wIU6XrMYY0yk1NTXk5+dTWdlya+3ExEQ2buw95cu2xBsVFcWAAQMIDw9v8+u2pdXN1c3MOtfLsgU4fZk3jP+C9nWSZIwxrcrPzyc+Pp6srCycCgPvjh49Snx8fDdG1jmtxauqFBcXk5+fz5AhQ9r8unZlrDGm16msrCQ5ObnFJB+IRITk5ORW/8k0ZYneGNMrBVuSb9CR/Q6oRP+n5Vt5Y8uBng7DGGP8SkAl+kdW7uCNzZbojTG+FxoaSnZ2NuPGjePiiy+mtLQUgF27diEi/OxnPzu+7MGDBwkPD+e73/0uAJs3b2batGlkZ2czevRo5s1z7tr55ptvkpiYSHZ29vHHa6+91ulYAyrRx0eFcbSypqfDMMYEgejoaNatW0deXh59+/bloYceOj5v6NChLFmy5Pj4M888w9ixY4+P33jjjdx0002sW7eOjRs38r3vnbg/0NSpU1m3bt3xx3nnndfpWAMw0bf7rnnGGNMpZ5xxBnv3nrgtbXR0NKNHj2b16tUAPP3003zpSyduCrdv3z4GDDhxV8zx48f7NL6O9nXjl+KjwjlaZSV6Y4LJnS9u4JOCI17n1dXVERoa2u7XHJORwC8uHtv6gu42li9fzty5cxtNv+qqq1i0aBFpaWmEhoaSkZFBQYHTr+NNN93EjBkzOPPMMzn//PO59tprSUpKApzqm+zs7OOv89xzzzFs2LB274MnK9EbY0wHHDt2jOzsbJKTkzl06BAzZ85sNH/WrFksW7aMhQsXcuWVVzaad+2117Jx40auuOIKVqxYwemnn05VVRVwctVNZ5M8BGCJftfB8p4OwxjTjVoqefvygqmGOvrDhw9z0UUX8dBDD3HjjTcenx8REUFOTg733XcfGzZs4MUXX2y0fkZGBtdddx3XXXcd48aNIy8vzydxgpXojTGmUxITE/nTn/7EvffeS01N46rjm2++md/+9rckJyc3mv7KK68cX7awsJDi4mIyMzN9FmNAlegTosIt0Rtjut3EiROZMGECixYtYurUqcenjx07tlFrmwZLly7l+9//PlFRUQD8/ve/Jy0tDTi5jv6OO+7g8ssv71R8AZXo46PCqK6rp7Kmjqjw9p+AMcaYtiorK2s07lk1460a5pprruGaa64B4P777+f+++8/aZmpU6dy+PDhrg2UAKu6SYhyfresVG+MMScEVKKPj3K67bSLpowx5oQAS/ROif6IleiNCXjO3fOCT0f2O8ASvZXojQkGUVFRFBcXB12yb+iPvuEkblsF1MnYhGirozcmGAwYMID8/HwOHGi5E8PKysp2J8We1JZ4G+4w1R4BleitRG9McAgPD2/THZZWrFjBxIkTuyGiruGreFutuhGRx0WkSETyPKb1FZFlIrLVfe7jZb2RIrLO43FERH7QxfE3Em+tbowx5iRtqaOfD8xqMu02YLmqjgCWu+ONqOpmVc1W1WwgB6gAnu9UtK2IiwhDxE7GGmOMp1YTvaquBA41mXwpsMAdXgBc1srLnAtsV9VP2xtge4SECHER1ie9McZ4kractRaRLGCJqo5zx0tVNcljfomqnlR94zH/cWCtqj7YwjLzgHkAqampOYsWLWrrPgDOVWpxcXH8cEUFo/uGcv2pke1avzs1xNpbWLy+ZfH6Tm+KFToX7/Tp09eo6mSvM1W11QeQBeR5jJc2mV/SwroRwEEgtS3bUlVycnK0vXJzc1VV9fz739DrF3zQ7vW7U0OsvYXF61sWr+/0plhVOxcvsFqbyakdbUe/X0TSAdznohaWvRCnNL+/g9tql4Ro68HSGGM8dTTRLwbmuMNzgBdaWPZqYGEHt9NudpcpY4xprC3NKxcC7wIjRSRfROYC9wAzRWQrMNMdR0QyROQlj3Vj3Pn/8UXw3lif9MYY01irF0yp6tXNzDrXy7IFwGyP8QoguelyvmSJ3hhjGguovm7ArbqprAm6PjCMMaY5AZjow6ipUypr6ns6FGOM8QsBmOitvxtjjPEUcIk+wfqkN8aYRgIw0VuJ3hhjPAVcorceLI0xprEATPQNJXpL9MYYAwGZ6Bvq6K3qxhhjIIATvdXRG2OMI+ASfax78xGrujHGGEfAJfqQECE+0rpBMMaYBgGX6ME5IWt19MYY4wjQRG8lemOMaRCQiT7B7djMGGNMgCb6+KgwjhyzEr0xxkAAJ3q7y5QxxjgCNNGHWx29Mca4AjLRN9wg3G4+YowxAZro46PCqatXjtXU9XQoxhjT49pyc/DHRaRIRPI8pvUVkWUistV97tPMukki8qyIbBKRjSJyRlcG3xzrwdIYY05oS4l+PjCrybTbgOWqOgJY7o5780fgFVUdBUwANnYwznZp6MHyyDE7IWuMMa0melVdCRxqMvlSYIE7vAC4rOl6IpIAnAM85r5OtaqWdiLWNou3u0wZY8xx0pYTliKSBSxR1XHueKmqJnnML1HVPk3WyQYeAT7BKc2vAb6vquXNbGMeMA8gNTU1Z9GiRe3akbKyMuLi4gDYWlLHr9+v5Ic5kZyaEtau1+kOnrH2Bhavb1m8vtObYoXOxTt9+vQ1qjrZ60xVbfUBZAF5HuOlTeaXeFlnMlALTHHH/wj8qi3by8nJ0fbKzc09Pry58IgOvnWJLl63t92v0x08Y+0NLF7fsnh9pzfFqtq5eIHV2kxO7Wirm/0ikg7gPhd5WSYfyFfV993xZ4FJHdxeuyTYXaaMMea4jib6xcAcd3gO8ELTBVS1ENgjIiPdSefiVOP4nN18xBhjTmhL88qFwLvASBHJF5G5wD3ATBHZCsx0xxGRDBF5yWP17wH/FJGPgGzgN10cv1cxEaGEhoiV6I0xBmj1TKWqXt3MrHO9LFsAzPYYX4dTV9+tRIS4yDDrk94YYwjQK2PB+qQ3xpgGAZzorU96Y4yBgE70YXbBlDHGEMCJPsG6KjbGGCCgE32YVd0YYwwBnOjtZKwxxjgCONE7J2PVbj5ijAlyAZzow6hXKK+2m48YY4JbACf6hv5urJ7eGBPcAjjR212mjDEGAjjRJ0Rbid4YYyCAE73dZcoYYxwBm+gTGhK93TfWGBPkAjbRx9vNR4wxBgjoRG8nY40xBgI40UeHN9x8xKpujDHBLWATvYhYNwjGGEMAJ3po6MHSSvTGmODWlnvGPi4iRSKS5zGtr4gsE5Gt7nOfZtbdJSIfi8g6EVndlYG3hZXojTGmbSX6+cCsJtNuA5ar6ghguTvenOmqmq2q3X7vWOfmI1aiN8YEt1YTvaquBA41mXwpsMAdXgBc1rVhdY14u/mIMcYgbenGV0SygCWqOs4dL1XVJI/5Jap6UvWNiOwESgAFHlbVR1rYxjxgHkBqamrOokWL2rUjZWVlxMXFNZr26EdVbDpUx33TYtr1Wr7mLVZ/ZvH6lsXrO70pVuhcvNOnT1/TbM2Jqrb6ALKAPI/x0ibzS5pZL8N97g+sB85py/ZycnK0vXJzc0+adufiDTr6Zy+3+7V8zVus/szi9S2L13d6U6yqnYsXWK3N5NSOtrrZLyLpAO5zUTM/IgXucxHwPHBaB7fXIakJkVRU11nLG2NMUOtool8MzHGH5wAvNF1ARGJFJL5hGDgfyGu6nC+lJUYBUHi4sjs3a4wxfqUtzSsXAu8CI0UkX0TmAvcAM0VkKzDTHUdEMkTkJXfVVOAtEVkPrAL+p6qv+GInmpOa4Cb6I5bojTHBK6y1BVT16mZmnetl2QJgtju8A5jQqeg6KS3BSvTGGBPQV8Y2VN3stxK9MSaIBXSijwoPJSkm3KpujDFBLaATPTjVN4WHq3o6DGOM6TEBn+hTE6Ks6sYYE9QCPtGnJURZ1Y0xJqgFfKJPTYziYFkVNXX1PR2KMcb0iIBP9GkJUajCgaNWT2+MCU6Bn+gTIwHYZ23pjTFBKuATfcPVsXZC1hgTrAI+0acnRgN2dawxJngFfKLvExNORFiIleiNMUEr4BO9iJCaEGlNLI0xQSvgEz00XB1rid4YE5yCItHb1bHGmGAWFIm+4epYbcP9cY0xJtAER6JPjKKypp4jx2p7OhRjjOl2QZHoG9rS7ztyrIcjMcaY7hcUid7uHWuMCWZtuWfs4yJSJCJ5HtP6isgyEdnqPvdpYf1QEflQRJZ0VdDtlWZXxxpjglhbSvTzgVlNpt0GLFfVEcByd7w53wc2dii6LtI/wenvxm5AYowJRq0melVdCRxqMvlSYIE7vAC4zNu6IjIA+Bzw946H2HmRYaEkx0bYRVPGmKDU0Tr6VFXdB+A+929muQeAW4Ae7wze2tIbY4KVtKVtuYhkAUtUdZw7XqqqSR7zS1S1T5N1LgJmq+p3RGQa8CNVvaiFbcwD5gGkpqbmLFq0qF07UlZWRlxcXLPz/7CmkpJK5a6zotv1ur7QWqz+xuL1LYvXd3pTrNC5eKdPn75GVSd7namqrT6ALCDPY3wzkO4OpwObvazzf0A+sAsoBCqAp9qyvZycHG2v3NzcFuff9txHOumupe1+XV9oLVZ/Y/H6lsXrO70pVtXOxQus1mZyakerbhYDc9zhOcALXn5AblfVAaqaBVwFvK6qX+3g9jotLSGK4vJqqmrreioEY4zpEW1pXrkQeBcYKSL5IjIXuAeYKSJbgZnuOCKSISIv+TLgjmq401TREWt5Y4wJLmGtLaCqVzcz61wvyxYAs71MXwGsaGdsXarh6tjCI5UM7BvTk6EYY0y3CoorY8GujjXGBK+gSfTpCU5rG2tiaYwJNkGT6BOiw4gKD7ESvTEm6ARNoheR4/3SG2NMMAmaRA92dawxJjgFVaJPS7QSvTEm+ARXok+IYv+RKruloDEmqARVok9NiKK6tp6SipqeDsUYY7pNUCV6a0tvjAlGQZXoM5KctvRrd5f0cCTGGNN9girRn5qZyOTBffjdK5sospOyxpggEVSJPiRE+N3lp1JVW89Pns+zk7LGmKAQVIkeYGhKHD86fySvbdzPqxv293Q4xhjjc0GX6AGuO3sI0eGhfLCr6a1wjTEm8ARlog8NEdIS7SpZY0xwCMpED9A/PtJuQmKMCQpBm+hTrYMzY0yQCNpE31B1Yy1vjDGBLmgTff/4SKpq6zlyrLanQzHGGJ9qy83BHxeRIhHJ85jWV0SWichW97mPl/WiRGSViKwXkQ0icmdXB98ZnveQNcaYQNaWEv18YFaTabcBy1V1BLDcHW+qCpihqhOAbGCWiJze8VC7VkO/N9byxhgT6FpN9Kq6Emja4PxSYIE7vAC4zMt6qqpl7mi4+/CbCvHUeEv0xpjgIG05GSkiWcASVR3njpeqapLH/BJV9VZ9EwqsAYYDD6nqrS1sYx4wDyA1NTVn0aJF7dqRsrIy4uLi2rx8dZ0yb1kFXxgRziXDItq1rc5qb6w9zeL1LYvXd3pTrNC5eKdPn75GVSd7namqrT6ALCDPY7y0yfySVtZPAnKBcW3ZXk5OjrZXbm5uu9eZcOeresfzH7d7vc7qSKw9yeL1LYvXd3pTrKqdixdYrc3k1I62utkvIukA7nNRSwuraimwgpPr+ntUarxdHWuMCXwdTfSLgTnu8BzghaYLiEiKiCS5w9HAecCmDm7PJ/onRFqiN8YEvLY0r1wIvAuMFJF8EZkL3APMFJGtwEx3HBHJEJGX3FXTgVwR+Qj4AFimqkt8sRMd1XAPWWOMCWRhrS2gqlc3M+tcL8sWALPd4Y+AiZ2KzsdSE6I4UFZFXb0SGiI9HY4xxvhE0F4ZC5CaEEldvVJcZqV6Y0zgCupE3z+hoS29JXpjTOAK6kSflmAXTRljAl9QJ3rr78YYEwyCOtH3i4sgRKDIEr0xJoAFdaIPCw2hX1yk1dEbYwJaUCd6gIykaHYfqujpMIwxxmeCPtGPTI1n8/6jdqcpY0zAskSfFs+h8moOWFt6Y0yACvpEPyotHoDNhUd7OBJjjPGNoE/0Iy3RG2MCXNAn+uS4SPrFRVqiN8YErKBP9AAj0+LYvP/kRL/HWuMYYwKAJXpgZGoCW/Yfpa7+RMubd7YfZOrvclm55UAPRmaMMZ1niR7nhGxlTX2j9vSvfeLcNOvfq/f0VFjGGNMlLNHjeUL2yPFpb2xxEv3ST/Zz+FhNj8RljDFdwRI9cEpqPCKwubAMgPySCrYfKOfS7Ayqa+t5+eN9PRyhMcZ0nCV6IDoilMF9Y9i83ynRr9xyEIDvzRjOsJRY/rN2b0+GZ4wxndKWe8Y+LiJFIpLnMa2viCwTka3ucx8v6w0UkVwR2SgiG0Tk+10dfFcalZbAB7tKKK2o5o0tRWQmRTMsJY4vTBrAql2H2F1sLXCMMb1TW0r084FZTabdBixX1RHAcne8qVrgZlUdDZwO3CAiYzoRq09df85QDlfU8K2n1vD2tmLOOSUFEeHzEzMJCxG+t3AthYetO2NjTO/TaqJX1ZXAoSaTLwUWuMMLgMu8rLdPVde6w0eBjUBmZ4L1pZzBfbjni+N5b8chyqpq+ewpKYDTu+Vfv5rDtqIyLnnwLf7x7i5e+2Q/RUct6RtjeoeO1tGnquo+cBI60L+lhUUkC5gIvN/B7XWLL0wawI0zhpMUE86Zw5OPT585JpXnvnMm0RGh/PyFDXzjH6u59okPejBSY4xpO2lL97xuol6iquPc8VJVTfKYX6KqJ9XTu/PigDeAX6vqf1rYxjxgHkBqamrOokWL2rEbUFZWRlxcXLvWaU5dvRIaIidNr1flcJXywvYaVubX8vDMGMK9LNearoy1O1i8vmXx+k5vihU6F+/06dPXqOpkrzNVtdUHkAXkeYxvBtLd4XRgczPrhQOvAj9sy3YaHjk5Odpeubm57V6noxav26uDb12ieXtLO7R+d8baFSxe37J4fac3xarauXiB1dpMTu1o1c1iYI47PAd4oekCIiLAY8BGVb2/g9vxS6PTnQusNu2zjtCMMf6vLc0rFwLvAiNFJF9E5gL3ADNFZCsw0x1HRDJE5CV31bOArwEzRGSd+5jtk73oZlnJsUSEhXjtCM0YY/xNWGsLqOrVzcw618uyBcBsd/gtoP0V2L1AWGgIp6TGsXHfkdYXNsaYHmZXxnbQqLQENlkf9saYXsASfQeNSovnwNEqDtq9Zo0xfs4SfQeNSksA7BaExhj/Z4m+g0a5LW+snt4Y4+8s0XdQP/des1ZPb4zxd5boO2F0ejybCq1Eb4zxb5boO2FUWjxb9pc1utesMcb4G0v0nTCifzzVtfXkl1hf9cYY/2WJvhOG9Y8FYPuBsh6OxBhjmmeJvhOG9nN6mdteVN7DkRhjTPMs0XdCn9gIkmMjrERvjPFrlug7aWhKLDsOWIneGOO/LNF30rCUOCvRG2P8miX6ThqWEkdxeTUl5dU9HYoxxnhlib6TGlre7DhopXpjjH+yRN9Jw1Ks5Y0xxr9Zou+kAX1iiAgNsXp6Y4zfskTfSaEhQla/GEv0xhi/ZYm+CwxLibMmlsYYv9WWm4M/LiJFIpLnMa2viCwTka3uc5+2rhuIhqXE8emhCqpr63s6FGOMOUlbSvTzgVlNpt0GLFfVEcByd7yt6wacYf1jqatXHlm5nZq6tiX7j/JLueaJVSxctdt6vzQmiFTV1nX7NsNaW0BVV4pIVpPJlwLT3OEFwArg1jauG3AuGJvG9JEF3Lt0C/9dV8ADV2YzLjPR67Kqyqu7anh22TuEiLBi8wHmv72LB788kRGp8d0cuTHGF4qOVFJ0tIrUhCgqqmt5bWMRb209wMZ9Ryk8Uslnsvpw9WmDCBFhfX4pidHhfHHSAJ/FI6qtlybdZL1EVce546WqmuQxv0RVm6u+abRuC9uYB8wDSE1NzVm0aFEbd8FRVlZGXFxcu9bpauuKavnHJ9VU1Cg3ZEcyPqXx72hNvfLYx1W8t6+OSf1DuW5cJBsP1fF4XhXj+4XyneyoHoq8Zf7w3raHxetbvSne7oq1XpWiCmXH4Xre21fLxwfqaJpZ02OFrMQQ+kSGsGZ/LfsrnCUiQqGmDhQ4JVG5ZUosYSHS7himT5++RlUne5vXaom+u6jqI8AjAJMnT9Zp06a1a/0VK1bQ3nW62jTg6iOVXPPEBzzw4VH+dNVYPndqOgCHj9XwzSdX896+Cr44Ipx7r5uJiHARUPTMel7ZUMhZU88hPNT/zo/7w3vbHt7iratXVu08xIjUOPrFRfZMYM0IhPfXX7UWa+HhSmrq6hnYNwaAd7Yd5LG3drKzuJyiI1WcNTyZ70wbzqGKah5/ayf7DlcyPjORsRkJZCXHEhUeyst5+3jp432UVNQAkJYQxQ3ThzAuM4H9R6oQgc+eksLg5Njj262vVz7cU0psZCgj+sez/0glz3+4l/fytnPejOld/j50NNHvF5F0Vd0nIulAUVcG1ZulJkTx72+ezrVPfMCPn13PqPR4MpOiufaJVXy89zB/vCqbxNKtiJz4xT53dCrPrMnng12HOHNYvx6MPjDtLT3GTU+vY9XOQ4QITBmSzPXnDGHGqFQADpVX8+bWA6zfc5iDZVVMHdGPGaP60zc2otHn1KCmrp6Fq3bzwa4Spgzpy7SRKQzoE9Pdu2W8UFVe3VDI9gPlfH5iZqN5JeXVvL6piC1FR9ldXMFH+YfZW3oMgPGZifSNjeCNLQdIT4xi4qAkpgzpy/8+2serG/YDTgIfl5nAO9sP8vyHe4+/bnR4KDPHpHLmsGTGZSYyOj2B0FZK5CEhQs7gE5UgGUnR3DB9OGMlv6veikY6mugXA3OAe9znF7osogAQHxXOg1+exOw/vckN/1zLkH6xfLinlL98eRIXjk9nxYqtjZafOqIfEaEhvL6xyG8SfWlFNYKQGBPu822pqteE2hXe21HMvH+spq5eufOSsRSXVfHfdQVcN381547qT1JMBC9+VEB1bT1R4SHER4WzeH0BAGEhQmxkGFn9YhmTnkC/uAgAXvp4H9sPlNMnJpwX3WVH9I9j+qj+fGXKoEYlN9N1aurqeSWvkBWbD5CaEElWv1iG9ItlcHIM5VV1bC8q45GVO1i16xAA9y/bwui+ISzas4bSY9Ws3lVCbb0SERrCgL7RTBiYyNyzh1BXr7z4UQHr80u5ZdZIrjtrCFHhoQD8ZPZo/rN2L0kx4Vw4Lp2IMOcf96HyavYcquBQRTWnZfUlNtJvKke8ajU6EVmIUyvRT0TygV/gJPh/i8hcYDdwhbtsBvB3VZ3d3Lqq+pgP9sPvpCVGce8Vp3Ld/NVsKjzKHZ8bzYXj070uGxsZxunDknl9UxF3XDSmmyM9WV29csXf3qW6rp5Xf3COz7ZTX6/86n+f8OL6ffz84jFcMiGjS1//k4IjXL9gNf0TInn8ms8cT8DfnTGC+e/s5I+vbUWBKycP5IrJAxidnkBYiLCh4AhvbTvIkWM1HKmsYXtROS/n7ePwsRpUna6pH/36ZM4b3Z8dB8vJ3VR0/KT642/t5CtTBnHh+HSSYyPI7BNNTIR/JwF/daSyhj++tpX1e0qJCAthx4FyCo9UkhQTTlllLbVeWqslx0bwm8+P56zhyfzr/d28uHYn1QfLiAoP5RtTh/K58emMyTi5xH39OUO9xhAfFc6cM7NOmt43NoK+sRFdsp/doS2tbq5uZta5XpYtAGa3Yd2gMGNUKr+8eAzHauqZe/aQFpc9d1R/frF4AzsOlDE0pXtOdB2rrmP5pv1cMDat0bmBV/IK2VrkXOn7tze2k+2DPFVTV8/N/17P4vUFZCZFc+PCD3l1QyH3Xj6B6IjQTr/+nkMVXPPEKuKiwnhy7hQykqKPz4sIC2HeOcP48pTBCJxUGhuXmdhsqylo/A9kWEocw1Li+MbUoRQdqeSPy7fy1Pu7WfDup4Dzr2D8gESmDu/HFZMHHq8LDhbVtfXsKi7naGUNVbX1jM9MJD6q8b/Ed7cX89zafCqqa6mrV1LiI0mOjWThqt0cKKti8uA+VNfWMy4zgV9/fhzTR/anTpW9JcfYWVzO7uIK559Xcgyj0xOOf563zx7NGTH7mTbtsz2x637Fiho+ds1ZLSf4BjPcRP/6pqJuS/S/fWUT89/ZxbenDePWWaMAJ4n9+fWtDE2JZXRaAn9ZsZ27z4xkz6EKNhQcJntgH9ISO9c66NPicn78zEes2nWIW2aNZN7UoTy8cgf3Lt1MXEQYv7381E69/paSOn70l7eprq3n2W+f2SjJe4rr4N/t5qqZ+idE8evPj+eG6cPZebCc4vJqNhce4b0dh3gwdxt/zt1GzqA+hIeGUFNXT0xkGInR4SRU1zD2aBUp8Y1PEtfVa6t1vXDiHMMZQ5Ppn9DzLbe2FR3lhXUFLN2wn+0HyhqVvIemxPLct86kT2wEmwuP8ovFeby34xBJMeH0i4skROC9HYc4fKyGsRkJPPr1yUwYmHTSNkIQsvrFktXPqsnawhK9nxjYN4ZRafG8sK6AuWcP8VmddYPNhUd58r1P6RsbwV9XbOf0ocl89pQUlm8sYlPhUe67YgJnDe/His1F3PnuMcpX5h5fd1hKLP3jo4iJCKVfXCTpSVFcOC6dkWnNXwew82A5eXsP88m+I8x/exdhocIfrpzA5yc6bYdvmD6cY9V1PJi7jSlD+/KFdrYprqiuZd3uUt7adpCHV1UyKDmWR76W0yPXJmQkRZ/4cXGrowpKj/H0B3t4Y8sBoJ6IsBAOH6thx4Ey8kuqWfib17hgbBo3nz+S+Kgw7nrxE3I3F3H77NF8dcqgZo+Hqto65i74gA93lzqbG5DIjy4YydQRKd2wp41V1tTxy8UbWPTBHkIETh+azHljhnJKajx9YiIoqajmx89+xNwFH/DlKYP52X/ziIkI5ecXjeHLUwYdrxcHKK+qJSYi1Offg2Bhid6PzDkzi9v/8zFvbyvm7BG+Oymrqtz54gbiIsP4341nc83jH/DDp9cxY1R/3tlezMC+0VySnUF4aAi/uGQsf1+ex+VnnMKkQX1Yu7uEVTudElfB4RrW5zstVea/s4sl3zubAX1i+Nsb23n4je1cc+YQvvSZAfz59W386/3dx7d/7qj+3P35caQnNi5p/+C8EazadYifPp/Hlv1lZCRFMX1k/0bVHZ8Wl/PE27t4Z/tBpgxJ5qzh/cjdVMQL6/dSWVOPCExMCeWJb59FYrTvTyS3VUZSNDfNPIWbZp5y0rx/vvg6e8IyePLdXby6oZDo8FBq65VR6Qn87L95LN1QyNiMROrq64mJCCMpJpzR6Ql8Jqsvv1y8gQ93l/KLi8dQUV3HM6v38LXHVvH5iZlMG5lCVHgotXVKWVUNw/vHkTO4b5fuV0V1LZ8WV7D/SCV/WLaF9fmH+eY5Q5k7dQj940/+dxERGsJ3/rWWtbtLmTy4D3/5yiSv/0L8/eRmb2Pvph/5wqRMHnhtCw/lbvNJoj9YVsUreYWs+bSEd7YXc9elY0lPjOahr0zk2vkfsHLrAdISovjRBSOP19l/afJA+pdtZ9o5wwCYnNWXee5wg50Hy7nkz2/xnX+u5YuTBnDPy5vISo7hD69t4Q+vbSE0RPjG2UO4fPKA422PvQkLDeHPV09k7oIPeOytHdTUKVHhG/nhzFMYn5nEE2/vZNnG/YSFCJMG9eGZNXt48r1PiQoP4dIJmcwan8akQX348P23/SrJtyYzPoSvTBvF9VOH8NcV2zlYVsUPzjuFQX1j+Me7u7hv6Rbe33GI0BDhWM2Jy+eTYsIprajh29OGca1bRTj37CE8lLuNv72xvVETwAZfO30wt88eRUxEGPX1ynNr83nsrZ2kxEcycVAfvjgp86RWQ3sOVbB4fQHlVc4J0LpDNfTbe5hX8gpZ8O4ujlbWAk5V2MNfy+GCsWnN7uuF49P5w5ey2XGwnO9OH368FYvxLUv0fiQyLJTrpw7l7v9tZO3uEpKiw3l9UxEf7i5l58Fy7r9yAqPSEjr02rV19Xzl0ffZvP8oCVFhXJadwZdPGwTA8P7xvHnLjA7HPaRfLPd9aQLznlzDR/mHmTYyhUe/Ppm8vYd56eN9XJqd2eLJTU+pCVEs+d5U6uuVPSUV3P2/jfzmpU0A9IkJ54Zpw/naGYOPX1q+bk8pYzMSe1Vib05yXORJra6uOWtIo/M8dfVKaUU17+4o5uWPC4mLDONH5488Pj8qPJSbzx/JN84eSnF5FRXVdYSHhhATEcqCd3bx2Ns7WfJRAcNS4iirqmVT4VHGpCdwsKyaB1/fyqMrd3DHRaO56jOD2HfYqW56eOUOqmvrCQ8VRITq2noey3sLEZg1No2LTs2gf0Ikw1Li2tQS5bIm7duN71mi9zNXnzaIB3O38fXHVlFW5ZSUMpOiOXyshl//byNPzp3S4vqfFpdzz8ubmDYyhUsmZB5vwfLMmnw27z/KA1dmc2l2RpfXfZ4/No3bLhzF2k9LeOCqbMJDQ5g4qA8TB3ntGaNVISHCYLeePXdzESXlNXzu1PRG/wZiIsL85rqD7hIaIiTHRXLRqRlcdGrzzVETY8JPugbijovGcN6YVJ5ZnU9+SQXhoSE8cGU2l0zIICRE2Hf4GLc8+xE/fT6PXy7eQE2dcxL1suwMbpk1ioykaOrrlYUv5RKTeQrjMxMZ3t/6Z+oNLNH7mdjIMG65YBTPf5jPBWPTmD0+nYykaP7+5g7u/t9G3tp68Hi1jqrylxXbyS+p4OcXjUUEvvnkGjbvP8rLeYX85qVN3HTeCL6YM4D7lm5h8uA+PknyDb712WGtL9ROInL8ClbTeacPTeb0ocle56UnRrPg2tN4Zs0ethU5zXwnDEhiTMaJf5EhIUJmXAjTJvquAy7T9SzR+6EvTxnEl6cMajTta2cM5om3d/F/L2/kxWFnU6fK7f/5mGfXOJdMf1JwhEHJsWzef5THr/kMMeGhPJi7jV+++AkP5jr1vo9+PcdaMZgWhYQIV35mUOsLml7FEn0vERkWyo8uOIWbnl7P9PtWUFVTT+GRSm467xRGp8dz46IPWZ9/mBtnDGf6yP4AnDakL8+uyeeuJZ/whUmZHa5GMcb0bpboe5FLJ2Sy80A52w+Ug8D5Y1K5NNs5sfXst87kza0HmedxKbeIcMXkgVw8IcMve8U0xnQPS/S9SEiI8EOPFhaeWrpsv7nmjMaY4GDFPGOMCXCW6I0xJsBZojfGmABnid4YYwKcJXpjjAlwluiNMSbAWaI3xpgAZ4neGGMCnKiefIPdniYiB4BP27laP+CgD8Lxhd4UK1i8vmbx+k5vihU6F+9gVfV6azG/TPQdISKrVXVyT8fRFr0pVrB4fc3i9Z3eFCv4Ll6rujHGmABnid4YYwJcICX6R3o6gHboTbGCxetrFq/v9KZYwUfxBkwdvTHGGO8CqURvjDHGC0v0xhgT4Hp9oheRWSKyWUS2ichtPR0PgIg8LiJFIpLnMa2viCwTka3ucx+Pebe78W8WkQu6OdaBIpIrIhtFZIOIfN/P440SkVUist6N905/jtcjhlAR+VBElvh7vCKyS0Q+FpF1IrK6F8SbJCLPisgm9zg+wx/jFZGR7nva8DgiIj/ollhVtdc+gFBgOzAUiADWA2P8IK5zgElAnse03wG3ucO3Ab91h8e4cUcCQ9z9Ce3GWNOBSe5wPLDFjclf4xUgzh0OB94HTvfXeD3i/iHwL2CJPx8Pbgy7gH5NpvlzvAuAb7jDEUCSP8frxhEKFAKDuyPWbt05H7xZZwCveozfDtze03G5sWTRONFvBtLd4XRgs7eYgVeBM3ow7heAmb0hXiAGWAtM8ed4gQHAcmCGR6L353i9JXq/jBdIAHbiNizx93g9tns+8HZ3xdrbq24ygT0e4/nuNH+Uqqr7ANzn/u50v9kHEckCJuKUkv02XrcaZB1QBCxTVb+OF3gAuAWo95jmz/EqsFRE1ojIPHeav8Y7FDgAPOFWjf1dRGL9ON4GVwEL3WGfx9rbE714mdbb2ov6xT6ISBzwHPADVT3S0qJepnVrvKpap6rZOCXl00RkXAuL92i8InIRUKSqa9q6ipdp3X08nKWqk4ALgRtE5JwWlu3peMNwqkn/qqoTgXKc6o/m9HS8iEgEcAnwTGuLepnWoVh7e6LPBwZ6jA8ACnooltbsF5F0APe5yJ3e4/sgIuE4Sf6fqvofd7LfxttAVUuBFcAs/Dfes4BLRGQXsAiYISJP4b/xoqoF7nMR8DxwGv4bbz6Q7/6rA3gWJ/H7a7zg/ICuVdX97rjPY+3tif4DYISIDHF/Ja8CFvdwTM1ZDMxxh+fg1IU3TL9KRCJFZAgwAljVXUGJiACPARtV9f5eEG+KiCS5w9HAecAmf41XVW9X1QGqmoVzfL6uql/113hFJFZE4huGceqS8/w1XlUtBPaIyEh30rnAJ/4ar+tqTlTbNMTk21i7+ySED05qzMZpKbId+GlPx+PGtBDYB9Tg/CrPBZJxTshtdZ/7eiz/Uzf+zcCF3Rzr2Th/Bz8C1rmP2X4c76nAh268ecDP3el+GW+T2Kdx4mSsX8aLU+e93n1saPhO+Wu87vazgdXuMfFfoI+/xovTgKAYSPSY5vNYrQsEY4wJcL296sYYY0wrLNEbY0yAs0RvjDEBzhK9McYEOEv0xhgT4CzRG9MGIvKQ2+PgJyJyzKMHwst7OjZjWmPNK41pB7c/oCWqOq7J9FBVreuZqIxpmZXoTcARkSy3X/JHxemzfqmIRIvIChGZ7C7Tz+2WABG5RkT+KyIvishOEfmuiPzQ7STrPRHp28x2ponTl/+/gI/dztZ+LyIfiMhHIvJNj2V/7DG9oQ/9WBH5nzh96+eJyJW+f3dMMArr6QCM8ZERwNWqer2I/Bv4YivLj8PpuTMK2AbcqqoTReQPwNdxeqD05jRgnKrudHt6PKyqnxGRSOBtEVnqxjLCXVaAxW5HYSlAgap+DkBEEjuxv8Y0yxK9CVQ7VXWdO7wG5/4ALclV1aPAURE5DLzoTv8Yp9uF5qxS1Z3u8PnAqR719ok4Cf589/GhOz3Onf4mcK+I/BanOujNtuyYMe1lid4EqiqP4TogGqjlRHVlVAvL13uM19Py96TcY1iA76nqq54LuLeA+z9VfbjpyiKSg9O30P+JyFJVvauFbRnTIVZHb4LJLiDHHfZFa5lXgW+73T4jIqe4PUC+Clzn9vmPiGSKSH8RyQAqVPUp4F6c7nWN6XJWojfB5F7g3yLyNeB1H7z+33GqiNa63T8fAC5T1aUiMhp415lMGfBVYDjwexGpx+np9Ns+iMkYa15pjDGBzqpujDEmwFmiN8aYAGeJ3hhjApwlemOMCXCW6I0xJsBZojfGmABnid4YYwLc/wMdlEinKbVZ8wAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "numTrees = ee.List.sequence(5, 700, 5)\n",
    "\n",
    "\n",
    "def trees(t):\n",
    "    rfRegression = ee.Classifier.smileRandomForest(numberOfTrees=t, seed=0).setOutputMode('REGRESSION').train(\n",
    "    features= training, \n",
    "    classProperty= 'Percent_Tree_Cover', \n",
    "    inputProperties= predictionBands)\n",
    "    \n",
    "    rfRegressionTesting = testing.classify(rfRegression, 'y_hat')\n",
    "\n",
    "    # // Calculate the RMSE for RF regression and print the result\n",
    "    regression=rfRegressionTesting.map(set_sqError)\n",
    "    regression = ee.FeatureCollection(regression)\n",
    "    MSE = regression.reduceColumns(ee.Reducer.mean(),['sqError'])\n",
    "    RMSE =  ee.Number(MSE.get('mean')).sqrt()        \n",
    "    return RMSE\n",
    "\n",
    "\n",
    "rmse_trees=numTrees.map(trees)\n",
    "value_info = rmse_trees.getInfo()\n",
    "\n",
    "# Extracting tree value for minimum RMSE\n",
    "min_rmse=min(value_info)\n",
    "print(\"Min RMSE =\", min_rmse)\n",
    "\n",
    "index=value_info.index(min(value_info))\n",
    "trees_range = numTrees.getInfo()\n",
    "min_trees =trees_range[index]\n",
    "print(\"Min Trees =\", min_trees)\n",
    "\n",
    "# print(rmse_gama.getInfo())\n",
    "\n",
    "# Plotting the interations\n",
    "\n",
    "import pandas as pd\n",
    "df =pd.DataFrame(value_info,columns=['RMSE'])\n",
    "df['numTrees'] = numTrees.getInfo() \n",
    "\n",
    "ax =df.plot.line(x='numTrees', \n",
    "             y='RMSE',\n",
    "             title= 'Impact of Number of Trees in Random Forest'\n",
    "             )\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "same for svr"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min RMSE = 9.936551386540803\n",
      "Min Gamma = 6\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAAuaklEQVR4nO3deZiU9Znv//fdVb130w0NNjvtggRtBQV3TUAlItFoMr+4ZBm3/DgnE5fxJFdiTiZxzmROTuYk4/yS35iZcRKXOeNIEnWiEuISBDGJSQRFaEQEBaRlh6bthd7v88fzdFs01fRCNVXU83ldV131bN967urlrm/dz/I1d0dERLJXTroDEBGR4aVELyKS5ZToRUSynBK9iEiWU6IXEclySvQiIllOiV6ykpl9ysy2mVmjmZ2V7niGk5l9zsyeT3cckrmU6AUz22Jml6c7jkRm5mZ2ylG8xA+A2929xN1fT/L6Zma3m9kaM2s2s51mttzMbjiKfaaFuz/q7h8fSlszKzezB8P332Bmb5vZ18N1b5nZrUna3GVmK8Pp5WbWEn6g7jWzJ81s3NG9I0k1JXrJVlOAdUdY/yPgL4GvABXABOCvgPnDHllm+QegBJgOlAGfBN4J1z0C/HmSNl8I13W73d1LgFPC1/rBsEUrQ+PuekT8AWwBLg+nbwZ+R5AADgDvAheGy7cBu4GbEto+DPwz8ALQALwETElY/8Ow3QfAKuCShHUx4L8TJJaGcP0kYAXgQBPQCFyfJOYcgsS8NYzp3wgSVX7Yprv9O0nangp0ArP7+bncAqwPY3sX+C8J6+YAtcDXwv3vAK4FFgBvA/uB/56w/V8DvwD+PXy9tWEc3wjbbwM+PpB9J4nzZuC3CfMO/FdgI1AH3A9YH21rgGv7WDcR6Oj1+5wOtAGjw/nlwBcT1v8FsC7df9N6HPpQj16SOQ9YQ9DT/Q9gEXAOQY/t88A/mllJwvafA74DjAZWA48mrHsVmAmMCl/rF2ZWEK77b8CNBMlxBHAr0OzuHw3Xz/Cg9PKzJDHeHD7mAicR9CT/0d1bPehddrc/OUnbS4Ft7r6yn5/DbuCqMLZbgH8ws7MT1o8FCgi+DXwb+FeCn88s4BLg22Z2UsL2VwP/BxgJvA48R/CBNQH4G+BfBrHv/lxF8DubAVwHXNHHdn8A/qeZ3WJmUxNXuHstsIygB9/tz4El7r639wuZWQXwaWDTIOKUYyHdnzR6pP/B4T36jQnrziDoIVYmLNsHzAynHwYWJawrIegtT+pjX3UECRhgA3BNH9s5cMoRYl4K/EXC/DSgHYj3157gm8Afei2rJfgG00JCD7bXNr8E7gqn5wAHgVg4Xxru87yE7VcR9pYJevQvJKy7muCbR+/25f3tO8m6mzm8R39xwvzPgXv6aFtI8K1qVfjz2wRcmbD+88CGcDoHeA/4VML65UAzUB/udzUwOd1/03oc+lCPXpLZlTB9EMDdey9L7NFv655w90aCssV4ADP7ipmtN7N6MztAUF4ZHW4+iQ/rwYM1nqBs020rEAcqB9B2H3DIAUN3nxjGlQ9YGPuVZvYHM9sfxr4gIXaAfe7eGU4fDJ+P9HPqvW5vkvYlA9x3f3YmTDf3iqOHux909++6+yyCb3A/J/jWNSrc5ElgnJmdT/DhVgT8qtfL3OnuZcCZBN9WJg4iTjkGlOglFSZ1T4QlnVHAdjO7BPg6QelgpLuXE/T8LNx8G5CstDIQ2wkOuHabTFBP3pV880O8CEw0s9l9bWBm+cATBAcWK8PYl/Bh7MMmXft29w+A7wLFwInhsmbgcYKSzRcIvr219dF+LfC3wP1mNuw/Jxk4JXpJhQVmdrGZ5RHU6v/o7tsIyhEdwB4gbmbfJqg5d/sJ8B0zmxqe7nhmWOeFIGEn1rd7ewy428xODD9cvgv8zN07+gvW3TcQ1MMXmdk8Mys0sxjBQedueQS9+z1Ah5ldCQzpFMYhOGb7NrNvmdk5ZpYXHju5i6CEtSFhs0eA64E/49CzbZJ5BDiB4OwdyRBK9JIK/wHcS1CymUVwcBaCg42/JjgLZStB/XtbQrv7CEoFzxOclfNTgpoxBDXtR8zsgJldl2SfDxIc2FwBbA5f+45BxPxlglMs7wvjriX4kLoeeM/dG4A7w/jqgM8CTw/i9YfsGO/bgYeAvQTfkuYBnwhLcN1WEHwTe9/dXz3iiwW9/R8B3xqecGUoLDygIjIkZvYwUOvuf5XuWEQkOfXoRUSynBK9iEiWU+lGRCTLqUcvIpLl4ukOIJnRo0d7VVXVoNo0NTVRXFw8PAGlkOJMLcWZWooztY5lnKtWrdrr7mOSrkz3pbnJHrNmzfLBWrZs2aDbpIPiTC3FmVqKM7WOZZzAStctEEREokmJXkQkyynRi4hkuYw8GCsiciRmxubNm2lpaUl3KEdUVlbG+vXrU/qaBQUFTJw4kdzc3AG3UaIXkeNOcXExpaWlVFVVkck3ymxoaKC0tDRlr+fu7Nu3j9raWk488cQBt1PpRkSOO7FYjIqKioxO8sPBzKioqBj0NxklehE5LkUtyXcbyvvOqkT/o6UbeentPekOQ0Qko2RVov+Xl97hpQ1K9CIy/GKxGDNnzqS6upqrr76aAwcOALBlyxbMjG9968Nb8u/du5fc3Fxuv/12ADZs2MCcOXOYOXMm06dPZ+HChQAsX76csrIyZs6c2fP4zW9+c9SxZlWiLymI09janu4wRCQCCgsLWb16NTU1NYwaNYr777+/Z91JJ53E4sWLe+Z/8YtfcPrpp/fM33nnndx9992sXr2a9evXc8cdH46Zc8kll7B69eqex+WXX37UsWZXos+P09Ta2f+GIiIpdMEFF/D+++/3zBcWFjJ9+nRee+01AH72s59x3XUfDpS2Y8cOJk78cAz1M844Y1jjy6rTK0vy4zS09jtkqIhkkf/xzDre3P5BSl/ztPEjuPfq0/vfEOjs7GTp0qXcdttthyy/4YYbeOKJJzj55JOJxWKMHz+e7du3A3D33Xdz6aWXcuGFF/Lxj3+cW265hfLycgBefvllZs6c2fM63a9xNLKrR18Qp0mJXkSOgYMHDzJz5kwqKirYv38/8+bNO2T9/PnzWbZsGY899hjXX3/9IetuueUW1q9fz2c+8xmWL1/O+eefT2trK3B46eZokzxkWY++OC/O3obmdIchIsfQQHveqdZdo6+vr+eqq67i/vvv58477+xZn5eXx8yZM/n7v/971q1bxzPPPHNI+/Hjx3Prrbdy6623Ul1dTU1NzbDFmnU9+kb16EXkGCorK+NHP/oRP/jBD2hvP/RkkDvuuIO/+7u/o6Ki4pDlzz77bM+2O3fuZN++fUyYMGHYYsyqRF+ar0QvIsfeWWedxYwZM1i0aNEhy6dPn85NN9102PbPP/881dXVzJgxgyuuuILvf//7jB07FviwRt/9ePzxx486vuwq3YSJ3t0je9WciBwbjY2Nh8wnlmaSlWFuvvlmbr75ZgDuu+8+7rvvvsO2mTNnDvX19akNlCzr0ZcUxOnsclo7utIdiohIxsiuRJ8ffEFpaFH5RkSkW1Ymep1iKZL9gmFSo2co7zsrE70OyIpkt87OTvbt2xe5ZN99P/qCgoJBtcuqg7Eq3YhEQ1NTEw0NDezZk9k3MWxpaRl0Uu5P9whTg5Fdib5ApRuRKHD3QY2wlC7Lly/nrLPOSncY/ZduzOxBM9ttZjUJyz5jZuvMrMvMZvfRbpKZLTOz9eG2d6Uy8GSKVboRETnMQGr0DwPzey2rAT4NrDhCuw7gK+4+HTgf+LKZnTaUIAeqVIleROQw/ZZu3H2FmVX1WrYejjyklbvvAHaE0w1mth6YALx5FPEekXr0IiKHs4EctQ4T/WJ3r+61fDnwVXdfOYD2K4Bqd096P1EzWwgsBKisrJzV+1Li/jQ2NlJcXMytzzVz1Um5/NmpeYNqf6w0NjZSUlKS7jD6pThTS3GmluI83Ny5c1e5e9JSOu7e7wOoAmqSLF8OzO6nbQmwCvj0QPbl7syaNcsHa9myZe7uXn3vs37vUzWDbn+sdMeZ6RRnainO1FKchwNWeh85dVjPozezXOAJ4FF3f3I499WtRDc2ExE5xLAlegsK+D8F1rv74XfvGSbBcIJK9CIi3QZyeuVjwCvANDOrNbPbzOxTZlYLXAD8ysyeC7cdb2ZLwqYXAV8ALjWz1eFjwTC9jx7F6tGLiBxiIGfd3NjHqv9Msu12YEE4/VvgmN8ruFSDj4iIHCKr7nUDYY1et0AQEemRdYlepRsRkUNlXaLXWTciIofKykTfFA4nKCIi2ZjoC+J0ORxs70x3KCIiGSH7En33/W50QFZEBMjiRN+gOr2ICJDFiV5Xx4qIBLIu0RerdCMicoisS/SlBbonvYhIoqxL9CUafERE5BBZl+g1ypSIyKGyLtGrdCMicqisS/T58RxiOaaDsSIioaxL9GamwUdERBJkXaKH4ICsLpgSEQlkbaJXj15EJJCdiV6jTImI9MjKRF+sUaZERHoMZHDwB81st5nVJCz7jJmtM7MuM5t9hLbzzWyDmW0ys3tSFXR/SjX4iIhIj4H06B8G5vdaVgN8GljRVyMziwH3A1cCpwE3mtlpQwtzcIrzY0r0IiKhfhO9u68A9vdatt7dN/TT9Fxgk7u/6+5twCLgmiFHOggl+bk0tWrgERERgPgwvvYEYFvCfC1wXl8bm9lCYCFAZWUly5cvH9TOGhsbe9rs29lGY2sHLy5bRo7Z4KIeZolxZjLFmVqKM7UU5+AMZ6JPlmH7HMjV3R8AHgCYPXu2z5kzZ1A7W758Od1t3s55h6feeYtzLriY0oLcQb3OcEuMM5MpztRSnKmlOAdnOM+6qQUmJcxPBLYP4/56lOQHyV3lGxGR4U30rwJTzexEM8sDbgCeHsb99SjOjwHQ2Np+LHYnIpLRBnJ65WPAK8A0M6s1s9vM7FNmVgtcAPzKzJ4Ltx1vZksA3L0DuB14DlgP/Nzd1w3XG0n04R0s1aMXEem3Ru/uN/ax6j+TbLsdWJAwvwRYMuTohqi7dKOLpkREsvbKWJVuRES6ZWWiL+3u0at0IyKSnYm+p0ffoh69iEhWJvqS8GBsU5t69CIiWZno8+Mx8mI5NOhgrIhIdiZ6CMo3GnxERCSLE70GHxERCWRtoi/Oi6t0IyJCFif60gKNGysiAlmc6Es0ypSICJDFib44Xz16ERHI4kRfWhCnQYleRCR7E31xXlw3NRMRIYsTfUlBnIPtnXR29TmolYhIJGRvos/vvie9evUiEm1Zn+h1QFZEoi57E32BevQiIpDFib447NHr6lgRibqsTfSlKt2IiABZnOiLdTBWRAQYQKI3swfNbLeZ1SQsG2VmL5jZxvB5ZB9t7zazdWZWY2aPmVlBKoM/Ep11IyISGEiP/mFgfq9l9wBL3X0qsDScP4SZTQDuBGa7ezUQA244qmgHobT7YKxq9CIScf0mendfAezvtfga4JFw+hHg2j6ax4FCM4sDRcD2oYU5eMWq0YuIAGDu/V85amZVwOKwZ46ZHXD38oT1de5+WPnGzO4C/idwEHje3T93hH0sBBYCVFZWzlq0aNGg3khjYyMlJSWHLPt/n2/i8im5XD8tb1CvNZySxZmJFGdqKc7UUpyHmzt37ip3n510pbv3+wCqgJqE+QO91tclaTMSeBEYA+QCvwQ+P5D9zZo1ywdr2bJlhy07+2+e9288uWbQrzWcksWZiRRnainO1FKchwNWeh85dahn3ewys3EA4fPuJNtcDmx29z3u3g48CVw4xP0NSYkGHxERGXKifxq4KZy+CXgqyTbvAeebWZGZGXAZsH6I+xuSknzdwVJEZCCnVz4GvAJMM7NaM7sN+B4wz8w2AvPCecxsvJktAXD3PwKPA68Ba8N9PTAs76IPxRplSkSEeH8buPuNfay6LMm224EFCfP3AvcOObqjVJofZ+cHLenavYhIRsjaK2NBPXoREcjyRK+DsSIiWZ7oS/PjunuliEReVif64vw4rR1dtHd2pTsUEZG0yepEr1GmREQikuhVvhGRKMvuRB/ewbKpTYleRKIruxN9vm5VLCKS1Yleo0yJiGR5ou8ZfESJXkQiLKsTvQYfERHJ8kSvs25ERLI80RfnxQCVbkQk2rI60cdjORTmxlS6EZFIy+pED8G59OrRi0iUZX+iz4/T2NqZ7jBERNImGom+pT3dYYiIpE3WJ/ri/JhKNyISaVmf6Evyc1W6EZFIG8jg4A+a2W4zq0lYNsrMXjCzjeHzyD7alpvZ42b2lpmtN7MLUhn8QJQWxGlsVelGRKJrID36h4H5vZbdAyx196nA0nA+mR8Cz7r7R4AZwPohxjlkxfkxmtSjF5EI6zfRu/sKYH+vxdcAj4TTjwDX9m5nZiOAjwI/DV+nzd0PHEWsQ1KSn6u7V4pIpJm797+RWRWw2N2rw/kD7l6esL7O3Uf2ajMTeAB4k6A3vwq4y92b+tjHQmAhQGVl5axFixYN6o00NjZSUlJy2PKn32njyY3t/OvHi8jNsUG95nDoK85MozhTS3GmluI83Ny5c1e5++ykK9293wdQBdQkzB/otb4uSZvZQAdwXjj/Q+A7A9nfrFmzfLCWLVuWdPlDv33Xp3x9se9rbB30aw6HvuLMNIoztRRnainOwwErvY+cOtSzbnaZ2TiA8Hl3km1qgVp3/2M4/zhw9hD3N2QlBbmABh8RkegaaqJ/GrgpnL4JeKr3Bu6+E9hmZtPCRZcRlHGOqZJ83dhMRKJtIKdXPga8Akwzs1ozuw34HjDPzDYC88J5zGy8mS1JaH4H8KiZrQFmAt9Ncfz9KskPe/RK9CISUfH+NnD3G/tYdVmSbbcDCxLmVxPU6tOmOOzR6w6WIhJVWX9lbPdwgg1K9CISUVmf6HtKNzoYKyIRlfWJXqUbEYm67E/0eSrdiEi0ZX2iz8kxivM0nKCIRFfWJ3oIhxNUjV5EIioaiT5f48aKSHQp0YuIZLloJPoCJXoRia5IJPrivLgOxopIZEUi0ZcUxGnQwVgRiahoJHrV6EUkwiKT6JtaO7oHRBERiZRoJPqCOB1dTmtHV7pDERE55qKR6POD2yCofCMiURStRK8DsiISQZFI9MXq0YtIhEUi0Zcq0YtIhEUi0ZcUqHQjItEViUTfXbppalOiF5Ho6TfRm9mDZrbbzGoSlo0ysxfMbGP4PPII7WNm9rqZLU5V0IPVXbrR1bEiEkUD6dE/DMzvteweYKm7TwWWhvN9uQtYP6ToUqSnR68avYhEUL+J3t1XAPt7Lb4GeCScfgS4NllbM5sIfAL4ydBDPHpFeTHMdDBWRKLJBnJbADOrAha7e3U4f8DdyxPW17n7YeUbM3sc+F9AKfBVd7/qCPtYCCwEqKysnLVo0aJBvZHGxkZKSkr6XP+l3zRx8YQ4n5ueP6jXTbX+4swUijO1FGdqKc7DzZ07d5W7z0660t37fQBVQE3C/IFe6+uStLkK+HE4PYfgg2JA+5s1a5YP1rJly464/vzv/sa/+vPVg37dVOsvzkyhOFNLcaaW4jwcsNL7yKlDPetml5mNAwifdyfZ5iLgk2a2BVgEXGpm/z7E/R013cFSRKJqqIn+aeCmcPom4KneG7j7N9x9ortXATcAL7r754e4v6NWrEQvIhE1kNMrHwNeAaaZWa2Z3QZ8D5hnZhuBeeE8ZjbezJYMZ8BDVarhBEUkouL9beDuN/ax6rIk224HFiRZvhxYPsjYUqokP87O+pZ0hiAikhaRuDIWgtKNzqMXkSiKTKIvyY/ToEQvIhEUqUSv4QRFJIqik+gL4nQ5HGzvTHcoIiLHVHQSve5JLyIRFb1ErztYikjERC/Rq0cvIhETmUSvcWNFJKoik+hLNZygiERUZBK9hhMUkaiKTKLXwVgRiarIJPru0o2ujhWRqIlMos+P5xDLMd3vRkQiJzKJ3syCwUdUuhGRiIlMoofuUaZ0CwQRiZYIJvr2dIchInJMRSvRa5QpEYmgSCX6YpVuRCSCIpXoS/PjNLaodCMi0TKQwcEfNLPdZlaTsGyUmb1gZhvD55FJ2k0ys2Vmtt7M1pnZXakOfrCK82M0qUcvIhEzkB79w8D8XsvuAZa6+1RgaTjfWwfwFXefDpwPfNnMTjuKWI9aSX6uavQiEjn9Jnp3XwHs77X4GuCRcPoR4Nok7Xa4+2vhdAOwHphwNMEerZKCOE1tHXR1aThBEYkOG8gYqmZWBSx29+pw/oC7lyesr3P3w8o3vdqvAKrd/YM+tlkILASorKyctWjRooG/C6CxsZGSkpIjbvPrze38bEMb/3R5EYVxG9Trp8pA4swEijO1FGdqKc7DzZ07d5W7z0660t37fQBVQE3C/IFe6+uO0LYEWAV8eiD7cndmzZrlg7Vs2bJ+t3n0D1t9ytcX+44DBwf9+qkykDgzgeJMLcWZWorzcMBK7yOnDvWsm11mNg4gfN6dbCMzywWeAB519yeHuK+UKc6PARp8RESiZaiJ/mngpnD6JuCp3huYmQE/Bda7+31D3E9K9Qw+okQvIhEykNMrHwNeAaaZWa2Z3QZ8D5hnZhuBeeE8ZjbezJaETS8CvgBcamarw8eCYXkXA1SSnwugO1iKSKTE+9vA3W/sY9VlSbbdDiwIp38LpOeIZx+6SzcNuoOliERIxK6MDXr0Kt2ISJREKtF39+hVuhGRKIlUoi/RwVgRiaBIJfr8eIy8WI4SvYhESqQSPQTlGw0nKCJR0u9ZN9lGg49krpb2Trbua2bz3ia27GuiubWD6eNGUD2hjIkjCwkuzRCRwYpcoi/OU6JPF3entaOL2rpmNu9tZsveJjbva2LL3uCx44MWEm+9lGPQff+5kUW5VE8o44zw0Z38RaR/kUv0pQXxrCzdrN52gN0ftFCYF6MoL0ZhbpyicLogL0ZRbox4LLWVutaOTjbsbGBNbT1b9jbR1NZJU2sHzW0dNLV20tzWQWNrB809yzvp6HXn0LLCXE4cXcx5J1VQVVFM1egiThxdTNXoYvJiOWzY2cDa9+upeb+eNbX1PLDi3Z7XGFmUy/jCTv5w8C3Ki3LJMTAMMzAzcgxy7MN5I/j9nz15pL4hSKRELtEX58fZ39SW7jBSxt35xxc38fcvvN3vtnmxHArzYowojDNlVDFTKoqoqgifRxczeVQRBbmxpG3bOrp4e1eQdNfU1rP2/QNs2NlAe2eQdAtzY5QUxCnJDz5givPijCzOY+LIomA+P05xfoyivDjjygqoGl3MiRXFjCzOO2LMMyaVM2NSec98S3vnIcn/lQ3v85OX3z3sA6Q/J5TmM7tqJLOmjGLWlJGcPn4EuSn+IBRpbHOefmM7expaaWoNOj4NLR09040tHTS0fjhfnB/j5a9dmvI4IpfoS/LjvLevOd1hpER7Zxff/M+1/HxlLZ8+awK3XnwizW1BT7qlvTOc7uRg+Nzc3kFLWyf7m9t5b18Ti9fsoP7goUMrjisrYEpFEXltrazzTeyoP8ja2nrW72igrbMLgBEFcc6cWM5tF5/EmRODUsqx6iEX5MYOSf7Ll+/noks+SntnF+7Q5R6Ue8JpJ3z24ENxT2Mrr22tY+XWOlZuqWPJ2p3h6+YwY2I5s6tGMnvKKM6ePJKyotxhfz+SffY1tvL8m7tYsnYHv9/UTKe/3rOuIDeHkvxcSsNOUXF+jAnlhT3zFSVH7vgMVeQSfWmWHIxtaGnnLx59jZc37uXOS0/h7nmnDinRHmhuY+u+Zrbsazrk+c3dHayo3UBpfpzqCWXcclEV1RPKOHNiGZNHFWVU2SM3ljPg3vgJIwo4fXwZX7igCoCd9S2s2lrHyq37WbW1jn9+6V06u94B6Cl9FeUF31IKw28qwXOMwnD5iIJcxpUVMK68gHFlhYwvL6AoL3L/WpG2+4MWnlu3k1/X7OQP7+6jy2FKRRFXVOWycMG5nFhRTHF+6sunAxW5v8ZsOBi7o/4gtzz0Kpt2N/K//+xMrjtn0pBfq7woj/KivEPKIwDLly/n3AsvpiAeIycnc5J6qo0tK+ATZ47jE2eOA6C5rYPV2w7w+nsHqGtqo7m9k+bw+MLB9uBYw97G1oRvSx00tR0+DnFZYZD8x5cXMq6sgIP729g/opZpY0s5tbJUZaIssKP+IM/W7OTXa3fy6tb9uMPJY4r58txTuLJ6HNPHlfLSSy8xs9f/VjpELtGXFMRpbuuks8uJHYcJbN32em59+FWaWjt56JZzuGTqmGHbVxR7pUV5cS48eTQXnjx6wG1aOzrZVd/K9vqD7Kg/yPYDLeysb+mZfv29Ouqa23ly4xsA5Mdzer4dzZgYlKGqKjLrW5Ik19XlLH97Nw/+dgu/3bQXgI+MLeUvLzuVBWeMZWplaZojTC5y/8kl+cFbbmrrYERB/zXYN7d/QCzHmHpCSdp7ti+9vYe/+PdVjCjM5fEvXcBHxo5IazwSyI/HmFxRxOSKoj63eW7pMk4+4xzWbQ8OZr+x7QCP/ek9HvrdFuDD4x4zJpVx5sRyzp48kjGl+cfoHUh/mlo7eOK1Wh763RY2721i7IgC/tu8U/nEmeM4eUzmD2kY2UTf2NJ/ot+0u4Fr7v8t7Z1OeVEus6eMZHbVKM6pGsUZE8rIix+7r9+P/ek9/uqXNUyrLOWhW86hckTBMdu3HL38mHHKCSWcckIJ18ycAEBHZxcbdzeypvYAq7fVs6b2AP/yUnAGkRmcWzWKq2aMZ0H1WCpKlPTTobaumX97ZSuP/ek9Glo6mDGpnB/deBZXVo89rspv0Uv04Y3N+ruDpbtz79PrKMyN8TfXTGf1ewd4dct+frM+GDWxIDeHmZPKOSdM/GdPGdnzIZJKXV3OD57fwI+Xv8PHTh3D/Z87e1j2I8dePJbD9HEjmD5uBNefEyxrae/kzR0fsOLtPTzzxna+9csa/vrpdVx4cgVXnzmeK04fq7OBhpm7s3JrHQ/9bjPP1uzEzLiyeiy3XnwiZ08eme7whiRyGaM4TJIN/ST6X63dwe827eM715zOjedO5sZzJwOwp6GVVVv386fNdby6ZT/3L9tElwdXcVZPKOPiU0ZzydQxzJoy8qh7/PXN7XzrqRqefmM7N547me9cc3rajtrLsVGQG+PsySM5e/JI7rpsKm/tbGDxmu0888YOvvbEGr75y7V8dOoYrp4xnstPq9SHfgo1tnbwXM1OHnllC2tq6ykrzGXhR0/mzy+Ywvjy4/sq7Mj9lZQmlG760tTawd8uXs/p40fw2fOmHLJuTGk+86vHMb86OEujsbWD19+r49XN+3nl3X08sOJdfrz8HYryYpx/UgWXTA0S/8ljio94sK2to4u3dn7A6m0Heh7v7mkC4Gvzp/Glj52sg3URY2Y9Pf6vfnwaa2rrWbxmO4vX7GDpW7vJj+cwd9oJXHRKBdUTypg+bkSfF7xJch+0tPObN3exZO1OVmzcQ1tHFyeNKeZvr63m02dPyJoTErLjXQxCd4/+SKWbH724kZ0ftHD/587u98yckvw4l0wd03P2S0NLO394dz8vb9zDyxv38uJbQalnfFkBl0wdQ0V7BzOa2mho6eD1bXU9SX3d9g9o6wguSBpdks/MSeV8+qwJXHTKaM46Tr8uSuqYWc+FYt+4cjqr3qtj8RvbWVKzk2fXBRd9xXOMUytLg/sBhReyfWRcKflxJf9EB5rbeOHNXfy6Zicvb9xDe6czrqyAz503mSurxzF7ysi0n3iRapFL9CX9lG427W7gpy9v5rrZE5k1ZfAJtrQgl3mnVTLvtEoAtu1v5uWNe3l54x6W1OygoaWDH7/xQs/2Bbk5nDGhjJsumMKMSeXMnFTOhHLdh0X6lpNjPceG/vqTp7O9voW1tQd6bk/x3Js7+dnKbQDkxoLkf+bEMmZOKufiqWOYcJyXIYZif1Mbz6/byZKanfx+0146upwJ5YXcfGEVV54xjpkTy7MuuSfqN9Gb2YPAVcBud68Ol40CfgZUAVuA69y9Lknb+cAPgRjwE3f/XsoiH6LSIxyM7T4AW5QX4+vzP5KS/U0aVcRnz5vMZ8+bTEdnFw8/s4zWsimUFeYyc1I508bq4hkZOjNjQnkhE8oLe8qJ7k5t3UHWvl8fPGrrWbJ2J4/9KUj+J48p5pKpY/joqaM578SKnm+52eadPY28uH43S9/axatb6ujscqZUFPHFS05iwRljOWNCWWQ6VAP5DT8M/CPwbwnL7gGWuvv3zOyecP7riY3MLAbcD8wDaoFXzexpd38zFYEPVfERavSJB2CH43S2eCyHU8pjzJlzSspfW6SbmTFpVBGTRhWx4IwPk//G3Y2seDsoKS569T0e/v0WcmPGrCkjuWTqGD526hi6fHA3h8skbR1drNyyn6Vv7ebFt3azeW9wjOsjY0v50sdO5sozxnLauBGRSe6J+k307r7CzKp6Lb4GmBNOPwIsp1eiB84FNrn7uwBmtihsl9ZEnxvLIT9++HCCRzoAK3K8MwtKOKdWlvLFS06ipb2TlVvqeHnjHlZs3Mv3n9vA95/bQGkuzNn5OqeNG8HkUUVMqQg+MMoKM/OUzv1NbSzfsJul63ez4u09NLR2kBfP4cKTK7j1oiounV4ZyVJVb+YD+AQPE/3ihNLNAXcvT1hf5+4je7X5f4D57v7FcP4LwHnufnsf+1gILASorKyctWjRokG9kcbGRkpKBnaF2h0vNjG7Ms5Np3/Ya//5hjaWbG7nr84v4JTy4Tt4NZg400lxplamx3mgtYs393WxekcLb3+Qw4HWQ/NCcS6cUJTDmEILnouMEwpzOKHIGFlg5ByDXrK7s/eg8259F2/taWFrU4zN9V04UJZvzBgTY+aYGKdXxMiPZ0av/Vj+3ufOnbvK3WcnWzecxblkP+k+P1Xc/QHgAYDZs2f7nDlzBrWz5cuXM9A2o15dxoiKcubMOQsIDsA+//zLXDd7Il+8dsag9jtYg4kznRRnah0PcV7Lh3E2tLTz3v5mtu1vZuu+Zt7bHzy27mtm1e6DdCbc/z+eY4wvL2TiyEImjSwKnkd9+DymJH9IBzo/aGlnzbZ6Viecnba3MRhLIjfHOHNSGVfPHs1l00+genxZRh5MzZTf+1AT/S4zG+fuO8xsHLA7yTa1QOJtFScC24e4v5QqyY/3HIx1d779VGoPwIoc70oLcjl9fBmnjy87bF1HZxfbD7T0JP/auma21R2ktq6ZFzfsZk9D6yHb58VyGF9eQElBnLxYDnnxHPLiMfLCMmpePCdheQ71B9tZve0A7+xp7Bla8qQxxXz01DGcNamcmZNGsuvt17j80guPxY8iKww10T8N3AR8L3x+Ksk2rwJTzexE4H3gBuCzQ9xfShXnx2kID8b+au0Ofv/OPr5zbbXuJyIyAPFYzhFv4tbS3klt3UG21TVTG34AvF93kOa2Tto6umjr6KL+YHs43Ul7pwfTncG6gtwYMyaWcc2M8cG1AxPLD7vtw75Nmdd7z2QDOb3yMYIDr6PNrBa4lyDB/9zMbgPeAz4Tbjue4DTKBe7eYWa3A88RnF75oLuvG563MTil+XF2NbTQ2NrBdxa/SfWEEXw2vMWBiBydgtxYzw3cJDMM5KybG/tYdVmSbbcDCxLmlwBLhhzdMCnOj9O4p4P/f+lGdn3Qyj99ftZxeW96EZGBiOSVOiUFcXbUt/DT327m+tmTjts70omIDEQ0E31+nNaOLorz43xt/rR0hyMiMqwim+gBvnrFNB2AFZGsl503uejHgjPG0tbRpQOwIhIJkUz0p5xQylevUMlGRKIhkqUbEZEoUaIXEclySvQiIllOiV5EJMsp0YuIZDklehGRLKdELyKS5ZToRUSy3ICGEjzWzGwPsHWQzUYDe4chnFRTnKmlOFNLcabWsYxziruPSbYiIxP9UJjZyr7GS8wkijO1FGdqKc7UypQ4VboREclySvQiIlkumxL9A+kOYIAUZ2opztRSnKmVEXFmTY1eRESSy6YevYiIJKFELyKS5Y77RG9m881sg5ltMrN70h1PNzObZGbLzGy9ma0zs7vC5aPM7AUz2xg+Z8TI5GYWM7PXzWxxOJ9xcZpZuZk9bmZvhT/XCzI0zrvD33mNmT1mZgWZEKeZPWhmu82sJmFZn3GZ2TfC/6sNZnZFmuP8fvh7X2Nm/2lm5ZkYZ8K6r5qZm9nodMcJx3miN7MYcD9wJXAacKOZnZbeqHp0AF9x9+nA+cCXw9juAZa6+1RgaTifCe4C1ifMZ2KcPwSedfePADMI4s2oOM1sAnAnMNvdq4EYcAOZEefDwPxey5LGFf6t3gCcHrb5cfj/lq44XwCq3f1M4G3gGxkaJ2Y2CZgHvJewLJ1xHt+JHjgX2OTu77p7G7AIuCbNMQHg7jvc/bVwuoEgKU0giO+RcLNHgGvTEmACM5sIfAL4ScLijIrTzEYAHwV+CuDube5+gAyLMxQHCs0sDhQB28mAON19BbC/1+K+4roGWOTure6+GdhE8P+Wljjd/Xl37whn/wBMzMQ4Q/8AfA1IPNMlbXHC8Z/oJwDbEuZrw2UZxcyqgLOAPwKV7r4Dgg8D4IQ0htbt/yP4w+xKWJZpcZ4E7AEeCktMPzGzYjIsTnd/H/gBQW9uB1Dv7s+TYXEm6CuuTP7fuhX4dTidUXGa2SeB9939jV6r0hrn8Z7oLcmyjDpf1MxKgCeAv3T3D9IdT29mdhWw291XpTuWfsSBs4F/cvezgCYyo5x0iLDGfQ1wIjAeKDazz6c3qiHJyP8tM/smQVn00e5FSTZLS5xmVgR8E/h2stVJlh2zOI/3RF8LTEqYn0jwNTkjmFkuQZJ/1N2fDBfvMrNx4fpxwO50xRe6CPikmW0hKH1damb/TubFWQvUuvsfw/nHCRJ/psV5ObDZ3fe4ezvwJHAhmRdnt77iyrj/LTO7CbgK+Jx/eAFQJsV5MsEH/Bvh/9NE4DUzG0ua4zzeE/2rwFQzO9HM8ggOdjyd5pgAMDMjqCevd/f7ElY9DdwUTt8EPHWsY0vk7t9w94nuXkXw83vR3T9P5sW5E9hmZtPCRZcBb5JhcRKUbM43s6Lwb+AyguMzmRZnt77iehq4wczyzexEYCrwpzTEBwRn1wFfBz7p7s0JqzImTndf6+4nuHtV+P9UC5wd/u2mN053P64fwAKCo/DvAN9MdzwJcV1M8NVsDbA6fCwAKgjObtgYPo9Kd6wJMc8BFofTGRcnMBNYGf5MfwmMzNA4/wfwFlAD/B8gPxPiBB4jOG7QTpCEbjtSXARliHeADcCVaY5zE0GNu/t/6Z8zMc5e67cAo9Mdp7vrFggiItnueC/diIhIP5ToRUSynBK9iEiWU6IXEclySvQiIllOiV4iwcwqzew/zOxdM1tlZq+Y2afSHZfIsaBEL1kvvHDpl8AKdz/J3WcRXBw28YgNRbKEzqOXrGdmlwHfdvePJVlXRXBRU3G46HZ3/72ZzSG48GkXwYVaTwJrCW7nXAhc6+7vmNnDwEHgI8AU4BaCK0wvAP7o7jeH+/kn4Jyw7ePufm/q36lIcvF0ByByDJwOvNbHut3APHdvMbOpBFc7zg7XzQCmE9yK9l3gJ+5+rgWDyNwB/GW43UjgUuCTwDME9w/6IvCqmc1099UEV23vD+9BvtTMznT3NSl+nyJJqXQjkWNm95vZG2b2KpAL/KuZrQV+QTCATbdXPRhXoJXg0vXnw+VrgaqE7Z7x4KvxWmCXB/c86QLWJWx3nZm9BrxO8MGTKQPkSASoRy9RsA74s+4Zd/9yOMTbSuBugvLMDIKOT0tCu9aE6a6E+S4O/d9pTbJNz3bhTay+Cpzj7nVhuafgKN+TyICpRy9R8CJQYGZfSlhWFD6XATvCHvgXCIb+S7URBPfPrzezSoKhL0WOGfXoJeu5u5vZtcA/mNnXCEaqaiK47e1rwBNm9hlgWbg81ft/w8xeJ/hm8S7wu1TvQ+RIdNaNiEiWU+lGRCTLKdGLiGQ5JXoRkSynRC8ikuWU6EVEspwSvYhIllOiFxHJcv8X0IRdBEKu1t4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "gamma_vals = ee.List.sequence(1, 150, 5)\n",
    "\n",
    "\n",
    "def gammas(t):\n",
    "    svrRegression = ee.Classifier.libsvm(svmType ='EPSILON_SVR',kernelType='RBF',gamma=t,cost=100).setOutputMode('REGRESSION').train(\n",
    "      features= training, \n",
    "      classProperty= 'Percent_Tree_Cover', \n",
    "      inputProperties= predictionBands\n",
    "    )\n",
    "    \n",
    "    svrRegressionTesting = testing.classify(svrRegression, 'y_hat');\n",
    "\n",
    "    # // Calculate the RMSE for SVM regression and print the result\n",
    "    regression=svrRegressionTesting.map(set_sqError)\n",
    "    regression = ee.FeatureCollection(regression)\n",
    "    MSE = regression.reduceColumns(ee.Reducer.mean(),['sqError'])\n",
    "    RMSE =  ee.Number(MSE.get('mean')).sqrt()        \n",
    "    return RMSE\n",
    "\n",
    "\n",
    "rmse_gama=gamma_vals.map(gammas)\n",
    "value_info = rmse_gama.getInfo()\n",
    "\n",
    "# Extracting gamma value for minimum RMSE\n",
    "min_rmse=min(value_info)\n",
    "print(\"Min RMSE =\", min_rmse)\n",
    "\n",
    "index=value_info.index(min(value_info))\n",
    "gamma_range = gamma_vals.getInfo()\n",
    "min_gamma =gamma_range[index]\n",
    "print(\"Min Gamma =\", min_gamma)\n",
    "\n",
    "# print(rmse_gama.getInfo())\n",
    "\n",
    "# Plotting the interations\n",
    "import pandas as pd\n",
    "df =pd.DataFrame(value_info,columns=['RMSE'])\n",
    "df['Gamma'] = gamma_vals.getInfo() \n",
    "\n",
    "ax =df.plot.line(x='Gamma', \n",
    "             y='RMSE',\n",
    "             title= 'Impact of Gamma in SVR'\n",
    "             )\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "same for cart regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Min RMSE = 10.564642576458905\n",
      "Min Leaf = 21\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXoAAAEWCAYAAABollyxAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5e0lEQVR4nO3deXiU5bn48e+dSSY7CSQkkAQI+w5BEBVEcaniUuni2lrrVrqpPbY91VbPac/p/rNqteVYrVpaF1Db2ipaRS0RBDeQfV8SkhAIhJB9nZnn98e8wRAmJJPMzDuZuT/XNVdm3vWe7Z4nz/ssYoxBKaVU5IqxOwCllFLBpYleKaUinCZ6pZSKcJrolVIqwmmiV0qpCKeJXimlIpwmehVWROTzIlIqIvUiMiNE58wWkVUiUiciD4binKeJZYmI/MyG835ZRFaE+rwqNDTR92MiUiwiF9sdR0ciYkRkTB8O8RvgDmNMijFmQxCO78sioBIYYIz5XoCP3S8YY54zxlzS2/1FZJyIvCQilSJSIyKbReS7IuLo5fFs+cGLVJroVbgZAWyz4ZzbTYT1Huxtku3FeUYDHwKlwFRjTBpwDTALSO3F8UISd1Qxxuitn96AYuBi6/7NwBrgYaAa2A/MsZaXAkeAr3bYdwnwB+AtoA54FxjRYf0j1n61wHpgXod1DuBHwD5r3/XAMGAVYIAGoB64zkfMMcD9wAErpr8AaUC8tU/7/vu6eM4GGONjeTze/wZKgArruSVa6wYCy4GjwHHrfl6H16ENaLXOf7GPYy8BFgOvWc/3Q2C0tS7fiim2w/aFwO1Bel8mWOuqgF3AtZ32fQx43XoNLwYuB7ZbxzoIfL+L1/Vm4L1Or/M3gD3Wa7YYkC72fRZ4rZvP6kvAYaDG+pxMPk3cizq9J6/a/V3r7zfbA9BbH968UxO9C7gFbyL+mZX0FltJ8BLry55ibb/Eenyetf6RTl/0G4EMIBb4nvUlTbDW/SewBRgPCDAdyLDW+UzEHY57K7AXGAWkAH8Hnumwvrv9u0r0vwVeAQbhLUW+CvzSWpcBfBFIsta9BPyjw75LgJ+d5pxLrMQ623o9ngOWWevy6T7RB+R9AZLx/jjcYsVxBt4qp8kd9q0B5uL9QU0ADmH9SOP9wTuji+d4M6cm+uVAOjAc74/kgi72PQzc0s1n9VbrtY+33quNnV7fznGf9j3Rm5+5wu4A9NaHN+/URL+nw7qp1pc1u8OyY0CBdX9Je7KyHqcAbmBYF+c6Dky37u8CFnaxXXeJ+h3gWx0ej8dbeovt4f6nrMf7Y9OAVcq2lp0DFHVxjALgeIfHp00q1vonOzy+HNhp3c+n+0QfkPcFuA5Y3Sm2x4Efd9j3L53WlwBfx3v94XSfpZs5NdGf2+Hxi8C9XezbRhc/Al1sn24dP+00cZ/2PdGbfzeto48sFR3uNwEYYzovS+nwuLT9jjGmHm+pNQdARL4nIjusC2vVeKtXMq3Nh+GttumNHLzVNu0O4C2dZvfyeACD8ZbW14tItRXvG9ZyRCRJRB4XkQMiUou36iDdz7rgwx3uN3Ly69idQL0vI4Cz2p+j9Ty/DAzxta/li3h/mA6IyLsico4fcff0OR8DhnZ1EBFxiMivRGSf9foXW6syO2zWOW4VQJroo9uw9jsikoK32qNcROYB9wDXAgONMel4/7UWa/NSYHQvz1mON2G1G463aqPC9+Y9Uok3WU42xqRbtzRjTHti+h7e/xzOMsYMwFstAp8+n75osP4mdVg2xNeGfvD5vuB93d/t8BzTjbd10jc77HvSBWVjzMfGmIVAFvAPvCXzQHsb7w9KV74ELMR7zSAN739BcPLr3/lCeERdGLebJvrodrmInCsiTuCnwIfGmFK8dakuvPWysSLy38CADvs9CfxURMaK1zQRybDWVeCtf+/KUuBuERlpJbFfAC8YY1x+xO0UkYT2G96E8UfgYRHJAhCRXBG51No+Fe8PQbWIDAJ+7Me5TssYcxTvRc4brZLrrfT+R7BdV+/LcmCciHxFROKs25kiMtHXQUTEabWPTzPGtOG9sO7uY2y+/BiYIyIPiMgQ69xjRORZEUnH+/q34C35J+F9z7vT3edI+UETfXR7Hu+XtAqYibcaAOBN4F/AbrxVK82c/K/1Q3hLhivwJo+ngERr3U+AP1tVC9f6OOfTwDN4q0+KrGPf6Wfc2/Am7vbbLXj/A9kLfGBVD7yNtxQP3ot/iXhL/h/grdYJpK/hvUB9DJgMrO3j8Xy+L8aYOrwXb6/HW8I/DPwa7wXOrnwFKLZek2/gvcgeUMaYfXivieQD20SkBvgbsA7vheW/4P0cHcTbAuiDHhz2KWCS9Tn6R6BjjjZiXfhQUUZElgBlxpj77Y5FKRVcWqJXSqkIp4leKaUinFbdKKVUhNMSvVJKRbhYuwPwJTMz0+Tn5/u1T0NDA8nJycEJKIA0zsDSOANL4wysUMa5fv36SmPMYJ8r7e6a6+s2c+ZM46+VK1f6vY8dNM7A0jgDS+MMrFDGCawzOgSCUkpFJ030SikV4TTRK6VUhNNEr5RSEU4TvVJKRThN9EopFeE00SulVISLmERvjOF37+zh3d1H7Q5FKaXCSsQkehHhidX7WbnziN2hKKVUWImYRA8wOCWeo/UtdoehlFJhJaISfUaKk2Oa6JVS6iQRlegzU+KprG+1OwyllAorEZjotUSvlFIdRVyir25so83tsTsUpZQKG5GV6FOdABzT6hullDohshJ9SjyAVt8opVQHEZnotYmlUkp9KqIS/eD2En2dJnqllGoXUYm+vY5em1gqpfxljOGNrYcisi9ORCX6JGcsSU6H1tErpfz2pzXFfOPZT/jSHz+kprHN7nACKqISPWhbeqWU/9YfqOIXr+9gxvB0iiobuGXJRzS2uuwOK2AiMNE7NdErpXqssr6Fbz33CbkDE1lyy2weub6AjaXVfP2Z9bS6IqNPTgQm+ngq67SOXinVPbfHcNfSDVQ3tvHYl2eSlhjHZVOH8ssvTGX1nkrufnEjbo8JSSzGGJpa3UE5dmxQjmqjzNR41h84bncYSql+4KG3drF23zEeuHoak3IGnFh+3ZnDqW5s45f/2smAhDh+8fkpiIjfx29xG/68tphDNc3Ut7TR0OKmvsVFg3Wrt24NLW4aWl1kpcbz4Y8uDuRTBCIx0afEU9XYisvtIdYRcf+wKKUC5J0dFSxeuY/rzxzGNbOGnbL+6+ePprqpjccK9zEwKY4fLJjQ42MbY1ixvYIfrW7iWPM2nI4YkuMdpCTEkuyMJSU+lvQkJ3kDk0iJjyU5PpaUeAcDk52BfIonRFyiH5zixBioamwlKzXB7nCUUmGotKqRu1/YyOScAfzkqsldbveDS8dT3djG/xXuIz0pjkXnje722AeONfCTV7axctdRclOEFxadzVmjMgIZvt8iLtGfGAahThO9UupUzW1uvvncegAe+/JMEuIcXW4rIvzsc1OobW7jF6/vJC0xjuvOHN7lcR8r3Mdj7+4jLka4/4qJ5LcdsD3JQyQm+lQd70Yp1bX/eXUbWw/W8uRNsxiekdTt9o4Y4eFrC6htauOHf99CWmIcC6YMPWmblTuP8ONXtlFS1ciV04Zy/xWTGJKWQGFhSbCehl8irhJbBzZTSnXlpXWlLP2olG/NH83Fk7J7vJ8zNobHvzKTgmHp3LV0I+/tqQSg7Hgji/6yjluWfEycQ3ju9rP4/ZfOYEhaeNUmRF6JPqV9GARN9EqpT20vr+X+f2zlnFEZfPcz4/zeP8kZy59uns11T7zPomfW8aXZw3n2wwMIwg8WjOf2c0fhjA3PsnN4RtUHKfGxxMfG6Hg3SqkTapra+OZz60lPiuPRG2b0ukVeWlIcf7l1Npkp8Tz5XhHzx2Xx9vfO51vzx4RtkocelOhF5GngSuCIMWaKteynwELAAxwBbjbGlPvYdwHwCOAAnjTG/CqAsXcVr9VpSkv0SilvU8f/fGkTB483sWzR2Qy2ruP1VtaABP72zTmUVDUyc8TAAEUZXD35CVoCLOi07AFjzDRjTAGwHPjvzjuJiANYDFwGTAJuEJFJfYq2hzJT43VMeqUUm0qruXPpBlZsr+DeyyYwK39QQI47ODW+3yR56EGJ3hizSkTyOy2r7fAwGfDVR3g2sNcYsx9ARJbh/S9ge6+j7aHBKU4OVjcH+zRKqTDU1Orm1U3lPPPBAbYcrCHZ6eCb80dz27kj7Q7NNmJM9+M4WIl+eXvVjbXs58BNQA1wgTHmaKd9rgYWGGNutx5/BTjLGHNHF+dYBCwCyM7Onrls2TK/nkh9fT0pKSkAPL21hU1H3TxyQfdNp0KtY5zhTOMMLI0zsHzFebjBw79L2njvoItGF+SmCBcOj2NOTiyJsf4PXxCsOIPlggsuWG+MmeVrXa9b3Rhj7gPuE5EfAncAP+60ia9XtstfFWPME8ATALNmzTLz58/3K57CwkLa9/m4ZSdryvdz3nnnExNjzxvclY5xhjONM7A0zsBqj9Pl9vD2jgqe/aCE9/ZWEucQFkzJ4cazhjN75KBejU8TjDjtFojmlc8Dr3Fqoi8DOg4gkQeccsE2GDJT4nF7DNVNbQwK0tgRSin7GGP446r9PPVeEYdrm8lJS+D7l4zj2jOHaY94H3qV6EVkrDFmj/XwKmCnj80+BsaKyEjgIHA98KVeRemnjp2mNNErFXl2Vnn49cc7OHvUIP534WQunJClgxieRk+aVy4F5gOZIlKGt+R+uYiMx9u88gDwDWvbHLzNKC83xrhE5A7gTbzNK582xmwLztM4WWaHScLHZaeG4pRKqRB6s7iNjGQnS26ZfdqxapRXT1rd3OBj8VNdbFsOXN7h8evA672OrpcGW5OEaxNLpSLP/qP1bDzq5jsXjdIk30MR+b/Op1U32jtWqUjz9JoiYmPgxrNH2B1KvxGRiT4tMY44h+h4N0pFmOrGVv66voxzhsb2uYdrNInIRC8iZCTrMAhKRZrnPiyhuc3DpflxdofSr0RkogfITHVqiV6pCNLq8vDntcXMG5tJXmrEpq6giNhXKzMlXuvolYogr20p50hdS1QPZdBbEZ7otUSvVCQwxvDk6iLGZqVw/rjBdofT70R0oj9W30pPxvJRSoW3D4uq2FZey63njrR9WIP+KIITvZNWt4faJpfdoSil+ujJ1UUMSnby+Rm5dofSL0Vsom9veqWdppTq34oqG3hnZwU3njVcO0j1UsQmep0kXKnI8Kc1RcTFxHDjOdpBqrc00SulwlZ1YysvrSvjqoIcHZWyDyI40XvHu9FOU0r1X89/VEJTm1ubVPZRxCb6gUlOHDGibemVslmLy82v39jJB/uP+bVfewepuWMymDh0QJCiiw4Rm+hjYoRBydo7Vik7GWP48T+38VjhPr70xw/47du7cXt61uT59S2HqKht4fZzRwU5ysgXiBmmwpZ2mlLKXs99WMKyj0v52ryRHKtv5bdv7+HD/VU8cn0BWQO6rnM3xvDke/sZPThZO0gFQMSW6MFbT39Uq26UssVHRVX85JVtXDghi3svm8hD1xXwwNXT2FhazWWPrGbV7qOn3XfrQW8HqXCb97k/iuhEPzhFR7BUyg7l1U1867n1DB+UxMPXFeCwkvU1s4bxyh1zyUhx8tU/fcQDb+7E5facsv+T7xUxMCmOL8zIC3XoESmiE31mqrfqRodBUCp0mtvcfP2Z9TS3eXjippmkJZ48pPDY7FT++e1zuW7WMBav3McNf/yAQzVNJ9YXVzbw9o4KvnzWCBKd2kEqECI70ac4aXF5qG/RYRCU6okml+HK363mf1/dTnOb2+/9jTH86OUtbDlYw8PXFTAmy/eczYlOB7/64jQeub6A7eW1XP7Iav69swLwdpCKjRFu0g5SARPhiV6nFFTKH6vKXGw9WMvTa4r47O/eY+vBGr/2/9OaYv7+yUHuvngcn5mU3e32CwtyefXOcxmSlsitS9bxk1e28eK6Mj47Pee0F2uVfyI60WdYif6YtrxRqlsut4cVxW3Mzh/En2+dTU1TG5//vzX8X+HeHjWJXLu3kp+/voNLJ2dz54VjenzeUYNTePlbc/jK2SNYsrZYO0gFQYQ3r7R6x2qiV6pbb26r4Fiz4RfzRnL+uMG8+R/ncd8/tvD/3thF4c6jPHjtdIYNSvK5b2lVI99+/hNGZSbz4LUFfreUSYhz8NPPTeHcsZmUHW9ick5aIJ6SskR0iX5wSvsIllp1o1R3nnxvP1lJwsUTvVUuA5OdLP7SGTx4zXS2H6rlskdW89f1Zac0bmhsdbHomfW4PYY/3jSLlPjelx8vnTxES/NBENGJflCyExEd70ap7qw/UMWGkmouGRF3oikkgIjwxZl5/Os785g0dADff2kT33ruE6oavIUnYww/+Otmdh6u5dEbZpCfmWzXU1CnEdFVN7GOGAYm6TAISnXnydVFpCXGMS/Xd0oYNiiJpYvO5o+r9/Pgil2sO3CcB66exs7DdSzffIh7Fkxg/visEEeteiqiEz146+k10SvVtZJjjby57TBfP3808bGHu9zOESN84/zRzBubyd0vbOTmP32MCFwxbSjfOF/Howln3VbdiMjTInJERLZ2WPaAiOwUkc0i8rKIpHexb7GIbBGRjSKyLoBx95h3vButo1eqK0+vKSJGhK+ek9+j7SfnpPHKHefytXkjuWB8Fg9cPU3ncQ1zPamjXwIs6LTsLWCKMWYasBv44Wn2v8AYU2CMmdW7EPtGBzZTqms1TW28uK6Uq6bnMCSt5+3WE+Ic3HfFJJ6++UySnBFfMdDvdZvojTGrgKpOy1YYY9q7m34AhO2AFJk63o1SXVr2UQmNrW5u1ZYuEU16Mg6MiOQDy40xU3ysexV4wRjzrI91RcBxwACPG2OeOM05FgGLALKzs2cuW7asp88BgPr6elJSUk5Zvnx/K3/d3cbjn0ki3mH/v5ddxRluNM7ACsc4XR7DD1Y1kZ0k3DM7EQjPOH3ROE91wQUXrO+y5sQY0+0NyAe2+lh+H/Ay1g+Gj/U51t8sYBNwXk/ON3PmTOOvlStX+lz+wsclZsQ9y03JsQa/jxkMXcUZbjTOwApWnCXHGkxLm7tX+/5jQ5kZcc9y886OwyeWRfvrGWihjBNYZ7rIqb1uRy8iXwWuBL5sncTXj0i59feI9YMwu7fn661PO01p9Y2KLDsO1TL/N4XcsuQjvwcgM8bwx9X7GTU4mfnjtFlkpOtVoheRBcA9wFXGmMYutkkWkdT2+8AlwFZf2wbTiYHNtJ5eRRBjDD9dvp342BjW7jvGN55dT4ur58n+Q2tij9t0Yo+o0JPmlUuB94HxIlImIrcBvwdSgbesppN/sLbNEZHXrV2zgfdEZBPwEfCaMeaNoDyL08hMbR/vRptYqsixYnsFa/cd497LJvDLz0+lcNdR7nh+A20+JvHw5cnV3ok9vnhG2LajUAHUbbsoY8wNPhY/1cW25cDl1v39wPQ+RRcAGcntQxVriV5FhhaXm1+8voOxWSl8afZwYh0xtLo9/Pc/t/EfyzbyyPUFxDq6LsPtP1rPOzsruPOCMSTE6cQe0SDiG8A6Y2NIS4zTRK8ixpI1xRw41shfbp19IqHfdE4+rS4PP3ttB3EO4cFrC04as6ajp9cUERcTw406sUfUiPhEDzoMgoocR+ta+N2/93LRhCzOGzf4pHW3zxtFi8vDA2/uwhkbw6++MO2U+vfjDa38dX0ZCwtyyErViT2iRZQk+ngq67SOXvV/D721i+Y2Nz+6YqLP9d++YAwtLg+PvrMHZ2wMP1045aThCZ7/qITmNg+3z9OxaaJJdCT61Hh2lNfaHYZSfbKtvIZlH5dy69yRjB7cdSecuy8eS4vLzePv7sfpcPBfV05ERGhxuVmytph5YzMZP8T3XK4qMkVFoh+cEs8qrbpR/Vh7c8r0xDjuunDsabcVEe5dMIFWl4en1xThjI3hngXjeXXTIY7WtfCba2xvI6FCLCoSfWaKk7pmF81tbm1loPqlN7cd5oP9Vfz0c1NIS4rrdnsR4b+vnESry8Mf3t1HfGwMb247zLjsFM4bmxmCiFU4iZJEb00S3tBKbnqizdEo5Z8Wl5ufv76D8dmp3HDmsB7vJyL8dOEUWl0eHnlnDwD/74s6pHA0iqpEX1nXoole9TtPv1dMaVUTz9521mnbx/sSEyP86ovTANhQWs1VBTnBCFGFuehI9KnaaUr1T0fqmlm8ci8XT8zm3F5WuThihAeumY4xRkvzUSqiJwdvl5nSPgyCJnrVvzz45m5aXG7u66I5pT80yUevKEn07SV6bUuv+o+tB2t4cX0pN8/JZ2Rmst3hqH4sKhJ9QpyDlPhYLdGrfsMYw/8u387AJCd3dNOcUqnuREWih/ZhELREr/qHN7Ye5qOiKr53yTjSErtvTqnU6URRote5Y1X/0NzmbU45YUgq183qeXNKpboSFa1uwJvo9x2ttzsMFaWMMdz/jy0UVzbS5vbQ5vbg8hhaXd6/LreHNrehze2huc1NbbOL52/3vzmlUr5ET6JPdfJhkZbolT0ONxie/aCE0YOTyUiJJ8kZS6xDiHPEEOcQYmNiPr3vEKbmpjFnjPZgVYERPYk+JZ7jjW20uT3EaSlJhdj+Gu80f4/dOJNx2TqgmAqtqMl47U0sqxr0gqwKvX3VHlLiY0876qRSwRJ1if6oXpBVNthX42H6sLQuZ31SKpiiJtEPTtXescoeTa1uyuo8FAxLtzsUFaWiJtFr71hll63lNbgNFAwbaHcoKkpFYaLXEr0KrY0l1QBaole2iZpEnxwfS2KcQztNqZDbWFpNZqIw2BpFValQi5pED9629FqiV6G2oeQ4o9Ki6qumwkxUffoyU+K1jl6F1JHaZsprmhmdrlNYKvtEYaLXEr0KnQ2l1QCM1hK9slG3nz4ReVpEjojI1g7LHhCRnSKyWUReFpH0LvZdICK7RGSviNwbwLh7RRO9CrUNJdXEOYThAzTRK/v05NO3BFjQadlbwBRjzDRgN/DDzjuJiANYDFwGTAJuEJFJfYq2jwanOKlqaMXtMXaGoaLIxtLjTBw6AKdDO0op+3Sb6I0xq4CqTstWGGNc1sMPgDwfu84G9hpj9htjWoFlwMI+xtsnmanxeIwOg6BCw+0xbCmrYYY2q1Q2C8SgZrcCL/hYnguUdnhcBpzV1UFEZBGwCCA7O5vCwkK/gqivr+92n4rD3t+mNwrXMCzVnn+lexJnONA4+660zkNDq5v4hsPUO5rDNs6Owvn17Ejj9E+fEr2I3Ae4gOd8rfaxrMs6E2PME8ATALNmzTLz58/3K5bCwkK62yepqIrFG98nf8JU5o0d7NfxA6UncYYDjbPvln5UAmzhhkvO4cDWj8M2zo7C+fXsSOP0T6+LtSLyVeBK4MvGGF8JvAzoOD1OHlDe2/MFQmaKjnejQmdjSTXpSXHkZyTZHYqKcr1K9CKyALgHuMoY09jFZh8DY0VkpIg4geuBV3oXZmBkWj0TK+u0jl4F38bSagqGpSOiF2KVvXrSvHIp8D4wXkTKROQ24PdAKvCWiGwUkT9Y2+aIyOsA1sXaO4A3gR3Ai8aYbUF6Hj2SGh+LMzZGS/Qq6Oqa29h9pE7Ht1Fhods6emPMDT4WP9XFtuXA5R0evw683uvoAkxEGJwSz1FN9CrItpTVYIwOZKbCQ9T14shIcXJMh0FQQdbeI1YTvQoHUZfotXesCoUNJdWMykwmPclpdyhKRWOi1xEsVXAZY05ciFUqHERhoo/nWH0rHh0GQQXJweomKutbKBiebncoSgFRmuhdHkNNU5vdoagItdGqn5+hUweqMBF9iT5VpxRUwbWhpJr42BgmDE21OxSlgGhM9FbvWG1iqYJlY2k1U3LTiHNE3ddLhamo+yQOPjFJuDaxVIHX6vKw9aCOWKnCS9Ql+sz2RK+ThKsg2Hm4lhaXRy/EqrASdYk+LTGO2BjROnoVFBu1o5QKQ1GX6GNihAxtS6+CZGNJNZkp8eSmJ9odilInRF2ih/besVpHrwJvQ2k1M4briJUqvERxotcSvQqs6sZWiiobtNpGhZ3oTfR6MVYF2KcdpdJtjUOpzqIz0ac6qaxvxffEWEr1zsbSakRgmiZ6FWaiMtHnpCXS6vZQUaul+khQ29wWFmMXbSipZlxWKinxfZqKWamAi8pEPyV3AABbDtbYHInqi1aXh5f3tHLG/77F02uKbI3FGMOmMh2xUoWnqEz0k4am4YgRtpRV2x2K6qXt5bUsXLyGf+5rwxEjrNheYWs8xccaqW5s045SKixF5f+YiU4HY7NS2FSmJfr+ps3tYfHKvfz+33tJT3Jy14x4WlJzeXpNEY2tLpKc9nykN5QcB2CGJnoVhqKyRA8wLS+NLQdr9IJsP7K9vJaFv1/Db9/ewxXThvLW3edxRnYsc8Zk0uY2fFx83LbYNpZWk+x0MDZLR6xU4SeKE306VQ2tHKxusjsU1Y02t4dH3t7DVb9/jyN1zTz+lZk8cv0MBiZ7RyI9M38gTkcMa/ZW2hbjxtJqpuZ5qwSVCjdRWXUD3hI9wOayGvIGJtkcjerKjkO1fP+lTWwrr+Wq6Tn8z1WTTyT4dknOWGYMT7ct0Te3udleXsvXzhtly/mV6k7UlujHD0klziFs1nr6sNTm9vDoO95SfEVtM3+48QwevWHGKUm+3dwxmWw/VMvxhtAPbbGtvAaXx2iLGxW2ojbRx8c6mDh0AJu15U1Y+tOaIh56azcLpgxlxd3ns2DK0NNuP3dMJsbA+/uPhSjCT20oqQa0R6wKX1Gb6AGm5novyIZDZxt1sk8OVDMyM5nf3TCDQV2U4juanpdGSnws79lQfbOxtJqctASyBiSE/NxK9URUJ/ppeWnUNbs4UNVodyiqk91H6hif3fMWLLGOGM4aOYi1NiT6DSXVzBiuE4Gr8BXliT4dQKtvwkxzm5viygbGDfGvqeLcMZkUH2uk7HjofriP1rVwsLpJ6+dVWOs20YvI0yJyRES2dlh2jYhsExGPiMw6zb7FIrJFRDaKyLpABR0oY7NSiI+N0QuyYWbf0Xo8BsZlp/i139wxmQCs3Ru6evoTM0ppRykVxnpSol8CLOi0bCvwBWBVD/a/wBhTYIzp8gfBLrGOGCbnDGCLJvqwsqeiHsCvqhvw/jBkpsSzZl/oqm/e33eM2BhhSk5ayM6plL+6TfTGmFVAVadlO4wxu4IWVQhNy0tna3kNbr0gGzZ2VdQR5xDyM5P92k9EmDsmgzV7j4Wkx3NTq5u/fVLGpZOHkOh0BP18SvVWsDtMGWCFiBjgcWPME11tKCKLgEUA2dnZFBYW+nWi+vp6v/cBiKtro7HVzbLXVpKbGvxLFr2NM9TsjHPttmayE2HN6u7/YewcZ4arjcr6Vp5bvpK8IL+f75a2UdPUxrTE492+Vvq+B5bG6SdjTLc3IB/Y6mN5ITDrNPvlWH+zgE3AeT0538yZM42/Vq5c6fc+xhizp6LOjLhnuXlpXWmv9vdXb+MMNTvjPPfX75g7nv+kR9t2jrO0qsGMuGe5eWr1/iBE9imPx2Muffhdc+nD7xqPx9Pt9vq+B5bGeSpgnekipwa1yGOMKbf+HgFeBmYH83y9MSozmWSnQ1vehImGFhelVU2M9/NCbLu8gUnkZySxNsj19B8XH2fn4TpunpOvE4GrsBe0RC8iySKS2n4fuATvRdywEhMjTMlN05Y3YWLPEe+F2LF+XojtaM6YTD7YX4XL7QlUWKf489pi0hLjWFiQG7RzKBUoPWleuRR4HxgvImUicpuIfF5EyoBzgNdE5E1r2xwRed3aNRt4T0Q2AR8Brxlj3gjO0+ib6cPS2X6olrYgJgbVM7sr6gD/W9x0NHd0JvUtrqDNN3C4ppk3th3mujOH6UVY1S90ezHWGHNDF6te9rFtOXC5dX8/ML1P0YXI1Nw0Wl0edh2uY0quNpOz0+7DdcTHxjBsUO9HFD1ndAYisHZvJTNHBL7H6nMfHsBjDDeeNSLgx1YqGKK6Z2y79iGLdQ5Z++2qqGNsdkqfxnUflOxk0tABQRn3psXlZulHJVw0IYvhGTq8teofNNEDwwclkZYYp/X0YWBPRT3j+lBt027umEw2lFTT1OoOQFSfen3LISrrW7npnPyAHlepYNJEj7ejzbS8NG15Y7OapjYO1zYHLNG3uj18XFzV/cZ+WLL2AKMGJ3OuNdyCUv2BJnrL1Nw0dh2uo7ktsCVA1XN7AnAhtt2Z+QOJc0hAZ53aWFrNptJqbjp7BDE6ZaDqRzTRW6blpePyGHYcqrU7lKi1y0r0/o5a6Yt3esGBAR335i9ri0l2OvjizLyAHVOpUNBEb9ELsvbbfbiOlPhYctICM4HHuWMy2VYemOkFK+tbWL75EFfPzCM1IS4A0SkVOproLUPTEshMceoFWRvtrqhnbHZKwHqazh2TEbDpBZd9VEKr28NX9CKs6oc00Vu8F2TT9YKsjXZX+DerVHem5aWT7HT0uZ7e5fbw7AclzBubyZis3g3NoJSdNNF3MDU3jb1H6mlocdkdStSprG/hWENrn4Y+6CzOEcPZozJYu69vJfoV2ys4XNusTSpVv6WJvoPpw9LwGNiuF2RDLhBDH/gyZ0wmRZUNHKxu6vUx/ry2mLyBiVw4ISuAkSkVOproO2gf/mCTNT2cCp3dh9tb3AS2amTumAyAXlff7DhUy4dFVXzl7BF96q2rlJ000XeQlZrA0LQEbXljg10V9aQnxTE4JT6gxx2fnUpmipO1vUz0f3m/mIS4GK47c1hA41IqlDTRdzItL03nkLXBnoo6xmWnBnxsdxFhzuhM1uzzf3rBmsY2Xt5wkM8V5JKe5AxoXEqFkib6TqblpbO/soGapja7Q4kaxhh2VdQxrpeTjXRn7pgMjta1nBjrvqdeXFdKc5tHL8Kqfk8TfSdTrXr6bVp9EzKHa5upa3YF/EJsu7nWuDT+1NO7PYZnPjjA7PxBTMoZEJS4lAoVTfSdtPeQ3ayJPmR2V3hL2oEYzMyXvIFJjMhI8ivRF+46QklVIzfN0THnVf+nib6T9CQnwwclacepEDrR4iZIiR5gzuhMPuzh9II1jW38cfV+sgfEc+nkIUGLSalQ6XaGqWg0NS9Nm1iG0K6KOganxjMwOXgXPM8dk8nSj0rYfLCGM4afPOuUMYb9lQ28s6OCd3YcYd2B47g9hvsun0icQ8tCqv/TRO/D9Lw0Xtt8iKqGVgYFMfkorz0BHvrAl3NGW+3p91RyxvCBtLq8Y9W/s+MI/95ZQfGxRgAmDEnlG+eP4sIJ2ZwxPD2oMSkVKprofZiamw7A5rJq5o/X3pDB5PEYdlfUc/3s4LZTb59e8J+bytl5uI5Vu49S1+LCGRvDnNEZ3DZvFBdOyCI3PTGocShlB030PkzJHYAIbCmr0UQfZGXHm2hqcwe9RA9wwYTBLF65j5qmNq6YNpSLJmYzd0wGSU79GqjIpp9wH1IT4hiVmcwm7TgVdLsDONlId+68cCyfn5HLqMwUnSFKRRW90tSFaXnpbDlYbXcYEa99VqmxIRj+NyHOwZisVE3yKupoou/CtLw0KmpbqKhttjuUiLanoo7c9ESdtUmpINJE34UTHae0+iaodlXUB23oA6WUlyb6LkwamoYjRtiiHaeCxuX2sO9IfVA7SimlepDoReRpETkiIls7LLtGRLaJiEdEZp1m3wUisktE9orIvYEKOhQSnQ7GZqXoUAhBdKCqkVa3RxO9UkHWkxL9EmBBp2VbgS8Aq7raSUQcwGLgMmAScIOITOpdmPaYlpfG5rIav4e3VT0TiqEPlFI9SPTGmFVAVadlO4wxu7rZdTaw1xiz3xjTCiwDFvY6UhtMzUunqqG1T9PQqa7tqqhDBJ1wW6kgC2Y7+lygtMPjMuCsrjYWkUXAIoDs7GwKCwv9Oll9fb3f+3THVeMGYOmbazlzSGBeqmDEGQyhiHPN1mayEoUP167u9TH09QwsjTOwwiXOYCZ6X42Vu6wDMcY8ATwBMGvWLDN//ny/TlZYWIi/+3TnHJebX368gn1tA/n++WcEZPajYMQZDKGI82efvMu0/GTmz+/yMk+39PUMLI0zsMIlzmC2uikDOg5gkgeUB/F8ARcf6+A7F43ljW2HWfZxafc7qB5rcbkpqmwIydAHSkW7YCb6j4GxIjJSRJzA9cArQTxfUHzz/NHMG5vJT17Zxs7DtXaHEzGKKhtwe0xIhj5QKtr1pHnlUuB9YLyIlInIbSLyeREpA84BXhORN61tc0TkdQBjjAu4A3gT2AG8aIzZFqwnEiwxMcJD1xaQmhDHHc9voLHVZXdIEWHXiRY3eiFWqWDrto7eGHNDF6te9rFtOXB5h8evA6/3OrowMTg1nkeuL+DGpz7kx//cxgPXTLc7pH5vd0UdsTHCqExN9EoFm/aM7aG5YzK544IxvLS+jJc3lNkdTr+3u6KekZnJOGP1I6hUsOm3zA/fuWgss/MHcd/LW9l/tN7ucPq13RV12lFKqRDRRO+HWEcMj9xQQHxsDN9+fgPNbW67Q+qXmlrdlFQ1aqJXKkQ00ftpaFoiD147nR2Havn5azvsDqdf2nukHmP0QqxSoaKJvhcunJDN1+aN5JkPDvCvLYfsDqff2RXCWaWUUproe+0/L53A9GHp/OBvmymtarQ7nH5ld0UdztgYRgxKsjsUpaKCJvpecsbG8PsbZgBwx9INtLo8NkfUf+yuqGP04BRiHfrxUyoU9JvWB8MGJfHrL05jU2k1v1nR3WCeqt3uw3WM1/p5pUJGE30fXT51KDeePZwnVu3n3zsr7A4n7NU2t1Fe06z180qFkCb6ALj/iklMGJLK917cxBGdTPy09lR4+x+My9JEr1SoaKIPgIQ4B7//0hk0trr5wd8264xUp7HbanEzXkv0SoWMJvoAGZOVwg8vm0DhrqM8/1GJ3eGErd0VdSQ5HeSmJ9odilJRI5gTj0Sdm87J5+0dR/jZ8h3MHZ1Jfmay3SGFnd0VdYzNSiEmpu+TuKjoJSIUFRXR3BzeVaVpaWns2BHYjpUJCQnk5eURFxfX43000QdQTIzwwDXTuPThVXz3xY28+PVztAlhB42tLnYequPCCVl2h6L6ueTkZFJTU8nPzw/IzG/BUldXR2pq4KopjTEcO3aMsrIyRo4c2eP9NAsF2NC0RH76uSl8UlLN46v22x1O2NhUWs2Vj77HsYZWLpqoiV71jcPhICMjI6yTfDCICBkZGX7/J6OJPgiump7DFdOG8vBbu9l6sMbucGzlcnt49J09fOGxtTS3uXn+9rNYMGWo3WGpCBBtSb5db563JvogEBF+/rkpDEp2cvcLG6N2lMviygauefx9HnprN1dOG8q//uM85ozJtDsspaKOJvogSU9y8sA109lzpJ7fvBldvWaNMSz9qITLH13NviP1PHrDDB65fgZpiT2/eKRUuHM4HBQUFDBlyhQ++9nPUl1dDUBxcTEiwn/913+d2LayspK4uDjuuOMOAHbt2sX8+fMpKChg4sSJLFq0CIDCwkLS0tIoKCg4cXv77bf7HKsm+iA6f9xgvnL2CJ5aU8TafZV2hxMSlfUtfO0v6/nh37dQMCydN+8+j6um59gdllIBl5iYyMaNG9m6dSuDBg1i8eLFJ9aNGjWK5cuXn3j80ksvMXny5BOP77rrLu6++242btzIjh07uPPOO0+smzdvHhs3bjxxu/jii/scq7a6CbIfXj6B9/ZW8v0XN/HG3efZHU5Qvb29gnv/vpnaZhf3XzGRW+eO1GaUKuj+59VtbC+vDegxJ+UM4Mefndz9hpZzzjmHzZs3n3icmJjIxIkT+eSTTzj//PN54YUXuPbaaykvLwfg0KFD5OXlndh+6tSpgQveB030QZbkjOWha6dz9R/e5yevbOOqHjQ4McawqayGF9eV8smB48Q5YkiIiyE+1nHib3ynx4lxDiYMTWXO6AxSE0JXRdLiclNa1chT7xWz9KMSJg4dwHO3F2jPVxU13G4377zzDrfddttJy6+//nr+9re/MXr0aBwOBzk5OScS/d13382FF17InDlzuOSSS7jllltIT08HYPXq1RQUFJw4Tvsx+kITfQjMGD6Qb88fzaP/3ktOQTzzu9juaF0L/9hwkBfXlbLnSD0JcTGcPSoDAZrbPDS2uqhq8NDictPc5qHF5b3f0uah1e0dJjk2Rjhj+EDOG5fJeeMGMyUnrc+l6ja3h7LjTRRXNlBU2UDxMe/fosoGyqub8BgQga+fP4rvfmYc8bGOPp1PKX/4U/IOpKamJgoKCiguLmbmzJl85jOfOWn9ggULuO+++xg2bBjXXXfdSetuueUWLr30Ut544w3++c9/8vjjj7Np0ybAW3XTsdonEDTRh8idF41l5a6jLNlWw011zWSlJgDeJLpy5xFeWl/Gyp1HcHkMZwxP51dfmMoV04b2uHTe6vLwSclxVu0+yqo9R/nNit38ZsVuBibFce7YwZw31pv4swcknLJvXXMbh2uaOVTTzOGaZsprmvhkZwtP7/+IkmMNlB5vwu35dPye1PhY8jOTmTF8IF+YkUt+ZjLT8tIYowOVqSjSXkdfU1PDlVdeyeLFi7nrrrtOrHc6nRQUFPDggw+ybds2Xn311ZP2z8nJ4dZbb+XWW29lypQpbN26NWixaqIPkThHDA9fN53LfruKe/+2hXsvm8BL60p5ecNBKutbGZwaz23zRnLNzGGMyfJ/rHZnrLf0f/aoDH6wYAKV9S2s2VvJu7uPsnpPJa9u8v7LOGFIKpOGDuBofcuJ5F7f4jrleAOcwojBrUzOSeOKaUPJz0hmZGYy+ZnJZCQ7o7YNs1KdpaWl8eijj7Jw4UK++c1vnrTuzjvv5OKLLyYjI+Ok5W+88QYXXXQRcXFxHD58mGPHjpGbm8vOnTuDEqMm+hAak5XKteOcPLfzCP/eeYTYGOGiiVlcO2sY548bHNDhEjJT4llYkMvCglyMMew4VMeqPUdZveco7+8/RtaABEYPTmHumExy0hMYkpbI0LQEhgxIIHtAAmvfW8X8+ecGLB6lItmMGTOYPn06y5YtY968eSeWT5w4kdmzZ5+y/YoVK/jOd75DQoL3P+wHHniAIUOGsHPnzlPq6O+//36uvvrqPsWniT7ELhoRy6Cc4aQlxvG5GblkpsQH/ZwiwqScAUzKGcA3zu/bRR2llFd9ff1JjztWzfiqhrn55pu5+eabAXjooYd46KGHTtlm/vz51NQEvjd9t4leRJ4GrgSOGGOmWMsGAS8A+UAxcK0x5riPfYuBOsANuIwxswIVeH8VI8L3LhlvdxhKqSjSk7qCJcCCTsvuBd4xxowF3rEed+UCY0yBJnmllLJHt4neGLMKqOq0eCHwZ+v+n4HPBTYspZQ6vWidya03z1t6spOI5APLO1TdVBtj0jusP26MGehjvyLgOGCAx40xT5zmHIuARQDZ2dkzly1b5tcTqa+vJyXF/9YqoaZxBpbGGVj9Jc6EhARyc3NJS0sL6xZgbrcbhyNw/UqMMdTU1FBRUXHKNYILLrhgfVc1J8G+GDvXGFMuIlnAWyKy0/oP4RTWj8ATALNmzTLz58/360SFhYX4u48dNM7A0jgDq7/E+e677+JyuTh48KDdoZxWc3PziZY1gZKQkMD06dNDMsNUhYgMNcYcEpGhwBFfGxljyq2/R0TkZWA24DPRK6VUTxlj/JphyS6FhYXMmDHD7jB6PXrlK8BXrftfBf7ZeQMRSRaR1Pb7wCVA8Lp+KaWU8qnbRC8iS4H3gfEiUiYitwG/Aj4jInuAz1iPEZEcEXnd2jUbeE9ENgEfAa8ZY94IxpNQSinVtW6rbowxN3Sx6iIf25YDl1v39wPT+xSdUkqpPutRq5tQE5GjwAE/d8sE+sPsHhpnYGmcgaVxBlYo4xxhjBnsa0VYJvreEJF1/aFTlsYZWBpnYGmcgRUucepUgkopFeE00SulVISLpETfZa/bMKNxBpbGGVgaZ2CFRZwRU0evlFLKt0gq0SullPJBE71SSkW4fp/oRWSBiOwSkb0icrpx8UNKRIaJyEoR2SEi20TkO9byQSLylojssf6eMuqnHUTEISIbRGS59Tjs4hSRdBH5q4jstF7Xc8I0zrut93yriCwVkYRwiFNEnhaRIyKytcOyLuMSkR9a36tdInKpzXE+YL3vm0XkZRFJD8c4O6z7vogYEcm0O07o54leRBzAYuAyYBJwg4hMsjeqE1zA94wxE4GzgW9bsfkzaUsofQfY0eFxOMb5CPCGMWYC3l7XOwizOEUkF7gLmGUN6+0Aric84lxCDycRsj6r1wOTrX3+z/q+2RXnW8AUY8w0YDfwwzCNExEZhndomJIOy+yMs38neryjYe41xuw3xrQCy/BOimI7Y8whY8wn1v06vEkplzCctEVE8oArgCc7LA6rOEVkAHAe8BSAMabVGFNNmMVpiQUSRSQWSALKCYM4/ZxEaCGwzBjTYowpAvbi/b7ZEqcxZoUxxmU9/ADIC8c4LQ8DP8A7D0c72+KE/p/oc4HSDo/LrGVhxZq4ZQbwIZBtjDkE3h8DIMvG0Nr9Fu8H09NhWbjFOQo4CvzJqmJ60hoVNaziNMYcBH6DtzR3CKgxxqwgzOLsoKu4wvm7dSvwL+t+WMUpIlcBB40xmzqtsjXO/p7ofU0tE1btRUUkBfgb8B/GmFq74+lMRNonfl9vdyzdiAXOAB4zxswAGgiP6qSTWHXcC4GRQA6QLCI32htVr4Tld0tE7sNbLfpc+yIfm9kSp4gkAfcB/+1rtY9lIYuzvyf6MmBYh8d5eP9NDgsiEoc3yT9njPm7tbjCmqyF003aEkJzgatEpBhv1deFIvIs4RdnGVBmjPnQevxXvIk/3OK8GCgyxhw1xrQBfwfmEH5xtusqrrD7bonIV4ErgS+bTzsAhVOco/H+wG+yvk95wCciMgSb4+zvif5jYKyIjBQRJ96LHa/YHBMAIiJ465N3GGMe6rCq20lbQskY80NjTJ4xJh/v6/dvY8yNhF+ch4FSERlvLboI2E6YxYm3yuZsEUmyPgMX4b0+E25xtusqrleA60UkXkRGAmPxzithCxFZANwDXGWMaeywKmziNMZsMcZkGWPyre9TGXCG9dm1N05jTL++4R3/fjewD7jP7ng6xHUu3n/NNgMbrdvlQAbe1g17rL+D7I61Q8zz8U4CTzjGCRQA66zX9B/AwDCN83+AnXhnVHsGiA+HOIGleK8btOFNQredLi681RD7gF3AZTbHuRdvHXf7d+kP4Rhnp/XFQKbdcRpjdAgEpZSKdP296kYppVQ3NNErpVSE00SvlFIRThO9UkpFOE30SikV4TTRKwWISH0AjhEvIm+LyEYRuS4QcSkVCLF2B6BUBJkBxBljCuwORKmOtESvVBdEZLSIvCEi60VktYhMsJZ/VkQ+tAZXe1tEskUkC3gWKLBK9KPtjV6pT2mHKaXwVt0YY1I6LXsH+IYxZo+InAX80hhzoTVwWbUxxojI7cBEY8z3RGQ+8H1jzJUhfwJKnYZW3SjlgzXq6BzgJe+QNYB3KAPwDkj1gjUImBMoCn2ESvWcJnqlfIvBW2ov8LHud8BDxphXrFL8T0IXllL+0zp6pXww3rkDikTkGvCORioi063VacBB6/5Xfe2vVDjRRK+UV5KIlHW4fRf4MnCbiGwCtvHpNJU/wVulsxqotCdcpXpOL8YqpVSE0xK9UkpFOE30SikV4TTRK6VUhNNEr5RSEU4TvVJKRThN9EopFeE00SulVIT7/6VZ/c/Lts87AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "leaf_vals = ee.List.sequence(1, 150, 5)\n",
    "\n",
    "def leafs(t):\n",
    "    cartRegression = ee.Classifier.smileCart(minLeafPopulation=t).setOutputMode('REGRESSION').train(\n",
    "    features= training, \n",
    "    classProperty= 'Percent_Tree_Cover', \n",
    "    inputProperties= predictionBands\n",
    "    )\n",
    "\n",
    "    cartRegressionTesting = testing.classify(cartRegression, 'y_hat');\n",
    "\n",
    "    # // Calculate the RMSE for CART regression and print the result\n",
    "    regression=cartRegressionTesting.map(set_sqError)\n",
    "    regression = ee.FeatureCollection(regression)\n",
    "    MSE = regression.reduceColumns(ee.Reducer.mean(),['sqError'])\n",
    "    RMSE =  ee.Number(MSE.get('mean')).sqrt()        \n",
    "    return RMSE\n",
    "\n",
    "\n",
    "rmse_leaf=leaf_vals.map(leafs)\n",
    "value_info = rmse_leaf.getInfo()\n",
    "\n",
    "# Extracting gamma value for minimum RMSE\n",
    "min_rmse=min(value_info)\n",
    "print(\"Min RMSE =\", min_rmse)\n",
    "\n",
    "index=value_info.index(min(value_info))\n",
    "leaf_range = leaf_vals.getInfo()\n",
    "min_leaf =leaf_range[index]\n",
    "print(\"Min Leaf =\", min_leaf)\n",
    "\n",
    "# print(rmse_gama.getInfo())\n",
    "\n",
    "# Plotting the interations\n",
    "import pandas as pd\n",
    "df =pd.DataFrame(value_info,columns=['RMSE'])\n",
    "df['Leaf'] = leaf_vals.getInfo() \n",
    "\n",
    "ax =df.plot.line(x='Leaf', \n",
    "             y='RMSE',\n",
    "             title= 'Impact of Leaf numbers in Cart'\n",
    "             )\n",
    "ax.grid()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Assignment\n",
    "\n",
    "Recreate one notebook for Logan, UT, a dry climate with a variety of forest (natural to agriculture, Path 38, Row 31) and another notebook for New Orleans, Path 22 Row 39). Make sure your Cart, SVR and Random Forest in each notebook uses the tuned hyperparameters when discussing your found RMSE . Are the 4 models behaving similarly? Which model was able to model adequately Tree Percentage at Landsat resolution? Discuss it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.9.7"
  },
  "latex_envs": {
   "LaTeX_envs_menu_present": true,
   "autoclose": false,
   "autocomplete": true,
   "bibliofile": "biblio.bib",
   "cite_by": "apalike",
   "current_citInitial": 1,
   "eqLabelWithNumbers": true,
   "eqNumInitial": 1,
   "hotkeys": {
    "equation": "Ctrl-E",
    "itemize": "Ctrl-I"
   },
   "labels_anchors": false,
   "latex_user_defs": false,
   "report_style_numbering": false,
   "user_envs_cfg": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}