mangrovedigital/tide-engine-api

1# ===== LEGACY / EXPERIMENTAL — NOT SHIPPED =====
2# Part of the superseded engine.core research stack (neighbour-transfer
3# nowcast). NOT imported by the product path (api_server / point_forecast /
4# engine.forecast). Kept for reference only. See LEGACY.md.
5"""Synoptic meteorological correction: the learned, direction-dependent
6wind/pressure surge.
7
8Design (validated against held-out data):
9  * NO intercept  -- an intercept encodes the training period's mean state and
10    injects bias when applied to another period.
11  * High-pass both target and forcing (subtract a centered running mean) to
12    isolate the synoptic weather band from the slow seasonal/runoff band.
13  * Wind enters as eastward/northward velocity (u, v) AND stress (speed*u,
14    speed*v ~ speed^2). The direction dependence is LEARNED, not assumed.
15  * Lags capture the bays' delayed response to wind.
16
17The centered high-pass uses neighbouring hours; this is operationally valid
18because the forcing comes from a met source that provides forecast hours on both
19sides of the target time.
20"""
21import math
22from . import linalg
23
24P_REF = 1013.25
25
26
27def running_mean(x, win):
28    n = len(x); half = win // 2; out = [None] * n
29    for i in range(n):
30        lo = max(0, i - half); hi = min(n, i + half + 1)
31        vals = [x[j] for j in range(lo, hi) if x[j] is not None]
32        if len(vals) >= max(3, win // 3):
33            out[i] = sum(vals) / len(vals)
34    return out
35
36
37def highpass(x, win):
38    if not win:
39        return list(x)
40    rm = running_mean(x, win)
41    return [None if (x[i] is None or rm[i] is None) else x[i] - rm[i]
42            for i in range(len(x))]
43
44
45def _wind_comp(speed, dirf):
46    r = math.radians(dirf)
47    u = -speed * math.sin(r)   # eastward velocity (toward)
48    v = -speed * math.cos(r)   # northward velocity (toward)
49    return u, v, speed * u, speed * v
50
51
52class SynopticModel:
53    def __init__(self, coeffs, lags, hp_win):
54        self.coeffs = coeffs
55        self.lags = lags
56        self.hp_win = hp_win
57
58    def _forcing_arrays(self, speed, dirf, pres):
59        n = len(speed)
60        U = [0.0] * n; V = [0.0] * n; SX = [0.0] * n; SY = [0.0] * n
61        for i in range(n):
62            U[i], V[i], SX[i], SY[i] = _wind_comp(speed[i], dirf[i])
63        dP = highpass(pres, self.hp_win) if self.hp_win else \
64            [pres[i] - P_REF for i in range(n)]
65        if self.hp_win:
66            U = highpass(U, self.hp_win); V = highpass(V, self.hp_win)
67            SX = highpass(SX, self.hp_win); SY = highpass(SY, self.hp_win)
68        return U, V, SX, SY, dP
69
70    def rows(self, speed, dirf, pres):
71        """Feature rows aligned to the hourly arrays; None where a lag underflows
72        or a high-pass value is missing (-> caller treats as HOLD)."""
73        U, V, SX, SY, dP = self._forcing_arrays(speed, dirf, pres)
74        out = []
75        for i in range(len(speed)):
76            row = []; ok = True
77            for L in self.lags:
78                j = i - L
79                if j < 0 or U[j] is None or dP[j] is None:
80                    ok = False; break
81                row += [U[j], V[j], SX[j], SY[j], dP[j]]
82            out.append(row if ok else None)
83        return out
84
85    def predict(self, speed, dirf, pres):
86        rows = self.rows(speed, dirf, pres)
87        return [None if r is None else sum(r[k] * self.coeffs[k]
88                                           for k in range(len(self.coeffs)))
89                for r in rows]
90
91
92def fit(target_resid, speed, dirf, pres, lags=(0, 1, 2, 3, 4, 5, 6),
93        hp_win=120, ridge=1e-6):
94    """Fit on hourly arrays. target_resid = observed - astronomical (the
95    non-tidal residual). It is high-passed to the synoptic band before fitting."""
96    model = SynopticModel(None, list(lags), hp_win)
97    tgt = highpass(target_resid, hp_win) if hp_win else list(target_resid)
98    rows = model.rows(speed, dirf, pres)
99    X, y = [], []
100    for i in range(len(tgt)):
101        if tgt[i] is None or rows[i] is None:
102            continue
103        X.append(rows[i]); y.append(tgt[i])
104    if len(X) < len(rows[0] or [1]) * 3:
105        raise RuntimeError("too few usable samples for synoptic fit")
106    model.coeffs = linalg.solve_lstsq(X, y, ridge=ridge)
107    return model, len(X)
108
109
110def direction_response(model, speed=10.0):
111    """Diagnostic: steady-state surge (m) the model predicts for a sustained wind
112    of `speed` m/s blowing FROM each compass direction. Reveals the LEARNED
113    geometry (which winds fill vs drain the bays)."""
114    out = {}
115    dirs = {"N": 0, "NE": 45, "E": 90, "SE": 135,
116            "S": 180, "SW": 225, "W": 270, "NW": 315}
117    # steady wind -> high-pass of a constant is ~0, so evaluate the LINEAR wind
118    # response directly from the per-lag coefficients (sum over lags), using the
119    # raw (un-high-passed) wind components, which is what a sustained anomaly is.
120    ncoef_per_lag = 5
121    for name, deg in dirs.items():
122        u, v, sx, sy = _wind_comp(speed, deg)
123        s = 0.0
124        for li in range(len(model.lags)):
125            base = li * ncoef_per_lag
126            s += (model.coeffs[base] * u + model.coeffs[base + 1] * v +
127                  model.coeffs[base + 2] * sx + model.coeffs[base + 3] * sy)
128        out[name] = s
129    return out
130