{"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":"code","source":"from pathlib import Path","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":"import pandas as pd\ndata_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\ntruths = (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":"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":"for i in tqdm(range(1,len(df_test))):\n    lat1 = df_test.loc[i-1,'latDeg']\n    lon1 = df_test.loc[i-1,'lngDeg']\n    lat2 = df_test.loc[i,'latDeg']\n    lon2 = df_test.loc[i,'lngDeg']\n    df_test.loc[i,'dist_pre'] = calc_haversine(lat1, lon1, lat2, lon2)\n\nlist_phone = df_test['phone'].unique()\nfor phone in list_phone:\n    ind = df_test[df_test['phone'] == phone].index[0]\n    df_test.loc[ind,'dist_pre'] = 0","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:46:25.729039Z","iopub.execute_input":"2021-05-28T10:46:25.72941Z","iopub.status.idle":"2021-05-28T10:47:26.716564Z","shell.execute_reply.started":"2021-05-28T10:46:25.72938Z","shell.execute_reply":"2021-05-28T10:47:26.715603Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in tqdm(range(0,len(df_test)-1)):\n    lat1 = df_test.loc[i,'latDeg']\n    lon1 = df_test.loc[i,'lngDeg']\n    lat2 = df_test.loc[i+1,'latDeg']\n    lon2 = df_test.loc[i+1,'lngDeg']\n    df_test.loc[i,'dist_pro'] = calc_haversine(lat1, lon1, lat2, lon2)\nlist_phone = df_test['phone'].unique()\nfor phone in list_phone:\n    ind = df_test[df_test['phone'] == phone].index[-1]\n    df_test.loc[ind,'dist_pro'] = 0","metadata":{"execution":{"iopub.status.busy":"2021-05-28T10:47:27.398396Z","iopub.execute_input":"2021-05-28T10:47:27.398807Z","iopub.status.idle":"2021-05-28T10:48:28.229393Z","shell.execute_reply.started":"2021-05-28T10:47:27.398764Z","shell.execute_reply":"2021-05-28T10:48:28.22818Z"},"trusted":true},"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, 1, T, 0, 0.5 * T ** 2, 0], [0, 1, 0, T, 0, 0.5 * T ** 2], [0, 1, 1, 0, T, 0],\n                             [0, 1, 0, 1, 0, T], [0, 1, 0, 0, 1, 0], [0, 1, 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-10, 5e-10]) + np.ones((2, 2)) * 1e-9\n#ecept 2nd all 6 has 2 on 2nd vlaue 5e-5.\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"}}}]}