{"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":"# 1.0 Introduction","metadata":{}},{"cell_type":"code","source":"## put something here to toggle between kg env and home env","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:43.538643Z","iopub.execute_input":"2021-06-03T00:24:43.538962Z","iopub.status.idle":"2021-06-03T00:24:43.543144Z","shell.execute_reply.started":"2021-06-03T00:24:43.538932Z","shell.execute_reply":"2021-06-03T00:24:43.542121Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 2.0 Installs and Imports","metadata":{}},{"cell_type":"markdown","source":"## 2.1 simdkalman installation","metadata":{"heading_collapsed":true}},{"cell_type":"code","source":"!pip install simdkalman","metadata":{"hidden":true,"execution":{"iopub.status.busy":"2021-06-03T00:24:43.54824Z","iopub.execute_input":"2021-06-03T00:24:43.548608Z","iopub.status.idle":"2021-06-03T00:24:49.170015Z","shell.execute_reply.started":"2021-06-03T00:24:43.548579Z","shell.execute_reply":"2021-06-03T00:24:49.169155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2.2 Imports","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\nimport numpy as np\nimport pandas as pd\nimport simdkalman\nfrom tqdm.notebook import tqdm\nimport itertools\nfrom skopt import gp_minimize\nfrom skopt.space import Real, Integer","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.171311Z","iopub.execute_input":"2021-06-03T00:24:49.171538Z","iopub.status.idle":"2021-06-03T00:24:49.17698Z","shell.execute_reply.started":"2021-06-03T00:24:49.171514Z","shell.execute_reply":"2021-06-03T00:24:49.176214Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 3.0 Classes and Functions","metadata":{}},{"cell_type":"markdown","source":"## 3.1 Kalman Filter Model","metadata":{}},{"cell_type":"markdown","source":"T, size, transitions noise and observation noise to be tuned","metadata":{}},{"cell_type":"code","source":"T=1.0\nsize = 4\nnoise = 1e-5\nobs_noise = 5e-5","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.178261Z","iopub.execute_input":"2021-06-03T00:24:49.178554Z","iopub.status.idle":"2021-06-03T00:24:49.18775Z","shell.execute_reply.started":"2021-06-03T00:24:49.178526Z","shell.execute_reply":"2021-06-03T00:24:49.18704Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_shifted_matrix(vector):\n    matrix =[]\n    vlen = len(vector)\n    for i in range(vlen):\n        matrix.append([0] * i + vector[:vlen-i])\n    \n    return np.array(matrix)","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.188528Z","iopub.execute_input":"2021-06-03T00:24:49.188777Z","iopub.status.idle":"2021-06-03T00:24:49.19893Z","shell.execute_reply.started":"2021-06-03T00:24:49.188755Z","shell.execute_reply":"2021-06-03T00:24:49.198219Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_state_vector(T, size):\n    vector = [1,0]\n    step = 2\n    for i in range (size - 2):\n        if i % 2 == 0:\n            vector.append(T)\n            T *= T / step\n            step += 1\n        else:\n            vector.append(0)\n    return vector","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.202261Z","iopub.execute_input":"2021-06-03T00:24:49.202786Z","iopub.status.idle":"2021-06-03T00:24:49.209706Z","shell.execute_reply.started":"2021-06-03T00:24:49.202752Z","shell.execute_reply":"2021-06-03T00:24:49.208926Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_noise_vector(noise, size):\n    noise_vector = []\n    for i in range(size):\n        if i>0 and i % 2 == 0:\n            noise *= 0.5\n        noise_vector.append(noise)\n    return noise_vector","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.21078Z","iopub.execute_input":"2021-06-03T00:24:49.21117Z","iopub.status.idle":"2021-06-03T00:24:49.224085Z","shell.execute_reply.started":"2021-06-03T00:24:49.211139Z","shell.execute_reply":"2021-06-03T00:24:49.223353Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def make_kalman_filter(T, size, noise, obs_noise):\n    vec = make_state_vector(T, size)\n    #print(vec)\n    state_transition = make_shifted_matrix(vec)\n    #print(state_transition)\n    \n    process_noise = (np.diag(make_noise_vector(noise,size))\n                     + np.ones(size) * 1e-9)\n    #print(process_noise)\n    \n    observation_model = np.array([[1]+[0]*(size-1), \n                                 [0,1] + [0]*(size-2)])\n    #print(observation_model)\n    \n    observation_noise = (np.diag([obs_noise]*2) +\n                         np.ones(2) * 1e-9)\n    \n    #print(observation_noise)\n    \n    kf = simdkalman.KalmanFilter(\n            state_transition = state_transition,\n            process_noise = process_noise,\n            observation_model = observation_model,\n            observation_noise = observation_noise)\n    \n    return kf","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.225239Z","iopub.execute_input":"2021-06-03T00:24:49.225709Z","iopub.status.idle":"2021-06-03T00:24:49.235435Z","shell.execute_reply.started":"2021-06-03T00:24:49.225668Z","shell.execute_reply":"2021-06-03T00:24:49.234685Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def apply_kf_smoothing(df, kf_):\n    unique_paths = df[['collectionName', 'phoneName']].drop_duplicates().to_numpy()\n    \n    for collection, phone in unique_paths:\n        #XXX this part feels unneccessary, examine later\n        cond = np.logical_and(df['collectionName'] == collection,\n                              df['phoneName'] == phone)\n        data = df[cond][['latDeg', 'lngDeg']].to_numpy()\n        data = data.reshape(1,len(data), 2)\n        \n        #so  think he is reshaping it to make the values rows\n        # then the \n        \n        #apply the kalman filter and overwrite\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    \n    return df","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.236418Z","iopub.execute_input":"2021-06-03T00:24:49.236643Z","iopub.status.idle":"2021-06-03T00:24:49.245734Z","shell.execute_reply.started":"2021-06-03T00:24:49.236621Z","shell.execute_reply":"2021-06-03T00:24:49.245067Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.2 Ingest Data and apply filter","metadata":{}},{"cell_type":"code","source":"#get a generator to run through all those nested files\n\ndata_path = Path('../input/google-smartphone-decimeter-challenge/')\ntruths = (data_path / 'train').rglob('ground_truth.csv')","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.246821Z","iopub.execute_input":"2021-06-03T00:24:49.247286Z","iopub.status.idle":"2021-06-03T00:24:49.259283Z","shell.execute_reply.started":"2021-06-03T00:24:49.247245Z","shell.execute_reply":"2021-06-03T00:24:49.258508Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_list = []\ncols = ['collectionName', \n        'phoneName',\n        'millisSinceGpsEpoch',\n        'latDeg',\n        'lngDeg']","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.260296Z","iopub.execute_input":"2021-06-03T00:24:49.260842Z","iopub.status.idle":"2021-06-03T00:24:49.268247Z","shell.execute_reply.started":"2021-06-03T00:24:49.260804Z","shell.execute_reply":"2021-06-03T00:24:49.267572Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"#takes and iterator of csv's and a kalman fitler\n\ndef calculate_locations(truths, kf):\n    \n    for t in truths:\n        df_phone = pd.read_csv(t, usecols=cols)\n        df_list.append(df_phone)\n    df_truth = pd.concat(df_list, ignore_index=True)\n    \n    baselines = '../input/google-smartphone-decimeter-challenge/baseline_locations_train.csv'\n    df_basepreds_kf = apply_kf_smoothing(pd.read_csv(baselines, \n                                                     usecols=cols),\n                                         kf_ = kf)\n    #don't like the on parameters, not clear XXX\n    df_all = df_truth.merge(df_basepreds_kf,\n                           how = 'inner',\n                           on = cols[:3], \n                           suffixes=('_truth','_basepred'))\n    \n    return df_all","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.26922Z","iopub.execute_input":"2021-06-03T00:24:49.26956Z","iopub.status.idle":"2021-06-03T00:24:49.277608Z","shell.execute_reply.started":"2021-06-03T00:24:49.269535Z","shell.execute_reply":"2021-06-03T00:24:49.276938Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.3 haversine distance","metadata":{}},{"cell_type":"code","source":"#great circle distance on earth \n#need to tear this a part later :)\ndef calc_haversine(lat1, lon1, lat2, lon2):\n    lat1, lon1, lat2, lon2 = map(np.radians, [lat1, lon1, lat2, lon2])\n    \n    dlat = lat2 - lat1\n    dlon = lon2 - lon1\n    \n    a = (np.sin(dlat/2.0)**2 + \n         np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2 )\n    \n    c = 2 * np.arcsin(a**0.5)\n    \n    return  6367000 * c\n    ","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.279713Z","iopub.execute_input":"2021-06-03T00:24:49.279988Z","iopub.status.idle":"2021-06-03T00:24:49.287612Z","shell.execute_reply.started":"2021-06-03T00:24:49.279967Z","shell.execute_reply":"2021-06-03T00:24:49.286989Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3.4 error calcs and optimization","metadata":{}},{"cell_type":"code","source":"def get_error(truths, kf):\n    df_all = calculate_locations(truths, kf)\n    df_all['dist'] = calc_haversine(df_all.latDeg_truth,\n                                    df_all.lngDeg_truth,\n                                    df_all.latDeg_basepred,\n                                    df_all.lngDeg_basepred)\n    error = df_all.dist.mean()\n    return error","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.288725Z","iopub.execute_input":"2021-06-03T00:24:49.289098Z","iopub.status.idle":"2021-06-03T00:24:49.300651Z","shell.execute_reply.started":"2021-06-03T00:24:49.28905Z","shell.execute_reply":"2021-06-03T00:24:49.300046Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def optimize(params):\n    T, half_size, noise, obs_noise = params\n    size = half_size * 2\n    kf = make_kalman_filter(T, size, noise, obs_noise)\n    error = get_error(truths, kf)\n    print(f'T = {T}, size = {size}, noise = {noise}, obs_noise = {obs_noise} => error = {error:.3f}m')\n    \n    return error","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.301574Z","iopub.execute_input":"2021-06-03T00:24:49.301784Z","iopub.status.idle":"2021-06-03T00:24:49.30946Z","shell.execute_reply.started":"2021-06-03T00:24:49.301764Z","shell.execute_reply":"2021-06-03T00:24:49.308867Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 4.0 Hyperparameter turning","metadata":{}},{"cell_type":"code","source":"space = [Real(0.5, 1.5, name='T'), \n         Integer(1, 4, name='half_size'), \n         Real(1e-7, 1e-4, \"log-uniform\", name='noise'), \n         Real(1e-7, 1e-4, \"log-uniform\", name='obs_noise')]\n\nresults = gp_minimize(optimize, space, n_calls=10 )","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:24:49.31015Z","iopub.execute_input":"2021-06-03T00:24:49.310355Z","iopub.status.idle":"2021-06-03T00:31:46.232737Z","shell.execute_reply.started":"2021-06-03T00:24:49.310335Z","shell.execute_reply":"2021-06-03T00:31:46.23174Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"num_models = 4\nbest_param_indices = sorted(range(len(results.func_vals)), \n                        key = lambda x: results.func_vals[x])[:num_models]\n\nbest_params = [results.x_iters[i] for i in best_param_indices]\n\nbest_params","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:31:46.233712Z","iopub.execute_input":"2021-06-03T00:31:46.233987Z","iopub.status.idle":"2021-06-03T00:31:46.242971Z","shell.execute_reply.started":"2021-06-03T00:31:46.233961Z","shell.execute_reply":"2021-06-03T00:31:46.242054Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 5.0 Submission preparation","metadata":{}},{"cell_type":"code","source":"test_base_path = '../input/google-smartphone-decimeter-challenge/baseline_locations_test.csv'\nsub_path = '../input/google-smartphone-decimeter-challenge/sample_submission.csv'\n\ntest_base = pd.read_csv(test_base_path)\nsub = pd.read_csv(sub_path)\n\n\nT, half_size, noise, obs_noise = best_params[0]\nsize = half_size * 2\nkf = make_kalman_filter(T, size, noise, obs_noise)\nkf_smoothed_baseline = apply_kf_smoothing(test_base, kf)\n\nlatDeg = kf_smoothed_baseline.latDeg\nlngDeg = kf_smoothed_baseline.lngDeg\n\n\nprint(latDeg[0], len(latDeg))\n\nprint(lngDeg[0], len(lngDeg))\n\nprint(len(sub))\nsub = sub.assign(\n    latDeg = latDeg,\n    lngDeg = lngDeg\n)\n\nsub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:33:24.704793Z","iopub.execute_input":"2021-06-03T00:33:24.705159Z","iopub.status.idle":"2021-06-03T00:33:55.125303Z","shell.execute_reply.started":"2021-06-03T00:33:24.70513Z","shell.execute_reply":"2021-06-03T00:33:55.124617Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(len(sub))","metadata":{"execution":{"iopub.status.busy":"2021-06-03T00:31:46.282691Z","iopub.status.idle":"2021-06-03T00:31:46.283149Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"","metadata":{},"execution_count":null,"outputs":[]}]}