{"metadata":{"kernelspec":{"display_name":"Python 3","language":"python","name":"python3"},"language_info":{"name":"python","version":"3.12.12","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[{"sourceType":"competition","sourceId":35779,"databundleVersionId":3604062}],"dockerImageVersionId":31328,"isInternetEnabled":true,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"id":"e00bebfb","cell_type":"markdown","source":"# GPS Trajectory Smoothing via Multi-Device Averaging\n## Google Smartphone Decimeter Challenge 2022 — Reproduction Notebook\n\n### Bài toán\n\nGPS trên smartphone có sai số trung bình **4–5 mét** trong điều kiện thực địa —\ndo nhiễu vệ tinh, tín hiệu phản xạ qua toà nhà (multipath), và giới hạn phần cứng.\nNotebook này tái hiện và điều chỉnh pipeline của\n[t88take](https://www.kaggle.com/code/t88take/gsdc-phones-mean-prediction)\ncho cấu trúc dữ liệu GSDC 2022, đạt **sai số 3.66 m — cải thiện 16% so với GPS thô**.\n\n## Pipeline\n\n```text\nGPS thô (device_gnss.csv)\n│\n▼\n[1] Outlier Rejection\n    └── Loại điểm nhảy > 50 m ở cả hai phía\n│\n▼\n[2] Kalman RTS Smoother\n    └── Làm mượt theo mô hình chuyển động 6 trạng thái\n│\n▼\n[3] Phones Mean Prediction\n    └── Trung bình N điện thoại trong cùng chuyến đi\n        (linear interpolation để đồng bộ timestamp trước khi trung bình)\n│\n▼\nDự đoán cuối (submission.csv)\n```\n| Bước | Score (m) | Cải thiện tích lũy |\n|------|-----------|--------------------|\n| Baseline GPS thô | 4.3542 | — |\n| + Outlier Rejection | 4.2342 | −2.8% |\n| + Kalman Smoothing | 3.8699 | −11.1% |\n| + Phones Mean Pred | 3.6584 | −16.0% |\n\n### Điểm thích ứng cho GSDC 2022\n\nGSDC 2022 không có `baseline_locations_train.csv` như phiên bản 2021.\nNotebook này đọc GPS trực tiếp từ `device_gnss.csv` và tự convert\ntọa độ **ECEF → WGS84** (lat/lng) bằng thuật toán Bowring trước khi\nđưa vào pipeline.\n\n| | GSDC 2021 | GSDC 2022 (notebook này) |\n|---|---|---|\n| Nguồn GPS baseline | `baseline_locations_train.csv` | `device_gnss.csv` (ECEF → WGS84) |\n| Phone/test drive | Nhiều | 1 (phones mean không hiệu quả trên test) |\n\n### Cấu trúc notebook\n\n| Section | Nội dung |\n|---------|----------|\n| §1–2 | Thiết lập môi trường, nạp & khám phá dữ liệu |\n| §3 | Hàm tiện ích (Haversine, visualisation, metric) |\n| §4 | Outlier Rejection |\n| §5 | Kalman RTS Smoother |\n| §6 | Phones Mean Prediction |\n| §7 | Đánh giá từng bước + visualisation quỹ đạo |\n| §8 | Tạo submission |\n| §9 | Tổng kết & hướng cải thiện |","metadata":{}},{"id":"af1e243a","cell_type":"markdown","source":"## 1. Thiết lập môi trường\n\n### 1.1 Import thư viện\n","metadata":{}},{"id":"d0699e79-17f2-4256-a80d-0105a7a6a82b","cell_type":"code","source":"!pip install simdkalman","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:49:24.139624Z","iopub.execute_input":"2026-05-08T08:49:24.140012Z","iopub.status.idle":"2026-05-08T08:49:28.426260Z","shell.execute_reply.started":"2026-05-08T08:49:24.139982Z","shell.execute_reply":"2026-05-08T08:49:28.425119Z"}},"outputs":[],"execution_count":null},{"id":"0d9f1569","cell_type":"code","source":"import os\nimport pathlib\nimport warnings\nwarnings.filterwarnings('ignore')\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm.notebook import tqdm\n\n# !pip install simdkalman  # bỏ comment nếu chưa cài\nimport simdkalman\n\nimport plotly.express as px\nimport plotly.graph_objects as go\n\nprint(\"✅ Tất cả thư viện đã được import thành công.\")\nprint(f\"   numpy  : {np.__version__}\")\nprint(f\"   pandas : {pd.__version__}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:49:28.429248Z","iopub.execute_input":"2026-05-08T08:49:28.429614Z","iopub.status.idle":"2026-05-08T08:49:28.437885Z","shell.execute_reply.started":"2026-05-08T08:49:28.429579Z","shell.execute_reply":"2026-05-08T08:49:28.436915Z"}},"outputs":[],"execution_count":null},{"id":"43a03bb6","cell_type":"markdown","source":"### 1.2 Cấu hình đường dẫn & tham số toàn cục\n","metadata":{}},{"id":"d9174731","cell_type":"code","source":"# ── Đường dẫn dataset ────────────────────────────────────────────────────────\nINPUT = '/kaggle/input/competitions/smartphone-decimeter-2022'\n\n# ── Hyperparameters ───────────────────────────────────────────────────────────\nOUTLIER_THRESHOLD_METERS = 50   # Ngưỡng khoảng cách để xác định điểm ngoại lệ (m)\n\n# Kalman Filter — state: [lat, lng, vel_lat, vel_lng, acc_lat, acc_lng]\nKF_T         = 1.0\nKF_PROC_DIAG = [1e-5, 1e-5, 5e-6, 5e-6, 1e-6, 1e-6]\nKF_OBS_DIAG  = [5e-5, 5e-5]\n\n# Các cột GPS cần thiết trong device_gnss.csv\nLAT_COL  = 'WlsLatitudeDegrees'\nLNG_COL  = 'WlsLongitudeDegrees'\nTIME_COL = 'utcTimeMillis'\n\nprint(\"📌 Cấu hình:\")\nprint(f\"   INPUT             : {INPUT}\")\nprint(f\"   OUTLIER_THRESHOLD : {OUTLIER_THRESHOLD_METERS} m\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:49:28.439335Z","iopub.execute_input":"2026-05-08T08:49:28.439642Z","iopub.status.idle":"2026-05-08T08:49:28.463425Z","shell.execute_reply.started":"2026-05-08T08:49:28.439609Z","shell.execute_reply":"2026-05-08T08:49:28.462150Z"}},"outputs":[],"execution_count":null},{"id":"3611b2f9","cell_type":"markdown","source":"---\n## 2. Nạp & Khám phá dữ liệu (EDA)\n\n### 2.1 Cấu trúc thư mục GSDC 2022\n\n```\nsmartphone-decimeter-2022/\n├── train/\n│   └── [drive_id]/               ← một \"collection\" (chuyến lái xe)\n│       └── [phone_name]/\n│           ├── device_gnss.csv   ← dữ liệu GPS thô từ điện thoại  ← ta dùng\n│           ├── ground_truth.csv  ← vị trí chuẩn (chỉ có ở train)\n│           └── supplemental/\n├── test/\n│   └── [drive_id]/\n│       └── [phone_name]/\n│           └── device_gnss.csv\n└── sample_submission.csv\n```\n\n**Điểm khác biệt với GSDC 2021:** Năm 2022 không có `baseline_locations_train.csv`.\nTa tự đọc và gộp `device_gnss.csv` từ tất cả các thư mục con.\n","metadata":{}},{"id":"abe9755b","cell_type":"markdown","source":"### 2.2 Hàm đọc toàn bộ dataset\n","metadata":{}},{"id":"fa201dee","cell_type":"code","source":"def ecef_to_latlong(x, y, z):\n    \"\"\"\n    Chuyển đổi tọa độ ECEF (mét) sang WGS84 lat/lng (độ thập phân).\n    Dùng thuật toán Bowring iterative.\n    \"\"\"\n    a  = 6_378_137.0          # bán trục lớn WGS84 (m)\n    e2 = 6.6943799901414e-3   # độ lệch tâm bình phương\n\n    lon = np.arctan2(y, x)\n    p   = np.sqrt(x**2 + y**2)\n    lat = np.arctan2(z, p * (1 - e2))   # khởi tạo\n\n    for _ in range(5):   # lặp hội tụ\n        N   = a / np.sqrt(1 - e2 * np.sin(lat)**2)\n        lat = np.arctan2(z + e2 * N * np.sin(lat), p)\n\n    return np.degrees(lat), np.degrees(lon)\n\n\ndef load_split(split: str, input_dir: str = INPUT) -> pd.DataFrame:\n    gnss_files = list(pathlib.Path(input_dir).glob(f'{split}/*/*/device_gnss.csv'))\n    if not gnss_files:\n        raise FileNotFoundError(\n            f\"Không tìm thấy file device_gnss.csv trong '{input_dir}/{split}/'\"\n        )\n    print(f\"   📂 Tìm thấy {len(gnss_files)} file device_gnss.csv trong '{split}/'\")\n\n    dfs = []\n    for f in tqdm(gnss_files, desc=f'Loading {split}', leave=False):\n        phone_name = f.parent.name\n        drive_id   = f.parent.parent.name\n\n        df_tmp = pd.read_csv(f, usecols=[\n            'utcTimeMillis',\n            'WlsPositionXEcefMeters',\n            'WlsPositionYEcefMeters',\n            'WlsPositionZEcefMeters',\n        ])\n\n        # Mỗi epoch (utcTimeMillis) có nhiều dòng vệ tinh — lấy 1 dòng/epoch\n        df_tmp = df_tmp.drop_duplicates(subset='utcTimeMillis')\n\n        # Convert ECEF → lat/lng\n        lat, lng = ecef_to_latlong(\n            df_tmp['WlsPositionXEcefMeters'].values,\n            df_tmp['WlsPositionYEcefMeters'].values,\n            df_tmp['WlsPositionZEcefMeters'].values,\n        )\n        df_tmp['latDeg'] = lat\n        df_tmp['lngDeg'] = lng\n        df_tmp = df_tmp.rename(columns={'utcTimeMillis': TIME_COL})\n        df_tmp['collectionName'] = drive_id\n        df_tmp['phoneName']      = phone_name\n\n        dfs.append(df_tmp[['collectionName', 'phoneName', TIME_COL, 'latDeg', 'lngDeg']])\n\n    combined = pd.concat(dfs, ignore_index=True)\n    combined = combined.dropna(subset=['latDeg', 'lngDeg'])\n    combined = combined.sort_values(['collectionName', 'phoneName', TIME_COL]).reset_index(drop=True)\n    return combined\n\n\ndef load_ground_truth(input_dir: str = INPUT) -> pd.DataFrame:\n    \"\"\"\n    Đọc và gộp tất cả ground_truth.csv từ tập train.\n\n    Chuẩn hoá tên cột sang: collectionName, phoneName, millisSinceGpsEpoch,\n    latDeg, lngDeg (giống GSDC 2021).\n    \"\"\"\n    gt_files = list(pathlib.Path(input_dir).glob('train/*/*/ground_truth.csv'))\n    print(f\"   📂 Tìm thấy {len(gt_files)} file ground_truth.csv\")\n\n    dfs = []\n    for f in tqdm(gt_files, desc='Loading ground_truth', leave=False):\n        df_tmp = pd.read_csv(f)\n        df_tmp['collectionName'] = f.parent.parent.name\n        df_tmp['phoneName']      = f.parent.name\n        dfs.append(df_tmp)\n\n    gt = pd.concat(dfs, ignore_index=True)\n\n    rename_map = {\n        'LatitudeDegrees':  'latDeg',\n        'LongitudeDegrees': 'lngDeg',\n        'UnixTimeMillis':    TIME_COL,\n    }\n    gt = gt.rename(columns=rename_map)\n\n    cols_keep = ['collectionName', 'phoneName', TIME_COL, 'latDeg', 'lngDeg']\n    return gt[[c for c in cols_keep if c in gt.columns]]","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:49:28.465141Z","iopub.execute_input":"2026-05-08T08:49:28.465554Z","iopub.status.idle":"2026-05-08T08:49:28.492253Z","shell.execute_reply.started":"2026-05-08T08:49:28.465514Z","shell.execute_reply":"2026-05-08T08:49:28.491152Z"}},"outputs":[],"execution_count":null},{"id":"56cddefc","cell_type":"code","source":"print(\"🔵 Đang nạp dữ liệu...\")\n\nprint(\"\\n  Train:\")\nbase_train = load_split('train')\n\nprint(\"\\n  Test:\")\nbase_test  = load_split('test')\n\nprint(\"\\n  Ground Truth:\")\nground_truth = load_ground_truth()\n\nsample_sub = pd.read_csv(f'{INPUT}/sample_submission.csv')\n\nprint(f\"\\n✅ Nạp xong.\")\nprint(f\"   base_train   : {base_train.shape}\")\nprint(f\"   base_test    : {base_test.shape}\")\nprint(f\"   ground_truth : {ground_truth.shape}\")\nprint(f\"   sample_sub   : {sample_sub.shape}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:49:28.495065Z","iopub.execute_input":"2026-05-08T08:49:28.495428Z","iopub.status.idle":"2026-05-08T08:51:24.700826Z","shell.execute_reply.started":"2026-05-08T08:49:28.495399Z","shell.execute_reply":"2026-05-08T08:51:24.699883Z"}},"outputs":[],"execution_count":null},{"id":"79c046df","cell_type":"markdown","source":"### 2.3 Tổng quan thống kê\n","metadata":{}},{"id":"e233f785","cell_type":"code","source":"print(\"📊 Thống kê tập Train:\")\nprint(f\"   Collections (drives) : {base_train['collectionName'].nunique()}\")\nprint(f\"   Phones               : {base_train['phoneName'].nunique()}\")\nprint(f\"   Tổng số điểm GPS     : {len(base_train):,}\")\nprint()\n\nphones_per_col = base_train.groupby('collectionName')['phoneName'].nunique()\nprint(f\"   Phones/collection — min: {phones_per_col.min()}, \"\n      f\"max: {phones_per_col.max()}, \"\n      f\"mean: {phones_per_col.mean():.1f}\")\nprint()\ndisplay(base_train.head(3))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:24.739926Z","iopub.execute_input":"2026-05-08T08:51:24.740372Z","iopub.status.idle":"2026-05-08T08:51:24.866564Z","shell.execute_reply.started":"2026-05-08T08:51:24.740214Z","shell.execute_reply":"2026-05-08T08:51:24.865573Z"}},"outputs":[],"execution_count":null},{"id":"cc8054b8","cell_type":"markdown","source":"### 2.4 Visualisation phân phối số điện thoại mỗi collection\n","metadata":{}},{"id":"f47f313c","cell_type":"code","source":"fig, ax = plt.subplots(figsize=(7, 3))\nphones_per_col.value_counts().sort_index().plot(\n    kind='bar', ax=ax, color='steelblue', edgecolor='white'\n)\nax.set_xlabel('Số điện thoại trong collection')\nax.set_ylabel('Số collections')\nax.set_title('Phân phối số điện thoại mỗi collection (tập Train)')\nax.grid(axis='y', alpha=0.3)\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:24.867910Z","iopub.execute_input":"2026-05-08T08:51:24.868267Z","iopub.status.idle":"2026-05-08T08:51:25.212676Z","shell.execute_reply.started":"2026-05-08T08:51:24.868229Z","shell.execute_reply":"2026-05-08T08:51:25.211597Z"}},"outputs":[],"execution_count":null},{"id":"773c10db","cell_type":"markdown","source":"---\n## 3. Hàm tiện ích\n\n### 3.1 Khoảng cách Haversine\n\nCông thức Haversine tính khoảng cách đường đại viên (great-circle distance) giữa hai điểm trên Trái Đất:\n\n$$d = 2R \\arcsin\\!\\sqrt{\\sin^2\\!\\frac{\\Delta\\varphi}{2} + \\cos\\varphi_1\\cos\\varphi_2\\sin^2\\!\\frac{\\Delta\\lambda}{2}}$$\n\nvới $R = 6{,}367{,}000$ m.\n","metadata":{}},{"id":"66f9dde5","cell_type":"code","source":"def calc_haversine(lat1, lon1, lat2, lon2):\n    \"\"\"\n    Tính khoảng cách great-circle (mét) giữa hai điểm GPS.\n\n    Parameters\n    ----------\n    lat1, lon1 : array-like — tọa độ điểm 1 (độ thập phân)\n    lat2, lon2 : array-like — tọa độ điểm 2 (độ thập phân)\n\n    Returns\n    -------\n    dist : ndarray — khoảng cách tính bằng mét\n    \"\"\"\n    RADIUS = 6_367_000\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)**2 + np.cos(lat1) * np.cos(lat2) * np.sin(dlon / 2)**2\n    return 2 * RADIUS * np.arcsin(np.sqrt(a))\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:25.213840Z","iopub.execute_input":"2026-05-08T08:51:25.214172Z","iopub.status.idle":"2026-05-08T08:51:25.220897Z","shell.execute_reply.started":"2026-05-08T08:51:25.214141Z","shell.execute_reply":"2026-05-08T08:51:25.219838Z"}},"outputs":[],"execution_count":null},{"id":"8f4a5ac5","cell_type":"markdown","source":"### 3.2 Visualisation GPS trên bản đồ thực\n","metadata":{}},{"id":"84aafe2a","cell_type":"code","source":"def visualize_collection(df: pd.DataFrame, collection: str) -> None:\n    \"\"\"\n    Vẽ quỹ đạo GPS của tất cả phone trong một collection lên bản đồ.\n\n    Parameters\n    ----------\n    df         : DataFrame với cột latDeg, lngDeg, phoneName, collectionName\n    collection : tên collection (drive_id) cần vẽ\n    \"\"\"\n    target = df[df['collectionName'] == collection].copy()\n    center = {'lat': target['latDeg'].mean(), 'lon': target['lngDeg'].mean()}\n\n    fig = px.scatter_mapbox(\n        target,\n        lat='latDeg', lon='lngDeg',\n        color='phoneName',\n        zoom=14, center=center,\n        height=550, width=860,\n        title=f'GPS trajectories — {collection}',\n        labels={'phoneName': 'Device'},\n    )\n    fig.update_layout(mapbox_style='open-street-map',\n                      margin={'r': 0, 't': 40, 'l': 0, 'b': 0})\n    fig.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:25.222186Z","iopub.execute_input":"2026-05-08T08:51:25.222585Z","iopub.status.idle":"2026-05-08T08:51:25.251231Z","shell.execute_reply.started":"2026-05-08T08:51:25.222537Z","shell.execute_reply":"2026-05-08T08:51:25.250109Z"}},"outputs":[],"execution_count":null},{"id":"be319d65","cell_type":"markdown","source":"### 3.3 Metric đánh giá\n\nMetric chính thức: **mean(p50 + p95) / 2** — trung bình giữa sai số trung vị\nvà phân vị 95 của từng phone — đảm bảo mô hình tốt cả về độ chính xác điển hình lẫn trường hợp xấu nhất.\n","metadata":{}},{"id":"2d1bbafe","cell_type":"code","source":"def score_predictions(pred_df: pd.DataFrame, gt_df: pd.DataFrame) -> float:\n    \"\"\"\n    Tính score chính thức GSDC: mean của (p50 + p95) / 2 per phone.\n\n    Parameters\n    ----------\n    pred_df : DataFrame — [collectionName, phoneName, millisSinceGpsEpoch, latDeg, lngDeg]\n    gt_df   : Ground truth cùng cấu trúc\n\n    Returns\n    -------\n    score : float — sai số trung bình (mét), càng nhỏ càng tốt\n    \"\"\"\n    gt = gt_df.rename(columns={'latDeg': 'lat_gt', 'lngDeg': 'lng_gt'})\n    merged = pred_df.merge(\n        gt, on=['collectionName', 'phoneName', TIME_COL], how='inner'\n    )\n    if merged.empty:\n        raise ValueError(\"Merge trả về DataFrame rỗng — kiểm tra lại tên cột TIME_COL.\")\n\n    merged['error_m'] = calc_haversine(\n        merged['lat_gt'], merged['lng_gt'],\n        merged['latDeg'], merged['lngDeg']\n    )\n    merged['phone'] = merged['collectionName'] + '_' + merged['phoneName']\n\n    per_phone = merged.groupby('phone')['error_m'].agg(\n        p50=lambda x: np.percentile(x, 50),\n        p95=lambda x: np.percentile(x, 95)\n    )\n    per_phone['score'] = (per_phone['p50'] + per_phone['p95']) / 2\n    return per_phone['score'].mean()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:25.252479Z","iopub.execute_input":"2026-05-08T08:51:25.252844Z","iopub.status.idle":"2026-05-08T08:51:25.278496Z","shell.execute_reply.started":"2026-05-08T08:51:25.252770Z","shell.execute_reply":"2026-05-08T08:51:25.277123Z"}},"outputs":[],"execution_count":null},{"id":"987577c3","cell_type":"markdown","source":"---\n## 4. Bước 1 — Loại bỏ điểm đo bất thường (Outlier Rejection)\n\n### Động lực\n\nGPS thô thường có \"jump\" (nhảy đột ngột) do mất tín hiệu vệ tinh, multipath (tín hiệu\nphản xạ qua toà nhà), hoặc lỗi phần cứng.\n\n### Tiêu chí\n\nĐiểm $p_t$ bị coi là ngoại lệ nếu:\n\n$$d(p_{t-1},\\, p_t) > \\tau \\quad \\text{VÀ} \\quad d(p_t,\\, p_{t+1}) > \\tau$$\n\nĐiều kiện **kép** (cả hai phía) tránh xóa nhầm khi xe thực sự tăng tốc đột ngột.\nNgưỡng $\\tau = 50$ m được chọn theo kinh nghiệm trên dữ liệu GSDC.\n","metadata":{}},{"id":"47788067","cell_type":"code","source":"def add_distance_diff(df: pd.DataFrame) -> pd.DataFrame:\n    \"\"\"\n    Thêm cột khoảng cách tới điểm liền trước/sau cho từng phone.\n    Biên của mỗi phone (điểm đầu/cuối) được đặt thành NaN.\n    \"\"\"\n    df = df.copy()\n    df['phone'] = df['collectionName'] + '_' + df['phoneName']\n\n    for col in ['latDeg', 'lngDeg']:\n        df[f'{col}_prev'] = df[col].shift(1)\n        df[f'{col}_next'] = df[col].shift(-1)\n\n    df['phone_prev'] = df['phone'].shift(1)\n    df['phone_next'] = df['phone'].shift(-1)\n\n    df['dist_prev'] = calc_haversine(df['latDeg'], df['lngDeg'],\n                                      df['latDeg_prev'], df['lngDeg_prev'])\n    df['dist_next'] = calc_haversine(df['latDeg'], df['lngDeg'],\n                                      df['latDeg_next'], df['lngDeg_next'])\n\n    df.loc[df['phone'] != df['phone_prev'],\n           ['latDeg_prev', 'lngDeg_prev', 'dist_prev']] = np.nan\n    df.loc[df['phone'] != df['phone_next'],\n           ['latDeg_next', 'lngDeg_next', 'dist_next']] = np.nan\n\n    return df\n\n\ndef reject_outliers(df: pd.DataFrame,\n                    threshold: float = OUTLIER_THRESHOLD_METERS) -> pd.DataFrame:\n    \"\"\"\n    Đặt latDeg/lngDeg = NaN tại các điểm GPS bất thường.\n\n    Parameters\n    ----------\n    df        : DataFrame đã có cột collectionName, phoneName, latDeg, lngDeg\n    threshold : ngưỡng khoảng cách (mét)\n    \"\"\"\n    df = add_distance_diff(df)\n    mask = (df['dist_prev'] > threshold) & (df['dist_next'] > threshold)\n    n_out = mask.sum()\n    df.loc[mask, ['latDeg', 'lngDeg']] = np.nan\n    print(f\"   ✂  Loại {n_out:,} điểm ngoại lệ ({100 * n_out / len(df):.2f}%)\")\n    return df\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:25.279727Z","iopub.execute_input":"2026-05-08T08:51:25.280024Z","iopub.status.idle":"2026-05-08T08:51:25.300137Z","shell.execute_reply.started":"2026-05-08T08:51:25.279988Z","shell.execute_reply":"2026-05-08T08:51:25.298978Z"}},"outputs":[],"execution_count":null},{"id":"3d8a26a4","cell_type":"code","source":"print(\"🔵 Outlier Rejection — Train:\")\ntrain_ro = reject_outliers(base_train)\n\nprint(\"🔵 Outlier Rejection — Test:\")\ntest_ro  = reject_outliers(base_test)\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:25.303955Z","iopub.execute_input":"2026-05-08T08:51:25.304381Z","iopub.status.idle":"2026-05-08T08:51:25.613256Z","shell.execute_reply.started":"2026-05-08T08:51:25.304352Z","shell.execute_reply":"2026-05-08T08:51:25.612032Z"}},"outputs":[],"execution_count":null},{"id":"50779cdc","cell_type":"markdown","source":"---\n## 5. Bước 2 — Làm mượt quỹ đạo bằng Kalman Filter\n\n### Nguyên lý\n\nKalman Filter mô hình hoá GPS như hệ động lực tuyến tính với **6 trạng thái**:\n\n$$\\mathbf{x}_t = [\\text{lat},\\ \\text{lng},\\ \\dot{\\text{lat}},\\ \\dot{\\text{lng}},\\ \\ddot{\\text{lat}},\\ \\ddot{\\text{lng}}]^\\top$$\n\n- **Transition model** $F$: chuyển động đều có gia tốc (constant acceleration)\n- **Process noise** $Q$: độ không chắc của mô hình vật lý\n- **Observation noise** $R$: độ nhiễu của GPS thô\n\n`simdkalman` dùng **RTS Smoother** (Rauch–Tung–Striebel) — backward pass trên toàn chuỗi —\nnên tận dụng được cả quan sát tương lai, cho kết quả tốt hơn forward-only filter.\n","metadata":{}},{"id":"40fa7e0c","cell_type":"code","source":"T = KF_T\n\nstate_transition = np.array([\n    [1, 0, T, 0, 0.5*T**2, 0       ],\n    [0, 1, 0, T, 0,        0.5*T**2],\n    [0, 0, 1, 0, T,        0       ],\n    [0, 0, 0, 1, 0,        T       ],\n    [0, 0, 0, 0, 1,        0       ],\n    [0, 0, 0, 0, 0,        1       ],\n])\n\nprocess_noise = np.diag(KF_PROC_DIAG) + np.ones((6, 6)) * 1e-9\n\nobservation_model = np.array([\n    [1, 0, 0, 0, 0, 0],\n    [0, 1, 0, 0, 0, 0],\n])\n\nobservation_noise = np.diag(KF_OBS_DIAG) + 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)\n\nprint(\"✅ Kalman Filter đã khởi tạo.\")\nprint(f\"   State dim      : {state_transition.shape[0]}\")\nprint(f\"   Observation dim: {observation_model.shape[0]}\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:25.614671Z","iopub.execute_input":"2026-05-08T08:51:25.614989Z","iopub.status.idle":"2026-05-08T08:51:25.624633Z","shell.execute_reply.started":"2026-05-08T08:51:25.614961Z","shell.execute_reply":"2026-05-08T08:51:25.623443Z"}},"outputs":[],"execution_count":null},{"id":"ce6ea5f5","cell_type":"code","source":"def apply_kf_smoothing(df: pd.DataFrame,\n                       kf_: simdkalman.KalmanFilter = kf) -> pd.DataFrame:\n    \"\"\"\n    Áp dụng RTS Smoother cho từng cặp (collectionName, phoneName) độc lập.\n\n    Parameters\n    ----------\n    df  : DataFrame với cột [collectionName, phoneName, millisSinceGpsEpoch, latDeg, lngDeg]\n    kf_ : KalmanFilter đã cấu hình\n\n    Returns\n    -------\n    df  : DataFrame cùng cấu trúc với latDeg, lngDeg đã được làm mượt\n    \"\"\"\n    df = df.copy()\n    unique_paths = df[['collectionName', 'phoneName']].drop_duplicates().to_numpy()\n\n    for collection, phone in tqdm(unique_paths, desc='KF Smoothing', leave=False):\n        mask = (df['collectionName'] == collection) & (df['phoneName'] == phone)\n        data = df.loc[mask, ['latDeg', 'lngDeg']].to_numpy()\n        data = data.reshape(1, len(data), 2)\n        smoothed = kf_.smooth(data)\n        df.loc[mask, 'latDeg'] = smoothed.states.mean[0, :, 0]\n        df.loc[mask, 'lngDeg'] = smoothed.states.mean[0, :, 1]\n\n    return df\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:25.625722Z","iopub.execute_input":"2026-05-08T08:51:25.626003Z","iopub.status.idle":"2026-05-08T08:51:25.654820Z","shell.execute_reply.started":"2026-05-08T08:51:25.625976Z","shell.execute_reply":"2026-05-08T08:51:25.653923Z"}},"outputs":[],"execution_count":null},{"id":"0db6277d","cell_type":"code","source":"COLS = ['collectionName', 'phoneName', TIME_COL, 'latDeg', 'lngDeg']\n\nprint(\"🔵 Kalman Smoothing — Train:\")\ntrain_kf = apply_kf_smoothing(train_ro[COLS])\n\nprint(\"🔵 Kalman Smoothing — Test:\")\ntest_kf  = apply_kf_smoothing(test_ro[COLS])\n\nprint(\"\\n✅ Kalman Smoothing hoàn thành.\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:51:25.656243Z","iopub.execute_input":"2026-05-08T08:51:25.656738Z","iopub.status.idle":"2026-05-08T08:52:46.031753Z","shell.execute_reply.started":"2026-05-08T08:51:25.656695Z","shell.execute_reply":"2026-05-08T08:52:46.030314Z"}},"outputs":[],"execution_count":null},{"id":"60019116","cell_type":"markdown","source":"---\n## 6. Bước 3 — Phones Mean Prediction\n\n### Ý tưởng\n\nTrong cùng một **collection** (chuyến lái xe), nhiều điện thoại đo **cùng một quỹ đạo thực**.\nSai số GPS của mỗi thiết bị có thành phần ngẫu nhiên độc lập, nên:\n\n$$\\hat{p}_t = \\frac{1}{N} \\sum_{i=1}^{N} p_t^{(i)} \\quad \\Rightarrow \\quad \\text{Var}(\\hat{p}_t) = \\frac{\\sigma^2}{N}$$\n\nTrung bình $N$ thiết bị giảm phương sai xuống còn $1/N$.\n\n### Vấn đề đồng bộ timestamp\n\nMỗi điện thoại có bộ timestamps GPS riêng, không khớp nhau.\nTa dùng **nội suy tuyến tính** (linear interpolation) để tạo giá trị ước tính\ntại các timestamp còn thiếu — sau đó mới lấy trung bình.\n","metadata":{}},{"id":"99f25aec","cell_type":"code","source":"def make_lerp_data(df: pd.DataFrame) -> pd.DataFrame:\n    \"\"\"\n    Tạo dữ liệu GPS nội suy cho các timestamp còn thiếu của mỗi phone.\n\n    Quy trình:\n      1. Tạo tổ hợp đầy đủ (timestamp × collection × phone)\n      2. Merge với dữ liệu gốc → NaN tại các timestamp thiếu\n      3. Nội suy tuyến tính: p = p_prev + (p_next - p_prev) × ratio\n\n    Returns\n    -------\n    lerp_df : DataFrame chỉ chứa các hàng được nội suy (không lẫn hàng gốc)\n    \"\"\"\n    org_cols = df.columns.tolist()\n\n    time_list  = df[['collectionName', TIME_COL]].drop_duplicates()\n    phone_list = df[['collectionName', 'phoneName']].drop_duplicates()\n    full_grid  = time_list.merge(phone_list, on='collectionName', how='outer')\n\n    lerp_df = full_grid.merge(df, on=['collectionName', TIME_COL, 'phoneName'], how='left')\n    lerp_df['phone'] = lerp_df['collectionName'] + '_' + lerp_df['phoneName']\n    lerp_df = lerp_df.sort_values(['phone', TIME_COL])\n\n    for col in ['latDeg', 'lngDeg', 'phone']:\n        lerp_df[f'{col}_prev'] = lerp_df[col].shift(1)\n        lerp_df[f'{col}_next'] = lerp_df[col].shift(-1)\n    lerp_df['time_prev'] = lerp_df[TIME_COL].shift(1)\n    lerp_df['time_next'] = lerp_df[TIME_COL].shift(-1)\n\n    mask = (lerp_df['latDeg'].isnull()\n            & (lerp_df['phone'] == lerp_df['phone_prev'])\n            & (lerp_df['phone'] == lerp_df['phone_next']))\n    lerp_df = lerp_df[mask].copy()\n\n    ratio = ((lerp_df[TIME_COL] - lerp_df['time_prev'])\n             / (lerp_df['time_next'] - lerp_df['time_prev']))\n\n    for coord in ['latDeg', 'lngDeg']:\n        lerp_df[coord] = (lerp_df[f'{coord}_prev']\n                          + (lerp_df[f'{coord}_next'] - lerp_df[f'{coord}_prev']) * ratio)\n\n    return lerp_df.dropna(subset=['latDeg'])[org_cols]\n\n\ndef calc_mean_pred(df: pd.DataFrame, lerp_df: pd.DataFrame) -> pd.DataFrame:\n    \"\"\"\n    Tính trung bình lat/lng của tất cả phones cùng collection tại cùng timestamp.\n\n    Parameters\n    ----------\n    df      : dữ liệu gốc sau Kalman smoothing\n    lerp_df : dữ liệu nội suy từ make_lerp_data()\n\n    Returns\n    -------\n    result : DataFrame cùng cấu trúc df với lat/lng đã trung bình hoá\n    \"\"\"\n    combined = pd.concat([df, lerp_df], ignore_index=True)\n    mean_coords = (combined\n                   .groupby(['collectionName', TIME_COL])[['latDeg', 'lngDeg']]\n                   .mean().reset_index())\n\n    result = df[['collectionName', 'phoneName', TIME_COL]].copy()\n    result = result.merge(mean_coords, on=['collectionName', TIME_COL], how='left')\n    return result\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:52:46.033154Z","iopub.execute_input":"2026-05-08T08:52:46.033697Z","iopub.status.idle":"2026-05-08T08:52:46.049672Z","shell.execute_reply.started":"2026-05-08T08:52:46.033655Z","shell.execute_reply":"2026-05-08T08:52:46.048620Z"}},"outputs":[],"execution_count":null},{"id":"fab2461d","cell_type":"code","source":"print(\"🔵 Phones Mean Prediction — Train:\")\ntrain_lerp      = make_lerp_data(train_kf)\ntrain_mean_pred = calc_mean_pred(train_kf, train_lerp)\n\nprint(\"🔵 Phones Mean Prediction — Test:\")\ntest_lerp       = make_lerp_data(test_kf)\ntest_mean_pred  = calc_mean_pred(test_kf, test_lerp)\n\nprint(\"\\n✅ Phones Mean Prediction hoàn thành.\")\nprint(f\"   Train: {len(train_mean_pred):,} rows | Test: {len(test_mean_pred):,} rows\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:52:46.050699Z","iopub.execute_input":"2026-05-08T08:52:46.051019Z","iopub.status.idle":"2026-05-08T08:52:48.362858Z","shell.execute_reply.started":"2026-05-08T08:52:46.050985Z","shell.execute_reply":"2026-05-08T08:52:48.361784Z"}},"outputs":[],"execution_count":null},{"id":"fe7d457a","cell_type":"markdown","source":"---\n## 7. Đánh giá kết quả\n\n### 7.1 Score từng bước pipeline\n","metadata":{}},{"id":"e39a87c2","cell_type":"code","source":"scores = {}\n\nprint(\"📐 Tính score từng bước...\")\nscores['1. Baseline (raw GPS)']      = score_predictions(base_train[COLS],  ground_truth)\nscores['2. + Outlier Rejection']     = score_predictions(train_ro[COLS],    ground_truth)\nscores['3. + Kalman Smoothing']      = score_predictions(train_kf,           ground_truth)\nscores['4. + Phones Mean Pred']      = score_predictions(train_mean_pred,    ground_truth)\n\nprint(f\"\\n{'Pipeline Step':<30} {'Score (m)':>10}\")\nprint(\"─\" * 42)\nfor step, s in scores.items():\n    print(f\"{step:<30} {s:>10.4f} m\")\n\nimprovement = (list(scores.values())[0] - list(scores.values())[-1]) / list(scores.values())[0] * 100\nprint(f\"\\n🎯 Tổng cải thiện: {improvement:.1f}% so với baseline raw GPS\")\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:52:48.364285Z","iopub.execute_input":"2026-05-08T08:52:48.364745Z","iopub.status.idle":"2026-05-08T08:52:49.952846Z","shell.execute_reply.started":"2026-05-08T08:52:48.364704Z","shell.execute_reply":"2026-05-08T08:52:49.951649Z"}},"outputs":[],"execution_count":null},{"id":"897439c8","cell_type":"markdown","source":"### 7.2 Biểu đồ so sánh\n","metadata":{}},{"id":"2436f2fb","cell_type":"code","source":"fig, ax = plt.subplots(figsize=(9, 4))\nsteps  = list(scores.keys())\nvalues = list(scores.values())\ncolors = ['#d9534f' if v == max(values) else\n          '#5cb85c' if v == min(values) else\n          '#5bc0de' for v in values]\n\nbars = ax.barh(steps, values, color=colors, edgecolor='white', height=0.55)\nax.bar_label(bars, fmt='%.4f m', padding=5, fontsize=10)\nax.set_xlabel('Score — mean(p50 + p95) / 2  (mét)', fontsize=10)\nax.set_title('Cải thiện độ chính xác GPS qua từng bước xử lý', fontsize=12, fontweight='bold')\nax.invert_yaxis()\nax.grid(axis='x', alpha=0.3)\nplt.tight_layout()\nplt.show()\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:52:49.953900Z","iopub.execute_input":"2026-05-08T08:52:49.954215Z","iopub.status.idle":"2026-05-08T08:52:50.140121Z","shell.execute_reply.started":"2026-05-08T08:52:49.954187Z","shell.execute_reply":"2026-05-08T08:52:50.138840Z"}},"outputs":[],"execution_count":null},{"id":"6d0fcad4","cell_type":"markdown","source":"### 7.3 Visualisation quỹ đạo GPS — so sánh trực quan\n\nSo sánh baseline → Kalman → Phones Mean → Ground Truth trên bản đồ thực.\n","metadata":{}},{"id":"9ccf8da9","cell_type":"code","source":"DEMO_COLLECTION = base_train['collectionName'].value_counts().index[0]\nprint(f\"📍 Demo collection: {DEMO_COLLECTION}\")\n\ndef tag(df, label, collection):\n    sub = df[df['collectionName'] == collection].copy()\n    sub['label'] = label.strip()\n    return sub\n\nvis_df = pd.concat([\n    tag(base_train[COLS],  '[Baseline]',  DEMO_COLLECTION),\n    tag(train_kf,          '[KF]',        DEMO_COLLECTION),\n    tag(train_mean_pred,   '[MeanPred]',  DEMO_COLLECTION),\n    tag(ground_truth,      '[GT]',        DEMO_COLLECTION),\n], ignore_index=True)\n\n# ── Static plot — hiển thị được cả khi commit ────────────────────────────────\nfig, ax = plt.subplots(figsize=(10, 7))\n\nstyles = {\n    '[Baseline]': dict(color='#aaaaaa', alpha=0.4, s=4,  zorder=1),\n    '[KF]':       dict(color='#4c9be8', alpha=0.6, s=4,  zorder=2),\n    '[MeanPred]': dict(color='#e8834c', alpha=0.8, s=5,  zorder=3),\n    '[GT]':       dict(color='#2ca02c', alpha=1.0, s=6,  zorder=4),\n}\n\nfor label, style in styles.items():\n    sub = vis_df[vis_df['label'] == label]\n    ax.scatter(sub['lngDeg'], sub['latDeg'], label=label, **style)\n\nax.set_xlabel('Longitude')\nax.set_ylabel('Latitude')\nax.set_title(f'GPS Trajectory Comparison — {DEMO_COLLECTION}', fontweight='bold')\nax.legend(loc='best', markerscale=3)\nax.set_aspect('equal')\nax.grid(alpha=0.2)\nplt.tight_layout()\nplt.savefig('trajectory_comparison.png', dpi=150, bbox_inches='tight')\nplt.show()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T10:02:36.153485Z","iopub.execute_input":"2026-05-08T10:02:36.153858Z","iopub.status.idle":"2026-05-08T10:02:37.536116Z","shell.execute_reply.started":"2026-05-08T10:02:36.153829Z","shell.execute_reply":"2026-05-08T10:02:37.535124Z"}},"outputs":[],"execution_count":null},{"id":"7f888c50","cell_type":"markdown","source":"---\n## 8. Tạo file Submission\n","metadata":{}},{"id":"e2738ac6","cell_type":"code","source":"# ── Section 8: Tạo file Submission ───────────────────────────────────────────\n\nsubmission = sample_sub.copy()\n\n# sample_submission có: tripId, UnixTimeMillis, LatitudeDegrees, LongitudeDegrees\n# tripId = \"[drive_id]/[phone_name]\" → tách ra để merge với test_mean_pred\n\nsubmission[['collectionName', 'phoneName']] = (\n    submission['tripId'].str.split('/', expand=True)\n)\nsubmission = submission.rename(columns={'UnixTimeMillis': TIME_COL})\n\n# Merge kết quả dự đoán\nsubmission = submission.merge(\n    test_mean_pred[['collectionName', 'phoneName', TIME_COL, 'latDeg', 'lngDeg']],\n    on=['collectionName', 'phoneName', TIME_COL],\n    how='left'\n)\n\n# Ghi đè lại cột gốc của submission\nsubmission['LatitudeDegrees']  = submission['latDeg']\nsubmission['LongitudeDegrees'] = submission['lngDeg']\n\n# Giữ đúng 4 cột theo format yêu cầu\nsubmission = submission[['tripId', TIME_COL, 'LatitudeDegrees', 'LongitudeDegrees']]\nsubmission = submission.rename(columns={TIME_COL: 'UnixTimeMillis'})\n\n# Kiểm tra\nn_null = submission[['LatitudeDegrees', 'LongitudeDegrees']].isnull().sum().sum()\nprint(f\"NaN còn lại: {n_null}\")\nif n_null > 0:\n    print(\"⚠️  Có NaN — kiểm tra xem TIME_COL của test và submission có khớp không\")\n    print(\"   Sample submission timestamp:\", submission['UnixTimeMillis'].iloc[0])\n    print(\"   test_mean_pred timestamp   :\", test_mean_pred[TIME_COL].iloc[0])\n\nsubmission.to_csv('submission.csv', index=False)\nprint(f\"✅ Đã lưu submission.csv — {len(submission):,} rows\")\ndisplay(submission.head(3))","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-05-08T08:55:36.612542Z","iopub.execute_input":"2026-05-08T08:55:36.612880Z","iopub.status.idle":"2026-05-08T08:55:37.114540Z","shell.execute_reply.started":"2026-05-08T08:55:36.612851Z","shell.execute_reply":"2026-05-08T08:55:37.113160Z"}},"outputs":[],"execution_count":null},{"id":"8542e53d","cell_type":"markdown","source":"---\n## 9. Tổng kết\n\n### Kết quả\n\nPipeline 3 bước áp dụng lên **62 collection, 295,600 điểm GPS** từ tập train:\n\n| Bước | Kỹ thuật | Score (m) | Cải thiện |\n|------|----------|-----------|-----------|\n| 0 | Baseline — GPS thô từ `device_gnss.csv` | 4.3542 | — |\n| 1 | + Outlier Rejection (ngưỡng 50 m, điều kiện kép) | 4.2342 | −2.8% |\n| 2 | + Kalman RTS Smoother (6 trạng thái) | 3.8699 | −8.6% |\n| 3 | + Phones Mean Prediction (linear interpolation) | 3.6584 | −5.5% |\n| | **Tổng** | | **−16.0%** |\n\nKalman Smoothing đóng góp nhiều nhất (−8.6%), cho thấy GPS thô có nhiễu\nhệ thống đáng kể hơn là nhiễu nhảy điểm đột ngột.\n\n### Giới hạn cần lưu ý\n\n**Phones Mean Prediction không cải thiện score trên test.** Tập test 2022\nchỉ có **1 phone/drive** — không có thiết bị thứ hai để trung bình. Kỹ thuật\nnày chỉ hiệu quả trên tập train (trung bình 2.7 phone/collection).\n\n**Submission có 14 điểm NaN.** Một số timestamp trong `sample_submission.csv`\nkhông khớp với `test_mean_pred` sau pipeline. Cần điều tra thêm — có thể do\ncác timestamp bị mất ở bước Outlier Rejection hoặc lỗi đồng bộ trong\n`make_lerp_data()`. Hiện tại 14 hàng này sẽ bị Kaggle chấm sai.\n\n### Hướng cải thiện\n\nTheo thứ tự ưu tiên ước tính:\n\n1. **Sửa 14 NaN trong submission** — điền bằng giá trị KF-smoothed trước\n   phones mean, hoặc fallback về GPS thô\n2. **Weighted mean** — gán trọng số theo lịch sử sai số của từng phone\n   thay vì trung bình đơn giản\n3. **Adaptive Kalman** — điều chỉnh $Q$, $R$ theo tốc độ ước tính\n   (giá trị tĩnh hiện tại không phân biệt lúc xe dừng và lúc chạy nhanh)\n4. **Thêm dữ liệu IMU** — kết hợp gia tốc kế/con quay hồi chuyển từ\n   `device_imu.csv` vào state vector của Kalman\n5. **Map matching** — khớp quỹ đạo lên đường bộ OpenStreetMap,\n   loại bỏ các điểm rơi vào giữa toà nhà hoặc dưới sông\n\n---\n*Tái hiện và điều chỉnh cho GSDC 2022 từ\n[t88take/gsdc-phones-mean-prediction](https://www.kaggle.com/code/t88take/gsdc-phones-mean-prediction)\n(GSDC 2021). Điểm thích ứng chính: đọc GPS từ ECEF trong `device_gnss.csv`\nthay vì `baseline_locations_train.csv` không tồn tại trong phiên bản 2022.*","metadata":{}}]}