426 lines
15 KiB
Python
426 lines
15 KiB
Python
"""Run the complete DLL-quality/WLS and code-smoothing comparison."""
|
||
|
||
from __future__ import annotations
|
||
|
||
import argparse
|
||
import csv
|
||
import math
|
||
from concurrent.futures import ProcessPoolExecutor
|
||
from dataclasses import dataclass
|
||
from pathlib import Path
|
||
|
||
import matplotlib
|
||
|
||
matplotlib.use("Agg")
|
||
import matplotlib.pyplot as plt
|
||
import numpy as np
|
||
|
||
from .code_smoothing_model import run_smoothing_trial
|
||
from .firmware_wls_model import SigmaMethod
|
||
from .measurement_scenarios import EpochSample, run_scenario, scenario_matrix
|
||
|
||
|
||
METHOD_ORDER = (
|
||
SigmaMethod.CURRENT,
|
||
SigmaMethod.NAIVE_ADD_DLL,
|
||
SigmaMethod.DLL_EXCESS,
|
||
SigmaMethod.DLL_EXCESS_GUARD,
|
||
SigmaMethod.ORACLE,
|
||
)
|
||
METHOD_LABEL = {
|
||
SigmaMethod.CURRENT: "Mevcut C/N0",
|
||
SigmaMethod.NAIVE_ADD_DLL: "Naif DLL ekle",
|
||
SigmaMethod.DLL_EXCESS: "DLL excess",
|
||
SigmaMethod.DLL_EXCESS_GUARD: "DLL excess+mean",
|
||
SigmaMethod.ORACLE: "Oracle variance",
|
||
}
|
||
|
||
|
||
@dataclass(frozen=True)
|
||
class Summary:
|
||
scenario: str
|
||
method: SigmaMethod
|
||
epochs: int
|
||
horizontal_rms_m: float
|
||
position_3d_rms_m: float
|
||
position_3d_p95_m: float
|
||
valid_percent: float
|
||
valid_3d_rms_m: float
|
||
horizontal_2sigma_coverage_percent: float
|
||
false_valid_over_25m_percent: float
|
||
fde_percent: float
|
||
mean_measurement_sigma_m: float
|
||
|
||
|
||
def _worker(task: tuple[int, int]) -> tuple[int, int, list[EpochSample]]:
|
||
scenario_index, seed = task
|
||
scenario = scenario_matrix()[scenario_index]
|
||
return scenario_index, seed, run_scenario(scenario, seed, METHOD_ORDER)
|
||
|
||
|
||
def _rms(values: np.ndarray) -> float:
|
||
return float(np.sqrt(np.mean(np.square(values)))) if len(values) else math.nan
|
||
|
||
|
||
def summarize(scenario_name: str, method: SigmaMethod, samples: list[EpochSample]) -> Summary:
|
||
selected = [
|
||
item
|
||
for item in samples
|
||
if item.method == method and item.time_s >= 5.0
|
||
]
|
||
horizontal = np.asarray([item.horizontal_error_m for item in selected])
|
||
error_3d = np.asarray([item.position_error_3d_m for item in selected])
|
||
valid = np.asarray([item.valid for item in selected], dtype=bool)
|
||
horizontal_sigma = np.asarray([item.horizontal_sigma_m for item in selected])
|
||
fde = np.asarray([item.fde_applied for item in selected], dtype=bool)
|
||
mean_sigma = np.asarray([item.mean_measurement_sigma_m for item in selected])
|
||
false_valid = valid & (error_3d > 25.0)
|
||
return Summary(
|
||
scenario=scenario_name,
|
||
method=method,
|
||
epochs=len(selected),
|
||
horizontal_rms_m=_rms(horizontal),
|
||
position_3d_rms_m=_rms(error_3d),
|
||
position_3d_p95_m=float(np.percentile(error_3d, 95.0)),
|
||
valid_percent=100.0 * float(np.mean(valid)),
|
||
valid_3d_rms_m=_rms(error_3d[valid]),
|
||
horizontal_2sigma_coverage_percent=100.0
|
||
* float(np.mean(horizontal <= 2.0 * horizontal_sigma)),
|
||
false_valid_over_25m_percent=100.0 * float(np.mean(false_valid)),
|
||
fde_percent=100.0 * float(np.mean(fde)),
|
||
mean_measurement_sigma_m=float(np.mean(mean_sigma)),
|
||
)
|
||
|
||
|
||
def write_csv(path: Path, summaries: list[Summary]) -> None:
|
||
with path.open("w", newline="", encoding="utf-8") as output:
|
||
writer = csv.DictWriter(output, fieldnames=list(Summary.__dataclass_fields__))
|
||
writer.writeheader()
|
||
for item in summaries:
|
||
row = dict(item.__dict__)
|
||
row["method"] = item.method.value
|
||
writer.writerow(row)
|
||
|
||
|
||
def _comparison_plot(
|
||
path: Path,
|
||
summaries: list[Summary],
|
||
field: str,
|
||
ylabel: str,
|
||
*,
|
||
logarithmic: bool = False,
|
||
) -> None:
|
||
scenarios = scenario_matrix()
|
||
x = np.arange(len(scenarios))
|
||
width = 0.16
|
||
figure, axis = plt.subplots(figsize=(15, 6.5))
|
||
for method_index, method in enumerate(METHOD_ORDER):
|
||
values = [
|
||
getattr(
|
||
next(
|
||
item
|
||
for item in summaries
|
||
if item.scenario == scenario.name and item.method == method
|
||
),
|
||
field,
|
||
)
|
||
for scenario in scenarios
|
||
]
|
||
axis.bar(
|
||
x + (method_index - 2) * width,
|
||
values,
|
||
width,
|
||
label=METHOD_LABEL[method],
|
||
)
|
||
axis.set_ylabel(ylabel)
|
||
if logarithmic:
|
||
axis.set_yscale("log")
|
||
axis.set_xticks(x)
|
||
axis.set_xticklabels([item.name for item in scenarios], rotation=28, ha="right")
|
||
axis.grid(axis="y", alpha=0.25)
|
||
axis.legend(ncol=3)
|
||
figure.tight_layout()
|
||
figure.savefig(path, dpi=150)
|
||
plt.close(figure)
|
||
|
||
|
||
def _transition_plot(path: Path, samples: list[EpochSample]) -> None:
|
||
figure, axes = plt.subplots(2, 1, figsize=(12, 7), sharex=True)
|
||
for method in (
|
||
SigmaMethod.CURRENT,
|
||
SigmaMethod.DLL_EXCESS,
|
||
SigmaMethod.DLL_EXCESS_GUARD,
|
||
):
|
||
selected = [item for item in samples if item.method == method]
|
||
axes[0].plot(
|
||
[item.time_s for item in selected],
|
||
[item.position_error_3d_m for item in selected],
|
||
label=METHOD_LABEL[method],
|
||
linewidth=1.2,
|
||
)
|
||
reference = [item for item in samples if item.method == SigmaMethod.CURRENT]
|
||
axes[1].plot(
|
||
[item.time_s for item in reference],
|
||
[item.maximum_dll_ratio for item in reference],
|
||
color="tab:red",
|
||
)
|
||
axes[0].axvspan(22.0, 45.0, color="tab:orange", alpha=0.15)
|
||
axes[1].axvspan(22.0, 45.0, color="tab:orange", alpha=0.15)
|
||
axes[0].set_ylabel("3D konum hatası (m)")
|
||
axes[1].set_ylabel("Maks. DLL variance ratio")
|
||
axes[1].set_xlabel("Zaman (s)")
|
||
axes[0].legend()
|
||
for axis in axes:
|
||
axis.grid(alpha=0.25)
|
||
figure.tight_layout()
|
||
figure.savefig(path, dpi=150)
|
||
plt.close(figure)
|
||
|
||
|
||
def smoothing_summary(seed_count: int) -> dict[str, dict[str, float]]:
|
||
result: dict[str, dict[str, float]] = {}
|
||
for label, cfo_mps in (
|
||
("cfo_compensated", 0.0),
|
||
("cfo_residual_0p5_mps", 0.5),
|
||
("cfo_residual_3_mps", 3.0),
|
||
):
|
||
trials = [
|
||
run_smoothing_trial(seed, doppler_cfo_mps=cfo_mps)
|
||
for seed in range(1, seed_count + 1)
|
||
]
|
||
result[label] = {
|
||
field: float(np.mean([getattr(item, field) for item in trials]))
|
||
for field in trials[0].__dataclass_fields__
|
||
}
|
||
return result
|
||
|
||
|
||
def _smoothing_plot(path: Path, smoothing: dict[str, dict[str, float]]) -> None:
|
||
labels = list(smoothing)
|
||
x = np.arange(len(labels))
|
||
width = 0.20
|
||
figure, axis = plt.subplots(figsize=(10, 5.5))
|
||
for offset, (field, title) in enumerate(
|
||
(
|
||
("raw_rmse_m", "Ham code"),
|
||
("scalar_iir_rmse_m", "Kör IIR"),
|
||
("doppler_smoothed_rmse_m", "Raw Doppler-predicted"),
|
||
(
|
||
"firmware_iir_doppler_smoothed_rmse_m",
|
||
"FW-IIR Doppler-predicted",
|
||
),
|
||
)
|
||
):
|
||
axis.bar(
|
||
x + (offset - 1.5) * width,
|
||
[smoothing[label][field] for label in labels],
|
||
width,
|
||
label=title,
|
||
)
|
||
axis.set_xticks(x)
|
||
axis.set_xticklabels(labels)
|
||
axis.set_ylabel("Range RMSE (m)")
|
||
axis.set_yscale("log")
|
||
axis.grid(axis="y", alpha=0.25)
|
||
axis.legend()
|
||
figure.tight_layout()
|
||
figure.savefig(path, dpi=150)
|
||
plt.close(figure)
|
||
|
||
|
||
def write_report(
|
||
path: Path,
|
||
summaries: list[Summary],
|
||
focused_rms: dict[str, dict[SigmaMethod, float]],
|
||
smoothing: dict[str, dict[str, float]],
|
||
seed_count: int,
|
||
) -> None:
|
||
current = {
|
||
item.scenario: item
|
||
for item in summaries
|
||
if item.method == SigmaMethod.CURRENT
|
||
}
|
||
lines = [
|
||
"# GPS code measurement model sonuçları",
|
||
"",
|
||
f"Monte Carlo seed sayısı: `{seed_count}`. İlk 5 saniye metriklere alınmadı.",
|
||
"",
|
||
"## Yöntem karşılaştırması",
|
||
"",
|
||
"| Senaryo | Yöntem | 3D RMS (m) | Mevcuda göre | P95 (m) | Valid (%) | >25 m false-valid (%) |",
|
||
"|---|---|---:|---:|---:|---:|---:|",
|
||
]
|
||
for scenario in scenario_matrix():
|
||
baseline = current[scenario.name].position_3d_rms_m
|
||
for method in METHOD_ORDER:
|
||
item = next(
|
||
value
|
||
for value in summaries
|
||
if value.scenario == scenario.name and value.method == method
|
||
)
|
||
improvement = 100.0 * (baseline - item.position_3d_rms_m) / baseline
|
||
lines.append(
|
||
f"| {scenario.name} | {METHOD_LABEL[method]} | "
|
||
f"{item.position_3d_rms_m:.2f} | {improvement:+.1f}% | "
|
||
f"{item.position_3d_p95_m:.2f} | {item.valid_percent:.1f} | "
|
||
f"{item.false_valid_over_25m_percent:.2f} |"
|
||
)
|
||
|
||
lines.extend(
|
||
[
|
||
"",
|
||
"## Modelden çıkan karar özeti",
|
||
"",
|
||
"- C/N0 ile uyumlu nominal DLL jitter'da `dll_excess`, mevcut çözümü değiştirmemelidir.",
|
||
"- Bütün measurement sigmaları aynı sabitle çarpılırsa WLS konum noktası değişmez; yalnız covariance/integrity değişir. Fakat ortak DLL jitter çarpanı, C/N0/elevation bileşimi kanaldan kanala farklı olduğu için toplam sigmaları eşit oranda büyütmez.",
|
||
"- Kanal-bağımlı ve C/N0'da görünmeyen excess jitter varsa dinamik sigma o kanalı down-weight ederek konumu iyileştirebilir.",
|
||
"- Stable multipath bias düşük DLL variance bırakabildiği için variance-only yöntem çözüm değildir.",
|
||
"- 50 Hz'de DLL ratio alpha=0.10 yaklaşık 0.19 s zaman sabiti ve 0.57 s yüzde-95 geçiş verir.",
|
||
"- Aşağıdaki yüzdeler seçilen sentetik hata büyüklüklerine koşulludur; kart için kazanç tahmini değildir.",
|
||
"",
|
||
"## Bozulma penceresi odaklı 3D RMS",
|
||
"",
|
||
"| Senaryo/pencere | Mevcut (m) | DLL excess (m) | DLL excess+mean (m) | Oracle variance (m) |",
|
||
"|---|---:|---:|---:|---:|",
|
||
]
|
||
)
|
||
for label, values in focused_rms.items():
|
||
lines.append(
|
||
f"| {label} | {values[SigmaMethod.CURRENT]:.2f} | "
|
||
f"{values[SigmaMethod.DLL_EXCESS]:.2f} | "
|
||
f"{values[SigmaMethod.DLL_EXCESS_GUARD]:.2f} | "
|
||
f"{values[SigmaMethod.ORACLE]:.2f} |"
|
||
)
|
||
lines.extend(
|
||
[
|
||
"",
|
||
"## Code smoothing ve carrier CFO",
|
||
"",
|
||
"| Carrier-rate koşulu | Ham RMSE (m) | Kör IIR RMSE (m) | Raw Doppler RMSE (m) | FW-IIR Doppler RMSE (m) | Kör IIR peak (m) | Raw Doppler peak (m) | FW-IIR Doppler peak (m) |",
|
||
"|---|---:|---:|---:|---:|---:|---:|---:|",
|
||
]
|
||
)
|
||
for label, values in smoothing.items():
|
||
lines.append(
|
||
f"| {label} | {values['raw_rmse_m']:.2f} | "
|
||
f"{values['scalar_iir_rmse_m']:.2f} | "
|
||
f"{values['doppler_smoothed_rmse_m']:.2f} | "
|
||
f"{values['firmware_iir_doppler_smoothed_rmse_m']:.2f} | "
|
||
f"{values['scalar_iir_peak_dynamic_error_m']:.2f} | "
|
||
f"{values['doppler_peak_dynamic_error_m']:.2f} | "
|
||
f"{values['firmware_iir_doppler_peak_dynamic_error_m']:.2f} |"
|
||
)
|
||
|
||
lines.extend(
|
||
[
|
||
"",
|
||
"## Yorum sınırları",
|
||
"",
|
||
"- `current_cn0`, firmware'deki URA/elevation/C/N0 sigma modelidir.",
|
||
"- `naive_add_dll`, aynı thermal jitter'ı iki kere saymanın negatif kontrolüdür.",
|
||
"- `dll_excess`, yalnız ölçülen DLL variance'ın C/N0 beklentisini aşan kısmını code variance'a ekler.",
|
||
"- `dll_excess_guard`, variance'a ek olarak sürekli discriminator mean stress'ini kullanır.",
|
||
"- `oracle_variance`, sentetik gerçek jitter çarpanını bilir; yalnız variance kaynakları için erişilebilir üst sınırdır, bias'ı bilmez.",
|
||
"- Stable multipath DLL'nin yeni denge noktasında düşük residual variance bırakabilir. Bu nedenle variance-only yönteminin onu çözmesi beklenmez.",
|
||
"- Discriminator gürültü ölçeği kart correlator kayıtları gelene kadar kalibrasyon varsayımıdır; oran tabanlı yöntem mutlak ölçeğe göre daha az hassastır.",
|
||
"- Sonuçlar bitstream veya firmware değişikliği değildir; karar öncesi executable tasarım keşfidir.",
|
||
"",
|
||
]
|
||
)
|
||
path.write_text("\n".join(lines), encoding="utf-8")
|
||
|
||
|
||
def main() -> None:
|
||
parser = argparse.ArgumentParser()
|
||
parser.add_argument("--seeds", type=int, default=12)
|
||
parser.add_argument("--jobs", type=int, default=4)
|
||
parser.add_argument(
|
||
"--output",
|
||
type=Path,
|
||
default=Path(__file__).resolve().parent / "output",
|
||
)
|
||
arguments = parser.parse_args()
|
||
arguments.output.mkdir(parents=True, exist_ok=True)
|
||
|
||
scenarios = scenario_matrix()
|
||
tasks = [
|
||
(scenario_index, seed)
|
||
for scenario_index in range(len(scenarios))
|
||
for seed in range(1, arguments.seeds + 1)
|
||
]
|
||
grouped: dict[int, list[EpochSample]] = {
|
||
index: [] for index in range(len(scenarios))
|
||
}
|
||
transition_samples: list[EpochSample] = []
|
||
with ProcessPoolExecutor(max_workers=max(1, arguments.jobs)) as executor:
|
||
for scenario_index, seed, samples in executor.map(_worker, tasks):
|
||
grouped[scenario_index].extend(samples)
|
||
if scenarios[scenario_index].name == "noise_burst_recovery" and seed == 1:
|
||
transition_samples = samples
|
||
print(
|
||
f"completed {scenarios[scenario_index].name} seed={seed}",
|
||
flush=True,
|
||
)
|
||
|
||
summaries = [
|
||
summarize(scenario.name, method, grouped[index])
|
||
for index, scenario in enumerate(scenarios)
|
||
for method in METHOD_ORDER
|
||
]
|
||
focused_windows = {
|
||
"noise_burst_recovery [22,45) s": ("noise_burst_recovery", 22.0, 45.0),
|
||
"dll_transient_stress [25,34) s": ("dll_transient_stress", 25.0, 34.0),
|
||
"cn0_step_recovery [5,30) s": ("cn0_step_recovery", 5.0, 30.0),
|
||
}
|
||
focused_rms: dict[str, dict[SigmaMethod, float]] = {}
|
||
for label, (scenario_name, start_s, stop_s) in focused_windows.items():
|
||
scenario_index = next(
|
||
index
|
||
for index, scenario in enumerate(scenarios)
|
||
if scenario.name == scenario_name
|
||
)
|
||
focused_rms[label] = {}
|
||
for method in METHOD_ORDER:
|
||
values = np.asarray(
|
||
[
|
||
item.position_error_3d_m
|
||
for item in grouped[scenario_index]
|
||
if item.method == method
|
||
and start_s <= item.time_s < stop_s
|
||
]
|
||
)
|
||
focused_rms[label][method] = _rms(values)
|
||
smoothing = smoothing_summary(max(arguments.seeds, 20))
|
||
|
||
write_csv(arguments.output / "scenario_summary.csv", summaries)
|
||
write_report(
|
||
arguments.output / "ANALYSIS.md",
|
||
summaries,
|
||
focused_rms,
|
||
smoothing,
|
||
arguments.seeds,
|
||
)
|
||
_comparison_plot(
|
||
arguments.output / "01_position_3d_rms.png",
|
||
summaries,
|
||
"position_3d_rms_m",
|
||
"3D RMS (m)",
|
||
logarithmic=True,
|
||
)
|
||
_comparison_plot(
|
||
arguments.output / "02_solution_availability.png",
|
||
summaries,
|
||
"valid_percent",
|
||
"Integrity-valid çözüm (%)",
|
||
)
|
||
_transition_plot(
|
||
arguments.output / "03_noise_burst_transition.png",
|
||
transition_samples,
|
||
)
|
||
_smoothing_plot(arguments.output / "04_code_smoothing_cfo.png", smoothing)
|
||
print(f"wrote analysis to {arguments.output}", flush=True)
|
||
|
||
|
||
if __name__ == "__main__":
|
||
main()
|