#!/usr/bin/env python3 """Embedded-oriented PPS closed-loop scenario and Monte-Carlo analysis. This is intentionally deterministic and dependency-free. It exercises the same Q32/Q46 controller, quality hysteresis, Doppler-aided software pre-filter, runtime profile changes, lock detector and 125 MHz edge quantization used by the implementation plan. """ from __future__ import annotations from dataclasses import dataclass, replace import math import random import statistics from typing import Callable try: from .adaptive_pps import ( Observation, Profile, QualityMetrics, run_observations, ) except ImportError: from adaptive_pps import ( Observation, Profile, QualityMetrics, run_observations, ) @dataclass(frozen=True) class Scenario: name: str duration_s: float seeds: int oscillator_ppm: float noise_sigma: Callable[[float], float] quality_profile: Callable[[float], Profile] gaps: tuple[tuple[float, float], ...] = () timing_step: Callable[[float], float] = lambda _time_s: 0.0 outlier_period: int = 0 outlier_size_s: float = 0.0 drift_sigma_mps: float = 0.25 drift_valid: bool = True drift_valid_profile: Callable[[float], bool] | None = None drop_period: int = 0 invalid_intervals: tuple[tuple[float, float], ...] = () reference_range_rate: Callable[[float], float] = lambda _time_s: 0.0 reference_id: Callable[[float], int] = lambda _time_s: 1 reference_datum_offset: Callable[[float], float] = ( lambda _time_s: 0.0 ) @dataclass(frozen=True) class ScenarioResult: name: str lock_success: int runs: int first_lock_p95_s: float settle_p95_s: float phase_rms_median_ns: float edge_jitter_median_ns: float edge_jitter_p95_ns: float unlock_total: int coarse_total: int rejected_total: int def quality(profile: Profile) -> QualityMetrics: if profile == Profile.PRECISION: return QualityMetrics( 1.7, 3.2, 18.0, 38.0, 27.0, 0.35, 0.25, 7 ) if profile == Profile.ROBUST: return QualityMetrics( 5.0, 8.0, 65.0, 27.0, 18.0, 0.85, 1.8, 5 ) return QualityMetrics( 12.0, 18.0, 180.0, 18.0, 10.0, 2.0, 4.0, 4 ) def _percentile(values: list[float], percentile: float) -> float: if not values: return math.nan ordered = sorted(values) position = (len(ordered) - 1) * percentile lower = math.floor(position) upper = math.ceil(position) if lower == upper: return ordered[lower] fraction = position - lower return ordered[lower] * (1.0 - fraction) + ordered[upper] * fraction def observations(scenario: Scenario, seed: int) -> list[Observation]: rng = random.Random(seed) result: list[Observation] = [] updates = round(scenario.duration_s / 0.2) for update in range(updates): time_s = update * 0.2 if any(start <= time_s < stop for start, stop in scenario.gaps): continue if scenario.drop_period and update % scenario.drop_period == 0: continue timing_error_s = ( rng.gauss(0.0, scenario.noise_sigma(time_s)) + scenario.timing_step(time_s) ) if scenario.outlier_period and update % scenario.outlier_period == 0: timing_error_s += rng.choice((-1.0, 1.0)) * scenario.outlier_size_s quality_metrics = quality(scenario.quality_profile(time_s)) solution_valid = not any( start <= time_s < stop for start, stop in scenario.invalid_intervals ) if not solution_valid: quality_metrics = replace( quality_metrics, solution_valid=False, used_satellites=2, velocity_valid=False, ) drift_valid = ( scenario.drift_valid_profile(time_s) if scenario.drift_valid_profile is not None else scenario.drift_valid ) if not drift_valid: quality_metrics = replace( quality_metrics, velocity_valid=False, velocity_residual_rms_mps=0.0, ) result.append( Observation( time_s=time_s, timing_error_s=timing_error_s, clock_drift_error_mps=rng.gauss( 0.0, scenario.drift_sigma_mps ), quality=quality_metrics, clock_drift_valid=drift_valid, reference_range_rate_mps=scenario.reference_range_rate( time_s ), reference_id=scenario.reference_id(time_s), reference_datum_offset_m=scenario.reference_datum_offset( time_s ), ) ) return result def run_scenario(scenario: Scenario) -> ScenarioResult: traces = [ run_observations( observations(scenario, seed), finish_time_s=scenario.duration_s, oscillator_ppm=scenario.oscillator_ppm, ) for seed in range(1, scenario.seeds + 1) ] first_locks = [ value for trace in traces if (value := trace.first_lock_time()) is not None ] settle_times = [ value for trace in traces if ( value := trace.settling_time( rms_limit_s=100.0e-9, window_s=10.0 ) ) is not None ] analysis_start_s = max(30.0, scenario.duration_s * 0.70) phase_rms = [ value * 1.0e9 for trace in traces if math.isfinite( value := trace.phase_rms(after_s=analysis_start_s) ) ] edge_jitter = [ value * 1.0e9 for trace in traces if math.isfinite( value := trace.one_second_edge_jitter_rms( after_s=analysis_start_s ) ) ] return ScenarioResult( name=scenario.name, lock_success=len(first_locks), runs=len(traces), first_lock_p95_s=_percentile(first_locks, 0.95), settle_p95_s=_percentile(settle_times, 0.95), phase_rms_median_ns=( statistics.median(phase_rms) if phase_rms else math.nan ), edge_jitter_median_ns=( statistics.median(edge_jitter) if edge_jitter else math.nan ), edge_jitter_p95_ns=_percentile(edge_jitter, 0.95), unlock_total=sum(trace.unlock_events for trace in traces), coarse_total=sum(trace.coarse_commits for trace in traces), rejected_total=sum(trace.rejected_observations for trace in traces), ) def scenario_matrix() -> tuple[Scenario, ...]: precision = lambda _time_s: Profile.PRECISION robust = lambda _time_s: Profile.ROBUST return ( Scenario( "clean_precision", 120.0, 30, 1.65, lambda _time_s: 50.0e-9, precision, ), Scenario( "board_noise_outliers", 180.0, 50, 1.65, lambda _time_s: 320.0e-9, precision, outlier_period=41, outlier_size_s=900.0e-9, ), Scenario( "weak_geometry_robust", 240.0, 30, 1.65, lambda _time_s: 650.0e-9, robust, outlier_period=37, outlier_size_s=1.2e-6, drift_sigma_mps=0.5, ), Scenario( "short_dropouts", 180.0, 30, -3.0, lambda _time_s: 320.0e-9, precision, drop_period=7, ), Scenario( "long_wls_gap", 180.0, 30, 1.65, lambda _time_s: 320.0e-9, precision, gaps=((60.0, 75.0),), ), Scenario( "quality_transitions", 240.0, 30, 2.5, lambda time_s: ( 320.0e-9 if time_s < 70.0 or time_s >= 150.0 else 650.0e-9 ), lambda time_s: ( Profile.PRECISION if time_s < 70.0 or time_s >= 150.0 else ( Profile.ROBUST if time_s < 120.0 else Profile.HOLDOVER ) ), outlier_period=43, outlier_size_s=900.0e-9, ), Scenario( "startup_noisy_then_clean", 240.0, 30, 1.65, lambda time_s: 1.1e-6 if time_s < 55.0 else 120.0e-9, lambda time_s: ( Profile.ROBUST if time_s < 55.0 else Profile.PRECISION ), outlier_period=47, outlier_size_s=1.4e-6, drift_sigma_mps=0.5, ), Scenario( "noise_bursts_then_recover", 300.0, 30, -2.5, lambda time_s: ( 900.0e-9 if ( 55.0 <= time_s < 78.0 or 135.0 <= time_s < 151.0 or 205.0 <= time_s < 214.0 ) else 150.0e-9 ), lambda time_s: ( Profile.ROBUST if ( 55.0 <= time_s < 78.0 or 135.0 <= time_s < 151.0 or 205.0 <= time_s < 214.0 ) else Profile.PRECISION ), outlier_period=53, outlier_size_s=1.2e-6, ), Scenario( "satellite_loss_recovery", 260.0, 30, 1.65, lambda _time_s: 320.0e-9, precision, invalid_intervals=((55.0, 78.0), (145.0, 151.0)), ), Scenario( "doppler_loss_recovery", 240.0, 30, 1.65, lambda _time_s: 320.0e-9, precision, drift_valid_profile=lambda time_s: not ( 65.0 <= time_s < 105.0 ), ), Scenario( "reference_range_rate", 180.0, 30, 1.65, lambda _time_s: 320.0e-9, precision, reference_range_rate=lambda _time_s: 350.0, ), Scenario( "reference_switch", 220.0, 30, 1.65, lambda _time_s: 320.0e-9, precision, reference_range_rate=lambda time_s: ( 180.0 if time_s < 70.0 else -90.0 ), reference_id=lambda time_s: 3 if time_s < 70.0 else 8, reference_datum_offset=lambda time_s: ( 0.0 if time_s < 70.0 else 38000.0 ), ), Scenario( "multipath_step", 240.0, 30, 1.65, lambda _time_s: 320.0e-9, precision, timing_step=lambda time_s: ( 300.0e-9 if 70.0 <= time_s < 120.0 else 0.0 ), ), Scenario( "no_doppler_drift", 180.0, 30, 1.65, lambda _time_s: 320.0e-9, precision, drift_valid=False, ), Scenario( "near_slew_limit", 240.0, 30, 18.0, lambda _time_s: 320.0e-9, precision, ), Scenario( "beyond_slew_limit", 120.0, 10, 25.0, lambda _time_s: 50.0e-9, precision, ), ) def validate(results: list[ScenarioResult]) -> None: by_name = {result.name: result for result in results} for name in ( "clean_precision", "board_noise_outliers", "weak_geometry_robust", "short_dropouts", "long_wls_gap", "quality_transitions", "startup_noisy_then_clean", "noise_bursts_then_recover", "satellite_loss_recovery", "doppler_loss_recovery", "reference_range_rate", "reference_switch", "multipath_step", "no_doppler_drift", "near_slew_limit", ): result = by_name[name] if result.lock_success != result.runs: raise AssertionError(f"{name}: not every run acquired lock") if by_name["clean_precision"].first_lock_p95_s >= 15.0: raise AssertionError("clean acquisition exceeded 15 s") if by_name["board_noise_outliers"].edge_jitter_p95_ns >= 20.0: raise AssertionError("board-noise p95 edge jitter exceeded 20 ns") if by_name["near_slew_limit"].settle_p95_s >= 120.0: raise AssertionError("near-limit oscillator did not settle in 120 s") if by_name["beyond_slew_limit"].lock_success != 0: raise AssertionError("out-of-range oscillator falsely claimed lock") def main() -> int: results = [run_scenario(scenario) for scenario in scenario_matrix()] print( "scenario lock lock_p95 settle_p95 " "phase_rms edge_jit_med edge_jit_p95 unlock coarse reject" ) for result in results: print( f"{result.name:27s} " f"{result.lock_success:2d}/{result.runs:<2d} " f"{result.first_lock_p95_s:8.2f}s " f"{result.settle_p95_s:8.2f}s " f"{result.phase_rms_median_ns:8.2f}ns " f"{result.edge_jitter_median_ns:8.2f}ns " f"{result.edge_jitter_p95_ns:8.2f}ns " f"{result.unlock_total:6d} " f"{result.coarse_total:6d} " f"{result.rejected_total:6d}" ) validate(results) print("PPS scenario matrix passed") return 0 if __name__ == "__main__": raise SystemExit(main())