219 lines
10 KiB
Python
219 lines
10 KiB
Python
|
#!/usr/bin/env python
|
||
|
|
"""Profile the archive's real availability process into a sampler-ready model.
|
||
|
|
|
||
|
|
profile_archive.py --out /media/shared/NOAA_GOES/corpus_v8
|
||
|
|
|
||
|
|
Classifies every slot of the requested era, per satellite and band, as GOOD (a
|
||
|
|
non-stub frame judged good), BAD (judged bad) or MISSING (no usable frame), and
|
||
|
|
writes:
|
||
|
|
|
||
|
|
* ``availability.npz`` -- the state sequences themselves (:mod:`suvi.availability`).
|
||
|
|
This is the statistical model: an empirical joint process whose sampled windows
|
||
|
|
reproduce real run lengths, cross-satellite outage correlation, and the diurnal /
|
||
|
|
eclipse-season structure with no parametric assumptions. The corpus builder
|
||
|
|
selects pristine segments from it, and training draws simulated outages from it.
|
||
|
|
* ``availability_report.json`` -- the human-readable statistics: state fractions,
|
||
|
|
per-band state fractions and unusable-run quantiles, how many bands slots lose,
|
||
|
|
the commonest damage patterns, per-band dual-outage rates and runs, outage-start
|
||
|
|
time-of-day and month histograms, and per-band bad rates.
|
||
|
|
|
||
|
|
"Judged" means the stored production judgment (``filter_FITS.py judge``: header,
|
||
|
|
lunar, disc, field and clean-neighbour checks), the same verdicts every other
|
||
|
|
consumer reads -- not a header-only approximation, which passed 41% of the real
|
||
|
|
defects in a hand-confirmed audit. Refuses to profile an era with unjudged frames,
|
||
|
|
because "unjudged" silently counted as "good" would understate exactly the badness
|
||
|
|
this model exists to capture. Read-only against the index.
|
||
|
|
"""
|
||
|
|
|
||
|
|
import argparse
|
||
|
|
import collections
|
||
|
|
import datetime
|
||
|
|
import json
|
||
|
|
import os
|
||
|
|
import sys
|
||
|
|
|
||
|
|
import numpy as np
|
||
|
|
|
||
|
|
sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
|
||
|
|
|
||
|
|
from suvi import availability, dataset, db, detectors, paths
|
||
|
|
|
||
|
|
#: Default era. Starts at 2023: the 2022 directories were never fully fetched
|
||
|
|
#: (and are excluded from the puller), so their absences are download artifacts
|
||
|
|
#: that would masquerade as spacecraft outages in the model.
|
||
|
|
DEFAULT_FROM = "2023-01-01"
|
||
|
|
DEFAULT_TO = "2025-03-25"
|
||
|
|
|
||
|
|
|
||
|
|
def unix(day):
|
||
|
|
stamp = datetime.datetime.fromisoformat(day).replace(tzinfo=datetime.timezone.utc)
|
||
|
|
return int(stamp.timestamp())
|
||
|
|
|
||
|
|
|
||
|
|
def classify_era(conn, t_from, t_to, satellites, run_id):
|
||
|
|
"""Build the (satellites, bands, slots) state array from a stored judgment run.
|
||
|
|
|
||
|
|
A band is GOOD where a frame was judged good, BAD where the only frames were
|
||
|
|
judged bad, MISSING where there is no frame (or only a stub). Returns
|
||
|
|
(states, bad rate per band, number of frames the run has not judged).
|
||
|
|
"""
|
||
|
|
slots = (t_to - t_from) // paths.CADENCE
|
||
|
|
bands = availability.BANDS
|
||
|
|
states = np.full((len(satellites), len(bands), slots), availability.MISSING,
|
||
|
|
dtype=np.uint8)
|
||
|
|
row_of = {satellite: index for index, satellite in enumerate(satellites)}
|
||
|
|
column_of = {wavelength: index for index, wavelength in enumerate(bands)}
|
||
|
|
band_bad = collections.Counter()
|
||
|
|
band_total = collections.Counter()
|
||
|
|
unjudged = 0
|
||
|
|
|
||
|
|
cursor = conn.execute(
|
||
|
|
"""
|
||
|
|
SELECT f.satellite, f.wavelength, f.t_start, f.size_bytes, d.verdict
|
||
|
|
FROM frame f LEFT JOIN detection d ON d.frame_id = f.id AND d.run_id = ?
|
||
|
|
WHERE f.t_start >= ? AND f.t_start < ?
|
||
|
|
""",
|
||
|
|
(run_id, t_from, t_to),
|
||
|
|
)
|
||
|
|
for row in cursor:
|
||
|
|
if row["satellite"] not in row_of or row["wavelength"] not in column_of:
|
||
|
|
continue
|
||
|
|
index = (row["t_start"] // paths.CADENCE * paths.CADENCE - t_from) \
|
||
|
|
// paths.CADENCE
|
||
|
|
if not 0 <= index < slots:
|
||
|
|
continue
|
||
|
|
if row["size_bytes"] is not None and row["size_bytes"] < dataset.MIN_FRAME_BYTES:
|
||
|
|
continue # a stub is not a frame
|
||
|
|
if row["verdict"] is None:
|
||
|
|
unjudged += 1
|
||
|
|
continue
|
||
|
|
band_total[row["wavelength"]] += 1
|
||
|
|
state = availability.GOOD
|
||
|
|
if row["verdict"] == detectors.BAD:
|
||
|
|
state = availability.BAD
|
||
|
|
band_bad[row["wavelength"]] += 1
|
||
|
|
cell = (row_of[row["satellite"]], column_of[row["wavelength"]], index)
|
||
|
|
# GOOD < BAD < MISSING: a good frame in the slot wins over a bad one.
|
||
|
|
states[cell] = min(states[cell], state)
|
||
|
|
|
||
|
|
rates = {
|
||
|
|
int(wavelength): band_bad[wavelength] / band_total[wavelength]
|
||
|
|
for wavelength in sorted(band_total)
|
||
|
|
}
|
||
|
|
return states, rates, unjudged
|
||
|
|
|
||
|
|
|
||
|
|
def quantiles(lengths):
|
||
|
|
if lengths.size == 0:
|
||
|
|
return {}
|
||
|
|
marks = (0.5, 0.9, 0.99, 1.0)
|
||
|
|
values = np.quantile(lengths, marks)
|
||
|
|
return {f"p{int(mark * 100)}": float(value)
|
||
|
|
for mark, value in zip(marks, values)} | {"n": int(lengths.size)}
|
||
|
|
|
||
|
|
|
||
|
|
def summarise(model, t_from):
|
||
|
|
"""The report: every statistic a human (or a test) needs to trust the model."""
|
||
|
|
report = {"fractions": {str(s): f for s, f in model.fractions().items()}}
|
||
|
|
names = [str(b) for b in availability.BANDS]
|
||
|
|
|
||
|
|
for satellite in model.satellites:
|
||
|
|
row = model.states[model.satellites.index(satellite)] # (bands, slots)
|
||
|
|
report[f"bands_g{satellite}"] = {
|
||
|
|
name: {state: float((row[i] == value).mean())
|
||
|
|
for value, state in enumerate(availability.STATE_NAMES)}
|
||
|
|
for i, name in enumerate(names)}
|
||
|
|
report[f"unusable_runs_g{satellite}"] = {
|
||
|
|
name: quantiles(availability.run_lengths((row[i] != availability.GOOD)
|
||
|
|
.astype(np.uint8), 1))
|
||
|
|
for i, name in enumerate(names)}
|
||
|
|
lost = (row != availability.GOOD).sum(axis=0)
|
||
|
|
report[f"bands_lost_per_slot_g{satellite}"] = np.bincount(
|
||
|
|
lost, minlength=len(names) + 1).tolist()
|
||
|
|
|
||
|
|
masks, probabilities = model.damage_patterns()
|
||
|
|
order = np.argsort(-probabilities)[:12]
|
||
|
|
report["damage_patterns"] = {
|
||
|
|
"+".join(n for n, lost in zip(names, masks[k]) if lost): float(probabilities[k])
|
||
|
|
for k in order}
|
||
|
|
|
||
|
|
# Per band, both satellites without it -- what the fill policy is applied to.
|
||
|
|
dual = (model.states != availability.GOOD).all(axis=0) # (bands, slots)
|
||
|
|
report["dual_unusable_fraction"] = {name: float(dual[i].mean())
|
||
|
|
for i, name in enumerate(names)}
|
||
|
|
report["dual_unusable_runs"] = quantiles(
|
||
|
|
availability.run_lengths(dual.any(axis=0).astype(np.uint8), 1))
|
||
|
|
|
||
|
|
# When do whole-slot outages *start*? Eclipse seasons put daily structure here.
|
||
|
|
dark = (model.states[0] == availability.MISSING).all(axis=0)
|
||
|
|
starts = np.flatnonzero(np.diff(dark.astype(np.int8)) == 1)
|
||
|
|
when = t_from + (starts + 1) * paths.CADENCE
|
||
|
|
hours = ((when % 86400) // 3600).astype(np.int64)
|
||
|
|
months = np.array([datetime.datetime.fromtimestamp(
|
||
|
|
int(t), datetime.timezone.utc).month for t in when], dtype=np.int64)
|
||
|
|
report["outage_start_by_utc_hour"] = np.bincount(hours, minlength=24).tolist()
|
||
|
|
report["outage_start_by_month"] = np.bincount(months, minlength=13).tolist()[1:]
|
||
|
|
return report
|
||
|
|
|
||
|
|
|
||
|
|
def main(argv=None):
|
||
|
|
parser = argparse.ArgumentParser(description=__doc__,
|
||
|
|
formatter_class=argparse.RawDescriptionHelpFormatter)
|
||
|
|
parser.add_argument("--out", required=True,
|
||
|
|
help="directory for availability.npz + report")
|
||
|
|
parser.add_argument("--from", dest="from_day", default=DEFAULT_FROM)
|
||
|
|
parser.add_argument("--to", dest="to_day", default=DEFAULT_TO)
|
||
|
|
parser.add_argument("--satellites", default="16,18")
|
||
|
|
parser.add_argument("--run", type=int, default=None,
|
||
|
|
help="judgment run id (default: the run made under the current rules)")
|
||
|
|
parser.add_argument("--db", default=None)
|
||
|
|
parser.add_argument("--allow-partial", action="store_true",
|
||
|
|
help="profile despite unjudged frames (development only)")
|
||
|
|
args = parser.parse_args(argv)
|
||
|
|
|
||
|
|
satellites = tuple(int(s) for s in args.satellites.split(","))
|
||
|
|
t_from, t_to = unix(args.from_day), unix(args.to_day)
|
||
|
|
conn = db.connect(args.db, readonly=True)
|
||
|
|
run_id = args.run or detectors.judgment_run(conn)
|
||
|
|
if run_id is None:
|
||
|
|
raise SystemExit("No judgment run under the current rules; run "
|
||
|
|
"`filter_FITS.py judge` over the era first.")
|
||
|
|
|
||
|
|
states, band_rates, unjudged = classify_era(conn, t_from, t_to, satellites, run_id)
|
||
|
|
if unjudged:
|
||
|
|
message = (f"{unjudged:,} frames in the era have no verdict in judgment run "
|
||
|
|
f"{run_id} (run filter_FITS.py judge to completion)")
|
||
|
|
if not args.allow_partial:
|
||
|
|
raise SystemExit(f"REFUSING to profile: {message}. An unjudged frame "
|
||
|
|
"silently counted as good would understate the badness "
|
||
|
|
"this model exists to capture.")
|
||
|
|
print(f"WARNING: {message}; proceeding under --allow-partial")
|
||
|
|
|
||
|
|
model = availability.Availability(states=states, t0=t_from,
|
||
|
|
cadence=paths.CADENCE, satellites=satellites)
|
||
|
|
report = summarise(model, t_from)
|
||
|
|
report["era"] = {"from": args.from_day, "to": args.to_day,
|
||
|
|
"slots": int(states.shape[2]), "unjudged_frames": int(unjudged)}
|
||
|
|
report["bad_rate_by_band"] = band_rates
|
||
|
|
report["judgment_run"] = int(run_id)
|
||
|
|
|
||
|
|
os.makedirs(args.out, exist_ok=True)
|
||
|
|
model.save(os.path.join(args.out, "availability.npz"))
|
||
|
|
with open(os.path.join(args.out, "availability_report.json"), "w") as handle:
|
||
|
|
json.dump(report, handle, indent=2)
|
||
|
|
|
||
|
|
print(f"era {args.from_day}..{args.to_day}: {states.shape[2]:,} slots "
|
||
|
|
f"(judgment run {run_id})")
|
||
|
|
for satellite, fractions in model.fractions().items():
|
||
|
|
print(f" G{satellite}: " + " ".join(f"{k} {v:.2%}"
|
||
|
|
for k, v in fractions.items()))
|
||
|
|
worst = max(report["dual_unusable_fraction"].values())
|
||
|
|
print(f" a band unusable on both at once: up to {worst:.2%} of slots; "
|
||
|
|
f"longest run {report['dual_unusable_runs'].get('p100', 0):.0f} slots")
|
||
|
|
print(f" wrote availability.npz + availability_report.json under {args.out}")
|
||
|
|
return 0
|
||
|
|
|
||
|
|
|
||
|
|
if __name__ == "__main__":
|
||
|
|
sys.exit(main())
|