"""Stage 5 - calibration points CSV and HydroBayesCal config emit.
Turns the hydraulic measurements into the calibration-points CSV that
HydroBayesCal reads (``id, x, y, z, <QTY>_DATA, <QTY>_ERROR`` columns) and emits
a ``config_Telemac.py`` wiring the produced case to the calibration parameters
and BAL settings declared in the axqua config.
"""
from __future__ import annotations
import logging
from collections.abc import Callable, Mapping
from pathlib import Path
import numpy as np
import pandas as pd
from axqua.config import Config
from axqua import ground_truth
log = logging.getLogger("axqua")
def _scalar_velocity_with_error(df: pd.DataFrame) -> tuple[pd.Series, pd.Series | None]:
"""Speed from the velocity components, with the component errors propagated.
Shared by ``SCALAR VELOCITY`` (TELEMAC) and ``U_MAG`` (OpenFOAM) so the two
cannot drift apart: they are the same quantity under two solvers' names, and
the error propagation is the part that would be easy to get subtly different.
sigma_V = (1/V) * sqrt(sum_i (c_i * sigma_i)^2)
"""
vel = ground_truth.scalar_velocity(df)
terms = [(df[c].astype(float) * df[f"{c}_err"].astype(float)) ** 2
for c in ("u", "v", "w") if c in df and f"{c}_err" in df]
err = None
if terms:
safe_v = vel.replace(0.0, np.nan)
err = (np.sqrt(sum(terms)) / safe_v).fillna(0.0)
return vel, err
def _resolve_quantity(df: pd.DataFrame, qty: str) -> tuple[pd.Series, pd.Series | None] | None:
"""Map a SELAFIN quantity name to (values, measured_error) from a tidy table.
Returns ``None`` when the tidy table lacks the columns to produce *qty*, so
callers can report missing quantities. ``measured_error`` is ``None`` when no
per-point error was measured (the caller falls back to the configured
fraction).
"""
q = qty.upper()
if q == "WATER DEPTH" and "h" in df:
return df["h"].astype(float), None
if q in ("TURBULENT ENERG.", "TURBULENT ENERGY") and "tke" in df:
return df["tke"].astype(float), None
if q in ("CUMUL BED EVOL", "BED EVOLUTION") and "dz" in df:
return df["dz"].astype(float), None
if q == "SCALAR VELOCITY" and "u_mag" in df:
err = df["u_mag_std"].astype(float) if "u_mag_std" in df else None
return df["u_mag"].astype(float), err
if q == "SCALAR VELOCITY" and any(c in df for c in ("u", "v", "w")):
return _scalar_velocity_with_error(df)
# --- OpenFOAM field names (axqua.solvers.openfoam.calibration) ------------
# HydroBayesCal's OpenFOAM binding extracts U_x/U_y/U_z/U_MAG/TKE, not the
# SELAFIN names, and compares them case-SENSITIVELY - so the spelling the
# caller passes is the spelling that must reach the CSV header.
for name, comp in (("U_X", "u"), ("U_Y", "v"), ("U_Z", "w")):
if q == name and comp in df:
err = df[f"{comp}_err"].astype(float) if f"{comp}_err" in df else None
return df[comp].astype(float), err
if q == "U_MAG" and "u_mag" in df:
err = df["u_mag_std"].astype(float) if "u_mag_std" in df else None
return df["u_mag"].astype(float), err
if q == "U_MAG" and any(c in df for c in ("u", "v", "w")):
return _scalar_velocity_with_error(df)
if q == "TKE" and "tke" in df:
return df["tke"].astype(float), None
if q == "VELOCITY U" and "u" in df:
return df["u"].astype(float), df["u_err"].astype(float) if "u_err" in df else None
if q == "VELOCITY V" and "v" in df:
return df["v"].astype(float), df["v_err"].astype(float) if "v_err" in df else None
return None
def _pick_table(tables: dict[str, pd.DataFrame], quantities: list[str]) -> pd.DataFrame:
"""Choose the tidy category whose columns satisfy all calibration quantities."""
for df in (*([tables["hydraulics"]] if "hydraulics" in tables else []), *tables.values()):
if all(_resolve_quantity(df, q) is not None for q in quantities):
return df
available = {cat: [c for c in df.columns if c not in ("x", "y", "z")]
for cat, df in tables.items()}
raise ValueError(
f"no ground-truth category provides all calibration quantities {quantities}. "
f"Available columns per category: {available}"
)
def _positional_errors(cfg, df, quantities) -> dict:
"""Per-quantity uncertainty induced by not knowing where the point is.
Returns ``{}`` unless the case declares a survey position quality, so every
existing config produces a byte-identical CSV. Never raises: a positional
error bar is an improvement to the error budget, and failing to compute one
must not cost the calibration its inputs.
"""
from axqua import positional_error as pe
# resolve_sigma's ValueError (an unknown quality name) is deliberately NOT
# caught: that is a typo in the config, and silently calibrating without the
# error bar the user asked for is the wrong response to it.
sigma, quality = pe.resolve_sigma(cfg)
if not sigma:
return {}
geometry = Path(cfg.model_path(cfg.geometry_slf))
results = Path(cfg.model_path("r2d.slf"))
if not (geometry.is_file() and results.is_file()):
log.warning(
"%s, but there is no converged 2D result to measure it against yet "
"(%s). The calibration CSV carries measurement error only; rebuild it "
"after the steady run to include the positional term.",
pe.describe(sigma, quality), results.name)
return {}
x = df["x"].to_numpy(float)
y = df["y"].to_numpy(float)
out = {}
for qty in quantities:
variable = _SELAFIN_VARIABLE.get(qty, qty)
try:
sampler = pe.telemac_sampler(geometry, results, variable)
except (KeyError, OSError) as exc:
log.debug("no positional error for %s: %s", qty, exc)
continue
spread = pe.sample_spread(sampler, x, y, sigma)
if np.any(spread > 0):
out[qty] = spread
log.info("%s -> %s positional error: median %.3f, max %.3f "
"(added in quadrature to the measurement error)",
pe.describe(sigma, quality), qty,
float(np.median(spread)), float(np.max(spread)))
return out
#: Calibration quantity -> the SELAFIN variable whose spatial variation carries
#: the positional uncertainty. A component velocity varies like the speed does.
_SELAFIN_VARIABLE = {
"U_x": "VELOCITY U", "U_y": "VELOCITY V", "U_z": "SCALAR VELOCITY",
"U_MAG": "SCALAR VELOCITY", "TKE": "SCALAR VELOCITY",
}
[docs]
def build_calibration_csv(cfg: Config, *,
floors: Mapping[str, float] | None = None,
path: Path | None = None,
z_resolver: "Callable[[pd.DataFrame], np.ndarray] | None" = None,
category: str | None = None,
) -> Path | None:
"""Write the HydroBayesCal calibration-points CSV from the tidy ground truth.
The tidy table is compiled from ``ground_truth.sources`` when given, otherwise
read from a user-authored ``ground_truth.measurements``. Returns ``None`` if neither
is configured.
*floors* clamps the per-quantity ``_ERROR`` column from below. This is a
correctness guard, not a nicety: HydroBayesCal builds each observation's
variance from this column, so a measured error that is zero (or a value near
zero multiplied by the relative fraction) gives that point a near-zero variance
and lets the least trustworthy measurement dominate the whole likelihood. The
default ``None`` leaves the existing behavior exactly as it was.
*path* overrides the output location, for a second CSV alongside the default
one (the OpenFOAM calibration targets the same points with different columns).
*z_resolver* replaces the ``z`` column. The default (``None``) writes the tidy
table's own ``z``, which is the **surveyed bed** - fine for a 2D calibration,
where ``z`` is ignored, and wrong for a 3D one, where it places the extraction
point on the bed and couples it to the survey's vertical accuracy. A solver
adapter supplies a resolver that places targets in the model's own column
instead (see :mod:`axqua.model_column`). It is a callback rather than a flag
because the two HydroBayesCal bindings disagree on what ``z`` means - OpenFOAM
wants an absolute elevation, TELEMAC-3D a height above the model bed - and that
knowledge belongs to the adapter, not here.
*category* selects one tidy tab by name instead of letting the quantities
decide. This is what lets a 2D and a 3D calibration read the *same* survey
under **different profile rules** - the 2D one depth-averaged, the 3D one
per-reading - by compiling them as separate categories (``hydraulics`` and
``hydraulics-3d``). They must not share a category: ``compile_ground_truth``
concatenates sources that do, which would silently hand both calibrations the
union of the two row sets.
"""
compile_ground_truth(cfg) # no-op when the table is user-supplied
if not cfg.ground_truth.sources and cfg.ground_truth.measurements is None \
and cfg.ground_truth.targets is None:
return None
if not Path(cfg.ground_truth_path).exists():
# e.g. a calibration-target template with only the parameters tab filled:
# no measurement points yet, but the parameter ranges still reach the HBC
# config via merged_parameters()
log.info("no ground-truth measurements compiled (yet); skipping the "
"calibration-points CSV")
return None
tables = ground_truth.read_tidy(cfg.ground_truth_path)
quantities = cfg.calibration.calibration_quantities
if category is not None:
if category not in tables:
raise ValueError(
f"ground-truth category {category!r} is not in "
f"{Path(cfg.ground_truth_path).name} (has: {', '.join(tables)}). "
"Add a ground_truth source with that category, or drop the "
"category= argument to pick the tab by its quantities.")
df = tables[category].reset_index(drop=True)
else:
df = _pick_table(tables, quantities).reset_index(drop=True)
if z_resolver is not None:
z_column = np.asarray(z_resolver(df), dtype=float)
elif "z" in df:
z_column = df["z"].astype(float)
else:
z_column = 0.0
out = pd.DataFrame({
"id": np.arange(1, len(df) + 1),
"x": df["x"].astype(float),
"y": df["y"].astype(float),
"z": z_column,
})
frac = cfg.calibration.measurement_error
positional = _positional_errors(cfg, df, quantities)
for qty in quantities:
values, measured_err = _resolve_quantity(df, qty)
err = measured_err if measured_err is not None else values.abs() * frac
floor = (floors or {}).get(qty)
if floor is not None:
err = err.clip(lower=float(floor))
extra = positional.get(qty)
if extra is not None:
# In quadrature: the position error and the instrument error are
# independent. Applied AFTER the floor, because the floor is a
# statement about the instrument and this is not.
err = np.sqrt(err.astype(float) ** 2 + np.asarray(extra) ** 2)
out[f"{qty}_DATA"] = values.round(6)
out[f"{qty}_ERROR"] = np.round(err, 6)
path = Path(path) if path is not None else cfg.calibration_path(cfg.calibration_csv)
path.parent.mkdir(parents=True, exist_ok=True)
out.to_csv(path, index=False)
return path
[docs]
def compile_ground_truth(cfg: Config) -> Path | None:
"""Compile configured raw ground-truth sources into the tidy table (if any)."""
return ground_truth.compile_ground_truth(cfg)
[docs]
def roughness_zone_parameters(cfg: Config) -> list:
"""``zone<N>`` calibration parameters declared by the roughness table.
A zone is a calibration parameter when it has a non-degenerate
``[ks_min, ks_max]`` prior and is not flagged ``calibration = False``. Writing
a genuine range is itself the declaration, so a table with bounds but no
``calibration`` column calibrates all of them rather than silently none; the
column is how you *opt a zone out*. Excluded zones are held at their nominal ks
- they still get a ``.tbl`` row and still influence the flow, they are just not
inferred.
Returns an empty list for the legacy single-``ks`` schema, which carries no
bounds; there the config's ``calibration.parameters`` remains the only source.
"""
from axqua.config import CalibrationParameter
from axqua.core.geodata import read_roughness_zones
if cfg.geodata.roughness_table is None:
return []
out = []
for zid, z in sorted(read_roughness_zones(cfg.geodata.roughness_table).items()):
if z.calibrate is False or z.ks_min is None or z.ks_max is None:
continue
if z.ks_max <= z.ks_min:
log.warning(" roughness zone %d flagged for calibration but its prior "
"is degenerate (ks_min=%g >= ks_max=%g) - held fixed at %g m",
zid, z.ks_min, z.ks_max, z.ks)
continue
out.append(CalibrationParameter(
name=f"zone{zid}", min=z.ks_min, max=z.ks_max,
comment=f"roughness zone {zid} ks [m] (Nikuradse)"))
return out
[docs]
def merged_parameters(cfg: Config) -> list:
"""Calibration parameters, merged from the three places they can be declared.
Precedence, later wins on a name collision: config ``calibration.parameters``
< the roughness table's ``zone<N>`` rows (``geodata.roughness_table`` with the
bounds schema) < the ``parameters`` tab of the filled calibration-target
template. The roughness table outranks the config because it is where the
zonation itself is defined - adding a zone there must not require editing the
parameter list in two files - while the spreadsheet stays the live override.
"""
from axqua import targets as targets_mod
by_name = {p.name: p for p in cfg.calibration.parameters}
for p in roughness_zone_parameters(cfg):
by_name[p.name] = p
for p in targets_mod.read_target_parameters(cfg):
by_name[p.name] = p
return list(by_name.values())
[docs]
def telemac_parameters(cfg: Config) -> list:
"""The declared parameters minus the ones that belong to another solver.
A case that calibrates both codes declares both sets in one
``calibration.parameters`` list, so this filter is the mirror image of
:func:`axqua.solvers.openfoam.calibration.openfoam_parameters`. Without it an
OpenFOAM ``ks`` / ``Cmu`` would be emitted into ``config_Telemac.py`` as though
it were a TELEMAC keyword, and HydroBayesCal would spend the whole campaign
rewriting a ``.cas`` line that does not exist.
"""
from axqua.solvers.openfoam.spec import OPENFOAM_PARAMETERS
keep, drop = [], []
for p in merged_parameters(cfg):
# a zone<N> roughness parameter is TELEMAC's and must survive: only the
# names that are exclusively OpenFOAM's are removed
if str(p.name).strip().lower() in OPENFOAM_PARAMETERS:
drop.append(p.name)
else:
keep.append(p)
if drop:
log.info("not TELEMAC calibration parameters, skipped: %s (they belong to "
"the OpenFOAM calibration)", ", ".join(sorted(drop)))
return keep
[docs]
def emit_hbc_config(cfg: Config, calibration_csv: Path | None) -> Path:
"""Emit a HydroBayesCal ``config_Telemac.py`` for the produced case."""
c = cfg.calibration
parameters = telemac_parameters(cfg)
params = [p.name for p in parameters]
ranges = [[p.min, p.max] for p in parameters]
# which case HydroBayesCal drives: the built steady 2D case by default, or the
# 3D steering file when calibrating vertical profiles (calibration.solver_name)
control_file = c.control_file or cfg.cas_file
if c.results_base:
results_base = c.results_base
elif c.solver_name == "Telemac3d":
# HydroBayesCal sets '3D RESULT FILE' from this base and derives the
# companion 2D file as <base>_<id>_2d.slf, so the base must be the 3D one.
results_base = Path(cfg.results3d_slf).stem
else:
results_base = Path(cfg.results_slf).stem
# HydroBayesCal resolves this against <model_dir>/user_fortran/, and rewrites a
# line in it for any calibration parameter prefixed `f.` - which is how the
# gain-lose conductivity (`f.HMP_KF`) is calibrated. Only meaningful when a
# gain-lose reach actually generated the routine.
fortran_file = ("user_rain.f"
if (cfg.gain_lose.active
and cfg.gain_lose.implementation == "fortran") else None)
def pylist(items, q=True):
if q:
return "[" + ", ".join(repr(str(i)) for i in items) + "]"
return "[" + ", ".join(str(i) for i in items) + "]"
# bal_telemac.py reads config.morphodynamic_simulation['gaia_cas'] /
# ['gaia_results_filename_base'] unconditionally, so always emit the block;
# the GAIA values are None when morphodynamics is off (hydraulic-only case).
if cfg.morphodynamics.enabled:
gaia_cas = repr(cfg.gaia_cas)
gaia_base = repr(cfg.morphodynamics.gaia_results_base)
else:
gaia_cas = gaia_base = "None"
morph = (
"morphodynamic_simulation = {\n"
f" 'gaia_cas': {gaia_cas},\n"
f" 'gaia_results_filename_base': {gaia_base},\n"
"}\n\n"
)
text = f'''"""HydroBayesCal TELEMAC config generated by axqua for case '{cfg.name}'.
Run a surrogate-assisted Bayesian calibration with:
python bal_telemac.py --config {cfg.hbc_config}
"""
import os
BASE_DIR = {str(cfg.model_dir)!r}
paths = {{
'case_template_dir': BASE_DIR,
'model_dir': BASE_DIR,
'res_dir': {str(cfg.calibration_dir)!r},
'calibration_pts_file_path': {(str(calibration_csv) if calibration_csv else "")!r},
}}
hydrodynamic_simulation = {{
'solver_name': {c.solver_name!r},
'n_processors': {cfg.telemac.n_processors},
'results_filename_base': {results_base!r},
'control_file': {control_file!r},
'friction_file': {cfg.friction_tbl!r},
'fortran_file': {fortran_file!r},
}}
{morph}calibration = {{
'parameters': {pylist(params)},
'param_values': {ranges},
'extraction_quantities': {pylist(c.extraction_quantities)},
'calibration_quantities': {pylist(c.calibration_quantities)},
'dict_output_name': 'extraction-data',
}}
# WHEN to read the model outputs. "last" takes the converged frame; the stock
# default "mean_last" averages the final n_last frames, which on a run marching to
# steady state from a dry or pre-wetted start folds the residual transient into the
# values the surrogate is trained on.
# Emitting it here rather than relying on axqua's source patch of the staged
# driver: the patch rewrites a string literal in someone else's file, so it silently
# no-ops if that literal ever changes. This says it in the config, which is the
# interface the drivers actually document. (The patch is kept as well - hydrobayescal
# 1.4.3's multiflow driver ignores this block, and only honors it from 1.4.4 on.)
# n_last is unused while output_extraction_time is 'last'; 1 is the value that keeps
# the intent if anyone flips this back to 'mean_last'.
extraction = {{
'output_extraction_time': 'last',
'n_last': 1,
}}
sampling = {{
'init_runs': {c.init_runs},
'max_runs': {c.max_runs},
'parameter_distribution': 'uniform',
'parameter_sampling_method': {c.parameter_sampling_method!r},
# HydroBayesCal >= 1.5: grow the initial design in Sobol blocks and stop when
# it is measurably sufficient. init_runs is the ceiling, so this can only save
# runs; the saved ones become BAL iterations.
'adaptive_init_runs': {bool(c.adaptive_init_runs)},
'init_runs_min': {c.init_runs_min!r},
# exploit the posterior until >1 well-separated mode is found, then explore
'bal_exploration_tradeoff': {c.bal_exploration_tradeoff!r},
'tp_selection_criteria': 'dkl',
'eval_steps': 1,
# prior_samples drives a (prior_samples x prior_samples) dense covariance PER
# location in gpytorch; kept modest (calibration.prior_samples) to stay in RAM.
'prior_samples': {c.prior_samples},
'mc_samples_al': 2000,
'mc_exploration': 1000,
'gp_library': {c.gp_library!r},
}}
execution = {{
'complete_bal_mode': True,
'only_bal_mode': False,
'delete_complex_outputs': True,
'validation': False,
'user_param_values': False,
}}
'''
path = cfg.calibration_path(cfg.hbc_config)
path.write_text(text)
return path