{"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_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"**Note**: all my cell outputs were gone when I uploaded my notebook to Kaggle...","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\nimport time\n\nimport numpy as np # linear algebra\nfrom scipy.stats import uniform\nnp.random.seed(42)\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nimport geopandas as gpd\nimport seaborn as sns\n\nfrom sklearn.utils.fixes import loguniform\nfrom sklearn.compose import ColumnTransformer\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.impute import SimpleImputer, KNNImputer\nfrom sklearn.decomposition import PCA\nfrom sklearn.feature_selection import SelectFromModel\nfrom feature_engine.creation import CyclicalTransformer\n\nfrom sklearn.pipeline import make_pipeline\nfrom sklearn.linear_model import \\\n        Ridge, Lasso, SGDRegressor\nfrom sklearn.tree import DecisionTreeRegressor\nfrom sklearn.svm import SVR\n\nimport lightgbm as lgb\nfrom lightgbm import LGBMRegressor\n\nfrom sklearn.model_selection import cross_val_score, RandomizedSearchCV, \\\n        PredefinedSplit\nfrom sklearn.metrics import mean_squared_error","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# data_dir = Path(\"/kaggle/input/eurecom-aml-2022-challenge-1\")\ndata_dir = Path(\"data/\")\n\ndf = pd.read_csv(data_dir / \"public/train.csv\", index_col=0)\nfeats_test = pd.read_csv(data_dir / \"public/test_feat.csv\", index_col=0)\n\nfeats = df.iloc[:, :-1]\nlabels = df.iloc[:, -1]\n\nfeats.shape, labels.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data exploration","metadata":{}},{"cell_type":"code","source":"# Find missing vals in both train & test set\nna_cols = set()\nfor df_tmp in [feats, feats_test]:\n    na_counts = df_tmp.isna().sum()\n    na_cols.update(list(na_counts[na_counts > 0].index))\nna_cols = sorted(list(na_cols))\n \nmissing_cnts = []\nfor df_tmp in [feats, feats_test]:\n    na_counts = df_tmp.isna().sum()\n    missing_cnts.append(na_counts[na_cols])\npd.concat(missing_cnts,axis=1).rename(columns={0: \"train\", 1: \"test\"}).T","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Double check missing values in gfs_soil_temperature\nnum_soil_temp_na = (feats[\"gfs_soil_temperature_available\"] == 0).sum()\nprint(f\"Number of missing vals in gfs_soil_temp (train): {num_soil_temp_na} ({num_soil_temp_na / len(feats)})\")\nnum_soil_temp_na_test = (feats_test[\"gfs_soil_temperature_available\"] == 0).sum()\nprint(f\"Number of missing vals in gfs_soil_temp (test): {num_soil_temp_na_test} ({num_soil_temp_na_test / len(feats)})\")\nfeats[\"gfs_soil_temperature\"].value_counts()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feats = feats.reindex(sorted(feats.columns), axis=1)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Look for correlated features \nplt.figure(figsize=(20,10))\nsns.heatmap(feats.drop([\n        \"gfs_soil_temperature_available\",\n        \"cmc_available\",\n        \"gfs_available\",\n        \"wrf_available\"\n    ], axis=1).corr().abs())","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"corr_temp_feats = [\n    \"gfs_temperature_25000\",\n    \"gfs_temperature_30000\",\n    \"gfs_temperature_35000\",\n    \"gfs_temperature_40000\",\n    \"gfs_temperature_45000\",\n    \"gfs_temperature_50000\",\n    \"gfs_temperature_55000\",\n    \"gfs_temperature_60000\",\n    \"gfs_temperature_65000\",\n    \"gfs_temperature_70000\",\n    \"gfs_temperature_75000\",\n    \"gfs_temperature_80000\",\n    \"gfs_temperature_85000\",\n    \"gfs_temperature_90000\",\n    \"gfs_temperature_92500\",\n    \"gfs_temperature_95000\",\n    \"gfs_temperature_97500\",\n]\ncorr_cmc_feats = [\n 'cmc_0_1_65_0',\n 'cmc_0_1_65_0_next',\n 'cmc_0_1_66_0',\n 'cmc_0_1_66_0_next',\n 'cmc_0_1_67_0',\n 'cmc_0_1_67_0_next',\n 'cmc_0_1_68_0',\n 'cmc_0_1_68_0_next',\n 'cmc_0_3_0_0',\n 'cmc_0_3_0_0_next',\n 'cmc_0_3_1_0',\n 'cmc_0_3_5_1000',\n 'wrf_t2_interpolated',\n 'wrf_t2_next',\n]\ncorr_temp_mat = feats[corr_temp_feats].corr().abs()\ncorr_cmc_mat = feats[corr_cmc_feats].corr().abs()","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, axs = plt.subplots(ncols=2, figsize=(30, 10))\n\nsns.heatmap(corr_cmc_mat, vmin=0, annot=True, ax=axs[0])\nsns.heatmap(corr_temp_mat, vmin=0, annot=True, ax=axs[1])","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"skewed_cols = ['cmc_0_1_67_0','cmc_0_1_68_0', 'gfs_clouds_sea', \n               'wrf_rain', 'wrf_snow', 'wrf_graupel', \n               'cmc_0_1_66_0_grad', 'cmc_0_1_67_0_grad', 'cmc_0_1_67_0_next', \n               'cmc_0_1_68_0_grad', 'cmc_0_1_68_0_next']\nfeats[skewed_cols[:10]].hist(bins=50, figsize=(25,10), layout=(2,5))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessing","metadata":{}},{"cell_type":"code","source":"def cyclical_transform(var, max_value):\n    var_sin = np.sin(var * (2. * np.pi / max_value))\n    var_cos = np.cos(var * (2. * np.pi / max_value))\n    return var_sin, var_cos\n\ndef preprocess(feats, labels=None):\n    feats_new = feats.copy()\n    if labels is not None:\n        labels_new = labels.copy()\n    # -9999 were used as NA for this feature\n    # Set these missing values to actual nan in order to impute\n    feats_new.loc[feats_new[\"gfs_soil_temperature_available\"] == 0.0, \"gfs_soil_temperature\"] = np.nan\n    \n    # Drop rows with missing CMC or WRF measurements\n    for col in [\"cmc_available\", \"wrf_available\"]:\n        dropped_idxs = (feats_new[feats_new[col] == 0]).index\n        feats_new.drop(dropped_idxs, axis=0, inplace=True)\n        if labels is not None:\n            labels_new.drop(dropped_idxs, axis=0, inplace=True)\n             \n    # Drop some columns\n    feats_new.drop([\n        # variables that indicates missing values for other feats\n        \"gfs_soil_temperature_available\",\n        \"cmc_available\",\n        \"gfs_available\",\n        \"wrf_available\",\n        # These have ~1.0 correlation coefs\n        'cmc_0_1_65_0_next',\n        'cmc_0_1_66_0_next',\n        'cmc_0_1_67_0_next',\n        'cmc_0_1_68_0_next',\n        'cmc_0_3_0_0_next',\n        'cmc_0_3_5_1000',\n        'wrf_t2_next',\n    ], axis=1, inplace=True)\n    \n    # feats_new.drop(skewed_cols, axis=1, inplace=True)\n    \n    # Encode datetime using cyclical transformations\n    dts = pd.to_datetime(feats_new[\"fact_time\"], unit='s')\n    feats_new[\"sin_hour\"], feats_new[\"cos_hour\"] = \\\n            cyclical_transform(dts.dt.hour, 23)\n    feats_new[\"sin_day\"], feats_new[\"cos_day\"] = \\\n            cyclical_transform(dts.dt.day, dts.dt.daysinmonth)\n    feats_new[\"sin_month\"], feats_new[\"cos_month\"] = \\\n            cyclical_transform(dts.dt.month, 12)\n    feats_new.drop(\"fact_time\", axis=1, inplace=True)\n    \n    if labels is not None:\n        return feats_new, labels_new\n    else:\n        return feats_new\n\nX, y = preprocess(feats, labels)\nN = len(X)\n\n# Train/val split\nidxs = np.arange(N)\nnp.random.shuffle(idxs)\nn_trains = int(0.7 * N)\nX_train = X.iloc[idxs[:n_trains]]\ny_train = y.iloc[idxs[:n_trains]]\nX_val = X.iloc[idxs[n_trains:]]\ny_val = y.iloc[idxs[n_trains:]]\nX.shape, y.shape","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Modeling","metadata":{}},{"cell_type":"code","source":"num_iters = 1500\n\ndef build_model():\n    model = make_pipeline(\n        LGBMRegressor(boosting_type='goss', \n                      num_iterations=num_iters,\n                      verbose=-1,\n                      # verbose_eval =False\n                      learning_rate=0.2,\n                      reg_alpha=0.126,      \n                      num_leaves=200,\n                      max_depth=9\n        )\n    )\n    return model\n\nmodel = build_model()\ncallbacks = [lgb.early_stopping(int(0.1 * num_iters), verbose=0), \n             lgb.log_evaluation(period=10)]\nmodel.fit(X_train, y_train,\n          lgbmregressor__eval_set=[(X_val, y_val)],\n          lgbmregressor__eval_metric=\"rmse\",\n          lgbmregressor__callbacks=callbacks)\nprint(f\"Train RMSE: {mean_squared_error(y_train, model.predict(X_train), squared=False)}\")\nprint(f\"Val RMSE: {mean_squared_error(y_val, model.predict(X_val), squared=False)}\")","metadata":{"tags":[]},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_fold = np.zeros(N)\ntest_fold[idxs[:n_trains]] = -1\nps = PredefinedSplit(test_fold=test_fold)\n\nparam_dist = {\n    \"lgbmregressor__learning_rate\": uniform(0.1, 0.5),\n    \"lgbmregressor__reg_alpha\": loguniform(1e-2, 10),\n    \"lgbmregressor__max_depth\": range(2, 10),\n    \"lgbmregressor__num_leaves\": range(200, 500, 20)\n}\n\nn_iter_search = 20\nrandom_search = RandomizedSearchCV(\n    model,\n    param_distributions=param_dist,\n    n_iter=n_iter_search,\n    cv=ps,\n    refit=True,\n    scoring=\"neg_root_mean_squared_error\",\n    verbose=10\n)\nrandom_search.fit(X, y,\n                  lgbmregressor__eval_set=[(X_val, y_val)],\n                  lgbmregressor__eval_metric=\"rmse\",\n                  lgbmregressor__callbacks=callbacks);","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"best_model = random_search.best_estimator_","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"random_search.best_params_","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Prediction","metadata":{}},{"cell_type":"code","source":"X_test = preprocess(feats_test)\nprint(X_test.shape)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_test = best_model.predict(X_test)\ndf_pred = pd.DataFrame(data={\"fact_temperature\": y_test}, index=feats_test.index)\ndf_pred.to_csv(\"submission.csv\")","metadata":{},"execution_count":null,"outputs":[]}]}