"""Sequential event-by-event unit hydrograph fitting.
Uses :class:`~sparsehydro.calibration.CalibrationProblem` for data preparation
and objective evaluation, and constructs a
:class:`~sparsehydro.calibration.CalibrationResult` per event.
Warm starting from the previous event's fitted parameters is handled by calling
``scipy.optimize.minimize`` directly (bypassing ``ScipySolver``, which always
restarts from the bounds midpoint).
"""
from __future__ import annotations
from dataclasses import dataclass
from typing import Callable
import numpy as np
import pandas as pd
from scipy.optimize import Bounds, minimize as scipy_minimize
from scipy.optimize import differential_evolution
from scipy.signal import savgol_filter
from ...calibration.objectives import RMSE, KGE, NashSutcliffe, IObjective
from ...calibration.problem import CalibrationProblem
from ...calibration.result import CalibrationResult
from ..base import IModel
from ...events import EventRecord
[docs]
@dataclass
class SequentialFitSummary:
"""Results of a sequential unit hydrograph fitting run.
:ivar events: Events that were fitted, in chronological order.
:vartype events: list[EventRecord]
:ivar calibration_results: One
:class:`~sparsehydro.calibration.CalibrationResult` per event.
``pareto_X[0]`` holds fitted parameters;
``objective_display_values()[0]`` returns [RMSE, NSE, KGE].
:vartype calibration_results: list[CalibrationResult]
:ivar global_predicted: Full-domain accumulated predicted flow.
Columns: ``datetime``, ``Q_pred``.
:vartype global_predicted: pandas.DataFrame
:ivar global_observed: Full-domain stormflow.
Columns: ``datetime``, ``stormflow``.
:vartype global_observed: pandas.DataFrame
:ivar model_class_name: Name of the model class produced by ``model_factory``.
:vartype model_class_name: str
:ivar fitted_effective_areas: Fitted effective area per event:
``sum(Q_pred) / sum(rain)`` over the event window. Parallel to
*events* and *calibration_results*.
:vartype fitted_effective_areas: list[float]
"""
events: list[EventRecord]
calibration_results: list[CalibrationResult]
global_predicted: pd.DataFrame
global_observed: pd.DataFrame
model_class_name: str
fitted_effective_areas: list[float]
[docs]
def parameter_evolution(self) -> pd.DataFrame:
"""Return fitted parameters indexed by event.
Columns: ``event_id``, ``start_datetime``, ``end_datetime``,
then one column per calibrated parameter.
:returns: One row per fitted event with its calibrated parameters.
:rtype: pandas.DataFrame
"""
rows = []
for event, result in zip(self.events, self.calibration_results):
row: dict = {
"event_id": event.event_id,
"start_datetime": event.start_datetime,
"end_datetime": event.end_datetime,
}
row.update(dict(zip(result.param_names, result.pareto_X[0])))
rows.append(row)
return pd.DataFrame(rows)
[docs]
def metrics_summary(self) -> pd.DataFrame:
"""Return per-event objective values in display form.
Columns: ``event_id``, then one column per objective name.
:returns: One row per fitted event with display-form objective values.
:rtype: pandas.DataFrame
"""
rows = []
for event, result in zip(self.events, self.calibration_results):
display_F = result.objective_display_values()[0]
row: dict = {"event_id": event.event_id}
for name, val in zip(result.objective_names, display_F):
row[name] = val
rows.append(row)
return pd.DataFrame(rows)
[docs]
def calibrated_areas(self) -> list[float]:
"""Return the calibrated UH area parameter per event.
For a single UH model this is the ``A`` parameter.
For an ensemble (aliases like ``fast_A``, ``slow_A``) it is the sum of all
``*_A`` parameters — consistent with ``mode="sum"`` and fixed weights of 1.0.
:returns: Calibrated total ``A`` value per fitted event.
:rtype: list[float]
:raises ValueError: If no ``A``-bearing parameter is found in a
calibration result.
"""
result: list[float] = []
for cal in self.calibration_results:
params = dict(zip(cal.param_names, cal.pareto_X[0]))
a_params = {k: v for k, v in params.items()
if k == "A" or k.endswith("_A")}
if not a_params:
raise ValueError(
f"No 'A' parameter found among calibrated params: {list(params.keys())}"
)
result.append(float(sum(a_params.values())))
return result
[docs]
def effective_area_summary(self) -> pd.DataFrame:
"""Return observed and fitted effective area per event.
Columns: ``event_id``, ``start_datetime``, ``observed_ae``, ``fitted_ae``.
:returns: One row per fitted event comparing observed vs fitted area.
:rtype: pandas.DataFrame
"""
rows = [
{
"event_id": e.event_id,
"start_datetime": e.start_datetime,
"observed_ae": e.effective_area,
"fitted_ae": ae,
}
for e, ae in zip(self.events, self.fitted_effective_areas)
]
return pd.DataFrame(rows)
[docs]
class SequentialFitter:
"""Fit an :class:`~sparsehydro.models.IModel` to storm events one by one.
:param model_factory: Zero-argument callable returning a fresh
:class:`~sparsehydro.models.IModel` in CREATED state. Works with single
UH models and :class:`~sparsehydro.models.EnsembleModel`::
SequentialFitter(lambda: GammaUH(), data, events)
SequentialFitter(lambda: make_ensemble(), data, events)
:type model_factory: Callable[[], IModel]
:param rain_stormflow: Full time series with columns ``datetime``,
``rain``, ``stormflow``.
:type rain_stormflow: pandas.DataFrame
:param events: Events to fit (sorted internally by ``start_datetime``).
:type events: list[EventRecord]
:param output_column: Column name produced by ``model.predict()``.
Default ``"Q_pred"``.
:type output_column: str
"""
def __init__(
self,
model_factory: Callable[[], IModel],
rain_stormflow: pd.DataFrame,
events: list[EventRecord],
output_column: str = "Q_pred",
) -> None:
self._factory = model_factory
self._data = rain_stormflow.copy()
self._data["datetime"] = pd.to_datetime(self._data["datetime"])
self._data["stormflow"] = np.maximum(self._data["stormflow"].values, 0.0)
self._events = sorted(events, key=lambda e: e.start_datetime)
self._output_column = output_column
[docs]
def fit(
self,
time_range: tuple[str, str] | None = None,
objectives: list[IObjective] | None = None,
calibration_objective_index: int = 0,
method: str = "Nelder-Mead",
smooth_obs: bool = True,
verbose: bool = True,
) -> SequentialFitSummary:
"""Run the sequential fitting loop.
:param time_range: Restrict fitting to this ``(start, end)`` date range.
:type time_range: tuple[str, str] | None
:param objectives: Objectives to evaluate; defaults to
``[RMSE(), NashSutcliffe(), KGE()]``.
:type objectives: list[IObjective] | None
:param calibration_objective_index: Index of *objectives* to minimise
(default 0 = RMSE).
:type calibration_objective_index: int
:param method: ``scipy.optimize.minimize`` method, or
``"differential_evolution"`` for a global search.
:type method: str
:param smooth_obs: Apply a light Savitzky-Golay smooth to observed flow
within each event.
:type smooth_obs: bool
:param verbose: Print per-event progress.
:type verbose: bool
:returns: Summary of the sequential fitting run.
:rtype: SequentialFitSummary
"""
if objectives is None:
objectives = [RMSE(), NashSutcliffe(), KGE()]
data = self._filter_data(time_range)
events = self._filter_events(time_range)
if not events:
return SequentialFitSummary(
events=[],
calibration_results=[],
global_predicted=pd.DataFrame({"datetime": data["datetime"], "Q_pred": np.zeros(len(data))}),
global_observed=data[["datetime", "stormflow"]].copy(),
model_class_name=self._factory().__class__.__name__,
fitted_effective_areas=[],
)
datetime_arr = data["datetime"].values
Q_residual = data["stormflow"].values.copy()
Q_global = np.zeros(len(data))
cal_results: list[CalibrationResult] = []
fitted_events: list[EventRecord] = []
fitted_areas: list[float] = []
last_x: np.ndarray | None = None
last_param_names: list[str] = []
n_events = len(events)
end_dt = pd.Timestamp(datetime_arr[-1])
for i, event in enumerate(events):
if verbose:
print(f"[{i+1}/{n_events}] Event {event.event_id}: {event.start_datetime} → {event.end_datetime}")
fit_mask = (data["datetime"] >= event.start_datetime) & (data["datetime"] <= event.end_datetime)
if not fit_mask.any():
if verbose:
print(" Skipped: no data in event window.")
continue
fit_indices = np.where(fit_mask.values)[0]
idx_s, idx_e = int(fit_indices[0]), int(fit_indices[-1]) + 1
app_end = events[i + 1].start_datetime if i < n_events - 1 else end_dt
app_mask = (data["datetime"] >= event.start_datetime) & (data["datetime"] < app_end)
app_indices = np.where(app_mask.values)[0]
event_data = data[fit_mask].copy().reset_index(drop=True)
obs = Q_residual[idx_s:idx_e].copy()
if smooth_obs and len(obs) >= 5:
wlen = 5
obs = savgol_filter(obs, window_length=wlen, polyorder=2, mode="interp")
obs = np.maximum(obs, 0.0)
event_data["stormflow"] = obs
model_i = self._factory()
model_i.initialize()
if last_x is not None and last_param_names:
for name, val in zip(last_param_names, last_x):
if name in model_i.scalar_parameter_names:
p = model_i.get_scalar_parameter(name)
p.update(value=float(np.clip(val, p.lower_bound, p.upper_bound)))
model_i.validate()
try:
problem = CalibrationProblem(
model=model_i,
data=event_data,
objectives=objectives,
column_map={"observed": "stormflow", "predicted": self._output_column},
)
except Exception as exc:
if verbose:
print(f" CalibrationProblem failed: {exc!r} — skipped.")
continue
xl, xu = problem.bounds
x0 = np.clip(last_x, xl, xu) if (last_x is not None and last_param_names == problem.param_names) else (xl + xu) / 2.0
best_x, success = self._optimise(problem, x0, xl, xu, method, calibration_objective_index)
all_F = problem.evaluate(best_x).reshape(1, -1)
cal_result = CalibrationResult(
history=[],
pareto_X=best_x.reshape(1, -1),
pareto_F=all_F,
param_names=problem.param_names,
objective_names=problem.objective_names,
minimize_flags=problem.minimize_flags,
)
if verbose:
disp = cal_result.objective_display_values()[0]
obj_strs = ", ".join(f"{n}={v:.4f}" for n, v in zip(cal_result.objective_names, disp))
print(f" {obj_strs}, success={success}")
for name, val in zip(problem.param_names, best_x):
if name in model_i.scalar_parameter_names:
model_i.get_scalar_parameter(name).update(value=float(val))
orig_rain = float(data[fit_mask]["rain"].values.sum())
orig_flow = float(np.maximum(data[fit_mask]["stormflow"].values, 0.0).sum())
event.effective_area = orig_flow / orig_rain if orig_rain > 0 else 0.0
app_data = data[app_mask].copy().reset_index(drop=True)
fitted_ae = 0.0
if len(app_data) > 0:
try:
model_i.prepare(app_data)
app_pred = model_i.predict()
q_app = app_pred[self._output_column].values
pred_len = min(len(q_app), len(app_indices))
Q_global[app_indices[:pred_len]] += q_app[:pred_len]
Q_residual[app_indices[:pred_len]] -= q_app[:pred_len]
Q_residual = np.maximum(Q_residual, 0.0)
n_event = len(event_data)
rain_event_sum = float(event_data["rain"].values.sum())
if rain_event_sum > 0 and n_event <= len(q_app):
fitted_ae = float(q_app[:n_event].sum()) / rain_event_sum
except Exception as exc:
if verbose:
print(f" Application window prediction failed: {exc!r}")
last_x = best_x
last_param_names = problem.param_names
cal_results.append(cal_result)
fitted_events.append(event)
fitted_areas.append(fitted_ae)
return SequentialFitSummary(
events=fitted_events,
calibration_results=cal_results,
global_predicted=pd.DataFrame({"datetime": datetime_arr, "Q_pred": Q_global}),
global_observed=data[["datetime", "stormflow"]].copy(),
model_class_name=self._factory().__class__.__name__,
fitted_effective_areas=fitted_areas,
)
def _optimise(self, problem, x0, xl, xu, method, obj_idx):
"""Minimise a single objective of *problem* from start point *x0*.
:param problem: Calibration problem exposing ``evaluate``.
:type problem: CalibrationProblem
:param x0: Initial parameter vector.
:type x0: numpy.ndarray
:param xl: Lower parameter bounds.
:type xl: numpy.ndarray
:param xu: Upper parameter bounds.
:type xu: numpy.ndarray
:param method: ``scipy.optimize.minimize`` method or
``"differential_evolution"``.
:type method: str
:param obj_idx: Index of the objective to minimise.
:type obj_idx: int
:returns: Tuple ``(best_x, success)``.
:rtype: tuple[numpy.ndarray, bool]
"""
bounds_list = list(zip(xl.tolist(), xu.tolist()))
def scalar_fn(x):
return float(problem.evaluate(x)[obj_idx])
if method == "differential_evolution":
res = differential_evolution(scalar_fn, bounds=bounds_list, seed=42, polish=True)
return res.x, bool(res.success)
try:
res = scipy_minimize(
scalar_fn, x0=x0, bounds=Bounds(lb=xl, ub=xu), method=method,
options={"maxiter": 2000, "xatol": 1e-6, "fatol": 1e-8},
)
return res.x, bool(res.success)
except Exception:
return x0.copy(), False
def _filter_data(self, time_range):
"""Return the forcing data optionally restricted to *time_range*.
:param time_range: ``(start, end)`` date range, or ``None`` for all data.
:type time_range: tuple[str, str] | None
:returns: Filtered forcing DataFrame with a reset index.
:rtype: pandas.DataFrame
"""
data = self._data.copy()
if time_range is not None:
s, e = pd.to_datetime(time_range[0]), pd.to_datetime(time_range[1])
data = data[(data["datetime"] >= s) & (data["datetime"] <= e)]
return data.reset_index(drop=True)
def _filter_events(self, time_range):
"""Return the events optionally restricted to *time_range*.
:param time_range: ``(start, end)`` date range, or ``None`` for all events.
:type time_range: tuple[str, str] | None
:returns: Events fully contained within the range.
:rtype: list[EventRecord]
"""
evts = self._events
if time_range is not None:
s, e = pd.to_datetime(time_range[0]), pd.to_datetime(time_range[1])
evts = [ev for ev in evts if ev.start_datetime >= s and ev.end_datetime <= e]
return evts
__all__ = ["SequentialFitter", "SequentialFitSummary"]