{"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":"!pip install simdkalman","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true},"execution_count":null,"outputs":[]},{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"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\ndef 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":"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\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    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\n\ndef 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\n\ndef 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":"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\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":"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":"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":[]},{"cell_type":"code","source":"!","metadata":{},"execution_count":null,"outputs":[]}]}