Merged.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import os\n",
    "import pandas as pd\n",
    "import glob\n",
    "import h5py\n",
    "import tables\n",
    "import pickle\n",
    "from ast import literal_eval\n",
    "from pandas.tseries.holiday import USFederalHolidayCalendar as calendar\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import seaborn as sns\n",
    "from scipy.optimize import curve_fit\n",
    "\n",
    "from ast import literal_eval\n",
    "from sklearn.ensemble import RandomForestRegressor\n",
    "from sklearn.tree import DecisionTreeRegressor\n",
    "from sklearn.model_selection import train_test_split\n",
    "from sklearn import metrics\n",
    "from sklearn.model_selection import KFold\n",
    "from sklearn.model_selection import TimeSeriesSplit\n",
    "from sklearn.model_selection import GridSearchCV\n",
    "\n",
    "from sklearn import datasets, linear_model\n",
    "from sklearn.metrics import mean_squared_error, r2_score\n",
    "\n",
    "\n",
    "\n",
    "import psutil\n",
    "import sys\n",
    "from multiprocessing import Pool\n",
    "from keras.models import Sequential\n",
    "from keras.layers import Conv2D, MaxPooling2D, Dropout, Flatten, Dense\n",
    "from keras.layers import Convolution1D, MaxPooling1D, AveragePooling1D\n",
    "from keras.layers import Flatten\n",
    "from keras.layers import LSTM\n",
    "from keras.layers.convolutional import Conv1D\n",
    "from keras.layers.convolutional import MaxPooling1D\n",
    "from sklearn.model_selection import train_test_split\n",
    "from keras.callbacks import EarlyStopping\n",
    "from keras.optimizers import *\n",
    "import math\n",
    "from sklearn.metrics import mean_squared_error\n",
    "from sklearn.model_selection import train_test_split\n",
    "import matplotlib.pyplot as plt\n",
    "import h5py\n",
    "from keras.callbacks import ModelCheckpoint\n",
    "from keras.models import load_model\n",
    "from matplotlib import pyplot as plt\n",
    "from keras.layers import TimeDistributed"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Preprocessing Data\n",
    "\n",
    "This script takes raw data, create a pandas df and saves it to disk as a hdf5 file called `data.h5`. Before running this code, make sure your file directory looks like this...\n",
    "\n",
    "```\n",
    ".\n",
    "│   Data Merged.ipynb\n",
    "│   ...    \n",
    "│\n",
    "└───Data\n",
    "    │   search_grid.csv\n",
    "    │   impact_score_2017.csv\n",
    "    │   impact_score_2017.csv\n",
    "    |   impact_score_2019_public_test.csv\n",
    "    |\n",
    "    └───GFS_2017\n",
    "    |   │   gfs_4_20170101_0000_000.csv\n",
    "    |   │   ...\n",
    "    |\n",
    "    └───GFS_2018\n",
    "    |   │   gfs_4_20180101_0000_000.csv\n",
    "    |   │   ...     \n",
    "    | \n",
    "    └───GFS_2019\n",
    "        │   gfs_4_20190101_0000_000.csv\n",
    "        │   ...     \n",
    "   \n",
    "\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Loading Provided Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "#### Path to all data files\n",
    "\n",
    "# All weather\n",
    "weather17 = [\"./Data/GFS_2017/\" + f for f in os.listdir(\"./Data/GFS_2017\")]\n",
    "weather18 = [\"./Data/GFS_2018/\" + f for f in os.listdir(\"./Data/GFS_2018\")]\n",
    "weather19 = [\"./Data/GFS_2019/\" + f for f in os.listdir(\"./Data/GFS_2019\")]\n",
    "allweather = weather17 + weather18 + weather19\n",
    "\n",
    "# All impact\n",
    "allimpact = [\"./Data/\" + f for f in os.listdir(\"./Data\") if f.startswith('impact')]\n",
    "\n",
    "# Loc2zip map\n",
    "grid = \"./Data/search_grid.csv\""
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load weather into 1 df\n",
    "# Concat into a giant dataframe\n",
    "weather = pd.concat(map(pd.read_csv, allweather), ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Load and flatten loc2zip map\n",
    "loc2zip = pd.read_csv(grid, converters={\"mapped_zipcodes\": literal_eval})\n",
    "loc2zip = loc2zip.explode(\"mapped_zipcodes\").reset_index()\n",
    "loc2zip['mapped_zipcodes'] = pd.to_numeric(loc2zip['mapped_zipcodes'])\n",
    "#loc2zip.rename(columns={\"mapped_zipcodes\":\"zip5\"}, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "### Load impact into 1 df\n",
    "impact = pd.concat(map(pd.read_csv, allimpact), ignore_index=True).drop([\"Unnamed: 0\"], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Join impact and loc2zip\n",
    "impact_loc = impact.join(loc2zip.set_index('mapped_zipcodes'), on='zip5').drop([\"index\"], axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Join impact and weather in loc and date\n",
    "# Preprocessing on join columns \n",
    "impact_loc[\"date_key\"] = pd.to_datetime(impact_loc[\"date_key\"])\n",
    "weather[\"Date\"] = pd.to_datetime(weather[\"Date\"], format=\"%Y%m%d\")\n",
    "\n",
    "merged = impact_loc.merge(weather, left_on=[\"date_key\", \"grid_lat\", \"grid_lon\"], right_on=[\"Date\", \"lat\", \"lng\"], how='left')\n",
    "merged = merged.sort_values(by=\"Date\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "#impact_loc.to_hdf(\"impact_scores.h5\", key='df', mode='a')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Add Empty Rows for Missing Weather Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "left = merged.loc[pd.isna(merged[\"Date\"])].sort_values(\"date_key\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add empty rows for impact scores that don't have any weather data\n",
    "missing = left.loc[left.index.repeat(4)]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add empty rows for impact scores that have less than 4 weather timestamps\n",
    "complete = merged.loc[pd.notna(merged[\"Date\"])]\n",
    "g = complete.groupby([\"zip5\", \"date_key\"])\n",
    "miss3 = g.filter(lambda x: len(x) == 1)\n",
    "miss2 = g.filter(lambda x: len(x) == 2)\n",
    "miss1 = g.filter(lambda x: len(x) == 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(637, 123)\n",
      "(1911, 123)\n"
     ]
    }
   ],
   "source": [
    "# missing 3 timestamps\n",
    "miss3.iloc[:, 5:] = np.NaN\n",
    "add3 = miss3.loc[miss3.index.repeat(3)]\n",
    "print(miss3.shape)\n",
    "print(add3.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {
    "scrolled": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(182, 123)\n",
      "(182, 123)\n"
     ]
    }
   ],
   "source": [
    "# missing 2 timestamps\n",
    "add2 = miss2\n",
    "add2.iloc[:, 5:] = np.NaN\n",
    "print(miss2.shape)\n",
    "print(add2.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(1092, 123)\n",
      "(364, 123)\n"
     ]
    }
   ],
   "source": [
    "# missing 2 timestamp\n",
    "miss1_ = miss1.iloc[:, 0:5].drop_duplicates([\"date_key\", \"zip5\"])\n",
    "add1 = miss1.loc[miss1_.index.repeat(1)]\n",
    "add1.iloc[:, 5:] = np.NaN\n",
    "print(miss1.shape)\n",
    "print(add1.shape)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add missing value rows to dataframe\n",
    "missing = missing.append(add1)\n",
    "missing = missing.append(add2)\n",
    "missing = missing.append(add3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0, 123)\n",
      "(394700, 123)\n",
      "(0, 123)\n"
     ]
    }
   ],
   "source": [
    "# Verify that each \n",
    "data = complete.append(missing)\n",
    "gg = data.groupby([\"zip5\", \"date_key\"])\n",
    "print((gg.filter(lambda x: len(x) < 4)).shape)\n",
    "print((gg.filter(lambda x: len(x) == 4)).shape)\n",
    "print((gg.filter(lambda x: len(x) > 4).shape))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Impute Missing Weather Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# use haversine formula to calculate shortest distance over the earth's surface\n",
    "# https://en.wikipedia.org/wiki/Haversine_formula\n",
    "\n",
    "def latlong_distance(lat1, long1, lat2, long2):\n",
    "    lat1 = np.radians(lat1)\n",
    "    lat2 = np.radians(lat2)\n",
    "    lat_delta = lat2 - lat1\n",
    "    long_delta = np.radians(long2 - long1)\n",
    "    #a = sin^2(latitude delta / 2) + cos latitude 1 * cos latitude 2\n",
    "    #* sin^2(longitude delta / 2)\n",
    "    a = np.sin(lat_delta / 2) ** 2 + np.cos(lat1) * np.cos(lat2) * np.sin(long_delta / 2) ** 2\n",
    "\n",
    "    #c = 2 * arcsin(a^0.5)\n",
    "    c = 2 * np.arcsin(np.sqrt(a))\n",
    "    \n",
    "    #d = R * c\n",
    "    R = 6371 #earth's radius = 6371km\n",
    "    d = R * c\n",
    "    return d"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def closest_n_search(latlong, latlongs, **kwargs):\n",
    "    # find n closest latlongs to the given coordinates and corresponding weights\n",
    "    # input = tuple of target coordinates, list of tuples [(lat, long)], n\n",
    "    # output = dict, key = weights, values = tuple(lat, long)\n",
    "    target_lat, target_long = latlong\n",
    "    min_latlong = tuple()\n",
    "    min_distance = float('inf')\n",
    "    dist = {}\n",
    "    \n",
    "    # evaluate each lat long pairs\n",
    "    for (lat, long) in latlongs:\n",
    "        # only evaluate lat long coordinates within a certain range\n",
    "        if np.abs(lat - target_lat) <= 10 and np.abs(long - target_long) <= 10:\n",
    "            dist = latlong_distance(target_lat, target_long, lat, long)\n",
    "\n",
    "            if len(min_latlong) < 1 or dist <= min_distance:\n",
    "                min_latlong = (lat, long)\n",
    "                min_distance = dist\n",
    "    return min_latlong, min_distance"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Pre-calculate closest n neighbors and corresponding weights of all coordinates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Pre-compute n nearest neighbor based on distance\n",
    "# Select unique lat long combinations from the data\n",
    "latlongs = data[['grid_lat', 'grid_lon']].values\n",
    "latlongs = [(ll[0], ll[1]) for ll in latlongs]\n",
    "latlongs = list(set(latlongs))\n",
    "len(latlongs)\n",
    "\n",
    "latlong_matches = {}\n",
    "for i in range(len(latlongs)):\n",
    "    latlong_match = closest_n_search(latlongs[i], latlongs[:i] + latlongs[i + 1:])\n",
    "    latlong_matches[latlongs[i]] = latlong_match\n",
    "    # return dict\n",
    "    # key = (target_lat, target_long), value = (n_lat, n_long, dist)\n",
    "        \n",
    "print(\"example of coordinates at a similar distance to multiple points, target point :(36.0, -86.5), closest neighbor:\",\\\n",
    "latlong_matches[(36.0, -86.5)])\n",
    "\n",
    "dist_closest_n = [i[1] for i in list(latlong_matches.values())]\n",
    "\n",
    "print(f\"Average distance with geo-imputing ref (km): {np.mean(dist_closest_n)}\")\n",
    "print(f\"Min. distance with geo-imputing ref (km): {np.min(dist_closest_n)}\")\n",
    "print(f\"Max. distance with geo-imputing ref (km): {np.max(dist_closest_n)}\")\n",
    "\n",
    "#\"Distribution of distance to closest neighbor (km)\", plt.hist(dist_closest_n, bins=[0,100,200,300,400,500,600,700,800])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sns.set(style=\"darkgrid\")\n",
    "fig, ax = plt.subplots(1, 1, figsize=(9, 5))\n",
    "ax = sns.distplot(dist_closest_n,bins=9,kde=False)\n",
    "ax.set(xlabel='Distance to closest FC (km)', ylabel='Count')\n",
    "fig.suptitle(\"Histogram of distance to closest neighbor of FCs\", fontsize=15)\n",
    "#plt.show()\n",
    "fig.savefig('Figure_3.png')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Resize data to ensure there are data records for all dates, zip5 and time (4 times a day)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for index, row in data.iterrows():\n",
    "    f = (data['date_key'] == row['date_key']) & (data['zip5'] == row['zip5'])\n",
    "    if pd.isna(row['Date']):\n",
    "        data.at[index,'Date'] = row['date_key']\n",
    "    if pd.isna(row['Time']):\n",
    "        missing_times = list(set(time) - set(data[f]['Time']))\n",
    "        data.at[index,'Time'] = missing_times[0]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Drop columns with > 80% missing data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = data.drop(columns = [\"ForecastRange\",\"Categorical_Snow_surface\",\\\n",
    "                              \"Composite_reflectivity_entire_atmosphere\",\n",
    "                              \"Graupel_snow_pellets_hybrid\",\\\n",
    "                              \"Graupel_snow_pellets_isobaric\",\\\n",
    "                              \"Snow_mixing_ratio_hybrid\",\\\n",
    "                              \"Snow_mixing_ratio_isobaric\",\\\n",
    "                              'Geopotential_height_potential_vorticity_surface', \\\n",
    "                              'u_component_of_wind_potential_vorticity_surface',\\\n",
    "                              'v_component_of_wind_potential_vorticity_surface',\\\n",
    "                              'Categorical_Ice_Pellets_surface',\\\n",
    "                              'Ice_water_mixing_ratio_hybrid',\\\n",
    "                              'Ice_water_mixing_ratio_isobaric',\\\n",
    "                              'Precipitation_rate_surface',\\\n",
    "                              'Vertical_velocity_geometric_isobaric',\\\n",
    "                              'Ice_growth_rate_altitude_above_msl',\\\n",
    "                              'Land_sea_coverage_nearest_neighbor_land1sea0_surface',\\\n",
    "                              'Pressure_potential_vorticity_surface',\\\n",
    "                              'Temperature_potential_vorticity_surface',\\\n",
    "                              'Vertical_Speed_Shear_potential_vorticity_surface',\\\n",
    "                              'Rain_mixing_ratio_hybrid',\\\n",
    "                              'Total_cloud_cover_isobaric',\\\n",
    "                              'Cloud_mixing_ratio_hybrid',\\\n",
    "                              'Categorical_Freezing_Rain_surface',\\\n",
    "                              'Categorical_Rain_surface', 'Visibility_surface'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Spatial Imputation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def spatial_impute(lat, long, date, time, colname, data):\n",
    "    n_latlong, dist = latlong_matches[(lat, long)]\n",
    "    #print('\\n=================')\n",
    "    #print(\"target point:\" , \"(\" ,lat,long, \"), closest neighbor:\", n_latlong, \", distance:\", dist)\n",
    "    b = data[(data['grid_lat'] == n_latlong[0]) & (data['grid_lon'] == n_latlong[1])\\\n",
    "        & (data['date_key'] == date) & (data['Time'] == time)]\n",
    "    #print('\\nData of closet neighbor\\n',b,'\\n')\n",
    "    #print(np.sum(b))\n",
    "    if b.shape[0] == 0:\n",
    "        output = np.nan\n",
    "    else: \n",
    "        output = b[colname].values[0]\n",
    "    #print('imputed value:', output, 'distance',dist)\n",
    "    #print('=================\\n')\n",
    "    return output, dist"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Code for all features to impute\n",
    "\n",
    "cols_to_impute = ['Snow_depth_surface', 'Haines_Index_surface',\\\n",
    "'u_component_of_wind_altitude_above_msl', 'v_component_of_wind_altitude_above_msl',\\\n",
    "'Soil_temperature_depth_below_surface_layer', 'Temperature_altitude_above_msl',\\\n",
    "'Volumetric_Soil_Moisture_Content_depth_below_surface_layer', 'Field_Capacity_surface', \\\n",
    "'Water_equivalent_of_accumulated_snow_depth_surface', 'Wilting_Point_surface', \\\n",
    "'Field_Capacity_surface', \"Water_equivalent_of_accumulated_snow_depth_surface\",\"Wilting_Point_surface\"]\n",
    "\n",
    "#Find index of column to impute\n",
    "for colname in cols_to_impute:\n",
    "    #colname = cols_to_impute[0]\n",
    "    colind = list(data.columns).index(colname)\n",
    "    print('Impute for column', colname)\n",
    "\n",
    "    for i in range(data.shape[0]):\n",
    "        data_target = data.iloc[i,:]\n",
    "        #print(data_target)\n",
    "        if np.isnan(data_target[colname]): #If data is missing, we need to impute\n",
    "            print(i, '=>impute')\n",
    "            lat, long, date, time = data_target['grid_lat'], data_target['grid_lon'], data_target['date_key'], data_target['Time']\n",
    "            i_value, dist = spatial_impute(lat, long, date, time, colname, data)\n",
    "\n",
    "            #Sucessful spatial imputation is when the imputed value is non-nan and the distance to the closest neighbor is < 200\n",
    "            #The code currently set imputed value to Nan if the closest neighbor does not have value. \n",
    "            #However, it still returns a real value even when the distance to the closest neighbor is > 200. Code below is to fix it:\n",
    "            if dist > 200:\n",
    "                i_value = np.nan\n",
    "\n",
    "            #Put the imputed values and flags inside the original dataframe:\n",
    "            data.iloc[i,-3] = 1 #Update the column \"any_impute_col\" value to 1\n",
    "            data.iloc[i,colind] = i_value\n",
    "        #else:\n",
    "        #    print(i, '=> no impute')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Temporal Interpolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#data = pd.read_hdf('imputeddata.h5')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Time-series interpolation\n",
    "pd.to_datetime(data['datetime'])\n",
    "data.set_index('datetime')\n",
    "data.head()\n",
    "data_interpolate = data.groupby(['grid_lat', 'grid_lon']).apply(lambda x: x.interpolate(method='linear', limit_direction = 'both'))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Temporal Extrapolation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "extrapcols = []\n",
    "cols = data_interpolate.columns\n",
    "for col in cols:\n",
    "    if (sum(data_interpolate[col].isna()) != 0) & (col != \"impact_score\"): \n",
    "        extrapcols.append(col)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "testdat = data_interpolate.copy()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in extrapcols:\n",
    "    zips = data_interpolate[pd.isna(data_interpolate[col])][\"zip5\"].unique()\n",
    "    for zipcode in zips:\n",
    "        print(\"for \" + str(zipcode) + \" column: \" + col)\n",
    "        datum = data_interpolate[data_interpolate[\"zip5\"] == zipcode][col].copy().reset_index()\n",
    "           \n",
    "        # Function to curve fit to the data\n",
    "        def func(x, a, b, c, d):\n",
    "            return a * (x ** 3) + b * (x ** 2) + c * x + d\n",
    "\n",
    "        # Initial parameter guess, just to kick off the optimization\n",
    "        guess = (0.5, 0.5, 0.5, 0.5)\n",
    "\n",
    "        # Create copy of data to remove NaNs for curve fitting\n",
    "        fit_datum = datum.dropna()\n",
    "\n",
    "        # Place to store function parameters for each column\n",
    "        col_params = {}\n",
    "\n",
    "        # Get x & y and fit curve\n",
    "        x = fit_datum.index.astype(float).values\n",
    "        y = fit_datum[col].values\n",
    "        # Curve fit column and get curve parameters\n",
    "        params = curve_fit(func, x, y, guess)\n",
    "        # Store optimized parameters\n",
    "        col_params[col] = params[0]\n",
    "\n",
    "        # Extrapolate for the columns\n",
    "        # Get the index values for NaNs in the column\n",
    "        x = datum[pd.isnull(datum[col])].index.astype(float).values\n",
    "        # Extrapolate those points with the fitted function\n",
    "        datum[col][x] = func(x, *col_params[col])\n",
    "        print(\"testdat before assignment\", testdat.loc[testdat[\"zip5\"] == zipcode, col])\n",
    "        print(\"datum extrap\", datum[col])\n",
    "        idx = datum[\"index\"]\n",
    "        testdat[col][idx] = datum[col]\n",
    "        print(\"testdat after assignment\", testdat.loc[testdat[\"zip5\"] == zipcode, col])\n",
    "        print(\"\\n\")\n",
    "        #data_interpolate.loc[data_interpolate[\"zip5\"] == zipcode][col] = datum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# clipping data to the right range\n",
    "cols_to_clip_0= ['Snow_depth_surface',\\\n",
    "'Soil_temperature_depth_below_surface_layer', 'Temperature_altitude_above_msl',\\\n",
    "'Volumetric_Soil_Moisture_Content_depth_below_surface_layer','Field_Capacity_surface', \\\n",
    "\"Water_equivalent_of_accumulated_snow_depth_surface\",\"Wilting_Point_surface\"]\n",
    "cols_to_clip_2 = ['Haines_Index_surface']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in cols_to_clip_0:\n",
    "    testdat[col].clip(lower=0, inplace=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for col in cols_to_clip_2:\n",
    "    testdat[col].clip(lower=2, upper=6, inplace=True)\n",
    "    testdat[col] = testdat[col].round()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "data = testdat"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Adding Features from External Data Scources"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "### Add state, region, and urban/suburban information from FAR codes dataset\n",
    "data = data.sort_values([\"zip5\", 'date_key', 'Time'])\n",
    "x = data['zip5']\n",
    "# Read in Far codes\n",
    "far = pd.read_excel('FARcodesZIPdata2010WithAKandHI.xlsx', sheet_name = 1)\n",
    "cols = ['ZIP', 'density', 'state']\n",
    "far = far[cols]\n",
    "new_row = {'ZIP': 2722, 'density': 0, 'state': 'MA'}\n",
    "far = far.append(new_row, ignore_index=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add density and state\n",
    "data = data.merge(far, left_on = 'zip5', right_on = 'ZIP')\n",
    "y = data['zip5']"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "# State to region dictionary\n",
    "New_England_Northeast = ['CT', 'ME', 'MA', 'NH', 'RI', 'VT']\n",
    "Mid_Atlantic_Northeast = ['NJ', 'NY', 'PA']\n",
    "East_North_Central_Midwest = ['IL', 'IN', 'MI', 'OH', 'WI'] \n",
    "West_North_Central_Midwest = ['IA', 'KS', 'MN', 'MO', 'NE', 'ND', 'SD']\n",
    "South_Atlantic_South = ['DE', 'FL', 'GA', 'MD', 'NC', 'SC', 'VA', 'DC', 'WV']\n",
    "East_South_Central_South = ['AL', 'KY', 'MS', 'TN']\n",
    "West_South_Central_South = ['AR', 'LA', 'OK', 'TX']\n",
    "Mountain_West = ['AZ', 'CO', 'ID', 'MT', 'NV', 'NM', 'UT', 'WY']\n",
    "Pacific_West = ['AK', 'CA', 'HI', 'OR', 'WA']\n",
    "\n",
    "# Add region\n",
    "data['Region'] = ['New England Northeast' if x in New_England_Northeast else\n",
    "                 'Mid Atlantic Northeast' if x in Mid_Atlantic_Northeast else\n",
    "                 'East North Central Midwest' if x in East_North_Central_Midwest else\n",
    "                 'West North Central Midwest' if x in West_North_Central_Midwest else\n",
    "                 'South Atlantic South' if x in South_Atlantic_South else\n",
    "                 'East South Central South' if x in East_South_Central_South else\n",
    "                 'West South Central South' if x in West_South_Central_South else\n",
    "                 'Mountain West' if x in Mountain_West else\n",
    "                 'Pacific West' for x in data['state']]\n",
    "data = data.drop([\"ZIP\"], axis=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add suburban or urban\n",
    "data['USR'] = ['Urban' if x > 3000 else 'Rural' if x < 1000 else 'Suburban' for x in data['density']]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add day of the week\n",
    "data['Weekday'] = data['date_key'].dt.dayofweek\n",
    "data['is_weekend'] = (data['date_key'].dt.dayofweek < 5).astype(int)\n",
    "\n",
    "# Add date, month and year\n",
    "data['day'] = data['date_key'].dt.day\n",
    "data['month'] = data['date_key'].dt.month\n",
    "data['year'] = data['date_key'].dt.year"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add federal holidays\n",
    "cal = calendar()\n",
    "holidays = cal.holidays(start=\"2017-01-01\", end=\"2019-12-31\")\n",
    "\n",
    "data['is_holiday'] = (data[\"date_key\"].isin(holidays)).astype(int)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add categorical variables for date and zip\n",
    "data['year'] = pd.Categorical(data['year'])\n",
    "data['month'] = pd.Categorical(data['month'])\n",
    "data['day'] = pd.Categorical(data['day'])\n",
    "data['zip5'] = pd.Categorical(data['zip5'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Add binary variables for state, rural/urban, region, weekday, year, month, and day of week\n",
    "cols = ['state', 'USR', 'Region', 'Weekday', 'year', 'month', 'day']\n",
    "df = data[cols]\n",
    "dataDummies = pd.get_dummies(df,drop_first=True)\n",
    "data = pd.concat([data, dataDummies], axis = 1)\n",
    "data = data.drop(['state', 'USR', 'Region', 'year', 'month', 'day'], axis = 1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>date_key</th>\n",
       "      <th>zip5</th>\n",
       "      <th>impact_score</th>\n",
       "      <th>grid_lat</th>\n",
       "      <th>grid_lon</th>\n",
       "      <th>Date</th>\n",
       "      <th>Time</th>\n",
       "      <th>ForecastRange</th>\n",
       "      <th>x</th>\n",
       "      <th>y</th>\n",
       "      <th>...</th>\n",
       "      <th>density</th>\n",
       "      <th>state</th>\n",
       "      <th>Region</th>\n",
       "      <th>USR</th>\n",
       "      <th>Weekday</th>\n",
       "      <th>is_weekend</th>\n",
       "      <th>day</th>\n",
       "      <th>month</th>\n",
       "      <th>year</th>\n",
       "      <th>is_holiday</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>2017-01-01</td>\n",
       "      <td>2722</td>\n",
       "      <td>20.268081</td>\n",
       "      <td>41.5</td>\n",
       "      <td>-71.0</td>\n",
       "      <td>2017-01-01</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>578.0</td>\n",
       "      <td>97.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>MA</td>\n",
       "      <td>New England Northeast</td>\n",
       "      <td>Rural</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2017</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>2017-01-01</td>\n",
       "      <td>2722</td>\n",
       "      <td>20.268081</td>\n",
       "      <td>41.5</td>\n",
       "      <td>-71.0</td>\n",
       "      <td>2017-01-01</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>578.0</td>\n",
       "      <td>97.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>MA</td>\n",
       "      <td>New England Northeast</td>\n",
       "      <td>Rural</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2017</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>2017-01-01</td>\n",
       "      <td>2722</td>\n",
       "      <td>20.268081</td>\n",
       "      <td>41.5</td>\n",
       "      <td>-71.0</td>\n",
       "      <td>2017-01-01</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>578.0</td>\n",
       "      <td>97.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>MA</td>\n",
       "      <td>New England Northeast</td>\n",
       "      <td>Rural</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2017</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>2017-01-01</td>\n",
       "      <td>2722</td>\n",
       "      <td>20.268081</td>\n",
       "      <td>41.5</td>\n",
       "      <td>-71.0</td>\n",
       "      <td>2017-01-01</td>\n",
       "      <td>18.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>578.0</td>\n",
       "      <td>97.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>MA</td>\n",
       "      <td>New England Northeast</td>\n",
       "      <td>Rural</td>\n",
       "      <td>6</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>2017</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>2017-01-02</td>\n",
       "      <td>2722</td>\n",
       "      <td>16.868994</td>\n",
       "      <td>41.5</td>\n",
       "      <td>-71.0</td>\n",
       "      <td>2017-01-02</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>578.0</td>\n",
       "      <td>97.0</td>\n",
       "      <td>...</td>\n",
       "      <td>0.000000</td>\n",
       "      <td>MA</td>\n",
       "      <td>New England Northeast</td>\n",
       "      <td>Rural</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>2</td>\n",
       "      <td>1</td>\n",
       "      <td>2017</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>...</th>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "      <td>...</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>394695</th>\n",
       "      <td>2019-12-30</td>\n",
       "      <td>98421</td>\n",
       "      <td>NaN</td>\n",
       "      <td>47.5</td>\n",
       "      <td>-122.5</td>\n",
       "      <td>2019-12-30</td>\n",
       "      <td>18.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>950.0</td>\n",
       "      <td>170.0</td>\n",
       "      <td>...</td>\n",
       "      <td>208.237785</td>\n",
       "      <td>WA</td>\n",
       "      <td>Pacific West</td>\n",
       "      <td>Rural</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "      <td>30</td>\n",
       "      <td>12</td>\n",
       "      <td>2019</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>394696</th>\n",
       "      <td>2019-12-31</td>\n",
       "      <td>98421</td>\n",
       "      <td>NaN</td>\n",
       "      <td>47.5</td>\n",
       "      <td>-122.5</td>\n",
       "      <td>2019-12-31</td>\n",
       "      <td>0.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>950.0</td>\n",
       "      <td>170.0</td>\n",
       "      <td>...</td>\n",
       "      <td>208.237785</td>\n",
       "      <td>WA</td>\n",
       "      <td>Pacific West</td>\n",
       "      <td>Rural</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>31</td>\n",
       "      <td>12</td>\n",
       "      <td>2019</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>394697</th>\n",
       "      <td>2019-12-31</td>\n",
       "      <td>98421</td>\n",
       "      <td>NaN</td>\n",
       "      <td>47.5</td>\n",
       "      <td>-122.5</td>\n",
       "      <td>2019-12-31</td>\n",
       "      <td>6.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>950.0</td>\n",
       "      <td>170.0</td>\n",
       "      <td>...</td>\n",
       "      <td>208.237785</td>\n",
       "      <td>WA</td>\n",
       "      <td>Pacific West</td>\n",
       "      <td>Rural</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>31</td>\n",
       "      <td>12</td>\n",
       "      <td>2019</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>394698</th>\n",
       "      <td>2019-12-31</td>\n",
       "      <td>98421</td>\n",
       "      <td>NaN</td>\n",
       "      <td>47.5</td>\n",
       "      <td>-122.5</td>\n",
       "      <td>2019-12-31</td>\n",
       "      <td>12.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>950.0</td>\n",
       "      <td>170.0</td>\n",
       "      <td>...</td>\n",
       "      <td>208.237785</td>\n",
       "      <td>WA</td>\n",
       "      <td>Pacific West</td>\n",
       "      <td>Rural</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>31</td>\n",
       "      <td>12</td>\n",
       "      <td>2019</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>394699</th>\n",
       "      <td>2019-12-31</td>\n",
       "      <td>98421</td>\n",
       "      <td>NaN</td>\n",
       "      <td>47.5</td>\n",
       "      <td>-122.5</td>\n",
       "      <td>2019-12-31</td>\n",
       "      <td>18.0</td>\n",
       "      <td>0.0</td>\n",
       "      <td>950.0</td>\n",
       "      <td>170.0</td>\n",
       "      <td>...</td>\n",
       "      <td>208.237785</td>\n",
       "      <td>WA</td>\n",
       "      <td>Pacific West</td>\n",
       "      <td>Rural</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "      <td>31</td>\n",
       "      <td>12</td>\n",
       "      <td>2019</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "<p>394700 rows × 133 columns</p>\n",
       "</div>"
      ],
      "text/plain": [
       "         date_key   zip5  impact_score  grid_lat  grid_lon       Date  Time  \\\n",
       "0      2017-01-01   2722     20.268081      41.5     -71.0 2017-01-01   0.0   \n",
       "1      2017-01-01   2722     20.268081      41.5     -71.0 2017-01-01   6.0   \n",
       "2      2017-01-01   2722     20.268081      41.5     -71.0 2017-01-01  12.0   \n",
       "3      2017-01-01   2722     20.268081      41.5     -71.0 2017-01-01  18.0   \n",
       "4      2017-01-02   2722     16.868994      41.5     -71.0 2017-01-02   0.0   \n",
       "...           ...    ...           ...       ...       ...        ...   ...   \n",
       "394695 2019-12-30  98421           NaN      47.5    -122.5 2019-12-30  18.0   \n",
       "394696 2019-12-31  98421           NaN      47.5    -122.5 2019-12-31   0.0   \n",
       "394697 2019-12-31  98421           NaN      47.5    -122.5 2019-12-31   6.0   \n",
       "394698 2019-12-31  98421           NaN      47.5    -122.5 2019-12-31  12.0   \n",
       "394699 2019-12-31  98421           NaN      47.5    -122.5 2019-12-31  18.0   \n",
       "\n",
       "        ForecastRange      x      y  ...     density  state  \\\n",
       "0                 0.0  578.0   97.0  ...    0.000000     MA   \n",
       "1                 0.0  578.0   97.0  ...    0.000000     MA   \n",
       "2                 0.0  578.0   97.0  ...    0.000000     MA   \n",
       "3                 0.0  578.0   97.0  ...    0.000000     MA   \n",
       "4                 0.0  578.0   97.0  ...    0.000000     MA   \n",
       "...               ...    ...    ...  ...         ...    ...   \n",
       "394695            0.0  950.0  170.0  ...  208.237785     WA   \n",
       "394696            0.0  950.0  170.0  ...  208.237785     WA   \n",
       "394697            0.0  950.0  170.0  ...  208.237785     WA   \n",
       "394698            0.0  950.0  170.0  ...  208.237785     WA   \n",
       "394699            0.0  950.0  170.0  ...  208.237785     WA   \n",
       "\n",
       "                       Region    USR  Weekday  is_weekend  day  month  year  \\\n",
       "0       New England Northeast  Rural        6           0    1      1  2017   \n",
       "1       New England Northeast  Rural        6           0    1      1  2017   \n",
       "2       New England Northeast  Rural        6           0    1      1  2017   \n",
       "3       New England Northeast  Rural        6           0    1      1  2017   \n",
       "4       New England Northeast  Rural        0           1    2      1  2017   \n",
       "...                       ...    ...      ...         ...  ...    ...   ...   \n",
       "394695           Pacific West  Rural        0           1   30     12  2019   \n",
       "394696           Pacific West  Rural        1           1   31     12  2019   \n",
       "394697           Pacific West  Rural        1           1   31     12  2019   \n",
       "394698           Pacific West  Rural        1           1   31     12  2019   \n",
       "394699           Pacific West  Rural        1           1   31     12  2019   \n",
       "\n",
       "        is_holiday  \n",
       "0                0  \n",
       "1                0  \n",
       "2                0  \n",
       "3                0  \n",
       "4                1  \n",
       "...            ...  \n",
       "394695           0  \n",
       "394696           0  \n",
       "394697           0  \n",
       "394698           0  \n",
       "394699           0  \n",
       "\n",
       "[394700 rows x 133 columns]"
      ]
     },
     "execution_count": 60,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Store full dataset"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "data.to_hdf(\"full.h5\", key='df', mode='a')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load data\n",
    "data = pd.read_hdf('full.h5')\n",
    "data.head()\n",
    "# Preprocessing\n",
    "data = data.drop(['zip5'], axis = 1)\n",
    "data = data[data['datetime'].dt.year < 2019]\n",
    "data = data.drop(['x', 'y', 'any_impute_col', 'impute_row', 'datetime'], axis = 1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Training and Test Split"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Split data to training and validation (based on percentage)\n",
    "\n",
    "end_ind = int(0.7*data.shape[0])\n",
    "trainData, testData = data.loc[0:end_ind,:], data.loc[end_ind+1:data.shape[0],:]\n",
    "print(trainData.shape, testData.shape)\n",
    "\n",
    "X_train = trainData.drop([\"impact_score\"], axis = 1)\n",
    "X_test = testData.drop([\"impact_score\"], axis = 1)\n",
    "\n",
    "Y_train = trainData[\"impact_score\"]\n",
    "Y_test = testData[\"impact_score\"]\n",
    "print('train data: ', X_train.shape, Y_train.shape, '\\ntest data: ', X_test.shape, Y_test.shape)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Regression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Create linear regression object\n",
    "regr = linear_model.LinearRegression()\n",
    "\n",
    "# Train the model using the training sets\n",
    "regr.fit(X_train.values, Y_train.values)\n",
    "\n",
    "# Make predictions using the testing set\n",
    "y_pred_MLR_ori = regr.predict(X_test.values)\n",
    "\n",
    "# The coefficients\n",
    "#print('Coefficients: \\n', regr.coef_)\n",
    "# The mean squared error\n",
    "\n",
    "y_pred_MLR = y_pred_MLR_ori.flatten()\n",
    "\n",
    "# Clipping values\n",
    "y_pred_MLR[y_pred_MLR<-1] = 0\n",
    "y_pred_MLR[y_pred_MLR>35] = 35\n",
    "\n",
    "print('MLR Error: %.2f'\n",
    "      % mean_squared_error(Y_test.values, y_pred_MLR))\n",
    "# The coefficient of determination: 1 is perfect prediction\n",
    "print('Coefficient of determination: %.2f'\n",
    "      % r2_score(Y_test.values, y_pred_MLR))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random Forest Regression modeling - Parameter Tuning: n_est\n",
    "n_est = 40 #The number of trees in the forest.\n",
    "max_d = 10 #Max depth of the tree\n",
    "min_samples = 100 #The minimum number of samples required to be at a leaf node. \n",
    "\n",
    "n_est_dict_y_e = {}\n",
    "n_est_dict_mse_e = {}\n",
    "\n",
    "for n_est in range(40,220,40):\n",
    "    m1 = RandomForestRegressor(n_estimators = n_est, max_depth = max_d, min_samples_leaf = min_samples)\n",
    "    m1.fit(X_train,Y_train)\n",
    "    y_pred_RF = m1.predict(X_test)\n",
    "    n_est_dict_mse_e[n_est] = metrics.mean_squared_error(Y_test, y_pred_RF)\n",
    "    n_est_dict_y_e[n_est] = y_pred_RF\n",
    "    print(\"RF Error:\", metrics.mean_squared_error(Y_test, y_pred_RF))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random Forest Regression modeling - Parameter Tuning: max_d\n",
    "n_est = 40 #The number of trees in the forest.\n",
    "max_d = 10 #Max depth of the tree\n",
    "min_samples = 100 #The minimum number of samples required to be at a leaf node. \n",
    "\n",
    "n_est_dict_y_d = {}\n",
    "n_est_dict_mse_d = {}\n",
    "\n",
    "for max_d in range(10,22,2):\n",
    "    m1 = RandomForestRegressor(n_estimators = n_est, max_depth = max_d, min_samples_leaf = min_samples)\n",
    "    m1.fit(X_train,Y_train)\n",
    "    y_pred_RF = m1.predict(X_test)\n",
    "    n_est_dict_mse_d[max_d] = metrics.mean_squared_error(Y_test, y_pred_RF)\n",
    "    n_est_dict_y_d[max_d] = y_pred_RF\n",
    "    print(\"RF Error:\", metrics.mean_squared_error(Y_test, y_pred_RF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Random Forest Hyperparameter Tuning\n",
    "\n",
    "Varying number of trees (n_estimators)\n",
    "\n",
    "n_estimators|max_d|MSE|\n",
    "---|---|---|\n",
    "40|10|43.69|\n",
    "80|10|43.69|\n",
    "120|10|43.79|\n",
    "160|10|43.46|\n",
    "200|10|43.75|\n",
    "\n",
    "\n",
    "Varying minimum number of samples per leaf node (n_estimators)\n",
    "\n",
    "n_estimators|max_d|MSE|\n",
    "---|---|---|\n",
    "40|10|43.70|\n",
    "40|12|43.83|\n",
    "40|14|42.69|\n",
    "40|16|41.94|\n",
    "40|18|43.46|\n",
    "40|20|42.43|"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Random Forest Regression modeling - Final modeling\n",
    "n_est = 160 #The number of trees in the forest.\n",
    "max_d = 16 #Max depth of the tree\n",
    "min_samples = 100 #The minimum number of samples required to be at a leaf node. \n",
    "\n",
    "m1 = RandomForestRegressor(n_estimators = n_est, max_depth = max_d, min_samples_leaf = min_samples)\n",
    "m1.fit(X_train,Y_train)\n",
    "y_pred_RF = m1.predict(X_test)\n",
    "print(\"RF Error:\", metrics.mean_squared_error(Y_test, y_pred_RF))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Result Visualization of Regression and Random Forest"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Join testing data: X_test, y_test, y_pred to a single dataframe for visualization\n",
    "y_vals = X_test.copy()\n",
    "y_vals['true_IS'], y_vals['pred_RF'], y_vals['pred_MLR'] = y_test, y_pred_RF, y_pred_MLR\n",
    "\n",
    "y_vals = y_vals.reset_index()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Plot a certain zip code data. i.e. lat = 47, long = -122\n",
    "zip_lat, zip_lon = 47, -122\n",
    "data_zip = y_vals[(y_vals['lat']==zip_lat) & (y_vals['lng']==zip_lon)]\n",
    "\n",
    "print(data_zip.shape)\n",
    "\n",
    "fig, ax = plt.subplots(1,1,figsize=(15,5))\n",
    "ax.plot(data_zip['true_IS'], label = 'True Value')\n",
    "ax.plot(data_zip['pred_RF'], label = 'Random Forest')\n",
    "ax.plot(data_zip['pred_MLR'], label = 'Regression')\n",
    "ax.set(xlabel='Index', ylabel='Impact Score')\n",
    "ax.legend()\n",
    "fig.suptitle('lat '+str(zip_lat)+ ' & long ' + str(zip_lon), fontsize=15)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## CNN"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#reshaping functions Chitwan made\n",
    "def rolling_window(a, window):\n",
    "    a_copy = a.copy()\n",
    "    shape = a_copy.shape[:-1] + (a_copy.shape[-1] - window + 1, window)\n",
    "    strides = a_copy.strides + (a_copy.strides[-1],)\n",
    "    return np.lib.stride_tricks.as_strided(a_copy, shape=shape, strides=strides)\n",
    "\n",
    "def reshape_X(X, time_interval):\n",
    "    XT = np.transpose(X.values)\n",
    "    reshaped_X = np.transpose(rolling_window(XT, time_interval), (1,2,0))\n",
    "    return reshaped_X"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# list of zips\n",
    "zips = list(data.zip5.unique())\n",
    "zips1 = zips[:23]\n",
    "zips2 = zips[23:46]\n",
    "zips3 = zips[46:69]\n",
    "zips4 = zips[69:]\n",
    "zipsx = [29172]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#dictionary of predictions\n",
    "predictions = {}\n",
    "\n",
    "#list of rsme\n",
    "rsme = []\n",
    "\n",
    "#change zips to your set of zips\n",
    "for z in zips2:\n",
    "    lag = 28\n",
    "    # Save model path\n",
    "    checkpointer_path = f\"./models1/{z}.hdf5\"\n",
    "    \n",
    "    #subsetting data for model\n",
    "    test = data[data['zip5'] == z]\n",
    "    test = test.drop(['zip5', 'datetime'], axis =1)\n",
    "    test_17_18 = test[test['year_2019'] != 1]\n",
    "    \n",
    "    #19 test\n",
    "    test_19 = test[test['year_2019'] == 1]\n",
    "    test_19 = test_17_18[-lag+1:].append(test_19, ignore_index = True)\n",
    "    test_19 = test_19.drop(['impact_score'], axis = 1)\n",
    "    test_19 = reshape_X(test_19, lag)\n",
    "    \n",
    "    #final y list\n",
    "    y = test_17_18['impact_score']\n",
    "    y = y[lag-1:]\n",
    "    \n",
    "    #final X\n",
    "    X = test_17_18.drop(['impact_score'], axis = 1)\n",
    "    X = reshape_X(X, lag)\n",
    "    \n",
    "    #data split\n",
    "    test_size = int(0.33 * X.shape[0])\n",
    "    X_train, X_test, y_train, y_test = X[:-test_size], X[-test_size:], y[:-test_size], y[-test_size:]\n",
    "    \n",
    "    checkpointer = ModelCheckpoint(checkpointer_path, monitor='val_loss', verbose=1, save_best_only=True)\n",
    "    \n",
    "    #model architecture\n",
    "    model = Sequential((\n",
    "        Convolution1D(input_shape=(lag, 173), \n",
    "                      kernel_size=4, activation=\"relu\", filters=64),\n",
    "        MaxPooling1D(),    \n",
    "        Dropout(0.2),\n",
    "        Convolution1D(input_shape=(lag-4+1, 173), \n",
    "                      kernel_size=4, activation=\"relu\", filters=64),\n",
    "        MaxPooling1D(),    \n",
    "        Dropout(0.3),\n",
    "        Flatten(),\n",
    "        Dense(50, activation='relu'),\n",
    "        Dense(1),\n",
    "    ))\n",
    "    opt = Adam(lr=0.0005)\n",
    "    \n",
    "    #output of model while running\n",
    "    model.compile(loss='mean_squared_error', optimizer=opt, metrics=['mse'])\n",
    "    model.summary()\n",
    "    \n",
    "    #stop to prevent overfitting\n",
    "    #es = EarlyStopping(monitor='val_loss', mode='min', verbose=1, patience=20)\n",
    "    \n",
    "    #model fit and test\n",
    "    val = model.fit(X_train, y_train, epochs=100, batch_size=5, validation_data=(X_test, y_test), callbacks=[checkpointer])\n",
    "    pred = model.predict(X_test)\n",
    "    testScore = math.sqrt(mean_squared_error(y_test,pred))\n",
    "    rsme.append(testScore)\n",
    "    \n",
    "    #model predictions\n",
    "    model = load_model(f\"./models1/{z}.hdf5\")\n",
    "    pred_19 = model.predict(test_19)\n",
    "    predictions[z] = pred_19"
   ]
  }
 ],
 "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.7.0"
  },
  "nbdime-conflicts": {
   "local_diff": [
    {
     "key": "language_info",
     "op": "add",
     "value": {
      "codemirror_mode": {
       "name": "ipython",
       "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.0"
     }
    }
   ],
   "remote_diff": [
    {
     "key": "language_info",
     "op": "add",
     "value": {
      "codemirror_mode": {
       "name": "ipython",
       "version": 3
      },
      "file_extension": ".py",
      "mimetype": "text/x-python",
      "name": "python",
      "nbconvert_exporter": "python",
      "pygments_lexer": "ipython3",
      "version": "3.7.3"
     }
    }
   ]
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}