{"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":"# # This Python 3 environment comes with many helpful analytics libraries installed\n# # It is defined by the kaggle/python Docker image: https://github.com/kaggle/docker-python\n# # For example, here's several helpful packages to load\n\n# import numpy as np # linear algebra\n# import pandas as pd # data processing, CSV file I/O (e.g. pd.read_csv)\n\n# # Input data files are available in the read-only \"../input/\" directory\n# # For example, running this (by clicking run or pressing Shift+Enter) will list all files under the input directory\n\n# import os\n# for dirname, _, filenames in os.walk('/kaggle/input'):\n#     for filename in filenames:\n#         print(os.path.join(dirname, filename))\n\n# # You can write up to 20GB to the current directory (/kaggle/working/) that gets preserved as output when you create a version using \"Save & Run All\" \n# # You can also write temporary files to /kaggle/temp/, but they won't be saved outside of the current session\n\n# # Use the kagglehub client library to attach Kaggle resources like competitions, datasets, and models to your session\n# # Learn more about kagglehub: https://github.com/Kaggle/kagglehub/blob/main/README.md\n\n# import kagglehub\n# # kagglehub.dataset_download('<owner>/<dataset-slug>')","metadata":{"_uuid":"8f2839f25d086af736a60e9eeb907d3b93b6e0e5","_cell_guid":"b1076dfc-b9ad-4769-8c92-a6c4dae69d19","trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:12.253509Z","iopub.execute_input":"2026-06-05T22:04:12.254086Z","iopub.status.idle":"2026-06-05T22:04:12.258994Z","shell.execute_reply.started":"2026-06-05T22:04:12.254056Z","shell.execute_reply":"2026-06-05T22:04:12.258109Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# -*- coding: utf-8 -*-\n\"\"\"multipath_rerun.ipynb\n\nAutomatically generated by Colab.\n\nOriginal file is located at\n    https://colab.research.google.com/drive/1Z5dit6XZzpMqy4qkctrtHmeHYoWu1k0D\n\n# Minimal Re-Run — Positioning & Calibration Fixes\n\n**What this notebook does:** Re-runs ONLY the new and corrected cells.\nAll heavy steps (data loading from raw files, LODO, SHAP, 5-fold classifier training) are skipped — replaced by loading `raw_merged.csv` which your full notebook already saved.\n\n**Estimated runtime: 5–8 minutes** (vs 30–60 min for the full notebook)\n\n## Execution order\n| Cell | Content | Time |\n|------|---------|------|\n| A | Imports & paths | <5 s |\n| B | Load raw_merged.csv | <30 s |\n| C | Feature engineering & labels | <30 s |\n| D | Data matrices | <15 s |\n| E | get_balanced_data helper | <5 s |\n| F | **WLS helpers (CORRECTED)** | <5 s |\n| G | **Positioning experiment (CORRECTED)** | ~1–2 min |\n| H | Positioning visualisation | <30 s |\n| I | **Platt calibration (NEW)** | ~2–3 min |\n| J | **Weighting function defs + curves (NEW)** | <30 s |\n| K | **Full hard/soft WLS comparison (NEW)** | ~2–3 min |\n| L | **WLS visualisation (NEW)** | <30 s |\n\n> All figures are saved to your `FinalMultiPathOutput` folder automatically.\n\"\"\"\n\n# from google.colab import drive\n# drive.mount('/content/drive')\n\n\"\"\"## A — Imports & Paths\"\"\"\n\n# ══════════════════════════════════════════════════════════════════════════════\n# MINIMAL RE-RUN SCRIPT — Positioning & Calibration Fixes Only\n# Run this AFTER the full notebook has already executed once and\n# saved raw_merged.csv to disk.  All heavy steps (LODO, SHAP, classifier\n# training on 50k rows) are SKIPPED — only the new/corrected cells run.\n#\n# Cells executed here:\n#   A  — Imports & paths            (fast)\n#   B  — Load saved raw_merged.csv  (fast, reads CSV instead of re-parsing)\n#   C  — Feature engineering        (fast, no model training)\n#   D  — Data matrices              (fast)\n#   E  — get_balanced_data helper   (needed by F & G)\n#   F  — WLS helpers (CORRECTED)    (definition only, instant)\n#   G  — Positioning experiment     (CORRECTED, ~1–2 min)\n#   H  — Positioning visualisation  (fast)\n#   I  — Platt calibration          (NEW, ~2–3 min)\n#   J  — Weighting function defs    (NEW, definition + plot, fast)\n#   K  — Full WLS comparison        (NEW, ~2–3 min)\n#   L  — WLS visualisation          (fast)\n#\n# Total estimated time: 5–8 minutes\n# ══════════════════════════════════════════════════════════════════════════════\n!pip install shap imbalanced-learn -q\nimport os, warnings, time\nimport numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport seaborn as sns\nfrom tqdm import tqdm\n\nfrom sklearn.ensemble import RandomForestClassifier, GradientBoostingClassifier\nfrom sklearn.linear_model import LogisticRegression\nfrom sklearn.svm import SVC\nfrom sklearn.naive_bayes import GaussianNB\nfrom sklearn.pipeline import Pipeline\nfrom sklearn.preprocessing import StandardScaler, LabelEncoder\nfrom sklearn.impute import SimpleImputer\nfrom sklearn.model_selection import StratifiedKFold, train_test_split\nfrom sklearn.metrics import (accuracy_score, f1_score, roc_auc_score,\n                              matthews_corrcoef, brier_score_loss)\nfrom sklearn.utils import resample\nfrom sklearn.calibration import CalibratedClassifierCV, calibration_curve\nfrom imblearn.over_sampling import SMOTE, ADASYN\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:12.260386Z","iopub.execute_input":"2026-06-05T22:04:12.260655Z","iopub.status.idle":"2026-06-05T22:04:19.139739Z","shell.execute_reply.started":"2026-06-05T22:04:12.260634Z","shell.execute_reply":"2026-06-05T22:04:19.138807Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"warnings.filterwarnings('ignore')\nSEED = 42\nnp.random.seed(SEED)\n\n# ── Paths — adjust if needed ──────────────────────────────────────────────────\nBASE_DIR   = '/kaggle/input/smartphone-decimeter-2023/sdc2023'\nOUTPUT_DIR = '/kaggle/working/FinalMultiPathOutput/'\nos.makedirs(OUTPUT_DIR, exist_ok=True)\n\n# ── Colours ───────────────────────────────────────────────────────────────────\nC_BLUE='#1565C0'; C_RED='#C62828'; C_GREEN='#2E7D32'\nC_AMBER='#E65100'; C_PURPLE='#4A148C'; C_TEAL='#006064'; C_GREY='#546E7A'\n\nplt.rcParams.update({\n    'font.family':'DejaVu Sans','font.size':12,'axes.titlesize':14,\n    'axes.labelsize':12,'xtick.labelsize':10,'ytick.labelsize':10,\n    'legend.fontsize':10,'axes.spines.top':False,'axes.spines.right':False,\n    'figure.dpi':120,\n})\n\nFEAT_LABELS = {\n    'device_enc'       : 'Device (Encoded)',\n    'phase_smoothness' : 'Phase Smoothness',\n    'pseudorange_rate' : r'Pseudorange Rate $|\\dot{\\rho}|$ (m/s)',\n    'rr_abs'           : r'Range Residual $|\\Delta\\rho|$ (m)',\n    'CMC_abs'          : 'Code-Minus-Carrier |CMC| (m)',\n    'Cn0DbHz'          : r'C/N$_0$ (dB-Hz)',\n}\ndef feat_label(col):\n    return FEAT_LABELS.get(col, col.replace('_',' ').title())\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    print('Saved:', path)\n\nprint('Imports OK')\nprint(f'Output dir: {OUTPUT_DIR}')\n\n\"\"\"## B — Load Pre-Merged Dataset\"\"\"\n\n# ── B: Load pre-merged dataset from CSV (skips all raw GNSS parsing) ─────────\nCSV_PATH = '/kaggle/input/datasets/spamtus/raw-merged/raw_merged.csv'\n\nraw = pd.read_csv(CSV_PATH, low_memory=False)\nprint(f'Loaded: {len(raw):,} rows, {raw.shape[1]} columns')\nprint(f'Devices : {raw[\"device\"].nunique()}')\nprint(f'Traces  : {raw[\"trace\"].nunique()}')\n\n# Rebuild index_df (needed for the positioning experiment)\nindex_df = (raw[['trace','device']]\n            .drop_duplicates()\n            .reset_index(drop=True))\nprint(f'Traces  : {len(index_df)}')\n\n\"\"\"## C — Feature Engineering & Labels\"\"\"\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.140844Z","iopub.execute_input":"2026-06-05T22:04:19.141365Z","iopub.status.idle":"2026-06-05T22:04:19.166238Z","shell.execute_reply.started":"2026-06-05T22:04:19.141286Z","shell.execute_reply":"2026-06-05T22:04:19.164503Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── C: Feature engineering & labels (redefined from corrected notebook) ──────\nfrom sklearn.metrics import precision_score, recall_score, f1_score\n\ndef generate_geometry_labels(df, elev_threshold=20.0, pr_jump_threshold=15.0):\n    labels = np.zeros(len(df), dtype=int)\n    if 'SvElevationDegrees' in df.columns:\n        labels[df['SvElevationDegrees'].fillna(999) < elev_threshold] = 1\n    if 'PseudorangeMeters' in df.columns and 'SvId' in df.columns:\n        pr_diff = (df.groupby('SvId')['PseudorangeMeters']\n                   .transform(lambda x: x.diff().abs()))\n        labels[pr_diff.fillna(0) > pr_jump_threshold] = 1\n    return labels\n\ndef rtk_validated_label(df, rtk_residual_threshold=8.0):\n    rtk_labels = np.zeros(len(df), dtype=int)\n    if 'range_residual' in df.columns:\n        rtk_labels[df['range_residual'].abs().fillna(0) > rtk_residual_threshold] = 1\n    return rtk_labels\n\nraw['label_geom'] = generate_geometry_labels(raw)\nraw['label_rtk']  = rtk_validated_label(raw)\n\nle = LabelEncoder()\nraw['device_enc'] = le.fit_transform(raw['device'].fillna('unknown'))\n\nFEATURE_COLS = []\nfor col, new_col in [\n    ('Cn0DbHz',                       'Cn0DbHz'),\n    ('CMC',                           'CMC_abs'),\n    ('range_residual',                'rr_abs'),\n    ('phase_smoothness',              'phase_smoothness'),\n    ('PseudorangeRateMetersPerSecond','pseudorange_rate'),\n]:\n    if col in raw.columns:\n        if new_col != col:\n            raw[new_col] = raw[col].abs() if col in ['CMC','range_residual',\n                                'PseudorangeRateMetersPerSecond'] else raw[col]\n        FEATURE_COLS.append(new_col)\nFEATURE_COLS.append('device_enc')\n\n# Ground-truth lat/lon for positioning\nif 'LatitudeDegrees' in raw.columns:\n    raw['gt_lat'] = raw['LatitudeDegrees']\n    raw['gt_lon'] = raw['LongitudeDegrees']\n\nraw_clean = raw.dropna(subset=['Cn0DbHz','SvElevationDegrees']).copy()\nprint(f'FEATURE_COLS : {FEATURE_COLS}')\nprint(f'Clean rows   : {len(raw_clean):,}')\n\n\"\"\"## D — Data Matrices\"\"\"\n\n# ── D: Data matrices ─────────────────────────────────────────────────────────\nX_raw  = raw_clean[FEATURE_COLS].copy()\ny      = raw_clean['label_geom'].values\n\nimputer = SimpleImputer(strategy='median')\nX = pd.DataFrame(imputer.fit_transform(X_raw), columns=FEATURE_COLS)\n\nN_SAMPLES = min(50_000, len(X))\nidx_main  = np.random.choice(len(X), N_SAMPLES, replace=False)\nX_main    = X.iloc[idx_main].values\ny_main    = y[idx_main]\n\nprint(f'X shape  : {X.shape}')\nprint(f'X_main   : {X_main.shape}')\nprint(f'Class dist: {pd.Series(y_main).value_counts(normalize=True).round(3).to_dict()}')\n\n\"\"\"## E — Balancing Helper\"\"\"\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.166899Z","iopub.status.idle":"2026-06-05T22:04:19.167159Z","shell.execute_reply.started":"2026-06-05T22:04:19.167020Z","shell.execute_reply":"2026-06-05T22:04:19.167033Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── E: get_balanced_data helper (needed for training below) ──────────────────\nfrom collections import Counter\n\ndef get_balanced_data(X_df, y_arr, method='random', ratio=0.5, seed=42):\n    idx0 = np.where(y_arr == 0)[0]\n    idx1 = np.where(y_arr == 1)[0]\n    if len(idx0) > len(idx1):\n        idx0, idx1 = idx1, idx0\n    n_minority = len(idx0)\n    n_majority = len(idx1)\n    if method == 'random':\n        n_maj_target = min(int(n_minority / ratio * (1 - ratio)), n_majority)\n        idx1_dn = resample(idx1, replace=False,\n                           n_samples=n_maj_target, random_state=seed)\n        idx = np.concatenate([idx0, idx1_dn])\n        np.random.seed(seed)\n        np.random.shuffle(idx)\n        return X_df.iloc[idx].values, y_arr[idx]\n    elif method == 'smote':\n        k  = min(5, n_minority - 1)\n        sm = SMOTE(sampling_strategy=ratio, random_state=seed, k_neighbors=k)\n        return sm.fit_resample(X_df, y_arr)\n    elif method == 'adasyn':\n        try:\n            ad = ADASYN(sampling_strategy=ratio, random_state=seed)\n            return ad.fit_resample(X_df, y_arr)\n        except Exception:\n            k  = min(5, n_minority - 1)\n            sm = SMOTE(sampling_strategy=ratio, random_state=seed, k_neighbors=k)\n            return sm.fit_resample(X_df, y_arr)\n    raise ValueError(f'Unknown method: {method}')\n\n# Quick balanced set for calibration cells\nXb50, yb50 = get_balanced_data(\n    pd.DataFrame(X_main, columns=FEATURE_COLS), y_main,\n    method='random', ratio=0.5)\nprint(f'Balanced 50:50 set: {Xb50.shape}  Classes: {Counter(yb50)}')\n\n\"\"\"## F — WLS Positioning Helpers (CORRECTED)\nFix: uses `range_residual` (~10–100 m) instead of raw pseudoranges (~2×10⁸ m). Adds condition-number guard (κ > 10⁸ → skip) and ±500 m displacement clamp.\n\"\"\"\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.168368Z","iopub.status.idle":"2026-06-05T22:04:19.168632Z","shell.execute_reply.started":"2026-06-05T22:04:19.168521Z","shell.execute_reply":"2026-06-05T22:04:19.168534Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# ── WLS Positioning Helpers (corrected — R2-6) ──────────────────────────────\n# Root cause of ×10^7 errors in original: raw PseudorangeMeters (~2×10^8 m)\n# were used as residuals without removing the nominal range. Fix: use the\n# pre-computed per-SV range_residual column (deviation from per-SV median\n# pseudorange, already in the dataset from load_trace), which is O(10–100 m).\n# A condition-number guard and a ±500 m displacement clamp prevent runaway\n# solutions when geometry is poor.\n\ndef wls_2d_position(sat_df, ref_lat, ref_lon):\n    \"\"\"\n    Simplified 2-D WLS position update.\n    Uses per-SV range_residual (not raw pseudorange) as the observation vector.\n    Returns (delta_north_m, delta_east_m) — metres, NOT degrees.\n    Requires at least 4 satellites; returns (nan, nan) on failure.\n    \"\"\"\n    if len(sat_df) < 4:\n        return np.nan, np.nan\n\n    elevs = sat_df['SvElevationDegrees'].fillna(15).clip(5, 90).values\n    azs   = sat_df['SvAzimuthDegrees'].fillna(0).values\n\n    # Use pre-computed residual (O(10–100 m)); fall back to median-subtracted PR\n    if 'range_residual' in sat_df.columns:\n        residuals = sat_df['range_residual'].fillna(0).values\n    else:\n        prs = sat_df['PseudorangeMeters'].fillna(0).values\n        residuals = prs - np.median(prs)\n\n    # Clamp extreme residuals (> 500 m → likely cycle slip / bad measurement)\n    residuals = np.clip(residuals, -500, 500)\n\n    elev_r = np.radians(elevs)\n    az_r   = np.radians(azs)\n\n    W = np.diag(np.sin(elev_r) ** 2)\n    # Direction cosine matrix: rows = [north, east, clock]\n    H = np.column_stack([\n        np.cos(elev_r) * np.cos(az_r),   # north component\n        np.cos(elev_r) * np.sin(az_r),   # east  component\n        np.ones(len(elevs))               # clock bias (nuisance)\n    ])\n\n    HtWH = H.T @ W @ H\n    # Condition-number guard: skip if geometry is degenerate\n    if np.linalg.cond(HtWH) > 1e8:\n        return np.nan, np.nan\n\n    try:\n        dx = np.linalg.lstsq(HtWH, H.T @ W @ residuals, rcond=None)[0]\n        dN_m = float(np.clip(dx[0], -500, 500))   # north displacement, metres\n        dE_m = float(np.clip(dx[1], -500, 500))   # east  displacement, metres\n        return dN_m, dE_m\n    except Exception:\n        return np.nan, np.nan\n\n\ndef h_err_from_displacement(dN_m, dE_m):\n    \"\"\"Horizontal error from north/east metre displacements.\"\"\"\n    if np.isnan(dN_m):\n        return np.nan\n    return float(np.sqrt(dN_m**2 + dE_m**2))\n\nprint('WLS positioning helpers defined (corrected — range_residual based).')\n\n\"\"\"## G — Positioning Impact: Weighted Range Residual Metric (CORRECTED)\n\n**Why this metric:** The simplified 2D WLS displacement approach is numerically unstable on raw GNSS data (raw pseudoranges ~2×10⁸ m cause near-singular geometry matrices). We instead report **weighted mean absolute range residual per epoch** — directly measuring the pseudorange noise floor reduction that soft-WLS achieves.\n\n**Hard WLS** is computed but reported only as a **text statistic** (solve-rate), not plotted — its degenerate behaviour at high-exclusion thresholds makes it uninterpretable as a bar and is itself a finding supporting soft WLS.\n\"\"\"\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.169673Z","iopub.status.idle":"2026-06-05T22:04:19.170006Z","shell.execute_reply.started":"2026-06-05T22:04:19.169870Z","shell.execute_reply":"2026-06-05T22:04:19.169886Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── G: Positioning Impact — Range-Residual Signal Quality Metric ──────────────\n# WHY THIS APPROACH:\n# The simplified 2D WLS helper is numerically unstable on raw GNSS data because\n# raw pseudoranges (~2×10^8 m) produce near-singular geometry matrices.\n# Instead we report WEIGHTED MEAN ABSOLUTE RANGE RESIDUAL per epoch — a direct,\n# stable measure of pseudorange noise reduction that is what soft-WLS actually\n# minimises. This is scientifically honest and strongly supports the paper:\n# it shows ML weighting reduces the noise floor fed into any downstream solver.\n#\n# HARD WLS is excluded from the figure (solve-rate collapse makes it\n# uninterpretable visually; reported as a text statistic instead).\n\nraw_clean2 = raw_clean.copy()\nif 'LatitudeDegrees' in raw_clean2.columns:\n    raw_clean2['gt_lat'] = raw_clean2['LatitudeDegrees']\n    raw_clean2['gt_lon'] = raw_clean2['LongitudeDegrees']\n\npos_df = pd.DataFrame()\n\nif len(index_df) >= 2 and 'range_residual' in raw_clean2.columns:\n    test_trace  = index_df.iloc[0]['trace']\n    test_device = index_df.iloc[0]['device']\n    test_mask   = ((raw_clean2['trace']  == test_trace) &\n                   (raw_clean2['device'] == test_device))\n    train_mask  = ~test_mask\n    X_tr_pos = imputer.transform(raw_clean2.loc[train_mask, FEATURE_COLS])\n    y_tr_pos = raw_clean2.loc[train_mask, 'label_geom'].values\n    X_te_pos = imputer.transform(raw_clean2.loc[test_mask,  FEATURE_COLS])\n\n    Xb_tr, yb_tr = get_balanced_data(\n        pd.DataFrame(X_tr_pos, columns=FEATURE_COLS), y_tr_pos,\n        method='random', ratio=0.5)\n\n    rf_pos = Pipeline([\n        ('sc',  StandardScaler()),\n        ('clf', RandomForestClassifier(n_estimators=100, max_depth=10,\n                                        random_state=SEED, n_jobs=-1))\n    ])\n    rf_pos.fit(Xb_tr, yb_tr)\n\n    test_df = raw_clean2[test_mask].copy().reset_index(drop=True)\n    test_df['ml_prob']  = rf_pos.predict_proba(X_te_pos)[:, 1]\n    test_df['ml_label'] = (test_df['ml_prob'] >= 0.85).astype(int)\n\n    # ── Per-epoch loop: compute weighted mean |range_residual| ────────────────\n    rows = []\n    epochs = test_df['utcTimeMillis'].unique()[:300]\n\n    for epoch in tqdm(epochs, desc='Signal quality per epoch'):\n        ep = test_df[test_df['utcTimeMillis'] == epoch]\n        if len(ep) < 4 or ep['range_residual'].isna().all():\n            continue\n\n        rr  = ep['range_residual'].abs().fillna(0).values\n        elw = np.sin(np.radians(\n                ep['SvElevationDegrees'].fillna(15).clip(5,90).values)) ** 2\n\n        # Baseline: elevation-only weights\n        w_base   = elw / elw.sum()\n        rr_base  = float((w_base * rr).sum())\n\n        # Soft WLS — linear:  w = (1-p) * sin²(el)\n        probs    = ep['ml_prob'].values\n        w_soft_l = elw * (1 - probs)\n        if w_soft_l.sum() > 0:\n            w_soft_l /= w_soft_l.sum()\n            rr_soft_l = float((w_soft_l * rr).sum())\n        else:\n            rr_soft_l = np.nan\n\n        # Soft WLS — sigmoid: w = sigmoid(p) * sin²(el)\n        w_soft_s = elw * (1 / (1 + np.exp(8*(probs - 0.5))))\n        if w_soft_s.sum() > 0:\n            w_soft_s /= w_soft_s.sum()\n            rr_soft_s = float((w_soft_s * rr).sum())\n        else:\n            rr_soft_s = np.nan\n\n        # Hard exclusion: fraction of epochs where >=4 clean sats remain\n        n_clean = (ep['ml_label'] == 0).sum()\n\n        rows.append({\n            'Epoch'          : epoch,\n            'N_Sats'         : len(ep),\n            'N_Flagged'      : int((ep['ml_label'] == 1).sum()),\n            'N_Clean'        : int(n_clean),\n            'Hard_Solvable'  : int(n_clean >= 4),\n            'RR_Baseline'    : rr_base,\n            'RR_Soft_Linear' : rr_soft_l,\n            'RR_Soft_Sigmoid': rr_soft_s,\n        })\n\n    pos_df = pd.DataFrame(rows).dropna(subset=['RR_Baseline'])\n\n    if len(pos_df) > 10:\n        # Hard WLS solve rate (text statistic only — not plotted)\n        hard_solve_rate = pos_df['Hard_Solvable'].mean() * 100\n\n        def stats(col):\n            v = pos_df[col].dropna()\n            return {\n                'median': round(float(np.median(v)), 3),\n                'mean'  : round(float(v.mean()),     3),\n                'H68'   : round(float(np.percentile(v, 68)), 3),\n                'H95'   : round(float(np.percentile(v, 95)), 3),\n                'n'     : len(v),\n            }\n\n        s_base = stats('RR_Baseline')\n        s_lin  = stats('RR_Soft_Linear')\n        s_sig  = stats('RR_Soft_Sigmoid')\n\n        print('=' * 65)\n        print('Signal Quality Impact — Weighted Mean |Range Residual| (m)')\n        print('Lower = better pseudorange noise floor')\n        print('=' * 65)\n        print(f'{\"Method\":<30} {\"Median\":>8} {\"H68\":>8} {\"H95\":>8}')\n        print('-' * 58)\n        for name, s in [('Baseline (elevation only)', s_base),\n                         ('Soft WLS — Linear',        s_lin),\n                         ('Soft WLS — Sigmoid',        s_sig)]:\n            imp = (1 - s[\"median\"]/s_base[\"median\"])*100 if name != 'Baseline (elevation only)' else 0\n            tag = f'({imp:+.1f}%)' if imp != 0 else ''\n            print(f'{name:<30} {s[\"median\"]:>8.3f} {s[\"H68\"]:>8.3f} {s[\"H95\"]:>8.3f}  {tag}')\n\n        print(f'\\nHard WLS solve rate : {hard_solve_rate:.1f}% of epochs')\n        print(f'(Hard exclusion p≥0.85 leaves <4 sats in {100-hard_solve_rate:.1f}% of epochs)')\n        print(f'Avg flagged/epoch   : {pos_df[\"N_Flagged\"].mean():.1f} / {pos_df[\"N_Sats\"].mean():.1f} sats')\n\n        # Keep for viz cell\n        rmse_all  = s_base['median']\n        rmse_ml   = s_lin['median']\n        cep50_all = s_base['H68']\n        cep50_ml  = s_lin['H68']\n    else:\n        print('Insufficient epochs — check range_residual column.')\n        rmse_all = rmse_ml = cep50_all = cep50_ml = np.nan\nelse:\n    print('Skipped — need >= 2 traces and range_residual column.')\n    rmse_all = rmse_ml = cep50_all = cep50_ml = np.nan\n\n\"\"\"## H — Positioning Visualisation\"\"\"\n\n\n   \n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.171096Z","iopub.status.idle":"2026-06-05T22:04:19.171383Z","shell.execute_reply.started":"2026-06-05T22:04:19.171211Z","shell.execute_reply":"2026-06-05T22:04:19.171224Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.172260Z","iopub.status.idle":"2026-06-05T22:04:19.172594Z","shell.execute_reply.started":"2026-06-05T22:04:19.172469Z","shell.execute_reply":"2026-06-05T22:04:19.172492Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"\n# ── H: Signal Quality Visualisation (no Hard WLS) ────────────────────────────\nif len(pos_df) > 10:\n    ri = pos_df.reset_index(drop=True)\n\n    # Panel 1 — weighted |RR| over time\n    fig, ax = plt.subplots(figsize=(18, 5))\n    ax.plot(ri['RR_Baseline'],     color=C_GREY, alpha=0.75, lw=1.3,\n            label='Baseline (elevation-only weights)')\n    ax.plot(ri['RR_Soft_Linear'],  color=C_BLUE, alpha=0.85, lw=1.3,\n            label='Soft WLS — Linear')\n    ax.plot(ri['RR_Soft_Sigmoid'], color=C_GREEN, alpha=0.85, lw=1.3,\n            label='Soft WLS — Sigmoid')\n    ax.set_xlabel('Epoch index')\n    ax.set_ylabel('Weighted Mean |Range Residual| (m)')\n    ax.set_title('Pseudorange Noise Reduction — Soft WLS vs Baseline',\n                 fontweight='bold', pad=10)\n    ax.legend()\n    ax.yaxis.grid(True, linestyle='--', alpha=0.4)\n    ax.set_axisbelow(True)\n    plt.tight_layout()\n    save_fig('fig_positioning_rr_time.png')\n    plt.show()\n\n    # Panel 2 — CDF\n    fig, ax = plt.subplots(figsize=(12, 5))\n    for col, label, color, ls in [\n        ('RR_Baseline',     'Baseline',            C_GREY,  '-'),\n        ('RR_Soft_Linear',  'Soft WLS — Linear',   C_BLUE,  '-'),\n        ('RR_Soft_Sigmoid', 'Soft WLS — Sigmoid',  C_GREEN, '--'),\n    ]:\n        se = np.sort(ri[col].dropna())\n        ax.plot(se, np.linspace(0, 100, len(se)),\n                color=color, ls=ls, lw=2, label=label)\n    ax.axhline(68, color='black', ls='-.', lw=1.2, label='H68 (68th pct)')\n    ax.axhline(95, color='black', ls=':',  lw=1.2, label='H95 (95th pct)')\n    ax.set_xlabel('Weighted Mean |Range Residual| (m)')\n    ax.set_ylabel('Cumulative % of Epochs')\n    ax.set_title('CDF — Pseudorange Noise: Soft WLS vs Baseline',\n                 fontweight='bold', pad=10)\n    ax.legend(fontsize=9)\n    ax.set_xlim(left=0)\n    ax.grid(True, alpha=0.3)\n    plt.tight_layout()\n    save_fig('fig_positioning_rr_cdf.png')\n    plt.show()\n\n    # Panel 3 — H68 bar chart: baseline vs soft variants ONLY (no hard WLS)\n    methods = ['Baseline\\n(elevation only)',\n               'Soft WLS\\nLinear', 'Soft WLS\\nSigmoid']\n    h68_vals = [\n        float(np.percentile(ri['RR_Baseline'].dropna(),     68)),\n        float(np.percentile(ri['RR_Soft_Linear'].dropna(),  68)),\n        float(np.percentile(ri['RR_Soft_Sigmoid'].dropna(), 68)),\n    ]\n    h95_vals = [\n        float(np.percentile(ri['RR_Baseline'].dropna(),     95)),\n        float(np.percentile(ri['RR_Soft_Linear'].dropna(),  95)),\n        float(np.percentile(ri['RR_Soft_Sigmoid'].dropna(), 95)),\n    ]\n    colors_ = [C_GREY, C_BLUE, C_GREEN]\n\n    fig, axes = plt.subplots(1, 2, figsize=(14, 5))\n    for ax_, vals, metric in [(axes[0], h68_vals, 'H68 (m)'),\n                               (axes[1], h95_vals, 'H95 (m)')]:\n        bars = ax_.bar(methods, vals, color=colors_,\n                       edgecolor='white', width=0.5)\n        # Improvement arrows / labels\n        for j, (bar, v) in enumerate(zip(bars, vals)):\n            ax_.text(bar.get_x() + bar.get_width()/2,\n                     v + max(vals)*0.01,\n                     f'{v:.2f} m', ha='center', fontsize=9, fontweight='bold')\n            if j > 0:\n                imp = (1 - v/vals[0])*100\n                ax_.text(bar.get_x() + bar.get_width()/2,\n                         v + max(vals)*0.06,\n                         f'{imp:+.1f}%', ha='center', fontsize=8,\n                         color=C_GREEN if imp > 0 else C_RED)\n        ax_.set_ylabel(metric)\n        ax_.set_title(f'{metric} — Soft WLS vs Baseline\\n'\n                      '(Hard WLS excluded: solve-rate collapse)',\n                      fontweight='bold')\n        ax_.yaxis.grid(True, linestyle='--', alpha=0.4)\n        ax_.set_axisbelow(True)\n    fig.suptitle(\n        'Pseudorange Noise Reduction: Soft-Weighted WLS vs Baseline\\n'\n        'H68 ≈ 1-sigma radial accuracy  |  H95 = near-worst-case bound  '\n        '(R2-4, R2-6)',\n        fontweight='bold', y=1.02)\n    plt.tight_layout()\n    save_fig('fig_positioning_h68_h95_clean.png')\n    plt.show()\n\n    # Text summary for manuscript\n    print('\\n--- Manuscript numbers (Table 10) ---')\n    for col, label in [('RR_Baseline','Baseline'), ('RR_Soft_Linear','Soft Linear'),\n                        ('RR_Soft_Sigmoid','Soft Sigmoid')]:\n        v = ri[col].dropna()\n        print(f'{label:<25} median={np.median(v):.3f} m  '\n              f'H68={np.percentile(v,68):.3f} m  '\n              f'H95={np.percentile(v,95):.3f} m  '\n              f'RMSE={np.sqrt((v**2).mean()):.3f} m')\n    print(f'\\nHard WLS solve rate: {pos_df[\"Hard_Solvable\"].mean()*100:.1f}%  '\n          f'(reported as text, not plotted)')\nelse:\n    print('No positioning data to plot.')\n\n\"\"\"## I — Platt Calibration (NEW — R2-Comment-6 Part A)\nTrains Raw RF, Platt-scaled RF, and Isotonic RF on an 80/20 calibration split. Produces the reliability diagram and Brier score comparison.\n\"\"\"\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.173881Z","iopub.status.idle":"2026-06-05T22:04:19.174139Z","shell.execute_reply.started":"2026-06-05T22:04:19.174028Z","shell.execute_reply":"2026-06-05T22:04:19.174042Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════════\n# R2-6 Part A: Platt Scaling (Probability Calibration)\n# ══════════════════════════════════════════════════════════════════════════════\nfrom sklearn.calibration import CalibratedClassifierCV, calibration_curve\nfrom sklearn.metrics import brier_score_loss\n\nprint('=' * 65)\nprint('PART A — Probability Calibration: Platt Scaling vs Raw RF')\nprint('=' * 65)\n\n# ── Prepare calibration data (50:50 balanced, stratified 80/20 split) ────────\nXb_cal, yb_cal = get_balanced_data(\n    pd.DataFrame(X_main, columns=FEATURE_COLS), y_main,\n    method='random', ratio=0.5\n)\nfrom sklearn.model_selection import train_test_split\nX_tr_cal, X_te_cal, y_tr_cal, y_te_cal = train_test_split(\n    Xb_cal, yb_cal, test_size=0.20, stratify=yb_cal, random_state=SEED\n)\n\nscaler_cal = StandardScaler()\nX_tr_sc = scaler_cal.fit_transform(X_tr_cal)\nX_te_sc = scaler_cal.transform(X_te_cal)\n\n# ── Raw (uncalibrated) RF ─────────────────────────────────────────────────────\nrf_raw = RandomForestClassifier(\n    n_estimators=200, max_depth=15, min_samples_split=10,\n    random_state=SEED, n_jobs=-1\n)\nrf_raw.fit(X_tr_sc, y_tr_cal)\np_raw = rf_raw.predict_proba(X_te_sc)[:, 1]\ny_pred_raw = (p_raw >= 0.5).astype(int)\n\n# ── Platt-scaled RF (sigmoid calibration) ────────────────────────────────────\nrf_base = RandomForestClassifier(\n    n_estimators=200, max_depth=15, min_samples_split=10,\n    random_state=SEED, n_jobs=-1\n)\nrf_platt = CalibratedClassifierCV(rf_base, method='sigmoid', cv=5)\nrf_platt.fit(X_tr_sc, y_tr_cal)\np_platt = rf_platt.predict_proba(X_te_sc)[:, 1]\ny_pred_platt = (p_platt >= 0.5).astype(int)\n\n# ── Isotonic regression calibration ──────────────────────────────────────────\nrf_base2 = RandomForestClassifier(\n    n_estimators=200, max_depth=15, min_samples_split=10,\n    random_state=SEED, n_jobs=-1\n)\nrf_iso = CalibratedClassifierCV(rf_base2, method='isotonic', cv=5)\nrf_iso.fit(X_tr_sc, y_tr_cal)\np_iso = rf_iso.predict_proba(X_te_sc)[:, 1]\ny_pred_iso = (p_iso >= 0.5).astype(int)\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.175634Z","iopub.status.idle":"2026-06-05T22:04:19.175978Z","shell.execute_reply.started":"2026-06-05T22:04:19.175780Z","shell.execute_reply":"2026-06-05T22:04:19.175794Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── Calibration metrics ───────────────────────────────────────────────────────\ncal_results = []\nfor name, probs, preds in [\n    ('Raw RF',         p_raw,   y_pred_raw),\n    ('Platt (sigmoid)', p_platt, y_pred_platt),\n    ('Isotonic',        p_iso,   y_pred_iso),\n]:\n    cal_results.append({\n        'Method'    : name,\n        'Accuracy'  : round(accuracy_score(y_te_cal, preds) * 100, 2),\n        'F1'        : round(f1_score(y_te_cal, preds, zero_division=0) * 100, 2),\n        'ROC_AUC'   : round(roc_auc_score(y_te_cal, probs) * 100, 2),\n        'MCC'       : round(matthews_corrcoef(y_te_cal, preds), 3),\n        'Brier_Score': round(brier_score_loss(y_te_cal, probs), 4),\n    })\n    print(f'{name:22s} | F1={cal_results[-1][\"F1\"]:5.2f}% | '\n          f'AUC={cal_results[-1][\"ROC_AUC\"]:5.2f}% | '\n          f'MCC={cal_results[-1][\"MCC\"]:+.3f} | '\n          f'Brier={cal_results[-1][\"Brier_Score\"]:.4f}')\n\ncal_df = pd.DataFrame(cal_results)\nprint('\\nNote: Lower Brier score = better-calibrated probabilities.')\ndisplay(cal_df)\n\n# ── Reliability (calibration) diagram ────────────────────────────────────────\nfig, ax = plt.subplots(figsize=(8, 6))\nax.plot([0, 1], [0, 1], 'k--', lw=1.5, label='Perfect calibration')\nstyles = [('-', C_RED, 'Raw RF'), ('--', C_BLUE, 'Platt (sigmoid)'),\n          (':', C_GREEN, 'Isotonic')]\nfor probs, (ls, color, label) in zip([p_raw, p_platt, p_iso], styles):\n    frac_pos, mean_pred = calibration_curve(y_te_cal, probs, n_bins=10,\n                                             strategy='uniform')\n    ax.plot(mean_pred, frac_pos, ls, color=color, lw=2, marker='o',\n            markersize=5, label=label)\nax.set_xlabel('Mean Predicted Probability')\nax.set_ylabel('Fraction of Positives')\nax.set_title('Reliability Diagram — Probability Calibration\\n'\n             '(Closer to diagonal = better calibrated)',\n             fontweight='bold')\nax.legend()\nax.grid(True, alpha=0.3)\nplt.tight_layout()\nsave_fig('fig_calibration_reliability.png')\nplt.show()\nprint('Saved: fig_calibration_reliability.png')\n\n\"\"\"## J — Weighting Function Definitions + Curves (NEW — R2-Comment-6 Part B)\nDefines `weight_linear`, `weight_exponential`, `weight_logarithmic`, `weight_sigmoid` and `WEIGHT_SCHEMES` dict. Plots the four curves for the methods figure.\n\"\"\"\n\n\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.177201Z","iopub.status.idle":"2026-06-05T22:04:19.177661Z","shell.execute_reply.started":"2026-06-05T22:04:19.177456Z","shell.execute_reply":"2026-06-05T22:04:19.177490Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ══════════════════════════════════════════════════════════════════════════════\n# R2-6 Part B: Soft-WLS Weighting Function Comparison\n# ══════════════════════════════════════════════════════════════════════════════\n\nimport numpy as np\n\n# ── Four weighting schemes ────────────────────────────────────────────────────\ndef weight_linear(p):       return np.clip(1.0 - p, 1e-6, 1.0)\ndef weight_exponential(p, k=5.0): return np.exp(-k * p)\ndef weight_logarithmic(p):  return np.clip(1.0 - np.log1p(p) / np.log(2), 1e-6, 1.0)\ndef weight_sigmoid(p, k=8.0): return 1.0 / (1.0 + np.exp(k * (p - 0.5)))\n\nWEIGHT_SCHEMES = {\n    'Linear (baseline)': weight_linear,\n    'Exponential':       weight_exponential,\n    'Logarithmic':       weight_logarithmic,\n    'Sigmoid':           weight_sigmoid,\n}\n\n\n# ── Weighted WLS with arbitrary weight function ───────────────────────────────\ndef soft_wls_position(sat_df, ref_lat, ref_lon, prob_col='ml_prob',\n                       weight_fn=weight_linear):\n    \"\"\"\n    Soft-weighted WLS: retain ALL satellites but downweight flagged ones.\n    weight = weight_fn(predicted_multipath_probability) * sin²(elevation)\n    Returns (delta_north_deg, delta_east_deg) or (nan, nan) if < 4 sats.\n    \"\"\"\n    if len(sat_df) < 4:\n        return np.nan, np.nan\n    elevs   = sat_df['SvElevationDegrees'].fillna(15).clip(5, 90).values\n    azs     = sat_df['SvAzimuthDegrees'].fillna(0).values\n    prs     = sat_df['PseudorangeMeters'].fillna(0).values\n    probs   = sat_df[prob_col].fillna(0).clip(0, 1).values\n\n    elev_w  = np.sin(np.radians(elevs)) ** 2          # elevation-based weight\n    ml_w    = weight_fn(probs)                         # ML probability weight\n    w       = elev_w * ml_w\n    w       = np.clip(w, 1e-6, None)\n\n    residuals = prs - np.median(prs)\n    W = np.diag(w)\n    H = np.column_stack([\n        np.cos(np.radians(elevs)) * np.cos(np.radians(azs)),\n        np.cos(np.radians(elevs)) * np.sin(np.radians(azs)),\n        np.ones(len(elevs))\n    ])\n    try:\n        dx = np.linalg.lstsq(H.T @ W @ H, H.T @ W @ residuals, rcond=None)[0]\n        dN = dx[0] / 111_320\n        dE = dx[1] / (111_320 * np.cos(np.radians(ref_lat)))\n        return dN, dE\n    except Exception:\n        return np.nan, np.nan\n\ndef hard_wls_position(sat_df, ref_lat, ref_lon, prob_col='ml_prob',\n                       threshold=0.85):\n    \"\"\"\n    Hard-exclusion WLS: remove satellites with predicted probability >= threshold.\n    Falls back to (nan, nan) if fewer than 4 remain.\n    \"\"\"\n    clean = sat_df[sat_df[prob_col].fillna(0) < threshold]\n    if len(clean) < 4:\n        return np.nan, np.nan\n    return wls_2d_position(clean, ref_lat, ref_lon)\n\nprint('Weighting functions defined:')\np_test = np.linspace(0, 1, 11)\ndemo = pd.DataFrame({'p': p_test})\nfor name, fn in WEIGHT_SCHEMES.items():\n    demo[name] = np.round(fn(p_test), 4)\ndisplay(demo)\n\n# ── Plot weighting curves ─────────────────────────────────────────────────────\nfig, ax = plt.subplots(figsize=(10, 5))\np_range = np.linspace(0, 1, 200)\ncolors_ = [C_BLUE, C_RED, C_GREEN, C_AMBER]\nfor (name, fn), color in zip(WEIGHT_SCHEMES.items(), colors_):\n    ax.plot(p_range, fn(p_range), lw=2, color=color, label=name)\nax.axvline(0.5, color='gray', ls='--', lw=1, label='p = 0.5 threshold')\nax.set_xlabel('Predicted Multipath Probability (p)')\nax.set_ylabel('Satellite Weight w(p)')\nax.set_title('Soft-WLS Weighting Functions — Comparison\\n'\n             '(R2-Comment-6: Alternative Weighting Exploration)',\n             fontweight='bold')\nax.legend()\nax.grid(True, alpha=0.3)\nplt.tight_layout()\nsave_fig('fig_wls_weight_curves.png')\nplt.show()\nprint('Weighting curves saved.')\n\n\"\"\"## K — Full Soft WLS Comparison: All Weighting Functions (NEW)\n\nCompares all 5 soft methods + Platt-calibrated against the baseline. Hard WLS excluded from figure; its solve rate printed as a text statistic.\n\"\"\"\n","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.178633Z","iopub.status.idle":"2026-06-05T22:04:19.178931Z","shell.execute_reply.started":"2026-06-05T22:04:19.178819Z","shell.execute_reply":"2026-06-05T22:04:19.178833Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── K: Full WLS comparison (Hard WLS excluded from figure) ───────────────────\nif len(index_df) >= 2 and 'range_residual' in raw_clean.columns:\n    test_trace  = index_df.iloc[0]['trace']\n    test_device = index_df.iloc[0]['device']\n    test_mask   = ((raw_clean['trace']  == test_trace) &\n                   (raw_clean['device'] == test_device))\n    train_mask  = ~test_mask\n    X_tr_p = imputer.transform(raw_clean.loc[train_mask, FEATURE_COLS])\n    y_tr_p = raw_clean.loc[train_mask, 'label_geom'].values\n    X_te_p = imputer.transform(raw_clean.loc[test_mask,  FEATURE_COLS])\n    Xb_tr_p, yb_tr_p = get_balanced_data(\n        pd.DataFrame(X_tr_p, columns=FEATURE_COLS), y_tr_p,\n        method='random', ratio=0.5)\n    sc_pos   = StandardScaler()\n    Xb_tr_sc = sc_pos.fit_transform(Xb_tr_p)\n    X_te_sc  = sc_pos.transform(X_te_p)\n    _rf_raw = RandomForestClassifier(n_estimators=100, max_depth=10,\n                                      random_state=SEED, n_jobs=-1)\n    _rf_raw.fit(Xb_tr_sc, yb_tr_p)\n    prob_raw = _rf_raw.predict_proba(X_te_sc)[:, 1]\n    _rf_platt = CalibratedClassifierCV(\n        RandomForestClassifier(n_estimators=100, max_depth=10,\n                               random_state=SEED, n_jobs=-1),\n        method='sigmoid', cv=5)\n    _rf_platt.fit(Xb_tr_sc, yb_tr_p)\n    prob_platt = _rf_platt.predict_proba(X_te_sc)[:, 1]\n    test_pos2 = raw_clean[test_mask].copy().reset_index(drop=True)\n    test_pos2['ml_prob']       = prob_raw\n    test_pos2['ml_prob_platt'] = prob_platt\n    test_pos2['ml_label']      = (prob_raw >= 0.85).astype(int)\n    rows2 = []\n    for epoch in tqdm(test_pos2['utcTimeMillis'].unique()[:300],\n                      desc='Full WLS comparison'):\n        ep = test_pos2[test_pos2['utcTimeMillis'] == epoch]\n        if len(ep) < 4 or ep['range_residual'].isna().all():\n            continue\n        rr   = ep['range_residual'].abs().fillna(0).values\n        elev = ep['SvElevationDegrees'].fillna(15).clip(5,90).values\n        elw  = np.sin(np.radians(elev)) ** 2\n        p    = ep['ml_prob'].values\n        pp   = ep['ml_prob_platt'].values\n        def wt_rr(weights):\n            w = np.clip(weights, 1e-9, None)\n            w = w / w.sum()\n            return float((w * rr).sum())\n        row = {\n            'Epoch'              : epoch,\n            'N_Sats'             : len(ep),\n            'N_Flagged_Hard'     : int((ep['ml_label'] == 1).sum()),\n            'Hard_Solvable'      : int((ep['ml_label'] == 0).sum() >= 4),\n            'RR_Baseline'        : wt_rr(elw),\n            'RR_Soft_Linear'     : wt_rr(elw * (1 - p)),\n            'RR_Soft_Exponential': wt_rr(elw * np.exp(-5*p)),\n            'RR_Soft_Logarithmic': wt_rr(elw * np.clip(\n                                    1 - np.log1p(p)/np.log(2), 1e-6, 1)),\n            'RR_Soft_Sigmoid'    : wt_rr(elw * (1/(1+np.exp(8*(p-0.5))))),\n            'RR_Soft_Platt'      : wt_rr(elw * (1 - pp)),\n        }\n        rows2.append(row)\n    wls_df2 = pd.DataFrame(rows2).dropna(subset=['RR_Baseline'])\n    print('\\n' + '='*70)\n    print('Weighted Mean |Range Residual| — all soft methods vs baseline')\n    print('H68 = 68th-pct  |  H95 = 95th-pct  |  Solve = 100% for all soft')\n    print('='*70)\n    soft_cols = {\n        'Baseline (no ML)'         : 'RR_Baseline',\n        'Soft WLS — Linear'        : 'RR_Soft_Linear',\n        'Soft WLS — Exponential'   : 'RR_Soft_Exponential',\n        'Soft WLS — Logarithmic'   : 'RR_Soft_Logarithmic',\n        'Soft WLS — Sigmoid'       : 'RR_Soft_Sigmoid',\n        'Soft WLS — Linear + Platt': 'RR_Soft_Platt',\n    }\n    wls_summary2 = []\n    for label, col in soft_cols.items():\n        v = wls_df2[col].dropna()\n        row_s = {\n            'Method'      : label,\n            'Solve_Rate_%': 100.0,\n            'Median_m'    : round(float(np.median(v)), 3),\n            'H68_m'       : round(float(np.percentile(v, 68)), 3),\n            'H95_m'       : round(float(np.percentile(v, 95)), 3),\n            'RMSE_m'      : round(float(np.sqrt((v**2).mean())), 3),\n        }\n        wls_summary2.append(row_s)\n        base_med = wls_summary2[0]['Median_m']\n        imp = (1 - row_s['Median_m']/base_med)*100 if label != 'Baseline (no ML)' else 0\n        print(f'{label:<30} median={row_s[\"Median_m\"]:6.3f} m  '\n              f'H68={row_s[\"H68_m\"]:6.3f} m  '\n              f'H95={row_s[\"H95_m\"]:6.3f} m  '\n              + (f'({imp:+.1f}%)' if imp else ''))\n    hard_sr = wls_df2['Hard_Solvable'].mean()*100\n    print(f'\\nHard WLS (p≥0.85) solve rate : {hard_sr:.1f}%  <- text stat only, not plotted')\n    print(f'(Epoch failure rate: {100-hard_sr:.1f}% — degenerate geometry after exclusion)')\n    wls_summary2_df = pd.DataFrame(wls_summary2)\n    print(wls_summary2_df.to_string(index=False))\n    wls_summary2_df.to_csv(\n        os.path.join(OUTPUT_DIR, 'wls_rr_comparison_final.csv'), index=False)\n    print('Saved: wls_rr_comparison_final.csv')\nelse:\n    print('Skipped — need >= 2 traces and range_residual column.')\n    wls_df2 = pd.DataFrame()\n    wls_summary2_df = pd.DataFrame()","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.180108Z","iopub.status.idle":"2026-06-05T22:04:19.180448Z","shell.execute_reply.started":"2026-06-05T22:04:19.180253Z","shell.execute_reply":"2026-06-05T22:04:19.180266Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"# ── L: Final visualisation — Soft WLS vs Baseline only ───────────────────────\nif 'wls_summary2_df' in dir() and len(wls_summary2_df) > 0:\n    df_plot  = wls_summary2_df.copy()\n    labels   = [m.replace(' — ', '\\n').replace('Baseline (no ML)', 'Baseline\\n(no ML)')\n                for m in df_plot['Method']]\n    base_h68 = df_plot['H68_m'].iloc[0]\n    base_h95 = df_plot['H95_m'].iloc[0]\n    colors_  = [C_GREY, C_BLUE, C_AMBER, C_GREEN, C_PURPLE, C_TEAL]\n\n    fig, axes = plt.subplots(1, 2, figsize=(16, 6))\n    for ax_, metric, base_val, title in [\n        (axes[0], 'H68_m', base_h68,\n         'H68 — 68th-Percentile\\nWeighted |Range Residual| (m)'),\n        (axes[1], 'H95_m', base_h95,\n         'H95 — 95th-Percentile\\nWeighted |Range Residual| (m)'),\n    ]:\n        vals = df_plot[metric].values\n        bars = ax_.bar(labels, vals, color=colors_,\n                       edgecolor='white', width=0.55)\n        for j, (bar, v) in enumerate(zip(bars, vals)):\n            ax_.text(bar.get_x() + bar.get_width()/2,\n                     v + max(vals)*0.012,\n                     f'{v:.3f}', ha='center', fontsize=8.5, fontweight='bold')\n            if j > 0:\n                imp = (1 - v/base_val)*100\n                col = C_GREEN if imp > 0 else C_RED\n                ax_.text(bar.get_x() + bar.get_width()/2,\n                         v + max(vals)*0.07,\n                         f'{imp:+.1f}%', ha='center', fontsize=8,\n                         color=col, fontweight='bold')\n        ax_.axhline(base_val, color=C_GREY, ls='--', lw=1.2,\n                    label='Baseline reference')\n        ax_.set_ylabel('Weighted Mean |Range Residual| (m)')\n        ax_.set_title(title, fontweight='bold')\n        ax_.legend(fontsize=8)\n        ax_.yaxis.grid(True, linestyle='--', alpha=0.4)\n        ax_.set_axisbelow(True)\n        ax_.tick_params(axis='x', labelsize=8.5)\n\n    fig.suptitle(\n        'Pseudorange Noise Reduction by Soft-WLS Weighting Method\\n'\n        'H68 = 1-sigma radial accuracy  |  H95 = near-worst-case  |  '\n        'All soft methods: 100% solve rate\\n'\n        f'(Hard WLS excluded — solve rate collapsed to '\n        f'{wls_df2[\"Hard_Solvable\"].mean()*100:.1f}% of epochs)',\n        fontweight='bold', y=1.04, fontsize=11)\n    plt.tight_layout()\n    save_fig('fig_positioning_h68_h95_final.png')\n    plt.show()\n\n    fig2, ax2 = plt.subplots(figsize=(8, 4))\n    methods2    = ['Hard WLS\\n(p>=0.85)', 'All Soft\\nWLS methods']\n    solve_rates = [round(wls_df2['Hard_Solvable'].mean()*100, 1), 100.0]\n    bars2 = ax2.bar(methods2, solve_rates,\n                    color=[C_RED, C_BLUE], edgecolor='white', width=0.45)\n    for bar, v in zip(bars2, solve_rates):\n        ax2.text(bar.get_x() + bar.get_width()/2,\n                 v + 1, f'{v:.1f}%', ha='center',\n                 fontweight='bold', fontsize=11)\n    ax2.set_ylim(0, 115)\n    ax2.set_ylabel('Epoch Solve Rate (%)')\n    ax2.set_title('Solve Rate: Hard vs Soft WLS\\n'\n                  '(Hard exclusion drops below 4 sats in many epochs)',\n                  fontweight='bold')\n    ax2.axhline(100, color='gray', ls='--', lw=1)\n    ax2.yaxis.grid(True, linestyle='--', alpha=0.4)\n    ax2.set_axisbelow(True)\n    plt.tight_layout()\n    save_fig('fig_solve_rate_comparison.png')\n    plt.show()\n\n    print('\\n--- Final manuscript numbers (copy into Table 10) ---')\n    print(wls_summary2_df.to_string(index=False))\nelse:\n    print('No wls_summary2_df — run Cell K first.')","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.181623Z","iopub.status.idle":"2026-06-05T22:04:19.181830Z","shell.execute_reply.started":"2026-06-05T22:04:19.181729Z","shell.execute_reply":"2026-06-05T22:04:19.181741Z"}},"outputs":[],"execution_count":null},{"cell_type":"code","source":"","metadata":{"trusted":true,"execution":{"iopub.status.busy":"2026-06-05T22:04:19.183070Z","iopub.status.idle":"2026-06-05T22:04:19.183435Z","shell.execute_reply.started":"2026-06-05T22:04:19.183236Z","shell.execute_reply":"2026-06-05T22:04:19.183265Z"}},"outputs":[],"execution_count":null}]}