{"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":"# Overview\n\nThis notebook is designed to give you an introduction on how to approach this competition and use gnss data. I perform outlier correction and apply a savgol filter after hyperparameter tuning with bayesian optimization.\n\nThis notebook is broken down into a few sections. \n1. Standard Functions and Constants - this code is mostly helper functions borrowed from the notebook by @saitodevel01. It is used to generate the baseline which I have included as a datasource to save time. This section also contains my imports and evaluation function.\n2. Outlier Correction - here I detect outliers by comparing the lat and lon at each timestep to the timestep before and after. If the haversine distance between the points is greater than a threshold, it is flagged as an outlier. I then replace outliers with the mean of the lat and lon at the previous and future timestep.\n3. Savgol Filter - here I have defined a function to apply scipy’s savgol filter algorithm to the lat and lon columns. The function is set up to hyperparameter tune the window length and poly order. \n4. Bayesian Optimization - here I use skopt’s gp_minimize function in order to apply Bayesian optimization using Gaussian Processes. I optimize the outlier correction threshold, savgol filter window length, and savgol filter poly order\n5. Submit - uses optimal parameters to generate submission file\n\n**References**\n\nhttps://www.kaggle.com/code/saitodevel01/gsdc2-baseline-submission by @saitodevel01 - used for baseline generation\n\nhttps://www.kaggle.com/code/dehokanta/baseline-post-processing-by-outlier-correction by @dehokanta - notebook from last year’s competition inspired outlier correction technique\n\nhttps://www.kaggle.com/code/tqa236/kalman-filter-hyperparameter-search-with-bo by @tqa236 - notebook from last year’s competition inspiration for bayesian optimization","metadata":{}},{"cell_type":"markdown","source":"# Standard Functions and Constants","metadata":{}},{"cell_type":"code","source":"import glob\nfrom dataclasses import dataclass\nimport numpy as np\nimport pandas as pd\nfrom tqdm.notebook import tqdm\nfrom pathlib import Path\n\nfrom scipy.interpolate import InterpolatedUnivariateSpline\nfrom scipy.signal import savgol_filter\n\nfrom skopt import gp_minimize\nfrom skopt.space import Real, Integer\n\nimport warnings\nwarnings.filterwarnings('ignore')\n\nINPUT_PATH = '../input/smartphone-decimeter-2022'\nbl_path = '../input/gsdc2-baseline-submission'\nbl_train = pd.read_csv(f'{bl_path}/baseline_train.csv')\nbl_test = pd.read_csv(f'{bl_path}/baseline_test.csv')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-05-07T01:34:17.876943Z","iopub.execute_input":"2022-05-07T01:34:17.877314Z","iopub.status.idle":"2022-05-07T01:34:20.286036Z","shell.execute_reply.started":"2022-05-07T01:34:17.877223Z","shell.execute_reply":"2022-05-07T01:34:20.285265Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"WGS84_SEMI_MAJOR_AXIS = 6378137.0\nWGS84_SEMI_MINOR_AXIS = 6356752.314245\nWGS84_SQUARED_FIRST_ECCENTRICITY  = 6.69437999013e-3\nWGS84_SQUARED_SECOND_ECCENTRICITY = 6.73949674226e-3\n\nHAVERSINE_RADIUS = 6_371_000","metadata":{"execution":{"iopub.status.busy":"2022-05-07T01:34:20.287604Z","iopub.execute_input":"2022-05-07T01:34:20.288235Z","iopub.status.idle":"2022-05-07T01:34:20.292620Z","shell.execute_reply.started":"2022-05-07T01:34:20.288203Z","shell.execute_reply":"2022-05-07T01:34:20.291830Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# reference https://www.kaggle.com/code/saitodevel01/gsdc2-baseline-submission\n\n@dataclass\nclass ECEF:\n    x: np.array\n    y: np.array\n    z: np.array\n\n    def to_numpy(self):\n        return np.stack([self.x, self.y, self.z], axis=0)\n\n    @staticmethod\n    def from_numpy(pos):\n        x, y, z = [np.squeeze(w) for w in np.split(pos, 3, axis=-1)]\n        return ECEF(x=x, y=y, z=z)\n\n@dataclass\nclass BLH:\n    lat : np.array\n    lng : np.array\n    hgt : np.array\n\ndef ECEF_to_BLH(ecef):\n    a = WGS84_SEMI_MAJOR_AXIS\n    b = WGS84_SEMI_MINOR_AXIS\n    e2  = WGS84_SQUARED_FIRST_ECCENTRICITY\n    e2_ = WGS84_SQUARED_SECOND_ECCENTRICITY\n    x = ecef.x\n    y = ecef.y\n    z = ecef.z\n    r = np.sqrt(x**2 + y**2)\n    t = np.arctan2(z * (a/b), r)\n    B = np.arctan2(z + (e2_*b)*np.sin(t)**3, r - (e2*a)*np.cos(t)**3)\n    L = np.arctan2(y, x)\n    n = a / np.sqrt(1 - e2*np.sin(B)**2)\n    H = (r / np.cos(B)) - n\n    return BLH(lat=B, lng=L, hgt=H)\n\ndef haversine_distance(blh_1, blh_2):\n    dlat = blh_2.lat - blh_1.lat\n    dlng = blh_2.lng - blh_1.lng\n    a = np.sin(dlat/2)**2 + np.cos(blh_1.lat) * np.cos(blh_2.lat) * np.sin(dlng/2)**2\n    dist = 2 * HAVERSINE_RADIUS * np.arcsin(np.sqrt(a))\n    return dist\n\ndef pandas_haversine_distance(df1, df2):\n    blh1 = BLH(\n        lat=np.deg2rad(df1['LatitudeDegrees'].to_numpy()),\n        lng=np.deg2rad(df1['LongitudeDegrees'].to_numpy()),\n        hgt=0,\n    )\n    blh2 = BLH(\n        lat=np.deg2rad(df2['LatitudeDegrees'].to_numpy()),\n        lng=np.deg2rad(df2['LongitudeDegrees'].to_numpy()),\n        hgt=0,\n    )\n    return haversine_distance(blh1, blh2)","metadata":{"execution":{"iopub.status.busy":"2022-05-07T01:34:20.293950Z","iopub.execute_input":"2022-05-07T01:34:20.294312Z","iopub.status.idle":"2022-05-07T01:34:20.314511Z","shell.execute_reply.started":"2022-05-07T01:34:20.294263Z","shell.execute_reply":"2022-05-07T01:34:20.313711Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def calc_score(tripID, pred_df, gt_df):\n    d = pandas_haversine_distance(pred_df, gt_df)\n    score = np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])    \n    return score","metadata":{"execution":{"iopub.status.busy":"2022-05-07T01:34:20.317274Z","iopub.execute_input":"2022-05-07T01:34:20.318253Z","iopub.status.idle":"2022-05-07T01:34:20.331235Z","shell.execute_reply.started":"2022-05-07T01:34:20.318207Z","shell.execute_reply":"2022-05-07T01:34:20.330618Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Outlier Correction","metadata":{}},{"cell_type":"code","source":"def calc_haversine(lat1, lon1, lat2, lon2):\n    \"\"\"Calculates the great circle distance between two points\n    on the earth. Inputs are array-like and specified in decimal degrees.\n    \"\"\"\n    RADIUS = 6_367_000\n    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])\n    dlat = lat2 - lat1\n    dlon = lon2 - lon1\n    a = np.sin(dlat/2)**2 + \\\n        np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2)**2\n    dist = 2 * RADIUS * np.arcsin(a**0.5)\n    return dist\n\ndef correct_outliers(df, th=2):\n    df['dist_pre'] = 0\n    df['dist_pro'] = 0\n\n    df['latDeg_pre'] = df['LatitudeDegrees'].shift(periods=1,fill_value=0)\n    df['lngDeg_pre'] = df['LongitudeDegrees'].shift(periods=1,fill_value=0)\n    df['latDeg_pro'] = df['LatitudeDegrees'].shift(periods=-1,fill_value=0)\n    df['lngDeg_pro'] = df['LongitudeDegrees'].shift(periods=-1,fill_value=0)\n    df['dist_pre'] = calc_haversine(df.latDeg_pre, df.lngDeg_pre, df.LatitudeDegrees, df.LongitudeDegrees)\n    df['dist_pro'] = calc_haversine(df.LatitudeDegrees, df.LongitudeDegrees, df.latDeg_pro, df.lngDeg_pro)\n\n    df.loc[df.index.min(), 'dist_pre'] = 0\n    df.loc[df.index.max(), 'dist_pro'] = 0\n    \n    pro_95 = df['dist_pro'].mean() + (df['dist_pro'].std() * th)\n    pre_95 = df['dist_pre'].mean() + (df['dist_pre'].std() * th)\n\n    ind = df[(df['dist_pro'] > pro_95)&(df['dist_pre'] > pre_95)][['dist_pre','dist_pro']].index\n\n    for i in ind:\n        df.loc[i,'LatitudeDegrees'] = (df.loc[i-1,'LatitudeDegrees'] + df.loc[i+1,'LatitudeDegrees'])/2\n        df.loc[i,'LongitudeDegrees'] = (df.loc[i-1,'LongitudeDegrees'] + df.loc[i+1,'LongitudeDegrees'])/2\n    \n    return df","metadata":{"execution":{"iopub.status.busy":"2022-05-07T01:34:20.332490Z","iopub.execute_input":"2022-05-07T01:34:20.332880Z","iopub.status.idle":"2022-05-07T01:34:20.348001Z","shell.execute_reply.started":"2022-05-07T01:34:20.332834Z","shell.execute_reply":"2022-05-07T01:34:20.347155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Savgol Filter","metadata":{}},{"cell_type":"code","source":"def apply_savgol_filter(df, wl, poly):\n    df.LatitudeDegrees = savgol_filter(df.LatitudeDegrees, wl, poly)\n    df.LongitudeDegrees = savgol_filter(df.LongitudeDegrees, wl, poly)\n    return df","metadata":{"execution":{"iopub.status.busy":"2022-05-07T01:34:20.349263Z","iopub.execute_input":"2022-05-07T01:34:20.349725Z","iopub.status.idle":"2022-05-07T01:34:20.363577Z","shell.execute_reply.started":"2022-05-07T01:34:20.349660Z","shell.execute_reply":"2022-05-07T01:34:20.362862Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Bayesian Optimization","metadata":{}},{"cell_type":"code","source":"def optimize(params):\n    th, wl, poly = params\n    if wl%2==0:\n        wl+=1\n    \n    score_list = []\n\n    for tripID in sorted(bl_train.tripId.unique()):\n\n        gt_df   = pd.read_csv(f'{INPUT_PATH}/train/{tripID}/ground_truth.csv')\n        pred_df = bl_train[bl_train.tripId == tripID]\n\n        pred_df = correct_outliers(pred_df, th)\n        pred_df = apply_savgol_filter(pred_df, wl, poly)\n\n        score = calc_score(tripID, pred_df, gt_df)\n        score_list.append(score)\n\n    mean_score = np.mean(score_list)\n    return mean_score","metadata":{"execution":{"iopub.status.busy":"2022-05-07T01:34:20.364924Z","iopub.execute_input":"2022-05-07T01:34:20.365783Z","iopub.status.idle":"2022-05-07T01:34:20.376070Z","shell.execute_reply.started":"2022-05-07T01:34:20.365743Z","shell.execute_reply":"2022-05-07T01:34:20.375480Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"space = [Real(1.5, 2.5, name='threshhold'), \n         Integer(7, 31, name='window_len'), \n         Integer(2, 6, name='poly_order')]\n\nresult = gp_minimize(optimize, space, n_calls=100)","metadata":{"execution":{"iopub.status.busy":"2022-05-07T01:34:20.377436Z","iopub.execute_input":"2022-05-07T01:34:20.378266Z","iopub.status.idle":"2022-05-07T01:57:29.832025Z","shell.execute_reply.started":"2022-05-07T01:34:20.378233Z","shell.execute_reply":"2022-05-07T01:57:29.831168Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'best train score: {result.fun}')","metadata":{"execution":{"iopub.status.busy":"2022-05-07T01:57:29.833418Z","iopub.execute_input":"2022-05-07T01:57:29.833739Z","iopub.status.idle":"2022-05-07T01:57:29.839935Z","shell.execute_reply.started":"2022-05-07T01:57:29.833699Z","shell.execute_reply":"2022-05-07T01:57:29.839159Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"if result.x[1]%2==0:\n    result.x[1]+=1\n\nprint(f'best params:\\noutlier threshhold: {result.x[0]}\\nsavgol filter window length: {result.x[1]}\\nsavgol filter poly order: {result.x[2]}')","metadata":{"execution":{"iopub.status.busy":"2022-05-07T02:17:05.293045Z","iopub.execute_input":"2022-05-07T02:17:05.293366Z","iopub.status.idle":"2022-05-07T02:17:05.299061Z","shell.execute_reply.started":"2022-05-07T02:17:05.293333Z","shell.execute_reply":"2022-05-07T02:17:05.298364Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Submit","metadata":{}},{"cell_type":"code","source":"preds = list()\n\nfor tripID in sorted(bl_test.tripId.unique()):\n    pred_df = bl_test[bl_test.tripId == tripID]\n\n    pred_df = correct_outliers(pred_df, result.x[0])\n    pred_df = apply_savgol_filter(pred_df, result.x[1], result.x[2])\n\n    preds.append(pred_df)\n    \nsub = pd.concat(preds)\nsub = sub[[\"tripId\", \"UnixTimeMillis\", \"LatitudeDegrees\", \"LongitudeDegrees\"]]\nsub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2022-05-07T02:19:01.806408Z","iopub.execute_input":"2022-05-07T02:19:01.806977Z","iopub.status.idle":"2022-05-07T02:19:04.150180Z","shell.execute_reply.started":"2022-05-07T02:19:01.806932Z","shell.execute_reply":"2022-05-07T02:19:04.149269Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}