{"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":"# 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\n\ndef 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 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)\n    \n# directory setting\nINPUT = '../input/google-smartphone-decimeter-challenge'\nbase_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')\n\n# 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())\n\n\ndef 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\n\n\n# 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\n\n!pip install simdkalman\nimport simdkalman\n\nT = 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\n\ncols = ['collectionName', 'phoneName', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg']\ntrain_ro_kf = apply_kf_smoothing(train_ro[cols])\n\ndef 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]\n\ndef 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\n\ntrain_lerp = make_lerp_data(train_ro_kf)\ntrain_mean_pred = calc_mean_pred(train_ro_kf, train_lerp)\n\n\ntmp1 = 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')\n\ndef percentile50(x):\n    return np.percentile(x, 50)\ndef percentile95(x):\n    return np.percentile(x, 95)\n\ndef 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\n\nprint('kf + reject_outlier : ', get_train_score(train_ro_kf, ground_truth))\nprint('+ phones_mean_pred : ', get_train_score(train_mean_pred, ground_truth))\n\nbase_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)\n\n\n\n\n\n\n\n    \n    \n    \n    \n    \n    \n    ","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-01-27T18:54:02.194148Z","iopub.execute_input":"2022-01-27T18:54:02.194537Z","iopub.status.idle":"2022-01-27T18:55:30.745139Z","shell.execute_reply.started":"2022-01-27T18:54:02.194418Z","shell.execute_reply":"2022-01-27T18:55:30.744201Z"},"trusted":true},"execution_count":null,"outputs":[]}]}