464 lines
16 KiB
Go
464 lines
16 KiB
Go
// Command compare-tawhiri runs identical predictions against a local predictor
|
|
// and a hosted Tawhiri instance and reports how closely they agree.
|
|
//
|
|
// To make the comparison test the engine rather than data drift, it discovers
|
|
// the local predictor's loaded GFS run via /ready and asks Tawhiri to use the
|
|
// same run (the `dataset` parameter), so both integrate identical wind data.
|
|
// It compares the burst apex (terrain-independent) and the landing point
|
|
// (terrain-dependent) separately, since without the ruaumoko elevation dataset
|
|
// the local predictor terminates descent at sea level while Tawhiri uses
|
|
// ground elevation.
|
|
//
|
|
// Usage:
|
|
//
|
|
// compare-tawhiri --server http://localhost:8080 # built-in suite
|
|
// compare-tawhiri --lat 52.2 --lng 0.1 --burst 30000 # single site
|
|
package main
|
|
|
|
import (
|
|
"encoding/json"
|
|
"flag"
|
|
"fmt"
|
|
"io"
|
|
"math"
|
|
"net/http"
|
|
"net/url"
|
|
"os"
|
|
"strings"
|
|
"text/tabwriter"
|
|
"time"
|
|
)
|
|
|
|
func main() {
|
|
var (
|
|
server = flag.String("server", "http://localhost:8080", "local predictor base URL")
|
|
tawhiri = flag.String("tawhiri", "https://api.v2.sondehub.org/tawhiri", "hosted Tawhiri base URL")
|
|
lat = flag.Float64("lat", math.NaN(), "launch latitude (single-site mode)")
|
|
lng = flag.Float64("lng", math.NaN(), "launch longitude (single-site mode)")
|
|
alt = flag.Float64("alt", 0, "launch altitude m")
|
|
ascent = flag.Float64("ascent-rate", 5, "ascent rate m/s")
|
|
burst = flag.Float64("burst", 30000, "burst altitude m")
|
|
descent = flag.Float64("descent-rate", 5, "descent rate m/s")
|
|
launch = flag.String("launch", "", "launch time RFC3339 (default: epoch + 3h)")
|
|
align = flag.Bool("align-dataset", true, "ask Tawhiri to use the local predictor's GFS run")
|
|
)
|
|
flag.Parse()
|
|
|
|
epoch, err := fetchActiveEpoch(*server)
|
|
if err != nil {
|
|
fmt.Fprintln(os.Stderr, "local /ready:", err)
|
|
os.Exit(1)
|
|
}
|
|
fmt.Printf("local dataset epoch: %s\n", epoch.Format(time.RFC3339))
|
|
|
|
launchTime := epoch.Add(3 * time.Hour)
|
|
if *launch != "" {
|
|
launchTime, err = time.Parse(time.RFC3339, *launch)
|
|
if err != nil {
|
|
fmt.Fprintln(os.Stderr, "invalid --launch:", err)
|
|
os.Exit(1)
|
|
}
|
|
}
|
|
datasetParam := ""
|
|
if *align {
|
|
datasetParam = epoch.Format(time.RFC3339)
|
|
}
|
|
|
|
sites := suite()
|
|
if !math.IsNaN(*lat) && !math.IsNaN(*lng) {
|
|
sites = []site{{name: "custom", lat: *lat, lng: *lng}}
|
|
}
|
|
|
|
tw := tabwriter.NewWriter(os.Stdout, 0, 0, 2, ' ', 0)
|
|
fmt.Fprintln(tw, "\nsite\tburst Δ\tlanding Δ\tapex alt Δ\tland alt Δ\tasc pts\tdesc pts\tnotes")
|
|
fmt.Fprintln(tw, "----\t-------\t---------\t----------\t----------\t-------\t--------\t-----")
|
|
|
|
var worst float64
|
|
compared := 0
|
|
for _, s := range sites {
|
|
// The site suite compares trajectories, not notations, so fold western
|
|
// longitudes into [0, 360) here: that is the only range Tawhiri accepts.
|
|
// The longitude suite below sends raw values precisely to test that edge.
|
|
lng := s.lng
|
|
if lng < 0 {
|
|
lng += 360
|
|
}
|
|
p := params{lat: s.lat, lng: lng, alt: *alt, launch: launchTime,
|
|
ascent: *ascent, burst: *burst, descent: *descent}
|
|
|
|
ours, _, err := predict(*server+"/api/v1/prediction", p, "")
|
|
if err != nil {
|
|
fmt.Fprintf(tw, "%s\tlocal error: %v\n", s.name, err)
|
|
continue
|
|
}
|
|
theirs, _, err := predict(*tawhiri, p, datasetParam)
|
|
if err != nil {
|
|
fmt.Fprintf(tw, "%s\ttawhiri error: %v\n", s.name, err)
|
|
continue
|
|
}
|
|
compared++
|
|
|
|
burstD := haversine(ours.apexLat, ours.apexLng, theirs.apexLat, theirs.apexLng)
|
|
landD := haversine(ours.landLat, ours.landLng, theirs.landLat, theirs.landLng)
|
|
if landD > worst {
|
|
worst = landD
|
|
}
|
|
note := ""
|
|
if theirs.dataset != "" && ours.dataset != "" && theirs.dataset != ours.dataset {
|
|
note = fmt.Sprintf("dataset mismatch (theirs=%s)", theirs.dataset)
|
|
}
|
|
fmt.Fprintf(tw, "%s\t%.0f m\t%.2f km\t%.0f m\t%.0f m\t%d/%d\t%d/%d\t%s\n",
|
|
s.name, burstD, landD/1000,
|
|
math.Abs(ours.apexAlt-theirs.apexAlt), math.Abs(ours.landAlt-theirs.landAlt),
|
|
ours.ascPts, theirs.ascPts, ours.descPts, theirs.descPts, note)
|
|
}
|
|
tw.Flush()
|
|
|
|
breaks := comparePairs(*server, *tawhiri, longitudeParams(launchTime, *alt, *ascent, *burst, *descent), datasetParam)
|
|
breaks += compareLongitudes(*server, *tawhiri, longitudeParams(launchTime, *alt, *ascent, *burst, *descent), datasetParam)
|
|
|
|
if compared == 0 {
|
|
fmt.Println("\nVERDICT: NO COMPARISONS (every site errored — see rows above)")
|
|
os.Exit(1)
|
|
}
|
|
if breaks > 0 {
|
|
fmt.Printf("\nVERDICT: INCOMPATIBLE (%d longitude(s) Tawhiri accepts and we reject)\n", breaks)
|
|
os.Exit(3)
|
|
}
|
|
fmt.Printf("\ncompared %d/%d sites; worst landing distance: %.2f km\n", compared, len(sites), worst/1000)
|
|
switch {
|
|
case worst < 1000:
|
|
fmt.Println("VERDICT: MATCH (all landings < 1 km — engine agrees with Tawhiri)")
|
|
case worst < 50000:
|
|
fmt.Println("VERDICT: CLOSE (< 50 km — consistent with elevation/dataset differences)")
|
|
default:
|
|
fmt.Println("VERDICT: DIVERGENT (> 50 km — investigate)")
|
|
os.Exit(2)
|
|
}
|
|
}
|
|
|
|
type site struct {
|
|
name string
|
|
lat, lng float64
|
|
}
|
|
|
|
// suite is a small set of diverse launch points: UK (lands on land/sea
|
|
// depending on winds), mid-Atlantic and mid-Pacific (ocean landings, so the
|
|
// sea-level-vs-terrain difference vanishes), and southern hemisphere.
|
|
func suite() []site {
|
|
return []site{
|
|
{"cambridge-uk", 52.2135, 0.0964},
|
|
{"mid-atlantic", 35.0, -40.0},
|
|
{"mid-pacific", 0.0, -160.0},
|
|
{"new-zealand", -41.3, 174.8},
|
|
{"colorado-us", 39.0, -105.5},
|
|
}
|
|
}
|
|
|
|
type params struct {
|
|
lat, lng, alt float64
|
|
launch time.Time
|
|
ascent, burst, descent float64
|
|
}
|
|
|
|
type result struct {
|
|
apexLat, apexLng, apexAlt float64
|
|
landLat, landLng, landAlt float64
|
|
ascPts, descPts int
|
|
dataset string
|
|
}
|
|
|
|
// predict sends p verbatim and returns the parsed result plus the HTTP status.
|
|
//
|
|
// Nothing is normalised here on purpose. This function used to fold negative
|
|
// longitudes into [0, 360) before sending, which made every longitude notation
|
|
// look identical to both endpoints — so the tool could not see a disagreement
|
|
// about notation even in principle, and ours drifted to [-180, 360) with an
|
|
// asymmetry nobody had chosen. Callers that want a fold do it themselves.
|
|
//
|
|
// The status is returned separately from the error so a caller can tell a
|
|
// deliberate rejection from a transport failure.
|
|
func predict(endpoint string, p params, dataset string) (result, int, error) {
|
|
q := url.Values{}
|
|
q.Set("launch_latitude", fmt.Sprintf("%.4f", p.lat))
|
|
q.Set("launch_longitude", fmt.Sprintf("%.4f", p.lng))
|
|
q.Set("launch_altitude", fmt.Sprintf("%.0f", p.alt))
|
|
q.Set("launch_datetime", p.launch.Format(time.RFC3339))
|
|
q.Set("ascent_rate", fmt.Sprintf("%.2f", p.ascent))
|
|
q.Set("burst_altitude", fmt.Sprintf("%.0f", p.burst))
|
|
q.Set("descent_rate", fmt.Sprintf("%.2f", p.descent))
|
|
if dataset != "" {
|
|
q.Set("dataset", dataset)
|
|
}
|
|
|
|
full := endpoint + "?" + q.Encode()
|
|
var body []byte
|
|
var status int
|
|
var lastErr error
|
|
for range 3 {
|
|
resp, err := http.Get(full)
|
|
if err != nil {
|
|
lastErr = err
|
|
time.Sleep(time.Second)
|
|
continue
|
|
}
|
|
body, _ = io.ReadAll(resp.Body)
|
|
status = resp.StatusCode
|
|
resp.Body.Close()
|
|
lastErr = nil
|
|
break
|
|
}
|
|
if lastErr != nil {
|
|
return result{}, 0, lastErr
|
|
}
|
|
if status != 200 {
|
|
return result{}, status, fmt.Errorf("HTTP %d: %s", status, truncate(string(body), 160))
|
|
}
|
|
|
|
var doc struct {
|
|
Prediction []struct {
|
|
Stage string `json:"stage"`
|
|
Trajectory []struct {
|
|
Latitude float64 `json:"latitude"`
|
|
Longitude float64 `json:"longitude"`
|
|
Altitude float64 `json:"altitude"`
|
|
} `json:"trajectory"`
|
|
} `json:"prediction"`
|
|
Request struct {
|
|
Dataset string `json:"dataset"`
|
|
} `json:"request"`
|
|
}
|
|
if err := json.Unmarshal(body, &doc); err != nil {
|
|
return result{}, status, err
|
|
}
|
|
|
|
var r result
|
|
r.dataset = doc.Request.Dataset
|
|
for _, st := range doc.Prediction {
|
|
if len(st.Trajectory) == 0 {
|
|
continue
|
|
}
|
|
last := st.Trajectory[len(st.Trajectory)-1]
|
|
switch st.Stage {
|
|
case "ascent":
|
|
r.ascPts = len(st.Trajectory)
|
|
r.apexLat, r.apexLng, r.apexAlt = last.Latitude, last.Longitude, last.Altitude
|
|
case "descent":
|
|
r.descPts = len(st.Trajectory)
|
|
r.landLat, r.landLng, r.landAlt = last.Latitude, last.Longitude, last.Altitude
|
|
}
|
|
}
|
|
return r, status, nil
|
|
}
|
|
|
|
type readinessResp struct {
|
|
Status string `json:"status"`
|
|
DatasetTime string `json:"dataset_time"`
|
|
}
|
|
|
|
func fetchActiveEpoch(base string) (time.Time, error) {
|
|
resp, err := http.Get(base + "/ready")
|
|
if err != nil {
|
|
return time.Time{}, err
|
|
}
|
|
defer resp.Body.Close()
|
|
body, _ := io.ReadAll(resp.Body)
|
|
if resp.StatusCode != 200 {
|
|
return time.Time{}, fmt.Errorf("HTTP %d: %s", resp.StatusCode, string(body))
|
|
}
|
|
var r readinessResp
|
|
if err := json.Unmarshal(body, &r); err != nil {
|
|
return time.Time{}, err
|
|
}
|
|
if r.Status != "ok" {
|
|
return time.Time{}, fmt.Errorf("server status %q (no dataset loaded yet)", r.Status)
|
|
}
|
|
return time.Parse(time.RFC3339, r.DatasetTime)
|
|
}
|
|
|
|
func haversine(lat1, lng1, lat2, lng2 float64) float64 {
|
|
const R = 6371000.0
|
|
phi1 := lat1 * math.Pi / 180
|
|
phi2 := lat2 * math.Pi / 180
|
|
dphi := (lat2 - lat1) * math.Pi / 180
|
|
dlam := (lng2 - lng1) * math.Pi / 180
|
|
a := math.Sin(dphi/2)*math.Sin(dphi/2) + math.Cos(phi1)*math.Cos(phi2)*math.Sin(dlam/2)*math.Sin(dlam/2)
|
|
return R * 2 * math.Atan2(math.Sqrt(a), math.Sqrt(1-a))
|
|
}
|
|
|
|
func truncate(s string, n int) string {
|
|
if len(s) <= n {
|
|
return s
|
|
}
|
|
return s[:n] + "…"
|
|
}
|
|
|
|
// longitudeCases are the notations a client can put in launch_longitude, sent
|
|
// verbatim to both endpoints.
|
|
//
|
|
// Longitude is on a circle, so every one of these names a real meridian and no
|
|
// value here is geometrically wrong; what differs is which notations each service
|
|
// agrees to read. Tawhiri implements a half-open [0, 360) — one clean fundamental
|
|
// domain, every meridian spelled exactly once. We accept a superset deliberately,
|
|
// because the frontend holds longitudes in [-180, 180] and a western launch
|
|
// therefore arrives negative.
|
|
//
|
|
// The pairs below name the same meridian in the two notations, so an accepted
|
|
// pair must produce the same trajectory: (-51.7, 308.3), (-180, 180), (0, 360),
|
|
// (-90, 270). That is what makes this a comparison rather than a status table.
|
|
func longitudeCases() []float64 {
|
|
return []float64{
|
|
0, 180, 270, 308.3, 359.999, // inside Tawhiri's domain
|
|
-0.0001, -51.7, -90, -180, -200, // signed convention: our extension
|
|
360, // 0 spelled redundantly; Tawhiri refuses it
|
|
5170, // a plain typo, refused by both
|
|
}
|
|
}
|
|
|
|
func longitudeParams(launch time.Time, alt, ascent, burst, descent float64) params {
|
|
// 64.1 N is Nuuk's latitude: high enough to be a realistic Arctic site, far
|
|
// enough from the pole that the comparison is about longitude alone.
|
|
return params{lat: 64.1, alt: alt, launch: launch, ascent: ascent, burst: burst, descent: descent}
|
|
}
|
|
|
|
// compareLongitudes reports how the two services read each notation and returns
|
|
// the number of compatibility breaks — longitudes Tawhiri accepts and we refuse.
|
|
// Those are the only rows that are a defect: refusing what a Tawhiri client is
|
|
// entitled to send is what a drop-in replacement must never do. The reverse, us
|
|
// accepting what Tawhiri refuses, is the extension the product depends on.
|
|
func compareLongitudes(server, tawhiri string, base params, datasetParam string) int {
|
|
tw := tabwriter.NewWriter(os.Stdout, 0, 0, 2, ' ', 0)
|
|
fmt.Fprintln(tw, "\nlongitude\tours\ttawhiri\tlanding Δ\tverdict")
|
|
fmt.Fprintln(tw, "---------\t----\t-------\t---------\t-------")
|
|
|
|
breaks := 0
|
|
for _, lng := range longitudeCases() {
|
|
p := base
|
|
p.lng = lng
|
|
|
|
ours, oursStatus, oursErr := predict(server+"/api/v1/prediction", p, "")
|
|
theirs, theirsStatus, theirsErr := predict(tawhiri, p, datasetParam)
|
|
|
|
code := func(status int, err error) string {
|
|
if err == nil {
|
|
return "200"
|
|
}
|
|
if status == 0 {
|
|
return "transport"
|
|
}
|
|
return fmt.Sprintf("%d", status)
|
|
}
|
|
|
|
delta := ""
|
|
verdict := ""
|
|
switch {
|
|
case oursErr == nil && theirsErr == nil:
|
|
d := haversine(ours.landLat, ours.landLng, theirs.landLat, theirs.landLng)
|
|
delta = fmt.Sprintf("%.2f km", d/1000)
|
|
verdict = "agree"
|
|
if d > 50000 {
|
|
verdict = "ACCEPTED BY BOTH, TRAJECTORIES DIVERGE"
|
|
}
|
|
// v1 must publish the convention Tawhiri publishes, not merely the same
|
|
// place in a different notation. It used to answer -51.2115 where
|
|
// Tawhiri answers 308.7884; haversine cannot see that, because a 360
|
|
// degree difference vanishes from sin(dlam/2).
|
|
if (ours.landLng < 0) != (theirs.landLng < 0) {
|
|
verdict = fmt.Sprintf("OUTPUT FORMAT DIFFERS (ours %.4f, theirs %.4f)", ours.landLng, theirs.landLng)
|
|
breaks++
|
|
}
|
|
case !aboutLongitude(oursErr) || !aboutLongitude(theirsErr):
|
|
// A rejection for some other reason says nothing about notation. Before
|
|
// this check, an upstream "No matching dataset found" — which happens as
|
|
// soon as sondehub rotates past our GFS run — was reported as "Tawhiri
|
|
// refuses this longitude", i.e. the tool inventing compatibility news.
|
|
verdict = "inconclusive: " + firstReason(oursErr, theirsErr)
|
|
case oursErr != nil && theirsErr != nil:
|
|
verdict = "agree (both refuse)"
|
|
case oursErr != nil:
|
|
verdict = "BREAK: Tawhiri accepts, we refuse"
|
|
breaks++
|
|
default:
|
|
verdict = "extension: we accept, Tawhiri refuses"
|
|
}
|
|
fmt.Fprintf(tw, "%g\t%s\t%s\t%s\t%s\n",
|
|
lng, code(oursStatus, oursErr), code(theirsStatus, theirsErr), delta, verdict)
|
|
}
|
|
tw.Flush()
|
|
return breaks
|
|
}
|
|
|
|
// aboutLongitude reports whether err is a rejection of the longitude itself
|
|
// rather than of something else in the request. A nil error passes: a successful
|
|
// call is always informative.
|
|
func aboutLongitude(err error) bool {
|
|
if err == nil {
|
|
return true
|
|
}
|
|
return strings.Contains(strings.ToLower(err.Error()), "longitude")
|
|
}
|
|
|
|
// firstReason summarises whichever side failed for an unrelated reason.
|
|
func firstReason(errs ...error) string {
|
|
for _, err := range errs {
|
|
if err != nil && !aboutLongitude(err) {
|
|
return truncate(err.Error(), 90)
|
|
}
|
|
}
|
|
return "unknown"
|
|
}
|
|
|
|
// longitudePairs are two notations naming one meridian: signed, then unsigned.
|
|
//
|
|
// Tawhiri reads only the unsigned one, so this is the only way the signed form —
|
|
// the extension the product actually depends on — can be checked against the
|
|
// reference at all. Our answer for the signed value must match Tawhiri's answer
|
|
// for its unsigned twin. Without this the "extension" rows in the table above are
|
|
// merely unrefuted, not verified.
|
|
//
|
|
// (0, 360) is absent on purpose: Tawhiri refuses 360, so that pair has no
|
|
// reference side.
|
|
func longitudePairs() [][2]float64 {
|
|
return [][2]float64{{-51.7, 308.3}, {-180, 180}, {-90, 270}, {-0.0001, 359.9999}}
|
|
}
|
|
|
|
// comparePairs checks each signed notation we accept against Tawhiri's answer for
|
|
// the same meridian written the way Tawhiri accepts it. Returns the number of
|
|
// pairs that disagree by more than the site suite's own tolerance.
|
|
func comparePairs(server, tawhiri string, base params, datasetParam string) int {
|
|
tw := tabwriter.NewWriter(os.Stdout, 0, 0, 2, ' ', 0)
|
|
fmt.Fprintln(tw, "\nours (signed)\tvs tawhiri (unsigned)\tlanding Δ\tverdict")
|
|
fmt.Fprintln(tw, "-------------\t---------------------\t---------\t-------")
|
|
|
|
bad := 0
|
|
for _, pair := range longitudePairs() {
|
|
signed, unsigned := pair[0], pair[1]
|
|
|
|
p := base
|
|
p.lng = signed
|
|
ours, _, oursErr := predict(server+"/api/v1/prediction", p, "")
|
|
|
|
q := base
|
|
q.lng = unsigned
|
|
theirs, _, theirsErr := predict(tawhiri, q, datasetParam)
|
|
|
|
switch {
|
|
case oursErr != nil:
|
|
fmt.Fprintf(tw, "%g\t%g\t\tours refused: %v\n", signed, unsigned, oursErr)
|
|
bad++
|
|
case theirsErr != nil:
|
|
fmt.Fprintf(tw, "%g\t%g\t\ttawhiri refused the unsigned twin: %v\n", signed, unsigned, theirsErr)
|
|
bad++
|
|
default:
|
|
d := haversine(ours.landLat, ours.landLng, theirs.landLat, theirs.landLng)
|
|
verdict := "agree — the extension reads the same meridian"
|
|
if d > 50000 {
|
|
verdict = "DIVERGENT — signed notation is not the same meridian"
|
|
bad++
|
|
}
|
|
fmt.Fprintf(tw, "%g\t%g\t%.2f km\t%s\n", signed, unsigned, d/1000, verdict)
|
|
}
|
|
}
|
|
tw.Flush()
|
|
return bad
|
|
}
|