Metadata-Version: 2.4
Name: hydromend
Version: 0.1.0
Summary: Learnable lag operators for post-processing hydrodynamic model output against observations.
Project-URL: Homepage, https://github.com/thomasmonahan/hydromend
Project-URL: Documentation, https://github.com/thomasmonahan/hydromend/tree/main/docs
Project-URL: Issues, https://github.com/thomasmonahan/hydromend/issues
Author: Thomas Monahan
License: MIT License
        
        Copyright (c) 2026 Thomas Monahan
        
        Permission is hereby granted, free of charge, to any person obtaining a copy
        of this software and associated documentation files (the "Software"), to deal
        in the Software without restriction, including without limitation the rights
        to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
        copies of the Software, and to permit persons to whom the Software is
        furnished to do so, subject to the following conditions:
        
        The above copyright notice and this permission notice shall be included in all
        copies or substantial portions of the Software.
        
        THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
        IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
        FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
        AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
        LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
        OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
        SOFTWARE.
License-File: LICENSE
Keywords: GESLA,GTSM,Volterra,oceanography,post-processing,regression,storm surge,tide gauge
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Atmospheric Science
Classifier: Topic :: Scientific/Engineering :: Hydrology
Requires-Python: >=3.9
Requires-Dist: matplotlib>=3.5
Requires-Dist: netcdf4>=1.5
Requires-Dist: numpy>=1.22
Requires-Dist: pandas>=1.4
Requires-Dist: pyarrow>=8
Requires-Dist: scikit-learn>=1.1
Requires-Dist: scipy>=1.8
Requires-Dist: seaborn>=0.12
Requires-Dist: xarray>=2022.3
Provides-Extra: dev
Requires-Dist: pytest>=7; extra == 'dev'
Requires-Dist: ruff>=0.1; extra == 'dev'
Requires-Dist: statsmodels>=0.13; extra == 'dev'
Provides-Extra: examples
Requires-Dist: jupyter; extra == 'examples'
Requires-Dist: numba>=0.57; extra == 'examples'
Provides-Extra: gam
Requires-Dist: statsmodels>=0.13; extra == 'gam'
Description-Content-Type: text/markdown

# hydromend

**Learnable lag operators for post-processing hydrodynamic model output against observations.**

`hydromend` turns a co-located pair of a **model** series (e.g. a GTSM / tide–surge
reanalysis point) and an **observed** series (e.g. a tide gauge) into a small,
interpretable *lag operator* that maps a window of recent model values onto the
observation. Operators can be

- **linear** — a lag kernel *w(τ)*, or
- **bilinear** — a second-order Volterra kernel *w(τ₁,τ₂)* that captures tide–surge
  interaction and shallow-water non-linearity.

The default estimator is a **variational-Bayes ARD** regressor (`vb_ard`): it
sparsifies irrelevant lags and returns coefficient uncertainties. Ordinary least
squares, a GAM, and Gaussian-process-prior variants are available as options.

```bash
pip install hydromend            # fit + apply
pip install "hydromend[gam]"     # + statsmodels GAM backend
```

---

## Fit your own operator — train, then test

Give `hydromend` a per-site `DataFrame` indexed by time with an `observations`
column and a `model` column. Build the feature set, hold out a period, fit on the
train part with the default VB-ARD estimator, and evaluate on the unseen test part.

```python
import hydromend as hm
import pandas as pd

# --- feature recipe: 0–24 h memory, full bilinear (Volterra) kernel -----------
recipe = dict(lags_hours=range(1, 25), feature_set="bilinear",
              include_current=True, include_squares=True)

X = hm.build_feature_set(site_df, **recipe)
d = pd.concat([X, site_df["observations"]], axis=1).dropna()

# --- temporal train/test split ------------------------------------------------
split = int(len(d) * 0.7)
Xtr, ytr = d.iloc[:split][X.columns], d.iloc[:split]["observations"]
Xte, yte = d.iloc[split:][X.columns], d.iloc[split:]["observations"]

# --- fit (VB-ARD by default) and evaluate on the held-out period --------------
reg  = hm.make_regressor().fit(Xtr, ytr)              # make_regressor("ols"|"gp_volterra"|... to switch
pred = pd.Series(reg.predict(Xte), index=Xte.index)

mae_raw = (yte - site_df["model"].reindex(Xte.index)).abs().mean()
mae_op  = (yte - pred).abs().mean()
print(f"held-out MAE: {mae_raw:.3f} -> {mae_op:.3f} m")
```

### Reuse the fitted operator anywhere

Wrap the fitted regressor in a portable `Operator` and apply it to any future
model series for that site:

```python
op = hm.Operator.from_regressor(reg, feature_names=list(X.columns), **recipe,
                                metadata={"site": "MySite", "lat": 51.5, "lon": 1.4})
op.predict(new_model_series)          # -> corrected Series
```

### See what it learned

```python
w = hm.extract_linear_lag_weights(reg, X.columns)
hm.plot_lag_weight_kernel(w)                          # linear kernel w(τ), with 95% bands

surf = hm.extract_bilinear_lag_surface(reg, X.columns)
hm.plot_bilinear_kernel_surface(surf)                 # bilinear kernel w(τ1,τ2)
```

Or in the frequency domain — the admittance and the input-weighted quadratic
transfer function:

```python
f, W = hm.weighted_quadratic_transfer(reg, model_series, X.columns)  # |H2|·|X(f1)||X(f2)|
hm.plot_quadratic_transfer(f, W)                       # which frequency pairs interact
```

