{"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 Keras\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/jeongyoonlee/google-smartphone-decimeter-eda-keras-tpu\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","metadata":{}},{"cell_type":"markdown","source":"Install additional libraries\n- `simdkalman` is for Kalman Filters","metadata":{}},{"cell_type":"code","source":"!pip install simdkalman","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-07-16T03:30:38.442647Z","iopub.execute_input":"2021-07-16T03:30:38.443239Z","iopub.status.idle":"2021-07-16T03:30:49.161892Z","shell.execute_reply.started":"2021-07-16T03:30:38.443196Z","shell.execute_reply":"2021-07-16T03:30:49.160520Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"%matplotlib inline\nfrom matplotlib import pyplot as plt\nimport numpy as np # linear algebra\nfrom pathlib import Path\nimport pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\nfrom scipy import sparse\nfrom sklearn.metrics import mean_squared_error\nfrom sklearn.model_selection import KFold\nfrom sklearn.preprocessing import StandardScaler\nimport seaborn as sns\nimport simdkalman\nimport tensorflow as tf\nfrom tensorflow import keras\nfrom tqdm.notebook import tqdm\nfrom warnings import simplefilter\n\nsimplefilter('ignore')\nplt.style.use('fivethirtyeight')\npd.set_option('max_columns', 100)\npd.set_option('max_rows', 100)","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","execution":{"iopub.status.busy":"2021-07-16T03:30:49.164576Z","iopub.execute_input":"2021-07-16T03:30:49.165132Z","iopub.status.idle":"2021-07-16T03:30:49.195294Z","shell.execute_reply.started":"2021-07-16T03:30:49.165068Z","shell.execute_reply":"2021-07-16T03:30:49.188770Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model config","metadata":{}},{"cell_type":"code","source":"model_name = 'nn_v2'\n\ndata_dir = Path('../input/google-smartphone-decimeter-challenge')\ntrain_file = data_dir / 'baseline_locations_train.csv'\ntest_file = data_dir / 'baseline_locations_test.csv'\nsample_file = data_dir / 'sample_submission.csv'\n\nbuild_dir = Path('./build')\nbuild_dir.mkdir(parents=True, exist_ok=True)\npredict_val_file = build_dir / f'{model_name}.val.txt'\npredict_tst_file = build_dir / f'{model_name}.tst.txt'\nsubmission_file = 'submission.csv'\n\ncname_col = 'collectionName'\npname_col = 'phoneName'\nphone_col = 'phone'\nts_col = 'millisSinceGpsEpoch'\ndt_col = 'datetime'\nlat_col = 'latDeg'\nlon_col = 'lngDeg'\n\nlrate = 0.1\nbatch_size = 32 #1024\nepochs = 2000 #100\nn_stop = 10\nn_fold = 5\nseed = 42","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:49.199150Z","iopub.execute_input":"2021-07-16T03:30:49.199556Z","iopub.status.idle":"2021-07-16T03:30:49.219338Z","shell.execute_reply.started":"2021-07-16T03:30:49.199516Z","shell.execute_reply":"2021-07-16T03:30:49.217740Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://www.kaggle.com/sohier/loading-gnss-logs\ndef gnss_log_to_dataframes(path):\n    print('Loading ' + path, flush=True)\n    gnss_section_names = {'Raw','UncalAccel', 'UncalGyro', 'UncalMag', 'Fix', 'Status', 'OrientationDeg'}\n    with open(path) as f_open:\n        datalines = f_open.readlines()\n\n    datas = {k: [] for k in gnss_section_names}\n    gnss_map = {k: [] for k in gnss_section_names}\n    for dataline in datalines:\n        is_header = dataline.startswith('#')\n        dataline = dataline.strip('#').strip().split(',')\n        # skip over notes, version numbers, etc\n        if is_header and dataline[0] in gnss_section_names:\n            gnss_map[dataline[0]] = dataline[1:]\n        elif not is_header:\n            datas[dataline[0]].append(dataline[1:])\n\n    results = dict()\n    for k, v in datas.items():\n        results[k] = pd.DataFrame(v, columns=gnss_map[k])\n    # pandas doesn't properly infer types from these lists by default\n    for k, df in results.items():\n        for col in df.columns:\n            if col == 'CodeType':\n                continue\n            results[k][col] = pd.to_numeric(results[k][col])\n\n    return results","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-16T03:30:49.223938Z","iopub.execute_input":"2021-07-16T03:30:49.224634Z","iopub.status.idle":"2021-07-16T03:30:49.270853Z","shell.execute_reply.started":"2021-07-16T03:30:49.224573Z","shell.execute_reply":"2021-07-16T03:30:49.267316Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://www.kaggle.com/dannellyz/start-here-simple-folium-heatmap-for-geo-data\nimport folium\nfrom folium import plugins\n\n\ndef simple_folium(df:pd.DataFrame, lat_col:str, lon_col:str):\n    \"\"\"\n    Descrption\n    ----------\n        Returns a simple Folium HeatMap with Markers\n    ----------\n    Parameters\n    ----------\n        df : padnas DataFrame, required\n            The DataFrane with the data to map\n        lat_col : str, required\n            The name of the column with latitude\n        lon_col : str, required\n            The name of the column with longitude\n    \"\"\"\n    #Preprocess\n    #Drop rows that do not have lat/lon\n    df = df[df[lat_col].notnull() & df[lon_col].notnull()]\n\n    # Convert lat/lon to (n, 2) nd-array format for heatmap\n    # Then send to list\n    df_locs = list(df[[lat_col, lon_col]].values)\n\n    #Set up folium map\n    fol_map = folium.Map([df[lat_col].median(), df[lon_col].median()])\n\n    # plot heatmap\n    heat_map = plugins.HeatMap(df_locs)\n    fol_map.add_child(heat_map)\n\n    # plot markers\n    markers = plugins.MarkerCluster(locations = df_locs)\n    fol_map.add_child(markers)\n\n    #Add Layer Control\n    folium.LayerControl().add_to(fol_map)\n\n    return fol_map","metadata":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-16T03:30:49.274112Z","iopub.execute_input":"2021-07-16T03:30:49.274686Z","iopub.status.idle":"2021-07-16T03:30:49.295542Z","shell.execute_reply.started":"2021-07-16T03:30:49.274646Z","shell.execute_reply":"2021-07-16T03:30:49.293585Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://www.kaggle.com/jpmiller/baseline-from-host-data\n# 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":{"_kg_hide-input":true,"execution":{"iopub.status.busy":"2021-07-16T03:30:49.297978Z","iopub.execute_input":"2021-07-16T03:30:49.298484Z","iopub.status.idle":"2021-07-16T03:30:49.326023Z","shell.execute_reply.started":"2021-07-16T03:30:49.298433Z","shell.execute_reply":"2021-07-16T03:30:49.324595Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://www.kaggle.com/emaerthin/demonstration-of-the-kalman-filter\nT = 1.0\nstate_transition = np.array([[1, 0, T, 0, 0.5 * T ** 2, 0], [0, 1, 0, T, 0, 0.5 * T ** 2], [0, 0, 1, 0, T, 0],\n                             [0, 0, 0, 1, 0, T], [0, 0, 0, 0, 1, 0], [0, 0, 0, 0, 0, 1]])\nprocess_noise = np.diag([1e-5, 1e-5, 5e-6, 5e-6, 1e-6, 1e-6]) + np.ones((6, 6)) * 1e-9\nobservation_model = np.array([[1, 0, 0, 0, 0, 0], [0, 1, 0, 0, 0, 0]])\nobservation_noise = np.diag([5e-5, 5e-5]) + np.ones((2, 2)) * 1e-9\n\nkf = simdkalman.KalmanFilter(\n        state_transition = state_transition,\n        process_noise = process_noise,\n        observation_model = observation_model,\n        observation_noise = observation_noise)\n\ndef apply_kf_smoothing(df, kf_=kf):\n    unique_paths = df[phone_col].unique()\n    for phone in tqdm(unique_paths):\n        data = df.loc[df[phone_col] == phone][[lat_col, lon_col]].values\n        data = data.reshape(1, len(data), 2)\n        smoothed = kf_.smooth(data)\n        df.loc[df[phone_col] == phone, lat_col] = smoothed.states.mean[0, :, 0]\n        df.loc[df[phone_col] == phone, lon_col] = smoothed.states.mean[0, :, 1]\n    return df","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:49.328495Z","iopub.execute_input":"2021-07-16T03:30:49.329226Z","iopub.status.idle":"2021-07-16T03:30:49.356398Z","shell.execute_reply.started":"2021-07-16T03:30:49.329162Z","shell.execute_reply":"2021-07-16T03:30:49.353931Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn = pd.read_csv(train_file)\nprint(trn.shape)\ntrn.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:49.360216Z","iopub.execute_input":"2021-07-16T03:30:49.360640Z","iopub.status.idle":"2021-07-16T03:30:49.923810Z","shell.execute_reply.started":"2021-07-16T03:30:49.360605Z","shell.execute_reply":"2021-07-16T03:30:49.922028Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tst = pd.read_csv(test_file)\nprint(tst.shape)\ntst.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:49.926688Z","iopub.execute_input":"2021-07-16T03:30:49.927233Z","iopub.status.idle":"2021-07-16T03:30:50.285193Z","shell.execute_reply.started":"2021-07-16T03:30:49.927177Z","shell.execute_reply":"2021-07-16T03:30:50.284130Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub = pd.read_csv(sample_file)\nprint(sub.shape)\nsub.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:50.286828Z","iopub.execute_input":"2021-07-16T03:30:50.287524Z","iopub.status.idle":"2021-07-16T03:30:50.565289Z","shell.execute_reply.started":"2021-07-16T03:30:50.287479Z","shell.execute_reply":"2021-07-16T03:30:50.563957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Exploratory Data Analysis (EDA)\n\n### Aggregator Data\n\nLets see if each phone has given good amount of readings","metadata":{}},{"cell_type":"code","source":"for col in [cname_col, pname_col]:\n    print(f'# of unique {col:>14s} in training: {trn[col].nunique():4d}')\n    print(f'# of unique {col:>14s}     in test: {tst[col].nunique():4d}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:50.566420Z","iopub.execute_input":"2021-07-16T03:30:50.566739Z","iopub.status.idle":"2021-07-16T03:30:50.725918Z","shell.execute_reply.started":"2021-07-16T03:30:50.566705Z","shell.execute_reply":"2021-07-16T03:30:50.724661Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn[pname_col].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:50.730018Z","iopub.execute_input":"2021-07-16T03:30:50.730469Z","iopub.status.idle":"2021-07-16T03:30:50.782204Z","shell.execute_reply.started":"2021-07-16T03:30:50.730429Z","shell.execute_reply":"2021-07-16T03:30:50.781276Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tst[pname_col].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:50.783265Z","iopub.execute_input":"2021-07-16T03:30:50.783601Z","iopub.status.idle":"2021-07-16T03:30:50.826573Z","shell.execute_reply.started":"2021-07-16T03:30:50.783565Z","shell.execute_reply":"2021-07-16T03:30:50.824953Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(f'# of unique phone in training: {trn[phone_col].nunique():4d}')\nprint(f'    # of unique phone in test: {tst[phone_col].nunique():4d}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:50.829608Z","iopub.execute_input":"2021-07-16T03:30:50.830152Z","iopub.status.idle":"2021-07-16T03:30:50.909696Z","shell.execute_reply.started":"2021-07-16T03:30:50.830105Z","shell.execute_reply":"2021-07-16T03:30:50.907604Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn[phone_col].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:50.910972Z","iopub.execute_input":"2021-07-16T03:30:50.911336Z","iopub.status.idle":"2021-07-16T03:30:50.959809Z","shell.execute_reply.started":"2021-07-16T03:30:50.911298Z","shell.execute_reply":"2021-07-16T03:30:50.958303Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tst[phone_col].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:50.961389Z","iopub.execute_input":"2021-07-16T03:30:50.961737Z","iopub.status.idle":"2021-07-16T03:30:51.003700Z","shell.execute_reply.started":"2021-07-16T03:30:50.961703Z","shell.execute_reply":"2021-07-16T03:30:51.002351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Each phone has fair amount of data points ranging between 577 and 3,517.","metadata":{}},{"cell_type":"code","source":"overlapping_phones = [x for x in tst[phone_col] if x in trn[phone_col]]\nprint(len(overlapping_phones))","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:51.005503Z","iopub.execute_input":"2021-07-16T03:30:51.005947Z","iopub.status.idle":"2021-07-16T03:30:51.869078Z","shell.execute_reply.started":"2021-07-16T03:30:51.005906Z","shell.execute_reply":"2021-07-16T03:30:51.867863Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"There's **no** overlapping phone between the training and test data.","metadata":{}},{"cell_type":"markdown","source":"`millisSinceGpsEpoch` is used as a timestamp. Its \"an integer number of milliseconds since the GPS epoch (1980/1/6 midnight UTC). Its value equals\".","metadata":{}},{"cell_type":"code","source":"tst[ts_col].min(), tst[ts_col].max()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:51.870707Z","iopub.execute_input":"2021-07-16T03:30:51.871091Z","iopub.status.idle":"2021-07-16T03:30:51.881810Z","shell.execute_reply.started":"2021-07-16T03:30:51.871038Z","shell.execute_reply":"2021-07-16T03:30:51.878379Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can convert them to `datatime64` using `pd.to_datetime()` as follows:","metadata":{}},{"cell_type":"code","source":"dt_offset = pd.to_datetime('1980-01-06 00:00:00')\nprint(dt_offset)\ndt_offset_in_ms = int(dt_offset.value / 1e6)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:51.884421Z","iopub.execute_input":"2021-07-16T03:30:51.884860Z","iopub.status.idle":"2021-07-16T03:30:51.903794Z","shell.execute_reply.started":"2021-07-16T03:30:51.884818Z","shell.execute_reply":"2021-07-16T03:30:51.902249Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn[dt_col] = pd.to_datetime(trn[ts_col] + dt_offset_in_ms, unit='ms')\ntst[dt_col] = pd.to_datetime(tst[ts_col] + dt_offset_in_ms, unit='ms')\nprint(f'Training data range: {trn[dt_col].min()} - {trn[dt_col].max()}')\nprint(f'    Test data range: {tst[dt_col].min()} - {tst[dt_col].max()}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:51.905442Z","iopub.execute_input":"2021-07-16T03:30:51.905825Z","iopub.status.idle":"2021-07-16T03:30:51.954271Z","shell.execute_reply.started":"2021-07-16T03:30:51.905787Z","shell.execute_reply":"2021-07-16T03:30:51.952639Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"First, let's see how estimated locations between the training and test data look like. The ground truth for training data is available per `phone` in `{collectionName}/{phoneName}/ground_truth.csv`.","metadata":{}},{"cell_type":"code","source":"latlon_trn = trn[[lat_col, lon_col]].round(3)\nlatlon_trn['counts'] = 1\nlatlon_trn = latlon_trn.groupby([lat_col, lon_col]).sum().reset_index()\nlatlon_trn.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:51.956648Z","iopub.execute_input":"2021-07-16T03:30:51.958392Z","iopub.status.idle":"2021-07-16T03:30:52.001282Z","shell.execute_reply.started":"2021-07-16T03:30:51.958325Z","shell.execute_reply":"2021-07-16T03:30:51.999260Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see the heatmap for the training data.","metadata":{}},{"cell_type":"code","source":"simple_folium(latlon_trn, lat_col, lon_col)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:52.003507Z","iopub.execute_input":"2021-07-16T03:30:52.004044Z","iopub.status.idle":"2021-07-16T03:30:54.026881Z","shell.execute_reply.started":"2021-07-16T03:30:52.003985Z","shell.execute_reply":"2021-07-16T03:30:54.024806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see the heatmap for the test data too.","metadata":{}},{"cell_type":"code","source":"latlon_tst = tst[[lat_col, lon_col]].round(3)\nlatlon_tst['counts'] = 1\nlatlon_tst = latlon_tst.groupby([lat_col, lon_col]).sum().reset_index()\n\nsimple_folium(latlon_tst, lat_col, lon_col)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:54.035065Z","iopub.execute_input":"2021-07-16T03:30:54.035578Z","iopub.status.idle":"2021-07-16T03:30:55.856874Z","shell.execute_reply.started":"2021-07-16T03:30:54.035532Z","shell.execute_reply":"2021-07-16T03:30:55.855365Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Phone Level Data","metadata":{}},{"cell_type":"markdown","source":"Read GNSS Logs","metadata":{}},{"cell_type":"code","source":"cname = trn[cname_col][0]\npname = trn[pname_col][0]\ndfs = gnss = gnss_log_to_dataframes(str(data_dir / 'train' / cname / pname / f'{pname}_GnssLog.txt'))\nprint(dfs.keys())","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:30:55.859925Z","iopub.execute_input":"2021-07-16T03:30:55.860330Z","iopub.status.idle":"2021-07-16T03:31:20.251423Z","shell.execute_reply.started":"2021-07-16T03:30:55.860292Z","shell.execute_reply":"2021-07-16T03:31:20.249808Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_raw = dfs['Raw']\nprint(df_raw.shape)\ndf_raw.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:20.253309Z","iopub.execute_input":"2021-07-16T03:31:20.253805Z","iopub.status.idle":"2021-07-16T03:31:20.295032Z","shell.execute_reply.started":"2021-07-16T03:31:20.253753Z","shell.execute_reply":"2021-07-16T03:31:20.293457Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_raw.info()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:20.297338Z","iopub.execute_input":"2021-07-16T03:31:20.297725Z","iopub.status.idle":"2021-07-16T03:31:20.331212Z","shell.execute_reply.started":"2021-07-16T03:31:20.297687Z","shell.execute_reply":"2021-07-16T03:31:20.329957Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From the [post](https://www.kaggle.com/c/google-smartphone-decimeter-challenge/discussion/238583) by @sohier and [slides](https://www.kaggle.com/google/android-smartphones-high-accuracy-datasets?select=ION+GNSS+2020+Slides+Android+Raw+GNSS+Measurement+Datasets+for+Precise+Positioning.pdf) by the data provider: \n\nMeasurements from GNSS chipsets of mobile phones are often noisier and more erroneous. Example of filters your can apply (to exclude) are:\n1. `FullBiasNanos` (GNSS Raw) is zero or invalid\n2. `BiasUncertaintyNanos` (GNSS Raw) too large (> 1e6)\n3. Arrival time is negative or unrealistically large - can be calculated from `rawPrM` (Derived)\n4. Unknown constellation (`constellationType == 0`) (Derived, GNSS Raw)\n5. `TimeNanos` is empty (GNSS Raw)\n6. `State` is not in (`STATE_TOW_DECODED`, `STATE_TOW_KNOWN`, `STATE_GLO_TOD_DECODED`, `STATE_GLO_TOD_KNOWN`) (GNSS Raw)\n7. `ReceivedSvTimeUncertaintyNanos` is high (500 ns) (GNSS Raw)\n8. `AccumulatedDeltaRangeState` violating this condition: `ADR_STATE_VALID == 1 & ADR_STATE_RESET == 0 & ADR_STATE_CYCLE_SLIP == 0` (GNSS Raw)\n9. `AccumulatedDeltaRangeUncertaintyMeters` is high (GNSS Raw)\n10. `Cn0DbHz` is less than 20 db-Hz (GNSS Raw)","metadata":{}},{"cell_type":"code","source":"df_raw['ArrivalTime'] = df_raw['TimeNanos'] - df_raw['FullBiasNanos'] - df_raw['BiasNanos']\nprint(df_raw['ArrivalTime'].describe())\ndf_raw['ArrivalTime'].hist(bins=20)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:20.332691Z","iopub.execute_input":"2021-07-16T03:31:20.333035Z","iopub.status.idle":"2021-07-16T03:31:20.615987Z","shell.execute_reply.started":"2021-07-16T03:31:20.332997Z","shell.execute_reply":"2021-07-16T03:31:20.614368Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_raw['BiasUncertaintyNanos'].describe())\ndf_raw['BiasUncertaintyNanos'].hist(bins=20)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:20.617542Z","iopub.execute_input":"2021-07-16T03:31:20.618003Z","iopub.status.idle":"2021-07-16T03:31:20.863464Z","shell.execute_reply.started":"2021-07-16T03:31:20.617963Z","shell.execute_reply":"2021-07-16T03:31:20.861510Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_raw['ReceivedSvTimeUncertaintyNanos'].describe())\ndf_raw['ReceivedSvTimeUncertaintyNanos'].hist(bins=20)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:20.864867Z","iopub.execute_input":"2021-07-16T03:31:20.865209Z","iopub.status.idle":"2021-07-16T03:31:21.125878Z","shell.execute_reply.started":"2021-07-16T03:31:20.865176Z","shell.execute_reply":"2021-07-16T03:31:21.124488Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_raw.AccumulatedDeltaRangeUncertaintyMeters.describe())\ndf_raw.AccumulatedDeltaRangeUncertaintyMeters.hist(bins=20)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:21.127521Z","iopub.execute_input":"2021-07-16T03:31:21.127859Z","iopub.status.idle":"2021-07-16T03:31:21.396651Z","shell.execute_reply.started":"2021-07-16T03:31:21.127826Z","shell.execute_reply":"2021-07-16T03:31:21.393497Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_raw.Cn0DbHz.describe())\ndf_raw.Cn0DbHz.hist(bins=20)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:21.397958Z","iopub.execute_input":"2021-07-16T03:31:21.398459Z","iopub.status.idle":"2021-07-16T03:31:21.644139Z","shell.execute_reply.started":"2021-07-16T03:31:21.398407Z","shell.execute_reply":"2021-07-16T03:31:21.642636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_raw = df_raw.loc[\n    ~pd.isnull(df_raw.FullBiasNanos) &\n    (df_raw.BiasUncertaintyNanos < 100) &\n    (df_raw.ArrivalTime > 0) &\n    (df_raw.ConstellationType != 0) &\n    ~pd.isnull(df_raw.TimeNanos) &\n    (df_raw.State != 3) & (df_raw.State != 14) & (df_raw.State != 7) & (df_raw.State != 15) &\n    (df_raw.ReceivedSvTimeUncertaintyNanos < 100) &\n    (df_raw.AccumulatedDeltaRangeUncertaintyMeters < 0.3) &\n    (df_raw.Cn0DbHz > 20)\n]\nprint(df_raw.shape)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:21.645407Z","iopub.execute_input":"2021-07-16T03:31:21.645730Z","iopub.status.idle":"2021-07-16T03:31:21.699958Z","shell.execute_reply.started":"2021-07-16T03:31:21.645699Z","shell.execute_reply":"2021-07-16T03:31:21.698872Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"See organizer's [Loading GNSS logs](https://www.kaggle.com/sohier/loading-gnss-logs) notebook for more details.","metadata":{}},{"cell_type":"markdown","source":"Derived values are used to generate baseline location estimates in `baseline_locations_{train|test}.csv`.","metadata":{}},{"cell_type":"code","source":"derived = pd.read_csv(data_dir / 'train' / cname / pname / f'{pname}_derived.csv')\nprint(derived.shape)\nderived.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:21.701368Z","iopub.execute_input":"2021-07-16T03:31:21.701839Z","iopub.status.idle":"2021-07-16T03:31:22.002509Z","shell.execute_reply.started":"2021-07-16T03:31:21.701788Z","shell.execute_reply":"2021-07-16T03:31:22.000940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived.info()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:22.004838Z","iopub.execute_input":"2021-07-16T03:31:22.005388Z","iopub.status.idle":"2021-07-16T03:31:22.042808Z","shell.execute_reply.started":"2021-07-16T03:31:22.005332Z","shell.execute_reply":"2021-07-16T03:31:22.041843Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived = derived.loc[derived.constellationType != 0]\nprint(derived.shape)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:22.043894Z","iopub.execute_input":"2021-07-16T03:31:22.044309Z","iopub.status.idle":"2021-07-16T03:31:22.059241Z","shell.execute_reply.started":"2021-07-16T03:31:22.044261Z","shell.execute_reply":"2021-07-16T03:31:22.057950Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's calculate `correctedPrM` as described in the data description:\n```\ncorrectedPrM = rawPrM + satClkBiasM - isrbM - ionoDelayM - tropoDelayM\n```\n\"The baseline locations are computed using correctedPrM and the satellite positions, using a standard Weighted Least Squares (WLS) solver, with the phone's position (x, y, z), clock bias (t), and isrbM for each unique signal type as states for each epoch.\"","metadata":{}},{"cell_type":"code","source":"derived['correctedPrM'] = (derived['rawPrM'] + derived['satClkBiasM'] - derived['isrbM'] - \n                           derived['ionoDelayM'] - derived['tropoDelayM'])\nsns.pairplot(data=derived, vars=['correctedPrM', 'rawPrM'], size=3)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:22.060964Z","iopub.execute_input":"2021-07-16T03:31:22.061369Z","iopub.status.idle":"2021-07-16T03:31:24.230986Z","shell.execute_reply.started":"2021-07-16T03:31:22.061330Z","shell.execute_reply":"2021-07-16T03:31:24.228806Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived[dt_col] = pd.to_datetime(derived[ts_col] + dt_offset_in_ms, unit='ms')\nprint(f'Data range for {cname}/{pname}: {derived[dt_col].min()} - {derived[dt_col].max()}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.232474Z","iopub.execute_input":"2021-07-16T03:31:24.232806Z","iopub.status.idle":"2021-07-16T03:31:24.250044Z","shell.execute_reply.started":"2021-07-16T03:31:24.232774Z","shell.execute_reply":"2021-07-16T03:31:24.248636Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The data is for 30 minutes or 1,800 seconds. However, we have a lot more samples (55K). This is because, for each second, there are multiple samples with different `constellationType`, `svid`, and `signalType`.","metadata":{}},{"cell_type":"code","source":"derived[['constellationType', 'svid', 'signalType']].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.251285Z","iopub.execute_input":"2021-07-16T03:31:24.251631Z","iopub.status.idle":"2021-07-16T03:31:24.282451Z","shell.execute_reply.started":"2021-07-16T03:31:24.251595Z","shell.execute_reply":"2021-07-16T03:31:24.280812Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived[[ts_col, 'constellationType', 'correctedPrM']].groupby([ts_col, 'constellationType']).agg(['mean', 'std', 'count']).describe()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.284140Z","iopub.execute_input":"2021-07-16T03:31:24.284476Z","iopub.status.idle":"2021-07-16T03:31:24.338785Z","shell.execute_reply.started":"2021-07-16T03:31:24.284444Z","shell.execute_reply":"2021-07-16T03:31:24.337528Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived.loc[derived.constellationType == 1][[ts_col, 'svid', 'correctedPrM']].groupby([ts_col, 'svid']).agg(['mean', 'std', 'count']).describe()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.340006Z","iopub.execute_input":"2021-07-16T03:31:24.340509Z","iopub.status.idle":"2021-07-16T03:31:24.414938Z","shell.execute_reply.started":"2021-07-16T03:31:24.340458Z","shell.execute_reply":"2021-07-16T03:31:24.412749Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Each epoch, given the constellation type of `1` (or GPS), from the same satellite, `coorectedPrM` can be different - because of different signal types.","metadata":{}},{"cell_type":"code","source":"derived.loc[derived.signalType == 'GPS_L1'][[ts_col, 'svid', 'correctedPrM']].groupby([ts_col, 'svid']).agg(['mean', 'std', 'count'])","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.418007Z","iopub.execute_input":"2021-07-16T03:31:24.418517Z","iopub.status.idle":"2021-07-16T03:31:24.479407Z","shell.execute_reply.started":"2021-07-16T03:31:24.418476Z","shell.execute_reply":"2021-07-16T03:31:24.477382Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived.loc[derived.signalType == 'GPS_L1'][[ts_col, 'svid', 'correctedPrM']].groupby([ts_col, 'svid']).agg(['mean', 'std', 'count']).describe()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.480902Z","iopub.execute_input":"2021-07-16T03:31:24.481435Z","iopub.status.idle":"2021-07-16T03:31:24.556158Z","shell.execute_reply.started":"2021-07-16T03:31:24.481381Z","shell.execute_reply":"2021-07-16T03:31:24.554247Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Each epoch, given the signal type of `GPS_L1`, from the same satellite, `correctedPrM` is unique.","metadata":{}},{"cell_type":"code","source":"derived.loc[derived.signalType == 'GPS_L1'][[ts_col, 'svid']].drop_duplicates().groupby([ts_col]).agg(['mean', 'std', 'count']).describe()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.557771Z","iopub.execute_input":"2021-07-16T03:31:24.558627Z","iopub.status.idle":"2021-07-16T03:31:24.655957Z","shell.execute_reply.started":"2021-07-16T03:31:24.558572Z","shell.execute_reply":"2021-07-16T03:31:24.647833Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Each epoch, given the signal type of `GPS_L1`, there are signals from at least 3 satellites.","metadata":{}},{"cell_type":"code","source":"gps_l1 = derived.loc[derived.signalType == 'GPS_L1'][[ts_col, 'svid', 'correctedPrM']].drop_duplicates([ts_col, 'svid'])\nprint(gps_l1.shape)\ngps_l1.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.657242Z","iopub.execute_input":"2021-07-16T03:31:24.657669Z","iopub.status.idle":"2021-07-16T03:31:24.706561Z","shell.execute_reply.started":"2021-07-16T03:31:24.657636Z","shell.execute_reply":"2021-07-16T03:31:24.705111Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Read Ground Truth","metadata":{}},{"cell_type":"code","source":"label = pd.read_csv(data_dir / 'train' / cname / pname / 'ground_truth.csv')\nprint(label.shape)\nlabel.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.708545Z","iopub.execute_input":"2021-07-16T03:31:24.708933Z","iopub.status.idle":"2021-07-16T03:31:24.762594Z","shell.execute_reply.started":"2021-07-16T03:31:24.708894Z","shell.execute_reply":"2021-07-16T03:31:24.761308Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"In the `*derived.csv`, we have 55K rows, but in the `ground_truth.csv`, we only have 1,740 rows.","metadata":{}},{"cell_type":"code","source":"label[dt_col] = pd.to_datetime(label[ts_col] + dt_offset_in_ms, unit='ms')\nprint(f'Labels range for {cname}/{pname}: {label[dt_col].min()} - {label[dt_col].max()}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.764201Z","iopub.execute_input":"2021-07-16T03:31:24.764552Z","iopub.status.idle":"2021-07-16T03:31:24.782252Z","shell.execute_reply.started":"2021-07-16T03:31:24.764516Z","shell.execute_reply":"2021-07-16T03:31:24.779675Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Hmm, this is weird. The label data starts 1 second earlier than the derived data. This means that if we join the derived and label data, the first second will have NaNs for derived columns. Let's check another phone data.","metadata":{}},{"cell_type":"code","source":"cname = trn[cname_col][10]\npname = trn[pname_col][10]\nderived2 = pd.read_csv(data_dir / 'train' / cname / pname / f'{pname}_derived.csv')\nlabel2 = pd.read_csv(data_dir / 'train' / cname / pname / 'ground_truth.csv')\nprint(f\"Derived data starts at: {pd.to_datetime(derived2[ts_col].min() + dt_offset_in_ms, unit='ms')}\")\nprint(f\"  Label data starts at: {pd.to_datetime(label2[ts_col].min() + dt_offset_in_ms, unit='ms')}\")","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:24.784473Z","iopub.execute_input":"2021-07-16T03:31:24.785321Z","iopub.status.idle":"2021-07-16T03:31:25.121713Z","shell.execute_reply.started":"2021-07-16T03:31:24.785256Z","shell.execute_reply":"2021-07-16T03:31:25.119492Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"It's the same. We don't have the first second data in the derived data. Let's take a note and move on.","metadata":{}},{"cell_type":"markdown","source":"## Feature Generation\n\nAdd features columns useful in predictions.","metadata":{}},{"cell_type":"markdown","source":"### Label Data Aggregation","metadata":{}},{"cell_type":"markdown","source":"First, let's add previous latitude and longitude estimates as features.","metadata":{}},{"cell_type":"code","source":"trn.sort_values([phone_col, ts_col], inplace=True)\ntrn[['prev_lat']] = trn[lat_col].shift().where(trn[phone_col].eq(trn[phone_col].shift()))\ntrn[['prev_lon']] = trn[lon_col].shift().where(trn[phone_col].eq(trn[phone_col].shift()))\n\ntst.sort_values([phone_col, ts_col], inplace=True)\ntst[['prev_lat']] = tst[lat_col].shift().where(tst[phone_col].eq(tst[phone_col].shift()))\ntst[['prev_lon']] = tst[lon_col].shift().where(tst[phone_col].eq(tst[phone_col].shift()))\ntrn.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:25.124190Z","iopub.execute_input":"2021-07-16T03:31:25.125228Z","iopub.status.idle":"2021-07-16T03:31:25.318186Z","shell.execute_reply.started":"2021-07-16T03:31:25.125172Z","shell.execute_reply":"2021-07-16T03:31:25.317286Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://www.kaggle.com/jpmiller/baseline-from-host-data\nlabel_files = (data_dir / 'train').rglob('ground_truth.csv')\ncols = [phone_col, ts_col, lat_col, lon_col]\n\ndf_list = []\nfor t in tqdm(label_files, total=73):\n    label = pd.read_csv(t, usecols=[cname_col, pname_col, ts_col, lat_col, lon_col])\n    df_list.append(label)\n\ndf_label = pd.concat(df_list, ignore_index=True)\ndf_label[phone_col] = df_label[cname_col] + '_' + df_label[pname_col]\n\ndf = df_label.merge(trn[cols + ['prev_lat', 'prev_lon']], how='inner', on=[phone_col, ts_col], \n                    suffixes=('_gt', '')).drop([cname_col, pname_col], axis=1)\ndf['sSinceGpsEpoch'] = df[ts_col] // 1000\nprint(df.shape)\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:25.319340Z","iopub.execute_input":"2021-07-16T03:31:25.319811Z","iopub.status.idle":"2021-07-16T03:31:27.029777Z","shell.execute_reply.started":"2021-07-16T03:31:25.319776Z","shell.execute_reply":"2021-07-16T03:31:27.027837Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_tst = sub[[phone_col, ts_col]].merge(tst[[phone_col, ts_col, lat_col, lon_col, 'prev_lat', 'prev_lon']], \n                                        how='left', on=[phone_col, ts_col], suffixes=('', '_basepred'))\ndf_tst['sSinceGpsEpoch'] = df_tst[ts_col] // 1000\nprint(df_tst.shape)\ndf_tst.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:27.032034Z","iopub.execute_input":"2021-07-16T03:31:27.032399Z","iopub.status.idle":"2021-07-16T03:31:27.142366Z","shell.execute_reply.started":"2021-07-16T03:31:27.032365Z","shell.execute_reply":"2021-07-16T03:31:27.141239Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"### Derived Data Aggregation","metadata":{}},{"cell_type":"code","source":"derived_files = (data_dir / 'train').rglob('*_derived.csv')\ncols = [ts_col, 'svid', 'correctedPrM']\n\ndf_list = []\nfor t in tqdm(derived_files, total=73):\n    derived = pd.read_csv(t).drop_duplicates([ts_col, 'svid'])\n    derived['correctedPrM'] = (derived['rawPrM'] + derived['satClkBiasM'] - derived['isrbM'] - \n                               derived['ionoDelayM'] - derived['tropoDelayM'])\n    df_list.append(derived[[cname_col, pname_col, ts_col, 'svid', 'correctedPrM']])\n    \ndf_derived = pd.concat(df_list, ignore_index=True)\ndf_derived[phone_col] = df_derived[cname_col] + '_' + df_derived[pname_col]\ndf_derived.drop([cname_col, pname_col], axis=1, inplace=True)\n\nprint(df_derived.shape)\ndf_derived.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:27.143839Z","iopub.execute_input":"2021-07-16T03:31:27.144225Z","iopub.status.idle":"2021-07-16T03:31:48.526768Z","shell.execute_reply.started":"2021-07-16T03:31:27.144190Z","shell.execute_reply":"2021-07-16T03:31:48.525507Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_derived_pivot = pd.pivot_table(df_derived, \n                                  values='correctedPrM', \n                                  index=[phone_col, ts_col],\n                                  columns=['svid'],\n                                  aggfunc=np.mean)\ndf_derived_pivot.columns = [f'svid_{x}' for x in df_derived_pivot.columns]\ndf_derived_pivot.reset_index(inplace=True)\ndf_derived_pivot['sSinceGpsEpoch'] = df_derived_pivot[ts_col] // 1000\n\nprint(df_derived_pivot.shape)\ndf_derived_pivot.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:48.528194Z","iopub.execute_input":"2021-07-16T03:31:48.528553Z","iopub.status.idle":"2021-07-16T03:31:51.135598Z","shell.execute_reply.started":"2021-07-16T03:31:48.528515Z","shell.execute_reply":"2021-07-16T03:31:51.134451Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = df.merge(df_derived_pivot, how='left', on=[phone_col, 'sSinceGpsEpoch'], suffixes=['', '_2'])\ndf.drop(['sSinceGpsEpoch', ts_col + '_2'], axis=1, inplace=True)\nprint(df.shape)\ndf.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:51.137041Z","iopub.execute_input":"2021-07-16T03:31:51.137394Z","iopub.status.idle":"2021-07-16T03:31:51.389495Z","shell.execute_reply.started":"2021-07-16T03:31:51.137358Z","shell.execute_reply":"2021-07-16T03:31:51.388220Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df['d_lat'] = df['latDeg_gt'] - df[lat_col]\ndf['d_lon'] = df['lngDeg_gt'] - df[lon_col]\ndf[['d_lat', 'd_lon']].describe()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:51.392156Z","iopub.execute_input":"2021-07-16T03:31:51.392845Z","iopub.status.idle":"2021-07-16T03:31:51.459617Z","shell.execute_reply.started":"2021-07-16T03:31:51.392788Z","shell.execute_reply":"2021-07-16T03:31:51.458342Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived_files = (data_dir / 'test').rglob('*_derived.csv')\ncols = [ts_col, 'svid', 'correctedPrM']\n\ndf_list = []\nfor t in tqdm(derived_files, total=48):\n    derived = pd.read_csv(t)\n    derived['sSinceGpsEpoch'] = derived[ts_col] // 1000\n    derived.drop_duplicates(['sSinceGpsEpoch', 'svid'], inplace=True)\n    derived['correctedPrM'] = (derived['rawPrM'] + derived['satClkBiasM'] - derived['isrbM'] - \n                               derived['ionoDelayM'] - derived['tropoDelayM'])\n    df_list.append(derived[[cname_col, pname_col, 'sSinceGpsEpoch', 'svid', 'correctedPrM']])\n    \ndf_derived = pd.concat(df_list, ignore_index=True)\ndf_derived[phone_col] = df_derived[cname_col] + '_' + df_derived[pname_col]\ndf_derived.drop([cname_col, pname_col], axis=1, inplace=True)\n\ndf_derived_pivot = pd.pivot_table(df_derived, \n                                  values='correctedPrM', \n                                  index=[phone_col, 'sSinceGpsEpoch'],\n                                  columns=['svid'],\n                                  aggfunc=np.mean)\ndf_derived_pivot.columns = [f'svid_{x}' for x in df_derived_pivot.columns]\ndf_derived_pivot.reset_index(inplace=True)\n\ndf_tst = df_tst.merge(df_derived_pivot, how='left', \n                      on=[phone_col, 'sSinceGpsEpoch']).drop(['sSinceGpsEpoch'], axis=1)\nprint(df_tst.shape)\ndf_tst.head()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:31:51.461114Z","iopub.execute_input":"2021-07-16T03:31:51.461451Z","iopub.status.idle":"2021-07-16T03:32:05.557164Z","shell.execute_reply.started":"2021-07-16T03:31:51.461418Z","shell.execute_reply":"2021-07-16T03:32:05.555856Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_tst.describe()","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:05.559107Z","iopub.execute_input":"2021-07-16T03:32:05.559461Z","iopub.status.idle":"2021-07-16T03:32:05.822071Z","shell.execute_reply.started":"2021-07-16T03:32:05.559425Z","shell.execute_reply":"2021-07-16T03:32:05.820208Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Model Training","metadata":{}},{"cell_type":"code","source":"tpu = tf.distribute.cluster_resolver.TPUClusterResolver.connect()\ntpu_strategy = tf.distribute.experimental.TPUStrategy(tpu)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:05.824223Z","iopub.execute_input":"2021-07-16T03:32:05.824571Z","iopub.status.idle":"2021-07-16T03:32:11.385827Z","shell.execute_reply.started":"2021-07-16T03:32:05.824536Z","shell.execute_reply":"2021-07-16T03:32:11.384847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"feature_cols = [x for x in df_tst.columns if x not in [phone_col, ts_col]]\ntarget_cols = ['d_lat', 'd_lon']\ninput_dim = len(feature_cols)\noutput_dim = len(target_cols)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:11.387290Z","iopub.execute_input":"2021-07-16T03:32:11.387922Z","iopub.status.idle":"2021-07-16T03:32:11.394156Z","shell.execute_reply.started":"2021-07-16T03:32:11.387882Z","shell.execute_reply":"2021-07-16T03:32:11.392688Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"scaler = StandardScaler()\nlabel_scaler = StandardScaler()\nscaler.fit(pd.concat([df[feature_cols], df_tst[feature_cols]], axis=0).fillna(0).values)\nX = scaler.transform(df[feature_cols].fillna(0).values)\nX_tst = scaler.transform(df_tst[feature_cols].fillna(0).values)\nY = label_scaler.fit_transform(df[target_cols].values)\nprint(X.shape, Y.shape, X_tst.shape)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:11.395763Z","iopub.execute_input":"2021-07-16T03:32:11.396233Z","iopub.status.idle":"2021-07-16T03:32:11.804123Z","shell.execute_reply.started":"2021-07-16T03:32:11.396162Z","shell.execute_reply":"2021-07-16T03:32:11.802735Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def build_model():\n    inputs = keras.layers.Input((input_dim,))\n    x = keras.layers.Dense(128, activation='relu')(inputs)\n    x = keras.layers.BatchNormalization()(x)\n    x = keras.layers.Dense(128, activation='relu')(x)\n    x = keras.layers.Dropout(.3)(x)\n    \n#     ox = x\n        \n#     x = keras.layers.Dense(128, activation='relu')(x)\n#     x = keras.layers.BatchNormalization()(x)\n#     x = keras.layers.Dense(128, activation='relu')(x)\n#     x = keras.layers.Dropout(.3)(x)    \n    \n#     x = keras.layers.Add()([x, ox])\n    x = keras.layers.Add()([x])\n    \n    x = keras.layers.Dense(128, activation='relu')(x)\n    x = keras.layers.BatchNormalization()(x)\n    x = keras.layers.Dense(128, activation='relu')(x)\n    x = keras.layers.Dropout(.3)(x)\n    \n    outputs = keras.layers.Dense(output_dim, activation='linear')(x)\n    \n    model = keras.Model(inputs, outputs)\n    model.compile(optimizer=keras.optimizers.Adam(lrate), loss='mean_squared_error')\n    return model","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:11.805312Z","iopub.execute_input":"2021-07-16T03:32:11.805738Z","iopub.status.idle":"2021-07-16T03:32:11.817312Z","shell.execute_reply.started":"2021-07-16T03:32:11.805688Z","shell.execute_reply":"2021-07-16T03:32:11.815702Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"with tpu_strategy.scope():\n    model = build_model()\n    model.summary()","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-07-16T03:32:11.818672Z","iopub.execute_input":"2021-07-16T03:32:11.819159Z","iopub.status.idle":"2021-07-16T03:32:12.069892Z","shell.execute_reply.started":"2021-07-16T03:32:11.819121Z","shell.execute_reply":"2021-07-16T03:32:12.066352Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def scheduler(epoch, lr, warmup=5):\n    if epoch < warmup:\n        return lr * 1.5\n    else:\n        return lr * tf.math.exp(-.1)\n\nes = keras.callbacks.EarlyStopping(patience=n_stop, restore_best_weights=True)\nlr = keras.callbacks.LearningRateScheduler(scheduler)\n\ncv = KFold(n_splits=n_fold, shuffle=True, random_state=seed)\n\nP = np.zeros_like(Y, dtype=float)\nP_tst = np.zeros((X_tst.shape[0], output_dim), dtype=float)\nfor i, (i_trn, i_val) in enumerate(cv.split(X), 1):\n    print(f'Training for CV #{i}')\n    model = build_model()\n    history = model.fit(X[i_trn], Y[i_trn], validation_data=(X[i_val], Y[i_val]), \n                        epochs=epochs, batch_size=batch_size, callbacks=[es, lr], verbose=0)\n    P[i_val] = label_scaler.inverse_transform(model.predict(X[i_val]))\n    P_tst += label_scaler.inverse_transform(model.predict(X_tst)) / n_fold\n    \n    distance_i = calc_haversine(df.latDeg_gt.values[i_val], \n                                df.lngDeg_gt.values[i_val], \n                                P[i_val, 0] + df.latDeg.values[i_val], \n                                P[i_val, 1] + df.lngDeg.values[i_val]).mean()\n    print(f'CV #{i}: {np.percentile(distance_i, [50, 95])}')","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-07-16T03:32:12.071799Z","iopub.status.idle":"2021-07-16T03:32:12.072600Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(P.mean(axis=0), P_tst.mean(axis=0))\nnp.savetxt(predict_val_file, P, delimiter=',', fmt='%.6f')\nnp.savetxt(predict_tst_file, P_tst, delimiter=',', fmt='%.6f')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:12.074558Z","iopub.status.idle":"2021-07-16T03:32:12.075332Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"distance = calc_haversine(df.latDeg_gt, df.lngDeg_gt, P[:, 0] + df.latDeg, P[:, 1] + df.lngDeg)\nprint(f'CV All: {np.percentile(distance, [50, 95])}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:12.077019Z","iopub.status.idle":"2021-07-16T03:32:12.077790Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df.sort_values([phone_col, ts_col], inplace=True)\ndf_smoothed = df.copy()\ndf_smoothed[lat_col] = df[lat_col] + P[:, 0]\ndf_smoothed[lon_col] = df[lon_col] + P[:, 1]\ndf_smoothed = apply_kf_smoothing(df_smoothed)\ndistance = calc_haversine(df_smoothed.latDeg_gt, df_smoothed.lngDeg_gt, df_smoothed.latDeg, df_smoothed.lngDeg)\nprint(f'CV All (smoothed): {np.percentile(distance, [50, 95])}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:12.079547Z","iopub.status.idle":"2021-07-16T03:32:12.080298Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(history.history['lr'])","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:12.082291Z","iopub.status.idle":"2021-07-16T03:32:12.083037Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Submission\n\nPrepare submission file","metadata":{}},{"cell_type":"code","source":"distance_tst = calc_haversine(df_tst.latDeg, df_tst.lngDeg, P_tst[:, 0] + df_tst.latDeg, P_tst[:, 1] + df_tst.lngDeg)\nprint(f'CV All: {np.percentile(distance_tst, [50, 95])}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:12.084658Z","iopub.status.idle":"2021-07-16T03:32:12.085162Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_tst.sort_values([phone_col, ts_col], inplace=True)\ndf_tst_smoothed = df_tst.copy()\ndf_tst_smoothed[lat_col] = df_tst_smoothed[lat_col] + P_tst[:, 0]\ndf_tst_smoothed[lon_col] = df_tst_smoothed[lon_col] + P_tst[:, 1]\ndf_tst_smoothed = apply_kf_smoothing(df_tst_smoothed)\ndistance_tst = calc_haversine(df_tst.latDeg, df_tst.lngDeg, df_tst_smoothed.latDeg, df_tst_smoothed.lngDeg)\nprint(f'CV All (smoothed): {np.percentile(distance_tst, [50, 95])}')","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:12.086578Z","iopub.status.idle":"2021-07-16T03:32:12.087122Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_tst_smoothed[[phone_col, ts_col, lat_col, lon_col]].to_csv(submission_file, index=False)","metadata":{"execution":{"iopub.status.busy":"2021-07-16T03:32:12.088666Z","iopub.status.idle":"2021-07-16T03:32:12.089197Z"},"trusted":true},"execution_count":null,"outputs":[]}]}