{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.7.10","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceId":26933,"databundleVersionId":2264643,"sourceType":"competition"}],"dockerImageVersionId":30096,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"markdown","source":"# Predict Next Point with the IMU Data\nAuthor: **Alvin.ai**<br>\nIMU is Inertial Measurement Unit, whicn involves accelerometer, gyroscope asn so on. In Microsoft Research's Indoor Location & Navigation Competition, many kagglers used the IMU Data. In recent days, I attempted to build a model to utilize the given sensors' data and I found that it works for me as below:<br>\n1. Reject Outlier + KF Smooth + Phone Mean: LB-**5.653**\n2. **IMU Prediction** + Reject Outlier + KF Smooth + Phone Mean: LB-**5.476**\n\nWhat I do is to simple, **use historical points with sensors' dataset to predict where the next point is.**<br>\nTalk is cheap, show you the code. Hope it is helpful to you.<br>\nBy the way, **I am looking forward to find a teammates, if you're interested, pls see the end of this kernal.**<br>\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19"}},{"cell_type":"code","source":"import math\nimport numpy as np\nfrom cv2 import Rodrigues\nfrom math import sin, cos, atan2, sqrt\nimport pandas as pd\nfrom pathlib import Path\nimport pyproj\nfrom pyproj import Proj, transform\nimport matplotlib.pyplot as plt\nimport plotly.express as px\nfrom sklearn.preprocessing import LabelEncoder\nfrom sklearn.model_selection import KFold, TimeSeriesSplit\nfrom sklearn.metrics import accuracy_score\nimport lightgbm as lgb\nfrom tqdm import tqdm\nimport warnings\nwarnings.filterwarnings(\"ignore\", category=Warning)","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:40.894042Z","iopub.execute_input":"2024-05-19T13:25:40.894567Z","iopub.status.idle":"2024-05-19T13:25:44.022396Z","shell.execute_reply.started":"2024-05-19T13:25:40.894464Z","shell.execute_reply":"2024-05-19T13:25:44.021160Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"data_dir = Path(\"../input/google-smartphone-decimeter-challenge\")\nbl_trn_fname = 'baseline_locations_train.csv'\nbl_tst_fname = 'baseline_locations_test.csv'\nsample_fname = 'sample_submission.csv'","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:44.024309Z","iopub.execute_input":"2024-05-19T13:25:44.024790Z","iopub.status.idle":"2024-05-19T13:25:44.030232Z","shell.execute_reply.started":"2024-05-19T13:25:44.024750Z","shell.execute_reply":"2024-05-19T13:25:44.028898Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"bl_trn_df = pd.read_csv(data_dir / bl_trn_fname)\nbl_tst_df = pd.read_csv(data_dir / bl_tst_fname)\nsample_df = pd.read_csv(data_dir / sample_fname)","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:44.032321Z","iopub.execute_input":"2024-05-19T13:25:44.032690Z","iopub.status.idle":"2024-05-19T13:25:45.007250Z","shell.execute_reply.started":"2024-05-19T13:25:44.032618Z","shell.execute_reply":"2024-05-19T13:25:45.006057Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"print('Baseline Train shape:', bl_trn_df.shape)\nprint('Baseline Test shape:', bl_tst_df.shape)\nprint('Test shape:', sample_df.shape)","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.009515Z","iopub.execute_input":"2024-05-19T13:25:45.010071Z","iopub.status.idle":"2024-05-19T13:25:45.018715Z","shell.execute_reply.started":"2024-05-19T13:25:45.010013Z","shell.execute_reply":"2024-05-19T13:25:45.017404Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 1. Euler Angles to Rotation Vector\nEuler Angles <-> Rotation Matrix <-> Rotation Vector<br>\nMore info.:[About IMU: OrientationDeg to Rotation Vector](https://www.kaggle.com/c/google-smartphone-decimeter-challenge/discussion/247834)","metadata":{"execution":{"iopub.status.busy":"2021-07-03T08:49:08.128564Z","iopub.execute_input":"2021-07-03T08:49:08.129017Z","iopub.status.idle":"2021-07-03T08:49:08.134112Z","shell.execute_reply.started":"2021-07-03T08:49:08.128971Z","shell.execute_reply":"2021-07-03T08:49:08.132931Z"}}},{"cell_type":"code","source":"# pitch:y\n# yaw:z\n# roll:x\ndef an2v(y_delta, z_delta, x_delta):\n    '''\n    Euler Angles ->Rotation Matrix -> Rotation Vector\n\n    Input：\n        1. y_delta          (float): the angle with rotateing around y-axis.\n        2. z_delta         (float): the angle with rotateing around z-axis. \n        3. x_delta         (float): the angle with rotateing around x-axis. \n    Output：\n        rx/ry/rz             (float): the rotation vector with rotateing \n    \n    Code Ref.: https://www.zacobria.com/universal-robots-knowledge-base-tech-support-forum-hints-tips/python-code-example-of-converting-rpyeuler-angles-to-rotation-vectorangle-axis-for-universal-robots/\n    (Note：In Code Ref: pitch=y,yaw=z,roll=x. But Google is pitch=x,yaw=z,roll=y)\n    '''\n    # yaw: z\n    Rz_Matrix = np.matrix([\n    [math.cos(z_delta), -math.sin(z_delta), 0],\n    [math.sin(z_delta), math.cos(z_delta), 0],\n    [0, 0, 1]\n    ])\n    \n    # pitch: y\n    Ry_Matrix = np.matrix([\n    [math.cos(y_delta), 0, math.sin(y_delta)],\n    [0, 1, 0],\n    [-math.sin(y_delta), 0, math.cos(y_delta)]\n    ])\n    \n    # roll: x\n    Rx_Matrix = np.matrix([\n    [1, 0, 0],\n    [0, math.cos(x_delta), -math.sin(x_delta)],\n    [0, math.sin(x_delta), math.cos(x_delta)]\n    ])\n\n    R = Rz_Matrix * Ry_Matrix * Rx_Matrix\n\n    theta = math.acos(((R[0, 0] + R[1, 1] + R[2, 2]) - 1) / 2)\n    multi = 1 / (2 * math.sin(theta))\n\n    rx = multi * (R[2, 1] - R[1, 2]) * theta\n    ry = multi * (R[0, 2] - R[2, 0]) * theta\n    rz = multi * (R[1, 0] - R[0, 1]) * theta\n\n    return rx, ry, rz","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.020381Z","iopub.execute_input":"2024-05-19T13:25:45.020790Z","iopub.status.idle":"2024-05-19T13:25:45.035071Z","shell.execute_reply.started":"2024-05-19T13:25:45.020753Z","shell.execute_reply":"2024-05-19T13:25:45.034039Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def v2a(rotation_v):\n    '''\n    Rotation Vector -> Rotation Matrix -> Euler Angles\n\n    Input：\n        rx/ry/rz             (float): the rotation vector with rotateing around x/y/z-axis.\n    Output：\n        1. y_delta          (float): the angle with rotateing around y-axis.\n        2. z_delta         (float): the angle with rotateing around z-axis. \n        3. x_delta         (float): the angle with rotateing around x-axis.  \n    '''\n    # Rotation Vector -> Rotation Matrix\n    R = Rodrigues(rotation_v)[0]\n\n    sq = sqrt(R[2,1] ** 2 +  R[2,2] ** 2)\n\n    if  not (sq < 1e-6) :\n        x_delta = atan2(R[2,1] , R[2,2])\n        y_delta = atan2(-R[2,0], sq)\n        z_delta = atan2(R[1,0], R[0,0])\n    else :\n        x_delta = atan2(-R[1,2], R[1,1])\n        y_delta = atan2(-R[2,0], sq)\n        z_delta = 0\n\n    return y_delta, z_delta, x_delta","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.036390Z","iopub.execute_input":"2024-05-19T13:25:45.036876Z","iopub.status.idle":"2024-05-19T13:25:45.052824Z","shell.execute_reply.started":"2024-05-19T13:25:45.036841Z","shell.execute_reply":"2024-05-19T13:25:45.051756Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 2. Prepare IMU Dataset\nThis part is to prepare the dataset for the model. I divided this part into the following steps:<br>\n(1) **Load GNSS Log**<br>\n(2) **Merge sub-dataset** (Status/UncalAccel/UncalGyro/UncalMag/OrientationDeg)<br>\n(3) **UTC to GpsEpoch**<br>\n(4) **OrientationDeg to Rotation Vecto**r<br>\n(5) **Calibrate Sensors' data**<br>\n(6) **LatDeg&lngDeg to x/y/z**<br>\n(7) **Orgainze Data** (eg. t1 t2 t3 t4 t5 -> t6)<br>\n(8) **Clean Data** (unrelated-aixs features and uncalibrated features)<br>\n(9) **Add Statistic Features**","metadata":{"execution":{"iopub.status.busy":"2021-07-03T09:01:39.921887Z","iopub.execute_input":"2021-07-03T09:01:39.922246Z","iopub.status.idle":"2021-07-03T09:01:39.930463Z","shell.execute_reply.started":"2021-07-03T09:01:39.922217Z","shell.execute_reply":"2021-07-03T09:01:39.92891Z"}}},{"cell_type":"code","source":"def gnss_log_to_dataframes(path):\n    '''Load GNSS Log'''\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":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.054288Z","iopub.execute_input":"2024-05-19T13:25:45.055126Z","iopub.status.idle":"2024-05-19T13:25:45.069641Z","shell.execute_reply.started":"2024-05-19T13:25:45.055083Z","shell.execute_reply":"2024-05-19T13:25:45.068686Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def UTC2GpsEpoch(df):\n    '''UTC to GpsEpoch\n    \n    utcTimeMillis         : UTC epoch (1970/1/1)\n    millisSinceGpsEpoch   : GPS epoch(1980/1/6 midnight 12:00 UTC)\n    \n    Ref: https://www.kaggle.com/c/google-smartphone-decimeter-challenge/discussion/239187\n    '''\n    dt_offset = pd.to_datetime('1980-01-06 00:00:00') \n    dt_offset_in_ms = int(dt_offset.value / 1e6)\n    df['millisSinceGpsEpoch'] = df['utcTimeMillis'] - dt_offset_in_ms + 18000\n    return df","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.072354Z","iopub.execute_input":"2024-05-19T13:25:45.072851Z","iopub.status.idle":"2024-05-19T13:25:45.079826Z","shell.execute_reply.started":"2024-05-19T13:25:45.072814Z","shell.execute_reply":"2024-05-19T13:25:45.078940Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def prepare_imu_data(data_dir, dataset_name, cname, pname, bl_df):\n    '''Prepare IMU Dataset (For Train: IMU+GT+BL; For Test: IMU+BL)\n    Input：\n        1. data_dir: data_dir\n        2. dataset_name: dataset name（'train'/'test'）\n        3. cname: CollectionName\n        4. pname: phoneName\n        5. bl_df: baseline's dataframe\n    Output：df_all\n    '''\n    # load GNSS log\n    gnss_df = gnss_log_to_dataframes(str(data_dir / dataset_name / cname / pname / f'{pname}_GnssLog.txt'))\n    print('sub-dataset shape：')\n    print('Raw:', gnss_df['Raw'].shape)\n    print('Status:', gnss_df['Status'].shape)\n    print('UncalAccel:', gnss_df['UncalAccel'].shape)\n    print('UncalGyro:', gnss_df['UncalGyro'].shape)\n    print('UncalMag:', gnss_df['UncalMag'].shape)\n    print('OrientationDeg:', gnss_df['OrientationDeg'].shape)\n    print('Fix:', gnss_df['Fix'].shape)\n\n    # merge sub-datasets\n    # accel + gyro\n    imu_df = pd.merge_asof(gnss_df['UncalAccel'].sort_values('utcTimeMillis'),\n                           gnss_df['UncalGyro'].drop('elapsedRealtimeNanos', axis=1).sort_values('utcTimeMillis'),\n                           on = 'utcTimeMillis',\n                           direction='nearest')\n    # (accel + gyro) + mag\n    imu_df = pd.merge_asof(imu_df.sort_values('utcTimeMillis'),\n                           gnss_df['UncalMag'].drop('elapsedRealtimeNanos', axis=1).sort_values('utcTimeMillis'),\n                           on = 'utcTimeMillis',\n                           direction='nearest')\n    # ((accel + gyro) + mag) + OrientationDeg\n    imu_df = pd.merge_asof(imu_df.sort_values('utcTimeMillis'),\n                           gnss_df['OrientationDeg'].drop('elapsedRealtimeNanos', axis=1).sort_values('utcTimeMillis'),\n                           on = 'utcTimeMillis',\n                           direction='nearest')\n   \n    # UTC->GpsEpoch\n    imu_df = UTC2GpsEpoch(imu_df)\n\n    # print IMU time\n    dt_offset = pd.to_datetime('1980-01-06 00:00:00')\n    dt_offset_in_ms = int(dt_offset.value / 1e6)\n    tmp_datetime = pd.to_datetime(imu_df['millisSinceGpsEpoch'] + dt_offset_in_ms, unit='ms')\n    print(f\"imu_df time scope: {tmp_datetime.min()} - {tmp_datetime.max()}\")\n\n\n    if dataset_name == 'train':\n        # read GT dataset\n        gt_path = data_dir / dataset_name / cname / pname / 'ground_truth.csv'\n        gt_df = pd.read_csv(gt_path, usecols = ['collectionName', 'phoneName', 'millisSinceGpsEpoch', 'latDeg', 'lngDeg'])\n\n        # print GT time\n        tmp_datetime = pd.to_datetime(gt_df['millisSinceGpsEpoch'] + dt_offset_in_ms, unit='ms')\n        print(f\"gt_df time scope: {tmp_datetime.min()} - {tmp_datetime.max()}\")\n\n        # merge GT dataset\n        imu_df = pd.merge_asof(gt_df.sort_values('millisSinceGpsEpoch'),\n                               imu_df.drop(['elapsedRealtimeNanos'], axis=1).sort_values('millisSinceGpsEpoch'),\n                               on = 'millisSinceGpsEpoch',\n                               direction='nearest')\n    elif dataset_name == 'test':\n        # merge smaple_df\n        imu_df = pd.merge_asof(sample_df.sort_values('millisSinceGpsEpoch'),\n                           imu_df.drop(['elapsedRealtimeNanos'], axis=1).sort_values('millisSinceGpsEpoch'),\n                           on = 'millisSinceGpsEpoch',\n                           direction='nearest')\n\n    # OrientationDeg -> Rotation Vector\n    rxs = []\n    rys = []\n    rzs = []\n    for i in range(len(imu_df)):\n        y_delta = imu_df['rollDeg'].iloc[i]\n        z_delta = imu_df['yawDeg'].iloc[i]\n        x_delta = imu_df['pitchDeg'].iloc[i]\n        rx, ry, rz = an2v(y_delta, z_delta, x_delta)\n        rxs.append(rx)\n        rys.append(ry)\n        rzs.append(rz)\n\n    imu_df['ahrsX'] = rxs\n    imu_df['ahrsY'] = rys\n    imu_df['ahrsZ'] = rzs\n\n    # calibrate sensors' reading\n    for axis in ['X', 'Y', 'Z']:\n        imu_df['Accel{}Mps2'.format(axis)] = imu_df['UncalAccel{}Mps2'.format(axis)] - imu_df['Bias{}Mps2'.format(axis)]\n        imu_df['Gyro{}RadPerSec'.format(axis)] = imu_df['UncalGyro{}RadPerSec'.format(axis)] - imu_df['Drift{}RadPerSec'.format(axis)]\n        imu_df['Mag{}MicroT'.format(axis)] = imu_df['UncalMag{}MicroT'.format(axis)] - imu_df['Bias{}MicroT'.format(axis)]\n\n        # clearn bias features\n        imu_df.drop(['Bias{}Mps2'.format(axis), 'Drift{}RadPerSec'.format(axis), 'Bias{}MicroT'.format(axis)], axis = 1, inplace = True) \n\n    if dataset_name == 'train':\n        # merge Baseline dataset：imu_df + bl_df = (GT + IMU) + Baseline\n        df_all = pd.merge(imu_df.rename(columns={'latDeg':'latDeg_gt', 'lngDeg':'lngDeg_gt'}),\n                      bl_df.drop(['phone'], axis=1).rename(columns={'latDeg':'latDeg_bl','lngDeg':'lngDeg_bl'}),\n                      on = ['collectionName', 'phoneName', 'millisSinceGpsEpoch'])\n    elif dataset_name == 'test':\n        df_all = pd.merge(imu_df,\n              bl_df[(bl_df['collectionName']==cname) & (bl_df['phoneName']==pname)].drop(['phone'], axis=1).rename(columns={'latDeg':'latDeg_bl','lngDeg':'lngDeg_bl'}),\n              on = ['millisSinceGpsEpoch'])\n        df_all.drop(['phone'], axis=1, inplace=True)\n        \n    return df_all","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.081636Z","iopub.execute_input":"2024-05-19T13:25:45.082137Z","iopub.status.idle":"2024-05-19T13:25:45.111815Z","shell.execute_reply.started":"2024-05-19T13:25:45.082100Z","shell.execute_reply":"2024-05-19T13:25:45.110737Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def WGS84_to_ECEF(lat, lon, alt):\n    # convert to radians\n    rad_lat = lat * (np.pi / 180.0)\n    rad_lon = lon * (np.pi / 180.0)\n    a    = 6378137.0\n    # f is the flattening factor\n    finv = 298.257223563\n    f = 1 / finv   \n    # e is the eccentricity\n    e2 = 1 - (1 - f) * (1 - f)    \n    # N is the radius of curvature in the prime vertical\n    N = a / np.sqrt(1 - e2 * np.sin(rad_lat) * np.sin(rad_lat))\n    x = (N + alt) * np.cos(rad_lat) * np.cos(rad_lon)\n    y = (N + alt) * np.cos(rad_lat) * np.sin(rad_lon)\n    z = (N * (1 - e2) + alt)        * np.sin(rad_lat)\n    return x, y, z\n\ntransformer = pyproj.Transformer.from_crs(\n    {\"proj\":'geocent', \"ellps\":'WGS84', \"datum\":'WGS84'},\n    {\"proj\":'latlong', \"ellps\":'WGS84', \"datum\":'WGS84'},)\ndef ECEF_to_WGS84(x,y,z):\n    lon, lat, alt = transformer.transform(x,y,z,radians=False)\n    return lon, lat, alt","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.113245Z","iopub.execute_input":"2024-05-19T13:25:45.113623Z","iopub.status.idle":"2024-05-19T13:25:45.184506Z","shell.execute_reply.started":"2024-05-19T13:25:45.113589Z","shell.execute_reply":"2024-05-19T13:25:45.183086Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def get_xyz(df_all, dataset_name):\n    # baseline: lat/lngDeg -> x/y/z\n    df_all['Xbl'], df_all['Ybl'], df_all['Zbl'] = zip(*df_all.apply(lambda x: WGS84_to_ECEF(x.latDeg_bl, x.lngDeg_bl, x.heightAboveWgs84EllipsoidM), axis=1))\n    \n    if dataset_name == 'train':\n        # gt: lat/lngDeg -> x/y/z\n        df_all['Xgt'], df_all['Ygt'], df_all['Zgt'] = zip(*df_all.apply(lambda x: WGS84_to_ECEF(x.latDeg_gt, x.lngDeg_gt, x.heightAboveWgs84EllipsoidM), axis=1))\n        # copy lat/lngDeg\n        lat_lng_df = df_all[['latDeg_gt','lngDeg_gt', 'latDeg_bl', 'lngDeg_bl']]\n        df_all.drop(['latDeg_gt','lngDeg_gt', 'latDeg_bl', 'lngDeg_bl'], axis = 1, inplace = True)\n    elif dataset_name == 'test':\n        # copy lat/lngDeg\n        lat_lng_df = df_all[['latDeg_bl', 'lngDeg_bl']]\n        df_all.drop(['latDeg_bl', 'lngDeg_bl', 'latDeg','lngDeg',], axis = 1, inplace = True)     \n        \n    return lat_lng_df, df_all","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.186329Z","iopub.execute_input":"2024-05-19T13:25:45.186783Z","iopub.status.idle":"2024-05-19T13:25:45.200482Z","shell.execute_reply.started":"2024-05-19T13:25:45.186712Z","shell.execute_reply":"2024-05-19T13:25:45.199080Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def prepare_df_train(df_all_train, window_size):\n    '''prepare training dataset with all aixses'''\n    tgt_df = df_all_train.copy()\n    total_len = len(tgt_df) \n    moving_times = total_len - window_size \n    \n    tgt_df.rename(columns = {'yawDeg':'yawZDeg', 'rollDeg':'rollYDeg', 'pitchDeg':'pitchXDeg'}, inplace = True)\n\n    feature_cols = [f for f in list(tgt_df) if f not in ['Xgt', 'Ygt', 'Zgt']]\n\n    # Historical Feature names\n    hist_feats = []\n    for time_flag in range(1, window_size + 1):\n        for fn in feature_cols:\n            hist_feats.append(fn + '_' + str(time_flag))\n\n    # Window Sliding\n    # t1 t2 t3 t4 t5 -> t6\n    # t2 t3 t4 t5 t6 -> t7\n\n    # Add historical data \n    df_train = pd.DataFrame()\n    features = []\n    xs = []\n    ys = []\n    zs = []\n\n    for start_idx in range(moving_times):\n        feature_list = list()\n        x_list = list()\n        y_list = list()\n        z_list = list()\n\n        for window_idx in range(window_size):\n            feature_list.extend(tgt_df[feature_cols].iloc[start_idx + window_idx,:].to_list())\n        x_list.append(tgt_df['Xgt'].iloc[start_idx + window_size])\n        y_list.append(tgt_df['Ygt'].iloc[start_idx + window_size])\n        z_list.append(tgt_df['Zgt'].iloc[start_idx + window_size])\n\n        features.append(feature_list)\n        xs.extend(x_list)\n        ys.extend(y_list)\n        zs.extend(z_list)\n\n    df_train = pd.DataFrame(features, columns = hist_feats)\n    df_train['Xgt'] = xs\n    df_train['Ygt'] = ys\n    df_train['Zgt'] = zs\n    \n    # clean single-value feature: collectionName_[1-5]\\phoneName_[1-5]\n    tmp_feats = []\n    for fn in list(df_train):\n        if (fn.startswith('collectionName_') == False) and (fn.startswith('phoneName_') == False):\n            tmp_feats.append(fn)\n    df_train = df_train[tmp_feats]\n\n    # clean time feature\n    tmp_drop_feats = []\n    for f in list(df_train):\n        if (f.startswith('millisSinceGpsEpoch') == True) or (f.startswith('timeSinceFirstFixSeconds') == True) or (f.startswith('utcTimeMillis') == True):\n            tmp_drop_feats.append(f)\n    df_train.drop(tmp_drop_feats, axis = 1, inplace = True)\n    \n    return df_train","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.201873Z","iopub.execute_input":"2024-05-19T13:25:45.202374Z","iopub.status.idle":"2024-05-19T13:25:45.219998Z","shell.execute_reply.started":"2024-05-19T13:25:45.202335Z","shell.execute_reply":"2024-05-19T13:25:45.219015Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def prepare_df_test(df_all_test, window_size):\n    '''prepare testing dataset with all aixses'''\n    tgt_df = df_all_test.copy()\n    total_len = len(tgt_df) \n    moving_times = total_len - window_size \n    \n    tgt_df.rename(columns = {'yawDeg':'yawZDeg', 'rollDeg':'rollYDeg', 'pitchDeg':'pitchXDeg'}, inplace = True)\n\n    feature_cols = [f for f in list(tgt_df) if f not in ['Xgt', 'Ygt', 'Zgt']] \n    \n    hist_feats = []\n    for time_flag in range(1, window_size + 1):\n        for fn in feature_cols:\n            hist_feats.append(fn + '_' + str(time_flag))\n\n    # t1 t2 t3 t4 t5 -> t6\n    # t2 t3 t4 t5 t6 -> t7\n    df_test = pd.DataFrame()\n    features = []\n\n    for start_idx in range(moving_times):\n        feature_list = list()\n\n        for window_idx in range(window_size):\n            feature_list.extend(tgt_df[feature_cols].iloc[start_idx + window_idx,:].to_list())\n        features.append(feature_list)\n\n    df_test = pd.DataFrame(features, columns = hist_feats)\n\n    tmp_feats = []\n    for fn in list(df_test):\n        if (fn.startswith('collectionName_') == False) and (fn.startswith('phoneName_') == False):\n            tmp_feats.append(fn)\n    df_test = df_test[tmp_feats]\n\n    tmp_drop_feats = []\n    for f in list(df_test):\n        if (f.startswith('millisSinceGpsEpoch') == True) or (f.startswith('timeSinceFirstFixSeconds') == True) or (f.startswith('utcTimeMillis') == True) or (f.startswith('elapsedRealtimeNanos') == True):\n            tmp_drop_feats.append(f)\n    df_test.drop(tmp_drop_feats, axis = 1, inplace = True)\n    \n    return df_test","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.221807Z","iopub.execute_input":"2024-05-19T13:25:45.222630Z","iopub.status.idle":"2024-05-19T13:25:45.238331Z","shell.execute_reply.started":"2024-05-19T13:25:45.222578Z","shell.execute_reply":"2024-05-19T13:25:45.237360Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def remove_other_axis_feats(df_all, tgt_axis):\n    '''unrelated-aixs features and uncalibrated features'''\n    # Clean unrelated-aixs features\n    all_imu_feats = ['UncalAccelXMps2', 'UncalAccelYMps2', 'UncalAccelZMps2',\n                     'UncalGyroXRadPerSec', 'UncalGyroYRadPerSec', 'UncalGyroZRadPerSec',\n                     'UncalMagXMicroT', 'UncalMagYMicroT', 'UncalMagZMicroT',\n                     'ahrsX', 'ahrsY', 'ahrsZ',\n                     'AccelXMps2', 'AccelYMps2', 'AccelZMps2',\n                     'GyroXRadPerSec', 'GyroZRadPerSec', 'GyroYRadPerSec',\n                     'MagXMicroT', 'MagYMicroT', 'MagZMicroT',\n                     'yawZDeg', 'rollYDeg', 'pitchXDeg',\n                     'Xbl', 'Ybl', 'Zbl']\n    tgt_imu_feats = []\n    for axis in ['X', 'Y', 'Z']:\n        if axis != tgt_axis:\n            for f in all_imu_feats:\n                if f.find(axis) >= 0:\n                    tgt_imu_feats.append(f)\n            \n    tmp_drop_feats = []\n    for f in list(df_all):\n        if f.split('_')[0] in tgt_imu_feats:\n            tmp_drop_feats.append(f)\n\n    tgt_df = df_all.drop(tmp_drop_feats, axis = 1)\n    \n    # Clean uncalibrated features\n    uncal_feats = [f for f in list(tgt_df) if f.startswith('Uncal') == True]\n    tgt_df = tgt_df.drop(uncal_feats, axis = 1)\n    \n    return tgt_df","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.240135Z","iopub.execute_input":"2024-05-19T13:25:45.240950Z","iopub.status.idle":"2024-05-19T13:25:45.252945Z","shell.execute_reply.started":"2024-05-19T13:25:45.240900Z","shell.execute_reply":"2024-05-19T13:25:45.251775Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def add_stat_feats(data, tgt_axis):\n    for f in ['yawZDeg', 'rollYDeg', 'pitchXDeg']:\n        if f.find(tgt_axis) >= 0:\n            ori_feat = f\n            break\n            \n    cont_feats = ['heightAboveWgs84EllipsoidM', 'ahrs{}'.format(tgt_axis),\n           'Accel{}Mps2'.format(tgt_axis), 'Gyro{}RadPerSec'.format(tgt_axis), 'Mag{}MicroT'.format(tgt_axis),\n            '{}bl'.format(tgt_axis)] + [ori_feat]\n    \n    for f in cont_feats:\n        data[f + '_' + str(window_size) + '_mean'] = data[[f + f'_{i}' for i in range(1,window_size)]].mean(axis=1)\n        data[f + '_' + str(window_size) + '_std'] = data[[f + f'_{i}' for i in range(1,window_size)]].std(axis=1)\n        data[f + '_' + str(window_size) + '_max'] = data[[f + f'_{i}' for i in range(1,window_size)]].max(axis=1)\n        data[f + '_' + str(window_size) + '_min'] = data[[f + f'_{i}' for i in range(1,window_size)]].min(axis=1)\n        data[f + '_' + str(window_size) + '_median'] = data[[f + f'_{i}' for i in range(1,window_size)]].median(axis=1)\n    return data","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.254715Z","iopub.execute_input":"2024-05-19T13:25:45.255209Z","iopub.status.idle":"2024-05-19T13:25:45.269357Z","shell.execute_reply.started":"2024-05-19T13:25:45.255158Z","shell.execute_reply":"2024-05-19T13:25:45.267750Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 3. Modeling\nNote: I only use the given axis features for predict the target axis location.<br>\nFor example, use features contains x-axis to predict the next x location.<br>\nMore, I used LGBM here.","metadata":{"execution":{"iopub.status.busy":"2021-07-03T09:14:32.286714Z","iopub.execute_input":"2021-07-03T09:14:32.287331Z","iopub.status.idle":"2021-07-03T09:14:32.293464Z","shell.execute_reply.started":"2021-07-03T09:14:32.287291Z","shell.execute_reply":"2021-07-03T09:14:32.292014Z"}}},{"cell_type":"code","source":"# LightGBM\nparams = {\n    'metric':'mse',\n    'objective':'regression',\n    'seed':2021,\n    'boosting_type':'gbdt',\n    'early_stopping_rounds':10,\n    'subsample':0.7,\n    'feature_fraction':0.7,\n    'bagging_fraction': 0.7,\n    'reg_lambda': 10\n}\nwindow_size = 30\nverbose_flag = True\nfolds = 5","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.270931Z","iopub.execute_input":"2024-05-19T13:25:45.271296Z","iopub.status.idle":"2024-05-19T13:25:45.284067Z","shell.execute_reply.started":"2024-05-19T13:25:45.271239Z","shell.execute_reply":"2024-05-19T13:25:45.283178Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Example: I use SJC's dataset for training \ntgt_cns = ['2021-04-22-US-SJC-1', '2021-04-28-US-SJC-1', '2021-04-29-US-SJC-2']\ncn2pn_df = bl_trn_df[['collectionName', 'phoneName']].drop_duplicates()\n\ndf_trains = []\nlat_lng_df_trains = []\nfor tgt_cn in tqdm(tgt_cns):\n    pns = list(cn2pn_df[cn2pn_df['collectionName'] == tgt_cn]['phoneName'].values)\n    for tgt_pn in pns:\n        print('Prepare Training Dataset：', tgt_cn + '_' + tgt_pn)  \n        df_all_train = prepare_imu_data(data_dir, 'train', tgt_cn, tgt_pn, bl_trn_df)\n        lat_lng_df_train, df_all_train = get_xyz(df_all_train, 'train')\n        df_train = prepare_df_train(df_all_train,  window_size) # 所有轴的数据\n        df_trains.append(df_train)\n        lat_lng_df_trains.append(lat_lng_df_train)\n        print('_'*20)\n        \ndf_train = pd.concat(df_trains, axis = 0)\nlat_lng_df_train = pd.concat(lat_lng_df_trains, axis = 0)\nprint('Final Dataset shape：', df_train.shape)","metadata":{"execution":{"iopub.status.busy":"2024-05-19T13:25:45.285695Z","iopub.execute_input":"2024-05-19T13:25:45.286097Z"},"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Example: I choose one of SJC collection from the test dataset as my test dataset, you can choose what as you like\ncname_test = '2021-04-29-US-SJC-3'\npname_test = 'SamsungS20Ultra'\ndf_all_test = prepare_imu_data(data_dir, 'test', cname_test, pname_test, bl_tst_df)\nlat_lng_df_test, df_all_test = get_xyz(df_all_test, 'test')\ndf_test = prepare_df_test(df_all_test,  window_size)\nprint('df_test:', df_test.shape)\nprint('df_test.columns:', df_test.columns)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def training(df_train, df_test, tgt_axis, window_size):\n    '''For the given axis target to train the model. Also, it has validation and prediciton.'''\n    df_train = remove_other_axis_feats(df_train, tgt_axis)\n    df_train = add_stat_feats(df_train, tgt_axis)\n    df_test = remove_other_axis_feats(df_test, tgt_axis)\n    df_test = add_stat_feats(df_test, tgt_axis)\n    \n    feature_names = [f for f in list(df_train) if f not in ['Xgt', 'Ygt', 'Zgt']]\n    target = '{}gt'.format(tgt_axis)\n\n    kfold = KFold(n_splits=folds, shuffle=True, random_state=params['seed'])\n\n    pred_valid = np.zeros((len(df_train),)) \n    pred_test = np.zeros((len(df_test),)) \n    scores = []\n    for fold_id, (trn_idx, val_idx) in enumerate(kfold.split(df_train, df_train[target])):\n        X_train = df_train.iloc[trn_idx][feature_names]\n        Y_train = df_train.iloc[trn_idx][target]\n        X_val = df_train.iloc[val_idx][feature_names]\n        Y_val = df_train.iloc[val_idx][target]\n\n        model = lgb.LGBMRegressor(**params)\n        lgb_model = model.fit(X_train, \n                              Y_train,\n                              eval_names=['train', 'valid'],\n                              eval_set=[(X_train, Y_train), (X_val, Y_val)],\n                              verbose=0,\n                              eval_metric=params['metric'],\n                              early_stopping_rounds=params['early_stopping_rounds'])\n\n        pred_valid[val_idx] = lgb_model.predict(X_val, num_iteration =  lgb_model.best_iteration_)\n        pred_test += lgb_model.predict(df_test[feature_names], num_iteration =  lgb_model.best_iteration_)\n\n        scores.append(lgb_model.best_score_['valid']['l2'])\n    \n    pred_test = pred_test /  kfold.n_splits\n    \n    if verbose_flag == True:\n        print(\"Each Fold's MSE：{}, Average MSE：{:.4f}\".format([np.round(v,2) for v in scores], np.mean(scores)))\n        print(\"-\"*60)\n    return df_train, df_test, pred_valid, pred_test","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"df_train_x, df_test_x, pred_valid_x, pred_test_x = training(df_train, df_test, 'X', window_size)\ndf_train_y, df_test_y, pred_valid_y, pred_test_y = training(df_train, df_test, 'Y', window_size)\ndf_train_z, df_test_z, pred_valid_z, pred_test_z = training(df_train, df_test, 'Z', window_size)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val_compare_df = pd.DataFrame({'Xgt':df_train_x['Xgt'].values, 'Xpred':pred_valid_x,\n                               'Ygt':df_train_y['Ygt'].values, 'Ypred':pred_valid_y,\n                                'Zgt':df_train_z['Zgt'].values, 'Zpred':pred_valid_z})","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val_compare_df[['Zgt', 'Zpred']].plot(figsize=(16,8))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"val_compare_df[['Ygt', 'Ypred']].plot(figsize=(16,8))","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"We can see that, gt and pre is close to each other.","metadata":{}},{"cell_type":"code","source":"# xyz -> lng, lat\nlng_gt, lat_gt, _ = ECEF_to_WGS84(val_compare_df['Xgt'].values,val_compare_df['Ygt'].values,val_compare_df['Zgt'].values)\nlng_pred, lat_pred, _ = ECEF_to_WGS84(val_compare_df['Xpred'].values,val_compare_df['Ypred'].values,val_compare_df['Zpred'].values)\nlng_test_pred, lat_test_pred, _ = ECEF_to_WGS84(pred_test_x, pred_test_y, pred_test_z)\n\n    \nval_compare_df['latDeg_gt'] = lat_gt\nval_compare_df['lngDeg_gt'] = lng_gt\nval_compare_df['latDeg_pred'] = lat_pred\nval_compare_df['lngDeg_pred'] = lng_pred\ntest_pred_df = pd.DataFrame({'latDeg':lat_test_pred, 'lngDeg':lng_test_pred})","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# From：https://www.kaggle.com/emaerthin/demonstration-of-the-kalman-filter\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":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Baseline vs. GT\nlat_lng_df_train['dist'] = calc_haversine(lat_lng_df_train.latDeg_gt, lat_lng_df_train.lngDeg_gt, \n                                lat_lng_df_train.latDeg_bl, lat_lng_df_train.lngDeg_bl)\nprint('dist_50:',np.percentile(lat_lng_df_train['dist'],50) )\nprint('dist_95:',np.percentile(lat_lng_df_train['dist'],95) )\nprint('avg_dist_50_95:',(np.percentile(lat_lng_df_train['dist'],50) + np.percentile(lat_lng_df_train['dist'],95))/2)\nprint('avg_dist:', lat_lng_df_train['dist'].mean())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# IMU Prediction vs. GT\nval_compare_df['dist'] = calc_haversine(val_compare_df.latDeg_gt, val_compare_df.lngDeg_gt, \n                                val_compare_df.latDeg_pred, val_compare_df.lngDeg_pred)\n# IMU预测vsGT（多collection）\nprint('dist_50:',np.percentile(val_compare_df['dist'],50) )\nprint('dist_95:',np.percentile(val_compare_df['dist'],95) )\nprint('avg_dist_50_95:',(np.percentile(val_compare_df['dist'],50) + np.percentile(val_compare_df['dist'],95))/2)\nprint('avg_dist:', val_compare_df['dist'].mean())","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"def visualize_trafic(df, center, zoom=15):\n    fig = px.scatter_mapbox(df,\n                            \n                            # Here, plotly gets, (x,y) coordinates\n                            lat=\"latDeg\",\n                            lon=\"lngDeg\",\n                            \n                            #Here, plotly detects color of series\n                            color=\"phoneName\",\n                            labels=\"phoneName\",\n                            \n                            zoom=zoom,\n                            center=center,\n                            height=600,\n                            width=800)\n    fig.update_layout(mapbox_style='stamen-terrain')\n    fig.update_layout(margin={\"r\": 0, \"t\": 0, \"l\": 0, \"b\": 0})\n    fig.update_layout(title_text=\"GPS trafic\")\n    fig.show()\n    \ndef visualize_collection(df):\n    target_df = df\n    lat_center = target_df['latDeg'].mean()\n    lng_center = target_df['lngDeg'].mean()\n    center = {\"lat\":lat_center, \"lon\":lng_center}\n    \n    visualize_trafic(target_df, center)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualization: Train dataset\ncname = '2021-04-29-US-SJC-2'\npname = 'SamsungS20Ultra'\n# IMU Prediciton\ntmp0 = val_compare_df.copy()\ntmp0.rename(columns={'latDeg_pred':'latDeg', 'lngDeg_pred':'lngDeg'}, inplace=True)\ntmp0['phoneName'] = [cname + '_' + pname + '_imu_pred' for i in range(len(tmp0))]\n# GT\ntmp1 = val_compare_df.copy()\ntmp1.rename(columns={'latDeg_gt':'latDeg', 'lngDeg_gt':'lngDeg'}, inplace=True)\ntmp1['phoneName'] = [cname + '_' + pname + '_gt' for i in range(len(tmp1))]\n# Baseline\ntmp2 = lat_lng_df_train.copy()\ntmp2.rename(columns={'latDeg_bl':'latDeg', 'lngDeg_bl':'lngDeg'}, inplace=True)\ntmp2['phoneName'] = [cname + '_' + pname + '_bl_pred' for i in range(len(tmp2))]\n\ntmp = pd.concat([tmp0, tmp1, tmp2])\nvisualize_collection(tmp)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"code","source":"# Visualization:: Test dataset\ncname = '2021-04-29-US-SJC-3'\npname = 'SamsungS20Ultra'\ntmp3 = test_pred_df.copy()\ntmp3['phoneName'] = cname_test + '_' + pname_test + '_imu_pred' \n\ntmp4 = bl_tst_df.iloc[bl_tst_df[bl_tst_df['phone']==cname_test + '_' + pname_test].index[window_size:],3:5].copy()\ntmp4['phoneName'] = cname_test + '_' + pname_test + '_bl_pred' \n\ntmp5 = pd.concat([tmp3, tmp4])\nvisualize_collection(tmp5)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"## 4. Submission","metadata":{}},{"cell_type":"code","source":"bl_tst_df.iloc[bl_tst_df[bl_tst_df['phone']==cname_test + '_' + pname_test].index[window_size:],3] = test_pred_df['latDeg'].values\nbl_tst_df.iloc[bl_tst_df[bl_tst_df['phone']==cname_test + '_' + pname_test].index[window_size:],4] = test_pred_df['lngDeg'].values\n\n# bl_tst_df.to_csv('../submit/imu_baseline_locations_test.csv', index = False)","metadata":{"trusted":true},"execution_count":null,"outputs":[]},{"cell_type":"markdown","source":"The above code is just a template. I cannot promise it could improve a lot for you because i didn't publish my entire code (which contains some tricks). However, i hope my code could inspire you how to use GNSS log data. <br>\n\nOn the other hand, i am seeking for a nice teammate, here is my basic infomation:\n1. I am fresh of Kaggle, but attend some competitions in China before.\n2. I am Data Analyst from OPPO for more than one year (well, indoor game uses our phone's data).\n3. I am graduated from University of Manchester, UK. I guess my English can handle the normal communication.\n4. I have another teammate from china, too. (he is nice)\n5. I can use 1-3 hrs/day for playing this competition but the weekend i got one day for fee.\n\nWhat I hope:\n1. You are Top 50 or related experience about phone locating.\n2. You have time and energy to keep working on the game.\n3. You are like to publish the code when we finished this game.\n4. Most importance!!! you have great ideas!\n\nIf you want to join us, feel free to contact me. My email: **alvinai9603@outlook.com**<br>\n\nThanks.","metadata":{}}]}