Cette page a été traduite automatiquement de l’anglais et peut contenir des erreurs. Original en anglais
Code source de axqua.bayescal
"""HydroBayesCal integration: build inputs, emit config, launch calibration.
Keeps the per-case ``run_Bayes_cal.py`` / ``run_Bayes_cal_multiflow.py`` scripts
thin - they only declare case-specific choices (calibration target, and for the
multi-flow case the list of :class:`FlowSpec`) and call
:func:`run_single_flow_calibration` / :func:`run_multiflow_calibration`.
Two calibration modes:
* **single-flow** - one steady discharge, ground truth compiled from the
configured ``ground_truth.sources`` (see :mod:`axqua.ground_truth`); the
case's built ``.cas`` is calibrated directly with HydroBayesCal's stock
``bal_telemac.py``.
* **multi-flow** - several steady discharges, each with its own field campaign
(:mod:`axqua.campaigns`) and its own steering file derived from the built
base case; run with the additive ``bal_telemac_multiflow.py`` +
``MultiflowTelemacModel``. One shared roughness is calibrated against all
campaigns jointly.
**HydroBayesCal is a pip dependency** (``pip install 'axqua[calibration]'``),
not an external checkout: the drivers ship with the package from 1.4.2 on
(``hydroBayesCal.copy_driver``), so there is no directory to point at and no
second conda environment to name. It is an *optional* extra because its surrogate
stack - bayesvalidrox, gpytorch, scikit-learn - is heavy and irrelevant to
building or running a TELEMAC case, so :func:`require_hbc` imports it lazily and
explains the install if it is missing.
A staged driver is still run as a **subprocess**, but in this interpreter's own
environment: the driver is a script with an ``argparse`` ``main()``, and it must
run with TELEMAC sourced (every surrogate iteration launches the solver), which
mutates the environment in ways not to inflict on the caller's process. Pass
``checkout=`` to develop against an unreleased driver from a source tree.
Lessons baked in (see the module functions):
* the built ``.cas`` must output ``SCALAR VELOCITY`` (graphic printout ``M``);
* the staged driver extracts the **last** (converged) SELAFIN frame, not the
``mean_last`` window that would average a dry-start transient;
* HydroBayesCal's prior-sample count drives a per-location dense covariance in
gpytorch, so it is kept modest (``calibration.prior_samples``, default 5000) -
ample for a low-dimensional roughness posterior and safe on a 16 GB box;
* per-flow steering files differ only in Q, the outflow stage (from the
synthesized rating at that Q), ``DURATION`` and the results name.
"""
from __future__ import annotations
import logging
import shutil
import sys
from contextlib import contextmanager
from dataclasses import dataclass
from pathlib import Path
import numpy as np
import pandas as pd
from axqua import campaigns, hbc
from axqua.config import Config
from axqua.ground_truth import compile_ground_truth, read_tidy
from axqua.rating import synthesize_outflow_rating
log = logging.getLogger("axqua")
# --------------------------------------------------------------------------- #
# install location
# --------------------------------------------------------------------------- #
# require_hbc / _hbc_version / stage_driver now live in axqua.hbc (they are the
# same whichever solver is calibrated). Re-exported, not reimplemented, so every
# existing import of them from axqua.bayescal keeps working.
require_hbc = hbc.require_hbc
_hbc_version = hbc.hbc_version
stage_driver = hbc.stage_driver
# --------------------------------------------------------------------------- #
# shared HydroBayesCal plumbing
# --------------------------------------------------------------------------- #
[docs]
def ensure_graphic_output(cas_path: Path, code: str = "M") -> bool:
"""Ensure the ``.cas`` writes a graphic-printout variable (default ``M`` =
SCALAR VELOCITY). Returns True if the file was patched. HydroBayesCal reads
the quantity straight from the results SELAFIN, so it must be printed.
"""
cas_path = Path(cas_path)
if not cas_path.exists():
raise SystemExit(f"no built case at {cas_path}; run preprocessing.py first")
import re
text = cas_path.read_text()
m = re.search(r"(VARIABLES FOR GRAPHIC PRINTOUTS\s*:\s*')([^']*)(')", text)
if not m:
print(f" warning: no VARIABLES FOR GRAPHIC PRINTOUTS in {cas_path.name}")
return False
codes = [c.strip() for c in m.group(2).split(",") if c.strip()]
if code in codes:
return False
codes.append(code)
cas_path.write_text(text[:m.start()] + m.group(1) + ",".join(codes)
+ m.group(3) + text[m.end():])
print(f" added graphic printout {code!r} to {cas_path.name}")
return True
[docs]
def other_calibration_running(model_dir: Path) -> bool:
"""A TELEMAC run from this model dir already active (would fight over the
shared friction.tbl)."""
return hbc.solver_running("telemac2d", model_dir, extra=("steady2d",))
def _telemac_environment(cfg: Config):
"""The TELEMAC environment a staged driver is launched inside.
``environment`` is read defensively: the previous implementation needed nothing
but ``telemac.pysource``, so requiring the newer structured block would narrow
what a caller may pass for no gain (``from_config(None, ...)`` already falls
back to a POSIX environment against the legacy script).
"""
return hbc.solver_environment(getattr(cfg.telemac, "environment", None),
cfg.telemac.pysource)
[docs]
def launch(cfg: Config, driver: Path, config_path: Path, *, env: str | None = None,
note: str = "HydroBayesCal calibration") -> int:
"""Run a staged TELEMAC driver with TELEMAC sourced. See :func:`axqua.hbc.launch_driver`.
*env* is accepted and ignored; it named the separate conda environment this used
to need, before HydroBayesCal became a dependency of axqua itself.
"""
if env:
print(f"note: env={env!r} ignored - HydroBayesCal now runs in axqua's own "
"environment (it is a dependency), so no separate env is needed")
return hbc.launch_driver(_telemac_environment(cfg), driver, config_path, note=note)
# --------------------------------------------------------------------------- #
# single-flow calibration
# --------------------------------------------------------------------------- #
[docs]
def build_velocity_csv(cfg: Config, *, vel_err_floor: float = campaigns.VELOCITY_ERROR_FLOOR,
depth_err_floor: float = campaigns.DEPTH_ERROR_FLOOR
) -> tuple[Path, pd.DataFrame]:
"""Compile the configured ground truth and write the calibration-points CSV.
The velocity target is the horizontal FlowTracker speed (0.6-depth proxy for
the depth-averaged velocity) with propagated, floored error; water depth is
carried alongside. Reads ``ground_truth.sources`` via
:func:`axqua.ground_truth.compile_ground_truth`.
"""
compile_ground_truth(cfg)
tables = read_tidy(cfg.ground_truth_path)
if "hydraulics" not in tables:
raise SystemExit(f"no 'hydraulics' tab in {cfg.ground_truth_path}; check "
"ground_truth.sources in case-config.yml")
df = tables["hydraulics"].reset_index(drop=True)
for col in ("x", "y", "u", "v", "h"):
if col not in df.columns:
raise SystemExit(f"hydraulics ground truth missing '{col}' (has {list(df.columns)})")
out = campaigns.assemble_velocity_csv(
df["x"], df["y"], df["z"] if "z" in df.columns else 0.0,
df["u"], df["v"],
df["u_err"] if "u_err" in df.columns else 0.0,
df["v_err"] if "v_err" in df.columns else 0.0,
df["h"], [f"pt-{i + 1}" for i in range(len(df))],
vel_err_floor, depth_err_floor).drop(columns="label")
out = _inside_the_model(cfg, out)
path = cfg.calibration_path(cfg.calibration_csv)
path.parent.mkdir(parents=True, exist_ok=True)
out.to_csv(path, index=False)
return path, out
def _inside_the_model(cfg: Config, points: pd.DataFrame) -> pd.DataFrame:
"""Leave out the measurement points that lie outside the model.
A campaign often includes the cross sections at which the inflow was gauged,
upstream of where the model begins. HydroBayesCal reads the model at the node
nearest to a point with no limit on the distance, so such a point is compared with
a node on the model boundary: a velocity measured in mid-channel against one at a
wall. Of the 56 verticals of the Isar campaign of November 2025, 23 lie upstream of
the model.
"""
from shapely.geometry import Point
outline = cfg.geodata.boundary
if outline is None or not Path(outline).exists() or points.empty:
return points
try:
import geopandas as gpd
layer = gpd.read_file(outline)
if layer.crs is not None and layer.crs.to_epsg() != cfg.crs_epsg:
layer = layer.to_crs(epsg=cfg.crs_epsg)
shape = layer.geometry.union_all() if hasattr(layer.geometry, "union_all") \
else layer.geometry.unary_union
if shape.geom_type in ("LineString", "MultiLineString"):
from shapely.ops import polygonize, unary_union
shape = unary_union(list(polygonize(shape)))
inside = [shape.covers(Point(float(x), float(y)))
for x, y in zip(points["x"], points["y"])]
except Exception as exc: # noqa: BLE001 - never lose the data
log.warning("could not tell which measurement points lie inside the model "
"(%s); all of them are used", exc)
return points
outside = len(inside) - sum(inside)
if outside:
log.warning("%d of %d measurement points lie outside the model outline and "
"are left out", outside, len(inside))
points = points.loc[inside].reset_index(drop=True)
if "id" in points.columns:
points["id"] = range(1, len(points) + 1)
return points
[docs]
def run_single_flow_calibration(
cfg: Config, *, calibration_quantities=None,
extraction_quantities=None,
vel_err_floor: float = campaigns.VELOCITY_ERROR_FLOOR,
depth_err_floor: float = campaigns.DEPTH_ERROR_FLOOR,
prepare_only: bool = False, checkout: Path | None = None,
env: str | None = None) -> int:
"""End-to-end single-discharge calibration of the built case.
Compiles the calibration CSV, ensures the ``.cas`` prints the target,
emits the HydroBayesCal config, stages the driver and (unless
``prepare_only``) launches it. Returns the driver exit code (0 when
``prepare_only``).
"""
from axqua.calibration import emit_hbc_config
cfg.ensure_dirs()
# The case file decides what is calibrated against. Until 2026-10-17 a built-in
# ("SCALAR VELOCITY",) replaced it whenever the caller named nothing, which is
# what every job did: a case asking for depth and velocity was calibrated
# against velocity alone, without a word.
if calibration_quantities is None:
calibration_quantities = (cfg.calibration.calibration_quantities
or ["SCALAR VELOCITY"])
if extraction_quantities is None:
extraction_quantities = (cfg.calibration.extraction_quantities
or ["WATER DEPTH", "SCALAR VELOCITY"])
extraction_quantities = list(dict.fromkeys([*extraction_quantities,
*calibration_quantities]))
cfg.calibration.calibration_quantities = list(calibration_quantities)
cfg.calibration.extraction_quantities = list(extraction_quantities)
print(f"calibration target(s): {list(calibration_quantities)} "
f"parameters: {[p.name for p in cfg.calibration.parameters]}")
csv, table = build_velocity_csv(cfg, vel_err_floor=vel_err_floor,
depth_err_floor=depth_err_floor)
print(f"calibration points ({len(table)}) -> {csv}")
if "SCALAR VELOCITY_DATA" in table:
print(f" scalar velocity [m/s]: {table['SCALAR VELOCITY_DATA'].min():.3f} .. "
f"{table['SCALAR VELOCITY_DATA'].max():.3f}")
ensure_graphic_output(cfg.model_path(cfg.cas_file))
config_path = emit_hbc_config(cfg, csv)
print(f"HydroBayesCal config -> {config_path}")
driver = stage_driver(cfg.calibration_dir, "bal_telemac.py", checkout=checkout)
if prepare_only:
print("\n--prepare-only: not launched. Run it with:")
print(f" cd {driver.parent} && {sys.executable} {driver.name} "
f"--config {config_path}")
return 0
with case_as_built(cfg):
return launch(cfg, driver, config_path, env=env,
note="HydroBayesCal calibration")
#: Folder in the calibration directory that holds the files of the built case while a
#: calibration rewrites them.
AS_BUILT = "case-as-built"
def _rewritten_by_the_driver(cfg: Config) -> list[str]:
"""The files of the built case a HydroBayesCal run writes into, relative to it."""
names = [cfg.calibration.control_file or cfg.cas_file, cfg.friction_tbl]
if cfg.gain_lose.active and cfg.gain_lose.implementation == "fortran":
names.append("user_fortran/user_rain.f")
if cfg.morphodynamics.enabled:
names.append(cfg.gaia_cas)
return [str(name) for name in names if name]
[docs]
@contextmanager
def case_as_built(cfg: Config):
"""Hand the built case back the way the calibration found it.
HydroBayesCal runs TELEMAC in the folder of the built case and writes every
parameter set it tests into the files there: the roughness values into the
friction table, the name of the results file into the steering file. When it has
finished, the case holds the **last set it tried**, which is neither the case
that was built nor the calibrated one, and the next steady run uses it without a
word. On the Isar example the friction table read ks = 0.195 m in zone 4 after a
calibration of a case built with 0.05 m.
The files are copied into the calibration folder before the driver starts and
copied back when it returns, whether it succeeded or not. A calibration that is
killed restores nothing; see :func:`restore_as_built` for that case.
"""
model = Path(cfg.model_dir)
keep = cfg.calibration_path(AS_BUILT)
restore_as_built(cfg)
for name in _rewritten_by_the_driver(cfg):
if (model / name).is_file():
(keep / name).parent.mkdir(parents=True, exist_ok=True)
shutil.copy2(model / name, keep / name)
try:
yield
finally:
restore_as_built(cfg)
[docs]
def restore_as_built(cfg: Config) -> bool:
"""Copy the files a calibration set aside back into the built case.
Returns whether there was anything to restore. A copy is only ever found here
outside a running calibration when that calibration was killed, by a cancelled job
or a power cut, before it could restore the case itself. The copy is then the case
as built, and what is in the model folder is the last parameter set the driver had
reached. Every TELEMAC run of the case calls this first, so a cancelled calibration
cannot leak into the next simulation.
"""
keep = cfg.calibration_path(AS_BUILT)
if not keep.is_dir():
return False
model = Path(cfg.model_dir)
for source in sorted(path for path in keep.rglob("*") if path.is_file()):
target = model / source.relative_to(keep)
target.parent.mkdir(parents=True, exist_ok=True)
shutil.copy2(source, target)
shutil.rmtree(keep, ignore_errors=True)
return True
# --------------------------------------------------------------------------- #
# multi-flow calibration
# --------------------------------------------------------------------------- #
[docs]
@dataclass
class FlowSpec:
"""One steady discharge in a multi-flow calibration.
``kind`` + ``values`` (+ ``positions``) feed :func:`axqua.campaigns`; the
steering file is derived from the built base case, differing only in Q, the
outflow stage, ``DURATION`` and the results name. ``vel_err_floor`` raises
the per-point velocity error for a campaign that should inform but not
dominate the joint likelihood (e.g. a single cross-section).
"""
name: str # short tag (steering + results names)
discharge: float # m3/s
kind: str # adapter | transect | inline | points | csv
values: Path # FlowTracker xlsx, point layer ("points") or a ready CSV ("csv")
positions: Path | None = None # DGPS point layer (adapter/transect)
sheet: str = "KB15"
duration: float | None = None # steady sim seconds (default: cfg duration)
vel_err_floor: float = campaigns.VELOCITY_ERROR_FLOOR
# Continue from an already-converged result instead of re-filling the reach
# from the pre-wet seed on every calibration run. The flow field then only has
# to re-equilibrate to the perturbed roughness, which is a small fraction of
# the time a fill takes - the difference between a feasible and an infeasible
# calibration on a large mesh. See write_flow_cas for the caveat on DURATION.
hotstart_from: str | None = None # e.g. "r2d.slf"
[docs]
def compile(self, crs_epsg: int) -> pd.DataFrame:
if self.kind == "csv":
# pre-compiled calibration points (already in the HydroBayesCal
# schema: id, x, y, z, <QTY>_DATA, <QTY>_ERROR); passed through so
# case-specific target preparation (e.g. bathymetry-corrected
# depths, profile-averaged velocities) survives the launcher's
# recompilation step.
return pd.read_csv(self.values)
if self.kind == "adapter":
return campaigns.compile_adapter(self.values, self.positions, crs_epsg,
self.vel_err_floor, self.name)
if self.kind == "transect":
return campaigns.compile_transect(self.values, self.positions,
self.sheet, self.vel_err_floor,
label=self.name)
if self.kind == "inline":
return campaigns.compile_inline(self.values, self.sheet,
self.vel_err_floor, label=self.name)
if self.kind == "points":
return campaigns.compile_points(self.values, crs_epsg,
self.vel_err_floor, label=self.name)
raise ValueError(f"flow {self.name!r}: unknown kind {self.kind!r}")
def _patch_cas(text: str, key: str, value: str) -> str:
import re
new, n = re.subn(rf"^{re.escape(key)}\s*:.*$", f"{key} : {value}", text,
count=1, flags=re.M)
if n != 1:
raise ValueError(f"keyword {key!r} not found in the base steering file")
return new
def _set_cas(text: str, key: str, value: str) -> str:
"""Like :func:`_patch_cas` but appends the keyword when the base case has no
such line (a hotstart adds keywords a pre-wet start never carried)."""
import re
new, n = re.subn(rf"^{re.escape(key)}\s*:.*$", f"{key} : {value}", text,
count=1, flags=re.M)
if n == 1:
return new
return text.rstrip("\n") + f"\n{key} : {value}\n"
[docs]
def outflow_stage(cfg: Config, discharge: float, out_dir: Path) -> float:
"""Outflow water-surface elevation at ``discharge`` from the synthetic
normal-flow rating (extrapolation warning: unreliable far above the geometry
calibration discharge / overbank)."""
rating = synthesize_outflow_rating(cfg, discharge,
out=out_dir / f"rating-q{discharge:g}.csv")
return float(pd.read_csv(rating)["WSE"].iloc[0])
[docs]
def write_flow_cas(cfg: Config, base_cas: Path, flow: FlowSpec, rating_dir: Path,
*, duration: float | None = None, dest_dir: Path | None = None,
listing: int | None = None, graphic: int | None = None) -> Path:
"""Derive one flow's steering file from the built base ``.cas``.
With ``flow.hotstart_from`` the run continues from that result file instead of
the base case's pre-wet/dry seed. ``DURATION`` then measures *re-equilibration*
after the roughness change, not the fill - but it must still be long enough for
the field to forget the seed, or every run is pulled toward the roughness the
seed was converged at and the calibration under-reports sensitivity. Verify it
at the prior extremes before trusting a short hotstart duration.
"""
text = base_cas.read_text()
wse = outflow_stage(cfg, flow.discharge, rating_dir)
dur = duration if duration is not None else (flow.duration
or cfg.hydrodynamics.duration)
text = _patch_cas(text, "TITLE", f"'{cfg.name} steady {flow.name}'")
if flow.hotstart_from:
text = _set_cas(text, "PREVIOUS COMPUTATION FILE", flow.hotstart_from)
text = _set_cas(text, "PREVIOUS COMPUTATION FILE FORMAT", "'SERAFIN'")
text = _set_cas(text, "INITIAL TIME SET TO ZERO", "YES")
text = _patch_cas(text, "PRESCRIBED FLOWRATES", f"0.;{flow.discharge:.4f}")
text = _patch_cas(text, "PRESCRIBED ELEVATIONS", f"{wse:.4f};0.")
text = _patch_cas(text, "DURATION", f"{dur}")
text = _patch_cas(text, "RESULTS FILE", f"r2d-{flow.name}.slf")
if listing is not None:
text = _patch_cas(text, "LISTING PRINTOUT PERIOD", str(listing))
if graphic is not None:
text = _patch_cas(text, "GRAPHIC PRINTOUT PERIOD", str(graphic))
dest = (dest_dir or base_cas.parent) / f"steady2d-{flow.name}.cas"
dest.write_text(text)
print(f" {dest.name}: Q={flow.discharge} m3/s, outflow WSE={wse:.3f} m, "
f"DURATION={dur} s")
return dest
[docs]
def emit_multiflow_config(cfg: Config, out_dir: Path, model_dir: Path,
flows: list[FlowSpec], csv_paths: dict[str, Path], *,
n_cpus: int, init_runs: int, max_runs: int,
complete_bal_mode: bool = True,
only_bal_mode: bool = False,
prior_samples: int | None = None,
adaptive_init_runs: bool | None = None) -> Path:
"""Write ``config_Telemac_multiflow.py`` (standard blocks + a ``multiflow``
block of per-flow steering/CSV entries)."""
from axqua.calibration import merged_parameters
c = cfg.calibration
# merged_parameters, not c.parameters: the roughness table's calibrated zones
# and the target template's rows must reach the multiflow config too
parameters = merged_parameters(cfg)
params = [p.name for p in parameters]
ranges = [[p.min, p.max] for p in parameters]
prior = prior_samples if prior_samples is not None else getattr(
c, "prior_samples", 5000)
adaptive = bool(getattr(c, "adaptive_init_runs", True)
if adaptive_init_runs is None else adaptive_init_runs)
flow_entries = ",\n".join(
" {{'name': {n!r}, 'control_file': {cf!r},\n"
" 'results_filename_base': {rb!r},\n"
" 'calibration_pts_file_path': {csv!r}}}".format(
n=f.name, cf=f"steady2d-{f.name}.cas", rb=f"r2d-{f.name}",
csv=str(csv_paths[f.name]))
for f in flows)
text = f'''"""HydroBayesCal multi-discharge config generated by axqua for '{cfg.name}'.
Run: python bal_telemac_multiflow.py --config config_Telemac_multiflow.py
"""
import os
BASE_DIR = {str(model_dir)!r}
paths = {{
'case_template_dir': BASE_DIR,
'model_dir': BASE_DIR,
'res_dir': {str(out_dir)!r},
'calibration_pts_file_path': {str(csv_paths[flows[0].name])!r},
}}
multiflow = {{
'flows': [
{flow_entries},
],
}}
hydrodynamic_simulation = {{
'solver_name': 'Telemac2d',
'n_processors': {n_cpus},
'results_filename_base': {f"r2d-{flows[0].name}"!r},
'control_file': {f"steady2d-{flows[0].name}.cas"!r},
'friction_file': {cfg.friction_tbl!r},
'fortran_file': None,
}}
morphodynamic_simulation = {{
'gaia_cas': None,
'gaia_results_filename_base': None,
}}
calibration = {{
'parameters': {params!r},
'param_values': {ranges!r},
'extraction_quantities': {list(c.extraction_quantities)!r},
'calibration_quantities': {list(c.calibration_quantities)!r},
'dict_output_name': 'extraction-data',
}}
sampling = {{
'init_runs': {init_runs},
'max_runs': {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 - and here one run means one TELEMAC simulation PER FLOW.
'adaptive_init_runs': {adaptive},
'init_runs_min': {getattr(c, "init_runs_min", None)!r},
# exploit the posterior until >1 well-separated mode is found, then explore
'bal_exploration_tradeoff': {getattr(c, "bal_exploration_tradeoff", "auto")!r},
'tp_selection_criteria': 'dkl',
'eval_steps': 1,
# prior_samples drives a (prior_samples x prior_samples) dense covariance PER
# location in gpytorch; keep modest (a low-dim roughness posterior needs few).
'prior_samples': {prior},
'mc_samples_al': 2000,
'mc_exploration': 1000,
'gp_library': {c.gp_library!r},
}}
execution = {{
'complete_bal_mode': {complete_bal_mode},
'only_bal_mode': {only_bal_mode},
'delete_complex_outputs': True,
'validation': False,
'user_param_values': False,
}}
extraction = {{
# Read the CONVERGED final frame, not a mean over the last n frames: these runs
# march to steady state from a pre-wetted start, so a mean_last window folds the
# residual transient into the values the surrogate is trained on. The multiflow
# driver silently defaulted to 'mean_last' until HydroBayesCal 1.5 made it read
# this block (single-flow runs on the same config already honored it), so the
# two modes were fitted to different data - emit it explicitly.
'output_extraction_time': 'last',
'n_last': 80,
}}
'''
path = Path(out_dir) / "config_Telemac_multiflow.py"
path.parent.mkdir(parents=True, exist_ok=True)
path.write_text(text)
return path
def _restart_ready(multiflow_dir: Path, flows: list[FlowSpec]) -> list[str]:
"""Names of flows whose completed initial design is NOT on disk."""
return [f.name for f in flows if not (
multiflow_dir / f"flow-{f.name}" / "auto-saved-results-HydroBayesCal"
/ "restart_data" / "initial-model-outputs.json").exists()]
[docs]
def compile_flow_csvs(cfg: Config, flows: list[FlowSpec],
multiflow_dir: Path) -> dict[str, Path]:
"""Compile each flow's calibration-points CSV; returns {flow name: path}."""
multiflow_dir.mkdir(parents=True, exist_ok=True)
paths: dict[str, Path] = {}
print("compiling per-flow calibration points:")
for f in flows:
df = f.compile(cfg.crs_epsg)
out = multiflow_dir / f"measurements-calibration-{f.name}.csv"
df.to_csv(out, index=False)
paths[f.name] = out
print(f" {f.name} (Q={f.discharge} m3/s): {len(df)} points, "
f"U {df['SCALAR VELOCITY_DATA'].min():.3f}..{df['SCALAR VELOCITY_DATA'].max():.3f} "
f"m/s -> {out.name}")
return paths
[docs]
def run_multiflow_smoke(cfg: Config, flows: list[FlowSpec], multiflow_dir: Path,
csv_paths: dict[str, Path], *, checkout: Path | None = None,
env: str | None = None, init_runs: int = 3) -> int:
"""Isolated end-to-end plumbing test: tiny 30 s runs in a copied model dir.
Verifies steering/friction updates, extraction at every point, per-run JSON
accumulation (rows must differ across ``init_runs`` distinct ks levels) and
the combined matrix shape - without the full solver cost.
"""
smoke = multiflow_dir / "smoke"
model = smoke / "model"
if smoke.exists():
shutil.rmtree(smoke)
model.mkdir(parents=True)
src = cfg.model_dir
inputs = ["geometry.slf", "boundaries.cli", "initial-conditions.slf"]
# a hotstarted flow reads its seed instead of (or besides) the pre-wet initial
# conditions, so it has to exist in the smoke model dir too
inputs += sorted({f.hotstart_from for f in flows if f.hotstart_from})
for name in inputs:
if (src / name).exists():
(model / name).symlink_to(src / name)
elif name not in ("initial-conditions.slf",):
raise SystemExit(f"smoke: {name} missing from {src}")
shutil.copy2(src / cfg.friction_tbl, model / cfg.friction_tbl)
base_cas = src / cfg.cas_file
print("smoke steering files (DURATION 30 s):")
for f in flows:
write_flow_cas(cfg, base_cas, f, smoke, duration=30.0, dest_dir=model,
listing=50, graphic=50)
config_path = emit_multiflow_config(
cfg, smoke, model, flows, csv_paths, n_cpus=2, init_runs=init_runs,
max_runs=init_runs, complete_bal_mode=False,
# the plumbing checks below assert exactly init_runs rows; the adaptive
# ladder may stop at a different size, so pin the design for the smoke test
adaptive_init_runs=False)
driver = stage_driver(smoke, "bal_telemac_multiflow.py", checkout=checkout,
also=("bal_telemac.py",))
rc = launch(cfg, driver, config_path, env=env,
note=f"multiflow SMOKE ({init_runs * len(flows)} tiny runs)")
if rc != 0:
return rc
combined = (smoke / "auto-saved-results-HydroBayesCal" / "calibration-data"
/ "_".join(cfg.calibration.calibration_quantities)
/ "model-results-multiflow.csv")
if not combined.exists():
print(f"missing combined results {combined}")
return 3
arr = np.loadtxt(combined, delimiter=",", skiprows=1, ndmin=2)
n_pts = (sum(len(pd.read_csv(p)) for p in csv_paths.values())
* len(cfg.calibration.calibration_quantities))
checks = {
"shape": arr.shape == (init_runs, n_pts),
"no all-zero run": bool(np.all(arr.any(axis=1))),
"runs differ (ks varies)":
len({tuple(np.round(r, 6)) for r in arr}) == init_runs,
}
for name, ok in checks.items():
print(f" smoke check - {name}: {'OK' if ok else 'FAIL'}")
return 0 if all(checks.values()) else 3
[docs]
def run_multiflow_calibration(
cfg: Config, flows: list[FlowSpec], *, launch_mode: str = "prepare",
checkout: Path | None = None, env: str | None = None,
force: bool = False) -> int:
"""Compile CSVs + per-flow steering, emit config, stage drivers, then act by
``launch_mode``:
* ``prepare`` - write everything, do not launch (default);
* ``smoke`` - run the isolated plumbing test;
* ``run`` - launch the full calibration (initial design + BAL);
* ``resume`` - reuse completed per-flow initial designs (only_bal_mode) and
go straight to surrogate + BAL.
"""
cfg.ensure_dirs()
multiflow_dir = cfg.calibration_path("multiflow")
csv_paths = compile_flow_csvs(cfg, flows, multiflow_dir)
if launch_mode == "smoke":
return run_multiflow_smoke(cfg, flows, multiflow_dir, csv_paths,
checkout=checkout, env=env)
base_cas = cfg.model_path(cfg.cas_file)
if not base_cas.exists():
raise SystemExit(f"no built case at {base_cas}; run preprocessing.py first")
print("per-flow steering files:")
for f in flows:
write_flow_cas(cfg, base_cas, f, multiflow_dir)
resume = launch_mode == "resume"
if resume:
missing = _restart_ready(multiflow_dir, flows)
if missing:
raise SystemExit(f"resume: no completed initial design for {missing}; "
"run the initial design (launch_mode='run') first")
config_path = emit_multiflow_config(
cfg, multiflow_dir, cfg.model_dir, flows, csv_paths,
n_cpus=cfg.telemac.n_processors, init_runs=cfg.calibration.init_runs,
max_runs=cfg.calibration.max_runs, only_bal_mode=resume)
print(f"multiflow config -> {config_path}"
+ (" (only_bal_mode: reusing completed initial designs)" if resume else ""))
driver = stage_driver(multiflow_dir, "bal_telemac_multiflow.py",
checkout=checkout, also=("bal_telemac.py",))
if launch_mode == "prepare":
print("\nprepared (not launched). Launch with launch_mode='run' "
"(or 'resume' to reuse completed initial designs).")
print(f" cd {driver.parent} && {sys.executable} {driver.name} "
f"--config {config_path}")
return 0
if other_calibration_running(cfg.model_dir) and not force:
raise SystemExit("another TELEMAC calibration appears active on this case "
"(shared friction.tbl); wait or pass force=True")
note = ("multi-discharge calibration - BAL only (reusing initial designs)"
if resume else "multi-discharge calibration")
return launch(cfg, driver, config_path, env=env, note=note)