{"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":"# This Python 3 environment comes with many helpful analytics libraries installed\n# It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# For example, here's several helpful packages to load\n\nimport numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport matplotlib.pyplot as plt\nfrom matplotlib_venn import venn2, venn2_circles\nimport seaborn as sns\nfrom tqdm.notebook import tqdm\nfrom pathlib import Path\nimport plotly\nimport plotly.express as px\n\n# Input data files are available in the read-only \"../input/\" directory\n# For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\nimport os\nfor dirname, _, filenames in os.walk('/kaggle/input'):\n    for filename in filenames:\n        print(os.path.join(dirname, filename))\n\n# You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-07-19T22:22:56.280328Z","iopub.execute_input":"2021-07-19T22:22:56.280852Z","iopub.status.idle":"2021-07-19T22:22:59.220669Z","shell.execute_reply.started":"2021-07-19T22:22:56.280806Z","shell.execute_reply":"2021-07-19T22:22:59.219632Z"},"collapsed":true,"jupyter":{"outputs_hidden":true},"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-07-19T22:22:59.222397Z","iopub.execute_input":"2021-07-19T22:22:59.222700Z","iopub.status.idle":"2021-07-19T22:22:59.587893Z","shell.execute_reply.started":"2021-07-19T22:22:59.222670Z","shell.execute_reply":"2021-07-19T22:22:59.586571Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"truths = (data_path / 'train').rglob('ground_truth.csv')\ndf_list = []\ncols = ['collectionName', 'phoneName', 'millisSinceGpsEpoch', 'latDeg',\n       'lngDeg']\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)\ndf_phone.head()\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'))\n#df_all.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-19T22:22:59.589826Z","iopub.execute_input":"2021-07-19T22:22:59.590175Z","iopub.status.idle":"2021-07-19T22:23:01.391811Z","shell.execute_reply.started":"2021-07-19T22:22:59.590142Z","shell.execute_reply":"2021-07-19T22:23:01.390786Z"},"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-07-19T22:23:01.393699Z","iopub.execute_input":"2021-07-19T22:23:01.394158Z","iopub.status.idle":"2021-07-19T22:23:01.401631Z","shell.execute_reply.started":"2021-07-19T22:23:01.394111Z","shell.execute_reply":"2021-07-19T22:23:01.400561Z"},"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)\ndf_all.dist.describe()\ndf_all.sort_values(by = 'dist',ascending = False)[['collectionName','dist']].head(10)","metadata":{"execution":{"iopub.status.busy":"2021-07-19T22:23:01.403179Z","iopub.execute_input":"2021-07-19T22:23:01.403469Z","iopub.status.idle":"2021-07-19T22:23:01.535640Z","shell.execute_reply.started":"2021-07-19T22:23:01.403441Z","shell.execute_reply":"2021-07-19T22:23:01.534906Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test['dist_pre'] = 0\ndf_test['dist_pro'] = 0\ndf_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\ndf_test.dist_pre.describe()","metadata":{"execution":{"iopub.status.busy":"2021-07-19T22:23:01.536735Z","iopub.execute_input":"2021-07-19T22:23:01.537146Z","iopub.status.idle":"2021-07-19T22:23:03.025270Z","shell.execute_reply.started":"2021-07-19T22:23:01.537117Z","shell.execute_reply":"2021-07-19T22:23:03.024539Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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-07-19T22:23:03.026225Z","iopub.execute_input":"2021-07-19T22:23:03.026610Z","iopub.status.idle":"2021-07-19T22:23:03.184939Z","shell.execute_reply.started":"2021-07-19T22:23:03.026570Z","shell.execute_reply":"2021-07-19T22:23:03.184154Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"!pip install simdkalman\nimport numpy as np\nimport pandas as pd\nimport simdkalman\nfrom tqdm.notebook import tqdm","metadata":{"execution":{"iopub.status.busy":"2021-07-19T22:23:03.186611Z","iopub.execute_input":"2021-07-19T22:23:03.187005Z","iopub.status.idle":"2021-07-19T22:23:12.374579Z","shell.execute_reply.started":"2021-07-19T22:23:03.186976Z","shell.execute_reply":"2021-07-19T22:23:12.373531Z"},"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-07-19T22:23:12.376902Z","iopub.execute_input":"2021-07-19T22:23:12.377376Z","iopub.status.idle":"2021-07-19T22:23:12.387818Z","shell.execute_reply.started":"2021-07-19T22:23:12.377310Z","shell.execute_reply":"2021-07-19T22:23:12.386645Z"},"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-07-19T22:23:12.389552Z","iopub.execute_input":"2021-07-19T22:23:12.390028Z","iopub.status.idle":"2021-07-19T22:23:12.401058Z","shell.execute_reply.started":"2021-07-19T22:23:12.389982Z","shell.execute_reply":"2021-07-19T22:23:12.400194Z"},"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-07-19T22:23:12.402176Z","iopub.execute_input":"2021-07-19T22:23:12.402461Z","iopub.status.idle":"2021-07-19T22:23:46.901957Z","shell.execute_reply.started":"2021-07-19T22:23:12.402434Z","shell.execute_reply":"2021-07-19T22:23:46.900845Z"},"trusted":true},"execution_count":null,"outputs":[]}]}