fuzzyworld

shadow-zone / files / src / validator.ts

validator.ts

5.5 KB · history · blame

1import type { EarthquakeDataset, PhaseArrival } from './earthquakeData';
2import type { EarthModel } from './earth';
3import { buildTravelTimeCurve } from './rayTracer';
4
5export interface StationResidual {
6 station: string;
7 network: string;
8 distanceDeg: number;
9 observedSec: number;
10 predictedSec: number;
11 residualSec: number; // predicted − observed
12}
13
14export interface ValidationResult {
15 datasetId: string;
16 datasetName: string;
17 stationCount: number;
18 /** Root-mean-square of residuals in seconds. Lower = better fit to real Earth. */
19 rmsSec: number;
20 meanResidualSec: number;
21 /** Fraction of arrivals within ±5s of observed (rough "hit rate"). */
22 hitRate5s: number;
23 perStation: StationResidual[];
24 /** True if the dataset has no arrivals (fetch script not yet run). */
25 noData: boolean;
26}
27
28/**
29 * Interpolate predicted travel time from the curve at a given epicentral distance.
30 * Returns null if the distance is inside the model's shadow zone (no ray arrives there).
31 */
32function interpolateTravelTime(
33 curve: Array<{ distanceDeg: number; travelTimeSec: number }>,
34 distanceDeg: number,
35): number | null {
36 if (curve.length === 0) return null;
37
38 // Find bracketing points
39 let lo = 0;
40 let hi = curve.length - 1;
41
42 // Extrapolation guard
43 if (distanceDeg < curve[0].distanceDeg - 2 || distanceDeg > curve[hi].distanceDeg + 2) {
44 return null;
45 }
46
47 // Binary search for the closest bracketing interval
48 while (lo < hi - 1) {
49 const mid = (lo + hi) >> 1;
50 if (curve[mid].distanceDeg < distanceDeg) lo = mid;
51 else hi = mid;
52 }
53
54 // Check for gap (shadow zone): if the gap between lo and hi is large, no ray arrives
55 const gap = curve[hi].distanceDeg - curve[lo].distanceDeg;
56 if (gap > 5 && distanceDeg > curve[lo].distanceDeg + 0.5 && distanceDeg < curve[hi].distanceDeg - 0.5) {
57 return null; // in shadow zone
58 }
59
60 const t = (distanceDeg - curve[lo].distanceDeg) / (curve[hi].distanceDeg - curve[lo].distanceDeg + 1e-10);
61 return curve[lo].travelTimeSec + t * (curve[hi].travelTimeSec - curve[lo].travelTimeSec);
62}
63
64export function validateModel(
65 model: EarthModel,
66 dataset: EarthquakeDataset,
67): ValidationResult {
68 if (dataset.arrivals.length === 0) {
69 return {
70 datasetId: dataset.id,
71 datasetName: dataset.name,
72 stationCount: 0,
73 rmsSec: NaN,
74 meanResidualSec: NaN,
75 hitRate5s: NaN,
76 perStation: [],
77 noData: true,
78 };
79 }
80
81 const curve = buildTravelTimeCurve(model, 'P', dataset.depthKm);
82 const perStation: StationResidual[] = [];
83
84 for (const arrival of dataset.arrivals) {
85 if (arrival.pTravelTimeSec === null) continue;
86 const predicted = interpolateTravelTime(curve, arrival.distanceDeg);
87 if (predicted === null) continue; // model says shadow zone here
88
89 perStation.push({
90 station: arrival.station,
91 network: arrival.network,
92 distanceDeg: arrival.distanceDeg,
93 observedSec: arrival.pTravelTimeSec,
94 predictedSec: predicted,
95 residualSec: predicted - arrival.pTravelTimeSec,
96 });
97 }
98
99 if (perStation.length === 0) {
100 return {
101 datasetId: dataset.id,
102 datasetName: dataset.name,
103 stationCount: 0,
104 rmsSec: NaN,
105 meanResidualSec: NaN,
106 hitRate5s: NaN,
107 perStation: [],
108 noData: false,
109 };
110 }
111
112 const residuals = perStation.map(s => s.residualSec);
113 const mean = residuals.reduce((a, b) => a + b, 0) / residuals.length;
114 const rms = Math.sqrt(residuals.reduce((a, r) => a + r * r, 0) / residuals.length);
115 const hits = residuals.filter(r => Math.abs(r) <= 5).length;
116
117 return {
118 datasetId: dataset.id,
119 datasetName: dataset.name,
120 stationCount: perStation.length,
121 rmsSec: rms,
122 meanResidualSec: mean,
123 hitRate5s: hits / perStation.length,
124 perStation,
125 noData: false,
126 };
127}
128
129/**
130 * Classify a model fitness score as a human-readable label.
131 * Thresholds are generous — even PREM won't be perfect against raw ISC picks.
132 */
133export function fitnessLabel(result: ValidationResult): string {
134 if (result.noData) return '(no data — run fetch script)';
135 if (result.stationCount === 0) return 'No matched arrivals';
136 const rms = result.rmsSec;
137 if (rms < 3) return '🟢 Excellent fit';
138 if (rms < 8) return '🟡 Good fit';
139 if (rms < 20) return '🟠 Fair fit';
140 return '🔴 Poor fit';
141}
142
143/**
144 * Count observed arrivals that fall in the model's shadow zone — if the model
145 * is correct, this should be near zero. A water-core model would have many
146 * shadow-zone mismatches.
147 */
148export function countShadowZoneMismatches(
149 model: EarthModel,
150 dataset: EarthquakeDataset,
151): { inShadow: number; total: number } {
152 const curve = buildTravelTimeCurve(model, 'P', dataset.depthKm);
153 let inShadow = 0;
154 let total = 0;
155
156 for (const arrival of dataset.arrivals) {
157 if (arrival.pTravelTimeSec === null) continue;
158 total++;
159 const predicted = interpolateTravelTime(curve, arrival.distanceDeg);
160 if (predicted === null) inShadow++;
161 }
162
163 return { inShadow, total };
164}
165
166/**
167 * Build a distance → observed travel time map from dataset arrivals,
168 * for use in the renderer's scatter plot overlay.
169 */
170export function getObservedPoints(
171 dataset: EarthquakeDataset,
172 wave: 'P' | 'S',
173): Array<{ distanceDeg: number; travelTimeSec: number }> {
174 return dataset.arrivals
175 .filter((a: PhaseArrival) => {
176 const t = wave === 'P' ? a.pTravelTimeSec : a.sTravelTimeSec;
177 return t !== null;
178 })
179 .map((a: PhaseArrival) => ({
180 distanceDeg: a.distanceDeg,
181 travelTimeSec: (wave === 'P' ? a.pTravelTimeSec : a.sTravelTimeSec)!,
182 }));
183}