the-final-kernel.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Collecting meteocalc\r\n",
      "  Downloading https://files.pythonhosted.org/packages/6c/f7/95473a929f0a02547461fa3698b7f8082ff40445ba5e21601f5d9a5e48ec/meteocalc-1.1.0.tar.gz\r\n",
      "Building wheels for collected packages: meteocalc\r\n",
      "  Building wheel for meteocalc (setup.py) ... \u001b[?25l-\b \b\\\b \bdone\r\n",
      "\u001b[?25h  Created wheel for meteocalc: filename=meteocalc-1.1.0-cp36-none-any.whl size=8196 sha256=52e4e258dbb85797d3667fa4c977e14403b403b48f5b25ead51668f906f68b4e\r\n",
      "  Stored in directory: /tmp/.cache/pip/wheels/9e/34/13/83d36ecc28837e3c2a5b696542e697538e7c1025382f4ded55\r\n",
      "Successfully built meteocalc\r\n",
      "Installing collected packages: meteocalc\r\n",
      "Successfully installed meteocalc-1.1.0\r\n"
     ]
    }
   ],
   "source": [
    "!pip install meteocalc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Requirement already satisfied: holidays in /opt/conda/lib/python3.6/site-packages (0.9.11)\r\n",
      "Requirement already satisfied: python-dateutil in /opt/conda/lib/python3.6/site-packages (from holidays) (2.8.0)\r\n",
      "Requirement already satisfied: six in /opt/conda/lib/python3.6/site-packages (from holidays) (1.12.0)\r\n"
     ]
    }
   ],
   "source": [
    "!pip install holidays"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "_cell_guid": "d665b504-3089-4cda-9e2c-adb138822fdd",
    "_uuid": "d77b2edf-0bce-4f96-a789-d52769133c05"
   },
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.6/site-packages/sklearn/externals/joblib/__init__.py:15: DeprecationWarning: sklearn.externals.joblib is deprecated in 0.21 and will be removed in 0.23. Please import this functionality directly from joblib, which can be installed with: pip install joblib. If this warning is raised when loading pickled models, you may need to re-serialize those models with scikit-learn 0.21+.\n",
      "  warnings.warn(msg, category=DeprecationWarning)\n"
     ]
    }
   ],
   "source": [
    "import gc\n",
    "import os\n",
    "import random\n",
    "\n",
    "import holidays\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "import lightgbm as lgb\n",
    "from catboost import Pool, CatBoostRegressor\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from tqdm import tqdm_notebook as tqdm\n",
    "\n",
    "from sklearn.metrics import mean_squared_error, mean_squared_log_error\n",
    "from sklearn.preprocessing import LabelEncoder\n",
    "from sklearn.model_selection import RandomizedSearchCV\n",
    "from sklearn.model_selection import KFold, GroupKFold,StratifiedKFold\n",
    "from sklearn.externals import joblib\n",
    "\n",
    "# matplotlib and seaborn for plotting\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sns\n",
    "\n",
    "from pandas.plotting import register_matplotlib_converters\n",
    "register_matplotlib_converters()\n",
    "\n",
    "import calendar\n",
    "import datetime\n",
    "\n",
    "path_data = \"/kaggle/input/ashrae-energy-prediction/\"\n",
    "path_train = path_data + \"train.csv\"\n",
    "path_test = path_data + \"test.csv\"\n",
    "path_building = path_data + \"building_metadata.csv\"\n",
    "path_weather_train = path_data + \"weather_train.csv\"\n",
    "path_weather_test = path_data + \"weather_test.csv\"\n",
    "path_dataset_region_sites = \"/kaggle/input/ashrae-cities-and-timezones/\"\n",
    "\n",
    "myfavouritenumber = 666\n",
    "seed = myfavouritenumber\n",
    "random.seed(seed)\n",
    "\n",
    "import xgboost as xgb\n",
    "from xgboost import plot_importance, plot_tree\n",
    "\n",
    "from meteocalc import Temp, dew_point, heat_index, wind_chill, feels_like\n",
    "\n",
    "\n",
    "DO_TEST = False\n",
    "DO_TRAIN = True"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['ashrae-cities-and-timezones', 'lgb-models-ashrae', 'ashrae-leak-data-station', 'ashrae-energy-prediction', 'outliers-energydemand']\n"
     ]
    }
   ],
   "source": [
    "print(os.listdir('/kaggle/input/'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "3a3bc974-e6c0-49a0-8f1d-e1cd8f9ca518",
    "_uuid": "58b023be-56ab-4549-bbbf-42ce385d70f0"
   },
   "source": [
    "## Reading train data\n",
    "\n",
    "Reading train data along with building and weather metadata."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "_cell_guid": "f01f65b1-a6ce-4bda-97bd-ec7141e64f1d",
    "_uuid": "901bac8f-44a3-4d57-9deb-12770f6a90d5"
   },
   "outputs": [],
   "source": [
    "#https://www.mytecbits.com/internet/python/week-number-of-month\n",
    "def get_week_of_month(date_value):\n",
    "    # isocalendar -- > year, week_del_año, dia_de_la_semana\n",
    "    # empeizan ene le 1\n",
    "    # .dt.date.values[0]\n",
    "     return (date_value.isocalendar()[1] - date_value.replace(day=1).isocalendar()[1] + 1)\n",
    "# 11 dummies para los meses\n",
    "def set_dummies(df):\n",
    "    meses_dummie = 11\n",
    "    for mes_dummie in np.arange(1, meses_dummie + 1):\n",
    "        df['mes' + str(mes_dummie)] = np.full(df.shape[0], 0)\n",
    "        df.loc[df['timestamp'].dt.month == mes_dummie, 'mes' + str(mes_dummie)] = 1\n",
    "\n",
    "# estacionalidades\n",
    "def get_summer(df):\n",
    "    summer = [6, 7, 8, 9]\n",
    "    df['is_summer'] = np.zeros(df.shape[0], dtype=int)\n",
    "    df.loc[df.timestamp.dt.month.isin(summer), 'is_summer'] = 1\n",
    "def get_winter(df):\n",
    "    winter = [1, 2, 11, 12]\n",
    "    df['is_winter'] = np.zeros(df.shape[0], dtype=int)\n",
    "    df.loc[df.timestamp.dt.month.isin(winter), 'is_winter'] = 1\n",
    "\n",
    "# centrarse en fines de semana\n",
    "def is_weekend(df):\n",
    "    weekend = [5,6] # empieza en 0\n",
    "    df['is_weekend'] = np.zeros(df.shape[0], dtype=int)\n",
    "    df.loc[df.timestamp.dt.dayofweek.isin(weekend), 'is_weekend'] = 1\n",
    "\n",
    "#RMSLE\n",
    "def get_rmsle(y_true, y_pred):\n",
    "    return np.sqrt(mean_squared_log_error( np.clip(y_true, 0, np.max(y_true)),\n",
    "                                                   np.clip(y_pred, 0,np.max(y_pred) )))\n",
    "# importancia de las carac\n",
    "def plot_fi(model):\n",
    "    importance_df = pd.DataFrame(model.feature_importance(),\n",
    "                                 index=model.feature_name(),\n",
    "                                 columns=['importance']).sort_values('importance')\n",
    "    fig, ax = plt.subplots(figsize=(8, 8))\n",
    "    importance_df.plot.barh(ax=ax)\n",
    "    fig.show()\n",
    "    \n",
    "# importancia para CATBOOST\n",
    "def plot_fi_cat(model):\n",
    "    importance_df = pd.DataFrame(model.get_feature_importance(),\n",
    "                                 index=model.feature_names_,\n",
    "                                 columns=['importance']).sort_values('importance')\n",
    "    fig, ax = plt.subplots(figsize=(8, 8))\n",
    "    importance_df.plot.barh(ax=ax)\n",
    "    fig.show()\n",
    "# TEMPERATURA caracs\n",
    "def humidity(df):\n",
    "    df['humidity'] = 100 - 5 * (df['air_temperature'] - df['dew_temperature'])\n",
    "\n",
    "# feels like\n",
    "def feels_like_df(df):\n",
    "    #lambda row: my_test(row['a'], row['c'])\n",
    "    df['feels'] = df.apply(lambda row: feels_like(temperature=row['air_temperature'],humidity=row['humidity'],\n",
    "                                                      wind_speed=row['wind_speed']), axis=1)\n",
    "    df['feels'] = df['feels'].astype(np.float16)\n",
    "    df['feels'] = np.round(df['feels'], 2)\n",
    "    \n",
    "#https://www.kaggle.com/caesarlupum/ashrae-ligthgbm-simple-fe/comments\n",
    "def beaufort_scale(df):\n",
    "    beaufort = [(0, 0, 0.3), (1, 0.3, 1.6), (2, 1.6, 3.4), (3, 3.4, 5.5), (4, 5.5, 8), (5, 8, 10.8), (6, 10.8, 13.9), \n",
    "          (7, 13.9, 17.2), (8, 17.2, 20.8), (9, 20.8, 24.5), (10, 24.5, 28.5), (11, 28.5, 33), (12, 33, 200)]\n",
    "\n",
    "    for item in beaufort:\n",
    "        df.loc[(df['wind_speed']>=item[1])\n",
    "               & (df['wind_speed']<item[2]),\n",
    "               'beaufort_scale'] = item[0]\n",
    "        \n",
    "    df['beaufort_scale'] = df['beaufort_scale'].fillna(-999)\n",
    "    df['beaufort_scale'] = df['beaufort_scale'].astype(np.int16) \n",
    "        \n",
    "# Convert wind direction into categorical feature. We can split 360 degrees into 16-wind compass rose. \n",
    "# See this: https://en.wikipedia.org/wiki/Points_of_the_compass#16-wind_compass_rose\n",
    "def degToCompass(num):\n",
    "    val=int((num/22.5))\n",
    "    arr=[i for i in range(0,16)]\n",
    "    return arr[(val % 16)]\n",
    "\n",
    "def wind_direction_to_categorical(df):\n",
    "    df['wind_direction'] = df['wind_direction'].apply(degToCompass)\n",
    "    "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**VACACIONES**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# VACACIONES\n",
    "def get_holidays(site):\n",
    "    \n",
    "    df_sites_regions = pd.read_csv(path_dataset_region_sites + 'site_ids.csv')\n",
    "    #print(df_sites_regions[['site_id', 'location', 'country']])\n",
    "    \n",
    "    not_USA_holidays = [1, 5, 7, 11, 12]\n",
    "    \n",
    "    # es americano\n",
    "    if site not in not_USA_holidays:\n",
    "        # limpiamos\n",
    "        state = df_sites_regions[df_sites_regions.site_id == site]\n",
    "        state = state.location.values[0]\n",
    "        state = state.replace(\" \", \"\")\n",
    "        state = state.split(',')[1]\n",
    "\n",
    "        list_holidays = holidays.UnitedStates(state=state)\n",
    "        \n",
    "    # UK\n",
    "    elif site in [1, 5]:\n",
    "        list_holidays = holidays.UnitedKingdom()\n",
    "        \n",
    "    # CANADA\n",
    "    elif site in [7, 11]:\n",
    "        if site == 7:\n",
    "            list_holidays = holidays.Canada(prov='QC')\n",
    "        else:\n",
    "            list_holidays = holidays.Canada(prov='ON')\n",
    "    \n",
    "    # IRLANDA\n",
    "    else:\n",
    "        list_holidays = holidays.Ireland()\n",
    "        \n",
    "    return list_holidays['2016-01-01':'2017-12-31'] +\\\n",
    "    list_holidays['2018-01-01':'2018-12-31'] #esto por el bug que tiene la libreria"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "#https://www.kaggle.com/purist1024/ashrae-simple-data-cleanup-lb-1-08-no-leaks\n",
    "# 3000, 7500\n",
    "def make_is_bad_zero(Xy_subset, min_interval=48,\n",
    "                     summer_start=datetime.datetime(2016,4,1),\n",
    "                     summer_end=datetime.datetime(2016,10,1)):\n",
    "    \"\"\"Helper routine for 'find_bad_zeros'.\n",
    "    \n",
    "    This operates upon a single dataframe produced by 'groupby'. We expect an \n",
    "    additional column 'meter_id' which is a duplicate of 'meter' because groupby \n",
    "    eliminates the original one.\"\"\"\n",
    "    meter = Xy_subset.meter_id.iloc[0]\n",
    "    is_zero = Xy_subset.meter_reading == 0\n",
    "    if meter == 0:\n",
    "        # Electrical meters should never be zero. Keep all zero-readings in this table so that\n",
    "        # they will all be dropped in the train set.\n",
    "        return is_zero\n",
    "\n",
    "    transitions = (is_zero != is_zero.shift(1))\n",
    "    all_sequence_ids = transitions.cumsum()\n",
    "    ids = all_sequence_ids[is_zero].rename(\"ids\")\n",
    "    if meter in [2, 3]:\n",
    "        # It's normal for steam and hotwater to be turned off during the summer\n",
    "        # ++ .timestamp()\n",
    "        keep = set(ids[(Xy_subset.timestamp < summer_start) |\n",
    "                       (Xy_subset.timestamp > summer_end)].unique())\n",
    "        is_bad = ids.isin(keep) & (ids.map(ids.value_counts()) >= min_interval)\n",
    "    elif meter == 1:\n",
    "        time_ids = ids.to_frame().join(Xy_subset.timestamp).set_index(\"timestamp\").ids\n",
    "        is_bad = ids.map(ids.value_counts()) >= min_interval\n",
    "\n",
    "        # Cold water may be turned off during the winter\n",
    "        jan_id = time_ids.get(0, False)\n",
    "        dec_id = time_ids.get(8283, False)\n",
    "        if (jan_id and dec_id and jan_id == time_ids.get(500, False) and\n",
    "                dec_id == time_ids.get(8783, False)):\n",
    "            is_bad = is_bad & (~(ids.isin(set([jan_id, dec_id]))))\n",
    "    else:\n",
    "        raise Exception(f\"Unexpected meter type: {meter}\")\n",
    "\n",
    "    result = is_zero.copy()\n",
    "    result.update(is_bad)\n",
    "    return result\n",
    "\n",
    "def find_bad_zeros(X, y):\n",
    "    \"\"\"Returns an Index object containing only the rows which should be deleted.\"\"\"\n",
    "    Xy = X.assign(meter_reading=y, meter_id=X.meter)\n",
    "    is_bad_zero = Xy.groupby([\"building_id\", \"meter\"]).apply(make_is_bad_zero)\n",
    "    return is_bad_zero[is_bad_zero].index.droplevel([0, 1])\n",
    "\n",
    "def find_bad_sitezero(X):\n",
    "    \"\"\"Returns indices of bad rows from the early days of Site 0 (UCF).\"\"\"\n",
    "    veinte_del_cinco = datetime.datetime(2016,5,20)\n",
    "    return X[(X.timestamp < veinte_del_cinco) & (X.site_id == 0) & (X.meter == 0)].index\n",
    "\n",
    "def find_bad_building1099(X, y):\n",
    "    \"\"\"Returns indices of bad rows (with absurdly high readings) from building 1099.\"\"\"\n",
    "    return X[(X.building_id == 1099) & (X.meter == 2) & (y > 3e4)].index\n",
    "\n",
    "def find_bad_rows(X, y):\n",
    "    return find_bad_zeros(X, y).union(find_bad_sitezero(X)).union(find_bad_building1099(X, y))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "#https://www.kaggle.com/rohanrao/ashrae-divide-and-conquer\n",
    "def create_lag_features(df, window):\n",
    "    \"\"\"\n",
    "    Creating lag-based features looking back in time.\n",
    "    \"\"\"\n",
    "    \n",
    "    feature_cols = [\"air_temperature\", \"humidity\"]\n",
    "    df_site = df.groupby(\"site_id\")\n",
    "    \n",
    "    df_rolled = df_site[feature_cols].rolling(window=window, min_periods=0)\n",
    "    \n",
    "    df_mean = df_rolled.mean().reset_index().astype(np.float16)\n",
    "    df_median = df_rolled.median().reset_index().astype(np.float16)\n",
    "    df_min = df_rolled.min().reset_index().astype(np.float16)\n",
    "    df_max = df_rolled.max().reset_index().astype(np.float16)\n",
    "    df_std = df_rolled.std().reset_index().astype(np.float16)\n",
    "    df_skew = df_rolled.skew().reset_index().astype(np.float16)\n",
    "    \n",
    "    for feature in feature_cols:\n",
    "        df[f\"{feature}_mean_lag{window}\"] = df_mean[feature]\n",
    "        #df[f\"{feature}_median_lag{window}\"] = df_median[feature]\n",
    "        df[f\"{feature}_min_lag{window}\"] = df_min[feature]\n",
    "        df[f\"{feature}_max_lag{window}\"] = df_max[feature]\n",
    "        #df[f\"{feature}_std_lag{window}\"] = df_std[feature]\n",
    "        #df[f\"{feature}_skew_lag{window}\"] = df_skew[feature]\n",
    "        \n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "3b823785-5d08-4f12-8037-6f1473d7f624",
    "_uuid": "8d4f21ef-b0ba-4578-a8b9-a950e5339c9e"
   },
   "source": [
    "Para ajustar la temperatura\n",
    "# # ---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "_cell_guid": "6546cbfe-d4ff-49c7-ba6c-25afe5f12f1a",
    "_uuid": "4c202726-7a35-4ac9-957f-619716ff9cbe"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 768 ms, sys: 100 ms, total: 868 ms\n",
      "Wall time: 1.18 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "weather_dtypes = {\n",
    "    'site_id': np.uint8,\n",
    "    'air_temperature': np.float32,\n",
    "    'cloud_coverage': np.float32,\n",
    "    'dew_temperature': np.float32,\n",
    "    'precip_depth_1_hr': np.float32,\n",
    "    'sea_level_pressure': np.float32,\n",
    "    'wind_direction': np.float32,\n",
    "    'wind_speed': np.float32,\n",
    "}\n",
    "\n",
    "weather_train = pd.read_csv(path_weather_train,dtype=weather_dtypes,\n",
    "    parse_dates=['timestamp'])\n",
    "weather_test = pd.read_csv(path_weather_test,dtype=weather_dtypes,\n",
    "    parse_dates=['timestamp'])\n",
    "\n",
    "weather = pd.concat([weather_train,weather_test],ignore_index=True)\n",
    "#del weather_train, weather_test\n",
    "weather_key = ['site_id', 'timestamp']\n",
    "temp_skeleton = weather[weather_key + ['air_temperature']].drop_duplicates(subset=weather_key).sort_values(by=weather_key).copy()\n",
    "del weather, weather_train, weather_test"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "_cell_guid": "7a01e052-46f3-4234-a9e1-8c99ab1351f9",
    "_uuid": "92f7a0bb-03a8-4f47-80dd-4d05c65c7c65"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 620 ms, sys: 16 ms, total: 636 ms\n",
      "Wall time: 585 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "# calculate ranks of hourly temperatures within date/site_id chunks\n",
    "temp_skeleton['temp_rank'] = temp_skeleton.groupby(['site_id', temp_skeleton.timestamp.dt.date])['air_temperature'].rank('average')\n",
    "\n",
    "# create a dataframe of site_ids (0-16) x mean hour rank of temperature within day (0-23)\n",
    "df_2d = temp_skeleton.groupby(['site_id', temp_skeleton.timestamp.dt.hour])['temp_rank'].mean().unstack(level=1)\n",
    "\n",
    "# Subtract the columnID of temperature peak by 14, getting the timestamp alignment gap.\n",
    "# La hora pico de temperatura suele estar sobre las 8-9 de la tarde en el dataset (deberia estar), por lo que se le resta 14\n",
    "#(a las dos suele ser la hora pico), dando un array de offsets/desplazamientos --> cuantas veces hay que moverse\n",
    "site_ids_offsets = pd.Series(df_2d.values.argmax(axis=1) - 14)\n",
    "site_ids_offsets.index.name = 'site_id'\n",
    "\n",
    "del df_2d\n",
    "gc.collect()\n",
    "\n",
    "def timestamp_align(df):\n",
    "    df['offset'] = df.site_id.map(site_ids_offsets)\n",
    "    df['timestamp_aligned'] = (df.timestamp - pd.to_timedelta(df.offset, unit='H'))\n",
    "    df['timestamp'] = df['timestamp_aligned']\n",
    "    del df['timestamp_aligned']\n",
    "    return df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "20fcdb98-f933-416b-8cd9-ae927f51e604",
    "_uuid": "fd89ed16-8ad5-4c66-bb56-1493800a64b4"
   },
   "source": [
    "El resto\n",
    " ---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "_cell_guid": "20b2f8fc-b3d2-4806-b72b-29267eced38a",
    "_uuid": "97730aab-8af2-4579-8e96-5ed1492194a7"
   },
   "outputs": [],
   "source": [
    "CATEGORY_COLS = [  \"building_id\",\"site_id\", \"meter\", \"primary_use\",\n",
    "                'is_weekend',\"hour\",\"weekday\" ,\"year\"] #, 'weekofmonth',\n",
    "# ,'is_holiday'\"hour\", \"weekday\", , 'year_built'\n",
    "#,'beaufort_scale', 'wind_direction''is_summer',\n",
    "\n",
    "FEATURE_COLS = ['square_feet', 'floor_count'] + [\n",
    "    'air_temperature', 'cloud_coverage',\n",
    "    'dew_temperature',  'wind_speed']  #+windspeed --'precip_depth_1_hr',\n",
    "\n",
    "USAR_LEAKED_DATA = False"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "_cell_guid": "5fd71733-7979-43fd-bc76-8c2615ab5f0f",
    "_kg_hide-input": true,
    "_kg_hide-output": true,
    "_uuid": "463a9374-2e79-49a4-907b-6775f9ca7d43"
   },
   "outputs": [],
   "source": [
    "## Memory optimization\n",
    "\n",
    "# Original code from https://www.kaggle.com/gemartin/load-data-reduce-memory-usage by @gemartin\n",
    "# Modified to support timestamp type, categorical type\n",
    "# Modified to add option to use float16\n",
    "\n",
    "from pandas.api.types import is_datetime64_any_dtype as is_datetime\n",
    "from pandas.api.types import is_categorical_dtype\n",
    "\n",
    "def reduce_mem_usage(df, use_float16=False):\n",
    "    \"\"\"\n",
    "    Iterate through all the columns of a dataframe and modify the data type to reduce memory usage.        \n",
    "    \"\"\"\n",
    "    \n",
    "    start_mem = df.memory_usage().sum() / 1024**2\n",
    "    print(\"Memory usage of dataframe is {:.2f} MB\".format(start_mem))\n",
    "    \n",
    "    for col in df.columns:\n",
    "        if is_datetime(df[col]) or is_categorical_dtype(df[col]):\n",
    "            continue\n",
    "        col_type = df[col].dtype\n",
    "        \n",
    "        if col_type != object:\n",
    "            c_min = df[col].min()\n",
    "            c_max = df[col].max()\n",
    "            if str(col_type)[:3] == \"int\":\n",
    "                if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:\n",
    "                    df[col] = df[col].astype(np.int8)\n",
    "                elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:\n",
    "                    df[col] = df[col].astype(np.int16)\n",
    "                elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:\n",
    "                    df[col] = df[col].astype(np.int32)\n",
    "                elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:\n",
    "                    df[col] = df[col].astype(np.int64)  \n",
    "            else:\n",
    "                if use_float16 and c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:\n",
    "                    df[col] = df[col].astype(np.float16)\n",
    "                elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:\n",
    "                    df[col] = df[col].astype(np.float32)\n",
    "                else:\n",
    "                    df[col] = df[col].astype(np.float64)\n",
    "        else:\n",
    "            df[col] = df[col].astype(\"category\")\n",
    "\n",
    "    end_mem = df.memory_usage().sum() / 1024**2\n",
    "    print(\"Memory usage after optimization is: {:.2f} MB\".format(end_mem))\n",
    "    print(\"Decreased by {:.1f}%\".format(100 * (start_mem - end_mem) / start_mem))\n",
    "    \n",
    "    return df"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Original code from https://www.kaggle.com/aitude/ashrae-missing-weather-data-handling by @aitude\n",
    "def fill_weather_dataset(weather_df):\n",
    "    \n",
    "    # Find Missing Dates\n",
    "    time_format = \"%Y-%m-%d %H:%M:%S\"\n",
    "    start_date = datetime.datetime.strptime(weather_df['timestamp'].min(),time_format)\n",
    "    end_date = datetime.datetime.strptime(weather_df['timestamp'].max(),time_format)\n",
    "    total_hours = int(((end_date - start_date).total_seconds() + 3600) / 3600)\n",
    "    hours_list = [(end_date - datetime.timedelta(hours=x)).strftime(time_format) for x in range(total_hours)]\n",
    "\n",
    "    missing_hours = []\n",
    "    for site_id in range(16):\n",
    "        site_hours = np.array(weather_df[weather_df['site_id'] == site_id]['timestamp'])\n",
    "        new_rows = pd.DataFrame(np.setdiff1d(hours_list,site_hours),columns=['timestamp'])\n",
    "        new_rows['site_id'] = site_id\n",
    "        weather_df = pd.concat([weather_df,new_rows])\n",
    "\n",
    "        weather_df = weather_df.reset_index(drop=True)     \n",
    "    \n",
    "        # Add new Features\n",
    "    weather_df[\"datetime\"] = pd.to_datetime(weather_df[\"timestamp\"])\n",
    "    weather_df[\"day\"] = weather_df[\"datetime\"].dt.day\n",
    "    weather_df[\"week\"] = weather_df[\"datetime\"].dt.week\n",
    "    weather_df[\"month\"] = weather_df[\"datetime\"].dt.month\n",
    "    \n",
    "    # Reset Index for Fast Update\n",
    "    weather_df = weather_df.set_index(['site_id','day','month'])\n",
    "\n",
    "    air_temperature_filler = pd.DataFrame(weather_df.groupby(['site_id','day','month'])['air_temperature'].mean(),columns=[\"air_temperature\"])\n",
    "    weather_df.update(air_temperature_filler,overwrite=False)\n",
    "\n",
    "    # Step 1\n",
    "    cloud_coverage_filler = weather_df.groupby(['site_id','day','month'])['cloud_coverage'].mean()\n",
    "    # Step 2\n",
    "    cloud_coverage_filler = pd.DataFrame(cloud_coverage_filler.fillna(method='ffill'),columns=[\"cloud_coverage\"])\n",
    "\n",
    "    weather_df.update(cloud_coverage_filler,overwrite=False)\n",
    "\n",
    "    due_temperature_filler = pd.DataFrame(weather_df.groupby(['site_id','day','month'])['dew_temperature'].mean(),columns=[\"dew_temperature\"])\n",
    "    weather_df.update(due_temperature_filler,overwrite=False)\n",
    "\n",
    "    # Step 1\n",
    "    sea_level_filler = weather_df.groupby(['site_id','day','month'])['sea_level_pressure'].mean()\n",
    "    # Step 2\n",
    "    sea_level_filler = pd.DataFrame(sea_level_filler.fillna(method='ffill'),columns=['sea_level_pressure'])\n",
    "\n",
    "    weather_df.update(sea_level_filler,overwrite=False)\n",
    "\n",
    "    wind_direction_filler =  pd.DataFrame(weather_df.groupby(['site_id','day','month'])['wind_direction'].mean(),columns=['wind_direction'])\n",
    "    weather_df.update(wind_direction_filler,overwrite=False)\n",
    "\n",
    "    wind_speed_filler =  pd.DataFrame(weather_df.groupby(['site_id','day','month'])['wind_speed'].mean(),columns=['wind_speed'])\n",
    "    weather_df.update(wind_speed_filler,overwrite=False)\n",
    "\n",
    "    # Step 1\n",
    "    precip_depth_filler = weather_df.groupby(['site_id','day','month'])['precip_depth_1_hr'].mean()\n",
    "    # Step 2\n",
    "    precip_depth_filler = pd.DataFrame(precip_depth_filler.fillna(method='ffill'),columns=['precip_depth_1_hr'])\n",
    "\n",
    "    weather_df.update(precip_depth_filler,overwrite=False)\n",
    "\n",
    "    weather_df = weather_df.reset_index()\n",
    "    weather_df = weather_df.drop(['datetime','day','week','month'],axis=1)\n",
    "        \n",
    "    return weather_df"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Outliers**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "# outlierrs\n",
    "def detect_outlier(data_1):\n",
    "    \n",
    "    # threshold=8\n",
    "    return data_1[np.abs((data_1- np.mean(data_1))/np.std(data_1)) > 8].values\n",
    "\n",
    "def remove_outliers(df):\n",
    "    \n",
    "    outliers_tuples = []\n",
    "    for building in df.building_id.unique():\n",
    "        print(building)\n",
    "        for meter in df[df['building_id'] == building].meter.unique():\n",
    "            outliers = detect_outlier(df[(df['building_id'] == building) & (df['meter'] == meter)].meter_reading)\n",
    "            if len(outliers) == 1:\n",
    "                outliers_tuples.append([building, meter, outliers[0]])\n",
    "                \n",
    "    return outliers_tuples\n",
    "            \n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "candidatos_out = [\n",
    "    [60,1], [101,0],[103,1],[112,3], [138,3], [144,3], [745,1], [755,1], [762,2], [765,1], [769,1],[769,2], [770,1], [771,1],\n",
    "    [776,1], [777,1],[787,1], [1000,0], [1000,3],[1003,1], [1017,1], [1017,3], [1018,1], [1276, 1], [1426, 2]\n",
    "]\n",
    "#candidatos_out = []"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "_cell_guid": "d4639a4f-461a-47a3-b7e0-8246fbab4140",
    "_uuid": "a5f4e519-ed08-40cf-8dfa-50b98a820628"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[112.0, 3.0]\n",
      "0\n",
      "[138.0, 3.0]\n",
      "1\n",
      "[144.0, 3.0]\n",
      "2\n",
      "[745.0, 1.0]\n",
      "3\n",
      "[755.0, 1.0]\n",
      "4\n",
      "[762.0, 2.0]\n",
      "5\n",
      "[765.0, 1.0]\n",
      "6\n",
      "[769.0, 2.0]\n",
      "7\n",
      "[769.0, 1.0]\n",
      "8\n",
      "[770.0, 1.0]\n",
      "9\n",
      "[771.0, 1.0]\n",
      "10\n",
      "[776.0, 1.0]\n",
      "11\n",
      "[777.0, 1.0]\n",
      "12\n",
      "[787.0, 1.0]\n",
      "13\n",
      "[1000.0, 0.0]\n",
      "14\n",
      "[1000.0, 3.0]\n",
      "15\n",
      "[1003.0, 1.0]\n",
      "16\n",
      "[1017.0, 1.0]\n",
      "17\n",
      "[1017.0, 3.0]\n",
      "18\n",
      "[1018.0, 1.0]\n",
      "19\n",
      "[1276.0, 1.0]\n",
      "20\n",
      "[1426.0, 2.0]\n",
      "21\n",
      "[60.0, 1.0]\n",
      "22\n",
      "[103.0, 1.0]\n",
      "23\n",
      "[101.0, 0.0]\n",
      "24\n",
      "CPU times: user 37.7 s, sys: 20.4 s, total: 58.1 s\n",
      "Wall time: 51.5 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "if DO_TRAIN == True:\n",
    "    df_train = pd.read_csv(path_train)#.sample(10000)\n",
    "    df_train['timestamp'] = pd.to_datetime(df_train['timestamp'])\n",
    "    df_train = df_train.query('not (building_id <= 104 & meter == 0 & timestamp <= \"2016-05-20\")') #++\n",
    "    df_train = df_train.query('not (building_id == 1099 & meter == 2 & meter_reading > 3e4)') #outlier\n",
    "\n",
    "    # csv outliers\n",
    "    outliers = pd.read_csv('/kaggle/input/outliers-energydemand/outliers.csv')\n",
    "    count = 0\n",
    "    for index, fila in outliers.iterrows():\n",
    "        if [fila.building_id, fila.meter] in candidatos_out:\n",
    "            print([fila.building_id, fila.meter])\n",
    "            print(count)\n",
    "            count += 1\n",
    "            df_train = df_train.query('not (building_id == ' + str(int(fila.building_id)) +\\\n",
    "                                      ' & meter == ' + str(int(fila.meter)) +\\\n",
    "                                      ' & meter_reading >= ' + str(fila.meter_reading) + ')')\n",
    "    del outliers\n",
    "    gc.collect()\n",
    "\n",
    "    #get_summer(df_train)\n",
    "    # fin de semana\n",
    "    is_weekend(df_train) #++\n",
    "\n",
    "# BUILDING META\n",
    "building = pd.read_csv(path_building)\n",
    "le = LabelEncoder()\n",
    "building.primary_use = le.fit_transform(building.primary_use)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.6/site-packages/ipykernel_launcher.py:16: FutureWarning: Sorting because non-concatenation axis is not aligned. A future version\n",
      "of pandas will change to not sort by default.\n",
      "\n",
      "To accept the future behavior, pass 'sort=False'.\n",
      "\n",
      "To retain the current behavior and silence the warning, pass 'sort=True'.\n",
      "\n",
      "  app.launch_new_instance()\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 19.7 s, sys: 236 ms, total: 19.9 s\n",
      "Wall time: 19.3 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "if DO_TRAIN == True:\n",
    "    weather_train = pd.read_csv(path_weather_train)\n",
    "    #rellenar nans \n",
    "    weather_train = fill_weather_dataset(weather_train)\n",
    "    # a datetime\n",
    "    weather_train['timestamp'] = pd.to_datetime(weather_train['timestamp'])\n",
    "    # ajustar horas de temperaturas\n",
    "    weather_train = timestamp_align(weather_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Humedad\n",
      "CPU times: user 8 ms, sys: 0 ns, total: 8 ms\n",
      "Wall time: 4.13 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "if DO_TRAIN == True:\n",
    "    print('Humedad')\n",
    "    humidity(weather_train)\n",
    "    #print('Feels like')\n",
    "    #feels_like_df(weather_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 9.78 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "#if DO_TRAIN == True:\n",
    "    #weather_train = create_lag_features(weather_train, 24)\n",
    "    #weather_train = create_lag_features(weather_train, 48)\n",
    "    # TODO probar +48"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Calculando semana del mes...\n",
      "(19866510, 5)\n",
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 138 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "if DO_TRAIN == True:\n",
    "    print('Calculando semana del mes...')\n",
    "    #df_train['weekofmonth'] = df_train['timestamp'].apply(get_week_of_month) #++\n",
    "\n",
    "    # para los meses negativos \n",
    "    # 53 porque es bisiesto\n",
    "    #df_train.loc[df_train.weekofmonth < 0, 'weekofmonth'] = 53 + df_train.loc[df_train.weekofmonth < 0, 'weekofmonth'] #++\n",
    "    print(df_train.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "_cell_guid": "73922343-28da-4e46-b6df-87ae9eb6bd9d",
    "_uuid": "199f5e6f-65c9-4839-b26f-7d173d5b1718"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Memory usage of dataframe is 909.42 MB\n",
      "Memory usage after optimization is: 454.71 MB\n",
      "Decreased by 50.0%\n",
      "Memory usage of dataframe is 0.07 MB\n",
      "Memory usage after optimization is: 0.02 MB\n",
      "Decreased by 66.5%\n",
      "Memory usage of dataframe is 11.80 MB\n",
      "Memory usage after optimization is: 5.63 MB\n",
      "Decreased by 52.3%\n"
     ]
    }
   ],
   "source": [
    "if DO_TRAIN == True:\n",
    "    df_train = reduce_mem_usage(df_train)\n",
    "    building = reduce_mem_usage(building)\n",
    "    weather_train = reduce_mem_usage(weather_train)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "b2343bc5-14b5-4ee5-8f55-3c2a375cda74",
    "_uuid": "dba843bd-10a1-4c8b-99ba-cc98c741b855"
   },
   "source": [
    "## Preparing data\n",
    "There are two files with features that need to be merged with the data. One is building metadata that has information on the buildings and the other is weather data that has information on the weather.   \n",
    "\n",
    "Note that the only features created are hour, weekday and is_holiday!"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {
    "_cell_guid": "5cad6753-3c2a-41ef-bfe8-ae96eb684d28",
    "_uuid": "fa0ee6c6-c5d4-407c-8165-b3cd7b5a8d53"
   },
   "outputs": [],
   "source": [
    "def prepare_data(X, building_data, weather_data, test=False):\n",
    "    \"\"\"\n",
    "    Preparing final dataset with all features.\n",
    "    \"\"\"\n",
    "    \n",
    "    X = X.merge(building_data, on=\"building_id\", how=\"left\")\n",
    "    X = X.merge(weather_data, on=[\"site_id\", \"timestamp\"], how=\"left\")\n",
    "    \n",
    "    X.sort_values(\"timestamp\")\n",
    "    X.reset_index(drop=True)\n",
    "    \n",
    "    gc.collect()\n",
    "    \n",
    "    #X.timestamp = pd.to_datetime(X.timestamp, format=\"%Y-%m-%d %H:%M:%S\")\n",
    "    X.square_feet = np.log1p(X.square_feet)\n",
    "    \n",
    "    X['year'] = np.uint8(X['timestamp'].dt.year) #++\n",
    "    #X['week'] = np.uint8(X['timestamp'].dt.week) #++\n",
    "    #X['day'] = X['timestamp'].dt.day #++\n",
    "    X[\"hour\"] = np.uint8(X.timestamp.dt.hour) #--\n",
    "    X[\"weekday\"] = np.uint8(X.timestamp.dt.weekday) #--\n",
    "    #X['year_built'] = np.uint8(X['year_built']-1900) #++\n",
    "    #X['hourofyear']  = np.uint16(np.uint16(X.timestamp.dt.dayofyear) * (np.uint8(X.timestamp.dt.hour) + 1))\n",
    "    \n",
    "    #X[\"is_holiday\"] = (X.timestamp.dt.date.astype(\"str\").isin(holidays)).astype(int) #--\n",
    "    '''\n",
    "    X[\"is_holiday\"] = np.zeros(X.shape[0], dtype=int)\n",
    "    for site in range(0, 16):\n",
    "        print('Llenando vacaciones de sitio ', site)\n",
    "        X.loc[X.site_id == site, \"is_holiday\"] = (X.timestamp.dt.date.isin(get_holidays(site))).astype(int)\n",
    "        gc.collect()\n",
    "\n",
    "    # https://www.kaggle.com/ragnar123/another-1-08-lb-no-leak\n",
    "    X['group'] = X['timestamp'].dt.month\n",
    "    X['group'].replace((1, 2, 3, 4), 1, inplace = True)\n",
    "    X['group'].replace((5, 6, 7, 8), 2, inplace = True)\n",
    "    X['group'].replace((9, 10, 11, 12), 3, inplace = True)\n",
    "    '''\n",
    "    # nuevas\n",
    "    #beaufort_scale(X)\n",
    "    #wind_direction_to_categorical(X)\n",
    "    \n",
    "    drop_features = [\"sea_level_pressure\", 'floor_count', 'offset',\n",
    "                    'precip_depth_1_hr','year_built', \"wind_direction\",\n",
    "                     'dew_temperature']  #'timestamp', 'feels', 'beaufort_scale',\n",
    "    \n",
    "    if 'floor_count' not in drop_features:\n",
    "        # para tratar el enorme numero de floor_count a Nan\n",
    "        X['floor_count'] = X['floor_count'].fillna(1)\n",
    "\n",
    "    X.drop(drop_features, axis=1, inplace=True)\n",
    "\n",
    "    if test:\n",
    "        row_ids = X.row_id\n",
    "        X.drop(\"row_id\", axis=1, inplace=True)\n",
    "        X.drop(['timestamp'], axis=1, inplace=True) #++\n",
    "        return X, row_ids\n",
    "    else:\n",
    "        y = np.log1p(X.meter_reading)\n",
    "        X.drop(\"meter_reading\", axis=1, inplace=True)\n",
    "        '''\n",
    "        # bad rows\n",
    "        print('shape', X.shape)\n",
    "        bad_rows = find_bad_rows(X, y)\n",
    "        X = X.drop(bad_rows)\n",
    "        print('bad rows', bad_rows.shape)\n",
    "        print('shape', X.shape)\n",
    "        y = y.drop(bad_rows)\n",
    "        '''\n",
    "        X.drop(['timestamp'], axis=1, inplace=True)\n",
    "\n",
    "        return X, y\n",
    "   \n",
    "#----------------------------------------------------------------------------------------------------------------------------------------\n",
    "def get_leak_validation():\n",
    "    \n",
    "    leak_df = pd.read_feather('../input/ashrae-leak-data-station/leak.feather')#.sample(10000)\n",
    "    leak_df.fillna(0, inplace=True)\n",
    "    leak_df = leak_df[(leak_df.timestamp.dt.year > 2016) & (leak_df.timestamp.dt.year < 2019)]\n",
    "    leak_df.loc[leak_df.meter_reading < 0, 'meter_reading'] = 0 # remove large negative values\n",
    "    leak_df = leak_df[leak_df.building_id!=245]\n",
    "    print(leak_df.shape)\n",
    "    \n",
    "    # lo que se necesita para estar ok\n",
    "    #print('Calculando semana del mes...')\n",
    "    leak_df['timestamp'] = pd.to_datetime(leak_df['timestamp'])\n",
    "    #leak_df['weekofmonth'] = leak_df['timestamp'].apply(get_week_of_month)\n",
    "    # para los meses negativos\n",
    "    #leak_df.loc[leak_df.weekofmonth < 0, 'weekofmonth'] = 52 + leak_df.loc[leak_df.weekofmonth < 0, 'weekofmonth'] \n",
    "    \n",
    "    is_weekend(leak_df)\n",
    "\n",
    "    print('temperatura leak')\n",
    "    weather_test = pd.read_csv(path_weather_test)\n",
    "    weather_test = fill_weather_dataset(weather_test)\n",
    "    weather_test['timestamp'] = pd.to_datetime(weather_test['timestamp'])\n",
    "    weather_test = timestamp_align(weather_test)\n",
    "    humidity(weather_test)\n",
    "    #feels_like_df(weather_test)\n",
    "\n",
    "    leak_df = reduce_mem_usage(leak_df)\n",
    "    weather_test = reduce_mem_usage(weather_test)\n",
    "    gc.collect()\n",
    "    \n",
    "    print('Juntando leak')\n",
    "    X_train_leak, y_train_leak = prepare_data(leak_df, building, weather_test)\n",
    "    del leak_df, weather_test\n",
    "    gc.collect()\n",
    "    \n",
    "    return X_train_leak, y_train_leak"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {
    "_cell_guid": "731b6804-a594-481b-9921-fe8b4e8c2395",
    "_kg_hide-output": false,
    "_uuid": "0e10bd0c-1ef6-427a-bbb0-c467f632c7d2"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 17.3 s, sys: 10.5 s, total: 27.8 s\n",
      "Wall time: 25 s\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "if DO_TRAIN == True:\n",
    "    X_train, y_train = prepare_data(df_train, building, weather_train)\n",
    "    #X_train = X_train.drop(columns=[\"index\"])\n",
    "    del df_train, weather_train\n",
    "    gc.collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "33c30319-6bbe-469e-a8b3-1461011da7f4",
    "_uuid": "e0889554-d431-4044-864d-d43c6108591a"
   },
   "source": [
    "AQUI LOS ALGORITMOS \n",
    "---"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {
    "_cell_guid": "a4e648f6-93aa-4902-a74b-419e587e61c4",
    "_uuid": "b8f445e5-21ad-4a5b-a84d-72cb7715447b"
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'\\nparams = {\"objective\": \"regression\",\\n                  \"num_leaves\": 600,\\n                  \"learning_rate\": 0.005,\\n                  \"bagging_freq\": 5,\\n                  \"bagging_fraction\": 0.51,\\n                  \"feature_fraction\": 0.81,\\n                  \"metric\": \"rmse\"\\n                  }\\n'"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "metric = 'rmse' \n",
    "tecnica_split = 'normal'\n",
    "params = {\n",
    "    \"objective\": \"regression\",\n",
    "    \"boosting\": \"gbdt\", #gbdt\n",
    "    \"num_leaves\": 600, #600 numero hojas por arbol\n",
    "    \"learning_rate\": 0.005, #0.005\n",
    "    \"feature_fraction\": 0.85,\n",
    "    #'bagging_fraction': 0.85,\n",
    "    #'max_bin': 512,\n",
    "    #'reg_alpha': 5,#lambda1\n",
    "    \"reg_lambda\": 2, #lambda2\n",
    "    #'max_depth': 15, #++\n",
    "    #'min_data_in_leaf': 5, #++\n",
    "    #'min_child_weight': 1,\n",
    "    #'min_split_gain': 0.001,\n",
    "    \"metric\": metric\n",
    "}\n",
    "'''\n",
    "params = {\"objective\": \"regression\",\n",
    "                  \"num_leaves\": 600,\n",
    "                  \"learning_rate\": 0.005,\n",
    "                  \"bagging_freq\": 5,\n",
    "                  \"bagging_fraction\": 0.51,\n",
    "                  \"feature_fraction\": 0.81,\n",
    "                  \"metric\": \"rmse\"\n",
    "                  }\n",
    "'''"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {
    "_cell_guid": "f05544d5-bbc6-40a6-b578-195a6a182aa4",
    "_uuid": "c995fab0-c470-4664-b126-271bfc86ebac"
   },
   "outputs": [],
   "source": [
    "#---------------------------------------------------------------------------------------------\n",
    "def LIGHTGBM_model(metric, params, train_set, valid_sets):\n",
    "    \n",
    "    print('Lightgbm model')\n",
    "    model_lgbm = lgb.train(params,\n",
    "                           train_set=train_set, num_boost_round=6000, valid_sets=valid_sets,\n",
    "                           verbose_eval=200, early_stopping_rounds=500)\n",
    "    print('best_score', model_lgbm.best_score)\n",
    "    log = {'train/mae': model_lgbm.best_score['training'][metric],\n",
    "           'valid/mae': model_lgbm.best_score['valid_1'][metric]}\n",
    "    \n",
    "    return model_lgbm, log\n",
    "\n",
    "#---------------------------------------------------------------------------------------------  "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {
    "_cell_guid": "c58429b8-ef5a-4e0b-88d2-179883f60d4e",
    "_uuid": "c9a25548-ff86-4a82-a2e6-3907b6a29583"
   },
   "outputs": [],
   "source": [
    "\n",
    "folds =3\n",
    "seed = 666\n",
    "shuffle = False\n",
    "group_folds_label = ['building_id'] # ,building_id, weekofmonth\n",
    "DO_FOLDS = True\n",
    "kf = KFold(n_splits=folds, shuffle=shuffle, random_state=seed)\n",
    "skf = StratifiedKFold(n_splits=folds, shuffle=shuffle, random_state=seed)\n",
    "LIBRERIA = 'lgb' #cat, lgb, xgb"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {
    "_cell_guid": "60dceda4-1f42-4800-ae5a-aec25941e3fc",
    "_uuid": "61fff533-738b-4cc1-8cea-6edd2eb461e3"
   },
   "outputs": [],
   "source": [
    "#------------------------------------------------------------------------------------------------------------\n",
    "def validate_with_leak(X_train, y_train):\n",
    "    \n",
    "    leak_valid_x, leak_valid_y = get_leak_validation()\n",
    "    y_pred_valid = np.zeros(leak_valid_x.shape[0])\n",
    "    \n",
    "    print('Preparando datasets')\n",
    "    d_train = lgb.Dataset(X_train, label=y_train, categorical_feature=CATEGORY_COLS)\n",
    "    d_valid = lgb.Dataset(leak_valid_x, label=leak_valid_y, categorical_feature=CATEGORY_COLS)\n",
    "    watchlist = [d_train, d_valid]\n",
    "\n",
    "    model, log = LIGHTGBM_model(metric, params, d_train, watchlist)\n",
    "\n",
    "    print('Prediccion...')\n",
    "    y_pred_valid = model.predict(leak_valid_x, num_iteration=model.best_iteration)\n",
    "\n",
    "    plot_fi(model)\n",
    "    # save model\n",
    "    joblib.dump(model, LIBRERIA + '_' + 'leak' + '.pkl')\n",
    "    \n",
    "    print('RMSE:', np.sqrt(mean_squared_error( np.clip(leak_valid_y, 0, np.max(leak_valid_y)),\n",
    "                                              np.clip(y_pred_valid, 0,np.max(y_pred_valid) ))))\n",
    "    return [model]\n",
    "#-------------------------------------------------------------------------------------------------------------\n",
    "def leak_kfolds(X_train, y_train, kf):\n",
    "    \n",
    "    n_folds = 0\n",
    "    leak_valid_x, leak_valid_y = get_leak_validation()\n",
    "    y_valid_pred_total = np.zeros(leak_valid_x.shape[0])\n",
    "    \n",
    "    modelos = []\n",
    "    for train_idx, valid_idx in kf.split(X_train, y_train, groups=X_train[group_folds_label]): #kf.split(...)\n",
    "        print('Folds --> ', n_folds)\n",
    "        n_folds += 1\n",
    "        X_train_fold, y_train_fold = X_train.iloc[train_idx,:], y_train[train_idx]\n",
    "        print('train', len(train_idx), 'valid', len(leak_valid_y))\n",
    "        #--------------------------------------------------------------------------------------------\n",
    "        d_train = lgb.Dataset(X_train_fold, label=y_train_fold, categorical_feature=CATEGORY_COLS)\n",
    "        d_valid = lgb.Dataset(leak_valid_x, label=leak_valid_y, categorical_feature=CATEGORY_COLS)\n",
    "        watchlist = [d_train, d_valid]\n",
    "\n",
    "        model, log = LIGHTGBM_model(metric, params, d_train, watchlist)\n",
    "\n",
    "        y_pred_valid = np.zeros(leak_valid_x.shape[0])\n",
    "        y_pred_valid = model.predict(leak_valid_x, num_iteration=model.best_iteration)\n",
    "        print('RMSE SUB MODELO: ', np.sqrt(mean_squared_error( np.clip(leak_valid_y, 0, np.max(leak_valid_y)),\n",
    "                                              np.clip(y_pred_valid, 0,np.max(y_pred_valid) ))))\n",
    "\n",
    "        plot_fi(model)\n",
    "        #--------------------------------------------------------------------------------------------\n",
    "        # save model\n",
    "        joblib.dump(model, LIBRERIA + '_' + 'leak_' +  str(n_folds) + '.pkl')\n",
    "        \n",
    "        y_valid_pred_total += y_pred_valid\n",
    "        \n",
    "        modelos.append(model)\n",
    "        gc.collect()\n",
    "     \n",
    "    #++++++++++++++\n",
    "    y_valid_pred_total = y_valid_pred_total/len(modelos)\n",
    "    print('RMSE:', np.sqrt(mean_squared_error( np.clip(leak_valid_y, 0, np.max(leak_valid_y)),\n",
    "                                              np.clip(y_valid_pred_total, 0,np.max(y_valid_pred_total) ))))\n",
    "    \n",
    "    return modelos\n",
    "    \n",
    "#-------------------------------------------------------------------------------------------------------------\n",
    "def kfolds(X_train, y_train, kf, libreria='cat'):\n",
    "\n",
    "    n_folds = 0\n",
    "    y_valid_pred_total = np.zeros(X_train.shape[0])\n",
    "    modelos = []\n",
    "    for train_idx, valid_idx in kf.split(X_train, y_train, groups=X_train[group_folds_label]): #kf.split(...)\n",
    "        print('Folds --> ', n_folds)\n",
    "        n_folds += 1\n",
    "        X_train_fold, y_train_fold = X_train.iloc[train_idx,:], y_train[train_idx]\n",
    "        X_valid_fold, y_valid_fold = X_train.iloc[valid_idx,:], y_train[valid_idx]\n",
    "        print('train', len(train_idx), 'valid', len(valid_idx))\n",
    "        #--------------------------------------------------------------------------------------------\n",
    "        d_train = lgb.Dataset(X_train_fold, label=y_train_fold, categorical_feature=CATEGORY_COLS)\n",
    "        d_valid = lgb.Dataset(X_valid_fold, label=y_valid_fold, categorical_feature=CATEGORY_COLS)\n",
    "        watchlist = [d_train, d_valid]\n",
    "\n",
    "        model, log = LIGHTGBM_model(metric, params, d_train, watchlist)\n",
    "\n",
    "        y_pred_valid = model.predict(X_valid_fold, num_iteration=model.best_iteration)\n",
    "\n",
    "        plot_fi(model)\n",
    "        #--------------------------------------------------------------------------------------------\n",
    "        # save model\n",
    "        joblib.dump(model, LIBRERIA + '_' + str(n_folds) + '.pkl')\n",
    "        #model.save_model(LIBRERIA + '_' + str(n_folds) + '.txt')\n",
    "        \n",
    "        y_valid_pred_total[valid_idx] = y_pred_valid\n",
    "        \n",
    "        modelos.append(model)\n",
    "        gc.collect()\n",
    "     \n",
    "    #++++++++++++++\n",
    "    print('RMSE:', np.sqrt(mean_squared_error( np.clip(y_train, 0, np.max(y_train)),\n",
    "                                              np.clip(y_valid_pred_total, 0,np.max(y_valid_pred_total) ))))\n",
    "    return modelos\n",
    "#-------------------------------------------------------------------------------------------------------------\n",
    "#TODO\n",
    "def Stratify_kfolds(X_train, y_train, kf):\n",
    "\n",
    "    n_folds = 0\n",
    "    y_valid_pred_total = np.zeros(X_train.shape[0])\n",
    "    modelos = []\n",
    "    for train_idx, valid_idx in tqdm(kf.split(X_train, X_train['group']), total = 3):\n",
    "        print('Folds --> ', n_folds)\n",
    "        n_folds += 1\n",
    "        X_train_fold, y_train_fold = X_train.iloc[train_idx,:], y_train[train_idx]\n",
    "        X_valid_fold, y_valid_fold = X_train.iloc[valid_idx,:], y_train[valid_idx]\n",
    "        print('train', len(train_idx), 'valid', len(valid_idx))\n",
    "        #--------------------------------------------------------------------------------------------\n",
    "        d_train = lgb.Dataset(X_train_fold, label=y_train_fold, categorical_feature=CATEGORY_COLS)\n",
    "        d_valid = lgb.Dataset(X_valid_fold, label=y_valid_fold, categorical_feature=CATEGORY_COLS)\n",
    "        watchlist = [d_train, d_valid]\n",
    "\n",
    "        model, log = LIGHTGBM_model(metric, params, d_train, watchlist)\n",
    "\n",
    "        y_pred_valid = model.predict(X_valid_fold, num_iteration=model.best_iteration)\n",
    "\n",
    "        plot_fi(model)\n",
    "        #--------------------------------------------------------------------------------------------\n",
    "        # save model\n",
    "        joblib.dump(model, LIBRERIA + '_' + str(n_folds) + '.pkl')\n",
    "        #model.save_model(LIBRERIA + '_' + str(n_folds) + '.txt')\n",
    "        \n",
    "        y_valid_pred_total[valid_idx] = y_pred_valid\n",
    "        \n",
    "        modelos.append(model)\n",
    "        gc.collect()\n",
    "     \n",
    "    #++++++++++++++\n",
    "    print('RMSE:', np.sqrt(mean_squared_error( np.clip(y_train, 0, np.max(y_train)),\n",
    "                                              np.clip(y_valid_pred_total, 0,np.max(y_valid_pred_total) ))))\n",
    "    return modelos\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "e5bfc513-0e06-46fa-abd3-24958679a1b3",
    "_uuid": "99a3ff01-62df-403c-b2fd-50b7c9b1e4a7"
   },
   "source": [
    "Aqui entrenamos\n",
    "----"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "if DO_TRAIN == True:\n",
    "    X_train.head()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {
    "_cell_guid": "3a6ed1da-57a1-495b-a990-8396373cba8b",
    "_uuid": "85c105e1-9dc5-4ed6-9664-fd502911e11d"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Folds -->  0\n",
      "train 13244340 valid 6622170\n",
      "Lightgbm model\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/opt/conda/lib/python3.6/site-packages/lightgbm/basic.py:1243: UserWarning: Using categorical_feature in Dataset.\n",
      "  warnings.warn('Using categorical_feature in Dataset.')\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Training until validation scores don't improve for 500 rounds\n",
      "[200]\ttraining's rmse: 1.25224\tvalid_1's rmse: 1.36931\n",
      "[400]\ttraining's rmse: 0.999468\tvalid_1's rmse: 1.17745\n",
      "[600]\ttraining's rmse: 0.910382\tvalid_1's rmse: 1.13241\n",
      "[800]\ttraining's rmse: 0.862201\tvalid_1's rmse: 1.11935\n",
      "[1000]\ttraining's rmse: 0.830308\tvalid_1's rmse: 1.11546\n",
      "[1200]\ttraining's rmse: 0.807685\tvalid_1's rmse: 1.11542\n",
      "[1400]\ttraining's rmse: 0.791478\tvalid_1's rmse: 1.11583\n",
      "Early stopping, best iteration is:\n",
      "[1081]\ttraining's rmse: 0.820159\tvalid_1's rmse: 1.1151\n",
      "best_score defaultdict(<class 'collections.OrderedDict'>, {'training': OrderedDict([('rmse', 0.820159388638609)]), 'valid_1': OrderedDict([('rmse', 1.1151045483197137)])})\n",
      "Folds -->  1\n",
      "train 13244340 valid 6622170\n",
      "Lightgbm model\n",
      "Training until validation scores don't improve for 500 rounds\n",
      "[200]\ttraining's rmse: 1.24787\tvalid_1's rmse: 1.34059\n",
      "[400]\ttraining's rmse: 0.994204\tvalid_1's rmse: 1.13162\n",
      "[600]\ttraining's rmse: 0.914843\tvalid_1's rmse: 1.06781\n",
      "[800]\ttraining's rmse: 0.87776\tvalid_1's rmse: 1.04249\n",
      "[1000]\ttraining's rmse: 0.85338\tvalid_1's rmse: 1.0315\n",
      "[1200]\ttraining's rmse: 0.833863\tvalid_1's rmse: 1.02684\n",
      "[1400]\ttraining's rmse: 0.819975\tvalid_1's rmse: 1.02511\n",
      "[1600]\ttraining's rmse: 0.809799\tvalid_1's rmse: 1.02393\n",
      "[1800]\ttraining's rmse: 0.801892\tvalid_1's rmse: 1.02296\n",
      "[2000]\ttraining's rmse: 0.79523\tvalid_1's rmse: 1.0224\n",
      "[2200]\ttraining's rmse: 0.789778\tvalid_1's rmse: 1.02225\n",
      "[2400]\ttraining's rmse: 0.783902\tvalid_1's rmse: 1.02205\n",
      "[2600]\ttraining's rmse: 0.778988\tvalid_1's rmse: 1.02186\n",
      "[2800]\ttraining's rmse: 0.774586\tvalid_1's rmse: 1.02176\n",
      "[3000]\ttraining's rmse: 0.770349\tvalid_1's rmse: 1.02188\n",
      "[3200]\ttraining's rmse: 0.766774\tvalid_1's rmse: 1.02174\n",
      "[3400]\ttraining's rmse: 0.76394\tvalid_1's rmse: 1.02179\n",
      "[3600]\ttraining's rmse: 0.761555\tvalid_1's rmse: 1.02174\n",
      "[3800]\ttraining's rmse: 0.758943\tvalid_1's rmse: 1.02165\n",
      "[4000]\ttraining's rmse: 0.756676\tvalid_1's rmse: 1.0216\n",
      "[4200]\ttraining's rmse: 0.754248\tvalid_1's rmse: 1.02164\n",
      "[4400]\ttraining's rmse: 0.752006\tvalid_1's rmse: 1.0216\n",
      "[4600]\ttraining's rmse: 0.749647\tvalid_1's rmse: 1.02157\n",
      "[4800]\ttraining's rmse: 0.747237\tvalid_1's rmse: 1.0215\n",
      "[5000]\ttraining's rmse: 0.744798\tvalid_1's rmse: 1.02138\n",
      "[5200]\ttraining's rmse: 0.742617\tvalid_1's rmse: 1.02143\n",
      "[5400]\ttraining's rmse: 0.740413\tvalid_1's rmse: 1.02134\n",
      "[5600]\ttraining's rmse: 0.738306\tvalid_1's rmse: 1.0214\n",
      "[5800]\ttraining's rmse: 0.735859\tvalid_1's rmse: 1.02149\n",
      "Early stopping, best iteration is:\n",
      "[5387]\ttraining's rmse: 0.740555\tvalid_1's rmse: 1.02133\n",
      "best_score defaultdict(<class 'collections.OrderedDict'>, {'training': OrderedDict([('rmse', 0.7405554943253191)]), 'valid_1': OrderedDict([('rmse', 1.0213273905502769)])})\n",
      "Folds -->  2\n",
      "train 13244340 valid 6622170\n",
      "Lightgbm model\n",
      "Training until validation scores don't improve for 500 rounds\n",
      "[200]\ttraining's rmse: 1.23542\tvalid_1's rmse: 1.35737\n",
      "[400]\ttraining's rmse: 0.96694\tvalid_1's rmse: 1.18001\n",
      "[600]\ttraining's rmse: 0.873503\tvalid_1's rmse: 1.14293\n",
      "[800]\ttraining's rmse: 0.821707\tvalid_1's rmse: 1.13448\n",
      "[1000]\ttraining's rmse: 0.78944\tvalid_1's rmse: 1.13331\n",
      "[1200]\ttraining's rmse: 0.765399\tvalid_1's rmse: 1.1348\n",
      "[1400]\ttraining's rmse: 0.74844\tvalid_1's rmse: 1.13677\n",
      "Early stopping, best iteration is:\n",
      "[981]\ttraining's rmse: 0.79234\tvalid_1's rmse: 1.13322\n",
      "best_score defaultdict(<class 'collections.OrderedDict'>, {'training': OrderedDict([('rmse', 0.7923398749189288)]), 'valid_1': OrderedDict([('rmse', 1.1332238452914014)])})\n",
      "RMSE: 1.0908135764337845\n",
      "CPU times: user 22h 49min 33s, sys: 5min 29s, total: 22h 55min 2s\n",
      "Wall time: 5h 53min 57s\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAHVCAYAAADb6QDfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucXVV9///XmwQTbga5aCOoQRpAJRDIgAKC4AWteC1BVGoB9ZvijfbXL/WbWqtUa4vf+hMQ8RIVQUVFvCJYAZVwv2QCIYFyh7GCFhUlolyE8Pn+cfbgMMxMJskkZ2fm9Xw85jH7rLP2Wmsvz5g3a+1zTqoKSZKkNtmg2wOQJEkazIAiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJaZ3K3BzCRbbXVVjVjxoxuD0OSpHVi8eLFv66qrUdT14DSRTNmzKC3t7fbw5AkaZ1I8tPR1nWLR5IktY4BRZIktY4BRZIktY4BRZIktY4BRZIktY4BRZIktY5vM+6iZXctZ8b8c7o9DEmSnqDvuIO62r8rKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXVaE1CSzEhy3SrUf02S+c3xsUmOGanNJD1JPjF2I4YkH0ry0iHK909y9lj2JUnSRLLevs24qs4CzlqF+r3AmH51cFV9YCzbkyRJHa1ZQWlMTnJakqVJvplk4yR9SbaCx1ZBFjbHRyT55OAGksxJcm2Sy4F3DSh/bFWjWXE5JcnCJLcnOXpAvX9OcmOS85N8baiVmQF1T00ytzl+RXPeJcBfjtF8SJI0IbUtoOwILKiqXYDfAe9cjTa+CBxdVXutpN5OwMuBPYEPJtkwSQ9wMLAbnZDRM5oOk0wFPge8GtgX+LMR6s5L0pukd8X9y0fTvCRJE07bAsrPqurS5vgrwAtX5eQk04DNq+rCpujLI1Q/p6oeqqpfA78Entb0972qeqCq7gO+P8qudwLuqKpbqqqasQ+pqhZUVU9V9UzaeNoom5ckaWJpW0CpIR4/wp/GOXUl52eINobz0IDjFXTux8kozx3KaPuVJEkr0baA8swk/VszbwIuAfqAOU3ZwSOdXFX3AsuT9K+8HLaK/V8CvDrJ1CSbAqP9IoIbge2SbN88ftMq9itJkgZoW0C5ATg8yVJgC+DTwL8AJya5mM5Kx8ocCZzc3CT7wKp0XlWL6Lwz6Frg23Te9bPSG0Wq6kFgHnBOc5PsT1elX0mS9Hjp3DKhfkk2rarfJ9kYuAiYV1VXr42+pkyfWdMPP2FtNC1J0hpZG99mnGRxVY3qDSjr7eegrEULkjyXzv0up62tcCJJkoZnQBmkqt48uCzJycA+g4pPrKovrptRSZI0sbjF00U9PT3V2zumH24rSVJrrcoWT9tukpUkSTKgSJKk9jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1pnc7QFMZMvuWs6M+ed0exjSkPqOO6jbQ5A0gbmCIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWme1A0qSHyTZfBXqv291+1qXkhyR5OndHockSRPZageUqnplVd07sCwdw7XZmoCSZNIITx8BrFJASeLbtSVJGkOjCihJvptkcZLrk8xryvqSbJVkRpIbknwKuBp4xhDnHwdslGRJktObsr9KclVT9tn+0JDk90k+2vT3oyR7JlmY5PYkr2nqHJHke0l+mOSmJB8c0NdI7X4oyZXAXkk+kGRRkuuSLGjC1VygBzi9OX+j/uts2uhJsrA5PrY57zzgS0kmJfmPps2lSf5mmLmcl6Q3Se+K+5ePZvolSZpwRruC8taqmkPnH++jk2w56PkdgS9V1W5V9dPBJ1fVfOCBqppdVYcleQ5wKLBPVc0GVgCHNdU3ARY2/d0H/CvwMuD1wIcGNLtnc85s4JAmPKys3euq6vlVdQnwyarao6p2BjYCXlVV3wR6gcOasT6wknmZA7y2qt4MvA1YXlV7AHsA/yvJdkPMxYKq6qmqnkkbT1tJ85IkTUyj3Zo4Osnrm+NnADMHPf/TqrpiFfp9CZ1/3BclgU5A+GXz3B+BHzbHy4CHqurhJMuAGQPaOL+q7gFI8m3ghcAjI7S7AvjWgPMPSPJeYGNgC+B64PurcA0AZw0IMQcCuzSrMADT6MzTHavYpiRJE95KA0qS/YGXAntV1f3NFsfUQdX+sIr9Bjitqv5xiOcerqpqjh8FHgKoqkcH3etRg86rlbT7YFWtAEgyFfgU0FNVP0tyLE+8pn6P8KeVppGuO8B7qurcYdqRJEmjNJotnmnAb5twshPwgtXs6+EkGzbHPwbmJnkqQJItkjxrFdt7WXPeRsDrgEtXod3+oPHrJJsCcwc8dx+w2YDHfXRWZQAOHmE85wLv6L/GJDsk2WQVr0mSJDG6LZ4fAkclWQrcBKzKVs5AC4ClSa5u7kN5P3Be866fh4F3AU+4f2UElwBfBv4c+GpV9QKMpt2qujfJ5+hsIfUBiwY8fSrwmSQPAHsB/wJ8oXmb9JUjjOfzdLagrk5nf+lXdIKTJElaRfnTbsr6I8kRdLZn3t3tsayJKdNn1vTDT+j2MKQh+W3GksZaksVV1TOaun6SrCRJap0xX0FpPmdkyqDit1TVsjHtaBzo6emp3t7ebg9DkqR1YlVWUMb8E1Cr6vlj3aYkSZpY3OKRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtM7nbA5jIlt21nBnzz+n2MDQKfccd1O0hSNKE4gqKJElqHQOKJElqHQOKJElqHQOKJElqHQPKKCWZneSV3R6HJEkTgQFl9GYDqxRQkvguKUmSVsOECihJZiS5Mcnnk1yX5PQkL01yaZJbkuyZZJMkpyRZlOSaJK9N8iTgQ8ChSZYkOXSoek0fRyQ5M8n3gfO6esGSJK2nJuJ/4f85cAgwD1gEvBl4IfAa4H3AfwE/qaq3JtkcuAr4EfABoKeq3g2Q5N8G10vyo6aPvYBdquo36/C6JEkaNyZiQLmjqpYBJLke+HFVVZJlwAxgW+A1SY5p6k8FnjlEOweOUO/84cJJknl0whGTnrz1GFyOJEnjz0QMKA8NOH50wONH6czHCuDgqrpp4ElJnj+onYxQ7w/DdV5VC4AFAFOmz6zVuQBJksa7CXUPyiidC7wnSQCS7NaU3wdsNop6kiRpDRlQnujDwIbA0iTXNY8BLgCe23+T7Aj1JEnSGkqVuwzdMmX6zJp++AndHoZGwS8LlKQ1l2RxVfWMpq4rKJIkqXUMKJIkqXUMKJIkqXUm4tuMW2PWNtPo9d4GSZKewBUUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOpO7PYCJbNldy5kx/5xuD2Pc6TvuoG4PQZK0hlxBkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAaSSZkeS6bo9DkiQZUNaqJL5LSpKk1WBAebxJST6X5Pok5yXZKMnsJFckWZrkO0meApBkYZKe5nirJH3N8RFJzkzyfeC87l2KJEnrLwPK480ETq6q5wH3AgcDXwL+T1XtAiwDPjiKdvYCDq+qFw9+Ism8JL1Jelfcv3wMhy5J0vhhQHm8O6pqSXO8GNge2LyqLmzKTgP2G0U751fVb4Z6oqoWVFVPVfVM2njamo9YkqRxyIDyeA8NOF4BbD5C3Uf40/xNHfTcH8ZyUJIkTTQGlJEtB36bZN/m8VuA/tWUPmBOczx3HY9LkqRxzXeZrNzhwGeSbAzcDhzZlH8M+EaStwA/6dbgJEkaj1JV3R7DhDVl+syafvgJ3R7GuOOXBUpSOyVZXFU9o6nrFo8kSWodA4okSWod70HpolnbTKPX7QhJkp7AFRRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6k7s9gIls2V3LmTH/nG4Po5X6jjuo20OQJHWRKyiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DChDSLJTkiVJrkmy/Wqc/3dJNl4bY5MkaSKYEAElyaRVPOV1wPeqarequm01uvw7wIAiSdJq6lpASbJJknOSXJvkuiSHJnlFkhuTXJLkE0nObuoem+SYAedel2RGc/zdJIuTXJ9k3oA6v0/yoSRXAnslmZPkwqbuuUmmDzOuV9IJGG9PckFT9ldJrmpWVT7bH3iSHJjk8iRXJzkzyaZJjgaeDlzQf/6g9ucl6U3Su+L+5WM0m5IkjS/dXEF5BfDzqtq1qnYGfgh8Dng1sC/wZ6Ns561VNQfoAY5OsmVTvglwXVU9H7gSOAmY29Q9BfjIUI1V1Q+AzwDHV9UBSZ4DHArsU1WzgRXAYUm2At4PvLSqdgd6gb+vqk8APwcOqKoDhmh/QVX1VFXPpI2njfISJUmaWLr5SbLLgI8l+ShwNnAfcEdV3QKQ5CvAvBHO73d0ktc3x88AZgL30AkS32rKdwR2Bs5PAjAJ+MUox/kSYA6wqDl3I+CXwAuA5wKXNuVPAi4fZZuSJGkEXQsoVXVzkjnAK4F/B84Dapjqj/D41Z6pAEn2B14K7FVV9ydZ2P8c8GBVrWiOA1xfVXutxlADnFZV//i4wuTVwPlV9abVaFOSJI2gm/egPB24v6q+AnwM2BvYbsC7Zgb+w98H7N6ctzuwXVM+DfhtE052orOqMZSbgK2T7NW0sWGS541yqD8G5iZ5anPuFkmeBVwB7JPkz5vyjZPs0JxzH7DZKNuXJEmDdHOLZxbwH0keBR4G3gFsBZyT5NfAJXS2ZaCzVfPXSZYAi4Cbm/IfAkclWUonhFwxVEdV9cckc4FPJJlG57pPAK5f2SCr6r+SvB84L8kGzVjfVVVXJDkC+FqSKU319zdjWwD8Z5JfDHUfiiRJGlmqhttV6a5m++aYqnpVt8eytkyZPrOmH35Ct4fRSn6bsSSNP0kWV1XPaOpOiM9BkSRJ65dubvGMqKoWAgvXZh9JTgb2GVR8YlV9cW32K0mSRtbaLZ6JoKenp3p7e7s9DEmS1gm3eCRJ0nrNgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklpncrcHMJEtu2s5M+af0+1hrFN9xx3U7SFIktYDrqBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWWS8DSpIZSa5bC+1+KMlLhyjfP8nZzfFrksxvjl+X5LljPQ5JkiY638UzQFV9YBR1zgLOah6+Djgb+K+1OS5Jkiaa9XIFpTEpyeeSXJ/kvCQbJVmYpAcgyVZJ+prjI5J8N8n3k9yR5N1J/j7JNUmuSLJFU+/UJHOb41ckuTHJJcBf9nfatPXJJHsDrwH+I8mSJNsnuXpAvZlJFq+76ZAkafxYnwPKTODkqnoecC9w8Erq7wy8GdgT+Ahwf1XtBlwO/PXAikmmAp8DXg3sC/zZ4Maq6jI6Kyn/UFWzq+o2YHmS2U2VI4FTB5+XZF6S3iS9K+5fPtprlSRpQlmfA8odVbWkOV4MzFhJ/Quq6r6q+hWwHPh+U75siHN3atq/paoK+Moox/R54Mgkk4BDga8OrlBVC6qqp6p6Jm08bZTNSpI0sazPAeWhAccr6NxP8wh/uqapI9R/dMDjRxn6XpxajTF9C/gL4FXA4qq6ZzXakCRpwlufA8pQ+oA5zfHcNWjnRmC7JNs3j980TL37gM36H1TVg8C5wKeBL65B/5IkTWjjLaB8DHhHksuArVa3kSZozAPOaW6S/ekwVb8O/ENzs21/mDmdzurLeavbvyRJE106t1horCQ5BphWVf+8srpTps+s6YefsA5G1R5+WaAkTVxJFldVz2jq+jkoYyjJd4DtgRd3eyySJK3PDChjqKpe3+0xSJI0HhhQumjWNtPodctDkqQnGG83yUqSpHHAgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklpncrcHMJEtu2s5M+af0+1hjFrfcQd1ewiSpAnCFRRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6axRQkhyb5JixGEiShUl6xqItSZK0fnMFZZAkvvVakqQuW6WAkuSvkyxNcm2SLw96bnaSK5rnv5PkKU35YysjSbZK0tccb5Tk6039M4CNVtL3K5Jc3fT946ZsiyTfbdq4IskuSTZI0pdk8wHn3prkaUm2TvKtJIuan32a549NsiDJecCXksxIcnHT39VJ9m7qbZDkU0muT3J2kh8kmds8NyfJhUkWJzk3yfRhrmNekt4kvSvuX74q0y9J0oQx6oCS5HnAPwEvrqpdgb8dVOVLwP+pql2AZcAHV9LkO4D7m/ofAeaM0PfWwOeAg5u+D2me+hfgmqaN9wFfqqpHge8Br2/OfT7QV1V3AycCx1fVHsDBwOcHdDMHeG1VvRn4JfCyqtodOBT4RFPnL4EZwCzg7cBeTR8bAicBc6tqDnBKc01PUFULqqqnqnombTxtJVMkSdLEtCrbGS8GvllVvwaoqt8kASDJNGDzqrqwqXsacOZK2tuP5h/+qlqaZOkIdV8AXFRVd/T33ZS/kE7QoKp+kmTLZixnAB8Avgi8sXkM8FLguf3jBp6cZLPm+KyqeqA53hD4ZJLZwApghwH9ndmEoP9JckFTviOwM3B+0/Yk4BcruX5JkjSMVQkoAWo1+niEP63UTB303GjbG67vDFFWwOXAnzcrL68D/rV5bgNgrwFBpNNIJ1T8YUDR/wfcDezanPPgCP31l19fVXut9EokSdJKrco9KD8G3pBkS+jc/9H/RFUtB36bZN+m6C1A/2pKH3/avpk7oL2LgMOatnYGdhmh78uBFyXZblDfA9vYH/h1Vf2uqgr4DvBx4Iaquqepfx7w7v5GmxWSoUwDftGslLyFzooIwCXAwc29KE8D9m/KbwK2TvLYlk+zJSZJklbDqFdQqur6JB8BLkyyAriGTvjodzjwmSQbA7cDRzblHwO+keQtwE8G1P808MVma2cJcNUIff8qyTzg20k2oLlHBDh2QBv3N2PodwawCDhiQNnRwMlN/cl0As5RQ3T5KeBbSQ4BLuBPqyvfAl4CXAfcDFwJLK+qPzY3y36i2WKaDJwAXD/cNUmSpOGls9ig0UqyaVX9vllJugrYp6r+Z3XamjJ9Zk0//ISxHeBa5LcZS5LWRJLFVTWqzzzzMz9W3dnNW5ifBHx4dcOJJEkaXusCSpIrgSmDit9SVcu6MZ7Bqmr/bo9BkqTxrnUBpaqe3+0xrCuztplGr9smkiQ9gR91L0mSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWmdytwcwkS27azkz5p+zTvvsO+6gddqfJEmrwxUUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOq0NKEl+kGTzVag/I8l1a3NMqyLJ77s9BkmS1letfZtxVb2y22OQJEnd0bUVlCTvTXJ0c3x8kp80xy9J8pUkfUm2alZGbkjyuSTXJzkvyUZN3TlJrk1yOfCulfT3vCRXJVmSZGmSmU3bNyY5rSn7ZpKNB7R9YZLFSc5NMr0p3z7JD5vyi5Ps1JRvl+TyJIuSfHgtTp0kSeNeN7d4LgL2bY57gE2TbAi8ELh4UN2ZwMlV9TzgXuDgpvyLwNFVtdco+jsKOLGqZjf93dmU7wgsqKpdgN8B72zGcRIwt6rmAKcAH2nqLwDe05QfA3yqKT8R+HRV7QH8z3CDSDIvSW+S3hX3Lx/FsCVJmni6GVAWA3OSbAY8BFxOJzjsyxMDyh1VtWTAeTOSTAM2r6oLm/Ivr6S/y4H3Jfk/wLOq6oGm/GdVdWlz/BU6AWlHYGfg/CRLgPcD2ybZFNgbOLMp/ywwvTl3H+BrKxtLVS2oqp6q6pm08bSVDFmSpImpa/egVNXDSfqAI4HLgKXAAcD2wA2Dqj804HgFsBEQoFahv68muRI4CDg3yduB24doo5q2rx+8MpPkycC9zSrMkN2MdjySJGl43X4Xz0V0tkkuorNqchSwpKpW+g99Vd0LLE/ywqbosJHqJ3k2cHtVfQI4C9ileeqZSfqDyJuAS4CbgK37y5NsmOR5VfU74I4khzTlSbJrc+6lwBtHMxZJkjSybgeUi+lskVxeVXcDD/LE7Z2RHAmc3Nwk+8BK6h4KXNdszewEfKkpvwE4PMlSYAs695H8EZgLfDTJtcASOls70Akfb2vKrwde25T/LfCuJIsA924kSVoDGcVixbiVZAZwdlXt3I3+p0yfWdMPP2Gd9um3GUuSuiXJ4qrqGU3dbq+gSJIkPUFrP6htdSV5OfDRQcV3VNXrB9etqj4679aRJEktMqG3eLqtp6enent7uz0MSZLWCbd4JEnSes2AIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWmdytwcwkS27azkz5p+zTvrqO+6gddKPJEljwRUUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOgYUSZLUOhMyoCT5fJLnNsfvW4N2fpBk8yHKj01yzJqMUZKkiWxCBpSqentV/VfzcLUDSlW9sqruHaNhSZKkxrgPKEk2SXJOkmuTXJfk0CQLk/QkOQ7YKMmSJKc39f8qyVVN2WeTTBqh7b4kWzXH/5TkpiQ/AnYc4Zx5SXqT9K64f/lYX64kSePCuA8owCuAn1fVrlW1M/DD/ieqaj7wQFXNrqrDkjwHOBTYp6pmAyuAw1bWQZI5wBuB3YC/BPYYrm5VLaiqnqrqmbTxtDW6MEmSxquJ8Emyy4CPJfkocHZVXZxkuLovAeYAi5o6GwG/HEUf+wLfqar7AZKctcajliRpAhv3AaWqbm5WOF4J/HuS80aoHuC0qvrH1elqtQYoSZKeYNxv8SR5OnB/VX0F+Biw+6AqDyfZsDn+MTA3yVObc7dI8qxRdHMR8PokGyXZDHj1GA1fkqQJadyvoACzgP9I8ijwMPAOOkGl3wJgaZKrm/tQ3g+cl2SDpv67gJ+O1EFVXZ3kDGBJU/fitXAdkiRNGKlyZ6JbpkyfWdMPP2Gd9OW3GUuSui3J4qrqGU3dcb/FI0mS1j8TYYtnjSW5EpgyqPgtVbWsG+ORJGm8M6CMQlU9f220O2ubafS69SJJ0hO4xSNJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklrHgCJJklpncrcHMJEtu2s5M+afMyZt9R130Ji0I0lSG7iCIkmSWseAIkmSWseAIkmSWseAIkmSWmfCB5QkC5P0rKTOEUk+ua7GJEnSRDfhA4okSWqf9S6gJHlvkqOb4+OT/KQ5fkmSryQ5MMnlSa5OcmaSTZvn5yS5MMniJOcmmT6o3Q2SnJbkX5vHRya5OcmFwD4D6r06yZVJrknyoyRPa869JcnWA9q6NclW62haJEkaV9a7gAJcBOzbHPcAmybZEHghsAx4P/DSqtod6AX+vnn+JGBuVc0BTgE+MqDNycDpwM1V9f4mvPwLnWDyMuC5A+peArygqnYDvg68t6oeBb4CHNbUeSlwbVX9emwvXZKkiWF9/KC2xcCcJJsBDwFX0wkq+wJn0QkTlyYBeBJwObAjsDNwflM+CfjFgDY/C3yjqvpDy/OBhVX1K4AkZwA7NM9tC5zRhJgnAXc05acA3wNOAN4KfHGowSeZB8wDmPTkrVd3DiRJGtfWuxWUqnoY6AOOBC4DLgYOALanExbOr6rZzc9zq+ptQIDrB5TPqqoDBzR7GXBAkqkDuxpmCCcBn6yqWcDfAFObcf0MuDvJi+kEnP8cZvwLqqqnqnombTxtteZAkqTxbr0LKI2LgGOa3xcDRwFLgCuAfZL8OUCSjZPsANwEbJ1kr6Z8wyTPG9DeF4AfAGcmmQxcCeyfZMtme+iQAXWnAXc1x4cPGtfn6Wz1fKOqVozZ1UqSNMGsrwHlYmA6cHlV3Q08CFzcbMkcAXwtyVI6gWWnqvojMBf4aJJr6YSZvQc2WFUfp7Nd9GXgbuBYOttDP2rK+x1LJ8hcDAy+x+QsYFOG2d6RJEmjk6rhdjK0qprPUzm+qvZdaWVgyvSZNf3wE8akb78sUJLUdkkWV9WInz3Wb328SbaVkswH3sGf3skjSZJW0/q6xdM6VXVcVT2rqi7p9lgkSVrfGVAkSVLruMXTRbO2mUav945IkvQErqBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWMaBIkqTWmdztAUxky+5azoz556xxO33HHTQGo5EkqT1cQZEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa3TyoCS5ENJXtrtcUiSpO5o3duMk0yqqg+shXYDpKoeHeu2JUnS2FqnKyhJZiS5MclpSZYm+WaSjZP0JflAkkuAQ5KpHtAsAAARoUlEQVScmmRuc05fkn9LcnmS3iS7Jzk3yW1JjmrqbJrkx0muTrIsyWsH9HdDkk8BVwP/nOT4AeP5X0k+PsJYrxvw+JgkxzbHRyf5r+Yavt6UbZLklCSLklzTP4Yh2p3XXEfvivuXj8W0SpI07nRjBWVH4G1VdWmSU4B3NuUPVtULAZK8YtA5P6uqvZpwcSqwDzAVuB74DPAg8Pqq+l2SrYArkpw1oL8jq+qdSTYBliZ5b1U9DBwJ/M1qXMN8YLuqeijJ5k3ZPwE/qaq3NmVXJflRVf1h4IlVtQBYADBl+sxajb4lSRr3unEPys+q6tLm+CvAC5vjM0Y4pz9sLAOurKr7qupXwINNGAjwb0mWAj8CtgGe1pzz06q6AqAJCz8BXpVkJ2DDqlq2GtewFDg9yV8BjzRlBwLzkywBFtIJUM9cjbYlSZrwurGCMnjVoP/xHwZXHOCh5vejA477H08GDgO2BuZU1cNJ+ugEhKHa/TzwPuBG4Isj9PkIjw9wUwccHwTsB7yGzrbR8+iEpIOr6qYR2pQkSaPQjRWUZybZqzl+E3DJGLQ5DfhlE04OAJ41XMWquhJ4BvBm4GsjtHk38NQkWyaZArwKIMkGwDOq6gLgvcDmwKbAucB7mptxSbLbml+WJEkTUzcCyg3A4c12zBbAp8egzdOBniS9dFZTblxJ/W8Al1bVb4er0Nyj8iHgSuDsAW1OAr6SZBlwDXB8Vd0LfBjYkM49Ltc1jyVJ0mpI1bq7TzPJDODsqtp5nXU69DjOphMsftzNcUyZPrOmH37CGrfjtxlLktYHSRZXVc9o6rbyg9rWliSbJ7kZeKDb4USSJA1vnd4kW1V9QNdWT5qtmB0GliXZEhgqrLykqu5ZJwOTJEmP07pPkl3XmhAyuxt9z9pmGr1uz0iS9AQTaotHkiStHwwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdQwokiSpdSZ3ewAT2bK7ljNj/jlr1EbfcQeN0WgkSWoPV1AkSVLrGFAkSVLrGFAkSVLrGFAkSVLrGFAkSVLrdDWgJLmsm/03Yzg1ydy10O7CJD1j3a4kSRNBVwNKVe3dzf4lSVI7dXsF5ffN7+lJLkqyJMl1SfYdpv4bkny8Of7bJLc3x9snuaQ5npPkwiSLk5ybZPqAOj9syi9OstMQ7X+4WVHZYIR2Fib5aJKrktzcP9YkGyX5epKlSc4ANloLUyZJ0oTQlntQ3gycW1WzgV2BJcPUuwjoDy/7Avck2QZ4IXBxkg2Bk4C5VTUHOAX4SFN/AfCepvwY4FMDG07yf4GnAkcCk0ZoB2ByVe0J/B3wwabsHcD9VbVLU3fOUBeQZF6S3iS9K+5fvpJpkSRpYmrLJ8kuAk5pAsZ3q2rIgFJV/5Nk0ySbAc8AvgrsRyesfBvYEdgZOD8JdILGL5JsCuwNnNmUA0wZ0PQ/A1dW1TyAJEO2M6D+t5vfi4EZzfF+wCeacS5NsnSYa1hAJywxZfrMGnFWJEmaoFoRUKrqoiT7AQcBX07yH1X1pWGqX05nleMm4GLgrcBewP8GnglcX1V7DTwhyZOBe5sVmqEsAuYk2aKqfgNkqHYGeKj5vYLHz6GBQ5KkMdCKLZ4kzwJ+WVWfA74A7D5C9YvobNFcBFwDHAA8VFXL6YSWrZPs1bS7YZLnVdXvgDuSHNKUJ8muA9r8IXAccE6zOjNkOyu5jIuAw5r6OwO7jH4GJEnSQK0IKMD+wJIk1wAHAyeOUPdiOts7F1XVCuBnwCUAVfVHYC7w0STX0rmXpf+dQocBb2vKrwdeO7DRqjoT+BxwFp0tneHaGc6ngU2brZ33Alet/LIlSdJQUuWuRLdMmT6zph9+whq14bcZS5LWF0kWV9WoPiOsLSsokiRJj2nFTbJDSXIlj3+nDcBbqmpZN8YjSZLWndYGlKp6frfHsLbN2mYavW7RSJL0BG7xSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1jGgSJKk1pnc7QHo8R5++GHuvPNOHnzwwW4PZUKYOnUq2267LRtuuGG3hyJJGsCA0jJ33nknm222GTNmzCBJt4czrlUV99xzD3feeSfbbbddt4cjSRrALZ6WefDBB9lyyy0NJ+tAErbccktXqySphQwoLWQ4WXeca0lqJwPKWpRkUrfHIEnS+sh7UBpJPgz8uqpObB5/BLgbmAK8ofn9nar6YPP8d4FnAFOBE6tqQVP+e+DjwMuB/w1csibjmjH/nDU5/Qn6jjtopXX23ntvLrvssjHtdyR9fX1cdtllvPnNb15nfUqS2s0VlD/5AnA4QJINgDfSCSgzgT2B2cCcJPs19d9aVXOAHuDoJFs25ZsA11XV86vqCeEkybwkvUl6f/WrX63dK1pN6zKcPPLII/T19fHVr351nfUpSWo/A0qjqvqAe5LsBhwIXAPsMeD4amAnOoEFOqHkWuAKOisp/eUrgG+N0M+Cquqpqp6tt956bVzKGtt0000BWLhwIS960Yt4wxvewA477MD8+fM5/fTT2XPPPZk1axa33XYbAEcccQRHHXUU++67LzvssANnn3020Lnh98gjj2TWrFnstttuXHDBBQCceuqpHHLIIbz61a/mwAMPZP78+Vx88cXMnj2b448/nr6+Pvbdd1923313dt9998cC08KFC9l///2ZO3cuO+20E4cddhhVBcCiRYvYe++92XXXXdlzzz257777WLFiBf/wD//AHnvswS677MJnP/vZdT2VkqTV5BbP430eOAL4M+AU4CXAv1fV4/5lS7I/8FJgr6q6P8lCOls9AA9W1Yp1NeC17dprr+WGG25giy224NnPfjZvf/vbueqqqzjxxBM56aSTOOGEE4DONs2FF17IbbfdxgEHHMCtt97KySefDMCyZcu48cYbOfDAA7n55psBuPzyy1m6dClbbLEFCxcu5GMf+9hjweb+++/n/PPPZ+rUqdxyyy286U1vore3F4BrrrmG66+/nqc//enss88+XHrppey5554ceuihnHHGGeyxxx787ne/Y6ONNuILX/gC06ZNY9GiRTz00EPss88+HHjggb6lWJLWAwaUx/sO8CFgQ+DNwCPAh5OcXlW/T7IN8DAwDfhtE052Al7QtRGvZXvssQfTp08HYPvtt+fAAw8EYNasWY+tiAC84Q1vYIMNNmDmzJk8+9nP5sYbb+SSSy7hPe95DwA77bQTz3rWsx4LKC972cvYYosthuzz4Ycf5t3vfjdLlixh0qRJj50DsOeee7LtttsCMHv2bPr6+pg2bRrTp09njz32AODJT34yAOeddx5Lly7lm9/8JgDLly/nlltuMaBI0nrAgDJAVf0xyQXAvc0qyHlJngNc3rwd9ffAXwE/BI5KshS4ic42z7g0ZcqUx4432GCDxx5vsMEGPPLII489N/jtukke234ZyiabbDLsc8cffzxPe9rTuPbaa3n00UeZOnXqY88NHM+kSZN45JFHqKoh3y5cVZx00km8/OUvH+EKJUlt5D0oAzQ3x76Azg2zAFTViVU1q/nZq6puq6qHquovqmqXqjqkqvavqoVN/U27NPyuOvPMM3n00Ue57bbbuP3229lxxx3Zb7/9OP300wG4+eab+e///m923HHHJ5y72Wabcd999z32ePny5UyfPp0NNtiAL3/5y6xYMfKO2U477cTPf/5zFi1aBMB9993HI488wstf/nI+/elP8/DDDz82hj/84Q9jdcmSpLXIFZRGkucCZ9N5K/Et3R5Pv9G8LbgNdtxxR170ohdx991385nPfIapU6fyzne+k6OOOopZs2YxefJkTj311MetgPTbZZddmDx5MrvuuitHHHEE73znOzn44IM588wzOeCAA0ZcbQF40pOexBlnnMF73vMeHnjgATbaaCN+9KMf8fa3v52+vj523313qoqtt96a7373u2trCiRJYygjLcNr7erp6an+mz/73XDDDTznOc/p0ohWzxFHHMGrXvUq5s6d2+2hrJb1cc4laX2UZHFV9Yymrls8kiSpddzi0Ro79dRTuz0ESdI44wpKC7nttu4415LUTgaUlpk6dSr33HOP/3CuA1XFPffc87i3MUuS2sEtnpbZdtttufPOO2nr9/SMN1OnTn3sg98kSe1hQGmZDTfc0E86lSRNeG7xSJKk1jGgSJKk1jGgSJKk1vGTZLsoyX10vmxQa24r4NfdHsQ44DyOHedy7DiXY6fbc/msqtp6NBW9Sba7bhrtR/5qZEl6ncs15zyOHedy7DiXY2d9mku3eCRJUusYUCRJUusYULprQbcHMI44l2PDeRw7zuXYcS7Hznozl94kK0mSWscVFEmS1DoGFEmS1DoGlC5J8ookNyW5Ncn8bo+nLZL0JVmWZEmS3qZsiyTnJ7ml+f2UpjxJPtHM4dIkuw9o5/Cm/i1JDh9QPqdp/9bm3Kz7q1w7kpyS5JdJrhtQttbnbrg+1lfDzOOxSe5qXpdLkrxywHP/2MzJTUlePqB8yL/xJNslubKZrzOSPKkpn9I8vrV5fsa6ueK1J8kzklyQ5IYk1yf526bc1+UqGGEex/frsqr8Wcc/wCTgNuDZwJOAa4HndntcbfgB+oCtBpX9X2B+czwf+Ghz/ErgP4EALwCubMq3AG5vfj+lOX5K89xVwF7NOf8J/EW3r3kM524/YHfgunU5d8P1sb7+DDOPxwLHDFH3uc3f7xRgu+bvetJIf+PAN4A3NsefAd7RHL8T+Exz/EbgjG7PxRjM5XRg9+Z4M+DmZs58XY7NPI7r16UrKN2xJ3BrVd1eVX8Evg68tstjarPXAqc1x6cBrxtQ/qXquALYPMl04OXA+VX1m6r6LXA+8IrmuSdX1eXV+Wv70oC21ntVdRHwm0HF62LuhutjvTTMPA7ntcDXq+qhqroDuJXO3/eQf+PNf92/GPhmc/7g/0365/GbwEvW9xW+qvpFVV3dHN8H3ABsg6/LVTLCPA5nXLwuDSjdsQ3wswGP72TkF9tEUsB5SRYnmdeUPa2qfgGdP1TgqU35cPM4UvmdQ5SPZ+ti7obrY7x5d7PtcMqA7YJVncctgXur6pFB5Y9rq3l+eVN/XGi2BnYDrsTX5WobNI8wjl+XBpTuGCp9+n7vjn2qanfgL4B3JdlvhLrDzeOqlk9Ezt2q+TSwPTAb+AXw/zflYzmP43aOk2wKfAv4u6r63UhVhyjzddkYYh7H9evSgNIddwLPGPB4W+DnXRpLq1TVz5vfvwS+Q2dJ8u5mKZfm9y+b6sPN40jl2w5RPp6ti7kbro9xo6rurqoVVfUo8Dk6r0tY9Xn8NZ1ti8mDyh/XVvP8NEa/1dRaSTak84/q6VX17abY1+UqGmoex/vr0oDSHYuAmc1d00+ic+PRWV0eU9cl2STJZv3HwIHAdXTmpv+u/cOB7zXHZwF/3dz5/wJgebOUey5wYJKnNEueBwLnNs/dl+QFzR7qXw9oa7xaF3M3XB/jRv8/dI3X03ldQufa39i802E7YCadmzaH/Btv7pO4AJjbnD/4f5P+eZwL/KSpv95qXitfAG6oqo8PeMrX5SoYbh7H/etyXd2N688T7rJ+JZ07sW8D/qnb42nDD507y69tfq7vnxc6+50/Bm5pfm/RlAc4uZnDZUDPgLbeSufGsFuBIweU99D5I74N+CTNpymPhx/ga3SWeR+m8189b1sXczdcH+vrzzDz+OVmnpbS+T/s6QPq/1MzJzcx4F1hw/2NN6/zq5r5PROY0pRPbR7f2jz/7G7PxRjM5QvpbAcsBZY0P6/0dTlm8ziuX5d+1L0kSWodt3gkSVLrGFAkSVLrGFAkSVLrGFAkSVLrGFAkSVLrGFAkSVLrGFAkSVLr/D+0l6OU8tFtjwAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAigAAAHVCAYAAADb6QDfAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucXVV9///XmyQk3Axy0Y5gDdIAKoFABjQgCopoxUstQURqAes3X7zR/vq1NlqrVmuL3/oTEPESFfGCinhFsAIq4X7JBEJC5CYyVtCCUokoFyF8vn+cPTBMJpNJmOTsZF7Px2Mes886a6+19tJj3q6195lUFZIkSW2ySbcHIEmSNJQBRZIktY4BRZIktY4BRZIktY4BRZIktY4BRZIktY4BRZIktY4BRZIktY4BRZIktc7Ebg9gPNtuu+1q2rRp3R6GJEnrxaJFi35TVduPpq4BpYumTZtGX19ft4chSdJ6keTno63rFo8kSWodA4okSWodA4okSWodA4okSWodA4okSWodA4okSWodHzPuoqV3LGfavHO7PQxJklbSf8KhXe3fFRRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6rQkoSaYluX4N6r8qybzm+P1J3jFSm0l6k3xs7EYMST6Q5OBhyg9Mcs5Y9iVJ0niywT5mXFVnA2evQf0+YEz/dHBVvXcs25MkSR2tWUFpTEzyhSRLknwjyeZJ+pNsB4+ugixojo9J8vGhDSSZleS6JFcAbx1U/uiqRrPiclqSBUl+luT4QfX+OcmNSS5I8tXhVmYG1T09yZzm+GXNeZcCfzlG8yFJ0rjUtoCyKzC/qvYAfge8ZS3a+DxwfFXNXk293YCXAvsC70syKUkvcBiwF52Q0TuaDpNMAT4DvBI4APiTEerOTdKXpG/FfctH07wkSeNO2wLKL6rqsub4y8Dz1+TkJFOBravqoqboSyNUP7eqHqyq3wB3AU9t+vtuVd1fVfcC3xtl17sBt1XVLVVVzdiHVVXzq6q3qnonbD51lM1LkjS+tC2g1DCvH+axcU5ZzfkZpo1VeXDQ8Qo69+NklOcOZ7T9SpKk1WhbQPnTJANbM0cClwL9wKym7LCRTq6qe4DlSQZWXo5aw/4vBV6ZZEqSLYHR/iGCG4GdkuzcvD5yDfuVJEmDtC2g3AAcnWQJsA3wSeBfgJOTXEJnpWN1jgVObW6SvX9NOq+qhXSeDLoO+Badp35We6NIVT0AzAXObW6S/fma9CtJkh4vnVsmNCDJllX1+ySbAxcDc6vqmnXR1+Se6dVz9EnromlJkp6QdfHXjJMsqqpRPYCywX4Pyjo0P8mz6dzv8oV1FU4kSdKqGVCGqKrXDy1Lciqw/5Dik6vq8+tnVJIkjS9u8XRRb29v9fWN6ZfbSpLUWmuyxdO2m2QlSZIMKJIkqX0MKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUmdnsA49nSO5Yzbd653R6GJGkD0n/Cod0ewnrhCookSWodA4okSWodA4okSWodA4okSWodA4okSWqdtQ4oSb6fZOs1qP/ute1rfUpyTJKndXsckiSNZ2sdUKrq5VV1z+CydKyqzdYElCQTRnj7GGCNAkoSH9eWJGkMjSqgJPlOkkVJliWZ25T1J9kuybQkNyT5BHAN8PRhzj8B2CzJ4iRnNGV/leTqpuzTA6Ehye+TfLjp74dJ9k2yIMnPkryqqXNMku8m+UGSm5K8b1BfI7X7gSRXAbOTvDfJwiTXJ5nfhKs5QC9wRnP+ZgPX2bTRm2RBc/z+5rzzgS8mmZDkP5o2lyT536uYy7lJ+pL0rbhv+WimX5KkcWe0KyhvrKpZdP7xPj7JtkPe3xX4YlXtVVU/H3pyVc0D7q+qmVV1VJJnAUcA+1fVTGAFcFRTfQtgQdPfvcC/Ai8BXgN8YFCz+zbnzAQOb8LD6tq9vqqeW1WXAh+vqn2qandgM+AVVfUNoA84qhnr/auZl1nAq6vq9cDfAMurah9gH+B/JdlpmLmYX1W9VdU7YfOpq2lekqTxabRbE8cneU1z/HRg+pD3f15VV65Bvy+m84/7wiTQCQh3Ne/9EfhBc7wUeLCqHkqyFJg2qI0LqupugCTfAp4PPDxCuyuAbw46/6Ak7wQ2B7YBlgHfW4NrADh7UIg5BNijWYUBmEpnnm5bwzYlSRr3VhtQkhwIHAzMrqr7mi2OKUOq/WEN+w3whap61zDvPVRV1Rw/AjwIUFWPDLnXo4acV6tp94GqWgGQZArwCaC3qn6R5P2sfE0DHuaxlaaRrjvA26vqvFW0I0mSRmk0WzxTgd824WQ34Hlr2ddDSSY1xz8C5iR5CkCSbZI8Yw3be0lz3mbAXwCXrUG7A0HjN0m2BOYMeu9eYKtBr/vprMoAHDbCeM4D3jxwjUl2SbLFGl6TJElidFs8PwCOS7IEuAlYk62cweYDS5Jc09yH8h7g/Oapn4eAtwIr3b8ygkuBLwF/BnylqvoARtNuVd2T5DN0tpD6gYWD3j4d+FSS+4HZwL8An2sek75qhPF8ls4W1DXp7C/9mk5wkiRJayiP7aZsOJIcQ2d75m3dHssTMblnevUcfVK3hyFJ2oBsyH/NOMmiquodTV2/SVaSJLXOmK+gNN8zMnlI8RuqaumYdrQR6O3trb6+vm4PQ5Kk9WJNVlDG/BtQq+q5Y92mJEkaX9zikSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrTOx2wMYz5besZxp887t9jAkaY30n3Bot4egccAVFEmS1DoGFEmS1DoGFEmS1DoGFEmS1DobZEBJMi3J9eug3Q8kOXiY8gOTnNMcvyrJvOb4L5I8e6zHIUnSeOdTPINU1XtHUeds4Ozm5V8A5wA/WZfjkiRpvNkgV1AaE5J8JsmyJOcn2SzJgiS9AEm2S9LfHB+T5DtJvpfktiRvS/L3Sa5NcmWSbZp6pyeZ0xy/LMmNSS4F/nKg06atjyfZD3gV8B9JFifZOck1g+pNT7Jo/U2HJEkbjw05oEwHTq2q5wD3AIetpv7uwOuBfYEPAfdV1V7AFcBfD66YZArwGeCVwAHAnwxtrKoup7OS8g9VNbOqbgWWJ5nZVDkWOH3tLk2SpPFtQw4ot1XV4uZ4ETBtNfUvrKp7q+rXwHLge0350mHO3a1p/5aqKuDLoxzTZ4Fjk0wAjgC+MrRCkrlJ+pL0rbhv+SiblSRpfNmQA8qDg45X0Lmf5mEeu6YpI9R/ZNDrRxj+XpxaizF9E/hz4BXAoqq6e6VGq+ZXVW9V9U7YfOpadCFJ0sZvQw4ow+kHZjXHc55AOzcCOyXZuXl95Crq3QtsNfCiqh4AzgM+CXz+CfQvSdK4trEFlI8Ab05yObDd2jbSBI25wLnNTbI/X0XVrwH/0NxsOxBmzqCz+nL+2vYvSdJ4l84tFhorSd4BTK2qf15d3ck906vn6JPWw6gkaez4xwK1tpIsqqre0dT1e1DGUJJvAzsDL+r2WCRJ2pAZUMZQVb2m22OQJGljsLHdgyJJkjYCrqB00YwdptLnXq4kSStxBUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLWOAUWSJLXOxG4PYDxbesdyps07t9vDkAT0n3Bot4cgaRBXUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUus8oYCS5P1J3jEWA0myIEnvWLQlSZI2bK6gDJHER68lSeqyNQooSf46yZIk1yX50pD3Zia5snn/20me3JQ/ujKSZLsk/c3xZkm+1tQ/E9hsNX2/LMk1Td8/asq2SfKdpo0rk+yRZJMk/Um2HnTuT5M8Ncn2Sb6ZZGHzs3/z/vuTzE9yPvDFJNOSXNL0d02S/Zp6myT5RJJlSc5J8v0kc5r3ZiW5KMmiJOcl6VnFdcxN0pekb8V9y9dk+iVJGjdGHVCSPAf4J+BFVbUn8LdDqnwR+Meq2gNYCrxvNU2+Gbivqf8hYNYIfW8PfAY4rOn78OatfwGubdp4N/DFqnoE+C7wmubc5wL9VXUncDJwYlXtAxwGfHZQN7OAV1fV64G7gJdU1d7AEcDHmjp/CUwDZgBvAmY3fUwCTgHmVNUs4LTmmlZSVfOrqreqeidsPnU1UyRJ0vi0JtsZLwK+UVW/Aaiq/0kCQJKpwNZVdVFT9wvAWatp7wU0//BX1ZIkS0ao+zzg4qq6baDvpvz5dIIGVfXjJNs2YzkTeC/weeB1zWuAg4FnD4wbeFKSrZrjs6vq/uZ4EvDxJDOBFcAug/o7qwlB/53kwqZ8V2B34IKm7QnAr1Zz/ZIkaRXWJKAEqLXo42EeW6mZMuS90ba3qr4zTFkBVwB/1qy8/AXwr817mwCzBwWRTiOdUPGHQUX/H3AnsGdzzgMj9DdQvqyqZq/2SiRJ0mqtyT0oPwJem2Rb6Nz/MfBGVS0HfpvkgKboDcDAako/j23fzBnU3sXAUU1buwN7jND3FcALk+w0pO/BbRwI/KaqfldVBXwb+ChwQ1Xd3dQ/H3jbQKPNCslwpgK/alZK3kBnRQTgUuCw5l6UpwIHNuU3AdsneXTLp9kSkyRJa2HUKyhVtSzJh4CLkqwArqUTPgYcDXwqyebAz4Bjm/KPAF9P8gbgx4PqfxL4fLO1sxi4eoS+f51kLvCtJJvQ3CMCvH9QG/c1YxhwJrAQOGZQ2fHAqU39iXQCznHDdPkJ4JtJDgcu5LHVlW8CLwauB24GrgKWV9Ufm5tlP9ZsMU0ETgKWreqaJEnSqqWz2KDRSrJlVf2+WUm6Gti/qv57bdqa3DO9eo4+aWwHKGmt+NeMpXUvyaKqGtV3nvmdH2vunOYR5k2BD65tOJEkSavWuoCS5Cpg8pDiN1TV0m6MZ6iqOnCs2pqxw1T6/H9tkiStpHUBpaqe2+0xSJKk7vKr7iVJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUusYUCRJUutM7PYAxrOldyxn2rxzuz0MaUz0n3Bot4cgaSPiCookSWodA4okSWodA4okSWodA4okSWodA4okSWqd1gaUJN9PsvUa1J+W5Pp1OaY1keT33R6DJEkbqtY+ZlxVL+/2GCRJUnd0bQUlyTuTHN8cn5jkx83xi5N8OUl/ku2alZEbknwmybIk5yfZrKk7K8l1Sa4A3rqa/p6T5Ooki5MsSTK9afvGJF9oyr6RZPNBbV+UZFGS85L0NOU7J/lBU35Jkt2a8p2SXJFkYZIPjjCOuUn6kvStuG/5mMylJEkbm25u8VwMHNAc9wJbJpkEPB+4ZEjd6cCpVfUc4B7gsKb888DxVTV7FP0dB5xcVTOb/m5vyncF5lfVHsDvgLc04zgFmFNVs4DTgA819ecDb2/K3wF8oik/GfhkVe0D/PeqBlFV86uqt6p6J2w+dRTDliRp/OlmQFkEzEqyFfAgcAWd4HAAKweU26pq8aDzpiWZCmxdVRc15V9aTX9XAO9O8o/AM6rq/qb8F1V1WXP8ZToBaVdgd+CCJIuB9wA7JtkS2A84qyn/NNDTnLs/8NVRjkWSJI2ga/egVNVDSfqBY4HLgSXAQcDOwA1Dqj846HgFsBkQoNagv68kuQo4FDgvyZuAnw3TRjVtLxu6MpPkScA9zSrMsN2MdjySJGnVuv0Uz8V0tkkuprNqchywuKpW+w99Vd0DLE/y/KboqJHqJ3km8LOq+hhwNrBH89afJhkIIkcClwI3AdsPlCeZlOQ5VfU74LYkhzflSbJnc+5lwOtGMxZJkjSybgeUS+hskVxRVXcCD7Dy9s5IjgVObW6SvX81dY8Arm+2ZnYDvtiU3wAcnWQJsA2d+0j+CMwBPpzkOmAxna0d6ISPv2nKlwGvbsr/FnhrkoWAN5dIkvQEZBSLFRutJNOAc6pq9270P7lnevUcfVI3upbGnH/NWNLqJFlUVb2jqdvtFRRJkqSVtPaL2tZWkpcCHx5SfFtVvWZo3arqp/O0jiRJapFxvcXTbb29vdXX19ftYUiStF64xSNJkjZoBhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6E7s9gPFs6R3LmTbv3G4PQ1oj/Scc2u0hSBoHXEGRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0BpJJmW5Ppuj0OSJBlQ1qkkPiUlSdJaMKA83oQkn0myLMn5STZLMjPJlUmWJPl2kicDJFmQpLc53i5Jf3N8TJKzknwPOL97lyJJ0obLgPJ404FTq+o5wD3AYcAXgX+sqj2ApcD7RtHObODoqnrR0DeSzE3Sl6RvxX3Lx3DokiRtPAwoj3dbVS1ujhcBOwNbV9VFTdkXgBeMop0Lqup/hnujquZXVW9V9U7YfOoTH7EkSRshA8rjPTjoeAWw9Qh1H+ax+Zsy5L0/jOWgJEkabwwoI1sO/DbJAc3rNwADqyn9wKzmeM56HpckSRs1nzJZvaOBTyXZHPgZcGxT/hHg60neAPy4W4OTJGljlKrq9hjGrck906vn6JO6PQxpjfjHAiWtrSSLqqp3NHXd4pEkSa1jQJEkSa3jPShdNGOHqfS5XC5J0kpcQZEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa0zsdsDGM+W3rGcafPO7fYwpNXqP+HQbg9B0jjjCookSWodA4okSWodA4okSWodA4okSWodA4okSWodA8owkuyWZHGSa5PsvBbn/12SzdfF2CRJGg/GRUBJMmENT/kL4LtVtVdV3boWXf4dYECRJGktdS2gJNkiyblJrktyfZIjkrwsyY1JLk3ysSTnNHXfn+Qdg869Psm05vg7SRYlWZZk7qA6v0/ygSRXAbOTzEpyUVP3vCQ9qxjXy+kEjDclubAp+6skVzerKp8eCDxJDklyRZJrkpyVZMskxwNPAy4cOH9I+3OT9CXpW3Hf8jGaTUmSNi7dXEF5GfDLqtqzqnYHfgB8BnglcADwJ6Ns541VNQvoBY5Psm1TvgVwfVU9F7gKOAWY09Q9DfjQcI1V1feBTwEnVtVBSZ4FHAHsX1UzgRXAUUm2A94DHFxVewN9wN9X1ceAXwIHVdVBw7Q/v6p6q6p3wuZTR3mJkiSNL938JtmlwEeSfBg4B7gXuK2qbgFI8mVg7gjnDzg+yWua46cD04G76QSJbzbluwK7AxckAZgA/GqU43wxMAtY2Jy7GXAX8Dzg2cBlTfmmwBWjbFOSJI2gawGlqm5OMgt4OfDvwPlAraL6wzx+tWcKQJIDgYOB2VV1X5IFA+8BD1TViuY4wLKqmr0WQw3whap61+MKk1cCF1TVkWvRpiRJGkE370F5GnBfVX0Z+AiwH7DToKdmBv/D3w/s3Zy3N7BTUz4V+G0TTnajs6oxnJuA7ZPMbtqYlOQ5oxzqj4A5SZ7SnLtNkmcAVwL7J/mzpnzzJLs059wLbDXK9iVJ0hDd3OKZAfxHkkeAh4A3A9sB5yb5DXApnW0Z6GzV/HWSxcBC4Oam/AfAcUmW0AkhVw7XUVX9Mckc4GNJptK57pOAZasbZFX9JMl7gPOTbNKM9a1VdWWSY4CvJpncVH9PM7b5wH8m+dVw96FIkqSRpWpVuyrd1WzfvKOqXtHtsawrk3umV8/RJ3V7GNJq+deMJY2FJIuqqnc0dcfF96BIkqQNSze3eEZUVQuABeuyjySnAvsPKT65qj6/LvuVJEkja+0Wz3jQ29tbfX193R6GJEnrhVs8kiRpg2ZAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrTOx2wMYz5besZxp887t9jC0keo/4dBuD0GS1porKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKJIkqXUMKKOUZGaSl3d7HJIkjQcGlNGbCaxRQEniY9ySJK2FcRVQkkxLcmOSzya5PskZSQ5OclmSW5Lsm2SLJKclWZjk2iSvTrIp8AHgiCSLkxwxXL2mj2OSnJXke8D5Xb1gSZI2UOPx/+H/GXA4MBdYCLweeD7wKuDdwE+AH1fVG5NsDVwN/BB4L9BbVW8DSPJvQ+sl+WHTx2xgj6r6n6GdJ5nb9M2EJ22/7q5SkqQN2HgMKLdV1VKAJMuAH1VVJVkKTAN2BF6V5B1N/SnAnw7TziEj1LtguHACUFXzgfkAk3um1xhcjyRJG53xGFAeHHT8yKDXj9CZjxXAYVV10+CTkjx3SDsZod4fxnTEkiSNM+PqHpRROg94e5IAJNmrKb8X2GoU9SRJ0hNkQFnZB4FJwJIk1zevAS4Enj1wk+wI9SRJ0hOUKm+D6JbJPdOr5+iTuj0MbaT8a8aS2ibJoqrqHU1dV1AkSVLrGFAkSVLrjMeneFpjxg5T6XMZXpKklbiCIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWseAIkmSWmditwcwni29YznT5p3b7WFoiP4TDu32ECRp3HMFRZIktY4BRZIktY4BRZIktY4BRZIktc64DyhJFiTpXU2dY5J8fH2NSZKk8W7cBxRJktQ+G1xASfLOJMc3xycm+XFz/OIkX05ySJIrklyT5KwkWzbvz0pyUZJFSc5L0jOk3U2SfCHJvzavj01yc5KLgP0H1XtlkquSXJvkh0me2px7S5LtB7X10yTbradpkSRpo7LBBRTgYuCA5rgX2DLJJOD5wFLgPcDBVbU30Af8ffP+KcCcqpoFnAZ8aFCbE4EzgJur6j1NePkXOsHkJcCzB9W9FHheVe0FfA14Z1U9AnwZOKqpczBwXVX9Zujgk8xN0pekb8V9y5/oXEiStFHaEL+obREwK8lWwIPANXSCygHA2XTCxGVJADYFrgB2BXYHLmjKJwC/GtTmp4GvV9VAaHkusKCqfg2Q5Exgl+a9HYEzmxCzKXBbU34a8F3gJOCNwOeHG3xVzQfmA0zumV5rOwmSJG3MNriAUlUPJekHjgUuB5YABwE70wkLF1TVkYPPSTIDWFZVs1fR7OXAQUn+/6p6YKCrVdQ9BfhoVZ2d5EDg/c24fpHkziQvohNwjlrF+ZIkaTU2xC0e6GzzvKP5fQlwHLAYuBLYP8mfASTZPMkuwE3A9klmN+WTkjxnUHufA74PnJVkInAVcGCSbZvtocMH1Z0K3NEcHz1kXJ+ls9Xz9apaMWZXK0nSOLOhBpRLgB7giqq6E3gAuKTZkjkG+GqSJXQCy25V9UdgDvDhJNfRCTP7DW6wqj5KZ7voS8CddFZGrgB+2JQPeD+dIHMJMPQek7OBLVnF9o4kSRqdVHkbxFhpvk/lxKo6YLWV6dyD0nP0Set4VFpT/rFASVo3kiyqqhG/e2zABncPSlslmQe8Ge89kSTpCdtQt3hap6pOqKpnVNWl3R6LJEkbOgOKJElqHbd4umjGDlPp834HSZJW4gqKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqHQOKJElqnYndHsB4tvSO5Uybd263hzFm+k84tNtDkCRtJFxBkSRJrWNAkSRJrWNAkSRJrWNAkSRJrWNAkSRJrdPVgJLk8m7234zh9CRz1kG7C5L0jnW7kiSNB10NKFW1Xzf7lyRJ7dTtFZTfN797klycZHGS65McsIr6r03y0eb4b5P8rDneOcmlzfGsJBclWZTkvCQ9g+r8oCm/JMluw7T/wWZFZZMR2lmQ5MNJrk5y88BYk2yW5GtJliQ5E9hsHUyZJEnjQlvuQXk9cF5VzQT2BBavot7FwEB4OQC4O8kOwPOBS5JMAk4B5lTVLOA04ENN/fnA25vydwCfGNxwkv8LPAU4FpgwQjsAE6tqX+DvgPc1ZW8G7quqPZq6s4a7gCRzk/Ql6Vtx3/LVTIskSeNTW75JdiFwWhMwvlNVwwaUqvrvJFsm2Qp4OvAV4AV0wsq3gF2B3YELkkAnaPwqyZbAfsBZTTnA5EFN/zNwVVXNBUgybDuD6n+r+b0ImNYcvwD4WDPOJUmWrOIa5tMJS0zumV4jzookSeNUKwJKVV2c5AXAocCXkvxHVX1xFdWvoLPKcRNwCfBGYDbwf4A/BZZV1ezBJyR5EnBPs0IznIXArCTbVNX/ABmunUEebH6v4PFzaOCQJGkMtGKLJ8kzgLuq6jPA54C9R6h+MZ0tmouBa4GDgAerajmd0LJ9ktlNu5OSPKeqfgfcluTwpjxJ9hzU5g+AE4Bzm9WZYdtZzWVcDBzV1N8d2GP0MyBJkgZrRUABDgQWJ7kWOAw4eYS6l9DZ3rm4qlYAvwAuBaiqPwJzgA8nuY7OvSwDTwodBfxNU74MePXgRqvqLOAzwNl0tnRW1c6qfBLYstnaeSdw9eovW5IkDSdV7kp0y+Se6dVz9EndHsaY8a8ZS5JGkmRRVY3qO8LasoIiSZL0qFbcJDucJFfx+CdtAN5QVUu7MR5JkrT+tDagVNVzuz2GdW3GDlPpc1tEkqSVuMUjSZJax4AiSZJax4CwBZdHAAASQUlEQVQiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJax4AiSZJaZ2K3BzCeLb1jOdPmnbte++w/4dD12p8kSWvDFRRJktQ6BhRJktQ6BhRJktQ6BhRJktQ6BhRJktQ64zKgJPlskmc3x+9+Au18P8nWw5S/P8k7nsgYJUkaz8ZlQKmqN1XVT5qXax1QqurlVXXPGA1LkiQ1NvqAkmSLJOcmuS7J9UmOSLIgSW+SE4DNkixOckZT/6+SXN2UfTrJhBHa7k+yXXP8T0luSvJDYNcRzpmbpC9J34r7lo/15UqStFHY6AMK8DLgl1W1Z1XtDvxg4I2qmgfcX1Uzq+qoJM8CjgD2r6qZwArgqNV1kGQW8DpgL+AvgX1WVbeq5ldVb1X1Tth86hO6MEmSNlbj4ZtklwIfSfJh4JyquiTJquq+GJgFLGzqbAbcNYo+DgC+XVX3ASQ5+wmPWpKkcWyjDyhVdXOzwvFy4N+TnD9C9QBfqKp3rU1XazVASZK0ko1+iyfJ04D7qurLwEeAvYdUeSjJpOb4R8CcJE9pzt0myTNG0c3FwGuSbJZkK+CVYzR8SZLGpY1+BQWYAfxHkkeAh4A30wkqA+YDS5Jc09yH8h7g/CSbNPXfCvx8pA6q6pokZwKLm7qXrIPrkCRp3EiVOxPdMrlnevUcfdJ67dO/ZixJ6pYki6qqdzR1N/otHkmStOEZD1s8T1iSq4DJQ4rfUFVLuzEeSZI2dgaUUaiq566LdmfsMJU+t1wkSVqJWzySJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1DCiSJKl1JnZ7AOPZ0juWM23eueu8n/4TDl3nfUiSNJZcQZEkSa1jQJEkSa1jQJEkSa1jQJEkSa1jQJEkSa3TyoCS5ANJDu72OCRJUne07jHjJBOq6r3roN0AqapHxrptSZI0ttbrCkqSaUluTPKFJEuSfCPJ5kn6k7w3yaXA4UlOTzKnOac/yb8luSJJX5K9k5yX5NYkxzV1tkzyoyTXJFma5NWD+rshySeAa4B/TnLioPH8ryQfHWGs1w96/Y4k72+Oj0/yk+YavtaUbZHktCQLk1w7MAZJkrTmurGCsivwN1V1WZLTgLc05Q9U1fMBkrxsyDm/qKrZTbg4HdgfmAIsAz4FPAC8pqp+l2Q74MokZw/q79iqekuSLYAlSd5ZVQ8BxwL/ey2uYR6wU1U9mGTrpuyfgB9X1RubsquT/LCq/jD4xCRzgbkAE560/Vp0LUnSxq8b96D8oqoua46/DDy/OT5zhHMGwsZS4Kqqureqfg080ISBAP+WZAnwQ2AH4KnNOT+vqisBmrDwY+AVSXYDJlXV0rW4hiXAGUn+Cni4KTsEmJdkMbCAToD606EnVtX8quqtqt4Jm09di64lSdr4dWMFpVbx+g9DKw7yYPP7kUHHA68nAkcB2wOzquqhJP10AsJw7X4WeDdwI/D5Efp8mMcHuCmDjg8FXgC8is620XPohKTDquqmEdqUJEmj0I0VlD9NMrs5PhK4dAzanArc1YSTg4BnrKpiVV0FPB14PfDVEdq8E3hKkm2TTAZeAZBkE+DpVXUh8E5ga2BL4Dzg7c3NuCTZ64lfliRJ41M3AsoNwNHNdsw2wCfHoM0zgN4kfXRWU25cTf2vA5dV1W9XVaG5R+UDwFXAOYPanAB8OclS4FrgxKq6B/ggMInOPS7XN68lSdJaSNXQHZd12FkyDTinqnZfb50OP45z6ASLH3VzHJN7plfP0Set8378a8aSpDZIsqiqekdTt5Vf1LauJNk6yc3A/d0OJ5IkadXW602yVdUPdG31pNmK2WVwWZJtgeHCyour6u71MjBJkvQ4rfsm2fWtCSEzu9H3jB2m0uf2iyRJKxlXWzySJGnDYECRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtY0CRJEmtM7HbA9DjPfTQQ9x+++088MAD3R7KuDBlyhR23HFHJk2a1O2hSJIGMaC0zO23385WW23FtGnTSNLt4WzUqoq7776b22+/nZ122qnbw5EkDeIWT8s88MADbLvttoaT9SAJ2267ratVktRCBpQWMpysP861JLWTAWUdSjKh22OQJGlD5D0ojSQfBH5TVSc3rz8E3AlMBl7b/P52Vb2vef87wNOBKcDJVTW/Kf898FHgpcD/AS59IuOaNu/cJ3L6SvpPOHS1dfbbbz8uv/zyMe13JP39/Vx++eW8/vWvX299SpLazRWUx3wOOBogySbA6+gElOnAvsBMYFaSFzT131hVs4Be4Pgk2zblWwDXV9Vzq2qlcJJkbpK+JH2//vWv1+0VraX1GU4efvhh+vv7+cpXvrLe+pQktZ8BpVFV/cDdSfYCDgGuBfYZdHwNsBudwAKdUHIdcCWdlZSB8hXAN0foZ35V9VZV7/bbb78uLuUJ23LLLQFYsGABL3zhC3nta1/LLrvswrx58zjjjDPYd999mTFjBrfeeisAxxxzDMcddxwHHHAAu+yyC+eccw7QueH32GOPZcaMGey1115ceOGFAJx++ukcfvjhvPKVr+SQQw5h3rx5XHLJJcycOZMTTzyR/v5+DjjgAPbee2/23nvvRwPTggULOPDAA5kzZw677bYbRx11FFUFwMKFC9lvv/3Yc8892Xfffbn33ntZsWIF//AP/8A+++zDHnvswac//en1PZWSpLXkFs/jfRY4BvgT4DTgxcC/V9Xj/mVLciBwMDC7qu5LsoDOVg/AA1W1Yn0NeF277rrruOGGG9hmm2145jOfyZve9CauvvpqTj75ZE455RROOukkoLNNc9FFF3Hrrbdy0EEH8dOf/pRTTz0VgKVLl3LjjTdyyCGHcPPNNwNwxRVXsGTJErbZZhsWLFjARz7ykUeDzX333ccFF1zAlClTuOWWWzjyyCPp6+sD4Nprr2XZsmU87WlPY//99+eyyy5j33335YgjjuDMM89kn3324Xe/+x2bbbYZn/vc55g6dSoLFy7kwQcfZP/99+eQQw7xkWJJ2gAYUB7v28AHgEnA64GHgQ8mOaOqfp9kB+AhYCrw2yac7AY8r2sjXsf22Wcfenp6ANh555055JBDAJgxY8ajKyIAr33ta9lkk02YPn06z3zmM7nxxhu59NJLefvb3w7AbrvtxjOe8YxHA8pLXvISttlmm2H7fOihh3jb297G4sWLmTBhwqPnAOy7777suOOOAMycOZP+/n6mTp1KT08P++yzDwBPetKTADj//PNZsmQJ3/jGNwBYvnw5t9xyiwFFkjYABpRBquqPSS4E7mlWQc5P8izgiuZx1N8DfwX8ADguyRLgJjrbPBulyZMnP3q8ySabPPp6k0024eGHH370vaGP6yZ5dPtlOFtsscUq3zvxxBN56lOfynXXXccjjzzClClTHn1v8HgmTJjAww8/TFUN+7hwVXHKKafw0pe+dIQrlCS1kfegDNLcHPs8OjfMAlBVJ1fVjOZndlXdWlUPVtWfV9UeVXV4VR1YVQua+lt2afhdddZZZ/HII49w66238rOf/Yxdd92VF7zgBZxxxhkA3HzzzfzXf/0Xu+6660rnbrXVVtx7772Pvl6+fDk9PT1ssskmfOlLX2LFipF3zHbbbTd++ctfsnDhQgDuvfdeHn74YV760pfyyU9+koceeujRMfzhD38Yq0uWJK1DrqA0kjwbOIfOo8S3dHs8A0bzWHAb7LrrrrzwhS/kzjvv5FOf+hRTpkzhLW95C8cddxwzZsxg4sSJnH766Y9bARmwxx57MHHiRPbcc0+OOeYY3vKWt3DYYYdx1llncdBBB4242gKw6aabcuaZZ/L2t7+d+++/n80224wf/vCHvOlNb6K/v5+9996bqmL77bfnO9/5zrqaAknSGMpIy/Bat3p7e2vg5s8BN9xwA8961rO6NKK1c8wxx/CKV7yCOXPmdHsoa2VDnHNJ2hAlWVRVvaOp6xaPJElqHbd49ISdfvrp3R6CJGkj4wpKC7nttv4415LUTgaUlpkyZQp33323/3CuB1XF3Xff/bjHmCVJ7eAWT8vsuOOO3H777bT17/RsbKZMmfLoF79JktrDgNIykyZN8ptOJUnjnls8kiSpdQwokiSpdQwokiSpdfwm2S5Kci+dPzaox2wH/Kbbg2gh52VlzsnwnJeVOScr69acPKOqth9NRW+S7a6bRvuVv+NFkj7nZGXOy8qck+E5LytzTla2IcyJWzySJKl1DCiSJKl1DCjdNb/bA2gh52R4zsvKnJPhOS8rc05W1vo58SZZSZLUOq6gSJKk1jGgSJKk1jGgdEmSlyW5KclPk8zr9njGQpLTktyV5PpBZdskuSDJLc3vJzflSfKx5vqXJNl70DlHN/VvSXL0oPJZSZY253wsSUbqow2SPD3JhUluSLIsyd825eN2XpJMSXJ1kuuaOfmXpnynJFc14z0zyaZN+eTm9U+b96cNautdTflNSV46qHzYz9eq+miLJBOSXJvknOa1c5L0N//9Xpykrykbt58fgCRbJ/lGkhub/22ZvVHOSVX5s55/gAnArcAzgU2B64Bnd3tcY3BdLwD2Bq4fVPZ/gXnN8Tzgw83xy4H/BAI8D7iqKd8G+Fnz+8nN8ZOb964GZjfn/Cfw5yP10YYfoAfYuzneCrgZePZ4npdmnFs2x5OAq5pr/Trwuqb8U8Cbm+O3AJ9qjl8HnNkcP7v57EwGdmo+UxNG+nytqo+2/AB/D3wFOGek8Y6zOekHthtSNm4/P814vgC8qTneFNh6Y5yTrk/0ePxp/oM/b9DrdwHv6va4xujapvH4gHIT0NMc99D5cjqATwNHDq0HHAl8elD5p5uyHuDGQeWP1ltVH238Ab4LvMR5eXS8mwPXAM+l862WE5vyRz8jwHnA7OZ4YlMvQz83A/VW9flqzhm2jzb8ADsCPwJeBJwz0njHy5w0Y+pn5YAybj8/wJOA22gectmY58Qtnu7YAfjFoNe3N2Ubo6dW1a8Amt9PacpXNQcjld8+TPlIfbRKswy/F50Vg3E9L81WxmLgLuACOv/v/p6qeripMvg6Hr325v3lwLas+VxtO0IfbXAS8E7gkeb1SOMdL3MCUMD5SRYlmduUjefPzzOBXwOfb7YDP5tkCzbCOTGgdEeGKRtvz3uvag7WtHyDkGRL4JvA31XV70aqOkzZRjcvVbWiqmbSWTXYF3jWcNWa32M1J62dqySvAO6qqkWDi4epOm7mZJD9q2pv4M+BtyZ5wQh1N8brH2oina30T1bVXsAf6Gy3rMoGOycGlO64HXj6oNc7Ar/s0ljWtTuT9AA0v+9qylc1ByOV7zhM+Uh9tEKSSXTCyRlV9a2meNzPC0BV3QMsoLM3vnWSgb8PNvg6Hr325v2pwP+w5nP1mxH66Lb9gVcl6Qe+Rmeb5yTG95wAUFW/bH7fBXybTqAdz5+f24Hbq+qq5vU36ASWjW5ODCjdsRCY3tw9vymdm9zO7vKY1pWzgYG7w4+mcw/GQPlfN3eYPw9Y3iwZngcckuTJzR3ih9DZE/8VcG+S5zV3lP/1kLaG66PrmrF+Drihqj466K1xOy9Jtk+ydXO8GXAwcANwITCnqTZ0TgauYw7w4+psgp8NvC6dJ1p2AqbTublv2M9Xc86q+uiqqnpXVe1YVdPojPfHVXUU43hOAJJskWSrgWM6/72/nnH8+amq/wZ+kWTXpujFwE/YGOek2zf8jNcfOndW30xn7/2fuj2eMbqmrwK/Ah6ik8L/hs4e94+AW5rf2zR1A5zaXP9SoHdQO28Eftr8HDuovJfO/zjdCnycx74Jedg+2vADPJ/O8ugSYHHz8/LxPC/AHsC1zZxcD7y3KX8mnX9MfwqcBUxuyqc0r3/avP/MQW39U3PdN9E8adCUD/v5WlUfbfoBDuSxp3jG9Zw0Y7uu+Vk2MO7x/PlpxjYT6Gs+Q9+h8xTORjcnftW9JElqHbd4JElS6xhQJElS6xhQJElS6xhQJElS6xhQJElS6xhQJElS6xhQJElS6/w/bQirQtD+YUQAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAHVCAYAAAA5NRumAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4zLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvnQurowAAIABJREFUeJzt3XucXVV9///XmwQTbga5aCNYgzSCSrhlQAFR8IJWvJYoKrWA9csXb7S/fqlNrVWqtcVv/Qp4NyriBRXximAFVML9kgmEBModxgpaVJSIchHC5/fH2YPDODOZSWbmzM68no/HPGafddZea+3FSfJmrX3OSVUhSZLUFht1ewCSJEljYXiRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtMrPbA5jOttlmm5o3b163hyFJ0qRYvnz5L6tq2/Vtx/DSRfPmzaO3t7fbw5AkaVIk+fF4tOO2kSRJahXDiyRJahXDiyRJahXDiyRJahXDiyRJahXDiyRJahXfKt1Fq+5YzbzFZ3V7GJIk/ZG+4w/u9hCG5cqLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqlSkTXpLMS3LNGOq/PMni5vi4JMeO1GaSniQfHr8RQ5L3JnnBEOUHJDlzPPuSJEkdrX2rdFWdAZwxhvq9wLh+hXNVvXs825MkSWs3ZVZeGjOTfD7JyiRfT7Jpkr4k28AjqydLm+Mjknx0cANJFia5OsmlwFsHlD+yGtKs1JycZGmSW5McM6DePye5Psm5Sb4y1IrOgLqnJFnUHL+4Oe8i4C/GaT4kSdIgUy287AQsqapdgd8Ab1mHNj4HHFNV+6yl3s7Ai4C9gfck2ThJD3AIsAedANIzmg6TzAY+DbwM2B/4kxHqHpWkN0nvmntXj6Z5SZI0wFQLLz+pqoub4y8Bzx7LyUnmAFtW1flN0RdHqH5WVT1QVb8Efg48oenvO1V1X1XdA3x3lF3vDNxWVTdVVTVjH1JVLamqnqrqmbHpnFE2L0mS+k218FJDPH6IP4xz9lrOzxBtDOeBAcdr6Nz/k1GeO5TR9itJktbDVAsvf5qkf7vndcBFQB+wsCk7ZKSTq+puYHWS/hWbw8bY/0XAy5LMTrI5MNovdrge2CHJjs3j142xX0mSNEpTLbxcBxyeZCWwFfAJ4F+Ak5JcSGeFZG2OBD7W3LB731g6r6pldN7BdDXwTTrvTlrrjSlVdT9wFHBWc8Puj8fSryRJGr10btFQvySbV9Vvk2wKXAAcVVVXTkRfs+bOr7mHnzgRTUuStF4m4lulkyyvqlG9GWYkrf2clwm0JMnT6dxf8/mJCi6SJGndGF4GqarXDy5L8jFgv0HFJ1XV5yZnVJIkqZ/bRl3U09NTvb3j+qG/kiRNWeO1bTTVbtiVJEkakeFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1iuFFkiS1ysxuD2A6W3XHauYtPqvbw5Aepe/4g7s9BEkakSsvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVdY5vCT5XpItx1D/neva12RKckSSJ3Z7HJIkaWjrHF6q6iVVdffAsnQM1+aUCS9JZozw9BHAmMJLEt9yLknSJBlVeEny7STLk1yb5KimrC/JNknmJbkuyceBK4EnDXH+8cAmSVYkObUp+8skVzRln+oPFEl+m+QDTX8/SLJ3kqVJbk3y8qbOEUm+k+T7SW5I8p4BfY3U7nuTXA7sk+TdSZYluSbJkiZ4LQJ6gFOb8zfpv86mjZ4kS5vj45rzzgG+kGRGkv9o2lyZ5H8PM5dHJelN0rvm3tWjmX5JkjTAaFde3lhVC+n8w35Mkq0HPb8T8IWq2qOqfjz45KpaDNxXVbtX1WFJngYcCuxXVbsDa4DDmuqbAUub/u4B/hV4IfAq4L0Dmt27OWd34NVNsFhbu9dU1TOr6iLgo1W1V1XtAmwCvLSqvg70Aoc1Y71vLfOyEHhFVb0e+GtgdVXtBewF/K8kOwwxF0uqqqeqemZsOmctzUuSpMFGu91xTJJXNcdPAuYPev7HVXXZGPp9Pp1/+JclgU54+Hnz3O+B7zfHq4AHqurBJKuAeQPaOLeq7gJI8k3g2cBDI7S7BvjGgPMPTPIOYFNgK+Ba4LtjuAaAMwYEnIOAXZvVG4A5dObptjG2KUmSRrDW8JLkAOAFwD5VdW+zbTJ7ULXfjbHfAJ+vqn8c4rkHq6qa44eBBwCq6uFB95bUoPNqLe3eX1VrAJLMBj4O9FTVT5Icxx9fU7+H+MMK1UjXHeDtVXX2MO1IkqRxMJptoznAr5vgsjPwrHXs68EkGzfHPwQWJXk8QJKtkjx5jO29sDlvE+CVwMVjaLc/hPwyyebAogHP3QNsMeBxH53VHIBDRhjP2cCb+68xyVOTbDbGa5IkSWsxmm2j7wNHJ1kJ3ACMZXtooCXAyiRXNve9vAs4p3l30oPAW4E/ul9mBBcBXwT+DPhyVfUCjKbdqro7yafpbEv1AcsGPH0K8Mkk9wH7AP8CfLZ5q/flI4znM3S2ta5MZ8/qF3RClSRJGkf5ww5NeyQ5gs6Wz9u6PZb1MWvu/Jp7+IndHob0KH6rtKSJkmR5VfWsbzt+wq4kSWqVcV95aT5HZdag4jdU1apx7WgD0NPTU729vd0ehiRJk2K8Vl7G/ZNhq+qZ492mJElSP7eNJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSqxheJElSq8zs9gCms1V3rGbe4rO6PQyto77jD+72ECRpWnLlRZIktYrhRZIktYrhRZIktYrhRZIktYrhZZSS7J7kJd0ehyRJ053hZfR2B8YUXpL4bi5JksbZtAovSeYluT7JZ5Jck+TUJC9IcnGSm5LsnWSzJCcnWZbkqiSvSPIY4L3AoUlWJDl0qHpNH0ckOT3Jd4FzunrBkiRtgKbjysCfAa8GjgKWAa8Hng28HHgn8F/Aj6rqjUm2BK4AfgC8G+ipqrcBJPm3wfWS/KDpYx9g16r61SRelyRJ08J0DC+3VdUqgCTXAj+sqkqyCpgHbA+8PMmxTf3ZwJ8O0c5BI9Q7d7jgkuQoOsGJGY/ddhwuR5Kk6WU6hpcHBhw/PODxw3TmYw1wSFXdMPCkJM8c1E5GqPe74TqvqiXAEoBZc+fXulyAJEnT2bS652WUzgbeniQASfZoyu8BthhFPUmSNIEML3/sfcDGwMok1zSPAc4Dnt5/w+4I9SRJ0gRKlTsX3TJr7vyae/iJ3R6G1pFfzChJY5NkeVX1rG87rrxIkqRWMbxIkqRWMbxIkqRWmY5vlZ4yFmw3h17vm5AkaUxceZEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa1ieJEkSa0ys9sDmM5W3bGaeYvP6vYwWqPv+IO7PQRJ0hTgyoskSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw8sQkuycZEWSq5LsuA7n/22STSdibJIkTXfTIrwkmTHGU14JfKeq9qiqW9ahy78FDC+SJE2AroWXJJslOSvJ1UmuSXJokhcnuT7JRUk+nOTMpu5xSY4dcO41SeY1x99OsjzJtUmOGlDnt0nem+RyYJ8kC5Oc39Q9O8ncYcb1Ejrh401JzmvK/jLJFc1qzKf6w1CSg5JcmuTKJKcn2TzJMcATgfP6zx/U/lFJepP0rrl39TjNpiRJ00c3V15eDPy0qnarql2A7wOfBl4G7A/8ySjbeWNVLQR6gGOSbN2UbwZcU1XPBC4HPgIsauqeDLx/qMaq6nvAJ4ETqurAJE8DDgX2q6rdgTXAYUm2Ad4FvKCq9gR6gb+rqg8DPwUOrKoDh2h/SVX1VFXPjE3njPISJUlSv25+wu4q4INJPgCcCdwD3FZVNwEk+RJw1Ajn9zsmyaua4ycB84G76ISMbzTlOwG7AOcmAZgB/GyU43w+sBBY1py7CfBz4FnA04GLm/LHAJeOsk1JkrSOuhZequrGJAuBlwD/DpwD1DDVH+LRq0SzAZIcALwA2Keq7k2ytP854P6qWtMcB7i2qvZZh6EG+HxV/eOjCpOXAedW1evWoU1JkrSOunnPyxOBe6vqS8AHgX2BHQa8u2dgKOgD9mzO2xPYoSmfA/y6CS4701kNGcoNwLZJ9mna2DjJM0Y51B8Ci5I8vjl3qyRPBi4D9kvyZ035pkme2pxzD7DFKNuXJElj0M1towXAfyR5GHgQeDOwDXBWkl8CF9HZ6oHO9s9fJVkBLANubMq/DxydZCWdgHLZUB1V1e+TLAI+nGQOnes+Ebh2bYOsqv9K8i7gnCQbNWN9a1VdluQI4CtJZjXV39WMbQnwn0l+NtR9L5Ikad2laridmu5qtoSOraqXdnssE2XW3Pk19/ATuz2M1vBbpSWp3ZIsr6qe9W1nWnzOiyRJ2nBM2ZWXyZDkY8B+g4pPqqrPTUb/PT091dvbOxldSZLUdeO18tLNe166rqre2u0xSJKksXHbSJIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktYrhRZIktcrMbg9gOlt1x2rmLT6r28OYsvqOP7jbQ5AkTUGuvEiSpFYxvEiSpFYxvEiSpFYxvEiSpFYxvDSSzEtyTbfHIUmSRmZ4mUBJfDeXJEnjzPDyaDOSfDrJtUnOSbJJkt2TXJZkZZJvJXkcQJKlSXqa422S9DXHRyQ5Pcl3gXO6dymSJG2YDC+PNh/4WFU9A7gbOAT4AvAPVbUrsAp4zyja2Qc4vKqeN/iJJEcl6U3Su+be1eM4dEmSpgfDy6PdVlUrmuPlwI7AllV1flP2eeA5o2jn3Kr61VBPVNWSquqpqp4Zm85Z/xFLkjTNGF4e7YEBx2uALUeo+xB/mL/Zg5773XgOSpIk/YHhZWSrgV8n2b95/AagfxWmD1jYHC+a5HFJkjRt+W6YtTsc+GSSTYFbgSOb8g8CX0vyBuBH3RqcJEnTTaqq22OYtmbNnV9zDz+x28OYsvxiRknasCRZXlU969uO20aSJKlVDC+SJKlVDC+SJKlVvGG3ixZsN4de7+uQJGlMXHmRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtMrPbA5jOVt2xmnmLz+r2MCZE3/EHd3sIkqQNlCsvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVVoZXpLMS3LNBLT73iQvGKL8gCRnNscvT7K4OX5lkqeP9zgkSdLwfLfRAFX17lHUOQM4o3n4SuBM4L8mclySJOkPWrny0piR5NNJrk1yTpJNkixN0gOQZJskfc3xEUm+neS7SW5L8rYkf5fkqiSXJdmqqXdKkkXN8YuTXJ/kIuAv+jtt2vpokn2BlwP/kWRFkh2TXDmg3vwkyydvOiRJmh7aHF7mAx+rqmcAdwOHrKX+LsDrgb2B9wP3VtUewKXAXw2smGQ28GngZcD+wJ8MbqyqLqGzAvP3VbV7Vd0CrE6ye1PlSOCUweclOSpJb5LeNfeuHu21SpKkRpvDy21VtaI5Xg7MW0v986rqnqr6BbAa+G5TvmqIc3du2r+pqgr40ijH9BngyCQzgEOBLw+uUFVLqqqnqnpmbDpnlM1KkqR+bQ4vDww4XkPn/p2H+MM1zR6h/sMDHj/M0Pf+1DqM6RvAnwMvBZZX1V3r0IYkSRpBm8PLUPqAhc3xovVo53pghyQ7No9fN0y9e4At+h9U1f3A2cAngM+tR/+SJGkYG1p4+SDw5iSXANusayNNCDkKOKu5YffHw1T9KvD3zY2//UHnVDqrNuesa/+SJGl46dzSofGS5FhgTlX989rqzpo7v+YefuIkjGry+cWMkqTBkiyvqp71bcfPeRlHSb4F7Ag8r9tjkSRpQ2V4GUdV9apuj0GSpA2d4aWLFmw3h163VyRJGpMN7YZdSZK0gTO8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVpnZ7QFMZ6vuWM28xWdNeD99xx884X1IkjRZXHmRJEmtYniRJEmtYniRJEmtYniRJEmtYniRJEmtsl7hJclxSY4dj4EkWZqkZzzakiRJGy5XXgZJ4tvHJUmawsYUXpL8VZKVSa5O8sVBz+2e5LLm+W8leVxT/siKSpJtkvQ1x5sk+WpT/zRgk7X0/eIkVzZ9/7Ap2yrJt5s2Lkuya5KNkvQl2XLAuTcneUKSbZN8I8my5me/5vnjkixJcg7whSTzklzY9Hdlkn2behsl+XiSa5OcmeR7SRY1zy1Mcn6S5UnOTjJ3mOs4Kklvkt41964ey/RLkiTGEF6SPAP4J+B5VbUb8DeDqnwB+Ieq2hVYBbxnLU2+Gbi3qf9+YOEIfW8LfBo4pOn71c1T/wJc1bTxTuALVfUw8B3gVc25zwT6qupO4CTghKraCzgE+MyAbhYCr6iq1wM/B15YVXsChwIfbur8BTAPWAC8Cdin6WNj4CPAoqpaCJzcXNMfqaolVdVTVT0zNp2zlimSJEmDjWWL5HnA16vqlwBV9askACSZA2xZVec3dT8PnL6W9p5DEwqqamWSlSPUfRZwQVXd1t93U/5sOiGEqvpRkq2bsZwGvBv4HPDa5jHAC4Cn948beGySLZrjM6rqvuZ4Y+CjSXYH1gBPHdDf6U1A+p8k5zXlOwG7AOc2bc8AfraW65ckSetgLOElQK1DHw/xhxWe2YOeG217w/WdIcoKuBT4s2bF5pXAvzbPbQTsMyCkdBrpBI7fDSj6/4A7gd2ac+4fob/+8murap+1XokkSVovY7nn5YfAa5JsDZ37TfqfqKrVwK+T7N8UvQHoX4Xp4w9bQosGtHcBcFjT1i7AriP0fSnw3CQ7DOp7YBsHAL+sqt9UVQHfAj4EXFdVdzX1zwHe1t9os7IylDnAz5oVljfQWUkBuAg4pLn35QnAAU35DcC2SR7ZRmq22SRJ0jgb9cpLVV2b5P3A+UnWAFfRCSb9Dgc+mWRT4FbgyKb8g8DXkrwB+NGA+p8APtdsF60Arhih718kOQr4ZpKNaO5JAY4b0Ma9zRj6nQYsA44YUHYM8LGm/kw64efoIbr8OPCNJK8GzuMPqzLfAJ4PXAPcCFwOrK6q3zc37n642baaCZwIXDvcNUmSpHWTziKFRivJ5lX122YF6gpgv6r6n3Vpa9bc+TX38BPHd4BD8FulJUlTQZLlVbXen+nmZ5qM3ZnN27AfA7xvXYOLJElaN1MuvCS5HJg1qPgNVbWqG+MZrKoO6PYYJEmazqZceKmqZ3Z7DJNlwXZz6HVLR5KkMfHrASRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqsYXiRJUqvM7PYAprNVd6xm3uKzJryfvuMPnvA+JEmaLK68SJKkVjG8SJKkVjG8SJKkVjG8SJKkVjG8SJKkVpmW4SXJZ5I8vTl+53q0870kWw5RflySY9dnjJIkaWjTMrxU1Zuq6r+ah+scXqrqJVV19zgNS5IkjcIGH16SbJbkrCRXJ7kmyaFJlibpSXI8sEmSFUlOber/ZZIrmrJPJZkxQtt9SbZpjv8pyQ1JfgDsNDlXJ0nS9LPBhxfgxcBPq2q3qtoF+H7/E1W1GLivqnavqsOSPA04FNivqnYH1gCHra2DJAuB1wJ7AH8B7DVC3aOS9CbpXXPv6vW6MEmSpqPpEF5WAS9I8oEk+1fVSInh+cBCYFmSFc3jp4yij/2Bb1XVvVX1G+CM4SpW1ZKq6qmqnhmbzhnDZUiSJJgGXw9QVTc2KyMvAf49yTkjVA/w+ar6x3Xpap0GKEmSxmSDX3lJ8kTg3qr6EvBBYM9BVR5MsnFz/ENgUZLHN+duleTJo+jmAuBVSTZJsgXwsnEaviRJGmSDX3kBFgD/keRh4EHgzXRCTL8lwMokVzb3vbwLOCfJRk39twI/HqmDqroyyWnAiqbuhRNwHZIkCUiVux3dMmvu/Jp7+IkT3o/fKi1JmgqSLK+qnvVtZ4PfNpIkSRuW6bBttN6SXA7MGlT8hqpa1Y3xSJI0nRleRqGqnjkR7S7Ybg69bulIkjQmbhtJkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWmdntAUxnq+5YzbzFZ417u33HHzzubUqSNFW48iJJklrF8CJJklrF8CJJklrF8CJJklrF8CJJklplyoaXJN9LsuUY6s9Lcs1Ejmkskvy222OQJGlDNGXfKl1VL+n2GCRJ0tTTtZWXJO9IckxzfEKSHzXHz0/ypSR9SbZpVlSuS/LpJNcmOSfJJk3dhUmuTnIp8Na19PeMJFckWZFkZZL5TdvXJ/l8U/b1JJsOaPv8JMuTnJ1kblO+Y5LvN+UXJtm5Kd8hyaVJliV53wjjOCpJb5LeNfeuHpe5lCRpOunmttEFwP7NcQ+weZKNgWcDFw6qOx/4WFU9A7gbOKQp/xxwTFXtM4r+jgZOqqrdm/5ub8p3ApZU1a7Ab4C3NOP4CLCoqhYCJwPvb+ovAd7elB8LfLwpPwn4RFXtBfzPcIOoqiVV1VNVPTM2nTOKYUuSpIG6GV6WAwuTbAE8AFxKJ1Tszx+Hl9uqasWA8+YlmQNsWVXnN+VfXEt/lwLvTPIPwJOr6r6m/CdVdXFz/CU64WknYBfg3CQrgHcB2yfZHNgXOL0p/xQwtzl3P+AroxyLJElaR12756WqHkzSBxwJXAKsBA4EdgSuG1T9gQHHa4BNgAA1hv6+nORy4GDg7CRvAm4doo1q2r528IpOkscCdzerN0N2M9rxSJKkddPtdxtdQGfr5QI6qy1HAyuqaq0hoKruBlYneXZTdNhI9ZM8Bbi1qj4MnAHs2jz1p0n6Q8rrgIuAG4Bt+8uTbJzkGVX1G+C2JK9uypNkt+bci4HXjmYskiRp3XU7vFxIZ9vl0qq6E7ifP94yGsmRwMeaG3bvW0vdQ4Frmu2enYEvNOXXAYcnWQlsRee+ld8Di4APJLkaWEFnuwg6weSvm/JrgVc05X8DvDXJMsCbWSRJmiAZxSLHBivJPODMqtqlG/3Pmju/5h5+4ri367dKS5KmoiTLq6pnfdvp9sqLJEnSmEzZD6lbV0leBHxgUPFtVfWqwXWrqo/Ou4okSVJLTOtto27r6emp3t7ebg9DkqRJ4baRJEmalgwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVQwvkiSpVWZ2ewDT2ao7VjNv8Vnj1l7f8QePW1uSJE1VrrxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWMbxIkqRWmZLhJcl7k7yg2+OQJElTz5R7q3SSGVX17gloN0Cq6uHxbluSJE2eSV15STIvyfVJPp9kZZKvJ9k0SV+Sdye5CHh1klOSLGrO6Uvyb0kuTdKbZM8kZye5JcnRTZ3Nk/wwyZVJViV5xYD+rkvyceBK4J+TnDBgPP8ryYdGGOs1Ax4fm+S45viYJP/VXMNXm7LNkpycZFmSq/rHIEmSxlc3Vl52Av66qi5OcjLwlqb8/qp6NkCSFw865ydVtU8TPE4B9gNmA9cCnwTuB15VVb9Jsg1wWZIzBvR3ZFW9JclmwMok76iqB4Ejgf+9DtewGNihqh5IsmVT9k/Aj6rqjU3ZFUl+UFW/G3hikqOAowBmPHbbdehakqTprRv3vPykqi5ujr8EPLs5Pm2Ec/qDyCrg8qq6p6p+AdzfBIUA/5ZkJfADYDvgCc05P66qywCaIPEj4KVJdgY2rqpV63ANK4FTk/wl8FBTdhCwOMkKYCmdcPWng0+sqiVV1VNVPTM2nbMOXUuSNL11Y+Wlhnn8u8EVB3ig+f3wgOP+xzOBw4BtgYVV9WCSPjrhYah2PwO8E7ge+NwIfT7Eo8Pd7AHHBwPPAV5OZyvqGXQC1CFVdcMIbUqSpPXUjZWXP02yT3P8OuCicWhzDvDzJrgcCDx5uIpVdTnwJOD1wFdGaPNO4PFJtk4yC3gpQJKNgCdV1XnAO4Atgc2Bs4G3NzcGk2SP9b8sSZI0WDfCy3XA4c0Wz1bAJ8ahzVOBniS9dFZhrl9L/a8BF1fVr4er0NwT817gcuDMAW3OAL6UZBVwFXBCVd0NvA/YmM49Ndc0jyVJ0jhL1eBdnAnsLJkHnFlVu0xap0OP40w6oeOH3RzHrLnza+7hJ45be36rtCRpKkuyvKp61redKfkhdRMlyZZJbgTu63ZwkSRJ62ZSb9itqj6ga6suzfbOUweWJdkaGCrIPL+q7pqUgUmSpFGbcp+wO9magLJ7N/pesN0cet3qkSRpTKbVtpEkSWo/w4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWqVmd0ewHS26o7VzFt81jqf33f8weM4GkmS2sGVF0mS1CqGF0mS1CqGF0k2ta1dAAARsElEQVSS1CqGF0mS1CrTPrwkWZqkZy11jkjy0ckakyRJGt60Dy+SJKldWhdekrwjyTHN8QlJftQcPz/Jl5IclOTSJFcmOT3J5s3zC5Ocn2R5krOTzB3U7kZJPp/kX5vHRya5Mcn5wH4D6r0syeVJrkrygyRPaM69Kcm2A9q6Ock2kzQtkiRNG60LL8AFwP7NcQ+weZKNgWcDq4B3AS+oqj2BXuDvmuc/AiyqqoXAycD7B7Q5EzgVuLGq3tUEm3+hE1peCDx9QN2LgGdV1R7AV4F3VNXDwJeAw5o6LwCurqpfDh58kqOS9CbpXXPv6vWdC0mSpp02fkjdcmBhki2AB4Ar6YSY/YEz6ASNi5MAPAa4FNgJ2AU4tymfAfxsQJufAr5WVf2B5pnA0qr6BUCS04CnNs9tD5zWBJzHALc15ScD3wFOBN4IfG6owVfVEmAJwKy582tdJ0GSpOmqdeGlqh5M0gccCVwCrAQOBHakEyTOrarXDTwnyQLg2qraZ5hmLwEOTPL/qur+/q6GqfsR4ENVdUaSA4DjmnH9JMmdSZ5HJ/wcNsz5kiRpPbRx2wg6W0fHNr8vBI4GVgCXAfsl+TOAJJsmeSpwA7Btkn2a8o2TPGNAe58FvgecnmQmcDlwQJKtmy2nVw+oOwe4ozk+fNC4PkNn++hrVbVm3K5WkiQ9oq3h5UJgLnBpVd0J3A9c2GzzHAF8JclKOmFm56r6PbAI+ECSq+kEnX0HNlhVH6KzBfVF4E46KyqXAj9oyvsdRyfkXAgMvqflDGBzhtkykiRJ6y9V3nYxXprPizmhqvZfa2U697zMPfzEde7PL2aUJLVJkuVVNeJnq41G6+55maqSLAbejPe6SJI0odq6bTTlVNXxVfXkqrqo22ORJGlDZniRJEmt4rZRFy3Ybg693rciSdKYuPIiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaxfAiSZJaZWa3BzCdrbpjNfMWnzXm8/qOP3gCRiNJUju48iJJklrF8CJJklrF8CJJklrF8CJJklrF8CJJklqlq+ElySXd7L8ZwylJFk1Au0uT9Ix3u5IkTXddDS9VtW83+5ckSe3T7ZWX3za/5ya5IMmKJNck2X+Y+q9J8qHm+G+S3Noc75jkouZ4YZLzkyxPcnaSuQPqfL8pvzDJzkO0/75mJWajEdpZmuQDSa5IcmP/WJNskuSrSVYmOQ3YZAKmTJKkaW+q3PPyeuDsqtod2A1YMUy9C4D+YLM/cFeS7YBnAxcm2Rj4CLCoqhYCJwPvb+ovAd7elB8LfHxgw0n+L/B44EhgxgjtAMysqr2BvwXe05S9Gbi3qnZt6i4c6gKSHJWkN0nvmntXr2VaJEnSYFPlE3aXASc34ePbVTVkeKmq/0myeZItgCcBXwaeQyfIfBPYCdgFODcJdELIz5JsDuwLnN6UA8wa0PQ/A5dX1VEASYZsZ0D9bza/lwPzmuPnAB9uxrkyycphrmEJnSDFrLnza8RZkSRJf2RKhJequiDJc4CDgS8m+Y+q+sIw1S+lszpyA3Ah8EZgH+D/AH8KXFtV+ww8IcljgbublZ2hLAMWJtmqqn4FZKh2Bnig+b2GR8+hYUSSpAk2JbaNkjwZ+HlVfRr4LLDnCNUvoLPtcwFwFXAg8EBVraYTaLZNsk/T7sZJnlFVvwFuS/LqpjxJdhvQ5veB44GzmlWdIdtZy2VcABzW1N8F2HX0MyBJkkZrSoQX4ABgRZKrgEOAk0aoeyGdLaMLqmoN8BPgIoCq+j2wCPhAkqvp3DvT/46mw4C/bsqvBV4xsNGqOh34NHAGnW2i4doZzieAzZvtoncAV6z9siVJ0lilyp2Obpk1d37NPfzEMZ/nt0pLktooyfKqWu/PQJsqKy+SJEmjMiVu2B1Kkst59DuCAN5QVau6MR5JkjQ1TNnwUlXP7PYYJtqC7ebQ6xaQJElj4raRJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqFcOLJElqlZndHoAe7cEHH+T222/n/vvv7/ZQpoXZs2ez/fbbs/HGG3d7KJKkUTK8TDG33347W2yxBfPmzSNJt4ezQasq7rrrLm6//XZ22GGHbg9HkjRKbhtNMffffz9bb721wWUSJGHrrbd2lUuSWsbwMgUZXCaPcy1J7WN4mUBJZnR7DJIkbWi856WR5H3AL6vqpObx+4E7gVnAa5rf36qq9zTPfxt4EjAbOKmqljTlvwU+BLwI+D/AReszrnmLz1qf0/9I3/EHr7XOvvvuyyWXXDKu/Y6kr6+PSy65hNe//vWT1qckqb1cefmDzwKHAyTZCHgtnfAyH9gb2B1YmOQ5Tf03VtVCoAc4JsnWTflmwDVV9cyq+qPgkuSoJL1Jen/xi19M7BWto8kMLg899BB9fX18+ctfnrQ+JUntZnhpVFUfcFeSPYCDgKuAvQYcXwnsTCfMQCewXA1cRmcFpr98DfCNEfpZUlU9VdWz7bbbTsSlrLfNN98cgKVLl/Lc5z6X17zmNTz1qU9l8eLFnHrqqey9994sWLCAW265BYAjjjiCo48+mv3335+nPvWpnHnmmUDn5uMjjzySBQsWsMcee3DeeecBcMopp/DqV7+al73sZRx00EEsXryYCy+8kN13350TTjiBvr4+9t9/f/bcc0/23HPPR8LU0qVLOeCAA1i0aBE777wzhx12GFUFwLJly9h3333Zbbfd2HvvvbnnnntYs2YNf//3f89ee+3Frrvuyqc+9anJnkpJ0gRw2+jRPgMcAfwJcDLwfODfq+pR/+olOQB4AbBPVd2bZCmd7SOA+6tqzWQNeKJdffXVXHfddWy11VY85SlP4U1vehNXXHEFJ510Eh/5yEc48cQTgc7Wz/nnn88tt9zCgQceyM0338zHPvYxAFatWsX111/PQQcdxI033gjApZdeysqVK9lqq61YunQpH/zgBx8JPffeey/nnnsus2fP5qabbuJ1r3sdvb29AFx11VVce+21PPGJT2S//fbj4osvZu+99+bQQw/ltNNOY6+99uI3v/kNm2yyCZ/97GeZM2cOy5Yt44EHHmC//fbjoIMO8m3RktRyhpdH+xbwXmBj4PXAQ8D7kpxaVb9Nsh3wIDAH+HUTXHYGntW1EU+wvfbai7lz5wKw4447ctBBBwGwYMGCR1ZSAF7zmtew0UYbMX/+fJ7ylKdw/fXXc9FFF/H2t78dgJ133pknP/nJj4SXF77whWy11VZD9vnggw/ytre9jRUrVjBjxoxHzgHYe++92X777QHYfffd6evrY86cOcydO5e99toLgMc+9rEAnHPOOaxcuZKvf/3rAKxevZqbbrrJ8CJJLWd4GaCqfp/kPODuZvXknCRPAy5t3lL7W+Avge8DRydZCdxAZ+togzRr1qxHjjfaaKNHHm+00UY89NBDjzw3+C3HSR7Z0hnKZpttNuxzJ5xwAk94whO4+uqrefjhh5k9e/Yjzw0cz4wZM3jooYeoqiHf8lxVfOQjH+FFL3rRCFcoSWob73kZoLlR91l0bt4FoKpOqqoFzc8+VXVLVT1QVX9eVbtW1aur6oCqWtrU37xLw++q008/nYcffphbbrmFW2+9lZ122onnPOc5nHrqqQDceOON/Pd//zc77bTTH527xRZbcM899zzyePXq1cydO5eNNtqIL37xi6xZM/Iu3M4778xPf/pTli1bBsA999zDQw89xIte9CI+8YlP8OCDDz4yht/97nfjdcmSpC5x5aWR5OnAmXTeDn1Tt8fTbzRvbZ4KdtppJ5773Ody55138slPfpLZs2fzlre8haOPPpoFCxYwc+ZMTjnllEetnPTbddddmTlzJrvtthtHHHEEb3nLWzjkkEM4/fTTOfDAA0dcpQF4zGMew2mnncbb3/527rvvPjbZZBN+8IMf8KY3vYm+vj723HNPqoptt92Wb3/72xM1BZKkSZKRlvY1sXp6eqr/RtR+1113HU972tO6NKJ1c8QRR/DSl76URYsWdXso66SNcy5JbZRkeVX1rG87bhtJkqRWcdtI6+2UU07p9hAkSdOIKy9TkFt5k8e5lqT2MbxMMbNnz+auu+7yH9VJUFXcddddj3ortiRp6nPbaIrZfvvtuf3225mq33u0oZk9e/YjH3onSWoHw8sUs/HGG/sJsJIkjcBtI0mS1CqGF0mS1CqGF0mS1Cp+wm4XJbmHzhc7amJsA/yy24PYQDm3E8v5nVjO78RZ29w+uaq2Xd9OvGG3u24Yj49J1tCS9Dq/E8O5nVjO78RyfifOZM2t20aSJKlVDC+SJKlVDC/dtaTbA9jAOb8Tx7mdWM7vxHJ+J86kzK037EqSpFZx5UWSJLWK4UWSJLWK4aVLkrw4yQ1Jbk6yuNvjmcqS9CVZlWRFkt6mbKsk5ya5qfn9uKY8ST7czOvKJHsOaOfwpv5NSQ4fUL6waf/m5txM/lVOniQnJ/l5kmsGlE34fA7Xx4ZkmLk9Lskdzet3RZKXDHjuH5t5uiHJiwaUD/n3Q5IdklzezOFpSR7TlM9qHt/cPD9vcq54ciV5UpLzklyX5Nokf9OU+/pdTyPM7dR8/VaVP5P8A8wAbgGeAjwGuBp4erfHNVV/gD5gm0Fl/xdY3BwvBj7QHL8E+E8gwLOAy5vyrYBbm9+Pa44f1zx3BbBPc85/An/e7Wue4Pl8DrAncM1kzudwfWxIP8PM7XHAsUPUfXrzZ38WsEPzd8KMkf5+AL4GvLY5/iTw5ub4LcAnm+PXAqd1ey4maH7nAns2x1sANzbz6Ot34uZ2Sr5+XXnpjr2Bm6vq1qr6PfBV4BVdHlPbvAL4fHP8eeCVA8q/UB2XAVsmmQu8CDi3qn5VVb8GzgVe3Dz32Kq6tDp/cr4woK0NUlVdAPxqUPFkzOdwfWwwhpnb4bwC+GpVPVBVtwE30/m7Yci/H5oVgOcBX2/OH/zfqX9uvw48f0NcQayqn1XVlc3xPcB1wHb4+l1vI8ztcLr6+jW8dMd2wE8GPL6dkV8k010B5yRZnuSopuwJVfUz6PyhAx7flA83tyOV3z5E+XQzGfM5XB/TwduabYuTB2w3jHVutwburqqHBpU/qq3m+dVN/Q1Ws7WwB3A5vn7H1aC5hSn4+jW8dMdQidL3rA9vv6raE/hz4K1JnjNC3eHmdqzl6nA+198ngB2B3YGfAf+vKR/PuZ1W855kc+AbwN9W1W9GqjpEma/fEQwxt1Py9Wt46Y7bgScNeLw98NMujWXKq6qfNr9/DnyLzrLknc0SL83vnzfVh5vbkcq3H6J8upmM+Ryujw1aVd1ZVWuq6mHg03RevzD2uf0lnW2PmYPKH9VW8/wcRr991SpJNqbzj+upVfXNptjX7zgYam6n6uvX8NIdy4D5zZ3Xj6Fzg9IZXR7TlJRksyRb9B8DBwHX0Jmv/ncIHA58pzk+A/ir5l0GzwJWN0u8ZwMHJXlcs+x5EHB289w9SZ7V7LH+1YC2ppPJmM/h+tig9f+D13gVndcvdObjtc07LXYA5tO5WXTIvx+aezDOAxY15w/+79Q/t4uAHzX1NyjNa+qzwHVV9aEBT/n6XU/Dze2Uff2O9x3L/oz6zu6X0Lmb+xbgn7o9nqn6Q+eO9aubn2v754rOfugPgZua31s15QE+1szrKqBnQFtvpHNT2c3AkQPKe5o/kLcAH6X55OkN9Qf4Cp3l3wfp/B/PX0/GfA7Xx4b0M8zcfrGZu5XNX9JzB9T/p2aebmDAu9yG+/uh+fNwRTPnpwOzmvLZzeObm+ef0u25mKD5fTad7YSVwIrm5yW+fid0bqfk69evB5AkSa3itpEkSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWoVw4skSWqV/x+0msILP+AGSQAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 576x576 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "%%time\n",
    "if DO_TRAIN == True:\n",
    "    modelos = kfolds(X_train, y_train, kf)\n",
    "    #modelos = Stratify_kfolds(X_train, y_train, skf)\n",
    "    #modelos = validate_with_leak(X_train, y_train)\n",
    "    #modelos =  leak_kfolds(X_train, y_train, kf)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "if DO_TRAIN == True:\n",
    "    del X_train, y_train\n",
    "    gc.collect()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "271d15be-39eb-480e-beee-55a907f48f35",
    "_uuid": "3344c390-9730-4009-8372-7268f05d8128"
   },
   "source": [
    "## Preparing test data\n",
    "# # Preparing test data with same features as train data."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {
    "_cell_guid": "82c93e89-4adb-41c7-9e96-130a71588863",
    "_uuid": "76bdab25-7a4a-4a01-b465-ff60530f0a54"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 11.7 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "if DO_TEST == True:\n",
    "    df_test = pd.read_csv(path_test)\n",
    "    df_test['timestamp'] = pd.to_datetime(df_test['timestamp'])\n",
    "    #print('Calculando semana del mes...')\n",
    "    #df_test['weekofmonth'] = df_test['timestamp'].apply(get_week_of_month) #++\n",
    "    #df_test.loc[df_test.weekofmonth < 0, 'weekofmonth'] = 52 + df_test.loc[df_test.weekofmonth < 0, 'weekofmonth']\n",
    "\n",
    "    gc.collect()\n",
    "\n",
    "    #get_summer(df_test)\n",
    "    #get_winter(df_test)\n",
    "    is_weekend(df_test) #++\n",
    "\n",
    "    print('Temperatura test')\n",
    "    weather_test = pd.read_csv(path_weather_test)\n",
    "    #rellenar nans \n",
    "    weather_test = fill_weather_dataset(weather_test)\n",
    "    # a datetime\n",
    "    weather_test['timestamp'] = pd.to_datetime(weather_test['timestamp'])\n",
    "    # ajustar horas de temperaturas\n",
    "    weather_test = timestamp_align(weather_test)\n",
    "    \n",
    "    # temp features\n",
    "    humidity(weather_test)\n",
    "    #feels_like_df(weather_test)\n",
    "    \n",
    "    weather_test = create_lag_features(weather_test, 24)\n",
    "    # TODO probar +48\n",
    "\n",
    "    df_test = reduce_mem_usage(df_test)\n",
    "    weather_test = reduce_mem_usage(weather_test)\n",
    "    gc.collect()\n",
    "\n",
    "    print('Preparando datos')\n",
    "    X_test, row_ids = prepare_data(df_test, building, weather_test, test=True)\n",
    "\n",
    "    print(\"Test memory size: \",  X_test.memory_usage().sum() / 1024 / 1024)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {
    "_cell_guid": "572c3c78-b3e0-4dab-8aa5-e1e3bb7d081e",
    "_kg_hide-input": true,
    "_kg_hide-output": true,
    "_uuid": "0aceb6bf-efbd-481a-9bb4-f8be04127165"
   },
   "outputs": [],
   "source": [
    "if DO_TEST == True:\n",
    "    del  df_test, building, weather_test\n",
    "    gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "if DO_TEST == True:\n",
    "    print(X_test.isna().sum())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "97fc2d64-1baa-4e61-8dc8-85d7f1e52d47",
    "_uuid": "e900c191-cda9-4af9-ba41-06590b55da85"
   },
   "source": [
    "## Scoring test data\n",
    "# # Averaging predictions from the two half train data models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {
    "_cell_guid": "de4cb5c9-7c78-4ca4-a650-911f21fb438b",
    "_uuid": "11f55a56-7515-40a0-9816-6872bef0091b"
   },
   "outputs": [],
   "source": [
    "#---------------------------------------------------------------------------------------------------------------\n",
    "def pred_kfolds(X_test, models, batch_size=1000000):\n",
    "    iterations = (X_test.shape[0] + batch_size -1) // batch_size\n",
    "    print('iterations', iterations)\n",
    "\n",
    "    y_test_pred_total = np.zeros(X_test.shape[0])\n",
    "    print(models)\n",
    "    for i, model in enumerate(models):\n",
    "        print(f'predicting {i}-th model')\n",
    "        for k in tqdm(range(iterations)):\n",
    "            y_pred_test = model.predict(X_test[k*batch_size:(k+1)*batch_size], num_iteration=model.best_iteration)\n",
    "            y_test_pred_total[k*batch_size:(k+1)*batch_size] += y_pred_test\n",
    "            \n",
    "        del model\n",
    "        gc.collect()\n",
    "    \n",
    "\n",
    "    y_test_pred_total /= 3 # prediccion con media de modelos\n",
    "    return y_test_pred_total\n",
    "#---------------------------------------------------------------------------------------------------------------\n",
    "# predecir todo el vector\n",
    "def pred_todo(X_test, models):\n",
    "    \n",
    "    y_test_pred_total = np.zeros(X_test.shape[0])\n",
    "    for i, model in enumerate(models):\n",
    "        print(f'predicting {i}-th model')\n",
    "        y_test_pred_total += model.predict(X_test, num_iteration=model.best_iteration)\n",
    "        \n",
    "        del model\n",
    "        gc.collect()\n",
    "        \n",
    "    y_test_pred_total /= 3 # prediccion con media de modelos\n",
    "    return y_test_pred_total\n",
    "#---------------------------------------------------------------------------------------------------------------\n",
    "def pred_kfolds_from_dict(X_test, models, batch_size=1000000):\n",
    "    iterations = (X_test.shape[0] + batch_size -1) // batch_size\n",
    "    print('iterations', iterations)\n",
    "\n",
    "    y_test_pred_total = np.zeros(X_test.shape[0])\n",
    "    print(modelos)\n",
    "    for key, model in models.items():\n",
    "        #if key == 'lgb_2':\n",
    "        print('predicting '+ key + ' model')\n",
    "        print(model)\n",
    "        for k in tqdm(range(iterations)):\n",
    "            y_pred_test = model.predict(X_test[k*batch_size:(k+1)*batch_size], num_iteration=model.best_iteration)\n",
    "            y_test_pred_total[k*batch_size:(k+1)*batch_size] += y_pred_test\n",
    "\n",
    "        del model\n",
    "        gc.collect()\n",
    "    \n",
    "\n",
    "    y_test_pred_total /= 3 # prediccion con media de modelos\n",
    "    return y_test_pred_total"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {
    "_cell_guid": "92f0cf7e-3b7d-4ffc-8328-b82859d495c2",
    "_uuid": "9675cb6c-8160-48c7-b606-e2f32dd88dc0"
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 20.5 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "from sklearn.externals import joblib\n",
    "\n",
    "PATH_dataset_modelos = '/kaggle/input/lgb-models-ashrae/modelos_v251/' #_v251\n",
    "\n",
    "if DO_TEST == True and DO_TRAIN == False:\n",
    "    modelos = {}\n",
    "    for file in os.listdir(PATH_dataset_modelos):\n",
    "        if file.endswith('.pkl'):\n",
    "            print(file)\n",
    "            lgb_pickle = joblib.load(PATH_dataset_modelos + file)                                                                                                                                                                                                                                                                           \n",
    "            # load model\n",
    "            #lgb_pickle = lgb.Booster(model_file='/kaggle/input/lgb-models-ashrae/' + file)\n",
    "            #modelos.append(lgb_pickle)\n",
    "            key = file.split('.')[0]\n",
    "            print(key)\n",
    "            modelos[key] = lgb_pickle\n",
    "            \n",
    "    del lgb_pickle\n",
    "    gc.collect()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 0 ns, sys: 0 ns, total: 0 ns\n",
      "Wall time: 11 µs\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "if DO_TEST == True and DO_TRAIN == False:\n",
    "    pred = np.expm1(pred_kfolds_from_dict(X_test, modelos)) \n",
    "\n",
    "if DO_TEST == True and DO_TRAIN == True:\n",
    "    pred = np.expm1(pred_kfolds(X_test, modelos)) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "_cell_guid": "4fff1d8f-c508-4a37-ae5a-37c86bff1410",
    "_uuid": "a986bfe6-6756-4e1d-9fab-8314947d2296"
   },
   "source": [
    "## Submission\n",
    "# # Preparing final file for submission."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {
    "_cell_guid": "8bb62888-2077-4ac5-a3a3-c4d35ebc8ad7",
    "_uuid": "879525f3-7201-4ef1-b259-e565ca2eb350"
   },
   "outputs": [],
   "source": [
    "if DO_TEST == True:\n",
    "    submission = pd.DataFrame({\"row_id\": row_ids, \"meter_reading\": np.clip(pred, 0, a_max=None)})\n",
    "    submission['meter_reading'] = submission['meter_reading'].clip(lower=0)\n",
    "    print(submission.head(15))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {
    "_cell_guid": "eb31ad6b-42dc-418b-afe2-65a60561e920",
    "_uuid": "73206f38-5283-4d59-bfc0-704808f154a8"
   },
   "outputs": [],
   "source": [
    "if DO_TEST == True:\n",
    "    submission.to_csv(\"submission_lgb2_v251_2.csv\", index=False)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {
    "_cell_guid": "",
    "_uuid": ""
   },
   "outputs": [],
   "source": [
    "import pandas as pd\n",
    "site_ids = pd.read_csv(\"../input/ashrae-cities-and-timezones/site_ids.csv\")"
   ]
  }
 ],
 "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.6.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}