{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This example demonstrates how to normalize the color contrast between 2 images from Landsat 8 imagery covering Queenstown, New Zealand. It runs through the steps of 1) acquiring imagery scenes from the EarthAI Catalog, 2) using RasterFrames to read and normalize imagery, 3) viewing imagery on a slippy map, and 4) showing histograms of cell values.\n",
    "\n",
    "Image normalization is the process of scaling a range of pixel values. Images of the same area acquired at different times of day or taken from different satellites may have high color contrast, which makes it hard to compare. Normalizing the color contrast between images is important for a variety of image applications, such as mosaicking and machine learning. \n",
    "\n",
    "In this example, we will:\n",
    "\n",
    "1. acquire two Landsat 8 scenes that cover a specified area of interest;\n",
    "1. read in the imagery using the EarthAI chip reader;\n",
    "1. display the imagery covering our area of interest on a slippy map;\n",
    "1. display histograms of pixel values;\n",
    "1. normalize the imagery using z-score normalization; and\n",
    "1. normalize the imagery using min-max normalization.\n",
    "\n",
    "*Note: If you are following along by running the attached companion notebook in EarthAI Notebook, this example takes about 4 minutes to run end-to-end on a Dedicated Instance type.*\n",
    "\n",
    "# Import Libraries\n",
    "\n",
    "We will start by importing all of the Python libraries used in this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n",
      "Importing EarthAI libraries.\n",
      "EarthAI version 1.5.0; RasterFrames version 0.9.0; PySpark version 2.4.6\n",
      "\n",
      "Creating SparkSession...\n",
      "  SparkSession is available as `spark`.\n"
     ]
    }
   ],
   "source": [
    "from earthai.init import *\n",
    "import earthai.chipping.strategy\n",
    "import pyspark.sql.functions as F\n",
    "\n",
    "from shapely.geometry import Polygon\n",
    "import geopandas\n",
    "import ipyleaflet\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import math"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Query the EarthAI Catalog\n",
    "\n",
    "The polygon below is a bounding box surrounding Queenstown, New Zealand in the standard \"EPSG:4326\" projection. This projection is required for querying the [EarthAI Catalog API](https://docs.astraea.earth/hc/en-us/articles/360043639052-The-EarthAI-Catalog-Query-API).  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "poly = Polygon([(168.59011352062228,-45.05876119993573),\n",
    "                (168.71783092617991,-45.05876119993573),\n",
    "                (168.71783092617991,-44.982437864876836),\n",
    "                (168.59011352062228,-44.982437864876836),\n",
    "                (168.59011352062228,-45.05876119993573)\n",
    "               ])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can now input this geometry into `earth_ondemand.read_catalog`. We use [Landsat 8 Collection 1 Tier 1](https://www.usgs.gov/media/files/landsat-collection-1-level-1-product-definition), and set `max_cloud_cover = 2` to select scenes with very low cloud coverage.\n",
    "\n",
    "It is possible to perform image normalization with cloudy scenes, but the additional step of masking out clouds would be necessary as clouds will skew the statistics used to calculate normalization parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "0cbfd85ffacd4a65818f72acac123fd9",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "HBox(children=(FloatProgress(value=0.0, max=11.0), HTML(value='')))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "\n"
     ]
    }
   ],
   "source": [
    "cat = earth_ondemand.read_catalog(poly,\n",
    "                                  max_cloud_cover = 2,\n",
    "                                  collections = 'landsat8_l1tp')\n",
    "\n",
    "cat['queenstown_geometry'] = poly"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We filter the catalog down to just two Landsat 8 scenes to speed up calculation in this example."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "cat = cat[cat['id'].isin(['LO08_L1TP_075091_20150210_20170413_01_T1_L1TP','LC08_L1TP_076091_20170222_20170301_01_T1_L1TP' ])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The code block below uses [ipyleaflet](https://ipyleaflet.readthedocs.io/en/latest/) to display the Queenstown boundary and the extents of the intersecting Landsat 8 scenes. If we wanted to create a mosaic of these Landsat 8 scenes, we would first want to normalize the bands in each scene to prevent the mosaic from looking choppy in the overlapping area. \n",
    "\n",
    "In this example, we will just look at imagery located within the Queenstown boundary to illustrate 2 different normalization techniques: z-score normalization and min-max normalization."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "53cc640c207a4610ba9e7c9cb2648d94",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[-45.02059953240628, 168.6539722234011], controls=(ZoomControl(options=['position', 'zoom_in_text',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m = ipyleaflet.Map(center=(poly.centroid.y, poly.centroid.x), zoom=8)\n",
    "\n",
    "queenstown_layer = ipyleaflet.GeoData(geo_dataframe=geopandas.GeoDataFrame([{'geometry':poly}]), \n",
    "                                      name='Queenstown', \n",
    "                                      style={'fillColor': '#f003fc', 'color': '#f003fc'})\n",
    "\n",
    "scenes_layer = ipyleaflet.GeoData(geo_dataframe=cat[['geometry']], \n",
    "                                  name='Landsat 8 Scenes', \n",
    "                                  style={'fillColor': '#32a852', 'color': '#32a852'})\n",
    "\n",
    "m.add_layer(queenstown_layer)    \n",
    "m.add_layer(scenes_layer)    \n",
    "m.add_control(ipyleaflet.LayersControl())\n",
    "m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Read in Landsat 8 Imagery\n",
    "\n",
    "Since the Landsat 8 scenes cover a much larger geographic area than Queenstown, we will make use of EarthAI's chip reader to load in only the data that intersects with the Queenstown bounding box. We specify `geometry_col_name = 'queenstown_geometry'` to specify the region of the scenes to load in. The data will be read in as a RasterFrame with uniform tiles with dimensions of 256 x 256. We use the `FeatureAlignedGrid` chipping strategy since we want the mosaic to be gridded in reference to the Queenstown geometry. \n",
    "\n",
    "We read in bands B4 (red), B3 (green), and B2 (blue) and convert the cell type for each band to `uint16`, whose NoData value is 0. This will cause any zero-valued cells in the bands to be considered NoData. In Landsat 8, these areas correspond to the band quality assurance fill areas and we don't want any zeros to skew the normalization statistics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf = spark.read.chip(cat.reset_index(),\n",
    "                     catalog_col_names=[\"B4\", \"B3\", \"B2\"],\n",
    "                     geometry_col_name = 'queenstown_geometry',\n",
    "                     chipping_strategy=earthai.chipping.strategy.FeatureAlignedGrid(256, 256)) \\\n",
    "            .select('datetime',\n",
    "                   rf_convert_cell_type('B4', 'uint16').alias('red'),\n",
    "                   rf_convert_cell_type('B3', 'uint16').alias('grn'), \n",
    "                   rf_convert_cell_type('B2', 'uint16').alias('blu'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We split the Landsat 8 scenes into 2 separate RasterFrames."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_scene1 = rf.filter(F.year('datetime') == 2015)\n",
    "rf_scene2 = rf.filter(F.year('datetime') == 2017)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To view the imagery, we use [ipyleaflet](https://ipyleaflet.readthedocs.io/en/latest/) and the RasterFrames slippy map writer to show our imagery overlaid on a basemap. Using the layer control on the top left of the map, you can toggle between the scene from 2015 and the scene from 2017 to see the color contrast between images."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "application/vnd.jupyter.widget-view+json": {
       "model_id": "19481f2fabc748b8beae3a6785a9380e",
       "version_major": 2,
       "version_minor": 0
      },
      "text/plain": [
       "Map(center=[-45.02059953240628, 168.6539722234011], controls=(ZoomControl(options=['position', 'zoom_in_text',…"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "m1 = ipyleaflet.Map(center=(poly.centroid.y, poly.centroid.x), zoom=12)\n",
    "\n",
    "rf_scene1_map = rf_scene1.select('red', 'grn', 'blu').to_map('Scene_2015', rendering_mode=earthai.RenderingMode.UNIFORM)\n",
    "rf_scene2_map = rf_scene2.select('red', 'grn', 'blu').to_map('Scene_2017', rendering_mode=earthai.RenderingMode.UNIFORM)\n",
    "\n",
    "m1.add_layer(rf_scene1_map.layers[-1]) \n",
    "m1.add_layer(rf_scene2_map.layers[-1]) \n",
    "\n",
    "m1.add_control(ipyleaflet.LayersControl())\n",
    "m1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `rf_agg_approx_histogram` function computes a count of cell values across all of the rows of tile in a RasterFrame. We use this function to compute histograms for each band and Landsat 8 scene. A histogram will show us if the distribution and range of cell values for each band are similar across scenes. If the histograms don't look similar, then image normalization may be required to scale the values to a similar data distribution. This is important for speeding up algorithm convergence and making sure one feature doesn't dominate another in machine learning applications. It can also be used to remove color contrast when creating a mosaic of images.\n",
    "\n",
    "As you can see in the histograms below, the cell values are shifted to the right for the scene from 2015."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABJEAAAEVCAYAAABQV3TbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/QUlEQVR4nO3de7iV5Xng/+8dUMQYPKIFNlRSwAPW0Ih0m/bXkFjEpDVoqxFjKmO0dNSYxEyNMs5MbKdc6m86MdEYDYmJmlQJQzORcRQFgybNDxVQooJRiRDYQpQQteSgAXL//ljvJovN2kf2Ouy9v5/rWtd61/0e9vO+wPOw7/UcIjORJEmSJEmSOvK2ehdAkiRJkiRJjc8kkiRJkiRJkjplEkmSJEmSJEmdMokkSZIkSZKkTplEkiRJkiRJUqdMIkmSJEmSJKlTJpEkSVKfFRFTI6JlH6/xi4h4Z2+VSZIkqb8yiSRJknpVRGyIiF8XyZnXIuL/RsToOpUlI2Jcm9i1EfHN1s+ZeVBmvtTJdfY5WSVJktTXmUSSJEnVcEZmHgSMAF4Bbq5zeRpaRAyqdxkkSZI6YxJJkiRVTWa+CSwEjm+NRcRfRMRTEfHvEbEpIq4t23d00XtoVkRsjIifRcQ1ZfuHRsQdRQ+ntcDJ+1rG8t5KEfHBiFgbEdsj4uWI+PuIeDvwADCy6F31i4gYGRFDIuLzEbG5eH0+IoaUXfczEbGl2Hdxm59zR0TcGhH3R8Qvgfd18blcWOx7LSL+Y0ScHBFPR8TrEfHFsuPHRcSjEfFG8Qy/ta/PSZIkaXC9CyBJkvqviDgQOBd4rCz8S+ACYA1wArAkIlZn5nfKjvlT4BhgAvBERHw7M58DPgv8QfFqTe70ptuBD2fm9yPiUGBsZv4yIj4AfDMzm8ru7R+BZmASkMC9wH8B/mtEnA58GjgVWA98ucLP+gjwQeAvgf2La3X2XP4YGA/8GbAIWAz8ObAf8FRE/K/MfBT478BDwPuKa0/ex+ciSZJkTyRJklQV34mI14F/B6YB/6N1R2Y+kpnPZOZvM/Np4B7gvW3O/4fM/HVm/hD4IfCuIv5hYG5m/jwzNwE3daEsTxY9dV4vynR1B8fuAI6PiGGZ+VpmPtnBsecD/5iZr2bmVuAfgL8pK+fXM3NNZv6q2NfWvZn5g+I5vNnF5/Lfi2MfopSMu6f4+S8D3wf+qOw+fh8YWRz/bx3chyRJUpeYRJIkSdVwZmYeAgwBPg48GhG/BxARfxwRyyJia0S8AfxH4Ig25/+0bPtXwEHF9khgU9m+n3ShLO/OzENaX8D1HRz715R6B/2kGA52SgfHjmzz839SxCqVs3y7YqyLz+WVsu1fV/jc+pw+AwSlXlxrIuJjHdyHJElSl5hEkiRJVZOZuzLz28AuSkPUAO6mNBRrdGYeDNxGKeHRFVuA8pXexvRWWQEyc0VmzgCOBL4DLGjdVeHwzZR6+5SXZXNZOZvK9lVana7tNffluex54cyfZubfZuZI4O+AL7VdpU6SJKm7TCJJkqSqiZIZwKHAc0X4HcDPM/PNiJhCaW6grloAzImIQyOiCbi8F8u6f0ScHxEHZ+YOSkPxdhW7XwEOj4iDy065B/gvETE8Io4A/hvwzbJyXhgRxxXzQv23LhRhX55L23s5p3g+AK9RSljt6uAUSZKkTplEkiRJ1fB/IuIXlBIxc4FZmbmm2Hcp8I8RsZ1ScmVBO9eo5B8oDRtbT2ni6G/0XpGB0pxGGyLi3ykNJ/soQGb+iFLS6KVibqWRwD8BK4GngWeAJ4sYmfkApfmalgHrgOXF9d/q4Gfvy3Np62Tg8eLPYBHwycxcvw/XkyRJIjIr9c6WJElSb4mI44BngSGZubPe5ZEkSeoJeyJJkiRVQUScVQyROxS4Afg/JpAkSVJfZhJJkiSpOv4O2Ar8mNJ8RJfUtziSJEn7xiSSVEURcW1EfLPzIyVJ/U1mnp6ZB2fmYZl5VmZuqXeZJEn1EREbIuLPK8SnRkRLPcok9YRJJEmSJEmSJHXKJJLUDRExuN5lkCT1XbYjkiSpLzOJJHWi6Hp6VUQ8DfwyIv40Iv6/YonnH0bE1LJjx0bEoxGxPSKWAEfUq9ySpNqJiHdHxFNF/f+/IuJbEfFPrcMUinbkp8DXi6HOCyLiruL4NRExud73IEmqupMjYm1EvBYRX4+IA9oeEBEZEePKPt8REf9U22JK7TOJJHXNecBfAO8E7gX+CTgM+HvgXyNieHHc3cAqSsmj/w7Mqn1RJUm1FBH7A/8buINS23APcFbZIb9XxH8fmF3EPgTMBw4BFgFfrE1pJUl1dD4wHfgDYALwX+pbHKn7TCJJXXNTZm4CPgrcn5n3Z+ZvM3MJsBL4YESMAU4G/mtmvpWZ3wP+Tx3LLEmqjWZgMKW2Ykdmfht4omz/b4HPFm3Dr4vYvxVtyS7gG8C7altkSVIdfDEzN2Xmz4G5lL6olvoUk0hS12wq3n8fOKcYyvZ6RLwO/CkwAhgJvJaZvyw77ye1LaYkqQ5GAi9nZpbFNpVtb83MN9uc89Oy7V8BBzhfkiT1e+Vtw08otR9Sn2ISSeqa1l8MNgHfyMxDyl5vz8zrgS3AoRHx9rLzxtS8pJKkWtsCjIqIKIuNLttOJEnas20YA2yucMyvgAPLPv9eVUskdZNJJKl7vgmcERHTI2JQRBxQTJralJk/oTS07R8iYv+I+FPgjPoWV5JUA8uBXcDHI2JwRMwAptS5TJKkxnNZRDRFxGHAfwa+VeGY1cBHit81TgfeW8sCSp0xiSR1QzEv0gxKlf5WSj2TruR3/5Y+Avwx8HPgs8BddSimJKmGMvM3wF8BFwGvU5o/7z7grToWS5LUeO4GHgJeKl6VVl37JKUvol+nNBH3d2pUNqlLYs/h+5IkSdpXEfE4cFtmfr3eZZEkSeot9kSSJEnaRxHx3oj4vWI42yzgRGBxvcslSZLUm1wFRJIkad8dAywADgJ+DJydmVvqWyRJkqTe5XA2SZIkSZIkdcrhbJIkSZIkSepUnx3OdsQRR+TRRx9d72JIUsNZtWrVzzJzeL3LUW+2E5JUme1Eie2EJFXWUTvRZ5NIRx99NCtXrqx3MSSp4UTET+pdhkZgOyFJldlOlNhOSFJlHbUTDmeTJEmSJElSp0wiSZIkSZIkqVMmkSRJkiRJktSpPjsnkqT627FjBy0tLbz55pv1LsqAdMABB9DU1MR+++1X76JIUkW2E/VlOyGp0dlO1FdP2gmTSJJ6rKWlhXe84x0cffTRRES9izOgZCbbtm2jpaWFsWPH1rs4klSR7UT92E5I6gtsJ+qnp+2Ew9kk9dibb77J4YcfboVfBxHB4Ycf7rc2khqa7UT92E5I6gtsJ+qnp+2ESSRJ+8QKv3589pL6Auuq+vHZS+oLrKvqpyfP3iSSJEmSJEmSOuWcSJJ6zY1LXujV610xbUKvXk+SVF+2E5KkjthOND57InVm2XWVX5Iaxty5c5k4cSInnngikyZN4vHHH69LOa688kqOPfZYTjzxRM466yxef/313fuuu+46xo0bxzHHHMODDz64O37NNdcwevRoDjrooD2udccddzB8+HAmTZrEpEmT+OpXv1qr21CDu3HJC3u9JHXMdkKqj0ptlu2WGpHtRNeZRJLUpy1fvpz77ruPJ598kqeffpqlS5cyevToupRl2rRpPPvsszz99NNMmDCB664rJZzXrl3L/PnzWbNmDYsXL+bSSy9l165dAJxxxhk88cQTFa937rnnsnr1alavXs3FF19cs/uQpP7EdkKS1BHbie4xiSSpT9uyZQtHHHEEQ4YMAeCII45g5MiRrFixgve85z28613vYsqUKWzfvp1du3Zx5ZVXcvLJJ3PiiSfy5S9/GYBHHnmEqVOncvbZZ3Psscdy/vnnk5kArFq1ive+972cdNJJTJ8+nS1btrRbltNOO43Bg0ujhJubm2lpaQHg3nvvZebMmQwZMoSxY8cybty43RV9c3MzI0aMqNrzUd/lt7dS77CdkCR1xHaie0wiSerTTjvtNDZt2sSECRO49NJLefTRR/nNb37Dueeeyxe+8AV++MMfsnTpUoYOHcrtt9/OwQcfzIoVK1ixYgVf+cpXWL9+PQBPPfUUn//851m7di0vvfQSP/jBD9ixYweXX345CxcuZNWqVXzsYx/jmmuu6VK5vva1r/GBD3wAgJdffnmPbzOampp4+eWXO73Gv/7rv3LiiSdy9tlns2nTph48HUmS7UT9RcTXIuLViHi2LPY/IuJHEfF0RPzviDikbN+ciFgXEc9HxPSy+EkR8Uyx76YolhWKiCER8a0i/nhEHF3L+5PUt9lOdI8Ta0vq0w466CBWrVrF97//fZYtW8a5557LNddcw4gRIzj55JMBGDZsGAAPPfQQTz/9NAsXLgTgjTfe4MUXX2T//fdnypQpNDU1ATBp0iQ2bNjAIYccwrPPPsu0adMA2LVrV5ey/HPnzmXw4MGcf/75ALu/hSjX2XKaZ5xxBueddx5DhgzhtttuY9asWXz3u9/t4lORJLWynWgIdwBfBO4qiy0B5mTmzoi4AZgDXBURxwMzgYnASGBpREzIzF3ArcBs4DHgfuB04AHgIuC1zBwXETOBG4Bza3Jnkvo824nuMYkkqc8bNGgQU6dOZerUqfzhH/4ht9xyS8VKNTO5+eabmT59+h7xRx55ZHf31dbr7dy5k8xk4sSJLF++vMtlufPOO7nvvvt4+OGHd5ehqalpj8x/S0sLI0eO7PA6hx9++O7tv/3bv+Wqq67qchkkSXuynaivzPxe295BmflQ2cfHgLOL7RnA/Mx8C1gfEeuAKRGxARiWmcsBIuIu4ExKSaQZwLXF+QuBL0ZEZKXfuiSpAtuJrjOJJKnX1GMJzeeff563ve1tjB8/HoDVq1dz3HHHsXjxYlasWMHJJ5/M9u3bGTp0KNOnT+fWW2/l/e9/P/vttx8vvPACo0aNavfaxxxzDFu3bmX58uWccsop7NixgxdeeIGJEydWPH7x4sXccMMNPProoxx44IG74x/60If4yEc+wqc//Wk2b97Miy++yJQpUzq8ry1btuz+lmLRokUcd9xx3X00ktRwbCdsJ9rxMeBbxfYoSkmlVi1FbEex3Tbees4mgKJn0xvA4cDPqlhmSVVgO9H47YRJJEl92i9+8Qsuv/xyXn/9dQYPHsy4ceOYN28eF154IZdffjm//vWvGTp0KEuXLuXiiy9mw4YNvPvd7yYzGT58ON/5znfavfb+++/PwoUL+cQnPsEbb7zBzp07+dSnPtVupf/xj3+ct956a3d31ebmZm677TYmTpzIhz/8YY4//ngGDx7MLbfcwqBBgwD4zGc+w913382vfvUrmpqauPjii7n22mu56aabWLRoEYMHD+awww7jjjvu6O1HJ0kDgu1EY4uIa4CdwL+0hioclh3EOzqn0s+bTWlIHGPGjOlWWSX1T7YT3RN9tZfn5MmTc+XKldX/Qcuuqxx/35zq/2ypwT333HN9/ZvPPq/Sn0FErMrMyXUqUsOoWTtRJd1Zia0e39pJXWE7UX+N0k4Uw9nuy8wTymKzgP8InJqZvypicwAy87ri84OUhqptAJZl5rFF/Dxgamb+Xesxmbk8IgYDPwWGdzacra+3E31Be22Z7ZZa2U7UX3fbCVdnkyRJklRTEXE6cBXwodYEUmERMLNYcW0sMB54IjO3ANsjorlYle0C4N6yc2YV22cD33U+JEmqDpNIktRNl112GZMmTdrj9fWvf73exaqbdpZuPiwilkTEi8X7oWX7XLpZUr9mO7GniLgHWA4cExEtEXERpdXa3gEsiYjVEXEbQGauARYAa4HFwGXFymwAlwBfBdYBP6Y0qTbA7cDhxSTcnwaurs2dSVLP9OV2wjmRJKmbbrnllnoXodHcwd5LN18NPJyZ10fE1cVnl26WNCDYTuwpM8+rEL69g+PnAnMrxFcCJ1SIvwmcsy9llKRa6svthD2RJEn7JDO/B/y8TXgGcGexfSelZZhb4/Mz863MXE/p2+QpETGCYunmYgjCXW3Oab3WQuDU1l5KkiRJkmrHJJIkqRqOKuavoHg/sojvXoa50LpE8yi6uHQz0Lp0814iYnZErIyIlVu3bu2lW5EkSZIEJpEkSbVV1aWbM3NeZk7OzMnDhw/vYRElSZIkVeKcSJJ6z7Lrevd675vTu9dTLb0SESMyc0sxVO3VIt4CjC47rgnYXMSbKsTLz2kplm4+mL2Hz0nqC2wnJEkdsZ1oePZEktTnzZ07l4kTJ3LiiScyadIkHn/88bqU48orr+TYY4/lxBNP5KyzzuL111/fve+6665j3LhxHHPMMTz44IO749dccw2jR4/moIMO2uNaV1xxxe6VGiZMmMAhhxxSo7voNeXLLc9iz2WYXbpZUk3ZTkiSOmI70XUmkST1acuXL+e+++7jySef5Omnn2bp0qWMHj268xOrYNq0aTz77LM8/fTTTJgwgeuuK32TsnbtWubPn8+aNWtYvHgxl156Kbt2lVYrPuOMM3jiiSf2utaNN97I6tWrWb16NZdffjl/9Vd/VdN76Y52lm6+HpgWES8C04rPLt0sqeZsJyRJHbGd6B6TSJL6tC1btnDEEUcwZMgQAI444ghGjhzJihUreM973sO73vUupkyZwvbt29m1axdXXnklJ598MieeeCJf/vKXAXjkkUeYOnUqZ599Nsceeyznn38+rR1dVq1axXvf+15OOukkpk+fzpYtW9oty2mnncbgwaVRws3NzbS0lOaJvvfee5k5cyZDhgxh7NixjBs3bndF39zczIgRIzq8x3vuuYfzzqu0OnJjyMzzMnNEZu6XmU2ZeXtmbsvMUzNzfPH+87Lj52bmH2TmMZn5QFl8ZWaeUOz7eGtvo8x8MzPPycxxmTklM1+qx31K6ptsJyRJHbGd6B6TSJL6tNNOO41NmzYxYcIELr30Uh599FF+85vfcO655/KFL3yBH/7whyxdupShQ4dy++23c/DBB7NixQpWrFjBV77yFdavXw/AU089xec//3nWrl3LSy+9xA9+8AN27NjB5ZdfzsKFC1m1ahUf+9jHuOaaa7pUrq997Wt84AMfAODll1/e49uMpqYmXn755S5d5yc/+Qnr16/n/e9/fzefjCQJbCckSR2znegeJ9aW1KcddNBBrFq1iu9///ssW7aMc889l2uuuYYRI0Zw8sknAzBs2DAAHnroIZ5++mkWLlwIwBtvvMGLL77I/vvvz5QpU2hqKs3rPGnSJDZs2MAhhxzCs88+y7Rp0wDYtWtXp1l+KI2pHjx4MOeffz4AlabvKU3707n58+dz9tlnM2jQoC4dL0nak+2EJKkjthPdYxJJUp83aNAgpk6dytSpU/nDP/xDbrnlloqVamZy8803M3369D3ijzzyyO7uq63X27lzJ5nJxIkTWb58eZfLcuedd3Lffffx8MMP7y5DU1MTmzZt2n1MS0sLI0eO7NL15s+fzy233NLlny9J2pvthCSpI7YTXddpEikiRgN3Ab8H/BaYl5lfiIjDgG8BRwMbgA9n5mvFOXOAi4BdwCcy88EifhJwBzAUuB/4ZGZmRAwpfsZJwDbg3Mzc0Gt3Kak26rCE5vPPP8/b3vY2xo8fD8Dq1as57rjjWLx4MStWrODkk09m+/btDB06lOnTp3Prrbfy/ve/n/32248XXniBUaNGtXvtY445hq1bt7J8+XJOOeUUduzYwQsvvMDEiRMrHr948WJuuOEGHn30UQ488MDd8Q996EN85CMf4dOf/jSbN2/mxRdfZMqUKV26t9dee41TTjmlm09FkhqU7YTthCR1xHai4duJrvRE2gn8p8x8MiLeAayKiCXAfwAezszrI+JqSqvlXBURxwMzgYnASGBpREwoVt+5FZgNPEYpiXQ6pdV3LgJey8xxETETuAE4t9fuUlK/9Ytf/ILLL7+c119/ncGDBzNu3DjmzZvHhRdeyOWXX86vf/1rhg4dytKlS7n44ovZsGED7373u8lMhg8fzne+8512r73//vuzcOFCPvGJT/DGG2+wc+dOPvWpT7Vb6X/84x/nrbfe2t1dtbm5mdtuu42JEyfy4Q9/mOOPP57Bgwdzyy237O5O+pnPfIa7776bX/3qVzQ1NXHxxRdz7bXXAqUJ8GbOnNnlrqqSpL3ZTkhlll23d6wOv7RLjcR2onui0ti6Dk+IuBf4YvGamplbImIE8EhmHlP0QiIzryuOfxC4llJvpWWZeWwRP684/+9aj8nM5RExGPgpMDw7KNzkyZNz5cqV3bvbnqhU0YKVrQQ899xzHHfccfUuxoBW6c8gIlZl5uQ6Falh1KydqJIbl7zQ5WOvmDahiiWRes52ov5sJ9rX19uJHqlxEqm9tsx2S61sJ+qvu+1Et1Zni4ijgT8CHgeOyswtAMX7kcVho4BNZae1FLFRxXbb+B7nZOZO4A3g8O6UTZIkSZIkSdXT5Ym1I+Ig4F+BT2Xmv3fQHarSjuwg3tE5bcswm9JwOMaMGdNZkSWpKi677DJ+8IMf7BH75Cc/yYUXXlinEkmSGonthCSpI325nehSEiki9qOUQPqXzPx2EX4lIkaUDWd7tYi3AKPLTm8CNhfxpgrx8nNaiuFsBwM/b1uOzJwHzINS99OulF1SdWXmgJuLoVFWwenucGRJqgfbifqxnZDUF9hO1E9P2olOh7NF6U/zduC5zPxc2a5FwKxiexZwb1l8ZkQMiYixwHjgiWLI2/aIaC6ueUGbc1qvdTbw3Y7mQ5LUGA444AC2bdvmf1LrIDPZtm0bBxxwQL2LIkntsp2oH9sJSX2B7UT99LSd6EpPpD8B/gZ4JiJWF7H/DFwPLIiIi4CNwDlFQdZExAJgLaWV3S4rVmYDuAS4AxhKaVW2B4r47cA3ImIdpR5IM7t1F5LqoqmpiZaWFrZu3VrvogxIBxxwAE1NTZ0fKEl1YjtRX7YTkhqd7UR99aSd6DSJlJn/RuU5iwBObeecucDcCvGVwAkV4m9SJKEk9R377bcfY8eOrXcxJEkNynZCktQR24m+p8sTa0uSpH1Q42WVJUmSpN5mEkmSJEmSBiq/5JDUDZ1OrC1JkiRJkiSZRJIkSZIkSVKnTCJJkiRJkiSpUyaRJEmSJEmS1CmTSJIkSZKqJiK+FhGvRsSzZbHDImJJRLxYvB9atm9ORKyLiOcjYnpZ/KSIeKbYd1NERBEfEhHfKuKPR8TRNb1BSRpATCJJkiRJqqY7gNPbxK4GHs7M8cDDxWci4nhgJjCxOOdLETGoOOdWYDYwvni1XvMi4LXMHAfcCNxQtTuRpAHOJJIkSZKkqsnM7wE/bxOeAdxZbN8JnFkWn5+Zb2XmemAdMCUiRgDDMnN5ZiZwV5tzWq+1EDi1tZeSJKl3mUSSJEmSVGtHZeYWgOL9yCI+CthUdlxLERtVbLeN73FOZu4E3gAOr/RDI2J2RKyMiJVbt27tpVuRpIHDJJIkSZKkRlGpB1F2EO/onL2DmfMyc3JmTh4+fHgPiyhJA5dJJEmSJEm19koxRI3i/dUi3gKMLjuuCdhcxJsqxPc4JyIGAwez9/A5SVIvGFzvAvQHNy55Ya/YFdMm1KEkkiRJUp+wCJgFXF+831sWvzsiPgeMpDSB9hOZuSsitkdEM/A4cAFwc5trLQfOBr5bzJskSeplJpEkSZIkVU1E3ANMBY6IiBbgs5SSRwsi4iJgI3AOQGauiYgFwFpgJ3BZZu4qLnUJpZXehgIPFC+A24FvRMQ6Sj2QZtbgtiRpQDKJJElSvSy7bu/Y++bUvhySVEWZeV47u05t5/i5wNwK8ZXACRXib1IkoSRJ1eWcSJIkSZIkSeqUSSRJkiRJkiR1yiSSJEmSJEmSOmUSSZJUNRFxRUSsiYhnI+KeiDggIg6LiCUR8WLxfmjZ8XMiYl1EPB8R08viJ0XEM8W+myIi6nNHkiRJ0sBlEkmSVBURMQr4BDA5M08ABlFaMedq4OHMHA88XHwmIo4v9k8ETge+FBGDisvdCsymtNTz+GK/JEmSpBoyiSRJqqbBwNCIGAwcCGwGZgB3FvvvBM4stmcA8zPzrcxcD6wDpkTECGBYZi7PzATuKjtHkiRJUo2YRJIkVUVmvgz8M7AR2AK8kZkPAUdl5pbimC3AkcUpo4BNZZdoKWKjiu228b1ExOyIWBkRK7du3dqbtyNJkiQNeCaRJElVUcx1NAMYC4wE3h4RH+3olAqx7CC+dzBzXmZOzszJw4cP726RJUmSJHXAJJIkqVr+HFifmVszcwfwbeA9wCvFEDWK91eL41uA0WXnN1Ea/tZSbLeNS5IkSaohk0iSpGrZCDRHxIHFamqnAs8Bi4BZxTGzgHuL7UXAzIgYEhFjKU2g/UQx5G17RDQX17mg7BxJkiRJNTK43gWQJPVPmfl4RCwEngR2Ak8B84CDgAURcRGlRNM5xfFrImIBsLY4/rLM3FVc7hLgDmAo8EDxkiRJDeLGJS/UuwiSasAkkiSpajLzs8Bn24TfotQrqdLxc4G5FeIrgRN6vYCSJEmSuszhbJIkSZIkSeqUSSRJkiRJkiR1yiSSJEmSJEmSOmUSSZIkSZIkSZ0yiSRJkiRJkqROmUSSJEmSJElSp0wiSZIkSZIkqVMmkSRJkiRJktQpk0iSJEmSJEnqlEkkSZIkSZIkdarTJFJEfC0iXo2IZ8ti10bEyxGxunh9sGzfnIhYFxHPR8T0svhJEfFMse+miIgiPiQivlXEH4+Io3v5HiVJkiRJkrSPutIT6Q7g9ArxGzNzUvG6HyAijgdmAhOLc74UEYOK428FZgPji1frNS8CXsvMccCNwA09vBdJkiRJkiRVSadJpMz8HvDzLl5vBjA/M9/KzPXAOmBKRIwAhmXm8sxM4C7gzLJz7iy2FwKntvZSkiRJktR/RcQVEbEmIp6NiHsi4oCIOCwilkTEi8X7oWXHd2vUgySpdw3eh3M/HhEXACuB/5SZrwGjgMfKjmkpYjuK7bZxivdNAJm5MyLeAA4Hftb2B0bEbEq9mRgzZsw+FF2SJElSPUXEKOATwPGZ+euIWEBpVMPxwMOZeX1EXA1cDVzVZtTDSGBpREzIzF38btTDY8D9lEY9PFDzm1KX3Ljkhb1iV0ybUIeSSOqunk6sfSvwB8AkYAvwP4t4pYx/dhDv6Jy9g5nzMnNyZk4ePnx4twosSZIkqeEMBoZGxGDgQGAze45UuJM9RzB0d9SDJKkX9SiJlJmvZOauzPwt8BVgSrGrBRhddmgTpYagpdhuG9/jnKLxOJiuD5+TJEmS1Adl5svAPwMbKX0x/UZmPgQclZlbimO2AEcWp+wewVBoHd0wivZHPewhImZHxMqIWLl169bevB1JGhB6lEQqsv2tzgJaV25bBMwsVlwbS2kC7SeKyn97RDQX45MvAO4tO2dWsX028N3iGwRJkiRJ/VQx19EMYCyl4Wlvj4iPdnRKhVhnox72DDqyQZL2SadzIkXEPcBU4IiIaAE+C0yNiEmUKucNwN8BZOaaYizzWmAncFkxRhngEkorvQ2lND65dYzy7cA3ImIdpR5IM3vhviRJkiQ1tj8H1mfmVoCI+DbwHuCViBiRmVuKL69fLY7vyagHSVIv6jSJlJnnVQjf3sHxc4G5FeIrgRMqxN8EzumsHJIkSZL6lY1Ac0QcCPwaOJXSoj2/pDRS4frivXwEw90R8TlKPZdaRz3siojtEdEMPE5p1MPNNb0TSRog9mV1NkmSJEnqkcx8PCIWAk9SGsXwFDAPOAhYEBEXUUo0nVMc35NRD5KkXmQSSZIkSVJdZOZnKU2XUe4tSr2SKh3frVEPkqTe1aOJtSVJkiRJkjSwmESSJEmSJElSpxzOJkmSJEnq2LLr9o69b07tyyGpruyJJEmSJEmSpE6ZRJIkSZIkSVKnTCJJkiRJkiSpU86JJEmSJEmNptIcRFCbeYja+9lVdOOSFyrGr5g2ocYlkdQReyJJkiRJkiSpUyaRJElVExGHRMTCiPhRRDwXEadExGERsSQiXizeDy07fk5ErIuI5yNieln8pIh4pth3U0REfe5IkiRJGrhMIkmSqukLwOLMPBZ4F/AccDXwcGaOBx4uPhMRxwMzgYnA6cCXImJQcZ1bgdnA+OJ1ei1vQpIkSZJJJElSlUTEMODPgNsBMvM3mfk6MAO4szjsTuDMYnsGMD8z38rM9cA6YEpEjACGZebyzEzgrrJzJEmSJNWISSRJUrW8E9gKfD0inoqIr0bE24GjMnMLQPF+ZHH8KGBT2fktRWxUsd02LkmSJKmGTCJJkqplMPBu4NbM/CPglxRD19pRaZ6j7CC+9wUiZkfEyohYuXXr1u6WV5IkSVIHTCJJkqqlBWjJzMeLzwspJZVeKYaoUby/Wnb86LLzm4DNRbypQnwvmTkvMydn5uThw4f32o1IkiRJMokkSaqSzPwpsCkijilCpwJrgUXArCI2C7i32F4EzIyIIRExltIE2k8UQ962R0RzsSrbBWXnSJIkSaqRwfUugCSpX7sc+JeI2B94CbiQ0hcYCyLiImAjcA5AZq6JiAWUEk07gcsyc1dxnUuAO4ChwAPFS5IkSVINmUSSJFVNZq4GJlfYdWo7x88F5laIrwRO6NXCSZIkSeoWk0iSJKlh3LjkhYrxK6ZNqHFJJEmS1JZzIkmSJEmSJKlTJpEkSZIkSZLUKYezSZLURzVvnFfaWHb474Lvm1OfwkiSJKnfsyeSJEmSJEmSOmUSSZIkSZIkSZ0yiSRJkiRJkqROmUSSJEmSVBcRcUhELIyIH0XEcxFxSkQcFhFLIuLF4v3QsuPnRMS6iHg+IqaXxU+KiGeKfTdFRNTnjiSpfzOJJEmSJKlevgAszsxjgXcBzwFXAw9n5njg4eIzEXE8MBOYCJwOfCkiBhXXuRWYDYwvXqfX8iYkaaAwiSRJkiSp5iJiGPBnwO0AmfmbzHwdmAHcWRx2J3BmsT0DmJ+Zb2XmemAdMCUiRgDDMnN5ZiZwV9k5kqReZBJJkiRJUj28E9gKfD0inoqIr0bE24GjMnMLQPF+ZHH8KGBT2fktRWxUsd02LknqZSaRJEmSJNXDYODdwK2Z+UfALymGrrWj0jxH2UF87wtEzI6IlRGxcuvWrd0tryQNeCaRJEmSJNVDC9CSmY8XnxdSSiq9UgxRo3h/tez40WXnNwGbi3hThfheMnNeZk7OzMnDhw/vtRuRpIFicL0LIEmSOte8cV69iyBJvSozfxoRmyLimMx8HjgVWFu8ZgHXF+/3FqcsAu6OiM8BIylNoP1EZu6KiO0R0Qw8DlwA3Fzj25GkAcEkkiRJkqR6uRz4l4jYH3gJuJDSaIkFEXERsBE4ByAz10TEAkpJpp3AZZm5q7jOJcAdwFDggeIlSeplJpEkSZIk1UVmrgYmV9h1ajvHzwXmVoivBE7o1cJJkvZiEkmSJGkf3LjkhYrxK6ZNqHFJJKk2lr+0DYDHdlau/yT1X50mkSLia8BfAq9m5glF7DDgW8DRwAbgw5n5WrFvDnARsAv4RGY+WMRP4nddTO8HPpmZGRFDgLuAk4BtwLmZuaHX7lCSJEmSGsWy6/aOvW9O7cvRC9qbr++xMbNrXBJJtdKV1dnuAE5vE7saeDgzxwMPF5+JiOOBmcDE4pwvRcSg4pxbgdmUJsAbX3bNi4DXMnMccCNwQ09vRpIkSZIkSdXRaU+kzPxeRBzdJjwDmFps3wk8AlxVxOdn5lvA+ohYB0yJiA3AsMxcDhARdwFnUprwbgZwbXGthcAXIyIyM3t6U5IkqQv60bfhkqTGUamHkr2TpP6hp3MiHZWZWwAyc0tEHFnERwGPlR3XUsR2FNtt463nbCqutTMi3gAOB37W9odGxGxKvZkYM2ZMD4suSZIkSf1f2znbmjdu45R3Hl6n0kjqD7oynK07okIsO4h3dM7ewcx5mTk5MycPHz68h0WUJEmSJElSd/W0J9IrETGi6IU0Ani1iLcAo8uOawI2F/GmCvHyc1oiYjBwMPDzHpZLkqS+rRhi1rxxW50LIknqj1pXVmvLHkqSuqKnSaRFwCzg+uL93rL43RHxOWAkpQm0n8jMXRGxPSKagceBC4Cb21xrOXA28F3nQ5IkSR2qNJ8TOKeTJElSFXWaRIqIeyhNon1ERLQAn6WUPFoQERcBG4FzADJzTUQsANYCO4HLMnNXcalLKK30NpTShNoPFPHbgW8Uk3D/nNLqbpIkSZKkOmqv15Kkgasrq7Od186uU9s5fi4wt0J8JXBChfibFEkoSZIkSVIHXFlTUh319sTakiRJkiRJ6odMIkmSJEmSJKlTJpEkSZIkSZLUqZ6uziZJktTn3bjkhYrxK6ZNqHFJJEmSGp89kSRJVRURgyLiqYi4r/h8WEQsiYgXi/dDy46dExHrIuL5iJheFj8pIp4p9t0UEVGPe5EkSZIGMnsiSZKq7ZPAc8Cw4vPVwMOZeX1EXF18vioijgdmAhOBkcDSiJiQmbuAW4HZwGPA/cDpwAO1vQ31de31OpIkSVLXmESSJFVNRDQBfwHMBT5dhGcAU4vtO4FHgKuK+PzMfAtYHxHrgCkRsQEYlpnLi2veBZyJSSR1VbEcdvPGbbtDj42ZXa/SSJIk9VkmkSRJ1fR54DPAO8piR2XmFoDM3BIRRxbxUZR6GrVqKWI7iu228b1ExGxKPZYYM2ZMLxRfkqQ+atl1eyTPJak3OCeSJKkqIuIvgVczc1VXT6kQyw7iewcz52Xm5MycPHz48C7+WEmSJEldYU8kSVK1/AnwoYj4IHAAMCwivgm8EhEjil5II4BXi+NbgNFl5zcBm4t4U4W4JEmSpBqyJ5IkqSoyc05mNmXm0ZQmzP5uZn4UWATMKg6bBdxbbC8CZkbEkIgYC4wHniiGvm2PiOZiVbYLys6RJEmSVCP2RJIk1dr1wIKIuAjYCJwDkJlrImIBsBbYCVxWrMwGcAlwBzCU0oTaTqotSZIk1ZhJJElS1WXmI5RWYSMztwGntnPcXEorubWNrwROqF4JJUmSJHXG4WySJEmSJEnqlEkkSZIkSZIkdcokkiRJkqS6iYhBEfFURNxXfD4sIpZExIvF+6Flx86JiHUR8XxETC+LnxQRzxT7bioWYpAk9TLnRJIkSbW37Lp2dvx1TYshqSF8EngOGFZ8vhp4ODOvj4iri89XRcTxlFb7nAiMBJZGxIRiEYZbgdnAY8D9wOm4CIMk9Tp7IkmSJEmqi4hoAv4C+GpZeAZwZ7F9J3BmWXx+Zr6VmeuBdcCUiBgBDMvM5ZmZwF1l50iSepFJJEmSJEn18nngM8Bvy2JHZeYWgOL9yCI+CthUdlxLERtVbLeN7yUiZkfEyohYuXXr1l65AUkaSEwiSZIkSaq5iPhL4NXMXNXVUyrEsoP43sHMeZk5OTMnDx8+vIs/VpLUyjmRJElS76k019H75tS+HJL6gj8BPhQRHwQOAIZFxDeBVyJiRGZuKYaqvVoc3wKMLju/CdhcxJsqxOuv3fnfenhcFS1/aVu9iyCpDzCJJEmSqqtKvxxV+oXnsZ0vcMW0CVX5eZJ6V2bOAeYARMRU4O8z86MR8T+AWcD1xfu9xSmLgLsj4nOUJtYeDzyRmbsiYntENAOPAxcAN9fyXiRpoDCJJEmSJKmRXA8siIiLgI3AOQCZuSYiFgBrgZ3AZcXKbACXAHcAQymtyubKbJJUBSaRJEmSJNVVZj4CPFJsbwNObee4ucDcCvGVwAnVK6Hq5cYlL+wVs8epVD8mkSRJUlVVGnZ2yjsPr0NJuq7SLy3gLy6S1Aiso6X6MYkkSVIfV56keWzn7/5j7X+mJUmS1JtMIkmSNMC0foPbvPF3yae+1jOovOw90bxx3l6xx8bM7vzYZcVz6sKKc73xTbnDOCRJUiN5W70LIEmSJEmSpMZnTyRJkrR7SFz5cDiw10tb7T0nSZKkgcCeSJIkSZIkSeqUPZEkSeqn2puTZ1+vsbt30rLr9vn6bXVnriJJkiTVlkkkSZIkSdKA1xsLIkj9nUkkSZIkSVKf54qWUvWZRJIkSd3S+p/05o3b6lwSdYW/VElqVLUYwmzvIql3mUSSJEmSJFVVpYRRX2EyXvodV2eTJEmSJElSp/apJ1JEbAC2A7uAnZk5OSIOA74FHA1sAD6cma8Vx88BLiqO/0RmPljETwLuAIYC9wOfzMzcl7LVm90mJWkAq8KqZQNZX/72uqt6YyU9SZKkauuN4Wzvy8yflX2+Gng4M6+PiKuLz1dFxPHATGAiMBJYGhETMnMXcCswG3iMUhLpdOCBXiibJElqQMtfcj4lSf2MXyD0KSbvpZ6pxpxIM4CpxfadwCPAVUV8fma+BayPiHXAlKI307DMXA4QEXcBZ9LoSaSyRqJ8YtHenghOkqTuaK/Xju2TJPVflRIiLn4gqRr2dU6kBB6KiFUR0fq/06MycwtA8X5kER8FbCo7t6WIjSq228b3EhGzI2JlRKzcunXrPhZdkiRJkiRJXbWvPZH+JDM3R8SRwJKI+FEHx0aFWHYQ3zuYOQ+YBzB58uQ+PWeSJEm1NBDmFaq3Lj/jZYfD++ZUtzCSJElVsE9JpMzcXLy/GhH/G5gCvBIRIzJzS0SMAF4tDm8BRped3gRsLuJNFeKSJKnGKiVCHAonSaon2yapcfR4OFtEvD0i3tG6DZwGPAssAmYVh80C7i22FwEzI2JIRIwFxgNPFEPetkdEc0QEcEHZOZIkSZIkSWoA+zIn0lHAv0XED4EngP+bmYuB64FpEfEiMK34TGauARYAa4HFwGXFymwAlwBfBdYBP6bRJ9WWJHUqIkZHxLKIeC4i1kTEJ4v4YRGxJCJeLN4PLTtnTkSsi4jnI2J6WfykiHim2HdT8aWDJEmSpBrq8XC2zHwJeFeF+Dbg1HbOmQvMrRBfCZzQ07JIkhrSTuA/ZeaTRc/VVRGxBPgPwMOZeX1EXA1cDVwVEccDM4GJwEhgaURMKL5wuBWYDTwG3A+cjl841IzzKUmS+gKHvUnVt68Ta0uSVFExXLl1tc7tEfEcpdU3ZwBTi8PuBB4Briri8zPzLWB9RKwDpkTEBmBYZi4HiIi7gDMxiaQ+avlL23hs597LcUuSJDU6k0g1duOSyv9pvGLahBqXpI9bdt3eMVe6kRpWRBwN/BHwOHBUkWCiWIThyOKwUZR6GrVqKWI7iu228Uo/ZzalHkuMGTOmF+9AkiRJkkkkSVJVRcRBwL8Cn8rMf+9gOqNKO7KD+N7BzHnAPIDJkydXPEZqBA65kKSBzc4F6qv2ZWJtSZI6FBH7UUog/UtmfrsIvxIRI4r9I4BXi3gLMLrs9CZgcxFvqhCXJPVhLsAgSX2PPZEkSVVR/Af+duC5zPxc2a5FwCxKq3fOAu4ti98dEZ+jNLH2eOCJzNwVEdsjopnScLgLgJtrdBuSpOpxAQb1WZV6EtmLSAOBSSSpSmxYJP4E+BvgmYhYXcT+M6Xk0YKIuAjYCJwDkJlrImIBsJbSLxaXFb8YAFwC3AEMpfRLgb8YqNc14ip0Fcu07PB9ngfQYRRqBC7AIEl9j0kkSVJVZOa/UXk+I4BT2zlnLjC3QnwlcELvlU6S1EhcgEGS+gaTSGoMlVZbg15dcc2eQZIkSY3HBRgkqe9wYm1JkiRJdeECDJLUt9gTSb3K3j6SJNVJpV69vdijV+ptLsAgSX2PSSSpp2owBE+SJKkfcwEGSepjTCL1skqrqDw2ZnYdSlJDZcmU5o3bdm/3+/uWJKmBLH9p216xx3ZWXoVNagQuwCBJfY9JJO1Ws+V+2+vBI0l9jfWZJEmSBhCTSDVQqXcS9P+eOrvve9nh9S1Ib/AXRUmS9lCzL58kSVLDcHU2SZIkSZIkdcqeSA2iV1c168MTPredz6F1LocuPYs+fN+SJEmSJDU6k0iqaI8heOXD0UzISJIkSaqz9qYMkVRdJpGkCtr2DGveuI1T3lm9uZ261BPNnlaSJElSw3KuOA0EJpGkLqq0dDJQ1eSSJEkVueCDJHXJQF3kSKoWk0gNrCEz2b34n9b2kjK11N4z7o725nGSJEmSJKk/MYmkquuNZFGl4WXQ815AfiMhSZIkDVyVfh/wdwGpcyaR6siKa99VSlDtlVgqek+1Jp76uko9nxxnLUmqpt7ouStJA5V1qPoTk0h9UJcmYa5kgEzMvC89n6q9ysPePapKP6+5zXH1SCY25PBJSZIkSVLDMImkHvXUaYT5jCRJkiRJUu2YRGowPe0Js/z2rh/b3XmETBhJkiRJkiSTSFKDam8y8X25RiuHqEk94JLqkiRJGuBMIg1A9iyqrj658tuy6yomqRq6zJIkSZKkmjKJJPUDlXocNW+c1xATdkuSJEmS+geTSFKNVHvlN0l9k71DVU2V2p56faHQ49VlJakP6K361rpSjc4kkqR2NW+cB8vaTMT+vjm9cm3na5IkSVJf5JfDGshMIknqUNteEo/tfMFEjyRJkgaERkgYtfflayX+P13VZhJJ6uMaoWEDKq9c1Uu9liRJvadPLgAhSZIagkkkSd1WeSLvCqu72WtJkiRJkvoNk0iSqqo73W8lSfVTzXnw2uP8eJJUP9bB6gmTSJIktVVpeKYkSQ2g4qqeY2pfDnWuHitk9sYXuK4Qp46YRJJUNT2Zd8NGS5IaSJFQLR+yXKkO36O+L+/N1MOeTE4iK2kg6c4cp85fp3prmCRSRJwOfAEYBHw1M6+vc5EkNTgTTgNLf2gnKn57LPUx3VrQoaxXX2siqtMkVMFflNRd/aGdUP+3r4viNFJ96XC4gakhkkgRMQi4BZgGtAArImJRZq6tb8kkVUN3eihVnKMD2v1228asf7KdkPqP3l5V1HpfYDshtdXVZFNXj+tOD9F9HVLXnfrbNqD2GiKJBEwB1mXmSwARMR+YAVjpSwNId7/drrgiXDvfxDg0os/rU+2EPY6kfbevwzt6o97vTo9X25m661PthFQPXa1Xu9PbqRo9o6q1MI/1dO9olCTSKGBT2ecW4I/rVBZJfVhvfMO9/PauH9udRrI3fkkZwGwnJLVrX+v+9ur95jafHxszu6qrjvoLzj6xnZCqqDv1bG/3OO2JSvV62zq9Q5VGQlD6orCvDLeuVjsRmVmVC3erEBHnANMz8+Li898AUzLz8jbHzQZa/8SOAZ6vaUG75wjgZ/UuRJ147wOT9944fj8zh9e7EL2pn7YTtdZof08bgc9kbz6TvfXHZ2I7UdKI7UQj/31r1LJZru6xXN3TqOWC6pat3XaiUXoitQCjyz43AZvbHpSZ84D6pzW7ICJWZubkepejHrx3732gGcj3XkP9rp2oNf+e7s1nsjefyd58Jn1Gv2gnGvnvW6OWzXJ1j+XqnkYtF9SvbG+r9Q9sxwpgfESMjYj9gZnAojqXSZLUOGwnJEkdsZ2QpBpoiJ5ImbkzIj4OPEhpSc6vZeaaOhdLktQgbCckSR2xnZCk2miIJBJAZt4P3F/vcvSihu0mWwPe+8Dkvauq+mE7UWv+Pd2bz2RvPpO9+Uz6iH7STjTy37dGLZvl6h7L1T2NWi6oU9kaYmJtSZIkSZIkNbZGmRNJkiRJkiRJDcwkUjdExCERsTAifhQRz0XEKRFxWEQsiYgXi/dDy46fExHrIuL5iJheFj8pIp4p9t0UEVGfO+qaiDgmIlaXvf49Ij41EO4dICKuiIg1EfFsRNwTEQcMoHv/ZHHfayLiU0WsX957RHwtIl6NiGfLYr12rxExJCK+VcQfj4ija3qD6ncGet1cyUCur9szkOrxjljHq146qKuvjYiXy+IfLDunJv8OG7nObKdsjfDMGrJObadcNX9ejVzXdqdsEXF0RPy67NndVq2ytVOuc4o/y99GxOQ2x9e/fcpMX118AXcCFxfb+wOHAP8vcHURuxq4odg+HvghMAQYC/wYGFTsewI4BQjgAeAD9b63bjyDQcBPgd8fCPcOjALWA0OLzwuA/zBA7v0E4FngQErzpy0FxvfXewf+DHg38GxZrNfuFbgUuK3Yngl8q9737Kv/vAZa3dzOMxiw9XUHz2RA1eOdPAvreF91f7Wpq68F/r7CMTX5d9jIdWYHZav3M2vIOrWDctX8edHAdW03y3Z0+XFtrtOrZWunXMcBxwCPAJP35c9uX55Zey97InVRRAyj9Ad8O0Bm/iYzXwdmUEouUbyfWWzPAOZn5luZuR5YB0yJiBHAsMxcnqU/ybvKzukLTgV+nJk/YeDc+2BgaEQMplQ5b2Zg3PtxwGOZ+avM3Ak8CpxFP733zPwe8PM24d681/JrLQRO3ddvnaQyA7FurmSg1tftGVD1eEes49Ugyuvq9tTy32Ej15mVytaeWpWtUevU9srVnqqVq5Hr2m6WraJqlK1SuTLzucx8vsLhDdE+mUTquncCW4GvR8RTEfHViHg7cFRmbgEo3o8sjh8FbCo7v6WIjSq228b7ipnAPcV2v7/3zHwZ+GdgI7AFeCMzH2IA3DulbzT+LCIOj4gDgQ8CoxkY996qN+919zlFA/8GcHjVSq6BZkDVzZUM8Pq6PdbjHbOOV62V19UAH4+Ip4vhLK1DfGry77CR68wOygZ1fGY0bp3aXrmgvs+rVSPXte2VDWBs8Xv/oxHx/5T9/Hq2A43wzEwidcNgSt3Mbs3MPwJ+SanLW3sqZfeyg3jDi4j9gQ8B/6uzQyvE+uS9F5XtDErdBUcCb4+Ij3Z0SoVYn7z3zHwOuAFYAiym1HVyZwen9Jt774Ke3Gt/fA5qAAOxbq5kINfX7bEe7zHrePW6CnX1rcAfAJMoJUr+Z+uhFU7v9X+HjVxndlC2uj6zRq1TOyhXXZ9XFzRyXbsFGFP83v9p4O5iZFK9y9YQz8wkUte1AC2Z+XjxeSGlpNIrRfex1u5tr5YdP7rs/CZK3TBbiu228b7gA8CTmflK8Xkg3PufA+szc2tm7gC+DbyHgXHvZObtmfnuzPwzSt0sX2SA3HuhN+919zlF1+yD2btLrdQTA7FurmRA19ftsR7vkHW8ammPujozX8nMXZn5W+ArwJTiuFr9O2zkOrNi2RrgmTVsnVqpXI3wvAqNXNdWLFsxXGxbsb2K0txDE2pctkoa4ZmZROqqzPwpsCkijilCpwJrgUXArCI2C7i32F4EzCxmQx9LaXKzJ4puctsjorkYi3hB2TmN7jz27II7EO59I9AcEQcWZT4VeI6Bce9ExJHF+xjgryj9+Q+Iey/05r2WX+ts4LvFmGVpXw3EurmSAV1ft8d6vEPW8aqlPerq1l9cC2dRGpIEtft32Mh1ZsWyNcAza9g6tVK5GuF5lf28Rq1rK5YtIoZHxKBi+51F2V5qgHagEZ6Zq7N150WpK+BK4GngO8ChlMYTPkwpC/0wcFjZ8ddQylo+T9nM9sBkSv+Ifwx8EYh631sX7v1AYBtwcFlsoNz7PwA/Ksr9DUqz4Q+Ue/8+pWTpD4FT+/OfO6X/BGwBdlDK2F/Um/cKHECpC/s6SqsnvLPe9+yr778Gct3czvMYsPV1B89kwNTjnTwH63hfdXu1U1d/A3iG0u8Vi4ARZftq8u+wkevMdsrWCM+sIevUdspV8+dFA9e13Skb8NfAmuJ5PgmcUa2ytVOus4rtt4BXgAfr8czae7VeWJIkSZIkSWqXw9kkSZIkSZLUKZNIkiRJkiRJ6pRJJEmSJEmSJHXKJJIkSZIkSZI6ZRJJkiRJkiRJnTKJJEmSJEmSpE6ZRJIkSZIkSVKnTCJJkiRJkiSpU/8/jV2cAHKwaT4AAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 1440x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = 50\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(20,4))\n",
    "fig.suptitle('Band Histograms')\n",
    "\n",
    "for i, band in enumerate(['red', 'grn', 'blu']):\n",
    "    bins_list1 = rf_scene1.agg(\n",
    "        rf_agg_approx_histogram(band)['bins'].alias('bins')\n",
    "        ).collect()\n",
    "    values1 = [row['value'] for row in bins_list1[0].bins]\n",
    "    counts1 = [row['count'] for row in bins_list1[0].bins]\n",
    "\n",
    "    bins_list2 = rf_scene2.agg(\n",
    "        rf_agg_approx_histogram(band)['bins'].alias('bins')\n",
    "        ).collect()\n",
    "    values2 = [row['value'] for row in bins_list2[0].bins]\n",
    "    counts2 = [row['count'] for row in bins_list2[0].bins]\n",
    "\n",
    "    axs[i].hist(values1[:bins], weights=counts1[:bins], bins=bins, alpha = 0.5, label='Scene_2015')\n",
    "    axs[i].hist(values2[:bins], weights=counts2[:bins], bins=bins, alpha = 0.5, label='Scene_2017')\n",
    "    axs[i].title.set_text(band)\n",
    "    axs[i].legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we look at the average cell value for each band and compare across scenes, we can see that the average cell value is higher for the scene from 2015 for each band. The mean and distribution of cell values can be transformed using z-score normalization, which we talk about in the next section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>scene_2015_red_mean</th><th>scene_2015_grn_mean</th><th>scene_2015_blu_mean</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>7193.0</td><td>7709.0</td><td>8379.0</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/markdown": [
       "| scene_2015_red_mean | scene_2015_grn_mean | scene_2015_blu_mean |\n",
       "|---|---|---|\n",
       "| 7193.0 | 7709.0 | 8379.0 |"
      ],
      "text/plain": [
       "DataFrame[scene_2015_red_mean: double, scene_2015_grn_mean: double, scene_2015_blu_mean: double]"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_scene1.agg(rf_agg_stats('red').alias('red_stats'),\n",
    "              rf_agg_stats('grn').alias('grn_stats'),\n",
    "              rf_agg_stats('blu').alias('blu_stats')) \\\n",
    "         .select(F.round('red_stats.mean', 0).alias('scene_2015_red_mean'), \n",
    "                 F.round('grn_stats.mean', 0).alias('scene_2015_grn_mean'),\n",
    "                 F.round('blu_stats.mean', 0).alias('scene_2015_blu_mean'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>scene_2017_red_mean</th><th>scene_2017_grn_mean</th><th>scene_2017_blu_mean</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>6938.0</td><td>7397.0</td><td>7970.0</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/markdown": [
       "| scene_2017_red_mean | scene_2017_grn_mean | scene_2017_blu_mean |\n",
       "|---|---|---|\n",
       "| 6938.0 | 7397.0 | 7970.0 |"
      ],
      "text/plain": [
       "DataFrame[scene_2017_red_mean: double, scene_2017_grn_mean: double, scene_2017_blu_mean: double]"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_scene2.agg(rf_agg_stats('red').alias('red_stats'),\n",
    "              rf_agg_stats('grn').alias('grn_stats'),\n",
    "              rf_agg_stats('blu').alias('blu_stats')) \\\n",
    "         .select(F.round('red_stats.mean', 0).alias('scene_2017_red_mean'), \n",
    "                 F.round('grn_stats.mean', 0).alias('scene_2017_grn_mean'),\n",
    "                 F.round('blu_stats.mean', 0).alias('scene_2017_blu_mean'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also see that the range of cell values for each band is wider for the scene from 2017. The range of cell values can be transformed using min-max normalization, which we talk about in the last section."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>scene_2015_red_min</th><th>scene_2015_red_max</th><th>scene_2015_grn_min</th><th>scene_2015_grn_max</th><th>scene_2015_blu_min</th><th>scene_2015_blu_max</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>5752.0</td><td>14701.0</td><td>6393.0</td><td>15156.0</td><td>7532.0</td><td>14820.0</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/markdown": [
       "| scene_2015_red_min | scene_2015_red_max | scene_2015_grn_min | scene_2015_grn_max | scene_2015_blu_min | scene_2015_blu_max |\n",
       "|---|---|---|---|---|---|\n",
       "| 5752.0 | 14701.0 | 6393.0 | 15156.0 | 7532.0 | 14820.0 |"
      ],
      "text/plain": [
       "DataFrame[scene_2015_red_min: double, scene_2015_red_max: double, scene_2015_grn_min: double, scene_2015_grn_max: double, scene_2015_blu_min: double, scene_2015_blu_max: double]"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_scene1.agg(rf_agg_stats('red').alias('red_stats'),\n",
    "              rf_agg_stats('grn').alias('grn_stats'),\n",
    "              rf_agg_stats('blu').alias('blu_stats')) \\\n",
    "         .select(F.round('red_stats.min', 0).alias('scene_2015_red_min'), \n",
    "                 F.round('red_stats.max', 0).alias('scene_2015_red_max'), \n",
    "                 F.round('grn_stats.min', 0).alias('scene_2015_grn_min'),\n",
    "                 F.round('grn_stats.max', 0).alias('scene_2015_grn_max'),\n",
    "                 F.round('blu_stats.min', 0).alias('scene_2015_blu_min'),\n",
    "                 F.round('blu_stats.max', 0).alias('scene_2015_blu_max'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>scene_2017_red_min</th><th>scene_2017_red_max</th><th>scene_2017_grn_min</th><th>scene_2017_grn_max</th><th>scene_2017_blu_min</th><th>scene_2017_blu_max</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>5609.0</td><td>18498.0</td><td>6164.0</td><td>16269.0</td><td>7160.0</td><td>16785.0</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/markdown": [
       "| scene_2017_red_min | scene_2017_red_max | scene_2017_grn_min | scene_2017_grn_max | scene_2017_blu_min | scene_2017_blu_max |\n",
       "|---|---|---|---|---|---|\n",
       "| 5609.0 | 18498.0 | 6164.0 | 16269.0 | 7160.0 | 16785.0 |"
      ],
      "text/plain": [
       "DataFrame[scene_2017_red_min: double, scene_2017_red_max: double, scene_2017_grn_min: double, scene_2017_grn_max: double, scene_2017_blu_min: double, scene_2017_blu_max: double]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_scene2.agg(rf_agg_stats('red').alias('red_stats'),\n",
    "              rf_agg_stats('grn').alias('grn_stats'),\n",
    "              rf_agg_stats('blu').alias('blu_stats')) \\\n",
    "         .select(F.round('red_stats.min', 0).alias('scene_2017_red_min'), \n",
    "                 F.round('red_stats.max', 0).alias('scene_2017_red_max'),\n",
    "                 F.round('grn_stats.min', 0).alias('scene_2017_grn_min'),\n",
    "                 F.round('grn_stats.max', 0).alias('scene_2017_grn_max'),\n",
    "                 F.round('blu_stats.min', 0).alias('scene_2017_blu_min'),\n",
    "                 F.round('blu_stats.max', 0).alias('scene_2017_blu_max'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Z-Score Normalization\n",
    "\n",
    "Z-score normalization is a normalization technique that accounts for the mean offset and variance of cell values. EarthAI has a built-in function to perform z-score normalization called `rf_standardize`. It calculates the mean and variance of cell values, then subtracts the mean from each cell value and divides by the standard deviation. This technique centers the cell values at 0 and puts the range of cell values on roughly the same scale for each scene. It is good for transforming image scenes to similar distributions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "rf_scene1 = rf_scene1.withColumn('red_z_normalized', rf_standardize('red')) \\\n",
    "                     .withColumn('grn_z_normalized', rf_standardize('grn')) \\\n",
    "                     .withColumn('blu_z_normalized', rf_standardize('blu'))\n",
    "\n",
    "rf_scene2 = rf_scene2.withColumn('red_z_normalized', rf_standardize('red')) \\\n",
    "                     .withColumn('grn_z_normalized', rf_standardize('grn')) \\\n",
    "                     .withColumn('blu_z_normalized', rf_standardize('blu'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After performing z-score normalization, we see that the distribution of cell values in the band histograms are similar. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABJEAAAEVCAYAAABQV3TbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAABCjUlEQVR4nO3de5iddXno/e9NAiGIiISAIZMIlIRDEKKEdLDdNYohaEXUzSGKlSI0u8hJbBFSuivdrynQ+spZNAoGVIg0WuGlcggY0LITSAKIHCREiMlAhBgOjcohiff7x3omrMysOc+sw8z3c13rWmv9nsO6n2dmnt+s+/kdIjORJEmSJEmSOrNNrQOQJEmSJElS/TOJJEmSJEmSpC6ZRJIkSZIkSVKXTCJJkiRJkiSpSyaRJEmSJEmS1CWTSJIkSZIkSeqSSSRJktTQImJaRLSUvX8sIqb182fMi4gv9+c+eysi7omIU/qw/T9ExLf6MyZJkjQ0mESSJGmARMQJEfG7Co+MiH/qYJumiPhBRPw2Il6JiF9ExF9XOfS2Mf11EfM5bcpb+jtZ0x8yc1Jm3lOtzyvOz+ayn+/TEXFqtT6/TSwXRMR3K5RnROwDkJn/kpldJqH6mqySJEmDj0kkSZIGSGZ+LzN3LH8AnweeB77ZwWbfAdYA7wRGAZ8p1u83ETG8F5u9CJwbETvV6PPr3eKyn/ExwL9GxLtrHVS9GqS/A5IkDXomkSRJqpIiqXAJMDMz13aw2qHAvMz8fWZuysyHMvO2sn38eUT834h4OSLWtLZSioi3RcT1EbEuIn4dEf8YEdsUy/46Iu6LiEsi4kXggogYERFfiYjVEfF8RHw9IkZ2Ev4TwGLg7A6ObUREXBoRzxWPSyNiRLFsWtFq6dyI+A3w7aLFzL9HxHcjYkPR4mpiRMyOiBeKYzuibP8nRcQTxbpPR8T/6uQ8r4qIDxavXy5rIfT7okXOnsWyj0TEw8U6/zciDir/WUXEg8XnfR/YvpNzs5XMfLA4X/uX7e/fI+I3Reuyn0bEpLJl8yLiqoj4z+Lz7o+IPylbPj0ifllseyUQ3Y2lkvLWShGxffEzWF+ch6URsXtEzAH+B3Blce6uLNZ/b7HOK8Xze8v2u1dxbBsi4q7imFo/Z8/i3J8cEauBn3TzvHwtIm4rYrgvIt5R/G69VJyTd5etf25EPFt8/pMRcXhfzpMkSWrPJJIkSVUQETsDC4Avd9HVaglwVUTMjIjxbfYxHrgNuAIYDUwGHi4WXwG8DdgbeB+lFkwnlW3+p8DTwG7AHOBiYGKxj32AsUDFLnZl/jdwdkTsUmHZ+UBzsb+DganAP5YtfwewC6UWVrOKsqMotbx6O/AQcAel/03GAv8H+EbZ9i8AHwF2Ko7rkoh4Txfxkpk7l7UQugz4GfBsse21wP+i1OLrG8AtRTJsO+BHRWy7AP8O/M+uPqtVRBxK6dwuKyu+DZhA6fw/CHyvzWafBP6Z0rlYSelnRETsCvyA0rncFfgV8GfdjaUbTqT0ezOO0nn4W+DVzDyf0rk6vTh/pxc/9/8ELi/W/SrwnxExqtjXDcADxbILgL+q8Hnvo5Rcm1G87+q8HMebx/46pUTmg8X7BUUMRMS+wOnAoZn51mL/q3pzQiRJUsdMIkmSNMAiIoDrgEeBf+1i9WMpfXn/38AzRUuZQ4tlJwB3ZeaNmbkxM9dn5sMRMQw4HpidmRsycxXw/7L1l/jnMvOKzNwEvAb8DXB2Zr6YmRuAfwFmdhZYZj4M3AmcW2HxCcD/ycwXMnMdpYRI+ef/EfhSZr6ema8WZT/LzDuKmP6dUmLsoszcCMwH9iySb2Tmf2bmr7Lk3iKO/9FZvOUi4njgU8D/LPb/N8A3MvP+zNycmddRSlI0F49tgUuL87wAWNrFRzS3tnqilEj5DvBU68LMvLb42bxOKcFycES8rWz7H2bmA8W5+B6lZBzAh4HHM3NBEfelwG+6iOW4IpYtj07W3Ugp6bNPcR6WZ+Z/d7DuXwJPZeZ3ilZyNwK/BI4qEpyHAv+UmW9k5n8Bt1TYxwVFK7tXu3le/qOI6TXgP4DXMvP6zNwMfB9obYm0GRgBHBAR22bmqsz8VRfnSZIk9ZBJJEmSBt65wIHAiZmZrYVRmiWrtavV1wEy86XMPC8zJwG7U2pp9KMiETWOUkuUtnYFtgN+XVb2a0otelqtKXs9GtgBWF6WZLi9KO/KPwGnRsQ72pTvUeHz9yh7v65IBJQrH+vpVeC3RXKg9T3AjgAR8aGIWBIRLxbxfpjScXep6PJ0JfDxIsEFpRZRf9cm0TKuiHkP4Nnyn1WbY6tkSWurJ0qtriZRSswREcMi4qKI+FVE/DdvtpApj788MfSH1uMuYtnysytiKv9ZVnJTEcuWRyfrfodSC7D5UeqG+K8RsW0H67b9GcObv2d7AC9m5h/KllWKc0tZN89L29+Rtu93BMjMlZTGG7sAeCEi5kdE+e+fJEnqByaRJEkaQFGavex84JjMfLl8WTFLVuug23/bdtvM/C3wFUpf0Heh9AX8T9quB/yWUouSd5aVjQeeLd9dm/VfBSaVJRreViRAOpWZvwR+CPxDm0XPVfj85zr4/B6J0thKP6B0LnYvkiI/phtjA0XEaEotWE7PzIfKFq0B5rRJtuxQtK5ZC4wtEnflx9Mtmfl8Ee9RRdGngKOBD1LqOrZna3jd2N1aSsmt1uOJ8vd9VbS0+ufMPAB4L6Uug59pXdxm9bY/Y3jz92wtsEtE7FC2rFKc5fvsy3lpv+PMGzLzz4sYk1KXTUmS1I9MIkmSNEAiYgylblmfb5PA6GybiyPiwIgYHhFvBU4FVmbmekrdnD4YEccVy0dFxOSi9c5NwJyIeGtEvBP4AtBuqneAzPwjpdnhLomI3YrPHRsRMyqtX8E/UxqXaOeyshuBf4yI0cU4Pv/U0ef3wnaUuiqtAzZFxIeAIzrfZMsMYD8AvpeZ32+z+JvA30bEn0bJWyLiL4tzvhjYBJxZnOdPUBrjqVuKMYI+DjxWFL2VUle59ZRagP1Ld/dFaQyiSRHxieJ4zqTU0qlfRMT7I+JdRZfI/6aUjGxtDfY8pTG2Wv0YmBgRnyrOy/HAAcCtmflrSmNAXRAR20XEYbyZROtIX85L2+PYNyI+UCQcX6OUJN3cxWaSJKmHTCJJkjRw/oZSl7TLyrqtbdV9rYIdKLWceZnSQNjvBD4KkJmrKXXj+jvgRUpd3Q4utjsD+H2xzX9RGuT42k5iO5fSAM5Liq5EdwH7duegMvMZSt2g3lJW/GVKSYRHgF9QGvz4y93ZXzc+bwOl5MlNwEuUWrBUGm+nrSZK4yZ9vs25H5+Zyyj9fK4s9rkS+Ovi894APlG8f4nSeFM/7OKzDmvdP6WZ2dZR+pkAXE+p29ezwOOUBk/vlqI12rHARZSSLROA+7q7fTe8g9IA1f9NKe57eTP5dxlwTJRmQru8SGR+hNLv33rgi8BHihihNC7WYcWyL1Mas+j1Tj671+elghGUztFvKXUN3I32reUkSVIfxdbd/SVJkqS+i4jvA7/MzC/VOhZJktQ/bIkkSZKkPouIQyPiTyJim4g4ktJ4Rz+qcViSJKkfDa91AJIkSRoU3kGp298ooAU4tbtjgUmSpMZgdzZpAEXEBcA+mfnpWsciSZIkqTYiYhVwSmbe1aZ8GvDdzGyqQVhSj9mdTZIkSZIkSV0yiST1QDG9siRJvWI9IkmSGplJJKkLEbEqIs6NiEeA30fEn0fE/42IlyPi50UT1NZ194qIeyNiQ0QsBHatVdySpOqJiPdExEPF9f/fI+L7EfHliJgWES1FPfIb4NsRcUFE3BQR1xfrPxYRU2p9DJKkAXdoRDweES9FxLcjYvu2K0RERsQ+Ze/nRcSXqxum1DGTSFL3fBL4S2Bv4Gbgy8AuwN8DP4iI0cV6NwDLKSWP/h/gxOqHKkmqpojYDvgPYB6luuFG4ONlq7yjKH8nMKso+ygwH9gZuAW4sjrRSpJq6ARgBvAnwETgH2sbjtRzJpGk7rk8M9cAnwZ+nJk/zsw/ZuZCYBnw4YgYDxwK/O/MfD0zfwr8fzWMWZJUHc2UZry9PDM3ZuYPgQfKlv8R+FJRN7xalP1XUZdsBr4DHFzdkCVJNXBlZq7JzBeBOZRuVEsNxSSS1D1riud3AscWXdlejoiXgT8HxgB7AC9l5u/Ltvt1dcOUJNXAHsCzufWUt2vKXq/LzNfabPObstd/ALZ3vCRJGvTK64ZfU6o/pIZiEknqntYvBmuA72TmzmWPt2TmRcBa4O0R8Zay7cZXPVJJUrWtBcZGRJSVjSt7nUiStHXdMB54rsI6fwB2KHv/jgGNSOohk0hSz3wXOCoiZkTEsIjYvhg0tSkzf02pa9s/R8R2EfHnwFG1DVeSVAWLgc3A6RExPCKOBqbWOCZJUv05LSKaImIX4B+A71dY52HgU8V3jSOB91UzQKkrJpGkHijGRTqa0kV/HaWWSefw5t/Sp4A/BV4EvgRcX4MwJUlVlJlvAJ8ATgZepjR+3q3A6zUMS5JUf24A7gSeLh6VZl07i9KN6JcpDcT9oyrFJnVLbN19X5IkSX0VEfcDX8/Mb9c6FkmSpP5iSyRJkqQ+ioj3RcQ7iu5sJwIHAbfXOi5JkqT+5CwgkiRJfbcvcBOwI/Ar4JjMXFvbkCRJkvqX3dkkSZIkSZLUJbuzSZIkSZIkqUsN251t1113zT333LPWYUhS3Vm+fPlvM3N0reOoNesJSarMeqLEekKSKuusnmjYJNKee+7JsmXLah2GJNWdiPh1rWOoB9YTklSZ9USJ9YQkVdZZPWF3NkmSJEmSJHXJJJIkqU8i4tqIeCEiHi0r2yUiFkbEU8Xz28uWzY6IlRHxZETMKCs/JCJ+USy7PCKiKB8REd8vyu+PiD2reoCSJEmSAJNIkqS+mwcc2absPODuzJwA3F28JyIOAGYCk4ptvhYRw4ptrgZmAROKR+s+TwZeysx9gEuAiwfsSCRJkiR1qGHHRJJUexs3bqSlpYXXXnut1qEMSdtvvz1NTU1su+22NY0jM39aoXXQ0cC04vV1wD3AuUX5/Mx8HXgmIlYCUyNiFbBTZi4GiIjrgY8BtxXbXFDsawFwZUREZubAHJGk/mI9UVv1Uk9ExLXAR4AXMvPAouzfgKOAN4BfASdl5svFstmUbiBsBs7MzDuK8kMo3bgYCfwYOCszMyJGANcDhwDrgeMzc1W1jk9S71lP1FZv6gmTSJJ6raWlhbe+9a3sueeeFD2PVCWZyfr162lpaWGvvfaqdTiV7J6ZawEyc21E7FaUjwWWlK3XUpRtLF63LW/dZk2xr00R8QowCvht2w+NiFmUWjMxfvz4fjsYSb1jPVE7dVZPzAOupJToabUQmF1c1y8GZgPntmmxugdwV0RMzMzNvNlidQmlJNKRlG42bGmxGhEzKbVYPb4qRyapT6wnaqe39YTd2ST12muvvcaoUaO84NdARDBq1KhGvGtT6ZclOynvbJv2hZlzM3NKZk4ZPXrIz14t1Zz1RO3UUz2RmT8FXmxTdmdmbireLgGaitdbWqxm5jNAa4vVMRQtVouWqK0tVlu3ua54vQA4PPylkxqC9UTt9LaeMIkkqU+84NdOnZ/754t/+CmeXyjKW4BxZes1Ac8V5U0VyrfaJiKGA2+jzZcRSfWrzq9Vg1oDnfvPUmpRBGWtTwutLVPH0s0Wq0Bri9V2ImJWRCyLiGXr1q3rtwOQ1HsNdK0adHpz7k0iSZIGwi3AicXrE4Gby8pnFjOu7UVpAO0Hiq5vGyKiubh7/Jk227Tu6xjgJ46HJEmDQ0ScD2wCvtdaVGE1W6xKUp1wTCRJ/eaShSv6dX9nT5/Yr/vTwIiIGykNor1rRLQAXwIuAm6KiJOB1cCxAJn5WETcBDxO6UvDacU4FwCn8uaAqbfx5l3pa4DvFINwv0hprAxJDch6QuUi4kRKA24fXnZzoC8tVltssSo1NuuJ+mdLpL5YdGHlh6SqmjNnDpMmTeKggw5i8uTJ3H///TWJ45xzzmG//fbjoIMO4uMf/zgvv/zylmUXXngh++yzD/vuuy933HHHlvLzzz+fcePGseOOO261r3nz5jF69GgmT57M5MmT+da3vlWtw+ixzPxkZo7JzG0zsykzr8nM9Zl5eGZOKJ5fLFt/Tmb+SWbum5m3lZUvy8wDi2Wnt36hyMzXMvPYzNwnM6dm5tO1OM5BwXpLQ5T1RP2JiCMpzdr50cz8Q9kiW6yKSxauaPeQBpL1RPfZEmko6OgLwvtnVzcOaQAsXryYW2+9lQcffJARI0bw29/+ljfeeKMmsUyfPp0LL7yQ4cOHc+6553LhhRdy8cUX8/jjjzN//nwee+wxnnvuOT74wQ+yYsUKhg0bxlFHHcXpp5/OhAkT2u3v+OOP58orr6zBkUjS4GE9UXsdtFidDYwAFhZjcizJzL+1xaqkarOe6BlbIklqaGvXrmXXXXdlxIgRAOy6667sscceLF26lPe+970cfPDBTJ06lQ0bNrB582bOOeccDj30UA466CC+8Y1vAHDPPfcwbdo0jjnmGPbbbz9OOOEEWm9gLl++nPe9730ccsghzJgxg7Vr13YYyxFHHMHw4aXcfHNzMy0tpfE/b775ZmbOnMmIESPYa6+92GeffXjggQe2rDdmzJgBOz+SNNRZT9ReBy1W98nMcZk5uXj8bdn6tliVVDXWEz1jEklSQzviiCNYs2YNEydO5HOf+xz33nsvb7zxBscffzyXXXYZP//5z7nrrrsYOXIk11xzDW9729tYunQpS5cu5Zvf/CbPPPMMAA899BCXXnopjz/+OE8//TT33XcfGzdu5IwzzmDBggUsX76cz372s5x//vndiuvaa6/lQx/6EADPPvss48a9ObxDU1MTzz77bJf7+MEPfsBBBx3EMcccw5o1a7pcX5LUnvWEJKkz1hM9Y3c2SQ1txx13ZPny5fzsZz9j0aJFHH/88Zx//vmMGTOGQw89FICddtoJgDvvvJNHHnmEBQsWAPDKK6/w1FNPsd122zF16lSamkrjdU6ePJlVq1ax88478+ijjzJ9+nQANm/e3K0s/5w5cxg+fDgnnHACAJWGZehqOs2jjjqKT37yk4wYMYKvf/3rnHjiifzkJz/p5lmRJLWynpAkdcZ6omdMIg1irQPQNa9ev1X5YXuPqkU40oAZNmwY06ZNY9q0abzrXe/iqquuqnhRzUyuuOIKZsyYsVX5Pffcs6X5auv+Nm3aRGYyadIkFi9e3O1YrrvuOm699VbuvvvuLTE0NTVtlflvaWlhjz326HQ/o0a9+Xf6N3/zN5x77rndjkGStDXrCUlSZ6wnus8kkqR+U4spNJ988km22WabLQPJPfzww+y///7cfvvtLF26lEMPPZQNGzYwcuRIZsyYwdVXX80HPvABtt12W1asWMHYsWM73Pe+++7LunXrWLx4MYcddhgbN25kxYoVTJo0qeL6t99+OxdffDH33nsvO+yww5byj370o3zqU5/iC1/4As899xxPPfUUU6dO7fS41q5du+UuxS233ML+++/f01MjSXXHesJ6QpI6Yz1R//VEl0mkiLgW+AjwQmYeWJTtAnwf2BNYBRyXmS8Vy2YDJwObgTMz846i/BDenE3hx8BZmZkRMQK4HjgEWA8cn5mr+uXoJA16v/vd7zjjjDN4+eWXGT58OPvssw9z587lpJNO4owzzuDVV19l5MiR3HXXXZxyyimsWrWK97znPWQmo0eP5kc/+lGH+95uu+1YsGABZ555Jq+88gqbNm3i85//fIcX/dNPP53XX399S3PV5uZmvv71rzNp0iSOO+44DjjgAIYPH85VV13FsGHDAPjiF7/IDTfcwB/+8Aeampo45ZRTuOCCC7j88su55ZZbGD58OLvssgvz5s3r71MnSUOC9YQkqTPWEz0TlfrWbbVCxF8AvwOuL0si/SvwYmZeFBHnAW/PzHMj4gDgRmAqsAdwFzAxMzdHxAPAWcASSkmkyzPztoj4HHBQZv5tRMwEPp6Zx3cV+JQpU3LZsmW9Pe7+sejCyuXvn13dODrwZne2uVuVb+nOVidxqnE98cQT3vmssUo/g4hYnplTahRS3aiLeqLe1Hm9pcHHeqL2rCc6Zj1Rv1q/x5SrRQsVDTzridrraT3R5exsmflT4MU2xUcD1xWvrwM+VlY+PzNfz8xngJXA1IgYA+yUmYuLqTivb7NN674WAIdHVyNESZIkSZIkqap6OybS7pm5FiAz10bEbkX5WEotjVq1FGUbi9dty1u3WVPsa1NEvAKMAn7b9kMjYhYwC2D8+PG9DF2S+ua0007jvvvu26rsrLPO4qSTTqpRRNIQZusq1SHrCal+2cpJ9aCR64n+Hli7Ugui7KS8s23aF2bOBeZCqflpbwKUpL666qqrah2CJKmOWU9IkjrTyPVEl93ZOvB80UWN4vmForwFGFe2XhPwXFHeVKF8q20iYjjwNtp3n5MkSZIkSVIN9TaJdAtwYvH6RODmsvKZETEiIvYCJgAPFF3fNkREczHe0WfabNO6r2OAn2RXo31LkiRJkiSpqrrszhYRNwLTgF0jogX4EnARcFNEnAysBo4FyMzHIuIm4HFgE3BaZm4udnUqMA8YCdxWPACuAb4TESsptUCa2S9HJkmSJEmSpH7TZRIpMz/ZwaLDO1h/DjCnQvky4MAK5a9RJKEkNbiOBrjtLQfGlaTBxXpCktQZ64m619vubJJUN+bMmcOkSZM46KCDmDx5Mvfff39N4jjnnHPYb7/9OOigg/j4xz/Oyy+/vGXZhRdeyD777MO+++7LHXfcsaX8/PPPZ9y4cey4445b7evss89m8uTJTJ48mYkTJ7LzzjtX6SgkafCxnpAGr0sWrqj4kHrCeqL7TCJJamiLFy/m1ltv5cEHH+SRRx7hrrvuYty4cV1vOACmT5/Oo48+yiOPPMLEiRO58MLSnZTHH3+c+fPn89hjj3H77bfzuc99js2bSz19jzrqKB544IF2+7rkkkt4+OGHefjhhznjjDP4xCc+UdVjkaTBwnpCktQZ64meMYkkqaGtXbuWXXfdlREjRgCw6667sscee7B06VLe+973cvDBBzN16lQ2bNjA5s2bOeecczj00EM56KCD+MY3vgHAPffcw7Rp0zjmmGPYb7/9OOGEE2gd33/58uW8733v45BDDmHGjBmsXbu2w1iOOOIIhg8v9RJubm6mpaUFgJtvvpmZM2cyYsQI9tprL/bZZ58tF/rm5mbGjBnT6THeeOONfPKTHfUsliR1xnpCktQZ64meMYkkqaEdccQRrFmzhokTJ/K5z32Oe++9lzfeeIPjjz+eyy67jJ///OfcddddjBw5kmuuuYa3ve1tLF26lKVLl/LNb36TZ555BoCHHnqISy+9lMcff5ynn36a++67j40bN3LGGWewYMECli9fzmc/+1nOP//8bsV17bXX8qEPfQiAZ599dqu7GU1NTTz77LPd2s+vf/1rnnnmGT7wgQ/08MxIksB6QpLUOeuJnulyYG1Jqmc77rgjy5cv52c/+xmLFi3i+OOP5/zzz2fMmDEceuihAOy0004A3HnnnTzyyCMsWLAAgFdeeYWnnnqK7bbbjqlTp9LU1ATA5MmTWbVqFTvvvDOPPvoo06dPB2Dz5s1dZvmh1Kd6+PDhnHDCCQBb7kKUi4huHd/8+fM55phjGDZsWLfWlyRtzXpCGjqaV899882iUaVnB1ZWF6wnesYkkqSGN2zYMKZNm8a0adN417vexVVXXVXxopqZXHHFFcyYMWOr8nvuuWdL89XW/W3atInMZNKkSSxevLjbsVx33XXceuut3H333VtiaGpqYs2aNVvWaWlpYY899ujW/ubPn89VV13V7c+XaqHSAKZnT59Yg0ikyqwnJEmdsZ7oPpNIkvpPDe70PPnkk2yzzTZMmDABgIcffpj999+f22+/naVLl3LooYeyYcMGRo4cyYwZM7j66qv5wAc+wLbbbsuKFSsYO3Zsh/ved999WbduHYsXL+awww5j48aNrFixgkmTJlVc//bbb+fiiy/m3nvvZYcddthS/tGPfpRPfepTfOELX+C5557jqaeeYurUqd06tpdeeonDDjush2dFkuqU9YT1hCR1xnqi7usJk0iSGtrvfvc7zjjjDF5++WWGDx/OPvvsw9y5cznppJM444wzePXVVxk5ciR33XUXp5xyCqtWreI973kPmcno0aP50Y9+1OG+t9tuOxYsWMCZZ57JK6+8wqZNm/j85z/f4UX/9NNP5/XXX9/SXLW5uZmvf/3rTJo0ieOOO44DDjiA4cOHc9VVV21pTvrFL36RG264gT/84Q80NTVxyimncMEFFwClAfBmzpzZ7aaqkqT2rCckSZ2xnuiZqNS3rhFMmTIlly1bVtsgFl1YubxO+t22di/Yqm8wcNje9g9W/3jiiSfYf//9ax3GkFbpZxARyzNzSo1Cqht1UU/UmwGqt2rana3O6+Khznqi9qwnOmY9Ub96Uq+Ur1v+vaej7zx2wa4v1hO119N6wtnZJEmSJEmS1CW7s0lSD5122mncd999W5WdddZZnHTSSTWKSJJUT6wnJEmdaeR6wiSSpD7JzCE3FkO9zILTqN2RJQ0t1hO1Yz0hqRFYT9ROb+oJu7NJ6rXtt9+e9evX+09qDWQm69evZ/vtt691KJLUIeuJ2rGekNQIrCdqp7f1hC2RJPVaU1MTLS0trFu3rtahDEnbb789TU1NtQ5DkjpkPVFb9VJPRMS1wEeAFzLzwKJsF+D7wJ7AKuC4zHypWDYbOBnYDJyZmXcU5YcA84CRwI+BszIzI2IEcD1wCLAeOD4zV1Xp8DTAtgyWvWjUm4VOnjBoWE/UVm/qCZNIknpt2223Za+99qp1GJKkOmU9ocI84EpKiZ5W5wF3Z+ZFEXFe8f7ciDgAmAlMAvYA7oqIiZm5GbgamAUsoZREOhK4jVLC6aXM3CciZgIXA8dX5cgk9Yn1ROOxO5skSZKkAZOZPwVebFN8NHBd8fo64GNl5fMz8/XMfAZYCUyNiDHATpm5OEv9Xq5vs03rvhYAh8dQG2BFkqrElkiSJEmSqm33zFwLkJlrI2K3onwspZZGrVqKso3F67blrdusKfa1KSJeAUYBv237oRExi1JrJsaPH99vB6P6s/jp9QAs2bSixpFIg4stkSRJkiTVi0otiLKT8s62aV+YOTczp2TmlNGjR/cyREkaukwiSZIGTEScHRGPRcSjEXFjRGwfEbtExMKIeKp4fnvZ+rMjYmVEPBkRM8rKD4mIXxTLLrebgiQ1vOeLLmoUzy8U5S3AuLL1moDnivKmCuVbbRMRw4G30b77nCSpH5hEkiQNiIgYC5wJTClm4xlGabDU1sFUJwB3F+9pM5jqkcDXImJYsbvWwVQnFI8jq3gokqT+dwtwYvH6RODmsvKZETEiIvaidM1/oOj6tiEimosbCZ9ps03rvo4BfpLOFy5JA8IkkiRpIA0HRhZ3hnegdNe4PwdTlSTVuYi4EVgM7BsRLRFxMnARMD0ingKmF+/JzMeAm4DHgduB04qZ2QBOBb5FqX74FaWZ2QCuAUZFxErgCxQ3JyRJ/c+BtSVJAyIzn42IrwCrgVeBOzPzzojoz8FUJUl1LjM/2cGiwztYfw4wp0L5MuDACuWvAcf2JUZJUveYRJIkDYhirKOjgb2Al4F/j4hPd7ZJhbKuBlNt+5lDctadSxZWnnnm7OkTqxyJJEmSBjO7s0mSBsoHgWcyc11mbgR+CLyX/h1MdSvOuiNJkiQNHJNIkqSBshpojogdikFQDweeoH8HU5UkSZJUJXZnkyQNiMy8PyIWAA8Cm4CHgLnAjsBNxcCqqynGscjMxyKidTDVTbQfTHUeMJLSQKq3IUmSJKmqTCJJkgZMZn4J+FKb4tfpp8FUJUmSJFWP3dkkSZIkSZLUJZNIkiRJkiRJ6pJJJEmSJEmSJHWpT0mkiDg7Ih6LiEcj4saI2D4idomIhRHxVPH89rL1Z0fEyoh4MiJmlJUfEhG/KJZdXsy+I0mSJEmSpDrR6yRSRIwFzgSmZOaBwDBgJnAecHdmTgDuLt4TEQcUyycBRwJfi4hhxe6uBmZRms55QrFckiRJkiRJdaKv3dmGAyMjYjiwA/AccDRwXbH8OuBjxeujgfmZ+XpmPgOsBKZGxBhgp8xcnJkJXF+2jSRJkiRJkurA8N5umJnPRsRXgNXAq8CdmXlnROyemWuLddZGxG7FJmOBJWW7aCnKNhav25a3ExGzKLVYYvz48b0NXZIkdceiC9uXvX929eOQJDU26xNp0OhLd7a3U2pdtBewB/CWiPh0Z5tUKMtOytsXZs7NzCmZOWX06NE9DVmSJEmSJEm91JfubB8EnsnMdZm5Efgh8F7g+aKLGsXzC8X6LcC4su2bKHV/aylety2XJEmSJElSneh1dzZK3diaI2IHSt3ZDgeWAb8HTgQuKp5vLta/BbghIr5KqeXSBOCBzNwcERsiohm4H/gMcEUf4qpLlyxcUbH87OkTqxyJJEmSJElSz/VlTKT7I2IB8CCwCXgImAvsCNwUESdTSjQdW6z/WETcBDxerH9aZm4udncqMA8YCdxWPCRJkiRJklQn+tISicz8EvClNsWvU2qVVGn9OcCcCuXLgAP7EkstLH56fcXyw95f5UAkSZIkSZIGWF/GRJIkSZIkSdIQYRJJkiRJkiRJXTKJJEmSJEmSpC71aUwkSZKkSirNSuqMpJIkSY3NJJIkSZIkacBUmpBoyaYVPbq5UL6PJZva36iQVB12Z5MkSZIkSVKXTCJJkiRJkiSpSyaRJEmSJEmS1CWTSJIkSZIkSeqSSSRJkiRJNRERZ0fEYxHxaETcGBHbR8QuEbEwIp4qnt9etv7siFgZEU9GxIyy8kMi4hfFsssjImpzRJI0uJlEkiRJklR1ETEWOBOYkpkHAsOAmcB5wN2ZOQG4u3hPRBxQLJ8EHAl8LSKGFbu7GpgFTCgeR1bxUCRpyDCJJEmSJKlWhgMjI2I4sAPwHHA0cF2x/DrgY8Xro4H5mfl6Zj4DrASmRsQYYKfMXJyZCVxfto0kqR+ZRJIkSZJUdZn5LPAVYDWwFnglM+8Eds/MtcU6a4Hdik3GAmvKdtFSlI0tXrctbyciZkXEsohYtm7duv48HEkaEkwiSZIkSaq6Yqyjo4G9gD2At0TEpzvbpEJZdlLevjBzbmZOycwpo0eP7mnIkjTkmUSSJEmSVAsfBJ7JzHWZuRH4IfBe4PmiixrF8wvF+i3AuLLtmyh1f2spXrctlyT1M5NIkiRJkmphNdAcETsUs6kdDjwB3AKcWKxzInBz8foWYGZEjIiIvSgNoP1A0eVtQ0Q0F/v5TNk2kqR+NLzWAUiSJEkaejLz/ohYADwIbAIeAuYCOwI3RcTJlBJNxxbrPxYRNwGPF+uflpmbi92dCswDRgK3FQ9JUj8ziSRJkiSpJjLzS8CX2hS/TqlVUqX15wBzKpQvAw7s9wAlSVsxiSRJkga9SxauaFd29vSJNYhEkiSpcZlEkiRJkiT1i0pJ++YaxCFpYDiwtiRJkiRJkrpkEkmSNGAiYueIWBARv4yIJyLisIjYJSIWRsRTxfPby9afHRErI+LJiJhRVn5IRPyiWHZ5MfuOJEmSpCoyiSRJGkiXAbdn5n7AwZSmbj4PuDszJwB3F++JiAOAmcAk4EjgaxExrNjP1cAsStM5TyiWS5IkSaoik0iSpAERETsBfwFcA5CZb2Tmy8DRwHXFatcBHyteHw3Mz8zXM/MZYCUwNSLGADtl5uLMTOD6sm0kSZIkVYkDa0uSBsrewDrg2xFxMLAcOAvYPTPXAmTm2ojYrVh/LLCkbPuWomxj8bptuSRJUqeaV89tV7Zk/KwaRCINDrZEkiQNlOHAe4CrM/PdwO8puq51oNI4R9lJefsdRMyKiGURsWzdunU9jVeSJElSJ0wiSZIGSgvQkpn3F+8XUEoqPV90UaN4fqFs/XFl2zcBzxXlTRXK28nMuZk5JTOnjB49ut8ORJIkSZLd2SRJAyQzfxMRayJi38x8EjgceLx4nAhcVDzfXGxyC3BDRHwV2IPSANoPZObmiNgQEc3A/cBngCuqfDiSJGmIuWThiorlZ0+fWOVIpPphEkmSNJDOAL4XEdsBTwMnUWoFe1NEnAysBo4FyMzHIuImSkmmTcBpmbm52M+pwDxgJHBb8Rg6Fl1Yufz9s6sbhyRJkoY0k0iSpAGTmQ8DUyosOryD9ecAcyqULwMO7NfgJEmSJPVIn8ZEioidI2JBRPwyIp6IiMMiYpeIWBgRTxXPby9bf3ZErIyIJyNiRln5IRHxi2LZ5RFRaRBVSZIkSZIk1UhfB9a+DLg9M/cDDgaeoDTzzt2ZOQG4u3hPRBwAzAQmAUcCX4uIYcV+rgZmURr/YkKxXJIkSZIkSXWi10mkiNgJ+AvgGoDMfCMzXwaOBq4rVrsO+Fjx+mhgfma+npnPACuBqcXMPDtl5uLMTOD6sm0kSZIkSZJUB/oyJtLewDrg2xFxMLAcOAvYPTPXAmTm2ojYrVh/LLCkbPuWomxj8bptuSRJUu+0GYy8efV6AJaMn1WLaCRJkgaFvnRnGw68B7g6M98N/J6i61oHKo1zlJ2Ut99BxKyIWBYRy9atW9fTeCVJkiRJktRLfUkitQAtmXl/8X4BpaTS80UXNYrnF8rWH1e2fRPwXFHeVKG8ncycm5lTMnPK6NGj+xC6JEmSJEmSeqLXSaTM/A2wJiL2LYoOBx4HbgFOLMpOBG4uXt8CzIyIERGxF6UBtB8our5tiIjmYla2z5RtI0mSJEmSpDrQlzGRAM4AvhcR2wFPAydRSkzdFBEnA6uBYwEy87GIuIlSomkTcFpmbi72cyowDxgJ3FY8JElSvWkz1hCUxhtyrCFJkqTBr09JpMx8GJhSYdHhHaw/B5hToXwZcGBfYpEkSZIkSdLA6WtLJDWyCneTAXj/7OrGIUmSJEnF95PWGTUl1Z++DKwtSZIkSZKkIcIkkiRJkiRJkrpkdzYNnErd5ewqJ0mSJElSQ7IlkiRJkiRJkrpkSyRJkgapSxau2Op98+r1HLb3qBpFI0ntRcTOwLcozdScwGeBJ4HvA3sCq4DjMvOlYv3ZwMnAZuDMzLyjKD8EmAeMBH4MnJWZWb0jkaShwZZIkiRJkmrlMuD2zNwPOBh4AjgPuDszJwB3F++JiAOAmcAk4EjgaxExrNjP1cAsYELxOLKaByFJQ4UtkdQQ2t5Nb3X29IlVjkSSJEn9ISJ2Av4C+GuAzHwDeCMijgamFatdB9wDnAscDczPzNeBZyJiJTA1IlYBO2Xm4mK/1wMfA26r0qFI0pBhSyRJkiRJtbA3sA74dkQ8FBHfioi3ALtn5lqA4nm3Yv2xwJqy7VuKsrHF67bl7UTErIhYFhHL1q1b179HI0lDgC2RJEmSJNXCcOA9wBmZeX9EXEbRda0DUaEsOylvX5g5F5gLMGXKFMdMUq9U6iVhDwkNFSaRJElS9S26sHL5+2dXNw5JtdQCtGTm/cX7BZSSSM9HxJjMXBsRY4AXytYfV7Z9E/BcUd5UoVyqqHn13K0LFo2y/pG6ye5skiRJkqouM38DrImIfYuiw4HHgVuAE4uyE4Gbi9e3ADMjYkRE7EVpAO0Hii5vGyKiOSIC+EzZNpKkfmRLJDWkLXcPFrWZqto7CJIkSY3kDOB7EbEd8DRwEqUb3TdFxMnAauBYgMx8LCJuopRo2gSclpmbi/2cCswDRlIaUNtBtdUzFVvI/s+qhyHVO5NIkiRJkmoiMx8GplRYdHgH688B5lQoXwYc2K/BSZLasTubJEmSJEmSumQSSZIkSZIkSV2yO5skSVIZp26WJEmqzCSSJEkaVBY/vb7WIUiSJA1KdmeTJEmSJElSl0wiSZIkSZIkqUt2Z1O/aTuGRPPqUneCw/Ye1fudLrpwq3211bbLwpJNpRgcu0KSJEmSpP5lEkmSJNW/4qaCJEmSasckkiRJUlmSqrz165Lxs2oRjSRJUl1yTCRJkiRJkiR1ySSSJEmSJEmSumR3NknSgIqIYcAy4NnM/EhE7AJ8H9gTWAUcl5kvFevOBk4GNgNnZuYdRfkhwDxgJPBj4KzMzOoeiSRJGkqaV89tV2Y3Zw11tkSSJA20s4Anyt6fB9ydmROAu4v3RMQBwExgEnAk8LUiAQVwNTALmFA8jqxO6JIkSZJa2RJJA27x02UDlG5aseX12dMn1iIcSVUUEU3AXwJzgC8UxUcD04rX1wH3AOcW5fMz83XgmYhYCUyNiFXATpm5uNjn9cDHgNuqchCSJEmSAFsiSZIG1qXAF4E/lpXtnplrAYrn3YryscCasvVairKxxeu25e1ExKyIWBYRy9atW9cvByBJkiSpxCSSJGlARMRHgBcyc3l3N6lQlp2Uty/MnJuZUzJzyujRo7v5sZIkSZK6w+5s6rtFFwLQvHp9FytKGmL+DPhoRHwY2B7YKSK+CzwfEWMyc21EjAFeKNZvAcaVbd8EPFeUN1UolyRJklRFfU4iOeuOJKmSzJwNzAaIiGnA32fmpyPi34ATgYuK55uLTW4BboiIrwJ7UBpA+4HM3BwRGyKiGbgf+AxwRTWPRV2rNIMNOIuNJEnSYNIf3dmcdUeS1BMXAdMj4ilgevGezHwMuAl4HLgdOC0zNxfbnAp8C1gJ/AoH1ZYkSZKqrk8tkZx1pwNF965yzavXezdW0pCVmfdQqg/IzPXA4R2sN4dSndK2fBlw4MBFOMRVqLckSVIfVKpb3z+7+nFI/ayvLZEuxVl3JEmSJEmSBr1eJ5GcdUeSJEmSJGno6Et3NmfdqROXLFxR6xD6VfnxVGPGt47O39nTJw74Z0uS2rBrnSRJUt3qdUukzJydmU2ZuSelAbN/kpmfpjS7zonFam1n3ZkZESMiYi/enHVnLbAhIpojIijNunMzkiRJkiRJqht9Gli7AxcBN0XEycBq4FgozboTEa2z7myi/aw784CRlAbUbtxBtSVJUp8sfrpyK9TD9h5V5Ug6sejCinEedvJXahCMJElSdfRLEslZd7bW0T+/kiRJkqS+aV49t9YhSEPWQLREUp2reOe0Snd3t7rgLyr7zIGe7rKjMTacZlOSJKmmImIYsAx4NjM/EhG7AN8H9gRWAcdl5kvFurOBk4HNwJmZeUdRfghv9mz4MXBWZlacrEeS1HsmkTRkVEqeLdk0uAYllyRJakBnAU8AOxXvzwPuzsyLIuK84v25EXEApbFYJwF7AHdFxMRiiIyrgVnAEkpJpCNxiAxJ6ncmkcrZWkWSJEmqmohoAv6S0pAXXyiKjwamFa+vozRsxrlF+fzMfB14JiJWAlMjYhWwU2YuLvZ5PfAxTCJJUr8ziSRJkiSpVi4Fvgi8taxs92IGZzJzbUTsVpSPpdTSqFVLUbaxeN22vJ2ImEWpxRLjx4/vh/DVW82r5249vIWkhmASScDWXb3Ku3idPX1iLcKRJEnSIBcRHwFeyMzlETGtO5tUKMtOytsXZs4F5gJMmTLFMZPUL0yIaSgxiSRJktSBLRNC+OVAGgh/Bnw0Ij4MbA/sFBHfBZ6PiDFFK6QxwAvF+i3AuLLtm4DnivKmCuWSpH5mEkmSJDWsSpMmSGoMmTkbmA1QtET6+8z8dET8G3AicFHxfHOxyS3ADRHxVUoDa08AHsjMzRGxISKagfuBzwBXVPNYNHRsNdu0NASZRKpTlyxsP2uYXcskSbVkwkZSlVwE3BQRJwOrgWMBMvOxiLgJeBzYBJxWzMwGcCowDxhJaUBtB9WWpAFgEkmSJDUEk1jS4JWZ91CahY3MXA8c3sF6cyjN5Na2fBlw4MBFKEkCk0h1rV1TydbxGN4/u/rBNAjHrpCk+lXeyrZ5dSkhdNjetbteV4pHktQ3dveSBjeTSJIkacD4ZUKSJGnw2KbWAUiSJEmSJKn+2RJJkiS14/hDkiRJasskUiNadOFWb1vHcVgyflYtopEk1UpRH7Qdz8f6oDF0dyZWZ2yVpMZVflNmyaY3r+dex9WoTCJJkqSascWTJElS4zCJJPVEm1ZgWzhjnqQGZiJHkiRJ3WESSZIkSZKkPvCGjIYKk0hqZ6vpmBeNevO1rW0kSZIkSRqyTCLVgwpdpNoOklorDgQnSZIkSZIAtql1AJIkSZIkSap/tkSSJEkaaN1tddzajdwu5JIkqQ6ZRJIkaQhx4E9JkiT1lkmkKtpqwOpW5QNXS5IkwVYtl8pbLC0ZP6sW0UiSJAGOiSRJkiRJkqRusCWSVEHz6rm2EpMkSZIkqYxJJEmSJEmSBthWw5s4kYIalEkk9colC1dseV1xdhlJkiRJUreVf8cqd/b0iVWOROqYYyJJkgZERIyLiEUR8UREPBYRZxXlu0TEwoh4qnh+e9k2syNiZUQ8GREzysoPiYhfFMsuj4ioxTFJkiRJQ5ktkSRJA2UT8HeZ+WBEvBVYHhELgb8G7s7MiyLiPOA84NyIOACYCUwC9gDuioiJmbkZuBqYBSwBfgwcCdxW9SOSeuGShStstStJkgYFk0jqtq368EpSFzJzLbC2eL0hIp4AxgJHA9OK1a4D7gHOLcrnZ+brwDMRsRKYGhGrgJ0yczFARFwPfAyTSJIkSVJV9TqJFBHjgOuBdwB/BOZm5mURsQvwfWBPYBVwXGa+VGwzGzgZ2AycmZl3FOWHAPOAkZTuMJ+Vmdnb2CRJ9SUi9gTeDdwP7F4kmMjMtRGxW7HaWEotjVq1FGUbi9dtyyt9zixKLZYYP358Px7BwKs0DoJjIDSejsazkCRJGgz60hLJbgqSpC5FxI7AD4DPZ+Z/dzKcUaUF2Ul5+8LMucBcgClTpngzQkPHogsrlzvrjyRJ6ke9TiI1ejeFind87dyn3vKfd6miiNiWUgLpe5n5w6L4+YgYU7RCGgO8UJS3AOPKNm8CnivKmyqUS5IkSaqifkmb2E1BQ8nip9sPjnrY3qN6thOTThoCihnUrgGeyMyvli26BTgRuKh4vrms/IaI+CqlFqsTgAcyc3NEbIiIZkr1zGeAK6p0GFJVVapj+sT6RpIk9aM+J5HspiD1XE8SUR2Nr+FYKWoAfwb8FfCLiHi4KPsHSsmjmyLiZGA1cCxAZj4WETcBj1PqMn1a0eUZ4FTeHDvvNuzyLEkNzzFWJanx9CmJZDcFSVJHMvO/qHyjAODwDraZA8ypUL4MOLD/opMa05aZUhf1sAWsVJ8cY1WSGkxfZmezm4IkSVKZLUmeAdZRt7ced6+WaqjRx1iVBkq7uqT1xoFdkVUH+tISyW4KkiRJkvqsWmOsSpL6pi+zs9lNQZIkSVKfVHOMVSfqkaS+cVJ7qQP9PkOOJEmStlLtMVadqEeS+sYkkjSAOppZrbnKcUiSJNUbx1iVpMZjEkkaQNUaYFWSJKkBOcaqJDUYk0hlOpzp5P1VDkQNx65vkgZCa2vG5tVvXmOcfUvSYOEYq5LUeEwiSQ2gYoumRaOc5lOS1DuLLqxcbr0iSXWn9Yb1kk1bD5Vx9vSJtQhHQ5xJJKmeVPinvrwFgiRJkiRJtWISSZIkSZLUI04gIw1N29Q6AEmSJEmSJNU/k0iSJEmSJEnqkt3ZJEmSJEmqhbIxUR0LVY3AJJJUJ1pnXegRZ9eRJEmSJFWJSSRpkFn89Hqn/5SkIarSDYnD9h5Vg0gkSdJgZBJJGgI6mj3D5JIkSZLUGJpXz926YNEoeyCo6hxYW5IkSZIkSV2yJZIkSZIkSYOAPRA00EwiSZIkDWKOkyRJkvqLSSRpEGrXX7qwZPysKkciSZIkacC0ma25efV6/+fXgDKJJDWoSneWJUkaEG2+pGzhgK6SJA0pJpEkSZKGGLu4SeoxW7xIwiSSJEmSKCWWlmzaekBWB2KVJEnltql1AJIkSZIkSap/tkSShpB2A24vKrouOKaFJKk3HCtJkmrGMVJVCyaRJEmSJEk91tGMwGoQiy5sl4haMn6WXZnVKZNI3dHRXTZJkiRJkupI5eTeV6oehwYnk0iSes7uC5IkSZI05JhEkmRSSJJU0SULS7O1Na/eurvDYXuPqkU4kqQqaL32d5fd34YWk0iSJNUJx5ZQo6g0mOtWiSVvTkhSfennIVq2+p/FyXqGFJNI0hDW0YwO3mGWJPWHdgO2blrhHWtJkhqYSSRJkiRVTaVuEiaWpD6o1MKkJy1C+rq9BpVKraKXjJ9Vg0hUr0wiSWqn9c7xkk1v/qPfp3/w7dYgSeqM9YTUv7roumRrdPVEV93tK313gJ59f+hoHCZvMtSfukkiRcSRwGXAMOBbmXlRjUOSVKb8wl4+wKr/bKharCekxtfRF5HFHax/2PsHLhYNPtYTfdfvNxI1pLS9xi++ZoBaMdl6rqbqIokUEcOAq4DpQAuwNCJuyczHaxuZpK6U38lqe/ehVYdJp6ICaDdmRlHZ+E+LWllPSENUhS8Ki59e3+5LSaX6wm5zQ4v1hFSf+j5pyFfalVRqSdfR9xDw2t/f6iKJBEwFVmbm0wARMR84GvCiLw0yHTWfLtda2Sy+pvLydq2fOrrz0E/dI2xeWxcGVT2x+Jq/b1fWXIM4pHrXUZ1R6W53u3Uq7a/NepXukHf0hWfJ+Fm1ue7b1a+7BlU90aFOuql1538sqdFU+j+80vW9s2RV+bW/4s3qbrRs8sbEm+oliTQWWFP2vgX40xrFIqnOtfsn6en2X8g71cP1O/xyv6hyV74up77uDr8ctFXf9UQ/T5srqTp6coe8efXciskqu3XXjfquJ8p10MJOUvV0NFTHFm2+LzTT/sZDRzea68VAJbkiMwdkxz0KIuJYYEZmnlK8/ytgamae0Wa9WUDrT25f4MmqBtr/dgV+W+sg+pnH1Bg8psbQ22N6Z2aO7u9gaqnO64lG/t1r1NiNu7qMu7qqEbf1RMlg+D7RG436t9Ffhvrxg+cAPAddHX+H9US9tERqAcaVvW8Cnmu7UmbOBfraqbJuRMSyzJxS6zj6k8fUGDymxjAYj6kP6raeaOSfU6PGbtzVZdzV1ahx14G6rSfqzVD/HRvqxw+eA/Ac9OX4t+nvYHppKTAhIvaKiO2AmcAtNY5JklQ/rCckSZ2xnpCkKqiLlkiZuSkiTgfuoDQl57WZ+ViNw5Ik1QnrCUlSZ6wnJKk66iKJBJCZPwZ+XOs4qmwwNqX1mBqDx9QYBuMx9Vod1xON/HNq1NiNu7qMu7oaNe6aq+N6ot4M9d+xoX784DkAz0Gvj78uBtaWJEmSJElSfauXMZEkSZIkSZJUx0wi1VhEHBsRj0XEHyOioUeHj4gjI+LJiFgZEefVOp6+iohrI+KFiHi01rH0l4gYFxGLIuKJ4vfurFrH1FcRsX1EPBARPy+O6Z9rHVN/iIhhEfFQRNxa61jUtYj4t4j4ZUQ8EhH/ERE71zqm7mi0OqhR65lGrE8atb5o9DrBa7+qoVHrrP7QqPVIf2jU6/pAGOrX2ojYOSIWFNeBJyLisJ5sbxKp9h4FPgH8tNaB9EVEDAOuAj4EHAB8MiIOqG1UfTYPOLLWQfSzTcDfZeb+QDNw2iD4Ob0OfCAzDwYmA0dGRHNtQ+oXZwFP1DoIddtC4MDMPAhYAcyucTzd1TB1UIPXM/NovPqkUeuLRq8TvParGhq1zuqTBq9H+kOjXtcHwlC/1l4G3J6Z+wEH08NzYRKpxjLzicx8stZx9IOpwMrMfDoz3wDmA0fXOKY+ycyfAi/WOo7+lJlrM/PB4vUGSheMsbWNqm+y5HfF222LR0MP9hYRTcBfAt+qdSzqnsy8MzM3FW+XAE21jKe7GqwOath6phHrk0atLxq5TvDar2pp1DqrHzRsPdIfGvW63t+G+rU2InYC/gK4BiAz38jMl3uyD5NI6i9jgTVl71sYghelRhIRewLvBu6vcSh9VjRJfRh4AViYmY1+TJcCXwT+WOM41DufBW6rdRCDkPVMjTRafdHAdcKleO1X9Q2lOst6pNBo1/V+dilD+1q7N7AO+HbRpe9bEfGWnuzAJFIVRMRdEfFohcdgynxHhbKGuPM3FEXEjsAPgM9n5n/XOp6+yszNmTmZ0p20qRFxYI1D6rWI+AjwQmYur3Us2lp3ruURcT6l5uLfq12kWxtEdZD1TA00Yn3RiHWC1371t0atswaY9QiNeV3vL15rARgOvAe4OjPfDfwe6NH4YMMHIiptLTM/WOsYqqAFGFf2vgl4rkaxqBMRsS2liuN7mfnDWsfTnzLz5Yi4h9LYIw0zgG0bfwZ8NCI+DGwP7BQR383MT9c4riGvq2t5RJwIfAQ4PDPr5p/SQVQHWc9UWaPXFw1WJ3jtV79q1DprgA35eqTRr+v9wGtt6e+gpayV7gJ6mESyJZL6y1JgQkTsFRHbATOBW2ock9qIiKDU//WJzPxqrePpDxExunVWkYgYCXwQ+GVNg+qDzJydmU2ZuSelv6OfDLGKrSFFxJHAucBHM/MPtY5nkLKeqaJGrS8atU7w2q9qGsJ11pCuRxr1ut6fvNZCZv4GWBMR+xZFhwOP92QfJpFqLCI+HhEtwGHAf0bEHbWOqTeKwflOB+6gNEjbTZn5WG2j6puIuBFYDOwbES0RcXKtY+oHfwb8FfCBiHi4eHy41kH10RhgUUQ8Qumfg4WZOSSn61RNXQm8FVhY/F19vdYBdUcj1UGNXM80aH3SqPWFdYLUtYass/qqkeuRftKo13X1vzOA7xV15WTgX3qycQyd1ouSJEmSJEnqLVsiSZIkSZIkqUsmkSRJkiRJktQlk0iSJEmSJEnqkkkkSZIkSZIkdckkkiRJkiRJkrpkEkmSJEmSJEldMokkSZIkSZKkLplEkiRJkiRJUpf+f5GASGDYj+MuAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "bins = 50\n",
    "\n",
    "fig, axs = plt.subplots(1, 3, figsize=(20,4))\n",
    "fig.suptitle('Z-Score Normalized Band Histograms')\n",
    "\n",
    "for i, band in enumerate(['red', 'grn', 'blu']):\n",
    "    bins_list1 = rf_scene1.agg(\n",
    "        rf_agg_approx_histogram(band + '_z_normalized')['bins'].alias('bins')\n",
    "        ).collect()\n",
    "    values1 = [row['value'] for row in bins_list1[0].bins]\n",
    "    counts1 = [row['count'] for row in bins_list1[0].bins]\n",
    "\n",
    "    bins_list2 = rf_scene2.agg(\n",
    "        rf_agg_approx_histogram(band + '_z_normalized')['bins'].alias('bins')\n",
    "        ).collect()\n",
    "    values2 = [row['value'] for row in bins_list2[0].bins]\n",
    "    counts2 = [row['count'] for row in bins_list2[0].bins]\n",
    "\n",
    "    axs[i].hist(values1[:bins], weights=counts1[:bins], bins = bins, alpha = 0.5, label='Scene_2015')\n",
    "    axs[i].hist(values2[:bins], weights=counts2[:bins], bins = bins, alpha = 0.5, label='Scene_2017')\n",
    "    axs[i].title.set_text(band)\n",
    "    axs[i].legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The mean cell value of each band is zero for both scenes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>scene_2015_red_mean</th><th>scene_2015_grn_mean</th><th>scene_2015_blu_mean</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/markdown": [
       "| scene_2015_red_mean | scene_2015_grn_mean | scene_2015_blu_mean |\n",
       "|---|---|---|\n",
       "| 0.0 | 0.0 | 0.0 |"
      ],
      "text/plain": [
       "DataFrame[scene_2015_red_mean: double, scene_2015_grn_mean: double, scene_2015_blu_mean: double]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_scene1.agg(rf_agg_stats('red_z_normalized').alias('red_stats'),\n",
    "              rf_agg_stats('grn_z_normalized').alias('grn_stats'),\n",
    "              rf_agg_stats('blu_z_normalized').alias('blu_stats')) \\\n",
    "         .select(F.round('red_stats.mean', 0).alias('scene_2015_red_mean'), \n",
    "                 F.round('grn_stats.mean', 0).alias('scene_2015_grn_mean'),\n",
    "                 F.round('blu_stats.mean', 0).alias('scene_2015_blu_mean'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>scene_2017_red_mean</th><th>scene_2017_grn_mean</th><th>scene_2017_blu_mean</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>0.0</td><td>0.0</td><td>0.0</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/markdown": [
       "| scene_2017_red_mean | scene_2017_grn_mean | scene_2017_blu_mean |\n",
       "|---|---|---|\n",
       "| 0.0 | 0.0 | 0.0 |"
      ],
      "text/plain": [
       "DataFrame[scene_2017_red_mean: double, scene_2017_grn_mean: double, scene_2017_blu_mean: double]"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_scene2.agg(rf_agg_stats('red_z_normalized').alias('red_stats'),\n",
    "              rf_agg_stats('grn_z_normalized').alias('grn_stats'),\n",
    "              rf_agg_stats('blu_z_normalized').alias('blu_stats')) \\\n",
    "         .select(F.round('red_stats.mean', 0).alias('scene_2017_red_mean'), \n",
    "                 F.round('grn_stats.mean', 0).alias('scene_2017_grn_mean'),\n",
    "                 F.round('blu_stats.mean', 0).alias('scene_2017_blu_mean'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Min-Max Normalization\n",
    "\n",
    "Min-max normalization is a normalization technique that puts all cell values on the same scale by subtracting each value by the minimum cell value and dividing this result by the difference between the maximum and minimum cell values. The minimum value of each band is transformed to 0 and the maximum value of each band is transformed to 1. Unfortunately, min-max normalization does not handle outliers well since the calculation is based on the minimum and maximum cell values, but it does guarantee that all normalized data will be on the exact same scale.\n",
    "\n",
    "After performing min-max normalization, you can multiply all cell values by the same number to change the range of values. For example, you can multiply the cell values by 5, so that the minimum cell value is 0 and the maximum cell value is 5. \n",
    "\n",
    "In this example, we perform min-max normalization and then multiply all cell values by 100 because it's easier to view the distribution of cell values in a histogram with a wider range."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "def min_max_normalization(dataframe, col, factor=100):\n",
    "    stats = dataframe.agg(rf_agg_stats(col).alias('stats')).select('stats.min', 'stats.max').first()\n",
    "    \n",
    "    mn = stats.min\n",
    "    mx = stats.max\n",
    "    \n",
    "    dataframe = dataframe.withColumn(col + '_min_max_normalized', rf_local_multiply(rf_local_divide(rf_local_subtract(col, mn), (mx - mn)), factor))\n",
    "    return dataframe\n",
    "\n",
    "rf_scene1 = min_max_normalization(rf_scene1, 'red')\n",
    "rf_scene1 = min_max_normalization(rf_scene1, 'grn')\n",
    "rf_scene1 = min_max_normalization(rf_scene1, 'blu')\n",
    "\n",
    "rf_scene2 = min_max_normalization(rf_scene2, 'red')\n",
    "rf_scene2 = min_max_normalization(rf_scene2, 'grn')\n",
    "rf_scene2 = min_max_normalization(rf_scene2, 'blu')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After performing image normalization, you can see that the range of cell values in the band histograms are the same."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABJEAAAEVCAYAAABQV3TbAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy86wFpkAAAACXBIWXMAAAsTAAALEwEAmpwYAAA/bElEQVR4nO3dfbzVZZno/88lKFKmJlAhGxIH8AEzJpDBmpkoD2pNpjammJOOaVQ+Zo2V2fnlNMMxz5mihzGN0tAmRQcrPR4fUlMrBxUocxQTTUm2MEqkRmMa0PX7Y92bFpu199rPaz983q/Xfu21rvX9ftd9L/ZeF/ta90NkJpIkSZIkSVJ7dmh0AyRJkiRJktT/WUSSJEmSJElSXRaRJEmSJEmSVJdFJEmSJEmSJNVlEUmSJEmSJEl1WUSSJEmSJElSXRaRJElqoIi4NCL+Z6PbMVhExN9HxE+q7v8uIvbu4ee4KyJO7clrdlVErI6I/9GN8/35kyRJHWYRSZKkXlD+uP9DRIxuFX8gIjIi9gLIzA9n5j918Tlml2t9t1X8jSV+V1fb385zXlCu/d6q2PDqPvUnmblLZj7RV89XXp9NpXj1u4h4JCL+tq+ev1VbFkXEP7eK7VX+rYZDx3/+uluskiRJg4NFJEmSes+TwPEtdyLiDcDIHn6O9cCbI2JUVewkYFUPP0+13wCfi4hh3b1QSzFjkLmmFK92AT4K/FtEvLbBbeq3BunPgCRJg5JFJEmSes+3gROr7p8EXFl9QPVokTKyqDkiPh4Rz0bEuog4uc5z/AH4PjC3XGMYcCzwnVbP8+WIWBMRv42IFRHxV1WP3RQRX6i6f01EXN7Oc95Snvfvaj0YEbtFxJURsT4ifhURn4mIHcpjfx8R90TEgoj4DXBBeQ2+FhE3l9E790TE6yLiSxHxXET8IiL+vOr6n4qIX0bExohYGRFHt9XQMupmUkTsWTU66HcR8WJEZNVxHyijhp6LiFsj4vVVj80pbXghIv4ViHZem21k5q3ARuDPyrVeHRE3ltfmuXK7qeq57oqIfyqvwcaI+EH1aLaIeH95TTdExPkdbUdbWv38jS7teT4ifhMRP46IHSLi28AE4P+W1+4T5fh3R8TD5fi7ImK/quu+KSJ+Vvrw7+VnqvXP+Scj4r+Ab3XwdfnniPiP0ob/GxGjIuI75Wd6WZSRcFGxoPwOvRARD0bEAd19rSRJkkUkSZJ6073ArhGxXynuHAf8W51zXgfsBowDTgEujohX1znnSv5UrDoMeBhY2+qYZcA0YA/gKuDfI2Ln8tgHgPdHxNsj4gTgIODsdp4vgf8JfDYidqzx+FdLH/YG3lraVl0M+wvgCeA1wPwSOxb4DDAaeBlYCvy03F8CfLHq/F8Cf1We4x+pjPQZ2057ycy1LaODygih7wGLASLiKODTwHuAMcCPgavLY6OB66ra9kvgLe09V4tSzPgbYCdgZQnvAHwLeD2VwszvgX9tder7qLxerynn/kO53v7AJcD7gT2BUUATPefjQDOV1+C1VF6TzMz3A08BR5TX739HxBQqr9FHy/E3USky7RQRO1F5fRdR+Xm7Gmhd6Htdeez1wDw69rrMpdL3cVSKckvLOXsAjwCfLccdCvw1MAXYncrv3YYuvyqSJGkri0iSJPWultFIc4BfAE/XOX4T8LnM3JSZNwG/A/Zp74TM/A9gj4jYpzzXlTWO+bfM3JCZmzPzC8CIlutm5n8BHwauAL4MnJiZG+s85w1UptJts8B0VbHsvMzcmJmrgS9Q+eO/xdrM/Gppy+9L7HuZuSIzX6JSgHgpM6/MzC3ANcDWkUiZ+e+lKPTHzLwGeAyY2V57W7Xxk8C+VIpnAB8CLszMRzJzM/C/gGllNNI7gZWZuSQzNwFfAv6rzlMcGxHPA/8N3AD8r8x8vrR9Q2Zel5kvltd4PpVCW7VvZeaq8tpcS6X4B3AMcGNm/igzX6ZSyPtjnbb8Qxkp9Hxp04PtHLsJGAu8vvz8/Tgzs41jjwP+X2beVl6Xf6EyVfPNwCxgOPCVcp3vAve3Ov+PwGcz8+XM/H0nXpdfZuYLwM3ALzPz9vJv9u/86WdkE/AqKv/GUf5d19V5nSRJUgdYRJIkqXd9m8rIkr+nRnGnhg3lj+IWLwK7RMSE6ulYbTzPGcDbqBRhthGVKXKPlOk9z1MZxVO96PeNwDDg0cz8Sevz2/AZ4Hxg56rYaCqjZ35VFfsVldEjLdbUuNYzVbd/X+P+LlV9OTEqC5S3FEYOaNWXNkXEO6iMsjqqqoD1euDLVdf7DZUpa+OojPjZ2t5SVKnV/mrXZubumfkKKiNmToyID5Xnf0VEfL1MSfst8CNg99h2fanqItWLVX1v3Zb/pv4Im38pbdk9M3cHDmzn2P8DPA78ICKeiIhPtXPsnlT9G2fmH0vbWl6zp1sVoFq/ZutLwRDo8OvSoZ+RzPwhlVFMFwPPRMTCiNi1nb5IkqQOsogkSVIvysxfUVlg+53Ad+sc3t51nmo1Hau1bwOnATdl5ovVD0Rl/aNPUpky9upSTHiBbdf2mU9lStDYiDieDsjM26gUHU6rCv+aykiQ11fFJrDtCKy2RrfUVUYHfYNKwWxU6ctDdGCdojJS6wrg2MysLmqsAT5UXWzJzJFlhNc6YHzVNaL6fj1lJNbNwBEl9HEqI8D+IjN3pTLtio60v0ZbXkFlSluPKCPHPp6Ze5f2fiwiDml5uNXha6n6N656XZ4u7RxXYi1av2atr9ed16VWX76SmdOBqVSmtZ3bletIkqRtWUSSJKn3nQK8vYwc6RWZ+SSV6T+1Flt+FbCZyvSz4RHx/wFbR2ZExF9TWYPnxPL11YgYV+M6tZwPfKKqHVuoTMGaHxGvKkWfj1F/LaiOeiWVAsT60vaTqYxEalcZiXI98JkaI60uBc6LiKnl2N0i4r3lsf8HTI2I90RlF7GzqKzn0yFlcejDqaxTBZV/i98Dz0fEHvxpHZ+OWAK8KyL+sqw79Dl68P9yEfGuqCxCHsBvgS3lCyqjfvauOvxa4G8i4pCyLtbHqaxl9R9U1iraApwREcMj4kjqTzfszuvSuh8HRcRflHb9N/BSVT8kSVI3WESSJKmXlXVclvfB8/wkM1svqA1wK5XRMKuoTEF6iTK9qBRXrgTOyMynS4HlMio7ZtUdBZKZ97D9ejdnUvnj/QngJ1QW8m5vt7cOy8yVVNZYWkqlsPEG4J4OnPomKiNdvth6WmBmfg+4CFhcplI9BLyjPPZr4L3A56lMHZvcgec7rur6y8rx/1ge+xKVtYN+TWXh9Vs60HZKWx4GTqfyeq4DnqOyEHZPmQzcTmUdrqXA1zLzrvLYhcBnypS/f8jMR6nszvdVKn05gsrC23/IzD9QWaT8FOD5ctyNVIpMbfkSXXxdatiVymi156j8vG+gsmaTJEnqpmh7vURJkiSp+yLiPuDSzPxWo9siSZK6zpFIkiRJ6lER8daIeF2ZznYSlQW9uzO6SJIk9QPDG90ASZIkDTr7UFk3aRfgl8AxmbmusU2SJEnd5XQ2qRdFxAXApMz8u0a3RZIkSVJjRMRq4NTMvL1VfDbwb5nZ1IBmSZ3mdDZJkiRJkiTVZRFJ6oSyvbMkSV1iHpEkSQOZRSSpjohYHRGfjIgHgf+OiL+MiP8o2xz/vAxBbTl2YkTcHREbI+I2YHSj2i1J6jsR8aaI+Fl5///3iLgmIv45ImZHRHPJI/8FfCsiLoiIayPiynL8wxExo9F9kCT1uoMiYmVEPBcR34qInVsfEBEZEZOq7i+KiH/u22ZKbbOIJHXM8cDfAHsD1wP/DOwB/ANwXUSMKcddBaygUjz6J+Ckvm+qJKkvRcROwPeARVRyw9XA0VWHvK7EXw/MK7F3A4uB3YEbgH/tm9ZKkhroBOAw4M+AKcBnGtscqfMsIkkd85XMXAP8HXBTZt6UmX/MzNuA5cA7I2ICcBDwPzPz5cz8EfB/G9hmSVLfmEVlx9uvZOamzPwucH/V438EPltyw+9L7Ccll2wBvg28sW+bLElqgH/NzDWZ+RtgPpUPqqUBxSKS1DFryvfXA+8tU9mej4jngb8ExgJ7As9l5n9Xnfervm2mJKkB9gSezm23vF1TdXt9Zr7U6pz/qrr9IrCz6yVJ0qBXnRt+RSV/SAOKRSSpY1r+MFgDfDszd6/6emVmfh5YB7w6Il5Zdd6EPm+pJKmvrQPGRURUxcZX3U4kSdo2N0wA1tY45kXgFVX3X9erLZI6ySKS1Dn/BhwREYdFxLCI2LksmtqUmb+iMrXtHyNip4j4S+CIxjZXktQHlgJbgDMiYnhEHAnMbHCbJEn9z+kR0RQRewCfBq6pccwDwPvK3xqHA2/tywZK9VhEkjqhrIt0JJU3/fVURiady59+l94H/AXwG+CzwJUNaKYkqQ9l5h+A9wCnAM9TWT/vRuDlBjZLktT/XAX8AHiifNXade1sKh9EP09lIe7v91HbpA6JbafvS5Ikqbsi4j7g0sz8VqPbIkmS1FMciSRJktRNEfHWiHhdmc52EnAgcEuj2yVJktST3AVEkiSp+/YBrgV2AX4JHJOZ6xrbJEmSpJ7ldDZJkiRJkiTV5XQ2SZIkSZIk1TVgp7ONHj0699prr0Y3Q5L6nRUrVvw6M8c0uh2NZp6QpNrMExXmCUmqrb08MWCLSHvttRfLly9vdDMkqd+JiF81ug39gXlCkmozT1SYJySptvbyhNPZJEmSJEmSVJdFJEmSJEmSJNVlEUmSJEmSJEl1Ddg1kSQ13qZNm2hubuall15qdFOGpJ133pmmpiZ23HHHRjdFkmoyTzSWeUJSf2eeaKyu5AmLSJK6rLm5mVe96lXstddeRESjmzOkZCYbNmygubmZiRMnNro5klSTeaJxzBOSBgLzRON0NU84nU1Sl7300kuMGjXKN/wGiAhGjRrlpzaS+jXzROOYJyQNBOaJxulqnrCIJKlbfMNvHF97SQOB71WN42svaSDwvapxuvLaW0SSJEmSJElSXa6JJKnHLLhtVY9e75w5U3r0epKkxjJPSJLaY57o/xyJVGXBbatqfknq3+bPn8/UqVM58MADmTZtGvfdd19D2nHuueey7777cuCBB3L00Ufz/PPPb33swgsvZNKkSeyzzz7ceuutW+Pnn38+48ePZ5dddtnmWosWLWLMmDFMmzaNadOm8c1vfrOvuqF2mCekgck8oX7tzgtrf0nqM+aJjrOIJGlAW7p0KTfeeCM//elPefDBB7n99tsZP358Q9oyZ84cHnroIR588EGmTJnChRdW/gO4cuVKFi9ezMMPP8wtt9zCaaedxpYtWwA44ogjuP/++2te77jjjuOBBx7ggQce4NRTT+2zfkjSYGKekCS1xzzRORaRJA1o69atY/To0YwYMQKA0aNHs+eee7Js2TLe/OY388Y3vpGZM2eyceNGtmzZwrnnnstBBx3EgQceyNe//nUA7rrrLmbPns0xxxzDvvvuywknnEBmArBixQre+ta3Mn36dA477DDWrVvXZlsOPfRQhg+vzBKeNWsWzc3NAFx//fXMnTuXESNGMHHiRCZNmrT1jX7WrFmMHTu2114fSRrqzBOSpPaYJzrHIpKkAe3QQw9lzZo1TJkyhdNOO427776bP/zhDxx33HF8+ctf5uc//zm33347I0eO5LLLLmO33XZj2bJlLFu2jG984xs8+eSTAPzsZz/jS1/6EitXruSJJ57gnnvuYdOmTZx55pksWbKEFStW8IEPfIDzzz+/Q+26/PLLecc73gHA008/vc2nGU1NTTz99NN1r3Hddddx4IEHcswxx7BmzZouvDqSJPOEJKk95onOcWFtSQPaLrvswooVK/jxj3/MnXfeyXHHHcf555/P2LFjOeiggwDYddddAfjBD37Agw8+yJIlSwB44YUXeOyxx9hpp52YOXMmTU1NAEybNo3Vq1ez++6789BDDzFnzhwAtmzZ0qEq//z58xk+fDgnnHACwNZPIarV207ziCOO4Pjjj2fEiBFceumlnHTSSfzwhz/s4KsiSWphntBAtPSJDdy7efs191wkWOp55onOGZpFpDYXqvvbPm2GpJ4xbNgwZs+ezezZs3nDG97AxRdfXPNNNTP56le/ymGHHbZN/K677to6fLXleps3byYzmTp1KkuXLu1wW6644gpuvPFG7rjjjq1taGpq2qby39zczJ577tnudUaNGrX19gc/+EE++clPdrgNkqRtmSckSe0xT3Tc0CwiSeoVjfh07NFHH2WHHXZg8uTJADzwwAPst99+3HLLLSxbtoyDDjqIjRs3MnLkSA477DAuueQS3v72t7PjjjuyatUqxo0b1+a199lnH9avX8/SpUs5+OCD2bRpE6tWrWLq1Kk1j7/lllu46KKLuPvuu3nFK16xNf7ud7+b973vfXzsYx9j7dq1PPbYY8ycObPdfq1bt27rpxQ33HAD++23X2dfGknqd8wT5glJao95ov/nCYtIkga03/3ud5x55pk8//zzDB8+nEmTJrFw4UJOPvlkzjzzTH7/+98zcuRIbr/9dk499VRWr17Nm970JjKTMWPG8P3vf7/Na++0004sWbKEs846ixdeeIHNmzfz0Y9+tM03/TPOOIOXX35563DVWbNmcemllzJ16lSOPfZY9t9/f4YPH87FF1/MsGHDAPjEJz7BVVddxYsvvkhTUxOnnnoqF1xwAV/5yle44YYbGD58OHvssQeLFi3q6Zeux0TE5cC7gGcz84ASuwbYpxyyO/B8Zk6LiL2AR4BHy2P3ZuaHyznTgUXASOAm4OzMzIgYAVwJTAc2AMdl5ure75mkwcA8IUlqj3mic6LW3LqBYMaMGbl8+fKundzGdLYFm2tPZ3PusVTbI4884iefDVbr3yAiVmTmjL5qQ0T8NfA74MqWIlKrx78AvJCZnytFpBvbOO5+4GzgXipFpK9k5s0RcRpwYGZ+OCLmAkdn5nH12tWtPNGGBbdtvz4FmCektpgnGq8/5In+qjfyRJfU+Ntk6RMbah568N6j4G3n9XaLpD5jnmi8zuYJd2eTJHVLZv4I+E2tx6IykftY4Or2rhERY4FdM3NpVj7duBI4qjx8JHBFub0EOCTqrSQoSZIkqcc5nU2SOun000/nnnvu2SZ29tlnc/LJJzeoRf3aXwHPZOZjVbGJEfEz4LfAZzLzx8A4oLnqmOYSo3xfA5CZmyPiBWAU8OvWTxYR84B5ABMmTOjhrrSjrQ0b/LRYGpLME5Kk9gzkPGERSZI66eKLL250EwaS49l2FNI6YEJmbihrIH0/IqYCtUYWtcy3bu+xbYOZC4GFUJmm0OVWS1I3mCckSe0ZyHnC6WySpF4REcOB9wDXtMQy8+XM3FBurwB+CUyhMvKoqer0JmBtud0MjK+65m60MX1OktT/RMT4iLgzIh6JiIcj4uwS3yMibouIx8r3V1edc15EPB4Rj0bEYVXx6RHxn+Wxr7RMb46IERFxTYnfV9bgkyT1MItIkqTe8j+AX2Tm1mlqETEmIoaV23sDk4EnMnMdsDEiZpU/CE4Eri+n3QCcVG4fA/wwB+quEJI0NG0GPp6Z+wGzgNMjYn/gU8AdmTkZuKPcpzw2F5gKHA58rSV3AJdQmbY8uXwdXuKnAM9l5iRgAXBRX3RMkoYai0iSpG6JiKuBpcA+EdEcEaeUh+ay/YLafw08GBE/p7JI9oczs2VU0UeAbwKPUxmhdHOJXwaMiojHgY9R/siQJA0MmbkuM39abm8EHqGy3l31xglXsO2GCovL6NUnqeSFmW7CIEmN55pIknpOW4sLd5WLEg8ImXl8G/G/rxG7DriujeOXAwfUiL8EvLd7rZTUL5gnhrwyzezPgfuA15aRqGTmuoh4TTlsHHBv1Wktmy1sopubMDRsAwZJHWOe6PcciSRpwJs/fz5Tp07lwAMPZNq0adx3330Nace5557Lvvvuy4EHHsjRRx/N888/v/WxCy+8kEmTJrHPPvtw6623bo2ff/75jB8/nl122WWba51zzjlMmzaNadOmMWXKFHbfffc+6oUkDT7mif4hInah8kHCRzPzt+0dWiOW7cTbO2fbQObCzJyRmTPGjBlTr8mShgjzRMdZRJI0oC1dupQbb7yRn/70pzz44IPcfvvtjB8/viFtmTNnDg899BAPPvggU6ZM4cILK5+krFy5ksWLF/Pwww9zyy23cNppp7FlyxYAjjjiCO6///7trrVgwQIeeOABHnjgAc4880ze85739GlfJGmwME/0DxGxI5UC0ncy87sl/EyZokb5/myJb91QoWjZbMFNGCT1OPNE51hEkjSgrVu3jtGjRzNixAgARo8ezZ577smyZct485vfzBvf+EZmzpzJxo0b2bJlC+eeey4HHXQQBx54IF//+tcBuOuuu5g9ezbHHHMM++67LyeccAIt6zavWLGCt771rUyfPp3DDjuMdevWtdmWQw89lOHDK7OEZ82aRXNzZcT99ddfz9y5cxkxYgQTJ05k0qRJW9/oZ82axdixY9vt49VXX83xx9ecMSZJqsM80XhlbaLLgEcy84tVD1VvnHAS226oMLfsuDaRygLa97sJg6TeYJ7oHItIkga0Qw89lDVr1jBlyhROO+007r77bv7whz9w3HHH8eUvf5mf//zn3H777YwcOZLLLruM3XbbjWXLlrFs2TK+8Y1v8OSTTwLws5/9jC996UusXLmSJ554gnvuuYdNmzZx5plnsmTJElasWMEHPvABzj///A616/LLL+cd73gHAE8//fQ2n2Y0NTXx9NNPd+g6v/rVr3jyySd5+9vf3slXRpIE5ol+4i3A+4G3R8QD5eudwOeBORHxGDCn3CczHwauBVYCtwCnZ+aWci03YZDUo8wTnePC2pIGtF122YUVK1bw4x//mDvvvJPjjjuO888/n7Fjx3LQQQcBsOuuuwLwgx/8gAcffJAlS5YA8MILL/DYY4+x0047MXPmTJqaKiPkp02bxurVq9l999156KGHmDNnDgBbtmypW+WHypzq4cOHc8IJJwBQ64PQjm4Ys3jxYo455hiGDRtW/2BJ0nbME42XmT+h9ppFAIe0cc58YH6NuJswSOpR5onOsYgkacAbNmwYs2fPZvbs2bzhDW/g4osvrvmmmpl89atf5bDDDtsmftddd20dvtpyvc2bN5OZTJ06laVLl3a4LVdccQU33ngjd9xxx9Y2NDU1sWbNmq3HNDc3s+eee3boeosXL+biiy/u8PNLkrZnnpAktcc80XEWkST1nAZsofnoo4+yww47MHnyZAAeeOAB9ttvP2655RaWLVvGQQcdxMaNGxk5ciSHHXYYl1xyCW9/+9vZcccdWbVqFePGjWvz2vvssw/r169n6dKlHHzwwWzatIlVq1YxderUmsffcsstXHTRRdx999284hWv2Bp/97vfzfve9z4+9rGPsXbtWh577DFmzpzZob4999xzHHzwwZ18VSSpnzJPmCckqT3miX6fJywiSRrQfve733HmmWfy/PPPM3z4cCZNmsTChQs5+eSTOfPMM/n973/PyJEjuf322zn11FNZvXo1b3rTm8hMxowZw/e///02r73TTjuxZMkSzjrrLF544QU2b97MRz/60Tbf9M844wxefvnlrcNVZ82axaWXXsrUqVM59thj2X///Rk+fDgXX3zx1uGkn/jEJ7jqqqt48cUXaWpq4tRTT+WCCy4AKgvgzZ07t8NDVSVJ2zNPSJLaY57onBiomxbMmDEjly9f3rWT77ywZnjB5r+tGT9nzpSuPY80yD3yyCPst99+jW7GkFbr3yAiVmTmjAY1qd/oVp5ow4LbVtWMnzP8utonNODTNKk/MU80nnmibb2RJ7qkxt8mS5/YUPPQg/ceZW7RoGKeaLzO5gl3Z5MkSZIkSVJdTmeTpE46/fTTueeee7aJnX322Zx88skNapEkqT8xT0iS2jOQ80TdIlJEjAeuBF4H/BFYmJlfjog9gGuAvYDVwLGZ+Vw55zzgFGALcFZm3lri04FFwEjgJuDszMyIGFGeYzqwATguM1f3WC8l9ZrMHHJrMfSXXXAG6nRkSUOLeaJxzBOSBgLzRON0JU90ZDrbZuDjmbkfMAs4PSL2Bz4F3JGZk4E7yn3KY3OBqcDhwNciYli51iXAPGBy+Tq8xE8BnsvMScAC4KJO90RSn9t5553ZsGGD/0ltgMxkw4YN7Lzzzo1uiiS1yTzROOYJSQOBeaJxupon6o5Eysx1wLpye2NEPAKMA44EZpfDrgDuAj5Z4osz82XgyYh4HJgZEauBXTNzKUBEXAkcBdxczrmgXGsJ8K8REelPktSvNTU10dzczPr16xvdlCFp5513pqmpqdHNkKQ2mScayzwhqb8zTzRWV/JEp9ZEioi9gD8H7gNeWwpMZOa6iHhNOWwccG/Vac0ltqncbh1vOWdNudbmiHgBGAX8utXzz6MykokJEyZ0pumSesGOO+7IxIkTG90MSVI/ZZ6QJLXHPDHwdHh3tojYBbgO+Ghm/ra9Q2vEsp14e+dsG8hcmJkzMnPGmDFj6jVZkiRJkiRJPaRDRaSI2JFKAek7mfndEn4mIsaWx8cCz5Z4MzC+6vQmYG2JN9WIb3NORAwHdgN+09nOSJIkSZIkqXfULSJFZZn0y4BHMvOLVQ/dAJxUbp8EXF8VnxsRIyJiIpUFtO8vU982RsSscs0TW53Tcq1jgB+6HpIkSZIkSVL/0ZE1kd4CvB/4z4h4oMQ+DXweuDYiTgGeAt4LkJkPR8S1wEoqO7udnplbynkfARYBI6ksqH1ziV8GfLsswv0bKru7SZIkSZIkqZ/oyO5sP6H2mkUAh7Rxznxgfo34cuCAGvGXKEUoSdLAEhGXA+8Cns3MA0rsAuCDQMtWG5/OzJvKY+cBpwBbgLMy89YSn86fPmi4CTg7MzMiRgBXAtOBDcBxmbm6TzonSZIkaasOL6wtSVIbFgGH14gvyMxp5aulgLQ/ldGmU8s5X4uIYeX4S6jswDm5fLVc8xTgucycBCwALuqtjkiSJElqm0UkSVK3ZOaP6PhmCEcCizPz5cx8EngcmFk2aNg1M5eWNfGuBI6qOueKcnsJcEhZW0+SJElSH7KIJEnqLWdExIMRcXlEvLrExgFrqo5pLrFx5Xbr+DbnZOZm4AVgVK0njIh5EbE8IpavX7++1iGSJEmSusgikiSpN1wC/BkwDVgHfKHEa40gynbi7Z2zfTBzYWbOyMwZY8aM6VSDJUmSJLXPIpIkqcdl5jOZuSUz/wh8A5hZHmoGxlcd2gSsLfGmGvFtzomI4cBudHz6nCRJkqQeYhFJktTjyhpHLY4GHiq3bwDmRsSIiJhIZQHt+zNzHbAxImaV9Y5OBK6vOuekcvsY4Idl3SRJkiRJfWh4oxsgSRrYIuJqYDYwOiKagc8CsyNiGpVpZ6uBDwFk5sMRcS2wEtgMnJ6ZW8qlPkJlp7eRwM3lC+Ay4NsR8TiVEUhze71TkiRJkrZjEUmS1C2ZeXyN8GXtHD8fmF8jvhw4oEb8JeC93WmjJEmSpO5zOpskSZIkSZLqsogkSZIkSZKkuiwiSZIkSZIkqS6LSJIkSZIkSarLIpIkSZIkSZLqsogkSZIkSZKkuiwiSZIkSZIkqS6LSJIkSZIkSarLIpIkSZIkSZLqsogkSZIkSZKkuiwiSZIkSZIkqS6LSJIkSZIkSarLIpIkSZIkSZLqsogkSZIkSZKkuiwiSZIkSZIkqS6LSJIkSZIkSarLIpIkSZKkXhMRl0fEsxHxUFXsgoh4OiIeKF/vrHrsvIh4PCIejYjDquLTI+I/y2NfiYgo8RERcU2J3xcRe/VpByVpCLGIJEmSJKk3LQIOrxFfkJnTytdNABGxPzAXmFrO+VpEDCvHXwLMAyaXr5ZrngI8l5mTgAXARb3VEUka6iwiSZIkSeo1mfkj4DcdPPxIYHFmvpyZTwKPAzMjYiywa2YuzcwErgSOqjrninJ7CXBIyyglSVLPsogkSeqWNqYp/J+I+EVEPBgR34uI3Ut8r4j4fdX0hUurznGagiQNLWeUPHF5RLy6xMYBa6qOaS6xceV26/g252TmZuAFYFStJ4yIeRGxPCKWr1+/vud6IklDhEUkSVJ3LWL7aQq3AQdk5oHAKuC8qsd+WTV94cNVcacpSNLQcQnwZ8A0YB3whRKvNYIo24m3d872wcyFmTkjM2eMGTOmUw2WJFlEkiR1U61pCpn5g/JpMMC9QFN713CagiQNLZn5TGZuycw/At8AZpaHmoHxVYc2AWtLvKlGfJtzImI4sBsdnz4nSeoEi0iSpN72AeDmqvsTI+JnEXF3RPxViTlNQZKGkPLhQYujgZYp0TcAc8tU5olURqben5nrgI0RMat8kHAicH3VOSeV28cAPywfSEiSetjwRjdAkjR4RcT5wGbgOyW0DpiQmRsiYjrw/YiYSg9OUwAWAsyYMcM/ICSpH4iIq4HZwOiIaAY+C8yOiGlU3s9XAx8CyMyHI+JaYCWV/HF6Zm4pl/oIlSnUI6l8ONHyAcVlwLcj4nEqI5Dm9nqnJGmIsogkSeoVEXES8C7gkJZPhDPzZeDlcntFRPwSmELHpik0O01BkgaezDy+Rviydo6fD8yvEV8OHFAj/hLw3u60UZLUMU5nkyT1uIg4HPgk8O7MfLEqPiYihpXbe1OZpvCE0xQkSZKk/s+RSJKkbmljmsJ5wAjgtrIG9r1lJ7a/Bj4XEZuBLcCHM7NlVJHTFCRJkqR+rG4RKSIupzId4dnMPKDELgA+CLSsWvrpzLypPHYele2YtwBnZeatJT6dP/1xcBNwdmZmRIygsgvPdGADcFxmru6h/kmSellnpilk5nXAdW081u+nKcx6amHtB/auuc63JEmSNKh0ZDrbIuDwGvEFmTmtfLUUkPan8gnx1HLO11qmLQCXAPOoTF2YXHXNU4DnMnMSsAC4qIt9kSRJkiRJUi+pW0TKzB/R8QVMjwQWZ+bLmfkk8Dgws2zhuWtmLi3rWFwJHFV1zhXl9hLgkLIehiRJkiRJkvqJ7iysfUZEPBgRl0fEq0tsHLCm6pjmEhtXbreOb3NOZm4GXgBqzguIiHkRsTwilq9fv77WIZIkSZIkSeoFXS0iXQL8GTANWAd8ocRrjSDKduLtnbN9MHNhZs7IzBljxozpVIMlSZIkSZLUdV0qImXmM5m5JTP/CHwDmFkeagbGVx3aBKwt8aYa8W3OiYjhwG50fPqcJEmSJEmS+kCXikhljaMWRwMPlds3AHMjYkRETKSygPb9mbkO2BgRs8p6RycC11edc1K5fQzww7JukiRJkiRJkvqJ4fUOiIirgdnA6IhoBj4LzI6IaVSmna0GPgSQmQ9HxLXASmAzcHpmbimX+giVnd5GAjeXL6hsA/3tiHicygikuT3QL0mSJEmSJPWgukWkzDy+Rviydo6fD8yvEV8OHFAj/hLw3nrtkCRJkqQB684Lu32JpU9s4N7Nq7aLnzNnSrevLUkd0Z3d2SRJkiRJkjREWESSJEmSJElSXXWns0mSJEmSOqgHpq1JUn/lSCRJkiRJkiTVZRFJkiRJkiRJdVlEkiRJkiRJUl0WkSRJkiRJklSXRSRJkiRJkiTVZRFJkiRJkiRJdVlEkiRJkiRJUl0WkSRJkiRJklTX8EY3QJI0sEXE5cC7gGcz84AS2wO4BtgLWA0cm5nPlcfOA04BtgBnZeatJT4dWASMBG4Czs7MjIgRwJXAdGADcFxmru6j7nXPnRfWjr/tvL5tR7HgtlU14+fMmdLHLZEkSdJA5EgkSVJ3LQIObxX7FHBHZk4G7ij3iYj9gbnA1HLO1yJiWDnnEmAeMLl8tVzzFOC5zJwELAAu6rWeSJIkSWqTRSRJUrdk5o+A37QKHwlcUW5fARxVFV+cmS9n5pPA48DMiBgL7JqZSzMzqYw8OqrGtZYAh0RE9EZfJEmSJLXNIpIkqTe8NjPXAZTvrynxccCaquOaS2xcud06vs05mbkZeAEY1WstlyRJklSTRSRJUl+qNYIo24m3d872F4+YFxHLI2L5+vXru9hESZIkSbW4sLYkqTc8ExFjM3Ndmar2bIk3A+OrjmsC1pZ4U4149TnNETEc2I3tp88BkJkLgYUAM2bMqFlokiRpIJv11MLtg3eOatimDZKGFkciSZJ6ww3ASeX2ScD1VfG5ETEiIiZSWUD7/jLlbWNEzCrrHZ3Y6pyWax0D/LCsmyRJkiSpDzkSSZLULRFxNTAbGB0RzcBngc8D10bEKcBTwHsBMvPhiLgWWAlsBk7PzC3lUh+hstPbSODm8gVwGfDtiHicygikuX3QLUmSJEmtWESSJHVLZh7fxkOHtHH8fGB+jfhy4IAa8ZcoRShJkiRJjWMRSZIkSZIaYOkTGxrdBEnqFNdEkiRJkiRJUl0WkSRJkiRJklSX09kkSRps7rywjQf+tk+bIUmSpMHFkUiSJEmSJEmqyyKSJEmSJEmS6rKIJEmSJEmSpLosIkmSJEmSJKkuF9aWJKmXLH1iQ834wW/r44ZIUgNFxOXAu4BnM/OAEtsDuAbYC1gNHJuZz5XHzgNOAbYAZ2XmrSU+HVgEjARuAs7OzIyIEcCVwHRgA3BcZq7uo+5J0pDiSCRJkiRJvWkRcHir2KeAOzJzMnBHuU9E7A/MBaaWc74WEcPKOZcA84DJ5avlmqcAz2XmJGABcFGv9USShjiLSJIkSZJ6TWb+CPhNq/CRwBXl9hXAUVXxxZn5cmY+CTwOzIyIscCumbk0M5PKyKOjalxrCXBIRERv9EWShjqLSJIkSZL62mszcx1A+f6aEh8HrKk6rrnExpXbrePbnJOZm4EXgFG1njQi5kXE8ohYvn79+h7qiiQNHRaRJEmSJPUXtUYQZTvx9s7ZPpi5MDNnZOaMMWPGdLGJkjR0WUSSJEmS1NeeKVPUKN+fLfFmYHzVcU3A2hJvqhHf5pyIGA7sxvbT5yRJPcDd2SRJ6i/uvLB2/G3n9W07JKn33QCcBHy+fL++Kn5VRHwR2JPKAtr3Z+aWiNgYEbOA+4ATga+2utZS4Bjgh2XdpCFj6RMbuHfzqu3i58yZ0oDWSBrM6o5EiojLI+LZiHioKrZHRNwWEY+V76+ueuy8iHg8Ih6NiMOq4tMj4j/LY19pWewuIkZExDUlfl9E7NXDfZQkSZLUIBFxNZUCzz4R0RwRp1ApHs2JiMeAOeU+mfkwcC2wErgFOD0zt5RLfQT4JpXFtn8J3FzilwGjIuJx4GOUnd4kST2vIyORFgH/SmUHhBYtW3J+PiI+Ve5/stWWnHsCt0fElPLG37Il573ATVS25LyZqi05I2IulS05j+uJzkmSJElqrMw8vo2HDmnj+PnA/Brx5cABNeIvAe/tThslSR1TdySSW3JKkiRJkiSpqwtruyWnJEmSJEnSENLTC2v3+pacwEKAGTNmDKnF8iRpoImIfYBrqkJ7A/8fsDvwQaDl04BPZ+ZN5ZzzqExz3gKclZm3lvh0KtOrR1KZEn32UFs0VZI0cC19YkOjmyBJPaKrRaRnImJsZq7rwS05m92SU5IGj8x8FJgGEBHDgKeB7wEnAwsy81+qj+/iunpDW1u7uUmSBCy4bfsd28Bd2yR1XVens7Vsownbb8k5t+y4NpE/bcm5DtgYEbPKekcntjqn5VpDcktOSRoCDgF+mZm/aueYrqyrJ0mSJKmP1C0iuSWnJKkHzAWurrp/RkQ8GBGXR8SrS6wr6+pJkiRJ6iN1p7O5JackqTsiYifg3cB5JXQJ8E9U1r/7J+ALwAfo2rp6rZ9rHpVpb0yYMKFb7ZYkSZK0ra5OZ5MkqaPeAfw0M58ByMxnMnNLZv4R+AYwsxzXlXX1tpGZCzNzRmbOGDNmTA93Q5IkSRraLCJJknrb8VRNZStrHLU4Gnio3O7KunqSJEmS+khXd2dTT2prd523nVc7LkkDRES8gsraeR+qCv/viJhGZUra6pbHMvPhiGhZV28z26+rtwgYSWVNPXdmkyRJkvqYRSRJUq/JzBeBUa1i72/n+E6tqydJUn+y4LZVzHpqQ6ObIUm9xiKSJEl9ra0RqD1k6RNt/AHjWuOSJEnqBtdEkiRJkiRJUl2ORJIkqb/r5ZFLkiRJUkdYROoIF76WJEmSJElDnNPZJEmSJEmSVJcjkarMemph7Qf2HlU7LkkS7SxkLUmSJA0iFpEkSRqgLF5JkiSpLzmdTZIkSZIkSXVZRJIkSZIkSVJdFpEkSZIkSZJUl2siSZLUT7S1xtHBbvAgSZKkfsCRSJIkSZIkSarLIpIkSZIkSZLqsogkSZIkSZKkulwTSZIkSZI6684LtwvNeqr22naSNFg4EkmSJEmSJEl1ORJpIKrxqQcAbzuvb9shSZIkSZKGDItIHVBry2W3W5YkDVULbltVM37OnCl93BJJkiT1JYtIkiT1sVofTkiSJEn9nUWk3tBT083auo4kDRARsRrYCGwBNmfmjIjYA7gG2AtYDRybmc+V488DTinHn5WZt5b4dGARMBK4CTg7M7Mv+yJJ0mBRa0Spo0kldYQLa0uSetvbMnNaZs4o9z8F3JGZk4E7yn0iYn9gLjAVOBz4WkQMK+dcAswDJpevw/uw/ZIkSZJwJFL3OFJIkrriSGB2uX0FcBfwyRJfnJkvA09GxOPAzDKaadfMXAoQEVcCRwE392mrh6C21j6SJEnS0ORIJElSb0rgBxGxIiLmldhrM3MdQPn+mhIfB6ypOre5xMaV263jkiRJkvqQI5EkSb3pLZm5NiJeA9wWEb9o59ioEct24ttfoFKomgcwYcKEzrZVkiRJUjssIkmSek1mri3fn42I7wEzgWciYmxmrouIscCz5fBmYHzV6U3A2hJvqhGv9XwLgYUAM2bMGDQLb/f2bm5OW5MkSVJHWETqS729hlJP7QonST0gIl4J7JCZG8vtQ4HPATcAJwGfL9+vL6fcAFwVEV8E9qSygPb9mbklIjZGxCzgPuBE4Kt925vBYdZTC2vG750wr2ZckiRJqmYRSZLUW14LfC8ioJJvrsrMWyJiGXBtRJwCPAW8FyAzH46Ia4GVwGbg9MzcUq71EWARMJLKgtouqt2DLC5JkiSpIywiSZJ6RWY+AbyxRnwDcEgb58wH5teILwcO6Ok2SpIkSeo4i0iqq621Ms6ZM6WPWyJJkqTBJCJWAxuBLcDmzJwREXsA1wB7AauBYzPzuXL8ecAp5fizMvPWEp/On0as3gScnZmDZm08SeovLCL1U20tonrw3qM6fY17N29bBLL4I0mSpH7kbZn566r7nwLuyMzPR8Snyv1PRsT+wFxgKpW1826PiCll6vMlVHbnvJdKEelwnPosST1uh0Y3QJIkSZKqHAlcUW5fARxVFV+cmS9n5pPA48DMstPnrpm5tIw+urLqHElSD+pWESkiVkfEf0bEAxGxvMT2iIjbIuKx8v3VVcefFxGPR8SjEXFYVXx6uc7jEfGVKKuwSpIkSRrUEvhBRKyIiJbV/F+bmesAyvfXlPg4YE3Vuc0lNq7cbh3fTkTMi4jlEbF8/fr1PdgNSRoaemIk0tsyc1pmzij3W4afTgbuKPdpNfz0cOBrETGsnNMy/HRy+Tq8B9olSZIkqX97S2a+CXgHcHpE/HU7x9b6oDnbiW8fzFyYmTMyc8aYMWM631pJGuJ6Yzqbw08lSZIk1ZWZa8v3Z4HvATOBZ8rfCJTvz5bDm4HxVac3AWtLvKlGXJLUw7q7sHbL8NMEvp6ZC2k1/DQiqoef3lt1bssw0010cPipembBbUmSJKnRIuKVwA6ZubHcPhT4HHADcBLw+fL9+nLKDcBVEfFFKgtrTwbuz8wtEbExImYB9wEnAl/t295I0tDQ3SLSWzJzbSkU3RYRv2jn2G4PPy3zpOcBTJgwobNtlSRJktR/vBb4XlkOdThwVWbeEhHLgGsj4hTgKeC9AJn5cERcC6wENgOnl53ZAD4CLAJGUtmVzZ3ZJKkXdKuIVD38NCK2GX5aRiH16PDTMtJpIcCMGTNqFppU34LbVtWMnzNnSh+3RJLUn816auF2sXsnzKtxpCR1XmY+AbyxRnwDcEgb58wH5teILwcO6Ok2SpK21eUiksNPB5+2ikuSJEmSJEndGYnk8FNJkiRJ6gccPSqpL3S5iOTw085xQWxJkiRJkjSQ7dDoBkiSJEmSJKn/6+7ubENWfxtZ1FZ7JEmSJEmSeoJFJHWZu7xJkiRJkjR0OJ1NkiRJkiRJdVlEkiRJkiRJUl1OZ5MkST2irWnObXH6syT1rllPLawZv3fCvD5uiaTBwiKSJEmSJA1xrncqqSMsIjWYu6pJGqwiYjxwJfA64I/Awsz8ckRcAHwQWF8O/XRm3lTOOQ84BdgCnJWZt5b4dGARMBK4CTg7M7PveiNJkiTJItIQ4DBWSQ2yGfh4Zv40Il4FrIiI28pjCzLzX6oPjoj9gbnAVGBP4PaImJKZW4BLgHnAvVSKSIcDN/dRP1TFnCJJkjR0WUTqYY4skqSKzFwHrCu3N0bEI8C4dk45ElicmS8DT0bE48DMiFgN7JqZSwEi4krgKCwiSZIkSX3KItIQ5qfJkvpKROwF/DlwH/AW4IyIOBFYTmW00nNUCkz3Vp3WXGKbyu3WcUmSJEl9yCKSJKlXRcQuwHXARzPztxFxCfBPQJbvXwA+AESN07OdeK3nmkdl2hsTJkzoeqPvvLDr50qSJEmDlEUk9ZnObP3sLhDS4BARO1IpIH0nM78LkJnPVD3+DeDGcrcZGF91ehOwtsSbasS3k5kLgYUAM2bMcOHtfs6dgCRJkgaWHRrdAEnS4BQRAVwGPJKZX6yKj6067GjgoXL7BmBuRIyIiInAZOD+srbSxoiYVa55InB9n3RCkiRJ0laORJIk9Za3AO8H/jMiHiixTwPHR8Q0KlPSVgMfAsjMhyPiWmAllZ3dTi87swF8BFgEjKSyoLaLakuS1AccNSqpmkUkSVKvyMyfUHs9o5vaOWc+ML9GfDlwQM+1TpIkSVJnWUTSdty1TZIkSRq8/P++pK5yTSRJkiRJkiTVZRFJkiRJkiRJdVlEkiRJkiRJUl2uiaQe19YODpIkSZIkaeByJJIkSZIkSZLqsogkSZIkSZKkupzOpg7ry61A25oSd86cKT3+XJKk7nO7aEmSpMHPIpIkSZIkqVP80FcampzOJkmSJEmSpLociaQBxU88JEmSJElqDEciSZIkSZIkqS5HIkmSJEmSeoQzB6TBzSKSJEmSJMmdNiXVZRFJg4KfeEiSJEmS1LssIkmSpF7jp9qSNPD5Xi6phUUkDWqOUJKkgaet9+62+J4uSf2f7+3S4ODubJIkSZIkSarLkUjqNoe3SpIkSUOPfwdIQ0+/KSJFxOHAl4FhwDcz8/MNbpK6qVZSMaFI6irzxODSk394OHVZEpgnBhunv0n9U78oIkXEMOBiYA7QDCyLiBsyc2VjW6bByj84pIGlr/PE0ic29MZl1QF9UVxqizlAGrj6Ok8suG0Vs54yV7TFEUrS4NUvikjATODxzHwCICIWA0cCFpEGmYGYUCw4Sf2CeWKI64v8Uev93vd6acAwTwwAbb2Xt6Uz7/Gd/eCgs8wHUkV/KSKNA9ZU3W8G/qJBbVED9FRC6e4fGZ1JPj2VqNpKSJ25vklNQ4B5QjV1Nn901oLbeiff1NPZ3GAekMwTg1Fvv8e3pdZ7eaOKVJ193zdPqLf1lyJS1IjldgdFzANafqN/FxGPdvH5RgO/7uK5A9Eg7O8X2nqgjb62eXzDfax7p48Gft3Nawwkg/BnuU3d6evre7Ih/URf5Ymh9DMG9rcDOps/eibfdPZ9vY3j/fcd3MwT2+rLvyf8WRvcRsMX+ry/PfS+35Xjh9K/71DqK/RSnugvRaRmYHzV/SZgbeuDMnMh0O1ydEQsz8wZ3b3OQDGU+juU+gr2dzAbSn3toD7JE0Ptdbe/g5v9HdyGWn87oM/+nhhqr739HdyGUn+HUl+h9/q7Q09fsIuWAZMjYmJE7ATMBW5ocJskSf2HeUKS1B7zhCT1gX4xEikzN0fEGcCtVLbkvDwzH25wsyRJ/YR5QpLUHvOEJPWNflFEAsjMm4Cb+ujpGrNCW+MMpf4Opb6C/R3MhlJfO6SP8sRQe93t7+Bmfwe3odbfuvrw74mh9trb38FtKPV3KPUVeqm/kbndenOSJEmSJEnSNvrLmkiSJEmSJEnqx4ZUESkiDo+IRyPi8Yj4VKPb09MiYnxE3BkRj0TEwxFxdonvERG3RcRj5furG93WnhIRwyLiZxFxY7k/mPu6e0QsiYhflH/jgwd5f88pP8cPRcTVEbHzYOpvRFweEc9GxENVsTb7FxHnlfeuRyPisMa0evAyPwz836lazBGDur+DOkeAeaK/Ml8Mjt+vauaKQd3fQZ0rGpUnhkwRKSKGARcD7wD2B46PiP0b26oetxn4eGbuB8wCTi99/BRwR2ZOBu4o9weLs4FHqu4P5r5+GbglM/cF3kil34OyvxExDjgLmJGZB1BZIHMug6u/i4DDW8Vq9q/8Hs8FppZzvlbe09QDzA+D5neqFnPEIOzvEMkRYJ7od8wXg+r3q5q5YhD2d4jkikU0IE8MmSISMBN4PDOfyMw/AIuBIxvcph6Vmesy86fl9kYqbwrjqPTzinLYFcBRDWlgD4uIJuBvgG9WhQdrX3cF/hq4DCAz/5CZzzNI+1sMB0ZGxHDgFcBaBlF/M/NHwG9ahdvq35HA4sx8OTOfBB6n8p6mnmF+qBjQv1OtmSPMEQzw/pon+iXzRcWA//1qYa4wVzCA+9uoPDGUikjjgDVV95tLbFCKiL2APwfuA16bmeugkhiA1zSwaT3pS8AngD9WxQZrX/cG1gPfKsNtvxkRr2SQ9jcznwb+BXgKWAe8kJk/YJD2t0pb/RtS718NMKRe3yGSH8AcYY4YJP1txTzRWEPqdR4i+eJLmCsGZX+HcK7o9TwxlIpIUSM2KLemi4hdgOuAj2bmbxvdnt4QEe8Cns3MFY1uSx8ZDrwJuCQz/xz4bwb20Mt2lbm7RwITgT2BV0bE3zW2VQ01ZN6/GmTIvL5DIT+AOQJzxFA0ZN7HGmzIvM5DIV+YK8wVQ0yPvX8NpSJSMzC+6n4TleFsg0pE7EjlDf87mfndEn4mIsaWx8cCzzaqfT3oLcC7I2I1laHEb4+If2Nw9hUqP7/NmXlfub+EShIYrP39H8CTmbk+MzcB3wXezODtb4u2+jck3r8aaEi8vkMoP4A5whwxuPpbzTzRWEPidR5C+cJcYa4YTP1t0et5YigVkZYBkyNiYkTsRGVRqRsa3KYeFRFBZY7rI5n5xaqHbgBOKrdPAq7v67b1tMw8LzObMnMvKv+WP8zMv2MQ9hUgM/8LWBMR+5TQIcBKBml/qQw7nRURryg/14dQmZM/WPvboq3+3QDMjYgRETERmAzc34D2DVbmh4pB8ztljjBHMLj6W8080Vjmi4pB8ftlrjBXMLj626L380RmDpkv4J3AKuCXwPmNbk8v9O8vqQxJexB4oHy9ExhFZWX2x8r3PRrd1h7u92zgxnJ70PYVmAYsL/++3wdePcj7+4/AL4CHgG8DIwZTf4GrqczP3kTlk4FT2usfcH5573oUeEej2z/YvswPA/93qp2+myMGZ38HdY4ofTRP9MMv88Xg+P2q0W9zxeDs76DOFY3KE1EuJkmSJEmSJLVpKE1nkyRJkiRJUhdZRJIkSZIkSVJdFpEkSZIkSZJUl0UkSZIkSZIk1WURSZIkSZIkSXVZRJIkSZIkSVJdFpEkSZIkSZJUl0UkSZIkSZIk1fX/A3DkdcfvQmlcAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 1440x288 with 3 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "fig, axs = plt.subplots(1, 3, figsize=(20,4))\n",
    "fig.suptitle('Min-Max Normalized Band Histograms')\n",
    "\n",
    "for i, band in enumerate(['red', 'grn', 'blu']):\n",
    "    bins_list1 = rf_scene1.agg(\n",
    "        rf_agg_approx_histogram(band + '_min_max_normalized')['bins'].alias('bins')\n",
    "        ).collect()\n",
    "    values1 = [row['value'] for row in bins_list1[0].bins]\n",
    "    counts1 = [row['count'] for row in bins_list1[0].bins]\n",
    "\n",
    "    bins_list2 = rf_scene2.agg(\n",
    "        rf_agg_approx_histogram(band + '_min_max_normalized')['bins'].alias('bins')\n",
    "        ).collect()\n",
    "    values2 = [row['value'] for row in bins_list2[0].bins]\n",
    "    counts2 = [row['count'] for row in bins_list2[0].bins]\n",
    "\n",
    "    axs[i].hist(values1, weights=counts1, bins = bins, alpha = 0.5, label='Scene_2015')\n",
    "    axs[i].hist(values2, weights=counts2, bins = bins, alpha = 0.5, label='Scene_2017')\n",
    "    axs[i].title.set_text(band)\n",
    "    axs[i].legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The minimum cell value is 0 and the maximum cell value is 100 for each band and scene."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>scene_2015_red_min</th><th>scene_2015_red_max</th><th>scene_2015_grn_min</th><th>scene_2015_grn_max</th><th>scene_2015_blu_min</th><th>scene_2015_blu_max</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>0.0</td><td>100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td>100.0</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/markdown": [
       "| scene_2015_red_min | scene_2015_red_max | scene_2015_grn_min | scene_2015_grn_max | scene_2015_blu_min | scene_2015_blu_max |\n",
       "|---|---|---|---|---|---|\n",
       "| 0.0 | 100.0 | 0.0 | 100.0 | 0.0 | 100.0 |"
      ],
      "text/plain": [
       "DataFrame[scene_2015_red_min: double, scene_2015_red_max: double, scene_2015_grn_min: double, scene_2015_grn_max: double, scene_2015_blu_min: double, scene_2015_blu_max: double]"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_scene1.agg(rf_agg_stats('red_min_max_normalized').alias('red_stats'),\n",
    "              rf_agg_stats('grn_min_max_normalized').alias('grn_stats'),\n",
    "              rf_agg_stats('blu_min_max_normalized').alias('blu_stats')) \\\n",
    "         .select(F.round('red_stats.min', 0).alias('scene_2015_red_min'), \n",
    "                 F.round('red_stats.max', 0).alias('scene_2015_red_max'), \n",
    "                 F.round('grn_stats.min', 0).alias('scene_2015_grn_min'),\n",
    "                 F.round('grn_stats.max', 0).alias('scene_2015_grn_max'),\n",
    "                 F.round('blu_stats.min', 0).alias('scene_2015_blu_min'),\n",
    "                 F.round('blu_stats.max', 0).alias('scene_2015_blu_max'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<table>\n",
       "<thead>\n",
       "<tr><th>scene_2017_red_min</th><th>scene_2017_red_max</th><th>scene_2017_grn_min</th><th>scene_2017_grn_max</th><th>scene_2017_blu_min</th><th>scene_2017_blu_max</th></tr>\n",
       "</thead>\n",
       "<tbody>\n",
       "<tr><td>0.0</td><td>100.0</td><td>0.0</td><td>100.0</td><td>0.0</td><td>100.0</td></tr>\n",
       "</tbody>\n",
       "</table>"
      ],
      "text/markdown": [
       "| scene_2017_red_min | scene_2017_red_max | scene_2017_grn_min | scene_2017_grn_max | scene_2017_blu_min | scene_2017_blu_max |\n",
       "|---|---|---|---|---|---|\n",
       "| 0.0 | 100.0 | 0.0 | 100.0 | 0.0 | 100.0 |"
      ],
      "text/plain": [
       "DataFrame[scene_2017_red_min: double, scene_2017_red_max: double, scene_2017_grn_min: double, scene_2017_grn_max: double, scene_2017_blu_min: double, scene_2017_blu_max: double]"
      ]
     },
     "execution_count": 21,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "rf_scene2.agg(rf_agg_stats('red_min_max_normalized').alias('red_stats'),\n",
    "              rf_agg_stats('grn_min_max_normalized').alias('grn_stats'),\n",
    "              rf_agg_stats('blu_min_max_normalized').alias('blu_stats')) \\\n",
    "         .select(F.round('red_stats.min', 0).alias('scene_2017_red_min'), \n",
    "                 F.round('red_stats.max', 0).alias('scene_2017_red_max'),\n",
    "                 F.round('grn_stats.min', 0).alias('scene_2017_grn_min'),\n",
    "                 F.round('grn_stats.max', 0).alias('scene_2017_grn_max'),\n",
    "                 F.round('blu_stats.min', 0).alias('scene_2017_blu_min'),\n",
    "                 F.round('blu_stats.max', 0).alias('scene_2017_blu_max'))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "EarthAI Environment",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.9"
  },
  "zendesk": {
   "draft": true,
   "id": 360051360531,
   "section_id": 360008877912,
   "title": "How to Normalize the Color Contrast between Two Images"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}