{"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 numpy as np # linear algebra\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nimport os\nfrom pathlib import Path\nimport pyproj\nfrom pyproj import Proj, transform\n\ndef calc_haversine(lat1, lon1, lat2, lon2):\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 compute_dist(fname, fname2 = 'gt.csv'):\n    oof = pd.read_csv(fname)\n    gt = pd.read_csv(fname2)\n    df = oof.merge(gt, on = ['phone','millisSinceGpsEpoch'])\n    dst_oof = calc_haversine(df.latDeg_x,df.lngDeg_x, df.latDeg_y, df.lngDeg_y)\n    scores = pd.DataFrame({'phone': df.phone,'dst': dst_oof})\n    scores_grp = scores.groupby('phone')\n    d50 = scores_grp.quantile(.50).reset_index()\n    d50.columns = ['phone','q50']\n    d95 = scores_grp.quantile(.95).reset_index()\n    d95.columns = ['phone','q95']\n    return (scores_grp.quantile(.50).mean() + scores_grp.quantile(.95).mean())/2, d50.merge(d95)\n\ndef WGS84_to_ECEF(lat, lon, alt):\n    # convert to radians\n    rad_lat = lat * (np.pi / 180.0)\n    rad_lon = lon * (np.pi / 180.0)\n    a    = 6378137.0\n    # f is the flattening factor\n    finv = 298.257223563\n    f = 1 / finv   \n    # e is the eccentricity\n    e2 = 1 - (1 - f) * (1 - f)    \n    # N is the radius of curvature in the prime vertical\n    N = a / np.sqrt(1 - e2 * np.sin(rad_lat) * np.sin(rad_lat))\n    x = (N + alt) * np.cos(rad_lat) * np.cos(rad_lon)\n    y = (N + alt) * np.cos(rad_lat) * np.sin(rad_lon)\n    z = (N * (1 - e2) + alt)        * np.sin(rad_lat)\n    return x, y, z\n\ntransformer = pyproj.Transformer.from_crs(\n    {\"proj\":'geocent', \"ellps\":'WGS84', \"datum\":'WGS84'},\n    {\"proj\":'latlong', \"ellps\":'WGS84', \"datum\":'WGS84'},)\ndef ECEF_to_WGS84(x,y,z):\n    lon, lat, alt = transformer.transform(x,y,z,radians=False)\n    return lon, lat, alt\n\ndatadir = Path('/kaggle/input/google-smartphone-decimeter-challenge/')\ntestdir = datadir / 'test'\ntraindir = datadir / 'train'\n\nsample_sub = pd.read_csv(datadir/'sample_submission.csv')\nsub_columns = sample_sub.columns\n\nbaseline_train = pd.read_csv(datadir / 'baseline_locations_train.csv')\nbaseline_train[sub_columns].to_csv('btrain.csv',index = False)\nbaseline_test = pd.read_csv(datadir / 'baseline_locations_test.csv')\nbaseline_test[sub_columns].to_csv('btest.csv',index = False)\n\nmsge = 'millisSinceGpsEpoch'\n\ngt = pd.DataFrame()\nfor d in os.listdir(traindir):\n    for p in os.listdir(traindir/d):\n        gt = gt.append(pd.read_csv(traindir/d/p/'ground_truth.csv'))\n\ngt['phone'] = gt['collectionName'] + '_' + gt['phoneName']\ngt[sub_columns].to_csv('gt.csv', index = False)\ngt['heightAboveWgs84EllipsoidM'].describe()\n\n\nscore, scores = compute_dist('btrain.csv','gt.csv')\nprint(score)\nprint(scores)\n\n\nimport optuna\ndef position_shift(fname,a):\n\n    d = pd.read_csv(fname)\n    d['heightAboveWgs84EllipsoidM'] = 63.5\n    d['x'], d['y'], d['z'] = zip(*d.apply(lambda x: WGS84_to_ECEF(x.latDeg, x.lngDeg, x.heightAboveWgs84EllipsoidM), axis=1))\n\n    #a = -0.2\n    d.sort_values(['phone', msge], inplace=True)\n    for fi in ['x','y','z']:\n        d[[fi+'p']] = d[fi].shift().where(d['phone'].eq(d['phone'].shift()))\n        d[[fi+'diff']] = d[fi]-d[fi+'p']\n    #d[['yp']] = d['y'].shift().where(d['phone'].eq(d['phone'].shift()))\n    d[['dist']] = np.sqrt(d['xdiff']**2 + d['ydiff']**2+ d['zdiff']**2)\n    for fi in ['x','y','z']:\n        d[[fi+'new']] = d[fi+'p'] + d[fi+'diff']*(1-a/d['dist'])\n    lng, lat, alt = ECEF_to_WGS84(d['xnew'].values,d['ynew'].values,d['znew'].values)\n    \n    lng[np.isnan(lng)] = d.loc[np.isnan(lng),'lngDeg']\n    lat[np.isnan(lat)] = d.loc[np.isnan(lat),'latDeg']\n    d['latDeg'] = lat\n    d['lngDeg'] = lng\n    \n    d.sort_values(['phone',msge],inplace = True)\n    ffname = 'shifted_' + fname\n    d[sub_columns].to_csv(ffname, index = False)\n    return ffname \ndef objective(trial):\n    a = trial.suggest_uniform('a', -1, 1)\n    score, scores = compute_dist(position_shift('btrain.csv', a),'gt.csv')\n    return score\n\nstudy = optuna.create_study()\nstudy.optimize(objective, n_trials=10)\n\nprint(study.best_params)\n\n! cp /kaggle/input/baseline-post-processing-by-outlier-correction/submission.csv sub.csv\nposition_shift('sub.csv', a = study.best_params['a'])\n! cp shifted_sub.csv submission.csv\n\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2022-01-27T17:15:25.319901Z","iopub.execute_input":"2022-01-27T17:15:25.320301Z","iopub.status.idle":"2022-01-27T17:17:04.142811Z","shell.execute_reply.started":"2022-01-27T17:15:25.320268Z","shell.execute_reply":"2022-01-27T17:17:04.141677Z"},"trusted":true},"execution_count":null,"outputs":[]}]}