{"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":"code","source":"import pandas as pd\nimport numpy as np\nimport seaborn as sns\nimport matplotlib.pyplot as plt\nimport copy\nfrom sklearn.linear_model import Lasso\nfrom sklearn.model_selection import GridSearchCV\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.ensemble import RandomForestRegressor","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:19:22.433716Z","iopub.execute_input":"2022-05-10T21:19:22.434077Z","iopub.status.idle":"2022-05-10T21:19:24.115119Z","shell.execute_reply.started":"2022-05-10T21:19:22.434035Z","shell.execute_reply":"2022-05-10T21:19:24.114231Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = pd.read_csv(\"../input/eurecom-aml-2022-challenge-1/public/train.csv\", low_memory=True).iloc[:,1:]  # drop index column","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:19:24.117276Z","iopub.execute_input":"2022-05-10T21:19:24.117863Z","iopub.status.idle":"2022-05-10T21:21:01.459332Z","shell.execute_reply.started":"2022-05-10T21:19:24.117809Z","shell.execute_reply":"2022-05-10T21:21:01.458016Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.head()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:05:51.446541Z","iopub.execute_input":"2022-05-10T21:05:51.446804Z","iopub.status.idle":"2022-05-10T21:05:51.485424Z","shell.execute_reply.started":"2022-05-10T21:05:51.446774Z","shell.execute_reply":"2022-05-10T21:05:51.484594Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Exploratory Data Analysis","metadata":{}},{"cell_type":"markdown","source":"Let's inspect geographical data\nFirst, we look at temperature in different areas of the world\n","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots()\ntrain_df.plot(kind=\"scatter\", x=\"fact_longitude\", y=\"fact_latitude\", alpha=0.4, c=\"fact_temperature\",\n              cmap=plt.get_cmap(\"jet\"), colorbar=True, figsize=(15,7), ax=ax)\nplt.tight_layout()\nplt.grid()\nplt.xticks(fontsize=20)\nplt.yticks(fontsize=20)\nplt.xlabel(\"Longitude\", fontsize=20)\nplt.ylabel(\"Latitude\", fontsize=20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:05:51.487303Z","iopub.execute_input":"2022-05-10T21:05:51.487553Z","iopub.status.idle":"2022-05-10T21:06:21.891537Z","shell.execute_reply.started":"2022-05-10T21:05:51.487522Z","shell.execute_reply":"2022-05-10T21:06:21.880986Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Now, look at climate pressure in different areas\n","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots()\ntrain_df.plot(kind=\"scatter\", x=\"fact_longitude\", y=\"fact_latitude\", alpha=0.4, c=\"climate_pressure\",\n              cmap=plt.get_cmap(\"jet\"), colorbar=True, figsize=(15,7), ax=ax)\nplt.tight_layout()\nplt.grid()\nplt.xticks(fontsize=20)\nplt.yticks(fontsize=20)\nplt.xlabel(\"Longitude\", fontsize=20)\nplt.ylabel(\"Latitude\", fontsize=20)\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:06:21.893134Z","iopub.execute_input":"2022-05-10T21:06:21.893634Z","iopub.status.idle":"2022-05-10T21:06:51.786345Z","shell.execute_reply.started":"2022-05-10T21:06:21.893599Z","shell.execute_reply":"2022-05-10T21:06:51.785315Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see how training data is coming from all around the world, with some exceptions: center and north South America and center of Africa.","metadata":{}},{"cell_type":"markdown","source":"## Target Variable","metadata":{}},{"cell_type":"code","source":"import seaborn as sns\n\ndef plot_violin_plot(column):\n  sns.set_theme(style=\"whitegrid\")\n  ax = sns.violinplot(y=column)\n\nplt.figure(figsize=(8,6))\nplot_violin_plot(train_df[\"fact_temperature\"])\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:04:02.317642Z","iopub.execute_input":"2022-05-10T20:04:02.318413Z","iopub.status.idle":"2022-05-10T20:04:06.425651Z","shell.execute_reply.started":"2022-05-10T20:04:02.318375Z","shell.execute_reply":"2022-05-10T20:04:06.424911Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"percentile_99_99 = train_df[\"fact_temperature\"].quantile(0.9999)\npercentile_0_01 = train_df[\"fact_temperature\"].quantile(0.0001)\npercentile_99_9 = train_df[\"fact_temperature\"].quantile(0.999)\npercentile_0_1 = train_df[\"fact_temperature\"].quantile(0.001)\npercentile_99 = train_df[\"fact_temperature\"].quantile(0.99)\npercentile_1 = train_df[\"fact_temperature\"].quantile(0.1)\n\n\nprint(f\"Percentile 99.99 %: {percentile_99_99}\")\nprint(f\"Percentile 0.01 %: {percentile_0_01}\")\nprint(f\"Percentile 99.9 %: {percentile_99_9}\")\nprint(f\"Percentile 0.1 %: {percentile_0_1}\")\nprint(f\"Percentile 99 %: {percentile_99}\")\nprint(f\"Percentile 1 %: {percentile_1}\")\n\nprint(f\"Percentile 99.99 % data over: { ( train_df['fact_temperature'] >= percentile_99_99).sum()} \")\nprint(f\"Percentile 0.01 % data under: { ( train_df['fact_temperature'] <= percentile_0_01).sum()} \")\nprint(f\"Percentile 99.9 % data over: { ( train_df['fact_temperature'] >= percentile_99_9).sum()} \")\nprint(f\"Percentile 0.1 % data under: { ( train_df['fact_temperature'] <= percentile_0_1).sum()} \")\nprint(f\"Percentile 99 % data over: { ( train_df['fact_temperature'] >= percentile_99).sum()} \")\nprint(f\"Percentile 1 % data under: { ( train_df['fact_temperature'] <= percentile_1).sum()}\")","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:04:06.426826Z","iopub.execute_input":"2022-05-10T20:04:06.427040Z","iopub.status.idle":"2022-05-10T20:04:06.677875Z","shell.execute_reply.started":"2022-05-10T20:04:06.427014Z","shell.execute_reply":"2022-05-10T20:04:06.677022Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data cleaning","metadata":{}},{"cell_type":"markdown","source":"We check that all features are numerical -> yes, so we don't need encoding","metadata":{}},{"cell_type":"code","source":"train_df.info()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:04:06.679821Z","iopub.execute_input":"2022-05-10T20:04:06.680069Z","iopub.status.idle":"2022-05-10T20:04:06.697952Z","shell.execute_reply.started":"2022-05-10T20:04:06.680030Z","shell.execute_reply":"2022-05-10T20:04:06.697189Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"NaN values before dealing with them: {train_df.isna().sum().sum()}\")","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:04:14.949501Z","iopub.execute_input":"2022-05-10T20:04:14.950210Z","iopub.status.idle":"2022-05-10T20:04:15.401516Z","shell.execute_reply.started":"2022-05-10T20:04:14.950161Z","shell.execute_reply":"2022-05-10T20:04:15.400741Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"All NaN values are in the following records\n","metadata":{}},{"cell_type":"code","source":"print(f\"gfs records not available: {(train_df['gfs_available'] == 0).sum()}\")\nprint(f\"cmc records not available: {(train_df['cmc_available'] == 0).sum()}\")\nprint(f\"wrf records not available: {(train_df['wrf_available'] == 0).sum()}\")","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:04:20.498124Z","iopub.execute_input":"2022-05-10T20:04:20.498905Z","iopub.status.idle":"2022-05-10T20:04:20.522141Z","shell.execute_reply.started":"2022-05-10T20:04:20.498862Z","shell.execute_reply":"2022-05-10T20:04:20.519902Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We create a categorical variable called 'flag', which is equal to 0 if the record does not have any null features, 1 if cmc features are missing, 2 if wrf features are missing and 3 if both cmc and wrf features are missing. \n\nThis variable is created only for visualization purposes.\n","metadata":{}},{"cell_type":"code","source":"nan_points_cmc = (train_df['cmc_available'] == 0).replace({True: 1, False: 0})\nnan_points_wrf = (train_df['wrf_available'] == 0).replace({True: 2, False: 0})\nnan_points_flag = pd.concat([nan_points_cmc, nan_points_wrf]).groupby(level=0).sum()\ntrain_df[\"flag\"] = nan_points_flag\n","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:21:01.462850Z","iopub.execute_input":"2022-05-10T21:21:01.463316Z","iopub.status.idle":"2022-05-10T21:21:03.679983Z","shell.execute_reply.started":"2022-05-10T21:21:01.463231Z","shell.execute_reply":"2022-05-10T21:21:03.678809Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df.plot(kind=\"scatter\", x=\"fact_longitude\", y=\"fact_latitude\", alpha=0.4, c=\"flag\",\n              cmap=plt.get_cmap(\"jet\"), colorbar=True, figsize=(15,7))\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:08:45.964376Z","iopub.execute_input":"2022-05-10T20:08:45.964737Z","iopub.status.idle":"2022-05-10T20:09:21.622581Z","shell.execute_reply.started":"2022-05-10T20:08:45.964683Z","shell.execute_reply":"2022-05-10T20:09:21.621754Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Blue points represent regions with no null records. We can see that half of Australia points are records with NaN values (yellow points). ","metadata":{}},{"cell_type":"code","source":"# drop available data feature and its flag which was useful for visualization only\ntrain_df = train_df.drop(columns=[\"gfs_available\", \"cmc_available\", \"wrf_available\", \"flag\"])","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:21:03.681415Z","iopub.execute_input":"2022-05-10T21:21:03.681673Z","iopub.status.idle":"2022-05-10T21:21:04.594364Z","shell.execute_reply.started":"2022-05-10T21:21:03.681643Z","shell.execute_reply":"2022-05-10T21:21:04.593576Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Dealing with NaN values by deleting them will lead us to a loss of information about half of Oceania. Thus, we will sort by latitude and longitude, then use bfill and ffill approaches, then restore original order by index. With bfill and ffill we mean backward and forward fill. The former uses next valid observation to fill the gap, while the latter propagates last valid observation forward","metadata":{}},{"cell_type":"code","source":"train_df = train_df.sort_values([\"fact_latitude\", \"fact_longitude\"])\ntrain_df = train_df.fillna(method=\"bfill\")\ntrain_df = train_df.fillna(method=\"ffill\")\ntrain_df = train_df.sort_index()\ntrain_df","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:21:04.596228Z","iopub.execute_input":"2022-05-10T21:21:04.596633Z","iopub.status.idle":"2022-05-10T21:21:15.950671Z","shell.execute_reply.started":"2022-05-10T21:21:04.596600Z","shell.execute_reply":"2022-05-10T21:21:15.949839Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"As a result, we won't have any NaN record. ","metadata":{}},{"cell_type":"code","source":"train_df.isna().sum().sum()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:17:05.605397Z","iopub.execute_input":"2022-05-10T20:17:05.605775Z","iopub.status.idle":"2022-05-10T20:17:06.204458Z","shell.execute_reply.started":"2022-05-10T20:17:05.605739Z","shell.execute_reply":"2022-05-10T20:17:06.203810Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We remove the duplicates from our training set","metadata":{}},{"cell_type":"code","source":"train_df.drop_duplicates(inplace=True, ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:21:15.951892Z","iopub.execute_input":"2022-05-10T21:21:15.952273Z","iopub.status.idle":"2022-05-10T21:21:52.760378Z","shell.execute_reply.started":"2022-05-10T21:21:15.952242Z","shell.execute_reply":"2022-05-10T21:21:52.759512Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Outliers removal","metadata":{}},{"cell_type":"markdown","source":"IQR method was tried, but we chose not to include it since it removes very low and very high temperatures from the training, which is useful information","metadata":{}},{"cell_type":"code","source":"'''\npercentile25 = train_df['fact_temperature'].quantile(0.25)\npercentile75 = train_df['fact_temperature'].quantile(0.75)\niqr = percentile75 - percentile25\n\nupper_limit = percentile75 + 1.5 * iqr\nlower_limit = percentile25 - 1.5 * iqr\n\nprint(f\"Temperature upper limit: {upper_limit}\")\nprint(f\"Temperature lower limit: {lower_limit}\")\n\nprint(f\"Number of data that exceeds temperature upper limit: {train_df[train_df['fact_temperature'] > upper_limit].shape[0]}\")\nprint(f\"Number of data that exceeds temperature lower limit: {train_df[train_df['fact_temperature'] < lower_limit].shape[0]}\")\n\nprint(f\"Len of train_df: {train_df.shape[0]} \")\n\ntrain_df =  train_df.drop( train_df[train_df['fact_temperature'] > upper_limit].index) # drop record over the upper limit\ntrain_df =  train_df.drop( train_df[train_df['fact_temperature'] < lower_limit].index) # drop record under the upper limit\n\nprint(f\"Len of train_df: {train_df.shape[0]}\")\n'''\n\n","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:21:09.703520Z","iopub.execute_input":"2022-05-10T20:21:09.704270Z","iopub.status.idle":"2022-05-10T20:21:09.712113Z","shell.execute_reply.started":"2022-05-10T20:21:09.704214Z","shell.execute_reply":"2022-05-10T20:21:09.711425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Features Engineering","metadata":{}},{"cell_type":"markdown","source":"We create new features from the timestamp feature. In particular, we extract the month and the our from the UNIX timestamp, since we believe they are the most useful information to predict the temperature.\nThen, we one hot encoded them not to create a misleading relation between them and the target variable. For example an higher month does not translate to higher/lower temperatures.","metadata":{}},{"cell_type":"code","source":"from datetime import datetime\n\ndef get_hour(date, n_intervals = 6):\n    \n    hour = date.hour\n    hours_in_interval = 24 // n_intervals\n\n    return hour // hours_in_interval\n\n\ndef get_categorical_month_and_hour(df):\n\n    datetimes = df[\"fact_time\"].apply(datetime.utcfromtimestamp)\n\n    df[\"month\"] = datetimes.apply(lambda date: date.month)\n    df[\"hour\"] = datetimes.apply(get_hour) # get hour in intervals\n\n    return pd.get_dummies(data=df, columns=['month', 'hour'])","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:21:52.761713Z","iopub.execute_input":"2022-05-10T21:21:52.762100Z","iopub.status.idle":"2022-05-10T21:21:52.769419Z","shell.execute_reply.started":"2022-05-10T21:21:52.762068Z","shell.execute_reply":"2022-05-10T21:21:52.768270Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_df = get_categorical_month_and_hour(train_df)\n\ntrain_df = train_df.drop(columns=[\"fact_time\"])","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:21:52.770736Z","iopub.execute_input":"2022-05-10T21:21:52.771174Z","iopub.status.idle":"2022-05-10T21:22:24.425977Z","shell.execute_reply.started":"2022-05-10T21:21:52.771142Z","shell.execute_reply":"2022-05-10T21:22:24.425245Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Feature selection","metadata":{}},{"cell_type":"markdown","source":"Plotting the Covariance matrix : ","metadata":{}},{"cell_type":"code","source":"plt.figure(figsize = (13,8))\ncorr_matrix = train_df.corr().abs()\nsns.heatmap(corr_matrix)\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:28:18.292512Z","iopub.execute_input":"2022-05-10T20:28:18.292782Z","iopub.status.idle":"2022-05-10T20:29:44.249225Z","shell.execute_reply.started":"2022-05-10T20:28:18.292750Z","shell.execute_reply":"2022-05-10T20:29:44.248384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Plotting the covariance matrix between some gfs temperature features to vizualize the covariance between them.","metadata":{}},{"cell_type":"code","source":"columns_gfs_corr = ['gfs_temperature_7000','gfs_temperature_20000','gfs_temperature_35000','gfs_temperature_50000','gfs_temperature_65000','gfs_temperature_80000','gfs_temperature_92500']\ncorr_gfs_temp = train_df[columns_gfs_corr].corr()\nsns.heatmap(corr_gfs_temp)\nplt.xticks()\nplt.tight_layout()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:29:44.250393Z","iopub.execute_input":"2022-05-10T20:29:44.250636Z","iopub.status.idle":"2022-05-10T20:29:45.060320Z","shell.execute_reply.started":"2022-05-10T20:29:44.250606Z","shell.execute_reply":"2022-05-10T20:29:45.059364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We check the linear correlation between each pair of features. If two features are highly correlated, then we drop one of them.","metadata":{}},{"cell_type":"code","source":"def get_correlated_features(df, threshold, method='pearson'):\n    corr_matrix = df.corr(method=method).abs()\n    high_corr_var=np.where(corr_matrix>threshold)\n    high_corr_var=[(corr_matrix.columns[x],corr_matrix.columns[y]) for x,y in zip(*high_corr_var) if x!=y and x<y]\n    # features and how many times they appear in tuples of correlated features\n    features_count = pd.Series(sum(np.array(high_corr_var).tolist(), [])).value_counts()\n    to_del = set()\n    for pair in high_corr_var:\n        if features_count[pair[0]] > features_count[pair[1]]:\n            to_del.add(pair[0])\n        else:\n            to_del.add(pair[1])\n    return list(to_del)","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:22:24.427334Z","iopub.execute_input":"2022-05-10T21:22:24.428118Z","iopub.status.idle":"2022-05-10T21:22:24.437504Z","shell.execute_reply.started":"2022-05-10T21:22:24.428074Z","shell.execute_reply":"2022-05-10T21:22:24.436585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_cmc_gfs_wrf = train_df[train_df.columns.values[train_df.columns.str.contains(\"gfs\") | \n                                                  train_df.columns.str.contains(\"cmc\") |\n                                                  train_df.columns.str.contains(\"wrf\")]]\n\n# pearson: linear relation between two variables\nfeatures_to_remove_pearson = get_correlated_features(df_cmc_gfs_wrf, threshold=0.8, method='pearson')\n\nfeatures_to_remove = features_to_remove_pearson\nfeatures_to_remove","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:22:24.439100Z","iopub.execute_input":"2022-05-10T21:22:24.439517Z","iopub.status.idle":"2022-05-10T21:23:23.071052Z","shell.execute_reply.started":"2022-05-10T21:22:24.439474Z","shell.execute_reply":"2022-05-10T21:23:23.070017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We check correlation between target variable and other features\n","metadata":{}},{"cell_type":"code","source":"corr_target = train_df.corrwith(train_df[\"fact_temperature\"]).sort_values(ascending=False)\ncorr_target","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:32:15.264233Z","iopub.execute_input":"2022-05-10T20:32:15.265920Z","iopub.status.idle":"2022-05-10T20:32:19.192778Z","shell.execute_reply.started":"2022-05-10T20:32:15.265871Z","shell.execute_reply":"2022-05-10T20:32:19.191725Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We finally drop the features highly correlated with other features","metadata":{}},{"cell_type":"code","source":"final_train_df = train_df.drop(columns=features_to_remove)\nfinal_train_df","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:23:23.074018Z","iopub.execute_input":"2022-05-10T21:23:23.075549Z","iopub.status.idle":"2022-05-10T21:23:23.785010Z","shell.execute_reply.started":"2022-05-10T21:23:23.075491Z","shell.execute_reply":"2022-05-10T21:23:23.783949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Scaling of features","metadata":{}},{"cell_type":"markdown","source":"We used Standard Scaler in order to use a common scale for the features, without distorting differences in the range of values.","metadata":{}},{"cell_type":"code","source":"from sklearn.preprocessing import StandardScaler\n\n# Get X_train and y_train and normalize\nX_train = final_train_df.drop(columns=['fact_temperature'])\ny_train = final_train_df['fact_temperature']\n\nsc = StandardScaler()\nsc.fit(X_train)\nX_train = sc.transform(X_train)\n\nscy = StandardScaler()\nscy.fit(y_train.to_numpy().reshape(-1, 1))\ny_train = scy.transform(y_train.to_numpy().reshape(-1, 1))","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:23:23.786364Z","iopub.execute_input":"2022-05-10T21:23:23.786623Z","iopub.status.idle":"2022-05-10T21:23:29.741771Z","shell.execute_reply.started":"2022-05-10T21:23:23.786592Z","shell.execute_reply":"2022-05-10T21:23:29.740975Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# PCA feature selection","metadata":{}},{"cell_type":"code","source":"'''from sklearn.decomposition import PCA\nfrom sklearn.preprocessing import StandardScaler\n\n\nX_sf = X_train.copy()\n\n\npercentage_of_variance = 0.95\npca = PCA(percentage_of_variance)\nX_train_pca = pca.fit_transform(X_sf)\n'''","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:42:52.339813Z","iopub.execute_input":"2022-05-10T20:42:52.340141Z","iopub.status.idle":"2022-05-10T20:42:52.347149Z","shell.execute_reply.started":"2022-05-10T20:42:52.340105Z","shell.execute_reply":"2022-05-10T20:42:52.346173Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Forward selection","metadata":{}},{"cell_type":"code","source":"'''\n# Feature selection iterative model\nimport statsmodels.api as sm\nfrom statsmodels.formula.api import ols\n\n# print(X_train)\n# print(y_train)\n\nremaining_features = list(X_train.columns)\n\ndata = X_train.copy()\ndata[\"y\"] = y_train[:]\n\nprint(remaining_features)\n\nalpha = 0.025\n\n# provo tutte le colonne e se il p value è minore di 0.05 aggiungo a una lista, poi prendo la colonna con il p-value minore\n# e la aggiungo al mio modello. Poi riparto con il mio modello aggiornato\n\nselected_features_forward = []\nwhile remaining_features: \n  PF = []  #list of (P value, feature)\n  print(f\"SELECTED FEATURES: {selected_features_forward}\")\n  for f in remaining_features:\n    temp = selected_features_forward + [f]  #temporary list of features+\n\n    log_reg = sm.Logit(y_train, X_train[temp]).fit()\n\n    p_values = log_reg.pvalues\n    min_p_value = np.min(log_reg.pvalues)\n\n    if min_p_value < alpha:\n       PF.append((min_p_value,f))\n  if PF:  #if not empty\n     PF.sort(reverse=True)\n     (best_pval, best_f) = PF.pop()\n     remaining_features.remove(best_f)\n     print('selected feature {} with p-value = {:.2E}'.\n            format(best_f, best_pval))\n     selected_features_forward.append(best_f)\n  else:\n     break\n\nprint(\"---------------------------------------------------\")\nprint(selected_features_forward)\nprint(data)\n\n'''","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:42:53.684131Z","iopub.execute_input":"2022-05-10T20:42:53.684430Z","iopub.status.idle":"2022-05-10T20:42:53.691914Z","shell.execute_reply.started":"2022-05-10T20:42:53.684395Z","shell.execute_reply":"2022-05-10T20:42:53.691076Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Training","metadata":{}},{"cell_type":"code","source":"# grid search cv\n\nparam_grid = {\n            \"alpha\": [0.1, 0.01, 0.001]\n            }\ngs = GridSearchCV(estimator=Lasso(), param_grid=param_grid, \n                  scoring=\"neg_root_mean_squared_error\",cv=5, verbose=True)\ngs.fit(X_train,y_train)","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:44:43.548872Z","iopub.execute_input":"2022-05-10T20:44:43.549211Z","iopub.status.idle":"2022-05-10T20:48:25.287745Z","shell.execute_reply.started":"2022-05-10T20:44:43.549179Z","shell.execute_reply":"2022-05-10T20:48:25.286405Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f\"Best set of parameters: {gs.best_params_}\")\nprint(f\"Best score: {gs.best_score_}\")","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2022-05-10T20:48:25.294925Z","iopub.execute_input":"2022-05-10T20:48:25.295395Z","iopub.status.idle":"2022-05-10T20:48:25.310480Z","shell.execute_reply.started":"2022-05-10T20:48:25.295342Z","shell.execute_reply":"2022-05-10T20:48:25.309322Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Now train with full train dataset and best parameters\nmodel = Lasso(alpha=0.001)\nmodel.fit(X_train, y_train)","metadata":{"execution":{"iopub.status.busy":"2022-05-10T20:50:31.656925Z","iopub.execute_input":"2022-05-10T20:50:31.657262Z","iopub.status.idle":"2022-05-10T20:50:54.829742Z","shell.execute_reply.started":"2022-05-10T20:50:31.657231Z","shell.execute_reply":"2022-05-10T20:50:54.828623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Train using a neural network\nfrom keras.models import Sequential\nfrom keras.layers import Dense\nfrom keras.layers import Dropout\nfrom keras.layers import BatchNormalization\nfrom keras.callbacks import EarlyStopping\n\n\ndef model():\n    model = Sequential()\n    model.add(Dense(64, activation=\"relu\", kernel_initializer='he_normal'))\n    model.add(BatchNormalization())\n    model.add(Dense(32, activation=\"relu\",  kernel_initializer='he_normal'))\n    model.add(BatchNormalization())\n    model.add(Dense(32, activation=\"relu\",  kernel_initializer='he_normal'))\n    model.add(BatchNormalization())\n    model.add(Dense(16, activation=\"relu\",  kernel_initializer='he_normal'))\n    model.add(BatchNormalization())\n    model.add(Dense(1))\n    model.compile(loss=\"mean_squared_error\", optimizer=\"adam\")\n    return model\n\nmodel = model()\n\ncallback = EarlyStopping(monitor='loss', patience=3)\n\nhistory = model.fit(\n        X_train,\n        y_train,\n        epochs=50,\n        batch_size=1024,\n        validation_split=0.2,\n        verbose=True,\n        callbacks=[callback]\n)","metadata":{"scrolled":true,"execution":{"iopub.status.busy":"2022-05-10T21:23:29.743134Z","iopub.execute_input":"2022-05-10T21:23:29.743673Z","iopub.status.idle":"2022-05-10T21:35:40.791089Z","shell.execute_reply.started":"2022-05-10T21:23:29.743629Z","shell.execute_reply":"2022-05-10T21:35:40.789683Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_history(losses, val_losses):\n    plt.figure(figsize=(12,8))\n    plt.plot(losses)\n    plt.plot(val_losses)\n    plt.legend([\"train_loss\", \"val_loss\"], fontsize=20)\n    plt.xlabel(\"Epochs\", fontsize=20)\n    plt.ylabel(\"Loss\", fontsize=20)\n    plt.xticks(fontsize=20)\n    plt.yticks(fontsize=20)\n    plt.grid()\n    plt.savefig('train_val_loss.png')\n    plt.show()\n    \nplot_history(\n        history.history[\"loss\"],\n        history.history[\"val_loss\"]\n    )\n\n","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:36:34.565383Z","iopub.execute_input":"2022-05-10T21:36:34.566677Z","iopub.status.idle":"2022-05-10T21:36:34.987585Z","shell.execute_reply.started":"2022-05-10T21:36:34.566615Z","shell.execute_reply":"2022-05-10T21:36:34.986630Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Test","metadata":{}},{"cell_type":"markdown","source":"We load the test model and plot the distribution of the measurements around the globe","metadata":{}},{"cell_type":"code","source":"test_df = pd.read_csv(\"../input/eurecom-aml-2022-challenge-1/public/test_feat.csv\", low_memory=True)\ntest_df.plot(kind=\"scatter\", x=\"fact_longitude\", y=\"fact_latitude\", alpha=0.4,\n              cmap=plt.get_cmap(\"jet\"), figsize=(15,7))\nplt.show()\nindices = test_df.pop(\"index\")","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:36:40.690245Z","iopub.execute_input":"2022-05-10T21:36:40.690541Z","iopub.status.idle":"2022-05-10T21:37:05.077219Z","shell.execute_reply.started":"2022-05-10T21:36:40.690511Z","shell.execute_reply":"2022-05-10T21:37:05.076172Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Preprocessing of test data","metadata":{}},{"cell_type":"code","source":"# add info about month and hour\ntest_df = get_categorical_month_and_hour(test_df)\n# drop features\nfinal_test_df = test_df.drop(columns=features_to_remove)\nfinal_test_df = final_test_df.drop(columns=[\"fact_time\", \"gfs_available\", \"cmc_available\", \"wrf_available\"])\n\n# Get X_test and normalize\nX_test = final_test_df.to_numpy()\nX_test = sc.transform(X_test)","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:37:05.079072Z","iopub.execute_input":"2022-05-10T21:37:05.079323Z","iopub.status.idle":"2022-05-10T21:37:13.907445Z","shell.execute_reply.started":"2022-05-10T21:37:05.079292Z","shell.execute_reply":"2022-05-10T21:37:13.906440Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Prediction of test data","metadata":{}},{"cell_type":"code","source":"y_pred = model.predict(X_test)\ny_pred = scy.inverse_transform(y_pred)","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:37:13.908710Z","iopub.execute_input":"2022-05-10T21:37:13.908965Z","iopub.status.idle":"2022-05-10T21:37:42.571022Z","shell.execute_reply.started":"2022-05-10T21:37:13.908934Z","shell.execute_reply":"2022-05-10T21:37:42.569863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_df = pd.DataFrame(data={'index': indices.values,\n                                   'fact_temperature': y_pred.squeeze()})\n\n# Save the predictions into a csv file\n# Notice that this file should be saved under the directory `/kaggle/working` \n# so that you can download it later\nsubmission_df.to_csv(\"/kaggle/working/submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:37:42.573464Z","iopub.execute_input":"2022-05-10T21:37:42.573812Z","iopub.status.idle":"2022-05-10T21:37:44.055651Z","shell.execute_reply.started":"2022-05-10T21:37:42.573766Z","shell.execute_reply":"2022-05-10T21:37:44.054569Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2022-05-10T21:37:44.056991Z","iopub.execute_input":"2022-05-10T21:37:44.057250Z","iopub.status.idle":"2022-05-10T21:37:45.380769Z","shell.execute_reply.started":"2022-05-10T21:37:44.057212Z","shell.execute_reply":"2022-05-10T21:37:45.379623Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}