{"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":"# Google Smartphone Decimeter Challenge by Adaptive Gauss\n\nSolution is based on referring to various notebooks posted at Kaggle, along with my changes. Seeing this as a learning experience.\n\n## Notes\n\n If you’re outside, with open sky, the GPS accuracy from your phone is about five meters, and that’s been constant for a while. With raw GNSS measurements from the phones, this can now improve dramatically.\n \n### The GNSS problem description\nhttps://github.com/commaai/laika\n\nGNSS satellites orbit the earth broadcasting signals that allow the receiver to determine the distance to each satellite. These satellites have known orbits and so their positions are known. This makes determining the receiver's position a basic 3-dimensional trilateration problem. In practice observed distances to each satellite will be measured with some offset that is caused by the receiver's clock error. This offset also needs to be determined, making it a 4-dimensional trilateration problem. \n\n<img src= \"https://camo.githubusercontent.com/0d85f5131c63442f8e7b46de7dab8040a7d693effd5e611ebed25be0b7600a32/68747470733a2f2f75706c6f61642e77696b696d656469612e6f72672f77696b6970656469612f636f6d6d6f6e732f7468756d622f632f63332f33737068657265732e7376672f36323270782d33737068657265732e7376672e706e67\"  alt =\"GNSS\" style=\"width:400px;height:400px;\">\n\nSince this problem is generally overdetermined (more than 4 satellites to solve the 4d problem) there is a variety of methods to compute a position estimate from the measurements. One can use a basic weighted least squares solver for experimental purposes. This is far from optimal due to the dynamic nature of the system, this makes a Bayesian estimator like a Kalman filter the preferred estimator.\n\nHowever, the above description is over-simplified. Getting accurate distance estimates to satellites and the satellite's position from the receiver observations is not trivial. This is what we call processing of the GNSS observables and it is this procedure laika is designed to make easy.\n \n ### How positioning works?\n - Send a burst of these transactions and, as a consequence, the system can calculate ranging statistics, such as the mean and the variance.\n \n <img src= \"https://www.gpsworld.com/wp-content/uploads/2018/07/Android-Figure-3.jpg\" alt =\"Wifi distance\">\n \n Wi-Fi RTT principles, basic concept. Image by Frank van Diggelen, Roy Want and Wei Wang\n \n - Take these ranges of separate access points; if those ranges were accurate, they would define four circles that would intersect at a single point. In practice, because of error in each range, a maximum likelihood position is calculated using a least squares multilateration algorithm.\n \n - Further refine this position by repeating the process, particularly as the phone moves, and then calculate trajectory using filtering techniques, such as Kalman filtering, to optimize the estimate.\n \n  <img src= \"https://www.gpsworld.com/wp-content/uploads/2018/07/Android-Figure-4.jpg\" alt =\"Workflow\">\n  \n  Wi-Fi Workflow. Image by: Frank van Diggelen, Roy Want and Wei Wang.\n\n \n \n### Steps for the solution from Sohier Dane\n- Smoothing out the baseline estimates\n- Integrating readings from other phone instruments, like the accelerometer.\n- Satellite triangulation using the *derived.csv files.\n- Building triangulations directly from the raw gnss logs. \n- Incorporating external data for controls like satellite readings from base stations in the area.\n\n### References: \n- https://www.kaggle.com/t88take/gsdc-phones-mean-prediction\n- Discussions topics from Sohier Dane\n- Google I/O https://www.gpsworld.com/how-to-achieve-1-meter-accuracy-in-android/\n- GPS Survey Workshop video https://www.youtube.com/watch?v=vOJ3u7Zd_i0\n- Hardware: Centimeter Positioning with a Smartphone-Quality GNSS Antenna https://www.youtube.com/watch?v=rCOvklUB5vQ\n- https://www.kaggle.com/bpetrb/kalman-filter-hyperparameter-search-with-bo","metadata":{}},{"cell_type":"code","source":"import numpy as np\nimport pandas as pd\nfrom scipy.ndimage import gaussian_filter1d\nfrom scipy.interpolate import interp1d","metadata":{"execution":{"iopub.status.busy":"2021-07-24T08:16:05.116676Z","iopub.execute_input":"2021-07-24T08:16:05.117105Z","iopub.status.idle":"2021-07-24T08:16:05.547437Z","shell.execute_reply.started":"2021-07-24T08:16:05.117019Z","shell.execute_reply":"2021-07-24T08:16:05.546472Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def apply_gauss_smoothing(df, params):\n    SZ_1 = params['sz_1']\n    SZ_2 = params['sz_2']\n    SZ_CRIT = params['sz_crit']    \n    \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                \n        lat_g1 = gaussian_filter1d(data[:, 0], np.sqrt(SZ_1))\n        lon_g1 = gaussian_filter1d(data[:, 1], np.sqrt(SZ_1))\n        lat_g2 = gaussian_filter1d(data[:, 0], np.sqrt(SZ_2))\n        lon_g2 = gaussian_filter1d(data[:, 1], np.sqrt(SZ_2))\n\n        lat_dif = data[1:,0] - data[:-1,0]\n        lon_dif = data[1:,1] - data[:-1,1]\n\n        lat_crit = np.append(np.abs(gaussian_filter1d(lat_dif, np.sqrt(SZ_CRIT)) / (1e-9 + gaussian_filter1d(np.abs(lat_dif), np.sqrt(SZ_CRIT)))),[0])\n        lon_crit = np.append(np.abs(gaussian_filter1d(lon_dif, np.sqrt(SZ_CRIT)) / (1e-9 + gaussian_filter1d(np.abs(lon_dif), np.sqrt(SZ_CRIT)))),[0])           \n            \n        df.loc[cond, 'latDeg'] = lat_g1 * lat_crit + lat_g2 * (1.0 - lat_crit)\n        df.loc[cond, 'lngDeg'] = lon_g1 * lon_crit + lon_g2 * (1.0 - lon_crit)    \n                       \n    return df","metadata":{"execution":{"iopub.status.busy":"2021-07-24T08:16:05.548735Z","iopub.execute_input":"2021-07-24T08:16:05.548992Z","iopub.status.idle":"2021-07-24T08:16:05.560287Z","shell.execute_reply.started":"2021-07-24T08:16:05.548968Z","shell.execute_reply":"2021-07-24T08:16:05.559425Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def mean_with_other_phones(df):\n    collections_list = df[['collectionName']].drop_duplicates().to_numpy()\n\n    for collection in collections_list:\n        phone_list = df[df['collectionName'].to_list() == collection][['phoneName']].drop_duplicates().to_numpy()\n\n        phone_data = {}\n        corrections = {}\n        for phone in phone_list:\n            cond = np.logical_and(df['collectionName'] == collection[0], df['phoneName'] == phone[0]).to_list()\n            phone_data[phone[0]] = df[cond][['millisSinceGpsEpoch', 'latDeg', 'lngDeg']].to_numpy()\n\n        for current in phone_data:\n            correction = np.ones(phone_data[current].shape, dtype=np.float)\n            correction[:,1:] = phone_data[current][:,1:]\n            \n            # Telephones data don't complitely match by time, so - interpolate.\n            for other in phone_data:\n                if other == current:\n                    continue\n\n                loc = interp1d(phone_data[other][:,0], \n                               phone_data[other][:,1:], \n                               axis=0, \n                               kind='linear', \n                               copy=False, \n                               bounds_error=None, \n                               fill_value='extrapolate', \n                               assume_sorted=True)\n                \n                start_idx = 0\n                stop_idx = 0\n                for idx, val in enumerate(phone_data[current][:,0]):\n                    if val < phone_data[other][0,0]:\n                        start_idx = idx\n                    if val < phone_data[other][-1,0]:\n                        stop_idx = idx\n\n                if stop_idx - start_idx > 0:\n                    correction[start_idx:stop_idx,0] += 1\n                    correction[start_idx:stop_idx,1:] += loc(phone_data[current][start_idx:stop_idx,0])                    \n\n            correction[:,1] /= correction[:,0]\n            correction[:,2] /= correction[:,0]\n            \n            corrections[current] = correction.copy()\n        \n        for phone in phone_list:\n            cond = np.logical_and(df['collectionName'] == collection[0], df['phoneName'] == phone[0]).to_list()\n            \n            df.loc[cond, ['latDeg', 'lngDeg']] = corrections[phone[0]][:,1:]            \n            \n    return df","metadata":{"execution":{"iopub.status.busy":"2021-07-24T08:16:05.561893Z","iopub.execute_input":"2021-07-24T08:16:05.562153Z","iopub.status.idle":"2021-07-24T08:16:05.574568Z","shell.execute_reply.started":"2021-07-24T08:16:05.562129Z","shell.execute_reply":"2021-07-24T08:16:05.573652Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"test_base = pd.read_csv('../input/google-smartphone-decimeter-challenge/baseline_locations_test.csv')\nsub = pd.read_csv('../input/google-smartphone-decimeter-challenge/sample_submission.csv')\n\nsmoothed_baseline = apply_gauss_smoothing(test_base, {'sz_1' : 0.85, 'sz_2' : 5.65, 'sz_crit' : 0.9})\nsmoothed_baseline = mean_with_other_phones(smoothed_baseline)\n\nsub = sub.assign( latDeg=smoothed_baseline.latDeg, lngDeg=smoothed_baseline.lngDeg )\nsub.to_csv('submission.csv', index=False)","metadata":{"execution":{"iopub.status.busy":"2021-07-24T08:16:05.575751Z","iopub.execute_input":"2021-07-24T08:16:05.576014Z","iopub.status.idle":"2021-07-24T08:16:09.642429Z","shell.execute_reply.started":"2021-07-24T08:16:05.575989Z","shell.execute_reply":"2021-07-24T08:16:09.641530Z"},"trusted":true},"execution_count":null,"outputs":[]}]}