train_model.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Workflow of model training."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import s2spy.time\n",
    "import s2spy.rgdr\n",
    "from s2spy import RGDR"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prepare data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# configure path to the testing data\n",
    "file_path = '/home/yangliu/AI4S2S/s2spy/tests/test_rgdr/test_data'\n",
    "field = xr.open_dataset(f'{file_path}/sst_daily_1979-2018_5deg_Pacific_175_240E_25_50N.nc')\n",
    "target = xr.open_dataset(f'{file_path}/tf5_nc5_dendo_80d77.nc')\n",
    "\n",
    "cal = s2spy.time.AdventCalendar((8, 31), freq = \"30d\")\n",
    "cal = cal.map_to_data(field)\n",
    "field_resampled = s2spy.time.resample(cal, field)\n",
    "target_resampled = s2spy.time.resample(cal, target)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# select data of interest\n",
    "target_timeseries = target_resampled.sel(cluster=3).ts.isel(i_interval=0)\n",
    "precursor_field = field_resampled.sst.isel(i_interval=slice(1,5)) # Multiple lags: 1 through 4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Split data as training set (75%, 30 years) and testing set (25%, 9 years)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# training set\n",
    "train_year_start = 1980\n",
    "train_year_end = 2009\n",
    "precursor_field_train = precursor_field.sel(anchor_year=slice(train_year_start, train_year_end))\n",
    "target_timeseries_train = target_timeseries.sel(anchor_year=slice(train_year_start, train_year_end))\n",
    "\n",
    "# testing set\n",
    "test_year_start = 2010\n",
    "test_year_end = 2018\n",
    "precursor_field_test = precursor_field.sel(anchor_year=slice(test_year_start, test_year_end))\n",
    "target_timeseries_test = target_timeseries.sel(anchor_year=slice(test_year_start, test_year_end))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Perform dimensionality reduction for precursor field."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# configure RGDR\n",
    "rgdr = RGDR(eps_km=600, alpha=0.05, min_area_km2=3000**2)\n",
    "# perform RGDR\n",
    "rgdr.fit(precursor_field_train, target_timeseries_train)\n",
    "clustered_data_train = rgdr.transform(precursor_field_train)\n",
    "clustered_data_test = rgdr.transform(precursor_field_test)\n",
    "# transform dimensions for model training\n",
    "clustered_data_train = clustered_data_train.transpose(\"anchor_year\", \"cluster_labels\")\n",
    "clustered_data_test = clustered_data_test.transpose(\"anchor_year\", \"cluster_labels\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div><svg style=\"position: absolute; width: 0; height: 0; overflow: hidden\">\n",
       "<defs>\n",
       "<symbol id=\"icon-database\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M16 0c-8.837 0-16 2.239-16 5v4c0 2.761 7.163 5 16 5s16-2.239 16-5v-4c0-2.761-7.163-5-16-5z\"></path>\n",
       "<path d=\"M16 17c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "<path d=\"M16 26c-8.837 0-16-2.239-16-5v6c0 2.761 7.163 5 16 5s16-2.239 16-5v-6c0 2.761-7.163 5-16 5z\"></path>\n",
       "</symbol>\n",
       "<symbol id=\"icon-file-text2\" viewBox=\"0 0 32 32\">\n",
       "<path d=\"M28.681 7.159c-0.694-0.947-1.662-2.053-2.724-3.116s-2.169-2.030-3.116-2.724c-1.612-1.182-2.393-1.319-2.841-1.319h-15.5c-1.378 0-2.5 1.121-2.5 2.5v27c0 1.378 1.122 2.5 2.5 2.5h23c1.378 0 2.5-1.122 2.5-2.5v-19.5c0-0.448-0.137-1.23-1.319-2.841zM24.543 5.457c0.959 0.959 1.712 1.825 2.268 2.543h-4.811v-4.811c0.718 0.556 1.584 1.309 2.543 2.268zM28 29.5c0 0.271-0.229 0.5-0.5 0.5h-23c-0.271 0-0.5-0.229-0.5-0.5v-27c0-0.271 0.229-0.5 0.5-0.5 0 0 15.499-0 15.5 0v7c0 0.552 0.448 1 1 1h7v19.5z\"></path>\n",
       "<path d=\"M23 26h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 22h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "<path d=\"M23 18h-14c-0.552 0-1-0.448-1-1s0.448-1 1-1h14c0.552 0 1 0.448 1 1s-0.448 1-1 1z\"></path>\n",
       "</symbol>\n",
       "</defs>\n",
       "</svg>\n",
       "<style>/* CSS stylesheet for displaying xarray objects in jupyterlab.\n",
       " *\n",
       " */\n",
       "\n",
       ":root {\n",
       "  --xr-font-color0: var(--jp-content-font-color0, rgba(0, 0, 0, 1));\n",
       "  --xr-font-color2: var(--jp-content-font-color2, rgba(0, 0, 0, 0.54));\n",
       "  --xr-font-color3: var(--jp-content-font-color3, rgba(0, 0, 0, 0.38));\n",
       "  --xr-border-color: var(--jp-border-color2, #e0e0e0);\n",
       "  --xr-disabled-color: var(--jp-layout-color3, #bdbdbd);\n",
       "  --xr-background-color: var(--jp-layout-color0, white);\n",
       "  --xr-background-color-row-even: var(--jp-layout-color1, white);\n",
       "  --xr-background-color-row-odd: var(--jp-layout-color2, #eeeeee);\n",
       "}\n",
       "\n",
       "html[theme=dark],\n",
       "body.vscode-dark {\n",
       "  --xr-font-color0: rgba(255, 255, 255, 1);\n",
       "  --xr-font-color2: rgba(255, 255, 255, 0.54);\n",
       "  --xr-font-color3: rgba(255, 255, 255, 0.38);\n",
       "  --xr-border-color: #1F1F1F;\n",
       "  --xr-disabled-color: #515151;\n",
       "  --xr-background-color: #111111;\n",
       "  --xr-background-color-row-even: #111111;\n",
       "  --xr-background-color-row-odd: #313131;\n",
       "}\n",
       "\n",
       ".xr-wrap {\n",
       "  display: block;\n",
       "  min-width: 300px;\n",
       "  max-width: 700px;\n",
       "}\n",
       "\n",
       ".xr-text-repr-fallback {\n",
       "  /* fallback to plain text repr when CSS is not injected (untrusted notebook) */\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-header {\n",
       "  padding-top: 6px;\n",
       "  padding-bottom: 6px;\n",
       "  margin-bottom: 4px;\n",
       "  border-bottom: solid 1px var(--xr-border-color);\n",
       "}\n",
       "\n",
       ".xr-header > div,\n",
       ".xr-header > ul {\n",
       "  display: inline;\n",
       "  margin-top: 0;\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-obj-type,\n",
       ".xr-array-name {\n",
       "  margin-left: 2px;\n",
       "  margin-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-obj-type {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-sections {\n",
       "  padding-left: 0 !important;\n",
       "  display: grid;\n",
       "  grid-template-columns: 150px auto auto 1fr 20px 20px;\n",
       "}\n",
       "\n",
       ".xr-section-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-section-item input {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-item input + label {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label {\n",
       "  cursor: pointer;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-item input:enabled + label:hover {\n",
       "  color: var(--xr-font-color0);\n",
       "}\n",
       "\n",
       ".xr-section-summary {\n",
       "  grid-column: 1;\n",
       "  color: var(--xr-font-color2);\n",
       "  font-weight: 500;\n",
       "}\n",
       "\n",
       ".xr-section-summary > span {\n",
       "  display: inline-block;\n",
       "  padding-left: 0.5em;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label {\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in + label:before {\n",
       "  display: inline-block;\n",
       "  content: '►';\n",
       "  font-size: 11px;\n",
       "  width: 15px;\n",
       "  text-align: center;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:disabled + label:before {\n",
       "  color: var(--xr-disabled-color);\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label:before {\n",
       "  content: '▼';\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked + label > span {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-section-summary,\n",
       ".xr-section-inline-details {\n",
       "  padding-top: 4px;\n",
       "  padding-bottom: 4px;\n",
       "}\n",
       "\n",
       ".xr-section-inline-details {\n",
       "  grid-column: 2 / -1;\n",
       "}\n",
       "\n",
       ".xr-section-details {\n",
       "  display: none;\n",
       "  grid-column: 1 / -1;\n",
       "  margin-bottom: 5px;\n",
       "}\n",
       "\n",
       ".xr-section-summary-in:checked ~ .xr-section-details {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-array-wrap {\n",
       "  grid-column: 1 / -1;\n",
       "  display: grid;\n",
       "  grid-template-columns: 20px auto;\n",
       "}\n",
       "\n",
       ".xr-array-wrap > label {\n",
       "  grid-column: 1;\n",
       "  vertical-align: top;\n",
       "}\n",
       "\n",
       ".xr-preview {\n",
       "  color: var(--xr-font-color3);\n",
       "}\n",
       "\n",
       ".xr-array-preview,\n",
       ".xr-array-data {\n",
       "  padding: 0 5px !important;\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-array-data,\n",
       ".xr-array-in:checked ~ .xr-array-preview {\n",
       "  display: none;\n",
       "}\n",
       "\n",
       ".xr-array-in:checked ~ .xr-array-data,\n",
       ".xr-array-preview {\n",
       "  display: inline-block;\n",
       "}\n",
       "\n",
       ".xr-dim-list {\n",
       "  display: inline-block !important;\n",
       "  list-style: none;\n",
       "  padding: 0 !important;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list li {\n",
       "  display: inline-block;\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "}\n",
       "\n",
       ".xr-dim-list:before {\n",
       "  content: '(';\n",
       "}\n",
       "\n",
       ".xr-dim-list:after {\n",
       "  content: ')';\n",
       "}\n",
       "\n",
       ".xr-dim-list li:not(:last-child):after {\n",
       "  content: ',';\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-has-index {\n",
       "  font-weight: bold;\n",
       "}\n",
       "\n",
       ".xr-var-list,\n",
       ".xr-var-item {\n",
       "  display: contents;\n",
       "}\n",
       "\n",
       ".xr-var-item > div,\n",
       ".xr-var-item label,\n",
       ".xr-var-item > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-even);\n",
       "  margin-bottom: 0;\n",
       "}\n",
       "\n",
       ".xr-var-item > .xr-var-name:hover span {\n",
       "  padding-right: 5px;\n",
       "}\n",
       "\n",
       ".xr-var-list > li:nth-child(odd) > div,\n",
       ".xr-var-list > li:nth-child(odd) > label,\n",
       ".xr-var-list > li:nth-child(odd) > .xr-var-name span {\n",
       "  background-color: var(--xr-background-color-row-odd);\n",
       "}\n",
       "\n",
       ".xr-var-name {\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-var-dims {\n",
       "  grid-column: 2;\n",
       "}\n",
       "\n",
       ".xr-var-dtype {\n",
       "  grid-column: 3;\n",
       "  text-align: right;\n",
       "  color: var(--xr-font-color2);\n",
       "}\n",
       "\n",
       ".xr-var-preview {\n",
       "  grid-column: 4;\n",
       "}\n",
       "\n",
       ".xr-var-name,\n",
       ".xr-var-dims,\n",
       ".xr-var-dtype,\n",
       ".xr-preview,\n",
       ".xr-attrs dt {\n",
       "  white-space: nowrap;\n",
       "  overflow: hidden;\n",
       "  text-overflow: ellipsis;\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-var-name:hover,\n",
       ".xr-var-dims:hover,\n",
       ".xr-var-dtype:hover,\n",
       ".xr-attrs dt:hover {\n",
       "  overflow: visible;\n",
       "  width: auto;\n",
       "  z-index: 1;\n",
       "}\n",
       "\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  display: none;\n",
       "  background-color: var(--xr-background-color) !important;\n",
       "  padding-bottom: 5px !important;\n",
       "}\n",
       "\n",
       ".xr-var-attrs-in:checked ~ .xr-var-attrs,\n",
       ".xr-var-data-in:checked ~ .xr-var-data {\n",
       "  display: block;\n",
       "}\n",
       "\n",
       ".xr-var-data > table {\n",
       "  float: right;\n",
       "}\n",
       "\n",
       ".xr-var-name span,\n",
       ".xr-var-data,\n",
       ".xr-attrs {\n",
       "  padding-left: 25px !important;\n",
       "}\n",
       "\n",
       ".xr-attrs,\n",
       ".xr-var-attrs,\n",
       ".xr-var-data {\n",
       "  grid-column: 1 / -1;\n",
       "}\n",
       "\n",
       "dl.xr-attrs {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  display: grid;\n",
       "  grid-template-columns: 125px auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt,\n",
       ".xr-attrs dd {\n",
       "  padding: 0;\n",
       "  margin: 0;\n",
       "  float: left;\n",
       "  padding-right: 10px;\n",
       "  width: auto;\n",
       "}\n",
       "\n",
       ".xr-attrs dt {\n",
       "  font-weight: normal;\n",
       "  grid-column: 1;\n",
       "}\n",
       "\n",
       ".xr-attrs dt:hover span {\n",
       "  display: inline-block;\n",
       "  background: var(--xr-background-color);\n",
       "  padding-right: 10px;\n",
       "}\n",
       "\n",
       ".xr-attrs dd {\n",
       "  grid-column: 2;\n",
       "  white-space: pre-wrap;\n",
       "  word-break: break-all;\n",
       "}\n",
       "\n",
       ".xr-icon-database,\n",
       ".xr-icon-file-text2 {\n",
       "  display: inline-block;\n",
       "  vertical-align: middle;\n",
       "  width: 1em;\n",
       "  height: 1.5em !important;\n",
       "  stroke-width: 0;\n",
       "  stroke: currentColor;\n",
       "  fill: currentColor;\n",
       "}\n",
       "</style><pre class='xr-text-repr-fallback'>&lt;xarray.DataArray &#x27;sst&#x27; (anchor_year: 30, cluster_labels: 7)&gt;\n",
       "array([[286.91668162, 299.84783596, 289.93774506, 297.54843947,\n",
       "        278.24066238, 286.10327492, 284.25688101],\n",
       "       [287.01053665, 299.90999519, 290.99665391, 297.70846305,\n",
       "        278.28356247, 286.84083223, 284.66730471],\n",
       "       [286.54620289, 298.42320455, 289.64697831, 295.4500786 ,\n",
       "        279.14810486, 286.19085251, 284.32519283],\n",
       "       [285.76642198, 299.85415344, 289.79958178, 297.46877566,\n",
       "        278.55730413, 285.7835013 , 283.63555146],\n",
       "       [287.28029452, 300.53261312, 290.32029845, 298.61729337,\n",
       "        278.91383413, 286.93986272, 284.47967082],\n",
       "       [287.11752036, 298.77104119, 290.61151905, 297.72280291,\n",
       "        278.97152786, 286.79347391, 284.8633425 ],\n",
       "       [286.91728331, 299.26571418, 290.95184425, 296.86022156,\n",
       "        279.06490962, 286.69249845, 285.0038166 ],\n",
       "       [285.73728271, 299.58077257, 290.1544995 , 295.95557241,\n",
       "        278.3467481 , 286.25965795, 284.43501001],\n",
       "       [286.27027317, 300.35477566, 290.1053183 , 298.7395288 ,\n",
       "        278.16248144, 285.89035   , 284.09578366],\n",
       "       [287.12407151, 299.47631632, 290.59468686, 297.39734454,\n",
       "        279.18296738, 286.86859677, 284.88719095],\n",
       "...\n",
       "       [287.51545001, 299.90851373, 290.34061042, 298.11604053,\n",
       "        278.68831329, 285.86874303, 283.83763964],\n",
       "       [287.4697688 , 300.19267917, 290.31979773, 298.33634933,\n",
       "        278.50768102, 286.32706228, 284.13051603],\n",
       "       [287.33929671, 300.05178019, 290.21811635, 297.50691676,\n",
       "        279.52856318, 286.8266357 , 284.81265618],\n",
       "       [287.259754  , 299.29538133, 290.75069599, 297.42936244,\n",
       "        278.6126709 , 286.22619978, 284.37959048],\n",
       "       [288.19457254, 299.91475186, 291.53432816, 296.19082513,\n",
       "        279.2724528 , 287.59860329, 284.88366148],\n",
       "       [287.73373685, 299.87610711, 290.69421445, 298.0900667 ,\n",
       "        279.67929255, 287.21957964, 285.26486288],\n",
       "       [286.85064034, 299.8528795 , 290.42994901, 297.67371479,\n",
       "        278.86642609, 285.85201269, 283.9566866 ],\n",
       "       [286.16938397, 299.89861688, 290.5898762 , 297.87487282,\n",
       "        278.4716746 , 286.55275645, 284.31180943],\n",
       "       [287.67270923, 299.32189297, 290.13167005, 297.69282282,\n",
       "        278.98851878, 286.96470686, 284.75600036],\n",
       "       [288.75096715, 300.44966058, 291.05704933, 297.29816477,\n",
       "        279.8935379 , 287.42835661, 285.22563649]])\n",
       "Coordinates:\n",
       "  * anchor_year     (anchor_year) int64 1980 1981 1982 1983 ... 2007 2008 2009\n",
       "  * cluster_labels  (cluster_labels) &lt;U20 &#x27;lag:1_cluster:-1&#x27; ... &#x27;lag:4_clust...\n",
       "    latitude        (cluster_labels) float64 43.89 27.5 34.45 ... 38.79 40.26\n",
       "    longitude       (cluster_labels) float64 213.4 182.5 226.3 ... 220.7 217.6</pre><div class='xr-wrap' hidden><div class='xr-header'><div class='xr-obj-type'>xarray.DataArray</div><div class='xr-array-name'>'sst'</div><ul class='xr-dim-list'><li><span class='xr-has-index'>anchor_year</span>: 30</li><li><span class='xr-has-index'>cluster_labels</span>: 7</li></ul></div><ul class='xr-sections'><li class='xr-section-item'><div class='xr-array-wrap'><input id='section-5c0a044c-1561-4f82-a824-4c6b7f9d223e' class='xr-array-in' type='checkbox' checked><label for='section-5c0a044c-1561-4f82-a824-4c6b7f9d223e' title='Show/hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-array-preview xr-preview'><span>286.9 299.8 289.9 297.5 278.2 286.1 ... 291.1 297.3 279.9 287.4 285.2</span></div><div class='xr-array-data'><pre>array([[286.91668162, 299.84783596, 289.93774506, 297.54843947,\n",
       "        278.24066238, 286.10327492, 284.25688101],\n",
       "       [287.01053665, 299.90999519, 290.99665391, 297.70846305,\n",
       "        278.28356247, 286.84083223, 284.66730471],\n",
       "       [286.54620289, 298.42320455, 289.64697831, 295.4500786 ,\n",
       "        279.14810486, 286.19085251, 284.32519283],\n",
       "       [285.76642198, 299.85415344, 289.79958178, 297.46877566,\n",
       "        278.55730413, 285.7835013 , 283.63555146],\n",
       "       [287.28029452, 300.53261312, 290.32029845, 298.61729337,\n",
       "        278.91383413, 286.93986272, 284.47967082],\n",
       "       [287.11752036, 298.77104119, 290.61151905, 297.72280291,\n",
       "        278.97152786, 286.79347391, 284.8633425 ],\n",
       "       [286.91728331, 299.26571418, 290.95184425, 296.86022156,\n",
       "        279.06490962, 286.69249845, 285.0038166 ],\n",
       "       [285.73728271, 299.58077257, 290.1544995 , 295.95557241,\n",
       "        278.3467481 , 286.25965795, 284.43501001],\n",
       "       [286.27027317, 300.35477566, 290.1053183 , 298.7395288 ,\n",
       "        278.16248144, 285.89035   , 284.09578366],\n",
       "       [287.12407151, 299.47631632, 290.59468686, 297.39734454,\n",
       "        279.18296738, 286.86859677, 284.88719095],\n",
       "...\n",
       "       [287.51545001, 299.90851373, 290.34061042, 298.11604053,\n",
       "        278.68831329, 285.86874303, 283.83763964],\n",
       "       [287.4697688 , 300.19267917, 290.31979773, 298.33634933,\n",
       "        278.50768102, 286.32706228, 284.13051603],\n",
       "       [287.33929671, 300.05178019, 290.21811635, 297.50691676,\n",
       "        279.52856318, 286.8266357 , 284.81265618],\n",
       "       [287.259754  , 299.29538133, 290.75069599, 297.42936244,\n",
       "        278.6126709 , 286.22619978, 284.37959048],\n",
       "       [288.19457254, 299.91475186, 291.53432816, 296.19082513,\n",
       "        279.2724528 , 287.59860329, 284.88366148],\n",
       "       [287.73373685, 299.87610711, 290.69421445, 298.0900667 ,\n",
       "        279.67929255, 287.21957964, 285.26486288],\n",
       "       [286.85064034, 299.8528795 , 290.42994901, 297.67371479,\n",
       "        278.86642609, 285.85201269, 283.9566866 ],\n",
       "       [286.16938397, 299.89861688, 290.5898762 , 297.87487282,\n",
       "        278.4716746 , 286.55275645, 284.31180943],\n",
       "       [287.67270923, 299.32189297, 290.13167005, 297.69282282,\n",
       "        278.98851878, 286.96470686, 284.75600036],\n",
       "       [288.75096715, 300.44966058, 291.05704933, 297.29816477,\n",
       "        279.8935379 , 287.42835661, 285.22563649]])</pre></div></div></li><li class='xr-section-item'><input id='section-eee5bb3a-e036-4b1d-8804-972ebe4c42ec' class='xr-section-summary-in' type='checkbox'  checked><label for='section-eee5bb3a-e036-4b1d-8804-972ebe4c42ec' class='xr-section-summary' >Coordinates: <span>(4)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><ul class='xr-var-list'><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>anchor_year</span></div><div class='xr-var-dims'>(anchor_year)</div><div class='xr-var-dtype'>int64</div><div class='xr-var-preview xr-preview'>1980 1981 1982 ... 2007 2008 2009</div><input id='attrs-0d856a2a-64a9-4d0f-b86a-a5840b798182' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-0d856a2a-64a9-4d0f-b86a-a5840b798182' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-2154de21-be63-429d-8b7a-9f00af64cf74' class='xr-var-data-in' type='checkbox'><label for='data-2154de21-be63-429d-8b7a-9f00af64cf74' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([1980, 1981, 1982, 1983, 1984, 1985, 1986, 1987, 1988, 1989, 1990, 1991,\n",
       "       1992, 1993, 1994, 1995, 1996, 1997, 1998, 1999, 2000, 2001, 2002, 2003,\n",
       "       2004, 2005, 2006, 2007, 2008, 2009])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span class='xr-has-index'>cluster_labels</span></div><div class='xr-var-dims'>(cluster_labels)</div><div class='xr-var-dtype'>&lt;U20</div><div class='xr-var-preview xr-preview'>&#x27;lag:1_cluster:-1&#x27; ... &#x27;lag:4_cl...</div><input id='attrs-b5d35826-4bcf-4f24-b5c6-4ff79fcf6c17' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-b5d35826-4bcf-4f24-b5c6-4ff79fcf6c17' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-0c238645-19d0-4e89-82a9-e2042d3e3b68' class='xr-var-data-in' type='checkbox'><label for='data-0c238645-19d0-4e89-82a9-e2042d3e3b68' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([&#x27;lag:1_cluster:-1&#x27;, &#x27;lag:1_cluster:1&#x27;, &#x27;lag:2_cluster:-2&#x27;,\n",
       "       &#x27;lag:2_cluster:1&#x27;, &#x27;lag:3_cluster:-1&#x27;, &#x27;lag:3_cluster:-2&#x27;,\n",
       "       &#x27;lag:4_cluster:-2&#x27;], dtype=&#x27;&lt;U20&#x27;)</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>latitude</span></div><div class='xr-var-dims'>(cluster_labels)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>43.89 27.5 34.45 ... 38.79 40.26</div><input id='attrs-6a386480-e84a-4459-8e92-304d7a39603a' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-6a386480-e84a-4459-8e92-304d7a39603a' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-47fa4934-2316-457f-8ae5-90bdbde8fa0f' class='xr-var-data-in' type='checkbox'><label for='data-47fa4934-2316-457f-8ae5-90bdbde8fa0f' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([43.89088777, 27.5       , 34.44976166, 28.46026867, 47.5       ,\n",
       "       38.78952959, 40.26462677])</pre></div></li><li class='xr-var-item'><div class='xr-var-name'><span>longitude</span></div><div class='xr-var-dims'>(cluster_labels)</div><div class='xr-var-dtype'>float64</div><div class='xr-var-preview xr-preview'>213.4 182.5 226.3 ... 220.7 217.6</div><input id='attrs-2a1729cb-d722-43b9-ac42-42138b493a12' class='xr-var-attrs-in' type='checkbox' disabled><label for='attrs-2a1729cb-d722-43b9-ac42-42138b493a12' title='Show/Hide attributes'><svg class='icon xr-icon-file-text2'><use xlink:href='#icon-file-text2'></use></svg></label><input id='data-109965f2-6c1c-429f-80db-850cd3123ed7' class='xr-var-data-in' type='checkbox'><label for='data-109965f2-6c1c-429f-80db-850cd3123ed7' title='Show/Hide data repr'><svg class='icon xr-icon-database'><use xlink:href='#icon-database'></use></svg></label><div class='xr-var-attrs'><dl class='xr-attrs'></dl></div><div class='xr-var-data'><pre>array([213.44261605, 182.5       , 226.30005617, 186.440403  ,\n",
       "       195.        , 220.66657355, 217.61060762])</pre></div></li></ul></div></li><li class='xr-section-item'><input id='section-be5c0ec2-b2bc-43ed-9e5e-f499d84ec2a7' class='xr-section-summary-in' type='checkbox' disabled ><label for='section-be5c0ec2-b2bc-43ed-9e5e-f499d84ec2a7' class='xr-section-summary'  title='Expand/collapse section'>Attributes: <span>(0)</span></label><div class='xr-section-inline-details'></div><div class='xr-section-details'><dl class='xr-attrs'></dl></div></li></ul></div></div>"
      ],
      "text/plain": [
       "<xarray.DataArray 'sst' (anchor_year: 30, cluster_labels: 7)>\n",
       "array([[286.91668162, 299.84783596, 289.93774506, 297.54843947,\n",
       "        278.24066238, 286.10327492, 284.25688101],\n",
       "       [287.01053665, 299.90999519, 290.99665391, 297.70846305,\n",
       "        278.28356247, 286.84083223, 284.66730471],\n",
       "       [286.54620289, 298.42320455, 289.64697831, 295.4500786 ,\n",
       "        279.14810486, 286.19085251, 284.32519283],\n",
       "       [285.76642198, 299.85415344, 289.79958178, 297.46877566,\n",
       "        278.55730413, 285.7835013 , 283.63555146],\n",
       "       [287.28029452, 300.53261312, 290.32029845, 298.61729337,\n",
       "        278.91383413, 286.93986272, 284.47967082],\n",
       "       [287.11752036, 298.77104119, 290.61151905, 297.72280291,\n",
       "        278.97152786, 286.79347391, 284.8633425 ],\n",
       "       [286.91728331, 299.26571418, 290.95184425, 296.86022156,\n",
       "        279.06490962, 286.69249845, 285.0038166 ],\n",
       "       [285.73728271, 299.58077257, 290.1544995 , 295.95557241,\n",
       "        278.3467481 , 286.25965795, 284.43501001],\n",
       "       [286.27027317, 300.35477566, 290.1053183 , 298.7395288 ,\n",
       "        278.16248144, 285.89035   , 284.09578366],\n",
       "       [287.12407151, 299.47631632, 290.59468686, 297.39734454,\n",
       "        279.18296738, 286.86859677, 284.88719095],\n",
       "...\n",
       "       [287.51545001, 299.90851373, 290.34061042, 298.11604053,\n",
       "        278.68831329, 285.86874303, 283.83763964],\n",
       "       [287.4697688 , 300.19267917, 290.31979773, 298.33634933,\n",
       "        278.50768102, 286.32706228, 284.13051603],\n",
       "       [287.33929671, 300.05178019, 290.21811635, 297.50691676,\n",
       "        279.52856318, 286.8266357 , 284.81265618],\n",
       "       [287.259754  , 299.29538133, 290.75069599, 297.42936244,\n",
       "        278.6126709 , 286.22619978, 284.37959048],\n",
       "       [288.19457254, 299.91475186, 291.53432816, 296.19082513,\n",
       "        279.2724528 , 287.59860329, 284.88366148],\n",
       "       [287.73373685, 299.87610711, 290.69421445, 298.0900667 ,\n",
       "        279.67929255, 287.21957964, 285.26486288],\n",
       "       [286.85064034, 299.8528795 , 290.42994901, 297.67371479,\n",
       "        278.86642609, 285.85201269, 283.9566866 ],\n",
       "       [286.16938397, 299.89861688, 290.5898762 , 297.87487282,\n",
       "        278.4716746 , 286.55275645, 284.31180943],\n",
       "       [287.67270923, 299.32189297, 290.13167005, 297.69282282,\n",
       "        278.98851878, 286.96470686, 284.75600036],\n",
       "       [288.75096715, 300.44966058, 291.05704933, 297.29816477,\n",
       "        279.8935379 , 287.42835661, 285.22563649]])\n",
       "Coordinates:\n",
       "  * anchor_year     (anchor_year) int64 1980 1981 1982 1983 ... 2007 2008 2009\n",
       "  * cluster_labels  (cluster_labels) <U20 'lag:1_cluster:-1' ... 'lag:4_clust...\n",
       "    latitude        (cluster_labels) float64 43.89 27.5 34.45 ... 38.79 40.26\n",
       "    longitude       (cluster_labels) float64 213.4 182.5 226.3 ... 220.7 217.6"
      ]
     },
     "execution_count": 53,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "clustered_data_train"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Train a gradient boosting machine for regression."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hyperparameters for the boosting machine\n",
    "params = {\n",
    "    \"n_estimators\": 500,\n",
    "    \"max_depth\": 4,\n",
    "    \"min_samples_split\": 5,\n",
    "    \"learning_rate\": 0.01,\n",
    "    \"loss\": \"squared_error\",\n",
    "}"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import ensemble\n",
    "regressor = ensemble.GradientBoostingRegressor(**params)\n",
    "# the shape of input X must be [samples, features]\n",
    "regressor.fit(clustered_data_train, target_timeseries_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# score for testing\n",
    "from sklearn.metrics import mean_squared_error\n",
    "mse = mean_squared_error(target_timeseries_test, regressor.predict(clustered_data_test))\n",
    "print(\"The mean squared error (MSE) on test set: {:.2f}\".format(mse))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot training deviance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAGoCAYAAAATsnHAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABAPElEQVR4nO3deXhU1f3H8feXJIQlILsiyKYIgkDAgCi7sqkILlhxxRXUKkXrVv1p6WK1rVXrUi0qLnXDohVURISqWDcWBUSRfYsgRAUk7CHn98eZhCRMkglkcmcyn9fz3Gdm7ty5+c5V8sm599xzzDmHiIhIrKkSdAEiIiLhKKBERCQmKaBERCQmKaBERCQmKaBERCQmKaBERCQmKaBEYoCZNTOzbDNLCroWkVihgBIpAzNbbWY7zWybmW0xs0/M7BozO6R/S865tc65NOfcvvKqVSTeKaBEyu5M51wtoDlwH3Ab8HSwJYlUPgookYPknNvqnJsCnA+MNLPjzSzVzO43s7VmttHMnjCz6gBmttjMhuR93sySzewHM+tiZi3MzJlZcui9y0PbbzOzlWY2usDn+ppZppn92sw2mdkGM7u8wPvVzexvZrbGzLaa2f8K1NA91OrbYmYLzKxvxRwtkbJTQIkcIufcbCAT6AX8GTgWSAeOAZoAd4c2fRm4oMBHBwE/OOe+CLPbTcAQoDZwOfCgmXUp8P4RwGGh/V8JPGZmdUPv3Q+cAJwM1ANuBXLNrAnwNvDH0PqbgdfMrOHBfneRaFJAiZSP9fhf+lcDNzrnfnLObQP+BIwIbfMSMNTMaoReXxhadwDn3NvOuRXO+xCYjg/APHuB3zvn9jrnpgLZQJvQtbArgF85575zzu1zzn3inNsNXAxMdc5Ndc7lOufeA+YCp5fjcRApN8lBFyBSSTTB/3uqAcwzs7z1BiQBOOeWm9li4EwzexMYCnQOtzMzOw34Lb41ViW0368KbPKjcy6nwOsdQBrQAKgGrAiz2+bAeWZ2ZoF1KcD7kX9NkYqjgBI5RGbWFR9Qb+A7TLR3zn1XzOZ5p/mqAN8455aH2V8q8BpwKTDZObfXzN7Ah11pfgB2AUcDC4q8tw74l3Pu6gj2IxI4neITOUhmVjvU6eEV4AXn3ALgSfz1okahbZqY2aACH3sFGAhcSzGn94CqQCqQBeSEWlMDI6nJOZcLTAAeMLMjzSzJzE4Khd4L+NbboND6aqEOF03L/OVFKoACSqTs3jSzbfgWyZ3AA/iODOBbUMuBz8zsZ2AG0Cbvg865DcCn+A4ME8PtPHTtagzwKrAZf61qShnquxl/OnAO8BO+40YV59w6YBhwBz781gG3oN8DEqNMExaKiEgs0l9OIiISkxRQIiISkxRQIiISkxRQIiISk+LuPqgGDRq4Fi1aBF2GiIiUk3nz5v3gnDtgyK24C6gWLVowd+7coMsQEZFyYmZrwq3XKT4REYlJCigREYlJCigREYlJUbsGZWZHAc/j563JBcY75/5eZJu+wGRgVWjV686530erJhE5NHv37iUzM5Ndu3YFXYrEoWrVqtG0aVNSUlIi2j6anSRygF87574ws1r4KQjec859U2S7j5xzQ8J8XkRiTGZmJrVq1aJFixYUmFJEpFTOOX788UcyMzNp2bJlRJ+J2ik+59yGvJlCQ4NfLsZPSSAicWrXrl3Ur19f4SRlZmbUr1+/TK3vCrkGZWYt8BOzfR7m7ZPMbIGZvWNm7Yv5/Cgzm2tmc7OysqJZqoiUQuEkB6us/+9EPaDMLA0/+dpY59zPRd7+AmjunOsEPIKf8O0AzrnxzrkM51xGw4YH3MslIiKVUFQDysxS8OH0onPu9aLvO+d+ds5lh55PBVLMrEE0axKR+PXjjz+Snp5Oeno6RxxxBE2aNMl/vWfPnhI/O3fuXMaMGVPqzzj55JPLpdYdO3Zw0UUX0aFDB44//nh69uxJdnZ2iZ/505/+VOx7LVq0oEOHDnTo0IF27drxf//3f+zevfugalu/fj3Dhw8/qM9WpKjNB2W+Lfcc8JNzbmwx2xwBbHTOOTPrBkzCt6iKLSojI8NpJAmRYCxevJjjjjsu6DIAGDduHGlpadx8883563JyckhOjo0Bcu69916ysrJ44IEHAFiyZAktWrQgNTW12M+kpaUVG2J5o+g0aNCA7OxsRo0aRUpKCs8991xU6o+WcP8Pmdk851xG0W2j2YLqAVwCnGJm80PL6WZ2jZldE9pmOLDIzBYADwMjSgonEZGiLrvsMm666Sb69evHbbfdxuzZszn55JPp3LkzJ598MkuWLAHggw8+YMgQ32F43LhxXHHFFfTt25dWrVrx8MMP5+8vLS0tf/u+ffsyfPhw2rZty0UXXUTer6epU6fStm1bevbsyZgxY/L3W9CGDRto0mR/v7A2bdrkh9MLL7xAt27dSE9PZ/To0ezbt4/bb7+dnTt3kp6ezkUXXVTid05LS+OJJ57gjTfe4KeffgLgr3/9K127dqVjx4789re/BeC2227jH//4R/7nxo0bx9/+9jdWr17N8ccfD8Dq1avp1asXXbp0oUuXLnzyySelfv85c+Zw8skn06lTJ7p168a2bdvYt28ft9xyS34N//znPyP671eSqP2p4Zz7H1DiFTHn3KPAo9GqQUSiZ+xYmD+/fPeZng4PPVT2zy1dupQZM2aQlJTEzz//zKxZs0hOTmbGjBnccccdvPbaawd85ttvv+X9999n27ZttGnThmuvvfaA+3O+/PJLvv76a4488kh69OjBxx9/TEZGBqNHj2bWrFm0bNmSCy64IGxNV1xxBQMHDmTSpEmceuqpjBw5ktatW7N48WImTpzIxx9/TEpKCtdddx0vvvgi9913H48++ijzIzyotWvXpmXLlixbtoytW7eybNkyZs+ejXOOoUOHMmvWLEaMGMHYsWO57rrrAHj11VeZNm0aubm5+ftp1KgR7733HtWqVWPZsmVccMEF+eOdhvv+3bp14/zzz2fixIl07dqVn3/+merVq/P0009z2GGHMWfOHHbv3k2PHj0YOHBgxF3Kw4mNtrCIyCE477zzSEpKAmDr1q2MHDmSZcuWYWbs3bs37GfOOOMMUlNTSU1NpVGjRmzcuJGmTZsW2qZbt27569LT01m9ejVpaWm0atUq/xfvBRdcwPjx4w/Yf3p6OitXrmT69OnMmDGDrl278umnnzJz5kzmzZtH165dAdi5cyeNGjU6qO+d16KZPn0606dPp3PnzgBkZ2ezbNkyrrzySjZt2sT69evJysqibt26NGvWjNWrV+fvY+/evVx//fXMnz+fpKQkli5dWuL3P+yww2jcuHF+/bVr186vYeHChUyaNAkgPzQVUCJS4Q6mpRMtNWvWzH9+11130a9fP/7zn/+wevVq+vbtG/YzBa8FJSUlkZOTE9E2ZbkKkZaWxjnnnMM555xDlSpVmDp1KlWrVmXkyJHce++9Ee8nnG3btrF69WqOPfZYnHP85je/YfTo0QdsN3z4cCZNmsT333/PiBEjDnj/wQcf5PDDD2fBggXk5uZSrVq1/PeK+/7huos753jkkUcYNGjQIX2vghJuLL4//hF69Ai6ChGJlq1bt+Zf+3n22WfLff9t27Zl5cqV+a2QiRMnht3u448/ZvPmzQDs2bOHb775hubNm3PqqacyadIkNm3aBMBPP/3EmjV+tomUlJRiW3wFZWdnc91113HWWWdRt25dBg0axIQJE/I7WHz33Xf5+x8xYgSvvPIKkyZNCttzb+vWrTRu3JgqVarwr3/9i3379pX6/devX8+cOXMAH5Q5OTkMGjSIxx9/PL/+pUuXsn379lK/S0kSrgW1cyd8/jns2wehMwIiUonceuutjBw5kgceeIBTTjml3PdfvXp1/vGPfzB48GAaNGhAt27dwm63YsUKrr32Wpxz5ObmcsYZZ3DuuediZvzxj39k4MCB5ObmkpKSwmOPPUbz5s0ZNWoUHTt2pEuXLrz44osH7LNfv375+zv77LO56667ABg4cCCLFy/mpJNOAnzL7YUXXqBRo0a0b9+ebdu20aRJExo3bnzAPq+77jrOPfdc/v3vf9OvX79CrdFwqlatysSJE7nhhhvYuXMn1atXZ8aMGVx11VWsXr2aLl264JyjYcOGvPHGG2U8uoVFrZt5tBxqN/MnnoBrr4XMTGiigZdEyiSWupkHKTs7m7S0NJxz/PKXv6R169bceOONQZcVF2Klm3lMOuoo/7huXbB1iEj8evLJJ0lPT6d9+/Zs3bo17LUfOXQJd4ovr5POunXQvXuwtYhIfLrxxhvVYqoAakGJiEhMSriAqlsXatRQQImIxLqECygz34pSQImIxLaECyjwAZWZGXQVIiJSkoQMqKZN1YISiUeHMt0G+AFQ8wZDLWrjxo0MGTKETp060a5dO04//fQS97Vly5ZCA7EWlZSUlN/Tr1OnTjzwwAOFxsAri0inCqlsEq4XH/gW1IYNsHcvFBkbUkRiWP369fMHUw033UZpPvjgA9LS0sLO+XT33XczYMAAfvWrXwGwcOHCEveVF1B5A7EWVb169fxaN23axIUXXsjWrVv53e9+F3G9eTIyMsjIOOA2oUovIVtQRx0FzsH69UFXIiKHat68efTp04cTTjiBQYMGsWHDBgAefvhh2rVrR8eOHRkxYgSrV6/miSee4MEHHyQ9PZ2PPvqo0H42bNhQaLDYjh075j8PN5XF7bffzooVK0hPT+eWW24pscZGjRoxfvx4Hn30UZxzxU5Ncf755zN16tT8z1122WW89tprhaYKKW46kWeffZZzzjmHwYMH07p1a2699db8/UybNo0uXbrQqVMnTj31VAC2b9/OFVdcQdeuXencuTOTJ08u24GvAAnbggJ/mq9582BrEYlbMTDfhnOOG264gcmTJ9OwYUMmTpzInXfeyYQJE7jvvvtYtWoVqampbNmyhTp16nDNNdcU2+r65S9/yfnnn8+jjz5K//79ufzyyznyyCOZPn162Kks7rvvPhYtWhTx9BitWrUiNzeXTZs2MXny5LBTU4wYMYKJEydy+umns2fPHmbOnMnjjz/O559/nr+ftm3bFjudyPz58/nyyy9JTU2lTZs23HDDDVSrVo2rr746f3qQvPmj7rnnHk455RQmTJjAli1b6NatG/379y91qKOKlPABJSLxa/fu3SxatIgBAwYAsG/fvvzx5jp27MhFF13EWWedxVlnnVXqvgYNGsTKlSuZNm0a77zzDp07d2bRokXFTmXRrFmzMtdbcHqMcFNTnHbaaYwZM4bdu3czbdo0evfuTfXq1Qvto6TpRE499VQOO+wwANq1a8eaNWvYvHkzvXv3zp/2ol69evk1TJkyhfvvvx+AXbt2sXbt2pgayiqhA0o9+UQOQQzMt+Gco3379nz66acHvPf2228za9YspkyZwh/+8Ae+/vrrUvdXr149LrzwQi688EKGDBnCrFmzip3KouCcSpFYuXIlSUlJNGrUqMSpKfr27cu7777LxIkTw06GWNJ0ImWdHuO1116jTZs2ZfoeFSkhr0HVrg21aqkFJRLvUlNTycrKyg+ovXv38vXXX5Obm8u6devo168ff/nLX9iyZQvZ2dnUqlWLbdu2hd3Xf//7X3bs2AH4KSRWrFhBs2bNip3KoqR9FZWVlcU111zD9ddfj5mVODXFiBEjeOaZZ/joo4/CBlhZpxM56aST+PDDD1m1ahVA/im+QYMG8cgjj+S36r788suIvktFSsgWFOhmXZHKoEqVKkyaNIkxY8awdetWcnJyGDt2LMceeywXX3wxW7duxTnHjTfeSJ06dTjzzDMZPnw4kydP5pFHHqFXr175+5o3bx7XX389ycnJ5ObmctVVV+XPGhtuKoujjz6aHj16cPzxx3Paaafx17/+tVBtO3fuJD09nb1795KcnMwll1zCTTfdBFDi1BQDBw7k0ksvZejQoVStWvWA71zW6UQaNmzI+PHjOeecc8jNzc2f4v2uu+5i7NixdOzYEeccLVq04K233jqo/w7RknDTbeQZPBh++AHKYVciCUPTbcih0nQbEVALSkQktiV0QG3aBLt3B12JiIiEk9ABBfDdd8HWIRJv4u2ygMSOsv6/k7ABVXDiQhGJTLVq1fjxxx8VUlJmzjl+/PFHqlWrFvFnEroXHyigRMqiadOmZGZmkpWVFXQpEoeqVatWaDip0iigFFAiEUtJSckfkUAk2hL2FF/Nmn52XQWUiEhsStiAAmjWDNauDboKEREJJ6EDqnlzWLMm6CpERCQcBdQaPzeUiIjEloQPqG3bYMuWoCsREZGiEj6gAMo4ar6IiFSAhA6oFi38o65DiYjEnoQOqLwWlAJKRCT2JHRANWgA1asroEREYlFCB5SZupqLiMSqhA4o8Neh1ElCRCT2JHxAqQUlIhKbFFDN4ccfYfv2oCsREZGCFFDqySciEpMSPqDy7oXSdSgRkdiS8AGlFpSISGxK+IBq3BhSUhRQIiKxJuEDqkoVP7uuAkpEJLYkfECBvw6lgBIRiS0KKPx1KHWSEBGJLQoofEBt2AC7dwddiYiI5FFAsb8n37p1wdYhIiL7KaDQvFAiIrFIAcX+FtSqVcHWISIi+ymg8N3Mk5IUUCIisUQBBSQn+1aUAkpEJHYooEJatYKVK4OuQkRE8iigQlq2VECJiMQSBVRIq1aQlQXZ2UFXIiIioIDK16qVf9R1KBGR2KCACskLKJ3mExGJDQqoEAWUiEhsiVpAmdlRZva+mS02s6/N7FdhtjEze9jMlpvZQjPrEq16SlO3LtSurYASEYkVyVHcdw7wa+fcF2ZWC5hnZu85574psM1pQOvQciLweOixwpmpq7mISCyJWgvKObfBOfdF6Pk2YDHQpMhmw4DnnfcZUMfMGkerptK0aqVOEiIisaJCrkGZWQugM/B5kbeaAAXHEM/kwBCrMHkBlZsbVAUiIpIn6gFlZmnAa8BY59zPRd8O8xEXZh+jzGyumc3NysqKRpmAD6hdu+D776P2I0REJEJRDSgzS8GH04vOudfDbJIJHFXgdVNgfdGNnHPjnXMZzrmMhg0bRqdY/GgSoOtQIiKxIJq9+Ax4GljsnHugmM2mAJeGevN1B7Y65zZEq6bS6GZdEZHYEc1efD2AS4CvzGx+aN0dQDMA59wTwFTgdGA5sAO4PIr1lKp5c9+bTy0oEZHgRS2gnHP/I/w1poLbOOCX0aqhrFJToWlTBZSISCzQSBJFaFRzEZHYoIAqQjfriojEBgVUEa1awfr1sHNn0JWIiCQ2BVQRRx/tH9WTT0QkWAqoIo45xj8uWxZsHSIiiU4BVUTr1v5x+fJg6xARSXQKqCLq1oX69dWCEhEJmgIqjGOOUUCJiARNARVG69Y6xSciEjQFVBjHHAPr1vmRzUVEJBgKqDBatwbnYMWKoCsREUlcCqgw1JNPRCR4CqgwdC+UiEjwFFBh5HU1VwtKRCQ4CqhitG6tFpSISJAUUMXQvVAiIsFSQBWjdWvf1VyjmouIBEMBVYy8jhKaG0pEJBgKqGLkdTXXaT4RkWAooIqR14JSTz4RkWAooIqhUc1FRIKlgCqBupqLiARHAVWCY4+FpUuDrkJEJDEpoErQpg189x1s2xZ0JSIiiUcBVYK2bf2jWlEiIhVPAVWCvID69ttg6xARSUQKqBIcfTRUqQJLlgRdiYhI4lFAlSA1FVq1UgtKRCQICqhStGmjFpSISBAUUKVo29Z3kti3L+hKREQSiwKqFG3bwq5dsHZt0JWIiCQWBVQp2rTxjzrNJyJSsRIzoMpwvk5dzUVEgpF4AXX33ZCRAc5FtHmDBn7gWLWgREQqVuIF1JFHwvz5sHBhRJub+VaUWlAiIhUr8QJq+HBIToaXX474IwooEZGKl3gB1aABDBgAr7wCubkRfaRNG/j+e9i6Ncq1iYhIvsQLKIALLoA1a+DTTyPaPK+jhK5DiYhUnMQMqLPOgmrVIj7Np4ASEal4iRlQtWrBkCHw739DTk6pm7dq5S9bLV5cAbWJiAiQqAEFcOGFsGkT/Pe/pW6akuKnf//mmwqoS0REgEQOqNNOg9q14cUXI9q8fXv4+uso1yQiIvkSN6CqVYNf/AImTYpoTvf27WHFCti5swJqExGRBA4ogCuugB07/LWoUrRv7wef0P1QIiIVI7EDqnt3f5PThAmlbtqunX/UdSgRkYqR2AFlBpdfDh9/7Cd9KkHr1r4nn65DiYhUjMQOKIBLLoEqVeDZZ0vcrGpVOPZYBZSISEVRQB15JAweDM8/X+o0HOrJJyJScRRQ4DtLfPcdTJ9e4mbt28PKlb5fhYhIwsrNhe3bfQ/oPXui9mOSo7bneHLmmdCoEfzzn/7+qGIU7MnXpUsF1icilZtzflSbffsOXPbu3b+E2yYnxy952zi3f8nbd8HnBdflfTbv8zk5sHu3/yv855/9YAZZWfDTT37ZvNk/Zmfvr/2ZZ+Cyy6JyWBRQ4C8wXXEF/OUvkJkJTZuG3axgTz4FlEg5yMmBLVv8L8Nt2+CHH2DjRj99wMaN/nVuru/Q5Jy/EXH3bti1a//jrl1+G+cKP+bm+v3v3Fn4r3yz/Y8Fn0f6aOavWxdd8mosrpbcXB8oeY85Ob6uPXsiGnKtwplBvXrQsKGfBaJZM0hP9zO41qkDNWr4752REbUSFFB5rr4a/vxneOopGDcu7CatW/thj3QdShLSvn0+RNau9ee6d+70v7zq1/e/YFev9nPS5IXN5s0+RLKy/C/lTZv8X+Z5f6Fv3lzyHDZVq/r9JyXtD6nq1f1N9qmpfqlZ0//8pKT9wVHwMTnZf6Zq1f0BAsW3KiJ5DBc8eUvRn5/3aOZrzFuqVPG1Va26f0lOLrxN3pKSsn8Jt01y8v73UlL8vuHAAA4XtgU/m/e8alV/zGrW9OsCpIDK06oVDBoETz4Jd97p/0MVkZKinnwSR/bt82cEtmzxv8j27fPBsGcPNGnin1er5h9zcmDRIh8aK1f67Z2D9et94CxZ4lsz27dH/vPT0vz+GzTwvzQPP9z/9Z2a6n8B1q3r/0KvWxcOO8wP4tyggd/uiCP8urxfppKQFFAFXXstDBsGb70FZ58ddpP27WHu3AquSyon53znnM2boXFjHyR79vhf4GvX+jnL1qzxz/fs8WNH1qjhWyqbNvl1edcccnL8vXwNGviWzd69PlgORq1aPsySknwLpm5d6NHDh0nLlv4U+NFH+1o2bfL1V6kCLVr4Uz+1a/t95LUeRA6SAqqg00/3//gef7zEgPr3v/0fnTVqVHB9Ep927vQ9azZvhgUL4IMPfHgsWuRbJaU54gjfEtm2zbdgmjb1LaBatQqf0unVy++3Rg0fcnXr+m2aNPGnn5KS/GmbatV8y6hqVR9kVar45x06+M5CdetGHizHHXfQh0WkNAqogpKT/bWo3/7Wn9Jo0+aATfJ68i1eDCecEECNEttWrfJTuMyd6wNo1y5/TrjgKMMtWvjQGToUOnf2ofD99z5s6tTxn2nWDJo3h6OO8mEjkoAUUEWNHg333AMPPeRbUkV06OAfv/pKASX44Hn7bZgxAyZP9kED/vpJhw7+cfRo6NbNX8zv1MlfYxGRUimgijr8cLj4YnjuOfjDH/w5/QKOPtpf3124MKD6pOJt3gyffeZPl2Vnw7p18Mc/Fr4fJDkZzjrLd/UcOdI/5p16E5GDooAK56ab/Ajn//yn79FXQFISHH+8AqpS27ULPv3Ut4pmzPCn63JzC29z3HH+FF1qqh8Vv0cP39FBRMpN1ALKzCYAQ4BNzrnjw7zfF5gMrAqtet059/to1VMm7dv7LuePPgo333zANYBOneCNN/y1KHVSqgR27ID//AemTPE95L76yp+6S0ry4XPXXdCnjw+uqlX9daLOndVCEomyaLagngUeBZ4vYZuPnHNDoljDwbvpJh9Sr7ziT9kU0LGjv5/3++/1R3NccM6HzrRpMH++78GWm+u7aa9Zs39Yl6Qkf1f8qFHQvz/07u27TItIIKIWUM65WWbWIlr7j7oBA3xL6oEH4NJLCzWVOnb0jwsXKqBiUlaW777944/w5pvw7ru+ezf4/2D16vlWUErK/hAaMgROPNH3rhORmBD0NaiTzGwBsB642TkXO2M0mPlW1JVXwsyZ/i/qkLyAWrDAN7IkBsycCRMn+mtGq1YVfq9VK9/h5Yor/PQqIhIXzOWNLxWNnfsW1FvFXIOqDeQ657LN7HTg78651sXsZxQwCqBZs2YnrFmzJmo1F7J7t79n5bjj/L0tBTRr5v/4fuGFiilFiti4EVas8B1ZXnll/0gL/frBSSf5rpbNmvlWcOuw/1uJSIwws3nOuQNGnQ2sBeWc+7nA86lm9g8za+CcO+DWeufceGA8QEZGRvQStajUVLj1Vt+S+ugjf6d+SMeO+88aSZR9/70fU27NGt9h4aeffECB7/o9bJjvuRKmQ4uIxK/AAsrMjgA2OuecmXXDT574Y1D1FGv0aLjvPn+KqMCEhp06+Usbu3frd2K5ys319xa98II/wLNn77/5Ffx9aaec4u+SrlYNLrnED80jIpVONLuZvwz0BRqYWSbwWyAFwDn3BDAcuNbMcoCdwAgXzfONB6tGDbjlFr98+qk/fYRvQeXk+CHWOnUKuMZ4lncQlyyBF1+EWbN85waAY46Bvn19Z4YePfyEXF27+mASkUovqtegoiEjI8PNrejhxLdv99eiMjLgnXcAPxZfu3Z+wIlLL63YcuJaTo6/8fWjj3w3yDff3D8nUOPGvtdJkyY+iIYNC7ZWEakQMXcNKq7UrOmvb9x+e34rqnVr/4e8RpSIQG6uvwn2vff8VCZr1+5/r0cPuOoq32vyrLP82HUiIiigInf99fDgg3DbbfDhhyQnG+3bK6DyzZzpR2NYsQKGD/c3vtav7zs2PPKIvzepalXfgnr4YfjFL3wHh+rVNRyHiISlgIpUzZp+Kvhrr/WtgDPPpGNHP5B1Qtm1y4dQq1bwxBO+l0i9er4zCfhZVKdNK/yZ00/350HPPdd3gKhTp8LLFpH4o2tQZbF3rx8pNjkZFizgkceTGTPGT4paqe7/XLrUt37q1YMnn/TTSGzc6G/8evVV/4Xr1vWjfOfJyPBdwAcN8hfoGjeGZcvgm2/8KTyNWycixdA1qPKQkgL33utbAs89R5cuVwLwxRdxHFDr1vkQevllPxZdx45+yuDkZH86DvyNrkceCX//u+/V+Lvf+Vbk+ef7Eb3fesuP0pB3/Sg93T8efjj07BnI1xKR+KcWVFk5ByefDOvWkf3FUmofUYNx4+Duu4MrqcxeeAE+/th3Wlixwq/r0cP3plu0yI9L16KF723XsqXfvkoVWL7ch1bbtoGWLyKVi1pQ5cUM/vIX6N2btKceok2bO/jii6CLitDatfC3v/lOCgW98w4MHgz79vng6t7dd2go6phjKqZOERH86A1SVr16wdlnwz33MLDNmtgPKOf8NPbHHOPDqXdv33L64gvfBXzwYL9dUpJ/L1w4iYhUMAXUwXroIQBuWD6Gdet8L+qY9OijvrPD//2fv89o1Sr48EPfC69zZ3XxFpGYpYA6WM2awbhxHPP1FM5kCl9+GXRBRUybBh06wA03QJcu8Pjj8NJL/tqSiEgcUEAdirFj2Xfc8TzCDSz8dHvQ1XjO+etMZ57pOzSMG+evMV1zje+ZJyISJxRQhyIlhaTxj9OctRzz8h+CrsZ3cnjgAT8sU79+8Nln8Nvf6pqSiMSliAPKzKqbWZtoFhOXevbk/RaXM3TJX/y0HDt3BlPHU0/560o33+ynq582TePaiUhci+icj5mdCdwPVAVamlk68Hvn3NAo1hY3FlxyP8f9YSpH/OY3/qbVGjXgxht9i2baNH+T64ABfsDUM87wPeeOPdafcnv5Zdi0yU/jsXw5tGkD99/v3//rX0tv/TgHr70Go0b57uF/+5u/eVYjN4hInIvoRl0zmwecAnzgnOscWrfQOdcxyvUdIPAbdcOYPh2uHLSOz868hybTJvjx6PKGAUpJ8UMkFdW4sd9m164D38v7zODBMHasn7Rv+HA/okPB60g5OX6bmTP9zbOzZ/sBWEVE4sih3qib45zbauqSHFbnzpDJUbzS5wl+PflxP4L3M8/4+6V69PCjeTdq5Df8xz986+ipp/y9VFdc4Sfsy82Fbdv84KujR/sed6NH7x949e67fQtp8mQ/q+y0ab5n3syZvqU1ZoyuNYlIpRJpC+ppYCZwO3AuMAZIcc5dE93yDhSLLSjwvc579vS5Um5mzID16/3OX3/dh9rhh/s55pcs8QO6XnedHxtPfzyISJw61BbUDcCdwG7gJeBd4I/lV17869rVn2ErV/3773/ety9ccom/2bZBA3/t6txz/elAEZFKKKKAcs7twAfUndEtJ3517+4bOVlZ0LBhlH5I165+9HF1gBCRBBDRbzoze8/M6hR4XdfM3o1aVXHoxBP94+efR/kHKZxEJEFE+tuugXNuS94L59xmoFFUKopTJ5zgx1qNekCJiCSISAMq18ya5b0ws+ZAfE0kFWU1a/qh7z77LOhKREQqh0g7SdwJ/M/MPgy97g2Mik5J8at7d9+LLzdXZ+JERA5VRL9GnXPTgC7AROBV4ATnnK5BFXHiifDzz/62JhEROTRl+Ts/FfgJ2Aq0M7Pe0SkpfnXv7h91HUpE5NBFOhbfn4Hzga+B3NBqB8yKUl1x6dhj/fisn30Gl18edDUiIvEt0mtQZwFtnHO7o1hL3KtSxZ/mUwtKROTQRXqKbyWgIQsicOKJ8NVXkJ0ddCUiIvEt0hbUDmC+mc3ED3cEgHNuTFSqimPdu/tefPPmQZ8+QVcjIhK/Ig2oKaFFStGtm3/87DMFlIjIoYh0LL7nol1IZdGgARxzjK5DiYgcqkh78bUG7gXaAdXy1jvnWkWprrjWvbufpsk5zYIhInKwIu0k8QzwOJAD9AOeB/4VraLi3YknwoYNkJkZdCUiIvEr0oCq7pybiZ/gcI1zbhx+CngJI++GXY3LJyJy8CINqF1mVgVYZmbXm9nZaDTzYnXs6Ce91XUoEZGDF2lAjQVq4Kd6PwG4BBgZpZriXtWqfvqNTz4JuhIRkfgVaS++OaGn2YAG8YlAz57w4IOwYwfUqBF0NSIi8afEFpSZPRR6fNPMphRdKqTCONWnD+zdq+tQIiIHq7QWVF5PvfujXUhl07OnH5vvww/hFHUnEREpsxIDyjk3L/S0HjBVg8VGrnZt6NzZB5SIiJRdpJ0khgJLzexfZnaGmUU6RFJC69PHn+LbtSvoSkRE4k+kM+peDhwD/Bu4EFhhZk9Fs7DKoE8f2L0bZs8OuhIRkfgT8Yy6zrm9wDvAK8A8YFi0iqosevXyQx3pNJ+ISNlFFFBmNtjMngWWA8OBp4DGUayrUqhbFzp0UECJiByMSK8lXYZvOY1WR4my6dMHnnoK9uzxN/CKiEhkIr0GNQL4EugFYGbVzaxWNAurLPr0gZ07/QSGIiISuUhP8V0NTAL+GVrVFHgjSjVVKr17+0ed5hMRKZtIO0n8EugB/AzgnFuGBouNSMOG0K6dAkpEpKwiDajdzrk9eS9C90G56JRU+fTpA//7H+TkBF2JiEj8iDSgPjSzO4DqZjYAfz/Um9Erq3Lp0weys3UdSkSkLCINqNuBLOArYDQwFfi/aBVV2eSNxTdjRrB1iIjEk0h78eXiO0Vc55wb7px70jmnU3wRatgQ0tMVUCIiZVHadBtmZuPM7AfgW2CJmWWZ2d0VU17lMWAAfPwxbN8edCUiIvGhtBbUWHzvva7OufrOuXrAiUAPM7sx2sVVJv37+/mhPvoo6EpEROJDaQF1KXCBc25V3grn3Erg4tB7EqFevSA1Fd57L+hKRETiQ2kBleKc+6HoSudcFpASnZIqp+rVoUcPXYcSEYlUaQG15yDfkzAGDICFC2HjxqArERGJfaUFVCcz+znMsg3oUBEFVib9+/vHmTODrUNEJB6UGFDOuSTnXO0wSy3nnE7xlVHnzlCvnq5DiYhEIuIJC8vKzCaY2SYzW1TM+2ZmD5vZcjNbaGZdolVLrEhK8jftzpgBuotMRKRkUQso4FlgcAnvnwa0Di2jgMejWEvMGDAAMjNhyZKgKxERiW1RCyjn3CzgpxI2GQY877zPgDpmVuln6c27DqXefCIiJYtmC6o0TYB1BV5nhtZVaq1a+UXXoUREShZkQFmYdWGvzJjZKDOba2Zzs7KyolxW9PXvD++/r+k3RERKEmRAZQJHFXjdFFgfbkPn3HjnXIZzLqNhw4YVUlw0DRgA27bB558HXYmISOwKMqCmAJeGevN1B7Y65zYEWE+F6d8fkpPhTc2oJSJSrGh2M38Z+BRoY2aZZnalmV1jZteENpkKrASWA08C10WrllhTpw707QuTJwddiYhI7EqO1o6dcxeU8r4Dfhmtnx/rhg6FMWNg6VI49tigqxERiT1BnuJLaEOH+ked5hMRCU8BFZDmzaFTJ53mExEpjgIqQMOG+Vl2fzhgQhMREVFABWjoUMjNhbffDroSEZHYo4AKUJcu0KQJTJkSdCUiIrFHARUgM9+Kevdd2LUr6GpERGKLAipgQ4fC9u3w3/8GXYmISGxRQAWsXz9IS1NvPhGRohRQAUtNhcGD/f1QublBVyMiEjsUUDFg2DDYsAHmzg26EhGR2KGAigFnnOEHj33ttaArERGJHQqoGFC3rp+C49VXwYWdEUtEJPEooGLEeefB6tUwb17QlYiIxAYFVIw46yxISfGtKBERUUDFDJ3mExEpTAEVQ847D9asUW8+ERFQQMWUYcN0mk9EJI8CKobUrQsDB8K//63TfCIiCqgYk3eab86coCsREQmWAirG6DSfiIingIoxderoNJ+ICCigYtIvfgFr18Ls2UFXIiISHAVUDBo6FKpW1Wk+EUlsCqgYVPA0n6bgEJFEpYCKUeefD+vWwSefBF2JiEgwFFAx6qyzoEYN+Ne/gq5ERCQYCqgYlZYG55zjr0Pt2hV0NSIiFU8BFcMuuQS2bIG33w66EhGRiqeAimGnngpNmsDTTwddiYhIxVNAxbCkJLjiCpg2zd8XJSKSSBRQMe7KK/2jWlEikmgUUDGueXN/T9SECbBvX9DViIhUHAVUHBg1CjIz/ak+EZFEoYCKA2eeCYcfDk8+GXQlIiIVRwEVB1JS4LLL4K23YP36oKsREakYCqg4cdVV/hrUM88EXYmISMVQQMWJY46BU06Bp57SALIikhgUUHHk6qth9WqYMSPoSkREok8BFUfOPhvq11dnCRFJDAqoOJKaCiNHwhtvwMaNQVcjIhJdCqg4c9VVkJMDzz0XdCUiItGlgIozxx0HvXrB+PEaWUJEKjcFVBwaMwZWrIDXXw+6EhGR6FFAxaGzz4Zjj4V77wXngq5GRCQ6FFBxKCkJbrsNvvwSpk8PuhoRkehQQMWpiy+Gpk3hT38KuhIRkehQQMWpqlXh17+GWbPgk0+CrkZEpPwpoOLY1Vf7G3fvvTfoSkREyp8CKo7VrAm/+pUf5XzhwqCrEREpXwqoOHf99ZCWBvfdF3QlIiLlSwEV5+rWhWuugYkT/b1RIiKVhQKqErjpJkhOhj//OehKRETKjwKqEmjc2HeYeOYZWL486GpERMqHAqqSuPNOPzX8uHFBVyIiUj4UUJVE48Z+jL6XXoKvvgq6GhGRQ6eAqkRuvRVq1YK77gq6EhGRQ6eAqkTq1YNbboHJk+Gzz4KuRkTk0CigKplf/QoOP9z37NNI5yISzxRQlUytWn4A2U8/hVdeCboaEZGDF9WAMrPBZrbEzJab2e1h3u9rZlvNbH5ouTua9SSKkSOhc2d/TWrHjqCrERE5OFELKDNLAh4DTgPaAReYWbswm37knEsPLb+PVj2JJCkJHnoIMjPh/vuDrkZE5OBEswXVDVjunFvpnNsDvAIMi+LPkwJ694bhw/3oEpmZQVcjIlJ20QyoJsC6Aq8zQ+uKOsnMFpjZO2bWPtyOzGyUmc01s7lZWVnRqLVS+stfYN8+P/uuiEi8iWZAWZh1RfuVfQE0d851Ah4B3gi3I+fceOdchnMuo2HDhuVbZSXWsqXvdv7SS/D++0FXIyJSNtEMqEzgqAKvmwLrC27gnPvZOZcdej4VSDGzBlGsKeHccYcPquuugz17gq5GRCRy0QyoOUBrM2tpZlWBEcCUghuY2RFmZqHn3UL1/BjFmhJO9erw6KPw7bfwt78FXY2ISOSiFlDOuRzgeuBdYDHwqnPuazO7xsyuCW02HFhkZguAh4ERzun20vJ2+ulwzjnwhz/AqlVBVyMiEhmLtzzIyMhwc+fODbqMuJOZCccdBz17wtSpYOGuEIqIBMDM5jnnMoqu10gSCaJpU7j3Xpg2DV54IehqRERKp4BKINddBz16wNixsHFj0NWIiJRMAZVAqlSBp56C7Gw/d5SISCxTQCWYtm3h7rvh1VfhjTeCrkZEpHgKqAR0663QsSNcey1oYA4RiVUKqASUkgLPPw8//QRXXql5o0QkNimgElSnTn6svjffhMceC7oaEZEDKaAS2Jgx/ibem2+GhQuDrkZEpDAFVAIzg2eegTp14IILNLmhiMQWBVSCa9TIX4/65hv49a+DrkZEZD8FlDBwoD/N98QT8J//BF2NiIingBIA7rkHTjgBrrpKM/CKSGxQQAkAVavCyy/D7t1w8cWQkxN0RSKS6BRQkq91a3j8cfjwQ7jxxqCrEZFElxx0ARJbLrkEFizwkxu2a+dHmxARCYJaUHKAP/8ZzjgDbrgBZswIuhoRSVQKKDlAUhK89JIfWPa882Dp0qArEpFEpICSsGrX9sMgJSfDmWdqUFkRqXgKKClWy5Z+So61a2HQINiyJeiKRCSRKKCkRD16+Jt3Fy3y16Wys4OuSEQShQJKSjV4MLzyCnz+OQwbBrt2BV2RiCQCBZRE5Jxz4Nln4f33Yfhw2LMn6IpEpLJTQEnELr7Y38j79tv++b59QVckIpWZbtSVMhk92l+HuvlmqFkTnn4aqujPHBGJAgWUlNmvf+1Datw4P138U0/57ugiIuVJv1bkoNx9t2853X23D6sXX4TU1KCrEpHKRAElB8UM7roLatXyA8tu2QKTJvnZeUVEyoOuHsghGTvW9+6bNQu6d9ewSCJSfhRQcshGjvSDyv74I5x4IkyfHnRFIlIZKKCkXPTuDXPmwFFHwWmnwd//7jtQiIgcLAWUlJsWLeCTT/zgsmPHwqhRuqFXRA6eAkrKVVoavP463Hmn737ev79GQheRg6OAknJXpQr88Y9+Tqk5c6BrV1i4MOiqRCTeKKAkai64wPfu27sXTjoJHnpIwyOJSOQUUBJVXbvC3LnQt6+/X+qkk2DBgqCrEpF4oICSqGvcGN56y0/ZsWYNnHAC/OY3sHNn0JWJSCxTQEmFMIPzz4fFi/19U/fdBx06wMyZQVcmIrFKASUVql49PwL6f//rQ6t/f7jsMn+Tr4hIQQooCUS/fr5n3x13+IFm27aFF17Qzb0isp8CSgJTvTrccw/MmwetWsEll8DJJ8OnnwZdmYjEAgWUBK5jRz8CxTPP+E4UJ58MZ5/tg0tEEpcCSmJCUpK/FrV0qZ8I8YMPICMDBg/291KJSOJRQElMSUuD3/7Wt6Tuuw++/BL69IFeveCdd3SNSiSRKKAkJtWuDbfdBqtXwyOP+MA6/XR/D9WkSRqRQiQRKKAkplWvDtdfD8uXw4QJfnr5886DNm3gT3+C774LukIRiRYFlMSFqlXh8sv9jb4TJ/p5p+68E5o18y2rSZNg9+6gqxSR8qSAkriSlAS/+AW8/75vVd1xB3z1lW9VNWni56HSyOkilYMCSuLW0UfDH/7gr1NNmwanngqPPw6dOvkegI89Bt9/H3SVInKwFFAS95KSYNAgf+pv/Xo/3fzevf7aVePG0K0b/P73vkegegGKxA9zcfYvNiMjw82dOzfoMiTGOedP/b35pl9mz/brmjSBM86AIUN8i6tGjaArFREzm+ecyzhgvQJKEsHGjf4+qrfegnff9b0Bq1XzITV4sL/P6vjjfWtMRCqWAkokZPduPzrFW2/51tWqVX597dp+QsWePf3SrZtaWCIVQQElEoZz/ibg//3PLx9/DIsW+feSk6FLl/2B1aMHNGoUbL0ilZECSiRCP/3kR1T/+GMfWrNn77/HqnVrP5htp07+lGCHDnD44X5uKxE5OAookYO0ezd88cX+FtZnn/lrWnnq1/dBlRdYxx8P7dvDYYcFV7NIPFFAiZSjrCx/KvCrrwo/Zmfv3+bww/2QTMce65e85y1bQmpqcLWLxJriAio5iGJE4l3Dhn5W4H799q/LzYW1a31Qff21nzpk6VKYPNkHWkGNGvnhmo46Cpo2Lfx41FFw5JF+eCeRRKaAEiknVapAixZ+GTKk8HubN/uwWrLEj3yRmQnr1sGyZX7Ypq1bC29v5ltgxQVY06Y+xFJSKujLiQRAASVSAerWhRNP9Es427btD6116wo/X7IEZszw24Tbb4MG+5f69Qu/btDAb3PYYX6pUwdq1dL9XhIfohpQZjYY+DuQBDzlnLuvyPsWev90YAdwmXPui2jWJBKLatWC447zS3F+/rlweH33Hfzww/5l3TqYP9+fTty1q/SfVzC0DjvM3weWlgY1a/rHgs/DrUtL8/eJVavmr6lV0cBpUs6iFlBmlgQ8BgwAMoE5ZjbFOfdNgc1OA1qHlhOBx0OPIlJE7dq+d2D79qVvu2OHD62sLNiyxS9bt+5fCr7essUPqrt0KWzf7jt6bN/ur6mVRdWqPqxKWqpX3/88NdUvKSn+s4f6WPB5crJvJZa0KFBjXzRbUN2A5c65lQBm9gowDCgYUMOA553vSviZmdUxs8bOuQ1RrEuk0qtRw8+V1azZwX3eOd8KywusvNAq+Jid7YNw926/7a5dsHPn/udFl61bfff8gtvu3Qt79ux/rGilhVjBMCvLtlWq+OuI0XzMe150gfDri1sOdfszzoCOHaPz3yeaAdUEWFfgdSYHto7CbdMEKBRQZjYKGAXQ7GD/xYlIxMx8a6d6dX8dqyI4B/v2FQ6soo8lvVfwce9ev6/yWnJzy7Z9To7/jHPReSy6Lm/JO46RLiVtH6nGjeMzoMLdW1/0a0eyDc658cB48PdBHXppIhJrzPypuWR13YoZkQRaNP97RfN/hUzgqAKvmwLrD2IbEREJQMHTekGI5mXCOUBrM2tpZlWBEcCUIttMAS41rzuwVdefREQEotiCcs7lmNn1wLv4buYTnHNfm9k1ofefAKbiu5gvx3czvzxa9YiISHyJ6tle59xUfAgVXPdEgecO+GU0axARkfikOwFERCQmKaBERCQmKaBERCQmKaBERCQmKaBERCQmKaBERCQmKaBERCQmKaBERCQmKaBERCQmKaBERCQmKaBERCQmmSvLzFQxwMyygDWHsIsGwA/lVE5loWNSmI5HYToehel4FFYex6O5c65h0ZVxF1CHyszmOucygq4jluiYFKbjUZiOR2E6HoVF83joFJ+IiMQkBZSIiMSkRAyo8UEXEIN0TArT8ShMx6MwHY/ConY8Eu4alIiIxIdEbEGJiEgcUECJiEhMSqiAMrPBZrbEzJab2e1B11MRzGyCmW0ys0UF1tUzs/fMbFnosW6B934TOj5LzGxQMFVHj5kdZWbvm9liM/vazH4VWp+Qx8TMqpnZbDNbEDoevwutT8jjAWBmSWb2pZm9FXqdsMcCwMxWm9lXZjbfzOaG1lXMMXHOJcQCJAErgFZAVWAB0C7ouirge/cGugCLCqz7C3B76PntwJ9Dz9uFjksq0DJ0vJKC/g7lfDwaA11Cz2sBS0PfOyGPCWBAWuh5CvA50D1Rj0foO94EvAS8FXqdsMci9D1XAw2KrKuQY5JILahuwHLn3Ern3B7gFWBYwDVFnXNuFvBTkdXDgOdCz58Dziqw/hXn3G7n3CpgOf64VRrOuQ3OuS9Cz7cBi4EmJOgxcV526GVKaHEk6PEws6bAGcBTBVYn5LEoRYUck0QKqCbAugKvM0PrEtHhzrkN4H9hA41C6xPqGJlZC6AzvtWQsMckdEprPrAJeM85l8jH4yHgViC3wLpEPRZ5HDDdzOaZ2ajQugo5JskH+8E4ZGHWqY99YQlzjMwsDXgNGOuc+9ks3Ff3m4ZZV6mOiXNuH5BuZnWA/5jZ8SVsXmmPh5kNATY55+aZWd9IPhJmXaU4FkX0cM6tN7NGwHtm9m0J25brMUmkFlQmcFSB102B9QHVErSNZtYYIPS4KbQ+IY6RmaXgw+lF59zrodUJfUwAnHNbgA+AwSTm8egBDDWz1fhLAKeY2Qsk5rHI55xbH3rcBPwHf8quQo5JIgXUHKC1mbU0s6rACGBKwDUFZQowMvR8JDC5wPoRZpZqZi2B1sDsAOqLGvNNpaeBxc65Bwq8lZDHxMwahlpOmFl1oD/wLQl4PJxzv3HONXXOtcD/fvivc+5iEvBY5DGzmmZWK+85MBBYREUdk6B7iFRwb5TT8b22VgB3Bl1PBX3nl4ENwF78XzdXAvWBmcCy0GO9AtvfGTo+S4DTgq4/CsejJ/6Uw0Jgfmg5PVGPCdAR+DJ0PBYBd4fWJ+TxKPAd+7K/F1/CHgt8r+cFoeXrvN+bFXVMNNSRiIjEpEQ6xSciInFEASUiIjFJASUiIjFJASUiIjFJASUiIjFJASWVnpntC43EvMDMvjCzk8t5/3cUef1JOe23b4ERtfuWZ91m1sLMLizwOsPMHi6v/YuUBwWUJIKdzrl051wn4DfAveW8/0IB5Zwr1wAM6QuUab9mVtJQZi2A/IByzs11zo05qMpEokQBJYmmNrAZ/KgSZvZXM1sUmu/m/FLWNzazWaHW2CIz62Vm9wHVQ+teDG2XHXrsa2YfmNkkM/vWzF4MjWSBmZ0eWvc/M3s4r6UUTmhQ22uAG0M/p1doBIjXzGxOaOkR2nacmY03s+nA86GW0kehlmPB1uN9QK/Q/m4s0lqrZ2ZvmNlCM/vMzDoW2PeE0HdaaWZjQutrmtnboRbqorzjJXKoEmmwWElc1c2P1l0NPx/UKaH15wDpQCegATDHzGbhWyrh1l8IvOucu8fMkoAazrmPzOx651x6MT+7M9AePx7Zx0AP85O+/RPo7ZxbZWYvl1S8c261mT0BZDvn7gcws5eAB51z/zOzZsC7wHGhj5wA9HTO7TSzGsAA59wuM2uNH1kkAz+Hz83OuSGh/fUt8CN/B3zpnDvLzE4Bng8dD4C2QD/8XFpLzOxx/Nh9651zZ4T2dVhJ30ckUgooSQQ78wLEzE7CtyyOxw979LLzo3lvNLMPga4lrJ8DTDA/2Owbzrn5Efzs2c65zNDPno8/tZYNrHR+vhzwoTEq7KeL1x9oZ/tHYa+dN2YaMMU5tzP0PAV41MzSgX3AsRHsuydwLoBz7r9mVr9A6LztnNsN7DazTcDhwFfA/Wb2Z/zwQB+V8buIhKVTfJJQnHOf4ltFDQk/NQDFrXd+8sfewHfAv8zs0gh+5O4Cz/fh/ygsdm6PMqgCnBS6tpbunGvi/ASMANsLbHcjsBHfGszAzyZdmpKmTDjg+zjnluJbbV8B95rZ3WX4HiLFUkBJQjGztkAS8CMwCzjf/IR9DfHhM7u49WbWHD9f0JP4EdG7hHa7N9SqitS3QKvQtSWASK7ZbMOfVsszHbi+wPdKL+ZzhwEbnHO5wCX47x5ufwXNAi4K7bcv8INz7ufiCjOzI4EdzrkXgPvZf1xEDolO8UkiyLsGBb51MNI5t8/M/gOchB+p2QG3Oue+L2H9SOAWM9uLP02X14IaDyw0sy+ccxeVVkzo2tB1wDQz+4HIpiN4E5hkZsOAG4AxwGNmthD/73gWviNFUf8AXjOz84D32d+6WgjkmNkC4Fn8iOZ5xgHPhPa9g/3TKhSnA/BXM8vFj5p/bQTfR6RUGs1cJABmluacyw716nsMWOacezDoukRiiU7xiQTj6lCr7mv8abh/BluOSOxRC0pERGKSWlAiIhKTFFAiIhKTFFAiIhKTFFAiIhKTFFAiIhKT/h+GcLfuXHXyHQAAAABJRU5ErkJggg==",
      "text/plain": [
       "<Figure size 432x432 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "test_score = np.zeros((params[\"n_estimators\"],), dtype=np.float64)\n",
    "for i, target_timeseries_pred in enumerate(regressor.staged_predict(clustered_data_test)):\n",
    "    test_score[i] = regressor.loss_(target_timeseries_test, target_timeseries_pred)\n",
    "\n",
    "fig = plt.figure(figsize=(6, 6))\n",
    "plt.subplot(1, 1, 1)\n",
    "plt.title(\"Deviance\")\n",
    "plt.plot(\n",
    "    np.arange(params[\"n_estimators\"]) + 1,\n",
    "    regressor.train_score_,\n",
    "    \"b-\",\n",
    "    label=\"Training Set Deviance\",\n",
    ")\n",
    "plt.plot(\n",
    "    np.arange(params[\"n_estimators\"]) + 1, test_score, \"r-\", label=\"Test Set Deviance\"\n",
    ")\n",
    "plt.legend(loc=\"upper right\")\n",
    "plt.xlabel(\"Boosting Iterations\")\n",
    "plt.ylabel(\"Deviance\")\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot feature importance."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAA1gAAAGoCAYAAABbkkSYAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAABGZ0lEQVR4nO3de5hldX3n+/cnIHS0paEBDTZqj0oEQWylgjMnjAmYIQkcJYwOHEXxgmEIz5jMcSQyac20cRAcjIM64xAiM1GCZKKR2EQcYtooQUEp2uYWiChpUIkC0gKN4OHyPX+sVbC72HXZVatq1+X9ep719N5rrd9vfdfqqr3qu3+XlapCkiRJkjR7PzPsACRJkiRpqTDBkiRJkqSOmGBJkiRJUkdMsCRJkiSpIyZYkiRJktQREyxJkiRJ6ogJliRJ0jKQ5Nwk7xl2HMtVkn+b5JxhxzEbSX47yVnDjmOhM8GSJpBka5IHk2zvWZ7VQZ2/0lWM0zjehiR/Ol/Hm0ySNye5YthxSNJkxn32/zDJ/0qycgHE9SdJ/vMA+z/pM7eqTqmq981BbN5rppBkF+DdwNnt+7VJKsnOHdQ90M/GAPX+cpLvjVt9HvCGJM/o+nhLiQmWNLlXVdXKnuWOYQbTxQfxMCzWuCUtW6+qqpXAy4BfoPnDeNrS8G+sebbA7zXHADdX1feHHchsVNVDwBeAE4cdy0LmL780oCSrkpyf5J+SfD/Jf06yU7vt+Um+lORHSe5OcmGS3dttFwDPAS5pvxn93X7fDvW2crXfCn4myZ8muQ9482THn0bsleTUJLckuT/J+9qYr0xyX5I/b79le/ybqyS/157L1iQnjLsOn0xyV5Lbkrx77A+K9hvEryb5r0nuAf43cC7wL9pz/3G739FJvtke+7tJNvTUP/bt3puS3N7GsL5n+05tbN9pz+WaJM9ut+2f5ItJ7knyD0mOG+g/WZKA9o/hLwAHAST550m+luTHSa5N8stj+yb5cpIzknwV+AnwvAE/c5/U8tKWf0GSk4ETgN9tP0Mvabef3vMZ+PdJjm3XH0D/z9wdWjqS/GaSb7eflRvT00ujPfYpbezbkvz3JJnOdfNe09evA1/peX95+++P21j/RVvnW5Pc1F7zy5I8t12f9jzvTHJvkuuSHDTRz8a4/4++Zdttuyb5YHvuP0zTjfRnkzyN5mf/WXlyL54vA0dP9jOw7FWVi4tLnwXYCvxKn/V/CfwR8DTgGcA3gH/bbnsB8K+AXYG9aT5Az5moTuCXge9NdFxgA/Aw8Bs0X4j87GTH7xPrBuBPe94XsBHYDTgQ+CmwCXgesAr4e+BNPbE9AnyoPZ9fAh4AXthu/yTwOeDpwFrgW8BJ7bY3t2XfDuzcxv1m4Ipx8f0y8OL23A4Gfgj8RrttbRvvH7flX9LGe0C7/TTgeuCFQNrte7bX5bvAW9pjvwy4Gzhw2D9TLi4uC38Z9xn8bOBG4H3AGuBHwFHtZ9a/at/v3e77ZeD29rN1Z+ApA37m9vuMLOAF7es/Af7zuO3/BnhWG8/x7Wf0PpPU93gdwBHtZ+PL2s/4jwKXjzv2XwG703w5eBfwaxNcsw14r5n0XgNcDfybnvdjx925Z91vAN8GDmjrfDfwtXbbrwLXtP8faffZZ/z/6wTHnqzsOe3/1er2Gl8CnNlz3b7Xp76XAfcM+3d1IS+2YEmT+8s031T+OMlfJnkmzbdQ/76qHqiqO4H/Cvw/AFX17ar6YlX9tKruorlh/NIsY7iyqv6yqh6juVlNePxp+kBV3VdVNwI3AH9dVbdW1b0031a9dNz+72nP5yvA54Hj0rSYHQ/8x6q6v6q2An8IvLGn3B1V9dGqeqSqHuwXSFV9uaqur6rHquo64CKefL3eW1UPVtW1wLU0NzeAtwHvrqp/qMa1VfUj4P8GtlbV/2qPvRn4C+C1A1wjScvbX7atH1fQtDq8H3gDcGlVXdp+Zn0RGKVJuMb8SVXd2H72PNyuG/Qzd9qq6tNVdUcbz/8GbgEOnWbxE4D/WVWbq+qnwH+kaflZ27PPWVX146q6HfhbYN0A4Xmv2dHuwP2TXjH4tzTJzU1V9QjNz926thXrYZoEaH8g7T7/NEV9Y/qWbVskfxP4f6vqnqq6vz3mVH9T3E+TKGsCC7mvqrQQ/EZV/c3YmySH0nwr+U89PSV+huZbLNIM+vwI8C9pPsx+Btg2yxi+2/P6uZMdf5p+2PP6wT7vf67n/baqeqDn/W0035buBezSvu/dtmaCuPtK8nLgLJruN7vQfHv56XG7/aDn9U+AscHmzwa+06fa5wIvH+sa0toZuGCqeCSptcNnP0D7R+6/SfKqntVPoUk8xvT73BvkM3cgSU4E3kHTGgLN5+Ne0yz+LGDz2Juq2p7kRzSf41vb1RN9/k6H95odbaP5u2AyzwU+nOQPe8MH1lTVl5L8N+C/A89JcjHwzqq6b4o6magssAJ4KnBNz98UAaYadvB04N6pjruc2YIlDea7NF0H9qqq3dtlt6o6sN1+Jk2T/8FVtRvNN569fdZrXH0P0Hy4AU1fb5quhb16y0x1/K7t0fbDHvMc4A6abhAP09wMerf1Dt4df67j3wN8iqZrwrOrahVN3/lp9fGnuRbPn2D9V3quz+7VTFDyW9OsV5L6+S5wwbjPlqdVVe+U1f0+56Zr/P1gfOK1Q91twvfHwL8D9qyq3WlaitJv/z7uoOczvP2s35MdP8fny3K411wH/PwUcX6Xpst/b50/W1VfA6iqj1TVITTdLn+epvviRHXtYIKyd9Mkuwf2HG9VNRO8TFbvATStfJqACZY0gLY5/q+BP0yyW5KfaQfujnU1eDqwnWbQ6hqe+PAb80OaPuhjvgWsaAfgPoWmv/Wuszj+XHhvkl2S/EuaLhGfrqpHgT8Hzkjy9PZG/w5gsml6fwjsm3Zgc+vpNP24H2pbB18/QFwfB96XZL92AO/BSfakGTPw80nemOQp7fILaQZ9S9JM/SnwqiS/2k58sCLNBA37dlT/tcCBSdYlWUEzrqnX+PvH02j+AL4LIMlbaCfj6Nl//Gdur08Bb2mPtytN17Cvt93whmGp32suZcduiXcBj7Hj/+m5wH9MciA8PsHHv2lf/0KSl7d/KzwAPAQ82nPOvfXsYKKy1Qw9+GPgv7Y9cEiyJsmv9tS7Z5Lx3QF/iaabpyZggiUN7kSaLgZ/T9Pk/xlgn3bbe2kGf95L04f8s+PKngm8ux3T9c62L/qpNB/g36f54Bv/zIlBjt+1H7THuAO4EDilqm5ut72dJt5bacYpfAr4n5PU9SWaweI/SHJ3u+5U4A+S3A/8Ps2NdLo+1O7/18B9wPnAz7Z9yI+k6UN+R3sOH2CSxFWSplJV36WZavv3aP44/i7Nl2id/C1VVd8C/gD4G5qxVOOf5XQ+8KKxMcFV9fc045GupPlD+MXAV3v27/eZ23u8TcB7aMYN/RNNK80g43m7tBzuNZcA+4/NxFdVPwHOAL7a/p/+86q6uK3jz9LMHHwDzbhraMZg/zHNdbqNZoKVD7bbdvjZ6HPsycq+i2ZijavaY/4NzYQetP8HFwG3tnU/q03+jwI+MY1rt2ylajat2ZKWqjTTD/9pVXX17awkSTtYTveaNFOqv6iq/v2wY5mpJG+n6Wr5u8OOZSFzkgtJkiRpjlXVecOOYbaq6qPDjmExsIugJEmSJHXELoKSJEmS1BFbsCRJkiSpI47Bmkd77bVXrV27dthhSNKCd80119xdVeOfCad54L1KkqZnonuVCdY8Wrt2LaOjo8MOQ5IWvCS3DTuG5cp7lSRNz0T3KrsISpIkSVJHTLAkSZIkqSMmWJIkSZLUERMsSZIkSeqICZYkSZIkdcQES5IkSZI6YoIlSZIkSR0xwZIkSZKkjphgSZIkSVJHTLAkSZIkqSMmWJIkSZLUERMsSZIkSerIzsMOQJIkSRqzevVqtm3bNuwwZqX+027kvfcNO4yh2WOPPbjnnnuGHcbQmGBJkiRpwdi2bRtVNewwZmfDqsV/DrOQZNghDJVdBCVJkiSpIyZYkiRJktQRuwjOo+u/fy9rT//8sMOQpDm39ayjhx2CJElDYQuWJEmSJHXEBEuSJEmSOmKCJUmSNIXlPiuatFTNxe+2CZYkSZIkdcQES5IkSZI6YoIlSZIkSR0xwZIkLXpJtndUzwlJrmuXryV5yQzq2JDknTMot3uSUwctN0Wd65JcmeTG9pyO77J+SdKTmWBJkvSEfwR+qaoOBt4HnDePx94dGCjBSmOye/lPgBOr6kDg14Bzkuw+4wglSVMywZIkLRlJVibZlGRzkuuTHNOz7T1Jbk7yxSQX9WtlqqqvVdW29u1VwL5THO/EtmXo2iQX9Nn+5SQj7eu9kmxtXx+Y5BtJtrTl9wPOAp7frju73e+0JFe3+7y3Xbc2yU1JPgZsBp49UXxV9a2quqV9fQdwJ7B3nzhPTjKaZPSuu+6a7JSXtSQu87BoaRj2z9Ewf952npNaJUkajoeAY6vqviR7AVcl2QgcArwGeCnNvW8zcA1AklMAqurccXWdBHxhogMlORBYD/xiVd2dZPUAcZ4CfLiqLkyyC7ATcDpwUFWta+s/EtgPOBQIsDHJK4DbgRcCb6mqU9t9LwXe1iZRE8V7KLAL8J3x26rqPNrWupGRkRrgPJaVKi/NfDDJWhoWy+/LXPy8mWBJkpaSAO9vE5HHgDXAM4HDgM9V1YMASS4ZK9AnsSLJ4TQJ1mGTHOsI4DNVdXdbzz0DxHklsD7JvsBnq+qWPjf5I9vlm+37lTQJ1+3AbVV1Vc85HDXZwZLsA1wAvKmqHhsgTknSgOwiKElaSk6g6QJ3SNsS9ENgBU3iNS1JDgY+DhxTVT+abFdgqq9oH+GJe+2KsZVV9Sng1cCDwGVJjpig/jOral27vKCqzm+3PTBB7C9P08VwS5JXt+t2Az4PvLs3KZMkzQ0TLEnSUrIKuLOqHm5boZ7brr8CeFWSFUlWAkf3K5zkOcBngTdW1bemONYm4Lgke7Zl+3UR3ErTPRHgtT3HeR5wa1V9BNgIHAzcDzy9p+xlwFvbeEmyJskzJguoqr7ek5BtbLsfXgx8sqo+PcX5SJI6YIIlSVpKLgRGkozStGbdDFBVV9MkMtfSJFCjwL3QjMEaG4cF/D6wJ/CxthVodKIDVdWNwBnAV5JcC3yoz24fBH4rydeAvXrWHw/ckGQLsD9NAvQj4KtJbkhydlX9NfAp4Mok1wOfYccE7HFJLk3yrD6bjgNeAby5p2Vr3UTnJEmavSyWAWhLwa777Ff7vOmcYYchSXNu61l9G4imLck1VTXSUThjda6squ1JngpcDpxcVZu7PMZSMDIyUqOjE+aVy1aSRTNof7FbEtd6wyrYcO+woxiaxfR/OJtYJ7pXOcmFJGm5OC/Ji2jGQn3C5EqSNBemTLCSbK+qlbM9UJITgHe1b7cDv1VV1w5YxwZge1V9cMByuwOvr6qPDVJuijrXAf8D2A14FDijqv53V/VLkrpVVa+fSbl2jNWmPpteOcUkGFpCFsu38ZIGMxe/2/PZgvWPwC9V1bYkv07zvI2Xz9OxdwdOBaadYKWZLzeTTGf7E+DEdmrdZwHXJLmsqn4822AlSQtHm0StG3YckqTFYdqTXCRZmWRTks1Jrk9yTM+29yS5OckXk1yU5J3jy1fV16pqW/v2KmDfKY53Ypon11+b5II+27+cZKR9vVeSre3rA5N8ox3Ie12S/YCzgOe3685u9zstydXtPu9t161NclOSj9E8hPLZE8VXVd+qqlva13cAd9JMDSxJkiRpmRqkBesh4Niqui/JXsBVSTbSTD/7GuClbX2bgWugmZkJ+j7E8STgCxMdKMmBwHrgF6vq7gmmvp3IKcCHq+rCdnranYDTgYPaZ6KQ5EiahzUeSvOckY3tQylvB14IvKWqTm33vRR4W5tETRTvocAuwHf6bDsZOBlgp93MvyRJkqSlbJAEK8D720TkMWAN8Eyap9x/rqoeBEhyyViBPokV7XNJTmrLTeQI4DNVdXdbzz0DxHklsD7JvsBn2y584/c5sl2+2b5fSZNw3Q7c1vsgxqo6arKDJdkHuAB4U7/uhFV1Hk13SHbdZz87cEuSJElL2CDPwTqBpgvcIW1L0A9pZmJ6UvYykSQHAx8HjpliYHCAqZKRR3gi/hVjK6vqU8CrgQeBy5IcMUH9Z/Y8jPEFVXV+u+2BCWJ/ec8zRF7drtsN+Dzw7t6kTJIkSdLyNEiCtQq4s6oebluhntuuvwJ4VZIV7dPm+z78JMlzaB7u+Maq+tYUx9oEHNfO3MQEXQS30nRPBHhtz3GeB9xaVR+heajkwcD97PhwxsuAt7bxkmRNkmdMFlBVfb0nIdvYdj+8mObhkJ+e4nwkSZI0TUkW9bIUzmE2yx577DHkn6DhGqSL4IXAJWmear8FuBmgqq5OMxbrWuA2YBS4F540Buv3gT2Bj7U/eI9M9BDJqroxyRnAV5I8StOV783jdvsg8OdJ3gh8qWf98cAbkjwM/AD4g6q6J8lXk9wAfKGqTktyAHBlG8t24A00063vIBOPwToOeAWwZ5Kx2N5cVVv6nZMkSZKmtlSmxK8Nw45Aw5IufoiTrKyq7UmeClwOnOwDHJ9s1332q33edM6ww5CkObf1rL6dGaYtyTUTfQmnuTUyMlKjo6PDDkOSFryJ7lVdPQfrvCQvohkL9QmTK0mSJEnLUScJVlW9fibl0oyx2tRn0yunmARDkiRJkhacrlqwZqRNotYNMwZJkiRJ6sogswhKkiRJkiZhgiVJkiRJHTHBkiRJkqSOmGBJkiRJUkdMsCRJkiSpIyZYkiRJktQREyxJkiRJ6ogJliRJkiR1xARLkiRJkjpigiVJkiRJHTHBkiRJkqSO7DzsAJaTF69ZxehZRw87DEmSJElzxBYsSZIkSeqICZYkSZIkdcQES5IkSZI6YoIlSZIkSR0xwZIkSZKkjphgSZIkSVJHTLAkSZIkqSMmWJIkSZLUERMsSZIkSerIzsMOYDm5/vv3svb0zw87DGlObT3r6GGHIEkzsnr1arZt2zZQmfpPu5H33jflfnvssQf33HPPTEOTtIjYgiVJkgRs27aNqhpoAaa136CJm6TFywRLkiRJkjpigiVJkiRJHTHBkiRJkqSOmGBJkiRJUkdMsCRJkiSpIyZYkiRpWUmyrI8vaW6ZYEmSJElSR0ywJEmSJKkjJliSJEmS1BETLEnSkpJke0f17J/kyiQ/TfLOGdaxYSZlk+ye5NSZHHNcPf8uybeTVJK9ZlufJGlqJliSJPV3D/DbwAeHcOzdgYESrDTG39e/CvwKcFtHcUmSpmCCJUlakpKsTLIpyeYk1yc5pmfbe5LcnOSLSS7q18pUVXdW1dXAw9M83olJrktybZIL+mz/cpKR9vVeSba2rw9M8o0kW9ry+wFnAc9v153d7ndakqvbfd7brlub5KYkHwM2A88edw7frKqt04j95CSjSUbvuuuu6ZzuopfkSct8HVPS0rbzsAOQJGmOPAQcW1X3td3jrkqyETgEeA3wUpr74GbgGoAkpwBU1bmDHCjJgcB64Ber6u4kqwcofgrw4aq6MMkuwE7A6cBBVbWurf9IYD/gUCDAxiSvAG4HXgi8papObfe9FHhbVd0x3QCq6jzgPICRkZEaIPZFq+rJpznXyc/YMU2ypKXNBEuStFQFeH+biDwGrAGeCRwGfK6qHgRIcslYgUETqx5HAJ+pqrvbeu4ZoOyVwPok+wKfrapb+vwBfmS7fLN9v5Im4boduK2qrhrbsaqOmtkpSJK6YIIlSVqqTgD2Bg6pqofbLnkraBKvrgWYquXnEZ7omr9ibGVVfSrJ14GjgcuSvA24tU/9Z1bVH+2wMlkLPDCLuCVJHXMMliRpqVoF3NkmV4cDz23XXwG8KsmKJCtpEpvZ2gQcl2RPgAm6CG6l6Z4I8NqxlUmeB9xaVR8BNgIHA/cDT+8pexnw1jZekqxJ8owO4pYkdWzgBMvpb3eox+lvJWnhuhAYSTJK05p1M0A7ccVG4Frgs8AocC80Y7DGxmEl+bkk3wPeAbw7yfeS7NbvQFV1I3AG8JUk1wIf6rPbB4HfSvI1oPeecTxwQ5ItwP7AJ6vqR8BXk9yQ5Oyq+mvgU8CVSa4HPsOOCdjjklya5Fnt699uz2Ff4LokH5/6skmSZiP9BnlOWiDZXlUrZ33g5pu35wK/AWyrqoGnwU2yAdg+aNm2S8VfVdVBA5QJzfV6rGfdS4FtwJeBkbG+9xPZdZ/9ap83nTNIqNKis/WsLhoDtNwluaaqRuaw/pVVtT3JU4HLgZOravNcHW8xGRkZqdHR0WGHMaeSTDjJxaB/F7FhFWy4d6Bjzug4khacie5VM+4i6PS305v+tnfq20d/MvUHsCRpXpzXthhtBv7C5EqS1JXZTHLh9LfT0Dv17a777OfXVZK0AFTV62dSrh1jtanPple23fq0CAy79WjYx5c0t2aTYDn9rSRpWWmTqHXDjkOStHDNJsFy+ltJkiRJ6jGbadqd/laSJEmSeswmwXL6W6e/lSRJktRj4Gnap1Wp09/25TTtWg6cpl1dmOtp2jWx5TBN+0Tma5p2SUtD59O0T8HpbyVJ0qKTZKBlumX22GOPIZ+ZpPkym0kuJuT0t5IkabGZaQtTbeg2DkmL25wkWDPl9LeSJEmSFrO56iIoSZIkScuOCZYkSZIkdcQES5IkSZI6YoIlSZIkSR0xwZIkSZKkjphgSZIkSVJHTLAkSZIkqSMmWJIkSZLUERMsSZIkSeqICZYkSZIkdWTnYQewnLx4zSpGzzp62GFIkiRJmiO2YEmSJElSR0ywJEmSJKkjJliSJEmS1BETLEmSJEnqiAmWJEmSJHXEBEuSJEmSOmKCJUmSJEkdMcGSJEmSpI6YYEmSJElSR3YedgDLyfXfv5e1p39+2GFoCdh61tHDDkGSJEl92IIlSZIkSR0xwZIkSZKkjphgSZIkSVJHTLAkSZLGWb16NUlgwyqSTLisXr162KFKWmBMsCRJksbZtm0bVQVAVU24bNu2bciRSlpoTLAkSZIkqSMmWJIkSZLUERMsSZIkSeqICZYkSZIkdcQES5IkLWtJFnR9khYXEyxJkiRJ6ogJliRJkiR1xARLkiRJkjpigiVJkiRJHTHBkiQtekm2d1TP/kmuTPLTJO+cYR0bZlI2ye5JTp3JMaeo9/8k+XGSv+q6bknSk5lgSZL0hHuA3wY+OIRj7w4MlGClMdW9/GzgjTMNSpI0GBMsSdKSkWRlkk1JNie5PskxPdvek+TmJF9MclG/VqaqurOqrgYenubxTkxyXZJrk1zQZ/uXk4y0r/dKsrV9fWCSbyTZ0pbfDzgLeH677ux2v9OSXN3u89523dokNyX5GLAZePZkMVbVJuD+Kc7j5CSjSUbvuuuu6Zz6kpNkh2W+ykpaegZOsDrshnFCe8O4LsnXkrxkBnUMtRtGkn+X5NtJKsles61PkjRrDwHHVtXLgMOBP2xbeUaA1wAvBf41MDJWIMkpSU4Z9EBJDgTWA0dU1UuA3xmg+CnAh6tqXRvL94DTge9U1bqqOi3JkcB+wKHAOuCQJK9oy78Q+GRVvbSqbktyaZJnDXoOY6rqvKoaqaqRvffee6bVLGpVtcMyX2UlLT07D/HY/wj8UlVtS/LrwHnAy+fp2LvTdMP42HQLpPlKKlX1WM/qrwJ/BXy5y+AkSTMW4P1tIvIYsAZ4JnAY8LmqehAgySVjBarq3Bke6wjgM1V1d1vPPQOUvRJYn2Rf4LNVdUuflo8j2+Wb7fuVNAnX7cBtVXXV2I5VddTMTkGS1LUZdxHsoBvG16pqW/v2KmDfKY634LphVNU3q2rrdK+ZJGnOnQDsDRzStg79EFhBk3h1LcBUzRWP8MS9dsXYyqr6FPBq4EHgsiRHTFD/mW2L1rqqekFVnd9ue6BvQMnL23vbliSvHuRkJEndmM0YrC67YZwEfGGiAy3mbhi9/dof/cm9A4QtSZqBVcCdVfVwksOB57brrwBelWRFkpXA0R0caxNwXJI9AZKs7rPPVuCQ9vVrx1YmeR5wa1V9BNgIHEwzTurpPWUvA97axkuSNUmeMVlAVfX1noRs48xOS5I0G7PpIthJN4z2BnhSW24ii7YbRlWdR9P9kV332c+O2ZI0ty4ELkkyCmwBbgaoqquTbASuBW4DRoF7ofnyr93n3CQ/127bDXgsyb8HXlRV940/UFXdmOQM4CtJHqW5h7x53G4fBP48yRuBL/WsPx54Q5KHgR8Af1BV9yT5apIbgC+0XwAeAFzZ3re2A28AHh0fS5JLgbdV1R19tv0dsD+wMsn3gJOq6rLJLqIkaeZmk2D1dsN4uO2SN1A3jCQHAx8Hfr2qfjTZrsyiG0aSr9N8W3lZkrcBt/ap/8yq+qNx8a1lgm4YkqSFo6pWtv/eDfyLCXb7YFVtSPJU4HLgD9syj3/5V1U/YIou6+OO+wngE+PWbeh5fTNN69SYd7frzwTO7FPf68e9/zDw4T6HPmjcfhN++VdV/3LCE5AkdW42XQRn1Q0jyXOAzwJvrKpvTXGsBdcNQ5K06JyXZAvNmNq/qKrNQ45HkrQEzaYFa1bdMIDfB/YEPtZ2fXikqkboY6F2w0jy28DvAj8HXJfk0qp62/QunyRpPo1vHZqu9su9TX02vXKK3hdaJLqeWt2p2qXlLXPxIZBkZVVt7+mGcbLfFDZjsPZ50znDDkNLwNazuhifLy1cSa6Z6Es3za2RkZEaHR0ddhhDl6RJlDasgg0TT1L1+H6Slp2J7lVz9Rys85K8iGYs1CdMriRJkiQtB3OSYNkNQ5IkSdJyNFctWDPSJlHrhh2HJEmSJM3EbGYRlCRJkiT1MMGSJEnqo51ZmCQTLnvssceQo5S00CyoLoKSJEkLQe/MgE88OlqSpmYLliRJkiR1xARLkiRJkjpigiVJkiRJHTHBkiRJkqSOmGBJkiRJUkdMsCRJkiSpIyZYkiRJktQREyxJkiRJ6ogPGp5HL16zitGzjh52GJIkSZLmiC1YkiRJktQREyxJkiRJ6ogJliRJkiR1xARLkiRJkjpigiVJkiRJHTHBkiRJkqSOmGBJkiRJUkdMsCRJkiSpIyZYkiRJktSRnYcdwHJy/ffvZe3pnx92GFqEtp519LBDkCRJ0jTYgiVJkiRJHTHBkiRJkqSOmGBJkiRJUkdMsCRJkiSpIyZYkiRJktQREyxJkiRJ6ogJliRJWp42rGL16tXDjkLSEmOCJUmSlq1t27YNOwRJS4wJliRJkiR1xARLkiRJkjpigiVJkiRJHTHBkiRJkqSOmGBJkqRlJcmwQ5C0hJlgSZIkSVJHTLAkSZIkqSMmWJIkSZLUERMsSZK05F100UUcdNBB7LTTTo+/H7NixQqSsGLFCt7+9rcPK0RJS8SUCVaS7V0eMMkvJHk0yWtnUHZDknfOoNzuSU4dtNwUda5LcmWSG5Ncl+T4LuuXJE1fV/eqJMe0n+lbkowmOWwGdSyYe1Vb7/9J8uMkf9V13YvFRRddxPr16/noRz/KQw89BMD69esf3/7+97+fBx54gPe///2ce+65JlmSZmVeW7CS7AR8ALhsPo8L7A4MdNNKY7Lr8xPgxKo6EPg14Jwku884QknSQrAJeElVrQPeCnx8Ho+9O93fqwDOBt4406CWgjPOOIPzzz+fww8/nKc85SkAnH/++Y9vf8c73sFTn/pU3vGOd/CBD3yAP/7jPx5WqJKWgGknWElWJtmUZHOS65Mc07PtPUluTvLFJBdN8s3d24G/AO6cxvFObL9FvDbJBX22fznJSPt6ryRb29cHJvlG++3jdUn2A84Cnt+uO7vd77QkV7f7vLddtzbJTUk+BmwGnj1RfFX1raq6pX19R3tOe/eJ8+T2W9DRR39y71SnLUmahdneq6pqe1VV+/ZpQI3fZ9zxFvS9qj2nTcD9U5zH4/equ+66a7JdF6WbbrqJww7bsTHyiCOO6LvvKaecwk9/+tP5CEvSEjVIC9ZDwLFV9TLgcOAP22/ORoDXAC8F/jUwMlYgySlJTmlfrwGOBc6d6kBJDgTWA0dU1UuA3xkgzlOAD7ffPo4A3wNOB75TVeuq6rQkRwL7AYcC64BDkryiLf9C4JNV9dKqui3JpUmeNUW8hwK7AN8Zv62qzquqkaoa2empqwY4DUnSDMzqXtW+PzbJzcDnaVqx+lps96rJ9N6r9t77Sd8VLnoHHHAAV1xxxQ7rvvSlL/Xd99xzz2XXXXedj7AkLVGDJFgB3p/kOuBvgDXAM4HDgM9V1YNVdT9wyViBqjq3qsYSqnOAd1XVo9M41hHAZ6rq7raeewaI80rg95K8C3huVT3YZ58j2+WbNN/+7U9zEwO4raqu6jmHo9oWqr6S7ANcALylqh4bIE5JUvdme6+iqi6uqv2B3wDeN8mxFs29arlbv349J510En/7t3/Lww8/DMBJJ530+PYPfehD/OQnP+FDH/oQ73rXu/jN3/zNYYUqaQnYeYB9T6DpAndIVT3cdnNYQXMzm44R4M/SPD19L+CoJI9U1V/22TdM0S0DeIQnEsQVYyur6lNJvg4cDVyW5G3ArX3qP7Oq/miHlcla4IF+B0vycmBs/9+vqo1JdqP5hvPdvTc6SdLQzPZe9biqujzJ85PsNZZEjbMo7lVTxLcsvO51rwPg7W9/OzfddBPQjMviH5qGy9/7vd/jP/yH/8Cuu+7KKaecwkc/+tGhxSpp8RukBWsVcGd7wzoceG67/grgVUlWJFlJc7N4kqr6Z1W1tqrWAp8BTp0guYJmkPFxSfYESLK6zz5bgUPa14/PSJjkecCtVfURYCNwME3f86f3lL0MeGsbL0nWJHnGZCdfVV9vu22sa5OrXYCLabpofHqyspKkeTOre1WSF6T9JjDJy2i6f/9ogmMt+HvVZPsuN6973eu44YYbePTRRx9/P+ahhx6iqnjooYdMriTN2iAtWBcClyQZBbYANwNU1dVJNgLXArcBo8C90PRrb/eZctxVr6q6MckZwFeSPErTPeLN43b7IPDnSd4I9HakPh54Q5KHgR8Af1BV9yT5apIbgC+0fdsPAK5s76PbgTcAT+q+mORS4G19ul4cB7wC2DPJWGxvrqotg5yrJKlTs71XvQY4sb2HPAgc3zPpxQ4Wyb2KJH9H071wZZLvASdV1XzP5itJy0YmuG8MVkmysqq2J3kqcDlwclVtnnXFS8yu++xX+7zpnGGHoUVo61l9v2yXlqwk11TVyNR7DlSn96ppGBkZqdHR0WGHMaeSUFWwYRV573108beQpOVnonvVIC1YkzkvyYto+pd/whuWJGkB8l4lSZpznSRYVfX6mZRr+61v6rPplVU1UZ93SZIG5r1KY2yxkjSXumrBmpH2xrRumDFIkjQZ71WSpEEMMougJEmSJGkSJliSJEmS1BETLEmSJEnqiAmWJEmSJHXEBEuSJC1be+yxx7BDkLTEmGBJkqTlacO93HPPPcOOQtISY4IlSZIkSR0xwZIkSZKkjphgSZIkSVJHTLAkSZIkqSM7DzuA5eTFa1YxetbRww5DkiRJ0hyxBUuSJEmSOmKCJUmSJEkdMcGSJEmSpI6YYEmSJElSR0ywJEmSJKkjJliSJEmS1BETLEmSJEnqiAmWJEmSJHXEBEuSJEmSOrLzsANYTq7//r2sPf3zww5jwdt61tHDDkGSJEmaEVuwJEmSJKkjJliSJEmS1BETLEmSJEnqiAmWJEmSJHXEBEuSJEmSOmKCJUmSJEkdMcGSJEmSpI6YYEmSJElSR0ywJEmSJKkjJliSJEmS1BETLEmStOisXr2aJDssbFj1pHWDLqtXrx72qUla5EywJEnSorNt2zaqaocFeNK6QZdt27YN+cwkLXYmWJIkSZLUERMsSZIkSeqICZYkSZIkdcQES5IkSZI6YoIlSZIWjCTDDmGolvv5S0vBlAlWku1dHCjJMUmuS7IlyWiSw2ZQx4Yk75xBud2TnDpouSnqXJfkyiQ3tud1fJf1S5IkSVp85rMFaxPwkqpaB7wV+Pg8Hnt3YKAEK43Jrs9PgBOr6kDg14Bzkuw+4wglSZIkLXrTTrCSrEyyKcnmJNcnOaZn23uS3Jzki0ku6tfKVFXba+whFfA0oMbvM+54J7YtQ9cmuaDP9i8nGWlf75Vka/v6wCTfaFvKrkuyH3AW8Px23dntfqclubrd573turVJbkryMWAz8OyJ4quqb1XVLe3rO4A7gb0nOydJkiRJS9sgLVgPAcdW1cuAw4E/bFt5RoDXAC8F/jUwMlYgySlJTul5f2ySm4HP07Ri9ZXkQGA9cERVvQT4nQHiPAX4cNtSNgJ8Dzgd+E5Vrauq05IcCewHHAqsAw5J8oq2/AuBT1bVS6vqtiSXJnnWZAdMciiwC/CdPttObrtEjj76k3sHOA1JkrRcXXTRRRx00EHstNNOHHTQQVx00UXDDknSNA2SYAV4f5LrgL8B1gDPBA4DPldVD1bV/cAlYwWq6tyqOrfn/cVVtT/wG8D7JjnWEcBnquruttw9A8R5JfB7Sd4FPLeqHuyzz5Ht8k2alqr9aRIugNuq6qqemI9qW6j6SrIPcAHwlqp6bPz2qjqvqkaqamSnp64a4DQkSTPR4djhE9peDtcl+VqSl8ygjqGOHU7y75J8O0kl2Wu29Wl+XHTRRaxfv56PfvSjPPTQQ3z0ox9l/fr1JlnSIjFIgnUCTRe4Q9rWoR8CK2gSr4FU1eU0XfYm+rAPU3QhBB7hifhX9NT9KeDVwIPAZUmOmKD+M9sWrXVV9YKqOr/d9kDfgJKXt10MtyR5dbtuN5rWuHf3JmWSpCXhH4FfqqqDab4UPG8ej7073Ywd/irwK8BtHcWleXDGGWdw/vnnc/jhh/OUpzyFww8/nPPPP58zzjhj2KFJmoZBEqxVwJ1V9XCSw4HntuuvAF6VZEWSlcDR/QoneUHauUeTvIymS92PJjjWJuC4JHu2+6/us89W4JD29Wt7jvM84Naq+giwETgYuB94ek/Zy4C3tvGSZE2SZ0x28lX19Z6EbGOSXYCLaboTfnqyspKk+dfB2OGvVdW29u1VwL5THG/BjR2uqm9W1dZpXKvHu7PfddddU+0+55JMuQz7+HO1ANx0000cdtiOky0fdthh3HTTTXN63pK6sfMA+14IXJJkFNgC3AxQVVcn2QhcS/MN2ShwLzRjsNp9zqUZp3VikodpWpeO75n0YgdVdWOSM4CvJHmUpivfm8ft9kHgz5O8EfhSz/rjgTe0x/kB8AdVdU+Srya5AfhCOw7rAODK9sNsO/AG4NHxsSS5FHhbn26CxwGvAPZMMhbbm6tqS79zkiTNu7Gxw/e1PSauau9Xh/DE2OGdaRKTa+BJ961eJwFfmOhAeWLs8C9W1d0TfDE4kbGxwxe2X97tRDN2+KC2xwjZcexwgI1pxg7fTjN2+C1VdWq770T3rQlV1Xm0LXQjIyNT9SCZcxP8ebCDuUyypnP8uZKEAw44gCuuuILDDz/88fVXXHEFBxxwwNDikjR9UyZYVbWy/fdu4F9MsNsHq2pDkqcClwN/2JbpHX/1AeAD0w2sqj4BfGLcug09r2+maZ0a8+52/ZnAmX3qe/249x8GPtzn0AeN2++oCeL7U+BPJzsHSdJQjY0dfgXwGH3GDgMk2WHs8JMqaXptnNSWm8hsxw6vT7Iv8NmquqVP8tA7dhhgJU3CdTt9xg4PcGwtQOvXr+ekk07i/PPP57DDDuOKK67gpJNOsougtEgM0oI1mfOSvIhmLNQnqmpzR/VKkjRTvWOHH2675A00djjJwTTPbfz1qpqoWzvMcuxwkq/TdLG/LMnbgFv71H9mVf3RuPjWMsHYYS1er3vd6wB4+9vfzk033cQBBxzAGWec8fh6SQtbJwnW+Nah6WrHWG3qs+mVU9zIJEmaymRjh/8oyZk098GjgT8eXzjJc4DPAm+sqm9NcaxNwMVJ/mtV/SjJ6j6tWFtpuid+gwnGDrevD6bpdj9+7PD7klxYVduTrAEensY10CL1ute9zoRKWqS6asGakTaJWjfMGCRJS9Zsxw7/PrAn8LG2y94jVTVCHwt17HCS3wZ+F/g54Lokl1bV26Z3+SRJM5FhDuRcbnbdZ7/a503nDDuMBW/rWX0nopS0jCS5ZqJkpqP6V7YtQWNjh0+2e3tjZGSkRkdHh3b8JNOe5OJJ+21YBRvunZfjz5VhH1/S9E10rxpqC5YkSUPi2GFJ0pwwwZIkLTuOHV64lnvrzXI/f2kpMMGSJGmaHDssSZrKz0y9iyRJkiRpOkywJEmSJKkjJliSJEmS1BETLEmStCgl2WHpt27QZY899hjyWUla7JzkQpIkLToTzbZXG+Y3DkkazxYsSZIkSeqICZYkSZIkdcQugvPoxWtWMXrW0cMOQ5IkSdIcsQVLkiRJkjpigiVJkiRJHTHBkiRJkqSOmGBJkiRJUkdMsCRJkiSpIyZYkiRJktQREyxJkiRJ6ogJliRJkiR1xARLkiRJkjqy87ADWE6u//69rD3988MOY8HYetbRww5BkiRJ6pQtWJIkSZLUERMsSZIkSeqICZYkSZIkdcQES5IkSZI6YoIlSZIkSR0xwZIkSZKkjphgSZIkSVJHTLAkSZIkqSMmWJIkSZLUERMsSZIkSeqICZYkSZIkdcQES5IkSZI6YoIlSZIWrNWrV5Pk8YUNq3Z4P51l9erVwz4NScuICZYkSVqwtm3bRlU9vgA7vJ/Osm3btiGfhaTlxARLkiRJkjpigiVJkiRJHTHBkiRJkqSOTJlgJdnexYGS7J/kyiQ/TfLOGdaxYSZlk+ye5NSZHHOKev9Pkh8n+auu65YkSZK0+MxnC9Y9wG8DH5zHY47ZHRgowUpjqutzNvDGmQYlSZIkaWmZdoKVZGWSTUk2J7k+yTE9296T5OYkX0xyUb9Wpqq6s6quBh6e5vFOTHJdkmuTXNBn+5eTjLSv90qytX19YJJvJNnSlt8POAt4frvu7Ha/05Jc3e7z3nbd2iQ3JfkYsBl49mQxVtUm4P4pzuPkJKNJRh/9yb3TOXVJkha9JMMOYU4t9fOTNHM7D7DvQ8CxVXVfkr2Aq5JsBA4BXgO8tK1vM3ANQJJTAKrq3EGCSnIgsB74xaq6O8kgD7A4BfhwVV2YZBdgJ+B04KCqWtfWfySwH3AoEGBjklcAtwMvBN5SVae2+14KvK2q7hjkHMZU1XnAeQC77rNfzaQOSZIkSYvDIAlWgPe3ichjwBrgmcBhwOeq6kGAJJeMFRg0sepxBPCZqrq7reeeAcpeCaxPsi/w2aq6pc+3TEe2yzfb9ytpEq7bgduq6qqxHavqqJmdgiRJkqTlZpAxWCcAewOHtC1BPwRW0CReXQswVWvPIzwR/4qxlVX1KeDVwIPAZUmOmKD+M6tqXbu8oKrOb7c90Deg5OVtF8MtSV49yMlIkiRJWh4GSbBWAXdW1cNJDgee266/AnhVkhVJVgJHdxDXJuC4JHsCTNBFcCtN90SA146tTPI84Naq+giwETiYZpzU03vKXga8tY2XJGuSPGOygKrq6z0J2caZnZYkSZKkpWyQBOtCYCTJKE1r1s0A7cQVG4Frgc8Co8C90IzBGhuHleTnknwPeAfw7iTfS7JbvwNV1Y3AGcBXklwLfKjPbh8EfivJ14C9etYfD9yQZAuwP/DJqvoR8NUkNyQ5u6r+GvgUcGWS64HPsGMC9rgklyZ51gTb/g74NPDK9nx+td9+kqS51dUjRXrq+4UkjyZ57dR7P6nsgnmkSJJ17SNSbmwndTq+y/olSU+WqtnPu5BkZVVtT/JU4HLg5KraPOuKl5hd99mv9nnTOcMOY8HYelYXjZ2SlqIk11TVyAD7b6+qlR0deyfgizSTO/3PqvrMgOU3ANuraqDHkiRZC/xVVR00QJnQ3Msfm2D7zwPVjkd+Fs0kVAdU1Y8nqnNkZKRGR0cHCX2i2Ojob4wd69mwCjYMNitvV7HMdZ2SFpeJ7lVdPQfrvLbFaDPwFyZXkqRhmO0jRVpvB/4CuHMax1vQjxSpqm9V1S3t6zvac9q7T5yPP1Lkrrvumuq0py3JrJeFFMtcxCVp6RlkFsEJVdXrZ1KuHWO1qc+mV7bd+iRJGsSsHimSZA1wLM1str8w2YGyyB4pkuRQYBfgO+O39T5SZGRkpLNmma5asLowFy1YktRPJwnWTLVJ1LphxiBJWlJm+0iRc4B3VdWj0/gDetE8UiTJPsAFwJsm6k4oSerGUBMsSZI61vtIkYfbLnmDPFJkBPizNtnZCzgqySNV9Zd99p3VI0WSfJ1m5t3LkrwNuLVP/WdW1R/tsLIZqzXhI0WAsf1/v6o2pplQ6vPAu3uTMknS3OhqDJYkSQvBrB4pUlX/rKrWVtVamhlmT50guYJF8EiRtvvhxTQz6n56srKSpG6YYEmSlpJZPVJkEIvkkSLHAa8A3txOnrElybpBz1WSNH2dTNOu6XGa9h05TbukiUw09e0s6/SRItPgNO3Dq1PS4jLRvcoxWJKk5eK8JC+iGQv1CZMrSdJcMMGSJC0LPlJkfi311p2lfn6SZs4ES5KkSfhIEUnSIJzkQpIkSZI6YoIlSZIkSR0xwZIkSZKkjjgGax69eM0qRp2aXJKkgSR5/HX9p912eD8de+yxR9chSdKETLAkSdKC1W+2vtow/3FI0nTZRVCSJEmSOmKCJUmSJEkdMcGSJEmSpI6YYEmSJElSR0ywJEmSJKkjJliSJEmS1BETLEmSJEnqiAmWJEmSJHXEBw3Po+u/fy9rT/98p3VuPevoTuuTJEmSNHO2YEmSJElSR0ywJEmSJKkjJliSJEmS1BETLEmSJEnqiAmWJEmSJHXEBEuSJEmSOmKCJUmSJEkdMcGSJEmSpI6YYEmSJElSR0ywJEmSJKkjJliSJEmS1BETLEmSJEnqiAmWJEmSJHXEBEuSJEmSOmKCJUmSJEkdMcGSJEmSpI6YYEmSpDmzevVqksx4YcOqWZUfW1avXj3sSyFpmZgywUqyvYsDJTkmyXVJtiQZTXLYDOrYkOSdMyi3e5JTBy03jXr/T5IfJ/mrruuWJGkp2LZtG1U14wWYVfmxZdu2bUO+EpKWi/lswdoEvKSq1gFvBT4+j8feHRgowUpjqutzNvDGmQYlSZIkaWmZdoKVZGWSTUk2J7k+yTE9296T5OYkX0xyUb9WpqraXmNfRcHTgBq/z7jjndi2eF2b5II+27+cZKR9vVeSre3rA5N8o20puy7JfsBZwPPbdWe3+52W5Op2n/e269YmuSnJx4DNwLMni7GqNgH3T3EeJ7ctdqOP/uTeyXaVJEmStMjtPMC+DwHHVtV9SfYCrkqyETgEeA3w0ra+zcA1AElOAaiqc9v3xwJnAs8Ajp7oQEkOBNYDv1hVdycZpOP0KcCHq+rCJLsAOwGnAwe1rWckORLYDzgUCLAxySuA24EXAm+pqlPbfS8F3lZVdwwQw+Oq6jzgPIBd99lv0qRSkiRJ0uI2SIIV4P1tIvIYsAZ4JnAY8LmqehAgySVjBcYSq573FwMXt3W8D/iVCY51BPCZqrq7LXfPAHFeCaxPsi/w2aq6Jcn4fY5sl2+271fSJFy3A7dV1VU9MR81wLElSZIkLWODjME6AdgbOKRtCfohsIIm8RpIVV1O02Vvrwl2CVN0IQQe4Yn4V/TU/Sng1cCDwGVJjpig/jOral27vKCqzm+3PdA3oOTlbRfDLUlePUVskiRJkpahQRKsVcCdVfVwksOB57brrwBelWRFkpVM0PUvyQvSNiUleRmwC/CjCY61CTguyZ7t/v26CG6l6Z4I8Nqe4zwPuLWqPgJsBA6mGSf19J6ylwFvbeMlyZokz5js5Kvq6z0J2cbJ9pUkaTnq02NkQVts8UpaHAbpInghcEmSUWALcDNAVV3djsW6FrgNGAXuhSeNwXoNcGKSh2lal47vmfRiB1V1Y5IzgK8keZSmK9+bx+32QeDPk7wR+FLP+uOBN7TH+QHwB1V1T5KvJrkB+EJVnZbkAODK9sN1O/AG4NHxsUw2BivJ3wH7AyuTfA84qaoum/AKSpIkSVrSMkGOM1glycqq2p7kqcDlwMlVtXnWFS8xu+6zX+3zpnM6rXPrWRPOFSJJi1aSa6pqZNhxLEcjIyM1Ojo6o7JJGP93Rb91A9mwCjbMfhbeOYlN0rI20b1qkBasyZyX5EU0Y6E+YXIlSZIkaTnqJMGqqtfPpFw7xmpTn02vrKqJxmdJkiRJ0oI0yCQXnauqH/VMHNG7mFxJkqYtyfaO6jmmfQD9lvYh8YfNoI4NSd45g3K7Jzl10HJT1LkuyZVJbmzP6/gu65ckPdlQEyxJkhaYTcBL2seRvBX4+Dwee3dgoAQrjcnu5T8BTqyqA4FfA85JsvuMI5QkTckES5K0ZCRZmWRTks1Jrk9yTM+29yS5OckXk1zUr5Wpqrb3zHD7NKZ4JmOSE9uWoWuTXNBn+5eTjLSv90qytX19YJJvtC1l1yXZDziL5hmRW5Kc3e53WpKr233e265bm+SmJB8DNgPPnii+qvpWVd3Svr4DuJPmmZbj4zy5bbEbveuuuyY75Skl2WFZSBZybJKWjq4muZAkaSF4CDi2qu5rH2Z/VfsokUNoHhfyUpp732bgGnjSI0VIcixwJvAMJni2Y7vfgcB64Ber6u4Jntk4kVOAD1fVhUl2AXYCTgcOalvPSHIksB9wKBBgY5JXALcDLwTeUlWntvtO+EiRnngPpXkG5XfGb6uq84DzoJlFcIDzeJJ+M/UtFAs5NklLhwmWJGkpCfD+NhF5DFgDPBM4DPhcVT0IkOSSsQJjiVXP+4uBi9s63gf8ygTHOgL4TFXd3Za7Z4A4rwTWJ9kX+GxV3dLnj/0j2+Wb7fuVNAnX7cBtVXVVT8xHTXawJPsAFwBvqqrHBohTkjQguwhKkpaSE2i6wB3StgT9kOYRIgM3VVTV5TRd9vaaYJcwRRdC4BGeuNeu6Kn7U8CrgQeBy5IcMUH9Z/ZMAPWCqjq/3fZA34CSl7ddDLckeXW7bjfg88C7e5MySdLcMMGSJC0lq4A7q+rhJIcDz23XXwG8KsmKJCuZoOtfkhekbUpK8jKaLnUTzWy7CTiufeQIE3QR3ErTPRHgtT3HeR5wa1V9BNgIHAzcDzy9p+xlwFvbeEmyJskzJjv5qvp6T0K2se1+eDHwyar69GRlJUndsIugJGkpuRC4JMkosAW4GaCqrm7HYl0L3AaMAvfCk8ZgvQY4McnDNK1Lx9f4gTutqroxyRnAV5I8StOV783jdvsg8OdJ3gh8qWf98cAb2uP8APiDqronyVeT3AB8oapOS3IAcGWb820H3gA8Oj6WScZgHQe8AtgzyVhsb66qLf3OSZI0e5ngvqE5sOs++9U+bzqn0zq3njXh+GtJWrSSXFNVIx3XubKqtid5KnA5cHJVbe7yGEvByMhIjY6Ozqhskr4TSczqb40Nq2DDvTMvP0kcs45N0rI20b3KFixJ0nJxXpIX0YyF+oTJlSRpLphgSZKWhap6/UzKtWOsNvXZ9Mqqmmh81rK02FqDFlu8khYHEyxJkibRJlHrhh2HJGlxMMGaRy9es4pRx0xJkiRJS5bTtEuSJElSR0ywJEmSJKkjJliSJGlOJZnxMtvyY8see+wx5KsgablwDJYkSZozXczUVxtmH4ckzRdbsCRJkiSpIyZYkiRJktQREyxJkiRJ6ogJliRJkiR1xARLkiRJkjpigiVJkiRJHTHBkiRJkqSOmGBJkiRJUkdMsCRJkiSpIyZYkiRJktQREyxJkiRJ6ogJliRJkiR1JFU17BiWjST3A/8w7DhmYS/g7mEHMQuLOf7FHDsY/zAt1tifW1V7DzuI5SjJXcBtHVW3mH7+jHVuGOvcMNa5MWisfe9VJljzKMloVY0MO46ZMv7hWcyxg/EP02KOXYvfYvr5M9a5Yaxzw1jnRlex2kVQkiRJkjpigiVJkiRJHTHBml/nDTuAWTL+4VnMsYPxD9Nijl2L32L6+TPWuWGsc8NY50YnsToGS5IkSZI6YguWJEmSJHXEBEuSJEmSOmKC1ZEkv5bkH5J8O8npfbYnyUfa7dcledl0y86HWca/Ncn1SbYkGZ3fyKcV+/5Jrkzy0yTvHKTsfJhl/Av92p/Q/rxcl+RrSV4y3bLzYZbxD/XatzFMFf8xbexbkowmOWy6ZaWZSLI6yReT3NL+u0effZ6d5G+T3JTkxiS/M88xzvh+N99m8xk136b7mZLkF5I8muS18xnfuBimjDXJL7efnTcm+cp8x9gTx1Q/A6uSXJLk2jbWtwwpzv+Z5M4kN0ywfSH9Xk0V6+x/r6rKZZYLsBPwHeB5wC7AtcCLxu1zFPAFIMA/B74+3bILOf5221ZgrwV87Z8B/AJwBvDOQcou5PgXybX/v4A92te/vgh/7vvGP+xrP0D8K3lirO3BwM0L5fq7LM0F+C/A6e3r04EP9NlnH+Bl7eunA9+ar5+/2d7v5vlazuozaqHF2rPfl4BLgdcu1FiB3YG/B57Tvn/GAo7198Z+z4C9gXuAXYYQ6yuAlwE3TLB9QfxeTTPWWf9e2YLVjUOBb1fVrVX1/wF/Bhwzbp9jgE9W4ypg9yT7TLPsXJtN/MM2ZexVdWdVXQ08PGjZeTCb+IdtOrF/raq2tW+vAvadbtl5MJv4F4LpxL+92jsE8DSgpltWmqFjgE+0rz8B/Mb4Harqn6pqc/v6fuAmYM08xbeY7neL6TNqup8pbwf+ArhzPoMbZzqxvh74bFXdDs19eJ5jHDOdWAt4epLQfKl2D/DI/IYJVXV5e+yJLJTfqylj7eL3ygSrG2uA7/a8/x5PvllMtM90ys612cQPzS/3Xye5JsnJcxZlf7O5fovl2k9mMV37k2i+vZpJ2bkwm/hhuNcephl/kmOT3Ax8HnjrIGWlGXhmVf0TNIkUTQv8hJKsBV4KfH3uQwNmf7+bT7P9jJpPU8aaZA1wLHDuPMbVz3Su688DeyT5cvsZf+K8Rbej6cT634ADgDuA64HfqarH5ie8gSyU36tBzej3auc5CGQ5Sp914+e/n2if6ZSda7OJH+AXq+qOJM8Avpjk5vbbgfkwm+u3WK79ZBbFtU9yOM2H1NgYoEV17fvED8O99jDN+KvqYuDiJK8A3gf8ynTLSv0k+Rvg5/psWj9gPStpWjP+fVXd10Vs0zlsn3WD3O/m02w/o+bTdGI9B3hXVT3aNLYMzXRi3Rk4BHgl8LPAlUmuqqpvzXVw40wn1l8FtgBHAM+nuR/93Tz+Tk3XQvm9mrbZ/F6ZYHXje8Cze97vS/NNwnT22WUaZefabOKnqsb+vTPJxTRN2vP1h+Z0Yp+Lsl2ZVQyL4donORj4OPDrVfWjQcrOsdnEP+xrDwNew6q6PMnzk+w1aFmpV1X9ykTbkvwwyT5V9U9t95++XauSPIUmubqwqj47R6H2M6v73Tyb1WfUPJtOrCPAn7XJ1V7AUUkeqaq/nJcInzDdn4G7q+oB4IEklwMvoRkvOJ+mE+tbgLPa7uDfTvKPwP7AN+YnxGlbKL9X0zLr36tBB2259B0MtzNwK/DPeGIQ4oHj9jmaHQf3fWO6ZRd4/E8Dnt7z+mvAry2k2Hv23cCOk1wsims/SfwL/toDzwG+DfxfMz3vBRr/UK/9APG/gCcmuXgZ8P32d3jo199laS7A2ew4ycV/6bNPgE8C5wwhvhnf7xZorH0/oxZirOP2/xOGN8nFdK7rAcCmdt+nAjcABy3QWP8HsKF9/cz2c35Yk1+tZeKJIxbE79U0Y53179XQTmypLTSzo3yLZraX9e26U4BT2tcB/nu7/XpgZLKyiyV+mpltrm2XG4cR/zRi/zmab07uA37cvt5tEV37vvEvkmv/cWAbTfeFLcDoIvu57xv/Qrj204z/XW18W4ArgcMW0vV3WXoLsCfNH6a3tP+ubtc/C7i0fX0YTdeg63p+t46axxhnfL8ewvWc8WfsQot13L5/wpASrOnGCpxGM5PgDTTdWBdkrO3v1l+3P6s3AG8YUpwXAf9EMyHX92i61i3U36upYp3179XYN5uSJEmSpFlyFkFJkiRJ6ogJliRJkiR1xARLkiRJkjpigiVJkiRJHTHBkiRJkqSOmGBJkiRJUkdMsCRJkiSpI/8/kHY4RS7IbKEAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 864x432 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "feature_importance = regressor.feature_importances_\n",
    "sorted_idx = np.argsort(feature_importance)\n",
    "pos = np.arange(sorted_idx.shape[0]) + 0.5\n",
    "fig = plt.figure(figsize=(12, 6))\n",
    "plt.subplot(1, 2, 1)\n",
    "plt.barh(pos, feature_importance[sorted_idx], align=\"center\")\n",
    "plt.yticks(pos, np.array(clustered_data_test.cluster_labels)[sorted_idx])\n",
    "plt.title(\"Feature Importance\")\n",
    "\n",
    "# also include the permutation impotance score of features\n",
    "from sklearn.inspection import permutation_importance\n",
    "\n",
    "result = permutation_importance(\n",
    "    regressor, clustered_data_test, target_timeseries_test,\n",
    "    n_repeats=10, random_state=42, n_jobs=2\n",
    ")\n",
    "sorted_idx = result.importances_mean.argsort()\n",
    "plt.subplot(1, 2, 2)\n",
    "plt.boxplot(\n",
    "    result.importances[sorted_idx].T,\n",
    "    vert=False,\n",
    "    labels=np.array(clustered_data_test.cluster_labels)[sorted_idx],\n",
    ")\n",
    "plt.title(\"Permutation Importance (test set)\")\n",
    "fig.tight_layout()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# from sklearn.model_selection import train_test_split\n",
    "# X_train, X_test, y_train, y_test = train_test_split(\n",
    "#     X, y, test_size=0.1, random_state=13\n",
    "# )"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.10"
  },
  "vscode": {
   "interpreter": {
    "hash": "e7604e8ec5f09e490e10161e37a4725039efd3ab703d81b1b8a1e00d6741866c"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}