1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
// SPDX-License-Identifier: AGPL-3.0-only
//! Real-hardware numeric-parity check of the Allan-family estimators against a
//! **measured atomic clock**, not a synthetic test vector.
//!
//! Where `allan_reference.rs` / `allan_nist_sp1065_1000point.rs` validate the
//! estimators on the canonical *synthetic* NBS14 data, this island validates them
//! on **556 990 phase samples of a real 5071A caesium primary standard measured
//! against a hydrogen maser** (1 PPS into a 53230A time-interval counter, τ₀ = 1 s,
//! Feb 2014; collected by A. Wallin, distributed with `allantools`).
//!
//! The reference deviations are the ones **Stable32** (W. J. Riley; the de-facto
//! reference frequency-stability tool) computes for the same file, with its 1‑σ
//! (0.683) confidence band — committed under `tests/fixtures/cs5071a/`
//! (`oadev_decade.txt`, `ohdev_decade.txt`; see the `NOTICE.md` there).
//!
//! The raw phase file is third-party data without an explicit redistribution
//! licence, so it is **git-ignored**, not vendored. Fetch it with
//! `scripts/fetch_cs5071a.sh` (downloads into `realdata-cache/cs5071a/`); override
//! the location with `KSHANA_CS5071A_PATH`. When the raw file is absent the test
//! prints a skip notice and passes, so CI without the data stays green; when it is
//! present the estimators are checked against the Stable32 reference across the
//! whole decade τ ladder. Only the public reference numbers are committed; no
//! third-party code is used.
use kshana::allan::{hadamard_adev, overlapping_adev};
const OADEV_ORACLE: &str = include_str!("fixtures/cs5071a/oadev_decade.txt");
const OHDEV_ORACLE: &str = include_str!("fixtures/cs5071a/ohdev_decade.txt");
/// One Stable32 reference row: averaging factor, point count, and the
/// `[min, sigma, max]` deviation with its confidence band.
struct OracleRow {
af: usize,
npts: u64,
min: f64,
sigma: f64,
max: f64,
}
/// Parse a Stable32 `*_decade.txt` table: comment lines start with `#`; data rows
/// are `AF Tau # Alpha MinSigma Sigma MaxSigma`.
fn parse_oracle(text: &str) -> Vec<OracleRow> {
let mut rows = Vec::new();
for line in text.lines() {
let line = line.trim();
if line.is_empty() || line.starts_with('#') {
continue;
}
let c: Vec<&str> = line.split_whitespace().collect();
if c.len() < 7 {
continue;
}
rows.push(OracleRow {
af: c[0].parse().expect("AF"),
npts: c[2].parse().expect("# count"),
min: c[4].parse().expect("MinSigma"),
sigma: c[5].parse().expect("Sigma"),
max: c[6].parse().expect("MaxSigma"),
});
}
rows
}
/// Locate the git-ignored raw phase file, if it has been fetched.
fn phase_path() -> std::path::PathBuf {
if let Ok(p) = std::env::var("KSHANA_CS5071A_PATH") {
return std::path::PathBuf::from(p);
}
std::path::Path::new(env!("CARGO_MANIFEST_DIR")).join("realdata-cache/cs5071a/5071A_phase.txt")
}
/// Read the phase series (one dimensionless phase sample per line; `#` comments).
fn load_phase() -> Option<Vec<f64>> {
let path = phase_path();
let text = std::fs::read_to_string(&path).ok()?;
let phase: Vec<f64> = text
.lines()
.map(str::trim)
.filter(|l| !l.is_empty() && !l.starts_with('#'))
.map(|l| l.parse::<f64>().expect("phase sample"))
.collect();
Some(phase)
}
fn rel_err(got: f64, want: f64) -> f64 {
((got - want) / want).abs()
}
/// Relative tolerance for well-sampled averaging factors. On this dataset Kshana
/// reproduces every Stable32 deviation to < 3e-5; 1e-3 leaves a wide margin while
/// still being a tight parity bar.
const TOL: f64 = 1e-3;
/// Averaging factors with at least this many analysis points are checked to `TOL`;
/// sparser (long-τ) rows are only required to land inside Stable32's own
/// confidence band, where the point estimate itself is statistically loose.
const WELL_SAMPLED: u64 = 100;
/// Check one estimator against its Stable32 reference table over the decade ladder.
fn check_estimator(
name: &str,
oracle: &str,
phase: &[f64],
est: impl Fn(&[f64], f64, usize) -> f64,
) {
let rows = parse_oracle(oracle);
assert!(rows.len() >= 10, "{name}: oracle table looks truncated");
let mut checked = 0;
for r in &rows {
// Need at least a few independent differences for the estimator to be defined.
if r.af * 3 >= phase.len() {
continue;
}
let got = est(phase, 1.0, r.af);
if r.npts >= WELL_SAMPLED {
let rel = rel_err(got, r.sigma);
assert!(
rel < TOL,
"{name} AF={}: got {got:.6e}, Stable32 sigma {:.6e} (rel {rel:.2e} >= {TOL:.0e})",
r.af,
r.sigma
);
} else {
// Allow a hair of slack on the reported band edges (Stable32 prints 5 sig figs).
let lo = r.min * (1.0 - 1e-3);
let hi = r.max * (1.0 + 1e-3);
assert!(
got >= lo && got <= hi,
"{name} AF={}: got {got:.6e} outside Stable32 band [{:.6e}, {:.6e}]",
r.af,
r.min,
r.max
);
}
checked += 1;
}
assert!(
checked >= 10,
"{name}: only {checked} averaging factors checked"
);
eprintln!("[cs5071a] {name}: {checked} averaging factors match Stable32");
}
#[test]
fn cs5071a_real_caesium_estimators_match_stable32() {
let Some(phase) = load_phase() else {
eprintln!(
"[cs5071a] SKIP: raw phase data not found at {} \
(run scripts/fetch_cs5071a.sh or set KSHANA_CS5071A_PATH); CI stays green.",
phase_path().display()
);
return;
};
assert_eq!(
phase.len(),
556_990,
"expected the full 556 990-point 5071A series; got {}",
phase.len()
);
check_estimator("OADEV", OADEV_ORACLE, &phase, overlapping_adev);
check_estimator("OHDEV", OHDEV_ORACLE, &phase, hadamard_adev);
}
#[test]
fn cs5071a_oracle_tables_are_well_formed() {
// Hermetic guard: the committed Stable32 reference tables parse and are sane,
// so the island is meaningful even on a machine that never fetches the raw data.
for (name, txt) in [("oadev", OADEV_ORACLE), ("ohdev", OHDEV_ORACLE)] {
let rows = parse_oracle(txt);
assert!(rows.len() >= 15, "{name}: too few rows ({})", rows.len());
// τ=1 s overlapping ADEV of this caesium standard is ~3.3e-10.
let first = &rows[0];
assert_eq!(first.af, 1, "{name}: first row should be AF=1");
assert!(
first.min <= first.sigma && first.sigma <= first.max,
"{name}: AF=1 band is not ordered"
);
for r in &rows {
assert!(
r.sigma > 0.0 && r.min > 0.0 && r.max > 0.0,
"{name}: nonpositive deviation"
);
assert!(r.min <= r.max, "{name}: AF={} band inverted", r.af);
}
}
}