10_operational_forecast_example.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Operational forecast example\n",
    "\n",
    "This takes the XGBoost and Random Forest quantile regression models trained in the last notebook, [09_machine_learning_models.ipynb](09_machine_learning_models.ipynb), and ilustrates how forecasts could be produced on demand and fun through those models to produce an on-demand probabilistic forecast.\n",
    "\n",
    "To get all 30 GEFS members and the HRRR forecasts, this process should take 15-30 minutes. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "import numpy as np\n",
    "import pvlib\n",
    "import datetime\n",
    "from sffp_tools.herbie_gefs_tools import get_tcc_fcasts_fastherbie, get_tcc_fcast\n",
    "from sffp_tools.herbie_hrrr_tools import get_hrrr_fcasts_fastherbie, hrrr_ds_to_power\n",
    "import matplotlib.pyplot as plt\n",
    "import xgboost as xgb\n",
    "from quantile_forest import RandomForestQuantileRegressor\n",
    "import joblib"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load plant specs:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "plants_df = pd.read_csv('pv_fleet_metadata_rev3.csv') # read in as a Pandas DataFrame\n",
    "plants_df = plants_df[plants_df['ac_capacity']>=10].copy().reset_index()\n",
    "num_plants = len(plants_df.index) # number of plants"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate initialization time for GEFS:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2024-06-06 06:00:00\n"
     ]
    }
   ],
   "source": [
    "# get delay_to_first_forecast and update_freq from \n",
    "# https://github.com/SolarArbiter/solarforecastarbiter-core/blob/master/solarforecastarbiter/io/fetch/nwp.py\n",
    "delay_to_first_forecast = '280min'\n",
    "update_freq = '6h'\n",
    "\n",
    "time_now = pd.Timestamp.utcnow().tz_localize(None) # current time in UTC, but remove timezone awareness\n",
    "init_time = (time_now - pd.Timedelta(delay_to_first_forecast)).floor(update_freq)\n",
    "print(init_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get GEFS data. \n",
    "\n",
    "Note: The forecast range, `fxx_range`, has been expanded from 27-45 hours ahead to now cover 12-45 hours ahead, with the idea that this might be run in the morning _or_ afternoon/evening, and this allows it to still cover the full day ahead. You could adjust the range accordingly to speed things up a bit."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "downloading new data\n",
      "member 1 download complete\n",
      "member 2 download complete\n",
      "member 3 download complete\n",
      "member 4 download complete\n",
      "member 5 download complete\n",
      "member 6 download complete\n",
      "member 7 download complete\n",
      "member 8 download complete\n",
      "member 9 download complete\n",
      "member 10 download complete\n",
      "member 11 download complete\n",
      "member 12 download complete\n",
      "member 13 download complete\n",
      "member 14 download complete\n",
      "member 15 download complete\n",
      "member 16 download complete\n",
      "member 17 download complete\n",
      "member 18 download complete\n",
      "member 19 download complete\n",
      "member 20 download complete\n",
      "member 21 download complete\n",
      "member 22 download complete\n",
      "member 23 download complete\n",
      "member 24 download complete\n",
      "member 25 download complete\n",
      "member 26 download complete\n",
      "member 27 download complete\n",
      "member 28 download complete\n",
      "member 29 download complete\n",
      "member 30 download complete\n",
      "downloaded and processed \n"
     ]
    },
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>time_utc</th>\n",
       "      <th>step</th>\n",
       "      <th>tcc_std</th>\n",
       "      <th>tcc_p01</th>\n",
       "      <th>tcc_p02</th>\n",
       "      <th>tcc_p03</th>\n",
       "      <th>tcc_p04</th>\n",
       "      <th>tcc_p05</th>\n",
       "      <th>tcc_p06</th>\n",
       "      <th>tcc_p07</th>\n",
       "      <th>...</th>\n",
       "      <th>tcc_p21</th>\n",
       "      <th>tcc_p22</th>\n",
       "      <th>tcc_p23</th>\n",
       "      <th>tcc_p24</th>\n",
       "      <th>tcc_p25</th>\n",
       "      <th>tcc_p26</th>\n",
       "      <th>tcc_p27</th>\n",
       "      <th>tcc_p28</th>\n",
       "      <th>tcc_p29</th>\n",
       "      <th>tcc_p30</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>valid_time_utc_end_of_interval</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2024-06-06 16:00:00+00:00</th>\n",
       "      <td>2024-06-06 06:00:00</td>\n",
       "      <td>0 days 12:00:00</td>\n",
       "      <td>13.887529</td>\n",
       "      <td>20.078788</td>\n",
       "      <td>40.736216</td>\n",
       "      <td>31.561795</td>\n",
       "      <td>46.587923</td>\n",
       "      <td>34.282701</td>\n",
       "      <td>5.552620</td>\n",
       "      <td>15.997928</td>\n",
       "      <td>...</td>\n",
       "      <td>35.884132</td>\n",
       "      <td>42.630439</td>\n",
       "      <td>49.657024</td>\n",
       "      <td>30.067244</td>\n",
       "      <td>24.030164</td>\n",
       "      <td>53.952049</td>\n",
       "      <td>23.072680</td>\n",
       "      <td>26.706509</td>\n",
       "      <td>8.574980</td>\n",
       "      <td>18.012486</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2024-06-06 17:00:00+00:00</th>\n",
       "      <td>2024-06-06 06:00:00</td>\n",
       "      <td>0 days 12:00:00</td>\n",
       "      <td>13.887529</td>\n",
       "      <td>20.078788</td>\n",
       "      <td>40.736216</td>\n",
       "      <td>31.561795</td>\n",
       "      <td>46.587923</td>\n",
       "      <td>34.282701</td>\n",
       "      <td>5.552620</td>\n",
       "      <td>15.997928</td>\n",
       "      <td>...</td>\n",
       "      <td>35.884132</td>\n",
       "      <td>42.630439</td>\n",
       "      <td>49.657024</td>\n",
       "      <td>30.067244</td>\n",
       "      <td>24.030164</td>\n",
       "      <td>53.952049</td>\n",
       "      <td>23.072680</td>\n",
       "      <td>26.706509</td>\n",
       "      <td>8.574980</td>\n",
       "      <td>18.012486</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2024-06-06 18:00:00+00:00</th>\n",
       "      <td>2024-06-06 06:00:00</td>\n",
       "      <td>0 days 12:00:00</td>\n",
       "      <td>13.887529</td>\n",
       "      <td>20.078788</td>\n",
       "      <td>40.736216</td>\n",
       "      <td>31.561795</td>\n",
       "      <td>46.587923</td>\n",
       "      <td>34.282701</td>\n",
       "      <td>5.552620</td>\n",
       "      <td>15.997928</td>\n",
       "      <td>...</td>\n",
       "      <td>35.884132</td>\n",
       "      <td>42.630439</td>\n",
       "      <td>49.657024</td>\n",
       "      <td>30.067244</td>\n",
       "      <td>24.030164</td>\n",
       "      <td>53.952049</td>\n",
       "      <td>23.072680</td>\n",
       "      <td>26.706509</td>\n",
       "      <td>8.574980</td>\n",
       "      <td>18.012486</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2024-06-06 19:00:00+00:00</th>\n",
       "      <td>2024-06-06 06:00:00</td>\n",
       "      <td>0 days 15:00:00</td>\n",
       "      <td>16.598369</td>\n",
       "      <td>37.864650</td>\n",
       "      <td>75.886312</td>\n",
       "      <td>32.135484</td>\n",
       "      <td>74.036111</td>\n",
       "      <td>66.590157</td>\n",
       "      <td>17.947986</td>\n",
       "      <td>46.709980</td>\n",
       "      <td>...</td>\n",
       "      <td>46.458736</td>\n",
       "      <td>53.754433</td>\n",
       "      <td>69.790437</td>\n",
       "      <td>54.601647</td>\n",
       "      <td>50.219008</td>\n",
       "      <td>53.582972</td>\n",
       "      <td>53.429379</td>\n",
       "      <td>68.298119</td>\n",
       "      <td>24.553454</td>\n",
       "      <td>66.218685</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2024-06-06 20:00:00+00:00</th>\n",
       "      <td>2024-06-06 06:00:00</td>\n",
       "      <td>0 days 15:00:00</td>\n",
       "      <td>16.598369</td>\n",
       "      <td>37.864650</td>\n",
       "      <td>75.886312</td>\n",
       "      <td>32.135484</td>\n",
       "      <td>74.036111</td>\n",
       "      <td>66.590157</td>\n",
       "      <td>17.947986</td>\n",
       "      <td>46.709980</td>\n",
       "      <td>...</td>\n",
       "      <td>46.458736</td>\n",
       "      <td>53.754433</td>\n",
       "      <td>69.790437</td>\n",
       "      <td>54.601647</td>\n",
       "      <td>50.219008</td>\n",
       "      <td>53.582972</td>\n",
       "      <td>53.429379</td>\n",
       "      <td>68.298119</td>\n",
       "      <td>24.553454</td>\n",
       "      <td>66.218685</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 33 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                                          time_utc            step    tcc_std  \\\n",
       "valid_time_utc_end_of_interval                                                  \n",
       "2024-06-06 16:00:00+00:00      2024-06-06 06:00:00 0 days 12:00:00  13.887529   \n",
       "2024-06-06 17:00:00+00:00      2024-06-06 06:00:00 0 days 12:00:00  13.887529   \n",
       "2024-06-06 18:00:00+00:00      2024-06-06 06:00:00 0 days 12:00:00  13.887529   \n",
       "2024-06-06 19:00:00+00:00      2024-06-06 06:00:00 0 days 15:00:00  16.598369   \n",
       "2024-06-06 20:00:00+00:00      2024-06-06 06:00:00 0 days 15:00:00  16.598369   \n",
       "\n",
       "                                  tcc_p01    tcc_p02    tcc_p03    tcc_p04  \\\n",
       "valid_time_utc_end_of_interval                                               \n",
       "2024-06-06 16:00:00+00:00       20.078788  40.736216  31.561795  46.587923   \n",
       "2024-06-06 17:00:00+00:00       20.078788  40.736216  31.561795  46.587923   \n",
       "2024-06-06 18:00:00+00:00       20.078788  40.736216  31.561795  46.587923   \n",
       "2024-06-06 19:00:00+00:00       37.864650  75.886312  32.135484  74.036111   \n",
       "2024-06-06 20:00:00+00:00       37.864650  75.886312  32.135484  74.036111   \n",
       "\n",
       "                                  tcc_p05    tcc_p06    tcc_p07  ...  \\\n",
       "valid_time_utc_end_of_interval                                   ...   \n",
       "2024-06-06 16:00:00+00:00       34.282701   5.552620  15.997928  ...   \n",
       "2024-06-06 17:00:00+00:00       34.282701   5.552620  15.997928  ...   \n",
       "2024-06-06 18:00:00+00:00       34.282701   5.552620  15.997928  ...   \n",
       "2024-06-06 19:00:00+00:00       66.590157  17.947986  46.709980  ...   \n",
       "2024-06-06 20:00:00+00:00       66.590157  17.947986  46.709980  ...   \n",
       "\n",
       "                                  tcc_p21    tcc_p22    tcc_p23    tcc_p24  \\\n",
       "valid_time_utc_end_of_interval                                               \n",
       "2024-06-06 16:00:00+00:00       35.884132  42.630439  49.657024  30.067244   \n",
       "2024-06-06 17:00:00+00:00       35.884132  42.630439  49.657024  30.067244   \n",
       "2024-06-06 18:00:00+00:00       35.884132  42.630439  49.657024  30.067244   \n",
       "2024-06-06 19:00:00+00:00       46.458736  53.754433  69.790437  54.601647   \n",
       "2024-06-06 20:00:00+00:00       46.458736  53.754433  69.790437  54.601647   \n",
       "\n",
       "                                  tcc_p25    tcc_p26    tcc_p27    tcc_p28  \\\n",
       "valid_time_utc_end_of_interval                                               \n",
       "2024-06-06 16:00:00+00:00       24.030164  53.952049  23.072680  26.706509   \n",
       "2024-06-06 17:00:00+00:00       24.030164  53.952049  23.072680  26.706509   \n",
       "2024-06-06 18:00:00+00:00       24.030164  53.952049  23.072680  26.706509   \n",
       "2024-06-06 19:00:00+00:00       50.219008  53.582972  53.429379  68.298119   \n",
       "2024-06-06 20:00:00+00:00       50.219008  53.582972  53.429379  68.298119   \n",
       "\n",
       "                                  tcc_p29    tcc_p30  \n",
       "valid_time_utc_end_of_interval                        \n",
       "2024-06-06 16:00:00+00:00        8.574980  18.012486  \n",
       "2024-06-06 17:00:00+00:00        8.574980  18.012486  \n",
       "2024-06-06 18:00:00+00:00        8.574980  18.012486  \n",
       "2024-06-06 19:00:00+00:00       24.553454  66.218685  \n",
       "2024-06-06 20:00:00+00:00       24.553454  66.218685  \n",
       "\n",
       "[5 rows x 33 columns]"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "num_members = 30 # Number of ensemble members to pull, max of 30\n",
    "freq='1D' # how often to pull a new forecast\n",
    "fxx_interval = 3 # time step of the forecast, GEFS has 3 hour intervals\n",
    "fxx_range = range(12,46,fxx_interval)\n",
    "\n",
    "# lat, lon extent for region of interest\n",
    "# (left lon, right lon, lower lat, upper lat)\n",
    "region_extent = [-92, -80, 29, 36]\n",
    "\n",
    "start = init_time\n",
    "end = init_time\n",
    "\n",
    "save_dir = 'nc_files/gefs/'\n",
    "remove_grib=True\n",
    "w_avg_tcc = get_tcc_fcasts_fastherbie(start, end, freq, fxx_range, plants_df,\n",
    "                              region_extent=region_extent, num_members=num_members,\n",
    "                              save=False, save_dir='', attempts=5, remove_grib=True)\n",
    "\n",
    "# clean up indices and localize to UTC\n",
    "w_avg_tcc = w_avg_tcc.reset_index().set_index('valid_time_utc_end_of_interval').tz_localize('UTC')\n",
    "w_avg_tcc.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Calculate initialization time for HRRR:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "2024-06-06 06:00:00\n"
     ]
    }
   ],
   "source": [
    "# get delay_to_first_forecast and update_freq from \n",
    "# https://github.com/SolarArbiter/solarforecastarbiter-core/blob/master/solarforecastarbiter/io/fetch/nwp.py\n",
    "delay_to_first_forecast = '45min' \n",
    "update_freq = '6h' # HRRR updates hourly, but only goes out 48 hours at 00Z, 06Z, 12Z, and 18Z forecasts\n",
    "\n",
    "time_now = pd.Timestamp(datetime.datetime.now().isoformat())\n",
    "init_time = (time_now - pd.Timedelta(delay_to_first_forecast)).floor(update_freq)\n",
    "print(init_time)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Get HRRR data.\n",
    "\n",
    "Note: Just list the GEFS range above, the values of `fxx_range` have been modified to provide wider coverage, and could be adjust to get a narrower forecast horizon if needed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "downloading new data\n",
      "attempt 1, pausing for 1 min\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\core.py:1096: UserWarning: Will not remove GRIB file because it previously existed.\n",
      "  warnings.warn(\"Will not remove GRIB file because it previously existed.\")\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "variable group 1 download complete\n",
      "variable group 2 download complete\n",
      "variable group 3 download complete\n",
      "downloaded and processed \n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:258: UserWarning: More than one time coordinate present for variable  \"t2m_c\".\n",
      "  ds = ds.metpy.parse_cf()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:258: UserWarning: More than one time coordinate present for variable  \"wspd\".\n",
      "  ds = ds.metpy.parse_cf()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"dswrf\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"t2m\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"u10\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"v10\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"t2m_c\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"wspd\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:258: UserWarning: More than one time coordinate present for variable  \"t2m_c\".\n",
      "  ds = ds.metpy.parse_cf()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:258: UserWarning: More than one time coordinate present for variable  \"wspd\".\n",
      "  ds = ds.metpy.parse_cf()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"dswrf\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"t2m\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"u10\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"v10\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"t2m_c\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"wspd\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:258: UserWarning: More than one time coordinate present for variable  \"t2m_c\".\n",
      "  ds = ds.metpy.parse_cf()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:258: UserWarning: More than one time coordinate present for variable  \"wspd\".\n",
      "  ds = ds.metpy.parse_cf()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"dswrf\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"t2m\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"u10\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"v10\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"t2m_c\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"wspd\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:258: UserWarning: More than one time coordinate present for variable  \"t2m_c\".\n",
      "  ds = ds.metpy.parse_cf()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:258: UserWarning: More than one time coordinate present for variable  \"wspd\".\n",
      "  ds = ds.metpy.parse_cf()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"dswrf\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"t2m\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"u10\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"v10\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"t2m_c\".\n",
      "  ds = ds.metpy.assign_y_x()\n",
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\herbie\\accessors.py:262: UserWarning: More than one time coordinate present for variable  \"wspd\".\n",
      "  ds = ds.metpy.assign_y_x()\n"
     ]
    }
   ],
   "source": [
    "freq = '24h'\n",
    "fxx_range = range(12,45)\n",
    "query_len_days = 3\n",
    "\n",
    "save_dir = 'nc_files/hrrr/'\n",
    "remove_grib=True\n",
    "region_extent=[220,440,1080,1440]\n",
    "\n",
    "start = init_time\n",
    "end = init_time\n",
    "\n",
    "ds, _, _, _, _ = get_hrrr_fcasts_fastherbie(start, end, freq, fxx_range, plants_df,\n",
    "                              region_extent=region_extent, save=False, \n",
    "                              save_dir='', attempts=5, transposition_model='perez',\n",
    "                              decomposition_model='erbs', coarsen_type='none')\n",
    "\n",
    "window_size = 10\n",
    "\n",
    "# coarsened mean\n",
    "ds_mean = ds.coarsen(x=window_size, y=window_size, boundary='trim').mean()\n",
    "p_ac_all_fcast, _, p_ac_all_cs, _ = hrrr_ds_to_power(ds_mean,\n",
    "                                plants_df,\n",
    "                                transposition_model = 'perez',\n",
    "                                decomposition_model = 'erbs',\n",
    "                                eta_inv_nom = 0.98)\n",
    "\n",
    "# coarsened max\n",
    "ds_max = ds.coarsen(x=window_size, y=window_size, boundary='trim').max()\n",
    "p_ac_all_fcast_max, _, _, _ = hrrr_ds_to_power(ds_max,\n",
    "                                plants_df,\n",
    "                                transposition_model = 'perez',\n",
    "                                decomposition_model = 'erbs',\n",
    "                                eta_inv_nom = 0.98)\n",
    "\n",
    "# coarsened mmin\n",
    "ds_min = ds.coarsen(x=window_size, y=window_size, boundary='trim').min()\n",
    "p_ac_all_fcast_min, _, _, _ = hrrr_ds_to_power(ds_min,\n",
    "                                plants_df,\n",
    "                                transposition_model = 'perez',\n",
    "                                decomposition_model = 'erbs',\n",
    "                                eta_inv_nom = 0.98)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Load models that we trained previously:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "booster = xgb.Booster()\n",
    "booster.load_model(\"regression_models/xgboost_model.json\")\n",
    "qrf = joblib.load('regression_models/qrf_file.joblib')\n",
    "\n",
    "\n",
    "pi_list_xgb = [.9995, .999, .99, .96, .57] \n",
    "pi_list_rf = [.995, .99, .98, .95, .5]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Prep the data to get ready to feed into the models:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "p_ac_all_cs.rename(columns={'power_ac': 'p_ac_all_cs'}, inplace=True)\n",
    "p_ac_all_fcast.rename(columns={'power_ac': 'p_ac_all_fcast'}, inplace=True)\n",
    "p_ac_all_fcast_max.rename(columns={'power_ac': 'p_ac_all_fcast_max'}, inplace=True)\n",
    "p_ac_all_fcast_min.rename(columns={'power_ac': 'p_ac_all_fcast_min'}, inplace=True)\n",
    "\n",
    "# list of dataframes to combine first\n",
    "dfs = [\n",
    "    p_ac_all_fcast,\n",
    "    p_ac_all_fcast_max,\n",
    "    p_ac_all_fcast_min,\n",
    "    p_ac_all_cs,\n",
    "    w_avg_tcc.tcc_std,\n",
    "    ]\n",
    "\n",
    "# combine all dataframes and drop na rows\n",
    "df = pd.concat(dfs, axis=1).dropna()\n",
    "\n",
    "df.index.name = 'time_center_labeled'\n",
    "\n",
    "# calculate normalized max-min\n",
    "df['max_min_norm'] = (df['p_ac_all_fcast_max'] - df['p_ac_all_fcast_min'])/df['p_ac_all_cs']\n",
    "\n",
    "# normalize forecast to clear sky\n",
    "df['csi_fcast'] = df['p_ac_all_fcast']/df['p_ac_all_cs']\n",
    "\n",
    "latitude, longitude = 32.8 , -83.6\n",
    "solpos = pvlib.solarposition.get_solarposition(time=df.index,latitude=latitude, longitude=longitude)\n",
    "\n",
    "df['hour'] = df.index.hour\n",
    "df['month'] = df.index.month\n",
    "df['el'] = solpos.apparent_elevation\n",
    "df['zen'] = solpos.apparent_zenith\n",
    "df['az'] = solpos.azimuth\n",
    "\n",
    "filtered_df = df.dropna()\n",
    "approx_nameplate = plants_df['ac_capacity'].sum() \n",
    "filtered_df = filtered_df[filtered_df['p_ac_all_cs'] >= approx_nameplate*.001]\n",
    "\n",
    "X=filtered_df[['tcc_std']].copy()\n",
    "X['max_min_norm'] = filtered_df['max_min_norm']\n",
    "X['csi_fcast'] = filtered_df['csi_fcast']\n",
    "X['zen'] = filtered_df['zen']\n",
    "X['az'] = filtered_df['az']\n",
    "X['month'] = filtered_df['month']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the XGboost model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "# predict\n",
    "y_pred = booster.inplace_predict(X)\n",
    "\n",
    "# create percentile list:\n",
    "pi_list = pi_list_xgb\n",
    "percentile_list = [(1-x)/2 for x in pi_list] + [.5] + [1-(1-x)/2 for x in reversed(pi_list)]\n",
    "percentile_list = [round(elem, 6) for elem in percentile_list]\n",
    "num_intervals = len(pi_list)\n",
    "\n",
    "# create a results dataframe\n",
    "result_df = X.copy()\n",
    "\n",
    "# merge orginal DF, marking duplicate columns and then dropping them\n",
    "result_df = result_df.merge(df, left_index=True, right_index=True, suffixes=('', '_drop'))\n",
    "result_df.drop(result_df.filter(regex='_drop$').columns, axis=1, inplace=True)\n",
    "\n",
    "# adding CS-normalized quantiles/percentiles, and then un-normalizing them to get power\n",
    "for i in range(len(percentile_list)):\n",
    "    result_df[str(percentile_list[i])+'_csi']=y_pred[:,i]\n",
    "for i in range(len(percentile_list)):\n",
    "    result_df[str(percentile_list[i])+'_p_ac']=y_pred[:,i]*result_df['p_ac_all_cs']\n",
    "\n",
    "# fill in overnight \n",
    "result_df = result_df.reset_index().set_index('time_center_labeled').asfreq('1h')\n",
    "\n",
    "# save results into a new dataframe\n",
    "result_df_XGB = result_df.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Run the Random Forest model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "c:\\Users\\willh\\miniforge3\\envs\\pvlib_herbie_env\\lib\\site-packages\\sklearn\\base.py:486: UserWarning: X has feature names, but RandomForestQuantileRegressor was fitted without feature names\n",
      "  warnings.warn(\n"
     ]
    }
   ],
   "source": [
    "# predict\n",
    "y_pred = qrf.predict(X, quantiles=percentile_list)\n",
    "\n",
    "# create percentile list:\n",
    "pi_list = pi_list_rf\n",
    "percentile_list = [(1-x)/2 for x in pi_list] + [.5] + [1-(1-x)/2 for x in reversed(pi_list)]\n",
    "percentile_list = [round(elem, 6) for elem in percentile_list]\n",
    "num_intervals = len(pi_list)\n",
    "\n",
    "# create a results dataframe\n",
    "result_df = X.copy()\n",
    "\n",
    "# merge orginal DF, marking duplicate columns and then dropping them\n",
    "result_df = result_df.merge(df, left_index=True, right_index=True, suffixes=('', '_drop'))\n",
    "result_df.drop(result_df.filter(regex='_drop$').columns, axis=1, inplace=True)\n",
    "\n",
    "# adding CS-normalized quantiles/percentiles, and then un-normalizing them to get power\n",
    "for i in range(len(percentile_list)):\n",
    "    result_df[str(percentile_list[i])+'_csi']=y_pred[:,i]\n",
    "for i in range(len(percentile_list)):\n",
    "    result_df[str(percentile_list[i])+'_p_ac']=y_pred[:,i]*result_df['p_ac_all_cs']\n",
    "\n",
    "# fill in overnight \n",
    "result_df = result_df.reset_index().set_index('time_center_labeled').asfreq('1h')\n",
    "\n",
    "# save results into a new dataframe\n",
    "result_df_RF = result_df.copy()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Create the hybrid results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "extreme lower bound from rf 0.0025_p_ac\n",
      "lower bounds 0.005_p_ac 0.0005_p_ac\n",
      "lower bounds 0.01_p_ac 0.005_p_ac\n",
      "lower bounds 0.025_p_ac 0.02_p_ac\n",
      "lower bounds 0.25_p_ac 0.215_p_ac\n",
      "mid 0.5_p_ac 0.5_p_ac\n",
      "upper bounds 0.75_p_ac 0.785_p_ac\n",
      "upper bounds 0.975_p_ac 0.98_p_ac\n",
      "upper bounds 0.99_p_ac 0.995_p_ac\n",
      "upper bounds 0.995_p_ac 0.9995_p_ac\n",
      "extreme upper bound from rf 0.9975_p_ac\n"
     ]
    }
   ],
   "source": [
    "# create new hyrbid dataframe, starting with RF because PIs were accurate so we can keep column names\n",
    "result_df_hybrid = result_df_RF.copy()\n",
    "\n",
    "# find columns that end in `_p_ac` (e.g., '0.015_p_ac')\n",
    "rf_columns = result_df_RF.filter(regex='_p_ac$').columns\n",
    "xgb_columns = result_df_XGB.filter(regex='_p_ac$').columns\n",
    "\n",
    "# loop through, selecting max of lower bounds and min of upper bounds\n",
    "for i in range(len(rf_columns)):\n",
    "    if i == 0: # first\n",
    "        print('extreme lower bound from rf', rf_columns[i])\n",
    "        result_df_hybrid[rf_columns[i]] = result_df_RF[rf_columns[i]]\n",
    "    elif i == len(rf_columns)-1: # last\n",
    "        print('extreme upper bound from rf', rf_columns[i])\n",
    "        result_df_hybrid[rf_columns[i]] = result_df_RF[rf_columns[i]]\n",
    "    elif i < (len(rf_columns)-1)/2:\n",
    "        print('lower bounds', rf_columns[i], xgb_columns[i])\n",
    "        result_df_hybrid[rf_columns[i]] = np.maximum(result_df_XGB[xgb_columns[i]], result_df_RF[rf_columns[i]])\n",
    "        # result_df_hybrid_2023[rf_columns[i]] = pd.DataFrame([result_df_XGB_2023[xgb_columns[i]], result_df_RF_2023[rf_columns[i]]]).max(axis=0)\n",
    "    elif i == (len(rf_columns)-1)/2:\n",
    "        print('mid', rf_columns[i], xgb_columns[i])\n",
    "        result_df_hybrid[rf_columns[i]] = result_df_RF[rf_columns[i]]\n",
    "    else:\n",
    "        print('upper bounds', rf_columns[i], xgb_columns[i])\n",
    "        result_df_hybrid[rf_columns[i]] = np.minimum(result_df_XGB[xgb_columns[i]], result_df_RF[rf_columns[i]])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot the final model results:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAzoAAAGnCAYAAAB/1iIAAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguNCwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8fJSN1AAAACXBIWXMAAA9hAAAPYQGoP6dpAACr2ElEQVR4nOzdeXxU5dn/8c8kM5NMAgmb2TAssiqgshO0LMoiZRFcoA8KWhCsC4ogKq0oFAHF5xEqVGr9oSgoULW4tJatFRRZxEhaEESgQUAJQQkJJJPZMr8/jhkTkkDOZE++79frvMjMnOWekJlzrnPf93VZ/H6/HxERERERkVokpKobICIiIiIiUt4U6IiIiIiISK2jQEdERERERGodBToiIiIiIlLrKNAREREREZFaR4GOiIiIiIjUOgp0RERERESk1lGgIyIiIiIitY4CHRERERERqXUU6IiIiIiISK1TbQKdBQsWYLFYmDp1auA5v9/P7NmzSUhIwOFw0K9fP7766qtC27lcLqZMmUKTJk2IjIxkxIgRnDhxotA6GRkZjBs3jujoaKKjoxk3bhxnz56thHclIiIiIiJVoVoEOrt37+bPf/4zV199daHnFy5cyAsvvMDSpUvZvXs3cXFxDBw4kHPnzgXWmTp1KuvWrWPNmjVs27aN8+fPM2zYMHw+X2CdsWPHkpKSwvr161m/fj0pKSmMGzeu0t6fiIiIiIhULovf7/dXZQPOnz9Ply5deOmll3jmmWe49tprWbx4MX6/n4SEBKZOncrjjz8OGL03sbGxPPfcc9x7771kZmZy2WWXsXLlSsaMGQPA999/T2JiIh999BGDBw/mwIEDXHXVVezcuZOePXsCsHPnTpKSkvj6669p165dqdqZl5fH999/T/369bFYLBXzyxAREZFy5ff7OXfuHAkJCYSEVIv7uyJSSaxV3YAHHniAoUOHMmDAAJ555pnA86mpqaSlpTFo0KDAc2FhYfTt25ft27dz7733kpycjMfjKbROQkICHTt2ZPv27QwePJgdO3YQHR0dCHIAevXqRXR0NNu3by8x0HG5XLhcrsDj7777jquuuqo837qIiIhUkuPHj3P55ZdXdTNEpBJVaaCzZs0avvzyS3bv3l3ktbS0NABiY2MLPR8bG8u3334bWMdut9OwYcMi6+Rvn5aWRkxMTJH9x8TEBNYpzoIFC5gzZ06R548fP05UVNQl3pmIiIhUB1lZWSQmJlK/fv2qboqIVLIqC3SOHz/Oww8/zMaNGwkPDy9xvQuHifn9/ksOHbtwneLWv9R+Zs6cybRp0wKP878oo6KiFOiIiIjUMBp2LlL3VNlg1eTkZNLT0+natStWqxWr1crWrVt58cUXsVqtgZ6cC3td0tPTA6/FxcXhdrvJyMi46DqnTp0qcvzTp08X6S0qKCwsLBDUKLgREREREalZqizQufHGG9m7dy8pKSmBpVu3btxxxx2kpKRwxRVXEBcXx6ZNmwLbuN1utm7dSu/evQHo2rUrNput0DonT55k3759gXWSkpLIzMzk888/D6yza9cuMjMzA+uIiIiIiEjtUmVD1+rXr0/Hjh0LPRcZGUnjxo0Dz0+dOpX58+fTpk0b2rRpw/z584mIiGDs2LEAREdHM3HiRKZPn07jxo1p1KgRjz76KJ06dWLAgAEAXHnlldx0001MmjSJl19+GYDJkyczbNiwUmdcExERERGRmqXKs65dzGOPPYbT6eT+++8nIyODnj17snHjxkITChctWoTVamX06NE4nU5uvPFGVqxYQWhoaGCdN998k4ceeiiQnW3EiBEsXbq00t+PiIiI1E15eXm43e6qboZIjWaz2Qpd419KldfRqSmysrKIjo4mMzNT83VERERqiOpw/na73aSmppKXl1clxxepTRo0aEBcXFypEoxU6x4dERERkZrM7/dz8uRJQkNDSUxMVNFSkSD5/X5ycnJIT08HID4+/pLbKNARERERqSBer5ecnBwSEhKIiIio6uaI1GgOhwMwMizHxMRcchibbiuIiIiIVBCfzweA3W6v4paI1A75Nww8Hs8l11WgIyIiIlLBVLBUpHyY+Swp0BERERERkVpHc3REREREKpnPB5WZhC0kBExk5RWpFRToiIhIteD1QlUVPLBYwBrkGbEs7S7LcaXm8vng+HGozLI6djskJpoLds6dO8esWbNYt24d6enpdO7cmT/84Q90794dgFOnTvH444+zceNGzp49S58+fViyZAlt2rQpcZ8rVqzg17/+dZHnnU4n4eHhgFH/8IknniA7O5uJEyfy/PPPB9Y7evQogwYN4osvvlC5D7kkfb2KiEiV83qNC79SzC2tEDabcRFoNugoa7uDPa7UbHl5RpATGlo5//der3G8vDxzgc4999zDvn37WLlyJQkJCaxatYoBAwawf/9+EhISGDlyJDabjffff5+oqCheeOGFwOuRkZEl7jcqKoqDBw8Wei4/yPnhhx+45557WLFiBVdccQVDhw6lX79+DB06FID77ruPZ599VkGOlIq+WkVEpMr5/UawUFkXfgV5vcaxg+mVKUu7848b7LHVG1TzWa1GsFsZfkr+VmpOp5N3332X999/nz59+gAwe/Zs3nvvPZYtW8b48ePZuXMn+/bto0OHDgC89NJLxMTEsHr1au65554S922xWIiLiyv2tf/+979ER0czZswYAPr378/+/fsZOnQob731Fna7nVtuucXcm5E6S1+RIiJSbVTmhV9BZi8CLxRMuy0WcDqNHqFgqDdIKpLX68Xn8wV6WvI5HA62bdsWCEQKvh4aGordbmfbtm0XDXTOnz9P8+bN8fl8XHvttcydO5fOnTsD0KZNG3JyctizZw/Nmzdn9+7dTJgwgTNnzvDUU0/x8ccfV8C7ldpKWddEREQo3MNS2sXrDf54VitERUFYmPklNDT4niCR0qhfvz5JSUnMnTuX77//Hp/Px6pVq9i1axcnT56kffv2NG/enJkzZ5KRkYHb7ebZZ58lLS2NkydPlrjf9u3bs2LFCj744ANWr15NeHg41113HYcOHQKgYcOGvP7664wfP54ePXowfvx4Bg8ezKOPPsqUKVNITU2lc+fOdOzYkXfeeaeyfh1SQ1n8fn1NlkZWVhbR0dFkZmZqXKiISDnzeCA11biIr+weHa8XsrKMydrBcLuNgKUye1Y8HnC5oGXLqukBq0mq+vydm5tLamoqLVu2DPR+VPbfe7B/L0eOHGHChAl88sknhIaG0qVLF9q2bcuXX37J/v37SU5OZuLEifz73/8mNDSUAQMGEBJi3EP/6KOPSnWMvLw8unTpQp8+fXjxxReLXWfLli3MmDGDrVu30rp1a1avXk1cXBw9evTg0KFDxMTElP5NSY1X3GeqJOrwFhGROi2/Z6UsmdPye3jMCgkxLnZFqqNWrVqxdetWsrOzycrKIj4+njFjxtCyZUsAunbtSkpKCpmZmbjdbi677DJ69uxJt27dSn2MkJAQunfvHujRuZDL5eL+++9n1apVHD58GK/XS9++fQFo27Ytu3btYvjw4WV/s1IrKdAREZE6L9iaJm43HDhgzLUJhsMBnTsr2JHqLTIyksjISDIyMtiwYQMLFy4s9Hp0dDQAhw4d4osvvmDu3Lml3rff7yclJYVOnToV+/rcuXMZMmQIXbp0Yc+ePXgLjBf1eDz4yjrBTmo1BToiIlJugq0pU5a5LvlcrqoJVpxOaNDA/NA1j8fYtjKLRkr1Uh5/9xV5nA0bNuD3+2nXrh2HDx9mxowZtGvXLlAH5+233+ayyy6jWbNm7N27l4cffpiRI0cyaNCgwD7Gjx9P06ZNWbBgAQBz5syhV69etGnThqysLF588UVSUlL44x//WOT4X331FWvXriUlJQUw5veEhISwfPly4uLi+PrrrwM1fUSKo0BHRETKRVlrypw/bwRJwVyUVVWwAlC/ftnm9+Tmmt+urIkQpGqFhBh/M2532TP+lZbdbhzXjMzMTGbOnMmJEydo1KgRt956K/PmzcP200SfkydPMm3aNE6dOkV8fDzjx49n1qxZhfZx7NixwLwdgLNnzzJ58mTS0tKIjo6mc+fOfPLJJ/To0aPQdn6/n8mTJ7No0aJATR6Hw8GKFSt44IEHcLlcLF26lKZNmwbx25C6QskISqmqJzOKiFR3+ROsg6kp43LBv/9dtkrxJ09CeHhwwYrNZgQ6lcnthjNngjuu12u8z5tvhnr1yr1ptUpVn79Lmjgd7HDJYIWEmCsWKlJdKRlBDXOpoR4qCiciNUkwNWW8XkhLC/777tw5+OCDYHt0/ERF5NCqNdiCOLY7LwKwBHNgwsLgF78wAjQznE7Izq683gApf6GhCjxEKpoun6tYaYZ6qCiciNR2p07BW28ZPTtmh9cY/OS5c2jc2GyQ5ecPw35BmyYpwRwUgG9OX8t9f/0Uv99csOPzgccfwTXXWEwHhm638bsSEZGS6dK5ivn9RpBT0lCP/AJ2GmAoIrWZK9cIVKIiICLC7NZlD1bKou1lKfzz3vpBbbvnxHX88fNPcTjMBUm5ucZ5YdQo+CnhlYiIXECBTjXgdpfcW5M/4TTYSaca9iYilcntNuYdmPrO8vu59vT17Jq+vcLadSn7vruWW5aZ65WxWPz89b5f0LFpStDH7Xz5ZzhSc7BYIk1t5/MZyRuUkEBEpGS6BK5iLhfs2lVysOPzGcM4cnODq54cFgatWyvYEZGKl59QIH+ifGmF+rP5VUjZg5wzaeFsfqMFwXSAH3e1JAQLfouZnhULI5duwxLi54orLKYyr4Vbs3n7f2IBOHokm1yTmepcLjjnjCArK7i5QSIidYEuf6tYVhbs2GHcAS3uJOnz/TzZ1GxBOZ/P2Obyy5WVR0Qqns9nZD6zWs1Nrrd6nNDY+PndRW3xeoKapIPPYyHYpAAJIQd5ytHj0isW47+ezvzp4EpTx3YUuHH14d2xQR33y+PX4fV8auq4IiJ1iQKdKub1GsFOvXrFBzoej3HxEBlpLGYoK4+IVKaTJ2HtWuN7zUw2qXqhcM/Dxs/zMjbhCmlSMQ0shgU/vw4dT3zI10Hv4wrbHqKsZ/CFOExt98XRJLq12BH0cbskfsZ//DmAyZODiEgdoUCnimVlwaFDxrC04nps8oOUXr3Mpx9VVh4RqUy5ucYNFrM3ZuoVCIryQh14/KazEQDG92UwiVtedL1DiM/JFVeA3cQQYavfyf+c6gPAk5F9TB/3mzf9/JMx/L+c/4eZXpkIezb/ftroBUpLg+aZpg+N1Wr+5pmISE2jQKeK5WdVcziKzzTkdEJmJnz+ubGOGcrKIyKVyu/HYcshOhLM1GWMDMkO/Oz1gjeIYCUvz/jOC24+ooUQWwRWB1hNzLPB7+CUvTOx7j3BHBSwkMheQiwheCh9cFcw/8Cnn8Kx780fuX59GDZMwY6I1G4KdKoJq7X4oWtgnPgjgki5mpdnBEku18Xr9JREGdtEpNT8flofvZ7dM8qWVCAvD9xBfF+BcTOoRYvgvrdCQkr+Di6RxcJHjVdy+Gsndru549pxMsNm9AJFWJ24TRzWYc0J/Fy/PjRoYGJjjIDw3DllbKtqbnfl/h9c7DpDpLbSZWwNEBpqDFszO3TN5TKGcnz3XXCBjgqVikhp+b3Z1HOWLchJP+6gcWwETYMrSRNcsFJWFgseIrBAUNnegEDAU1qhtrzAzw6HP6hkM7m55reR8uN2GyM1zp+vvGPWqwc9epj7jJw7d45Zs2axbt060tPT6dy5M3/4wx/o3r07AKdOneLxxx9n48aNnD17lj59+rBkyRLatGlz0f2+++67zJo1iyNHjtCqVSvmzZvHqFGjAq+/+eabPPHEE2RnZzNx4kSef/75wGtHjx5l0KBBfPHFF0SZ6TqWOqlKL2GXLVvGsmXLOHr0KAAdOnTgqaeeYsiQIQDcfffdvP7664W26dmzJzt37gw8drlcPProo6xevRqn08mNN97ISy+9xOWXXx5YJyMjg4ceeogPPvgAgBEjRrBkyRIamL0NVoWCOSkVnJ+Tl1fyesXxeo0vYhUqFal7vF7zn33vOSf5o2uDzZx22NmF8MQI0zd1aiIPDo7ldaZZSLDD3n7aT7aT7GxzkU52tuZvVjWv1why7HbzGVWD4XL9XHfJTKBzzz33sG/fPlauXElCQgKrVq1iwIAB7N+/n4SEBEaOHInNZuP9998nKiqKF154IfB6ZAnjInfs2MGYMWOYO3cuo0aNYt26dYwePZpt27bRs2dPfvjhB+655x5WrFjBFVdcwdChQ+nXrx9Dhw4F4L777uPZZ59VkCOlUqWBzuWXX86zzz5L69atAXj99de5+eab2bNnDx06dADgpptu4rXXXgtsY7/gEzp16lQ+/PBD1qxZQ+PGjZk+fTrDhg0jOTmZ0J/S/owdO5YTJ06wfv16ACZPnsy4ceP48MMPK+NtlonFYgQ5e4I4F3o8xpdbcjI0MZnEKL8ORosWwdXvEZGayeuF48fN9wLn5UD7n36e++MmPFZzXzpeL2S7HVxpqo5N9WE+u6WFP3tXYrc4TWWoA4jiB8bQHIA9KbD/qLntXS7j3KL5m1UvLMz8aI1guc2MjwScTifvvvsu77//Pn36GL2Os2fP5r333mPZsmWMHz+enTt3sm/fvsA120svvURMTAyrV6/mnnvuKXa/ixcvZuDAgcycOROAmTNnsnXrVhYvXszq1av573//S3R0NGPGjAGgf//+7N+/n6FDh/LWW29ht9u55ZZbgvwtSF1TpYHO8OHDCz2eN28ey5YtY+fOnYEPTVhYGHFxccVun5mZyfLly1m5ciUDBgwAYNWqVSQmJrJ582YGDx7MgQMHWL9+PTt37qRnz54AvPLKKyQlJXHw4EHatWtXge+w7Gw2aNw4uJ6VnBzjiy0sTKmpRaR0/H4jyAkNNTdsNa/AXAO334HXxOR6KDzBviYJCTG+p/NLAZhjIccbQWSkuXTcHn7OTHP8OLhNJiNwuYxe/jNnICHB3LZSd3i9Xnw+H+EXRGIOh4Nt27YFApGCr4eGhmK329m2bVuJgc6OHTt45JFHCj03ePBgFi9eDECbNm3Iyclhz549NG/enN27dzNhwgTOnDnDU089xccff1yO71Jqu2oz+8Ln8/H222+TnZ1NUlJS4PktW7YQExNDgwYN6Nu3L/PmzSMmJgaA5ORkPB4PgwYNCqyfkJBAx44d2b59O4MHD2bHjh1ER0cHghyAXr16ER0dzfbt20sMdFwuF64CfftZWVnl/ZZLLdgeFY/HOAkHc8dIk1RF6i632/zcvLw6+p1ht0PbtuaHB4Pxez5ypGxDhOvVA08QZ/KsLH3Py8XVr1+fpKQk5s6dy5VXXklsbCyrV69m165dtGnThvbt29O8eXNmzpzJyy+/TGRkJC+88AJpaWmcPHmyxP2mpaURG1u4SG5sbCxpaWkANGzYkNdff53x48fjdDoZP348gwcPZsKECUyZMoXU1FRGjBiBx+Nh9uzZ3HbbbRX6e5CarcoDnb1795KUlERubi716tVj3bp1XHXVVQAMGTKE22+/nebNm5OamsqsWbO44YYbSE5OJiwsjLS0NOx2Ow0bNiy0z4IfmLS0tEBgVFBMTExgneIsWLCAOXPmlOM7FRGp/lwu2LXLfLBj9UCbBhXWrGqtKjNZWa1gMXl8JZiR0lq5ciUTJkygadOmhIaG0qVLF8aOHcuXX36JzWbj3XffZeLEiTRq1IjQ0FAGDBgQmGd9MZYLhqj6/f5Cz40aNapQcoItW7awd+9eli5dSuvWrVm9ejVxcXH06NGDPn36FHudJwLVINBp164dKSkpnD17lnfffZe77rqLrVu3ctVVVwW6RQE6duxIt27daN68OX//+98vOj7zwg/MhR+o4ta50MyZM5k2bVrgcVZWFomJiWbfnohIjZKVBTt2GD0UZi7gw/1AX+PnGjrNpkbyeMBtsjepstMaS83VqlUrtm7dSnZ2NllZWcTHxzNmzBhatmwJQNeuXUlJSSEzMxO3281ll11Gz5496datW4n7jIuLK3KjOT09vUgvTz6Xy8X999/PqlWrOHz4MF6vl759jS+btm3bsmvXriJTIUTyVXmgY7fbA8kIunXrxu7du/nDH/7Ayy+/XGTd+Ph4mjdvzqFDhwDjw+J2u8nIyCjUq5Oenk7v3r0D65w6darIvk6fPl3ihwqMuUFhlZEKRUSkGvF6jWCnXj1zgY69wPCrEPMJ1yRIqamQZTJxhMdj/D+fPVshTZJaKDIyksjISDIyMtiwYQMLFy4s9Hr0T1ktDh06xBdffMHcuXNL3FdSUhKbNm0qNE9n48aNgeu2C82dO5chQ4bQpUsX9uzZg7dAlO7xePBpMrFcRJUHOhfy+/2F5sYU9OOPP3L8+HHi4+MB406CzWZj06ZNjB49GoCTJ0+yb9++wIcwKSmJzMxMPv/8c3r06AHArl27yMzMLPFDVdvk5hqJBcwoa/rRYNLT5lOhUpGqk5UFhw4ZcwPN3OupFwr8NL0yLw+8JnsZdK0SHLsd7CZ70Hw+YwmmvpqUr8pK8x3scTZs2IDf76ddu3YcPnyYGTNm0K5dO379618D8Pbbb3PZZZfRrFkz9u7dy8MPP8zIkSMLzZ0eP348TZs2ZcGCBQA8/PDD9OnTh+eee46bb76Z999/n82bN7Nt27Yix//qq69Yu3YtKSkpALRv356QkBCWL19OXFwcX3/9daCmj0hxqvRy8re//S1DhgwhMTGRc+fOsWbNGrZs2cL69es5f/48s2fP5tZbbyU+Pp6jR4/y29/+liZNmgTGbUZHRzNx4kSmT59O48aNadSoEY8++iidOnUKZGG78soruemmm5g0aVKgl2jy5MkMGzas2mRcy69ZU1KtnGCL4FksxoksOdmoGG5Gbq4RqAwfDhHmkifh88HJk8GfRFWoVKTqeL3GZ9fhMPfZdxToxfH7wR1kkWL1BpljtYLV5E0ls+mspfxZrUav6fnz5tM+B6tePfPn1czMTGbOnMmJEydo1KgRt956K/PmzcP2U5akkydPMm3aNE6dOkV8fDzjx49n1qxZhfZx7NgxQgp8sHv37s2aNWt48sknmTVrFq1atWLt2rWFkkaBceN78uTJLFq0KFCTx+FwsGLFCh544AFcLhdLly6ladOmQfw2pK6o0kvJU6dOMW7cOE6ePEl0dDRXX30169evZ+DAgTidTvbu3csbb7zB2bNniY+Pp3///qxdu5b69X8um71o0SKsViujR48OFAxdsWJFoIYOGBV2H3roocAdhhEjRrB06dJKf7/FOXsWTp+GjIySs6vZbEZWH7PBjtVqpJWOiDAfrPj98OOPcPiw0TazfD5o0MD8CTX/IkuFSkWqjtdr9MqYySRWcNX4eGhqMqU9BH9Tp6Yz25vl1RyoGs9uhx49KneulNVq/vM1evTowIiZ4jz00EM89NBDF93Hli1bijx32223XTJbmsVi4bPPPivy/LBhwxg2bNhFtxXJV6WBzvLly0t8zeFwsGHDhkvuIzw8nCVLlrBkyZIS12nUqBGrVq0Kqo0VLb/2gt1e/BdQ/vCCYFKXghFohIcHl17aYoGDB80HSWBs061bcKmxNYRFpOqczfCTnZmDKxt+NHGGiLb/PD42PBxCKqkIYk0WbA0ed4EbQfVCfzR9XKsNztOYvLyQoHreNby4fJR03heR8qOvqmriYsX5quLC32o1KmZHRZkvNurxGF3xwQZnIlJ2Qc2T8/u55vT1HP+/7WU6ts0Gul9xacHW4IksMLxvyxNXBnXsY99expa0NFJTzY8V1PBiEakp9DVVB5Q09+dS2/h8wfUGQeWNORaRorxeOH7c/Dw5iy+HNt6yBTnpxx24vBFYdae6VIK5o+/zN+a7E5fR9PLTQR+3WfPTWHN+JC/vMlPb5c8p1fBiEakJFOjUYhaLEbDs2WN+2/xepOuvN9+jIyJVy+83gpyL9RQXx1KgG+aaOafAWvoPv40cZtj64PNYIK7mTSIpy9xAiyW4YbpBs4Rw98Y0vjvyI3a7uf/jKPuP/GuG0Qt09Fvwm0xU4/Uax2vRopLfs4hIEBTo1GI2GzRuHNzJOycHnE7NlxGpyaxWcxejlgKjmDxuCxZ/6QOWEELwUTNTpnk8RvKVYHqvwbih1Lhx5Qc7WZ7LsFvMZ13L5wg3fyPL6TTKD+jcICI1gQKdWi7YE6/HY5zQRKTu8Be4K/J0dB98npoVuATbK+PxGEFO587mg538XvOaOJTLag8uUY2ISE2hQEdERADIyy373Y1Ub2e8Fkeln1zK2isTHm5kiwx2exERqX4U6NQQwU7uL2tditxcY5iC2W2CSYAgItXHgqwNeKxNTG3j9UK228GVlsqfo+P3B98rA8Z3ZV0Lcjwu89/v2dngclVMe0REypsCnWrOYjEuHo4cCW77YIuNWizGHdLkZKNCuhkej3HR0KOH+W1FpHpw+x14MVdEqzqMasrvmZFL278fUk+Z2yY31wgqR40yShCIiFRnNWsAdh0UGmpMFs0vLGZmCQ0Nvtio1WocNyLC+NfMYrcrkYGISHUX5gju+z03V3N1yoPTCZmZlbcEM+/23LlzTJ06lebNm+NwOOjduze7d+8OvH7+/HkefPBBLr/8chwOB1deeSXLli275H7fffddrrrqKsLCwrjqqqtYt25dodfffPNNEhMTadSoETNmzCj02tGjR2nbti1ZWVnm35DUOerRqUJ+vx+/y4ndAjaMpVihDiC4oSBlCTZCQ4OroxNMpW0REalcYUEkI9CwtfLhdML770NGRuUds2FDuPlmcyMt7rnnHvbt28fKlStJSEhg1apVDBgwgP3799O0aVMeeeQRPv74Y1atWkWLFi3YuHEj999/PwkJCdx8883F7nPHjh2MGTOGuXPnMmrUKNatW8fo0aPZtm0bPXv25IcffuCee+5hxYoVXHHFFQwdOpR+/foxdOhQAO677z6effZZoqKiyuPXIrWcAp0q5HM6+eHx7ixOuPh6x/I686pvJcEGO1UlN9f8HSSPR3cKRcqD22305pr5POVpbp1IpXC7jSDH4aicuWG5ucbx3O7SBzpOp5N3332X999/nz59+gAwe/Zs3nvvPZYtW8YzzzzDjh07uOuuu+jXrx8AkydP5uWXX+aLL74oMdBZvHgxAwcOZObMmQDMnDmTrVu3snjxYlavXs1///tfoqOjGTNmDAD9+/dn//79DB06lLfeegu73c4tt9xStl+I1BkKdGqAZiF7sPmceEyOly8PwSQVyK+zsHev+bHyKkYnUnYuF/z73z9/nkor1A1twiquXVL+zPbaewvcL3O5jZppZjidwSfHkaLCg6hlFCyzNx69Xi8+n4/wCyIxh8PBtm3bALj++uv54IMPmDBhAgkJCWzZsoVvvvmGP/zhDyXud8eOHTzyyCOFnhs8eDCLFy8GoE2bNuTk5LBnzx6aN2/O7t27mTBhAmfOnOGpp57i448/NvdGpE5ToFOFQh0Omjy3mzvH/TwfpiA7TmbY+lRJ2yyWn+tDmOXz/TwvSMXoRCqfz2d8liIjzQ1TCbECP332qiBxmpgQEmLcDPJ4zH1fugvU+9m3F/JMBrYul3FMTY+o/erXr09SUhJz587lyiuvJDY2ltWrV7Nr1y7atGkDwIsvvsikSZO4/PLLsVqthISE8P/+3//j+uuvL3G/aWlpxMbGFnouNjaWtLQ0ABo2bMjrr7/O+PHjcTqdjB8/nsGDBzNhwgSmTJlCamoqI0aMwOPxMHv2bG677baK+yVIjadApwpZgFC7H6uNwFJQKHmE2vLweSr/isNmMyp9B1MELyfHuONntaoYnUhVcbmMz7GZjIt+3WCoMex2I6Om2WQzkQVSENntEGKy193nM77j9V1dN6xcuZIJEybQtGlTQkND6dKlC2PHjuXLL78EjEBn586dfPDBBzRv3pxPPvmE+++/n/j4eAYMGFDifi0X3Enx+/2Fnhs1ahSjRo0KPN6yZQt79+5l6dKltG7dmtWrVxMXF0ePHj3o06cPMTEx5fzOpbZQoFOVfDnc+EM9vl988dXSjzuYt7zyy24HO3TM4wkuu4uIlI/sbKM31mIxd7Mh3A+jehs/m53fA+qJrWzB1EgLK3B9abNBqMl9mBkKKTVfq1at2Lp1K9nZ2WRlZREfH8+YMWNo2bIlTqeT3/72t6xbty6QKODqq68mJSWF//3f/y0x0ImLiwv03uRLT08v0suTz+Vycf/997Nq1SoOHz6M1+ulb9++ALRt25Zdu3YxfPjwcnzXUpsovXQNEJPoxGEzOZBaROosr9cYemq3m0sd7Ljg7r7bbW7x+YyL5xCdWURqlcjISOLj48nIyGDDhg3cfPPNeDwePB4PIRd84ENDQ8m7SFdjUlISmzZtKvTcxo0b6d27d7Hrz507lyFDhtClSxd8Ph/eAndgPB4PPt1hkYvQvZmqFBrBP5ucZ9ydxc/RibKf5tMZLaumbSJS49lNpg+2F7heiImFy4PI3hoSElxPg4hUPxs2bMDv99OuXTsOHz7MjBkzaNeuHb/+9a+x2Wz07duXGTNm4HA4aN68OVu3buWNN97ghRdeCOxj/PjxNG3alAULFgDw8MMP06dPH5577jluvvlm3n//fTZv3hxIcFDQV199xdq1a0lJSQGgffv2hISEsHz5cuLi4vj666/p3r17pfwupGZSoFOVLBbyLJHkuCHEDpYL6s/YUC+OiATH7TaGkJoZguorcBPWbquctLcidVkwmU0r8ziZmZnMnDmTEydO0KhRI2699VbmzZuH7acvljVr1jBz5kzuuOMOzpw5Q/PmzZk3bx6/+c1vAvs4duxYoV6f3r17s2bNGp588klmzZpFq1atWLt2LT179ix0bL/fz+TJk1m0aBGRP2U2cjgcrFixggceeACXy8XSpUtp2rRpcG9O6gQFOiIitUxWFhw6ZAQ5YSayatULBboYPyvrWu3ncv2Uac8Et1uJCMqD3W4U8MzIqLw5rQ0bmu9tHT16NKNHjy7x9bi4OF577bWL7mPLli1FnrvtttsumS3NYrHw2WefFXl+2LBhDBs27KLbiuRToCMiUst4vUZSEIfDXC0rR4Gh9pp0XvudPAlZnkuvV1B+UeezZyukSXWGwwE331y5NYnsdnPp5kVqA53KaggbTtND2SyAG8dPP4lIXWO1mruDa9NXRZ0SzHwqn89YPCYDJCnK4VDgIVLRFOjUEFNtg/GZTZJng1RbZ7b4V6JgR0RECrJawWqyckFoaMW0RUSkIigJaDXmoey3elpa92D1q6iNiIiIiNQt6tGp1n7uhXne8wlOT2Spt7TjZIatT0U0qtRyc43ChWZkZxsTZEVEREREykKBTg3hIQIPJmYVVyGLxRi/nZxsfvxxbi74/TBqFERHV0z7ROoCr9eY6GwmrWxoFffxezzG5z/YbatSsMe3WMylABcRkdJToCPlzmr9uQCqmYxPAHl5kJmp9KUiZXH2LJw+baSuNXMRHVVg3cpOL+3xwI8/lq12T3i4McG+MoWEGMfNzQ0u2MnNhcaNFeyIiFQEBTpSIUJDjZO/2YsWDVsTKTuPx8iMZbeby6plL3CxbbOBr/ybViK/3/i+6Nw5+GAnP+ioTOHh0LWrcZPGrNxc2LMn+F4sERG5OAU6UmGCqcTsdFZuXQGR2iw01Fw9nOpQOyc83HxPcFWr7OBKRERKp0pHZC9btoyrr76aqKgooqKiSEpK4h//+Efgdb/fz+zZs0lISMDhcNCvXz+++uqrQvtwuVxMmTKFJk2aEBkZyYgRIzhx4kShdTIyMhg3bhzR0dFER0czbtw4ztahamdujxF0mF2CDTgslp/vVO7YYW758kujontWVvn+DkREpDD7T/XZzC6gLigRqRmq9P7d5ZdfzrPPPkvr1q0BeP3117n55pvZs2cPHTp0YOHChbzwwgusWLGCtm3b8swzzzBw4EAOHjxI/fr1AZg6dSoffvgha9asoXHjxkyfPp1hw4aRnJxM6E8J/8eOHcuJEydYv349AJMnT2bcuHF8+OGHVfPGK9l//wvBzJO12aBtW/MF5Ww2Y8x5MMMxfD7IydEcHRGRijbdMRifx+T9zmg4HN4Zv+qzldn588GNfAhWeDjUq2dum3PnzjFr1izWrVtHeno6nTt35g9/+APdu3cH4O677+b1118vtE3Pnj3ZuXPnRff77rvvMmvWLI4cOUKrVq2YN28eo0aNCrz+5ptv8sQTT5Cdnc3EiRN5/vnnA68dPXqUQYMG8cUXXxAVFWXuDUmdU6WBzvDhwws9njdvHsuWLWPnzp1cddVVLF68mN/97nfccsstgBEIxcbG8tZbb3HvvfeSmZnJ8uXLWblyJQMGDABg1apVJCYmsnnzZgYPHsyBAwdYv349O3fupGfPngC88sorJCUlcfDgQdq1a1e5b7oK2O3mT0f5la+DGXcOwU+srQ5DZ0RqPL+fCHsODhs4THwWHXaT+eCLEWzmtKrOmlZX+Cxlr8/WOmwPuJ1QQzKBVkfnz8Of/2wk4KgsjRvD5Mnmgp177rmHffv2sXLlShISEli1ahUDBgxg//79NG3aFICbbrqJ1157LbCN/RJ3R3fs2MGYMWOYO3cuo0aNYt26dYwePZpt27bRs2dPfvjhB+655x5WrFjBFVdcwdChQ+nXrx9Dhw4F4L777uPZZ59VkCOlUm0uK30+H2+//TbZ2dkkJSWRmppKWloagwYNCqwTFhZG37592b59O/feey/Jycl4PJ5C6yQkJNCxY0e2b9/O4MGD2bFjB9HR0YEgB6BXr15ER0ezffv2EgMdl8uFq8DM+KwqHktl9gLERg6htjx8HgtWa3ADDXyVORNZRMqH30/3s9fz/eLtlX7osmZOq4qsaXVOgXR6C5yf4PHXrPpstUVurvFZcTgqZ05aTo5xvNzc0gc6TqeTd999l/fff58+fYz/99mzZ/Pee++xbNkynnnmGcC4NouLiyt1WxYvXszAgQOZOXMmADNnzmTr1q0sXryY1atX89///pfo6GjGjBkDQP/+/dm/fz9Dhw7lrbfewm63B26Ai1xKlQc6e/fuJSkpidzcXOrVq8e6deu46qqr2L7dOEnHxsYWWj82NpZvv/0WgLS0NOx2Ow0bNiyyTlpaWmCdmJiYIseNiYkJrFOcBQsWMGfOnDK9NzN8vqLDtbwFumF2zij8eyit9OMO5i3XeGqROsOXQwNv2YKc9OMOXN4IrCaHrZY1c1pVZE2ry2x2P/680p8fQskL3ECT8hERYX44WbCcTnPre71efD4f4Rd8KB0OB9u2bQs83rJlCzExMTRo0IC+ffsyb968Yq+78u3YsYNHHnmk0HODBw9m8eLFALRp04acnBz27NlD8+bN2b17NxMmTODMmTM89dRTfPzxx+beiNRpVR7otGvXjpSUFM6ePcu7777LXXfdxdatWwOvWy4o5uD3+4s8d6EL1ylu/UvtZ+bMmUybNi3wOCsri8TExEu+H7NsNiMzUl5e0cn/bncEn6deR4+WnwW9/5hEJw5bDh5P6e/aiUjN5ff7A0NV313UFq/ZORhgXMjGBX8xWxMzp9VFu2cGfwMtRTmxa7369euTlJTE3LlzufLKK4mNjWX16tXs2rWLNm3aADBkyBBuv/12mjdvTmpqKrNmzeKGG24gOTmZsLCwYveblpZW7E3s/JvPDRs25PXXX2f8+PE4nU7Gjx/P4MGDmTBhAlOmTCE1NZURI0bg8XiYPXs2t912W8X+IqRGq/JAx263B5IRdOvWjd27d/OHP/yBxx9/HDA+EPHx8YH109PTAx+QuLg43G43GRkZhXp10tPT6d27d2CdU6dOFTnu6dOni3zQCgoLCyvxQ1qeGjSAyy4z/v0pv0IBFmZs/JTvj+Vgt5ubvxJlP82nM1qWX0NFpEbwOZ2BL3avJ8T8ZHMg1dsZr8VR9ScIKXcuXwT7Tl1Hx9iy3UALJQfQDbTabuXKlUyYMIGmTZsSGhpKly5dGDt2LF9++SVAYHgZQMeOHenWrRvNmzfn73//+0WHl13qJvaoUaMKJSfYsmULe/fuZenSpbRu3ZrVq1cTFxdHjx496NOnz0V7kKRuq3bnMb/fj8vlomXLlsTFxbFp0yY6d+4MgNvtZuvWrTz33HMAdO3aFZvNxqZNmxg9ejQAJ0+eZN++fSxcuBCApKQkMjMz+fzzz+nRowcAu3btIjMzMxAMVTWr1UgYUPyQDQtOTyQ+C1hN3EAzUoCKSF22IGsDHmsTU9t4vZDtdnDlJXrOpaay8MhHn3L0sG6gyaW1atWKrVu3kp2dTVZWFvHx8YwZM4aWLYv/O4iPj6d58+YcOnSoxH3GxcUVmTpQ8Cb2hVwuF/fffz+rVq3i8OHDeL1e+vbtC0Dbtm3ZtWtXkeRWIvmqNND57W9/y5AhQ0hMTOTcuXOsWbOGLVu2sH79eiwWC1OnTmX+/Pm0adOGNm3aMH/+fCIiIhg7diwA0dHRTJw4kenTp9O4cWMaNWrEo48+SqdOnQJZ2K688kpuuukmJk2axMsvvwwY6aWHDRtWJzKuiUjd5fY78JrMjKXM7nWBbqCJOZGRkURGRpKRkcGGDRsCN5Mv9OOPP3L8+PFCI3EulJSUxKZNmwrN09m4cWOJN5/nzp3LkCFD6NKlC3v27MFbYEKzx+PBp8xJchFVGuicOnWKcePGcfLkSaKjo7n66qtZv349AwcOBOCxxx7D6XRy//33k5GRQc+ePdm4cWOghg7AokWLsFqtjB49GqfTyY033siKFSsCNXTAyMf+0EMPBbKzjRgxgqVLl1bumxURERGpQTZs2IDf76ddu3YcPnyYGTNm0K5dO379619z/vx5Zs+eza233kp8fDxHjx7lt7/9LU2aNCk07Gz8+PE0bdqUBQsWAPDwww/Tp08fnnvuOW6++Wbef/99Nm/eXCjBQb6vvvqKtWvXkpKSAkD79u0JCQlh+fLlxMXF8fXXXwdq+ogUp0oDneXLl1/0dYvFwuzZs5k9e3aJ64SHh7NkyRKWLFlS4jqNGjVi1apVwTazTrswQUJphYSYLzQqIiJSl+RUUidZsMfJzMxk5syZnDhxgkaNGnHrrbcyb948bDYbXq+XvXv38sYbb3D27Fni4+Pp378/a9euLXRD+tixY4QUyBvfu3dv1qxZw5NPPsmsWbNo1aoVa9euLVQGBIypDJMnT2bRokVERhrzwRwOBytWrOCBBx7A5XKxdOnSQD0fkeJUuzk6Uj1YLMZY/SNHgtveZoO2bRXsiIiIXCg83Cjg+eOP5tM+B6txY/Pp20ePHh2YA30hh8PBhg0bLrmPLVu2FHnutttuu2S2NIvFwmefFU2aMWzYMIYNG3bJ44qAAh0pQWgoREYGV+Hc5zMKB+bllX+7REREarp69WDyZKOAZ2UJD6+8mj0i1YUCHSlRgWlOpmluoIiISMnq1VPgIVLRFOjUATacQWXL8eAAlGJWRGo/tzv4XmjNSRQRqZ4U6NQBU22D8WG+aOCxvM686luJgh0Rqc3cbjh1CiLMZeIOyMmB2FgFOyIi1Y0CnVrK6I0pm2Yhe7D5nHhM1uEQEalJ8vKMIKdnT3CY/Op0OmHXrrLNSfR4gtvOYjESv1SFvLzg2m2xmCtSKiJSFvq6qbV+7oV53vMJTk9kqbe042SGrU9FNEpEpNpyOIwkLJUlJMSYIJ6bG1zQkJtrZNKqimDn9GlITTW/nc0GiYkKdkSkcuirpg7wEKFeGRGp9YKdZxNsvbCyCg+Hrl2Da3NuLuzZE1xmzPJgtUNYmLltvF4joKuqNotI3aNAR0REaryyzrOJiChbpslgma1rUl1YrcH1JCkjp4hUJgU6IiJS45Vlng0YQU5NDTqCZTbo8BbIS+Nxmy906fEYvToiIpVFgY5US/lFR83SRFeRn+Xlmb+wLOsdd5/PGFYVYj7RY7mkaa7seTY1UUiI0Rvj8Zj7/3YXGHK2dy94TfboeL3G93OLFlWXREFE6hZdEkq14/fDyZNQv775bTXRVeRneXngDiJwsdmCC1S8XsjONtItBzMPQ2maK4fdDm3bmp8bFFngb8JqMx9QOp3G34eGrxkyM42/+coSEQHR0aVff/bs2cyZM6fQc7GxsaSlpQHg9/uZM2cOf/7zn8nIyKBnz5788Y9/pEOHDhfd77vvvsusWbM4cuQIrVq1Yt68eYwaNSrw+ptvvskTTzxBdnY2EydO5Pnnnw+8dvToUQYNGsQXX3xBVFRU6d+M1Fm6HKwhyjLEoCaxWIwhJGFhmugqUlZNmkBkE/PbBduz4vcbQVLXrhATY27b8kjTLKUXzP9vWIHzit1mfqifhq39LDMT5s6FH36ovGM2aQKzZpkLdjp06MDmzZsDj0MLTGRbuHAhL7zwAitWrKBt27Y888wzDBw4kIMHD1K/hDuVO3bsYMyYMcydO5dRo0axbt06Ro8ezbZt2+jZsyc//PAD99xzDytWrOCKK65g6NCh9OvXj6FDhwJw33338eyzzyrIkVJToFPNBTvEILTAxX5NvHumia4iZWcL4mK0PISHa/iYyMXk5BhBjsMRfAKNYI6Xk2Mu0LFarcTFxRV53u/3s3jxYn73u99xyy23APD6668TGxvLW2+9xb333lvs/hYvXszAgQOZOXMmADNnzmTr1q0sXryY1atX89///pfo6GjGjBkDQP/+/dm/fz9Dhw7lrbfewm63B44nUhoKdKq5YIcYhBf4n1UPh4iISPUTERHcMO1gmE0eAXDo0CESEhIICwujZ8+ezJ8/nyuuuILU1FTS0tIYNGhQYN2wsDD69u3L9u3bSwx0duzYwSOPPFLoucGDB7N48WIA2rRpQ05ODnv27KF58+bs3r2bCRMmcObMGZ566ik+/vhj829C6rQgRmFLZbPbjTukZhazw75ERERE8vXs2ZM33niDDRs28Morr5CWlkbv3r358ccfA/N0YmNjC21TcA5PcdLS0i66TcOGDXn99dcZP348PXr0YPz48QwePJhHH32UKVOmkJqaSufOnenYsSPvvPNOOb9jqY3UoyMXZcf8LSAL4Mbx00/BcbmUulRERKSqDBkyJPBzp06dSEpKolWrVrz++uv06tULAIul8Hne7/cXee5Cl9pm1KhRhZITbNmyhb1797J06VJat27N6tWriYuLo0ePHvTp04cYsxMCpU5RoFMHRNizcZsY+mYjh1BbHj6PhRm2PuYPaINUW2e2+FdiNtixWIxhevv3w48/mjusUpeKiIhUjMjISDp16sShQ4cYOXIkYPTQxMfHB9ZJT08v0mNTUFxcXJEen4tt43K5uP/++1m1ahWHDx/G6/XSt29fANq2bcuuXbsYPnx4Gd+Z1GYKdOqA3TNL/tK5mPTjDja/0YJgemZaWvewze8EzM2yzM/2FBGh1KUiIiLVhcvl4sCBA/ziF7+gZcuWxMXFsWnTJjp37gyA2+1m69atPPfccyXuIykpiU2bNhWap7Nx40Z69+5d7Ppz585lyJAhdOnShT179uAtMGzD4/Hg0wlfLkGBTi3l8kWw79R1dIz9LOh9xCQ6eYGtOD2ljzjsOIPrBSogJMSYY6TUpSIiIlXj0UcfZfjw4TRr1oz09HSeeeYZsrKyuOuuu7BYLEydOpX58+fTpk0b2rRpw/z584mIiGDs2LGBfYwfP56mTZuyYMECAB5++GH69OnDc889x80338z777/P5s2b2bZtW5Hjf/XVV6xdu5aUlBQA2rdvT0hICMuXLycuLo6vv/6a7t27V8rvQmouBTq1loVHPvqUo4dzsNvNFdB02LPZOcPoBfIQgcdkr0w+twe8uea2cbkUsIjUZW538Nt5PHD+fOVmmszONo4bbLuDrVkktUNlFQwN5jgnTpzgf/7nf/jhhx+47LLL6NWrFzt37qR58+YAPPbYYzidTu6///5AwdCNGzcWqqFz7NgxQgpUH+7duzdr1qzhySefZNasWbRq1Yq1a9fSs2fPQsf2+/1MnjyZRYsWEfnT8A6Hw8GKFSt44IEHcLlcLF26lKZNmwbx25C6RIFOrWbB6YnEZwFrFaSY/u9/wWNyG6/XmKNz/nyFNElEqqnQUGPIak6OETiY5fEYQcf588FtHyyXyzhufhvMysmB2FgFO3VNRIRRwPOHH4JL+xyMJk3M1exZs2bNRV+3WCzMnj2b2bNnl7jOli1bijx32223cdttt11y3599VnREyrBhwxg2bNhFtxUpSIGOVBi73fzsHr9f2dNE6qLwcEhKCn6O3fnzxjJyJNSrV65Nu+RxwTim2eM6nbBrl/k6aVLzRUfDrFmV16MDRpBjplioSG2gQEcqjNUKZjuSKvNOrIhUL2bn5RWUf5OkXr3Kv5gLCzOSp5hNoCJ1W3S0Ag+RiqZAR0REqgW3O/jejbp6kyTY922xKA2/iNR+CnRERKopr9f8xHqvt2Z+sbvdkJ4ODkdw23u9xrahoeXbruoqJMToAcvNDS7Yyc2Fxo0V7IhI7RbU+fD48eMcPXqUnJwcLrvsMjp06EBYWFh5t01EpM7yeuH4cfMXsXk50L5imlSh8vKMQOXaa4MbwpY/t6+unIrCw6Fr1+B6wHJzYc+eys1OJyJSFUIuvYrh22+/ZebMmbRo0YIWLVrQt29fhgwZQrdu3YiOjmbgwIG8/fbb5Jn41l2wYAHdu3enfv36xMTEMHLkSA4ePFhonbvvvhuLxVJo6dWrV6F1XC4XU6ZMoUmTJkRGRjJixAhOnDhRaJ2MjAzGjRtHdHQ00dHRjBs3jrNnz5a6rSIilSl/zkloqHHxXuqlBmfv8vmM92u1ml9strqXuSw83JhgbnYpy1yofC63kW3O7OJylf3YIiKlVaoenYcffpjXXnuNQYMG8fvf/54ePXrQtGlTHA4HZ86cYd++fXz66afMmjWLOXPm8Nprr5WqiNPWrVt54IEH6N69O16vl9/97ncMGjSI/fv3B/KmA9x000289tprgcf2C85mU6dO5cMPP2TNmjU0btyY6dOnM2zYMJKTkwn9aRzD2LFjOXHiBOvXrwdg8uTJjBs3jg8//LA0v4I6y2HPNrW+jRxCbXn4PGbzrYlIcfIv4ksrryaOW8MIcnJyjAvhkFLfgivMZjPmnkjFO/QNnDGZMSw31wjgR43SJHwRqRylOiXa7XaOHDnCZZddVuS1mJgYbrjhBm644QaefvppPvroI7799ttSBTr5QUe+1157jZiYGJKTk+nTp0/g+bCwMOLi4ordR2ZmJsuXL2flypUMGDAAgFWrVpGYmMjmzZsZPHgwBw4cYP369ezcuTNQlOqVV14hKSmJgwcP0q5du9L8Guqk/MKhZqUfdzBvucZFiEjp+P1GUHf55RAVFdw+LBZzxZEleOHh5rPM5eVBZqbKB4hI5SnVfbPnn3++2CCnOL/85S8vWQiqJJmZmQA0atSo0PNbtmwhJiaGtm3bMmnSJNLT0wOvJScn4/F4GDRoUOC5hIQEOnbsyPbt2wHYsWMH0dHRhSrv9urVi+jo6MA68jOnJ4LkY9eVaR8xiU4ctkosECAitUJ+D1Ywi4KcymO3G8GOmaWuDS0UkapX6tPCk08+yQ033EDv3r0JL48Bvhfw+/1MmzaN66+/no4dOwaeHzJkCLfffjvNmzcnNTWVWbNmccMNN5CcnExYWBhpaWnY7XYaNmxYaH+xsbGkpaUBkJaWRkxMTJFjxsTEBNa5kMvlwlVgMHFWVlZ5vM0awsL/vPppUIFKlP00n85oWQFtEhERkZrm6NGjtGzZkj179nDttdea2nb27Nm89957pKSklLjO3XffzdmzZ3nvvffK1M7apEWLFkydOpWpU6dWdVOqXKlHQq9evZoBAwbQoEED+vbty+zZs/nkk09wu93l0pAHH3yQ//znP6xevbrQ82PGjGHo0KF07NiR4cOH849//INvvvmGv//97xfdn9/vx1JgsLalmIHbF65T0IIFCwKJC6Kjo0lMTAziXdVkFpyeSPOLWxXzREREarq7776bkSNHFnl+y5YtWCyWSkno9Oijj/LPf/6zwo9THbRo0YLFixdXdTNqnVIHOkeOHOH48eO88sortG7dmpUrV9KvXz8aNmzIgAEDmDdvXtDDwKZMmcIHH3zAxx9/zOWXX37RdePj42nevDmHDh0CIC4uDrfbTUZGRqH10tPTiY2NDaxz6tSpIvs6ffp0YJ0LzZw5k8zMzMBy/PjxYN6aiEjQ3G5wOk0uuVXdahGRsvH7/Xi9XurVq0fjxo2rujlSg5nKbdO0aVPGjRvH8uXLOXLkCN9++y3Lli2jWbNmLFy4sFACgdLw+/08+OCD/PWvf+Vf//oXLVteesjTjz/+yPHjx4mPjwega9eu2Gw2Nm3aFFjn5MmT7Nu3j969ewOQlJREZmYmn3/+eWCdXbt2kZmZGVjnQmFhYURFRRVaREQqi8sF//437NoFO3eWfklOLrATZSATkQqUnZ1NVFQU77zzTqHnP/zwQyIjIzl37lzgua+//jow/aFDhw5s2bIl8Fp+L9GGDRvo1q0bYWFhfPrpp8yePbvQcDefz8e0adNo0KABjRs35rHHHsNfioJQ27dvp0+fPjgcDhITE3nooYfIzjayyr7xxhvUq1cvcAMdjBvwbdu2DazTokUL5s6dy9ixY6lXrx4JCQksWbKk0DEyMzOZPHkyMTExREVFccMNN/Dvf/+70DoffPAB3bp1Izw8nCZNmnDLLbcA0K9fP7799lseeeSRQCmV0rQdjBv7w4cPx+Fw0LJlS958881L/j7qkiCTeBo9PBs3bmTDhg1s2LABn89H//79Te3jgQceYNWqVbz11lvUr1+ftLQ00tLScDqdAJw/f55HH32UHTt2cPToUbZs2cLw4cNp0qQJo0aNAiA6OpqJEycyffp0/vnPf7Jnzx7uvPNOOnXqFMjCduWVV3LTTTcxadIkdu7cyc6dO5k0aRLDhg1TxrVqyuVSjQap23w+o4fGbjeyW5V2iYj4eR8hCnRqBLfbSL1sdimnkeNSyfx+P96cnCpZShMUmBEZGcmvfvWrQiVAwMiie9ttt1G/fv3AczNmzGD69Ons2bOH3r17M2LECH788cdC2z322GMsWLCAAwcOcPXVVxc53v/93//x6quvsnz5crZt28aZM2dYt27dRdu4d+9eBg8ezC233MJ//vMf1q5dy7Zt23jwwQcBGD9+PL/85S+544478Hq9rF+/npdffpk333yzUKmT559/nquvvpovv/ySmTNn8sgjjwRusvv9foYOHUpaWhofffQRycnJdOnShRtvvJEzZ84A8Pe//51bbrmFoUOHsmfPHv75z3/SrVs3AP76179y+eWX8/vf/56TJ09y8uTJUrUdjCGGR48e5V//+hfvvPMOL730UqGkXXVdqZMRpKam8vHHH/Pxxx+zZcsWMjMzue666+jbty8PPvgg3bt3x2oy5c2yZcsAI5It6LXXXuPuu+8mNDSUvXv38sYbb3D27Fni4+Pp378/a9euLfThWbRoEVarldGjR+N0OrnxxhtZsWJFoIYOwJtvvslDDz0UyM42YsQIli5daqq9Ujny8uA//4HvvjO3nWo0SG1ks5kr8BgCoISHNUJoqBGY5uQYxWHNysmB2FhlM6tpfE4nfylFCY6KMHr3bqwF74Zcwt/+9jfq1atX6Dmfz1fo8T333EPv3r35/vvvSUhI4IcffuBvf/tboZE2YMzFvvXWWwHj+m/9+vUsX76cxx57LLDO73//ewYOHFhiexYvXszMmTMD+/nTn/7Ehg0bLvoenn/+ecaOHRuYmN+mTRtefPFF+vbty7JlywgPD+fll1/m6quv5qGHHuKvf/0rTz/9dJEyKddddx1PPPEEAG3btuWzzz5j0aJFDBw4kI8//pi9e/eSnp5OWFgYAP/7v//Le++9xzvvvMPkyZOZN28ev/rVr5gzZ05gn9dccw1gZBsODQ2lfv36hcqpXKrtx44d4x//+Eeh8inLly/nyiuvvOjvpC4pdWTSqlUrmjVrxv33389DDz1Ely5dCgUSwbjUnQWHw3HJP2CA8PBwlixZUqQbsaBGjRqxatUq022UyhUSYqSIdThUo0FEarfwcEhKMnrvzHI6jWGNeXnl3y6RfP379w/clM63a9cu7rzzzsDjHj160KFDB9544w2eeOIJVq5cSbNmzYpMZ0hKSgr8bLVa6datGwcOHCi0Tn4PR3EyMzM5efJksfu52PVkcnIyhw8fLjSky+/3k5eXR2pqKldeeSUNGzZk+fLlDB48mN69ewcCmpLan/84P3lAcnIy58+fLzKfyOl0cuTIEQBSUlKYNGlSie0Mpu3ffPNN4HeQr3379jRo0MDUcWqzUgc6t99+O5988gkLFixg27Zt9O3bl/79+9O5c+cSM5eJBMNigbAwc3exQcPWRKTmqYBqDVLNhTocjN69u8qObUZkZCStW7cu9NyJEyeKrHfPPfewdOlSnnjiCV577TV+/etfl+ra8MJ1Is3e4SyFvLw87r33Xh566KEirzVr1izw8yeffEJoaCjff/99YO7RpeS3Py8vj/j4+ELzjvLlBx0Ok7/70rT94MGDhdohRZV6js7atWs5efIkO3bsYMiQIXz++ef88pe/pGHDhgwbNoznn3+e3VX0wRURERGpCSwWC9aIiCpZKuqC+M477+TYsWO8+OKLfPXVV9x1111F1tm5c2fgZ6/XS3JyMu3bty/1MaKjo4mPjy92PxfTpUsXvvrqK1q3bl1ksf807nP79u0sXLiQDz/8kKioKKZMmXLR9uc/zm9/ly5dSEtLw2q1FjlGkyZNALj66qsvmirbbrcXGRZ4qbZfeeWVeL1evvjii8A2Bw8erJTU3zWF6WQE7du357777mPt2rWkpaWxfft2rr32Wp555pki3XoiIiIiUrs1bNiQW265hRkzZjBo0KBiS4X88Y9/ZN26dXz99dc88MADZGRkMGHCBFPHefjhh3n22WcD+7n//vsveVH/+OOPs2PHDh544AFSUlI4dOgQH3zwQSCYOXfuHOPGjWPKlCkMGTKEt956i7/85S+8/fbbhfbz2WefsXDhQr755hv++Mc/8vbbb/Pwww8DMGDAAJKSkhg5ciQbNmzg6NGjbN++nSeffDIQhDz99NOsXr2ap59+mgMHDrB3714WLlwY2H+LFi345JNP+O677/jhhx9K1fZ27doFkm3t2rWL5ORk7rnnnqB6j2qroLKunTp1irVr13Lfffdxyy23MH/+fNxuN7/4xS/Ku30iIiIiUs1NnDgRt9tdYvDy7LPP8txzz3HNNdfw6aef8v777wd6O0pr+vTpjB8/nrvvvpukpCTq168fyMJbkquvvpqtW7dy6NAhfvGLX9C5c2dmzZoVKFPy8MMPExkZyfz58wHo0KEDzz33HL/5zW/4rkBWpOnTp5OcnEznzp2ZO3cu//d//8fgwYMBo5fuo48+ok+fPkyYMIG2bdvyq1/9iqNHjwbqNfbr14+3336bDz74gGuvvZYbbriBXbt2Bfb/+9//nqNHj9KqVSsuu+yyUrUdjAReiYmJ9O3bl1tuuSWQ4loMFn8pcw2+/fbbgYxrBw8exGq10qNHD/r370///v3p3bt3INNEbZSVlUV0dDSZmZnlWlMnJQXuvBMaNSr/TGG5uXDggJGRx2RCvKBF2U7zxe+MD1i3eelkeS4ztb3LZUyyve02MPs5zcyE06dh/nwoMOxWpEbKzIS334YGDeCCpEcXFZJ7mkE5xoen/4vp1LvM3GewLLKzISsLHngACpyHS+X8eTh7Fm6/vWZlTQz2/6mssrNhyxYjaYvZeT45ObBjh5HxzUzGtgjLadb+yvjb+s2GdOJamvvbqqrv6Io6f5dWbm4uqamptGzZkvBaPCnrzTff5OGHH+b7778PDAmrLVq0aMHUqVMD2c+kapn5TJX68veOO+6gW7dujBo1iv79+3Pdddepa6yGCCajDxhJAcqSWM+GE5vJPLd5gBMHqnQoIiJS/eXk5JCamsqCBQu49957a12QIzVbqQOdjIyMCsmGIRUnJMSoweHxBBfseL3G3cJgg52ptsH4zI6OtMF/wzrzvX8lCnZEqobHY9SkCmY7EalbFi5cyLx58+jTpw8zZ86s6uaIFFLqQEdBTs1jt0PbtsHVWXC74cgR8xc7Hsrey3eFbQ9peU6g9EXNRKR8eDzw44/BpT32+YzvnTKWWBORGmT27NnMnj27qptRoY4ePVrVTZAglTrQKW1x0AtT40nVqvwe5J97YZ73fILTU/oA2Y6TGbY+l15RpI5wuYy5GGYywoaUsZ6U328EOZ07mw92cnONRffFRESkOih1oOP3+2nevDl33XUXnTt3rsg2SS3hIQKPemVEgpKdDXv2GEGOmYAj3A+jehs/l6VkRni4MVndjJCQ4Ia8Sc3jdhtJDcxwOo3tREQqS6kDnV27dvHqq6/yhz/8gZYtWzJhwgTuuOMOGjZsWJHtExGpk7xeo3ckOtpcD4m9wFBVDSGTinL4MBw6YW4bl8sYGpmVVTFtEhG5UKkDne7du9O9e3cWLVrEO++8w2uvvcbjjz/O8OHDmThxIgMHDqzIdkoN5LBnm1rfRg6htjx8HiUhEMlnt5vr0bFr9LBUgrAwCDHZ4+fzGb1AXm/FtElE5EKmq6uEh4dz5513cuedd5KamsrEiRO56aabOH36NI0aNaqINkoNtXNGbFDbpR93sP2/Gv8iIlJd2WwQanIOaGXVcxMRyWcy96/hxIkTPPPMMwwcOJCDBw8yY8aMKinCJdWP0xNB8rHryrSPmEQn9hCTg79FRERERAoo9f0Vt9vNunXrWL58OZ9++ilDhgxh8eLF/PKXvyQkJKh4SWolC//z6qc4bOYDlSj7aT6d0bIC2iQiIiIidU2pI5T4+Hgef/xxkpKS2Lt3LytWrKBPnz6cP3+erKyswCICFpyeSPOLWzlpRQpyu41MVTk5pV+czqputYhI6bz33nu0bt2a0NBQpk6dWtXNKWT58uUMGjSoqptR5yxdupQRI0aU2/5KHehkZGRw7Ngx5s6dS7t27WjYsGGhpUGDBsrAJuXG5TZ3cZd/gafUpVJbZGXBoUPw5ZewY0fply+++HkfZUkvLXIxLtfPdZNKu7jdSkRQk9x9991YLBYsFgs2m43Y2FgGDhzIq6++Sl4JlcgHDRpEaGgoO3fuLNUx7r33Xm677TaOHz/O3Llzy7P5xbJYLLz33nuXXM/lcvHUU08xa9aswHOzZ8/m2muvLbLu0aNHsVgspKSkFHqcv0RHR9OrVy8+/PDDQtutWLGi0HqxsbEMHz6cr776qtB6Bf8frFYrzZo147777iMjI8P0+68JJk2axO7du9m2bVu57K/UQ9c+/vjjcjmgSGl8fQDcJlPjKnWp1CZer/H37HCYq2fjKHD7SpO/paKcPAlZHnPbeDzG3/XZsxXSJKkAN910E6+99ho+n49Tp06xfv16Hn74Yd555x0++OADrAW+ZI4dO8aOHTt48MEHWb58Ob169brovs+fP096ejqDBw8mISGhot+KKe+++y716tXjF7/4RdD72Lx5Mx06dODs2bO89NJL3HrrrXz55Zd07NgxsE5UVBQHDx7E7/fz3Xff8dhjjzF06FC++eYb7AUqvuf/P3i9Xvbv38+ECRM4e/Ysq1evLtP7LG9ut7tQu4MRFhbG2LFjWbJkCddff32Z21TqHp2+ffuWahEpD/kXd2YWq/XnE6lIbWG1GimmS7vYbFXd4rrJ4zHfw5Hfy1EThYSY+7u0241tfD7jd1Wn+f3gza6axWRF37CwMOLi4mjatCldunTht7/9Le+//z7/+Mc/WLFiRaF1X3vtNYYNG8Z9993H2rVryc4uucTEli1bqF+/PgA33HADFouFLVu2APDZZ5/Rt29fIiIiaNiwIYMHDw70Xqxfv57rr7+eBg0a0LhxY4YNG8aRI0cC+3W73Tz44IPEx8cTHh5OixYtWLBgAQAtWrQAYNSoUVgslsDj4qxZs6bMw6caN25MXFwc7du3Z968eXg8niKdBhaLhbi4OOLj4+nWrRuPPPII3377LQcPHiy0Xv7/w+WXX86gQYMYM2YMGzduvGQbXn31VTp06EBYWBjx8fE8+OCDgdeOHTvGzTffTL169YiKimL06NGcOnUq8PqRI0e4+eabiY2NpV69enTv3p3NmzcX2n+LFi145plnuPvuu4mOjmbSpEkX/T8AyMzMZPLkycTExBAVFcUNN9zAv//970L7HTFiBO+99x7OchiLXar7fdnZ2USaqFhndn2pvnxB1uSwWMpWrNBqBYtSl4pINRcaatyYcbuhhNE8F+V0QkyMEQjUJFYrWE1WAVAB25/4cuAv9arm2KPPg7Vs12c33HAD11xzDX/961+55557APD7/bz22mv88Y9/pH379rRt25a//OUv/PrXvy52H7179+bgwYO0a9eOd999l969e9OoUSNSUlK48cYbmTBhAi+++CJWq5WPP/4Y308XI9nZ2UybNo1OnTqRnZ3NU089xahRo0hJSSEkJIQXX3yRDz74gL/85S80a9aM48ePc/z4cQB2795NTEwMr732GjfddBOhF/mD/PTTT7njjjvK9HvK5/F4eOWVVwCwXeRO1NmzZ3nrrbcuud5///tf1q9ff9F1AJYtW8a0adN49tlnGTJkCJmZmXz22WeA8f81cuRIIiMj2bp1K16vl/vvv58xY8YEAs7z58/zy1/+kmeeeYbw8HBef/11hg8fzsGDB2nWrFngOM8//zyzZs3iySefBLjo/4Hf72fo0KE0atSIjz76iOjoaF5++WVuvPFGvvnmm0CZmm7duuHxePj888/L3IlSqkvD1q1bM2XKFO6+++4Suxf9fj+bN2/mhRdeoE+fPsycObNMDZOqFRJi3Bn2eIILdrxeo5q7TmwiUpuFhcE11xgX/mZ703JzISUluABJpCq1b9+e//znP4HHmzdvJicnh8GDBwNw5513snz58hIDHbvdTkxMDACNGjUiLi4OgIULF9KtWzdeeumlwLodOnQI/HzrrbcW2s/y5cuJiYlh//79dOzYkWPHjtGmTRuuv/56LBYLzZs3D6x72WWXAdCgQYPA8Ypz9uxZzp49W+z17t69e6lXr3CQ6i+hl6x3796EhITgdDrJy8ujRYsWjB49utA6mZmZ1KtXD7/fT06Oka12xIgRtG/fvtB6f/vb36hXrx4+n4/c3FwAXnjhhRLfA8AzzzzD9OnTefjhhwPPde/eHTD+v/7zn/+QmppKYmIiACtXrqRDhw7s3r2b7t27c80113DNNdcU2t+6dev44IMPCvUM3XDDDTz66KOBxxf7P/j444/Zu3cv6enphIWFAfC///u/vPfee7zzzjtMnjwZgMjISBo0aMDRo0crJ9DZsmULTz75JHPmzOHaa6+lW7duJCQkEB4eTkZGBvv372fHjh3YbDZmzpwZaKjUXHY7tG0b3AnY7YYjR0z3kIuI1Eh2uxHwaNiglEpohNGzUlXHLgd+vx9LgWwny5cvZ8yYMYE5O//zP//DjBkzAr02pZWSksLtt99e4utHjhxh1qxZ7Ny5kx9++CGQFOHYsWN07NiRu+++m4EDB9KuXTtuuukmhg0bZjpzWv5wqfDw8CKvtWvXjg8++KDQc9999x39+vUrsu7atWtp374933zzDVOnTuVPf/pToMciX/369fnyyy/xer1s3bqV559/nj/96U9F9tW/f3+WLVtGTk4O/+///T+++eYbpkyZUuJ7SE9P5/vvv+fGG28s9vUDBw6QmJgYCHIArrrqKho0aMCBAwfo3r072dnZzJkzh7/97W98//33eL1enE4nx44dK7Svbt26FXp8sf+D5ORkzp8/T+PGjQtt43Q6Cw1BBHA4HIHgryxKFei0a9eOt99+mxMnTvD222/zySefsH37dpxOJ02aNKFz58688sorqqlTy9S0oRQiIiLVnsVS5uFjVe3AgQO0bGnUvTtz5gzvvfceHo+HZcuWBdbx+Xy8+uqrPPfcc6Xer8PhuOjrw4cPJzExkVdeeYWEhATy8vLo2LEj7p8mu3Xp0oXU1FT+8Y9/sHnzZkaPHs2AAQN45513St2Gxo0bY7FYis1qZrfbad26daHnrCWMm09MTKRNmza0adOGevXqceutt7J///5ATxZASEhIYH/t27cnLS2NMWPG8MknnxTaV2RkZGC9F198kf79+zNnzpwSM9Vd6vd4YaBa3PMzZsxgw4YN/O///i+tW7fG4XBw2223BX7XBdtW0MX+D/Ly8oiPjw8MjyuoQYMGhR6fOXMm0AtXFqZmNVx++eU88sgjPPLII2U+sIiIXJzXa/SQ/jRSoVRCda9JKoHDno0/z1zvgN8G54kAlPe8JvvXv/7F3r17A9eCb775JpdffnmRtM3//Oc/WbBgAfPmzSsxGLjQ1VdfzT//+U/mzJlT5LUff/yRAwcO8PLLLweyoRWXgjgqKooxY8YwZswYbrvtNm666SbOnDlDo0aNsNlsgfk+JbHb7Vx11VXs37+/3Oro9O3bl44dOzJv3jz+8Ic/lLjeI488wgsvvMC6desYNWpUies9/fTTDBkyhPvuu6/YIXb169enRYsW/POf/6R///5FXr/qqqs4duwYx48fD/Tq7N+/n8zMTK688krAmKd09913B9px/vx5jh49Wqr3W9L/QZcuXUhLS8NqtV40GcSRI0fIzc2lc+fOpTrexWj6tohINXT2LJw+DRkZ5oZERRVYV3V0pKLsmNkyuO0OX8d5/6co2KkZXC4XaWlphdJLL1iwgGHDhjF+/HjAGLZ22223FUqbDNC8eXMef/xx/v73v3PzzTeX6ngzZ86kU6dO3H///fzmN7/Bbrfz8ccfc/vtt9OoUSMaN27Mn//8Z+Lj4zl27BhPPPFEoe0XLVpEfHw81157LSEhIbz99tvExcUFegvyL/6vu+46wsLCSqz/OHjwYLZt21auRUynT5/O7bffzmOPPUbTpk2LXScqKop77rmHp59+mpEjRxbb6wLQr18/OnTowPz581m6dGmx68yePZvf/OY3xMTEMGTIEM6dO8dnn33GlClTGDBgAFdffTV33HEHixcvDiQj6Nu3b2AoWuvWrfnrX//K8OHDsVgszJo1q8T6SQVd7P9gwIABJCUlMXLkSJ577jnatWvH999/z0cffcTIkSMDx/7000+54ooraNWqVWl+tRele38iItVQfiKQYNL45tOcESlPLl8E6ccvPiTmUpJaf0YIZR93L5Vj/fr1xMfH06JFC2666SY+/vhjXnzxRd5//31CQ0NJTk7m3//+d5EkAWD0KgwaNIjly5eX+nht27Zl48aN/Pvf/6ZHjx4kJSXx/vvvY7VaCQkJYc2aNSQnJ9OxY0ceeeQRnn/++ULb16tXj+eee45u3brRvXt3jh49ykcffRSYVvF///d/bNq0icTExIv2FkyaNImPPvqIzMzMUrf9UoYNG0aLFi2YN2/eRdd7+OGHOXDgAG+//fZF15s2bRqvvPJKIKPZhe666y4WL17MSy+9RIcOHRg2bBiHDh0Cfi6c2rBhQ/r06cOAAQO44oorWLt2bWD7RYsW0bBhQ3r37s3w4cMZPHgwXbp0ueT7vNj/gcVi4aOPPqJPnz5MmDCBtm3b8qtf/YqjR48SGxsb2Mfq1auZNGnSJY9VGhZ/SekipJCsrCyio6PJzMwkKiqq3PabkgJ33gmNGkF0dLnttkrl5sKBA8YFl5mUz1G203zxO2Ps6sjX0/HZzY3NzMyEM2dg1SoopnixSI2yaRPcfjvUq2euYGhZP0duN+TkQFKSueOC8dnPzoZ+/Yysi2acP2/0Yt1+e836LvR4IDU1uGQETifs3Gn8roqZ93xR2dmwZUtw2+bkwI4dxv+vmbmYoXk5jD7ajVCbnwXOT/BbS/8H4rBns3OGcSHzzybnuXFQ5c1Rqajzd2nl5uaSmppKy5Yti53gLtXT6NGj6dy5s7IIV7J9+/YF0k1Hl3AyMPOZ0tA1ERERKSULPo8FpycSv798MniJVEfPP/98kQxrUvG+//573njjjRKDHLNMDV3zer3MmTOnxG4ysxYsWED37t2pX78+MTExjBw5skg1WL/fz+zZs0lISMDhcNCvXz+++uqrQuu4XC6mTJlCkyZNiIyMZMSIEZw4caLQOhkZGYwbN47o6Giio6MZN24cZ8+eLZf3ISIidZvXa/TumF2CLcosIhWrefPmF03hLBVj0KBBgXpM5cFUoGO1Wnn++ecvmbGitLZu3coDDzzAzp072bRpE16vl0GDBpGdnR1YZ+HChbzwwgssXbqU3bt3ExcXx8CBAzl37lxgnalTp7Ju3TrWrFnDtm3bOH/+PMOGDSvUzrFjx5KSksL69etZv349KSkpjBs3rlzeh4iI1E0WizFkzecDl8v8kpNTM4MdO05s5JhaREQqm+mhawMGDGDLli3cfffdZT74+vXrCz1+7bXXiImJITk5mT59+uD3+1m8eDG/+93vuOWWWwB4/fXXiY2N5a233uLee+8lMzOT5cuXs3LlSgYMGADAqlWrSExMZPPmzQwePJgDBw6wfv16du7cSc+ePQF45ZVXSEpKMl3MSkREJJ/VComJwRVIzsoytq+JM2VnOvqYWj/UViBbU018wyJSI5kOdIYMGcLMmTPZt28fXbt2LVIoaMSIEUE3Jj+7RX7l2NTUVNLS0grlMQ8LC6Nv375s376de++9l+TkZDweT6F1EhIS6NixI9u3b2fw4MHs2LGD6OjoQJAD0KtXL6Kjo9m+fXuxgY7L5cLlcgUeZ2VlBf2+RESk9jKTdKU8tqsqPouDU/bOxLr3lG1HbidQr1zaVJMo95NI+TDzWTL9NXvfffcB8MILLxR5zWKxBD2sze/3M23aNK6//vpALva0tDSAQinn8h9/++23gXXsdnuRXOixsbGB7dPS0gpVos0XExMTWOdCCxYsKLZglYhIZfD7/dgtTmyAmWReNpxlPrbPZ2RQCzFZgOCCgtlS21gsfNR4JYe/dprPqskPjKF5xbWtGgsNDQXA7XZfsmK9iFxaTo4xFNZWilSXpgOd0hQLCsaDDz7If/7zn2Kr3F5YMMnv95dYRKmkdYpb/2L7mTlzJtOmTQs8zsrKClSPFRGpSH6/n7Mv3snihBTT2xYaIhQEr9dIXZyTE9wIo4gI+Om6TmojiwUPEVgAM38eHuruBb7VaiUiIoLTp09js9kCNV1ExBy/309OTg7p6ek0aNAgcBPhYsrUcZ6bm1suOeGnTJnCBx98wCeffMLll18eeD4uLg4wemTi4+MDz6enpwd6eeLi4nC73WRkZBTq1UlPT6d3796BdU6dOlXkuKdPny7SW5QvLCyMsLCwMr83ERGzfE4n3tSUsu/HYv7i0u83Jtd37QrFdIRfUmio+bouIrWZxWIhPj6e1NTUwGgUEQlegwYNAjHCpZgOdHw+H/Pnz+dPf/oTp06d4ptvvuGKK65g1qxZtGjRgokTJ5Z6X36/nylTprBu3Tq2bNlCy5YtC73esmVL4uLi2LRpU6CCrdvtZuvWrTz33HMAdO3aFZvNxqZNmxg9ejQAJ0+eZN++fSxcuBCApKQkMjMz+fzzz+nRowcAu3btIjMzMxAMiYhUR3OyPsFqYriLg2zG8FOEcome74sJDzdf9FNEime322nTpg1uje8UKRObzVaqnpx8pgOdefPm8frrr7Nw4UImTZoUeL5Tp04sWrTIVKDzwAMP8NZbb/H+++9Tv379wHyZ6OhoHA4HFouFqVOnMn/+fNq0aUObNm2YP38+ERERjB07NrDuxIkTmT59Oo0bN6ZRo0Y8+uijdOrUKZCF7corr+Smm25i0qRJvPzyywBMnjyZYcOGKeOaiFRrbr8DP6UvzGgtMKDI4wG/yZFsHo+59aVqBXPdnJtbM1Na13QhISHlMgpGRErPdKDzxhtv8Oc//5kbb7yR3/zmN4Hnr776ar7++mtT+1q2bBkA/fr1K/T8a6+9Fkhf/dhjj+F0Orn//vvJyMigZ8+ebNy4kfr16wfWX7RoEVarldGjR+N0OrnxxhtZsWJFoYjvzTff5KGHHgpkZxsxYgRLly411V4RkcrjJ9Tmx2HLxmYr/WwIh/3nOmROJ+SZnA7g84Hdrnk21V1oqDEfKifHfHDqdBrzsBwO4/9aRKS2Mh3ofPfdd7Ru3brI83l5eXhMftuWJj2cxWJh9uzZzJ49u8R1wsPDWbJkCUuWLClxnUaNGrFq1SpT7RMRqRJ+PwPHH+WyROfPw9CCcNVVUL/hpdcrKDfXWDRsrXoLD4ekpOB6ZtLTISVF5WxEpPYzHeh06NCBTz/9lObNC6eJfPvttwPzaEREpAx8OVyWWLY00bv+ex0WWwQRpR/1BhgppXUBXDMEOwpKo6dEpK4wHeg8/fTTjBs3ju+++468vDz++te/cvDgQd544w3+9re/VUQbpQ7yeMBtcm6B222kxhWpTbrNScVjvczUNi4XnMmK4Lbbgk9GECy3G4KpQqC5QSIiUt5MBzrDhw9n7dq1zJ8/H4vFwlNPPUWXLl348MMPGThwYEW0Ueqg1FTIMnnh4/EYgc7ZsxXSJJEqkeOOxOs3N47MVUVBg9ttDIsKpiai12tsp7lBIiJSXoKqozN48GAGDx5c3m2RWsbs2HFvgZvPdjvYTd6M9vmMRXeGRapGXp4RrFx7rfnhUfk3KlS+TEREyovpQOd3v/sd/fr147rrriPC7OBvqRNCQoyCgx6PuWDHXWBegMUCVpN/nboTLFL1fD7js2j28wvGd4eIiEh5MX0qSk5OZsmSJbhcLrp06UK/fv3o27cv119/PfXq1auINkoNY7dD27bmx+lHFrjI0WRokZrH5zPSHbtcwQUtNluZapyKiIgUYjrQWb9+PT6fj88//5ytW7eyZcsWXnrpJZxOJ126dGHnzp0V0U6pYYKpzRCmCxyRGs3vN3pyLr8coqLMbx9MT66IiEhJgjqlhIaGkpSURKNGjWjYsCH169fnvffe48iRI+XdPhERqWGsVqN3RkREpCqZHlywbNkyfvWrXxEfH88vfvELNm7cyC9+8QuSk5M5ffp0RbRRRERERETEFNM9Og888ACXXXYZ06dP5ze/+Q1RwYxPEBERkRqrLFk1RUQqi+kenb/+9a/ccccdrFmzhpiYGHr27Mnjjz/OP/7xD86fP18RbawTvF6jBkVxi9Ili4hIdZCfVdPnK/mcVdKSL/Nc1bVfROoW0z06I0eOZOTIkQBkZmby6aef8s4773DzzTdjsVhwuVzl3cZaLX8su9drZCsqTm4uNG6sMe8iIlK1yiOrpm7eiUhlCSoZwZkzZwIZ17Zs2cK+ffto3Lgxffv2Le/21XpRUdCmDTRqVHyWotxc2LNH6ZZFRKR6UFZNEakpTAc6V199Nfv376dRo0b06dOHSZMm0a9fPzp27FgR7asT7Hajmrjqr4qISGXxeAoPKSsti0UjDESkZjAd6EyePFmBjYiISA0VGmrcYMsv8GqWhlOLSE1hOtB58MEHAz/7fxpPZVEpaxERkRohMhI6doTwcGMxQ8OpRaQmMZ11DeCNN96gU6dOOBwOHA4HV199NStXrizvtomIiEgFCAszhkubXcwGRiIiVcl0j84LL7zArFmzePDBB7nuuuvw+/189tln/OY3v+GHH37gkUceqYh2iojUSXl5RlZGM8zWOBEREamNTAc6S5YsYdmyZYwfPz7w3M0330yHDh2YPXu2Ah0RkXKUlwduk4FLXp6Ruj40tGLaJCIiUhOYDnROnjxJ7969izzfu3dvTp48WS6NEhERQ5MmENnE3DYejzGXwuGomDaJiIjUBKbn6LRu3Zq//OUvRZ5fu3Ytbdq0KZdGiYiIwWb7edJ4aZewMGXEEhERMd2jM2fOHMaMGcMnn3zCddddh8ViYdu2bfzzn/8sNgASERERERGpbKZ7dG699VY+//xzmjRpwnvvvcdf//pXmjRpwueff86oUaMqoo0iIiIiIiKmmOrROXfuHDt37sTj8bB48WKaNDE5cFyC5vEEt50qWIuIiIhIXVTqQOc///kPQ4YMIS0tDb/fT1RUFO+88w4DBgyoyPbVeSEhxpj73Nzggp2aWsHajhMP5kp2G2/RAaiArYiIiEhdV+pA54knnqBZs2a8/fbbhIeHM2fOHB588EG+/vrrimxfnRceDl27GulizarJFaynOwbj85gcWRkNh8M74/evRMGOiIiISN1W6kDniy++4KOPPqJbt24AvPrqq8TExHD+/Hnq1atXYQ2UulOJ2mf5OReu1RZEZAe05ktwO4GIcmqViIiIiNREpQ50fvjhB5o1axZ43LhxYyIiIjh9+rQCHSkflp97YW595JugdpF+3EFKTezCEhEREZFyVeqxQRaLhXPnzpGVlUVWVhaZmZlFnsvKyjJ18E8++YThw4eTkJCAxWLhvffeK/T63XffjcViKbT06tWr0Doul4spU6bQpEkTIiMjGTFiBCdOnCi0TkZGBuPGjSM6Opro6GjGjRvH2bNnTbVVKp7LF8G+U9eVaR8xiU5CTc7tEREREZHap9Q9On6/n7Zt2xZ5rnPnzoGfLRYLPp+v1AfPzs7mmmuu4de//jW33nprsevcdNNNvPbaa4HHdru90OtTp07lww8/ZM2aNTRu3Jjp06czbNgwkpOTCQ0NBWDs2LGcOHGC9evXAzB58mTGjRvHhx9+WOq2SmWw8MhHn3L0cA52O1hN5ASMsp/m0xktK65pIiIiIlKjlPpS8uOPPy73gw8ZMoQhQ4ZcdJ2wsDDi4uKKfS0zM5Ply5ezcuXKQPa3VatWkZiYyObNmxk8eDAHDhxg/fr17Ny5k549ewLwyiuvkJSUxMGDB2nXrl35vikpIwtOTyQ+C1hNjECzqRdHRERERAoodaDTt2/fimxHibZs2UJMTAwNGjSgb9++zJs3j5iYGACSk5PxeDwMGjQosH5CQgIdO3Zk+/btDB48mB07dhAdHR0IcgB69epFdHQ027dvLzHQcblcuFyuwGOzw/JERGoqr9fI2vhTp3ipBVvvS0REpCKYKhha2YYMGcLtt99O8+bNSU1NZdasWdxwww0kJycTFhZGWloadrudhg0bFtouNjaWtLQ0ANLS0gKBUUExMTGBdYqzYMEC5syZU75v6CLcbuPCojghIXDBiD0RkYvy+YzvlBCTWdpzciAzE7KyjIDH7DEdDvMBkoiISEWo1oHOmDFjAj937NiRbt260bx5c/7+979zyy23lLhd/nyhfAV/LmmdC82cOZNp06YFHmdlZZGYmGj2LVyS1Wqkj3a7ITu7+HVyciA2VsGOiJSO12t8n+TkmK+j5XIZ30lXXQUNGpg/dlhY3UmJLyIi1Vu1DnQuFB8fT/PmzTl06BAAcXFxuN1uMjIyCvXqpKen07t378A6p06dKrKv06dPExsbW+KxwsLCCAsLK+d3UFRkJHTuDPXqGcuFnE7YtSu4gqEiUjf5/WCzGcWGi+nQvqjz542lY0eIijJ/bIvFXCIRERGRimJyUEPV+vHHHzl+/Djx8fEAdO3aFZvNxqZNmwLrnDx5kn379gUCnaSkJDIzM/n8888D6+zatYvMzMzAOlUtLMwIeIpbHI5Lby8iUpzw8JK/Wy62hIUZwYrNZn5RkCMiItVFqU9JPp+Pr776ijZt2uC44Oo7JyeHw4cP07FjR0JMDAg/f/48hw8fDjxOTU0lJSWFRo0a0ahRI2bPns2tt95KfHw8R48e5be//S1NmjRh1KhRAERHRzNx4kSmT59O48aNadSoEY8++iidOnUKZGG78soruemmm5g0aRIvv/wyYKSXHjZsmDKuiYiIiIjUUqWOSlauXMmECROK1LEBY5jXhAkTeOutt0wd/IsvvqBz586BWjzTpk2jc+fOPPXUU4SGhrJ3715uvvlm2rZty1133UXbtm3ZsWMH9evXD+xj0aJFjBw5ktGjR3PdddcRERHBhx9+GKihA/Dmm2/SqVMnBg0axKBBg7j66qtZuXKlqbaKiIiIiEjNUeoeneXLl/Poo48WCiDyhYaG8thjj7F06VLuvPPOUh+8X79++C8yU3bDhg2X3Ed4eDhLlixhyZIlJa7TqFEjVq1aVep2iYiIiIhIzVbqHp2DBw/Sq1evEl/v3r07Bw4cKJdGiYiIiIiIlEWpA53s7OyLFs08d+4cOTmqTi8iIiIiIlWv1IFOmzZt2L59e4mvb9u2jTZt2pRLo0RERERERMqi1IHO2LFjefLJJ/nPf/5T5LV///vfPPXUU4wdO7ZcGyciIiLFc7shNze4xeut6taLiFS8UicjeOSRR/jHP/5B165dGTBgAO3bt8disXDgwAE2b97MddddxyOPPFKRbRURERGMICc9Pbhaazk5cO6cUTNJRKQ2K3WgY7PZ2LhxI4sWLeKtt97ik08+we/307ZtW+bNm8fUqVOx2WwV2VYREREB8vKMIOfaa43CsGacPg3bt8NFkp6KiNQKpmpY22w2HnvsMR577LFiX09JSeHaa68tj3aJiIjIJYSHm+/VMRsYiYjUVKWeo1OSzMxMXnrpJbp06ULXrl3Lo00iIiIiIiJlEnSg869//Ys77riD+Ph4lixZwi9/+Uu++OKL8mybiIiIiIhIUEwNXTtx4gQrVqzg1VdfJTs7m9GjR+PxeHj33Xe56qqrKqqNIiIiIiIippS6R+eXv/wlV111Ffv372fJkiV8//33LFmypCLbJiIiIiIiEpRS9+hs3LiRhx56iPvuu0+FQUVEREREpFordY/Op59+yrlz5+jWrRs9e/Zk6dKlnD59uiLbJiIiIiIiEpRSBzpJSUm88sornDx5knvvvZc1a9bQtGlT8vLy2LRpE+fOnavIdkoZeDxGcTmzi8dT1S0XEREREQmO6axrERERTJgwgW3btrF3716mT5/Os88+S0xMDCNGjKiINkqQQkKMegkej1EJ2+zy448KdkRERESkZjKVde1C7dq1Y+HChSxYsIAPP/yQV199tbzaJeUgPBy6djUqaJuVmwt79qhytoiIiIjUTGUKdPKFhoYycuRIRo4cWR67k3KkCtgiIlKd+POKHy1gsYC1XK5KREQM+koRERGRSvPDj5CaWvR5mw0SExXsiEj50deJiIhIDeTzGT0jZgMDj8cY0uzxGMOUzcjNNY5bFnYbhIUVfs7rNdqj4dIiUp4U6MhFBZuMwGIx7s6JiEj58/mMpDEul5F4xuy2oaFG0GI2sHA6ITsbHA6w281tmy80tPjzQ1kDKBGRCynQkWLlZ2zLzQ0u2MnNhcaNgw92zJ7wvJbgjiMiUhP5/UZPzuWXQ1SUuW0bNzYS1dSrZyxmpKdDSop6XkSkZlCgI8WqqoxtISFGcOTxmAt23AWOlamSTiJSR1it5m8oWa3G0LHISGMxozwS3Hg8Rs/Qhc95vWXft4hIQQp0pERVkbHNboe2bc0HWJEFhm6o9o+ISPV17FvYubPwc16vEYC1aKFhzyJSfhToSLUTzLjvMA1dExGpEezhRXuS8uf+aJ6OiJQnBToiIiJSaey2oiMGNGxNRCqCAp0awu0ObruQkOAz44iIiIiI1FQKdKqJkuoZeDw/T84PZu5JTg7ExirYEREREZG6RYFOFQsNNeoRuN0lT8Bv1QoaNTIfrDidsGtXcJnTRERERERqMpNlxsrXJ598wvDhw0lISMBisfDee+8Vet3v9zN79mwSEhJwOBz069ePr776qtA6LpeLKVOm0KRJEyIjIxkxYgQnTpwotE5GRgbjxo0jOjqa6Ohoxo0bx9mzZyv43ZVOWBhcc42Ryrm4pWNHqF/fCIbyU4GWdnE4qvrdiYiIiIhUjSoNdLKzs7nmmmtYunRpsa8vXLiQF154gaVLl7J7927i4uIYOHAg5879XChl6tSprFu3jjVr1rBt2zbOnz/PsGHD8BVI3TJ27FhSUlJYv34969evJyUlhXHjxlX4+ysNi8UISux2Yz7NhQsYw8+UiUZEREREpPSqdOjakCFDGDJkSLGv+f1+Fi9ezO9+9ztuueUWAF5//XViY2N56623uPfee8nMzGT58uWsXLmSAQMGALBq1SoSExPZvHkzgwcP5sCBA6xfv56dO3fSs2dPAF555RWSkpI4ePAg7dq1q5w3WwKrFRITSy6umZVlrFMTq1AHW8/GYlEdBREREREpmyrt0bmY1NRU0tLSGDRoUOC5sLAw+vbty/bt2wFITk7G4/EUWichIYGOHTsG1tmxYwfR0dGBIAegV69eREdHB9YpjsvlIisrq9BSUfIrWxe3WGvgLKqQECN1qMdj9EaZXX78UUU/RURERKRsqu1ldFpaGgCxsbGFno+NjeXbb78NrGO322nYsGGRdfK3T0tLIyYmpsj+Y2JiAusUZ8GCBcyZM6dM76GuCg835hcFkwQhNxf27KmZPVgiIiIiUn1U20Ann8VSuOS93+8v8tyFLlynuPUvtZ+ZM2cybdq0wOOsrCwSExNL2+w678JicCIiIgAuN2RnF34uOxtcrqppj4jUXtU20ImLiwOMHpn4+PjA8+np6YFenri4ONxuNxkZGYV6ddLT0+ndu3dgnVOnThXZ/+nTp4v0FhUUFhZGWFhYubwXqVzuXMjMNL+d1WokhhARkYpz6Bs4k1P4udxcoyd/1CiIjq6adolI7VNtA52WLVsSFxfHpk2b6Ny5MwBut5utW7fy3HPPAdC1a1dsNhubNm1i9OjRAJw8eZJ9+/axcOFCAJKSksjMzOTzzz+nR48eAOzatYvMzMxAMCS1S8q/wRVqfrv69WHYMAU7IiIVKTy86PdsXp5xg8rrrZo2iUjtVKWBzvnz5zl8+HDgcWpqKikpKTRq1IhmzZoxdepU5s+fT5s2bWjTpg3z588nIiKCsWPHAhAdHc3EiROZPn06jRs3plGjRjz66KN06tQpkIXtyiuv5KabbmLSpEm8/PLLAEyePJlhw4ZVecY1qRiRkdCggbltcnPh3DmdZEVESsPjMQpdl5atQOoju73o8GYNWxORilClgc4XX3xB//79A4/z58TcddddrFixgsceewyn08n9999PRkYGPXv2ZOPGjdSvXz+wzaJFi7BarYwePRqn08mNN97IihUrCA39+Zb+m2++yUMPPRTIzjZixIgSa/dIzRceDvXqmd8uN7f82yIiUpuEhhqBis9nZMksrZACvexmAiQRkbKo0kCnX79++C+SXstisTB79mxmz55d4jrh4eEsWbKEJUuWlLhOo0aNWLVqVVmaWuU8HvMX4rm56qEQqem8XvMXhkrPLhUlMhI6djRuKJlJOhNWIAunsmqKSGWptnN0xBAaCg6HcaFjNl1zTo4xHEtzTkRqLl+euTvnYNxtt9uN7w+R8hYWBhER5gIdu6/i2iMiUhIFOtVcWBhcc83PRUXNOH0atm83gqRghmWFhBgXSyJSdeLjoM3V5rbJzTUW3eSQi9FIARGp7RTo1AB2uxHwmA10IiN/7g0KZqhATg7Exta8YKe4Gg2XohoOUl3ZbMbdczNCQjQ8SEqmkQIiUlco0KnFIiKga1fzY6kBnE7Ytcv8SbA6KK5Gw6WohoOIVJVgelbKMg+rPEYKKJAWkZpAgU4NEcxQAY/HOIlFRpoPdGqy4mo0XIpqOIhIZcvvWXE6zSec8PmMbYOdhxXsSIG6dC4RkZpPgU41Z7EYJyKPxzixmeFyGcMM6toQg+JqNFyKhq2JSGULDzd63YP9/gkLK1vgEewNtJrY0y8idZMCnWrOaoXExOCGCWRlGdtriIGISPVjtULr1sF/R1ssxj6C2a4sN9CCmdsjIlIVFOjUAMGcyMqynYiIVI6q+J4uyw00m01py0Wk5tClsIhIBfJ6zV9Qer36cpaKFWyApSBHRGoSnUtFRCqI1wuHDwcxByMHrqmQFomIiNQdCnRERCpIbi4kJxuBjpk74VY3XPNTEpGQkIppm4iISG2nQEdEpIL4fEbq4PzivaUV4gJ+miSuoUIiIiLBUaAjIlLBbDZzaYBDAEwWvRWpKfJLHxQUTC0hEZFLUaAjFxXsiSckxKhnIyIiUtCJE3D8h8LPuVxGuuusrKppk4jUTgp0pFihoRARYdx183jMb5+TA7GxCnZERKSwkBAIiyj8nM9nnDeCKWIqIlISBTpSrPBwSEoyX0wOjCEIu3apoJyIiBRlsxW9Caa6byJSEfTVIiUyM6dARESkNMJt2YSFFu7SCbeCwwb4IwBLlbRLRGofBToiIiJSaT6c1LLE1/7jPw9EVl5jRKRWU4UGERERqVAefwTpx03kWBcRKQfq0REREZEKZmHzGy0Itfn5f/U+IbpJ4aFrWVmQkQGvvBpRwvYiIuYp0JFax+0uWqPhUlTDQSqKywXZ2WAxMe0gxFVx7RGpOhZ8Hgu53kjCfYUDmlwvOD1oeo6IlCsFOlLrHD4Mh06Y20Y1HKQiZGfDnj1GkGMmuUe4H0b1Nn4O0QBjqWXcbsjNLfqcUkuLSHlToCO1TlgYhJgc/aAaDlIRvF7jgi46GiJNzK+2F0jNnpdX9KLwUtQ7KdXZqTT4Lr3wcx6P8Xk5e7ZKmiQitZQCnTrA7EVSPqu1ZqaYttkg1GShUtVwkIridhuBtJm6Uv4C67rcRs+QWRERRuFfkerGElK0jo7PZyzBFKgWESmJLu9qMasV6teHc+eCC3bOn4cWLWpmsCNSHWRlwaFDRvAdFlb67eqFAl2Mn9u0gm59zB87NFSfXameQkMhxFr0ORGR8qZApxaLjIRhw4IbjnX+PPzjH8YkfbNyczUETASMz4HHAw6H0cNSWo4C83LqR5kb9iYiIiIGBTq1XLAXSKGhxsWZ221uyA0Yc13OndPFmUg+q7XoUJ2LsSnzlIiISJlV63w+s2fPxmKxFFri4uICr/v9fmbPnk1CQgIOh4N+/frx1VdfFdqHy+ViypQpNGnShMjISEaMGMGJEyZTctVBYWFwzTXQsyf06mVu6dTJGKrj91f1uxARERGRuqpaBzoAHTp04OTJk4Fl7969gdcWLlzICy+8wNKlS9m9ezdxcXEMHDiQc+fOBdaZOnUq69atY82aNWzbto3z588zbNgwfD5fVbydGsVuNwIWq9XcEhpqvhdIRERERKQ8Vfuha1artVAvTj6/38/ixYv53e9+xy233ALA66+/TmxsLG+99Rb33nsvmZmZLF++nJUrVzJgwAAAVq1aRWJiIps3b2bw4MGV+l5qEovFCHI8HiMTjhkuV3BD3kRqI6+3+LohFxNaDreggv0MKuuVVDS7xVmkLqgNsFvA73egqqEiUl6qfaBz6NAhEhISCAsLo2fPnsyfP58rrriC1NRU0tLSGDRoUGDdsLAw+vbty/bt27n33ntJTk7G4/EUWichIYGOHTuyfft2BToXYbVCYmJww89sNmXQEQGjJsjp05CRYXwuSiuqwLrBfJbcbkhPN+bZmeX1GtvpMywV5cnIYtIIRv+0uHcDJguhiYiUoFoHOj179uSNN96gbdu2nDp1imeeeYbevXvz1VdfkZaWBkBsbGyhbWJjY/n2228BSEtLw26307BhwyLr5G9fEpfLhcvlCjzOysoqj7dUowRbW6Y8LpDM3lW2VftBmFIX5feI2u3mkhHYCwQ6EUEk9cjLM4KVa681n2I6v3CjmXTYIpfiD3WQGdWZ6Kw9Vd0UEalDqnWgM2TIkMDPnTp1IikpiVatWvH666/Tq1cvACyWwl3cfr+/yHMXKs06CxYsYM6cOUG2XIIVEmJcmOXmmgt2QgoEV14v6Ga0VCehoeZuHJRXAdvwcPO9OlarMfxUpFxZLOy9diUfvOsk3FE0kM7JMQrjrrIH0Q0pIlKCah3oXCgyMpJOnTpx6NAhRo4cCRi9NvHx8YF10tPTA708cXFxuN1uMjIyCvXqpKen07t374sea+bMmUybNi3wOCsri8TExHJ8N1Kc8HDo2tX83IKwAusr25uISDVkseAmglCKZkLyAG6/MT9URKS81KgBPy6XiwMHDhAfH0/Lli2Ji4tj06ZNgdfdbjdbt24NBDFdu3bFZrMVWufkyZPs27fvkoFOWFgYUVFRhRapHOHhRnFFM0vBu9bhtmzCQs0t4dZsHLZsRUlSbvx+P35XDnZLDjbMLkFU6hWpIXw+o+e94KJEqCJSEap1j86jjz7K8OHDadasGenp6TzzzDNkZWVx1113YbFYmDp1KvPnz6dNmza0adOG+fPnExERwdixYwGIjo5m4sSJTJ8+ncaNG9OoUSMeffRROnXqFMjCJrXPh5NaBr3tf/znAVU6lbLzOZ388Hh3FieY3zbUppSFUvvkD+HMyzMSZhSUl2e8biZph4jIpVTrQOfEiRP8z//8Dz/88AOXXXYZvXr1YufOnTRv3hyAxx57DKfTyf33309GRgY9e/Zk48aN1K9fP7CPRYsWYbVaGT16NE6nkxtvvJEVK1YQqpRCtYrHH0H6cQcxiboTLrWMVXMWpHZwOKBZM6PX/sKA5tw5I0thgwZV0TIRqa2qdaCzZs2ai75usViYPXs2s2fPLnGd8PBwlixZwpIlS8q5dVK9WNj8RgtCbX7+X71PiG5iLj1pVpaRAviVV5XWVMpHqMNBk+d2c+c4iIw0hliWloNsxhDz0yNNWpDaw2YzEhFcmIXQ5Sq/JBwiIvn0tSK1iAWfx0KuN5Jwn7mAJdcLTg+6ppRyY7FYsIRF4PaDHWOydWlZ0VwxERGRslKgI7WO2Sr0+dt4vRXTHhERERGpfAp0pMIEE3CAUUvHTHHFC51Kg+/SzW2TXyTx7NngjytSkvwsU6XlVc+iiIhImSnQkXJntRqTTd1uowCcWTk5EBsbfLBjCSJQ8vmMxUyRUpFLycw0MkkVl2XqYkILjFzTvAUREZHg6BQq5S4yEjp3hnr1jMUMpxN27TJfMLQgh9WJxeRftj8UfBbw+x1ooo6UF4/HCHRiYiAurvTbhRf4+zX7GRKpDMH02Ofmql6OiFQuBTpSIcLCjIAnsgpK0jwZ2cf8RtE/Le7dgDKvSfmy2YxeztIKU/Z7qabK0mPvdBrbOBxlG54sIlJaCnSkVvCHOsiM6kx01p6qboqISK1Vlh779HRISQG/kgqKSCVRoCO1g8XC3mtX8sG7TsIdRo+SGTk5xp3GVXYVZxQRuZhge+zN9GqKiJQHBTpSe1gsuIkgFAgxuakHcPvBouk5Ul78fkL9OUTYIdxmbjhauDWILB4iIiJSiAIdEZGK4MthNPUYvbCqGyIiIlI3mb3xLSIilSTlu+vwKjmGiIhIUNSjIyJSEUIj+AvneewxSGgKCfHmNs/JgfQzETzwgMZTioiIBEOBjtQ6ZqvQ528jUq4sFnyWSHLckOsBl8m/sVyvUU/K5TJfr0SFb0VERBToSC0SGmrUeDBbhR6MbUJDjUr2ItVB/t9xfkZAsxwOCNHgZKmGPJ6i39Fmb06JiJSGAh2pNRwOaNbMSGFqs5nbNi0NTp3SnXCpPvx+I1C56ipo08b89iEh5tOsi5SWx2O+p9HrNW5G+XxGAH/hazab8bqISHnRV4rUKjabcXFntuq22cBIpLKEhRlBvEh1EBpq/D263Uavo1nNmkGTJkWLjWZlwZkzEBVVPu0UEQEFOlKBgrnjl5urIQxS+3i95odT6nMg1VFYGFxzjdHzYvYG0enTsH27EShFXJBM0OMxf4NKRORSFOhIuSvLHb+cHDh3znzFbZHqKDvbGEJW3FCdS/F6jQs/DeWR6sZuNwIes4FOeHjFtEdEpCQ6hUq5K487fn5/xbRNpDJ5vUbg37QpXHutuW3PnzeSY9SvXyFNExERqfUU6EiF0B0/kZ/Z7UWH6lyKzwdOZ8W0R0REpC5QoCPVks9nPgOaxxPc5FiRihbMfDVlABQRESkbBTpSYYKZTO3zGUPecnLMb+90GosyVEl14XYbQ9d8PvO1cDwe4285NLRi2iYiIlLbKdCRcmexGEPWPB7jAs+M0FBo3twIcswWO8zKUtY2qV58PuPvOC4OevUyt21+prb8z5JZFosSGUj1VFyPvb63RaQi6DQo5c5qhcTE4BIK5J/sQkPNz+85cgQ2blQiA6k+8gP2YHplLBZjbo/XCy6X+e3dbqMmiYIdqQhl6bH3eosO5cwvGKoeTBEpTzoFSoUoy8VVRIRxsjO7j/waDKpZItVFXl7wPToej/F32aqV+c+C1wvHjyvol/JXlh57MIqFXnZZ0YKh589DRoaRxEZEpLwo0JFqpSwn0fxhQh6P+Zol+duanUchUhpWq/m5Y1ar0ZMTTJp2kYpSlh77H380Apm8vKI3l/Ly1JsjIuVPgY5UK2U5idps0LYtNGpkDNkxIyUFzp5Vz46IyKUE22MfHm705OTmFp2j4/MZr2mopYiUJ32lSLUT7IkuNNQYvuZwmK9Zkj/sTUREKkZkJAwYYPx8Ye9Nfg9+ZGTltklEajeTea1qtpdeeomWLVsSHh5O165d+fTTT6u6SVIBvF7jbqGZJZix5iLVXTCfBfVqSkWxWCA62rgRFRZWeImIMF6zWKq6lSJSm9SZHp21a9cydepUXnrpJa677jpefvllhgwZwv79+2nWrFlVN0/KQX6mtuIy+lxK/hwdXeRJbVDWCeM2my44pfxdamiyUqKLSHmrM18pL7zwAhMnTuSee+4BYPHixWzYsIFly5axYMGCKm6dlIewMLj8cuPOoNmha0ePGidZZamS2qAsc91AF5xScfR3JSKVqU585bjdbpKTk3niiScKPT9o0CC2b99eRa2S8ma1GkMfcnLg3Dlz2+Zn/AmmXonIpeTlmS/6WdbeRV1QiohIXVcnToU//PADPp+P2NjYQs/HxsaSlpZW7DYulwtXgaverKysCm2jlF29ejBkiPlhawA//AD79sG//mW+Bo9ISbZtM3pH8vKCC6I1hExERCR4dSLQyWe54IrB7/cXeS7fggULmDNnTmU0S8qJ1Qrt2gU3XKdTJ/joIzh82FhEyovFYqQ9b9kyuG3VMyMiIhKcOnEKbdKkCaGhoUV6b9LT04v08uSbOXMm06ZNCzzOysoiMTGxQtspZRfsReHEicbQte+/L9/2iCQkwN13V3UrRERE6p46EejY7Xa6du3Kpk2bGDVqVOD5TZs2cfPNNxe7TVhYGGFhYZXVRKkGdDEqIiIiUnvUiUAHYNq0aYwbN45u3bqRlJTEn//8Z44dO8ZvfvObqm6aiIiIiIiUszoT6IwZM4Yff/yR3//+95w8eZKOHTvy0Ucf0bx586pumoiIiIiIlDOL36/KIaWRlZVFdHQ0mZmZREVFVXVzREREpBR0/hapu0KqugEiIiIiIiLlTYGOiIiIiIjUOgp0RERERESk1lGgIyIiIiIitY4CHRERERERqXUU6IiIiIiISK2jQEdERERERGqdOlMwtKzyyw1lZWVVcUtERESktPLP2yobKFL3KNAppXPnzgGQmJhYxS0RERERs86dO0d0dHRVN0NEKpHFr1scpZKXl8f3339P/fr1sVgs5bbfrKwsEhMTOX78eI2p2FwT2ww1s91qc+VQmyuH2lw51ObC/H4/586dIyEhgZAQjdgXqUvUo1NKISEhXH755RW2/6ioqBpzQspXE9sMNbPdanPlUJsrh9pcOdTmn6knR6Ru0q0NERERERGpdRToiIiIiIhIraNAp4qFhYXx9NNPExYWVtVNKbWa2Gaome1WmyuH2lw51ObKoTaLiBiUjEBERERERGod9eiIiIiIiEito0BHRERERERqHQU6IiIiIiJS6yjQERERERGRWkeBjoiYUtPyl9S09kLNbLOIiEh1o0CnArlcLrZv3863335b1U0ptfwLLK/XW8UtKT2Px8N3330XeFwTLhI9Hg9Hjx7F4/EANaPNubm5/Pa3v+XPf/5zVTel1DweD+np6YHHNeH37PP5OHfuXFU3wxSn08mTTz7Jhg0bqroppVYTv59B39EiImYo0KkgixcvpkWLFtx7771cc801/OlPf+L/t3fncVHV+x/H3wMiIPu+6SgisriBkJKSYVzTwAVUROPhkl5K6+byQ9PM7bpkWRqmaZkieC0F1PSqLCapgJayRCBk4IZbikuCyerM5/cHd06OLELKmRn4PB8PHw85nDPzYhy/zHfmO2dkMpmqsxq1ceNG+Pn5AQDatWunEb+M1q5dC2dnZwQGBmL48OH48ccfIZFIVJ3VqHXr1qF3794ICQmBn58fzpw5A4lEArlcruq0Bm3btg22trY4c+YMdHR0UFFRoeqkp1q3bh3c3NwwevRojB07Fj/99JPa3zc+++wzDBw4EEFBQZg9ezYuXLgAAGp931izZg1MTEyQnp6O3r17qzqnSTRxfAZ4jGaMsWYj9twtWrSIXFxc6PDhw1RUVESLFy8mU1NTKi8vV3Vava5evUqTJk0iKysrkkgktHr1aiIievTokYrLGrd27Vrq0qUL7dmzh6KiomjUqFFkaWlJx44dU3VavXJycsjX15ecnZ0pNjaWdu3aRQEBAdSzZ09VpzXq/Pnz9NJLL9GXX36p6pQmkcvlNGfOHHJxcaH9+/fTrl27aPz48WRubk4nTpxQdV69CgsL6ZVXXiFnZ2eKjo6mBQsWkI+PDw0ZMkTVaQ368ccfqXv37uTg4EAHDhxQdU6Tadr4TMRjNGOM/V080XnO7ty5Q/3796dPP/1U2FZYWEju7u50+/ZtIqp9IKZOkpOTKSQkhOLi4mjp0qVkYWFBJSUlREQkk8lUXFeXTCajmpoaeu2112jGjBlK3/P19aVhw4bRzz//rJq4Rnz11Vc0YcIE4X5ARLR582YaPHgwVVRUEJH63TeIah+seHp6EhFRcXExLV68mLZu3UppaWkqLqvflStXyMPDg+Lj45W2d+zYkYYOHUo5OTkqKqufTCajtWvXUkBAAF2/fl3YHh8fT3369KFff/1VhXUNe//990lHR0e4H1y/fp0OHDhAubm5wvihbvdnTRyfiXiMZoyxv4uXrj1nBgYGyM3Nha6urrDtgw8+gJ2dHeLj43Hjxg0V1imj/y17eOWVVzBnzhyEhIRg/Pjx6NKlCyIiIgBALZcYaGlpQS6XIy8vD3379gVQ+/4RoHZJSmFhIZKSklBVVaXKTIHidp46dSo++OADWFpaCt9LSEiAs7MzMjMzQURqdXsrus+fP48hQ4YgMTER3t7eyMjIwJdffgl/f39s3rxZ7ZaxVVVVIT8/H66urgBq3/Ny//59mJqaIjs7G8nJyaiurlZx5V+ICM7OznjnnXdgb28vLFMzNDTEzZs3YWFhoeJCZYq+xYsXw83NDd9++y3efvtt9O/fHytXroSfnx+GDx+OO3fuqNX9GdCs8RngMZoxxp4VT3SeQVJSktIaaZlMBj09PcydOxdLlixBUFAQTExMUFRUhN69e+Pzzz/H6NGjcejQIZU179u3D5GRkfj+++9RWloKANDW1saLL74IAOjatStmzJiB/fv3C+9pUPWbXnfs2IE5c+Zgx44dwhvL27dvDx8fH+zcuRMAoKenB7lcDi8vL7z66qvYu3cv7ty5oxbNt2/fBlC7pr5Hjx4AgJSUFFhbW6OwsBCXL19GWFgYgoKC1Ka5pKREeABlamqKuLg4JCQkYOXKlTh06BAyMjLw7rvvIiYmBsePH1eL5lu3bgGovW/0798fK1asQHV1NbS1tbF9+3Z4enrC29sbcXFxqKmpUdn7Gx4fN4gI2traCAwMREBAAIC/HriWlZXB3t4e+vr6Kul83OPNWlpakMlk0NfXx6xZs/D111/j/Pnz+Prrr7Fr1y7s3bsXlZWVmDhxokofyGri+AzwGM0YY8+Val5I0mwFBQX08ssvk0QioRUrVhBR3eUDxcXFFB4eTuPHjxfWUcvlcurbty8tXLhQ9OUGWVlZ1KdPH+rWrRv5+fmRra0tjRgxQmkfxZKN4uJiGjFiBA0YMEDUxifdvHmT/P39ycHBgYKDg0kqlZKLiwv99NNPRES0c+dOcnBwoISEBCIiYfnX1atXSSKRCPuputnV1bVOy6FDh2jv3r1UVVVFVVVVVFxcTBKJhPbu3UtE4i6faaj51KlTRFR737G1tSUtLS06efKkcFxpaSk5OzvT+vXrRWttrNnFxYUyMjKIiGjfvn1kZGREnp6eZGtrS2ZmZnTmzBkqKioiiURCly9fFr25KeMG0V//9rNmzaKJEyc2uJ8YGmp+/P65fPlyOn36tNJxGRkZJJFIKDs7W9ReIs0cn4l4jGaMsZbAE51munTpEk2ePJnGjh1Ls2bNIiMjI7p16xYRKT84ffDgAXl4eFBiYiIR/TXA/+Mf/6CAgABRmzMyMsjHx4ciIiLowYMHVFpaSkePHiU9PT06evRonXYiooMHD5KFhQXt2LGDiIhOnjxJ165dE7U7Pj6e3N3dhdu3srKSXnjhBQoMDKS8vDy6desWBQUFka+vr3CMXC6nkpISkkqlFBUVJWpvY80jRoyo82DwcTKZjDp37kwLFy4UK1XQUPPw4cMpOzubqquracaMGaSvr08//vij0EtE1L9/f3rrrbfUpjkwMJCysrKIiOjnn3+m6Ohoio6OFh7MpqamUpcuXejs2bOi9jZ13Hicq6srbdu2TekyGtv/eXtas+I+UN+b+M+dO0fm5uYUExMjSmtTmxXUaXwm4jGaMcZaCi9dayZra2t4eHhg3rx5WLBgARwdHTFnzpw6++np6aGsrAxZWVnC14WFhXj48CFCQ0NFbSYi2NnZ4V//+hcMDQ1hbGwMd3d39OzZUzh9rWK5DP1vqceAAQMQFhaG999/H2PGjIGvry/y8/NFbU5LS4OVlRWMjIwgl8uhq6uLyMhIlJSUYOvWrbCyskJ4eDguX76MuXPnoqamBhKJBHl5edDV1RVOw6oOzbdu3cI333wjvDeEnlg2deTIEVhaWmLChAlq1RwTEwMtLS288847cHBwwIoVK1BUVAQtLS38+uuvkMvlatV8+/ZtxMTEoLq6Gh4eHpg0aRImT54MbW1tAEB6ejocHBzg5OQkanNTxw2F3Nxc3L9/H8OGDcO9e/cwbdo0dO3aFYWFhaK9J+NpzVpatb8+6ltad+rUKUilUgwdOlSUVgVNHJ8BHqMZY6zFqGqGpYkUz6hVVlYK2+Li4pReglc8y1ldXU2ffPIJSSQSCg0NpZkzZ5K1tTWNGTOG/vjjD1G7y8vL6f79+3W2d+zYkZKTk+s95u7duxQaGkoSiYQCAwNFPfOT4naeP38+ubm5EZHyaVQXLlxIPj4+wtmedu/eTR06dKC+ffvSlClTyNjYmN58802qqKgQ7dnvpjQPGDCAUlJShG03b96k4uJi+vjjj8ne3p7mzZtHlZWVatXs4+NDx48fJyKiU6dOUadOnUgqldK4cePIwsKCxo8fT2VlZaL0NrV5wIABwrPgRLXLZM6fP08ffvgh2dnZ0datW5UuS6zmpowbCrt376Z+/frR6tWrydjYmAYNGkQFBQWi9P7d5t9//52Ki4tpzZo15ODgQKtWraJHjx6p5e2sTuMzEY/RjDHWUnii8xQNDcKK7WVlZTR8+HDq169fvfutW7eO3nzzTQoODhbWKbe0pzUT1S6VkEqldOPGjTqfxXDlyhXy9/cnW1tbSkpKatHW6urqOm2Kv2dlZZGOjg59//33RPTXA5hLly6Ro6MjffHFF8IxKSkp9Omnn9Ibb7xBhw4dUtvmTZs2EVHtg8L169eTm5sbubi40MGDB9W6WbHv2bNn6ZtvvqH33ntPWPajrs0KCQkJ5O3tTc7Oziq5nR/X1HFj8uTJJJFIyNnZWXjfVkt5Hs03b96kNWvWUNeuXal79+4aczuLPT43tZtIfcbo+qj7GM0YYwo80XlCWVkZrVixgv773/82+ZhTp06Rrq6usFb60aNHVFpa2lKJdTSnWfEL6rPPPqP+/fs3+Gbox58RbwllZWU0YcIEmjRpklLX4+7cuUOjR4+mXr16CdsUv/ADAwMpNDS0RRuf9Dyba2pqKC8vr8UfELb121lxfGpqqsqbn1TfuKF4dSwqKoq2bNnScsH0/JoV79HJzs6mffv2tVwwPb9mMcdnouZ1q8sY/fDhQ8rNzaU7d+4QUf0nxFC3sYMxxp7EE53HfPbZZ2RgYECBgYEUFxdHVVVVTTquurqaZs+eTZ06daKCggJ6/fXXafHixfTnn3+2cHHzmxW/hIKCgmj27NlEVPvG3OnTp4v2zOCSJUtIV1eXdHR0qE+fPo0uFUlJSSEzMzNauXKlsK2qqor8/f1p3rx5ItTW4mZxPO9mMc6e1Zzmx9U3bnzwwQdNHneexfNsXrRoET18+LBlg+n5Nos1PhM1v1sdxuhVq1aRg4MD9erVi+zs7JSW3D5JXcYOxhirD090/ufEiRP0wgsv0K5du5S2N3X9cGZmJkkkEpJIJOTq6irKevq/21xVVUU9e/akEydOUHR0NFlYWJCjoyP99ttvLZlLsbGxZGNjQ927d6f09HTavHkzubm50Z07dxpsrqqqok2bNpFEIqGIiAg6duwYbdiwgezt7Rv95cvN3KyOzU8Se9zgZvHG52fpVtUYXV1dTeHh4eTu7k5Hjhyh9PR0Gj9+vPA+nIZaVfn/kDHGGtPmJzqKXzgzZ86k4OBgIiIqKiqijz76iPbt2/fUN3jK5XJKSUkhW1tbkkqloqw9ftbmpKQkkkgkZGJiQqamprRz584Wb87NzaVBgwbRRx99JGzLyckhiUQiPOho7Jf/J598Qr6+vuTm5kZSqZTi4+O5mZs1tlkV4wY3i9NM9OzdqhijiYgKCwvJzc2N9u/fL2zbs2cPDRo0SHj/TUOvkqri/yFjjD1Nm5/oENU+i+Xn50fbt2+nPXv2kJmZGfn7+5OTkxOZm5vTDz/80ODgXlNTQ6tXr6Zly5ZpTPMPP/xABgYGtHr1atF6Kyoq6izJycrKImdnZ2HtfH0e/xnkcrmoZ57iZnG0pWYFVYwb3CyeZ+1WxRhNRHTx4kWSSCR07NgxYduwYcNo4sSJFB8fL5xM4XGq/H/IGGNP007Vp7cWW3R0NK5evYp+/frB19cXBgYG0NHRgY2NDTZu3AgvLy9s27YNgYGBAIDRo0dj2bJl+Prrr9G9e3elyyIitGvXDu+9957wmRLq3gwAXl5euHXrFgwMDERtBmpvM8XnQfTo0QNlZWUoLy8HAMjl8jq3o+JrxXFubm7czM0a26w4RlXjBjerfzegujG6Y8eOGDduHMaMGYMhQ4bg8OHDcHFxgZWVFWbOnIktW7Zg/fr1Sv/XxPx/yBhjzSb61EpFLly4QF5eXuTo6Eh+fn5ka2tLr776Kt2+fZuIiGJiYsjCwoJsbGzo4sWLwnHFxcWkq6srvAlUzHP+t6ZmxZl7FBTPAgYGBlJISIhoffXhZnFwszi4WTya2F1f85AhQ4RmuVxOv/zyCwUEBAgnQyCq/dweIyMjiomJUVU6Y4w1W8s+zaVGDh48iHbt2iE/Px8JCQk4duwYsrKysGDBAty+fRuDBg2Cl5cXAMDc3BwAIJPJIJVKYW9vj7NnzwKAaJ9K3hqbr1y5AqD2mT8tLS3IZDJYWVmhoqICDx48EK2Rm7mZm7m5rXbX15ydnY0FCxaguLgYEokExsbGyM/PR3h4OACguroa5ubmsLe3R0ZGhsraGWOsudrERKempgYHDx5Er169oK+vj/bt28PV1RVffvkljh49iu+++w5dunTBW2+9hYqKCixduhQPHjyAtrY2Tp48CRMTE4wYMYKbn7E5ISEBQO3ESyaTQVtbGy4uLsjLy4ORkZGordzMzdzMzW2tu6nN7du3h0wmw7lz54SvT58+DV1dXQQHB6uknTHG/o5WP9EhIujo6MDExAQ3btwQtgHA2LFj4eXlhbi4OFy7dg1BQUHYuHEjtm3bhpdffhmTJk3CyJEjMXDgQHTu3Fk4jpub39y3b1/s27cPFy5cAPDXum5vb2/cu3cP2dnZonRyMzdzMze3xe6mNl+9ehV6enoICAjApEmTEBERgYiICAwfPhz9+vWDj4+P6O2MMfa3tfzqOPXwn//8hywsLOjs2bNERMKpMrOyskhXV5fS09OFfVNSUigyMpLeeustlX4OQGtuVqxZ/+6772jEiBF0/fp11QQTN4uFm8XBzeLRxO6nNaelpRERUUlJCUVERFBoaCi99tpr/Jk4jDGN1ComOvfu3VN60P84xekws7KyaPDgwRQaGip8T/HLp2fPnrRkyZKWD31MW28W4xPribiZmxvGzdzcGE3sfh7NixcvVjpOMRFijDFNpPFL1xYvXgwnJyds27YNJSUlwnaZTAYA0NHRAQD07dsX48ePR2pqKmJjYwHULie4efMmampq0KlTJ6XLpRZc8sXNyqck5WZu5mZuVpdmTe1+Xs1SqVSpVVdXt8WaGWOsxalogvXMZDIZTZs2jdzd3Sk5OZlKS0vp0aNHdfaLjY0lHR0d2rt3L927d4/mzJlDHTp0oE2bNtEvv/xCH3/8MXXr1o1yc3O5mZu5mZu5uY02a2q3JjYzxphYNHaiU1BQQJ6enpSVlUVEtZ8dU1RURH/88QcRET18+JBCQkLIysqKPv30U6VPqZ4+fTr16NGDnJycSCqV0uHDh7mZm7mZm7m5DTdrarcmNjPGmFg0bqKj+PDL5ORkcnV1pdLSUoqIiCAHBwfq06cPubm5UXx8PBERxcXF0Y0bN+ocK5fL6c8//6TMzExu5mZu5mZubsPNmtqtic2MMSa2dqpeOtcU6enpkEqlkEqlwodfXrp0Cba2ttiyZQt++eUX7Nq1C3p6etiwYQNWrFgBHR0dhISEKK2JVhxLRDAwMBA+bJObuZmbuZmb206zpnZrYjNjjKmUmLOq5kpJSSFHR0fq3Lkz2dnZ0cSJE4Vnni5evEjt2rUjExMTioyMFI65ePEihYaGUmhoaL3rlLmZm7mZm7m5bTZrarcmNjPGmDpQ27OuXbt2DYsWLcLrr7+O1NRUbNmyBXl5eZg/fz7OnTsHR0dHvPnmmygrK4O+vr5wnKOjI6ysrFBWViY8a8XN3MzN3MzNbbtZU7s1sZkxxtSGqmdaDTly5Ajp6elRYWGhsC05OZkGDx5M48aNIyKia9eukYmJCQUHB1NxcbGw3/Tp0yksLIybuZmbuZmbuVmgid2a2MwYY+pCbSc6u3fvJg8PDzp37pywTSaT0RdffEGdO3emhIQEIiKKiooiR0dHGjhwIMXGxtK8efPIzs6OEhMTuZmbuZmbuZmbNbpbE5sZY0xdqO1EJy8vj3R1dWn//v1K28+dO0ejRo2iyZMnC9uOHDlCwcHBNHjwYOrXrx+lpaWJXFuLm8XBzeLgZnFws3g0sVsTmxljTF1IiFr4I6afQUBAAMrLy3Ho0CEYGhoK22fNmoXCwkLs2bMHBgYGwva7d+/CwsJCFakCbhYHN4uDm8XBzeLRxG5NbGaMMXWgticjAIDVq1fj5MmT2LlzJ6qqqoTtVlZWKCgogK6urtL+6jCwc7M4uFkc3CwObhaPJnZrYjNjjKkDtf4cnT59+mD+/PlYvnw5tLW1MWHCBMjlcpw5cwZhYWFo10798rlZHNwsDm4WBzeLRxO7NbGZMcbUgVovXVN45513sHfvXkilUpSUlKBDhw6Ii4tDz549VZ3WIG4WBzeLg5vFwc3i0cRuTWxmjDFV0oiJTlVVFQoKCpCTk4P27dsjLCxM1UlPxc3i4GZxcLM4uFk8mtitic2MMaZKGjHRYYwxxhhjjLHmUOuTETDGGGOMMcbY38ETHcYYY4wxxlirwxMdxhhjjDHGWKvDEx3GGGOMMcZYq8MTHcYYY4wxxlirwxMdxhhjjDHGWKvDEx3GGGOMMcZYq8MTHcYYY4wxxlirwxMdxtqIZcuWwcPDQ9UZauPu3buwtrbG5cuXVZ3SoOjoaJiamja6jyr+XTdu3IiRI0eKep2MMcZYc/FEh7FWQCKRNPpnypQpmDt3LlJSUlSdCgDo0qULIiMj62yPjIxEly5dhH0a+5n8/PwAAOfPn8cbb7yBjh07QldXF46OjpgwYQIyMzMbbVi9ejVGjBghXF9rdfnyZUgkEuTk5NT5np+fH2bPnt3g1wpPTrjCw8ORkZGB9PT05x/MGGOMPSftVB3AGHt2v//+u/D32NhYLFmyBL/99puwTV9fH4aGhjA0NFRF3t+SkZEBmUwGADh16hTGjBmD3377DcbGxgCA9u3bIzMzE/7+/ujZsye++uoruLq64sGDBzhw4AAiIiJw4sSJei+7oqIC27ZtQ0JCgmg/T2uiq6uL119/HRs2bICvr6+qcxhjjLF68Ss6jLUCtra2wh8TExNIJJI6255c4jRlyhQEBQXhww8/hI2NDUxNTfHvf/8bjx49wrx582Bubo6OHTsiKipK6bquX7+O0NBQmJmZwcLCAqNGjWqR5V9WVlZCv7m5OQDA2tpa2GZmZoYpU6bA2dkZaWlpCAwMhJOTEzw8PLB06VIcOHCgwctOTExEu3bt8OKLLwrbjh8/DolEguTkZHh6ekJfXx+vvPIKSkpKkJiYCDc3NxgbG2PChAkoLy8XjvPz88O7776L2bNnw8zMDDY2NtiyZQsePnyIN954A0ZGRnByckJiYuJTf+bo6GhIpVJ06NABwcHBuHv3bp19PvroI9jY2MDIyAjTpk1DZWVlc27W52bkyJHYv38/KioqVHL9jDHG2NPwRIexNuyHH37AjRs3kJqainXr1mHZsmUYPnw4zMzMcPr0aUyfPh3Tp0/H1atXAQDl5eUYPHgwDA0NkZqaivT0dBgaGmLYsGGorq4WtT0nJwf5+fmIiIiAllbdoayx97akpqbC29u73u8tW7YMGzduxKlTp3D16lWMGzcOkZGR+Pbbb3H48GF8//332LBhg9IxMTExsLS0xJkzZ/Duu+9ixowZCAkJwYABA5CdnY2hQ4di4sSJShOkJ50+fRpTp07F22+/jZycHAwePBgrV65U2icuLg5Lly7FqlWrkJmZCTs7O2zatKmRW6nleHt7o6amBmfOnFHJ9TPGGGNPwxMdxtowc3NzfP7553BxccHUqVPh4uKC8vJyLFy4EM7Oznj//ffRvn17nDx5EgCwe/duaGlpYevWrejVqxfc3Nywfft2XLlyBcePHxe1vaioCADg6ura7GMvX74Me3v7er+3cuVKDBw4EJ6enpg2bRpOnDiBzZs3w9PTEy+99BLGjh2LY8eOKR3Tp08fLFq0SLjN9PX1YWlpifDwcDg7O2PJkiW4e/cucnNzG2xav349hg4digULFqB79+6YOXMmhg4dqrRPZGQkpk6din/+859wcXHBypUr4e7u3qSfecCAAcLyRcWftLS0Ovtt2rSpzn7Tp0+vs5+BgQFMTU3V+mQOjDHG2jae6DDWhvXo0UPp1RAbGxv06tVL+FpbWxsWFhYoKSkBAGRlZeH8+fMwMjISHgSbm5ujsrISFy5cELWdiADUnoihuSoqKqCnp1fv93r37i383cbGBh06dEDXrl2Vtiluj/qOUdxmj9+ONjY2ACAc16NHD+H2e+211wAAv/76q9JSOgB1vm7KPg2JjY1FTk6O0p/6XtUKCwurs9/y5cvrvUx9ff1GX6VijDHGVIlPRsBYG6ajo6P0tUQiqXebXC4HAMjlcnh5eeGbb76pc1lWVlZNvl5jY2OUlpbW2X7//n2YmJg06TK6d+8OoPbBf3NPr2xpaYk//vij3u89/vM/7fao75j6jlNMxhTHJSQkoKamBkDtZAH4a+LWUjp16oRu3bopbVNc9+NMTEzq7GdtbV3vZd67d69Z/+6MMcaYmPgVHcZYk/Xt2xdFRUWwtrZGt27dlP40dYIC1C43y8jIqLM9IyMDLi4uTboMDw8PuLu7Y+3atXUmHkDtpKkhnp6eKCgoaHLv89a5c2fhdnNwcAAAuLu746efflLa78mv3dzcnrqPWC5cuIDKykp4enqq5PoZY4yxp+GJDmOsycLCwmBpaYlRo0YhLS0Nly5dwokTJzBr1ixcu3atyZfzf//3f0hMTMTy5ctRUFCAgoICrFixAklJSYiIiGjSZUgkEmzfvh2FhYUYNGgQEhIScPHiReTm5mLVqlUYNWpUg8cOHToU+fn5Db6qowozZ85EUlIS1qxZg8LCQmzcuBFJSUlK+8yaNQtRUVGIiopCYWEhli5divz8fJX0pqWloWvXrnByclLJ9TPGGGNPwxMdxliTdejQAampqZBKpRg9ejTc3NwwdepUVFRUCJ9vozhNc2NvUvfx8UFycjKOHj0KX19f+Pr64siRI0hOTkb//v2b3NOvXz9kZmbCyckJ4eHhcHNzw8iRI5Gfn1/vB5Iq9OrVC97e3oiLi2vydbU0Hx8fbN26FRs2bICHhweOHDmCRYsWKe0TGhqKJUuWYP78+fDy8kJxcTFmzJihkt5du3YhPDxcJdfNGGOMNYWEWnphOGOsTYmOjsaqVatQUFBQ570r6iQhIQFz587F2bNn6z09NWvY2bNn4e/vj8LCwmYtWWSMMcbExCcjYIw9V0lJSfjwww/VepIDAAEBASgqKsL169fRqVMnVedolBs3bmDHjh08yWGMMabW+BUdxhhjjDHGWKvD6zUYY4wxxhhjrQ5PdBhjjDHGGGOtDk90GGOMMcYYY60OT3QYY4wxxhhjrQ5PdBhjjDHGGGOtDk90GGOMMcYYY60OT3QYY4wxxhhjrQ5PdBhjjDHGGGOtDk90GGOMMcYYY60OT3QYY4wxxhhjrc7/Ayxrn+vLJdYJAAAAAElFTkSuQmCC",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "result_df = result_df_hybrid.copy()\n",
    "\n",
    "for i in range(num_intervals):\n",
    "    alph=0.1+(i/10.)\n",
    "    plt.fill_between(result_df.index,# x\n",
    "                     result_df[str(percentile_list[i])+'_p_ac'].fillna(0), # lower\n",
    "                     result_df[str(percentile_list[len(percentile_list)-i-1])+'_p_ac'].fillna(0), # upper\n",
    "                     color='blue',\n",
    "                     alpha=alph,\n",
    "                     step='mid',\n",
    "                     label=(str(pi_list[i]*100)+'%'))\n",
    "    \n",
    "plt.plot(result_df.index,\n",
    "         result_df[str(percentile_list[int((len(percentile_list)-1)/2)])+'_p_ac'],\n",
    "         label='Hybrid expected',color='brown', drawstyle='steps-mid')\n",
    "# plt.plot(result_df.index,result_df['p_ac_all_act'],label='actual',color='black', drawstyle='steps-mid')\n",
    "plt.plot(result_df.index,result_df['p_ac_all_fcast'],label='DA fcast (HRRR coarse)',color='orange', drawstyle='steps-mid')\n",
    "plt.ylabel('AC Power (MW)')\n",
    "# plt.ylim([0,4500])\n",
    "plt.legend(bbox_to_anchor=(1.05, 1.0), loc='upper left')\n",
    "plt.gcf().autofmt_xdate()\n",
    "plt.xlabel('Time, UTC (mm-dd HH)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Preview the tabular results, which could be exported with `result_df.to_csv('results_file_name.csv')` or similar:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>tcc_std</th>\n",
       "      <th>max_min_norm</th>\n",
       "      <th>csi_fcast</th>\n",
       "      <th>zen</th>\n",
       "      <th>az</th>\n",
       "      <th>month</th>\n",
       "      <th>p_ac_all_fcast</th>\n",
       "      <th>p_ac_all_fcast_max</th>\n",
       "      <th>p_ac_all_fcast_min</th>\n",
       "      <th>p_ac_all_cs</th>\n",
       "      <th>...</th>\n",
       "      <th>0.005_p_ac</th>\n",
       "      <th>0.01_p_ac</th>\n",
       "      <th>0.025_p_ac</th>\n",
       "      <th>0.25_p_ac</th>\n",
       "      <th>0.5_p_ac</th>\n",
       "      <th>0.75_p_ac</th>\n",
       "      <th>0.975_p_ac</th>\n",
       "      <th>0.99_p_ac</th>\n",
       "      <th>0.995_p_ac</th>\n",
       "      <th>0.9975_p_ac</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>time_center_labeled</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>2024-06-06 18:00:00+00:00</th>\n",
       "      <td>13.887529</td>\n",
       "      <td>0.446797</td>\n",
       "      <td>0.870150</td>\n",
       "      <td>11.649303</td>\n",
       "      <td>212.139492</td>\n",
       "      <td>6.0</td>\n",
       "      <td>3221.029681</td>\n",
       "      <td>3667.311258</td>\n",
       "      <td>2013.404822</td>\n",
       "      <td>3701.694858</td>\n",
       "      <td>...</td>\n",
       "      <td>1061.858416</td>\n",
       "      <td>1487.555810</td>\n",
       "      <td>1778.045292</td>\n",
       "      <td>2983.331731</td>\n",
       "      <td>3268.192161</td>\n",
       "      <td>3432.998298</td>\n",
       "      <td>3601.465339</td>\n",
       "      <td>3626.209012</td>\n",
       "      <td>3649.251363</td>\n",
       "      <td>3678.076832</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2024-06-06 19:00:00+00:00</th>\n",
       "      <td>16.598369</td>\n",
       "      <td>0.458493</td>\n",
       "      <td>0.808130</td>\n",
       "      <td>21.589427</td>\n",
       "      <td>247.812447</td>\n",
       "      <td>6.0</td>\n",
       "      <td>2989.210341</td>\n",
       "      <td>3609.200097</td>\n",
       "      <td>1913.271904</td>\n",
       "      <td>3698.921102</td>\n",
       "      <td>...</td>\n",
       "      <td>906.110411</td>\n",
       "      <td>1071.142460</td>\n",
       "      <td>1609.218733</td>\n",
       "      <td>2738.077589</td>\n",
       "      <td>3126.532861</td>\n",
       "      <td>3318.119905</td>\n",
       "      <td>3546.131413</td>\n",
       "      <td>3586.561509</td>\n",
       "      <td>3624.854341</td>\n",
       "      <td>3642.646708</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2024-06-06 20:00:00+00:00</th>\n",
       "      <td>16.598369</td>\n",
       "      <td>0.508149</td>\n",
       "      <td>0.723995</td>\n",
       "      <td>33.772316</td>\n",
       "      <td>262.213148</td>\n",
       "      <td>6.0</td>\n",
       "      <td>2661.058360</td>\n",
       "      <td>3570.675807</td>\n",
       "      <td>1702.963161</td>\n",
       "      <td>3675.520056</td>\n",
       "      <td>...</td>\n",
       "      <td>595.725531</td>\n",
       "      <td>881.792957</td>\n",
       "      <td>1237.586023</td>\n",
       "      <td>2342.956508</td>\n",
       "      <td>2840.164572</td>\n",
       "      <td>3032.884559</td>\n",
       "      <td>3383.251950</td>\n",
       "      <td>3448.334786</td>\n",
       "      <td>3541.261568</td>\n",
       "      <td>3562.776366</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2024-06-06 21:00:00+00:00</th>\n",
       "      <td>16.598369</td>\n",
       "      <td>0.549489</td>\n",
       "      <td>0.702292</td>\n",
       "      <td>46.340846</td>\n",
       "      <td>271.235400</td>\n",
       "      <td>6.0</td>\n",
       "      <td>2516.617788</td>\n",
       "      <td>3457.253733</td>\n",
       "      <td>1488.196094</td>\n",
       "      <td>3583.435091</td>\n",
       "      <td>...</td>\n",
       "      <td>572.992246</td>\n",
       "      <td>730.489014</td>\n",
       "      <td>1155.484134</td>\n",
       "      <td>2154.910208</td>\n",
       "      <td>2670.120248</td>\n",
       "      <td>2902.766546</td>\n",
       "      <td>3260.280987</td>\n",
       "      <td>3322.461980</td>\n",
       "      <td>3427.355082</td>\n",
       "      <td>3440.054128</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2024-06-06 22:00:00+00:00</th>\n",
       "      <td>17.412338</td>\n",
       "      <td>0.673193</td>\n",
       "      <td>0.654215</td>\n",
       "      <td>58.880160</td>\n",
       "      <td>278.571800</td>\n",
       "      <td>6.0</td>\n",
       "      <td>2247.652475</td>\n",
       "      <td>3209.861779</td>\n",
       "      <td>897.006441</td>\n",
       "      <td>3435.648003</td>\n",
       "      <td>...</td>\n",
       "      <td>197.829622</td>\n",
       "      <td>379.747230</td>\n",
       "      <td>1016.378641</td>\n",
       "      <td>1860.679711</td>\n",
       "      <td>2323.433875</td>\n",
       "      <td>2643.200406</td>\n",
       "      <td>3048.085604</td>\n",
       "      <td>3106.795581</td>\n",
       "      <td>3270.737793</td>\n",
       "      <td>3453.222890</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>5 rows × 34 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "                             tcc_std  max_min_norm  csi_fcast        zen  \\\n",
       "time_center_labeled                                                        \n",
       "2024-06-06 18:00:00+00:00  13.887529      0.446797   0.870150  11.649303   \n",
       "2024-06-06 19:00:00+00:00  16.598369      0.458493   0.808130  21.589427   \n",
       "2024-06-06 20:00:00+00:00  16.598369      0.508149   0.723995  33.772316   \n",
       "2024-06-06 21:00:00+00:00  16.598369      0.549489   0.702292  46.340846   \n",
       "2024-06-06 22:00:00+00:00  17.412338      0.673193   0.654215  58.880160   \n",
       "\n",
       "                                   az  month  p_ac_all_fcast  \\\n",
       "time_center_labeled                                            \n",
       "2024-06-06 18:00:00+00:00  212.139492    6.0     3221.029681   \n",
       "2024-06-06 19:00:00+00:00  247.812447    6.0     2989.210341   \n",
       "2024-06-06 20:00:00+00:00  262.213148    6.0     2661.058360   \n",
       "2024-06-06 21:00:00+00:00  271.235400    6.0     2516.617788   \n",
       "2024-06-06 22:00:00+00:00  278.571800    6.0     2247.652475   \n",
       "\n",
       "                           p_ac_all_fcast_max  p_ac_all_fcast_min  \\\n",
       "time_center_labeled                                                 \n",
       "2024-06-06 18:00:00+00:00         3667.311258         2013.404822   \n",
       "2024-06-06 19:00:00+00:00         3609.200097         1913.271904   \n",
       "2024-06-06 20:00:00+00:00         3570.675807         1702.963161   \n",
       "2024-06-06 21:00:00+00:00         3457.253733         1488.196094   \n",
       "2024-06-06 22:00:00+00:00         3209.861779          897.006441   \n",
       "\n",
       "                           p_ac_all_cs  ...   0.005_p_ac    0.01_p_ac  \\\n",
       "time_center_labeled                     ...                             \n",
       "2024-06-06 18:00:00+00:00  3701.694858  ...  1061.858416  1487.555810   \n",
       "2024-06-06 19:00:00+00:00  3698.921102  ...   906.110411  1071.142460   \n",
       "2024-06-06 20:00:00+00:00  3675.520056  ...   595.725531   881.792957   \n",
       "2024-06-06 21:00:00+00:00  3583.435091  ...   572.992246   730.489014   \n",
       "2024-06-06 22:00:00+00:00  3435.648003  ...   197.829622   379.747230   \n",
       "\n",
       "                            0.025_p_ac    0.25_p_ac     0.5_p_ac    0.75_p_ac  \\\n",
       "time_center_labeled                                                             \n",
       "2024-06-06 18:00:00+00:00  1778.045292  2983.331731  3268.192161  3432.998298   \n",
       "2024-06-06 19:00:00+00:00  1609.218733  2738.077589  3126.532861  3318.119905   \n",
       "2024-06-06 20:00:00+00:00  1237.586023  2342.956508  2840.164572  3032.884559   \n",
       "2024-06-06 21:00:00+00:00  1155.484134  2154.910208  2670.120248  2902.766546   \n",
       "2024-06-06 22:00:00+00:00  1016.378641  1860.679711  2323.433875  2643.200406   \n",
       "\n",
       "                            0.975_p_ac    0.99_p_ac   0.995_p_ac  0.9975_p_ac  \n",
       "time_center_labeled                                                            \n",
       "2024-06-06 18:00:00+00:00  3601.465339  3626.209012  3649.251363  3678.076832  \n",
       "2024-06-06 19:00:00+00:00  3546.131413  3586.561509  3624.854341  3642.646708  \n",
       "2024-06-06 20:00:00+00:00  3383.251950  3448.334786  3541.261568  3562.776366  \n",
       "2024-06-06 21:00:00+00:00  3260.280987  3322.461980  3427.355082  3440.054128  \n",
       "2024-06-06 22:00:00+00:00  3048.085604  3106.795581  3270.737793  3453.222890  \n",
       "\n",
       "[5 rows x 34 columns]"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result_df.head()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For our next notebook, let's revisit our early forecast example time range in April, 2021: [11_revisting_early_forecast.ipynb](11_revisting_early_forecast.ipynb)."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}