{"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":"# EURECOM AML 2022:: Challenge 1\n### Group 14","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\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-11T16:44:28.930500Z","iopub.execute_input":"2022-05-11T16:44:28.932583Z","iopub.status.idle":"2022-05-11T16:44:28.966327Z","shell.execute_reply.started":"2022-05-11T16:44:28.932200Z","shell.execute_reply":"2022-05-11T16:44:28.965198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train = 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-11T16:44:28.969701Z","iopub.execute_input":"2022-05-11T16:44:28.970819Z","iopub.status.idle":"2022-05-11T16:45:56.507534Z","shell.execute_reply.started":"2022-05-11T16:44:28.970780Z","shell.execute_reply":"2022-05-11T16:45:56.506237Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Exploratory Data Analysis","metadata":{}},{"cell_type":"markdown","source":"### Dataset information","metadata":{}},{"cell_type":"code","source":"df_train.head()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:45:56.509326Z","iopub.execute_input":"2022-05-11T16:45:56.509582Z","iopub.status.idle":"2022-05-11T16:45:56.562908Z","shell.execute_reply.started":"2022-05-11T16:45:56.509549Z","shell.execute_reply":"2022-05-11T16:45:56.562017Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Check number of features in training and test data, as well as their data types:","metadata":{}},{"cell_type":"code","source":"len(df_train.columns), len(df_test.columns)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:45:56.564568Z","iopub.execute_input":"2022-05-11T16:45:56.565562Z","iopub.status.idle":"2022-05-11T16:45:56.572730Z","shell.execute_reply.started":"2022-05-11T16:45:56.565514Z","shell.execute_reply":"2022-05-11T16:45:56.571858Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train.dtypes.value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:45:56.574881Z","iopub.execute_input":"2022-05-11T16:45:56.575260Z","iopub.status.idle":"2022-05-11T16:45:56.592043Z","shell.execute_reply.started":"2022-05-11T16:45:56.575213Z","shell.execute_reply":"2022-05-11T16:45:56.591270Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We only have numerical features in the dataset. Lets see how many unique coordinates we have:","metadata":{}},{"cell_type":"code","source":"print(\"Number of data points:\", len(df_train[[\"fact_longitude\", \"fact_latitude\"]]))\nprint(\"Number of unique coordinates:\", len(df_train[[\"fact_longitude\", \"fact_latitude\"]].drop_duplicates()))","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:45:56.593657Z","iopub.execute_input":"2022-05-11T16:45:56.594190Z","iopub.status.idle":"2022-05-11T16:45:56.739223Z","shell.execute_reply.started":"2022-05-11T16:45:56.594157Z","shell.execute_reply":"2022-05-11T16:45:56.738312Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are a lot of data per location, which can be useful for imputing missing values later. Lastly, lets plot the coordinates to get a feel of the distribution of data points:","metadata":{}},{"cell_type":"code","source":"coords = df_train[[\"fact_longitude\", \"fact_latitude\"]]\nfig, ax = plt.subplots(1, 1, figsize=[6, 3], dpi=200)\nax.scatter(coords['fact_longitude'], coords['fact_latitude'], \n            s=1, c='tab:green', alpha=0.1)\nplt.title(\"All unique data point locations\")\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:45:56.740491Z","iopub.execute_input":"2022-05-11T16:45:56.740731Z","iopub.status.idle":"2022-05-11T16:46:00.044627Z","shell.execute_reply.started":"2022-05-11T16:45:56.740702Z","shell.execute_reply":"2022-05-11T16:46:00.043882Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It looks like most coordinates have roughly the same number of samples.","metadata":{}},{"cell_type":"markdown","source":"### Missing data","metadata":{}},{"cell_type":"markdown","source":"From Kaggle:\n > Weather forecast model predictions are values produced by the following weather forecast models: Global Forecast System (GFS), Global Deterministic Forecast System from the Canadian Meteorological Center (CMC), Weather Research and Forecasting (WRF) Model. Each model returns the following predicted values: wind, humidity, pressure, clouds, precipitation, dew point, snow depth, air and soil temperature characteristics.\n \nWe divide the columns into the different categories \"general\", \"gfs\", \"cmc\" and \"wrf\" for better readability when plotting:","metadata":{}},{"cell_type":"code","source":"column_categories = {\n    \"general\": [\n        \"fact_time\",\n        \"fact_latitude\",\n        \"fact_longitude\",\n        \"topography_bathymetry\",\n        \"sun_elevation\",\n        \"climate_temperature\",\n        \"climate_pressure\",\n        \"fact_temperature\",\n    ]\n}\n\nfor prefix in [\"gfs\", \"cmc\", \"wrf\"]:\n    column_categories[prefix] = [\n        col for col in df_train.columns\n        if col.startswith(prefix) and not col.endswith(\"available\")\n    ]\n    \nprint(column_categories)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:00.045660Z","iopub.execute_input":"2022-05-11T16:46:00.046403Z","iopub.status.idle":"2022-05-11T16:46:00.053066Z","shell.execute_reply.started":"2022-05-11T16:46:00.046365Z","shell.execute_reply":"2022-05-11T16:46:00.052377Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We then plot the sum of rows with missing data per feature:","metadata":{}},{"cell_type":"code","source":"for prefix, columns in column_categories.items():\n    df = df_train[columns]\n    fig, ax = plt.subplots(1, 1, figsize=(5,2))\n    ax.set_title(prefix)\n    ax.tick_params(axis='x', labelsize=6)\n    df.isna().sum().plot.bar()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:00.054403Z","iopub.execute_input":"2022-05-11T16:46:00.054861Z","iopub.status.idle":"2022-05-11T16:46:02.862155Z","shell.execute_reply.started":"2022-05-11T16:46:00.054828Z","shell.execute_reply":"2022-05-11T16:46:02.861254Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It looks like rows either contains all or none of the features from the specific forecast models. Also, the \"general\" features contains no missing values. We can therefore summarize the number of rows with missing features as such:","metadata":{}},{"cell_type":"code","source":"df_train[[\"gfs_available\", \"cmc_available\", \"wrf_available\"]].value_counts()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:02.863411Z","iopub.execute_input":"2022-05-11T16:46:02.863648Z","iopub.status.idle":"2022-05-11T16:46:03.038043Z","shell.execute_reply.started":"2022-05-11T16:46:02.863620Z","shell.execute_reply":"2022-05-11T16:46:03.037098Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now plot maps showing which datapoints are missing feature data from the different weather forecast models:","metadata":{}},{"cell_type":"code","source":"for col in [\"gfs_available\", \"cmc_available\", \"wrf_available\"]:\n    df_available = df_train[df_train[col] == 1][[\"fact_longitude\", \"fact_latitude\"]].drop_duplicates()\n    df_not_available = df_train[df_train[col] == 0][[\"fact_longitude\", \"fact_latitude\"]].drop_duplicates()\n\n    fig, ax = plt.subplots(1, 1, figsize=[6, 3], dpi=200)\n    ax.scatter(df_not_available['fact_longitude'], df_not_available['fact_latitude'], \n                s=1, c='tab:red', label=f'not {col}')\n    ax.scatter(df_available['fact_longitude'], df_available['fact_latitude'], \n                s=1, c='tab:green', label=f'{col}')\n\n    ax.legend()\n    plt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:03.039259Z","iopub.execute_input":"2022-05-11T16:46:03.039485Z","iopub.status.idle":"2022-05-11T16:46:07.943916Z","shell.execute_reply.started":"2022-05-11T16:46:03.039452Z","shell.execute_reply":"2022-05-11T16:46:07.943003Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The `wrf`-attributes are missing for all samples at some locations (see red dots). This is not the case for the `cmc`-attributes.","metadata":{}},{"cell_type":"markdown","source":"## Feature engineering","metadata":{}},{"cell_type":"markdown","source":"### Month","metadata":{}},{"cell_type":"code","source":"from datetime import datetime\n\ndef add_month(df):\n    df[\"month\"] = df[[\"fact_time\"]].apply(lambda row: datetime.utcfromtimestamp(row[\"fact_time\"]).month, axis=1)\n    return df\n\ndf_train = add_month(df_train)\ndf_test = add_month(df_test)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:07.945554Z","iopub.execute_input":"2022-05-11T16:46:07.946218Z","iopub.status.idle":"2022-05-11T16:46:44.106215Z","shell.execute_reply.started":"2022-05-11T16:46:07.946168Z","shell.execute_reply":"2022-05-11T16:46:44.105391Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"month_order = [9, 10, 11, 12, 1, 2]\ndf_train[\"month\"].value_counts().loc[month_order]","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:44.107911Z","iopub.execute_input":"2022-05-11T16:46:44.108453Z","iopub.status.idle":"2022-05-11T16:46:44.128786Z","shell.execute_reply.started":"2022-05-11T16:46:44.108406Z","shell.execute_reply":"2022-05-11T16:46:44.128114Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The samples are spread out evenly through the months. We can plot the mean temperature for all data points per month:","metadata":{}},{"cell_type":"code","source":"# group data by month and plot mean fact_temperature\ngrouped = df_train.groupby([\"month\"])[\"fact_temperature\"].mean().loc[month_order]\nfig, ax = plt.subplots(figsize=(10, 3))\nax.plot(np.arange(len(month_order)), grouped.values)\nax.set_xticks(np.arange(len(month_order)))\nax.set_xticklabels(grouped.index)\nax.set_title(\"Mean temperature by month\")\nax.set_xlabel(\"Month of the year\")\nax.set_ylabel(\"Temperature (°C)\")\nax.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:44.131345Z","iopub.execute_input":"2022-05-11T16:46:44.132000Z","iopub.status.idle":"2022-05-11T16:46:44.425707Z","shell.execute_reply.started":"2022-05-11T16:46:44.131966Z","shell.execute_reply":"2022-05-11T16:46:44.424745Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We now encode the month integer values as cyclic features instead, inspired by http://blog.davidkaleko.com/feature-engineering-cyclical-features.html. This preserves the cyclic nature of months (i.e. month 12 is closer to month 1 than to month 9).","metadata":{}},{"cell_type":"code","source":"def add_month_sin_cos(df):\n    df[\"month_sin\"] = np.sin((df[\"month\"] - 1) * (2.0 * np.pi / 12))\n    df[\"month_cos\"] = np.cos((df[\"month\"] - 1) * (2.0 * np.pi / 12))\n    return df\n\ndf_train = add_month_sin_cos(df_train)\ndf_test = add_month_sin_cos(df_test)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:44.427086Z","iopub.execute_input":"2022-05-11T16:46:44.427387Z","iopub.status.idle":"2022-05-11T16:46:44.578695Z","shell.execute_reply.started":"2022-05-11T16:46:44.427347Z","shell.execute_reply":"2022-05-11T16:46:44.577729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can plot the cyclic encoding:","metadata":{}},{"cell_type":"code","source":"fig, ax = plt.subplots(figsize=(10, 3))\nax.scatter(df_train[\"month_sin\"], df_train[\"month_cos\"])\nax.set_title(\"Month cyclic encoding\")\nax.set_xlabel(\"month_sin\")\nax.set_ylabel(\"month_cos\")\nax.margins(0.05)\nax.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:44.579843Z","iopub.execute_input":"2022-05-11T16:46:44.580105Z","iopub.status.idle":"2022-05-11T16:46:57.617005Z","shell.execute_reply.started":"2022-05-11T16:46:44.580075Z","shell.execute_reply":"2022-05-11T16:46:57.615884Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Hour","metadata":{}},{"cell_type":"markdown","source":"Since the timestamps were stored in the UTC timezone, we also need to add or subtract the timezone difference from the sample locations. This is not a big as a problem for months as it is for hour of day.","metadata":{}},{"cell_type":"code","source":"! pip install timezonefinder pytz","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:46:57.618561Z","iopub.execute_input":"2022-05-11T16:46:57.618796Z","iopub.status.idle":"2022-05-11T16:47:13.959562Z","shell.execute_reply.started":"2022-05-11T16:46:57.618769Z","shell.execute_reply":"2022-05-11T16:47:13.958373Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from datetime import datetime\nfrom timezonefinder import TimezoneFinder\nimport pytz\nimport functools\n\ntf = TimezoneFinder()\n\n@functools.lru_cache(maxsize=None)\ndef get_timezone_str(lng, lat):\n    return tf.timezone_at(lng=lng, lat=lat)\n\ndef add_hour(df):\n    def _add_hour(row):\n            longitude, latitude = row[\"fact_longitude\"], row[\"fact_latitude\"]\n            timezone = pytz.timezone(get_timezone_str(longitude, latitude))\n            utc_date = datetime.utcfromtimestamp(row[\"fact_time\"]).replace(tzinfo=pytz.utc)\n            local_date = utc_date.astimezone(timezone)\n            return local_date.hour\n            \n    df[\"hour\"] = df[[\"fact_time\", \"fact_longitude\", \"fact_latitude\"]].apply(_add_hour, axis=1)\n    return df\n\ndf_train = add_hour(df_train)\ndf_test = add_hour(df_test)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:47:13.961257Z","iopub.execute_input":"2022-05-11T16:47:13.961522Z","iopub.status.idle":"2022-05-11T16:48:58.414502Z","shell.execute_reply.started":"2022-05-11T16:47:13.961490Z","shell.execute_reply":"2022-05-11T16:48:58.413552Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"hour_order = list(range(24))\ndf_train[\"hour\"].value_counts().loc[hour_order]","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:48:58.416184Z","iopub.execute_input":"2022-05-11T16:48:58.416436Z","iopub.status.idle":"2022-05-11T16:48:58.437621Z","shell.execute_reply.started":"2022-05-11T16:48:58.416408Z","shell.execute_reply":"2022-05-11T16:48:58.436639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The samples are spread out evenly throughout the day. We can plot the mean temperature for all data points per hour:","metadata":{}},{"cell_type":"code","source":"grouped = df_train.groupby([\"hour\"])[\"fact_temperature\"].mean().loc[hour_order]\nfig, ax = plt.subplots(1, 1, figsize=(10, 3))\nax.plot(np.arange(len(hour_order)), grouped.values)\nax.set_xticks(np.arange(len(hour_order)))\nax.set_xticklabels(grouped.index)\nax.set_title(\"Mean temperature by hour\")\nax.set_xlabel(\"Hour of the day\")\nax.set_ylabel(\"Temperature (°C)\")\nax.tick_params(axis='x', labelsize=7)\nax.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:48:58.439419Z","iopub.execute_input":"2022-05-11T16:48:58.439780Z","iopub.status.idle":"2022-05-11T16:48:58.800059Z","shell.execute_reply.started":"2022-05-11T16:48:58.439731Z","shell.execute_reply":"2022-05-11T16:48:58.799019Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Looks reasonable! We convert these to cyclic values just as we did with months:","metadata":{}},{"cell_type":"code","source":"def add_hour_sin_cos(df):\n    df[\"hour_sin\"] = np.sin((df[\"hour\"]) * (2.0 * np.pi / 24))\n    df[\"hour_cos\"] = np.cos((df[\"hour\"]) * (2.0 * np.pi / 24))\n    return df\n\ndf_train = add_hour_sin_cos(df_train)\ndf_test = add_hour_sin_cos(df_test)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:48:58.801703Z","iopub.execute_input":"2022-05-11T16:48:58.802055Z","iopub.status.idle":"2022-05-11T16:48:58.959558Z","shell.execute_reply.started":"2022-05-11T16:48:58.802006Z","shell.execute_reply":"2022-05-11T16:48:58.958910Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"fig, ax = plt.subplots(1, 1, figsize=(10, 3))\nax.scatter(df_train[\"hour_sin\"], df_train[\"hour_cos\"])\nax.set_title(\"Hour cyclic encoding\")\nax.set_xlabel(\"hour_sin\")\nax.set_ylabel(\"hour_cos\")\nax.margins(0.05)\nax.grid()\nplt.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:48:58.960860Z","iopub.execute_input":"2022-05-11T16:48:58.961616Z","iopub.status.idle":"2022-05-11T16:49:12.191794Z","shell.execute_reply.started":"2022-05-11T16:48:58.961581Z","shell.execute_reply":"2022-05-11T16:49:12.190924Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Imputing missing values","metadata":{}},{"cell_type":"markdown","source":"We first create methods for getting all data points from the closest weather station of a given data point:","metadata":{}},{"cell_type":"code","source":"from sklearn.neighbors import BallTree\n\n@functools.lru_cache(maxsize=None)\ndef get_ball_tree(dataset):\n    df = globals()[f\"df_{dataset}\"]\n    return BallTree(\n        np.deg2rad(df[[\"fact_latitude\", \"fact_longitude\"]].values), metric=\"haversine\"\n    )\n\ndef get_nearest_neighbors(dataset, lat, lng):  ## dataset = \"test\" or \"train\"\n    df = globals()[f\"df_{dataset}\"]\n    tree = get_ball_tree(dataset)\n    indices = tree.query([[np.deg2rad(lat), np.deg2rad(lng)]], k=1, return_distance=False)\n    coords = df.iloc[indices.flatten()][[\"fact_latitude\", \"fact_longitude\"]]\n    return df[df[\"fact_longitude\"].isin(coords[\"fact_longitude\"]) & df[\"fact_latitude\"].isin(coords[\"fact_latitude\"])]","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:49:12.193457Z","iopub.execute_input":"2022-05-11T16:49:12.194007Z","iopub.status.idle":"2022-05-11T16:49:13.123588Z","shell.execute_reply.started":"2022-05-11T16:49:12.193957Z","shell.execute_reply":"2022-05-11T16:49:13.122627Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Testing for a datapoint in Australia:","metadata":{}},{"cell_type":"code","source":"get_nearest_neighbors(\"train\", -28.149833623103504, 139.14635543626596)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:49:13.125671Z","iopub.execute_input":"2022-05-11T16:49:13.125911Z","iopub.status.idle":"2022-05-11T16:49:21.194220Z","shell.execute_reply.started":"2022-05-11T16:49:13.125883Z","shell.execute_reply":"2022-05-11T16:49:21.193145Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"**NOTE: Kaggle notebook editor is not capable of running from here and onwards.**","metadata":{}},{"cell_type":"markdown","source":"We can then impute the missing values based on the methods above. For a given data sample missing some values, we find first find the closest data points in distance. We then filter these on month and hour, before taking the mean of the relevant columns to use for filling the gaps in the query data point:","metadata":{}},{"cell_type":"code","source":"def impute_missing_values(dataset):   ## dataset = \"test\" or \"train\"\n    df = globals()[f\"df_{dataset}\"]\n    for column_category in [\"gfs\", \"cmc\", \"wrf\"]:\n        rows_with_na = df[df[f\"{column_category}_available\"] == 0]\n        print(f\"Imputing missing values for {len(rows_with_na)} rows in {column_category}\")\n        i = 0\n        for index, row in rows_with_na.iterrows():\n            nearest_neighbors = get_nearest_neighbors(\n                dataset, row[\"fact_latitude\"], row[\"fact_longitude\"]\n            )\n            nearest_neighbors = nearest_neighbors[\n                (nearest_neighbors[\"month_sin\"] == row[\"month_sin\"])\n                & (nearest_neighbors[\"month_cos\"] == row[\"month_cos\"])\n            ]\n            nearest_neighbors = nearest_neighbors[\n                (np.abs(nearest_neighbors[\"hour_sin\"] - row[\"hour_sin\"]) <= 0.5) &\n                (np.abs(nearest_neighbors[\"hour_cos\"] - row[\"hour_cos\"]) <= 0.5)\n            ]\n            column_means = nearest_neighbors[column_categories[column_category]].mean()\n            for column in column_categories[column_category]:\n                df.iloc[index, df.columns.get_loc(column)] = column_means[column]\n            i += 1\n            if i % 1000 == 0:\n                print(f\"{i} rows processed\")\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# NOTE: Took ~20 minutes to run locally.\nimpute_missing_values(\"train\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# df_train.to_csv(\"/kaggle/working/train_imputed.csv\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"impute_missing_values(\"test\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# df_test.to_csv(\"/kaggle/working/test_imputed.csv\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Dropping features and rest of missing data rows","metadata":{}},{"cell_type":"markdown","source":"We drop features that are either not relevant for training or not expected to be useful for predicting the temperature:","metadata":{}},{"cell_type":"code","source":"columns_to_drop = [\n    \"index\", # not relevant for training\n    \"fact_time\", # we use cyclic features instead\n    \"cmc_available\", # not needed after imputation\n    \"gfs_available\", # not needed after imputation\n    \"gfs_soil_temperature_available\", # not needed after imputation\n    \"wrf_available\", # not needed after imputation\n    \"cmc_timedelta_s\", # we use cyclic features instead\n    \"gfs_timedelta_s\", # we use cyclic features instead\n    \"hour\",  # we use cyclic features instead\n    \"month\",  # we use cyclic features instead\n]\n\n\nselected_df_train = df_train.drop(columns_to_drop, axis=1)\nselected_df_test = df_test.drop(columns_to_drop, axis=1)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We still have some rows with missing values left due to the imputation above not finding any relevant neighbouring data to impute from. We drop these ~25 000 rows before continuing.","metadata":{}},{"cell_type":"code","source":"selected_df_train = selected_df_train.dropna()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Preparing data for training","metadata":{}},{"cell_type":"markdown","source":"We first separate the training features (`X`) from the target feature (`y`), then we split the training data set into two new sets: one for training and one for validation.","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import train_test_split\n\nX = selected_df_train.drop(\"fact_temperature\", axis=1)\ny = selected_df_train[\"fact_temperature\"]\n\nX_train, X_val, y_train, y_val = train_test_split(X, y, random_state=42, test_size=10000)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Hyperparameter optimization","metadata":{}},{"cell_type":"markdown","source":"We use random grid search to search through 10 random samples of hyperparameters for the LightGBM model. The values in the param grid is inspired by https://lightgbm.readthedocs.io/en/latest/Parameters-Tuning.html and https://towardsdatascience.com/kagglers-guide-to-lightgbm-hyperparameter-tuning-with-optuna-in-2021-ed048d9838b5.\n\nSome parameters affect the performance directly, some affects the training time and some are added to reduce overfitting.","metadata":{}},{"cell_type":"code","source":"from lightgbm import LGBMRegressor\nfrom sklearn.model_selection import RandomizedSearchCV\n\n\nmodel = LGBMRegressor(\n    random_state=42,\n    n_jobs=4,\n    metric=\"rmse\",\n)\n\nparam_grid = {\n    \"model__num_iterations\": [500, 1000],\n    \"model__num_leaves\": list(range(10, 3000, 10)),\n    \"model__learning_rate\": [0.01, 0.05, 0.1, 0.3],\n    \"model__min_data_in_leaf\": [100, 500, 1000],\n    \"model__max_depth\": [3, 7, 12],\n    \"model__lambda_l1\": list(range(0, 100, 5)),\n    \"model__lambda_l2\": list(range(0, 100, 5)),\n    \"model__min_gain_to_split\": list(range(15)),\n}\n\n\npipeline = Pipeline(\n    [\n        (\"standardization\", StandardScaler()),\n        (\"model\", model),\n    ]\n)\n\nsearch = RandomizedSearchCV(\n    pipeline,\n    param_grid,\n    n_iter=10,\n    cv=5,\n    verbose=2,\n    random_state=42,\n    scoring=\"neg_root_mean_squared_error\",\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The `search.fit` is commented out here due to the long running time (~3 hours when running locally).","metadata":{}},{"cell_type":"code","source":"# search.fit(X_train, y_train)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Commented out because of the same reason as above. The output is included in markdown.","metadata":{}},{"cell_type":"code","source":"# search.best_params_","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Output:\n```python\n{\n    'model__num_leaves': 2320,\n    'model__num_iterations': 1000,\n    'model__min_gain_to_split': 4,\n    'model__min_data_in_leaf': 1000,\n    'model__max_depth': 7,\n    'model__learning_rate': 0.3,\n    'model__lambda_l2': 80,\n    'model__lambda_l1': 85\n}\n```","metadata":{}},{"cell_type":"markdown","source":"## Model training and evaluation of predictions","metadata":{}},{"cell_type":"markdown","source":"We know train the model with the optimal parameters found above. Note that we double the number of iterations (this was at max 1000 during grid search to reduce fit time), and `num_leaves` were reduced to be less than `2^(max_depth)` as suggested by [LightGBM docs](https://lightgbm.readthedocs.io/en/latest/Parameters-Tuning.html).","metadata":{}},{"cell_type":"code","source":"model = LGBMRegressor(\n    random_state=42,\n    n_jobs=4,\n    metric=\"rmse\",\n    num_iterations=2000,\n    num_leaves=100,\n    max_depth=7,\n    min_gain_to_split=4,\n    learning_rate=0.3,\n    lambda_l1=85,\n    lambda_l2=80,\n    min_data_in_leaf=1000,\n)\n\npipeline = Pipeline(\n    [\n        (\"standardization\", StandardScaler()),\n        (\"model\", model),\n    ]\n)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pipeline.fit(X_train, y_train)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Training set prediction and error","metadata":{}},{"cell_type":"markdown","source":"We can now try to predict the training data:","metadata":{}},{"cell_type":"code","source":"y_train_pred = pipeline.predict(X_train)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from sklearn.metrics import mean_squared_error\n\n# training set error\nmean_squared_error(y_train, y_train_pred, squared=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Output when running locally: 1.651446211418595","metadata":{}},{"cell_type":"markdown","source":"### Validation set error","metadata":{}},{"cell_type":"markdown","source":"Then we predict the validation data. We expect the error to be a bit higher here.","metadata":{}},{"cell_type":"code","source":"y_val_pred = pipeline.predict(X_val)\nmean_squared_error(y_val, y_val_pred, squared=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Output when running locally: 1.7684266867801408","metadata":{}},{"cell_type":"markdown","source":"## Creating submission file","metadata":{}},{"cell_type":"markdown","source":"Finally, we predict the testing data and create a submission file for Kaggle:","metadata":{}},{"cell_type":"code","source":"y_test_pred = pipeline.predict(selected_df_test)\nsubmission_df = pd.DataFrame(df_test[\"index\"])\nsubmission_df[\"fact_temperature\"] = y_test_pred\nsubmission_df.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_df.to_csv(\"/kaggle/working/submission.csv\", index=False)\n! head -6 \"/kaggle/working/submission.csv\"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## SMOTE (not used)","metadata":{}},{"cell_type":"markdown","source":"The following code was used to generate data points. It is not used in our final model as the benchmark performance did not improve with the generated training data (sample = 10 000).","metadata":{}},{"cell_type":"code","source":"'''\nfrom global_land_mask import globe\nimport pandas as pd \nimport numpy as np \nfrom sklearn.model_selection import train_test_split\nimport random\n\n\"\"\"\n    function euclideanDistance: Finding the distance between two coordinate points.\n    param point1:               Point 1.\n    param point2:               Point 2.\n    return distance:            The euclidean distance between two points.\n\"\"\"\ndef euclideanDistance(point1,point2):\n\n    distance = np.linalg.norm(point1 - point2)\n    return distance\n\n\"\"\"\n    function generateNRandomHours:  Generate N random hours.\n    param N:                        Number of generated hours.\n    return hours:                   Numpy array of generated hours.\n\"\"\"\ndef generateNRandomHours(N=10):\n    hours = np.array([])\n    i=0\n    for i in range(N):\n        hours = np.append(hours, np.random.randint(24))\n    return hours\n\n\"\"\"\n    function generateNRandomMonths: Generate N random months.\n    param N:                        Number of generated months.\n    return months:                  Numpy array of generated months.\n\"\"\"\ndef generateNRandomMonths(N=10):\n    monthList = [9, 10, 11, 12, 1, 2]\n    months = np.array([])\n    i=0\n    for i in range(N):\n        months = np.append(months, monthList[np.random.randint(len(monthList))])\n        i+=1\n    return months\n\n\"\"\"\n    function generateNRandomCoordinates:    Generate N random valid coordinates.\n    param lon_min:                          Minimum value of longitude generated.\n    param lon_max:                          Maximum value of longitude generated.\n    param lat_min:                          Minimum value of latitude generated.\n    param lat_max:                          Maximum value of latitude generated.\n    param N:                                Number of generated coordinate points.    \n    return hours:                           Numpy array of generated coordinate points.\n\"\"\"\ndef generateNRandomCoordinates(lon_min = -100, lon_max = 150, lat_min=-20, lat_max=20, N=10):\n    coordinates = np.array([])\n    i=0\n    while i<N:\n        lon = round(random.uniform(lon_min, lon_max), 6)\n        lat = round(random.uniform(lat_min, lat_max), 6)\n        if globe.is_land(lat, lon):\n            coordinates = np.append(coordinates, np.array((lat, lon)))\n            i+=1\n            #print(\"Appending:\", np.array((lat, lon)))\n    coordinates = coordinates.reshape(-1, 2)\n    return coordinates     \n\n\"\"\"\n    function saveData:  Save dataframe as csv.\n    param data:         Dataframe to be saved as csv.\n    return 0:           No return.\n\"\"\"\ndef saveData(data):\n    data.to_csv('public/generated_trainV2.csv', index=False, mode='a', header=False) \n    print(\"Successfully saved data.\")\n    return 0   \n\n\"\"\"\n    function generateData:      Generate N data points.\n    param train:                Training data to be used for neighbour sampling.\n    param min_lat:              Minimum value of latitude of data point.\n    param max_lat:              Maximum value of latitude of data point.\n    param N:                    Number of generated data points.\n    param K:                    Number of neighbours.\n    return df_new_data:         Dataframe with N generated data points.\n\"\"\"\ndef generateData(train, min_lat, max_lat, N=10, K=3):\n    coordinates = generateNRandomCoordinates(-155, 155, min_lat, max_lat, N)\n    months = generateNRandomMonths(N)\n    hours = generateNRandomHours(N)\n    print(\"Created \", N, \" new point(s)\")\n    df_new_data = pd.DataFrame()\n    for i in range(N):\n        distCoord = []\n        for count in range(len(train)):\n            distCoord.append(euclideanDistance(coordinates[i], np.array((train.iloc[count]['fact_latitude'], train.iloc[count]['fact_longitude']))))\n            #print(\"Coordinate loop: \", count)\n            \n        sortedDistCoord = np.argsort(distCoord)[:(K*99)]\n        df_placeholder = train.iloc[sortedDistCoord]\n        distMonth = []\n        month = np.array((np.sin((months[i]-1) * (2.0 * np.pi / 12)), np.cos((months[i]-1) * (2.0 * np.pi / 12))))\n        for count2 in range(len(sortedDistCoord)):\n            distMonth.append(euclideanDistance(month, np.array((df_placeholder.iloc[count2]['month_sin'], df_placeholder.iloc[count2]['month_cos']))))\n            #print(\"Month loop: \", count2)\n\n        sortedDistMonth = np.argsort(distMonth)[:(K*49)]\n        df_placeholder = df_placeholder.iloc[sortedDistMonth]\n        distHour = []\n        hour = np.array((np.sin(hours[i] * (2.0 * np.pi / 24)), np.cos(hours[i] * (2.0 * np.pi / 24))))\n        for count3 in range(len(sortedDistMonth)):\n            distHour.append(euclideanDistance(hour, np.array((df_placeholder.iloc[count3]['hour_sin'], df_placeholder.iloc[count3]['hour_cos']))))\n            #print(\"Hour loop: \", count3)\n\n        sortedDistFinal = np.argsort(distHour)[:K]\n        df_placeholder = df_placeholder.iloc[sortedDistFinal].mean(axis=0).to_frame().T\n\n        df_placeholder['fact_latitude'] = coordinates[i][0]\n        df_placeholder['fact_longitude'] = coordinates[i][1]\n        df_placeholder['month_sin'] = np.sin((months[i] - 1) * (2.0 * np.pi / 12))\n        df_placeholder['month_cos'] = np.cos((months[i] - 1) * (2.0 * np.pi / 12))\n        df_placeholder['hour_sin'] = np.sin(hours[i] * (2.0 * np.pi / 24))\n        df_placeholder['hour_cos'] = np.cos(hours[i] * (2.0 * np.pi / 24))\n\n        df_new_data = pd.concat([df_new_data, df_placeholder])\n        if i%10 == 0:\n            print(\"New data generated: \", i+1)\n\n    return df_new_data\n\ndf_train = pd.read_csv('public/train_imputed.csv', low_memory=True)\ndf_test = pd.read_csv('public/test_imputed.csv', low_memory=True)\ngen_train = pd.read_csv('public/generated_trainV2.csv', low_memory=True)\n\ncolumns_to_drop = [\n    \"Unnamed: 0\",\n    \"index\",\n    \"fact_time\",\n    \"cmc_available\", \n    \"gfs_available\",\n    \"gfs_soil_temperature_available\", \n    \"wrf_available\",\n    \"cmc_timedelta_s\",\n    \"gfs_timedelta_s\",\n]\n\nnew_df_train = df_train.drop(columns_to_drop, axis=1)\nnew_df_train = new_df_train.dropna()\nnew_df_test = df_test.drop(columns_to_drop, axis=1)\n\nX = new_df_train.drop('fact_temperature', axis=1)\ny = new_df_train[\"fact_temperature\"]\n\nX_train, X_val, y_train, y_val = train_test_split(X, y, random_state=42, test_size=10000)\n\nSMOTE_train = pd.concat([new_df_train, gen_train])\n\n# Subsamples of training data in order to make the data generation faster.\nSMOTE_train1 = SMOTE_train.loc[(SMOTE_train['fact_latitude'] >= -15) & (SMOTE_train['fact_latitude'] <= 0)]\nSMOTE_train2 = SMOTE_train.loc[(SMOTE_train['fact_latitude'] >= 0) & (SMOTE_train['fact_latitude'] <= 15)]\nSMOTE_train3 = SMOTE_train.loc[(SMOTE_train['fact_latitude'] >= -25) & (SMOTE_train['fact_latitude'] <= 25)]\n\n# Generator of new data\ntry:\n    while True:\n        test = generateData(SMOTE_train1, -15, 0, 100, 5)\n        saveData(test)\n        test = generateData(SMOTE_train2, 0, 15, 100, 5)\n        saveData(test)\n        test = generateData(SMOTE_train3, -25, 25, 50, 5)\n        saveData(test)\nexcept KeyboardInterrupt:\n    pass\n'''","metadata":{},"execution_count":null,"outputs":[]}]}