{"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":"This notebook shares EDA on the aggregate and phone level data, then trains a neural network model to predict the residuals of base estimations provided with aggregated features.\n\nThe contents of the notebook are organized as follows:\n1. Aggregated Data EDA\n2. Phone Level Data EDA\n3. Feature Generation: generates aggregated features for training. Currently we only use previous lat/long and `correctedPrM` from derived files.\n4. Model Training: trains a neural network with a skip connection in Keras on TPU.\n\nCredits to other notebooks:\n* [Baseline from host data](https://www.kaggle.com/jpmiller/baseline-from-host-data) by @jpmiller: for the distance calculation with `calc_haversine()`\n* [Demonstration of the Kalman filter](https://www.kaggle.com/emaerthin/demonstration-of-the-kalman-filter) by @emaerthin: for Kalman filtering with `apply_kf_smoothing()`\n* [Loading GNSS logs](https://www.kaggle.com/sohier/loading-gnss-logs) by organizers: for GNSS log loading with `gnss_log_to_dataframes()`\n* [Ѫ Start Here: Simple Folium Heatmap for Geo-Data](https://www.kaggle.com/dannellyz/start-here-simple-folium-heatmap-for-geo-data) by @dannellyz: for geospatial heatmap with `simple_folium()`\n\nThanks for sharing.","metadata":{}},{"cell_type":"markdown","source":"# Load Libraries & Data","metadata":{}},{"cell_type":"code","source":"!pip install simdkalman","metadata":{"_kg_hide-output":true,"execution":{"iopub.status.busy":"2021-05-31T08:46:25.489371Z","iopub.execute_input":"2021-05-31T08:46:25.489920Z","iopub.status.idle":"2021-05-31T08:46:34.487426Z","shell.execute_reply.started":"2021-05-31T08:46:25.489834Z","shell.execute_reply":"2021-05-31T08:46:34.486212Z"},"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-05-31T08:46:34.489950Z","iopub.execute_input":"2021-05-31T08:46:34.490333Z","iopub.status.idle":"2021-05-31T08:46:41.928124Z","shell.execute_reply.started":"2021-05-31T08:46:34.490294Z","shell.execute_reply":"2021-05-31T08:46:41.926320Z"},"trusted":true},"execution_count":null,"outputs":[]},{"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 = .01\nbatch_size = 2048\nepochs = 150\nn_stop = 10\nn_fold = 25\nseed = 42","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:46:41.930878Z","iopub.execute_input":"2021-05-31T08:46:41.931375Z","iopub.status.idle":"2021-05-31T08:46:41.941483Z","shell.execute_reply.started":"2021-05-31T08:46:41.931329Z","shell.execute_reply":"2021-05-31T08:46:41.940268Z"},"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-05-31T08:46:41.943330Z","iopub.execute_input":"2021-05-31T08:46:41.943642Z","iopub.status.idle":"2021-05-31T08:46:41.958009Z","shell.execute_reply.started":"2021-05-31T08:46:41.943612Z","shell.execute_reply":"2021-05-31T08:46:41.956316Z"},"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-05-31T08:46:41.960014Z","iopub.execute_input":"2021-05-31T08:46:41.960458Z","iopub.status.idle":"2021-05-31T08:46:42.421684Z","shell.execute_reply.started":"2021-05-31T08:46:41.960424Z","shell.execute_reply":"2021-05-31T08:46:42.420120Z"},"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-05-31T08:46:42.423476Z","iopub.execute_input":"2021-05-31T08:46:42.423902Z","iopub.status.idle":"2021-05-31T08:46:42.431591Z","shell.execute_reply.started":"2021-05-31T08:46:42.423859Z","shell.execute_reply":"2021-05-31T08:46:42.430629Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# from https://www.kaggle.com/emaerthin/demonstration-of-the-kalman-filter\nT = 0.80\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-05-31T08:46:42.432928Z","iopub.execute_input":"2021-05-31T08:46:42.433539Z","iopub.status.idle":"2021-05-31T08:46:42.449521Z","shell.execute_reply.started":"2021-05-31T08:46:42.433504Z","shell.execute_reply":"2021-05-31T08:46:42.447669Z"},"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-05-31T08:46:42.454706Z","iopub.execute_input":"2021-05-31T08:46:42.455325Z","iopub.status.idle":"2021-05-31T08:46:42.860072Z","shell.execute_reply.started":"2021-05-31T08:46:42.455258Z","shell.execute_reply":"2021-05-31T08:46:42.858611Z"},"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-05-31T08:46:42.862368Z","iopub.execute_input":"2021-05-31T08:46:42.862685Z","iopub.status.idle":"2021-05-31T08:46:43.115657Z","shell.execute_reply.started":"2021-05-31T08:46:42.862653Z","shell.execute_reply":"2021-05-31T08:46:43.114252Z"},"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-05-31T08:46:43.116985Z","iopub.execute_input":"2021-05-31T08:46:43.117379Z","iopub.status.idle":"2021-05-31T08:46:43.290073Z","shell.execute_reply.started":"2021-05-31T08:46:43.117345Z","shell.execute_reply":"2021-05-31T08:46:43.288937Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# Aggregated Data EDA","metadata":{}},{"cell_type":"markdown","source":"## `collectionName`, `phoneName`","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-05-31T08:46:43.291370Z","iopub.execute_input":"2021-05-31T08:46:43.291645Z","iopub.status.idle":"2021-05-31T08:46:43.390025Z","shell.execute_reply.started":"2021-05-31T08:46:43.291619Z","shell.execute_reply":"2021-05-31T08:46:43.388974Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn[pname_col].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:46:43.391885Z","iopub.execute_input":"2021-05-31T08:46:43.392333Z","iopub.status.idle":"2021-05-31T08:46:43.431306Z","shell.execute_reply.started":"2021-05-31T08:46:43.392287Z","shell.execute_reply":"2021-05-31T08:46:43.429768Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tst[pname_col].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:46:43.432876Z","iopub.execute_input":"2021-05-31T08:46:43.433505Z","iopub.status.idle":"2021-05-31T08:46:43.642365Z","shell.execute_reply.started":"2021-05-31T08:46:43.433456Z","shell.execute_reply":"2021-05-31T08:46:43.641160Z"},"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-05-31T08:46:43.643897Z","iopub.execute_input":"2021-05-31T08:46:43.644432Z","iopub.status.idle":"2021-05-31T08:46:43.701650Z","shell.execute_reply.started":"2021-05-31T08:46:43.644397Z","shell.execute_reply":"2021-05-31T08:46:43.700093Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn[phone_col].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:46:43.703685Z","iopub.execute_input":"2021-05-31T08:46:43.704058Z","iopub.status.idle":"2021-05-31T08:46:43.742561Z","shell.execute_reply.started":"2021-05-31T08:46:43.704025Z","shell.execute_reply":"2021-05-31T08:46:43.741384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tst[phone_col].value_counts()","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:46:43.744260Z","iopub.execute_input":"2021-05-31T08:46:43.744563Z","iopub.status.idle":"2021-05-31T08:46:43.772338Z","shell.execute_reply.started":"2021-05-31T08:46:43.744532Z","shell.execute_reply":"2021-05-31T08:46:43.771594Z"},"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-05-31T08:46:43.773519Z","iopub.execute_input":"2021-05-31T08:46:43.773988Z","iopub.status.idle":"2021-05-31T08:46:44.556018Z","shell.execute_reply.started":"2021-05-31T08:46:43.773940Z","shell.execute_reply":"2021-05-31T08:46:44.554975Z"},"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`","metadata":{}},{"cell_type":"code","source":"tst[ts_col].min(), tst[ts_col].max()","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:46:44.557778Z","iopub.execute_input":"2021-05-31T08:46:44.558153Z","iopub.status.idle":"2021-05-31T08:46:44.566207Z","shell.execute_reply.started":"2021-05-31T08:46:44.558119Z","shell.execute_reply":"2021-05-31T08:46:44.565053Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"From the data description, `millisSinceGpsEpoch` is \"an integer number of milliseconds since the GPS epoch (1980/1/6 midnight UTC). Its value equals\". 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-05-31T08:46:44.567570Z","iopub.execute_input":"2021-05-31T08:46:44.567894Z","iopub.status.idle":"2021-05-31T08:46:44.583541Z","shell.execute_reply.started":"2021-05-31T08:46:44.567863Z","shell.execute_reply":"2021-05-31T08:46:44.582005Z"},"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-05-31T08:46:44.585180Z","iopub.execute_input":"2021-05-31T08:46:44.585494Z","iopub.status.idle":"2021-05-31T08:46:44.630839Z","shell.execute_reply.started":"2021-05-31T08:46:44.585465Z","shell.execute_reply":"2021-05-31T08:46:44.629329Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## `latDeg` and `lngDeg`","metadata":{}},{"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-05-31T08:46:44.632157Z","iopub.execute_input":"2021-05-31T08:46:44.632445Z","iopub.status.idle":"2021-05-31T08:46:44.674904Z","shell.execute_reply.started":"2021-05-31T08:46:44.632417Z","shell.execute_reply":"2021-05-31T08:46:44.673531Z"},"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-05-31T08:46:44.676575Z","iopub.execute_input":"2021-05-31T08:46:44.677025Z","iopub.status.idle":"2021-05-31T08:46:46.139338Z","shell.execute_reply.started":"2021-05-31T08:46:44.676984Z","shell.execute_reply":"2021-05-31T08:46:46.138491Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"Let's see the heatmap for the test data too.","metadata":{}},{"cell_type":"markdown","source":"# Phone Level Data EDA","metadata":{}},{"cell_type":"markdown","source":"## 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-05-31T08:46:46.144279Z","iopub.execute_input":"2021-05-31T08:46:46.144596Z","iopub.status.idle":"2021-05-31T08:47:07.931168Z","shell.execute_reply.started":"2021-05-31T08:46:46.144568Z","shell.execute_reply":"2021-05-31T08:47:07.930085Z"},"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-05-31T08:47:07.933798Z","iopub.execute_input":"2021-05-31T08:47:07.934201Z","iopub.status.idle":"2021-05-31T08:47:07.967945Z","shell.execute_reply.started":"2021-05-31T08:47:07.934160Z","shell.execute_reply":"2021-05-31T08:47:07.967232Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_raw.info()","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:47:07.968889Z","iopub.execute_input":"2021-05-31T08:47:07.969258Z","iopub.status.idle":"2021-05-31T08:47:07.994263Z","shell.execute_reply.started":"2021-05-31T08:47:07.969230Z","shell.execute_reply":"2021-05-31T08:47:07.993083Z"},"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-05-31T08:47:07.995647Z","iopub.execute_input":"2021-05-31T08:47:07.996034Z","iopub.status.idle":"2021-05-31T08:47:08.244209Z","shell.execute_reply.started":"2021-05-31T08:47:07.996003Z","shell.execute_reply":"2021-05-31T08:47:08.243112Z"},"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-05-31T08:47:08.245441Z","iopub.execute_input":"2021-05-31T08:47:08.245723Z","iopub.status.idle":"2021-05-31T08:47:08.418876Z","shell.execute_reply.started":"2021-05-31T08:47:08.245696Z","shell.execute_reply":"2021-05-31T08:47:08.417894Z"},"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-05-31T08:47:08.420293Z","iopub.execute_input":"2021-05-31T08:47:08.420710Z","iopub.status.idle":"2021-05-31T08:47:08.620686Z","shell.execute_reply.started":"2021-05-31T08:47:08.420668Z","shell.execute_reply":"2021-05-31T08:47:08.619630Z"},"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-05-31T08:47:08.622305Z","iopub.execute_input":"2021-05-31T08:47:08.622712Z","iopub.status.idle":"2021-05-31T08:47:08.807261Z","shell.execute_reply.started":"2021-05-31T08:47:08.622669Z","shell.execute_reply":"2021-05-31T08:47:08.806514Z"},"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-05-31T08:47:08.808316Z","iopub.execute_input":"2021-05-31T08:47:08.808750Z","iopub.status.idle":"2021-05-31T08:47:08.989743Z","shell.execute_reply.started":"2021-05-31T08:47:08.808698Z","shell.execute_reply":"2021-05-31T08:47:08.988952Z"},"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-05-31T08:47:08.990863Z","iopub.execute_input":"2021-05-31T08:47:08.991290Z","iopub.status.idle":"2021-05-31T08:47:09.020046Z","shell.execute_reply.started":"2021-05-31T08:47:08.991250Z","shell.execute_reply":"2021-05-31T08:47:09.018829Z"},"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","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-05-31T08:47:09.021637Z","iopub.execute_input":"2021-05-31T08:47:09.022087Z","iopub.status.idle":"2021-05-31T08:47:09.407176Z","shell.execute_reply.started":"2021-05-31T08:47:09.022045Z","shell.execute_reply":"2021-05-31T08:47:09.406150Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived.info()","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:47:09.410168Z","iopub.execute_input":"2021-05-31T08:47:09.410456Z","iopub.status.idle":"2021-05-31T08:47:09.444220Z","shell.execute_reply.started":"2021-05-31T08:47:09.410428Z","shell.execute_reply":"2021-05-31T08:47:09.442993Z"},"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-05-31T08:47:09.446268Z","iopub.execute_input":"2021-05-31T08:47:09.446794Z","iopub.status.idle":"2021-05-31T08:47:09.460538Z","shell.execute_reply.started":"2021-05-31T08:47:09.446748Z","shell.execute_reply":"2021-05-31T08:47:09.459503Z"},"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-05-31T08:47:09.461825Z","iopub.execute_input":"2021-05-31T08:47:09.462513Z","iopub.status.idle":"2021-05-31T08:47:11.002474Z","shell.execute_reply.started":"2021-05-31T08:47:09.462469Z","shell.execute_reply":"2021-05-31T08:47:11.001137Z"},"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-05-31T08:47:11.003869Z","iopub.execute_input":"2021-05-31T08:47:11.004174Z","iopub.status.idle":"2021-05-31T08:47:11.019560Z","shell.execute_reply.started":"2021-05-31T08:47:11.004144Z","shell.execute_reply":"2021-05-31T08:47:11.018241Z"},"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-05-31T08:47:11.021313Z","iopub.execute_input":"2021-05-31T08:47:11.021868Z","iopub.status.idle":"2021-05-31T08:47:11.047963Z","shell.execute_reply.started":"2021-05-31T08:47:11.021816Z","shell.execute_reply":"2021-05-31T08:47:11.047019Z"},"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-05-31T08:47:11.050725Z","iopub.execute_input":"2021-05-31T08:47:11.051198Z","iopub.status.idle":"2021-05-31T08:47:11.099509Z","shell.execute_reply.started":"2021-05-31T08:47:11.051154Z","shell.execute_reply":"2021-05-31T08:47:11.098634Z"},"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-05-31T08:47:11.100503Z","iopub.execute_input":"2021-05-31T08:47:11.100934Z","iopub.status.idle":"2021-05-31T08:47:11.150578Z","shell.execute_reply.started":"2021-05-31T08:47:11.100903Z","shell.execute_reply":"2021-05-31T08:47:11.149348Z"},"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-05-31T08:47:11.151900Z","iopub.execute_input":"2021-05-31T08:47:11.152199Z","iopub.status.idle":"2021-05-31T08:47:11.203050Z","shell.execute_reply.started":"2021-05-31T08:47:11.152169Z","shell.execute_reply":"2021-05-31T08:47:11.202243Z"},"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-05-31T08:47:11.204072Z","iopub.execute_input":"2021-05-31T08:47:11.204470Z","iopub.status.idle":"2021-05-31T08:47:11.260512Z","shell.execute_reply.started":"2021-05-31T08:47:11.204441Z","shell.execute_reply":"2021-05-31T08:47:11.259540Z"},"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-05-31T08:47:11.261968Z","iopub.execute_input":"2021-05-31T08:47:11.262271Z","iopub.status.idle":"2021-05-31T08:47:11.312755Z","shell.execute_reply.started":"2021-05-31T08:47:11.262242Z","shell.execute_reply":"2021-05-31T08:47:11.311755Z"},"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-05-31T08:47:11.313997Z","iopub.execute_input":"2021-05-31T08:47:11.314280Z","iopub.status.idle":"2021-05-31T08:47:11.346943Z","shell.execute_reply.started":"2021-05-31T08:47:11.314253Z","shell.execute_reply":"2021-05-31T08:47:11.345672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 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-05-31T08:47:11.348406Z","iopub.execute_input":"2021-05-31T08:47:11.348862Z","iopub.status.idle":"2021-05-31T08:47:11.387966Z","shell.execute_reply.started":"2021-05-31T08:47:11.348824Z","shell.execute_reply":"2021-05-31T08:47:11.386923Z"},"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-05-31T08:47:11.389828Z","iopub.execute_input":"2021-05-31T08:47:11.390218Z","iopub.status.idle":"2021-05-31T08:47:11.400208Z","shell.execute_reply.started":"2021-05-31T08:47:11.390183Z","shell.execute_reply":"2021-05-31T08:47:11.398360Z"},"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-05-31T08:47:11.401893Z","iopub.execute_input":"2021-05-31T08:47:11.402297Z","iopub.status.idle":"2021-05-31T08:47:11.651359Z","shell.execute_reply.started":"2021-05-31T08:47:11.402263Z","shell.execute_reply":"2021-05-31T08:47:11.650176Z"},"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","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-05-31T08:47:11.652475Z","iopub.execute_input":"2021-05-31T08:47:11.652772Z","iopub.status.idle":"2021-05-31T08:47:11.810992Z","shell.execute_reply.started":"2021-05-31T08:47:11.652745Z","shell.execute_reply":"2021-05-31T08:47:11.809642Z"},"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-05-31T08:47:11.812552Z","iopub.execute_input":"2021-05-31T08:47:11.813320Z","iopub.status.idle":"2021-05-31T08:47:13.932930Z","shell.execute_reply.started":"2021-05-31T08:47:11.813276Z","shell.execute_reply":"2021-05-31T08:47:13.931148Z"},"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-05-31T08:47:13.934217Z","iopub.execute_input":"2021-05-31T08:47:13.934502Z","iopub.status.idle":"2021-05-31T08:47:14.012148Z","shell.execute_reply.started":"2021-05-31T08:47:13.934475Z","shell.execute_reply":"2021-05-31T08:47:14.010840Z"},"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-05-31T08:47:14.013388Z","iopub.execute_input":"2021-05-31T08:47:14.013680Z","iopub.status.idle":"2021-05-31T08:47:37.614861Z","shell.execute_reply.started":"2021-05-31T08:47:14.013650Z","shell.execute_reply":"2021-05-31T08:47:37.613244Z"},"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-05-31T08:47:37.616327Z","iopub.execute_input":"2021-05-31T08:47:37.616740Z","iopub.status.idle":"2021-05-31T08:47:40.229927Z","shell.execute_reply.started":"2021-05-31T08:47:37.616684Z","shell.execute_reply":"2021-05-31T08:47:40.228804Z"},"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-05-31T08:47:40.230996Z","iopub.execute_input":"2021-05-31T08:47:40.231290Z","iopub.status.idle":"2021-05-31T08:47:40.502649Z","shell.execute_reply.started":"2021-05-31T08:47:40.231261Z","shell.execute_reply":"2021-05-31T08:47:40.501397Z"},"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-05-31T08:47:40.504334Z","iopub.execute_input":"2021-05-31T08:47:40.504639Z","iopub.status.idle":"2021-05-31T08:47:40.586897Z","shell.execute_reply.started":"2021-05-31T08:47:40.504609Z","shell.execute_reply":"2021-05-31T08:47:40.585285Z"},"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-05-31T08:47:40.588098Z","iopub.execute_input":"2021-05-31T08:47:40.588429Z","iopub.status.idle":"2021-05-31T08:47:58.235835Z","shell.execute_reply.started":"2021-05-31T08:47:40.588400Z","shell.execute_reply":"2021-05-31T08:47:58.234188Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_tst.describe()","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:47:58.237068Z","iopub.execute_input":"2021-05-31T08:47:58.237347Z","iopub.status.idle":"2021-05-31T08:47:58.547456Z","shell.execute_reply.started":"2021-05-31T08:47:58.237320Z","shell.execute_reply":"2021-05-31T08:47:58.546383Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Raw Data Aggregation - To Be Updated","metadata":{}},{"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-05-31T08:47:58.549021Z","iopub.execute_input":"2021-05-31T08:47:58.549360Z","iopub.status.idle":"2021-05-31T08:48:03.699543Z","shell.execute_reply.started":"2021-05-31T08:47:58.549330Z","shell.execute_reply":"2021-05-31T08:48:03.698600Z"},"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-05-31T08:48:03.701008Z","iopub.execute_input":"2021-05-31T08:48:03.701416Z","iopub.status.idle":"2021-05-31T08:48:03.707010Z","shell.execute_reply.started":"2021-05-31T08:48:03.701375Z","shell.execute_reply":"2021-05-31T08:48:03.705672Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test = df_tst\nfor i in tqdm(range(1,len(df_test))):\n    lat1 = df_test.loc[i-1,'latDeg']\n    lon1 = df_test.loc[i-1,'lngDeg']\n    lat2 = df_test.loc[i,'latDeg']\n    lon2 = df_test.loc[i,'lngDeg']\n    df_test.loc[i,'dist_pre'] = calc_haversine(lat1, lon1, lat2, lon2)\n\nlist_phone = df_test['phone'].unique()\nfor phone in list_phone:\n    ind = df_test[df_test['phone'] == phone].index[0]\n    df_test.loc[ind,'dist_pre'] = 0","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:48:03.708542Z","iopub.execute_input":"2021-05-31T08:48:03.709252Z","iopub.status.idle":"2021-05-31T08:49:08.149236Z","shell.execute_reply.started":"2021-05-31T08:48:03.709208Z","shell.execute_reply":"2021-05-31T08:49:08.148467Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in tqdm(range(0,len(df_test)-1)):\n    lat1 = df_test.loc[i,'latDeg']\n    lon1 = df_test.loc[i,'lngDeg']\n    lat2 = df_test.loc[i+1,'latDeg']\n    lon2 = df_test.loc[i+1,'lngDeg']\n    df_test.loc[i,'dist_pro'] = calc_haversine(lat1, lon1, lat2, lon2)\nlist_phone = df_test['phone'].unique()\nfor phone in list_phone:\n    ind = df_test[df_test['phone'] == phone].index[-1]\n    df_test.loc[ind,'dist_pro'] = 0","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:49:15.884136Z","iopub.execute_input":"2021-05-31T08:49:15.884671Z","iopub.status.idle":"2021-05-31T08:50:18.714994Z","shell.execute_reply.started":"2021-05-31T08:49:15.884624Z","shell.execute_reply":"2021-05-31T08:50:18.713657Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pro_95 = df_test['dist_pro'].mean() + (df_test['dist_pro'].std() * 2)\npre_95 = df_test['dist_pre'].mean() + (df_test['dist_pre'].std() * 2)\nind = df_test[(df_test['dist_pro'] > pro_95)&(df_test['dist_pre'] > pre_95)][['dist_pre','dist_pro']].index\n\nfor i in ind:\n    df_test.loc[i,'latDeg'] = (df_test.loc[i-1,'latDeg'] + df_test.loc[i+1,'latDeg'])/2\n    df_test.loc[i,'lngDeg'] = (df_test.loc[i-1,'lngDeg'] + df_test.loc[i+1,'lngDeg'])/2","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:50:43.159204Z","iopub.execute_input":"2021-05-31T08:50:43.159679Z","iopub.status.idle":"2021-05-31T08:50:43.315158Z","shell.execute_reply.started":"2021-05-31T08:50:43.159647Z","shell.execute_reply":"2021-05-31T08:50:43.314384Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_tst = df_test","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:51:32.798329Z","iopub.execute_input":"2021-05-31T08:51:32.798845Z","iopub.status.idle":"2021-05-31T08:51:32.802427Z","shell.execute_reply.started":"2021-05-31T08:51:32.798811Z","shell.execute_reply":"2021-05-31T08:51:32.801344Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_test = df\nfor i in tqdm(range(1,len(df_test))):\n    lat1 = df_test.loc[i-1,'latDeg']\n    lon1 = df_test.loc[i-1,'lngDeg']\n    lat2 = df_test.loc[i,'latDeg']\n    lon2 = df_test.loc[i,'lngDeg']\n    df_test.loc[i,'dist_pre'] = calc_haversine(lat1, lon1, lat2, lon2)\n\nlist_phone = df_test['phone'].unique()\nfor phone in list_phone:\n    ind = df_test[df_test['phone'] == phone].index[0]\n    df_test.loc[ind,'dist_pre'] = 0","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:52:10.730025Z","iopub.execute_input":"2021-05-31T08:52:10.730644Z","iopub.status.idle":"2021-05-31T08:53:53.338214Z","shell.execute_reply.started":"2021-05-31T08:52:10.730608Z","shell.execute_reply":"2021-05-31T08:53:53.336847Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"for i in tqdm(range(0,len(df_test)-1)):\n    lat1 = df_test.loc[i,'latDeg']\n    lon1 = df_test.loc[i,'lngDeg']\n    lat2 = df_test.loc[i+1,'latDeg']\n    lon2 = df_test.loc[i+1,'lngDeg']\n    df_test.loc[i,'dist_pro'] = calc_haversine(lat1, lon1, lat2, lon2)\nlist_phone = df_test['phone'].unique()\nfor phone in list_phone:\n    ind = df_test[df_test['phone'] == phone].index[-1]\n    df_test.loc[ind,'dist_pro'] = 0","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:53:53.433078Z","iopub.execute_input":"2021-05-31T08:53:53.433462Z","iopub.status.idle":"2021-05-31T08:55:38.799442Z","shell.execute_reply.started":"2021-05-31T08:53:53.433432Z","shell.execute_reply":"2021-05-31T08:55:38.798450Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"pro_95 = df_test['dist_pro'].mean() + (df_test['dist_pro'].std() * 2)\npre_95 = df_test['dist_pre'].mean() + (df_test['dist_pre'].std() * 2)\nind = df_test[(df_test['dist_pro'] > pro_95)&(df_test['dist_pre'] > pre_95)][['dist_pre','dist_pro']].index\n\nfor i in ind:\n    df_test.loc[i,'latDeg'] = (df_test.loc[i-1,'latDeg'] + df_test.loc[i+1,'latDeg'])/2\n    df_test.loc[i,'lngDeg'] = (df_test.loc[i-1,'lngDeg'] + df_test.loc[i+1,'lngDeg'])/2","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:55:43.339391Z","iopub.execute_input":"2021-05-31T08:55:43.339767Z","iopub.status.idle":"2021-05-31T08:55:43.399354Z","shell.execute_reply.started":"2021-05-31T08:55:43.339711Z","shell.execute_reply":"2021-05-31T08:55:43.398351Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df = df_test","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:55:45.168645Z","iopub.execute_input":"2021-05-31T08:55:45.169026Z","iopub.status.idle":"2021-05-31T08:55:45.173565Z","shell.execute_reply.started":"2021-05-31T08:55:45.168994Z","shell.execute_reply":"2021-05-31T08:55:45.172194Z"},"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-05-31T08:55:49.605457Z","iopub.execute_input":"2021-05-31T08:55:49.605993Z","iopub.status.idle":"2021-05-31T08:55:49.968260Z","shell.execute_reply.started":"2021-05-31T08:55:49.605918Z","shell.execute_reply":"2021-05-31T08:55:49.967155Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"from keras import backend as K","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:55:51.747725Z","iopub.execute_input":"2021-05-31T08:55:51.748082Z","iopub.status.idle":"2021-05-31T08:55:51.824701Z","shell.execute_reply.started":"2021-05-31T08:55:51.748054Z","shell.execute_reply":"2021-05-31T08:55:51.823284Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def nll1(y_true, y_pred):\n    \"\"\" Negative log likelihood. \"\"\"\n\n    # keras.losses.binary_crossentropy give the mean\n    # over the last axis. we require the sum\n    return K.sum(K.binary_crossentropy(y_true, y_pred), axis=-1)","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:55:52.073230Z","iopub.execute_input":"2021-05-31T08:55:52.073636Z","iopub.status.idle":"2021-05-31T08:55:52.079434Z","shell.execute_reply.started":"2021-05-31T08:55:52.073601Z","shell.execute_reply":"2021-05-31T08:55:52.078354Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def loss(y_true, y_pred):\n    PI_ON_180 = tf.constant(np.pi / 180, dtype=tf.float32)\n    RADIUS_M = tf.constant(6_377_000, dtype = tf.float32)\n    tf.dtypes.cast(y_true, tf.float32)\n    tf.dtypes.cast(y_pred, tf.float32)\n\n    yt_rad = y_true * PI_ON_180\n    yp_rad = y_pred * PI_ON_180\n\n    delta = yt_rad - yp_rad\n    v = delta / 2\n    v = tf.sin(v)\n    v = v**2\n\n    a = v[:,1] + tf.cos(yt_rad[:,1]) * tf.cos(yp_rad[:,1]) * v[:,0] \n    c = tf.sqrt(a)\n    c = 2* tf.math.asin(c)\n    c = c*RADIUS_M\n\n    final = tf.reduce_mean(c)\n    return final","metadata":{"execution":{"iopub.status.busy":"2021-05-31T08:55:52.509018Z","iopub.execute_input":"2021-05-31T08:55:52.509808Z","iopub.status.idle":"2021-05-31T08:55:52.519767Z","shell.execute_reply.started":"2021-05-31T08:55:52.509754Z","shell.execute_reply":"2021-05-31T08:55:52.518864Z"},"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    \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.Ftrl(lrate), loss=nll1)\n    return model\n","metadata":{"execution":{"iopub.status.busy":"2021-05-31T09:57:41.920757Z","iopub.execute_input":"2021-05-31T09:57:41.921088Z","iopub.status.idle":"2021-05-31T09:57:41.930950Z","shell.execute_reply.started":"2021-05-31T09:57:41.921060Z","shell.execute_reply":"2021-05-31T09:57:41.929849Z"},"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-05-31T09:57:42.503512Z","iopub.execute_input":"2021-05-31T09:57:42.503877Z","iopub.status.idle":"2021-05-31T09:57:42.964885Z","shell.execute_reply.started":"2021-05-31T09:57:42.503847Z","shell.execute_reply":"2021-05-31T09:57:42.963914Z"},"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-05-31T09:57:45.680271Z","iopub.execute_input":"2021-05-31T09:57:45.680618Z","iopub.status.idle":"2021-05-31T10:03:59.265278Z","shell.execute_reply.started":"2021-05-31T09:57:45.680590Z","shell.execute_reply":"2021-05-31T10:03:59.263850Z"},"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-05-31T10:03:59.267805Z","iopub.execute_input":"2021-05-31T10:03:59.268174Z","iopub.status.idle":"2021-05-31T10:03:59.985018Z","shell.execute_reply.started":"2021-05-31T10:03:59.268139Z","shell.execute_reply":"2021-05-31T10:03:59.984055Z"},"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-05-31T10:03:59.986600Z","iopub.execute_input":"2021-05-31T10:03:59.986933Z","iopub.status.idle":"2021-05-31T10:04:00.019916Z","shell.execute_reply.started":"2021-05-31T10:03:59.986904Z","shell.execute_reply":"2021-05-31T10:04:00.019208Z"},"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-05-31T10:04:00.021334Z","iopub.execute_input":"2021-05-31T10:04:00.021600Z","iopub.status.idle":"2021-05-31T10:04:50.940112Z","shell.execute_reply.started":"2021-05-31T10:04:00.021574Z","shell.execute_reply":"2021-05-31T10:04:50.939013Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(history.history['lr'])","metadata":{"execution":{"iopub.status.busy":"2021-05-31T10:04:50.941344Z","iopub.execute_input":"2021-05-31T10:04:50.941632Z","iopub.status.idle":"2021-05-31T10:04:51.079702Z","shell.execute_reply.started":"2021-05-31T10:04:50.941604Z","shell.execute_reply":"2021-05-31T10:04:51.078582Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"# 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-05-31T10:04:51.081303Z","iopub.execute_input":"2021-05-31T10:04:51.081721Z","iopub.status.idle":"2021-05-31T10:04:51.113122Z","shell.execute_reply.started":"2021-05-31T10:04:51.081678Z","shell.execute_reply":"2021-05-31T10:04:51.112165Z"},"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-05-31T10:04:51.114339Z","iopub.execute_input":"2021-05-31T10:04:51.114624Z","iopub.status.idle":"2021-05-31T10:05:26.437570Z","shell.execute_reply.started":"2021-05-31T10:04:51.114597Z","shell.execute_reply":"2021-05-31T10:05:26.436543Z"},"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-05-31T10:05:26.440130Z","iopub.execute_input":"2021-05-31T10:05:26.440834Z","iopub.status.idle":"2021-05-31T10:05:27.085322Z","shell.execute_reply.started":"2021-05-31T10:05:26.440789Z","shell.execute_reply":"2021-05-31T10:05:27.084332Z"},"trusted":true},"execution_count":null,"outputs":[]}]}