**Worked examples** ([`examples/`](examples/)):

- [`quickstart.ipynb`](examples/quickstart.ipynb) — end-to-end on **real GESLA +
  GTSM** data (load the pair → fit → inspect → evaluate → save/reload).
- [`synthetic_volterra.ipynb`](examples/synthetic_volterra.ipynb) — **arbitrary
  models on synthetic data**: impose a known Volterra operator, recover it,
  compare linear vs bilinear, see the overtides in the spectrum, and use several
  predictors at once.
- [`1-d_estuary.ipynb`](examples/1-d_estuary.ipynb) — recreate the paper's **1-D
  shallow-water mechanism figures**: a numba estuary solver generates baseline
  vs target series, and `hydromend` learns the linear impulse response, the
  bilinear Volterra kernel *K(τ₁,τ₂)*, and the weighted QTF for each mechanism.

---

## Apply a released operator (ERA5-GTSM)

Once a weight library exists (e.g. the ERA5-GTSM tide-gauge release on Zenodo),
correcting a model series needs no fitting at all:

```python
import hydromend as hm

lib = hm.OperatorLibrary.from_parquet("era5gtsm_operators.parquet")

op = lib.nearest(lat=51.5, lon=1.4)      # or lib.get("Sheerness")
corrected = op.predict(gtsm_series)       # hourly model Series -> corrected Series
```

`op.predict` rebuilds the exact lagged/bilinear features the operator was trained
on and evaluates `X · coef + intercept` (numpy + pandas only).

### The whole GESLA + GTSM compare, in a few lines

```python
from hydromend.datasets import load_pair

# reads the gauge, finds its nearest GTSM point, streams that station, aligns hourly
df, info = load_pair("gesla4/Sheerness.txt", "GTSM_ERA5_E/water_level_*/*.nc")

op = lib.get(info["site"])
df["corrected"] = op.predict(df[["model"]])
print((df["observations"] - df["corrected"]).abs().mean())
```

---

## Benchmark many sites / regressors at once (optional)

To compare feature sets or estimators across a whole network:

```python
groups = hm.load_model_groups({"gtsm": "pairs/*.nc"})           # {model: {site: df}}
results, preds = hm.run_benchmark(
    groups,
    lags_hours=range(1, 25),
    feature_sets=("linear", "bilinear"),
    regressors=("vb_ard", "ols"),        # default is ("vb_ard",); add others to compare
    split_config=hm.SplitConfig(test_size_hours=4 * 8760),
)
hm.summarise_results(results)          # mean improvement ratio by model/regressor/features
```

---

## What's in the box

| Module | Purpose |
|---|---|
| `hydromend.features` | lagged + bilinear (Volterra) feature construction |
| `hydromend.models` | regressors — `vb_ard` (default), `ols`, `gam`, `gp_lag`, `gp_volterra` |
| `hydromend.benchmark` | train/test splitting, single-case fit, multi-site benchmark |
| `hydromend.weights` | extract & plot the learned linear / bilinear kernels |
| `hydromend.spectral` | frequency-domain views: admittance *H(f)*, quadratic transfer *H₂(f₁,f₂)*, weighted QTF |
| `hydromend.data` | load observation/model pairs from NetCDF into per-site frames |
| `hydromend.plotting` | diagnostic time-series / residual / learning-curve plots |
| `hydromend.priors` | GP priors that couple neighbouring lag weights |
| `hydromend.pretrained` | `Operator` / `OperatorLibrary` — save, ship, and apply weights |
| `hydromend.gesla` | GESLA-4.1 tide-gauge reader (numpy/pandas only) |
| `hydromend.gtsm` | stream a single GTSM-ERA5 station out of the monthly NetCDFs |
| `hydromend.datasets` | `load_pair` — one gauge + its GTSM point, aligned hourly |

See [`docs/`](docs/) for the concepts guide and API reference.

## Estimators

`make_regressor()` returns the **default VB-ARD** operator; pass a name or a kwargs
dict to switch:

| name | class | what it is |
|---|---|---|
| **`vb_ard`** (default) | `VBARDRegressor` | variational-Bayes automatic relevance determination — sparsifies lags, gives coefficient uncertainties |
| `ols` | `LinearRegression` | least squares — fast baseline |
| `gam` | `GAMRegressor` | additive splines (needs `hydromend[gam]`) |
| `gp_lag` | `GPLagWeightRegressor` | linear lag model with a GP prior over `w(τ)` |
| `gp_volterra` | `GPVolterraRegressor` | linear **and** bilinear weights, each with a GP prior |

## Caveats

- An operator is only as good as the model↔observation pairing it was fit on;
  apply it at (or very near) the station it was trained for.
- Released `mae_*` fields describe *in-sample* fit quality unless stated
  otherwise — see the accompanying dataset / paper for cross-validated skill.
- `predict` returns values referenced to whatever datum the operator was trained
  against (released ERA5-GTSM operators use the most-recent gauge datum segment).

## Citation

If you use `hydromend` or the ERA5-GTSM operator library, please cite the
accompanying paper and dataset:
```
@article{monahan2026learning,
  title={Learning unresolved coastal dynamics in hydrodynamic models},
  author={Monahan, Thomas Carey and Polton, Jeff and Innocenti, Silvia and Matte, Pascal and Ayyad, Mahmoud and Saman, Krijn and Adcock, Thomas AA},
  year={2026},
  publisher={EarthArXiv}
}
```
