{"metadata":{"kernelspec":{"language":"python","display_name":"Python 3","name":"python3"},"language_info":{"name":"python","version":"3.12.13","mimetype":"text/x-python","codemirror_mode":{"name":"ipython","version":3},"pygments_lexer":"ipython3","nbconvert_exporter":"python","file_extension":".py"},"kaggle":{"accelerator":"none","dataSources":[],"dockerImageVersionId":28755,"isInternetEnabled":false,"language":"python","sourceType":"notebook","isGpuEnabled":false}},"nbformat_minor":4,"nbformat":4,"cells":[{"cell_type":"code","source":"\"\"\"\n════════════════════════════════════════════════════════════════════════════\nGNSS Multipath Detection (CNN) + ML-Weighted WLS Positioning\nAdapted for: Google Smartphone Decimeter Challenge 2023 (sdc2023) dataset\n════════════════════════════════════════════════════════════════════════════\n\nWHAT CHANGED FROM THE OLD (NavIC, single static benchmark) PIPELINE\n--------------------------------------------------------------------\nOld dataset : DATASET_ROOT/<trace>/<device>/*.txt  (raw GnssLogger text logs)\n              + one hand-surveyed, SINGLE static ground-truth point for the\n              whole dataset (the receiver never moved).\nNew dataset : COMPETITION_ROOT/{train,test}/<drive_id>/<phone_name>/\n                  device_gnss.csv   <- ALREADY PARSED raw+derived measurements\n                  ground_truth.csv  <- per-epoch, MOVING receiver trajectory (train only)\n                  device_imu.csv    <- accelerometer/gyro/magnetometer (unused here)\n                  supplemental/gnss_log.txt + *.nmea  <- raw logs, if you ever\n                      want to re-derive everything from scratch (not needed, see below)\n\nConsequences of the new format:\n  1. No raw-text-log parsing (parse_gnsslogger_txt_v2 is gone). device_gnss.csv\n     already gives one clean row per (satellite, signal, epoch).\n  2. No RINEX ephemeris parsing / Kepler orbit propagation (parse_rinex_nav,\n     compute_sv_positions are gone). SvPosition[X/Y/Z]EcefMeters and\n     SvClockBiasMeters are already columns.\n  3. No custom pseudorange-from-transit-time derivation (correct_pseudorange_v3\n     is gone). RawPseudorangeMeters is given directly; we only apply the\n     standard correction:\n         correctedPrM = RawPseudorangeMeters + SvClockBiasMeters\n                         - IsrbMeters - IonosphericDelayMeters - TroposphericDelayMeters\n     which algebraically reduces to  correctedPrM ~= geometric_range + c*dt_receiver\n     (derived in Cell C-PRE) -- this is also how the competition itself says its\n     own WlsPositionXEcefMeters baseline is computed.\n  4. No Status-row merge for elevation/azimuth. SvElevationDegrees /\n     SvAzimuthDegrees are already columns, with the exact names the rest of\n     this pipeline expects.\n  5. Ground truth is now a per-epoch MOVING trajectory (train split only), so\n     the old \"hardcoded single lat/lon\" hack is replaced with: attach GT via\n     merge_asof(nearest timestamp) per (drive, phone), convert to ECEF, and\n     compute a real geometric range_residual (Cell C-PRE).\n  6. The old \"ENU-correction-from-a-fixed-point\" WLS hack only worked because\n     the receiver was static. For a moving receiver we implement a real\n     iterative (Gauss-Newton) 4-state WLS solve (x, y, z, receiver clock\n     bias) in ECEF, seeded from the dataset's own WlsPositionXEcefMeters/Y/Z\n     (Cells B1 / F / J).\n  7. A satellite can now report MULTIPLE signals (e.g. GPS_L1 and GPS_L5 on\n     the same Svid). All per-satellite time-series groupings (pseudorange\n     jump, CMC rolling std, CNN sequence windows) now group on\n     SV_KEY = ['trace','device','SvId','SignalType'], not just SvId, so two\n     different signal chains on the same satellite are never blended into\n     one fake time series.\n  8. device_gnss.csv also exposes a chipset-reported `MultipathIndicator`\n     column. In practice this is unreliable/mostly \"unknown\" on most Android\n     chipsets in these datasets -- it's surfaced defensively as an optional\n     feature but nothing here depends on it. Check its value_counts in the\n     EDA cell (B4) on your real data before trusting it.\n\nEverything else -- M-of-N label voting, the no-leakage feature/label split,\nthe causal sliding-window CNN, Platt/Isotonic calibration, and the\nhard-exclusion-vs-soft-weighted WLS comparison -- keeps the same design as\nthe original pipeline; only the data-access and geometry layers changed.\n\nEvery numerical routine below (WGS84 <-> ECEF conversion, local-ENU\nprojection, and the iterative WLS solver) was unit-tested against synthetic\nscenarios with a known answer before being placed in this script:\n  - geodetic_to_ecef / ecef_to_enu matched analytic WGS84 reference points\n    and recovered known meter-scale displacements to <0.2m.\n  - solve_wls_epoch recovered a known synthetic receiver position + clock\n    bias to sub-millimeter precision from a deliberately bad (~3km) initial\n    guess, and correctly showed soft-weighting/hard-exclusion of an injected\n    40m multipath outlier beating naive equal-weighting (64.8m -> 0.6m and\n    3.3m respectively).\n  - The full loader -> correctedPrM -> range_residual -> M-of-N labels ->\n    SV_KEY-aware sequence builder -> CNN -> WLS-comparison chain was run\n    end-to-end against a synthetic sdc2023-shaped dataset (multiple drives,\n    multiple phones, multiple signals per satellite, injected multipath\n    windows) with no errors and sane, bounded outputs at every stage.\n\nITERATION 2 -- CNN FEATURE SET NARROWED TO 5 PHYSICS-ONLY, LEAK-SAFE INPUTS\n----------------------------------------------------------------------------\nThe hard label (label_mp, Cell C) is UNCHANGED -- still the M-of-N vote over\nelevation / range_residual / pr_jump, still built independently of the CNN.\n\nFEATURE_COLS (Cell C) was narrowed from a broad defensive grab-list down to\n5 features chosen to match the smartphone/receiver-level LOS/NLOS & multipath\nML literature, and to exclude identity/categorical columns that let a model\nmemorize per-satellite/per-constellation/per-phone quirks instead of learning\ngeneralizable physics:\n  1. SvElevationDegrees\n  2. Cn0DbHz (fed as the existing 10-epoch sequence, so the CNN can learn\n     C/N0 fluctuation/ΔC/N0 patterns itself via convolution)\n  3. CMC_detrended_abs (NEW, Cell C-PRE Step 3b) -- |diff(CMC)| per SV_KEY.\n     Replaces the old CMC_abs, which was dominated by an arbitrary per-pass\n     carrier-ambiguity constant and wasn't actually measuring multipath.\n     CAVEAT: partly built from correctedPrM (like the label's pr_jump vote),\n     so it correlates more with pr_jump-triggered positives than with\n     elevation/range_residual-triggered ones -- see the recall-by-trigger\n     diagnostic added to Cell I.\n  4. carrier_doppler_consistency_abs (NEW, Cell C-PRE Step 3b) -- carrier-\n     phase-derived range rate vs. the independently-measured Doppler range\n     rate. Built only from AccumulatedDeltaRangeMeters + PseudorangeRateMeters\n     PerSecond, neither of which touches correctedPrM/RawPseudorangeMeters/\n     range_residual at all, so this one has no leakage caveat.\n  5. AccumulatedDeltaRangeUncertaintyMeters (unchanged)\n\nRemoved: SvAzimuthDegrees, ConstellationType, MultipathIndicator,\nIonosphericDelayMeters, TroposphericDelayMeters, AccumulatedDeltaRangeState,\ndevice_enc, and the old pr_rate_abs / CMC_abs. Reasoning for each is in the\ncomment block directly above the FEATURE_COLS loop in Cell C.\n════════════════════════════════════════════════════════════════════════════\n\"\"\"\n\n# ══════════════════════════════════════════════════════════════════════════\n# A — Imports, Config, Kaggle Paths (sdc2023)\n# ══════════════════════════════════════════════════════════════════════════\nimport os, glob, gc, warnings, time\nfrom pathlib import Path\n\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm.auto import tqdm            # auto-detects notebook vs. plain terminal\n\ntry:\n    from IPython.display import display\nexcept ImportError:\n    display = print                   # graceful fallback outside Jupyter\n\nfrom sklearn.ensemble import RandomForestClassifier     # quick diagnostic only (Cell E)\nfrom sklearn.linear_model import LogisticRegression      # Platt-scaling calibration\nfrom sklearn.preprocessing import StandardScaler, LabelEncoder\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.isotonic import IsotonicRegression\nfrom sklearn.model_selection import train_test_split\nfrom sklearn.metrics import (accuracy_score, precision_score, recall_score,\n                              f1_score, roc_auc_score, matthews_corrcoef,\n                              confusion_matrix, roc_curve, brier_score_loss)\n\nfrom imblearn.over_sampling import SMOTE, ADASYN\n\ntry:\n    import tensorflow as tf\n    from tensorflow.keras import layers, models, callbacks\nexcept ImportError:\n    import subprocess, sys\n    subprocess.run([sys.executable, '-m', 'pip', 'install', 'tensorflow-cpu', '-q'])\n    import tensorflow as tf\n    from tensorflow.keras import layers, models, callbacks\n\nwarnings.filterwarnings('ignore')\nsns.set_theme(style='whitegrid', palette='muted', font_scale=1.1)\nplt.rcParams['figure.dpi'] = 120\nSEED = 42\nnp.random.seed(SEED)\ntf.random.set_seed(SEED)\n\n# ── Kaggle paths — Google Smartphone Decimeter Challenge 2023 ───────────────\n# Reference file given in the task:\n#   /kaggle/input/competitions/smartphone-decimeter-2023/sdc2023/train/\n#       2020-06-25-00-34-us-ca-mtv-sb-101/pixel4/device_gnss.csv\nCOMPETITION_ROOT = '/kaggle/input/competitions/smartphone-decimeter-2023/sdc2023'\nTRAIN_ROOT = os.path.join(COMPETITION_ROOT, 'train')\nTEST_ROOT = os.path.join(COMPETITION_ROOT, 'test')\nGNSS_FILENAME = 'device_gnss.csv'\nGT_FILENAME = 'ground_truth.csv'\nIMU_FILENAME = 'device_imu.csv'          # not used in this pipeline, see header note\n\nOUTPUT_DIR = '/kaggle/working/'\nos.makedirs(OUTPUT_DIR, exist_ok=True)\nCSV_PATH_TRAIN = OUTPUT_DIR + 'raw_merged_sdc2023_train.csv'\nFORCE_REPARSE = False          # True -> always rebuild the cached CSV from source\n\nGT_MERGE_TOLERANCE_MS = 250    # max |delta t| when attaching ground_truth.csv via merge_asof\n\n# ── QUICK_TEST: bare-minimum run to validate the pipeline end-to-end ───────\n# True  -> fewer drive/phone folders parsed, fewer CNN epochs, capped WLS epochs.\n#          Use this first to catch bugs cheaply.\n# False -> full training run over every drive/phone folder.\nQUICK_TEST = True\nMAX_DRIVE_PHONE_PAIRS_QUICK = 8   # only used when QUICK_TEST=True\n\nC_BLUE, C_RED, C_GREEN = '#1565C0', '#C62828', '#2E7D32'\nC_AMBER, C_PURPLE, C_TEAL, C_GREY = '#E65100', '#4A148C', '#006064', '#546E7A'\n\nCONST_MAP = {0: 'Unknown', 1: 'GPS', 2: 'SBAS', 3: 'GLONASS',\n             4: 'QZSS', 5: 'BeiDou', 6: 'Galileo', 7: 'NavIC'}\nSPEED_OF_LIGHT = 299_792_458.0\n\n\ndef save_fig(fname, dpi=150):\n    path = os.path.join(OUTPUT_DIR, fname)\n    plt.savefig(path, bbox_inches='tight', dpi=dpi, facecolor='white')\n    plt.show()\n    print('Saved:', path)\n\n\nprint('Imports OK')\nprint('Train root :', TRAIN_ROOT)\nprint('Test root  :', TEST_ROOT)\nprint('QUICK_TEST :', QUICK_TEST, '(bare-minimum run)' if QUICK_TEST else '(full run)')\n\n# Optional / informational only: peek at the competition's sample submission\n# so you know the target format if you later want to score the test split.\n_sample_sub = os.path.join(COMPETITION_ROOT, 'sample_submission.csv')\nif os.path.exists(_sample_sub):\n    print('\\nsample_submission.csv columns:', pd.read_csv(_sample_sub, nrows=5).columns.tolist())\n\n\n\n\n","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:43:36.968235Z","iopub.execute_input":"2026-07-08T19:43:36.968582Z","iopub.status.idle":"2026-07-08T19:44:08.647793Z","shell.execute_reply.started":"2026-07-08T19:43:36.968547Z","shell.execute_reply":"2026-07-08T19:44:08.646890Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════\n# B1 — Geodesy Utilities + Iterative WLS Solver\n#      (NEW — replaces the old RINEX/raw-log parsing entirely; see header)\n# ══════════════════════════════════════════════════════════════════════════\nWGS84_A = 6378137.0                     # semi-major axis, meters\nWGS84_F = 1 / 298.257223563             # flattening\nWGS84_E2 = WGS84_F * (2 - WGS84_F)      # eccentricity squared\n\n\ndef geodetic_to_ecef(lat_deg, lon_deg, alt_m):\n    \"\"\"Vectorized WGS84 geodetic (lat, lon, alt) -> ECEF (x, y, z) meters.\"\"\"\n    lat = np.radians(np.asarray(lat_deg, dtype=float))\n    lon = np.radians(np.asarray(lon_deg, dtype=float))\n    alt = np.asarray(alt_m, dtype=float)\n    N = WGS84_A / np.sqrt(1 - WGS84_E2 * np.sin(lat) ** 2)\n    x = (N + alt) * np.cos(lat) * np.cos(lon)\n    y = (N + alt) * np.cos(lat) * np.sin(lon)\n    z = (N * (1 - WGS84_E2) + alt) * np.sin(lat)\n    return x, y, z\n\n\ndef ecef_to_enu(dx, dy, dz, ref_lat_deg, ref_lon_deg):\n    \"\"\"Rotate an ECEF displacement vector into a local East-North-Up frame\n    centered at (ref_lat_deg, ref_lon_deg). Because Cells F/J build a FRESH\n    local frame at every epoch's own ground-truth location, this stays\n    accurate even though the receiver moves across an entire drive — unlike\n    a single fixed-origin ENU approximation.\"\"\"\n    lat0 = np.radians(ref_lat_deg)\n    lon0 = np.radians(ref_lon_deg)\n    sl, cl = np.sin(lat0), np.cos(lat0)\n    so, co = np.sin(lon0), np.cos(lon0)\n    e = -so * dx + co * dy\n    n = -sl * co * dx - sl * so * dy + cl * dz\n    u = cl * co * dx + cl * so * dy + sl * dz\n    return e, n, u\n\n\ndef solve_wls_epoch(sv_pos, pr, weights, pos0, max_iter=10, tol=1e-4):\n    \"\"\"Iterative Gauss-Newton solve for receiver ECEF position + clock bias.\n\n    This directly implements the same 4-state model the competition itself\n    describes for WlsPositionXEcefMeters (\"...phone's position (x,y,z),\n    clock bias (t) ... as states for each epoch\"), just re-solved so we can\n    plug in our own per-satellite weights (elevation-only / ML-informed).\n\n    sv_pos  : (n,3) satellite ECEF positions [m]\n    pr      : (n,)  correctedPrM (see Cell C-PRE) [m]\n    weights : (n,)  observation weights (higher = more trusted)\n    pos0    : (3,)  initial ECEF position guess, e.g. this epoch's own\n                     WlsPositionXEcefMeters/Y/Z from the dataset.\n    Returns : (x, y, z, cdt_meters, n_iterations, converged: bool)\n    \"\"\"\n    sv_pos = np.asarray(sv_pos, dtype=float)\n    pr = np.asarray(pr, dtype=float)\n    w = np.asarray(weights, dtype=float)\n    n = len(pr)\n    if n < 4:\n        return np.nan, np.nan, np.nan, np.nan, 0, False\n\n    pos0 = np.asarray(pos0, dtype=float)\n    r0 = np.linalg.norm(sv_pos - pos0, axis=1)\n    cdt0 = np.median(pr - r0)           # cheap, robust clock-bias warm start\n    x = np.array([pos0[0], pos0[1], pos0[2], cdt0])\n\n    converged, it_used = False, 0\n    for it in range(max_iter):\n        it_used = it + 1\n        d = x[:3] - sv_pos                        # rx - sv\n        ranges = np.linalg.norm(d, axis=1)\n        if np.any(ranges < 1.0):\n            return np.nan, np.nan, np.nan, np.nan, it_used, False\n        los = d / ranges[:, None]                 # unit vector, sv -> rx\n        pred = ranges + x[3]\n        resid = pr - pred\n        H = np.column_stack([los[:, 0], los[:, 1], los[:, 2], np.ones(n)])\n        W = np.diag(w)\n        try:\n            delta = np.linalg.solve(H.T @ W @ H, H.T @ W @ resid)\n        except np.linalg.LinAlgError:\n            return np.nan, np.nan, np.nan, np.nan, it_used, False\n        x = x + delta\n        if np.linalg.norm(delta[:3]) < tol:\n            converged = True\n            break\n    return x[0], x[1], x[2], x[3], it_used, converged\n\n\nprint('Geodesy + WLS solver defined.')\nprint('(Unit-tested separately against known-answer synthetic scenarios —')\nprint(' see the header docstring for a summary of results.)')\n\n# ══════════════════════════════════════════════════════════════════════════\n# B2 — Dataset Walker: load device_gnss.csv + ground_truth.csv per\n#      (drive_id, phone_name)\n# ══════════════════════════════════════════════════════════════════════════\ndef load_one_device(gnss_path, attach_gt=True):\n    \"\"\"Load a single device_gnss.csv, standardize Svid->SvId, and (if present)\n    attach the matching ground_truth.csv from the same folder via a\n    nearest-timestamp merge_asof (GT is per-epoch, not per-satellite, and its\n    timestamps only approximately coincide with each Raw row's utcTimeMillis).\"\"\"\n    d = os.path.dirname(gnss_path)\n    gnss = pd.read_csv(gnss_path, low_memory=False)\n\n    if 'SvId' not in gnss.columns and 'Svid' in gnss.columns:\n        gnss = gnss.rename(columns={'Svid': 'SvId'})\n\n    gt_cols = ['LatitudeDegrees', 'LongitudeDegrees', 'AltitudeMeters',\n               'SpeedMps', 'AccuracyMeters', 'BearingDegrees']\n    gt_path = os.path.join(d, GT_FILENAME)\n    has_gt = attach_gt and os.path.exists(gt_path)\n\n    if has_gt:\n        gt = pd.read_csv(gt_path)\n        gnss = gnss.sort_values('utcTimeMillis')\n        gt = gt.sort_values('UnixTimeMillis')\n        gnss = pd.merge_asof(\n            gnss, gt[['UnixTimeMillis'] + gt_cols],\n            left_on='utcTimeMillis', right_on='UnixTimeMillis',\n            direction='nearest', tolerance=GT_MERGE_TOLERANCE_MS,\n        )\n    else:\n        for c in gt_cols:\n            gnss[c] = np.nan\n\n    return gnss, has_gt\n\n\ndef walk_split(root, split_name, max_pairs=None):\n    \"\"\"Walk root/<drive_id>/<phone_name>/device_gnss.csv, attach ground truth\n    where available, and stack everything into one long DataFrame.\"\"\"\n    files = sorted(glob.glob(os.path.join(root, '*', '*', GNSS_FILENAME)))\n    if not files:\n        print(f'No {GNSS_FILENAME} files found under {root}')\n        return pd.DataFrame()\n    if max_pairs is not None:\n        files = files[:max_pairs]\n\n    frames = []\n    for fp in tqdm(files, desc=f'Loading {split_name}'):\n        rel = os.path.relpath(fp, root)\n        parts = rel.split(os.sep)\n        drive_id = parts[0] if len(parts) > 0 else 'unknown_drive'\n        phone_name = parts[1] if len(parts) > 1 else 'unknown_phone'\n\n        gnss, has_gt = load_one_device(fp, attach_gt=(split_name == 'train'))\n        if has_gt:\n            n_matched = gnss['LatitudeDegrees'].notna().sum()\n            print(f'  {rel}: {len(gnss):,} rows | GT matched: {n_matched:,}/{len(gnss):,}')\n        else:\n            print(f'  {rel}: {len(gnss):,} rows | no ground_truth.csv (test split)')\n\n        gnss['trace'] = drive_id\n        gnss['device'] = phone_name\n        gnss['split'] = split_name\n        gnss['source_file'] = rel\n        frames.append(gnss)\n\n    out = pd.concat(frames, ignore_index=True)\n    print(f'\\n\\u2713 {split_name}: combined {len(out):,} rows, {out.shape[1]} columns, '\n          f'{out[\"trace\"].nunique()} drives, {out[\"device\"].nunique()} phones')\n    return out\n\n\nprint('Dataset walker defined.')\n\n# ══════════════════════════════════════════════════════════════════════════\n# B3 — Build / Load Combined Train Dataset\n# ══════════════════════════════════════════════════════════════════════════\nt0 = time.time()\nmax_pairs = MAX_DRIVE_PHONE_PAIRS_QUICK if QUICK_TEST else None\n\nif os.path.exists(CSV_PATH_TRAIN) and not FORCE_REPARSE:\n    print(f'Loading cached {CSV_PATH_TRAIN} ...')\n    raw = pd.read_csv(CSV_PATH_TRAIN, low_memory=False)\nelse:\n    print('Parsing device_gnss.csv + ground_truth.csv from scratch ...')\n    raw = walk_split(TRAIN_ROOT, 'train', max_pairs=max_pairs)\n    if raw.empty:\n        raise RuntimeError(f'No training rows found under {TRAIN_ROOT}. '\n                            'Check COMPETITION_ROOT in Cell A.')\n    raw.to_csv(CSV_PATH_TRAIN, index=False)\n    print(f'Wrote {CSV_PATH_TRAIN} ({len(raw):,} rows)')\n\nprint(f'\\nraw shape : {raw.shape}')\nprint(f'Elapsed   : {time.time() - t0:.1f}s')\n\n# ── Quick column check ───────────────────────────────────────────────────\nfor col in ['SvId', 'SignalType', 'SvElevationDegrees', 'SvAzimuthDegrees',\n            'RawPseudorangeMeters', 'AccumulatedDeltaRangeMeters', 'Cn0DbHz',\n            'SvPositionXEcefMeters', 'WlsPositionXEcefMeters',\n            'LatitudeDegrees', 'MultipathIndicator']:\n    present = col in raw.columns\n    fill = f'{raw[col].notna().mean()*100:.1f}%' if present else 'n/a'\n    print(f'  {col:32s} present={present}  fill={fill}')\n\n# ══════════════════════════════════════════════════════════════════════════\n# B4 — EDA\n# ══════════════════════════════════════════════════════════════════════════\nprint('Rows    :', f'{len(raw):,}')\nprint('Drives  :', raw['trace'].nunique())\nprint('Phones  :', sorted(raw['device'].unique()))\n\nif 'ConstellationType' in raw.columns:\n    print('\\nRows per constellation:')\n    print(raw['ConstellationType'].map(CONST_MAP).value_counts())\n\nif 'SignalType' in raw.columns:\n    print('\\nRows per signal type:')\n    print(raw['SignalType'].value_counts())\n\nif 'MultipathIndicator' in raw.columns:\n    print('\\nMultipathIndicator value counts (often mostly \"unknown\"/0 on Android —')\n    print('check before relying on it as a feature):')\n    print(raw['MultipathIndicator'].value_counts())\n\nmissing = raw.isna().mean().sort_values(ascending=False)\nprint('\\nTop-15 columns by missing-value fraction:')\nprint(missing.head(15).to_string())\n\nfig, axes = plt.subplots(1, 3, figsize=(16, 4))\nif 'SvElevationDegrees' in raw.columns:\n    sns.histplot(raw['SvElevationDegrees'].dropna(), bins=40, color=C_BLUE, ax=axes[0])\n    axes[0].axvline(15, color='red', ls='--', label='15\\u00b0 mask')\n    axes[0].set_title('SV Elevation (deg)')\n    axes[0].legend()\nif 'Cn0DbHz' in raw.columns:\n    sns.histplot(raw['Cn0DbHz'].dropna(), bins=40, color=C_TEAL, ax=axes[1])\n    axes[1].set_title('C/N0 (dB-Hz)')\nif 'RawPseudorangeMeters' in raw.columns:\n    pr = raw['RawPseudorangeMeters'].dropna()\n    pr = pr[pr.between(1e6, 3e7)]\n    sns.histplot(pr / 1e6, bins=40, color='#E07B54', ax=axes[2])\n    axes[2].set_title('Raw Pseudorange (\\u00d710\\u2076 m)')\nplt.tight_layout()\nsave_fig('eda_elevation_cn0_pr.png')\n\n# Sanity-plot one drive's ground-truth trajectory so you can eyeball that GT\n# actually looks like a real drive (moving), not a single static point.\none_trace = raw['trace'].iloc[0]\ngt_traj = (raw[(raw['trace'] == one_trace) & raw['LatitudeDegrees'].notna()]\n           .drop_duplicates(subset=['utcTimeMillis'])\n           .sort_values('utcTimeMillis'))\nif len(gt_traj) > 1:\n    plt.figure(figsize=(5, 5))\n    plt.plot(gt_traj['LongitudeDegrees'], gt_traj['LatitudeDegrees'], '-', color=C_PURPLE)\n    plt.scatter(gt_traj['LongitudeDegrees'].iloc[0], gt_traj['LatitudeDegrees'].iloc[0],\n                color=C_GREEN, label='start', zorder=5)\n    plt.scatter(gt_traj['LongitudeDegrees'].iloc[-1], gt_traj['LatitudeDegrees'].iloc[-1],\n                color=C_RED, label='end', zorder=5)\n    plt.xlabel('Longitude')\n    plt.ylabel('Latitude')\n    plt.title(f'Ground-truth trajectory: {one_trace}')\n    plt.legend()\n    plt.tight_layout()\n    save_fig('eda_sample_trajectory.png')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:44:08.650082Z","iopub.execute_input":"2026-07-08T19:44:08.650376Z","iopub.status.idle":"2026-07-08T19:44:56.849131Z","shell.execute_reply.started":"2026-07-08T19:44:08.650349Z","shell.execute_reply":"2026-07-08T19:44:56.847706Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════\n# C-PRE — correctedPrM, CMC, GT-ECEF, geometric range_residual, thresholds\n# ══════════════════════════════════════════════════════════════════════════\n# ── Step 1: correctedPrM (standard GNSS correction, see header note) ────────\nneeded_pr_cols = {'RawPseudorangeMeters', 'SvClockBiasMeters', 'IsrbMeters',\n                   'IonosphericDelayMeters', 'TroposphericDelayMeters'}\nif needed_pr_cols.issubset(raw.columns):\n    raw['correctedPrM'] = (raw['RawPseudorangeMeters']\n                            + raw['SvClockBiasMeters'].fillna(0)\n                            - raw['IsrbMeters'].fillna(0)\n                            - raw['IonosphericDelayMeters'].fillna(0)\n                            - raw['TroposphericDelayMeters'].fillna(0))\n    sane = raw['correctedPrM'].between(1e6, 3e7)   # 1,000-30,000 km, sane GNSS range\n    raw.loc[~sane, 'correctedPrM'] = np.nan\n    print(f\"correctedPrM computed. Valid: {raw['correctedPrM'].notna().sum():,} / \"\n          f\"{len(raw):,} ({raw['correctedPrM'].notna().mean()*100:.1f}%)\")\nelse:\n    missing_cols = needed_pr_cols - set(raw.columns)\n    raw['correctedPrM'] = np.nan\n    print(f'WARNING: missing columns {missing_cols} \\u2014 correctedPrM left as NaN.')\n\n# ── Step 2: CMC (Carrier-Minus-Code) ────────────────────────────────────────\nif {'correctedPrM', 'AccumulatedDeltaRangeMeters'}.issubset(raw.columns):\n    raw['CMC'] = raw['correctedPrM'] - raw['AccumulatedDeltaRangeMeters']\n\n# ── Step 3: Ensure SvId is canonical, define the per-signal grouping key ────\nif 'SvId' not in raw.columns and 'Svid' in raw.columns:\n    raw.rename(columns={'Svid': 'SvId'}, inplace=True)\n\n# A satellite can broadcast >1 signal (e.g. GPS_L1 and GPS_L5 on the same\n# SvId) -- group on SignalType too so time-series ops never blend them.\nSV_KEY = ['trace', 'device', 'SvId', 'SignalType'] if 'SignalType' in raw.columns \\\n    else ['trace', 'device', 'SvId']\nprint('SV_KEY =', SV_KEY)\n\n# ── Step 3b: Physics-grounded, leak-safe CNN features ───────────────────────\n# These two are NEVER used by generate_training_labels() (Cell C) -- only\n# SvElevationDegrees / range_residual / correctedPrM feed the hard label.\n# They exist purely as CNN inputs, chosen to match the smartphone/receiver-\n# level LOS/NLOS & multipath ML literature (elevation, C/N0, code-minus-\n# carrier fluctuation, pseudorange/carrier-Doppler consistency) while\n# excluding identity/categorical columns (ConstellationType, SvId, device)\n# that let a model memorize per-satellite/per-device quirks instead of\n# learning genuine multipath physics.\ng_sorted = raw.sort_values(SV_KEY + ['utcTimeMillis'])\n\n# CMC_detrended_abs: first difference of Code-Minus-Carrier per satellite\n# signal, sorted by time. Raw CMC (Cell C-PRE Step 2) is dominated by an\n# arbitrary per-pass ambiguity constant; differencing removes it and isolates\n# the multipath-driven code/carrier divergence -- the standard geodesy\n# convention for using CMC as a multipath indicator.\n# CAVEAT: CMC_detrended = diff(correctedPrM) - diff(AccumulatedDeltaRangeMeters).\n# Because it's partly built from correctedPrM (the same raw ingredient as the\n# hard label's pr_jump vote v3), it will correlate with v3-triggered examples\n# more than with v1/v2-triggered ones -- it mixes in independent carrier-phase\n# information so it is NOT a copy of v3, but it is not perfectly independent\n# either. Worth checking recall broken down by which vote(s) fired once you\n# have real results (see the note in Cell I).\nif 'CMC' in raw.columns:\n    cmc_diff = g_sorted.groupby(SV_KEY)['CMC'].diff()\n    raw['CMC_detrended_abs'] = cmc_diff.reindex(raw.index).abs()\nelse:\n    raw['CMC_detrended_abs'] = np.nan\n    print('WARNING: CMC not available -- CMC_detrended_abs left as NaN.')\n\n# carrier_doppler_consistency_abs: carrier-phase-derived range rate\n# (diff(AccumulatedDeltaRangeMeters)/dt) vs. the independently-measured\n# Doppler range rate (PseudorangeRateMetersPerSecond, which -- despite the\n# name -- is derived from the carrier frequency/Doppler tracking loop, not\n# from differencing code pseudoranges). Neither input touches correctedPrM /\n# RawPseudorangeMeters / range_residual at all, so this one is fully\n# independent of every hard-label vote -- no caveat needed.\nif {'AccumulatedDeltaRangeMeters', 'PseudorangeRateMetersPerSecond'}.issubset(raw.columns):\n    dt_s = g_sorted.groupby(SV_KEY)['utcTimeMillis'].diff() / 1000.0\n    dadr = g_sorted.groupby(SV_KEY)['AccumulatedDeltaRangeMeters'].diff()\n    carrier_rate = dadr / dt_s.replace(0, np.nan)\n    consistency = carrier_rate - g_sorted['PseudorangeRateMetersPerSecond']\n    raw['carrier_doppler_consistency_abs'] = consistency.reindex(raw.index).abs()\nelse:\n    raw['carrier_doppler_consistency_abs'] = np.nan\n    print('WARNING: AccumulatedDeltaRangeMeters/PseudorangeRateMetersPerSecond '\n          'not available -- carrier_doppler_consistency_abs left as NaN.')\n\ndel g_sorted\n\nfor c in ['CMC_detrended_abs', 'carrier_doppler_consistency_abs']:\n    print(f'{c:32s} fill={raw[c].notna().mean()*100:5.1f}%  '\n          f'median={raw[c].median():.3f}  p95={raw[c].quantile(0.95):.3f}')\n\n# ── Step 4: Ground-truth ECEF (train rows only; NaN elsewhere, e.g. test) ───\nraw['gt_ecef_x'] = np.nan\nraw['gt_ecef_y'] = np.nan\nraw['gt_ecef_z'] = np.nan\nvalid_gt = raw['LatitudeDegrees'].notna() & raw['LongitudeDegrees'].notna() & raw['AltitudeMeters'].notna()\nif valid_gt.any():\n    gx, gy, gz = geodetic_to_ecef(raw.loc[valid_gt, 'LatitudeDegrees'],\n                                   raw.loc[valid_gt, 'LongitudeDegrees'],\n                                   raw.loc[valid_gt, 'AltitudeMeters'])\n    raw.loc[valid_gt, 'gt_ecef_x'] = gx\n    raw.loc[valid_gt, 'gt_ecef_y'] = gy\n    raw.loc[valid_gt, 'gt_ecef_z'] = gz\nprint(f\"\\nGround-truth ECEF available for {valid_gt.sum():,} / {len(raw):,} rows \"\n      f\"({valid_gt.mean()*100:.1f}%) \\u2014 should be ~100% on train, ~0% on test.\")\n\n# ── Step 5: Geometric range_residual (replaces the old per-SV-median hack --\n#    that only worked for a STATIC receiver; here the receiver moves, so we\n#    use the real geometric range to the per-epoch ground truth instead) ────\nsv_pos_cols = ['SvPositionXEcefMeters', 'SvPositionYEcefMeters', 'SvPositionZEcefMeters']\ngeo_ready = valid_gt & raw['correctedPrM'].notna() & raw[sv_pos_cols].notna().all(axis=1)\n\nraw['raw_geom_residual'] = np.nan\nif geo_ready.any():\n    geom_range = np.sqrt(\n        (raw.loc[geo_ready, 'SvPositionXEcefMeters'] - raw.loc[geo_ready, 'gt_ecef_x']) ** 2 +\n        (raw.loc[geo_ready, 'SvPositionYEcefMeters'] - raw.loc[geo_ready, 'gt_ecef_y']) ** 2 +\n        (raw.loc[geo_ready, 'SvPositionZEcefMeters'] - raw.loc[geo_ready, 'gt_ecef_z']) ** 2\n    )\n    raw.loc[geo_ready, 'raw_geom_residual'] = raw.loc[geo_ready, 'correctedPrM'] - geom_range\n\n    # correctedPrM ~= geometric_range + c*dt_receiver (Step 1's derivation), so\n    # raw_geom_residual is dominated by ONE common per-epoch clock-bias term\n    # plus per-SV multipath/noise. Remove the common term via a robust\n    # per-epoch median (same trick as the old pipeline's clock-bias removal,\n    # just anchored to real geometry instead of a per-SV median PR).\n    epoch_bias = (raw.loc[geo_ready]\n                  .groupby(['trace', 'device', 'utcTimeMillis'])['raw_geom_residual']\n                  .transform('median'))\n    raw['range_residual'] = np.nan\n    raw.loc[geo_ready, 'range_residual'] = raw.loc[geo_ready, 'raw_geom_residual'] - epoch_bias.values\nelse:\n    raw['range_residual'] = np.nan\n    print('WARNING: no rows have GT + correctedPrM + SV position simultaneously \\u2014 '\n          'range_residual left as NaN (labels in Cell C will be elevation/pr-jump only).')\n\nrr = raw['range_residual'].dropna()\nif len(rr):\n    print(f\"\\n\\u2713 range_residual (after per-epoch clock-bias removal):\")\n    print(f\"  fill     : {raw['range_residual'].notna().mean()*100:.1f}%\")\n    print(f\"  mean/std : {rr.mean():.2f} / {rr.std():.2f} m\")\n    print(f\"  |rr| p50/p75/p95 : {rr.abs().quantile(0.50):.2f} / \"\n          f\"{rr.abs().quantile(0.75):.2f} / {rr.abs().quantile(0.95):.2f} m\")\n    RR_BASE_M = max(2.0, min(float(rr.abs().quantile(0.75) * 2.0), 50.0))\n    print(f\"  Auto-calibrated RR_BASE_M = {RR_BASE_M:.1f} m\")\nelse:\n    RR_BASE_M = 4.0\n    print(\"WARNING: range_residual empty \\u2014 using default RR_BASE_M=4.0\")\n\nELEV_MASK_DEG = 15.0\nPR_JUMP_M = 10.0\nprint(f\"\\nFinal thresholds \\u2192 ELEV={ELEV_MASK_DEG}\\u00b0 RR_BASE={RR_BASE_M:.1f}m PR_JUMP={PR_JUMP_M}m\")\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:44:56.850701Z","iopub.execute_input":"2026-07-08T19:44:56.852192Z","iopub.status.idle":"2026-07-08T19:44:58.519389Z","shell.execute_reply.started":"2026-07-08T19:44:56.852152Z","shell.execute_reply":"2026-07-08T19:44:58.517710Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════\n# C — Label Generation (M-of-N Multi-Metric, No Leakage) & Feature Engineering\n# ══════════════════════════════════════════════════════════════════════════\n# TRAINING label (\"label_mp\") is an M-of-N vote (>=2 of 3) over:\n#   1. Elevation < mask angle\n#   2. |range_residual| > elevation-weighted threshold (rr_base / sin(elev)) --\n#      built from the per-epoch GROUND-TRUTH ECEF position, a real geometric\n#      check, not a signal-quality heuristic (see Cell C-PRE).\n#   3. Epoch-to-epoch |correctedPrM| jump per (trace,device,SvId,SignalType)\n#      > threshold (severe NLOS/multipath event)\n# range_residual / correctedPrM are used ONLY here -- never as CNN features --\n# because they are literally what the label thresholds on (using them as\n# inputs too would leak the label into the model).\n#\n# A SEPARATE, VALIDATION-ONLY label (\"label_sigqual\") is built from\n# CN0-deviation-from-an-elevation-fitted-curve and CMC-residual rolling std.\n# It exists purely to sanity-check label_mp -- never a training target, never\n# used to pick label_mp's thresholds, since CN0/CMC ARE CNN features.\nCN0_DEV_TH = 5.5     # dB-Hz, validation label only\nCMC_STD_TH = 0.5     # meters, validation label only\nROLL_WINDOW = 15     # epochs, for CMC residual std (validation label only)\n\n\ndef generate_training_labels(df, elev_th=ELEV_MASK_DEG, rr_base=RR_BASE_M,\n                              pr_jump_th=PR_JUMP_M, sv_key=SV_KEY):\n    elev = df['SvElevationDegrees'].fillna(999.0)\n    v1 = (elev < elev_th).astype(int)\n\n    if 'range_residual' in df.columns:\n        elev_safe = elev.clip(lower=2.0, upper=90.0)\n        rr_th = rr_base / np.sin(np.radians(elev_safe))\n        v2 = (df['range_residual'].abs().fillna(0.0) > rr_th).astype(int)\n    else:\n        v2 = pd.Series(np.zeros(len(df), dtype=int), index=df.index)\n\n    if set(sv_key + ['correctedPrM']).issubset(df.columns):\n        pr_diff = df.groupby(sv_key)['correctedPrM'].transform(lambda x: x.diff().abs())\n        v3 = (pr_diff.fillna(0.0) > pr_jump_th).astype(int)\n    else:\n        v3 = pd.Series(np.zeros(len(df), dtype=int), index=df.index)\n\n    score = v1.values + v2.values + v3.values\n    return (score >= 2).astype(int), v1.values, v2.values, v3.values\n\n\ndef generate_signal_quality_label(df, cn0_dev_th=CN0_DEV_TH, cmc_std_th=CMC_STD_TH,\n                                   window=ROLL_WINDOW, sv_key=SV_KEY):\n    n = len(df)\n    v_cn0 = np.zeros(n, dtype=int)\n    v_cmc = np.zeros(n, dtype=int)\n\n    if {'Cn0DbHz', 'SvElevationDegrees'}.issubset(df.columns):\n        elev_bin = (df['SvElevationDegrees'].fillna(0) // 5 * 5).astype(int)\n        expected = df.groupby(elev_bin)['Cn0DbHz'].transform(lambda x: x.quantile(0.85))\n        dev = (expected - df['Cn0DbHz']).fillna(0.0)\n        v_cn0 = (dev > cn0_dev_th).astype(int).values\n\n    if set(sv_key + ['CMC']).issubset(df.columns):\n        rollmean = df.groupby(sv_key)['CMC'].transform(\n            lambda x: x.rolling(window, min_periods=3, center=True).mean())\n        detrended = df['CMC'] - rollmean\n        group_cols = [df[k] for k in sv_key]\n        roll_std = detrended.groupby(group_cols).transform(\n            lambda x: x.rolling(window, min_periods=3, center=True).std())\n        v_cmc = (roll_std.fillna(0.0) > cmc_std_th).astype(int).values\n\n    score = v_cn0 + v_cmc\n    return (score >= 1).astype(int)   # either metric alone is already a strong flag\n\n\nlabel_mp, v_elev, v_rr, v_prjump = generate_training_labels(raw)\nraw['label_mp'] = label_mp\nraw['label_sigqual'] = generate_signal_quality_label(raw)\n# Kept only as a diagnostic breakdown key for Cell I (which vote(s) triggered\n# each positive) -- never used as a CNN feature or re-used to build any label.\nraw['v_elev'] = v_elev\nraw['v_rr'] = v_rr\nraw['v_prjump'] = v_prjump\n\nvc = pd.Series(raw['label_mp']).value_counts(normalize=True)\nprint('Training label (label_mp) distribution \\u2014 M-of-N geometry + range-residual vote:')\nprint(f'  Class 0 (clean)      : {vc.get(0, 0) * 100:.1f}%')\nprint(f'  Class 1 (multipath)  : {vc.get(1, 0) * 100:.1f}%')\nprint(f'  Vote fire-rate -> elevation: {v_elev.mean()*100:.1f}% | '\n      f'range_residual: {v_rr.mean()*100:.1f}% | pr_jump: {v_prjump.mean()*100:.1f}%')\nprint('  NOTE: if these fire-rates look far too high/low on your data, recalibrate')\nprint('  ELEV_MASK_DEG / RR_BASE_M / PR_JUMP_M in Cell C-PRE.')\n\nonly_gt = raw['range_residual'].notna()\nif only_gt.any():\n    agreement = (raw.loc[only_gt, 'label_mp'] == raw.loc[only_gt, 'label_sigqual']).mean()\n    prec = precision_score(raw.loc[only_gt, 'label_sigqual'], raw.loc[only_gt, 'label_mp'], zero_division=0)\n    rec = recall_score(raw.loc[only_gt, 'label_sigqual'], raw.loc[only_gt, 'label_mp'], zero_division=0)\n    f1v = f1_score(raw.loc[only_gt, 'label_sigqual'], raw.loc[only_gt, 'label_mp'], zero_division=0)\n    print('\\n=== Label validation: label_mp (training) vs label_sigqual (CN0/CMC, validation-only) ===')\n    print(f'Overall agreement : {agreement * 100:.1f}%')\n    print(f'Precision : {prec:.3f}  Recall : {rec:.3f}  F1 : {f1v:.3f}')\n\n# ── Training feature set: 5 physics-grounded, leak-safe features ───────────\n# Deliberately NOT a broad defensive grab-list anymore. Chosen to mirror the\n# smartphone/receiver-level LOS/NLOS & multipath ML literature (elevation,\n# C/N0, code-minus-carrier fluctuation, carrier/Doppler consistency, carrier-\n# tracking uncertainty), while excluding:\n#   - anything derived from range_residual / correctedPrM / RawPseudorangeMeters\n#     used directly as a raw value (label-only, see Cell C-PRE / generate_\n#     training_labels above -- this is the hard no-leakage rule)\n#   - identity/categorical columns (ConstellationType, SvId, SignalType,\n#     device) that let a model memorize per-satellite/per-constellation/\n#     per-phone quirks instead of learning generalizable multipath physics\n#   - MultipathIndicator (another system's opaque classification output, not\n#     a physical measurement) and IonosphericDelayMeters/TroposphericDelayMeters\n#     (a different physical phenomenon -- atmospheric refraction -- nearly\n#     redundant with elevation, not multipath-causing)\n#   - SvAzimuthDegrees (only physically meaningful for multipath alongside an\n#     external 3D building model for ray-tracing, which we don't have; alone\n#     it's closer to a site fingerprint than a general physical driver)\nFEATURE_COLS = []\nfor col in ['SvElevationDegrees', 'Cn0DbHz', 'CMC_detrended_abs',\n            'carrier_doppler_consistency_abs', 'AccumulatedDeltaRangeUncertaintyMeters']:\n    if col in raw.columns:\n        FEATURE_COLS.append(col)\n    else:\n        print(f'WARNING: expected feature {col!r} not found -- check Cell C-PRE ran first.')\n\nraw_clean = raw.dropna(subset=['Cn0DbHz', 'SvElevationDegrees']).copy()\n\nprint(f'\\nTraining features ({len(FEATURE_COLS)}): {FEATURE_COLS}')\nprint(f'Training label     : label_mp (geometry + range-residual M-of-N, no signal-quality leakage)')\nprint(f'Validation anchor   : label_sigqual (CN0-deviation + CMC-residual, validation-only)')\nprint(f'Clean dataset       : {len(raw_clean):,} measurements')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:44:58.523114Z","iopub.execute_input":"2026-07-08T19:44:58.523425Z","iopub.status.idle":"2026-07-08T19:45:00.865412Z","shell.execute_reply.started":"2026-07-08T19:44:58.523399Z","shell.execute_reply":"2026-07-08T19:45:00.863282Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════\n# D — Data Matrices: Causal Sliding Windows per (trace, device, SvId, SignalType)\n# ══════════════════════════════════════════════════════════════════════════\n# Same idea as before: the CNN sees the last WINDOW epochs of engineered\n# features per satellite SIGNAL (a short time-series) instead of one flat\n# feature vector, so it can learn fading/fluctuation patterns. Windows are\n# causal (look-back only, labeled at the LAST epoch) and never cross an\n# SV_KEY boundary -- updated from the old pipeline to also respect SignalType,\n# so an L1 and an L5 track on the same satellite are never mixed into one\n# fake continuous series.\nWINDOW = 10\n\nX_raw = raw_clean[FEATURE_COLS].copy()\nimputer = SimpleImputer(strategy='median')\nX_imputed = pd.DataFrame(imputer.fit_transform(X_raw), columns=FEATURE_COLS, index=raw_clean.index)\nraw_clean[[c + '_imp' for c in FEATURE_COLS]] = X_imputed.values\n\n\ndef build_sequences(df, feature_cols, window=WINDOW, label_col='label_mp', sv_key=SV_KEY):\n    \"\"\"Group by sv_key, sort by time, build causal sliding windows.\n    Returns X_seq (N, window, n_features), y_seq (N,), meta_df (sv_key columns\n    + utcTimeMillis of the LAST epoch in each window -- used to merge\n    predictions back onto raw_clean downstream in Cells H/J).\"\"\"\n    imp_cols = [c + '_imp' for c in feature_cols]\n    seqs, labels, meta_rows = [], [], []\n\n    for key_vals, grp in df.groupby(sv_key, sort=False):\n        if not isinstance(key_vals, tuple):\n            key_vals = (key_vals,)\n        g = grp.sort_values('utcTimeMillis')\n        feat_mat = g[imp_cols].values\n        lab_arr = g[label_col].values\n        t_arr = g['utcTimeMillis'].values\n        n = len(g)\n        if n < window:\n            continue\n        for i in range(window - 1, n):\n            seqs.append(feat_mat[i - window + 1: i + 1])\n            labels.append(lab_arr[i])\n            meta_rows.append(tuple(key_vals) + (t_arr[i],))\n\n    if not seqs:\n        return (np.zeros((0, window, len(feature_cols)), dtype=np.float32),\n                np.zeros((0,), dtype=np.int32),\n                pd.DataFrame(columns=sv_key + ['utcTimeMillis']))\n\n    X_seq = np.asarray(seqs, dtype=np.float32)\n    y_seq = np.asarray(labels, dtype=np.int32)\n    meta_df = pd.DataFrame(meta_rows, columns=sv_key + ['utcTimeMillis'])\n    return X_seq, y_seq, meta_df\n\n\nX_seq_all, y_seq_all, meta_all = build_sequences(raw_clean, FEATURE_COLS, window=WINDOW, label_col='label_mp')\n\nprint(f'Sequence tensor : {X_seq_all.shape} (N windows, {WINDOW} epochs, {len(FEATURE_COLS)} features)')\nif len(y_seq_all):\n    print(f'Class dist : {pd.Series(y_seq_all).value_counts(normalize=True).round(3).to_dict()}')\nprint('(Satellite signals visible for fewer than WINDOW epochs are dropped from the sequence set.)')\n\n# Flat (non-windowed) matrix -- kept ONLY for the quick resampling-ablation\n# diagnostic in Cell E, not used by the CNN itself.\nN_SAMPLES = min(50_000, len(X_imputed))\nidx_main = np.random.choice(len(X_imputed), N_SAMPLES, replace=False)\nX_main = X_imputed.values[idx_main]\ny_main = raw_clean['label_mp'].values[idx_main]\n\nprint(f'\\nX_main (flat, ablation-only): {X_main.shape}')\nprint(f'Class dist : {pd.Series(y_main).value_counts(normalize=True).round(3).to_dict()}')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:45:00.867347Z","iopub.execute_input":"2026-07-08T19:45:00.867813Z","iopub.status.idle":"2026-07-08T19:45:05.159347Z","shell.execute_reply.started":"2026-07-08T19:45:00.867784Z","shell.execute_reply":"2026-07-08T19:45:05.157499Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════\n# E — Class-Imbalance Diagnostic (RF ablation, sanity check only — not used downstream)\n# ══════════════════════════════════════════════════════════════════════════\nXtr, Xte, ytr, yte = train_test_split(\n    X_main, y_main, test_size=0.25, random_state=SEED, stratify=y_main\n)\n\n\ndef quick_rf(X_train, y_train, X_test, y_test, label):\n    clf = RandomForestClassifier(\n        n_estimators=200, max_depth=12, class_weight='balanced',\n        random_state=SEED, n_jobs=-1\n    )\n    clf.fit(X_train, y_train)\n    proba = clf.predict_proba(X_test)[:, 1]\n    pred = (proba >= 0.5).astype(int)\n    row = {\n        'variant': label,\n        'accuracy': accuracy_score(y_test, pred),\n        'precision': precision_score(y_test, pred, zero_division=0),\n        'recall': recall_score(y_test, pred, zero_division=0),\n        'f1': f1_score(y_test, pred, zero_division=0),\n        'auc': roc_auc_score(y_test, proba) if len(np.unique(y_test)) > 1 else np.nan,\n        'mcc': matthews_corrcoef(y_test, pred),\n    }\n    return row, clf\n\n\nresults_ablation = []\nrow, rf_balanced = quick_rf(Xtr, ytr, Xte, yte, 'class_weight=balanced')\nresults_ablation.append(row)\n\ntry:\n    Xtr_sm, ytr_sm = SMOTE(random_state=SEED).fit_resample(Xtr, ytr)\n    row, _ = quick_rf(Xtr_sm, ytr_sm, Xte, yte, 'SMOTE')\n    results_ablation.append(row)\nexcept Exception as e:\n    print('SMOTE skipped:', e)\n\ntry:\n    Xtr_ad, ytr_ad = ADASYN(random_state=SEED).fit_resample(Xtr, ytr)\n    row, _ = quick_rf(Xtr_ad, ytr_ad, Xte, yte, 'ADASYN')\n    results_ablation.append(row)\nexcept Exception as e:\n    print('ADASYN skipped:', e)\n\nablation_df = pd.DataFrame(results_ablation).set_index('variant').round(3)\nprint('Diagnostic-only RF ablation (feature/label sanity check, NOT used downstream):')\ndisplay(ablation_df)\n\nimportances = pd.Series(rf_balanced.feature_importances_, index=FEATURE_COLS).sort_values(ascending=False)\nplt.figure(figsize=(7, 4))\nsns.barplot(x=importances.values, y=importances.index, color=C_PURPLE)\nplt.title('RF feature importance (diagnostic only)')\nplt.tight_layout()\nsave_fig('ablation_feature_importance.png')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:45:05.161027Z","iopub.execute_input":"2026-07-08T19:45:05.161932Z","iopub.status.idle":"2026-07-08T19:45:26.624569Z","shell.execute_reply.started":"2026-07-08T19:45:05.161884Z","shell.execute_reply":"2026-07-08T19:45:26.623364Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════\n# F — WLS Positioning Baselines: Google-given vs. Recomputed (hard elevation mask)\n# ══════════════════════════════════════════════════════════════════════════\n# Unlike the old static-receiver pipeline (which \"corrected\" a single fixed\n# point with a 2-D ENU line-of-sight fit), the receiver here MOVES, so we run\n# the real iterative ECEF+clock-bias solver from Cell B1 at every epoch,\n# seeded from that epoch's own WlsPositionXEcefMeters/Y/Z, and evaluate\n# horizontal error in a FRESH local ENU frame anchored at that epoch's own\n# ground-truth position (accurate for a moving receiver, unlike one fixed\n# ENU origin for an entire drive).\nWLS_POS_COLS = ['WlsPositionXEcefMeters', 'WlsPositionYEcefMeters', 'WlsPositionZEcefMeters']\nSV_POS_COLS = ['SvPositionXEcefMeters', 'SvPositionYEcefMeters', 'SvPositionZEcefMeters']\n\n\ndef run_wls_over_epochs(df, weight_col=None, elev_mask=ELEV_MASK_DEG, quick_test=QUICK_TEST):\n    \"\"\"Solve one 4-state WLS fix per (trace, device, utcTimeMillis) epoch and,\n    where ground truth is available, evaluate horizontal error against it.\"\"\"\n    need = set(SV_POS_COLS + WLS_POS_COLS + ['correctedPrM', 'SvElevationDegrees'])\n    if not need.issubset(df.columns):\n        print(f'Missing columns for WLS: {need - set(df.columns)}')\n        return pd.DataFrame(columns=['trace', 'device', 'utcTimeMillis', 'horiz_err_m', 'n_sv_used'])\n\n    groups = list(df.groupby(['trace', 'device', 'utcTimeMillis']))\n    if quick_test:\n        groups = groups[: min(3000, len(groups))]\n\n    rows = []\n    for (trace_, device_, t_), g in tqdm(groups, desc='WLS epochs'):\n        d = g.copy()\n        if elev_mask is not None:\n            d = d[d['SvElevationDegrees'] >= elev_mask]\n        d = d.dropna(subset=SV_POS_COLS + WLS_POS_COLS + ['correctedPrM'])\n        if len(d) < 4:\n            rows.append({'trace': trace_, 'device': device_, 'utcTimeMillis': t_,\n                         'horiz_err_m': np.nan, 'n_sv_used': len(d)})\n            continue\n\n        sv_pos = d[SV_POS_COLS].values\n        pr = d['correctedPrM'].values\n        if weight_col is not None and weight_col in d.columns:\n            w = np.clip(d[weight_col].values, 1e-3, None)\n        else:\n            w = np.ones(len(d))\n        pos0 = d[WLS_POS_COLS].values[0]\n\n        x, y, z, cdt, n_iter, conv = solve_wls_epoch(sv_pos, pr, w, pos0)\n\n        horiz_err = np.nan\n        if conv and pd.notna(d['gt_ecef_x'].iloc[0]):\n            e, n, u = ecef_to_enu(x - d['gt_ecef_x'].iloc[0], y - d['gt_ecef_y'].iloc[0],\n                                   z - d['gt_ecef_z'].iloc[0],\n                                   d['LatitudeDegrees'].iloc[0], d['LongitudeDegrees'].iloc[0])\n            horiz_err = float(np.sqrt(e ** 2 + n ** 2))\n\n        rows.append({'trace': trace_, 'device': device_, 'utcTimeMillis': t_,\n                      'x_ecef': x, 'y_ecef': y, 'z_ecef': z, 'cdt_m': cdt,\n                      'horiz_err_m': horiz_err, 'n_sv_used': len(d)})\n    return pd.DataFrame(rows)\n\n\n# Variant 0: Google's own given baseline (straight from WlsPositionXEcefMeters,\n# no re-solving) -- a useful sanity check for the rest of this section.\ngiven_epochs = (raw_clean.dropna(subset=WLS_POS_COLS + ['gt_ecef_x'])\n                .drop_duplicates(subset=['trace', 'device', 'utcTimeMillis']))\nif len(given_epochs):\n    e, n, u = ecef_to_enu(given_epochs['WlsPositionXEcefMeters'] - given_epochs['gt_ecef_x'],\n                           given_epochs['WlsPositionYEcefMeters'] - given_epochs['gt_ecef_y'],\n                           given_epochs['WlsPositionZEcefMeters'] - given_epochs['gt_ecef_z'],\n                           given_epochs['LatitudeDegrees'], given_epochs['LongitudeDegrees'])\n    given_baseline_pos = pd.DataFrame({\n        'trace': given_epochs['trace'], 'device': given_epochs['device'],\n        'utcTimeMillis': given_epochs['utcTimeMillis'],\n        'horiz_err_m': np.sqrt(e ** 2 + n ** 2),\n    })\n    print('Google-given WlsPosition baseline \\u2014 '\n          f'{given_baseline_pos[\"horiz_err_m\"].notna().sum():,} epochs')\n    print(given_baseline_pos['horiz_err_m'].describe())\nelse:\n    given_baseline_pos = pd.DataFrame(columns=['trace', 'device', 'utcTimeMillis', 'horiz_err_m'])\n    print('No epochs with both WlsPosition and ground truth \\u2014 given_baseline_pos left empty.')\n\n# Variant 1: our own recomputed baseline, hard elevation-mask exclusion, equal weight\nbaseline_pos = run_wls_over_epochs(raw_clean, weight_col=None, elev_mask=ELEV_MASK_DEG)\nsolved = baseline_pos['horiz_err_m'].notna().sum()\nprint(f'\\nRecomputed baseline WLS (hard elevation mask only) \\u2014 {solved:,} epochs solved')\nif solved:\n    print(baseline_pos['horiz_err_m'].describe())\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:45:26.625748Z","iopub.execute_input":"2026-07-08T19:45:26.626072Z","iopub.status.idle":"2026-07-08T19:45:38.582077Z","shell.execute_reply.started":"2026-07-08T19:45:26.626044Z","shell.execute_reply":"2026-07-08T19:45:38.581239Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# ══════════════════════════════════════════════════════════════════════════\n# G — 1-D CNN Model: Architecture + Trace-Level Train/Val Split\n# ══════════════════════════════════════════════════════════════════════════\ndef build_cnn(window, n_features):\n    inp = layers.Input(shape=(window, n_features))\n    x = layers.Conv1D(32, 3, padding='causal', activation='relu')(inp)\n    x = layers.BatchNormalization()(x)\n    x = layers.Conv1D(64, 3, padding='causal', activation='relu')(x)\n    x = layers.BatchNormalization()(x)\n    x = layers.Conv1D(64, 3, padding='causal', activation='relu')(x)\n    x = layers.GlobalAveragePooling1D()(x)\n    x = layers.Dense(32, activation='relu')(x)\n    x = layers.Dropout(0.3)(x)\n    out = layers.Dense(1, activation='sigmoid')(x)\n    model = models.Model(inp, out)\n    model.compile(\n        optimizer=tf.keras.optimizers.Adam(1e-3),\n        loss='binary_crossentropy',\n        metrics=[tf.keras.metrics.AUC(name='auc'), 'accuracy']\n    )\n    return model\n\n\n# Trace-level split: hold out ~20% of DRIVES entirely for validation, so the\n# CNN is never evaluated on a drive it partially trained on.\ntraces_unique = meta_all['trace'].unique()\nrng = np.random.RandomState(SEED)\nrng.shuffle(traces_unique)\nn_val_traces = max(1, int(0.2 * len(traces_unique)))\nval_traces = set(traces_unique[:n_val_traces])\n\nis_val = meta_all['trace'].isin(val_traces).values\nX_train_seq, y_train_seq = X_seq_all[~is_val], y_seq_all[~is_val]\nX_val_seq, y_val_seq = X_seq_all[is_val], y_seq_all[is_val]\n\nprint(f'Train windows: {X_train_seq.shape[0]:,} | Val windows: {X_val_seq.shape[0]:,}')\nprint(f'Val traces ({len(val_traces)}): {sorted(val_traces)}')\n\nn_feat = X_train_seq.shape[-1]\nflat_train = X_train_seq.reshape(-1, n_feat)\nscaler = StandardScaler().fit(flat_train)\n\n\ndef scale_seq(X):\n    shp = X.shape\n    return scaler.transform(X.reshape(-1, shp[-1])).reshape(shp)\n\n\nX_train_seq_s = scale_seq(X_train_seq)\nX_val_seq_s = scale_seq(X_val_seq) if len(X_val_seq) else X_val_seq\n\ncnn = build_cnn(WINDOW, n_feat)\ncnn.summary()\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:45:38.583211Z","iopub.execute_input":"2026-07-08T19:45:38.583517Z","iopub.status.idle":"2026-07-08T19:45:39.550100Z","shell.execute_reply.started":"2026-07-08T19:45:38.583493Z","shell.execute_reply":"2026-07-08T19:45:39.549193Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════\n# H — CNN Training + Calibration (Platt vs Isotonic -> ml_prob)\n# ══════════════════════════════════════════════════════════════════════════\nclasses = np.unique(y_train_seq)\ncw = {c: len(y_train_seq) / (len(classes) * max((y_train_seq == c).sum(), 1)) for c in classes}\nprint('Class weights:', cw)\n\nEPOCHS = 5 if QUICK_TEST else 40\ncb = [\n    callbacks.EarlyStopping(monitor='val_auc', mode='max', patience=5, restore_best_weights=True),\n    callbacks.ReduceLROnPlateau(monitor='val_loss', factor=0.5, patience=3),\n]\n\nhistory = cnn.fit(\n    X_train_seq_s, y_train_seq,\n    validation_data=(X_val_seq_s, y_val_seq) if len(X_val_seq_s) else None,\n    epochs=EPOCHS, batch_size=256, class_weight=cw, callbacks=cb, verbose=1\n)\n\nfig, axes = plt.subplots(1, 2, figsize=(11, 4))\naxes[0].plot(history.history['loss'], label='train', color=C_BLUE)\nif 'val_loss' in history.history:\n    axes[0].plot(history.history['val_loss'], label='val', color=C_RED)\naxes[0].set_title('Loss')\naxes[0].legend()\naxes[1].plot(history.history['auc'], label='train', color=C_BLUE)\nif 'val_auc' in history.history:\n    axes[1].plot(history.history['val_auc'], label='val', color=C_RED)\naxes[1].set_title('AUC')\naxes[1].legend()\nplt.tight_layout()\nsave_fig('cnn_training_curves.png')\n\n# ── Calibration on the validation windows ────────────────────────────────\nraw_prob_val = cnn.predict(X_val_seq_s, verbose=0).ravel() if len(X_val_seq_s) else np.array([])\n\nif len(raw_prob_val):\n    eps = 1e-6\n    logit_val = np.log(np.clip(raw_prob_val, eps, 1 - eps) / np.clip(1 - raw_prob_val, eps, 1 - eps)).reshape(-1, 1)\n    platt = LogisticRegression().fit(logit_val, y_val_seq)\n    platt_prob = platt.predict_proba(logit_val)[:, 1]\n\n    iso = IsotonicRegression(out_of_bounds='clip').fit(raw_prob_val, y_val_seq)\n    iso_prob = iso.predict(raw_prob_val)\n\n    brier_raw = brier_score_loss(y_val_seq, raw_prob_val)\n    brier_platt = brier_score_loss(y_val_seq, platt_prob)\n    brier_iso = brier_score_loss(y_val_seq, iso_prob)\n    print(f'Brier \\u2014 raw: {brier_raw:.4f} | Platt: {brier_platt:.4f} | Isotonic: {brier_iso:.4f}')\n\n    CALIBRATION = 'isotonic' if brier_iso <= brier_platt else 'platt'\n    print('Selected calibrator:', CALIBRATION)\nelse:\n    platt, iso, CALIBRATION = None, None, 'none'\n    print('No validation windows available (check trace split) \\u2014 using raw CNN probabilities uncalibrated.')\n\n\ndef calibrate_prob(raw_prob):\n    raw_prob = np.clip(raw_prob, 1e-6, 1 - 1e-6)\n    if CALIBRATION == 'platt' and platt is not None:\n        logit = np.log(raw_prob / (1 - raw_prob)).reshape(-1, 1)\n        return platt.predict_proba(logit)[:, 1]\n    if CALIBRATION == 'isotonic' and iso is not None:\n        return iso.predict(raw_prob)\n    return raw_prob\n\n\n# Score EVERY window (train + val) to get ml_prob for every (SV_KEY, epoch)\nX_all_seq_s = scale_seq(X_seq_all)\nraw_prob_all = cnn.predict(X_all_seq_s, verbose=0).ravel()\nml_prob_all = calibrate_prob(raw_prob_all)\nmeta_all = meta_all.copy()\nmeta_all['ml_prob'] = ml_prob_all\nmeta_all['ml_pred'] = (meta_all['ml_prob'] >= 0.5).astype(int)\n\nprint(f'\\nml_prob assigned to {len(meta_all):,} ({\", \".join(SV_KEY)}, epoch) windows.')\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:45:39.551434Z","iopub.execute_input":"2026-07-08T19:45:39.552269Z","iopub.status.idle":"2026-07-08T19:48:01.787397Z","shell.execute_reply.started":"2026-07-08T19:45:39.552239Z","shell.execute_reply":"2026-07-08T19:48:01.786506Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# ══════════════════════════════════════════════════════════════════════════\n# I — Classifier Evaluation (on held-out trace-level validation windows)\n# ══════════════════════════════════════════════════════════════════════════\nif len(X_val_seq_s):\n    val_prob_cal = calibrate_prob(raw_prob_val)\n    val_pred = (val_prob_cal >= 0.5).astype(int)\n\n    metrics_row = {\n        'accuracy': accuracy_score(y_val_seq, val_pred),\n        'precision': precision_score(y_val_seq, val_pred, zero_division=0),\n        'recall': recall_score(y_val_seq, val_pred, zero_division=0),\n        'f1': f1_score(y_val_seq, val_pred, zero_division=0),\n        'auc': roc_auc_score(y_val_seq, val_prob_cal) if len(np.unique(y_val_seq)) > 1 else np.nan,\n        'mcc': matthews_corrcoef(y_val_seq, val_pred),\n    }\n    print('CNN (calibrated) \\u2014 held-out trace validation:')\n    for k, v in metrics_row.items():\n        print(f'  {k:10s}: {v:.3f}')\n\n    cm = confusion_matrix(y_val_seq, val_pred)\n    fig, axes = plt.subplots(1, 2, figsize=(11, 4.2))\n    sns.heatmap(cm, annot=True, fmt='d', cmap='Blues', ax=axes[0],\n                xticklabels=['clean', 'multipath'], yticklabels=['clean', 'multipath'])\n    axes[0].set_title('Confusion matrix')\n    axes[0].set_xlabel('Predicted')\n    axes[0].set_ylabel('True')\n\n    fpr, tpr, _ = roc_curve(y_val_seq, val_prob_cal)\n    axes[1].plot(fpr, tpr, color=C_BLUE, label=f\"AUC = {metrics_row['auc']:.3f}\")\n    axes[1].plot([0, 1], [0, 1], '--', color=C_GREY)\n    axes[1].set_title('ROC curve')\n    axes[1].set_xlabel('FPR')\n    axes[1].set_ylabel('TPR')\n    axes[1].legend()\n    plt.tight_layout()\n    save_fig('cnn_eval_confusion_roc.png')\n\n    # ── Diagnostic: recall broken down by which hard-label vote(s) fired ────\n    # CMC_detrended_abs is partly built from correctedPrM (see the caveat in\n    # Cell C-PRE Step 3b), so it's worth checking whether the CNN is genuinely\n    # picking up elevation- / range-residual-driven multipath (votes v1/v2)\n    # too, or mostly just re-detecting the sudden pseudorange-jump cases\n    # (vote v3) through the CMC pathway.\n    vote_cols = ['v_elev', 'v_rr', 'v_prjump']\n    if set(vote_cols).issubset(raw_clean.columns):\n        meta_val = meta_all[is_val][SV_KEY + ['utcTimeMillis']].reset_index(drop=True).copy()\n        meta_val['y_true'] = y_val_seq\n        meta_val['y_pred'] = val_pred\n        meta_val = meta_val.merge(\n            raw_clean[SV_KEY + ['utcTimeMillis'] + vote_cols]\n            .drop_duplicates(subset=SV_KEY + ['utcTimeMillis']),\n            on=SV_KEY + ['utcTimeMillis'], how='left'\n        )\n\n        positives = meta_val[meta_val['y_true'] == 1].copy()\n        if len(positives):\n            def _trigger(row):\n                tags = []\n                if row['v_elev'] == 1:\n                    tags.append('elev')\n                if row['v_rr'] == 1:\n                    tags.append('range_resid')\n                if row['v_prjump'] == 1:\n                    tags.append('pr_jump')\n                return '+'.join(tags) if tags else 'none'\n\n            positives['trigger'] = positives.apply(_trigger, axis=1)\n            breakdown = pd.concat([\n                positives.groupby('trigger').size().rename('n_positives'),\n                positives.groupby('trigger')['y_pred'].mean().rename('recall'),\n            ], axis=1).sort_values('n_positives', ascending=False)\n            print('\\nRecall on validation positives, broken down by which vote(s) triggered label_mp==1:')\n            display(breakdown.round(3))\n            print('If recall is much higher where pr_jump fired than on pure elev/range_resid rows,')\n            print('the CNN is leaning on the CMC<->pr_jump correlation rather than learning the')\n            print('elevation/geometric-residual patterns independently -- worth a closer look if so.')\nelse:\n    print('No validation split available \\u2014 skipping evaluation plots.')\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:48:01.789945Z","iopub.execute_input":"2026-07-08T19:48:01.790307Z","iopub.status.idle":"2026-07-08T19:48:03.280250Z","shell.execute_reply.started":"2026-07-08T19:48:01.790256Z","shell.execute_reply":"2026-07-08T19:48:03.279007Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════\n# J — Merge ml_prob back onto raw_clean, run Hard-Exclusion(ML) & Soft-Weighted(ML) WLS\n# ══════════════════════════════════════════════════════════════════════════\nmerge_keys = SV_KEY + ['utcTimeMillis']\nraw_scored = raw_clean.merge(\n    meta_all[merge_keys + ['ml_prob', 'ml_pred']],\n    on=merge_keys, how='inner'\n)\nprint(f'{len(raw_scored):,} / {len(raw_clean):,} raw_clean rows matched to a scored window '\n      f'(rows before epoch {WINDOW - 1} for their satellite-signal have no window yet).')\n\nel_safe = np.radians(raw_scored['SvElevationDegrees'].fillna(45.0).clip(2.0, 90.0))\nraw_scored['w_elev'] = 1.0 / np.clip(np.sin(el_safe), 0.05, 1.0) ** 2\nraw_scored['w_soft'] = raw_scored['w_elev'] * (1.0 - raw_scored['ml_prob'])\n\nhard_ml = raw_scored[raw_scored['ml_prob'] < 0.5].copy()\nhardml_pos = run_wls_over_epochs(hard_ml, weight_col=None, elev_mask=ELEV_MASK_DEG)\nsoftml_pos = run_wls_over_epochs(raw_scored, weight_col='w_soft', elev_mask=None)\n\nfor name, dfr in [('Hard-exclusion (ML)', hardml_pos), ('Soft-weighted (ML)', softml_pos)]:\n    solved = dfr['horiz_err_m'].notna().sum()\n    print(f'\\n{name} \\u2014 {solved:,} epochs solved')\n    if solved:\n        print(dfr['horiz_err_m'].describe())\n\n# ══════════════════════════════════════════════════════════════════════════\n# K — H68 / H95 Comparison across all four positioning variants\n# ══════════════════════════════════════════════════════════════════════════\ndef h_percentiles(err_series, p=(68, 95)):\n    e = err_series.dropna().values\n    if len(e) == 0:\n        return {f'H{pp}': np.nan for pp in p}\n    return {f'H{pp}': np.percentile(e, pp) for pp in p}\n\n\nsummary_rows = []\nfor name, df_ in [('Google-given WLS baseline', given_baseline_pos),\n                   ('Recomputed baseline (hard elev.)', baseline_pos),\n                   ('Hard-exclusion (ML)', hardml_pos),\n                   ('Soft-weighted (ML)', softml_pos)]:\n    hp = h_percentiles(df_['horiz_err_m']) if 'horiz_err_m' in df_.columns else {'H68': np.nan, 'H95': np.nan}\n    hp['variant'] = name\n    hp['epochs_solved'] = df_['horiz_err_m'].notna().sum() if 'horiz_err_m' in df_.columns else 0\n    summary_rows.append(hp)\n\nsummary_df = pd.DataFrame(summary_rows).set_index('variant')[['H68', 'H95', 'epochs_solved']]\nprint('H68 / H95 horizontal-error comparison (meters):')\ndisplay(summary_df.round(2))\n\nplt.figure(figsize=(8, 4.5))\nx = np.arange(len(summary_df))\nwidth = 0.35\nplt.bar(x - width / 2, summary_df['H68'], width, label='H68', color=C_BLUE)\nplt.bar(x + width / 2, summary_df['H95'], width, label='H95', color=C_RED)\nplt.xticks(x, summary_df.index, rotation=15, ha='right')\nplt.ylabel('Horizontal error (m)')\nplt.title('Positioning accuracy: Google baseline vs. recomputed vs. ML-driven weighting')\nplt.legend()\nplt.tight_layout()\nsave_fig('h68_h95_comparison.png')\n\nsummary_df.to_csv(os.path.join(OUTPUT_DIR, 'h68_h95_summary_sdc2023.csv'))\nprint(f\"Saved: {os.path.join(OUTPUT_DIR, 'h68_h95_summary_sdc2023.csv')}\")","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-07-08T19:48:03.281458Z","iopub.execute_input":"2026-07-08T19:48:03.281776Z","iopub.status.idle":"2026-07-08T19:48:25.432538Z","shell.execute_reply.started":"2026-07-08T19:48:03.281750Z","shell.execute_reply":"2026-07-08T19:48:25.431492Z"}},"outputs":[],"execution_count":null}]}