{"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":"Credit for the entire code for Kalman filter: https://www.kaggle.com/emaerthin/demonstration-of-the-kalman-filter\n\nI figure out that the hyperparameters have some effect on the result, so I refactor the code a little to parameterize all functions and search through it with Bayesian Optimization. The Kalman filter is treated as a black box that we can evaluate each separate data point.","metadata":{}},{"cell_type":"markdown","source":"## Install Kalman filter library","metadata":{}},{"cell_type":"code","source":"!pip install simdkalman","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Please read the documentation if you would like to learn more about this implementation of kf: https://simdkalman.readthedocs.io/en/latest/\n\nscikit-optimize documentation: https://scikit-optimize.github.io/stable/","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":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define kf model","metadata":{}},{"cell_type":"markdown","source":"T, size, noise (transition noise) and obs_noise (observation noise) will be the hyperparameters we try to tune.","metadata":{}},{"cell_type":"code","source":"T = 1.0\nsize = 4\nnoise = 1e-5\nobs_noise = 5e-5\n\ndef make_shifted_matrix(vec):\n    matrix = []\n    size = len(vec)\n    for i in range(size):\n        row = [0] * i + vec[:size-i]\n        matrix.append(row)\n    return np.array(matrix)\n\ndef 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\n\ndef 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\n\ndef make_kalman_filter(T, size, noise, obs_noise):\n    vec = make_state_vector(T, size)\n    state_transition = make_shifted_matrix(vec)\n    process_noise = np.diag(make_noise_vector(noise, size)) + np.ones(size) * 1e-9\n    observation_model = np.array([[1] + [0] * (size - 1), [0, 1] + [0] * (size - 2)])\n    observation_noise = np.diag([obs_noise] * 2) + np.ones(2) * 1e-9\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    return kf\n","metadata":{"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    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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Evaluate train error¶","metadata":{}},{"cell_type":"code","source":"data_path = Path(\"../input/google-smartphone-decimeter-challenge\")\n\ntruths = (data_path / 'train').rglob('ground_truth.csv')\n    # returns a generator\n\ndf_list = []\ncols = ['collectionName', 'phoneName', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg']\n\ndef calculate_location(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    df_basepreds_kf = apply_kf_smoothing(pd.read_csv('../input/google-smartphone-decimeter-challenge/baseline_locations_train.csv', usecols=cols), kf_=kf)\n    df_all = df_truth.merge(df_basepreds_kf, how='inner', on=cols[:3], suffixes=('_truth', '_basepred'))\n    return df_all","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# simplified haversine distance\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    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.0)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon/2.0)**2\n\n    c = 2 * np.arcsin(a**0.5)\n    dist = 6_367_000 * c\n    return dist","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_error(truths, kf):\n    df_all = calculate_location(truths, kf)\n    df_all['dist'] = calc_haversine(df_all.latDeg_truth, df_all.lngDeg_truth, \n        df_all.latDeg_basepred, df_all.lngDeg_basepred)\n    error = df_all.dist.mean()\n    return error","metadata":{"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    return error","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Hyperparameter tuning using Bayesian Optimization\n\nWe will use Gaussian Process of the package scikit-optimize for this. The objective is the mean error on the training set. The number of iterations is changed with the parameters n_calls. ","metadata":{}},{"cell_type":"code","source":"space = [Real(0.5, 1.5, name='T'), Integer(1, 4, name='half_size'), Real(1e-7, 1e-4, \"log-uniform\", name='noise'), Real(1e-7, 1e-4, \"log-uniform\", name='obs_noise')]\n\n# change \nresult = gp_minimize(optimize, space, n_calls=10)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"T, half_size, noise, obs_noise = result.x\nsize = half_size * 2\nprint(f'Best result: T = {T}, size = {size}, noise = {noise} => error = {result.fun:.3f}m')","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Print the best 5 models","metadata":{}},{"cell_type":"code","source":"num_models = 5\nbest_param_indices = sorted(range(len(result.func_vals)), key=lambda x: result.func_vals[x])[:num_models]\nbest_params = [result.x_iters[i] for i in best_param_indices]\nbest_params","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Prepare a submission based on the sample submission\n\nThe best hyperparameter below is the result of a search with 200 iterations. I also tried to ensemble several best models into one but the result on the test set is worse than only use the best model. Naive approach with mean or median ensembling doesn't yield better result.","metadata":{}},{"cell_type":"code","source":"best_params = [[1.5, 2, 2.1714133956113952e-06, 1.719317114286542e-05]]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_base = pd.read_csv(\n'../input/google-smartphone-decimeter-challenge/baseline_locations_test.csv')\nsub = pd.read_csv('../input/google-smartphone-decimeter-challenge/sample_submission.csv')\nlatDeg = []\nlngDeg = []\nfor i in range(len(best_params)):\n    T, half_size, noise, obs_noise = best_params[i]\n    size = half_size * 2\n    kf = make_kalman_filter(T, size, noise, obs_noise)\n    kf_smoothed_baseline = apply_kf_smoothing(test_base, kf)\n    latDeg.append(kf_smoothed_baseline.latDeg)\n    lngDeg.append(kf_smoothed_baseline.lngDeg)\nsub = sub.assign(\nlatDeg = np.mean(latDeg, axis =  0),\nlngDeg = np.mean(lngDeg, axis =  0)\n)\nsub.to_csv('submission.csv', index=False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]}]}