{"cells":[{"metadata":{},"cell_type":"markdown","source":"One of the key base models in my top 1% (25th place) solution. Greatly influenced by/borrowed from the following great kernels - please go and upvote the original work of the authors:\n\nhttps://www.kaggle.com/isaienkov/keras-nn-with-embeddings-for-cat-features-1-15  \nhttps://www.kaggle.com/purist1024/ashrae-simple-data-cleanup-lb-1-08-no-leaks"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load in \n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the \"../input/\" directory.\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# Any results you write to the current directory are saved as output.\n\nfrom keras.models import Model, load_model\nfrom keras.layers import Input, Dropout, Dense, Embedding, SpatialDropout1D, concatenate, BatchNormalization, Flatten\nfrom keras.preprocessing.sequence import pad_sequences\nfrom keras.preprocessing import text, sequence\nfrom keras.callbacks import Callback\nfrom keras import backend as K\nfrom keras.models import Model\nfrom keras.losses import mean_squared_error as mse_loss\nimport gc\n\nfrom keras import optimizers\nfrom keras.optimizers import RMSprop, Adam\nfrom keras.callbacks import EarlyStopping, ModelCheckpoint, ReduceLROnPlateau","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def average_imputation(df, column_name):\n    imputation = df.groupby(['timestamp'])[column_name].mean()\n    \n    df.loc[df[column_name].isnull(), column_name] = df[df[column_name].isnull()][[column_name]].apply(lambda x: imputation[df['timestamp'][x.index]].values)\n    del imputation\n    return df\n","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# cleansing routines ==========================================================================\ndef make_is_bad_zero(Xy_subset, min_interval=48, summer_start=3000, summer_end=7500):\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        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\ndef 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\ndef find_bad_sitezero(X):\n    \"\"\"Returns indices of bad rows from the early days of Site 0 (UCF).\"\"\"\n    return X[(X.timestamp < 3378) & (X.site_id == 0) & (X.meter == 0)].index\n\ndef 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# binds all together\ndef find_bad_rows(X, y):\n    return find_bad_zeros(X, y).union(find_bad_sitezero(X)).union(find_bad_building1099(X, y))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"categorical_columns = [\n    \"building_id\", \"meter\", \"site_id\", \"primary_use\", \"had_air_temperature\", \"had_cloud_coverage\",\n    \"had_dew_temperature\", \"had_precip_depth_1_hr\", \"had_sea_level_pressure\", \"had_wind_direction\",\n    \"had_wind_speed\", \"tm_day_of_week\", \"tm_hour_of_day\"\n]","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"pd.set_option(\"max_columns\", 500)\n\ndef input_file(file):\n    path = f\"../input/ashrae-energy-prediction/{file}\"\n    if not os.path.exists(path): return path + \".gz\"\n    return path\n\ndef compress_dataframe(df):\n    result = df.copy()\n    for col in result.columns:\n        col_data = result[col]\n        dn = col_data.dtype.name\n        if dn == \"object\":\n            result[col] = pd.to_numeric(col_data.astype(\"category\").cat.codes, downcast=\"integer\")\n        elif dn == \"bool\":\n            result[col] = col_data.astype(\"int8\")\n        elif dn.startswith(\"int\") or (col_data.round() == col_data).all():\n            result[col] = pd.to_numeric(col_data, downcast=\"integer\")\n        else:\n            result[col] = pd.to_numeric(col_data, downcast='float')\n    return result\n\ndef read_train():\n    df = pd.read_csv(input_file(\"train.csv\"), parse_dates=[\"timestamp\"])\n    df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n    return compress_dataframe(df)\n\ndef read_building_metadata():\n    return compress_dataframe(pd.read_csv(\n        input_file(\"building_metadata.csv\")).fillna(-1)).set_index(\"building_id\")\n\nsite_GMT_offsets = [-5, 0, -7, -5, -8, 0, -5, -5, -5, -6, -7, -5, 0, -6, -5, -5]\n\ndef read_weather_train(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    df = pd.read_csv(input_file(\"weather_train.csv\"), parse_dates=[\"timestamp\"])\n    df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n    if fix_timestamps:\n        GMT_offset_map = {site: offset for site, offset in enumerate(site_GMT_offsets)}\n        df.timestamp = df.timestamp + df.site_id.map(GMT_offset_map)\n    if interpolate_na:\n        site_dfs = []\n        for site_id in df.site_id.unique():\n            # Make sure that we include all possible hours so that we can interpolate evenly\n            site_df = df[df.site_id == site_id].set_index(\"timestamp\").reindex(range(8784))\n            site_df.site_id = site_id\n            for col in [c for c in site_df.columns if c != \"site_id\"]:\n                if add_na_indicators: site_df[f\"had_{col}\"] = ~site_df[col].isna()\n                site_df[col] = site_df[col].interpolate(limit_direction='both', method='linear')\n                # Some sites are completely missing some columns, so use this fallback\n                site_df[col] = site_df[col].fillna(df[col].median())\n            site_dfs.append(site_df)\n        df = pd.concat(site_dfs).reset_index()  # make timestamp back into a regular column\n    elif add_na_indicators:\n        for col in df.columns:\n            if df[col].isna().any(): df[f\"had_{col}\"] = ~df[col].isna()\n    return compress_dataframe(df).set_index([\"site_id\", \"timestamp\"])\n\ndef combined_train_data(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    Xy = compress_dataframe(read_train().join(read_building_metadata(), on=\"building_id\").join(\n        read_weather_train(fix_timestamps, interpolate_na, add_na_indicators),\n        on=[\"site_id\", \"timestamp\"]).fillna(-1))\n    return Xy.drop(columns=[\"meter_reading\"]), Xy.meter_reading\n\ndef _add_time_features(X):\n    return X.assign(tm_day_of_week=((X.timestamp // 24) % 7), tm_hour_of_day=(X.timestamp % 24))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X, y = combined_train_data()\n\nbad_rows = find_bad_rows(X, y)\npd.Series(bad_rows.sort_values()).to_csv(\"rows_to_drop.csv\", header=False, index=False)\n\nX = X.drop(index=bad_rows)\ny = y.reindex_like(X)\n\nbeaufort = [(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)]\nfor item in beaufort:\n    X.loc[(X['wind_speed']>=item[1]) & (X['wind_speed']<item[2]), 'beaufort_scale'] = item[0]\n\n# Additional preprocessing\nX = compress_dataframe(_add_time_features(X))\nX = X.drop(columns=\"timestamp\")  # Raw timestamp doesn't help when prediction\nX.columns = ['building_id', 'meter', 'site_id', 'primary_use', 'square_feet',\n       'year_built', 'floor_count', 'air_temperature', 'cloud_coverage',\n       'dew_temperature', 'precip_depth_1_hr', 'sea_level_pressure',\n       'wind_direction', 'wind_speed', 'had_air_temperature',\n       'had_cloud_coverage', 'had_dew_temperature', 'had_precip_depth_1_hr',\n       'had_sea_level_pressure', 'had_wind_direction', 'had_wind_speed',\n       'beaufort_scale', 'weekday', 'hour']\ny = np.log1p(y)\n\ntrain = X.copy()\ntarget = y\n\ndel X, y\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"categoricals = [\"site_id\", \"building_id\", \"primary_use\", \"hour\", \"weekday\",  \"meter\"]\n\ndrop_cols = [\"sea_level_pressure\",\"wind_speed\",\"wind_direction\",\"had_air_temperature\",\"had_cloud_coverage\",\"had_dew_temperature\",\"had_precip_depth_1_hr\",\"had_sea_level_pressure\",\"had_wind_direction\",\"had_wind_speed\"]\n\nnumericals = [\"square_feet\", \"year_built\", \"air_temperature\", \"cloud_coverage\",\n              \"dew_temperature\", \"precip_depth_1_hr\", \"floor_count\", 'beaufort_scale']\n\nfeat_cols = categoricals + numericals\n\ntrain.drop(drop_cols, axis = 1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def model(dense_dim_1=64, dense_dim_2=32, dense_dim_3=32, dense_dim_4=16, \ndropout1=0.2, dropout2=0.1, dropout3=0.1, dropout4=0.1, lr=0.001):\n\n    #Inputs\n    site_id = Input(shape=[1], name=\"site_id\")\n    building_id = Input(shape=[1], name=\"building_id\")\n    meter = Input(shape=[1], name=\"meter\")\n    primary_use = Input(shape=[1], name=\"primary_use\")\n    square_feet = Input(shape=[1], name=\"square_feet\")\n    year_built = Input(shape=[1], name=\"year_built\")\n    air_temperature = Input(shape=[1], name=\"air_temperature\")\n    cloud_coverage = Input(shape=[1], name=\"cloud_coverage\")\n    dew_temperature = Input(shape=[1], name=\"dew_temperature\")\n    hour = Input(shape=[1], name=\"hour\")\n    precip = Input(shape=[1], name=\"precip_depth_1_hr\")\n    weekday = Input(shape=[1], name=\"weekday\")\n    beaufort_scale = Input(shape=[1], name=\"beaufort_scale\")\n   \n    #Embeddings layers\n    emb_site_id = Embedding(16, 2)(site_id)\n    emb_building_id = Embedding(1449, 6)(building_id)\n    emb_meter = Embedding(4, 2)(meter)\n    emb_primary_use = Embedding(16, 2)(primary_use)\n    emb_hour = Embedding(24, 3)(hour)\n    emb_weekday = Embedding(7, 2)(weekday)\n\n    concat_emb = concatenate([\n           Flatten() (emb_site_id)\n         , Flatten() (emb_building_id)\n         , Flatten() (emb_meter)\n         , Flatten() (emb_primary_use)\n         , Flatten() (emb_hour)\n         , Flatten() (emb_weekday)\n    ])\n    \n    categ = Dropout(dropout1)(Dense(dense_dim_1,activation='relu') (concat_emb))\n    categ = BatchNormalization()(categ)\n    categ = Dropout(dropout2)(Dense(dense_dim_2,activation='relu') (categ))\n    \n    #main layer\n    main_l = concatenate([\n          categ\n        , square_feet\n        , year_built\n        , air_temperature\n        , cloud_coverage\n        , dew_temperature\n        , precip\n        , beaufort_scale\n    ])\n    \n    main_l = Dropout(dropout3)(Dense(dense_dim_3,activation='relu') (main_l))\n    main_l = BatchNormalization()(main_l)\n    main_l = Dropout(dropout4)(Dense(dense_dim_4,activation='relu') (main_l))\n    \n    #output\n    output = Dense(1) (main_l)\n\n    model = Model([ site_id,\n                    building_id, \n                    meter, \n                    primary_use, \n                    square_feet, \n                    year_built, \n                    air_temperature,\n                    cloud_coverage,\n                    dew_temperature, \n                    hour,\n                    weekday, \n                    precip,\n                    beaufort_scale], output)\n\n    model.compile(optimizer = Adam(lr=lr),\n                  loss= mse_loss,\n                  metrics=[root_mean_squared_error])\n    return model\n\ndef root_mean_squared_error(y_true, y_pred):\n    return K.sqrt(K.mean(K.square(y_pred - y_true), axis=0))","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def get_keras_data(df, num_cols, cat_cols):\n    cols = num_cols + cat_cols\n    X = {col: np.array(df[col]) for col in cols}\n    return X\n\ndef train_model(keras_model, X_t, y_train, batch_size, epochs, X_v, y_valid, fold, patience=3):\n    early_stopping = EarlyStopping(patience=patience, verbose=1)\n    rlrop = ReduceLROnPlateau(monitor='val_loss', mode='min', patience=3, factor=0.5, min_lr=1e-6, verbose=1)\n    model_checkpoint = ModelCheckpoint(\"model_\" + str(fold) + \".hdf5\",\n                                       save_best_only=True, verbose=1, monitor='val_root_mean_squared_error', mode='min')\n\n    hist = keras_model.fit(X_t, y_train, batch_size=batch_size, epochs=epochs,\n                            validation_data=(X_v, y_valid), verbose=1,\n                            callbacks=[early_stopping, rlrop, model_checkpoint])\n\n    keras_model = load_model(\"model_\" + str(fold) + \".hdf5\", custom_objects={'root_mean_squared_error': root_mean_squared_error})\n    \n    return keras_model","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from sklearn.model_selection import KFold, StratifiedKFold\n\noof = np.zeros(len(train))\nbatch_size = 1024\nepochs = 20\nmodels = []\n\nfolds = 4\nseed = 666\n\nkf = StratifiedKFold(n_splits=folds, shuffle=True, random_state=seed)\n\nfor fold_n, (train_index, valid_index) in enumerate(kf.split(train, train['building_id'])):\n    print('Fold:', fold_n)\n    X_train, X_valid = train.iloc[train_index], train.iloc[valid_index]\n    y_train, y_valid = target.iloc[train_index], target.iloc[valid_index]\n    X_t = get_keras_data(X_train, numericals, categoricals)\n    X_v = get_keras_data(X_valid, numericals, categoricals)\n    \n    keras_model = model(dense_dim_1=64, dense_dim_2=32, dense_dim_3=32, dense_dim_4=16, \n                        dropout1=0.2, dropout2=0.1, dropout3=0.1, dropout4=0.1, lr=0.001)\n    mod = train_model(keras_model, X_t, y_train, batch_size, epochs, X_v, y_valid, fold_n, patience=5)\n    models.append(mod)\n    print('*'* 50)\n    ","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"del train, target, X_train, X_valid, y_train, y_valid, X_t, X_v, kf\ngc.collect()","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"def read_test():\n    df = pd.read_csv(input_file(\"test.csv\"), parse_dates=[\"timestamp\"])\n    df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n    return compress_dataframe(df).set_index(\"row_id\")\n\ndef read_weather_test(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    df = pd.read_csv(input_file(\"weather_test.csv\"), parse_dates=[\"timestamp\"])\n    df.timestamp = (df.timestamp - pd.to_datetime(\"2016-01-01\")).dt.total_seconds() // 3600\n    if fix_timestamps:\n        GMT_offset_map = {site: offset for site, offset in enumerate(site_GMT_offsets)}\n        df.timestamp = df.timestamp + df.site_id.map(GMT_offset_map)\n    if interpolate_na:\n        site_dfs = []\n        for site_id in df.site_id.unique():\n            # Make sure that we include all possible hours so that we can interpolate evenly\n            site_df = df[df.site_id == site_id].set_index(\"timestamp\").reindex(range(8784, 26304))\n            site_df.site_id = site_id\n            for col in [c for c in site_df.columns if c != \"site_id\"]:\n                if add_na_indicators: site_df[f\"had_{col}\"] = ~site_df[col].isna()\n                site_df[col] = site_df[col].interpolate(limit_direction='both', method='linear')\n                # Some sites are completely missing some columns, so use this fallback\n                site_df[col] = site_df[col].fillna(df[col].median())\n            site_dfs.append(site_df)\n        df = pd.concat(site_dfs).reset_index()  # make timestamp back into a regular column\n    elif add_na_indicators:\n        for col in df.columns:\n            if df[col].isna().any(): df[f\"had_{col}\"] = ~df[col].isna()\n    return compress_dataframe(df).set_index([\"site_id\", \"timestamp\"])\n\ndef combined_test_data(fix_timestamps=True, interpolate_na=True, add_na_indicators=True):\n    X = compress_dataframe(read_test().join(read_building_metadata(), on=\"building_id\").join(\n        read_weather_test(fix_timestamps, interpolate_na, add_na_indicators),\n        on=[\"site_id\", \"timestamp\"]).fillna(-1))\n    return X","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"X = combined_test_data()\n\nfor item in beaufort:\n    X.loc[(X['wind_speed']>=item[1]) & (X['wind_speed']<item[2]), 'beaufort_scale'] = item[0]\n\n# Additional preprocessing\nX = compress_dataframe(_add_time_features(X))\nX = X.drop(columns=\"timestamp\")  # Raw timestamp doesn't help when prediction\nX.columns = ['building_id', 'meter', 'site_id', 'primary_use', 'square_feet',\n       'year_built', 'floor_count', 'air_temperature', 'cloud_coverage',\n       'dew_temperature', 'precip_depth_1_hr', 'sea_level_pressure',\n       'wind_direction', 'wind_speed', 'had_air_temperature',\n       'had_cloud_coverage', 'had_dew_temperature', 'had_precip_depth_1_hr',\n       'had_sea_level_pressure', 'had_wind_direction', 'had_wind_speed',\n       'beaufort_scale', 'weekday', 'hour']\n\ntest = X.copy()\n\ndel X\ngc.collect()\n\ntest.drop(drop_cols, axis = 1, inplace=True)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"from tqdm import tqdm\ni=0\nres = np.zeros((test.shape[0]),dtype=np.float32)\nstep_size = 50000\nfor j in tqdm(range(int(np.ceil(test.shape[0]/step_size)))):\n    for_prediction = get_keras_data(test.iloc[i:i+step_size], numericals, categoricals)\n    res[i:min(i+step_size,test.shape[0])] = \\\n       np.expm1(sum([model.predict(for_prediction, batch_size=1024)[:,0] for model in models])/folds)\n    i+=step_size","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"submission = pd.read_csv('/kaggle/input/ashrae-energy-prediction/sample_submission.csv')\nsubmission['meter_reading'] = res\nsubmission.loc[submission['meter_reading']<0, 'meter_reading'] = 0\nsubmission.to_csv('submission.csv', index=False)\nsubmission","execution_count":null,"outputs":[]}],"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"pygments_lexer":"ipython3","nbconvert_exporter":"python","version":"3.6.4","file_extension":".py","codemirror_mode":{"name":"ipython","version":3},"name":"python","mimetype":"text/x-python"}},"nbformat":4,"nbformat_minor":1}