{"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"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":26933,"databundleVersionId":2264643}],"dockerImageVersionId":30096,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"It's a very simple idea,  \nbut I thought it might be useful to use the average of the predictions of  \nseveral phones in the same collection as the final prediction.","metadata":{}},{"cell_type":"code","source":"# import library\nimport os\nimport pandas as pd\nimport numpy as np\nimport matplotlib.pyplot as plt\nfrom matplotlib_venn import venn2, venn2_circles\nimport seaborn as sns\nfrom tqdm.notebook import tqdm\nimport pathlib\nimport plotly\nimport plotly.express as px","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-06-02T22:21:03.341586Z","iopub.execute_input":"2021-06-02T22:21:03.341848Z","iopub.status.idle":"2021-06-02T22:21:03.34691Z","shell.execute_reply.started":"2021-06-02T22:21:03.341824Z","shell.execute_reply":"2021-06-02T22:21:03.345505Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# utils","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","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:03.348263Z","iopub.execute_input":"2021-06-02T22:21:03.348634Z","iopub.status.idle":"2021-06-02T22:21:03.367053Z","shell.execute_reply.started":"2021-06-02T22:21:03.348603Z","shell.execute_reply":"2021-06-02T22:21:03.366107Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def visualize_trafic(df, center, zoom=9):\n    fig = px.scatter_mapbox(df,\n                            \n                            # Here, plotly gets, (x,y) coordinates\n                            lat=\"latDeg\",\n                            lon=\"lngDeg\",\n                            \n                            #Here, plotly detects color of series\n                            color=\"phoneName\",\n                            labels=\"phoneName\",\n                            \n                            zoom=zoom,\n                            center=center,\n                            height=600,\n                            width=800)\n    fig.update_layout(mapbox_style='stamen-terrain')\n    fig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\n    fig.update_layout(title_text=\"GPS trafic\")\n    fig.show()\n    \ndef visualize_collection(df, collection):\n    target_df = df[df['collectionName']==collection].copy()\n    lat_center = target_df['latDeg'].mean()\n    lng_center = target_df['lngDeg'].mean()\n    center = {\"lat\":lat_center, \"lon\":lng_center}\n    \n    visualize_trafic(target_df, center)","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:03.368804Z","iopub.execute_input":"2021-06-02T22:21:03.369115Z","iopub.status.idle":"2021-06-02T22:21:03.379727Z","shell.execute_reply.started":"2021-06-02T22:21:03.369082Z","shell.execute_reply":"2021-06-02T22:21:03.378676Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# data prep","metadata":{}},{"cell_type":"code","source":"# directory setting\nINPUT = '../input/google-smartphone-decimeter-challenge'","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:03.38081Z","iopub.execute_input":"2021-06-02T22:21:03.381046Z","iopub.status.idle":"2021-06-02T22:21:03.392531Z","shell.execute_reply.started":"2021-06-02T22:21:03.381014Z","shell.execute_reply":"2021-06-02T22:21:03.391255Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"base_train = pd.read_csv(INPUT + '/' + 'baseline_locations_train.csv')\nbase_test = pd.read_csv(INPUT + '/' + 'baseline_locations_test.csv')\nsample_sub = pd.read_csv(INPUT + '/' + 'sample_submission.csv')","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:03.393894Z","iopub.execute_input":"2021-06-02T22:21:03.394216Z","iopub.status.idle":"2021-06-02T22:21:03.645342Z","shell.execute_reply.started":"2021-06-02T22:21:03.394183Z","shell.execute_reply":"2021-06-02T22:21:03.644379Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ground_truth\np = pathlib.Path(INPUT)\ngt_files = list(p.glob('train/*/*/ground_truth.csv'))\nprint('ground_truth.csv count : ', len(gt_files))\n\ngts = []\nfor gt_file in tqdm(gt_files):\n    gts.append(pd.read_csv(gt_file))\nground_truth = pd.concat(gts)\n\ndisplay(ground_truth.head())","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:03.646348Z","iopub.execute_input":"2021-06-02T22:21:03.646561Z","iopub.status.idle":"2021-06-02T22:21:04.192704Z","shell.execute_reply.started":"2021-06-02T22:21:03.646537Z","shell.execute_reply":"2021-06-02T22:21:04.191813Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# reject outlier","metadata":{}},{"cell_type":"code","source":"def add_distance_diff(df):\n    df['latDeg_prev'] = df['latDeg'].shift(1)\n    df['latDeg_next'] = df['latDeg'].shift(-1)\n    df['lngDeg_prev'] = df['lngDeg'].shift(1)\n    df['lngDeg_next'] = df['lngDeg'].shift(-1)\n    df['phone_prev'] = df['phone'].shift(1)\n    df['phone_next'] = df['phone'].shift(-1)\n    \n    df['dist_prev'] = calc_haversine(df['latDeg'], df['lngDeg'], df['latDeg_prev'], df['lngDeg_prev'])\n    df['dist_next'] = calc_haversine(df['latDeg'], df['lngDeg'], df['latDeg_next'], df['lngDeg_next'])\n    \n    df.loc[df['phone']!=df['phone_prev'], ['latDeg_prev', 'lngDeg_prev', 'dist_prev']] = np.nan\n    df.loc[df['phone']!=df['phone_next'], ['latDeg_next', 'lngDeg_next', 'dist_next']] = np.nan\n    \n    return df","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:04.193704Z","iopub.execute_input":"2021-06-02T22:21:04.193915Z","iopub.status.idle":"2021-06-02T22:21:04.201295Z","shell.execute_reply.started":"2021-06-02T22:21:04.193894Z","shell.execute_reply":"2021-06-02T22:21:04.200258Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# reject outlier\ntrain_ro = add_distance_diff(base_train)\nth = 50\ntrain_ro.loc[((train_ro['dist_prev'] > th) & (train_ro['dist_next'] > th)), ['latDeg', 'lngDeg']] = np.nan","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:04.30975Z","iopub.execute_input":"2021-06-02T22:21:04.310011Z","iopub.status.idle":"2021-06-02T22:21:04.320963Z","shell.execute_reply.started":"2021-06-02T22:21:04.309985Z","shell.execute_reply":"2021-06-02T22:21:04.319389Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# kalman filter\nhttps://www.kaggle.com/emaerthin/demonstration-of-the-kalman-filter","metadata":{}},{"cell_type":"code","source":"!pip install simdkalman","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:04.322341Z","iopub.execute_input":"2021-06-02T22:21:04.322625Z","iopub.status.idle":"2021-06-02T22:21:09.206364Z","shell.execute_reply.started":"2021-06-02T22:21:04.322594Z","shell.execute_reply":"2021-06-02T22:21:09.20527Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"import simdkalman","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:09.207792Z","iopub.execute_input":"2021-06-02T22:21:09.208079Z","iopub.status.idle":"2021-06-02T22:21:09.213118Z","shell.execute_reply.started":"2021-06-02T22:21:09.20805Z","shell.execute_reply":"2021-06-02T22:21:09.211351Z"},"trusted":true},"outputs":[],"execution_count":null},{"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)\n\ndef apply_kf_smoothing(df, kf_=kf):\n    unique_paths = df[['collectionName', 'phoneName']].drop_duplicates().to_numpy()\n    for collection, phone in 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-06-02T22:21:09.214434Z","iopub.execute_input":"2021-06-02T22:21:09.214736Z","iopub.status.idle":"2021-06-02T22:21:09.232292Z","shell.execute_reply.started":"2021-06-02T22:21:09.214706Z","shell.execute_reply":"2021-06-02T22:21:09.231335Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"cols = ['collectionName', 'phoneName', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg']\ntrain_ro_kf = apply_kf_smoothing(train_ro[cols])","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:21:09.233474Z","iopub.execute_input":"2021-06-02T22:21:09.233783Z","iopub.status.idle":"2021-06-02T22:22:24.135317Z","shell.execute_reply.started":"2021-06-02T22:21:09.233731Z","shell.execute_reply":"2021-06-02T22:22:24.134502Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# phones mean prediction","metadata":{}},{"cell_type":"code","source":"def make_lerp_data(df):\n    '''\n    Generate interpolated lat,lng values for different phone times in the same collection.\n    '''\n    org_columns = df.columns\n    \n    # Generate a combination of time x collection x phone and combine it with the original data (generate records to be interpolated)\n    time_list = df[['collectionName', 'millisSinceGpsEpoch']].drop_duplicates()\n    phone_list =df[['collectionName', 'phoneName']].drop_duplicates()\n    tmp = time_list.merge(phone_list, on='collectionName', how='outer')\n    \n    lerp_df = tmp.merge(df, on=['collectionName', 'millisSinceGpsEpoch', 'phoneName'], how='left')\n    lerp_df['phone'] = lerp_df['collectionName'] + '_' + lerp_df['phoneName']\n    lerp_df = lerp_df.sort_values(['phone', 'millisSinceGpsEpoch'])\n    \n    # linear interpolation\n    lerp_df['latDeg_prev'] = lerp_df['latDeg'].shift(1)\n    lerp_df['latDeg_next'] = lerp_df['latDeg'].shift(-1)\n    lerp_df['lngDeg_prev'] = lerp_df['lngDeg'].shift(1)\n    lerp_df['lngDeg_next'] = lerp_df['lngDeg'].shift(-1)\n    lerp_df['phone_prev'] = lerp_df['phone'].shift(1)\n    lerp_df['phone_next'] = lerp_df['phone'].shift(-1)\n    lerp_df['time_prev'] = lerp_df['millisSinceGpsEpoch'].shift(1)\n    lerp_df['time_next'] = lerp_df['millisSinceGpsEpoch'].shift(-1)\n    # Leave only records to be interpolated\n    lerp_df = lerp_df[(lerp_df['latDeg'].isnull())&(lerp_df['phone']==lerp_df['phone_prev'])&(lerp_df['phone']==lerp_df['phone_next'])].copy()\n    # calc lerp\n    lerp_df['latDeg'] = lerp_df['latDeg_prev'] + ((lerp_df['latDeg_next'] - lerp_df['latDeg_prev']) * ((lerp_df['millisSinceGpsEpoch'] - lerp_df['time_prev']) / (lerp_df['time_next'] - lerp_df['time_prev']))) \n    lerp_df['lngDeg'] = lerp_df['lngDeg_prev'] + ((lerp_df['lngDeg_next'] - lerp_df['lngDeg_prev']) * ((lerp_df['millisSinceGpsEpoch'] - lerp_df['time_prev']) / (lerp_df['time_next'] - lerp_df['time_prev']))) \n    \n    # Leave only the data that has a complete set of previous and next data.\n    lerp_df = lerp_df[~lerp_df['latDeg'].isnull()]\n    \n    return lerp_df[org_columns]","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:22:24.136377Z","iopub.execute_input":"2021-06-02T22:22:24.136924Z","iopub.status.idle":"2021-06-02T22:22:24.148521Z","shell.execute_reply.started":"2021-06-02T22:22:24.13689Z","shell.execute_reply":"2021-06-02T22:22:24.147578Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def calc_mean_pred(df, lerp_df):\n    '''\n    Make a prediction based on the average of the predictions of phones in the same collection.\n    '''\n    add_lerp = pd.concat([df, lerp_df])\n    mean_pred_result = add_lerp.groupby(['collectionName', 'millisSinceGpsEpoch'])[['latDeg', 'lngDeg']].mean().reset_index()\n    mean_pred_df = df[['collectionName', 'phoneName', 'millisSinceGpsEpoch']].copy()\n    mean_pred_df = mean_pred_df.merge(mean_pred_result[['collectionName', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg']], on=['collectionName', 'millisSinceGpsEpoch'], how='left')\n    return mean_pred_df","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:22:24.149982Z","iopub.execute_input":"2021-06-02T22:22:24.150513Z","iopub.status.idle":"2021-06-02T22:22:24.165724Z","shell.execute_reply.started":"2021-06-02T22:22:24.15048Z","shell.execute_reply":"2021-06-02T22:22:24.165028Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"train_lerp = make_lerp_data(train_ro_kf)\ntrain_mean_pred = calc_mean_pred(train_ro_kf, train_lerp)","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:22:24.167159Z","iopub.execute_input":"2021-06-02T22:22:24.167733Z","iopub.status.idle":"2021-06-02T22:22:24.974136Z","shell.execute_reply.started":"2021-06-02T22:22:24.167702Z","shell.execute_reply":"2021-06-02T22:22:24.973522Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"tmp1 = train_ro_kf.copy()\ntmp2 = train_mean_pred.copy()\ntmp2['phoneName'] = tmp2['phoneName'] + '_MEAN'\ntmp3 = ground_truth.copy()\ntmp3['phoneName'] = tmp3['phoneName'] + '_GT'\ntmp = pd.concat([tmp1, tmp2, tmp3])\nvisualize_collection(tmp, '2020-05-14-US-MTV-1')","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:22:24.975178Z","iopub.execute_input":"2021-06-02T22:22:24.975517Z","iopub.status.idle":"2021-06-02T22:22:25.153194Z","shell.execute_reply.started":"2021-06-02T22:22:24.975492Z","shell.execute_reply":"2021-06-02T22:22:25.152128Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# evaluate train score","metadata":{}},{"cell_type":"code","source":"def percentile50(x):\n    return np.percentile(x, 50)\ndef percentile95(x):\n    return np.percentile(x, 95)","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:22:25.154257Z","iopub.execute_input":"2021-06-02T22:22:25.154471Z","iopub.status.idle":"2021-06-02T22:22:25.15933Z","shell.execute_reply.started":"2021-06-02T22:22:25.154449Z","shell.execute_reply":"2021-06-02T22:22:25.158125Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"def get_train_score(df, gt):\n    gt = gt.rename(columns={'latDeg':'latDeg_gt', 'lngDeg':'lngDeg_gt'})\n    df = df.merge(gt, on=['collectionName', 'phoneName', 'millisSinceGpsEpoch'], how='inner')\n    # calc_distance_error\n    df['err'] = calc_haversine(df['latDeg_gt'], df['lngDeg_gt'], df['latDeg'], df['lngDeg'])\n    # calc_evaluate_score\n    df['phone'] = df['collectionName'] + '_' + df['phoneName']\n    res = df.groupby('phone')['err'].agg([percentile50, percentile95])\n    res['p50_p90_mean'] = (res['percentile50'] + res['percentile95']) / 2 \n    score = res['p50_p90_mean'].mean()\n    return score","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:22:25.160643Z","iopub.execute_input":"2021-06-02T22:22:25.160962Z","iopub.status.idle":"2021-06-02T22:22:25.175912Z","shell.execute_reply.started":"2021-06-02T22:22:25.160933Z","shell.execute_reply":"2021-06-02T22:22:25.175038Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"code","source":"print('kf + reject_outlier : ', get_train_score(train_ro_kf, ground_truth))\nprint('+ phones_mean_pred : ', get_train_score(train_mean_pred, ground_truth))","metadata":{"execution":{"iopub.status.busy":"2021-06-02T22:22:25.176968Z","iopub.execute_input":"2021-06-02T22:22:25.177192Z","iopub.status.idle":"2021-06-02T22:22:25.805584Z","shell.execute_reply.started":"2021-06-02T22:22:25.177168Z","shell.execute_reply":"2021-06-02T22:22:25.804537Z"},"trusted":true},"outputs":[],"execution_count":null},{"cell_type":"markdown","source":"# make submission","metadata":{}},{"cell_type":"code","source":"base_test = add_distance_diff(base_test)\nth = 50\nbase_test.loc[((base_test['dist_prev'] > th) & (base_test['dist_next'] > th)), ['latDeg', 'lngDeg']] = np.nan\n\ntest_kf = apply_kf_smoothing(base_test)\n\ntest_lerp = make_lerp_data(test_kf)\ntest_mean_pred = calc_mean_pred(test_kf, test_lerp)\n\nsample_sub['latDeg'] = test_mean_pred['latDeg']\nsample_sub['lngDeg'] = test_mean_pred['lngDeg']\nsample_sub.to_csv('submission.csv', index=False)","metadata":{},"outputs":[],"execution_count":null}]}