{"cells":[{"metadata":{"_kg_hide-input":true,"trusted":true},"cell_type":"code","source":"import numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os, gc\nimport random\nimport datetime\n\nfrom tqdm import tqdm_notebook as tqdm\n\n# matplotlib and seaborn for plotting\nimport matplotlib\nimport matplotlib.pyplot as plt\nimport seaborn as sns\n\n\n\nfrom pandas.plotting import register_matplotlib_converters\nregister_matplotlib_converters()\n\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import KFold, StratifiedKFold\nfrom sklearn import preprocessing\n\nimport lightgbm as lgb","execution_count":null,"outputs":[]},{"metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_kg_hide-input":true,"_kg_hide-output":true,"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"cell_type":"code","source":"path = '../input/ashrae-energy-prediction'\n# Input data files are available in the \"../input/\" directory.\nfor dirname, _, filenames in os.walk(path):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Load data and display samples"},{"metadata":{"trusted":true},"cell_type":"code","source":"%%time\n# unimportant features (see importance below)\nunimportant_cols = ['wind_direction', 'wind_speed', 'sea_level_pressure']\ntarget = 'meter_reading'\n\ndef load_data(source='train', path=path):\n    ''' load and merge all tables '''\n    assert source in ['train', 'test']\n    \n    building = pd.read_csv(f'{path}/building_metadata.csv', dtype={'building_id':np.uint16, 'site_id':np.uint8})\n    weather  = pd.read_csv(f'{path}/weather_{source}.csv', parse_dates=['timestamp'],\n                                                           dtype={'site_id':np.uint8, 'air_temperature':np.float16,\n                                                                  'cloud_coverage':np.float16, 'dew_temperature':np.float16,\n                                                                  'precip_depth_1_hr':np.float16},\n                                                           usecols=lambda c: c not in unimportant_cols)\n    df = pd.read_csv(f'{path}/{source}.csv', dtype={'building_id':np.uint16, 'meter':np.uint8}, parse_dates=['timestamp'])\n    df = df.merge(building, on='building_id', how='left')\n    df = df.merge(weather, on=['site_id', 'timestamp'], how='left')\n    return df\n\n# load and display some samples\ntrain = load_data('train')\ntrain.sample(7)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Test data"},{"metadata":{"trusted":true},"cell_type":"code","source":"test = load_data('test')\ntest.sample(7)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"## Function to reduce the DF size\ndef reduce_memory_usage(df, verbose=True):\n    numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']\n    start_mem = df.memory_usage().sum() / 1024**2    \n    for col in df.columns:\n        col_type = df[col].dtypes\n        if col_type in numerics:\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 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    end_mem = df.memory_usage().sum() / 1024**2\n    if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))\n    return df","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"reduce_memory_usage(test)\nreduce_memory_usage(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# target's log-log histogram:\nax = np.log1p(train.meter_reading).hist()\nax.set_yscale('log')\n\n# describe raw values first\ntrain.meter_reading.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# check the distribution in the types of meters\nmeters = train.groupby('building_id').meter.nunique()\nplt.title('Distribution of types of meters\\n{0:electricity, 1:water, 2:steam, 3:hotwater}') # from the official starter kernel\n_ = meters.hist()\n# from the graphs it looks like steam and hotwater are reversed (e.g.: 3:steam, 2:hotwater) but that shouldn't make any difference to the model","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### display a single time series (notice measurement errors and discontinuities)"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"building_id = 1258  # a building with all 4 meters\nmeters = train[train['building_id'] == building_id].meter.nunique()\n\nfor meter in range(meters):\n    fig, ax = plt.subplots()\n    plt.title(f'Building {building_id} Meter {meter}')\n    ax2 = ax.twinx()\n    # plot meter_reading\n    idx = (train['building_id'] == building_id) & (train['meter'] == meter)\n    dates = matplotlib.dates.date2num(train.loc[idx, 'timestamp'])\n    ax2.plot_date(dates, train.loc[idx, 'meter_reading'], '-', label='meter_reading', alpha=0.8)\n    # plot air_temperature\n    dates = matplotlib.dates.date2num(train.loc[train['building_id'] == building_id, 'timestamp'])\n    ax.plot_date(dates, train.loc[train['building_id'] == building_id, 'air_temperature'], '.', color='tab:cyan', label='air_temperature')\n    ax.set_ylabel('air_temperature'); ax2.set_ylabel('meter_reading')\n    ax.legend(loc='upper left'); ax2.legend(loc='upper right')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### now let's see what's the expected prediction in the test set for the same building"},{"metadata":{"trusted":true,"_kg_hide-input":true},"cell_type":"code","source":"meter = 1 # pick a meter\n\ntrain_sample = train[(train['building_id'] == building_id) & (train['meter'] == meter)]  # same train sample as above\n\ntest['meter_reading'] = 0.0\ntest_sample = test[(test['building_id'] == building_id) & (test['meter'] == meter)]  # and the same meter in the test set\n\nfig, ax = plt.subplots(figsize=(16,4))\nplt.title(f'Meter {meter}')\nax.xaxis.set_tick_params(rotation=30, labelsize=10)\nax2 = ax.twinx()\n\n# plot training sample\ndates = matplotlib.dates.date2num(train_sample['timestamp'])\nax2.plot_date(dates, train_sample['meter_reading'], '-', label='train', alpha=0.8)\nax.plot_date(dates, train_sample['air_temperature'], '.', color='tab:cyan', label='air_temperature')\n\n# plot test sample\ndates = matplotlib.dates.date2num(test_sample['timestamp'])\nax2.plot_date(dates, test_sample['meter_reading'], '*', label='test', alpha=0.8)\nax.plot_date(dates, test_sample['air_temperature'], '.', color='tab:cyan', label='air_temperature')\n\nax.set_ylabel('air_temperature'); ax2.set_ylabel('meter_reading')\nax.legend(loc='upper left'); ax2.legend(loc='upper right')\n\ndel train_sample; del test_sample; del dates","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# some feature stats\ntrain.describe()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# the counts above expose the missing data (Should we drop or refill the missing data?)\nprint(\"Ratio of available data (not NAN's):\")\ndata_ratios = train.count()/len(train)\ndata_ratios","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Is the same happening in the test set? Yes\nprint(\"Ratio of available data (not NAN's):\")\ntest.count()/len(test)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# we can refill with averages\ntrain.loc[:, data_ratios < 1.0].mean()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Preprocess data"},{"metadata":{"trusted":true},"cell_type":"code","source":"class ASHRAE3Preprocessor(object):\n    @classmethod\n    def fit(cls, df, data_ratios=data_ratios):\n        cls.avgs = df.loc[:,data_ratios < 1.0].mean()\n        cls.pu_le = LabelEncoder()\n        cls.pu_le.fit(df[\"primary_use\"])\n\n    @classmethod\n    def transform(cls, df):\n        df = df.fillna(cls.avgs) # refill NAN with averages\n        df['primary_use'] = np.uint8(cls.pu_le.transform(df['primary_use']))  # encode labels\n\n        # expand datetime into its components\n        df['hour'] = np.uint8(df['timestamp'].dt.hour)\n        df['day'] = np.uint8(df['timestamp'].dt.day)\n        df['weekday'] = np.uint8(df['timestamp'].dt.weekday)\n        df['month'] = np.uint8(df['timestamp'].dt.month)\n        df['year'] = np.uint8(df['timestamp'].dt.year-2000)\n        \n        # parse and cast columns to a smaller type\n        df.rename(columns={\"square_feet\": \"log_square_feet\"}, inplace=True)\n        df['log_square_feet'] = np.float16(np.log(df['log_square_feet']))\n        df['year_built'] = np.uint8(df['year_built']-1900)\n        df['floor_count'] = np.uint8(df['floor_count'])\n        \n        # remove redundant columns\n        for col in df.columns:\n            if col in ['timestamp', 'row_id']:\n                del df[col]\n    \n        # extract target column\n        if 'meter_reading' in df.columns:\n            df['meter_reading'] = np.log1p(df['meter_reading']).astype(np.float32) # comp metric uses log errors\n\n        return df\n        \nASHRAE3Preprocessor.fit(train)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train = ASHRAE3Preprocessor.transform(train)\ntrain.sample(7)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"train.dtypes","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"Percentage of missing values in the train dataset\")\ntrain.isna().sum()/len(train)*100","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def outlier_function(df, col_name):\n    ''' this function detects first and third quartile and interquartile range for a given column of a dataframe\n    then calculates upper and lower limits to determine outliers conservatively\n    returns the number of lower and uper limit and number of outliers respectively\n    '''\n    first_quartile = np.percentile(\n        np.array(df[col_name].tolist()), 25)\n    third_quartile = np.percentile(\n        np.array(df[col_name].tolist()), 75)\n    IQR = third_quartile - first_quartile\n                      \n    upper_limit = third_quartile+(3*IQR)\n    lower_limit = first_quartile-(3*IQR)\n    outlier_count = 0\n                      \n    for value in df[col_name].tolist():\n        if (value < lower_limit) | (value > upper_limit):\n            outlier_count +=1\n    return lower_limit, upper_limit, outlier_count","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"print(\"{} percent of {} are outliers.\"\n      .format((\n              (100 * outlier_function(train, 'meter_reading')[2])\n               / len(train['meter_reading'])),\n              'meter_reading'))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# Distribution of the meter reading in meters without zeros\nplt.figure(figsize=(12,10))\n\n#list of different meters\nmeters = sorted(train['meter'].unique().tolist()) # [0, 1, 2, 3]\n\n# plot meter_reading distribution for each meter\nfor meter_type in meters:\n    subset = train[train['meter'] == meter_type]\n    sns.kdeplot(np.log1p(subset[\"meter_reading\"]), \n                label=meter_type, linewidth=2)\n\n# set title, legends and labels\nplt.ylabel(\"Density\")\nplt.xlabel(\"Meter_reading\")\nplt.legend(['electricity', 'chilled water', 'steam', 'hot water'])\nplt.title(\"Density of Logartihm(Meter Reading + 1) Among Different Meters\", size=14)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Feature ranked correlation"},{"metadata":{"trusted":true,"_kg_hide-output":false,"_kg_hide-input":true},"cell_type":"code","source":"\nfig, ax = plt.subplots(figsize=(16,8))\n# use a ranked correlation to catch nonlinearities\ncorr = train[[col for col in train.columns if col != 'year']].sample(100100).corr(method='spearman')\n_ = sns.heatmap(corr, annot=True,\n                xticklabels=corr.columns.values,\n                yticklabels=corr.columns.values)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = train[['meter', 'building_id', 'primary_use', 'month', 'day','air_temperature','log_square_feet','floor_count', 'precip_depth_1_hr', 'cloud_coverage']]\ny = train['meter_reading']","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection  import train_test_split\nx_train, x_val, y_train, y_val = train_test_split(X,y, test_size = 0.2, random_state= 45)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn import preprocessing\nfrom keras import backend as K\nfrom keras.callbacks import ModelCheckpoint, EarlyStopping\nfrom keras.layers import Dense, LSTM, GRU, Dropout, BatchNormalization\nfrom keras.models import Sequential\nfrom keras.optimizers import RMSprop,Adam\nfrom keras import regularizers","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def root_mean_squared_error(y_true, y_pred):\n  return K.sqrt(K.mean(K.square(y_pred - y_true)))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def make_model(input_dim=10,metrics=root_mean_squared_error,loss='mse', optimizer=\"rmsprop\",drop_rate=0.5):\n\n  model = Sequential()\n  model.add(LSTM(128,return_sequences=True, input_shape=(None,input_dim)))\n  model.add(Dropout(drop_rate))\n  model.add(BatchNormalization())\n  model.add(LSTM(128,return_sequences=False))\n  model.add(BatchNormalization())\n  model.add(Dropout(drop_rate))\n  model.add(Dense(1))\n  model.compile(optimizer=optimizer, loss=loss, metrics=[metrics])\n  \n  return model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def run_model(model,x_train,y_train,epochs=50,batch_size=500,verbose=1,validation_data=(x_val,y_val),callbacks =None):\n  x_train = x_train.values[:]\n  x_train= x_train.reshape((x_train.shape[0],1,x_train.shape[-1]))\n  y_train = np.log1p(y_train)\n  if validation_data != None:\n    x_val = validation_data[0].values[:]\n    x_val = x_val.reshape((x_val.shape[0],1,x_val.shape[-1]))\n    y_val = np.log1p(validation_data[-1])\n      \n  return model.fit(x_train,y_train,epochs=epochs,batch_size=batch_size,verbose=verbose,validation_data=(x_val,y_val),callbacks=callbacks)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"em = EarlyStopping(monitor='val_root_mean_squared_error', min_delta=0.0001, patience=5, verbose=True, mode='auto')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model = make_model(input_dim=x_train.shape[-1],drop_rate=0.2)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"model.summary()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"history = run_model(model,x_train,y_train,epochs=1,batch_size=500,verbose=1,validation_data=(x_val,y_val), callbacks =[em]) ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"loss = history.history\nloss.keys()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"#rmse loss\nrmse_loss_train = loss['root_mean_squared_error']\nrmse_loss_val = loss['val_root_mean_squared_error']\nepochs_stops = es.stopped_epoch +1 # epochs number from early stopping\nepochs = range(1,epochs_stops + 1)  #len(loss_train)\nplt.figure(figsize=(12,6))\nplt.plot(epochs,rmse_loss_train,'r', label='RMSE train loss')\nplt.plot(epochs,rmse_loss_val,'b',label='RMSE val loss')\nplt.title(' root mean square error loss')\nplt.legend()\nplt.show()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Check prediction"},{"metadata":{"trusted":true,"_kg_hide-input":true,"_kg_hide-output":true},"cell_type":"code","source":"# load and pre-process test data\ntest = ASHRAE3Preprocessor.transform(test)\ntest.sample(7)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"# Test Inference and Submission"},{"metadata":{"trusted":true},"cell_type":"code","source":"# split test data into batches\nset_size = len(test)\niterations = 50\nbatch_size = set_size // iterations\n\nprint(set_size, iterations, batch_size)\nassert set_size == iterations * batch_size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submit = pd.read_csv('/kaggle/input/ashrae-energy-prediction/sample_submission.csv')","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":" x_test = test[['meter', 'building_id', 'primary_use', 'month', 'day','air_temperature','log_square_feet','floor_count', 'precip_depth_1_hr', 'cloud_coverage']]\n x_test = x_test.values[:]\n x_test = x_test.reshape((x_test.shape[0],1,x_test.shape[-1]))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"prediction = model.predict(x_test)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Save submission"}],"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.5"}},"nbformat":4,"nbformat_minor":1}