{"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 matplotlib.pyplot as plt \nimport numpy as np \nimport seaborn as sns\nfrom tqdm.notebook import tqdm \n\nimport matplotlib\nfrom matplotlib import rc\n\nrc(\"font\", **{\"family\": \"sans-serif\", \"sans-serif\": \"DejaVu Sans\"})\nrc(\"figure\", **{\"dpi\": 200})\nrc(\n    \"axes\",\n    **{\"spines.right\": False, \"spines.top\": False, \"xmargin\": 0.0, \"ymargin\": 0.05}\n)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-05-11T19:06:42.417218Z","iopub.execute_input":"2022-05-11T19:06:42.417643Z","iopub.status.idle":"2022-05-11T19:06:43.718493Z","shell.execute_reply.started":"2022-05-11T19:06:42.417521Z","shell.execute_reply":"2022-05-11T19:06:43.717443Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data import","metadata":{}},{"cell_type":"code","source":"df = pd.read_csv('../input/eurecom-aml-2022-challenge-1/public/train.csv', low_memory=True)\ndf_test = pd.read_csv('../input/eurecom-aml-2022-challenge-1/public/test_feat.csv', low_memory=True)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:06:43.720411Z","iopub.execute_input":"2022-05-11T19:06:43.720728Z","iopub.status.idle":"2022-05-11T19:08:29.697660Z","shell.execute_reply.started":"2022-05-11T19:06:43.720691Z","shell.execute_reply":"2022-05-11T19:08:29.696540Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.head()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:29.699408Z","iopub.execute_input":"2022-05-11T19:08:29.699733Z","iopub.status.idle":"2022-05-11T19:08:29.755571Z","shell.execute_reply.started":"2022-05-11T19:08:29.699690Z","shell.execute_reply":"2022-05-11T19:08:29.754663Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test.head()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:29.758156Z","iopub.execute_input":"2022-05-11T19:08:29.758670Z","iopub.status.idle":"2022-05-11T19:08:29.789589Z","shell.execute_reply.started":"2022-05-11T19:08:29.758622Z","shell.execute_reply":"2022-05-11T19:08:29.788359Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data analysis - PART 1\n\n* First, we look at the coordinates and samples present in the training and testing sets.","metadata":{}},{"cell_type":"code","source":"tr_coordinates = df[['fact_latitude','fact_longitude']].drop_duplicates()\nte_coordinates = df_test[['fact_latitude','fact_longitude']].drop_duplicates()\nprint (tr_coordinates)\nprint (te_coordinates)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:29.791288Z","iopub.execute_input":"2022-05-11T19:08:29.791743Z","iopub.status.idle":"2022-05-11T19:08:30.010221Z","shell.execute_reply.started":"2022-05-11T19:08:29.791702Z","shell.execute_reply":"2022-05-11T19:08:30.009190Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"merged = te_coordinates.merge(tr_coordinates, how='left', indicator=True)\ntropic_coordinates = merged[merged['_merge']=='left_only']\ncommon_coordinates = merged[merged['_merge']=='both']","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:30.011561Z","iopub.execute_input":"2022-05-11T19:08:30.011803Z","iopub.status.idle":"2022-05-11T19:08:30.060824Z","shell.execute_reply.started":"2022-05-11T19:08:30.011772Z","shell.execute_reply":"2022-05-11T19:08:30.060053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, (ax0, ax1) = plt.subplots(1, 2, figsize=[10, 3], dpi=200, sharex=True, sharey=True)\nax0.scatter(tr_coordinates['fact_longitude'], tr_coordinates['fact_latitude'], \n            s=1, c='tab:green', label='Train data')\nax1.scatter(common_coordinates['fact_longitude'], common_coordinates['fact_latitude'], \n            s=1, c='tab:orange', label='In-domain test data')\nax1.scatter(tropic_coordinates['fact_longitude'], tropic_coordinates['fact_latitude'], \n            s=1, c='tab:red', label='Out-domain test data')\nax0.legend(), ax1.legend()\nfig.suptitle('Distribution of train and test points in the dataset', y=1.01, fontsize=14)\nplt.savefig(\"Distribution\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:30.062009Z","iopub.execute_input":"2022-05-11T19:08:30.062810Z","iopub.status.idle":"2022-05-11T19:08:31.534562Z","shell.execute_reply.started":"2022-05-11T19:08:30.062771Z","shell.execute_reply":"2022-05-11T19:08:31.533595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* Before analyzing further, we wanted to proceed with feature extraction and the rest of the pre-processing, so that focusing on the most important features becomes easier.","metadata":{}},{"cell_type":"markdown","source":"# Data Pre-processing\n\nWe used general pre-processing steps to get a better data quality. First of all, we delete the constant columns, then the data is normalised and standardised. Also, based on the correlation between the temperature and all features, we keep only the features for which this value is above a certain threshold.\n\nHowever, the normalisation and standardisation are not critical since we decided to use a decision tree. It also means we do not necessarily have to worry about outliers, missing values, and whether the values are categorical or numerical (as long as booleans values are replaced by 0 and 1).\n\n\nBesides, the original method we used to extract features was Principal Component Analysis (PCA). Still, the correlation method achieved better results, and has a huge additional benefit: features are kept intact, providing better interpretability for the tree.","metadata":{}},{"cell_type":"code","source":"df.replace([np.inf, -np.inf], np.nan, inplace=True)\ndf.replace(True, 1, inplace=True)\ndf.replace(False, 0, inplace=True)\nX = df.dropna().iloc[:, 1:-1].values\ny = df.dropna().iloc[:, -1].values","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:31.536420Z","iopub.execute_input":"2022-05-11T19:08:31.537094Z","iopub.status.idle":"2022-05-11T19:08:35.871783Z","shell.execute_reply.started":"2022-05-11T19:08:31.537040Z","shell.execute_reply":"2022-05-11T19:08:35.870779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#PCA feature selection (unused here)\ndef PCA_feature_selection(X) :\n    covariance_matrix = np.cov(X.T)\n    eigen_values, eigen_vectors = np.linalg.eig(covariance_matrix)\n    variance_explained = []\n    args=np.argsort(eigen_values)\n    args=np.flip(args)\n    eigen_values=eigen_values[args]\n    eigen_vectors=eigen_vectors[args]\n    for i in eigen_values:\n        variance_explained.append((i/sum(eigen_values))*100)\n\n    cumulative_variance_explained = np.cumsum(variance_explained)\n\n    projection_matrix = (eigen_vectors.T[:][:22]).T #22 first components to represent 80% of the data\n    X_pca = X.dot(projection_matrix)\n    return X_pca, projection_matrix","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:35.873169Z","iopub.execute_input":"2022-05-11T19:08:35.873536Z","iopub.status.idle":"2022-05-11T19:08:35.881554Z","shell.execute_reply.started":"2022-05-11T19:08:35.873491Z","shell.execute_reply":"2022-05-11T19:08:35.880531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def corr_feature_selection(X,y, threshold = 0.7) :\n    cor=np.zeros(len(X[0]))\n    relevant_features=[]\n    for i in range (len(X[0])) :\n        if i==2 or i==3 :\n            relevant_features.append(i) #latitude and longitude must be kept to deal with different climates\n        cor[i]=np.corrcoef(X[:, i], y)[0,1]\n        if (abs(cor[i]))>threshold :\n            relevant_features.append(i)\n          \n    cor =cor[abs(cor)>threshold]\n    return cor, relevant_features","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:35.884658Z","iopub.execute_input":"2022-05-11T19:08:35.884943Z","iopub.status.idle":"2022-05-11T19:08:35.900297Z","shell.execute_reply.started":"2022-05-11T19:08:35.884910Z","shell.execute_reply":"2022-05-11T19:08:35.899219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Xmean, Xstd, ymean, ystd = X.mean(0), X.std(0), y.mean(), y.std()\n\n#Delete rows of X and values of y where the standard deviation is 0, they are useless because constant\nzeros=[]\nfor i in range(len(Xstd)) :\n    if Xstd[i]==0 :\n        zeros.append(i)\n        \nX=np.delete(X,zeros,1) \n\n#Xmean and Xstd are vectors, so we delete along the axis 0\nXmean=np.delete(Xmean,zeros,0)\nXstd=np.delete(Xstd,zeros,0)\n\n\nX = (X - Xmean) / Xstd\ny = (y - ymean) / ystd","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:35.901723Z","iopub.execute_input":"2022-05-11T19:08:35.901957Z","iopub.status.idle":"2022-05-11T19:08:40.634030Z","shell.execute_reply.started":"2022-05-11T19:08:35.901929Z","shell.execute_reply":"2022-05-11T19:08:40.633051Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"cor, relevant_features=corr_feature_selection(X,y)\n\n#Feature extraction\nXrf =X[:,relevant_features]\n\nfeat = df.columns[relevant_features].to_numpy()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:40.635468Z","iopub.execute_input":"2022-05-11T19:08:40.636446Z","iopub.status.idle":"2022-05-11T19:08:43.463877Z","shell.execute_reply.started":"2022-05-11T19:08:40.636377Z","shell.execute_reply":"2022-05-11T19:08:43.462759Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Features kept are:\")\nprint(feat)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:43.465207Z","iopub.execute_input":"2022-05-11T19:08:43.465700Z","iopub.status.idle":"2022-05-11T19:08:43.471712Z","shell.execute_reply.started":"2022-05-11T19:08:43.465647Z","shell.execute_reply":"2022-05-11T19:08:43.470802Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data analysis - PART 2\n\n* Now that the main features are kept and some statistics have been computed, we can resume some analytics about temperature.","metadata":{}},{"cell_type":"code","source":"#Denormalise after the feature selection\nX *= Xstd\nX += Xmean\ny *= ystd\ny += ymean\n\nrdmean = round(ymean*100)/100","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:43.474984Z","iopub.execute_input":"2022-05-11T19:08:43.475629Z","iopub.status.idle":"2022-05-11T19:08:43.816777Z","shell.execute_reply.started":"2022-05-11T19:08:43.475572Z","shell.execute_reply":"2022-05-11T19:08:43.815728Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* We can already see the relationship between the final features and temperature (graphs are not displayed as of now to avoid having too much information shown). Obviously, there is no feature with a perfect relationship to temperature, but they can provide useful information despite the variance.","metadata":{}},{"cell_type":"code","source":"def plot_temperature(feature_index):  \n    feature = X[:,feature_index]\n    funique, indices = np.unique(feature, return_index = True)\n    yunique = y[indices]\n    plt.plot(funique,yunique)\n    plt.ylabel('Temperature in °C')\n    plt.xlabel(feat[feature_index])\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:43.818926Z","iopub.execute_input":"2022-05-11T19:08:43.819310Z","iopub.status.idle":"2022-05-11T19:08:43.827216Z","shell.execute_reply.started":"2022-05-11T19:08:43.819258Z","shell.execute_reply":"2022-05-11T19:08:43.826199Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"* We can now see the distribution of temperature in the data, and look at how the different samples at each location behave with respect to temperature.","metadata":{}},{"cell_type":"code","source":"sns.set_theme(style=\"darkgrid\")\n\np = sns.displot(df['fact_temperature'], kde = False, height = 3, aspect = 2)\np.set_axis_labels(\"Temperature (°C)\", \"Number of samples\")\nplt.axvline(x=ymean,\n            color='green', label = \"Mean=\" + str(rdmean), lw = 0.5)\nplt.legend()\nplt.savefig(\"Temperature histogram\")","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:43.828672Z","iopub.execute_input":"2022-05-11T19:08:43.829255Z","iopub.status.idle":"2022-05-11T19:08:47.670824Z","shell.execute_reply.started":"2022-05-11T19:08:43.829202Z","shell.execute_reply":"2022-05-11T19:08:47.670071Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"avg_tmp = df.groupby(['fact_longitude','fact_latitude']).mean().reset_index().rename(columns = {\"fact_temperature\": \"temperature_mean\"})\nstd_tmp = df.groupby(['fact_longitude','fact_latitude']).std().reset_index().rename(columns = {\"fact_temperature\": \"temperature_std\"})","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:47.671948Z","iopub.execute_input":"2022-05-11T19:08:47.672655Z","iopub.status.idle":"2022-05-11T19:08:52.791180Z","shell.execute_reply.started":"2022-05-11T19:08:47.672610Z","shell.execute_reply":"2022-05-11T19:08:52.790470Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots()\navg_tmp.plot(x='fact_longitude', y = 'fact_latitude', c='temperature_mean', figsize = (14,5),colorbar=True, colormap = plt.get_cmap(\"rainbow\"), ax = ax, kind = 'scatter')\n\nax.set_xlabel(\"Longitude (°)\")\nax.set_ylabel(\"Latitude(°)\")\nfig.suptitle('Average temperature per location in the training samples', y=1.01, fontsize=14)\nplt.savefig(\"Temperature location\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:52.792402Z","iopub.execute_input":"2022-05-11T19:08:52.792761Z","iopub.status.idle":"2022-05-11T19:08:54.258266Z","shell.execute_reply.started":"2022-05-11T19:08:52.792731Z","shell.execute_reply":"2022-05-11T19:08:54.257540Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots()\nstd_tmp.plot(x='fact_longitude', y = 'fact_latitude', c='temperature_std', figsize = (14,5),colorbar=True, colormap = plt.get_cmap(\"nipy_spectral\"), ax = ax, kind = 'scatter')\n\nax.set_xlabel(\"Longitude (°)\")\nax.set_ylabel(\"Latitude(°)\")\nfig.suptitle('Standard deviation for temperature per location in the training samples', y=1.01, fontsize=14)\nplt.savefig(\"Std Temperature location\")\nplt.show()\n","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:54.259557Z","iopub.execute_input":"2022-05-11T19:08:54.259930Z","iopub.status.idle":"2022-05-11T19:08:55.802628Z","shell.execute_reply.started":"2022-05-11T19:08:54.259899Z","shell.execute_reply":"2022-05-11T19:08:55.800968Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# Model Selection\n\n* We decided to use a decision tree for the regression, because it has excellent benefits for this kind of task. Indeed, the need for pre-processing is reduced as explained before, it is a simple model that trains and predicts relatively fast, and above all, it is easily interpretable since its splits are determined upon the original features. We leave the maximum depth as a hyper-parameter to be tuned in the next steps, by initializing it to a small value.\n\n* Naturally, we wondered whether random forests would be a better choice. Yet, they lose many advantages of the decision tree, namely the interpretability, simplicity and are a lot more computationally intensive. Only a much better performance than with the tree could have made a difference, but that was not the case using standard random forests. That is why we stayed with the decision tree as the final model.","metadata":{}},{"cell_type":"code","source":"#We use y again, not X\ny = (y - ymean) / ystd","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:55.804127Z","iopub.execute_input":"2022-05-11T19:08:55.804383Z","iopub.status.idle":"2022-05-11T19:08:55.815127Z","shell.execute_reply.started":"2022-05-11T19:08:55.804353Z","shell.execute_reply":"2022-05-11T19:08:55.814221Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\nXtr, Xval, ytr, yval = train_test_split(Xrf, y, random_state=1, test_size=10000)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:55.816983Z","iopub.execute_input":"2022-05-11T19:08:55.818677Z","iopub.status.idle":"2022-05-11T19:08:57.525141Z","shell.execute_reply.started":"2022-05-11T19:08:55.818596Z","shell.execute_reply":"2022-05-11T19:08:57.523779Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.tree import DecisionTreeRegressor\n\ndef compute_rmse(y, ypred, ystd=1.):\n    return np.mean((y - ypred)**2)**0.5 * ystd\n\ntreemod = DecisionTreeRegressor(max_depth=3, criterion='squared_error')\ntreemod.fit(Xtr, ytr)\nypred_tr = treemod.predict(Xtr)\nypred_val = treemod.predict(Xval)\n\nprint(f'Train RMSE: {compute_rmse(ytr, ypred_tr, ystd):.3f}')\nprint(f'Valid RMSE: {compute_rmse(yval, ypred_val, ystd):.3f}')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:08:57.527321Z","iopub.execute_input":"2022-05-11T19:08:57.527724Z","iopub.status.idle":"2022-05-11T19:09:20.557375Z","shell.execute_reply.started":"2022-05-11T19:08:57.527666Z","shell.execute_reply":"2022-05-11T19:09:20.556531Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Parameter Optimisation\n\nWe kept a grid search with cross-validation to tune the maximum depth of the tree. Another try was to modify the impurity criterion by using the Friedman version of the mean square error, but it only gave small improvements. \n\nWe can only test few combinations for the parameters, because it takes too much time otherwise. The most important enhancement comes from the depth. The model almost does not overfit, as the cross-validation hints, which is highly desirable, even if it may means that we could make it more complex to achieve better performance. ","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import RandomizedSearchCV, GridSearchCV","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:09:20.558521Z","iopub.execute_input":"2022-05-11T19:09:20.558733Z","iopub.status.idle":"2022-05-11T19:09:20.563755Z","shell.execute_reply.started":"2022-05-11T19:09:20.558705Z","shell.execute_reply":"2022-05-11T19:09:20.562848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\n\ndef get_cv_idx(n, test_size=0.2, n_splits=2):\n    train_idx, test_idx = [], []\n    for _ in range(n_splits):\n        idx = np.random.permutation(n)\n        train_size = int(n * (1 - test_size)) if isinstance(test_size, float) else n - test_size\n        train_idx.append(idx[:train_size])\n        test_idx.append(idx[train_size:])\n    return train_idx, test_idx","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:09:20.565666Z","iopub.execute_input":"2022-05-11T19:09:20.566009Z","iopub.status.idle":"2022-05-11T19:09:20.578806Z","shell.execute_reply.started":"2022-05-11T19:09:20.565963Z","shell.execute_reply":"2022-05-11T19:09:20.577831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"train_idx, cv_idx = get_cv_idx(len(Xtr), test_size=10000, n_splits=10)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:09:20.580188Z","iopub.execute_input":"2022-05-11T19:09:20.580690Z","iopub.status.idle":"2022-05-11T19:09:21.523401Z","shell.execute_reply.started":"2022-05-11T19:09:20.580653Z","shell.execute_reply":"2022-05-11T19:09:21.522203Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"param_grid = {\n    \"max_depth\": [12, 15],\n    \"criterion\": ['squared_error', 'friedman_mse']\n}\n\nsearch = GridSearchCV(treemod, \n                      param_grid,\n                      n_jobs=-1, \n                      verbose=1,\n                      cv=zip(train_idx, cv_idx), \n                      scoring='neg_root_mean_squared_error').fit(Xtr, ytr)\nprint('Done!')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:09:21.525061Z","iopub.execute_input":"2022-05-11T19:09:21.525950Z","iopub.status.idle":"2022-05-11T19:31:08.340108Z","shell.execute_reply.started":"2022-05-11T19:09:21.525895Z","shell.execute_reply":"2022-05-11T19:31:08.339143Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Best parameters set found on cv set:\")\nprint(search.best_params_)\nprint()\nprint(\"Grid scores on cv set:\")\nmeans = search.cv_results_[\"mean_test_score\"]\nstds = search.cv_results_[\"std_test_score\"]\nfor mean, std, params in zip(means, stds, search.cv_results_[\"params\"]):\n    print(\"%0.3f (+/-%0.03f) for %r\" % (-mean * ystd, (std * ystd) * 2, params))\nprint()\nprint(\"Error on the validation set\")\nypred_val = search.predict(Xval)\nprint(f'Valid RMSE: {compute_rmse(yval, ypred_val, ystd):.3f}')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:31:08.342411Z","iopub.execute_input":"2022-05-11T19:31:08.342897Z","iopub.status.idle":"2022-05-11T19:31:08.358303Z","shell.execute_reply.started":"2022-05-11T19:31:08.342846Z","shell.execute_reply":"2022-05-11T19:31:08.357310Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Evaluation\n","metadata":{}},{"cell_type":"markdown","source":"The time complexity of the decision tree is only O(depth)=O(log(N)), with N being the number of training samples, since it is a binary tree. That means it provides a really quick prediction for testing points, as we can verify here (about a tenth of a second or less). We can also use the tree structure to understand the splits. Unexpectedly, the first three splits were done using the feature 'wrf_wind_v', before using other features. That means splitting using this feature based on the wind is more relevant than starting with a feature we would have considered more intuitive, like the latitude.\n\nWe used the R2 score to get another metric to understand how well the model can expect to perform. We obtained a score of above 0.90. Since a score of 1 is the best possible one, that means a lot of the variance has been captured thanks to the independent variables we used. Therefore, we can expect a fair generalization, which is what also seemed the case after cross-validation.","metadata":{}},{"cell_type":"code","source":"#To visualize the tree, we need to call best_estimator\nsbe = search.best_estimator_\nsbe.fit(Xtr, ytr)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:31:08.364453Z","iopub.execute_input":"2022-05-11T19:31:08.364736Z","iopub.status.idle":"2022-05-11T19:32:29.049152Z","shell.execute_reply.started":"2022-05-11T19:31:08.364703Z","shell.execute_reply":"2022-05-11T19:32:29.048115Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Whole structure of the three\n\nfrom sklearn import tree\nimport graphviz \ndot_data = tree.export_graphviz(sbe, out_file=None, \n                     feature_names=feat,  \n                     class_names=None,  \n                     filled=True, rounded=True,  \n                     impurity=True)  \ngraph = graphviz.Source(dot_data)  \ngraph.render(\"Decision tree\") ","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:32:29.050778Z","iopub.execute_input":"2022-05-11T19:32:29.051065Z","iopub.status.idle":"2022-05-11T19:33:38.031731Z","shell.execute_reply.started":"2022-05-11T19:32:29.051031Z","shell.execute_reply":"2022-05-11T19:33:38.030513Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"score = sbe.score(Xtr, ytr)\nprint(\"R^2 score on the training set:\", score)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:38.033796Z","iopub.execute_input":"2022-05-11T19:33:38.034170Z","iopub.status.idle":"2022-05-11T19:33:38.515048Z","shell.execute_reply.started":"2022-05-11T19:33:38.034122Z","shell.execute_reply":"2022-05-11T19:33:38.513949Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Visualise the first splits of the tree\n\nn_nodes = sbe.tree_.node_count\ncleft = sbe.tree_.children_left\ncright = sbe.tree_.children_right\nfeature = sbe.tree_.feature\nthreshold = sbe.tree_.threshold\n\nnode_depth = np.zeros(shape=n_nodes, dtype=np.int64)\nis_leaves = np.zeros(shape=n_nodes, dtype=bool)\nstack = [(0, 0)]  # start with the root node id (0) and its depth (0)\nwhile len(stack) > 0:\n    # Depile\n    node_id, depth = stack.pop()\n    node_depth[node_id] = depth\n\n    # The node is a split if its two children are different\n    is_split_node = cleft[node_id] != cright[node_id]\n    # Add children to the process\n    if is_split_node:\n        stack.append((cleft[node_id], depth + 1))\n        stack.append((cright[node_id], depth + 1))\n    else:\n        is_leaves[node_id] = True\n\nprint(\n    \"The binary tree structure has {n} nodes and has \"\n    \"the following splits for the three first layers:\\n\".format(n=n_nodes)\n)\nfor i in range(n_nodes):\n    if node_depth[i] <= 3:\n        findex = feature[i]\n        print(\n            \"{space}Node {node} test: \"\n            \"{feature} <= {threshold} \"\n            .format(\n                space=node_depth[i] * \"\\t\",\n                node=i,\n                feature=feat[findex],\n                threshold=threshold[i]\n            )\n        )\n  \n    ","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:38.516808Z","iopub.execute_input":"2022-05-11T19:33:38.517269Z","iopub.status.idle":"2022-05-11T19:33:38.563193Z","shell.execute_reply.started":"2022-05-11T19:33:38.517215Z","shell.execute_reply":"2022-05-11T19:33:38.561903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test = pd.read_csv('../input/eurecom-aml-2022-challenge-1/public/test_feat.csv', low_memory=True)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:38.564799Z","iopub.execute_input":"2022-05-11T19:33:38.565034Z","iopub.status.idle":"2022-05-11T19:33:53.746457Z","shell.execute_reply.started":"2022-05-11T19:33:38.565006Z","shell.execute_reply":"2022-05-11T19:33:53.745486Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#Feature extraction for the testing set\nXte = df_test.iloc[:, 1:].values\nXte=np.delete(Xte,zeros,1) \n\nXte = (Xte - Xmean) / Xstd\nXte = Xte[:,relevant_features]\n","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:53.747918Z","iopub.execute_input":"2022-05-11T19:33:53.748139Z","iopub.status.idle":"2022-05-11T19:33:54.545395Z","shell.execute_reply.started":"2022-05-11T19:33:53.748112Z","shell.execute_reply":"2022-05-11T19:33:54.544485Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from time import perf_counter","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:54.547580Z","iopub.execute_input":"2022-05-11T19:33:54.548349Z","iopub.status.idle":"2022-05-11T19:33:54.552590Z","shell.execute_reply.started":"2022-05-11T19:33:54.548311Z","shell.execute_reply":"2022-05-11T19:33:54.551494Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"start_pred = perf_counter()\n#De-normalize the predictions\nypred_te = search.predict(Xte) * ystd + ymean\nend_pred = perf_counter() ","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:54.554386Z","iopub.execute_input":"2022-05-11T19:33:54.555055Z","iopub.status.idle":"2022-05-11T19:33:54.691729Z","shell.execute_reply.started":"2022-05-11T19:33:54.555009Z","shell.execute_reply":"2022-05-11T19:33:54.690768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(\"Time to predict :\", end_pred - start_pred, \"seconds\")","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:54.693233Z","iopub.execute_input":"2022-05-11T19:33:54.693493Z","iopub.status.idle":"2022-05-11T19:33:54.699313Z","shell.execute_reply.started":"2022-05-11T19:33:54.693452Z","shell.execute_reply":"2022-05-11T19:33:54.698473Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"\n# Submission\n","metadata":{}},{"cell_type":"code","source":"submission_df = pd.DataFrame(data={'index': df_test['index'].values,\n                                   'fact_temperature': ypred_te.squeeze()})\n\n# Save the predictions into a csv file\nsubmission_df.to_csv(\"/kaggle/working/submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:54.700911Z","iopub.execute_input":"2022-05-11T19:33:54.701168Z","iopub.status.idle":"2022-05-11T19:33:56.576942Z","shell.execute_reply.started":"2022-05-11T19:33:54.701134Z","shell.execute_reply":"2022-05-11T19:33:56.575886Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Check the submission file\n! head -6 \"/kaggle/working/submission.csv\"","metadata":{"execution":{"iopub.status.busy":"2022-05-11T19:33:56.578227Z","iopub.execute_input":"2022-05-11T19:33:56.578491Z","iopub.status.idle":"2022-05-11T19:33:57.392749Z","shell.execute_reply.started":"2022-05-11T19:33:56.578460Z","shell.execute_reply":"2022-05-11T19:33:57.391441Z"},"trusted":true},"execution_count":null,"outputs":[]}]}