{"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)\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\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\nfrom pathlib import Path\ndatapath = Path(\"../input/google-smartphone-decimeter-challenge/\")\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-05-26T00:31:42.691966Z","iopub.execute_input":"2021-05-26T00:31:42.692558Z","iopub.status.idle":"2021-05-26T00:31:42.69792Z","shell.execute_reply.started":"2021-05-26T00:31:42.69247Z","shell.execute_reply":"2021-05-26T00:31:42.69699Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Background\n\nEvery satellite in space continuously broadcasts its current position coordinate information to the world through satellite signals. Any GNSS receiver can receive and read this information through an antenna.\n\nWhen the satellite sends position information, it also attaches the timestamp when the data was sent. After the GNSS receiver receives the data, it subtracts the time on the timestamp from the current time, which is the time taken for the data transmission. And the approximate distance between the satellite and the GPS receiver can be calculated by data transmission time multiplied by the transmission speed, which is the `rawPrM` in `_derived` table. \n\nThe most difficult problem in GPS positioning is error. There are many reasons for positioning errors, such as the ionosphere, the receiving device, and the blocking and multipath effects. Inorder to get a more accurate the approximate distance between the satellite and the GPS receiver with as less error as possible, we can use a corrected approximate distance as [home page](https://www.kaggle.com/c/google-smartphone-decimeter-challenge/data) suggested, which can be computed as: `correctedPrM = rawPrM + satClkBiasM - isrbM - ionoDelayM - tropoDelayM`\n\n<img src='https://www.researchgate.net/publication/324989880/figure/fig1/AS:623672356253699@1525706485241/Global-Navigation-Satellite-System-GNSS-spoofing-attack-illustration.png' width=400/>\n\nNormally these values are need to be calculated carefully from GNSSLog file by lots of domain knowledge. Thanks to the dataset provider, we can directly use the intermediate calculated value from derived files to built the initial version of the model.","metadata":{}},{"cell_type":"markdown","source":"### Basic Utils Functions","metadata":{"_kg_hide-input":true}},{"cell_type":"code","source":"def ecef2lla(x, y, z):\n  # x, y and z are scalars or vectors in meters\n  x = np.array([x]).reshape(np.array([x]).shape[-1], 1)\n  y = np.array([y]).reshape(np.array([y]).shape[-1], 1)\n  z = np.array([z]).reshape(np.array([z]).shape[-1], 1)\n\n  a=6378137\n  a_sq=a**2\n  e = 8.181919084261345e-2\n  e_sq = 6.69437999014e-3\n\n  f = 1/298.257223563\n  b = a*(1-f)\n\n  # calculations:\n  r = np.sqrt(x**2 + y**2)\n  ep_sq  = (a**2-b**2)/b**2\n  ee = (a**2-b**2)\n  f = (54*b**2)*(z**2)\n  g = r**2 + (1 - e_sq)*(z**2) - e_sq*ee*2\n  c = (e_sq**2)*f*r**2/(g**3)\n  s = (1 + c + np.sqrt(c**2 + 2*c))**(1/3.)\n  p = f/(3.*(g**2)*(s + (1./s) + 1)**2)\n  q = np.sqrt(1 + 2*p*e_sq**2)\n  r_0 = -(p*e_sq*r)/(1+q) + np.sqrt(0.5*(a**2)*(1+(1./q)) - p*(z**2)*(1-e_sq)/(q*(1+q)) - 0.5*p*(r**2))\n  u = np.sqrt((r - e_sq*r_0)**2 + z**2)\n  v = np.sqrt((r - e_sq*r_0)**2 + (1 - e_sq)*z**2)\n  z_0 = (b**2)*z/(a*v)\n  h = u*(1 - b**2/(a*v))\n  phi = np.arctan((z + ep_sq*z_0)/r)\n  lambd = np.arctan2(y, x)\n\n  return phi*180/np.pi, lambd*180/np.pi, h","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:42.699639Z","iopub.execute_input":"2021-05-26T00:31:42.700077Z","iopub.status.idle":"2021-05-26T00:31:42.71987Z","shell.execute_reply.started":"2021-05-26T00:31:42.700044Z","shell.execute_reply":"2021-05-26T00:31:42.718686Z"},"_kg_hide-input":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Load Data","metadata":{}},{"cell_type":"code","source":"df_baseline = pd.read_csv(datapath/\"baseline_locations_train.csv\")\ndf_sample_trail_gt = pd.read_csv(datapath/\"train/2020-05-14-US-MTV-1/Pixel4/ground_truth.csv\")\ndf_sample_trail = pd.read_csv(datapath/\"train/2020-05-14-US-MTV-1/Pixel4/Pixel4_derived.csv\")\ndf_sample_trail[\"correctedPrM\"] = df_sample_trail[\"rawPrM\"] + df_sample_trail[\"satClkBiasM\"] - df_sample_trail[\"isrbM\"] - df_sample_trail[\"ionoDelayM\"] - df_sample_trail[\"tropoDelayM\"] ","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:42.721304Z","iopub.execute_input":"2021-05-26T00:31:42.721614Z","iopub.status.idle":"2021-05-26T00:31:43.607485Z","shell.execute_reply.started":"2021-05-26T00:31:42.721585Z","shell.execute_reply":"2021-05-26T00:31:43.606285Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Simple Least Square Solution\n\nAn intuitive solution to solve this problem is utilizing the basic geometry information to infer the GNSS receiver location. Since we have multiple satellite, we can optimize the location to achieve the minimum projection error among all satellite based by least square.\n\nTo be more specific,\n\n$$r_u^k = ||x^k_{sat} - x_u||$$\n\nwhere,\n\n$r_u^k$: the actually distance between user device (GNSS Receiver) to Satellite;\n\n$x^k_{sat}$: the actually position of Satellite;\n\n$x_u$: the actually position of user device (GNSS Receiver).\n\n$$ Pr_u^k = r_u^k + b_u + \\sigma^k $$\n\nwhere,\n\n$Pr_u^k$: the corrected pseudorange (i.e. a closer approximation to the geometric range from the phone to the satellite) can be computed as: $correctedPrM = rawPrM + satClkBiasM - isrbM - ionoDelayM - tropoDelayM$;\n\n$b_u$: the user clock bias in meters equivalent, which $b_u = clock\\_bias_u * LIGHTSPEED$;\n\n$\\sigma^k$: the actually position of user device (GNSS Receiver).\n\nOur aim is to minimize the measured and the estimated pseudorange for each satellite:\n\n$$\\delta Pr_u^k = Pr_u^k - (r_u^k + b_u) = Pr_u^k - (||x^k_{sat} - x_u|| + b_u)$$\n\nwhich is our target function in least square solver.\n\n<!-- Which can be deducted to based equation (5) in [here](https://www.telesens.co/2017/07/17/calculating-position-from-raw-gps-data/):\n\n$$\\delta Pr_u^k = - \\frac{x^k_{sat} - x_u}{||x^k_{sat} - x_u||} \\delta x_u^k + \\delta b_u + \\sigma^k = - \\hat{x}_u^k \\delta x_u^k + \\delta b_u + \\sigma^k$$ -->\n\nLet's try it with 1 single epoch.","metadata":{}},{"cell_type":"code","source":"epoch_time = 1273529464442\ndf_sample_epoch = df_sample_trail[df_sample_trail.millisSinceGpsEpoch == epoch_time]\ndf_sample_epoch_gt = df_sample_trail_gt[df_sample_trail_gt.millisSinceGpsEpoch == epoch_time]\ndf_sample_epoch_baseline = df_baseline[(df_baseline.collectionName == \"2020-05-14-US-MTV-1\") & (df_baseline.phoneName == \"Pixel4\") & (df_baseline.millisSinceGpsEpoch == epoch_time)]\n\nsat_pos = df_sample_epoch[[\"xSatPosM\",\"ySatPosM\",\"zSatPosM\"]].to_numpy()\npseudoranges = np.squeeze(df_sample_epoch[[\"correctedPrM\"]].to_numpy())","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:43.608963Z","iopub.execute_input":"2021-05-26T00:31:43.609309Z","iopub.status.idle":"2021-05-26T00:31:43.677006Z","shell.execute_reply.started":"2021-05-26T00:31:43.609269Z","shell.execute_reply":"2021-05-26T00:31:43.676053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### self implement simple least square\n\nBuild the simplest version from [code](https://www.kaggle.com/realbigmitchell/simple-least-squares-solution) and [article](https://www.telesens.co/2017/07/17/calculating-position-from-raw-gps-data/).","metadata":{}},{"cell_type":"code","source":"def least_squares(sat_pos, pseudoranges, weights=1, x_hat=np.array([0, 0, 0, 0])):\n  \"\"\"\n  Args:\n    sat_pos: The satellite position (meters) in an ECEF coordinate frame\n    pseudoranges: The corrected pseudorange (i.e. a closer approximation to the geometric range from the phone to the satellite)\n    x_hat: the phone's initial/previous estimated position (x, y, z, b) and \n           b represent the user clock bias in units of distance = clock bias (t) * light speed (c)\n  \n  Returns:\n    x_hat: the phone's estimated position\n    norm_dp:\n  \"\"\"\n  dx = np.Inf*np.ones(3)\n  G = np.ones((pseudoranges.size, 4))\n  iterations = 0\n\n  if isinstance(weights, np.ndarray):\n    weights = np.diag(weights)\n  else:\n    weights = weights*np.eye(pseudoranges.size)\n\n  while np.linalg.norm(dx) > 1e-3:\n    norms = np.linalg.norm(sat_pos - x_hat[:3], axis=1)\n    dp = pseudoranges - norms - x_hat[3]\n    G[:, 0:3] = -(sat_pos - x_hat[:3])/norms[:, None]\n    # G_T = np.transpose(G)\n    # dx = np.linalg.inv(G_T@G) @ G_T @ dp\n    dx = np.linalg.pinv(weights@G) @ weights @ dp\n    x_hat = x_hat + dx\n    iterations += 1\n  return x_hat, np.linalg.norm(dp)","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:43.681318Z","iopub.execute_input":"2021-05-26T00:31:43.681682Z","iopub.status.idle":"2021-05-26T00:31:43.691381Z","shell.execute_reply.started":"2021-05-26T00:31:43.681651Z","shell.execute_reply":"2021-05-26T00:31:43.690448Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x, dp = least_squares(sat_pos, pseudoranges)\n\nprint(\"Ground truth:\", df_sample_epoch_gt[[\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"]].to_numpy())\nprint(\"Weighted Least Square Estimation (Baseline):\", df_sample_epoch_baseline[[\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"]].to_numpy())\nprint(\"Simple Least Square Estimation:\", *ecef2lla(*x[:3]))\n# print(b/LIGHTSPEED)\n# print(dp)","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:43.693312Z","iopub.execute_input":"2021-05-26T00:31:43.693699Z","iopub.status.idle":"2021-05-26T00:31:43.732853Z","shell.execute_reply.started":"2021-05-26T00:31:43.693602Z","shell.execute_reply":"2021-05-26T00:31:43.731188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### scipy optimize version\n\nBuild the simplest version from [this code](https://github.com/commaai/laika/blob/master/laika/raw_gnss.py#L254-L350)","metadata":{}},{"cell_type":"code","source":"import scipy.optimize as opt\n\ndef calc_pos_fix(sat_pos, pr, weights=1, x0=[0, 0, 0, 0]):\n  '''\n  Calculates gps fix with WLS optimizer\n  returns:\n  0 -> list with positions\n  1 -> pseudorange errs\n  '''\n  n = len(pr)\n  if n < 3:\n      return x0, []\n\n  Fx_pos = pr_residual(sat_pos, pr, weights=weights)\n  opt_pos = opt.least_squares(Fx_pos, x0).x\n  return opt_pos, Fx_pos(opt_pos, weights=1)\n\n\ndef pr_residual(sat_pos, pr, weights=1):\n  # solve for pos\n  def Fx_pos(x_hat, weights=weights):\n    rows = weights * (np.linalg.norm(sat_pos - x_hat[:3], axis=1) + x_hat[3] - pr)\n    return rows\n  return Fx_pos","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:43.738987Z","iopub.execute_input":"2021-05-26T00:31:43.739668Z","iopub.status.idle":"2021-05-26T00:31:43.933669Z","shell.execute_reply.started":"2021-05-26T00:31:43.739618Z","shell.execute_reply":"2021-05-26T00:31:43.932403Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"x, dp = calc_pos_fix(sat_pos, pseudoranges)\n\nprint(\"Ground truth:\", df_sample_epoch_gt[[\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"]].to_numpy())\nprint(\"Weighted Least Square Estimation (Baseline):\", df_sample_epoch_baseline[[\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"]].to_numpy())\nprint(\"Simple Least Square Estimation:\", *ecef2lla(*x[:3]))\n# print(x[3]/LIGHTSPEED)\n# print(dp)","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:43.935322Z","iopub.execute_input":"2021-05-26T00:31:43.935981Z","iopub.status.idle":"2021-05-26T00:31:43.984095Z","shell.execute_reply.started":"2021-05-26T00:31:43.935932Z","shell.execute_reply":"2021-05-26T00:31:43.982922Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Weighted Least Square Solution\n\nThe signal from each satellite may not that reliable. In this case, we can apply the reliability of the signal as weights for each signal to improve the least square result. Here, we use the `rawPrUncM` as the weights for signal, which represents the pseudorange uncertainty/deviation of each signal. The relative code can be found at [the google official repo](WLS shown in official repo](https://github.com/google/gps-measurement-tools/blob/master/opensource/WlsPvt.m#L62-L111)","metadata":{}},{"cell_type":"code","source":"pseudoranges_sigma = np.squeeze(df_sample_epoch[[\"rawPrUncM\"]].to_numpy())\n\nx, dp = least_squares(sat_pos, pseudoranges, 1/pseudoranges_sigma)\n\nprint(\"Ground truth:\", df_sample_epoch_gt[[\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"]].to_numpy())\nprint(\"Weighted Least Square Estimation (Baseline):\", df_sample_epoch_baseline[[\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"]].to_numpy())\nprint(\"Weighted Least Square Estimation:\", *ecef2lla(*x[:3]))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pseudoranges_sigma = np.squeeze(df_sample_epoch[[\"rawPrUncM\"]].to_numpy())\n\nx, dp = calc_pos_fix(sat_pos, pseudoranges, 1/pseudoranges_sigma)\n\nprint(\"Ground truth:\", df_sample_epoch_gt[[\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"]].to_numpy())\nprint(\"Weighted Least Square Estimation (Baseline):\", df_sample_epoch_baseline[[\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"]].to_numpy())\nprint(\"Weighted Least Square Estimation:\", *ecef2lla(*x[:3]))","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Evaluation","metadata":{}},{"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-05-26T00:31:43.985635Z","iopub.execute_input":"2021-05-26T00:31:43.98605Z","iopub.status.idle":"2021-05-26T00:31:43.993998Z","shell.execute_reply.started":"2021-05-26T00:31:43.985997Z","shell.execute_reply":"2021-05-26T00:31:43.992859Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Evaluate in one epoch","metadata":{}},{"cell_type":"code","source":"deg_gt = df_sample_epoch_gt[[\"latDeg\",\"lngDeg\"]].to_numpy().transpose()\ndeg_baseline = df_sample_epoch_baseline[[\"latDeg\",\"lngDeg\"]].to_numpy().transpose()\nprint(\"Weighted Least Square (baseline) haversine distance (M):\", calc_haversine(*deg_gt, *deg_baseline))\nprint(\"Weighted Least Square haversine distance (M):\", calc_haversine(*deg_gt, *ecef2lla(*x[:3])[:2]))","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:43.995301Z","iopub.execute_input":"2021-05-26T00:31:43.995645Z","iopub.status.idle":"2021-05-26T00:31:44.015323Z","shell.execute_reply.started":"2021-05-26T00:31:43.995617Z","shell.execute_reply":"2021-05-26T00:31:44.013982Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Evaluate in whole training","metadata":{}},{"cell_type":"code","source":"def simple_pipeline(df_trails):\n  \"\"\" simple pipeline to estimate the GNSS receiver location by least square\n  Args:\n    df_trails: the df read from derived file\n  \n  Returns:\n    result df with estimated degrees and heights\n  \"\"\"\n  df_trails[\"correctedPrM\"] = df_trails[\"rawPrM\"] + df_trails[\"satClkBiasM\"] - df_trails[\"isrbM\"] - df_trails[\"ionoDelayM\"] - df_trails[\"tropoDelayM\"]\n  \n  results = []\n  x = [0, 0, 0, 0]\n  df_epochs = df_trails.groupby([\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"])\n  for indices, df_epoch in tqdm(df_epochs, desc=\"Estimate location by LS for epoch\"):\n    sat_pos = df_epoch[[\"xSatPosM\",\"ySatPosM\",\"zSatPosM\"]].to_numpy()\n    pseudoranges = np.squeeze(df_epoch[[\"correctedPrM\"]].to_numpy())\n    pseudoranges_sigma = np.squeeze(df_epoch[[\"rawPrUncM\"]].to_numpy())\n    x, _ = calc_pos_fix(sat_pos, pseudoranges, 1/pseudoranges_sigma, x)\n    # x, _ = calc_pos_fix(sat_pos, pseudoranges, 1, x)\n    values = np.squeeze(ecef2lla(*x[:3]))\n    results.append([*indices, *values])\n  return pd.DataFrame(results,columns=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\", \"latDeg\", \"lngDeg\", \"heightAboveWgs84EllipsoidM\"])","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:44.016927Z","iopub.execute_input":"2021-05-26T00:31:44.017416Z","iopub.status.idle":"2021-05-26T00:31:44.028309Z","shell.execute_reply.started":"2021-05-26T00:31:44.017384Z","shell.execute_reply":"2021-05-26T00:31:44.027132Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from pathlib import Path\n# from tqdm import tqdm\nfrom tqdm.notebook import tqdm\n\n# datapath = Path(\"./data\")\nground_truths = (datapath / \"train\").rglob(\"ground_truth.csv\")\ndrived_files = (datapath / \"train\").rglob(\"*_derived.csv\")\n\ndf_sample_trails_baseline = pd.read_csv(datapath / 'baseline_locations_train.csv')\ndf_sample_trails_gt = pd.concat([pd.read_csv(filepath) for filepath in tqdm(ground_truths, total=73, desc=\"Reading ground truth data\")], ignore_index=True)\ndf_sample_trails = pd.concat([pd.read_csv(filepath) for filepath in tqdm(drived_files, total=73, desc=\"Reading drived data\")], ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:31:44.029673Z","iopub.execute_input":"2021-05-26T00:31:44.029994Z","iopub.status.idle":"2021-05-26T00:32:08.388653Z","shell.execute_reply.started":"2021-05-26T00:31:44.029965Z","shell.execute_reply":"2021-05-26T00:32:08.387039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trails_gt[\"receivedSvTimeInGpsNanos\"] = df_sample_trails_gt.millisSinceGpsEpoch*int(1e6)\ndf_sample_trails_raw = df_sample_trails.drop(\"millisSinceGpsEpoch\", axis=1)\n\ndf_merge = pd.merge_asof(df_sample_trails_raw.sort_values('receivedSvTimeInGpsNanos'), df_sample_trails_gt.sort_values('receivedSvTimeInGpsNanos'), \n                                           on=\"receivedSvTimeInGpsNanos\", by=[\"collectionName\", \"phoneName\"], direction='nearest',tolerance=int(1e9))\ndf_merge = df_merge.sort_values(by=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"], ignore_index=True)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trails_estimate = simple_pipeline(df_merge)","metadata":{"execution":{"iopub.status.busy":"2021-05-26T00:32:08.390143Z","iopub.execute_input":"2021-05-26T00:32:08.390514Z","iopub.status.idle":"2021-05-26T01:16:44.241904Z","shell.execute_reply.started":"2021-05-26T00:32:08.390479Z","shell.execute_reply":"2021-05-26T01:16:44.240659Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trails_merged_baseline = pd.merge_asof(df_sample_trails_gt.sort_values('millisSinceGpsEpoch'),\n                                                 df_sample_trails_baseline.sort_values('millisSinceGpsEpoch'), \n                                                 on=\"millisSinceGpsEpoch\", by=[\"collectionName\", \"phoneName\"], \n                                                 direction='nearest',tolerance=100000, suffixes=('_truth', '_pred'))\ndf_sample_trails_merged_baseline = df_sample_trails_merged_baseline.sort_values(by=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"], ignore_index=True)\n\ndf_sample_trails_merged_SL = pd.merge_asof(df_sample_trails_gt.sort_values('millisSinceGpsEpoch'), \n                                           df_sample_trails_estimate.sort_values('millisSinceGpsEpoch'), \n                                           on=\"millisSinceGpsEpoch\", by=[\"collectionName\", \"phoneName\"], \n                                           direction='nearest',tolerance=100000, suffixes=('_truth', '_pred'))\ndf_sample_trails_merged_SL = df_sample_trails_merged_SL.sort_values(by=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"], ignore_index=True)\n\ncompared_cols = [\"latDeg_truth\",\"lngDeg_truth\",\"latDeg_pred\",\"lngDeg_pred\"]\nprint(\"Weighted Least Square (baseline) haversine distance (M):\", calc_haversine(*df_sample_trails_merged_baseline[compared_cols].to_numpy().transpose()).mean())\nprint(\"Weighted Least Square haversine distance (M):\", calc_haversine(*df_sample_trails_merged_SL[compared_cols].to_numpy().transpose()).mean())","metadata":{"execution":{"iopub.status.busy":"2021-05-26T01:16:44.243426Z","iopub.execute_input":"2021-05-26T01:16:44.243723Z","iopub.status.idle":"2021-05-26T01:16:44.475114Z","shell.execute_reply.started":"2021-05-26T01:16:44.243694Z","shell.execute_reply":"2021-05-26T01:16:44.474069Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submission\n\nIt should be noticed that the epoch time in the drived file and submission/baseline file is not same, in this case, we need to merge it to submission format with epoch time tolerant and sort it into suitable output format (based on quick experiment, 100s tolerant can be merged without empty value).","metadata":{}},{"cell_type":"code","source":"from pathlib import Path\n# from tqdm import tqdm\nfrom tqdm.notebook import tqdm\n\n# datapath = Path(\"./data\")\ndrived_files = (datapath / \"test\").rglob(\"*_derived.csv\")\n\ndf_sample_trails_baseline = pd.read_csv(datapath / 'baseline_locations_test.csv')\ndf_sample_trails = pd.concat([pd.read_csv(filepath) for filepath in tqdm(drived_files, total=48, desc=\"Reading drived data\")], ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2021-05-26T01:16:44.47657Z","iopub.execute_input":"2021-05-26T01:16:44.476893Z","iopub.status.idle":"2021-05-26T01:17:01.581623Z","shell.execute_reply.started":"2021-05-26T01:16:44.476861Z","shell.execute_reply":"2021-05-26T01:17:01.580635Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trails_baseline[\"receivedSvTimeInGpsNanos\"] = df_sample_trails_baseline.millisSinceGpsEpoch*int(1e6)\ndf_sample_trails_raw = df_sample_trails.drop(\"millisSinceGpsEpoch\", axis=1)\n\ndf_merge = pd.merge_asof(df_sample_trails_raw.sort_values('receivedSvTimeInGpsNanos'), df_sample_trails_baseline.sort_values('receivedSvTimeInGpsNanos'), \n                                           on=\"receivedSvTimeInGpsNanos\", by=[\"collectionName\", \"phoneName\"], direction='nearest',tolerance=int(1e9))\ndf_merge = df_merge.sort_values(by=[\"collectionName\", \"phoneName\", \"millisSinceGpsEpoch\"], ignore_index=True)","metadata":{"execution":{"iopub.status.busy":"2021-05-26T01:17:01.58268Z","iopub.execute_input":"2021-05-26T01:17:01.583097Z","iopub.status.idle":"2021-05-26T01:47:45.432735Z","shell.execute_reply.started":"2021-05-26T01:17:01.583067Z","shell.execute_reply":"2021-05-26T01:47:45.431604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trails_estimate = simple_pipeline(df_merge)","metadata":{},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_sample_trails_baseline = df_sample_trails_baseline.drop([\"latDeg\",\"lngDeg\",\"heightAboveWgs84EllipsoidM\"], axis=1)\ndf_sample_trails_merged = pd.merge_asof(df_sample_trails_baseline.sort_values('millisSinceGpsEpoch'), \n                                        df_sample_trails_estimate.sort_values('millisSinceGpsEpoch'), \n                                        on=\"millisSinceGpsEpoch\", by=[\"collectionName\", \"phoneName\"], direction='nearest', tolerance=100000)\ndf_sample_trails_merged = df_sample_trails_merged.sort_values(by=[\"phone\", \"millisSinceGpsEpoch\"], ignore_index=True)\n\ndf_submission = df_sample_trails_merged[[\"phone\", \"millisSinceGpsEpoch\", \"latDeg\", \"lngDeg\"]].copy()\ndf_submission.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-05-26T01:47:45.434341Z","iopub.execute_input":"2021-05-26T01:47:45.434635Z","iopub.status.idle":"2021-05-26T01:47:46.037209Z","shell.execute_reply.started":"2021-05-26T01:47:45.434608Z","shell.execute_reply":"2021-05-26T01:47:46.035729Z"},"trusted":true},"execution_count":null,"outputs":[]}]}