{"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,"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","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 = .001\nbatch_size = 1024\nepochs = 100\nn_stop = 10\nn_fold = 5\nseed = 42","metadata":{"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,"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,"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,"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn = pd.read_csv(train_file)\nprint(trn.shape)\ntrn.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tst = pd.read_csv(test_file)\nprint(tst.shape)\ntst.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"sub = pd.read_csv(sample_file)\nprint(sub.shape)\nsub.head()","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn[pname_col].value_counts()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tst[pname_col].value_counts()","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"trn[phone_col].value_counts()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"tst[phone_col].value_counts()","metadata":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_raw = dfs['Raw']\nprint(df_raw.shape)\ndf_raw.head()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_raw.info()","metadata":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_raw['BiasUncertaintyNanos'].describe())\ndf_raw['BiasUncertaintyNanos'].hist(bins=20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_raw['ReceivedSvTimeUncertaintyNanos'].describe())\ndf_raw['ReceivedSvTimeUncertaintyNanos'].hist(bins=20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_raw.AccumulatedDeltaRangeUncertaintyMeters.describe())\ndf_raw.AccumulatedDeltaRangeUncertaintyMeters.hist(bins=20)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(df_raw.Cn0DbHz.describe())\ndf_raw.Cn0DbHz.hist(bins=20)","metadata":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived.info()","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"derived = derived.loc[derived.constellationType != 0]\nprint(derived.shape)","metadata":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_tst.describe()","metadata":{"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":{"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":{"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":{"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.Adam(lrate), loss='mean_squared_error')\n    return model","metadata":{"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,"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,"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":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"plt.plot(history.history['lr'])","metadata":{"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":{"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":{"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":{"trusted":true},"execution_count":null,"outputs":[]}]}