{"cells":[{"metadata":{},"cell_type":"markdown","source":"Hi guys!\nIn this notebook I want to share with you my baseline approaches of handling missing values in time-dependent weather data, in particular:\n* finding consecutive missing blocks and their position/length\n* correspondent (interactive) visualization example with `cufflinks` library\n* creation of simple feature-lightgbm-based imputer to handle longer missing chunks\n* **have not yet figured out what to write next :) **"},{"metadata":{},"cell_type":"markdown","source":"### I had to install latest `cufflinks` version due to internal `plotly` errors in original docker image"},{"metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"!pip install -U cufflinks","execution_count":null,"outputs":[]},{"metadata":{"_uuid":"d629ff2d2480ee46fbb7e2d37f6b5fab8052498a","_cell_guid":"79c7e3d0-c299-4dcb-8224-4455121ee9b0","trusted":true},"cell_type":"code","source":"import pandas as pd\nimport numpy as np\nimport gc\nimport cufflinks as cf\ncf.go_offline(connected=False)  # to make it works without plotly account\nfrom os.path import join as pjoin","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Load weather data"},{"metadata":{"trusted":true},"cell_type":"code","source":"RAW_DATA_DIR = '/kaggle/input/ashrae-energy-prediction/'\n\n# load and concatenate weather data\nweather_dtypes = {\n    'site_id': np.uint8,\n    'air_temperature': np.float32,\n    'cloud_coverage': np.float32,\n    'dew_temperature': np.float32,\n    'precip_depth_1_hr': np.float32,\n    'sea_level_pressure': np.float32,\n    'wind_direction': np.float32,\n    'wind_speed': np.float32,\n}\n\nweather_train = pd.read_csv(\n    pjoin(RAW_DATA_DIR, 'weather_train.csv'),\n    dtype=weather_dtypes,\n    parse_dates=['timestamp']\n)\nweather_test = pd.read_csv(\n    pjoin(RAW_DATA_DIR, 'weather_test.csv'),\n    dtype=weather_dtypes,\n    parse_dates=['timestamp']\n)\n\nweather = pd.concat(\n    [\n        weather_train,\n        weather_test\n    ],\n    ignore_index=True\n)\n\nunique_site_ids = sorted(np.unique(weather['site_id']))\nweather = weather.set_index(['site_id', 'timestamp'], drop=False).sort_index()\n\n# construct full index w/o missing dates\nfull_index = pd.MultiIndex.from_product(\n    [\n        unique_site_ids, \n        pd.date_range(start='2016-01-01 00:00:00', end='2018-12-31 23:00:00', freq='H')\n    ]\n)\n\nprint(f'init shape: {weather.shape}')\nweather = weather.reindex(full_index)\nprint(f'full shape: {weather.shape}')\n\nweather['site_id'] = weather.index.get_level_values(0).astype(np.uint8)\nweather['timestamp'] = weather.index.get_level_values(1)\n\n# drop redundant dfs\ndel weather_train, weather_test\ngc.collect()\n\nprint(weather.dtypes)\n\n# check missing values\nprint(weather.isnull().sum())\n\n# check data sample\nweather.head()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### define supportives to work with NaNs"},{"metadata":{"trusted":true,"_kg_hide-input":false,"_kg_hide-output":false},"cell_type":"code","source":"def get_nan_sequences(series: pd.Series, thld_nan: int = 2):\n    \"\"\"\n    Given sequence with missing data, builds joint index\n    from consecutive NaN blocks\n    1) of len  < thld_nan\n    2) of len >= thld_nan\n    and returns them as 1-D np.arrays\n\n    thld_nan >= 2\n    solution is based on:\n    https://stackoverflow.com/questions/42078259/indexing-a-numpy-array-using-a-numpy-array-of-slices\n    \"\"\"\n    b = series.values\n\n    idx0 = np.flatnonzero(np.r_[True, np.diff(np.isnan(b)) != 0, True])\n    count = np.diff(idx0)\n    idx = idx0[:-1]\n    # >=\n    valid_mask_gte = (count >= thld_nan) & np.isnan(b[idx])\n    out_idx = idx[valid_mask_gte]\n    out_count = count[valid_mask_gte]\n\n    if len(out_idx) == 0:\n        out_gte = np.empty(shape=0)\n    else:\n        out_gte = np.hstack([\n            np.array(range(series, series + n))\n            for (series, n) in zip(out_idx, out_count)\n        ])\n\n    # <\n    valid_mask_lt = (count < thld_nan) & np.isnan(b[idx])\n    out_idx = idx[valid_mask_lt]\n    out_count = count[valid_mask_lt]\n\n    if len(out_idx) == 0:\n        out_lt = np.empty(shape=0)\n    else:\n        out_lt = np.hstack([\n            np.array(range(st, st + n))\n            for (st, n) in zip(out_idx, out_count)\n        ])\n    # check if gte + lt = all NaNs\n    assert len(out_gte) + len(out_lt) == series.isnull().sum(), 'incorrect calculations'\n\n    return out_lt, out_gte\n\n\n# check distribution of consecutive NA parts in data\ndef plot_series_and_consequtive_nans(\n    df: pd.DataFrame, \n    column: str, \n    site_id: int, \n    clip: int = 24,\n    index_slice: slice = None\n):\n    \"\"\"\n    Estimates consecutive NA blocks and plots interactive timeseries with missing data\n    If slice is passed - perform that steps for selected data slice only\n    clips upper block length at 24 (hours)\n    \"\"\"\n    series = weather.loc[site_id][c].copy()\n    if index_slice:\n        series = series.loc[index_slice]\n    \n    # define consecutive nan intervals\n    nan_intervals = series.isnull().astype(int).groupby(\n        series.notnull().astype(int).cumsum()\n    ).sum().clip(0, 24).value_counts()\n    \n    nan_intervals = nan_intervals[nan_intervals.index != 0]\n    \n    nan_intervals.iplot(\n        kind='bar', \n        dimensions=(240*3, 240), \n        title=f'consecutive NaNs in site \"{site_id}\" for column \"{c}\": {nan_intervals.sum()}',\n        xTitle='block length, n points'\n    )\n    \n    # to show missing values as simple interpolations\n    interpolated = series.interpolate()\n    to_plot = pd.DataFrame({c: series, 'missing': interpolated})\n    to_plot.loc[~to_plot[c].isnull(), 'missing'] = np.nan\n    \n    to_plot.iplot(\n        dimensions=(240*3, 320),\n        title=f'site \"{site_id}\", timeseries for column \"{c}\"',\n        xTitle='timestamp',\n        yTitle=f'{c}'\n    )","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# let's see what percentage of missing values do we have in full 2016-2018 weather range\nnulls_by_site_id = (weather.groupby(level=[0]).apply(\n    lambda x: x.isnull().sum()) / len(weather))\nnulls_by_site_id = nulls_by_site_id.loc[:, nulls_by_site_id.any()]\nnulls_by_site_id.index.name = 'site_id'\nnulls_by_site_id.style.format(\"{:.2%}\").highlight_max(axis=0).highlight_min(axis=0, color='#11ff00')\n# btw, we see exact duplication of site_id pairs:\n# (0, 8) and (7, 11)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# let's explore some missing data\nc = 'air_temperature'\nsite_id = 7\nst  = '2017-10-01 00:00:00'\nend = '2018-04-01 00:00:00'\nindex_slice = slice(st, end)\n\nplot_series_and_consequtive_nans(weather, column=c, site_id=site_id, index_slice=index_slice)","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"As we can see, there are **two different types** of missing blocks:\n* single (or near-single) gaps of len 1-2\n* much longer gaps **that is poorly restored** by simple interpolation techniques (see blue lines)\n\nLet's treat them accordingly:\n* For simple missing points (even categorical) simple interpolation might be enough\n<br>this [page](https://docs.scipy.org/doc/scipy/reference/tutorial/interpolate.html) can guide you to different interpolation techniques under the hood of `pd.Series.interpolate`\n* For longer sequences feature-based imputer might come in handy\n\n"},{"metadata":{"trusted":true,"_kg_hide-output":true},"cell_type":"code","source":"# let's fill shorter blocks (of len 1-2) with handy interpolation\n# and then train imputer model to fit longer NaN sequences\n\nthld = 3\nindexes_gte = []\ncols_to_fill = [\n    'air_temperature',\n    'dew_temperature',\n    'precip_depth_1_hr',\n    'cloud_coverage',\n]\n\nint_cols = [\n    # it's ordered but lives in integer scale\n    'cloud_coverage'\n]\n\nnans_total = weather[cols_to_fill].isnull().sum().sum()\nnans_filled = 0\nfor col in cols_to_fill:\n    print(f'filling short NaN series in col \"{col}\"')\n    dtype = np.int8 if col in int_cols else np.float32\n    for sid in sorted(weather.site_id.unique()):\n        print(f'\\tfor site_id: \"{sid}\"')\n        s = weather.loc[sid, col].copy()\n        idx_lt, idx_gte = get_nan_sequences(s, thld)\n        interpolation = s.interpolate()\n        nans_before = weather.loc[sid, col].isnull().sum()\n        print(f'\\t\\tnans before: {nans_before}')\n        weather.loc[sid, col].iloc[idx_lt] = interpolation.iloc[idx_lt].values.astype(dtype)\n        nans_after = weather.loc[sid, col].isnull().sum()\n        print(f'\\t\\tnans  after: {nans_before}')\n        nans_filled += (nans_before - nans_after)\n\nprint(f'Nans filled: {nans_filled}/{nans_total}:   {np.round(nans_filled/nans_total*100, 2)}%')","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"### Let's create simple feature-based imputer for temperatures"},{"metadata":{"trusted":true},"cell_type":"code","source":"# define simple imputer for longer missing sequences\nimport lightgbm as lgb\nimport os\n\n\ndef nan_imputer(data: pd.DataFrame, tcol: str, window: int = 24):\n    \n    df = data.copy()\n    \n    reg = lgb.LGBMRegressor(\n        learning_rate=0.05,\n        objective='mae',\n        n_estimators=350,\n        num_threads=os.cpu_count(),\n        num_leaves=31,\n        max_depth=8,\n        subsample=0.8,\n        min_child_samples=50,\n        random_state=42,\n    )\n\n    init_cols = df.columns.tolist()\n\n    dtime_col = 'timestamp'\n    df['year'] = df[dtime_col].dt.year.astype(np.int16)\n    df['hour'] = df[dtime_col].dt.hour.astype(np.uint8)\n    df['month'] = df[dtime_col].dt.month.astype(np.uint8) - 1\n    df['weekday'] = df[dtime_col].dt.weekday.astype(np.uint8)\n    df['dayofyear'] = df[dtime_col].dt.dayofyear.astype(np.uint16) - 1\n    df['weekofyear'] = df[dtime_col].dt.weekofyear.astype(np.uint8) - 1\n    df['quarter'] = df[dtime_col].dt.quarter.astype(np.uint8) - 1\n    df['monthday'] = df[dtime_col].dt.day.astype(np.uint8) - 1\n\n    df['rolling_back'] = df.groupby(by='site_id')[tcol]\\\n        .rolling(window=window, min_periods=1).mean().interpolate().values\n\n    # reversed rolling\n    df['rolling_forw'] = df.iloc[::-1].groupby(by='site_id')[tcol]\\\n        .rolling(window=window, min_periods=1).mean().interpolate().values\n\n    # rolling mean for same hour of the day\n    df['rolling_back_h'] = df.groupby(by=['site_id', 'hour'])[tcol]\\\n        .rolling(window=3, min_periods=1).mean().interpolate().values\n\n    df['rolling_back_h_f'] = df.iloc[::-1].groupby(by=['site_id', 'hour'])[tcol]\\\n        .rolling(window=3, min_periods=1).mean().interpolate().values\n    \n#     sampler = np.random.RandomState(42)\n#     df['interpolation'] = df[tcol].interpolate() * (1 + sampler.randn(len(df)) * 0.25)\n#     df.iloc[\n#         np.random.choice(len(df), int(len(df)*3/5)),\n#         df.columns.tolist().index('interpolation')\n#     ] = np.nan\n\n    tr_idx, val_idx = ~df[tcol].isnull(), df[tcol].isnull()\n\n    features = [\n        'site_id', 'hour', 'month', 'dayofyear', 'weekofyear', 'year',\n        'rolling_back',\n        'rolling_forw',\n        'rolling_back_h',\n        'rolling_back_h_f',\n#         'interpolation'\n    ]\n    \n    print(f'training model for col \"{tcol}\"...')\n    reg.fit(\n        X=df.loc[tr_idx, features], \n        y=df.loc[tr_idx, tcol],\n        categorical_feature=['site_id', 'year'],\n    )\n\n    df[f'{tcol}_restored'] = np.nan\n    df.loc[val_idx, f'{tcol}_restored'] = reg.predict(df.loc[val_idx, features])\n    df.loc[val_idx, f'{tcol}'] = df.loc[val_idx, f'{tcol}_restored'].values\n    \n    # add simple rolling mean for comparison\n    df[f'{tcol}_rolling_mean'] = df.groupby(by='site_id')[tcol]\\\n    .rolling(window=24*3, min_periods=1).mean().values\n    \n    # check from what features our imputer learned the most\n    lgb.plot_importance(reg)\n    \n    return df.loc[:, init_cols + [f'{tcol}_restored', f'{tcol}_rolling_mean']]\n\ntcols = ['air_temperature', 'dew_temperature']\nrestored = weather\nfor tcol in tcols:\n    restored = nan_imputer(data=restored, tcol=tcol, window=24)","execution_count":null,"outputs":[]},{"metadata":{"trusted":true},"cell_type":"code","source":"# check the sample of imputation results\ntcol = 'air_temperature'\n\nst =  (11, '2017-07-01 00:00:00')\nend = (11, '2018-07-01 00:00:00')\n\nrestored.loc[st:end].set_index('timestamp')[[\n    f'{tcol}', \n    f'{tcol}_restored', \n    f'{tcol}_rolling_mean'\n]].iplot()","execution_count":null,"outputs":[]},{"metadata":{},"cell_type":"markdown","source":"At the end of the day, such (or similar) restoration techniques over ENTIRE possible time-grid allow us to **correctly** build some lag-based features, such as rolling means, shifts, etc.\nHowever, their importance for this particular competition remains unclear, we need to deep dive further\n\nHope you enjoyed this kernel and learned some hints & tricks!\n\n**P.s.** Comments, likes, new ideas are highly welcomed!\n<br>Happy kaggling!\n\n---\nCheck my latest notebooks:\n- [Aligning Temperature Timestamp](https://www.kaggle.com/frednavruzov/aligning-temperature-timestamp)\n- [Faster stratified cross-validation](https://www.kaggle.com/frednavruzov/faster-stratified-cross-validation-upd)"}],"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}