{"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\nI thought a simple and easy post-processing based on my experience of using google maps on a daily basis.  \nHave you ever had an experience where the google map location information was way off?  I thought that the data set of this competition might contain such outliers.\n\n\nIn this notebook, Correct the outliers in the baseline and apply the Kalman filter afterwards.\n\nScore(Public Leaderboard)  \nOriginal : 7.190  \noutlier correction : 7.180  \nKalman filter : 6.164  ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"markdown","source":"### ■ Reference\nI used these great notebooks as a guide in creating this notebook. Thank you for publishing the notebook.\n\n\n(Baseline :by JohnM)\nhttps://www.kaggle.com/jpmiller/baseline-from-host-data\n\n(Kalman filter : by Marcin Bodych)\nhttps://www.kaggle.com/emaerthin/demonstration-of-the-kalman-filter","metadata":{}},{"cell_type":"markdown","source":"### ■ Update (2021-06-02)\nI got some advice from @Kyosuke0924 .I changed some code to speed up the process.","metadata":{}},{"cell_type":"markdown","source":"### Step.1 Check for outliers in training data --------------------------------------------","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\nimport pandas as pd","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:00.760399Z","iopub.execute_input":"2021-05-28T10:46:00.760844Z","iopub.status.idle":"2021-05-28T10:46:00.764964Z","shell.execute_reply.started":"2021-05-28T10:46:00.760803Z","shell.execute_reply":"2021-05-28T10:46:00.764166Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_path = Path(\"../input/google-smartphone-decimeter-challenge\")\ndf_test = pd.read_csv(\n    data_path / 'baseline_locations_test.csv')\ndf_sub    = pd.read_csv(\n    data_path / 'sample_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:01.451614Z","iopub.execute_input":"2021-05-28T10:46:01.452027Z","iopub.status.idle":"2021-05-28T10:46:02.28517Z","shell.execute_reply.started":"2021-05-28T10:46:01.451993Z","shell.execute_reply":"2021-05-28T10:46:02.284036Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"import numpy as np\nfrom tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:03.304866Z","iopub.execute_input":"2021-05-28T10:46:03.305302Z","iopub.status.idle":"2021-05-28T10:46:03.309575Z","shell.execute_reply.started":"2021-05-28T10:46:03.305239Z","shell.execute_reply":"2021-05-28T10:46:03.308579Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"truths = (data_path / 'train').rglob('ground_truth.csv')\n    # returns a generator\n\ndf_list = []\ncols = ['collectionName', 'phoneName', 'millisSinceGpsEpoch', 'latDeg',\n       'lngDeg']\n\nfor t in tqdm(truths, total=73):\n    df_phone = pd.read_csv(t, usecols=cols)  \n    df_list.append(df_phone)\ndf_truth = pd.concat(df_list, ignore_index=True)\n\ndf_basepreds = pd.read_csv(data_path / 'baseline_locations_train.csv', usecols=cols)\ndf_all = df_truth.merge(df_basepreds, how='inner', on=cols[:3], suffixes=('_truth', '_basepred'))","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:03.986136Z","iopub.execute_input":"2021-05-28T10:46:03.986525Z","iopub.status.idle":"2021-05-28T10:46:05.003581Z","shell.execute_reply.started":"2021-05-28T10:46:03.986492Z","shell.execute_reply":"2021-05-28T10:46:05.002643Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:05.00484Z","iopub.execute_input":"2021-05-28T10:46:05.005255Z","iopub.status.idle":"2021-05-28T10:46:05.010958Z","shell.execute_reply.started":"2021-05-28T10:46:05.005224Z","shell.execute_reply":"2021-05-28T10:46:05.010261Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_all['dist'] = calc_haversine(df_all.latDeg_truth, df_all.lngDeg_truth, \n    df_all.latDeg_basepred, df_all.lngDeg_basepred)","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:05.333487Z","iopub.execute_input":"2021-05-28T10:46:05.333888Z","iopub.status.idle":"2021-05-28T10:46:05.366325Z","shell.execute_reply.started":"2021-05-28T10:46:05.333851Z","shell.execute_reply":"2021-05-28T10:46:05.3655Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_all.dist.describe()","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:05.967811Z","iopub.execute_input":"2021-05-28T10:46:05.968419Z","iopub.status.idle":"2021-05-28T10:46:05.987753Z","shell.execute_reply.started":"2021-05-28T10:46:05.968367Z","shell.execute_reply":"2021-05-28T10:46:05.98648Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_all.sort_values(by = 'dist',ascending = False)[['collectionName','dist']].head(10)","metadata":{"execution":{"iopub.status.busy":"2021-05-28T11:07:00.248104Z","iopub.execute_input":"2021-05-28T11:07:00.248511Z","iopub.status.idle":"2021-05-28T11:07:00.299834Z","shell.execute_reply.started":"2021-05-28T11:07:00.248476Z","shell.execute_reply":"2021-05-28T11:07:00.298766Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There are some points where the difference between baseline and ground truth is more than 1,000 meters.It is surprising.  \nDealing with these outliers is likely to improve the score.","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:23:31.014242Z","iopub.execute_input":"2021-05-28T10:23:31.014665Z","iopub.status.idle":"2021-05-28T10:23:31.022602Z","shell.execute_reply.started":"2021-05-28T10:23:31.014631Z","shell.execute_reply":"2021-05-28T10:23:31.020862Z"}}},{"cell_type":"markdown","source":"### Step.2 Correct outliers in the test data.----------------------------------------------","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:35:32.622185Z","iopub.execute_input":"2021-05-28T10:35:32.622722Z","iopub.status.idle":"2021-05-28T10:35:32.641461Z","shell.execute_reply.started":"2021-05-28T10:35:32.62268Z","shell.execute_reply":"2021-05-28T10:35:32.640763Z"}}},{"cell_type":"markdown","source":"Since there is no 'dist' in the test data, identify outliers based on their distance from the previous and next pass.","metadata":{}},{"cell_type":"code","source":"df_test['dist_pre'] = 0\ndf_test['dist_pro'] = 0","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:19.650898Z","iopub.execute_input":"2021-05-28T10:46:19.651314Z","iopub.status.idle":"2021-05-28T10:46:19.658459Z","shell.execute_reply.started":"2021-05-28T10:46:19.651265Z","shell.execute_reply":"2021-05-28T10:46:19.657192Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test['latDeg_pre'] = df_test['latDeg'].shift(periods=1,fill_value=0)\ndf_test['lngDeg_pre'] = df_test['lngDeg'].shift(periods=1,fill_value=0)\ndf_test['latDeg_pro'] = df_test['latDeg'].shift(periods=-1,fill_value=0)\ndf_test['lngDeg_pro'] = df_test['lngDeg'].shift(periods=-1,fill_value=0)\ndf_test['dist_pre'] = calc_haversine(df_test.latDeg_pre, df_test.lngDeg_pre, df_test.latDeg, df_test.lngDeg)\ndf_test['dist_pro'] = calc_haversine(df_test.latDeg, df_test.lngDeg, df_test.latDeg_pro, df_test.lngDeg_pro)\n\nlist_phone = df_test['phone'].unique()\nfor phone in list_phone:\n    ind_s = df_test[df_test['phone'] == phone].index[0]\n    ind_e = df_test[df_test['phone'] == phone].index[-1]\n    df_test.loc[ind_s,'dist_pre'] = 0\n    df_test.loc[ind_e,'dist_pro'] = 0","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test.dist_pre.describe()","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:48:51.583968Z","iopub.execute_input":"2021-05-28T10:48:51.584351Z","iopub.status.idle":"2021-05-28T10:48:51.597069Z","shell.execute_reply.started":"2021-05-28T10:48:51.58431Z","shell.execute_reply":"2021-05-28T10:48:51.596383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There seems to be an outlier, although it is smaller than the training data.","metadata":{}},{"cell_type":"markdown","source":"If the distance between the previous and next pass is large, it is either an outlier or a case of speed (straight line movement), so fix the value to the middle of the previous and next pass.","metadata":{}},{"cell_type":"code","source":"pro_95 = df_test['dist_pro'].mean() + (df_test['dist_pro'].std() * 2)\npre_95 = df_test['dist_pre'].mean() + (df_test['dist_pre'].std() * 2)\nind = df_test[(df_test['dist_pro'] > pro_95)&(df_test['dist_pre'] > pre_95)][['dist_pre','dist_pro']].index\n\nfor i in ind:\n    df_test.loc[i,'latDeg'] = (df_test.loc[i-1,'latDeg'] + df_test.loc[i+1,'latDeg'])/2\n    df_test.loc[i,'lngDeg'] = (df_test.loc[i-1,'lngDeg'] + df_test.loc[i+1,'lngDeg'])/2","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:53:18.483732Z","iopub.execute_input":"2021-05-28T10:53:18.484066Z","iopub.status.idle":"2021-05-28T10:53:18.648391Z","shell.execute_reply.started":"2021-05-28T10:53:18.484037Z","shell.execute_reply":"2021-05-28T10:53:18.64729Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Step.3 Kalman filter -----------------------------------------------------------------","metadata":{}},{"cell_type":"markdown","source":"Finally, apply the Kalman filter.","metadata":{}},{"cell_type":"code","source":"!pip install simdkalman","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:56:42.365687Z","iopub.execute_input":"2021-05-28T10:56:42.366127Z","iopub.status.idle":"2021-05-28T10:56:51.262943Z","shell.execute_reply.started":"2021-05-28T10:56:42.36609Z","shell.execute_reply":"2021-05-28T10:56:51.261929Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pathlib import Path\nimport numpy as np\nimport pandas as pd\nimport simdkalman\nfrom tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:56:51.265081Z","iopub.execute_input":"2021-05-28T10:56:51.265863Z","iopub.status.idle":"2021-05-28T10:56:51.276607Z","shell.execute_reply.started":"2021-05-28T10:56:51.265809Z","shell.execute_reply":"2021-05-28T10:56:51.275831Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"T = 1.0\nstate_transition = np.array([[1, 0, T, 0, 0.5 * T ** 2, 0], [0, 1, 0, T, 0, 0.5 * T ** 2], [0, 0, 1, 0, T, 0],\n                             [0, 0, 0, 1, 0, T], [0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1]])\nprocess_noise = np.diag([1e-5, 1e-5, 5e-6, 5e-6, 1e-6, 1e-6]) + np.ones((6, 6)) * 1e-9\nobservation_model = np.array([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]])\nobservation_noise = np.diag([5e-5, 5e-5]) + np.ones((2, 2)) * 1e-9\n\nkf = simdkalman.KalmanFilter(\n        state_transition = state_transition,\n        process_noise = process_noise,\n        observation_model = observation_model,\n        observation_noise = observation_noise)","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:56:53.364882Z","iopub.execute_input":"2021-05-28T10:56:53.365403Z","iopub.status.idle":"2021-05-28T10:56:53.374266Z","shell.execute_reply.started":"2021-05-28T10:56:53.36537Z","shell.execute_reply":"2021-05-28T10:56:53.373384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def apply_kf_smoothing(df, kf_=kf):\n    unique_paths = df[['collectionName', 'phoneName']].drop_duplicates().to_numpy()\n    for collection, phone in tqdm(unique_paths):\n        cond = np.logical_and(df['collectionName'] == collection, df['phoneName'] == phone)\n        data = df[cond][['latDeg', 'lngDeg']].to_numpy()\n        data = data.reshape(1, len(data), 2)\n        smoothed = kf_.smooth(data)\n        df.loc[cond, 'latDeg'] = smoothed.states.mean[0, :, 0]\n        df.loc[cond, 'lngDeg'] = smoothed.states.mean[0, :, 1]\n    return df","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:56:58.47818Z","iopub.execute_input":"2021-05-28T10:56:58.478944Z","iopub.status.idle":"2021-05-28T10:56:58.486798Z","shell.execute_reply.started":"2021-05-28T10:56:58.478897Z","shell.execute_reply":"2021-05-28T10:56:58.485848Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"kf_smoothed_baseline = apply_kf_smoothing(df_test)\ndf_sub = df_sub.assign(\n    latDeg = kf_smoothed_baseline.latDeg,\n    lngDeg = kf_smoothed_baseline.lngDeg\n)\ndf_sub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-05-28T11:10:19.987637Z","iopub.execute_input":"2021-05-28T11:10:19.988006Z","iopub.status.idle":"2021-05-28T11:10:54.757708Z","shell.execute_reply.started":"2021-05-28T11:10:19.987966Z","shell.execute_reply":"2021-05-28T11:10:54.756723Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"#### Thank you for reading to the end.I look forward to your comments.","metadata":{"execution":{"iopub.status.busy":"2021-05-28T11:10:54.759402Z","iopub.execute_input":"2021-05-28T11:10:54.759931Z","iopub.status.idle":"2021-05-28T11:10:54.765259Z","shell.execute_reply.started":"2021-05-28T11:10:54.759878Z","shell.execute_reply":"2021-05-28T11:10:54.764041Z"}}}]}