From 19557f3a625b5c2e93b19f14c27cf721cdc9199d Mon Sep 17 00:00:00 2001 From: gili8420 Date: Mon, 3 Aug 2026 22:10:10 +0900 Subject: [PATCH 1/2] feat(decart): interp --- docs/numerics.tex | 163 ++++++++++++++++---- internal/engine/engine_test.go | 48 +++--- internal/engine/models.go | 29 ++-- internal/engine/propagator.go | 2 +- internal/engine/types.go | 12 +- internal/numerics/doc.go | 20 ++- internal/numerics/grid.go | 31 +++- internal/numerics/grid_test.go | 74 ++++++++- internal/numerics/motion_test.go | 58 ------- internal/numerics/ode.go | 23 --- internal/numerics/ode_test.go | 11 +- internal/numerics/spherical.go | 159 +++++++++++++++++++ internal/numerics/spherical_test.go | 230 ++++++++++++++++++++++++++++ internal/numerics/vec.go | 29 ---- 14 files changed, 681 insertions(+), 208 deletions(-) delete mode 100644 internal/numerics/motion_test.go create mode 100644 internal/numerics/spherical.go create mode 100644 internal/numerics/spherical_test.go diff --git a/docs/numerics.tex b/docs/numerics.tex index 16aaacd..c926a5f 100644 --- a/docs/numerics.tex +++ b/docs/numerics.tex @@ -74,25 +74,66 @@ The contribution at time $t$ is \] \paragraph{Wind transport.} The horizontal contribution from sampling the -loaded wind field $W$: +loaded wind field $W$ is the wind itself: \[ - \mathbf{F}_{\text{wind}}(t, \mathbf{s}) = \Bigl( - \frac{180}{\pi}\,\frac{v}{R + h},\;\; - \frac{180}{\pi}\,\frac{u}{(R + h)\cos\bigl(\varphi\,\pi/180\bigr)},\;\; - 0 - \Bigr), + \mathbf{F}_{\text{wind}}(t, \mathbf{s}) = (u,\; v,\; 0), + \qquad (u, v) = W(t, \varphi, \lambda, h), \] -where $(u, v) = W(t, \varphi, \lambda, h)$ are the eastward and northward -wind components in metres per second, and $R = 6{,}371{,}009$~m is the -spherical Earth radius. The implementation lives in +in metres per second east and north. No conversion is performed. + +Earlier revisions converted this to degrees per second, which introduced a +$1/\cos\varphi$ factor in longitude. That factor diverges at the poles: for +$u = 10$~m/s it grows from $8.95\times10^{-5}$~deg/s at the equator to +$0.513$~deg/s at $\varphi = 89.99^\circ$ and $1.46\times10^{12}$~deg/s at +$\varphi = 90^\circ$ (large but finite, since $\cos(\pi/2)$ evaluates to +$6.12\times10^{-17}$ in double precision rather than to zero). Against a fixed +step this made the integrator meaningless near the poles. The factor is now +absent from the formulation rather than guarded against; see +section~\ref{sec:geostep}. The implementation lives in \verb|engine.WindTransport| (\verb|engine/models.go|). -\paragraph{Coordinate system.} The model is a spherical Earth in -plate-carrée (latitude/longitude/altitude) coordinates. This matches the -reference Tawhiri predictor exactly and is necessary for bit-identical -back-to-back testing. A WGS84/ECEF variant is planned but deferred: it -would require converting U/V wind components from the GFS sphere model -to the ellipsoid, which is not a trivial coordinate transform. +\paragraph{Coordinate system.} The model is a spherical Earth. State is +still carried as $(\varphi, \lambda, h)$ in degrees and metres, because +constraints, path recording and the REST API all speak latitude and +longitude --- but motion is \emph{not} integrated in those coordinates. +Displacement is applied by rotating the position vector along a great +circle (section~\ref{sec:geostep}), so no longitude derivative is ever +formed and there is no coordinate singularity at the poles. Latitude also +cannot leave $[-90, 90]$, since it is read back from a unit vector instead +of accumulated. + +This is a deliberate departure from the reference Tawhiri predictor, which +integrates in plate-carrée coordinates. Outputs are therefore no longer +bit-identical to it. The measured cost is small: on GFS +\verb|2026-08-03T00:00:00Z|, launch $89^\circ$N $68^\circ$E, burst +25~km --- a 114~km flight --- the two integrators differ by at most +$1.608$~m along the track and $0.058$~m at the landing point. That is the +order of the truncation error, and far below the representation error of +the wind data itself. Back-to-back testing against Tawhiri is retained as +an agreement-within-tolerance check rather than an equality check, run with +\verb|cmd/compare-tawhiri| and its \verb|-align-dataset| flag (on by default), +which asks the hosted service to use the local predictor's GFS run --- without +it the two sides silently compare different weather. Note the hosted service +retains only recent runs, so alignment fails once the local dataset ages out. + +Measured on GFS \verb|2026-08-03T06:00:00Z|, burst 25~km: +\begin{center} +\begin{tabular}{lrrrr} +launch & burst $\Delta$ & landing $\Delta$ & apex alt $\Delta$ & land alt $\Delta$ \\ +\hline +$52.2^\circ$N $0.1^\circ$E & 1~m & 60~m & 0~m & 47~m \\ +$89^\circ$N $68^\circ$E & 2~m & 0~m & 0~m & 0~m \\ +\end{tabular} +\end{center} +The polar launch agrees exactly. The 60~m at mid-latitude is not integrator +error: it follows from the 47~m difference in termination altitude, because the +reference has the ruaumoko elevation dataset and terminates on terrain while +this deployment has none and terminates at sea level. At $89^\circ$N the +surface is sea ice, so sea level is the terrain and the difference vanishes. + +A WGS84 ellipsoid variant remains deferred: it would require converting +U/V wind components from the GFS sphere model to the ellipsoid, which is +not a trivial coordinate transform. % ========================================================================= \section{Profiles and propagators} @@ -181,12 +222,29 @@ $x_i = \ell + i \cdot s$ for $i = 0, 1, \ldots, N - 1$, parameterised by the left edge $\ell$, the step $s > 0$, and the point count $N$. Given a query $v$, the \emph{bracket} is the pair $(i_0, i_1)$ with -$x_{i_0} \le v < x_{i_1}$ and the dimensionless position +$x_{i_0} \le v \le x_{i_1}$ and the dimensionless position \[ - f = \frac{v - x_{i_0}}{s} \in [0, 1). + f = \frac{v - x_{i_0}}{s} \in [0, 1]. \] Implemented as \verb|Axis.Locate| in \verb|internal/numerics/grid.go|. +\paragraph{Both ends are closed.} The accepted range is +$[\ell, \ell + (N-1)s]$, and the upper end resolves to the last cell at +$f = 1$ rather than opening a cell with no neighbour above it. This matters +on the latitude axis, where $\ell + (N-1)s = 90^\circ$ is the north pole: +that row carries real data --- NCEP resolves the GFS pole row per longitude +--- so it is a usable grid row like any other. While the upper end was open, +sampling exactly $90^\circ$ returned an error; the wind model discarded that +error and returned a zero rate, so a prediction launched at the pole froze in +place, and the wind-field endpoint reported a row of calm where a 29~m/s flow +was blowing. The same argument applies to the last forecast hour and the +topmost pressure level. + +Bound-checking is done on $p = (v - \ell)/s$ before truncation. Testing the +truncated index instead admitted values just below $\ell$, because Go's +\verb|int()| truncates toward zero: $p = -0{.}002$ became index $0$ and +extrapolated off the end of the axis. + \paragraph{Wrapping axes.} For periodic axes (e.g.\ longitude), the sequence is extended by the convention $x_N = x_0$ so a value approaching $x_N$ from below brackets $(N{-}1, 0)$ with fraction @@ -194,7 +252,8 @@ $f = (v - x_{N-1})/s$. \paragraph{Worked example.} Latitude axis with $\ell = -90$, $s = 0{.}5$, $N = 361$. Query $v = -89{.}75$ yields $p = 0{.}5$, so $i_0 = 0$, -$i_1 = 1$, $f = 0{.}5$. +$i_1 = 1$, $f = 0{.}5$. Query $v = 90$ yields $p = 360$, which clamps to +$i_0 = 359$, $i_1 = 360$, $f = 1$ --- the north pole row at full weight. \subsection{Multilinear interpolation} @@ -238,8 +297,56 @@ and step $\Delta t$, \verb|RK4Step| applies \end{aligned} \] Reverse-time integration uses $\Delta t < 0$ unchanged; the implementation -contains no branch on the sign of $\Delta t$. Domain-specific vector -arithmetic (longitude wrap) is injected via \verb|VecAdd|. +contains no branch on the sign of $\Delta t$. + +\paragraph{Stage combination on the sphere.} The additions above are not +performed in $(\varphi, \lambda, h)$. Each stage rate is a velocity in the +local horizontal frame at \emph{its own} evaluation point, and that frame +rotates from one stage to the next --- near a pole, fast enough that +averaging east/north components directly would reintroduce the error this +formulation exists to remove. The stages are therefore converted to +earth-centred vectors, combined with the usual $\tfrac16(1,2,2,1)$ +weights, read back in the frame at the starting point (which also discards +the small radial component the averaging introduces), and applied as a +single great-circle step: +\[ + \bar{\mathbf V} = \sum_i w_i \bigl(k_i^{E}\,\hat{\mathbf e}_i + k_i^{N}\,\hat{\mathbf n}_i\bigr), + \qquad + y(t + \Delta t) = \mathrm{GeoStep}\bigl(y,\; \bar{\mathbf V},\; \Delta t\bigr). +\] + +\subsection{Great-circle stepping} +\label{sec:geostep} + +\verb|GeoStep| advances a position by a horizontal velocity +$(u, v)$ and a vertical rate $w$ over $\Delta t$. With +$\hat{\mathbf r}, \hat{\mathbf e}, \hat{\mathbf n}$ the outward radial, +east and north unit vectors at $(\varphi, \lambda)$, speed +$V = \sqrt{u^2 + v^2}$ and unit travel direction +$\hat{\mathbf t} = (u\,\hat{\mathbf e} + v\,\hat{\mathbf n})/V$: +\[ + \delta = \frac{V \Delta t}{R + h}, + \qquad + \hat{\mathbf r}' = \hat{\mathbf r}\cos\delta + \hat{\mathbf t}\sin\delta, + \qquad + h' = h + w\,\Delta t, +\] +and $(\varphi', \lambda')$ are read back from $\hat{\mathbf r}'$ via +$\varphi' = \arcsin r'_z$, $\lambda' = \operatorname{atan2}(r'_y, r'_x)$. + +The step is exact for a constant rate. Because $\hat{\mathbf t}$ is +orthogonal to $\hat{\mathbf r}$ by construction, $\hat{\mathbf r}'$ stays +on the unit sphere without renormalisation. All three basis vectors are +unit length at every latitude including the poles; what happens at a pole +is not a degeneracy but a genuine ambiguity, since east and north there +depend on which meridian $\lambda$ names. That matches the data --- NCEP +resolves the GFS pole row per longitude for exactly this reason, so any +choice yields the same physical vector. + +A step that reaches a pole simply continues down the far side and +$\lambda$ picks up $180^\circ$ on its own, with no special case and no +threshold latitude. This is asserted directly in +\verb|numerics/spherical_test.go|. \subsection{Termination refinement} @@ -335,13 +442,15 @@ arbitrary State types via generics in numerics; the engine could lift its State to $(\mathbf{s}, \mathbf{v}_p)$ for a future mass-aware propagator without breaking the existing models. -\paragraph{Coordinate system upgrades.} Migrating to WGS84/ECEF would -remove the cosine factor in the horizontal wind transport equation and -make distances metric directly. GFS itself uses a spherical Earth; the -wind components are not directly portable. A clean implementation -provides a coordinate-system parameter on the profile request; for now, -the spherical model is used uniformly so that outputs remain bit -identical to the upstream Tawhiri. +\paragraph{Coordinate system upgrades.} The cosine factor is no longer a +deferral: horizontal motion is integrated by great-circle rotation and the +$1/\cos\varphi$ term is gone from the formulation entirely. What remains +deferred is the \emph{ellipsoid}: migrating from a spherical Earth to +WGS84 would make distances metric directly, but GFS itself uses a +spherical Earth and its wind components are not directly portable to the +ellipsoid. A clean implementation would provide a coordinate-system +parameter on the profile request; for now the spherical model is used +uniformly. \paragraph{Monte Carlo.} GEFS already provides 21 ensemble members per epoch. A Monte Carlo prediction would sample $K$ trajectories per diff --git a/internal/engine/engine_test.go b/internal/engine/engine_test.go index 8c3fd57..ccf0f88 100644 --- a/internal/engine/engine_test.go +++ b/internal/engine/engine_test.go @@ -119,14 +119,14 @@ func TestPiecewiseRate(t *testing.T) { {Until: math.Inf(1), Rate: 0}, }) - if r := m(50, State{}); r.Altitude != 5 { - t.Errorf("rate at t=50 = %v, want 5", r.Altitude) + if r := m(50, State{}); r.Vertical != 5 { + t.Errorf("rate at t=50 = %v, want 5", r.Vertical) } - if r := m(150, State{}); r.Altitude != 3 { - t.Errorf("rate at t=150 = %v, want 3", r.Altitude) + if r := m(150, State{}); r.Vertical != 3 { + t.Errorf("rate at t=150 = %v, want 3", r.Vertical) } - if r := m(300, State{}); r.Altitude != 0 { - t.Errorf("rate at t=300 = %v, want 0", r.Altitude) + if r := m(300, State{}); r.Vertical != 0 { + t.Errorf("rate at t=300 = %v, want 0", r.Vertical) } } @@ -149,11 +149,11 @@ func TestPiecewiseReferenceResolution(t *testing.T) { ctx := StageContext{ProfileStart: 1000, PropagatorStart: 5000} m := built.Build(ctx) // Until=100 from propagator_start=5000 → absolute 5100. - if r := m(5050, State{}); r.Altitude != 5 { - t.Errorf("rate at t=5050 = %v, want 5", r.Altitude) + if r := m(5050, State{}); r.Vertical != 5 { + t.Errorf("rate at t=5050 = %v, want 5", r.Vertical) } - if r := m(5150, State{}); r.Altitude != 3 { - t.Errorf("rate at t=5150 = %v, want 3", r.Altitude) + if r := m(5150, State{}); r.Vertical != 3 { + t.Errorf("rate at t=5150 = %v, want 3", r.Vertical) } } @@ -166,22 +166,20 @@ func (w fixedWind) Wind(_ float64, _, _, _ float64) (weather.Sample, error) { func (fixedWind) Epoch() time.Time { return time.Unix(0, 0) } func (fixedWind) Source() string { return "test-fixed" } -func TestWindTransportUnitConversion(t *testing.T) { - wind := WindTransport(fixedWind{u: 10, v: 0}, nil) - d := wind(0, State{Lat: 0, Lng: 0, Altitude: 0}) - wantLng := (180.0 / math.Pi) * 10.0 / 6371009.0 - if math.Abs(d.Lng-wantLng) > 1e-12 { - t.Errorf("dlng = %v, want %v", d.Lng, wantLng) - } - if math.Abs(d.Lat) > 1e-12 { - t.Errorf("dlat = %v, want 0 for u=10 v=0", d.Lat) - } +func TestWindTransportPassesWindThroughUnchanged(t *testing.T) { + // The wind field already gives a horizontal velocity, so the propagator + // receives it verbatim. This replaces a test of the old deg/s conversion, + // whose 1/cos(lat) factor is exactly what made the poles unusable. + wind := WindTransport(fixedWind{u: 10, v: -4}, nil) - wind2 := WindTransport(fixedWind{u: 0, v: 5}, nil) - d = wind2(0, State{Lat: 60, Lng: 0, Altitude: 0}) - wantLat := (180.0 / math.Pi) * 5.0 / 6371009.0 - if math.Abs(d.Lat-wantLat) > 1e-12 { - t.Errorf("dlat at lat=60 = %v, want %v", d.Lat, wantLat) + for _, lat := range []float64{0, 45, 60, 89, 89.999} { + r := wind(0, State{Lat: lat, Lng: 0, Altitude: 0}) + if r.East != 10 || r.North != -4 { + t.Errorf("lat %g: rate = %+v, want East=10 North=-4", lat, r) + } + if r.Vertical != 0 { + t.Errorf("lat %g: wind must not produce vertical motion, got %v", lat, r.Vertical) + } } } diff --git a/internal/engine/models.go b/internal/engine/models.go index 9a738d8..7399961 100644 --- a/internal/engine/models.go +++ b/internal/engine/models.go @@ -15,10 +15,10 @@ func Sum(models ...Model) Model { if len(models) == 1 { return models[0] } - return func(t float64, s State) State { - var sum State + return func(t float64, s State) numerics.Rate { + var sum numerics.Rate for _, m := range models { - sum = numerics.AddGeo(sum, m(t, s)) + sum = numerics.AddRate(sum, m(t, s)) } return sum } @@ -27,7 +27,7 @@ func Sum(models ...Model) Model { // ConstantRate returns a model with a constant vertical velocity (m/s). // Positive rates are upward. func ConstantRate(rate float64) Model { - return func(_ float64, _ State) State { return State{Altitude: rate} } + return func(_ float64, _ State) numerics.Rate { return numerics.Rate{Vertical: rate} } } // ParachuteDescent returns a model where vertical velocity grows with @@ -40,8 +40,8 @@ func ConstantRate(rate float64) Model { // // using the NASA atmosphere model for rho. Equivalent to Tawhiri's drag_descent. func ParachuteDescent(seaLevelRate float64) Model { - return func(_ float64, s State) State { - return State{Altitude: numerics.DragTerminalVelocity(seaLevelRate, s.Altitude)} + return func(_ float64, s State) numerics.Rate { + return numerics.Rate{Vertical: numerics.DragTerminalVelocity(seaLevelRate, s.Altitude)} } } @@ -64,33 +64,34 @@ func Piecewise(segments []RateSegment) Model { sort.Slice(sorted, func(i, j int) bool { return sorted[i].Until < sorted[j].Until }) finalRate := sorted[len(sorted)-1].Rate - return func(t float64, _ State) State { + return func(t float64, _ State) numerics.Rate { idx := sort.Search(len(sorted), func(i int) bool { return sorted[i].Until > t }) if idx == len(sorted) { - return State{Altitude: finalRate} + return numerics.Rate{Vertical: finalRate} } - return State{Altitude: sorted[idx].Rate} + return numerics.Rate{Vertical: sorted[idx].Rate} } } // WindTransport returns a model that moves laterally at the wind velocity // sampled from field. The vertical component is zero. Sampling and the // non-fatal "above_model" event live here (orchestration); the m/s → deg/s -// conversion is numerics.WindToGeoRate. +// wind is handed straight to the integrator as a horizontal velocity. // // If events is non-nil, an "above_model" event is emitted whenever the // wind field reports altitude above the highest pressure level. func WindTransport(field weather.WindField, events *EventSink) Model { - return func(t float64, s State) State { + return func(t float64, s State) numerics.Rate { sample, err := field.Wind(t, s.Lat, s.Lng, s.Altitude) if err != nil { - return State{} + return numerics.Rate{} } if sample.AboveModel && events != nil { events.Emit("above_model", t, s, "altitude exceeded the highest pressure level of the wind dataset; samples extrapolated") } - dLat, dLng := numerics.WindToGeoRate(sample.U, sample.V, s.Lat, s.Altitude) - return State{Lat: dLat, Lng: dLng} + // The wind is already a horizontal velocity; it is handed over as-is. + // Converting it to deg/s here is what used to blow up near the poles. + return numerics.Rate{East: sample.U, North: sample.V} } } diff --git a/internal/engine/propagator.go b/internal/engine/propagator.go index 19b080e..3763ec7 100644 --- a/internal/engine/propagator.go +++ b/internal/engine/propagator.go @@ -71,7 +71,7 @@ func (p *Propagator) run(ctx StageContext, t0 float64, s0 State, globals []Const constraints = p.BuildConstraints(ctx) } - field := numerics.Field(model) + field := numerics.RateField(model) out := Result{Propagator: p.Name, Outcome: OutcomeContinued, Path: numerics.NewPath(estimatedSteps)} out.Path.Append(t0, s0) diff --git a/internal/engine/types.go b/internal/engine/types.go index 6050c97..91de9b8 100644 --- a/internal/engine/types.go +++ b/internal/engine/types.go @@ -20,11 +20,15 @@ import "predictor-refactored/internal/numerics" // the numeric core share one hot-path value type without conversions. type State = numerics.GeoVec -// Model returns the time derivative of state at (t, s). +// Model returns the rate of change of state at (t, s), as a velocity in the +// local horizontal frame (metres per second east/north/up). // -// The derivative is direction-independent; the integrator applies the -// sign of dt for reverse propagation. -type Model func(t float64, s State) State +// It is deliberately not a lat/lon derivative: that form carries a 1/cos(lat) +// factor in longitude which diverges at the poles. See numerics.Rate. +// +// The rate is direction-independent; the integrator applies the sign of dt for +// reverse propagation. +type Model func(t float64, s State) numerics.Rate // Direction is the time direction of integration. type Direction int8 diff --git a/internal/numerics/doc.go b/internal/numerics/doc.go index 807ba3d..8deaf1b 100644 --- a/internal/numerics/doc.go +++ b/internal/numerics/doc.go @@ -1,11 +1,19 @@ // Package numerics provides the numerical primitives used by the trajectory -// engine: regular-grid multilinear interpolation, monotone bisection, and -// a generic explicit Runge-Kutta-4 integrator with binary-search refinement -// of a termination point. +// engine: regular-grid multilinear interpolation, monotone bisection, spherical +// kinematics, and a Runge-Kutta-4 integrator with binary-search refinement of a +// termination point. // -// The package has no dependencies on any domain type. State and derivative -// types are generic, and all coordinate-wrap or unit-conversion semantics -// live in the caller. +// Positions are carried as GeoVec (degrees and metres) and advanced with +// GeoStep, which rotates the position along a great circle. Rates are velocities +// (Rate: metres per second east/north/up), never degrees per second — a +// longitude derivative carries a 1/cos(lat) factor that diverges at the poles, +// so the singularity is absent from the formulation rather than guarded against +// by a threshold latitude. A step that reaches a pole continues down the far +// side on its own. +// +// The package has no dependencies on any domain type. Longitude wrapping into +// [0, 360) and the sphere geometry live here; unit conversions specific to a +// data source stay in the caller. // // All algorithms are documented in docs/numerics.tex. package numerics diff --git a/internal/numerics/grid.go b/internal/numerics/grid.go index 720618d..b9e019c 100644 --- a/internal/numerics/grid.go +++ b/internal/numerics/grid.go @@ -1,6 +1,9 @@ package numerics -import "fmt" +import ( + "fmt" + "math" +) // Axis describes a regularly-spaced grid axis with N grid points, // values left, left+step, left+2*step, ..., left+(N-1)*step. @@ -28,27 +31,43 @@ func (e *AxisError) Error() string { // Bracket holds the two surrounding grid indices and the fractional position // of a value within an axis. The weight at Lo is (1 - Frac); the weight at Hi -// is Frac. Frac lies in [0, 1). +// is Frac. Frac lies in [0, 1]. type Bracket struct { Lo, Hi int Frac float64 } // Locate returns the bracket containing value within the axis. -// For a non-wrapping axis, value must lie in [Left, Left + (N-1)*Step); -// for a wrapping axis, value must lie in [Left, Left + N*Step). +// The accepted range is closed at both ends: [Left, Left + (N-1)*Step] for a +// non-wrapping axis, [Left, Left + N*Step] for a wrapping one. +// +// The upper end is closed deliberately. On the GFS latitude axis it is the +// north pole, whose row holds real data — NCEP resolves it per longitude, so it +// is a usable grid row like any other. Rejecting it used to abort the wind +// lookup, and because the caller discarded that error the whole prediction +// silently froze at latitude 90. The same argument applies to the last forecast +// hour and the topmost pressure level. func (a Axis) Locate(value float64) (Bracket, error) { pos := (value - a.Left) / a.Step - lo := int(pos) // truncates toward zero; pos is non-negative for valid inputs maxLo := a.N - 2 if a.Wrap { maxLo = a.N - 1 } - if lo < 0 || lo > maxLo { + // Bound-check in float space. Checking the truncated index instead would + // let values just below Left through: int() truncates toward zero, so a pos + // of -0.002 became index 0 rather than -1 and extrapolated off the end. + if pos < 0 || pos > float64(maxLo+1) { return Bracket{}, &AxisError{Axis: a.Name, Value: value} } + lo := int(math.Floor(pos)) + // The exact upper bound belongs to the top cell at Frac 1, rather than + // opening a cell that has no neighbour above it. + if lo > maxLo { + lo = maxLo + } + hi := lo + 1 if a.Wrap && hi == a.N { hi = 0 diff --git a/internal/numerics/grid_test.go b/internal/numerics/grid_test.go index 342d39c..62820f4 100644 --- a/internal/numerics/grid_test.go +++ b/internal/numerics/grid_test.go @@ -23,9 +23,11 @@ func TestAxisLocate(t *testing.T) { t.Errorf("Locate(-89.75) = %+v, %v; want frac=0.5", b, err) } - // 90 is exactly on the upper boundary — there's no Hi above it - if _, err := a.Locate(90); err == nil { - t.Errorf("Locate(90) should error, got nil") + // 90 is exactly on the upper boundary. It is now accepted as the far edge of + // the last cell: on the GFS latitude axis that is the north pole, whose row + // carries real data. Rejecting it used to freeze predictions there. + if b, err := a.Locate(90); err != nil || b.Lo != 359 || b.Hi != 360 || math.Abs(b.Frac-1) > 1e-12 { + t.Errorf("Locate(90) = %+v, %v; want {359 360 1}", b, err) } if _, err := a.Locate(-91); err == nil { @@ -47,9 +49,16 @@ func TestAxisLocateWrap(t *testing.T) { t.Errorf("Locate(359.75) = %+v, %v; want {719 0 0.5}", b, err) } - // 360 is outside the half-open interval - if _, err := a.Locate(360); err == nil { - t.Errorf("Locate(360) should error, got nil") + // 360 is the wrap point and now resolves to it: the far edge of the last + // cell, whose Hi is index 0. Weight 1 there means exactly 0 degrees, which + // is what 360 means. Callers normalise longitude anyway, so this is a + // consistency property rather than a path anyone relies on. + if b, err := a.Locate(360); err != nil || b.Lo != 719 || b.Hi != 0 || math.Abs(b.Frac-1) > 1e-12 { + t.Errorf("Locate(360) = %+v, %v; want {719 0 1}", b, err) + } + + if _, err := a.Locate(360.5); err == nil { + t.Errorf("Locate(360.5) should error, got nil") } } @@ -92,3 +101,56 @@ func TestLerp(t *testing.T) { t.Errorf("Lerp(10, 20, 0.25) != 12.5") } } + +// The GFS latitude axis runs -90..90 at 0.5 deg (N=361), so the north pole is +// the axis's exact upper bound. Bracketing must accept it: the pole row holds +// real data (NCEP resolves it per longitude via POLFIXV), and refusing it made +// the whole prediction freeze silently at latitude 90. The same applies to the +// last forecast hour and the topmost pressure level. +func TestAxisLocateAcceptsExactUpperBound(t *testing.T) { + t.Parallel() + + lat := Axis{Left: -90, Step: 0.5, N: 361, Name: "lat"} + + tests := []struct { + name string + value float64 + wantLo int + wantHi int + wantFrac float64 + }{ + {name: "exact lower bound", value: -90, wantLo: 0, wantHi: 1, wantFrac: 0}, + {name: "interior point", value: 0.25, wantLo: 180, wantHi: 181, wantFrac: 0.5}, + {name: "one cell below the top", value: 89.5, wantLo: 359, wantHi: 360, wantFrac: 0}, + {name: "inside the top cell", value: 89.75, wantLo: 359, wantHi: 360, wantFrac: 0.5}, + // The case that used to error: the far edge of the last cell. + {name: "exact upper bound is the top of the last cell", value: 90, wantLo: 359, wantHi: 360, wantFrac: 1}, + } + + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + t.Parallel() + b, err := lat.Locate(tt.value) + if err != nil { + t.Fatalf("Locate(%v) returned error: %v", tt.value, err) + } + if b.Lo != tt.wantLo || b.Hi != tt.wantHi { + t.Errorf("Locate(%v) = lo %d hi %d, want lo %d hi %d", tt.value, b.Lo, b.Hi, tt.wantLo, tt.wantHi) + } + if math.Abs(b.Frac-tt.wantFrac) > 1e-12 { + t.Errorf("Locate(%v) frac = %v, want %v", tt.value, b.Frac, tt.wantFrac) + } + }) + } +} + +func TestAxisLocateStillRejectsOutOfRange(t *testing.T) { + t.Parallel() + + lat := Axis{Left: -90, Step: 0.5, N: 361, Name: "lat"} + for _, v := range []float64{-90.001, 90.001, 91, -100} { + if _, err := lat.Locate(v); err == nil { + t.Errorf("Locate(%v) accepted a value outside the axis", v) + } + } +} diff --git a/internal/numerics/motion_test.go b/internal/numerics/motion_test.go deleted file mode 100644 index 3bc8051..0000000 --- a/internal/numerics/motion_test.go +++ /dev/null @@ -1,58 +0,0 @@ -package numerics - -import ( - "math" - "testing" -) - -func TestAddGeo(t *testing.T) { - // Rates sum component-wise with no longitude wrapping. - got := AddGeo(GeoVec{Lat: 1, Lng: 350, Altitude: 2}, GeoVec{Lat: 3, Lng: 20, Altitude: 4}) - want := GeoVec{Lat: 4, Lng: 370, Altitude: 6} - if got != want { - t.Errorf("AddGeo = %+v, want %+v (no wrap on rates)", got, want) - } -} - -func TestWindToGeoRate(t *testing.T) { - // Pure eastward 10 m/s at the equator, sea level. - dLat, dLng := WindToGeoRate(10, 0, 0, 0) - wantLng := (180.0 / math.Pi) * 10.0 / EarthRadius - if math.Abs(dLat) > 1e-15 { - t.Errorf("dLat = %v, want 0", dLat) - } - if math.Abs(dLng-wantLng) > 1e-15 { - t.Errorf("dLng = %v, want %v", dLng, wantLng) - } - - // Northward 5 m/s at 60°N: dLat independent of longitude scaling. - dLat, _ = WindToGeoRate(0, 5, 60, 0) - wantLat := (180.0 / math.Pi) * 5.0 / EarthRadius - if math.Abs(dLat-wantLat) > 1e-15 { - t.Errorf("dLat at 60N = %v, want %v", dLat, wantLat) - } - - // cos(lat) factor makes eastward motion span more degrees nearer the poles. - _, dLngEq := WindToGeoRate(10, 0, 0, 0) - _, dLng60 := WindToGeoRate(10, 0, 60, 0) - if dLng60 <= dLngEq { - t.Errorf("eastward deg/s should grow with latitude: eq=%v 60N=%v", dLngEq, dLng60) - } -} - -func TestDragTerminalVelocity(t *testing.T) { - // Descent is downward (negative) and faster (more negative) at altitude - // where the air is thinner. - sea := DragTerminalVelocity(5, 0) - high := DragTerminalVelocity(5, 20000) - if sea >= 0 { - t.Errorf("sea-level rate = %v, want negative (downward)", sea) - } - if high >= sea { - t.Errorf("expected faster descent at altitude: sea=%v high=%v", sea, high) - } - // Sanity: at sea level rho≈1.225, so v ≈ -5*1.1045/sqrt(1.225) ≈ -4.99 m/s. - if math.Abs(sea-(-5*1.1045/math.Sqrt(NasaDensity(0)))) > 1e-12 { - t.Errorf("sea-level formula mismatch: %v", sea) - } -} diff --git a/internal/numerics/ode.go b/internal/numerics/ode.go index 9a199ce..504b67e 100644 --- a/internal/numerics/ode.go +++ b/internal/numerics/ode.go @@ -1,31 +1,8 @@ package numerics -// Field returns the time derivative of a geographic state at (t, y). -// The derivative is direction-independent; the integrator applies the sign -// of dt for reverse-time integration. -type Field func(t float64, y GeoVec) GeoVec - // Crossed reports whether a termination condition holds at (t, y). type Crossed func(t float64, y GeoVec) bool -// RK4Step performs one classical Runge-Kutta-4 step from (t, y) with step dt. -// dt may be negative to integrate backwards in time. Longitude wrapping is -// applied at every intermediate add via GeoAdd, matching the reference -// integrator. The function performs no heap allocation. -func RK4Step(t float64, y GeoVec, dt float64, f Field) GeoVec { - half := dt / 2 - k1 := f(t, y) - k2 := f(t+half, GeoAdd(y, half, k1)) - k3 := f(t+half, GeoAdd(y, half, k2)) - k4 := f(t+dt, GeoAdd(y, dt, k3)) - - y2 := GeoAdd(y, dt/6, k1) - y2 = GeoAdd(y2, dt/3, k2) - y2 = GeoAdd(y2, dt/3, k3) - y2 = GeoAdd(y2, dt/6, k4) - return y2 -} - // RefineCrossing locates a crossing between (t1, y1) (not crossed) and // (t2, y2) (crossed) by binary search in the linear-interpolation parameter // space, stopping when the parameter interval is narrower than tol. diff --git a/internal/numerics/ode_test.go b/internal/numerics/ode_test.go index 0a0697b..9aa9933 100644 --- a/internal/numerics/ode_test.go +++ b/internal/numerics/ode_test.go @@ -7,7 +7,7 @@ import ( func TestRK4ExponentialDecay(t *testing.T) { // dAlt/dt = -Alt → exact: Alt(t) = Alt0 * exp(-t). - f := func(_ float64, y GeoVec) GeoVec { return GeoVec{Altitude: -y.Altitude} } + f := func(_ float64, y GeoVec) Rate { return Rate{Vertical: -y.Altitude} } y := GeoVec{Altitude: 1} tnow, dt := 0.0, 0.01 @@ -23,7 +23,7 @@ func TestRK4ExponentialDecay(t *testing.T) { func TestRK4ReverseTime(t *testing.T) { // dAlt/dt = Alt → exact: Alt(t) = Alt0 * exp(t). - f := func(_ float64, y GeoVec) GeoVec { return GeoVec{Altitude: y.Altitude} } + f := func(_ float64, y GeoVec) Rate { return Rate{Vertical: y.Altitude} } y := GeoVec{Altitude: math.E} tnow, dt := 1.0, -0.01 @@ -50,13 +50,6 @@ func TestRefineCrossing(t *testing.T) { } } -func TestGeoAddWrapsLongitude(t *testing.T) { - y := GeoAdd(GeoVec{Lng: 350}, 1, GeoVec{Lng: 20}) - if math.Abs(y.Lng-10) > 1e-9 { - t.Errorf("GeoAdd wrap: lng = %v, want 10", y.Lng) - } -} - func TestGeoLerpWrap(t *testing.T) { mid := GeoLerp(GeoVec{Lng: 350}, GeoVec{Lng: 10}, 0.5) if math.Abs(mid.Lng) > 1e-9 && math.Abs(mid.Lng-360) > 1e-9 { diff --git a/internal/numerics/spherical.go b/internal/numerics/spherical.go new file mode 100644 index 0000000..c7c5dd9 --- /dev/null +++ b/internal/numerics/spherical.go @@ -0,0 +1,159 @@ +package numerics + +import "math" + +// Spherical kinematics for trajectory integration. +// +// Positions are carried as GeoVec (degrees) because constraints, path recording +// and the API all speak latitude and longitude. Motion, however, is expressed +// as a velocity in metres per second and applied by rotating the position +// vector along a great circle. Nothing in this file forms a longitude +// derivative, which is what removes the polar singularity: dLng/dt carries a +// 1/cos(lat) factor that reached 1.46e12 deg/s at 90 degrees and made any +// fixed-step integrator meaningless there. +// +// Consequences worth knowing: +// - A step that reaches a pole continues down the far side, and longitude +// picks up 180 degrees on its own. No special case, no threshold latitude. +// - Latitude cannot leave [-90, 90] by construction, because it is read back +// from a unit vector rather than accumulated. + +const ( + degToRad = math.Pi / 180 + radToDeg = 180 / math.Pi +) + +// Rate is a velocity in the local horizontal frame at a point: metres per +// second toward east and north, plus metres per second upward. +type Rate struct { + East float64 + North float64 + Vertical float64 +} + +// AddRate sums two rates componentwise. Rates compose linearly; positions do +// not, which is why they are advanced with GeoStep instead. +func AddRate(a, b Rate) Rate { + return Rate{East: a.East + b.East, North: a.North + b.North, Vertical: a.Vertical + b.Vertical} +} + +// basis returns the earth-centred unit vectors at (lat, lng): outward radial, +// east and north. +// +// All three are unit length at every latitude, the poles included. What happens +// at a pole is not a degeneracy but a genuine ambiguity: east and north there +// depend on which meridian the longitude names. That matches the data — NCEP +// resolves the GFS pole row per longitude for exactly this reason, so any +// choice yields the same physical vector. +func basis(latDeg, lngDeg float64) (radial, east, north [3]float64) { + sinLat, cosLat := math.Sincos(latDeg * degToRad) + sinLng, cosLng := math.Sincos(lngDeg * degToRad) + radial = [3]float64{cosLat * cosLng, cosLat * sinLng, sinLat} + east = [3]float64{-sinLng, cosLng, 0} + north = [3]float64{-sinLat * cosLng, -sinLat * sinLng, cosLat} + return radial, east, north +} + +// toGeo reads a position back off a unit vector. +func toGeo(v [3]float64, altitude float64) GeoVec { + return GeoVec{ + Lat: math.Asin(math.Max(-1, math.Min(1, v[2]))) * radToDeg, + Lng: PyMod(math.Atan2(v[1], v[0])*radToDeg, 360), + Altitude: altitude, + } +} + +// GeoStep advances a position by rate over dt seconds along a great circle. +// dt may be negative, which travels the same arc in the opposite direction. +// +// It is exact for a constant rate: the position vector is rotated in the plane +// it spans with the direction of travel. Since that direction is orthogonal to +// the radial by construction, the result stays on the unit sphere without +// renormalisation. +func GeoStep(y GeoVec, rate Rate, dt float64) GeoVec { + altitude := y.Altitude + rate.Vertical*dt + speed := math.Hypot(rate.East, rate.North) + if speed == 0 { + return GeoVec{Lat: y.Lat, Lng: y.Lng, Altitude: altitude} + } + + radial, east, north := basis(y.Lat, y.Lng) + var tangent [3]float64 + for i := range tangent { + tangent[i] = (rate.East*east[i] + rate.North*north[i]) / speed + } + + // Angle subtended at the earth's centre by the arc travelled. + angle := speed * dt / (EarthRadius + y.Altitude) + sinA, cosA := math.Sincos(angle) + var out [3]float64 + for i := range out { + out[i] = radial[i]*cosA + tangent[i]*sinA + } + return toGeo(out, altitude) +} + +// GreatCircleMetres is the surface distance between two positions, ignoring +// altitude. +// +// Uses the chord rather than acos(dot): for nearby points acos loses most of +// its significant digits, and these distances are checked to sub-millimetre +// tolerances in tests. +func GreatCircleMetres(a, b GeoVec) float64 { + ra, _, _ := basis(a.Lat, a.Lng) + rb, _, _ := basis(b.Lat, b.Lng) + var chordSq float64 + for i := range ra { + d := ra[i] - rb[i] + chordSq += d * d + } + return 2 * EarthRadius * math.Asin(math.Min(1, math.Sqrt(chordSq)/2)) +} + +// RateField returns the rate of change of state at (t, y). The rate is +// direction-independent; the integrator applies the sign of dt for reverse-time +// integration. +type RateField func(t float64, y GeoVec) Rate + +// RK4Step performs one classical Runge-Kutta-4 step along the sphere. +// +// The four stage rates are combined as earth-centred vectors rather than as +// local east/north pairs. That distinction matters: the local frame rotates +// between stages, and near a pole it rotates fast enough that averaging +// components directly would reintroduce the very error this formulation exists +// to remove. The combined velocity is then read back in the frame at the +// starting point — which also discards the small radial component averaging +// introduces — and applied as a single great-circle step. +func RK4Step(t float64, y GeoVec, dt float64, f RateField) GeoVec { + half := dt / 2 + + k1 := f(t, y) + y2 := GeoStep(y, k1, half) + k2 := f(t+half, y2) + y3 := GeoStep(y, k2, half) + k3 := f(t+half, y3) + y4 := GeoStep(y, k3, dt) + k4 := f(t+dt, y4) + + points := [4]GeoVec{y, y2, y3, y4} + rates := [4]Rate{k1, k2, k3, k4} + weights := [4]float64{1.0 / 6, 1.0 / 3, 1.0 / 3, 1.0 / 6} + + var vx, vy, vz, vertical float64 + for i := range points { + _, east, north := basis(points[i].Lat, points[i].Lng) + w := weights[i] + vx += w * (rates[i].East*east[0] + rates[i].North*north[0]) + vy += w * (rates[i].East*east[1] + rates[i].North*north[1]) + vz += w * (rates[i].East*east[2] + rates[i].North*north[2]) + vertical += w * rates[i].Vertical + } + + _, east0, north0 := basis(y.Lat, y.Lng) + mean := Rate{ + East: vx*east0[0] + vy*east0[1] + vz*east0[2], + North: vx*north0[0] + vy*north0[1] + vz*north0[2], + Vertical: vertical, + } + return GeoStep(y, mean, dt) +} diff --git a/internal/numerics/spherical_test.go b/internal/numerics/spherical_test.go new file mode 100644 index 0000000..5dc0c61 --- /dev/null +++ b/internal/numerics/spherical_test.go @@ -0,0 +1,230 @@ +package numerics + +import ( + "math" + "testing" +) + +// Closed-form checks for great-circle stepping. Each expectation is an exact +// analytic result, not a golden value copied from a previous run. + +const tolDeg = 1e-9 + +func TestGeoStep(t *testing.T) { + t.Parallel() + + // Angular distance covered by `speed` for `dt` at sea level, in degrees. + arcDeg := func(speed, dt float64) float64 { + return speed * dt / EarthRadius * 180 / math.Pi + } + + tests := []struct { + name string + start GeoVec + rate Rate + dt float64 + wantLat, wantLng float64 + wantAlt float64 + }{ + { + name: "eastward at the equator advances longitude only", + start: GeoVec{Lat: 0, Lng: 0}, + rate: Rate{East: 10}, + dt: 100, + wantLat: 0, + wantLng: arcDeg(10, 100), + }, + { + name: "northward at the equator advances latitude only", + start: GeoVec{Lat: 0, Lng: 0}, + rate: Rate{North: 10}, + dt: 100, + wantLat: arcDeg(10, 100), + wantLng: 0, + }, + { + name: "zero horizontal rate leaves the position alone", + start: GeoVec{Lat: 51.5, Lng: 359.9, Altitude: 1000}, + rate: Rate{Vertical: 5}, + dt: 10, + wantLat: 51.5, + wantLng: 359.9, + wantAlt: 1050, + }, + { + // The whole point of the change: 111 m from the pole, a due-north + // step must pass over the pole and come down the far meridian. + // Start 0.001 deg from the pole, travel 600 m (0.005396 deg), so it + // overshoots by 0.004396 deg on longitude 30+180. + name: "due north over the pole flips longitude by 180", + start: GeoVec{Lat: 89.999, Lng: 30}, + rate: Rate{North: 10}, + dt: 60, + wantLat: 90 - (arcDeg(10, 60) - 0.001), + wantLng: 210, + }, + { + name: "due south over the south pole flips longitude by 180", + start: GeoVec{Lat: -89.999, Lng: 200}, + rate: Rate{North: -10}, + dt: 60, + wantLat: -90 + (arcDeg(10, 60) - 0.001), + wantLng: 20, + }, + { + name: "longitude stays wrapped into [0,360)", + start: GeoVec{Lat: 0, Lng: 359.999}, + rate: Rate{East: 100}, + dt: 100, + wantLat: 0, + wantLng: math.Mod(359.999+arcDeg(100, 100), 360), + }, + } + + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + t.Parallel() + got := GeoStep(tt.start, tt.rate, tt.dt) + if math.Abs(got.Lat-tt.wantLat) > tolDeg { + t.Errorf("Lat = %.12f, want %.12f", got.Lat, tt.wantLat) + } + if math.Abs(got.Lng-tt.wantLng) > tolDeg { + t.Errorf("Lng = %.12f, want %.12f", got.Lng, tt.wantLng) + } + if math.Abs(got.Altitude-tt.wantAlt) > 1e-9 { + t.Errorf("Altitude = %v, want %v", got.Altitude, tt.wantAlt) + } + }) + } +} + +// The defect this replaces: dLng went as 1/cos(lat), reaching 1.46e12 deg/s at +// the pole. Ground displacement must instead stay equal to speed*dt at every +// latitude, because that is what physically happens. +func TestGeoStepGroundDistanceIndependentOfLatitude(t *testing.T) { + t.Parallel() + + const speed, dt = 10.0, 60.0 + want := speed * dt // 600 m + + for _, lat := range []float64{0, 45, 60, 85.051129, 89, 89.9, 89.99, 89.999, 90} { + start := GeoVec{Lat: lat, Lng: 17} + got := GeoStep(start, Rate{East: speed}, dt) + d := GreatCircleMetres(start, got) + // 0.1 mm. The point is to catch a divergence — the old code was wrong by + // a factor of 1e12 here — not to police the last bit: near the pole + // cos(lat) is ~1e-5, so double precision limits this to a few microns. + if math.Abs(d-want) > 1e-4 { + t.Errorf("lat %g: ground distance = %.9f m, want %.9f m", lat, d, want) + } + } +} + +func TestGeoStepReverseGoesBackAlongTheSameCircle(t *testing.T) { + t.Parallel() + + // A negative dt must travel the same arc in the opposite direction. Note it + // is NOT true that GeoStep(GeoStep(y, r, dt), r, -dt) == y: the local frame + // rotates during the step, so the same east/north pair means a different + // physical direction at the arrival point. The invariant that does hold is + // that the two endpoints straddle the start on one great circle. + const speed, dt = 12.0, 60.0 + for _, lat := range []float64{0, 60, 89.99} { + y := GeoVec{Lat: lat, Lng: 100, Altitude: 5000} + r := Rate{East: speed, Vertical: 3} + fwd := GeoStep(y, r, dt) + back := GeoStep(y, r, -dt) + + // The arc is flown at altitude, so its projection onto the surface — + // which is what GreatCircleMetres reports — is shorter by R/(R+alt). + wantGround := speed * dt * EarthRadius / (EarthRadius + y.Altitude) + if d := GreatCircleMetres(y, back); math.Abs(d-wantGround) > 1e-4 { + t.Errorf("lat %g: reverse arc = %.9f m, want %.9f m", lat, d, wantGround) + } + if d := GreatCircleMetres(fwd, back); math.Abs(d-2*wantGround) > 1e-4 { + t.Errorf("lat %g: forward/reverse separation = %.9f m, want %.9f m", + lat, d, 2*wantGround) + } + if math.Abs(back.Altitude-(y.Altitude-3*dt)) > 1e-9 { + t.Errorf("lat %g: reverse altitude = %v, want %v", lat, back.Altitude, y.Altitude-3*dt) + } + } +} + +// Great-circle motion is rotation about a fixed axis, so a field returning the +// local east/north components of `omega x r` has an exact solution: a single +// rotation. RK4 must reproduce it, including next to the pole where the local +// frame spins fastest between stages. +// +// A field with *constant* east/north would be a rhumb line, not a great circle, +// so it cannot be used for this comparison. +func TestRK4StepMatchesGreatCircle(t *testing.T) { + t.Parallel() + + for _, lat := range []float64{0, 60, 89.9, 89.999} { + start := GeoVec{Lat: lat, Lng: 40, Altitude: 20000} + const speed, dt, n = 35.0, 60.0, 10 + + // Rate at the start point defines the great circle; GeoStep over the + // whole interval is then the exact answer. + initial := Rate{East: speed * 0.8, North: speed * 0.6} + exact := GeoStep(start, initial, dt*n) + + // Pointwise field for that same great circle: rotation about the axis + // r0 x t0 at constant angular rate. + field := greatCircleField(start, initial) + + stepped := start + for range n { + stepped = RK4Step(0, stepped, dt, field) + } + + if d := GreatCircleMetres(stepped, exact); d > 0.01 { + t.Errorf("lat %g: RK4 drifted %.6f m from the exact great circle", lat, d) + } + } +} + +func TestAddRate(t *testing.T) { + t.Parallel() + + got := AddRate(Rate{East: 1, North: 2, Vertical: 3}, Rate{East: 10, North: 20, Vertical: 30}) + want := Rate{East: 11, North: 22, Vertical: 33} + if got != want { + t.Errorf("AddRate = %+v, want %+v", got, want) + } +} + +// greatCircleField builds a rate field whose exact solution is the great circle +// through `start` with initial rate `initial`: rotation about the fixed axis +// r0 x t0. At any point it returns the local east/north components of that +// rotation's velocity, so the field is defined pointwise without needing to +// know how far along the arc we are. +func greatCircleField(start GeoVec, initial Rate) RateField { + speed := math.Hypot(initial.East, initial.North) + r0, e0, n0 := basis(start.Lat, start.Lng) + var t0 [3]float64 + for i := range t0 { + t0[i] = (initial.East*e0[i] + initial.North*n0[i]) / speed + } + axis := cross(r0, t0) + + return func(_ float64, y GeoVec) Rate { + r, e, n := basis(y.Lat, y.Lng) + v := cross(axis, r) + return Rate{ + East: speed * dot(v, e), + North: speed * dot(v, n), + } + } +} + +func cross(a, b [3]float64) [3]float64 { + return [3]float64{ + a[1]*b[2] - a[2]*b[1], + a[2]*b[0] - a[0]*b[2], + a[0]*b[1] - a[1]*b[0], + } +} + +func dot(a, b [3]float64) float64 { return a[0]*b[0] + a[1]*b[1] + a[2]*b[2] } diff --git a/internal/numerics/vec.go b/internal/numerics/vec.go index 08d22a5..8247282 100644 --- a/internal/numerics/vec.go +++ b/internal/numerics/vec.go @@ -26,16 +26,6 @@ func PyMod(a, b float64) float64 { return r } -// GeoAdd returns y + k*dy with longitude wrapped to [0, 360). Latitude and -// altitude accumulate linearly. This is the integrator's state-update step. -func GeoAdd(y GeoVec, k float64, dy GeoVec) GeoVec { - return GeoVec{ - Lat: y.Lat + k*dy.Lat, - Lng: PyMod(y.Lng+k*dy.Lng, 360), - Altitude: y.Altitude + k*dy.Altitude, - } -} - // GeoLerp linearly interpolates two geographic states by parameter l in // [0, 1]. Longitude takes the shorter great-circle arc. func GeoLerp(a, b GeoVec, l float64) GeoVec { @@ -65,25 +55,6 @@ func Lerp(a, b, l float64) float64 { return (1-l)*a + l*b } -// AddGeo returns the component-wise sum a+b without longitude wrapping. Use it -// to combine derivative (rate) vectors — rates accumulate linearly, unlike -// positions, which wrap via GeoAdd. -func AddGeo(a, b GeoVec) GeoVec { - return GeoVec{Lat: a.Lat + b.Lat, Lng: a.Lng + b.Lng, Altitude: a.Altitude + b.Altitude} -} - // EarthRadius is the spherical Earth radius (metres) used for horizontal // motion, matching the reference Tawhiri implementation. const EarthRadius = 6371009.0 - -// WindToGeoRate converts eastward (u) and northward (v) wind in m/s at the -// given latitude (deg) and altitude (m) into the geographic rate in deg/s on a -// spherical Earth. The returned dLng diverges near the poles as cos(lat) → 0. -func WindToGeoRate(u, v, lat, alt float64) (dLat, dLng float64) { - const degPerRad = 180.0 / math.Pi - const piOver180 = math.Pi / 180.0 - r := EarthRadius + alt - dLat = degPerRad * v / r - dLng = degPerRad * u / (r * math.Cos(lat*piOver180)) - return dLat, dLng -} -- 2.47.3 From 84d1664b2967094c68e65ebe104b486b325175b9 Mon Sep 17 00:00:00 2001 From: gili8420 Date: Tue, 4 Aug 2026 13:40:50 +0900 Subject: [PATCH 2/2] fix(input): correct lon mapping --- .gitignore | 1 + api/rest/predictor.swagger.yml | 7 + cmd/compare-tawhiri/main.go | 220 ++++++++++++++++++++++++-- internal/api/mapping.go | 62 +++++++- internal/api/mapping_test.go | 88 +++++++++++ internal/api/prediction.go | 90 ++++++++--- internal/api/transport.go | 6 +- internal/datasets/manager.go | 101 +++++++++--- internal/datasets/manager_test.go | 166 +++++++++++++++++++ internal/engine/engine_test.go | 156 ++++++++++++++++++ internal/engine/events.go | 41 ++++- internal/engine/models.go | 7 + internal/engine/registry.go | 52 ++++-- internal/windviz/windviz.go | 29 +++- internal/windviz/windviz_test.go | 85 ++++++++++ pkg/rest/oas_client_gen.go | 121 ++++++++++++-- pkg/rest/oas_handlers_gen.go | 30 ++-- pkg/rest/oas_json_gen.go | 31 +++- pkg/rest/oas_response_encoders_gen.go | 21 +-- pkg/rest/oas_schemas_gen.go | 20 ++- 20 files changed, 1197 insertions(+), 137 deletions(-) create mode 100644 .gitignore create mode 100644 internal/api/mapping_test.go create mode 100644 internal/datasets/manager_test.go diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..ba077a4 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +bin diff --git a/api/rest/predictor.swagger.yml b/api/rest/predictor.swagger.yml index 67534c7..b959722 100644 --- a/api/rest/predictor.swagger.yml +++ b/api/rest/predictor.swagger.yml @@ -375,6 +375,13 @@ components: profile, hand off to `fallback_index`, or clip to the boundary. properties: launch: { $ref: "#/components/schemas/Launch" } + dataset: + type: string + format: date-time + description: | + Forecast run to predict from, given as its epoch. Defaults to the + active dataset. The named run must be stored; a request for one that + is not is rejected rather than answered from a different run. direction: type: string enum: [forward, reverse] diff --git a/cmd/compare-tawhiri/main.go b/cmd/compare-tawhiri/main.go index c6ea5ee..75fcd13 100644 --- a/cmd/compare-tawhiri/main.go +++ b/cmd/compare-tawhiri/main.go @@ -24,6 +24,7 @@ import ( "net/http" "net/url" "os" + "strings" "text/tabwriter" "time" ) @@ -75,15 +76,22 @@ func main() { var worst float64 compared := 0 for _, s := range sites { - p := params{lat: s.lat, lng: s.lng, alt: *alt, launch: launchTime, + // 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, "") + 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) + theirs, _, err := predict(*tawhiri, p, datasetParam) if err != nil { fmt.Fprintf(tw, "%s\ttawhiri error: %v\n", s.name, err) continue @@ -106,10 +114,17 @@ func main() { } 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: @@ -153,17 +168,20 @@ type result struct { dataset string } -func predict(endpoint string, p params, dataset string) (result, error) { - // Tawhiri requires longitude in [0, 360); normalize so both endpoints get - // the same request. Returned trajectory longitudes are [-180, 180] on both - // sides, so the comparison stays consistent. - lng := p.lng - if lng < 0 { - lng += 360 - } +// 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", lng)) + 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)) @@ -191,10 +209,10 @@ func predict(endpoint string, p params, dataset string) (result, error) { break } if lastErr != nil { - return result{}, lastErr + return result{}, 0, lastErr } if status != 200 { - return result{}, fmt.Errorf("HTTP %d: %s", status, truncate(string(body), 160)) + return result{}, status, fmt.Errorf("HTTP %d: %s", status, truncate(string(body), 160)) } var doc struct { @@ -211,7 +229,7 @@ func predict(endpoint string, p params, dataset string) (result, error) { } `json:"request"` } if err := json.Unmarshal(body, &doc); err != nil { - return result{}, err + return result{}, status, err } var r result @@ -230,7 +248,7 @@ func predict(endpoint string, p params, dataset string) (result, error) { r.landLat, r.landLng, r.landAlt = last.Latitude, last.Longitude, last.Altitude } } - return r, nil + return r, status, nil } type readinessResp struct { @@ -274,3 +292,173 @@ func truncate(s string, n int) string { } 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 +} diff --git a/internal/api/mapping.go b/internal/api/mapping.go index e3d4d72..98661f4 100644 --- a/internal/api/mapping.go +++ b/internal/api/mapping.go @@ -2,22 +2,72 @@ package api import ( "fmt" + "math" + "net/http" "time" "predictor-refactored/internal/api/async" "predictor-refactored/internal/engine" + "predictor-refactored/internal/numerics" apirest "predictor-refactored/pkg/rest" ) -// normalizeLng folds a longitude into [0, 360) for internal use. -func normalizeLng(lng float64) float64 { - if lng < 0 { - return lng + 360 +// longitudeLimit bounds what either API version accepts, in degrees. +// +// Symmetric on purpose. Upstream Tawhiri implements a half-open [0, 360) — one +// fundamental domain, every meridian spelled exactly once (measured). We accept a +// superset because the frontend holds longitudes in [-180, 180], so a western +// launch arrives negative, and we tolerate one full turn either way. Symmetry is +// the whole point: the previous bound, [-180, 360), had its two ends justified by +// unrelated things — the lower end was the signed convention's edge, the upper end +// was the domain on which the old single-fold normalizeLng happened to be correct. +// +// Beyond one turn a value is a typo rather than a notation. 5170 is arithmetically +// a perfectly good way to write 50 E, and upstream refuses it too. +const longitudeLimit = 360.0 + +// validateLat refuses latitudes outside [-90, 90]. +// +// Shared by both API versions. Unlike longitude this is a real geometric bound: +// the poles are the ends of the axis, so no notation makes 100 meaningful. v2 +// checked it and v1 did not. +func validateLat(lat float64) error { + if math.IsNaN(lat) || lat < -90 || lat > 90 { + return apiError(http.StatusBadRequest, fmt.Sprintf("launch latitude must be in [-90, 90], got %g", lat)) } - return lng + return nil } -// signedLng converts an internal [0, 360) longitude back to [-180, 180). +// validateLng refuses longitudes outside [-longitudeLimit, longitudeLimit]. +// +// Shared by both API versions. v1 used to apply no check at all while v2 enforced +// [-180, 360), so the same launch could be accepted by one and refused by the +// other — a difference in behaviour where only format may differ. +func validateLng(lng float64) error { + if math.IsNaN(lng) || lng < -longitudeLimit || lng > longitudeLimit { + return apiError(http.StatusBadRequest, + fmt.Sprintf("launch longitude must be in [-%g, %g], got %g", longitudeLimit, longitudeLimit, lng)) + } + return nil +} + +// normalizeLng folds a longitude into the [0, 360) the wind grid is indexed on. +// +// A true modulo, via the same helper the integrator applies on every step +// (numerics.PyMod in toGeo), so the request path and the engine agree by +// construction rather than by coincidence. It used to be a single +// `if lng < 0 { lng += 360 }`, correct only on [-360, 360): -400 became -40, +// reached the grid, and came back as "lng=-40 out of range" — a complaint about +// wind data for what was an input problem. +func normalizeLng(lng float64) float64 { + return numerics.PyMod(lng, 360) +} + +// signedLng converts an internal [0, 360) longitude to [-180, 180). +// +// This is v2's published format. v1 publishes the internal [0, 360) unchanged, +// because that is the convention upstream Tawhiri publishes and v1 exists to be a +// drop-in for it. func signedLng(lng float64) float64 { if lng > 180 { return lng - 360 diff --git a/internal/api/mapping_test.go b/internal/api/mapping_test.go new file mode 100644 index 0000000..37d5265 --- /dev/null +++ b/internal/api/mapping_test.go @@ -0,0 +1,88 @@ +package api + +import ( + "math" + "testing" +) + +// Longitude handling is shared by both API versions on purpose: the two may +// differ in format, names and features, never in what the maths does. v1 used to +// apply no range check at all while v2 enforced [-180, 360), so the same launch +// could be accepted by one and refused by the other. + +func TestValidateLngAcceptsBothConventionsAndOneTurnEitherWay(t *testing.T) { + // Upstream Tawhiri implements a half-open [0, 360) — measured: 0, 180, 270 and + // 308.3 accepted, everything negative and 360 refused. We accept a superset + // because the frontend holds longitudes in [-180, 180], so a western launch + // arrives negative. The bound is symmetric so there is nothing to explain + // about one end that does not hold for the other. + for _, lng := range []float64{0, 0.1, 180, 270, 308.3, 359.999, 360, -0.0001, -51.7, -180, -200, -360} { + if err := validateLng(lng); err != nil { + t.Errorf("validateLng(%g) = %v, want accepted", lng, err) + } + } +} + +func TestValidateLngRefusesTypos(t *testing.T) { + // Arithmetically every one of these names a real meridian — 5170 is 50 E. They + // are refused because they are far likelier to be a unit slip or a swapped + // field than an intended notation, and upstream refuses them too. + for _, lng := range []float64{360.0001, -360.0001, 400, -400, 5170, math.Inf(1), math.Inf(-1)} { + if err := validateLng(lng); err == nil { + t.Errorf("validateLng(%g) was accepted, want refused", lng) + } + } +} + +func TestNormalizeLngFoldsAnyInputIntoTheGridDomain(t *testing.T) { + // The wind grid is indexed from LonStart = 0 over 360 degrees, so every value + // handed to it must land in [0, 360). + // + // This was a single `if lng < 0 { lng += 360 }`, correct only on [-360, 360): + // -400 became -40, reached the grid, and came back as "lng=-40 out of range" — + // a complaint about wind data for what was an input problem. + cases := map[float64]float64{ + 0: 0, 0.5: 0.5, 180: 180, 359.5: 359.5, + 360: 0, 720: 0, -360: 0, + -0.5: 359.5, -51.7: 308.3, -180: 180, -200: 160, -400: 320, 5170: 130, + } + for in, want := range cases { + if got := normalizeLng(in); math.Abs(got-want) > 1e-9 { + t.Errorf("normalizeLng(%g) = %g, want %g", in, got, want) + } + if got := normalizeLng(in); got < 0 || got >= 360 { + t.Errorf("normalizeLng(%g) = %g, outside [0, 360)", in, got) + } + } +} + +func TestSignedLngIsTheV2Format(t *testing.T) { + // v2 publishes [-180, 180): its own convenient format. v1 publishes the + // internal [0, 360) unchanged, because that is what upstream Tawhiri publishes + // and v1 exists to be a drop-in for it. + cases := map[float64]float64{0: 0, 90: 90, 180: 180, 180.5: -179.5, 308.3: -51.7, 359.5: -0.5} + for in, want := range cases { + if got := signedLng(in); math.Abs(got-want) > 1e-9 { + t.Errorf("signedLng(%g) = %g, want %g", in, got, want) + } + } +} + +func TestValidateLatIsSharedByBothVersions(t *testing.T) { + // Unlike longitude, latitude is a genuine geometric bound: +-90 are the poles + // and there is nothing beyond them, so no notation makes 100 meaningful. + // + // v2 checked this and v1 did not, which is a difference in behaviour rather + // than format. On v1 an out-of-range latitude reached the wind grid and came + // back as "lat=... out of range" — a complaint about data for an input problem. + for _, lat := range []float64{-90, -89.999, 0, 52.2, 89.999, 90} { + if err := validateLat(lat); err != nil { + t.Errorf("validateLat(%g) = %v, want accepted", lat, err) + } + } + for _, lat := range []float64{-90.0001, 90.0001, 100, -100, math.NaN(), math.Inf(1)} { + if err := validateLat(lat); err == nil { + t.Errorf("validateLat(%g) was accepted, want refused", lat) + } + } +} diff --git a/internal/api/prediction.go b/internal/api/prediction.go index 5af501d..2399270 100644 --- a/internal/api/prediction.go +++ b/internal/api/prediction.go @@ -2,6 +2,7 @@ package api import ( "context" + "fmt" "net/http" "time" @@ -23,9 +24,45 @@ func (h *Handler) ReadinessCheck(_ context.Context) (*apirest.ReadinessResponse, return resp, nil } +// fieldFor resolves the dataset a request asked for, or the active one. +// +// A request naming a dataset that is not stored is refused. The parameter used +// to be declared in the spec and never read, so a client asking for a specific +// forecast run — an archival date, a reproducible comparison — was silently +// served from whatever happened to be loaded instead. +func (h *Handler) fieldFor(ctx context.Context, want apirest.OptDateTime) (weather.WindField, error) { + if epoch, ok := want.Get(); ok { + field, err := h.mgr.FieldFor(ctx, epoch.UTC()) + if err != nil { + return nil, apiError(http.StatusBadRequest, err.Error()) + } + return field, nil + } + field := h.mgr.Active() + if field == nil { + return nil, apiError(http.StatusServiceUnavailable, "no dataset loaded, service is starting up") + } + return field, nil +} + +// runFailure turns an integration failure into a 400. +// +// The sampler reporting an error means it had no data — most often a launch time +// past the dataset's horizon — so the trajectory is not a forecast. It used to be +// discarded, and the run continued with zero wind: the balloon took off and +// landed on the same spot, returned as a successful prediction. +func runFailure(events *engine.EventSink, field weather.WindField) error { + err := events.Err() + if err == nil { + return nil + } + return apiError(http.StatusBadRequest, + fmt.Sprintf("%s (dataset %s)", err, field.Epoch().UTC().Format(time.RFC3339))) +} + // PerformPredictionV2 implements POST /api/v2/prediction. -func (h *Handler) PerformPredictionV2(_ context.Context, req *apirest.PredictionV2Request) (*apirest.PredictionV2Response, error) { - resp, err := h.runPredictionV2(req) +func (h *Handler) PerformPredictionV2(ctx context.Context, req *apirest.PredictionV2Request) (*apirest.PredictionV2Response, error) { + resp, err := h.runPredictionV2(ctx, req) if err == nil { h.metrics.Prediction("v2", resp.CompletedAt.Sub(resp.StartedAt), nil) } @@ -60,23 +97,22 @@ func (h *Handler) CancelPredictionJob(_ context.Context, params apirest.CancelPr // runPredictionV2 is the synchronous prediction core, shared by the v2 // endpoint and the async worker pool. -func (h *Handler) runPredictionV2(req *apirest.PredictionV2Request) (*apirest.PredictionV2Response, error) { +func (h *Handler) runPredictionV2(ctx context.Context, req *apirest.PredictionV2Request) (*apirest.PredictionV2Response, error) { // Validate the request shape before checking dataset availability, so a // malformed request is a 400 regardless of startup state. lat := req.Launch.Latitude - rawLng := req.Launch.Longitude alt := req.Launch.Altitude.Or(0) - if lat < -90 || lat > 90 { - return nil, apiError(http.StatusBadRequest, "launch.latitude must be in [-90, 90]") + if err := validateLat(lat); err != nil { + return nil, err } - if rawLng < -180 || rawLng >= 360 { - return nil, apiError(http.StatusBadRequest, "launch.longitude must be in [-180, 360)") + if err := validateLng(req.Launch.Longitude); err != nil { + return nil, err } - lng := normalizeLng(rawLng) + lng := normalizeLng(req.Launch.Longitude) - field := h.mgr.Active() - if field == nil { - return nil, apiError(http.StatusServiceUnavailable, "no dataset loaded, service is starting up") + field, err := h.fieldFor(ctx, req.Dataset) + if err != nil { + return nil, err } events := engine.NewEventSink() @@ -90,6 +126,9 @@ func (h *Handler) runPredictionV2(req *apirest.PredictionV2Request) (*apirest.Pr started := time.Now().UTC() results := prof.Run(float64(req.Launch.Time.Unix()), engine.State{Lat: lat, Lng: lng, Altitude: alt}, events) completed := time.Now().UTC() + if err := runFailure(events, field); err != nil { + return nil, err + } resp := &apirest.PredictionV2Response{ Stages: make([]apirest.StageResult, 0, len(results)), @@ -105,10 +144,10 @@ func (h *Handler) runPredictionV2(req *apirest.PredictionV2Request) (*apirest.Pr } // PerformPrediction implements GET /api/v1/prediction (Tawhiri-compatible). -func (h *Handler) PerformPrediction(_ context.Context, params apirest.PerformPredictionParams) (*apirest.PredictionResponse, error) { - field := h.mgr.Active() - if field == nil { - return nil, apiError(http.StatusServiceUnavailable, "no dataset loaded, service is starting up") +func (h *Handler) PerformPrediction(ctx context.Context, params apirest.PerformPredictionParams) (*apirest.PredictionResponse, error) { + field, err := h.fieldFor(ctx, params.Dataset) + if err != nil { + return nil, err } profileKind := "standard_profile" @@ -118,6 +157,12 @@ func (h *Handler) PerformPrediction(_ context.Context, params apirest.PerformPre ascentRate := params.AscentRate.Or(5) descentRate := params.DescentRate.Or(5) launchAlt := params.LaunchAltitude.Or(0) + if err := validateLat(params.LaunchLatitude); err != nil { + return nil, err + } + if err := validateLng(params.LaunchLongitude); err != nil { + return nil, err + } lng := normalizeLng(params.LaunchLongitude) launchTime := float64(params.LaunchDatetime.Unix()) @@ -142,6 +187,9 @@ func (h *Handler) PerformPrediction(_ context.Context, params apirest.PerformPre started := time.Now().UTC() results := prof.Run(launchTime, engine.State{Lat: params.LaunchLatitude, Lng: lng, Altitude: launchAlt}, events) completed := time.Now().UTC() + if err := runFailure(events, field); err != nil { + return nil, err + } h.metrics.Prediction(profileKind, completed.Sub(started), nil) resp := &apirest.PredictionResponse{ @@ -229,9 +277,13 @@ func tawhiriItem(name string, r engine.Result) apirest.PredictionResponsePredict for i := range n { t, p := r.Path.At(i) traj = append(traj, apirest.TawhiriPoint{ - Datetime: time.Unix(int64(t), 0).UTC(), - Latitude: p.Lat, - Longitude: signedLng(p.Lng), + Datetime: time.Unix(int64(t), 0).UTC(), + Latitude: p.Lat, + // v1 publishes [0, 360), the convention upstream Tawhiri publishes — + // measured: it answers 308.7884 where we used to answer -51.2115 for + // the same meridian. v1 exists to be a drop-in, so it matches. v2 keeps + // the signed format (mapping.go). + Longitude: p.Lng, Altitude: p.Altitude, }) } diff --git a/internal/api/transport.go b/internal/api/transport.go index 75a9122..fe23987 100644 --- a/internal/api/transport.go +++ b/internal/api/transport.go @@ -74,7 +74,11 @@ func New(port int, d Deps) (*Server, error) { Workers: d.AsyncWorkers, QueueSize: d.AsyncQueueSize, ResultTTL: d.AsyncResultTTL, - }, h.runPredictionV2, d.Metrics, d.Log) + // The async queue outlives the HTTP request that enqueued the job, so a + // worker cannot inherit its context. + }, func(req *apirest.PredictionV2Request) (*apirest.PredictionV2Response, error) { + return h.runPredictionV2(context.Background(), req) + }, d.Metrics, d.Log) ogenSrv, err := apirest.NewServer(h, apirest.WithMiddleware(middleware.OgenLogging(d.Log))) if err != nil { diff --git a/internal/datasets/manager.go b/internal/datasets/manager.go index b56a265..2881545 100644 --- a/internal/datasets/manager.go +++ b/internal/datasets/manager.go @@ -161,6 +161,57 @@ func (m *Manager) SelectFor(t time.Time, lat, lng float64) weather.WindField { return nil } +// FieldFor returns the global dataset whose epoch is epoch, loading it from +// storage on demand when it is stored but not currently active. +// +// Being stored is enough. A caller asking for a specific run — an archival +// forecast, a reproducible comparison — must not have to make it the service's +// active dataset first, and asking for it must not repoint the service: Load +// appends, and Active keeps returning the first global it finds. +// +// Errors when no stored dataset carries that epoch, so a request naming a +// dataset the service cannot serve is refused rather than quietly answered from +// a different one. +// +// ponytail: loaded datasets are never evicted, so a long-lived service asked for +// many distinct epochs accumulates one mmap and fd each. Add LRU eviction here +// if that becomes real; the files are mmap-backed, so the cost is address space +// and descriptors, not resident memory. +func (m *Manager) FieldFor(ctx context.Context, epoch time.Time) (weather.WindField, error) { + if f := m.activeFieldAt(epoch); f != nil { + return f, nil + } + + stored, err := m.store.List() + if err != nil { + return nil, fmt.Errorf("list stored datasets: %w", err) + } + for _, id := range stored { + if !id.Subset.IsGlobal() || !id.Epoch.Equal(epoch) { + continue + } + if err := m.Load(ctx, id); err != nil { + return nil, fmt.Errorf("load %s: %w", id.Filename(), err) + } + if f := m.activeFieldAt(epoch); f != nil { + return f, nil + } + } + return nil, fmt.Errorf("dataset %s is not stored", epoch.UTC().Format(time.RFC3339)) +} + +// activeFieldAt returns the loaded global field with exactly this epoch, or nil. +func (m *Manager) activeFieldAt(epoch time.Time) weather.WindField { + m.activeMu.RLock() + defer m.activeMu.RUnlock() + for _, d := range m.active { + if d.ID.Subset.IsGlobal() && d.ID.Epoch.Equal(epoch) { + return d.Field + } + } + return nil +} + // LoadedDatasets returns snapshots of every currently-loaded dataset. func (m *Manager) LoadedDatasets() []LoadedDatasetInfo { m.activeMu.RLock() @@ -309,36 +360,48 @@ func (m *Manager) Load(ctx context.Context, id DatasetID) error { // // Returns the JobID started, or empty string when nothing was scheduled. func (m *Manager) Refresh(ctx context.Context, freshnessTTL time.Duration) (string, error) { - if a := m.activeGlobal(); a != nil && time.Since(a.ID.Epoch) < freshnessTTL { - return "", nil - } - - if datasets, err := m.store.List(); err == nil { - for _, id := range datasets { - if !id.Subset.IsGlobal() { - continue - } - if time.Since(id.Epoch) > freshnessTTL { - continue - } - if a := m.activeGlobal(); a != nil && a.ID.Equals(id) { - return "", nil - } - if err := m.Load(ctx, id); err == nil { - return "", nil + // Get something usable loaded first, whatever its age. + // + // freshnessTTL answers one question only: go fetch something newer? It must + // never decide whether data already on disk may be read. While one check + // served both, a stored run older than the TTL was skipped even when it was + // the only dataset present — active stayed empty, every prediction answered + // "no dataset loaded" with gigabytes of usable wind on disk, and with no + // reachable origin that state was permanent. + if m.activeGlobal() == nil { + if stored, err := m.store.List(); err == nil { + for _, id := range stored { // newest first + if !id.Subset.IsGlobal() { + continue + } + if err := m.Load(ctx, id); err == nil { + break + } } } } + active := m.activeGlobal() + if active != nil && time.Since(active.ID.Epoch) < freshnessTTL { + return "", nil + } + latest, err := m.src.LatestEpoch(ctx) if err != nil { return "", fmt.Errorf("latest epoch: %w", err) } - id := DatasetID{Epoch: latest} - if a := m.activeGlobal(); a != nil && !latest.After(a.ID.Epoch) { + if active != nil && !latest.After(active.ID.Epoch) { return "", nil } + // Another replica sharing this volume may have committed it already. + id := DatasetID{Epoch: latest} + if m.store.Exists(id) { + if err := m.Load(ctx, id); err == nil { + return "", nil + } + } + jobID := m.Download(id) go m.loadAfterCompletion(jobID, id) return jobID, nil diff --git a/internal/datasets/manager_test.go b/internal/datasets/manager_test.go new file mode 100644 index 0000000..c68f724 --- /dev/null +++ b/internal/datasets/manager_test.go @@ -0,0 +1,166 @@ +package datasets + +import ( + "context" + "errors" + "testing" + "time" + + "predictor-refactored/internal/weather" +) + +// stubField is a WindField that carries only its epoch. +type stubField struct{ epoch time.Time } + +func (f stubField) Wind(_, _, _, _ float64) (weather.Sample, error) { + return weather.Sample{}, nil +} +func (f stubField) Epoch() time.Time { return f.epoch } +func (f stubField) Source() string { return "fake" } + +// fakeSource records what it was asked to open and can simulate an unreachable +// origin via latestErr. +type fakeSource struct { + latest time.Time + latestErr error + opened []DatasetID +} + +func (s *fakeSource) ID() string { return "fake" } + +func (s *fakeSource) LatestEpoch(context.Context) (time.Time, error) { + if s.latestErr != nil { + return time.Time{}, s.latestErr + } + return s.latest, nil +} + +func (s *fakeSource) Download(context.Context, DatasetID, Storage, ProgressSink, Throttle) error { + return errors.New("fakeSource does not download") +} + +func (s *fakeSource) Open(_ context.Context, id DatasetID, _ Storage) (weather.WindField, error) { + s.opened = append(s.opened, id) + return stubField{epoch: id.Epoch}, nil +} + +func (s *fakeSource) Coverage(id DatasetID) Coverage { + return Coverage{ + Region: Region{MinLat: -90, MaxLat: 90, MinLng: 0, MaxLng: 360}, + StartTime: id.Epoch, + EndTime: id.Epoch.Add(192 * time.Hour), + } +} + +// fakeStore is an in-memory Storage holding a fixed set of committed datasets, +// newest first, as LocalStore.List promises. +type fakeStore struct{ ids []DatasetID } + +func (s *fakeStore) SourceID() string { return "fake" } +func (s *fakeStore) Path(id DatasetID) string { return "/dev/null/" + id.Filename() } + +func (s *fakeStore) Exists(id DatasetID) bool { + for _, have := range s.ids { + if have.Equals(id) { + return true + } + } + return false +} + +func (s *fakeStore) List() ([]DatasetID, error) { return s.ids, nil } +func (s *fakeStore) Remove(DatasetID) error { return nil } +func (s *fakeStore) BeginWrite(DatasetID) (TempHandle, error) { return nil, errors.New("not used") } +func (s *fakeStore) Lock(context.Context) (func(), error) { return func() {}, nil } + +// TestRefreshLoadsAStoredDatasetOlderThanTheFreshnessTTL pins the distinction +// the freshness TTL is allowed to make. +// +// A dataset on disk is usable wind data whatever its age. The TTL answers "go +// fetch something newer?", not "may I read what is already here?". While both +// questions shared one check, an archival run that was the only dataset present +// got skipped, active stayed empty, and every prediction answered "no dataset +// loaded" with gigabytes of usable data on disk — and, with no reachable origin, +// forever. +func TestRefreshLoadsAStoredDatasetOlderThanTheFreshnessTTL(t *testing.T) { + old := time.Now().UTC().Add(-30 * 24 * time.Hour).Truncate(time.Hour) + src := &fakeSource{latestErr: errors.New("origin unreachable")} + store := &fakeStore{ids: []DatasetID{{Epoch: old}}} + m := New(src, store, nil, nil) + + // The origin probe is expected to fail; loading what is on disk is not. + _, _ = m.Refresh(context.Background(), 48*time.Hour) + + if m.Active() == nil { + t.Fatal("Active() is nil: the stored dataset was never loaded") + } + if got := m.Active().Epoch(); !got.Equal(old) { + t.Errorf("loaded epoch = %s, want %s", got, old) + } +} + +// TestRefreshPrefersTheNewestStoredDataset guards the ordering the fix relies +// on: falling back to disk must not mean falling back to the oldest file there. +func TestRefreshPrefersTheNewestStoredDataset(t *testing.T) { + newer := time.Now().UTC().Add(-10 * 24 * time.Hour).Truncate(time.Hour) + older := newer.Add(-20 * 24 * time.Hour) + src := &fakeSource{latestErr: errors.New("origin unreachable")} + store := &fakeStore{ids: []DatasetID{{Epoch: newer}, {Epoch: older}}} + m := New(src, store, nil, nil) + + _, _ = m.Refresh(context.Background(), 48*time.Hour) + + if m.Active() == nil { + t.Fatal("Active() is nil") + } + if got := m.Active().Epoch(); !got.Equal(newer) { + t.Errorf("loaded epoch = %s, want the newest stored %s", got, newer) + } +} + +// TestFieldForLoadsAStoredDatasetThatIsNotActive covers asking for a specific +// run by epoch — the archival-forecast case. Being stored is enough; the caller +// should not have to make it the active dataset first. +func TestFieldForLoadsAStoredDatasetThatIsNotActive(t *testing.T) { + active := time.Now().UTC().Truncate(time.Hour) + archive := active.Add(-30 * 24 * time.Hour) + src := &fakeSource{latest: active} + store := &fakeStore{ids: []DatasetID{{Epoch: active}, {Epoch: archive}}} + m := New(src, store, nil, nil) + + if err := m.Load(context.Background(), DatasetID{Epoch: active}); err != nil { + t.Fatalf("Load active: %v", err) + } + + field, err := m.FieldFor(context.Background(), archive) + if err != nil { + t.Fatalf("FieldFor(%s): %v", archive, err) + } + if got := field.Epoch(); !got.Equal(archive) { + t.Errorf("epoch = %s, want %s", got, archive) + } + // The active dataset must be unaffected: asking for an archival run for one + // prediction should not repoint the service at it. + if got := m.Active().Epoch(); !got.Equal(active) { + t.Errorf("Active() moved to %s, want %s", got, active) + } +} + +// TestFieldForRejectsAnEpochThatIsNotStored is the half that matters most: a +// requested dataset that cannot be served must say so rather than quietly +// substituting a different one, which is what ignoring the parameter did. +func TestFieldForRejectsAnEpochThatIsNotStored(t *testing.T) { + active := time.Now().UTC().Truncate(time.Hour) + src := &fakeSource{latest: active} + store := &fakeStore{ids: []DatasetID{{Epoch: active}}} + m := New(src, store, nil, nil) + if err := m.Load(context.Background(), DatasetID{Epoch: active}); err != nil { + t.Fatalf("Load: %v", err) + } + + missing := active.Add(-365 * 24 * time.Hour) + field, err := m.FieldFor(context.Background(), missing) + if err == nil { + t.Fatalf("FieldFor(%s) returned field with epoch %s, want an error", missing, field.Epoch()) + } +} diff --git a/internal/engine/engine_test.go b/internal/engine/engine_test.go index ccf0f88..9344f7d 100644 --- a/internal/engine/engine_test.go +++ b/internal/engine/engine_test.go @@ -1,7 +1,9 @@ package engine import ( + "fmt" "math" + "strings" "testing" "time" @@ -261,3 +263,157 @@ func TestPolygonOutsideAntimeridian(t *testing.T) { t.Errorf("(0, 0) should be outside") } } + +// horizonWind has data up to horizon seconds and fails past it, the way a real +// dataset's time axis does (numerics.AxisError from Axis.Locate). +type horizonWind struct{ horizon float64 } + +func (w horizonWind) Wind(t, _, _, _ float64) (weather.Sample, error) { + if t > w.horizon { + return weather.Sample{}, fmt.Errorf("hour=%v out of range", t/3600) + } + return weather.Sample{U: 10, V: 0}, nil +} +func (w horizonWind) Epoch() time.Time { return time.Unix(0, 0).UTC() } +func (w horizonWind) Source() string { return "test" } + +// TestWindSamplingFailureIsRecordedNotSwallowed pins the difference between an +// event and a failure. +// +// WindTransport used to discard the sampler's error and return a zero rate. A +// zero rate is indistinguishable from calm air, so a launch past the dataset's +// horizon integrated with no wind at all and produced a balloon that took off +// and landed on the same spot — returned as a successful forecast. An error from +// the sampler means there was no data, so the trajectory past that point is not +// a forecast and must not be presented as one. +func TestWindSamplingFailureIsRecordedNotSwallowed(t *testing.T) { + sink := NewEventSink() + ascend := &Propagator{ + Name: "ascent", + Step: 60, + Model: Sum(ConstantRate(5), WindTransport(horizonWind{horizon: 600}, sink)), + Constraints: []Constraint{Altitude{Op: OpGreaterEqual, Limit: 30000, On: ActionStop}}, + } + prof := Profile{Stages: []*Propagator{ascend}, Direction: Forward} + prof.Run(0, State{Lat: 0, Lng: 0, Altitude: 0}, sink) + + err := sink.Err() + if err == nil { + t.Fatal("sink reports no failure; the sampling error was swallowed") + } + if !strings.Contains(err.Error(), "out of range") { + t.Errorf("error %q does not carry the sampler's reason", err) + } +} + +// TestWindSamplingSuccessLeavesNoFailure keeps the check above from firing on +// healthy runs, which would turn every prediction into a 400. +func TestWindSamplingSuccessLeavesNoFailure(t *testing.T) { + sink := NewEventSink() + ascend := &Propagator{ + Name: "ascent", + Step: 60, + Model: Sum(ConstantRate(5), WindTransport(horizonWind{horizon: 1e9}, sink)), + Constraints: []Constraint{Altitude{Op: OpGreaterEqual, Limit: 30000, On: ActionStop}}, + } + prof := Profile{Stages: []*Propagator{ascend}, Direction: Forward} + prof.Run(0, State{Lat: 0, Lng: 0, Altitude: 0}, sink) + + if err := sink.Err(); err != nil { + t.Fatalf("healthy run reported a failure: %v", err) + } +} + +// TestAboveModelStaysNonFatal separates the two paths explicitly: extrapolating +// above the highest pressure level is a warning the caller may ignore, and must +// not become a hard failure. +func TestAboveModelStaysNonFatal(t *testing.T) { + sink := NewEventSink() + ascend := &Propagator{ + Name: "ascent", + Step: 60, + Model: Sum(ConstantRate(5), WindTransport(aboveModelWind{}, sink)), + Constraints: []Constraint{Altitude{Op: OpGreaterEqual, Limit: 30000, On: ActionStop}}, + } + prof := Profile{Stages: []*Propagator{ascend}, Direction: Forward} + prof.Run(0, State{Lat: 0, Lng: 0, Altitude: 0}, sink) + + if err := sink.Err(); err != nil { + t.Fatalf("above_model became a failure: %v", err) + } + if len(sink.Snapshot()) == 0 { + t.Error("above_model event was not emitted") + } +} + +// TestEveryModelHonoursIncludeWind pins the promise ModelSpec.IncludeWind makes. +// +// buildConstantRate and buildParachuteDescent both took BuildDeps as `_` and +// returned a bare vertical model, so include_wind was accepted and discarded. +// Only buildPiecewise ever called maybeAddWind. POST /api/v2/prediction with the +// spec's own documented example therefore answered 200 with a trajectory that +// rose and fell on the spot — a wind-free forecast presented as a forecast. +// +// The loop over modelFactories is the part that matters going forward: a newly +// registered model type cannot be added without deciding here what include_wind +// does to it, which is exactly the step that was skipped before. +func TestEveryModelHonoursIncludeWind(t *testing.T) { + cases := map[string]struct { + spec ModelSpec + wantError bool + }{ + "constant_rate": {spec: ModelSpec{Type: "constant_rate", Rate: 5, IncludeWind: true}}, + "parachute_descent": {spec: ModelSpec{Type: "parachute_descent", SeaLevelRate: 5, IncludeWind: true}}, + "piecewise": {spec: ModelSpec{Type: "piecewise", IncludeWind: true, + Segments: []PiecewiseSegmentSpec{{Until: math.Inf(1), Rate: 5}}}}, + // include_wind on the wind model itself would sum the wind into itself. + // Refused rather than ignored: a silently doubled wind is a wrong forecast + // that looks right. + "wind": {spec: ModelSpec{Type: "wind", IncludeWind: true}, wantError: true}, + } + + for name := range modelFactories { + if _, ok := cases[name]; !ok { + t.Errorf("model %q is registered but has no include_wind case here", name) + } + } + + deps := BuildDeps{Wind: fixedWind{u: 10, v: -4}} + for name, tc := range cases { + t.Run(name, func(t *testing.T) { + built, err := BuildModel(tc.spec, deps) + if tc.wantError { + if err == nil { + t.Fatal("expected include_wind to be refused, got no error") + } + return + } + if err != nil { + t.Fatalf("BuildModel: %v", err) + } + m := built.Model + if built.Build != nil { + m = built.Build(StageContext{}) + } + rate := m(0, State{Lat: 0, Lng: 0, Altitude: 1000}) + if rate.East != 10 || rate.North != -4 { + t.Errorf("include_wind ignored: east/north = %v/%v, want 10/-4", rate.East, rate.North) + } + }) + } +} + +// TestIncludeWindWithoutAFieldIsRefused covers the other half of maybeAddWind, +// which returned the wind-free model when deps.Wind was nil. A request that asks +// for wind and cannot get it must fail, not quietly become a vertical drop. +func TestIncludeWindWithoutAFieldIsRefused(t *testing.T) { + for _, typ := range []string{"constant_rate", "parachute_descent", "piecewise"} { + t.Run(typ, func(t *testing.T) { + spec := ModelSpec{Type: typ, Rate: 5, SeaLevelRate: 5, IncludeWind: true, + Segments: []PiecewiseSegmentSpec{{Until: math.Inf(1), Rate: 5}}} + if _, err := BuildModel(spec, BuildDeps{}); err == nil { + t.Error("include_wind with no wind field was accepted") + } + }) + } +} diff --git a/internal/engine/events.go b/internal/engine/events.go index 7fde684..3032320 100644 --- a/internal/engine/events.go +++ b/internal/engine/events.go @@ -26,10 +26,13 @@ type EventSummary struct { } // EventSink collects events from models and the integrator, aggregating -// duplicate types into a single EventSummary. Safe for concurrent use. +// duplicate types into a single EventSummary. It also carries the run's first +// unrecoverable failure, which is a different thing from an event — see Fail. +// Safe for concurrent use. type EventSink struct { mu sync.Mutex summaries map[string]*EventSummary + err error } // NewEventSink returns an empty sink. @@ -61,6 +64,42 @@ func (s *EventSink) Emit(typ string, t float64, state State, message string) { } } +// Fail records the run's first unrecoverable error. Nil-safe, like Emit. +// +// Deliberately not an event. An event is an observation the caller may choose to +// ignore: "above_model" means samples above the highest pressure level were +// extrapolated, which degrades the answer without invalidating it. A failure +// means the sampler had no data at all, so the trajectory from that point on is +// not a forecast and must not be returned as one. +// +// Routing both through Emit is how the second came to be treated like the first. +// WindTransport discarded the sampler's error and returned a zero rate, which is +// indistinguishable from calm air, so any launch past the dataset's horizon +// integrated with no wind and produced a balloon that took off and landed on the +// same spot — reported as a successful prediction. +// +// Only the first failure is kept: it is the one that explains the run. +func (s *EventSink) Fail(err error) { + if s == nil || err == nil { + return + } + s.mu.Lock() + defer s.mu.Unlock() + if s.err == nil { + s.err = err + } +} + +// Err returns the first failure passed to Fail, or nil if the run was sound. +func (s *EventSink) Err() error { + if s == nil { + return nil + } + s.mu.Lock() + defer s.mu.Unlock() + return s.err +} + // Snapshot returns a stable copy of every summary in deterministic order // (sorted by Type). func (s *EventSink) Snapshot() []EventSummary { diff --git a/internal/engine/models.go b/internal/engine/models.go index 7399961..bb9e72b 100644 --- a/internal/engine/models.go +++ b/internal/engine/models.go @@ -1,7 +1,9 @@ package engine import ( + "fmt" "sort" + "time" "predictor-refactored/internal/numerics" "predictor-refactored/internal/weather" @@ -84,6 +86,11 @@ func WindTransport(field weather.WindField, events *EventSink) Model { return func(t float64, s State) numerics.Rate { sample, err := field.Wind(t, s.Lat, s.Lng, s.Altitude) if err != nil { + // Recorded, not swallowed. A zero rate here is indistinguishable + // from calm air, so returning one turns "we have no data for this + // time" into "the balloon did not move" — see EventSink.Fail. + events.Fail(fmt.Errorf("no wind data at %s: %w", + time.Unix(int64(t), 0).UTC().Format(time.RFC3339), err)) return numerics.Rate{} } if sample.AboveModel && events != nil { diff --git a/internal/engine/registry.go b/internal/engine/registry.go index 27a7a9f..9676ddc 100644 --- a/internal/engine/registry.go +++ b/internal/engine/registry.go @@ -205,21 +205,32 @@ func buildPolygon(spec ConstraintSpec, _ BuildDeps) (Constraint, error) { return NewPolygon(spec.Vertices, mode, act, spec.Label), nil } -func buildConstantRate(spec ModelSpec, _ BuildDeps) (BuiltModel, error) { - return BuiltModel{Model: ConstantRate(spec.Rate)}, nil +func buildConstantRate(spec ModelSpec, deps BuildDeps) (BuiltModel, error) { + if err := checkWind(spec.IncludeWind, deps); err != nil { + return BuiltModel{}, err + } + return BuiltModel{Model: addWind(ConstantRate(spec.Rate), spec.IncludeWind, deps)}, nil } -func buildParachuteDescent(spec ModelSpec, _ BuildDeps) (BuiltModel, error) { +func buildParachuteDescent(spec ModelSpec, deps BuildDeps) (BuiltModel, error) { if spec.SeaLevelRate <= 0 { return BuiltModel{}, fmt.Errorf("parachute_descent requires positive sea_level_rate") } - return BuiltModel{Model: ParachuteDescent(spec.SeaLevelRate)}, nil + if err := checkWind(spec.IncludeWind, deps); err != nil { + return BuiltModel{}, err + } + return BuiltModel{Model: addWind(ParachuteDescent(spec.SeaLevelRate), spec.IncludeWind, deps)}, nil } -func buildWind(_ ModelSpec, deps BuildDeps) (BuiltModel, error) { +func buildWind(spec ModelSpec, deps BuildDeps) (BuiltModel, error) { if deps.Wind == nil { return BuiltModel{}, fmt.Errorf("wind model requires a loaded wind field") } + // Refused rather than ignored: summing the wind into itself would double the + // drift, and a doubled wind is a wrong forecast that looks like a right one. + if spec.IncludeWind { + return BuiltModel{}, fmt.Errorf("wind model does not take include_wind: it would sum the wind into itself") + } return BuiltModel{Model: WindTransport(deps.Wind, deps.Events)}, nil } @@ -231,12 +242,17 @@ func buildPiecewise(spec ModelSpec, deps BuildDeps) (BuiltModel, error) { return BuiltModel{}, fmt.Errorf("piecewise: unknown segment reference %q", s.Reference) } } + // Checked here, not in the closure below: the profile runner calls Build once + // per stage and has nowhere to return an error to. + if err := checkWind(spec.IncludeWind, deps); err != nil { + return BuiltModel{}, err + } // Always build lazily: the profile runner supplies a StageContext before // each stage, which is what resolves absolute / profile-relative / // propagator-relative segment times uniformly. return BuiltModel{ Build: func(ctx StageContext) Model { - return maybeAddWind(Piecewise(resolveSegments(spec.Segments, ctx)), spec.IncludeWind, deps) + return addWind(Piecewise(resolveSegments(spec.Segments, ctx)), spec.IncludeWind, deps) }, }, nil } @@ -266,12 +282,26 @@ func segmentBase(reference string, ctx StageContext) float64 { } } -// maybeAddWind sums a WindTransport model into base when the spec asks for it. -func maybeAddWind(base Model, includeWind bool, deps BuildDeps) Model { - if !includeWind { - return base +// checkWind reports whether include_wind can be satisfied at all. +// +// A request that asks for wind and cannot have it is an error. This used to +// return the wind-free model instead, which turned a missing wind field into a +// balloon that rose and fell on the spot, reported as a successful forecast. +func checkWind(includeWind bool, deps BuildDeps) error { + if includeWind && deps.Wind == nil { + return fmt.Errorf("include_wind requires a loaded wind field") } - if deps.Wind == nil { + return nil +} + +// addWind sums a WindTransport model into base when the spec asks for it. +// +// Callers must have called checkWind first. The split exists because the +// piecewise builder combines inside a closure the profile runner invokes per +// stage, with no way to surface an error — so validation happens at build time +// and combination stays infallible, rather than an error being swallowed there. +func addWind(base Model, includeWind bool, deps BuildDeps) Model { + if !includeWind { return base } return Sum(base, WindTransport(deps.Wind, deps.Events)) diff --git a/internal/windviz/windviz.go b/internal/windviz/windviz.go index bec386a..fcf9117 100644 --- a/internal/windviz/windviz.go +++ b/internal/windviz/windviz.go @@ -68,9 +68,20 @@ const ( // Rasterize samples field over req and returns the U/V grid payload. // // Data is laid out in wind-js scan order: row 0 is the northernmost -// latitude (la1), each row runs west→east, longitudes increasing. Per-cell -// sampling errors (e.g. altitude outside the model) are written as 0 rather -// than failing the whole request; a time outside coverage is a hard error. +// latitude (la1), each row runs west→east, longitudes increasing. +// +// Cells that cannot be sampled individually — a regional dataset queried outside +// its region — are written as 0 and the request still succeeds, because a partly +// covered grid is worth drawing. A grid in which nothing could be sampled is an +// error instead. Time is one value for the whole request, so a time past the +// dataset's horizon fails every cell, and returning that as zeros made a velocity +// layer draw an atmosphere at perfect rest where the honest answer is "we have no +// data for this time". +// +// Known gap: an individual zero-filled cell is indistinguishable from genuine +// calm in the payload. The wind-js format expresses missing data as null, which +// would mean typing Data as []*float64; there is no consumer to validate that +// against since the browser wind layer was removed in the Cesium migration. func Rasterize(field weather.WindField, req Request) (Field, error) { step := req.Step if step <= 0 { @@ -125,6 +136,8 @@ func Rasterize(field weather.WindField, req Request) (Field, error) { v := make([]float64, nx*ny) // Row 0 = north (la1); rows descend in latitude. + var failed int + var firstErr error for j := range ny { lat := maxLat - float64(j)*step for i := range nx { @@ -132,12 +145,20 @@ func Rasterize(field weather.WindField, req Request) (Field, error) { s, err := field.Wind(req.Time, lat, normLng(lng), req.Altitude) idx := j*nx + i if err != nil { - continue // leave as 0 + failed++ + if firstErr == nil { + firstErr = err + } + continue // leave as 0; see the doc comment on this distinction } u[idx] = s.U v[idx] = s.V } } + if failed == nx*ny { + return nil, fmt.Errorf("no wind data anywhere in the requested grid at %s: %w", + time.Unix(int64(req.Time), 0).UTC().Format(time.RFC3339), firstErr) + } refTime := time.Unix(int64(req.Time), 0).UTC().Format("2006-01-02T15:04:05.000Z") mk := func(num int, name string, data []float64) Component { diff --git a/internal/windviz/windviz_test.go b/internal/windviz/windviz_test.go index 521cbdb..aace48e 100644 --- a/internal/windviz/windviz_test.go +++ b/internal/windviz/windviz_test.go @@ -1,6 +1,8 @@ package windviz import ( + "fmt" + "strings" "testing" "time" @@ -94,3 +96,86 @@ func TestCacheRoundTrip(t *testing.T) { t.Errorf("cache should hit after put") } } + +// horizonWind has data up to horizon seconds and fails past it, the way a real +// dataset's time axis does. +type horizonWind struct{ horizon float64 } + +func (w horizonWind) Wind(t, _, _, _ float64) (weather.Sample, error) { + if t > w.horizon { + return weather.Sample{}, fmt.Errorf("hour=%v out of range", t/3600) + } + return weather.Sample{U: 7, V: -3}, nil +} +func (w horizonWind) Epoch() time.Time { return time.Unix(0, 0).UTC() } +func (w horizonWind) Source() string { return "test" } + +// northOnlyWind has data only in the northern hemisphere, the way a regional +// subset does. Its failures are genuinely per-cell. +type northOnlyWind struct{} + +func (northOnlyWind) Wind(_, lat, _, _ float64) (weather.Sample, error) { + if lat < 0 { + return weather.Sample{}, fmt.Errorf("lat=%v outside region", lat) + } + return weather.Sample{U: 5, V: 1}, nil +} +func (northOnlyWind) Epoch() time.Time { return time.Unix(0, 0).UTC() } +func (northOnlyWind) Source() string { return "test" } + +// TestRasterizeRefusesWhenNoCellHasData is the case the package comment already +// promised and the code never implemented ("a time outside coverage is a hard +// error"). +// +// Time is one value for the whole request, so a time past the dataset's horizon +// fails every cell. Each failure was written as a zero, so the response was a +// complete grid of zero wind — a velocity layer draws that as an atmosphere at +// perfect rest. Missing data must not be rendered as calm air. +func TestRasterizeRefusesWhenNoCellHasData(t *testing.T) { + f := horizonWind{horizon: 600} + + out, err := Rasterize(f, Request{Time: 700, MinLng: 0, MaxLng: 360, Step: 90}) + + if err == nil { + t.Fatalf("Rasterize returned %d components instead of an error", len(out)) + } + if !strings.Contains(err.Error(), "out of range") { + t.Errorf("error %q does not carry the sampler's reason", err) + } +} + +// TestRasterizeStillSucceedsWithinCoverage keeps the check above from rejecting +// healthy requests. +func TestRasterizeStillSucceedsWithinCoverage(t *testing.T) { + f := horizonWind{horizon: 600} + + out, err := Rasterize(f, Request{Time: 300, MinLng: 0, MaxLng: 360, Step: 90}) + if err != nil { + t.Fatalf("Rasterize within coverage: %v", err) + } + if got := out[0].Data[0]; got != 7 { + t.Errorf("u = %v, want 7", got) + } +} + +// TestRasterizeZeroFillsIndividualGaps preserves the documented per-cell +// behaviour: a partly-covered grid is still worth drawing, so gaps stay zero and +// the request succeeds. Only a grid with nothing in it is refused. +func TestRasterizeZeroFillsIndividualGaps(t *testing.T) { + out, err := Rasterize(northOnlyWind{}, Request{ + MinLat: -30, MaxLat: 30, MinLng: 0, MaxLng: 30, Step: 30, + }) + if err != nil { + t.Fatalf("Rasterize partly-covered grid: %v", err) + } + + // Row 0 is the northernmost latitude (+30), the last row is -30. + nx := out[0].Header.Nx + ny := out[0].Header.Ny + if got := out[0].Data[0]; got != 5 { + t.Errorf("northern cell u = %v, want 5", got) + } + if got := out[0].Data[(ny-1)*nx]; got != 0 { + t.Errorf("southern gap u = %v, want 0", got) + } +} diff --git a/pkg/rest/oas_client_gen.go b/pkg/rest/oas_client_gen.go index b7dad2f..c74265d 100644 --- a/pkg/rest/oas_client_gen.go +++ b/pkg/rest/oas_client_gen.go @@ -4,6 +4,7 @@ package rest import ( "context" + "io" "net/url" "strings" "time" @@ -239,7 +240,13 @@ func (c *Client) sendCancelDatasetJob(ctx context.Context, params CancelDatasetJ return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeCancelDatasetJobResponse(resp) @@ -331,7 +338,13 @@ func (c *Client) sendCancelPredictionJob(ctx context.Context, params CancelPredi return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeCancelPredictionJobResponse(resp) @@ -408,7 +421,13 @@ func (c *Client) sendCreatePredictionJob(ctx context.Context, request *Predictio return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeCreatePredictionJobResponse(resp) @@ -500,7 +519,13 @@ func (c *Client) sendDeleteDataset(ctx context.Context, params DeleteDatasetPara return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeDeleteDatasetResponse(resp) @@ -592,7 +617,13 @@ func (c *Client) sendGetDatasetJob(ctx context.Context, params GetDatasetJobPara return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeGetDatasetJobResponse(resp) @@ -684,7 +715,13 @@ func (c *Client) sendGetPredictionJob(ctx context.Context, params GetPredictionJ return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeGetPredictionJobResponse(resp) @@ -758,7 +795,13 @@ func (c *Client) sendGetServiceStatus(ctx context.Context) (res *StatusResponse, return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeGetServiceStatusResponse(resp) @@ -955,7 +998,13 @@ func (c *Client) sendGetWindField(ctx context.Context, params GetWindFieldParams return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeGetWindFieldResponse(resp) @@ -1029,7 +1078,13 @@ func (c *Client) sendGetWindMeta(ctx context.Context) (res *WindMeta, err error) return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeGetWindMetaResponse(resp) @@ -1103,7 +1158,13 @@ func (c *Client) sendListDatasetJobs(ctx context.Context) (res []DownloadJob, er return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeListDatasetJobsResponse(resp) @@ -1177,7 +1238,13 @@ func (c *Client) sendListDatasets(ctx context.Context) (res *DatasetList, err er return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeListDatasetsResponse(resp) @@ -1433,7 +1500,13 @@ func (c *Client) sendPerformPrediction(ctx context.Context, params PerformPredic return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodePerformPredictionResponse(resp) @@ -1510,7 +1583,13 @@ func (c *Client) sendPerformPredictionV2(ctx context.Context, request *Predictio return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodePerformPredictionV2Response(resp) @@ -1584,7 +1663,13 @@ func (c *Client) sendReadinessCheck(ctx context.Context) (res *ReadinessResponse return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeReadinessCheckResponse(resp) @@ -1661,7 +1746,13 @@ func (c *Client) sendTriggerDatasetDownload(ctx context.Context, request *Downlo return res, errors.Wrap(err, "do request") } body := resp.Body - defer body.Close() + defer func() { + // Drain the body to EOF before closing, so the underlying + // connection can be reused by the Transport regardless of the + // response status code. See https://github.com/ogen-go/ogen/issues/1670. + _, _ = io.Copy(io.Discard, body) + _ = body.Close() + }() stage = "DecodeResponse" result, err := decodeTriggerDatasetDownloadResponse(resp) diff --git a/pkg/rest/oas_handlers_gen.go b/pkg/rest/oas_handlers_gen.go index 4867c13..9acf62c 100644 --- a/pkg/rest/oas_handlers_gen.go +++ b/pkg/rest/oas_handlers_gen.go @@ -71,7 +71,7 @@ func (s *Server) handleCancelDatasetJobRequest(args [1]string, argsEscaped bool, if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -225,7 +225,7 @@ func (s *Server) handleCancelPredictionJobRequest(args [1]string, argsEscaped bo if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -379,7 +379,7 @@ func (s *Server) handleCreatePredictionJobRequest(args [0]string, argsEscaped bo if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -533,7 +533,7 @@ func (s *Server) handleDeleteDatasetRequest(args [1]string, argsEscaped bool, w if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -687,7 +687,7 @@ func (s *Server) handleGetDatasetJobRequest(args [1]string, argsEscaped bool, w if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -841,7 +841,7 @@ func (s *Server) handleGetPredictionJobRequest(args [1]string, argsEscaped bool, if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -995,7 +995,7 @@ func (s *Server) handleGetServiceStatusRequest(args [0]string, argsEscaped bool, if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -1130,7 +1130,7 @@ func (s *Server) handleGetWindFieldRequest(args [0]string, argsEscaped bool, w h if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -1308,7 +1308,7 @@ func (s *Server) handleGetWindMetaRequest(args [0]string, argsEscaped bool, w ht if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -1443,7 +1443,7 @@ func (s *Server) handleListDatasetJobsRequest(args [0]string, argsEscaped bool, if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -1578,7 +1578,7 @@ func (s *Server) handleListDatasetsRequest(args [0]string, argsEscaped bool, w h if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -1713,7 +1713,7 @@ func (s *Server) handlePerformPredictionRequest(args [0]string, argsEscaped bool if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -1907,7 +1907,7 @@ func (s *Server) handlePerformPredictionV2Request(args [0]string, argsEscaped bo if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -2061,7 +2061,7 @@ func (s *Server) handleReadinessCheckRequest(args [0]string, argsEscaped bool, w if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) @@ -2196,7 +2196,7 @@ func (s *Server) handleTriggerDatasetDownloadRequest(args [0]string, argsEscaped if code != 0 { codeAttr := semconv.HTTPResponseStatusCode(code) attrs = append(attrs, codeAttr) - span.SetAttributes(codeAttr) + span.SetAttributes(attrs...) } attrOpt := metric.WithAttributes(attrs...) diff --git a/pkg/rest/oas_json_gen.go b/pkg/rest/oas_json_gen.go index 21c5d84..bbd3acf 100644 --- a/pkg/rest/oas_json_gen.go +++ b/pkg/rest/oas_json_gen.go @@ -4417,6 +4417,12 @@ func (s *PredictionV2Request) encodeFields(e *jx.Encoder) { e.FieldStart("launch") s.Launch.Encode(e) } + { + if s.Dataset.Set { + e.FieldStart("dataset") + s.Dataset.Encode(e, json.EncodeDateTime) + } + } { if s.Direction.Set { e.FieldStart("direction") @@ -4449,12 +4455,13 @@ func (s *PredictionV2Request) encodeFields(e *jx.Encoder) { } } -var jsonFieldsNameOfPredictionV2Request = [5]string{ +var jsonFieldsNameOfPredictionV2Request = [6]string{ 0: "launch", - 1: "direction", - 2: "profile", - 3: "globals", - 4: "options", + 1: "dataset", + 2: "direction", + 3: "profile", + 4: "globals", + 5: "options", } // Decode decodes PredictionV2Request from json. @@ -4477,6 +4484,16 @@ func (s *PredictionV2Request) Decode(d *jx.Decoder) error { }(); err != nil { return errors.Wrap(err, "decode field \"launch\"") } + case "dataset": + if err := func() error { + s.Dataset.Reset() + if err := s.Dataset.Decode(d, json.DecodeDateTime); err != nil { + return err + } + return nil + }(); err != nil { + return errors.Wrap(err, "decode field \"dataset\"") + } case "direction": if err := func() error { s.Direction.Reset() @@ -4488,7 +4505,7 @@ func (s *PredictionV2Request) Decode(d *jx.Decoder) error { return errors.Wrap(err, "decode field \"direction\"") } case "profile": - requiredBitSet[0] |= 1 << 2 + requiredBitSet[0] |= 1 << 3 if err := func() error { s.Profile = make([]StageSpec, 0) if err := d.Arr(func(d *jx.Decoder) error { @@ -4542,7 +4559,7 @@ func (s *PredictionV2Request) Decode(d *jx.Decoder) error { // Validate required fields. var failures []validate.FieldError for i, mask := range [1]uint8{ - 0b00000101, + 0b00001001, } { if result := (requiredBitSet[i] & mask) ^ mask; result != 0 { // Mask only required fields and check equality to mask using XOR. diff --git a/pkg/rest/oas_response_encoders_gen.go b/pkg/rest/oas_response_encoders_gen.go index 33a690c..341c7fd 100644 --- a/pkg/rest/oas_response_encoders_gen.go +++ b/pkg/rest/oas_response_encoders_gen.go @@ -14,14 +14,12 @@ import ( func encodeCancelDatasetJobResponse(response *CancelDatasetJobNoContent, w http.ResponseWriter, span trace.Span) error { w.WriteHeader(204) - span.SetStatus(codes.Ok, http.StatusText(204)) return nil } func encodeCancelPredictionJobResponse(response *CancelPredictionJobNoContent, w http.ResponseWriter, span trace.Span) error { w.WriteHeader(204) - span.SetStatus(codes.Ok, http.StatusText(204)) return nil } @@ -29,7 +27,6 @@ func encodeCancelPredictionJobResponse(response *CancelPredictionJobNoContent, w func encodeCreatePredictionJobResponse(response *PredictionJob, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(202) - span.SetStatus(codes.Ok, http.StatusText(202)) e := new(jx.Encoder) response.Encode(e) @@ -42,7 +39,6 @@ func encodeCreatePredictionJobResponse(response *PredictionJob, w http.ResponseW func encodeDeleteDatasetResponse(response *DeleteDatasetNoContent, w http.ResponseWriter, span trace.Span) error { w.WriteHeader(204) - span.SetStatus(codes.Ok, http.StatusText(204)) return nil } @@ -50,7 +46,6 @@ func encodeDeleteDatasetResponse(response *DeleteDatasetNoContent, w http.Respon func encodeGetDatasetJobResponse(response *DownloadJob, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) response.Encode(e) @@ -64,7 +59,6 @@ func encodeGetDatasetJobResponse(response *DownloadJob, w http.ResponseWriter, s func encodeGetPredictionJobResponse(response *PredictionJob, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) response.Encode(e) @@ -78,7 +72,6 @@ func encodeGetPredictionJobResponse(response *PredictionJob, w http.ResponseWrit func encodeGetServiceStatusResponse(response *StatusResponse, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) response.Encode(e) @@ -92,7 +85,6 @@ func encodeGetServiceStatusResponse(response *StatusResponse, w http.ResponseWri func encodeGetWindFieldResponse(response []WindComponent, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) e.ArrStart() @@ -110,7 +102,6 @@ func encodeGetWindFieldResponse(response []WindComponent, w http.ResponseWriter, func encodeGetWindMetaResponse(response *WindMeta, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) response.Encode(e) @@ -124,7 +115,6 @@ func encodeGetWindMetaResponse(response *WindMeta, w http.ResponseWriter, span t func encodeListDatasetJobsResponse(response []DownloadJob, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) e.ArrStart() @@ -142,7 +132,6 @@ func encodeListDatasetJobsResponse(response []DownloadJob, w http.ResponseWriter func encodeListDatasetsResponse(response *DatasetList, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) response.Encode(e) @@ -156,7 +145,6 @@ func encodeListDatasetsResponse(response *DatasetList, w http.ResponseWriter, sp func encodePerformPredictionResponse(response *PredictionResponse, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) response.Encode(e) @@ -170,7 +158,6 @@ func encodePerformPredictionResponse(response *PredictionResponse, w http.Respon func encodePerformPredictionV2Response(response *PredictionV2Response, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) response.Encode(e) @@ -184,7 +171,6 @@ func encodePerformPredictionV2Response(response *PredictionV2Response, w http.Re func encodeReadinessCheckResponse(response *ReadinessResponse, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(200) - span.SetStatus(codes.Ok, http.StatusText(200)) e := new(jx.Encoder) response.Encode(e) @@ -198,7 +184,6 @@ func encodeReadinessCheckResponse(response *ReadinessResponse, w http.ResponseWr func encodeTriggerDatasetDownloadResponse(response *DownloadAccepted, w http.ResponseWriter, span trace.Span) error { w.Header().Set("Content-Type", "application/json; charset=utf-8") w.WriteHeader(202) - span.SetStatus(codes.Ok, http.StatusText(202)) e := new(jx.Encoder) response.Encode(e) @@ -217,10 +202,8 @@ func encodeErrorResponse(response *DefaultErrorStatusCode, w http.ResponseWriter code = http.StatusOK } w.WriteHeader(code) - if st := http.StatusText(code); code >= http.StatusBadRequest { - span.SetStatus(codes.Error, st) - } else { - span.SetStatus(codes.Ok, st) + if code >= http.StatusInternalServerError { + span.SetStatus(codes.Error, http.StatusText(code)) } e := new(jx.Encoder) diff --git a/pkg/rest/oas_schemas_gen.go b/pkg/rest/oas_schemas_gen.go index bdcd84a..75fbfba 100644 --- a/pkg/rest/oas_schemas_gen.go +++ b/pkg/rest/oas_schemas_gen.go @@ -2615,13 +2615,15 @@ func (s *PredictionResponseWarnings) init() PredictionResponseWarnings { return m } -// A profile-driven prediction. `profile` is an ordered chain of -// propagators; each integrates from where the previous ended. A stage's -// `constraints` decide when it ends and what happens next: stop the -// profile, hand off to `fallback_index`, or clip to the boundary. +// A profile-driven prediction. `profile` is an ordered chain of propagators; each integrates from +// where the previous ended. A stage's `constraints` decide when it ends and what happens next: stop +// the profile, hand off to `fallback_index`, or clip to the boundary. // Ref: #/components/schemas/PredictionV2Request type PredictionV2Request struct { Launch Launch `json:"launch"` + // Forecast run to predict from, given as its epoch. Defaults to the active dataset. The named run must + // be stored; a request for one that is not is rejected rather than answered from a different run. + Dataset OptDateTime `json:"dataset"` // Forward integrates launch→landing; reverse integrates backward in time. Direction OptPredictionV2RequestDirection `json:"direction"` Profile []StageSpec `json:"profile"` @@ -2635,6 +2637,11 @@ func (s *PredictionV2Request) GetLaunch() Launch { return s.Launch } +// GetDataset returns the value of Dataset. +func (s *PredictionV2Request) GetDataset() OptDateTime { + return s.Dataset +} + // GetDirection returns the value of Direction. func (s *PredictionV2Request) GetDirection() OptPredictionV2RequestDirection { return s.Direction @@ -2660,6 +2667,11 @@ func (s *PredictionV2Request) SetLaunch(val Launch) { s.Launch = val } +// SetDataset sets the value of Dataset. +func (s *PredictionV2Request) SetDataset(val OptDateTime) { + s.Dataset = val +} + // SetDirection sets the value of Direction. func (s *PredictionV2Request) SetDirection(val OptPredictionV2RequestDirection) { s.Direction = val -- 2.47.3