#!/usr/bin/env python3 """airquality-evaluate.py - recompute the ozone design-value chain of 40 CFR Part 50 Appendix U from filed hourly readings, and recount the paper's headline census numbers from exact values. The paper: "More than the last digit" (airquality.pdf on the same site). Modes ----- default Rebuild the regulation's whole chain from the filed hourly readings in hourly-o3-exhibits.csv.gz (beside this script) for the three monitors worked in the paper, under four digit conventions, and check every value quoted on the page: law truncate 8-hour window means to 3 decimals (U 1.4), truncate the design value (U 1.7) all-round round-half-up at both digit sites 2010 exact window means carried through, one round-half-up at the end (EPA's proposed fix, 75 FR 3032) exact-trunc exact window means carried through, one truncation at the end Monitors: AQS site 17-031-4007 (Cook County IL, 2020-2022), 48-257-0005 (Kaufman County TX, 2022-2024), 09-009-0027 (New Haven County CT, 2020-2022 / 2021-2023 / 2023-2025). --census Recount every headline number of the paper's borderline census from the exact per-record values in census-terminal.json and census-wholechain.json (beside this script). The digit rules are re-applied to each record's exact three-year mean here, in this script; stored verdict labels are not shipped in those files and are never consulted. --check Run the default computation, then mutation controls that MUST fail (a control that passes is itself a failure, exit 1): (1) one filed hourly reading raised by 0.001 ppm inside the Cook County 2020 fourth-highest day must move at least one chain's exact mean; (2) one flipped byte in the gzipped data slice must break its sha256 comparison; (3) replacing the law chain's window-step truncation by rounding must break the Cook County law-chain check (0.070 ATTAIN); (4) moving one census record's exact law mean from 211/3000 to 214/3000 ppm must move the recounts exactly as the digit rules dictate. Every comparison is enforced (nonzero exit on any mismatch), not just printed. Python 3.8+ standard library only; no network, no other files. Measured runs: default and --census each finish in a few seconds on ordinary hardware, the full --check in well under a minute. Data provenance --------------- hourly-o3-exhibits.csv.gz was cut from EPA's AirData annual hourly files hourly_44201_{2020..2026}.zip (https://aqs.epa.gov/aqsweb/airdata/), archived and sha256-pinned at retrieval; the parent-file sha256 values are listed in the slice's own header, and the slice's sha256 is pinned below. Readings are verbatim as filed (ppm strings). Each span includes the first day of the following year because Appendix U's 8-hour windows cross midnight on Dec 31. For these three site-vintages every archived daily row carries Event Type "None" (no exceptional-event exclusions) and each site files a single monitor (POC 1), so the chain is exactly: filed hourly readings -> 8-hour window means -> daily maxima -> annual fourth-highest -> three-year mean -> design value, with no designation machinery. The zip-direct recompute of all five records reproduced the paper's certified values before this slice shipped, including Kaufman County's 857/12000 ppm. census-terminal.json and census-wholechain.json are the paper's certified census records (terminal-digit census generated 2026-07-18; whole-chain census generated 2026-08-13) with stored verdict and flip labels stripped; each file's meta block names its EPA inputs and their sha256 pins. """ import argparse import csv import gzip import hashlib import io import json import os import sys from datetime import date, timedelta from fractions import Fraction HERE = os.path.dirname(os.path.abspath(__file__)) SLICE_NAME = "hourly-o3-exhibits.csv.gz" SLICE_SHA256 = "91e177c6ef62ba862ce512cd1cf8c26de75a4a1c4abe0cdc217f0624c7d2d890" TERMINAL_NAME = "census-terminal.json" TERMINAL_SHA256 = "fc08fc548241398c5bfc17f4557944dd873eacd338d38437ea17dc49f1cf3ccd" WHOLECHAIN_NAME = "census-wholechain.json" WHOLECHAIN_SHA256 = "0980fb31c383bc5146da5d9d7d47f00aafbd09f59ce79958429b18683062f39a" # ------------------------------------------------------------------ small helpers def check(cond, msg): if not cond: print(f"FAIL: {msg}", flush=True) raise SystemExit(1) def sha256_bytes(b): return hashlib.sha256(b).hexdigest() def milli_from_str(s): """'0.052' -> 52 by string surgery (truncation onto the filing lattice, U 1.1).""" s = s.strip() neg = s.startswith("-") if s[0] in "+-": s = s[1:] if "e" in s.lower(): raise ValueError(s) ip, _, fp = s.partition(".") fp = (fp + "000")[:3] v = int(ip or "0") * 1000 + int(fp or "0") return -v if neg else v def fr(s): n, d = s.split("/") return Fraction(int(n), int(d)) def tz(f): """Truncate a nonnegative Fraction to an integer (toward zero).""" f = Fraction(f) return f.numerator // f.denominator def rhu(f): """Round a Fraction half-up to an integer.""" f = Fraction(f) return (2 * f.numerator + f.denominator) // (2 * f.denominator) def milli_to_ppm_str(f): """Exact milli-ppm Fraction -> 'num/den' in ppm units, keeping the milli fraction's own reduced terms over the 1000 (the paper's rendering: 211/3 milli-ppm prints as 211/3000, not 53/750-style re-reduction).""" f = Fraction(f) return f"{f.numerator}/{f.denominator * 1000}" def m4_str(x): """Fourth-highest value -> display string (integers as 0.072-style ppm).""" if isinstance(x, int): return f"0.{x:03d}" return milli_to_ppm_str(x) def det(dv_milli): return "ATTAIN" if dv_milli <= 70 else "VIOLATE" def iso_next(d): return (date.fromisoformat(d) + timedelta(days=1)).isoformat() def iso_prev(d): return (date.fromisoformat(d) - timedelta(days=1)).isoformat() # ------------------------------------------------------------------ the chain def day_chains(hours, dstr): """One day's maxima over the 17 begin-hours 07-23 LST (windows cross midnight). Returns (n_valid_windows, law_max, roundhalfup_max, exact_max_Fraction) with window validity 6-of-8 (available-hours divisor, U 1.2) or the zero-substitution sum threshold (U 1.3, divisor 8). All values on the integer milli-ppm lattice.""" nd = iso_next(dstr) nval, law, rh, ex = 0, None, None, None for h0 in range(7, 24): vals = [hours[k] for k in (((dstr, h) if h < 24 else (nd, h - 24)) for h in range(h0, h0 + 8)) if k in hours] n = len(vals) if n == 0: continue s = sum(vals) if n >= 6: pass elif s >= 568: n = 8 else: continue nval += 1 wl = s // n # law: truncate the window mean (U 1.4) wr = (2 * s + n) // (2 * n) # variant: round half up we = Fraction(s, n) # variant: exact, carried through if law is None or wl > law: law = wl if rh is None or wr > rh: rh = wr if ex is None or we > ex: ex = we return nval, law, rh, ex def year_m4s(hourmap, year): """Annual fourth-highest daily maxima for the three value chains. Daily validity: >=13/17 valid windows, waived when the day's own maximum exceeds the standard, evaluated in each chain's own value space (U 1.6b).""" ystr = str(year) dates = {d for (d, h) in hourmap if d[:4] == ystr} dates |= {iso_prev(d) for (d, h) in hourmap if h <= 6 and iso_prev(d)[:4] == ystr} law_v, rh_v, ex_v = [], [], [] for dstr in sorted(dates): nval, law, rh, ex = day_chains(hourmap, dstr) if law is None: continue if nval >= 13 or law >= 71: law_v.append(law) if nval >= 13 or rh >= 71: rh_v.append(rh) if nval >= 13 or ex > 70: ex_v.append(ex) for v in (law_v, rh_v, ex_v): v.sort(reverse=True) return (law_v[3] if len(law_v) >= 4 else None, rh_v[3] if len(rh_v) >= 4 else None, ex_v[3] if len(ex_v) >= 4 else None) def record_chains(hourmap, y0): """Three-year record -> per-chain (m4 list, exact pre-final mean, DV, verdict).""" m4l, m4r, m4e = [], [], [] for y in (y0, y0 + 1, y0 + 2): a, b, c = year_m4s(hourmap, y) check(a is not None and b is not None and c is not None, f"year {y}: fewer than four valid days in some chain") m4l.append(a) m4r.append(b) m4e.append(c) pre_a = Fraction(sum(m4l), 3) pre_b = Fraction(sum(m4r), 3) pre_e = sum((Fraction(x) for x in m4e), Fraction(0)) / 3 return { "law": {"m4": m4l, "pre": pre_a, "dv": tz(pre_a)}, "all-round": {"m4": m4r, "pre": pre_b, "dv": rhu(pre_b)}, "2010": {"m4": m4e, "pre": pre_e, "dv": rhu(pre_e)}, "exact-trunc": {"m4": m4e, "pre": pre_e, "dv": tz(pre_e)}, } # ------------------------------------------------------------------ slice loading def load_slice(data=None): """-> {site: {(date, hour): milli}}. Verifies the slice sha256 first unless raw bytes are handed in (the byte-flip control does its own comparison).""" if data is None: with open(os.path.join(HERE, SLICE_NAME), "rb") as f: data = f.read() check(sha256_bytes(data) == SLICE_SHA256, f"{SLICE_NAME}: sha256 mismatch (corrupted or altered file)") text = gzip.decompress(data).decode("utf-8") rows = [ln for ln in text.splitlines() if ln and not ln.startswith("#")] rdr = csv.reader(io.StringIO("\n".join(rows))) hdr = next(rdr) check(hdr == ["site", "poc", "date", "hour", "o3_ppm"], f"unexpected slice columns {hdr}") sites = {} npoc = {} for site, poc, dstr, hh, val in rdr: sites.setdefault(site, {})[(dstr, int(hh))] = milli_from_str(val) npoc.setdefault(site, set()).add(poc) for site, ps in npoc.items(): check(len(ps) == 1, f"site {site}: more than one monitor in the slice") return sites # ------------------------------------------------------------------ default mode RECORDS = [ # (site, first year, label) ("170314007", 2020, "Cook County IL (Chicago area)"), ("482570005", 2022, "Kaufman County TX"), ("090090027", 2020, "New Haven County CT"), ("090090027", 2021, "New Haven County CT"), ("090090027", 2023, "New Haven County CT"), ] # (record index, chain) -> (m4 list or None, pre-final milli, DV milli, verdict) EXPECT = { (0, "law"): ([72, 69, 70], Fraction(211, 3), 70, "ATTAIN"), (0, "all-round"): ([73, 70, 71], Fraction(214, 3), 71, "VIOLATE"), (0, "2010"): (None, Fraction(853, 12), 71, "VIOLATE"), (0, "exact-trunc"): (None, Fraction(853, 12), 71, "VIOLATE"), (1, "law"): ([69, 71, 72], Fraction(212, 3), 70, "ATTAIN"), (1, "all-round"): ([70, 72, 73], Fraction(215, 3), 72, "VIOLATE"), (1, "2010"): (None, Fraction(857, 12), 71, "VIOLATE"), (1, "exact-trunc"): (None, Fraction(857, 12), 71, "VIOLATE"), (2, "law"): ([68, 71, 72], Fraction(211, 3), 70, "ATTAIN"), (3, "law"): ([71, 72, 69], Fraction(212, 3), 70, "ATTAIN"), (4, "law"): ([69, 77, 65], Fraction(211, 3), 70, "ATTAIN"), } def run_default(sites, verbose=True): """Compute the five records' chains and enforce every page value. Returns {record index: chains dict} for reuse by the controls.""" out = {} for i, (site, y0, label) in enumerate(RECORDS): ch = record_chains(sites[site], y0) out[i] = ch if verbose: print(f"\n{site} @ {y0}-{y0+2} ({label})", flush=True) for cn in ("law", "all-round", "2010", "exact-trunc"): c = ch[cn] m4s = ", ".join(m4_str(x) for x in c["m4"]) print(f" {cn:11s} 4th-highest [{m4s}] ppm ; " f"3-yr mean {milli_to_ppm_str(c['pre'])} ppm ; " f"design value 0.{c['dv']:03d} ppm -> {det(c['dv'])}", flush=True) for (i, cn), (m4, pre, dv, verdict) in sorted(EXPECT.items()): site, y0, _ = RECORDS[i] rid = f"{site}@{y0}-{y0+2}" c = out[i][cn] if m4 is not None: check(c["m4"] == m4, f"{rid} {cn}: 4th-highest {c['m4']} != {m4}") check(c["pre"] == pre, f"{rid} {cn}: mean {c['pre']} != {pre} (milli-ppm)") check(c["dv"] == dv, f"{rid} {cn}: design value {c['dv']} != {dv}") check(det(c["dv"]) == verdict, f"{rid} {cn}: {det(c['dv'])} != {verdict}") # every attaining record's exact-chain value sits strictly below 71.875 ppb for i, ch in out.items(): if ch["law"]["dv"] <= 70: check(ch["2010"]["pre"] < Fraction(575, 8), f"record {i}: exact mean {ch['2010']['pre']} not below 575/8 milli") if verbose: nh = out[3]["law"]["pre"] print(f"\nNote: New Haven 2021-2023's exact law mean {milli_to_ppm_str(nh)} ppm " f"rounds half-up to 0.{rhu(nh):03d} ppm; the law truncates it to " f"0.{tz(nh):03d}.", flush=True) print("Every attaining record's exact-chain mean sits strictly below " "575/8000 ppm (71.875 ppb).", flush=True) print("\nAll page values check out (5 records x 4 chains).", flush=True) return out # ------------------------------------------------------------------ census mode def load_census(): out = [] for name, pin in ((TERMINAL_NAME, TERMINAL_SHA256), (WHOLECHAIN_NAME, WHOLECHAIN_SHA256)): with open(os.path.join(HERE, name), "rb") as f: b = f.read() check(sha256_bytes(b) == pin, f"{name}: sha256 mismatch (corrupted or altered file)") out.append(json.loads(b.decode("utf-8"))) return out # [terminal, wholechain] def census_counts(term, whole): """Recount the paper's headline numbers from exact values only. The digit rules are re-applied here to each record's exact mean; the shipped files carry no verdict labels to copy from.""" wrecs = whole["records"] ok = [r for r in wrecs if r.get("status") == "OK"] cls = [r for r in ok if r["workbook_valid"] and r["n_years"] == 3] def pre(r, chain): return fr(r["chains"][chain]["prefinal_milli"]) n = {} n["o3_headline"] = len(cls) n["terminal_only"] = sum(1 for r in cls if rhu(pre(r, "a_law")) >= 71) n["all_round"] = sum(1 for r in cls if rhu(pre(r, "b_allround")) >= 71) n["chain_2010"] = sum(1 for r in cls if rhu(pre(r, "c_2010")) >= 71) n["cell_shift"] = sum(1 for r in cls if tz(pre(r, "b_allround")) != tz(pre(r, "a_law"))) n["o3_at_standard"] = sum(1 for r in cls if tz(pre(r, "a_law")) == 70) attaining = [pre(r, "c_2010") for r in cls if tz(pre(r, "a_law")) <= 70] n["max_attaining_exact"] = max(attaining) n["all_attaining_below"] = all(v < Fraction(575, 8) for v in attaining) ann = term["records"]["pm25_annual"] annv = [r for r in ann if r["workbook_valid"] and r.get("n_years") == 3 and "exact_mean_ugm3" in r] n["pm_annual_valid"] = len(annv) n["pm_annual_flips"] = sum( 1 for r in annv if tz(10 * fr(r["exact_mean_ugm3"])) <= 90 < rhu(10 * fr(r["exact_mean_ugm3"]))) n["pm_annual_at_standard"] = sum( 1 for r in annv if rhu(10 * fr(r["exact_mean_ugm3"])) == 90) day = term["records"]["pm25_24hr"] dayv = [r for r in day if r["workbook_valid"] and r.get("n_years") == 3 and "exact_mean_ugm3" in r] n["pm_daily_valid"] = len(dayv) n["pm_daily_flips"] = sum( 1 for r in dayv if tz(fr(r["exact_mean_ugm3"])) <= 35 < rhu(fr(r["exact_mean_ugm3"]))) n["pm_daily_at_standard"] = sum( 1 for r in dayv if rhu(fr(r["exact_mean_ugm3"])) == 35) n["at_standard_total"] = (n["o3_at_standard"] + n["pm_annual_at_standard"] + n["pm_daily_at_standard"]) classes = ("o3", "pm25_annual", "pm25_24hr") n["match_exact"] = sum(1 for cl in classes for r in term["records"][cl] if r["match"] == "EXACT-MATCH") n["match_total"] = sum(len(term["records"][cl]) for cl in classes) return n, ok, cls CENSUS_EXPECT = { "o3_headline": 190, "terminal_only": 61, "all_round": 148, "chain_2010": 140, "cell_shift": 82, "o3_at_standard": 190, "pm_annual_valid": 75, "pm_annual_flips": 33, "pm_annual_at_standard": 42, "pm_daily_valid": 39, "pm_daily_flips": 14, "pm_daily_at_standard": 25, "at_standard_total": 257, "match_exact": 330, "match_total": 340, "max_attaining_exact": Fraction(857, 12), "all_attaining_below": True, } def run_census(verbose=True): term, whole = load_census() n, ok, cls = census_counts(term, whole) if verbose: print(f"\nheadline ozone class (recomputable, workbook-valid, 3 years): " f"{n['o3_headline']}", flush=True) print(f" terminal-digit-only flips (round-half-up of the exact law mean " f">= 0.071): {n['terminal_only']}/{n['o3_headline']}", flush=True) print(f" all-rounding chain flips: {n['all_round']}/{n['o3_headline']}", flush=True) print(f" 2010 single-cut chain flips: {n['chain_2010']}/{n['o3_headline']}", flush=True) print(f" integer-cell shifts (all-round vs law mean): " f"{n['cell_shift']}/{n['o3_headline']}", flush=True) print(f"PM2.5 annual: {n['pm_annual_flips']} of {n['pm_annual_valid']} valid " f"records flip on truncate-vs-round at 9.0", flush=True) print(f"PM2.5 24-hour: {n['pm_daily_flips']} of {n['pm_daily_valid']} valid " f"records flip on truncate-vs-round at 35", flush=True) print(f"records with design value exactly at a standard: " f"{n['at_standard_total']} = {n['o3_at_standard']} ozone + " f"{n['pm_annual_at_standard']} PM annual + {n['pm_daily_at_standard']} " f"PM 24-hour", flush=True) print(f"recomputation vs published record: {n['match_exact']} of " f"{n['match_total']} EXACT-MATCH", flush=True) print(f"largest attaining exact ozone mean: " f"{milli_to_ppm_str(n['max_attaining_exact'])} ppm " f"(all attaining means < 575/8000 ppm)", flush=True) for k, want in CENSUS_EXPECT.items(): check(n[k] == want, f"census recount {k}: {n[k]!r} != {want!r}") # cross-checks between the two shipped files (values only, no verdict labels) tmap = {(r["site"], r["period"]): r for r in term["records"]["o3"]} for r in ok: t = tmap[(r["site"], r["vintage"])] check(fr(t["exact_mean_ppm"]) * 1000 == fr(r["chains"]["a_law"]["prefinal_milli"]), f"{r['id']}: terminal exact mean disagrees with whole-chain law mean") check(t["recomputed_dv_milli"] == tz(fr(t["exact_mean_ppm"]) * 1000), f"{r['id']}: stored law DV is not the truncation of the exact mean") for r in cls: for cn, final in (("a_law", tz), ("b_allround", rhu), ("c_2010", rhu), ("d_exacttrunc", tz)): c = r["chains"][cn] check(c["dv_milli"] == final(fr(c["prefinal_milli"])), f"{r['id']} {cn}: stored DV is not the terminal rule of the mean") if verbose: print("\nAll census recounts and cross-checks pass.", flush=True) return n # ------------------------------------------------------------------ controls def find_bump_key(hourmap, year): """Locate one filed hour inside the year's fourth-highest day of the exact chain (the day whose maximum enters the exact three-year mean, 583/8000 ppm for Cook County 2020): the first present hour of the window achieving that day's maximum, so raising it by one filing step must raise the mean.""" ystr = str(year) dates = {d for (d, h) in hourmap if d[:4] == ystr} scored = [] for dstr in sorted(dates): nval, law, rh, ex = day_chains(hourmap, dstr) if law is None: continue if nval >= 13 or ex > 70: scored.append((ex, dstr)) scored.sort(key=lambda t: (-t[0], t[1])) ex4, d4 = scored[3] nd = iso_next(d4) for h0 in range(7, 24): keys = [((d4, h) if h < 24 else (nd, h - 24)) for h in range(h0, h0 + 8)] vals = [hourmap[k] for k in keys if k in hourmap] n = len(vals) if n == 0: continue s = sum(vals) if n >= 6: pass elif s >= 568: n = 8 else: continue if Fraction(s, n) == ex4: for k in keys: if k in hourmap: return d4, k raise SystemExit("control setup: no window found for the fourth-highest day") def run_checks(sites): print("\n--- mutation controls (each must FAIL) ---", flush=True) base = run_default(sites, verbose=False) n_controls = 0 # (1) +0.001 ppm to one filed hour inside Cook County's 2020 fourth-highest day site = "170314007" d4, key = find_bump_key(sites[site], 2020) mut = dict(sites[site]) mut[key] = mut[key] + 1 check(mut != sites[site], "control 1: mutant equals control (dead knob)") ch = record_chains(mut, 2020) moved = [cn for cn in ch if ch[cn]["pre"] != base[0][cn]["pre"]] check(moved, "control 1 PASSED (no chain mean moved) - must fail") print(f"control 1 failed as required: +0.001 ppm at {key[0]} hour {key[1]:02d} " f"(inside the 4th-highest day {d4}) moved chains {moved}", flush=True) n_controls += 1 # (2) flip one byte of the gzipped slice -> sha256 comparison must fail with open(os.path.join(HERE, SLICE_NAME), "rb") as f: data = f.read() flipped = bytearray(data) flipped[len(flipped) // 2] ^= 0xFF flipped = bytes(flipped) check(flipped != data, "control 2: mutant equals control (dead knob)") check(sha256_bytes(flipped) != SLICE_SHA256, "control 2 PASSED (sha256 unchanged after byte flip) - must fail") print("control 2 failed as required: one flipped byte breaks the slice's " "sha256 comparison", flush=True) n_controls += 1 # (3) window-step truncation replaced by rounding inside the law chain variant = record_chains(sites["170314007"], 2020)["all-round"] law = base[0]["law"] check(variant["m4"] != law["m4"], "control 3: mutant equals control (dead knob)") dv_variant = tz(variant["pre"]) # law's terminal truncation kept exp_pre, exp_dv, exp_det = Fraction(211, 3), 70, "ATTAIN" check(not (variant["pre"] == exp_pre and dv_variant == exp_dv and det(dv_variant) == exp_det), "control 3 PASSED (law-chain check survived rounded windows) - must fail") print(f"control 3 failed as required: rounded window means give mean " f"{milli_to_ppm_str(variant['pre'])} ppm, design value 0.{dv_variant:03d} " f"{det(dv_variant)} (page value 211/3000 ppm, 0.070 ATTAIN)", flush=True) n_controls += 1 # (4) census mutation: one exact law mean 211/3 -> 214/3 milli term, whole = load_census() n0, _, _ = census_counts(term, whole) target = None for r in whole["records"]: if (r.get("status") == "OK" and r["workbook_valid"] and r["n_years"] == 3 and r["chains"]["a_law"]["prefinal_milli"] == "211/3"): target = r break check(target is not None, "control 4: no census record with law mean 211/3") before = target["chains"]["a_law"]["prefinal_milli"] target["chains"]["a_law"]["prefinal_milli"] = "214/3" check(target["chains"]["a_law"]["prefinal_milli"] != before, "control 4: mutant equals control (dead knob)") n1, _, _ = census_counts(term, whole) pre_b = fr(target["chains"]["b_allround"]["prefinal_milli"]) cell_delta = int(tz(pre_b) != tz(Fraction(214, 3))) - int(tz(pre_b) != tz(Fraction(211, 3))) check(n1["terminal_only"] == n0["terminal_only"] + 1, f"control 4: terminal-only count {n1['terminal_only']} != " f"{n0['terminal_only'] + 1}") check(n1["o3_at_standard"] == n0["o3_at_standard"] - 1, f"control 4: at-standard ozone count {n1['o3_at_standard']} != " f"{n0['o3_at_standard'] - 1}") check(n1["all_round"] == n0["all_round"], f"control 4: all-round count moved ({n1['all_round']} != {n0['all_round']}) " f"though the all-round mean was untouched") check(n1["cell_shift"] == n0["cell_shift"] + cell_delta, f"control 4: cell-shift count {n1['cell_shift']} != " f"{n0['cell_shift'] + cell_delta}") print(f"control 4 failed as required: {target['id']} law mean 211/3 -> 214/3 " f"milli moves the recounts (terminal-only {n0['terminal_only']} -> " f"{n1['terminal_only']}, at-standard {n0['o3_at_standard']} -> " f"{n1['o3_at_standard']}, cell shifts {n0['cell_shift']} -> " f"{n1['cell_shift']}; all-round unchanged at {n1['all_round']})", flush=True) n_controls += 1 print(f"\nAll {n_controls} mutation controls failed as required.", flush=True) # ------------------------------------------------------------------ main def main(): ap = argparse.ArgumentParser( description="Recompute the ozone design-value chain for the three monitors " "of the air-quality page, and recount the paper's census.") ap.add_argument("--census", action="store_true", help="recount the paper's headline census numbers from the " "exact per-record values in the two shipped JSON files") ap.add_argument("--check", action="store_true", help="run the default computation plus mutation controls " "that must fail") args = ap.parse_args() mode = "--check" if args.check else ("--census" if args.census else "default") print(f"airquality-evaluate.py ({mode}): rebuilding 40 CFR Part 50 Appendix U " f"arithmetic, exact fractions only", flush=True) if args.census: run_census() else: sites = load_slice() print(f"slice OK: sha256 verified, " f"{sum(len(v) for v in sites.values())} hourly readings, " f"{len(sites)} monitors", flush=True) run_default(sites) if args.check: run_checks(sites) run_census() print("\nPASS", flush=True) if __name__ == "__main__": main()