{"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":"# Group4 AML Challenge 1 Notebook\n\nGeneral thoughts: (Thoughts that are not yet assigned to a section)\n- Model validation:\n    - Dataset is probably too large for k-fold-cross-validation (Philipp)\n    - We should put some regions only in the validation set and not the training set to have a good indicator on generalization performance. (Philipp)\n\n- Performance:\n    - We can quickly try out models with a small proxy subset (Philipp)\n- Other:\n    - We should balance the longitudes in terms of number samples\n","metadata":{}},{"cell_type":"code","source":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-05-11T15:54:20.57063Z","iopub.execute_input":"2022-05-11T15:54:20.57091Z","iopub.status.idle":"2022-05-11T15:54:20.584664Z","shell.execute_reply.started":"2022-05-11T15:54:20.570878Z","shell.execute_reply":"2022-05-11T15:54:20.583795Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Load data\ndf = pd.read_csv('../input/eurecom-aml-2022-challenge-1/public/train.csv',low_memory = True, dtype='float32')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:54:20.586853Z","iopub.execute_input":"2022-05-11T15:54:20.587422Z","iopub.status.idle":"2022-05-11T15:55:31.580954Z","shell.execute_reply.started":"2022-05-11T15:54:20.587379Z","shell.execute_reply":"2022-05-11T15:55:31.580167Z"},"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, dtype='float32')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:55:31.582476Z","iopub.execute_input":"2022-05-11T15:55:31.582767Z","iopub.status.idle":"2022-05-11T15:55:53.258506Z","shell.execute_reply.started":"2022-05-11T15:55:31.582727Z","shell.execute_reply":"2022-05-11T15:55:53.257664Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Create smaller proxy subsets 'df_small_[magnitude]' with 10^magnitude random samples\nfor mag in [4, 5]:\n    name = f'df_small_{mag}'\n    globals()[name] = df.sample(n=10**mag, random_state=1, ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:55:53.260551Z","iopub.execute_input":"2022-05-11T15:55:53.261357Z","iopub.status.idle":"2022-05-11T15:55:53.601007Z","shell.execute_reply.started":"2022-05-11T15:55:53.261314Z","shell.execute_reply":"2022-05-11T15:55:53.600239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data analysis\nThe purpose of the section is to explore the data:\n- Get insight into the data\n- Find potential problems for models\n- Find biases in the data\n\nTODOs:\n- Do we have outliers and in which features?\n- How many 'null' values are there and in which features?\n\nGiven description:\\\n\"\\\nIn this challenge, you are free (and encouraged) to explore in depth the data you have, you can run simple queries on the data, perform exploration and compute statistics.\n\n**NOTE**: finding the right question to ask is difficult! Don't be afraid to complement your analysis with your own questions. This can give you extra points!\n\n**NOTE 2**: the presentation quality is critical in any business-oriented data analysis. Take time to create few but informative plots, rather than endless tables!\\\n\"","metadata":{}},{"cell_type":"markdown","source":"### General Properties of the Data\nMotivation: Obvious\n\nResults:\n- ...","metadata":{}},{"cell_type":"code","source":"# df.describe() # takes a long time","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:55:53.602207Z","iopub.execute_input":"2022-05-11T15:55:53.604182Z","iopub.status.idle":"2022-05-11T15:55:53.607797Z","shell.execute_reply.started":"2022-05-11T15:55:53.604144Z","shell.execute_reply":"2022-05-11T15:55:53.60704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from matplotlib import pyplot as plt\nfig = plt.figure(figsize = (30,40))\nax = fig.gca()\ndf_small_5.hist(ax = ax)\npass","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:55:53.609387Z","iopub.execute_input":"2022-05-11T15:55:53.609924Z","iopub.status.idle":"2022-05-11T15:56:11.222031Z","shell.execute_reply.started":"2022-05-11T15:55:53.609881Z","shell.execute_reply":"2022-05-11T15:56:11.221259Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Analyse NaN","metadata":{}},{"cell_type":"code","source":"\n# Create a dataframe that includes the number of data points for each location\nnan_df = df[df.isnull().any(axis=1)]\nnan_location_count = nan_df.groupby(['fact_latitude', 'fact_longitude']).size()\n# Convert from series to dataframe and rename column to 'count'\nnan_location_count = nan_location_count.reset_index().rename({0:'count'}, axis=1)\n# Library reference for this code: https://plotly.com/python/mapbox-layers/\nimport plotly.express as px\n\nfig = px.scatter_mapbox(nan_location_count, lat=\"fact_latitude\", lon=\"fact_longitude\", hover_name='count',\n                        color='count', hover_data=[\"fact_latitude\", \"fact_longitude\"],\n                        color_discrete_sequence=[\"fuchsia\"], zoom=1, height=300)\nfig.update_layout(mapbox_style=\"open-street-map\")\nfig.update_layout(margin={\"r\":0,\"t\":0,\"l\":0,\"b\":0})\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:11.223148Z","iopub.execute_input":"2022-05-11T15:56:11.223409Z","iopub.status.idle":"2022-05-11T15:56:15.751521Z","shell.execute_reply.started":"2022-05-11T15:56:11.223375Z","shell.execute_reply":"2022-05-11T15:56:15.749656Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}},{"cell_type":"markdown","source":"### Plots Between fact_temp and other features \nMotivation : Make simple ploting functions to observe how our target feature (fact_temp) evolves wrt other features\n\n**functions**:  \n1.      **plot_tempeturs_in_function_of(column, df, color, label):**\n            Does scatter plots between the temp and any other feature.  \n2.      **plot_with_longitudes_steps(longitude_step, df, column):**\n            Iterats trough a df making plots for several longitudes segments. longitude_step indicates the size of each segments ie the numbers of plots done \n        \nResults: \n>        some features like cloud coverage give very noisy plots that show little corelation with fact_temp, whereas others like climate_temperature show easy to follow tendencies.\n\n> An idea could be to try and train simple perceptrons to define upper and lower bounds within which we can make more precise regressions (dont know if that makes good sense)","metadata":{}},{"cell_type":"code","source":"df_small_4.isnull().sum().sum()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:15.753245Z","iopub.execute_input":"2022-05-11T15:56:15.753837Z","iopub.status.idle":"2022-05-11T15:56:15.768874Z","shell.execute_reply.started":"2022-05-11T15:56:15.753783Z","shell.execute_reply":"2022-05-11T15:56:15.768077Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import matplotlib.pyplot as plt\n\ndef plot_tempeturs_in_function_of(column, df, color, label):\n    df_sort = df.sort_values(column)\n    plt.figure(figsize=(8, 6), dpi=80)\n    Y = []\n    X = []\n    for index, row in df.iterrows():\n        Y.append(np.float64(row.fact_temperature))\n        X.append(np.float64(row[column]))\n    # plot point clouds \n    plt.scatter(X,Y,c = color, label = label)\n    plt.legend()\n    plt.title(f\"tempeturs_in_function_of_{str(column)}\")\n    plt.xlabel(column)\n    plt.ylabel(\"fact_temperature\")","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:15.770506Z","iopub.execute_input":"2022-05-11T15:56:15.771051Z","iopub.status.idle":"2022-05-11T15:56:15.779919Z","shell.execute_reply.started":"2022-05-11T15:56:15.77101Z","shell.execute_reply":"2022-05-11T15:56:15.77904Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def plot_with_longitudes_steps(longitude_step, df, column):\n    df.sort_values('fact_longitude')\n    color = iter(plt.cm.rainbow(np.linspace(0, 1, longitude_step)))\n    longitude_start = -180\n    plt.figure(figsize=(8, 6), dpi=80)\n    Y = []\n    X = []\n    while(longitude_start < 180):\n        df_longitude_step = df.where(df.fact_longitude <= longitude_start + longitude_step)\n        df_longitude_step = df_longitude_step.where(df.fact_longitude >= longitude_start)\n        c = next(color)\n        plot_tempeturs_in_function_of(column, df_longitude_step, c, f'{str(longitude_start)}/{str(longitude_start + longitude_step)}')\n        longitude_start += longitude_step","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:15.786468Z","iopub.execute_input":"2022-05-11T15:56:15.786937Z","iopub.status.idle":"2022-05-11T15:56:15.794963Z","shell.execute_reply.started":"2022-05-11T15:56:15.786901Z","shell.execute_reply":"2022-05-11T15:56:15.794198Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_count_nulls = df.isnull().sum()\ndf_disp_null_count = df_count_nulls.loc[(df_count_nulls != 0)]\ndf_disp_null_count","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:15.79641Z","iopub.execute_input":"2022-05-11T15:56:15.797007Z","iopub.status.idle":"2022-05-11T15:56:16.279749Z","shell.execute_reply.started":"2022-05-11T15:56:15.796953Z","shell.execute_reply":"2022-05-11T15:56:16.279026Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#plot_with_longitudes_steps(30, df_small, 'climate_temperature')\n#plot_with_longitudes_steps(30,df_small_4,'wrf_t2_next')\n#plot_tempeturs_in_function_of('wrf_t2_next', df_small_4, 'r', '-180/180')\n#plot_tempeturs_in_function_of('climate_temperature', df_small, 'r', '-180/180')\n#plot_tempeturs_in_function_of('gfs_cloudness', df_small, 'r', '-180/180')\n#plot_tempeturs_in_function_of('gfs_humidity', df_small, 'r', '-180/180')\n#plot_tempeturs_in_function_of('gfs_precipitations', df_small, 'r', '-180/180')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:16.28112Z","iopub.execute_input":"2022-05-11T15:56:16.281591Z","iopub.status.idle":"2022-05-11T15:56:16.28609Z","shell.execute_reply.started":"2022-05-11T15:56:16.281549Z","shell.execute_reply":"2022-05-11T15:56:16.284993Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Correlation Between fact_temperature and other features\nMotivation: See how features correlate with target feature\\\nResults:\n- Features with the strongest correlation to the target/fact_temperature are:\n    1. wrf_t2_interpolated:      0.963\n    2. wrf_t2_next:              0.958\n    3. gfs_temperature_97500:    0.874\n    4. gfs_temperature_95000:    0.863","metadata":{}},{"cell_type":"code","source":"corr_small_5 = df_small_5.corr(method ='pearson')\nbest_features = list(corr_small_5['fact_temperature'].sort_values(ascending=False).head(60).index)\nbest_features","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:16.288008Z","iopub.execute_input":"2022-05-11T15:56:16.288552Z","iopub.status.idle":"2022-05-11T15:56:19.887686Z","shell.execute_reply.started":"2022-05-11T15:56:16.288509Z","shell.execute_reply":"2022-05-11T15:56:19.886888Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Analyse temperature distribution\nMotivation: Find representation of temperature ranges as well es outliers\nResults:\n- Majority of our temperature points are between 0 and 30 C,\n- Few strange outliers. Is t = 60 C real??\n\n\n**Note:** Consider plotting the temperatures during the day and night independently, e.g., 11:00-14:00 and 2:00-6:00 (Philipp)\n","metadata":{}},{"cell_type":"code","source":"temps = df['fact_temperature']\ntemps.hist(bins=60)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:19.889124Z","iopub.execute_input":"2022-05-11T15:56:19.890156Z","iopub.status.idle":"2022-05-11T15:56:20.293651Z","shell.execute_reply.started":"2022-05-11T15:56:19.89011Z","shell.execute_reply":"2022-05-11T15:56:20.292875Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Show data point density\nMotivation: See how dense the samples are in the different areas in the training and test data. Also see how many data points we have in the different regions.\n\nResults:\n- most data points in nothern hemisphere\n- few data points near equator and towards the poles","metadata":{"execution":{"iopub.status.busy":"2022-04-29T13:07:03.720198Z","iopub.execute_input":"2022-04-29T13:07:03.721134Z","iopub.status.idle":"2022-04-29T13:07:04.265567Z","shell.execute_reply.started":"2022-04-29T13:07:03.72107Z","shell.execute_reply":"2022-04-29T13:07:04.264574Z"}}},{"cell_type":"code","source":"# Create a dataframe that includes the number of data points for each location\ndf_location_count = df.groupby(['fact_latitude', 'fact_longitude']).size()\n# Convert from series to dataframe and rename column to 'count'\ndf_location_count = df_location_count.to_frame().reset_index().rename({0:'count'}, axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:20.295277Z","iopub.execute_input":"2022-05-11T15:56:20.295585Z","iopub.status.idle":"2022-05-11T15:56:20.46233Z","shell.execute_reply.started":"2022-05-11T15:56:20.295544Z","shell.execute_reply":"2022-05-11T15:56:20.461505Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import plotly.express as px\nfig = px.density_mapbox(df_location_count, lat='fact_latitude', lon='fact_longitude', z='count', radius=5,\n                        center=dict(lat=0, lon=180), zoom=0, title = \"Heatmap <br><sup>Heatmap of data point density (does not take number of data points per location into account)</sup>\",\n                        mapbox_style=\"stamen-terrain\")\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:20.463547Z","iopub.execute_input":"2022-05-11T15:56:20.464315Z","iopub.status.idle":"2022-05-11T15:56:20.543307Z","shell.execute_reply.started":"2022-05-11T15:56:20.464272Z","shell.execute_reply":"2022-05-11T15:56:20.5424Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Show data point count per location\nMotivation: See how many data points we have in the different measure locations\n\nResults:\n- Total of 8053 different locations\n- Most locations have below 20 data points\n- Locations in the US have lots of data points per location (around 600 for most)","metadata":{"execution":{"iopub.status.busy":"2022-04-29T13:07:03.720198Z","iopub.execute_input":"2022-04-29T13:07:03.721134Z","iopub.status.idle":"2022-04-29T13:07:04.265567Z","shell.execute_reply.started":"2022-04-29T13:07:03.72107Z","shell.execute_reply":"2022-04-29T13:07:04.264574Z"}}},{"cell_type":"code","source":"df_location_count['count'].hist(bins=60)\ndf_location_count['count'].describe()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:20.544784Z","iopub.execute_input":"2022-05-11T15:56:20.545559Z","iopub.status.idle":"2022-05-11T15:56:20.898693Z","shell.execute_reply.started":"2022-05-11T15:56:20.545516Z","shell.execute_reply":"2022-05-11T15:56:20.897897Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Library reference for this code: https://plotly.com/python/mapbox-layers/\nimport plotly.express as px\n\nfig = px.scatter_mapbox(df_location_count, lat=\"fact_latitude\", lon=\"fact_longitude\", hover_name='count',\n                        color='count', hover_data=[\"fact_latitude\", \"fact_longitude\"],\n                        color_discrete_sequence=[\"fuchsia\"], zoom=1, height=300)\nfig.update_layout(mapbox_style=\"open-street-map\")\nfig.update_layout(margin={\"r\":0,\"t\":0,\"l\":0,\"b\":0})\nfig.show()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:20.900233Z","iopub.execute_input":"2022-05-11T15:56:20.900677Z","iopub.status.idle":"2022-05-11T15:56:21.40164Z","shell.execute_reply.started":"2022-05-11T15:56:20.900633Z","shell.execute_reply":"2022-05-11T15:56:21.400903Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"###  Latitude Data point distribution\nMotivation: See which longitudes are over and under represented\n\nResults:\n- Most data points are between the latitude 40 - 60\n\nTODOs:\n- Integrate test data","metadata":{}},{"cell_type":"code","source":"df['fact_latitude'].hist(bins=60)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:21.403255Z","iopub.execute_input":"2022-05-11T15:56:21.403768Z","iopub.status.idle":"2022-05-11T15:56:21.790712Z","shell.execute_reply.started":"2022-05-11T15:56:21.403728Z","shell.execute_reply":"2022-05-11T15:56:21.789878Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Data Pre-processing/Feature Engineering\nIn this section we process the data by doing for example the following things:\n-    Normalising and standardising the given data;\n-    Removing outliers;\n-    Carrying out feature selection, possibly using metrics derived from information theory;\n-    Handling missing information in the dataset;\n-    Augmenting the dataset with external information;\n-    Combining existing features.\n\nIdeas:\n- [x] Hour of the day to distinguish between day and night\n- [ ] Temperature forecast for closest (& above & below) forecast height\n- [ ] Interpolated between closest two temperature forecast heights depending on distance","metadata":{"execution":{"iopub.status.busy":"2022-04-29T12:13:00.59631Z","iopub.execute_input":"2022-04-29T12:13:00.596602Z","iopub.status.idle":"2022-04-29T12:13:13.120938Z","shell.execute_reply.started":"2022-04-29T12:13:00.596572Z","shell.execute_reply":"2022-04-29T12:13:13.119755Z"}}},{"cell_type":"markdown","source":"### Feature: Elevation","metadata":{}},{"cell_type":"code","source":"def add_elevation(df_, filename='../input/elevation-files/elevation_train.csv'):\n    df_height = pd.read_csv(filename, index_col='Unnamed: 0')\n    df_['elevation'] = df_height\nadd_elevation(df)\nadd_elevation(df_test, filename='../input/elevation-files/elevation_test.csv')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:21.792335Z","iopub.execute_input":"2022-05-11T15:56:21.792876Z","iopub.status.idle":"2022-05-11T15:56:22.891123Z","shell.execute_reply.started":"2022-05-11T15:56:21.792832Z","shell.execute_reply":"2022-05-11T15:56:22.890275Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_elevation(df_, name='train'):\n    !pip install SRTM.py\n    import srtm\n    from tqdm import tqdm\n    location_count = df_.groupby(['fact_latitude', 'fact_longitude'])\n    elevation_data = srtm.get_data(local_cache_dir=\"/kaggle/working/geo_cache\", srtm1=False, srtm3=True)\n    for group_id in tqdm(location_count.groups):\n        height = elevation_data.get_elevation(*group_id, approximate=False)\n        df_.loc[location_count.groups[group_id], 'elevation'] = height\n    df_['elevation'].to_csv(f'elevation_{name}.csv')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:22.896241Z","iopub.execute_input":"2022-05-11T15:56:22.898728Z","iopub.status.idle":"2022-05-11T15:56:22.912471Z","shell.execute_reply.started":"2022-05-11T15:56:22.898674Z","shell.execute_reply":"2022-05-11T15:56:22.911148Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# get_elevation(df_test, name='test')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:22.913783Z","iopub.execute_input":"2022-05-11T15:56:22.914224Z","iopub.status.idle":"2022-05-11T15:56:22.925002Z","shell.execute_reply.started":"2022-05-11T15:56:22.914179Z","shell.execute_reply":"2022-05-11T15:56:22.924139Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Feature: Approximate Hour of the Day: 'approx_hour'\nMotivation: Help classifier to distinguish between day and night\n\nNot exact hour because then we would need take look up the time zone of every location. \nInstead, time zone offset is approximated based on the longitude.\\\nGreenwhich is longitude 0 and time zone UTC -> +-180 longitude is 12 hours of time offset with linear interpolation in between.","metadata":{}},{"cell_type":"code","source":"def add_hour(df_):\n    date_time = pd.to_datetime(df_['fact_time'], unit='s')\n    hour = date_time.dt.hour\n    minute = date_time.dt.minute\n    # Unix time does not take the time zone into account -> add offset based on longitude\n    longitude_time_offset = (df_['fact_longitude'] * 12 / 180)\n\n    df_['approx_hour'] = (hour + minute/60 + longitude_time_offset)%24\n    return df_","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:22.92904Z","iopub.execute_input":"2022-05-11T15:56:22.931151Z","iopub.status.idle":"2022-05-11T15:56:22.942175Z","shell.execute_reply.started":"2022-05-11T15:56:22.931097Z","shell.execute_reply":"2022-05-11T15:56:22.941316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_hour_floor(df_):\n    date_time = pd.to_datetime(df_['fact_time'], unit='s')\n    hour = date_time.dt.hour\n    minute = date_time.dt.minute\n    # Unix time does not take the time zone into account -> add offset based on longitude\n    longitude_time_offset = (df_['fact_longitude'] * 12 / 180)\n    df_['floor_hour'] = (hour + minute/60 + longitude_time_offset)%24\n    df_['floor_hour'] = df_['floor_hour'].apply(np.floor).astype('int8')\n    return df_\n\ndef add_month(df_):\n    date_time = pd.to_datetime(df_['fact_time'], unit='s')\n    month = date_time.dt.month\n    # Unix time does not take the time zone into account -> add offset based on longitude\n    df_['month'] = month\n    return df_\n\nfrom sklearn.preprocessing import OneHotEncoder\n\ndef OHE_df(df_, feature):\n    encoder = OneHotEncoder(handle_unknown='ignore')\n    encoder_df = pd.DataFrame(encoder.fit_transform(df_[[feature]]).toarray())\n    encoder_df.columns = encoder_df.columns.astype('str')\n    encoder_df.columns = [feature +' : ' + col_name for col_name in encoder_df.columns]\n    final_df = pd.concat([df_, encoder_df], axis=1, join='inner')\n    final_df = final_df.drop([feature], axis = 1)\n    return final_df","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:22.946526Z","iopub.execute_input":"2022-05-11T15:56:22.948023Z","iopub.status.idle":"2022-05-11T15:56:22.967707Z","shell.execute_reply.started":"2022-05-11T15:56:22.947952Z","shell.execute_reply":"2022-05-11T15:56:22.966704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Preprocessing","metadata":{}},{"cell_type":"code","source":"from sklearn.decomposition import PCA\nfrom sklearn.preprocessing import RobustScaler\nfrom sklearn.impute import SimpleImputer\n\nclass DataProcessor:\n    def __init__(self, config_:dict):\n        self.c = config_\n        self.mean_imputer = None\n        self.pca = None\n        self.norm_mean = None\n        self.norm_std = None\n        self.min_max_scaler = None\n        self.robust_scaler = None\n        self.features = self.c['white_features']\n\n    def fit_transform(self, df_):\n        \"\"\"For training data\"\"\"\n        # Feature filter\n        if self.c['feature_white_list']:\n            df_ = df_[self.features]\n        \n        df_ = df_[df_['fact_temperature'] < 42.0]\n        df_ = df_[df_['fact_temperature'] > -30.0]\n        #df_ = df_[df_['fact_latitude'] < 40]\n        # df_ = df_[df_['fact_latitude'] > -40.0]\n        \n\n        # Drop NaN Samples\n        if self.c['drop_nan_samples']:\n            len_before = len(df_)\n            df_ = df_.dropna()\n            print(f\"Dropped {len_before - len(df_)} NaN samples\")\n\n        \n        # Fill NaN\n        if self.c['fill_nan']:\n            df_.fillna(self.c['fill_nan_value'], inplace=True)\n\n        for feat in self.c['power2_features']:\n            df_[f'{feat}_p2'] = np.power((df_[feat]),2)\n\n        X_, y_ = df_.drop(['fact_temperature'], axis=1), df_['fact_temperature']\n\n        \n        # Impute with mean\n        if self.c['impute_with_mean']:\n            cols = X_.columns\n            indices = X_.index\n            self.mean_imputer = SimpleImputer(missing_values=np.nan, strategy='mean')\n            self.mean_imputer.fit(X_)\n            X_ = pd.DataFrame(self.mean_imputer.transform(X_), columns=cols, index=indices)\n        \n        # Normalize\n        if self.c['normalize']:\n            X_ = self.normalize(X_, fit=True)\n        \n        if self.c['c_feat_hour']:\n            X_ = add_hour(X_)\n        \n        if self.c['one_hot_encode']:\n            X_ = add_floor_hour(X_)\n            X_ = OHE_df(X_, 'floor_hour')\n            \n        # PCA\n        if self.c['pca']:\n            self.pca = PCA(n_components = self.c['pca_components'])\n            principal_components = self.pca.fit_transform(X_)\n            principal_df = pd.DataFrame(data=principal_components)\n            X_ = X_.join(principal_df)\n            Z_ = X_.join(y_)\n            len_before = len(Z_)\n            Z_.dropna(inplace=True)\n            print(f\"Dropped {len_before - len(Z_)} training NaN examples after PCA\")\n            \n            X_, y_ = Z_.iloc[:, :-1], Z_.iloc[:, -1]\n        feature_names = X_.columns\n        X_, y_ = X_.to_numpy(), y_.to_numpy()\n\n        if self.c['min_max_scaler']:\n            self.min_max_scaler = MinMaxScaler()\n            X_ = self.min_max_scaler.fit_transform(X_)\n        if self.c['standard_scaler']:\n            self.standard = StandardScaler()\n            X_ = self.standard.fit_transform(X_)\n        if self.c['robust_scaler']:\n            self.robust = RobustScaler()\n            X_ = self.robust.fit_transform(X_)\n        \n        return X_, y_, feature_names\n\n    def transform(self, df_, target=True):\n        \"\"\"For evaluation data\"\"\"\n        # Feature Filter: Keep only specified features. Ignore target column if specified.\n        if self.c['feature_white_list']:\n            df_ = df_[self.features] if target else df_[[x for x in self.features if x != \"fact_temperature\"]]\n        \n        # drop_high_nan(df_, threshold=self.c['drop_high_nan_cols_thres'])\n        \n            \n        # Drop NaN Samples\n        if self.c['drop_nan_samples']:\n            df_.dropna(inplace=True)\n\n        # Fill NaN\n        if self.c['fill_nan']:\n            df_.fillna(self.c['fill_nan_value'], inplace=True)\n        \n        for feat in self.c['power2_features']:\n            df_[f'{feat}_p2'] = np.power((df_[feat]),2)\n\n        if target:\n            X_, y_ = df_.drop(['fact_temperature'], axis=1), df_['fact_temperature']\n        else:\n            X_, y_ = df_, None\n\n        # Impute with mean\n        if self.c['impute_with_mean']:\n            cols = X_.columns\n            indices = X_.index\n            X_ = pd.DataFrame(self.mean_imputer.transform(X_), columns=cols, index=indices)\n\n        # Normalize\n        if self.c['normalize']:\n            X_ = self.normalize(X_, fit=False)\n            \n        if self.c['c_feat_hour']:\n            X_ = add_hour(X_)\n            \n        if self.c['one_hot_encode']:\n            X_ = add_floor_hour(X_)\n            X_ = OHE_df(X_, 'floor_hour')\n            \n        # PCA\n        if self.c['pca']:\n            principal_components = self.pca.transform(X_)\n            principal_df = pd.DataFrame(data=principal_components)\n            X_ = X_.join(principal_df)\n            # X_ = principal_df # ATTENTION\n            if target:\n                Z_ = X_.join(y_)\n                len_before = len(Z_)\n                Z_.dropna(inplace=True)\n                \n                print(f\"Transform: Dropped {len_before - len(Z_)} training NaN examples after PCA\")\n            \n                X_, y_ = Z_.drop(['fact_temperature'], axis=1), Z_['fact_temperature']\n            else:\n                print(f\"Test set NaN: {np.isnan(X_).sum().sum()}\")\n        X_ = X_.to_numpy()\n        y_ = y_.to_numpy() if target else y_\n\n        if self.c['min_max_scaler']:\n            X_ = self.min_max_scaler.transform(X_)\n        \n        if self.c['standard_scaler']:\n            X_ = self.standard.transform(X_)\n            \n        if self.c['robust_scaler']:\n            X_ = self.robust.transform(X_)\n\n        return X_, y_\n    \n        \n    \n    def normalize(self, df_, fit=True):\n        if fit:\n            self.norm_mean = df_.mean()\n            self.norm_std = df_.std()\n        return (df_-self.norm_mean)/self.norm_std","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:22.970208Z","iopub.execute_input":"2022-05-11T15:56:22.970461Z","iopub.status.idle":"2022-05-11T15:56:23.03443Z","shell.execute_reply.started":"2022-05-11T15:56:22.970409Z","shell.execute_reply":"2022-05-11T15:56:23.033398Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def identify_covariate_shift(df_, df_test_, threshold=0.7):\n    # importing libraries\n    import numpy as np\n    import pandas as pd\n    from pandas import Series, DataFrame\n    import os\n    import matplotlib.pyplot as plt\n    from sklearn.ensemble import RandomForestRegressor\n    from sklearn.ensemble import RandomForestClassifier\n    from sklearn.model_selection import cross_val_score\n    from sklearn.preprocessing import LabelEncoder\n    from sklearn.impute import SimpleImputer\n    from sklearn.model_selection import train_test_split\n\n    # df_train, df_val = train_test_split(df, test_size=0.1, random_state=21)\n    X_train, y_train = df.drop(['fact_temperature'], axis=1), df_['fact_temperature']\n    X_test = df_test_\n\n    X_train['origin'] = 0\n    X_test['origin'] = 1\n    train_cols = X_train.columns\n    test_cols = X_test.columns\n    imp_mean = SimpleImputer(missing_values=np.nan, strategy='mean')\n    imp_mean.fit(X_train)\n    X_train = imp_mean.transform(X_train)\n    X_test = imp_mean.transform(X_test)\n    X_train = pd.DataFrame(X_train, columns=train_cols)\n    X_test = pd.DataFrame(X_test, columns=test_cols)\n\n    training = X_train.sample(10000, random_state=12)\n    testing = X_test.sample(10000, random_state=11)\n\n    ## combining random samples\n    combi = training.append(testing)\n    y = combi['origin']\n    combi.drop('origin',axis=1,inplace=True)\n\n    ## modelling\n    model = RandomForestClassifier(n_estimators = 50,n_jobs=-1, max_depth = 5,min_samples_leaf = 5, verbose=False)\n    cov_shift_dict = {}\n    for i in combi.columns:\n        score = cross_val_score(model,pd.DataFrame(combi[i]),y,cv=2,scoring='roc_auc')\n        cov_shift_dict[i] = np.mean(score)\n        # print(i, \": \",np.mean(score))\n        \n    keep_features_list = []\n    #for name, value in cov_shift_dict.items():\n    #    if value<threshold and name or name in best_features[:10]:\n    #        keep_features_list.append(name)\n    #    else:\n    #        print(f'Dropped: {name}')\n    df_test_ = df_test.drop(['origin'], axis=1)\n    df_covariance = pd.Series(cov_shift_dict).drop('index')\n    ax = df_covariance.plot.hist(bins=18)\n    ax.set_xlabel('Estimated Covariant Shift')\n    ax.set_ylabel('Number of Features')\n    fig = ax.get_figure()\n    fig.savefig('cov_shift.png')\n    return keep_features_list, cov_shift_dict, df_covariance\n#keep_features_list, _, df_covariance = identify_covariate_shift(df, df_test)\n# df_test = df_test.drop(['origin'], axis=1)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:23.039682Z","iopub.execute_input":"2022-05-11T15:56:23.040021Z","iopub.status.idle":"2022-05-11T15:56:23.06069Z","shell.execute_reply.started":"2022-05-11T15:56:23.039968Z","shell.execute_reply":"2022-05-11T15:56:23.059851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Model Selection and Submission\n\n","metadata":{"execution":{"iopub.status.busy":"2022-04-29T13:02:25.21223Z","iopub.execute_input":"2022-04-29T13:02:25.214379Z","iopub.status.idle":"2022-04-29T13:02:25.239339Z","shell.execute_reply.started":"2022-04-29T13:02:25.214317Z","shell.execute_reply":"2022-04-29T13:02:25.238252Z"}}},{"cell_type":"markdown","source":"### Configuration","metadata":{}},{"cell_type":"code","source":"config = {\n    # Data processing\n    'feature_white_list': False,\n    'validation_ratio': 0.05,\n    'power2_features': ['climate_temperature', 'gfs_soil_temperature', 'wrf_t2_interpolated', 'wrf_t2_next'],\n    #'white_features': ['fact_temperature', 'gfs_soil_temperature', 'wrf_t2_interpolated','wrf_t2_next', 'fact_time', 'fact_latitude', 'gfs_temperature_97500', 'gfs_temperature_95000', 'climate_temperature','gfs_temperature_92500', 'gfs_temperature_90000','gfs_temperature_85000', 'gfs_temperature_80000', 'cmc_0_0_6_2','gfs_temperature_75000', 'gfs_temperature_70000','gfs_temperature_65000', 'cmc_0_1_0_0', 'gfs_temperature_60000','gfs_temperature_55000', 'gfs_temperature_50000','gfs_temperature_45000', 'gfs_temperature_40000','gfs_2m_dewpoint_next', 'gfs_2m_dewpoint', 'gfs_temperature_35000','gfs_temperature_30000', 'cmc_0_3_5_500', 'gfs_precipitable_water','gfs_temperature_25000', 'cmc_0_3_5_700', 'cmc_0_0_7_2','sun_elevation', 'cmc_0_3_5_850', 'gfs_temperature_20000','cmc_0_0_7_1000', 'fact_longitude', 'climate_pressure', 'cmc_0_3_0_0','gfs_pressure', 'cmc_0_3_0_0_next', 'cmc_0_0_7_500', 'cmc_0_2_3_500'],\n    #'white_features': ['fact_time', 'topography_bathymetry', 'sun_elevation', 'climate_pressure', 'cmc_0_0_7_1000', 'cmc_0_0_7_2', 'cmc_0_0_7_500', 'cmc_0_0_7_700', 'cmc_0_0_7_850', 'cmc_0_0_7_925', 'cmc_0_1_11_0', 'cmc_0_1_65_0', 'cmc_0_1_66_0', 'cmc_0_1_67_0', 'cmc_0_1_68_0', 'cmc_0_1_7_0', 'cmc_0_2_2_10', 'cmc_0_2_2_1000', 'cmc_0_2_2_500', 'cmc_0_2_2_700', 'cmc_0_2_2_850', 'cmc_0_2_2_925', 'cmc_0_2_3_10', 'cmc_0_2_3_1000', 'cmc_0_2_3_500', 'cmc_0_2_3_700', 'cmc_0_2_3_850', 'cmc_0_2_3_925', 'cmc_0_3_0_0', 'cmc_0_3_0_0_next', 'cmc_0_3_1_0', 'cmc_0_3_5_1000', 'cmc_0_3_5_700', 'cmc_0_3_5_850', 'cmc_0_3_5_925', 'cmc_0_6_1_0', 'cmc_available', 'cmc_horizon_h', 'cmc_precipitations', 'cmc_timedelta_s', 'gfs_a_vorticity', 'gfs_available', 'gfs_cloudness', 'gfs_clouds_sea', 'gfs_horizon_h', 'gfs_humidity', 'gfs_precipitations', 'gfs_pressure', 'gfs_r_velocity', 'gfs_soil_temperature', 'gfs_soil_temperature_available', 'gfs_temperature_5000', 'gfs_temperature_20000', 'gfs_timedelta_s', 'gfs_total_clouds_cover_high', 'gfs_total_clouds_cover_low', 'gfs_total_clouds_cover_middle', 'gfs_u_wind', 'gfs_v_wind', 'gfs_wind_speed', 'wrf_available', 'wrf_psfc', 'wrf_rh2', 'wrf_wind_u', 'wrf_wind_v', 'wrf_rain', 'wrf_snow', 'wrf_graupel', 'wrf_t2_grad', 'cmc_0_1_65_0_grad', 'cmc_0_1_65_0_next', 'cmc_0_1_66_0_grad', 'cmc_0_1_66_0_next', 'cmc_0_1_67_0_grad', 'cmc_0_1_67_0_next', 'cmc_0_1_68_0_grad', 'cmc_0_1_68_0_next', 'gfs_2m_dewpoint_grad', 'gfs_total_clouds_cover_low_grad', 'gfs_total_clouds_cover_low_next'],\n    # All features except those with high covariate shift\n    'white_features': ['fact_temperature', 'fact_longitude', 'fact_time', 'topography_bathymetry', 'sun_elevation', 'climate_temperature', 'climate_pressure', 'cmc_0_0_7_1000', 'cmc_0_0_7_2', 'cmc_0_0_7_500', 'cmc_0_0_7_700', 'cmc_0_0_7_850', 'cmc_0_0_7_925', 'cmc_0_1_11_0', 'cmc_0_1_65_0', 'cmc_0_1_66_0', 'cmc_0_1_67_0', 'cmc_0_1_68_0', 'cmc_0_1_7_0', 'cmc_0_2_2_10', 'cmc_0_2_2_1000', 'cmc_0_2_2_500', 'cmc_0_2_2_700', 'cmc_0_2_2_850', 'cmc_0_2_2_925', 'cmc_0_2_3_10', 'cmc_0_2_3_1000', 'cmc_0_2_3_500', 'cmc_0_2_3_700', 'cmc_0_2_3_850', 'cmc_0_2_3_925', 'cmc_0_3_0_0', 'cmc_0_3_0_0_next', 'cmc_0_3_1_0', 'cmc_0_3_5_1000', 'cmc_0_3_5_700', 'cmc_0_3_5_850', 'cmc_0_3_5_925', 'cmc_0_6_1_0', 'cmc_available', 'cmc_horizon_h', 'cmc_precipitations', 'cmc_timedelta_s', 'gfs_a_vorticity', 'gfs_available', 'gfs_cloudness', 'gfs_clouds_sea', 'gfs_horizon_h', 'gfs_humidity', 'gfs_precipitations', 'gfs_pressure', 'gfs_r_velocity', 'gfs_soil_temperature', 'gfs_soil_temperature_available', 'gfs_temperature_5000', 'gfs_temperature_20000', 'gfs_temperature_80000', 'gfs_temperature_85000', 'gfs_temperature_90000', 'gfs_temperature_92500', 'gfs_temperature_95000', 'gfs_temperature_97500', 'gfs_timedelta_s', 'gfs_total_clouds_cover_high', 'gfs_total_clouds_cover_low', 'gfs_total_clouds_cover_middle', 'gfs_u_wind', 'gfs_v_wind', 'gfs_wind_speed', 'wrf_available', 'wrf_t2_next', 'wrf_psfc', 'wrf_rh2', 'wrf_wind_u', 'wrf_wind_v', 'wrf_rain', 'wrf_snow', 'wrf_graupel', 'wrf_t2_interpolated', 'wrf_t2_grad', 'cmc_0_1_65_0_grad', 'cmc_0_1_65_0_next', 'cmc_0_1_66_0_grad', 'cmc_0_1_66_0_next', 'cmc_0_1_67_0_grad', 'cmc_0_1_67_0_next', 'cmc_0_1_68_0_grad', 'cmc_0_1_68_0_next', 'gfs_2m_dewpoint_grad', 'gfs_total_clouds_cover_low_grad', 'gfs_total_clouds_cover_low_next', 'elevation'],    \n    'drop_nan_samples': False,\n    'fill_nan': False,\n    'fill_nan_value': -1.0,\n    'impute_with_mean': True,\n    'normalize': False,\n    'c_feat_hour': True,\n    'one_hot_encode': False,\n    'pca': True,\n    'pca_components': 10,\n    'min_max_scaler': False,\n    'standard_scaler': False,\n    'robust_scaler': True,\n    \n    \n    # NN\n    'nn_structure': [128, 128, 128, 128, 64, 1],\n    'nn_activation': ['selu', 'selu', 'selu', 'selu', 'selu', 'linear'],\n    'nn_dropout': [0.0, 0.0, 0.0, 0.0, 0.0, 0.0],\n    'batch_norm': False,\n    'batch_size': 128,\n    'epoch_count': 50,\n    'lr': 1e-3,\n    # XGB\n    \n}","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:23.062908Z","iopub.execute_input":"2022-05-11T15:56:23.063538Z","iopub.status.idle":"2022-05-11T15:56:23.090271Z","shell.execute_reply.started":"2022-05-11T15:56:23.063496Z","shell.execute_reply":"2022-05-11T15:56:23.089304Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Data","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nimport seaborn as sns\nfrom tqdm import tqdm\nimport matplotlib.pyplot as plt\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\nfrom torch.utils.data import Dataset, DataLoader\nfrom sklearn.preprocessing import MinMaxScaler\nfrom sklearn.preprocessing import StandardScaler\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor\nfrom sklearn.linear_model import Ridge, RidgeCV\nfrom sklearn.svm import SVR\nfrom sklearn.ensemble import StackingRegressor\nfrom xgboost import XGBRegressor\nimport pickle\nimport time\n\ndef rmse(y, y_pred):\n    return np.sqrt(mean_squared_error(y, y_pred))\n\ndef fit_predict_eval(method, name:str, verbose=True, model_file=None):\n    start_time = time.time()\n    if verbose:\n        print(name)\n    if model_file:\n        model = pickle.load(open(model_file, 'rb'))\n    else:\n        model = method.fit(X_train, y_train)\n        pickle.dump(model, open(f\"/kaggle/working/{name}_model.sav\", 'wb'))\n    y_pred_train = model.predict(X_train)\n    if verbose:\n        print(f'{name} RSME training: {rmse(y_train, y_pred_train)}')\n    y_pred_val = model.predict(X_val)\n    val_score = rmse(y_val, y_pred_val)\n    if verbose:\n        print(f'{name} RSME validation: {val_score}')\n    y_pred_test = model.predict(X_test)\n    submission_df_ = pd.DataFrame(data={'index': df_test['index'].values, 'fact_temperature': y_pred_test.squeeze()})\n    if verbose:\n        print(f'Finished {name} in {time.time()-start_time} s')\n    return submission_df_, model, val_score\n    \nnn_test_pred = np.ndarray([0])\n\ndevice = torch.device(\"cuda:0\" if torch.cuda.is_available() else \"cpu\")\n\ndf_train, df_val = train_test_split(df, test_size=config['validation_ratio'], random_state=21)\n\nprocessor = DataProcessor(config)\nX_train, y_train, feature_names = processor.fit_transform(df_train)\nX_val, y_val = processor.transform(df_val)\nX_test, _ = processor.transform(df_test, target=False)\n\n\nX_train, y_train = np.array(X_train), np.array(y_train)\nX_val, y_val = np.array(X_val), np.array(y_val)\n\nprint(f\"Training samples: {len(X_train)} with {X_train.shape[1]} features\")","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:56:23.096568Z","iopub.execute_input":"2022-05-11T15:56:23.097502Z","iopub.status.idle":"2022-05-11T15:57:31.018901Z","shell.execute_reply.started":"2022-05-11T15:56:23.097453Z","shell.execute_reply":"2022-05-11T15:57:31.01717Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Models","metadata":{}},{"cell_type":"markdown","source":"Neural Network","metadata":{}},{"cell_type":"code","source":"class MultipleRegression(nn.Module):\n    def __init__(self, num_features, config_):\n        super(MultipleRegression, self).__init__()\n        self.config = config_\n        self.layers = []\n        self.structure = self.config['nn_structure']\n        for idx, neuron_count in enumerate(self.structure):\n            prev_neuron_count = num_features if idx==0 else self.structure[idx-1]\n            self.layers.append(nn.Linear(prev_neuron_count, neuron_count))\n            self.layers.append(nn.Dropout(p=self.config['nn_dropout'][idx]))\n            if self.config['batch_norm']:\n                self.layers.append(nn.BatchNorm1d(neuron_count))\n            if self.config['nn_activation'] == 'relu':\n                self.layers.append(nn.ReLU())\n            elif self.config['nn_activation'] == 'selu':\n                self.layers.append(nn.SELU())\n        self.classifier = nn.Sequential(*self.layers)\n\n    def forward(self, inputs):\n        return self.classifier(inputs)\n\n\ndef RMSELoss(output, target):\n    mse_loss = nn.MSELoss()\n    return torch.sqrt(mse_loss(output, target))\n\ndef run_nn():\n\n    EPOCHS = config['epoch_count']\n    BATCH_SIZE = config['batch_size']\n    LEARNING_RATE = config['lr']\n    NUM_FEATURES = X_train.shape[1]    \n\n    class RegressionDataset(Dataset):\n\n        def __init__(self, X_data, y_data):\n            self.X_data = X_data\n            self.y_data = y_data\n\n        def __getitem__(self, index):\n            return self.X_data[index], self.y_data[index]\n\n        def __len__ (self):\n            return len(self.X_data)\n    train_dataset = RegressionDataset(torch.from_numpy(X_train).float(), torch.from_numpy(y_train).float())\n    val_dataset = RegressionDataset(torch.from_numpy(X_val).float(), torch.from_numpy(y_val).float())\n    test_dataset = RegressionDataset(torch.from_numpy(X_test).float(), torch.zeros(len(X_test)))\n\n    train_loader = DataLoader(dataset=train_dataset, batch_size=BATCH_SIZE, shuffle=True)\n    val_loader = DataLoader(dataset=val_dataset, batch_size=1024)\n    test_loader = DataLoader(dataset=test_dataset, batch_size=1024)\n\n    model = MultipleRegression(NUM_FEATURES, config)\n    model.to(device)\n    print(model)\n    criterion = RMSELoss\n    optimizer = optim.Adam(model.parameters(), lr=LEARNING_RATE)    \n    y_test_pred_all = np.array([0])\n\n\n    loss_stats = {\n        'train': [],\n        \"val\": []\n    }   \n\n    print(\"Begin training.\")\n    for e in tqdm(range(1, EPOCHS+1)):\n\n        # TRAINING\n        train_epoch_loss = 0\n        model.train()\n        for X_train_batch, y_train_batch in (train_loader):\n            X_train_batch, y_train_batch = X_train_batch.to(device), y_train_batch.to(device)\n            optimizer.zero_grad()\n            y_train_pred = model(X_train_batch)\n            train_loss = criterion(y_train_pred, y_train_batch.unsqueeze(1))\n\n            train_loss.backward()\n            optimizer.step()\n\n            train_epoch_loss += train_loss.item()\n\n\n        #VALIDATION + TESTING\n        with torch.no_grad():\n\n            val_epoch_loss = 0\n\n            model.eval()\n            for X_val_batch, y_val_batch in val_loader:\n                X_val_batch, y_val_batch = X_val_batch.to(device), y_val_batch.to(device)\n\n                y_val_pred = model(X_val_batch)\n\n                val_loss = criterion(y_val_pred, y_val_batch.unsqueeze(1))\n\n                val_epoch_loss += val_loss.item()\n\n\n        loss_stats['train'].append(train_epoch_loss/len(train_loader))\n        loss_stats['val'].append(val_epoch_loss/len(val_loader))                              \n\n        print(f'Epoch {e+0:03}: | Train Loss: {train_epoch_loss/len(train_loader):.5f} | Val Loss: {val_epoch_loss/len(val_loader):.5f}')\n     #%%\n    with torch.no_grad():\n        model.eval()\n        for X_test_batch, _ in test_loader:\n            X_test_batch = X_test_batch.to(device)\n            y_test_pred = model(X_test_batch).cpu().numpy().squeeze()\n            nn_test_pred = np.concatenate((nn_test_pred, y_test_pred))\n    return pd.DataFrame(data={'index': df_test['index'].values.astype('int32'), 'fact_temperature': nn_test_pred})\n\n# nn_submission_df = run_nn()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:57:31.020779Z","iopub.execute_input":"2022-05-11T15:57:31.02109Z","iopub.status.idle":"2022-05-11T15:57:31.051605Z","shell.execute_reply.started":"2022-05-11T15:57:31.02105Z","shell.execute_reply":"2022-05-11T15:57:31.05078Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"CatBoost","metadata":{}},{"cell_type":"code","source":"from catboost import CatBoostRegressor\n\n \ncat = CatBoostRegressor(\n    loss_function ='RMSE',\n    iterations=50000,\n    learning_rate=0.12,\n    task_type=\"GPU\",\n    depth=8\n)\ncat_submission_df, cat_model, cat_val_score = fit_predict_eval(cat, name='cat')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:03:47.963433Z","iopub.execute_input":"2022-05-11T16:03:47.963723Z","iopub.status.idle":"2022-05-11T16:06:27.776207Z","shell.execute_reply.started":"2022-05-11T16:03:47.963692Z","shell.execute_reply":"2022-05-11T16:06:27.775316Z"},"_kg_hide-input":false,"collapsed":true,"jupyter":{"outputs_hidden":true},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"submission_df = cat_submission_df\nsubmission_df['index'] = submission_df['index'].astype('int64')\nsubmission_df.to_csv(\"/kaggle/working/submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:59:34.060939Z","iopub.execute_input":"2022-05-11T15:59:34.061399Z","iopub.status.idle":"2022-05-11T15:59:36.037759Z","shell.execute_reply.started":"2022-05-11T15:59:34.061359Z","shell.execute_reply":"2022-05-11T15:59:36.03693Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_imp = cat.get_feature_importance()","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:07:58.166176Z","iopub.execute_input":"2022-05-11T16:07:58.166467Z","iopub.status.idle":"2022-05-11T16:07:58.259907Z","shell.execute_reply.started":"2022-05-11T16:07:58.166436Z","shell.execute_reply":"2022-05-11T16:07:58.259158Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize=(5,5))\nplt.xlabel(\"CatBoost Feature Importance (top 15)\")\nax=(pd.Series(feature_imp, index=feature_names)/100).sort_values(ascending=True)[-15:].plot.barh()\nfig = ax.get_figure()\nfig.tight_layout()\nfig.savefig('cat_feat_imp.png')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:33:10.904971Z","iopub.execute_input":"2022-05-11T16:33:10.905581Z","iopub.status.idle":"2022-05-11T16:33:11.321905Z","shell.execute_reply.started":"2022-05-11T16:33:10.905541Z","shell.execute_reply":"2022-05-11T16:33:11.321034Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:17:07.414379Z","iopub.execute_input":"2022-05-11T16:17:07.414771Z","iopub.status.idle":"2022-05-11T16:17:07.428439Z","shell.execute_reply.started":"2022-05-11T16:17:07.414735Z","shell.execute_reply":"2022-05-11T16:17:07.427059Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"GB","metadata":{}},{"cell_type":"code","source":"gbr = GradientBoostingRegressor(\n    n_estimators=500, # 500\n    max_features='sqrt',\n    learning_rate=0.02,\n    max_depth=5,\n    loss='huber',\n    min_samples_split=8,\n    min_samples_leaf=13,\n    random_state=76\n)\n#gbr_submission_df, gbr_model, gbr_val_score = fit_predict_eval(gbr, name='gbr')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:59:36.043734Z","iopub.execute_input":"2022-05-11T15:59:36.044421Z","iopub.status.idle":"2022-05-11T15:59:36.053709Z","shell.execute_reply.started":"2022-05-11T15:59:36.044372Z","shell.execute_reply":"2022-05-11T15:59:36.052851Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"RF","metadata":{}},{"cell_type":"code","source":"rf = RandomForestRegressor(\n    n_estimators=180,\n    max_depth=28, # 20\n    min_samples_split=3, #3\n    min_samples_leaf=2, #2\n    oob_score=True,\n    n_jobs=-1,\n    # max_features=None,\n    random_state=86,\n    verbose=1\n)\n\n#rf_submission_df, rf_model, rf_val_score = fit_predict_eval(rf, name='rf')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:59:36.057736Z","iopub.execute_input":"2022-05-11T15:59:36.058082Z","iopub.status.idle":"2022-05-11T15:59:36.071569Z","shell.execute_reply.started":"2022-05-11T15:59:36.05804Z","shell.execute_reply":"2022-05-11T15:59:36.070435Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"XGB","metadata":{}},{"cell_type":"code","source":"xgb = XGBRegressor(\n    objective ='reg:squarederror',\n    eval_metric='mae',\n    learning_rate = 0.03,\n    max_depth = 4, \n    alpha = 6,\n    n_estimators = 500, # 600\n    nthread=-1,\n    random_state=35,\n    seed=64\n)\n# xgb_submission_df, xgb_model, xgb_val_score = fit_predict_eval(xgb, name='xgboost')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:59:36.073201Z","iopub.execute_input":"2022-05-11T15:59:36.07401Z","iopub.status.idle":"2022-05-11T15:59:36.082927Z","shell.execute_reply.started":"2022-05-11T15:59:36.073935Z","shell.execute_reply":"2022-05-11T15:59:36.081617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ridge Regressor","metadata":{}},{"cell_type":"code","source":"r_best_alpha = 0\nr_best_val_score = 0\nr_alphas = [1e-15, 1e-10, 1e-5, 1e-3, 0.1, 0.4, 1, 5, 10, 15, 18, 20, 30, 50, 75, 100]\nfor r_alpha in r_alphas:\n    ridge = Ridge(alpha=r_alpha)\n    ridge_submission_df, ridge_model, ridge_val_score = fit_predict_eval(ridge, name='ridge', verbose=False)\n    if ridge_val_score>r_best_val_score:\n        r_best_alpha = r_alpha\n        r_best_val_score = ridge_val_score\nridge = Ridge(alpha=r_best_alpha)\nridge_submission_df, ridge_model, ridge_val_score = fit_predict_eval(ridge, name='ridge')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T15:59:36.08522Z","iopub.execute_input":"2022-05-11T15:59:36.085573Z","iopub.status.idle":"2022-05-11T16:02:39.989376Z","shell.execute_reply.started":"2022-05-11T15:59:36.085528Z","shell.execute_reply":"2022-05-11T16:02:39.988375Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Ensemble","metadata":{}},{"cell_type":"code","source":"if False:\n    submission_df = gbr_submission_df.copy(deep=True)\n    submission_df['fact_temperature'] = 0.0\n    cat0 = []\n    cat1 = []\n    cat2 = []\n    cat3 = []\n    tuple_list = [\n        (gbr_val_score, gbr_submission_df),\n        (rf_val_score, rf_submission_df),\n        (xgb_val_score, xgb_submission_df),\n        (ridge_val_score, ridge_submission_df),\n        # (svr_val_score, svr_submission_df),\n        #(stack_val_score, stack_submission_df),\n    ]\n\n    for score, sub_df in tuple_list:\n        if score < 1.95:\n            print('cat0')\n            cat0.append(sub_df)\n        if score < 2.0:\n            print('cat1')\n            cat1.append(sub_df)\n        elif score < 2.15:\n            print('cat2')\n            cat2.append(sub_df)\n        elif score < 2.5:\n            print('cat3')\n            cat3.append(sub_df)\n        else:\n            print('no cat')\n    for sub_df in cat0:\n        submission_df['fact_temperature'] = submission_df['fact_temperature'] + sub_df['fact_temperature']*6\n    for sub_df in cat1:\n        submission_df['fact_temperature'] = submission_df['fact_temperature'] + sub_df['fact_temperature']*4\n    for sub_df in cat2:\n        submission_df['fact_temperature'] = submission_df['fact_temperature'] + sub_df['fact_temperature']*2\n    for sub_df in cat3:\n        submission_df['fact_temperature'] = submission_df['fact_temperature'] + sub_df['fact_temperature']\n\n    submission_df['fact_temperature'] = submission_df['fact_temperature'] / (len(cat0)*4 + len(cat1)*4 + len(cat2)*2 + len(cat3))\n    submission_df['index'] = submission_df['index'].astype('int64')\n    submission_df.to_csv(\"/kaggle/working/submission.csv\", index=False)","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:02:39.991132Z","iopub.execute_input":"2022-05-11T16:02:39.99163Z","iopub.status.idle":"2022-05-11T16:02:40.011625Z","shell.execute_reply.started":"2022-05-11T16:02:39.991583Z","shell.execute_reply":"2022-05-11T16:02:40.010397Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Stack","metadata":{}},{"cell_type":"code","source":"stack = StackingRegressor(\n    estimators=[('GBR', gbr), ('RF', rf), ('XGB', xgb), ('Ridge', ridge)],\n    final_estimator=xgb,\n    n_jobs=-1,\n    verbose=1\n)\n# stack_submission_df, stack_model, stack_val_score = fit_predict_eval(stack, name='stack')","metadata":{"execution":{"iopub.status.busy":"2022-05-11T16:02:40.013545Z","iopub.execute_input":"2022-05-11T16:02:40.013895Z","iopub.status.idle":"2022-05-11T16:02:40.02754Z","shell.execute_reply.started":"2022-05-11T16:02:40.013848Z","shell.execute_reply":"2022-05-11T16:02:40.026367Z"},"trusted":true},"execution_count":null,"outputs":[]}]}