{"metadata":{"kaggle":{"accelerator":"gpu","dataSources":[{"sourceId":60095,"databundleVersionId":6542333,"sourceType":"competition"}],"dockerImageVersionId":30646,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":true},"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"codemirror_mode":{"name":"ipython","version":3},"file_extension":".py","mimetype":"text/x-python","name":"python","nbconvert_exporter":"python","pygments_lexer":"ipython3","version":"3.10.13"}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"## Import and set constants","metadata":{}},{"cell_type":"code","source":"import glob\nfrom dataclasses import dataclass\nimport numpy as np\nimport pandas as pd\nfrom tqdm.notebook import tqdm\nfrom scipy.interpolate import InterpolatedUnivariateSpline\n\nimport math\n\nimport torch\n\nimport torch\nfrom torch import nn, Tensor\nfrom torch.nn import TransformerEncoder, TransformerEncoderLayer\nfrom torch.utils.data import dataset\n\nINPUT_PATH = '/kaggle/input/smartphone-decimeter-2023'\n\nWGS84_SEMI_MAJOR_AXIS = 6378137.0\nWGS84_SEMI_MINOR_AXIS = 6356752.314245\nWGS84_SQUARED_FIRST_ECCENTRICITY  = 6.69437999013e-3\nWGS84_SQUARED_SECOND_ECCENTRICITY = 6.73949674226e-3\n\nHAVERSINE_RADIUS = 6_371_000","metadata":{"_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Load Device","metadata":{}},{"cell_type":"code","source":"device = torch.device(\"cuda\" if torch.cuda.is_available() else \"cpu\")","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"device","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Class to manage converting between lat/lng BLH and GNSS","metadata":{}},{"cell_type":"code","source":"@dataclass\nclass ECEF:\n    x: np.array\n    y: np.array\n    z: np.array\n\n    def to_numpy(self):\n        return np.stack([self.x, self.y, self.z], axis=0)\n\n    @staticmethod\n    def from_numpy(pos):\n        x, y, z = [np.squeeze(w) for w in np.split(pos, 3, axis=-1)]\n        return ECEF(x=x, y=y, z=z)\n\n@dataclass\nclass BLH:\n    lat : np.array\n    lng : np.array\n    hgt : np.array = 0\n        \n\ndef ECEF_to_BLH(ecef):\n    a = WGS84_SEMI_MAJOR_AXIS\n    b = WGS84_SEMI_MINOR_AXIS\n    e2  = WGS84_SQUARED_FIRST_ECCENTRICITY\n    e2_ = WGS84_SQUARED_SECOND_ECCENTRICITY\n    x = ecef.x\n    y = ecef.y\n    z = ecef.z\n    r = np.sqrt(x**2 + y**2)\n    t = np.arctan2(z * (a/b), r)\n    B = np.arctan2(z + (e2_*b)*np.sin(t)**3, r - (e2*a)*np.cos(t)**3)\n    L = np.arctan2(y, x)\n    n = a / np.sqrt(1 - e2*np.sin(B)**2)\n    H = (r / np.cos(B)) - n\n    return BLH(lat=B, lng=L, hgt=H)\n\ndef haversine_distance(blh_1, blh_2):\n    dlat = blh_2.lat - blh_1.lat\n    dlng = blh_2.lng - blh_1.lng\n    a = np.sin(dlat/2)**2 + np.cos(blh_1.lat) * np.cos(blh_2.lat) * np.sin(dlng/2)**2\n    dist = 2 * HAVERSINE_RADIUS * np.arcsin(np.sqrt(a))\n    return dist\n\ndef pandas_haversine_distance(df1, df2):\n    blh1 = BLH(\n        lat=np.deg2rad(df1['LatitudeDegrees'].to_numpy()),\n        lng=np.deg2rad(df1['LongitudeDegrees'].to_numpy()),\n        hgt=0,\n    )\n    blh2 = BLH(\n        lat=np.deg2rad(df2['LatitudeDegrees'].to_numpy()),\n        lng=np.deg2rad(df2['LongitudeDegrees'].to_numpy()),\n        hgt=0,\n    )\n    return haversine_distance(blh1, blh2)\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def ecef_to_lat_lng(tripID, gnss_df, UnixTimeMillis):\n    ecef_columns = ['WlsPositionXEcefMeters', 'WlsPositionYEcefMeters', 'WlsPositionZEcefMeters']\n    columns = ['utcTimeMillis'] + ecef_columns\n    ecef_df = (gnss_df.drop_duplicates(subset='utcTimeMillis')[columns]\n               .dropna().reset_index(drop=True))\n    ecef = ECEF.from_numpy(ecef_df[ecef_columns].to_numpy())\n    blh  = ECEF_to_BLH(ecef)\n\n    TIME = ecef_df['utcTimeMillis'].to_numpy()\n    lat = InterpolatedUnivariateSpline(TIME, blh.lat, ext=3)(UnixTimeMillis)\n    lng = InterpolatedUnivariateSpline(TIME, blh.lng, ext=3)(UnixTimeMillis)\n    return pd.DataFrame({\n#         'tripId' : tripID,\n        'utcTimeMillis'   : UnixTimeMillis,\n        'LatitudeDegrees'  : np.degrees(lat),\n        'LongitudeDegrees' : np.degrees(lng),\n    })\n\ndef calc_score(pred_df, gt_df):\n    d = pandas_haversine_distance(pred_df, gt_df)\n    score = np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])    \n    return score\n","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def print_comparison(lat, lng, gt_lat, gt_lng):\n    for lat_val, lng_val, gt_lat_val, gt_lng_val in zip(lat, lng, gt_lat, gt_lng):\n        print(f'Pred: ({lat_val:<12.7f}, {lng_val:<12.7f}) Ground Truth: ({gt_lat_val:<12.7f}, {gt_lng_val:<12.7f})')\n        \ndef print_batch(amnt, lat_arr, lng_arr, gt_lat_arr, gt_lng_arr):\n    for batch in range(amnt):\n        print(f'Val data {batch}')\n        print_comparison(lat_arr[batch], lng_arr[batch], gt_lat_arr[batch], gt_lng_arr[batch])","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Loading Data","metadata":{}},{"cell_type":"code","source":"%%capture --no-stdout\n\npred_dfs  = []\ngt_dfs = []\n\n\n\nfor dirname in sorted(glob.glob(f'/kaggle/input/smartphone-decimeter-2023/sdc2023/train/*/*')):\n    drive, phone = dirname.split('/')[-2:]\n    tripID  = f'{drive}/{phone}'\n    gnss_df = pd.read_csv(f'{dirname}/device_gnss.csv')\n    gt_df   = pd.read_csv(f'{dirname}/ground_truth.csv')\n    \n    info_cols = ['IonosphericDelayMeters', 'TroposphericDelayMeters']\n    columns = ['utcTimeMillis'] + info_cols\n    info_df = (gnss_df.drop_duplicates(subset='utcTimeMillis')[columns].fillna(0).reset_index(drop=True))\n    \n    for col in info_cols:\n        info_df[col] = info_df[col].fillna((info_df[col].bfill() + info_df[col].ffill()) / 2)\n        \n    pred_df = ecef_to_lat_lng(tripID, gnss_df, gt_df['UnixTimeMillis'])\n    pred_df = pd.merge(pred_df, info_df, on='utcTimeMillis', how='left')\n    gt_df   = gt_df[['LatitudeDegrees', 'LongitudeDegrees']]\n    print(tripID)\n#     print(pred_df.shape)\n#     print(gt_df.shape)\n    pred_dfs.append(pred_df)\n    gt_dfs.append(gt_df)","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define Baseline Model","metadata":{}},{"cell_type":"code","source":"class PositionalEncoding(torch.nn.Module):\n\n    def __init__(self, d_model: int, dropout: float = 0.1, max_len: int = 5000):\n        super().__init__()\n        self.dropout = torch.nn.Dropout(p=dropout)\n\n        position = torch.arange(max_len).unsqueeze(1)\n        div_term = torch.exp(torch.arange(0, d_model, 2) * (-math.log(10000.0) / d_model))\n        pe = torch.zeros(max_len, 1, d_model)\n        pe[:, 0, 0::2] = torch.sin(position * div_term)\n        pe[:, 0, 1::2] = torch.cos(position * div_term)\n        self.register_buffer('pe', pe)\n\n    def forward(self, x: torch.Tensor) -> torch.Tensor:\n        \"\"\"\n        Arguments:\n            x: Tensor, shape ``[seq_len, batch_size, embedding_dim]``\n        \"\"\"\n        x = x + self.pe[:x.size(0)]\n        return self.dropout(x)\n    \n# time sequence encoder\nclass TransformerEncoder(torch.nn.Module):\n    def __init__(self, config):\n        super().__init__()\n        self.config = config\n        self.upscale = torch.nn.Linear(config.input_dim, config.d_model) \n        self.pos_encoder = PositionalEncoding(config.d_model, config.dropout, config.max_seq_len)\n        self.transformer = torch.nn.TransformerEncoder(\n            torch.nn.TransformerEncoderLayer(\n                d_model=config.d_model,\n                nhead=config.nhead,\n                dim_feedforward=config.dim_feedforward,\n                dropout=config.dropout,\n                activation=config.activation,\n                batch_first=True\n            ),\n            num_layers=config.num_layers\n        )\n        self.fc = torch.nn.Sequential(\n            torch.nn.Linear(config.d_model, config.d_model),\n            torch.nn.ReLU(),\n            torch.nn.Linear(config.d_model, config.d_model),\n            torch.nn.ReLU(),\n            torch.nn.Linear(config.d_model, config.output_dim)\n        )\n        \n        \n    def forward(self, x):\n        x = self.upscale(x)\n        x = self.pos_encoder(x)\n        x = self.transformer(x)\n        x = self.fc(x)\n        return x","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Declare Model","metadata":{}},{"cell_type":"code","source":"class Config:\n    def __init__(self, config_dict):\n        self.input_dim = config_dict.get('input_dim')\n        self.d_model = config_dict.get('d_model')\n        self.nhead = config_dict.get('nhead')\n        self.dim_feedforward = config_dict.get('dim_feedforward')\n        self.dropout = config_dict.get('dropout')\n        self.activation = config_dict.get('activation')\n        self.num_layers = config_dict.get('num_layers')\n        self.output_dim = config_dict.get('output_dim')\n        self.input_path = config_dict.get('input_path')\n        self.max_seq_len = config_dict.get('max_seq_len')\n        self.val_split = config_dict.get('val_split')\n\nconfig_dict = {\n    \"input_dim\": 4, \n    \"d_model\": 8,\n    \"nhead\": 4,\n    \"dim_feedforward\": 128,\n    \"dropout\": 0.1,\n    \"activation\": \"relu\",\n    \"num_layers\": 6,\n    \"output_dim\": 2, \n    \"input_path\": \"INPUT_PATH\",\n    \"max_seq_len\": 3500,\n    \"val_split\": 0.05,\n}\n\nconfig = Config(config_dict)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"model = TransformerEncoder(config)\nmodel.to(device)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Define Dataset Class","metadata":{}},{"cell_type":"code","source":"\nclass GNSSDataset(torch.utils.data.Dataset):\n\n    # dataset constructor.\n    def __init__(self, pred_dfs, gt_dfs):\n        '''\n        pred_dfs (list): list of dataframes containing the chosen columns\n        gt_df (pd.DataFrame): dataframe containing the ground truth lat, lng\n        '''\n        self.pred_dfs = pred_dfs\n        self.gt_dfs = gt_dfs\n        self.sequences = []\n        self.labels = []\n        \n        for pred_df in self.pred_dfs:\n            x_np = pred_df[['LatitudeDegrees', 'LongitudeDegrees', 'IonosphericDelayMeters', 'TroposphericDelayMeters']].to_numpy()\n            ## pad to max sequence length\n            pad = np.zeros((config.max_seq_len - x_np.shape[0], x_np.shape[1]))\n            x_np = np.vstack([x_np, pad])\n            #x_np = x_np/180\n            self.sequences.append(x_np)\n            \n        for gt_df in self.gt_dfs:\n            y_np = gt_df[['LatitudeDegrees', 'LongitudeDegrees']].to_numpy()\n            ## pad to max sequence length\n            pad = np.zeros((config.max_seq_len - y_np.shape[0], y_np.shape[1]))\n            y_np = np.vstack([y_np, pad])\n            #y_np = y_np/180\n            self.labels.append(y_np)\n\n        self.sequences = np.array(self.sequences, dtype=np.float32)\n        self.labels = np.array(self.labels,  dtype=np.float32)\n        \n        self.labels = self.labels - self.sequences[:, :, :2] # just the residuals\n        \n        print(\"seq and label shapes\")\n        print(self.sequences.shape)\n        print(self.labels.shape)\n\n    # return an instance from the dataset\n    def __getitem__(self, i):\n        # return tuple of ith sequence and label\n        return self.sequences[i], self.labels[i]\n\n    # return the size of the dataset\n    def __len__(self):\n        return self.sequences.shape[0]","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Initialize Dataset","metadata":{}},{"cell_type":"code","source":"import sklearn\nfrom sklearn.model_selection import train_test_split","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print(len(pred_dfs)*config.val_split)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# train test split into training and validation\nX_train, X_val, y_train, y_val = train_test_split(pred_dfs, gt_dfs, \n                            test_size=config.val_split, random_state=2)\n\n\ntrain_dataset = GNSSDataset(X_train, y_train)\nval_dataset = GNSSDataset(X_val, y_val)\n\ntrain_loader = torch.utils.data.DataLoader(train_dataset, batch_size=32, shuffle=True)\nval_loader = torch.utils.data.DataLoader(val_dataset, batch_size=32, shuffle=True)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## Train the Model","metadata":{}},{"cell_type":"code","source":"epochs = 1000\nn_eval = 10","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Initalize optimizer (for gradient descent) and loss function\noptimizer = torch.optim.Adam(params = model.parameters(), lr = 0.00001)\nloss_fn = torch.nn.MSELoss()\n\nstep = 0\n\nPATH = \"model.pt\"\nbest_loss = 99999999 #high number\n\n# Regularization strength\nlambda_l1 = 0.1\nlambda_l2 = 0.25\n\nfor epoch in range(epochs):\n    print(f\"Epoch {epoch + 1} of {epochs}\")\n\n    # Loop over each batch in the dataset\n    for batch in tqdm(train_loader):\n        \n        optimizer.zero_grad()\n\n        # Unpack the data and labels\n        features, labels = batch\n        \n        features = features.to(device)\n        labels = labels.to(device)\n\n        # Forward propagate\n\n        outputs = model(features)\n\n        # Backpropagation and gradient descent            \n\n        loss = loss_fn(outputs, labels)\n        # L1 regularization\n        l1_reg = sum(param.abs().sum() for param in model.parameters())\n        # L2 regularization\n        l2_reg = sum(param.pow(2).sum() for param in model.parameters())\n\n        # Add the regularization terms to the loss\n        loss += lambda_l1 * l1_reg + lambda_l2 * l2_reg\n\n        loss.backward()\n        optimizer.step()\n        \n        print(f'Loss/train {loss}')\n\n        # Periodically evaluate our model + log to Tensorboard\n        if step % n_eval == 0:\n            # Compute training loss and metrics\n            \n            model.eval()\n            mean_dist = 0\n            mean_score = 0\n            mean_raw_score = 3.218092441558838 # not passed thru model, precomputed for this seed\n            count = 0\n            val_loss = 0\n            for features, labels in val_loader:\n#                 print(features.shape)\n#                 print(features[0, :20])\n                features = features.to(device)\n                labels = labels.to(device)\n                \n                val_pred = model(features)\n                \n                val_loss = loss_fn(val_pred, labels)\n                \n                features = features.detach().cpu()# * 180\n                val_pred = val_pred.detach().cpu()# * 180\n                labels = labels.detach().cpu()# * 180\n                print(val_pred.shape)\n                pred_lats = val_pred[:, :, 0] + features[:, :, 0]\n                pred_lngs = val_pred[:, :, 1] + features[:, :, 1]\n                gt_lats = labels[:, :, 0] + features[:, :, 0]\n                gt_lns = labels[:, :, 1] + features[:, :, 1]\n                \n                start = 60\n                cutoff = 65\n                # print some samples\n                print_batch(1, pred_lats[:, start:cutoff], pred_lngs[:, start:cutoff], gt_lats[:, start:cutoff], gt_lns[:, start:cutoff]) \n                \n                blh1 = BLH(np.deg2rad(pred_lats), np.deg2rad(pred_lngs), hgt = 0)\n                blh2 = BLH(np.deg2rad(gt_lats), np.deg2rad(gt_lns), hgt = 0)\n                d = haversine_distance(blh1, blh2)\n                mean_score += np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])    \n                mean_dist += d\n#                 blh1 = BLH(np.deg2rad(features[:, :, 0]), np.deg2rad(features[:, :, 1]), hgt = 0)\n#                 d = haversine_distance(blh1, blh2)\n#                 mean_raw_score += np.mean([np.quantile(d, 0.50), np.quantile(d, 0.95)])    \n                \n                count += 1\n                \n                features = features.cpu()\n                labels = labels.cpu()\n                \n            if(val_loss < best_loss):\n                best_loss = val_loss\n                \n                torch.save({\n                            'epoch': epoch,\n                            'model_state_dict': model.state_dict(),\n                            'optimizer_state_dict': optimizer.state_dict(),\n                            'loss': val_loss,\n                            }, PATH)\n\n                \n#             print(f'Val mean dist {mean_dist.mean()}')\n            print(f'Val mean score {mean_score} vs baseline {mean_raw_score}')\n            print(f'Loss/val {val_loss}')\n\n            #turn on training, evaluate turns off training\n            model.train()\n\n        step += 1\n","metadata":{"scrolled":true,"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"","metadata":{}}]}