{"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":"# AML Flight delay challenge\n## Imports","metadata":{"id":"BdCEBnJzrWm3"}},{"cell_type":"code","source":"import sys\nimport time\nimport json\nimport pandas as pd\nimport numpy as np\nimport pickle\nfrom tqdm import tqdm\n\nfrom sklearn.metrics.pairwise import euclidean_distances\n\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.model_selection import train_test_split, GridSearchCV\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.linear_model import Ridge, BayesianRidge\nfrom sklearn.preprocessing import StandardScaler, PolynomialFeatures\nfrom sklearn.pipeline import make_pipeline\nfrom matplotlib import rc\nimport matplotlib.pyplot as plt\n\nimport torch\nimport torch.nn as nn\nimport torch.optim as optim\n\ntorch.set_default_tensor_type(torch.FloatTensor)\ndevice = \"cuda:0\" if torch.cuda.is_available() else \"cpu\"\nprint(device)\n\n!pip install gpflow\nimport gpflow\nfrom gpflow.utilities import print_summary\nfrom gpflow import set_trainable\nimport tensorflow as tf\n\nimport xgboost","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","id":"MmNSG9vvrWm6","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from distutils.dir_util import copy_tree\nfrom_dir = '../input/eurecom-aml-2021-challenge-1'\nto_dir = ''\ncopy_tree(from_dir, to_dir)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from_dir = '../input/weatherflights'\nto_dir = ''\ncopy_tree(from_dir, to_dir)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Preprocessing\nThis section reads all data and merges everthing in a single dataframe called `data`. That one must be used for inference and contains:\n- startting features\n- merged info about airports\n- merged info about weather","metadata":{"id":"4qEub2UQrWm7"}},{"cell_type":"code","source":"def funzione_del_tempo(weather_df, airports_df, train_df):\n    \"\"\"\n    A very english-friendly funzione_del_tempo that cleans the weather dataframe.\n    \n    Args:\n        The three dataframes, c'mon\n        \n    Returns:\n        A dataframe with the merged weather data for both origin and destination.\n        These columns will have suffixes '_origin' and '_destination' respectively.\n    \"\"\"\n    \n    # Expanding date format\n    date_cols = np.array(list(map(lambda date: [int(v) for v in date.split('-')], weather_df['DATE'].values)))\n    weather_df[[\"MONTH\", \"DAY\"]] = date_cols[:, 1:]\n    weather_df = weather_df.drop(['DATE'], axis=1)\n\n    # Evaluating the min distance between airports and weather station\n    distinct_latlon_flight_origin = train_df[['LATITUDE_ORIGIN', 'LONGITUDE_ORIGIN', 'ORIGIN_AIRPORT']].drop_duplicates()\n    distinct_latlon_flight_destination = train_df[['LATITUDE_DESTINATION', 'LONGITUDE_DESTINATION', 'DESTINATION_AIRPORT']].drop_duplicates()\n    rename_cols_origin = {'LATITUDE_ORIGIN': 'LATITUDE', 'LONGITUDE_ORIGIN': 'LONGITUDE', 'ORIGIN_AIRPORT': 'AIRPORT'}\n    rename_cols_destination = {'LATITUDE_DESTINATION': 'LATITUDE', 'LONGITUDE_DESTINATION': 'LONGITUDE', 'DESTINATION_AIRPORT': 'AIRPORT'}\n    distinct_latlon_flight = pd.concat([distinct_latlon_flight_origin.rename(rename_cols_origin, axis = 1), distinct_latlon_flight_destination.rename(rename_cols_destination, axis = 1)], axis = 0).drop_duplicates().values\n    \n    # Extract coords\n    iata_code = np.array(list(map(lambda v : v[2], distinct_latlon_flight)))\n    distinct_latlon_flight = np.array(list(map(lambda v : [v[0], v[1]], distinct_latlon_flight)))\n    \n    distinct_latlon_weather = weather_df[['LATITUDE', 'LONGITUDE', 'STATION']].value_counts().index.values\n    stations = np.array(list(map(lambda v : v[2], distinct_latlon_weather)))\n    distinct_latlon_weather = np.array(list(map(lambda v : [v[0], v[1]],  distinct_latlon_weather)))\n\n    distances = euclidean_distances(X = distinct_latlon_flight, Y = distinct_latlon_weather)\n    \n    # Creating a dict to map a wheather station to an airport\n    latlon_weather_airport_mapping = {}\n    for code, weather_latlon_id in zip(iata_code, distances.argmin(1)):\n        latlon_weather_airport_mapping[code] = stations[weather_latlon_id]\n\n    # attach the station to the airport\n    train_df['ORIGIN_STATION_WEATHER'] = list(map(lambda v: latlon_weather_airport_mapping[v] if v in latlon_weather_airport_mapping.keys() else 'NA' , train_df['ORIGIN_AIRPORT'].values))\n    train_df['DESTINATION_STATION_WEATHER'] = list(map(lambda v: latlon_weather_airport_mapping[v] if v in latlon_weather_airport_mapping.keys() else 'NA' , train_df['DESTINATION_AIRPORT'].values))\n    \n    # Finally merge the two things\n    origin_merge_df = pd.merge(\n      train_df, \n      weather_df,\n      how = 'left',\n      left_on = ['ORIGIN_STATION_WEATHER', 'MONTH', 'DAY'],\n      right_on = ['STATION', 'MONTH', 'DAY'],\n    )\n    origin_destination_merge_df = pd.merge(\n      origin_merge_df, \n      weather_df,\n      how = 'left',\n      suffixes=('_origin', '_destination'),\n      left_on = ['DESTINATION_STATION_WEATHER', 'MONTH', 'DAY'],\n      right_on = ['STATION', 'MONTH', 'DAY'],\n    )\n    #cols_to_drop = ['ORIGIN_STATION_WEATHER', 'DESTINATION_STATION_WEATHER', 'LATITUDE_origin', 'LONGITUDE_origin', 'LATITUDE_destination', 'LONGITUDE_destination'] \n    cols_to_drop = ['STATION_destination', 'LATITUDE_destination', 'LONGITUDE_destination', 'STATION_origin', 'LATITUDE_origin', 'LONGITUDE_origin', 'ORIGIN_STATION_WEATHER', 'DESTINATION_STATION_WEATHER']\n    \n    origin_destination_merge_df = origin_destination_merge_df.drop(cols_to_drop, axis = 1)\n\n    return origin_destination_merge_df","metadata":{"id":"oXNRsGjy5KUh","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def fill_nan_weather(weather):\n    \"\"\"\n    Fill NaNs in weather data. This data has a very peculiar format and NaNs take on different values according to the columns. So we need to clean them manually. We replace the NaNs with the mean value of each column.\n    \n    Returns:\n        the filled weather data\n    \"\"\"\n\n    nan_mapping = {\n        'TEMP': 9999.9,\n        'DEWP': 9999.9,\n        'SLP': 9999.9,\n        'STP': 999.9,\n        'VISIB': 999.9,\n        'WDSP': 999.9,\n        'MXSPD': 999.9,\n        'GUST': 999.9,\n        'MAX': 9999.9,\n        'MIN': 9999.9,\n        'PRCP': 99.99,\n        'SNDP': 999.9,\n    }\n\n    for ft in weather.columns:\n        if ft in nan_mapping.keys():\n            nan_value = nan_mapping[ft]\n            if ft == 'PRCP':\n                fill = 0\n            else:\n                fill = weather.loc[weather[ft] != nan_value][ft].mean()\n            print(f'Found {(weather[ft] == nan_value).sum()} missing values from {ft} replaced with {fill}')\n            weather[ft] = weather[ft].replace(nan_value, fill)\n    return weather","metadata":{"id":"eu3WPaVrZDX6","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def preprocess_data(features, airports, weather=None, targets=None, one_hot_time=False, realistic_fts=False):\n    \"\"\"\n    Preprocess data to a classifier-ready representation.\n    \n    Args:\n        features: the initial flights dataframe\n        airports: the airports dataframe\n        weather: dataframe with weather data\n        targets: the arrival delay column\n        one_hot_time: whether to encode time as one-hot columns. It's completely useless and maybe a bit stupid, so please don't do it\n        realistic_fts: if True, remove all features that can only be known after takeoff (e.g. departure delay). If False, all features are kept\n    \"\"\"\n    # 1) Keep only certain cols\n    known_in_advance_fts = ['ID','MONTH', 'DAY','DAY_OF_WEEK','AIRLINE',\n                            'ORIGIN_AIRPORT','DESTINATION_AIRPORT','SCHEDULED_DEPARTURE',\n                           'SCHEDULED_TIME','DISTANCE','SCHEDULED_ARRIVAL']\n    \n    known_on_the_spot_fts = ['DEPARTURE_TIME', 'DEPARTURE_DELAY', 'TAXI_OUT', 'AIR_TIME', 'WHEELS_OFF']\n    \n    if realistic_fts:\n        features_cols = known_in_advance_fts\n    else:\n        features_cols = known_in_advance_fts + known_on_the_spot_fts\n    \n    airports_cols = ['IATA_CODE','LATITUDE','LONGITUDE']\n    \n    features = features[features_cols]\n    airports = airports[airports_cols]\n    print(f'Initial shape after dropping cols: {features.shape}')\n    \n    # 2) Merge other information from the airport (mainly latitude/longitude)\n    # Origin lat/lon\n    features = pd.merge(features, airports.rename({'LATITUDE':'LATITUDE_ORIGIN', 'LONGITUDE':'LONGITUDE_ORIGIN'}, axis='columns', errors = 'raise'), how='left', left_on='ORIGIN_AIRPORT', right_on = 'IATA_CODE')\n    # Dest lat/lon\n    features = pd.merge(features, airports.rename({'LATITUDE':'LATITUDE_DESTINATION', 'LONGITUDE':'LONGITUDE_DESTINATION'}, axis='columns',errors = 'raise'), how='left', left_on='DESTINATION_AIRPORT', right_on = 'IATA_CODE')\n    # Merge the targets too, if we have them\n    if targets is not None:\n        features = pd.merge(features, targets, left_on='ID', right_on='ID')\n\n\n    # 2.something) Some values will have NaN lat/lon. We fill it with the mean lat and lon\n    ftstofill = ['LATITUDE_ORIGIN', 'LONGITUDE_ORIGIN', 'LATITUDE_DESTINATION', 'LONGITUDE_DESTINATION']\n        \n    for ft in ftstofill:\n        features[ft] = features[ft].fillna(value=features[ft].mean())\n\n\n    # 2.∞) facciamo cose magiche per aggiungere l'informazione del tempo\n    if weather is not None:\n      weather = weather.drop([\n          'Unnamed: 0', \n          'Unnamed: 0.1',\n          'NAME', \n          'TEMP_ATTRIBUTES', \n          'MAX_ATTRIBUTES', \n          'MIN_ATTRIBUTES',\n          'DEWP_ATTRIBUTES',\n          'STP_ATTRIBUTES',\n          'SLP_ATTRIBUTES',\n          'VISIB_ATTRIBUTES',\n          'WDSP_ATTRIBUTES',\n          'MAX_ATTRIBUTES',\n          'MIN_ATTRIBUTES',\n          'PRCP_ATTRIBUTES',\n          'ELEVATION',\n          'FRSHTT'], axis= 1)\n      weather = fill_nan_weather(weather)\n      features = funzione_del_tempo(weather, airports, features)\n    \n    # 3) We have created some rubbish features by merging things. We now remove them\n    # (actually we also remove some more features we initially chose to include since we change our mind every 6 seconds)\n    discarded_feat = ['ID','ORIGIN_AIRPORT', 'DESTINATION_AIRPORT', 'IATA_CODE_x', 'IATA_CODE_y']\n    if targets is None:\n        discarded_feat.remove('ID') # If we have no targets, we probably want to keep the IDS and discard them manually later on\n    \n    \n    categorical_feat = ['AIRLINE']\n    categorical_time = ['DAY','DAY_OF_WEEK','MONTH']\n    if one_hot_time:\n        categorical_feat += categorical_time\n    \n    # Discard features\n    features = features.drop(discarded_feat, axis=1)\n    \n    # Generation of categorical features\n    X_cat = features[categorical_feat].astype('category')\n    features = pd.concat([features, pd.get_dummies(X_cat)],axis=1)\n    features = features.drop(categorical_feat, axis=1)\n    print(f'Final shape: {features.shape}')\n    \n    return features","metadata":{"id":"Wqi5J3xKrWm7","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"folder = 'data'\ntrain_df = pd.read_csv(f'{folder}/train_features.csv', header=0)\ntargets_df = pd.read_csv(f'{folder}/train_targets.csv', header=0)\nairports_df = pd.read_csv(f'{folder}/airports.csv', header = 0)\ntest_df = pd.read_csv(f'{folder}/test_features.csv', header = 0)\nweather_df = pd.read_csv('weather365.csv', header = 0)","metadata":{"id":"-BuxOQTTrWm8","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Preprocess\ndata = preprocess_data(train_df, airports_df, weather=weather_df, targets=targets_df, realistic_fts=False)","metadata":{"id":"C2bxtQ6yrWm8","outputId":"150f4cab-9083-4a10-c7ad-f5f2948389e9","scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#add info on airport avg delays\ndelay_origin = data.groupby(['LATITUDE_ORIGIN','LONGITUDE_ORIGIN'])['DEPARTURE_DELAY'].aggregate(avg_dep_delay=np.mean)\ndelay_dest = data.groupby(['LATITUDE_DESTINATION','LONGITUDE_DESTINATION'])['ARRIVAL_DELAY'].aggregate(avg_arr_delay=np.mean)\n\ndata = pd.merge(data, delay_origin, how='left', left_on='LATITUDE_ORIGIN', right_on = 'LATITUDE_ORIGIN')\ndata = pd.merge(data, delay_dest, how='left', left_on='LATITUDE_DESTINATION', right_on = 'LATITUDE_DESTINATION')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Do this if your RAM explodes\ndel train_df","metadata":{"id":"9YyAzFGwLPa5","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data.columns","metadata":{"id":"_IoxkqjoLlAC","outputId":"684044a5-37a2-483e-88a0-7e25e7af3e87","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Final evaluation\nThe function used to produce the final submission.","metadata":{"id":"fYOLKSrIrWnH"}},{"cell_type":"code","source":"def produce_submission(ids, y_pred, filename):\n    out = pd.Series(data=y_pred, index=ids, name='ARRIVAL_DELAY')\n    out.to_csv(filename, header=True)","metadata":{"id":"JvQMIe4IrWnH","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Xtest = preprocess_data(test_df, airports_df, weather=weather_df, realistic_fts=True)\n\n#add info on airport avg delays\nXtest = pd.merge(Xtest, delay_origin, how='left', left_on='LATITUDE_ORIGIN', right_on = 'LATITUDE_ORIGIN')\nXtest = pd.merge(Xtest, delay_dest, how='left', left_on='LATITUDE_DESTINATION', right_on = 'LATITUDE_DESTINATION')","metadata":{"id":"bIhCIBEkrWnH","outputId":"80473ef0-a772-4e10-922d-ccc64f054047","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Models: Ridge, GP, GPFlow, XGBoost","metadata":{}},{"cell_type":"markdown","source":"## Ridge","metadata":{}},{"cell_type":"code","source":"val_split_size = 0.4\nestimator = Ridge(tol = 0.001)\nparam_grid = {'alpha': [0.001, 0.01, 0.1, 0., 1.0, 10.0]}","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X, y = data.drop('ARRIVAL_DELAY', axis=1), data['ARRIVAL_DELAY']\n\nXtrain, Xval, ytrain, yval = train_test_split(X, y, test_size = val_split_size, random_state = 42)\n\npipe = make_pipeline(StandardScaler(),\n                     SimpleImputer(),\n                     GridSearchCV(\n                        estimator = estimator,\n                        param_grid = param_grid,\n                        scoring = 'neg_root_mean_squared_error',\n                        verbose=4,\n                        n_jobs=1\n                        )\n                    )\n\npipe.fit(Xtrain,ytrain)\nprint(f'Best score in grid search: {round(pipe[\"gridsearchcv\"].best_score_, 3)}')\n\nval_score = mean_squared_error(yval, pipe.predict(Xval), squared=False)\nprint(f'Best score on validation set: {round(val_score, 3)}')\n\n# saving best classifier \nwith open(f'{int(time.time())}_model_ridge_mine.pickle', 'wb') as f:\n    pickle.dump(pipe, f)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Bayesian Ridge","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import cross_val_score\n\nestimator = BayesianRidge()\nX, y = data.drop('ARRIVAL_DELAY', axis=1), data['ARRIVAL_DELAY']\n\npipe = make_pipeline(StandardScaler(),\n                     SimpleImputer(),\n                     estimator\n                    )","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scores_cv = cross_val_score(pipe, X, y, verbose=10, scoring='neg_root_mean_squared_error')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'Averaged RMSE between folds: {np.array(scores_cv).mean()}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pipe.fit(X, y)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"y_pred = pipe.predict(Xtest.drop('ID', axis=1))\nids = Xtest['ID']\nproduce_submission(ids, y_pred, 'sub.csv')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Also producing results by trying to guess the overall interval of the delay. We count 1 if we guess, 0 otherwise","metadata":{}},{"cell_type":"code","source":"def conf_int_scorer(estimator, X, y, n_stds=1):\n    \"\"\"\n    Computes score by guessing the interval.\n    \"\"\"\n    means, stds = estimator.predict(X, return_std=True)\n    corrects = abs((n_stds * stds)) >= abs(means - y)\n    return corrects.mean()\n\ndef rich_conf_int_scorer(estimator, X, y, n_stds=1):\n    \"\"\"\n    Computes score by guessing the interval.\n    \"\"\"\n    means, stds = estimator.predict(X, return_std=True)\n    corrects = abs((n_stds * stds)) >= abs(means - y)\n    scores = {\n        'interval_accuracy': corrects.mean(),\n        'mean_interval_width': (stds * n_stds).mean()\n    }\n    return scores","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scores_cv_intervals = cross_val_score(pipe, X, y, verbose=10, scoring=conf_int_scorer)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'Confidence intervals accuracy: {np.array(scores_cv_intervals).mean()}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"A more refined version","metadata":{}},{"cell_type":"code","source":"from sklearn.model_selection import cross_validate\ncross_val_intervals = cross_validate(pipe, X, y, verbose=10, scoring=rich_conf_int_scorer)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'Mean accuracy: {cross_val_intervals[\"test_interval_accuracy\"].mean()}')\nprint(f'Mean interval width (std*2): {(cross_val_intervals[\"test_mean_interval_width\"]*2).mean()}')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## GPs","metadata":{}},{"cell_type":"code","source":"def args_as_tensors(*index):\n    def decorator(method):\n        def wrapper(*args, **kwargs):\n            converted_args = [torch.tensor(a).float() \n                              if i in index and type(a) is np.ndarray else a \n                              for i,a in enumerate(args)]\n            return method(*converted_args, **kwargs)\n        return wrapper  \n    return decorator\n\nclass RBF(nn.Module):\n    @args_as_tensors(0, 1)\n    def compute_squared_distance(X, Z):\n        return torch.cdist(X, Z) ** 2\n    \n    def __init__(self, lengthscale, variance):\n        super(RBF, self).__init__()\n        self.loglengthscale = nn.Parameter(torch.tensor(float(lengthscale)).log())\n        self.logvariance = nn.Parameter(torch.tensor(float(variance)).log())\n        \n    def __call__(self, X, Z=None):\n        if Z is None: \n            Z = X\n        return torch.exp(- 0.5 * torch.cdist(X, Z) ** 2 / (self.loglengthscale.exp() ** 2)) * self.logvariance.exp() \n    \n    def extra_repr(self):\n        return \"lengthscale = %.2f\\nvariance = %.2f\" % (self.loglengthscale.exp(), self.logvariance.exp())\n\n\nclass GPR(nn.Module):\n    @args_as_tensors(2, 3)\n    def __init__(self, kernel, noisevar=1e-3):\n        super(GPR, self).__init__() \n        self.kern = kernel\n        self.logvariance = torch.tensor(noisevar)#nn.Parameter(torch.tensor(noisevar).log()).to(device)\n    \n    def extra_repr(self):\n        return \"noise_variance = %1.2e\" % self.logvariance.exp()\n    \n    def compute_likelihood(self, X,y):        \n        K = self.kern(X) + torch.eye(X.shape[0]).to(device) * self.logvariance.exp()\n        #mask = K.diagonal() == 0\n        #K = torch.add(K, torch.diag(mask*1e-1))\n        L = torch.cholesky(K, upper=False)\n        V, _ = torch.solve(y.reshape(-1, 1), L)\n        ll = - 0.5 * (X.shape[0] * (np.log(2 * np.pi)) + torch.log(torch.diag(L)).sum() + (V**2).sum())\n        return ll\n   \n    @args_as_tensors(1)    \n    def predict_f_posterior(self, X,y,Xtest):\n        Kx = self.kern(X, Xtest)\n        K = self.kern(X) + torch.eye(X.shape[0]).to(device) * self.logvariance.exp()\n        #mask = K.diagonal() == 0\n        #K = torch.add(K, torch.diag(mask*1e-1))\n        L = torch.cholesky(K, upper=False)\n        A, _ = torch.solve(Kx, L)\n        V, _ = torch.solve(y.reshape(-1, 1), L)\n        fmean = A.T @ V\n        fvar = self.kern(Xtest) - A.T @ A\n        return fmean, fvar","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"class Dataset(torch.utils.data.Dataset):\n    \"\"\"\n    This hardcodes the mean and std values of the data\n    \"\"\"\n    def __init__(self, X, y, transform):\n        self.X = torch.tensor(transform.fit(X).transform(X)).type(torch.FloatTensor).to(device)\n        self.y = torch.tensor(y.values).type(torch.FloatTensor).to(device)\n\n    def __getitem__(self, index):\n        x,y = self.X[index], self.y[index]\n        return (x, y)\n  \n    def __len__(self):\n        return len(self.X)\n    \nclass TestDataset(torch.utils.data.Dataset):\n    def __init__(self, X, transform):\n        self.X = torch.tensor(transform.fit(X).transform(X)).type(torch.FloatTensor).to(device)\n\n    def __getitem__(self, index):\n        x = self.X[index]\n        return x\n  \n    def __len__(self):\n        return len(self.X)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"X, y = data.drop('ARRIVAL_DELAY', axis=1), data['ARRIVAL_DELAY']\nXtrain, Xval, ytrain, yval = train_test_split(X, y, test_size = 0.2, random_state = 42)\npipe = make_pipeline(StandardScaler(),\n                     SimpleImputer())\ntrain_loader = torch.utils.data.DataLoader(Dataset(Xtrain,ytrain, transform = pipe),batch_size = 256,shuffle = True, num_workers = 0)\nfit_loader = torch.utils.data.DataLoader(Dataset(Xtrain,ytrain, transform = pipe),batch_size = 5120,shuffle = True, num_workers = 0)\nval_loader = torch.utils.data.DataLoader(Dataset(Xval,yval, transform = pipe),batch_size = 1024,shuffle = True, num_workers = 0)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kernel = RBF(1, 1).to(device)\nmodel = GPR(kernel, noisevar = 1.).to(device)\nprint(model)\noptimizer = optim.Adam(model.parameters(), lr=.001)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"h_ll = []\nepochs = 1\nval_loss = [0]\nfor i in range(1,epochs+1):\n    epoch_loss = []\n    neg_rmse = np.mean(val_loss)\n    val_loss = []\n    for idx, (batch_train,batch_labels) in enumerate(train_loader,0):\n        optimizer.zero_grad()\n        loss = -model.compute_likelihood(batch_train,batch_labels)\n        loss.backward()\n        optimizer.step()\n        epoch_loss.append(loss.detach().cpu().numpy())\n        print('',end='\\r')\n        print(\"Epochs:[{}/{}] {}>{} train_loss: {} val_loss: {}\".format(i,epochs,\n                                                                        \"-\"*(idx//500),\"-\"*((len(train_loader)-idx-1)//500),\n                                                                        round(loss.detach().cpu().numpy()/(idx+1),3),round(neg_rmse,3)),end='')\n    mean_loss = np.mean(epoch_loss)\n    h_ll.append(mean_loss)\n    \n    fit_batch, fit_labels = next(iter(fit_loader))\n    for idx, (batch_val, val_labels) in enumerate(val_loader,0):\n        f_mean, f_cov = model.predict_f_posterior(fit_batch, fit_labels, batch_val)\n        f_mean, f_cov = f_mean.detach().cpu(), f_cov.detach().cpu().numpy()\n\n        val_loss.append(torch.sqrt(torch.mean(f_mean.flatten() - val_labels.cpu())**2))\n        \n    neg_rmse = np.mean(val_loss)\n    print('',end='\\r')\n    print(\"Epochs:[{}/{}] {}>{} train_loss: {} val_loss: {}\".format(i,epochs,\n                                                                        \"-\"*(idx//500),\"-\"*((len(train_loader)-idx-1)//500),\n                                                                        round(loss.detach().cpu().numpy()/(idx+1),3),round(neg_rmse,3)),end='')\ntorch.save(model,\"GP_model\")","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## GPFlow","metadata":{}},{"cell_type":"code","source":"X, y = data.drop('ARRIVAL_DELAY', axis=1), data['ARRIVAL_DELAY']\nXtrain, Xval, ytrain, yval = train_test_split(X, y, test_size = 0.2, random_state = 42)\npipe = make_pipeline(StandardScaler(),\n                     SimpleImputer())\nXtrain = pipe.fit(Xtrain).transform(Xtrain)\nXval = pipe.fit(Xval).transform(Xval)\n\nsc = StandardScaler()\nytrain = ytrain.values.reshape(-1,1)\nyval = yval.values.reshape(-1,1)\nytrain = sc.fit(ytrain).transform(ytrain)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\ndel X\ngc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"buffer_size = 1000\nbatch_size = 4000\n\n# Make dataset, do batching and shuffling\ndataset = tf.data.Dataset.from_tensor_slices((Xtrain, ytrain))\ndataset = dataset.shuffle(buffer_size)\ndataset = dataset.batch(batch_size, drop_remainder = True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"M = 1000 # Inducing points\nindexes = np.sort(np.random.randint(0, len(Xtrain),M))\nZ = Xtrain[indexes, :].copy() # Inducing points\nk = gpflow.kernels.RBF(lengthscales = np.random.uniform(0,20,Xtrain.shape[1]))#[]*Xtrain.shape[1])\nmean_f = None\nlikelihood = gpflow.likelihoods.Gaussian()\nm = gpflow.models.SVGP(kernel = k, mean_function = mean_f, inducing_variable = Z, likelihood = likelihood)\n#m = gpflow.models.SGPR(data = (Xtrain, ytrain), kernel = k, mean_function = mean_f, inducing_variable = Z, noise_variance = 0.1)\n#m.likelihood.variance.assign(0.3)\nset_trainable(likelihood.variance, True)\nset_trainable(k.lengthscales, True)\nset_trainable(m.inducing_variable, True)\nprint_summary(m)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def run_adam(model, lr, iterations):\n    \"\"\"\n    Utility function running the Adam optimizer\n\n    :param model: GPflow model\n    :param interations: number of iterations\n    \"\"\"\n    # Create an Adam Optimizer action\n    logf = []\n    train_iter = iter(dataset)\n    training_loss = model.training_loss_closure(train_iter, compile=True)\n    optimizer = tf.optimizers.Adam(learning_rate=lr)\n\n    @tf.function\n    def optimization_step():\n        optimizer.minimize(training_loss, model.trainable_variables)\n\n    for step in range(iterations):\n        print('',end='\\r')\n        print(\"Iteration: {}/{}\".format(step+1,iterations),end='')\n        optimization_step()\n        if step % 10 == 0:\n            elbo = -training_loss().numpy()\n            logf.append(elbo)\n    return logf\n\ndef optimize_model_with_scipy(model):\n    opt = gpflow.optimizers.Scipy()\n    training_loss = model.training_loss_closure(iter(dataset))\n    opt.minimize(\n        training_loss,\n        variables = model.trainable_variables,\n        #options=dict(maxiter=100),\n    )\n    \n#optimize_model_with_scipy(m)\nlogf = run_adam(m,0.03,960)\nprint_summary(m)\nplt.figure(figsize=(21,7))\nplt.plot(range(len(logf)), logf, marker = 'x')\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"indexes = np.sort(np.random.randint(0, len(Xval),200000))\nmean, var = m.predict_f(Xval[indexes,:])\nmean = sc.inverse_transform(mean)\nprint(mean_squared_error(mean, yval[indexes],squared=False))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"corrects = abs(1 * np.sqrt(var)) >= abs(mean - yval[indexes])\nprint(corrects.mean(), (1 * np.sqrt(var)).mean())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tf.saved_model.save(m, '')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"loaded_model = tf.saved_model.load('')","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import gc\ngc.collect()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Xtest = preprocess_data(test_df, airports_df, weather=weather_df, realistic_fts=False)\nids = Xtest['ID']\nXtest = Xtest.drop('ID', axis = 1)\n#add info on airport avg delays\nXtest = pd.merge(Xtest, delay_origin, how='left', left_on='LATITUDE_ORIGIN', right_on = 'LATITUDE_ORIGIN')\nXtest = pd.merge(Xtest, delay_dest, how='left', left_on='LATITUDE_DESTINATION', right_on = 'LATITUDE_DESTINATION')\nXtest = pipe.fit(Xtest).transform(Xtest)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with open(\"predsGPFin.txt\", \"a+\") as f:\n    f.write(\"ID,ARRIVAL_DELAY\\n\")\nfor i in range(0, Xtest.shape[0], 1000):\n    batch = Xtest[i:i+1000,:]\n    ids_p = ids[i:i+1000]\n    mean, var = m.predict_f(batch)\n    mean = sc.inverse_transform(mean)\n    with open(\"predsGPFin.txt\",\"a+\") as f:\n        for id_p,val in zip(ids_p,mean):\n            f.write(\"{},{}\\n\".format(id_p, val[0]))\n    print('',end='\\r')\n    print(\"Processed: {}\".format(i+1000),end='')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!kaggle competitions submit -c eurecom-aml-2021-challenge-1 -f predsGPFin.txt","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## XGBoost","metadata":{}},{"cell_type":"code","source":"params_cv = {\n    'n_estimators':[100,200,300,500],\n    'max_leaves':[1000,3000],\n    'min_child_weight': [10,15,20],\n    'gamma': [2, 2.5, 3],\n    'subsample': [0.8],\n    'colsample_bytree': [0.8],\n    'learning_rate':[0.1],  \n}\n\nmodel = xgboost.XGBRegressor(tree_method='gpu_hist', gpu_id=0,eval_metric='rmse',\n                             max_depth = 8, n_estimators = 500, #max_leaves = 1000,\n                             gamma = 2.5, learning_rate = 0.1, subsample = 0.8,\n                            colsample_bytree = 0.8)\npipe = make_pipeline(StandardScaler(),\n                     SimpleImputer())\n'''\npipe = make_pipeline(StandardScaler(),\n                     SimpleImputer(),\n                     GridSearchCV(\n                        estimator = model,\n                        param_grid = params_cv,\n                        scoring = 'neg_root_mean_squared_error',\n                        verbose=5,\n                        n_jobs=1\n                        )\n                    )\n'''\n\nX, y = data.drop('ARRIVAL_DELAY', axis=1).values, data['ARRIVAL_DELAY'].values\nXtrain, Xval, ytrain, yval = train_test_split(X, y, test_size = 0.2, random_state = 42)\nindexes = np.sort(np.random.randint(0, len(Xtrain), 3500000))\nXtrain, ytrain = Xtrain[indexes, :], ytrain[indexes]\nXtrain = pipe.fit(Xtrain).transform(Xtrain)\nXval = pipe.fit(Xval).transform(Xval)\nmodel.fit(Xtrain,ytrain)\n\n#pipe['gridsearchcv'].best_score_, pipe['gridsearchcv'].best_params_","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.figure(figsize = (30,5))\nplt.bar(data.drop('ARRIVAL_DELAY',axis=1).columns, model.feature_importances_)\nplt.xticks(rotation= 90)\nplt.savefig(\"priori_feature_importance.png\")\nplt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"indexes = np.sort(np.random.randint(0, len(Xval),200000))\npreds = model.predict(Xval[indexes,:])\nprint(mean_squared_error(preds, yval[indexes],squared=False))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"Xtest = preprocess_data(test_df, airports_df, weather=weather_df, realistic_fts=False)\nids = Xtest['ID']\nXtest = Xtest.drop('ID', axis = 1)\n#add info on airport avg delays\nXtest = pd.merge(Xtest, delay_origin, how='left', left_on='LATITUDE_ORIGIN', right_on = 'LATITUDE_ORIGIN')\nXtest = pd.merge(Xtest, delay_dest, how='left', left_on='LATITUDE_DESTINATION', right_on = 'LATITUDE_DESTINATION')\n#Xtest = pipe.fit(Xtest).transform(Xtest)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with open(\"xgb.txt\", \"w+\") as f:\n    f.write(\"ID,ARRIVAL_DELAY\\n\")\nfor i in range(0, Xtest.shape[0], 1000):\n    batch = Xtest.values[i:i+1000,:]\n    ids_p = ids[i:i+1000]\n    mean = model.predict(batch)\n    with open(\"xgb.txt\",\"a+\") as f:\n        for id_p,val in zip(ids_p,mean):\n            f.write(\"{},{}\\n\".format(id_p, val))\n    print('',end='\\r')\n    print(\"Processed: {}\".format(i+1000),end='')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Plot geographic data","metadata":{}},{"cell_type":"code","source":"from mpl_toolkits.basemap import Basemap\nimport matplotlib\n\ndef plot_geo_delay(df, latmin, latmax, longmin, longmax, title):\n    norm = matplotlib.colors.Normalize(vmin=df.values.min(), vmax=df.values.max())\n    cmap = matplotlib.cm.get_cmap('seismic')\n\n    fig, ax = plt.subplots(figsize=[6, 3], dpi=300)\n\n    m = Basemap(projection='cyl',\n                llcrnrlat=latmin, urcrnrlat=latmax,\n                llcrnrlon=longmin, urcrnrlon=longmax)\n    m.shadedrelief()\n    m.drawcoastlines()\n    m.drawcountries(linewidth = 0.3)\n    m.drawstates(color='0.3')\n\n\n    for index, delay in df.items():\n        ax.scatter(index[1],index[0], color = cmap(norm(delay)), s = 2.8,edgecolors = 'k',linewidths=0.3)\n\n\n    fig.colorbar(matplotlib.cm.ScalarMappable(norm=norm, cmap=cmap), ax=ax)\n    plt.savefig(title)\n    plt.show()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"delay_origin = data.groupby(['LATITUDE_ORIGIN','LONGITUDE_ORIGIN'])['ARRIVAL_DELAY'].aggregate(np.mean)\ndelay_dest = data.groupby(['LATITUDE_DESTINATION','LONGITUDE_DESTINATION'])['ARRIVAL_DELAY'].aggregate(np.mean)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plot_geo_delay(delay_origin, data['LATITUDE_ORIGIN'].min()-0.3, data['LATITUDE_ORIGIN'].max()+0.3, data['LONGITUDE_ORIGIN'].min()-0.3, data['LONGITUDE_ORIGIN'].max()+0.3, 'origindelay.png')\nplot_geo_delay(delay_dest, data['LATITUDE_DESTINATION'].min()-0.3, data['LATITUDE_DESTINATION'].max()+0.3, data['LONGITUDE_DESTINATION'].min()-0.3, data['LONGITUDE_DESTINATION'].max()+0.3,'destinationdelay.png')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}