From 2cab3a15a4e6c08b19a274eb16246eee215e488f Mon Sep 17 00:00:00 2001 From: BaptisteDE Date: Thu, 1 Oct 2026 15:55:00 +0200 Subject: [PATCH 1/2] Formula models, JSON candidates, training loop and control charts - FormulaModel: forward models declared as whitelisted text formulas (no eval), saved as their spec by PymcWrapper.save_model - models.py frozen; legacy_formulas maps the 23 legacy models to equivalent formulas (tested for exact equivalence) - candidates/training: candidate JSON config, fit/evaluate/compare and train_candidates loop (moved from sonergia) - control_charts: X-bar, EWMA, CUSUM stats and plots (moved from sonergia) - save_model/load_model accept str paths Co-Authored-By: Claude Opus 5.5 --- README.md | 46 +- bayesbuilding/candidates.py | 175 +++++++ bayesbuilding/control_charts.py | 567 +++++++++++++++++++++++ bayesbuilding/formula.py | 318 +++++++++++++ bayesbuilding/legacy_formulas.py | 152 ++++++ bayesbuilding/models.py | 8 + bayesbuilding/plotting.py | 98 ++++ bayesbuilding/training.py | 254 ++++++++++ bayesbuilding/wrapper.py | 26 +- tests/test_candidates_training.py | 266 +++++++++++ tests/test_control_charts.py | 266 +++++++++++ tests/test_formula.py | 184 ++++++++ tests/test_formula_legacy_equivalence.py | 140 ++++++ tests/test_models.py | 138 ++++++ 14 files changed, 2621 insertions(+), 17 deletions(-) create mode 100644 bayesbuilding/candidates.py create mode 100644 bayesbuilding/control_charts.py create mode 100644 bayesbuilding/formula.py create mode 100644 bayesbuilding/legacy_formulas.py create mode 100644 bayesbuilding/training.py create mode 100644 tests/test_candidates_training.py create mode 100644 tests/test_control_charts.py create mode 100644 tests/test_formula.py create mode 100644 tests/test_formula_legacy_equivalence.py diff --git a/README.md b/README.md index c535523..dc7d60e 100644 --- a/README.md +++ b/README.md @@ -23,10 +23,15 @@ project. - `PymcWrapper`: a thin wrapper around a PyMC model that handles prior/posterior sampling, scoring, LOO cross-validation, and saving/loading fitted models to disk. -- A library of ready-to-use `bayesbuilding.models` functions covering common building - energy patterns: seasonal change-point energy signatures, heating/cooling with - occupation change points, solar-radiation-augmented models, artificial lighting, and - PV panel production (constant efficiency and NOCT models). +- Forward models written as formulas (`bayesbuilding.formula.FormulaModel`), e.g. + `"g[occ]*max(tau - text, 0) - fs*rad + base[occ]"`: no Python function to write for + a new model variant. (`bayesbuilding.models` holds the older hand-written functions, + frozen and kept only to reload traces saved with them.) +- JSON candidate configs (`bayesbuilding.candidates`) and a training loop that fits, + scores (LOO, R2, NMBE, CV(RMSE), MAE at daily/weekly/monthly resolution) and ranks + candidates (`bayesbuilding.training.train_candidates`). +- Control charts on residuals (`bayesbuilding.control_charts`): X-bar, EWMA (with alarm + confirmation/reset) and CUSUM, with limits following the model's own sigma. - Plotting helpers (`bayesbuilding.plotting`) for prior/posterior comparison, HDI time series plots, and change-point diagnostic plots, with both `matplotlib` and `plotly` backends. @@ -47,7 +52,7 @@ import numpy as np import pandas as pd import pymc as pm -from bayesbuilding.models import season_cp_heating_es +from bayesbuilding.formula import FormulaModel from bayesbuilding.wrapper import PymcWrapper # Monthly external temperature and heating consumption @@ -57,10 +62,10 @@ data = pd.DataFrame( ) data["heating"] = 50 * np.maximum(14 - data["Text"], 0) + 50 + np.random.randn(12) * 5 -# Define the model and priors for a seasonal change-point energy signature: -# heating = g * max(tau - Text, 0) + base +# Define the model and priors for a seasonal change-point energy signature. +# Symbols that are not priors are inputs, i.e. columns of x. model = PymcWrapper( - model_function=season_cp_heating_es, + model_function=FormulaModel(mu="g*max(tau - Text, 0) + base", sigma="sigma"), priors_dict={ "g": (pm.Normal, dict(name="g", mu=40, sigma=5)), "tau": (pm.Normal, dict(name="tau", mu=12, sigma=1)), @@ -82,7 +87,28 @@ reloaded = PymcWrapper() reloaded.load_model("my_model") ``` -See `bayesbuilding/models.py` for the full list of built-in model functions and -`tests/test_wrapper.py` for a complete end-to-end example, including scoring on held-out +## Formula grammar + +- Operators `+ - * / **`, numeric constants, and the functions `max(a, b)`, + `min(a, b)`, `switch(cond, a, b)`, `sigmoid`, `sqrt`, `exp`, `log`, `abs`. +- `g[occ]` / `g[occ, 1]` index a vector/matrix prior by a categorical input. +- `sigma` is the likelihood scale and may be heteroscedastic, e.g. + `"sqrt(s0[occ]**2 + (s1*sigmoid((tau - text)/1.5))**2)"`; `lower` (default `0.0`) + is returned for a `TruncatedNormal` likelihood, `None` disables it. +- Input order (the columns of `x`): continuous drivers first, then categorical inputs, + each in order of first appearance. + +Expressions are parsed against a whitelist, never `eval`-ed. In a candidate config: + +```json +{ + "name": "dt_occ", + "model": {"mu": "g[occ]*dt - fs[occ]*rad", "sigma": "s0[occ]", "lower": 0.0}, + "inputs": {"dt": "dt__C", "rad": "GHI__W/m2", "occ": "occupation"}, + "priors": {"g": {"dist": "HalfNormal", "kwargs": {"sigma": 500, "shape": 2}}, "...": {}} +} +``` + +See `tests/test_wrapper.py` for a complete end-to-end example, including scoring on held-out data and plotting predictions with `bayesbuilding.plotting.time_series_hdi` and `changepoint_graph`. \ No newline at end of file diff --git a/bayesbuilding/candidates.py b/bayesbuilding/candidates.py new file mode 100644 index 0000000..67523df --- /dev/null +++ b/bayesbuilding/candidates.py @@ -0,0 +1,175 @@ +"""JSON configuration of candidate models to fit and compare. + +A config file looks like:: + + { + "config_id": 69370, + "target": "heating_energy", + "feature_pipe": {...}, + "candidates": [ + { + "name": "dt_occ", + "model": {"mu": "g[occ]*dt - fs[occ]*rad", "sigma": "s0[occ]", "lower": 0.0}, + "inputs": {"dt": "", "rad": "", "occ": ""}, + "priors": {"g": {"dist": "HalfNormal", "kwargs": {"sigma": 10, "shape": 2}}, ...}, + "draws": 4000, "tune": 1000 + } + ], + "unused_candidates": [] + } + +``model`` is a :class:`~bayesbuilding.formula.FormulaModel` spec. ``inputs`` maps +each input symbol of the formula to a real column name, since column names vary +from one site to the next. ``feature_pipe`` is an opaque dict left to the caller +(e.g. a Tide pipe computing the columns once for every candidate). +""" + +import json +import warnings +from dataclasses import dataclass, field +from pathlib import Path + +import pymc as pm + +from bayesbuilding.formula import FormulaModel +from bayesbuilding.legacy_formulas import legacy_to_formula + + +def build_priors_dict(priors: dict) -> dict: + """Turn a JSON-serializable priors spec into PymcWrapper's priors_dict. + + Input : {"g": {"dist": "Normal", "kwargs": {"mu": 100.0, "sigma": 10.0}}, ...} + Output : {"g": (pm.Normal, {"name": "g", "mu": 100.0, "sigma": 10.0}), ...} + """ + return { + name: (getattr(pm, spec["dist"]), {"name": name, **spec["kwargs"]}) + for name, spec in priors.items() + } + + +@dataclass +class CandidateConfig: + """One forward model to fit and compare, within a BayesConfig. + + ``model`` is given as a formula spec dict (or a FormulaModel) and is stored + as a FormulaModel bound to the prior names, so ``model.inputs`` is the + ordered list of input symbols. A legacy model name (a string) is still + accepted and translated with :func:`~bayesbuilding.legacy_formulas. + legacy_to_formula`, with a DeprecationWarning. + """ + + name: str + model: dict | FormulaModel + inputs: dict[str, str] + priors: dict + prior_draws: int = 4000 + draws: int = 4000 + tune: int = 1000 + chains: int | None = None + random_seed: int | None = None + + def __post_init__(self): + if isinstance(self.model, str): + warnings.warn( + f"Candidate {self.name!r}: legacy model name {self.model!r}, " + "replace it by its formula (see legacy_to_formula)", + DeprecationWarning, + stacklevel=3, + ) + self.model = legacy_to_formula(self.model) + if isinstance(self.model, dict): + self.model = FormulaModel.from_dict(self.model) + try: + self.model.bind(self.priors) + except ValueError as e: + raise ValueError(f"Candidate {self.name!r}: {e}") from e + + def to_dict(self) -> dict: + spec = self.model.to_dict() + spec.pop("inputs", None) # derived from the priors on load + return { + "name": self.name, + "model": spec, + "inputs": dict(self.inputs), + "priors": self.priors, + "prior_draws": self.prior_draws, + "draws": self.draws, + "tune": self.tune, + "chains": self.chains, + "random_seed": self.random_seed, + } + + +@dataclass +class BayesConfig: + """A set of candidate models sharing the same target and feature pipe.""" + + target: str + candidates: list[CandidateConfig] + config_id: int | str | None = None + feature_pipe: dict = field(default_factory=dict) + unused_candidates: list[CandidateConfig] = field(default_factory=list) + + @classmethod + def from_dict(cls, raw: dict) -> "BayesConfig": + raw = dict(raw) + candidates = [CandidateConfig(**c) for c in raw.pop("candidates")] + unused = [CandidateConfig(**c) for c in raw.pop("unused_candidates", [])] + return cls(candidates=candidates, unused_candidates=unused, **raw) + + @classmethod + def from_json(cls, path: Path | str) -> "BayesConfig": + return cls.from_dict(json.loads(Path(path).read_text(encoding="utf-8"))) + + def to_dict(self) -> dict: + return { + "config_id": self.config_id, + "target": self.target, + "feature_pipe": self.feature_pipe, + "candidates": [c.to_dict() for c in self.candidates], + "unused_candidates": [c.to_dict() for c in self.unused_candidates], + } + + def to_json(self, path: Path | str): + Path(path).write_text( + json.dumps(self.to_dict(), indent=2, ensure_ascii=False), encoding="utf-8" + ) + + def get_candidate(self, name: str) -> CandidateConfig: + for candidate in self.candidates + self.unused_candidates: + if candidate.name == name: + return candidate + raise KeyError(f"No candidate named {name!r}") + + +def posterior_as_priors(candidate: CandidateConfig, wrapper) -> dict: + """Turn a fitted wrapper's parameter posterior into a new priors spec. + + For each of ``candidate.priors``' variables, keeps the original prior's + distribution family (and any kwarg other than ``mu``/``sigma``, e.g. + ``TruncatedNormal``'s ``lower``) but replaces whichever of ``mu``/``sigma`` + the original spec already declared with the posterior sample's empirical + mean/std -- i.e. a moment-matched parametric refit of "today's posterior" + usable as "tomorrow's prior" (see :func:`~bayesbuilding.training. + fit_candidate`'s ``priors`` argument). This is a standard approximation for + sequential Bayesian updating: PyMC priors must be given as a parametric + family, so the posterior's exact (possibly non-Gaussian) shape isn't + preserved, only its first two moments. Only kwargs already present in the + original spec are touched -- a one-parameter family like ``HalfNormal`` + (kwarg ``sigma`` only, no ``mu``) doesn't get a spurious ``mu`` injected, + which would break the next ``fit_candidate(..., priors=updated)`` call. + """ + posterior = wrapper.traces["sampling"].posterior + updated = {} + for name, spec in candidate.priors.items(): + samples = posterior[name].to_numpy().reshape(-1) + moments = {} + if "mu" in spec["kwargs"]: + moments["mu"] = float(samples.mean()) + if "sigma" in spec["kwargs"]: + moments["sigma"] = float(samples.std()) + updated[name] = { + "dist": spec["dist"], + "kwargs": {**spec["kwargs"], **moments}, + } + return updated diff --git a/bayesbuilding/control_charts.py b/bayesbuilding/control_charts.py new file mode 100644 index 0000000..974e833 --- /dev/null +++ b/bayesbuilding/control_charts.py @@ -0,0 +1,567 @@ +"""Statistical process control charts on a fitted candidate's residuals. + +X-bar (Individuals), EWMA and CUSUM statistics computed from the posterior +predictive of ``mu`` and ``sigma`` (see :func:`bayesbuilding.training. +sample_mu_and_observations`), so control limits follow the model's own -- +possibly heteroscedastic -- noise scale. The ``*_control_stats`` functions take +plain numpy/pandas inputs (draws already flattened across chains) so they can be +unit-tested without running any MCMC; the ``plot_*_chart`` functions draw them +with plotly. +""" + +import numpy as np +import pandas as pd +import plotly.graph_objs as go + + +def flatten_predictive_samples(draws) -> np.ndarray: + """(chain, draw, time) or already-(samples, time), ndarray or + xarray.DataArray -> (samples, time) ndarray. + + Same positional convention as ``bayesbuilding.plotting._flatten_chains``: + never indexes by a dim name (PyMC/arviz auto-names an undeclared-``dims`` + axis like "observations_dim_0", which carries no semantic meaning), just + collapses every leading axis into one samples axis positionally. + """ + arr = np.asarray(draws) + if arr.ndim > 2: + arr = arr.reshape(-1, arr.shape[-1]) + return arr + + +def xbar_control_stats( + y_true: pd.Series, mu_draws, sigma_draws, L: float = 1.96 +) -> dict: + """Classical Individuals/X control-chart statistic and limits for the + residual ``mesure - modèle``, with the model's own (possibly + heteroscedastic) noise scale standing in for the usual historical + process sigma. + + ``mu_draws``/``sigma_draws`` are the posterior predictive of the + deterministic mean and the likelihood's noise scale (e.g. from + ``bayesbuilding.training.sample_mu_and_observations`` -- NOT the noisy + ``observations`` themselves, which would double-count noise: once in the + residual, once in the limits). Reduced across draws to one point + estimate per day: ``mu_hat = median(mu)``, ``sigma_hat = mean(sigma)``. + + ``residual[t] = y_true[t] - mu_hat[t]`` is then a single trajectory + (not one per draw) -- comparing a measurement to a *known*/calibrated + target and a *known* noise scale is exactly the classical control-chart + setup, unlike a credible band that would re-propagate the model's own + parameter uncertainty into the limits on every point. ``sigma_hat`` is + per-day (not pooled into one scalar for the whole period): for a model + with a heteroscedastic noise term (e.g. + ``bayesbuilding.models.heating_cp_occ_rad``'s ``sigma`` growing with + heating activity), the limits should genuinely narrow/widen with it + instead of being forced flat. + + Returns a dict with ``residual``, ``mu_hat``, ``sigma_hat``, ``ucl`` + (``L*sigma_hat``), ``lcl`` (``-ucl``), and ``out_of_control`` (bool array, + ``residual`` outside ``[lcl, ucl]``). + """ + mu_hat = np.median(flatten_predictive_samples(mu_draws), axis=0) + sigma_hat = np.mean(flatten_predictive_samples(sigma_draws), axis=0) + y = np.asarray(y_true, dtype=float) + residual = y - mu_hat + + ucl = L * sigma_hat + lcl = -ucl + out_of_control = (residual > ucl) | (residual < lcl) + + return { + "residual": residual, + "mu_hat": mu_hat, + "sigma_hat": sigma_hat, + "ucl": ucl, + "lcl": lcl, + "out_of_control": out_of_control, + } + + +def standardized_residuals(y_true: pd.Series, mu_draws, sigma_draws) -> np.ndarray: + """Per-day standardized residual (z-score) ``(y_true - mu_ref) / sig_ref``, + where ``mu_ref``/``sig_ref`` are the same per-day ``mu_hat``/``sigma_hat`` + point estimates computed by :func:`xbar_control_stats` -- reused rather + than recomputed, so both stay in sync -- but rescaled to a unitless z + instead of the raw ``residual``/``ucl``/``lcl`` that function returns. + + Meant as an M&V follow-up diagnostic once a candidate is calibrated on a + train period and applied forward to a verification/monitoring period: a + bimodal histogram of ``z`` signals a mixture of regimes rather than a + simple level shift; a scatter of ``z`` against an explanatory variable + (e.g. outdoor temperature) that clusters away from zero at one end + signals a localized, unmodeled effect; a run of consecutive out-of-range + days in ``z`` over time signals a dated event (drift, fault) rather than + i.i.d. noise. + """ + stats = xbar_control_stats(y_true, mu_draws, sigma_draws) + return stats["residual"] / stats["sigma_hat"] + + +def ewma_control_stats( + y_true: pd.Series, + mu_draws, + sigma_draws, + lam: float = 0.2, + L: float = 1.96, + n_confirm: int = 3, + n_reset: int = None, + reset_epsilon: float = 0.5, +) -> dict: + """Classical EWMA control-chart statistic and limits on the same + residual as :func:`xbar_control_stats` (one trajectory, not a per-draw + ensemble). + + The EWMA recursion ``Z[0] = lam*residual[0]``, ``Z[t] = lam*residual[t] + + (1-lam)*Z[t-1]`` is applied once, to the single residual trajectory. + Its limits use the standard Montgomery formula generalized to a + time-varying noise scale ``sigma_hat[t]`` (valid here because, unlike an + ensemble of posterior draws, day-to-day residuals against a *fixed* + ``mu_hat`` are independent given that estimate -- no shared per-draw + parameters correlating them across days): + ``Var(Z[t]) = lam**2 * sum_{j<=t} (1-lam)**(2*(t-j)) * sigma_hat[j]**2``, + computed via the equivalent recursion + ``Var(Z[t]) = lam**2*sigma_hat[t]**2 + (1-lam)**2*Var(Z[t-1])``. This + reproduces the classic control-chart look (limits starting tight near + the first sample and widening towards an asymptote) instead of a + per-day credible band. + + ``out_of_control`` is, by default (``n_reset=None``), the purely + instantaneous test ``ewma[t]`` outside ``[lcl[t], ucl[t]]`` -- no memory + of past points. That's a poor fit for a transient event (the error rises + then drops abruptly back to 0): the EWMA's own inertia keeps it pinned + outside the limits for several steps after the anomaly is already gone, + while a plain "N steps out -> reset" would risk instead masking a real, + sustained drift by prematurely re-labelling it "back to normal". + + Passing ``n_reset`` (not None) switches on a 3-state alternative instead, + layered on top of the same ``ewma`` trajectory: + + - NORMAL: default state. + - ALARME (confirmed): reached once ``ewma`` sits outside its limits for + ``n_confirm`` consecutive steps in a row (rules out a single noisy + point). + - Once in ALARME, the *raw* residual (not the EWMA) is watched instead: + ``|residual[t]| < reset_epsilon * sigma_hat[t]`` for ``n_reset`` + consecutive steps confirms the anomaly has genuinely cleared, and + resets the EWMA memory (``ewma[t] <- 0``, its variance recursion + restarted at ``lam**2*sigma_hat[t]**2`` as if ``t`` were a fresh first + point) before returning to NORMAL. If the residual never settles back + under the threshold, ALARME simply persists to the end of the series + -- a sustained drift never gets reset away. + + ``reset_epsilon`` is expressed in multiples of ``sigma_hat[t]`` (like + ``L`` here, and ``k``/``h`` for :func:`cusum_control_stats`), so it + stays meaningful across sites/candidates with different residual scales. + ``n_confirm`` is only used when ``n_reset`` is given. + + Returns a dict with ``ewma``, ``ucl``, ``lcl``, ``out_of_control`` (bool + array -- aliases ``alarm`` when the reset logic is active, otherwise the + instantaneous test described above), ``alarm`` (the persistent confirmed- + alarm state; equals ``out_of_control`` when ``n_reset`` is None), + ``reset_points`` (bool array, True at the exact step a reset fired -- + all False when ``n_reset`` is None), plus the underlying + ``residual``/``sigma_hat`` from :func:`xbar_control_stats`. + """ + xbar = xbar_control_stats(y_true, mu_draws, sigma_draws, L=L) + residual, sigma_hat = xbar["residual"], xbar["sigma_hat"] + n = residual.shape[0] + reset_active = n_reset is not None + + ewma = np.empty(n) + var_z = np.empty(n) + alarm = np.zeros(n, dtype=bool) + reset_points = np.zeros(n, dtype=bool) + + ewma[0] = lam * residual[0] + var_z[0] = lam**2 * sigma_hat[0] ** 2 + + state = "normal" + consec_out = 0 + consec_normal_err = 0 + + def _out_of_limits(z, v): + limit = L * np.sqrt(v) + return z > limit or z < -limit + + if reset_active: + consec_out = 1 if _out_of_limits(ewma[0], var_z[0]) else 0 + if consec_out >= n_confirm: + state = "alarm" + alarm[0] = state == "alarm" + + for t in range(1, n): + ewma[t] = lam * residual[t] + (1 - lam) * ewma[t - 1] + var_z[t] = lam**2 * sigma_hat[t] ** 2 + (1 - lam) ** 2 * var_z[t - 1] + + if reset_active: + if state == "normal": + out_t = _out_of_limits(ewma[t], var_z[t]) + consec_out = consec_out + 1 if out_t else 0 + if consec_out >= n_confirm: + state = "alarm" + consec_normal_err = 0 + else: # state == "alarm" + if abs(residual[t]) < reset_epsilon * sigma_hat[t]: + consec_normal_err += 1 + else: + consec_normal_err = 0 + if consec_normal_err >= n_reset: + state = "normal" + consec_out = 0 + ewma[t] = 0.0 + var_z[t] = lam**2 * sigma_hat[t] ** 2 + reset_points[t] = True + alarm[t] = state == "alarm" + + ucl = L * np.sqrt(var_z) + lcl = -ucl + + if reset_active: + out_of_control = alarm + else: + out_of_control = (ewma > ucl) | (ewma < lcl) + alarm = out_of_control + + return { + "ewma": ewma, + "ucl": ucl, + "lcl": lcl, + "out_of_control": out_of_control, + "alarm": alarm, + "reset_points": reset_points, + "residual": residual, + "sigma_hat": sigma_hat, + } + + +def cusum_control_stats( + y_true: pd.Series, mu_draws, sigma_draws, k: float = 0.5, h: float = 5.0, L: float = 1.96 +) -> dict: + """Classical two-sided tabular CUSUM (Page) on the standardized residual + ``z = residual / sigma_hat`` -- same ``mu_hat``/``sigma_hat`` point + estimates as :func:`xbar_control_stats`/:func:`ewma_control_stats` (one + trajectory, not a per-draw ensemble). + + ``k`` (the reference/slack value, standard default 0.5 sigma -- half the + smallest shift worth detecting) and ``h`` (the decision interval, + standard default 4-5 sigma) are both already in standardized units, so + the decision limits are the constant ``h``/``-h`` regardless of any + heteroscedasticity in ``sigma_hat`` -- unlike xbar/ewma's limits, which + track ``sigma_hat`` directly in the residual's own (unstandardized) + units. + + Running sums (reset at 0, never negative):: + + cusum_pos[t] = max(0, cusum_pos[t-1] + z[t] - k) + cusum_neg[t] = max(0, cusum_neg[t-1] - z[t] - k) + + Out of control the first time either sum exceeds ``h`` -- ``cusum_pos`` + signals a sustained upward shift, ``cusum_neg`` a sustained downward one. + A CUSUM is more sensitive than an X or EWMA chart to a small, sustained + drift (as opposed to a single large outlier), since it accumulates + evidence across days rather than resetting each step. + + Returns a dict with ``cusum_pos``, ``cusum_neg``, ``h``, + ``out_of_control`` (bool array, either sum > ``h``), plus the underlying + ``residual``/``sigma_hat`` from :func:`xbar_control_stats`. + """ + xbar = xbar_control_stats(y_true, mu_draws, sigma_draws, L=L) + residual, sigma_hat = xbar["residual"], xbar["sigma_hat"] + z = residual / sigma_hat + n = z.shape[0] + + cusum_pos = np.empty(n) + cusum_neg = np.empty(n) + prev_pos = prev_neg = 0.0 + for t in range(n): + prev_pos = max(0.0, prev_pos + z[t] - k) + prev_neg = max(0.0, prev_neg - z[t] - k) + cusum_pos[t] = prev_pos + cusum_neg[t] = prev_neg + + out_of_control = (cusum_pos > h) | (cusum_neg > h) + + return { + "cusum_pos": cusum_pos, + "cusum_neg": cusum_neg, + "h": h, + "out_of_control": out_of_control, + "residual": residual, + "sigma_hat": sigma_hat, + } + + + +def _contiguous_true_runs(mask: np.ndarray) -> list[tuple]: + """Positions ``(start_i, end_i)`` (inclusive) of each contiguous run of + ``True`` in ``mask``. Simpler than + ``bayesbuilding.plotting._boolean_blocks`` -- no single-point padding is + needed here, since a confirmed EWMA alarm run is, by construction, at + least ``n_confirm`` steps long (see + ``ewma_control_stats``). + """ + runs = [] + run_start = None + for i, val in enumerate(mask): + if val and run_start is None: + run_start = i + elif not val and run_start is not None: + runs.append((run_start, i - 1)) + run_start = None + if run_start is not None: + runs.append((run_start, len(mask) - 1)) + return runs + + +def _plot_control_chart( + index, + stat: np.ndarray, + ucl: np.ndarray, + lcl: np.ndarray, + out_of_control: np.ndarray, + y_label: str, + title: str, + series_name: str, + alarm: np.ndarray = None, + reset_points: np.ndarray = None, +) -> go.Figure: + """Shared renderer for :func:`plot_xbar_chart`/:func:`plot_ewma_chart`, + in the classical SPC style (Minitab et al.): a single trajectory (not a + per-draw ensemble/credible band), a solid center line at 0, solid red + UCL/LCL curves (flat for a homoscedastic sigma, following the model's + own heteroscedasticity otherwise -- see + ``xbar_control_stats``/``ewma_control_stats``), + and red markers directly on the trajectory's own out-of-limits points -- + the point genuinely sits outside ``[lcl, ucl]`` there, so it visibly + pokes past the red lines instead of a marker that's always contained + inside its own band by construction. + + ``alarm``/``reset_points`` are optional (only :func:`plot_ewma_chart`'s + 3-state reset logic passes them, see + ``ewma_control_stats``) -- ``alarm``'s ``True`` + runs are shaded pale red in the background (the *confirmed*, persistent + alarm period, as opposed to ``out_of_control``'s pointwise markers), and + ``reset_points`` gets its own green marker where the EWMA was reset back + to 0 after the error genuinely returned to normal. + """ + index = np.asarray(index) + fig = go.Figure() + + if alarm is not None: + for start_i, end_i in _contiguous_true_runs(alarm): + fig.add_vrect( + x0=index[start_i], + x1=index[end_i], + fillcolor="rgba(239,85,59,0.12)", + line_width=0, + layer="below", + ) + + fig.add_trace( + go.Scatter(x=index, y=ucl, mode="lines", line=dict(color="#EF553B", width=1.5)) + ) + fig.add_trace( + go.Scatter(x=index, y=lcl, mode="lines", line=dict(color="#EF553B", width=1.5)) + ) + fig.add_trace( + go.Scatter( + x=index, + y=stat, + mode="lines+markers", + line=dict(color="#2a78d6"), + marker=dict(size=5), + name=series_name, + ) + ) + if out_of_control.any(): + fig.add_trace( + go.Scatter( + x=index[out_of_control], + y=stat[out_of_control], + mode="markers", + marker=dict( + color="#EF553B", size=10, symbol="circle-open", line=dict(width=2) + ), + name="hors contrôle", + ) + ) + if reset_points is not None and reset_points.any(): + fig.add_trace( + go.Scatter( + x=index[reset_points], + y=stat[reset_points], + mode="markers", + marker=dict( + color="#2ca02c", size=10, symbol="circle-open", line=dict(width=2) + ), + name="reset EWMA", + ) + ) + + fig.add_hline(y=0, line_color="black") + fig.add_annotation( + x=index[-1], + y=ucl[-1], + xanchor="left", + yanchor="middle", + showarrow=False, + font=dict(color="#EF553B"), + text=f"LCS={ucl[-1]:.3g}", + ) + fig.add_annotation( + x=index[-1], + y=lcl[-1], + xanchor="left", + yanchor="middle", + showarrow=False, + font=dict(color="#EF553B"), + text=f"LCI={lcl[-1]:.3g}", + ) + + n_out = int(out_of_control.sum()) + fig.add_annotation( + x=0.02, + y=0.98, + xref="paper", + yref="paper", + showarrow=False, + align="left", + xanchor="left", + yanchor="top", + text=f"Points hors contrôle : {n_out}/{len(stat)} ({n_out / len(stat) * 100:.1f}%)", + ) + + fig.update_layout( + template="plotly_white", + title=title, + yaxis_title=y_label, + hovermode="x unified", + showlegend=False, + ) + return fig + + +def plot_xbar_chart( + measure_ts: pd.Series, + mu_prediction, + sigma_prediction, + L: float = 1.96, + y_label: str = None, + title: str = None, +) -> go.Figure: + """Carte de contrôle X (individus) classique sur le résidu (mesure - + modèle), centrée en 0. + + ``mu_prediction``/``sigma_prediction`` sont la posterior predictive de la + moyenne déterministe et de l'échelle de bruit du modèle (voir + ``bayesbuilding.training.sample_mu_and_observations``) -- PAS les observations + bruitées elles-mêmes, sinon le bruit serait compté deux fois (une fois + dans le résidu, une fois dans les limites). Les limites viennent du + sigma du modèle (voir ``xbar_control_stats``), + éventuellement variable dans le temps si le modèle est hétéroscédastique + -- pas une bande de crédibilité recalculée à chaque point. + """ + stats = xbar_control_stats(measure_ts, mu_prediction, sigma_prediction, L=L) + return _plot_control_chart( + measure_ts.index, + stats["residual"], + stats["ucl"], + stats["lcl"], + stats["out_of_control"], + y_label, + title or "Carte de contrôle X (résidu)", + "résidu (mesure - modèle)", + ) + + +def plot_ewma_chart( + measure_ts: pd.Series, + mu_prediction, + sigma_prediction, + lam: float = 0.2, + L: float = 1.96, + n_confirm: int = 3, + n_reset: int = None, + reset_epsilon: float = 0.5, + y_label: str = None, + title: str = None, +) -> go.Figure: + """Carte de contrôle EWMA classique sur le résidu (mesure - modèle), + centrée en 0 -- même résidu que :func:`plot_xbar_chart`, lissé par la + récursion EWMA, avec des limites qui s'élargissent depuis le premier + point jusqu'à une asymptote (voir + ``ewma_control_stats``), comme une carte EWMA + classique (Minitab et al.). + + ``n_confirm``/``n_reset``/``reset_epsilon`` sont transmis tels quels à + :func:`bayesbuilding.control_charts.ewma_control_stats` -- par défaut + (``n_reset=None``) rien ne change par rapport à la carte EWMA classique. + En passant ``n_reset``, la logique à 3 états (NORMAL / ALARME confirmée / + retour à la normale) s'active : la période d'alarme confirmée est alors + ombrée en rouge pâle en arrière-plan (au lieu des simples marqueurs + ponctuels de ``out_of_control``), et chaque reset de l'EWMA (l'erreur + brute est réellement revenue sous ``reset_epsilon`` sigma) est marqué + d'un cercle vert. + """ + stats = ewma_control_stats( + measure_ts, + mu_prediction, + sigma_prediction, + lam=lam, + L=L, + n_confirm=n_confirm, + n_reset=n_reset, + reset_epsilon=reset_epsilon, + ) + return _plot_control_chart( + measure_ts.index, + stats["ewma"], + stats["ucl"], + stats["lcl"], + stats["out_of_control"], + y_label, + title or f"Carte de contrôle EWMA (résidu, λ={lam})", + "EWMA résidu", + alarm=stats["alarm"] if n_reset is not None else None, + reset_points=stats["reset_points"] if n_reset is not None else None, + ) + + +def plot_cusum_chart( + measure_ts: pd.Series, + mu_prediction, + sigma_prediction, + k: float = 0.5, + h: float = 5.0, + L: float = 1.96, + y_label: str = None, + title: str = None, +) -> go.Figure: + """Carte de contrôle CUSUM tabulaire classique sur le résidu standardisé + (mesure - modèle)/sigma -- plus sensible qu'une carte X ou EWMA à une + petite dérive soutenue (voir + ``cusum_control_stats``). + + Tracée comme ``cusum_pos - cusum_neg`` : au plus une des deux sommes est + active à la fois en pratique (l'autre est retombée à 0 par construction), + donc cette différence signée reproduit la même trajectoire "centrée en 0, + limites symétriques" que les cartes X/EWMA, via le même renderer partagé + (:func:`_plot_control_chart`) -- ``h``/``-h`` sont des limites constantes + (déjà en unités sigma), contrairement aux limites de + :func:`plot_xbar_chart`/:func:`plot_ewma_chart` qui suivent + ``sigma_hat`` dans les unités brutes du résidu. + """ + stats = cusum_control_stats(measure_ts, mu_prediction, sigma_prediction, k=k, h=h, L=L) + signed = stats["cusum_pos"] - stats["cusum_neg"] + h_line = np.full_like(signed, stats["h"], dtype=float) + return _plot_control_chart( + measure_ts.index, + signed, + h_line, + -h_line, + stats["out_of_control"], + y_label, + title or f"Carte de contrôle CUSUM (résidu standardisé, k={k}, h={h})", + "CUSUM (C+ - C-)", + ) diff --git a/bayesbuilding/formula.py b/bayesbuilding/formula.py new file mode 100644 index 0000000..845d53d --- /dev/null +++ b/bayesbuilding/formula.py @@ -0,0 +1,318 @@ +"""Forward models declared as formulas instead of Python functions. + +A :class:`FormulaModel` is a ``model_function`` for :class:`~bayesbuilding.wrapper. +PymcWrapper` built from plain-text expressions, so that a new model variant is +a line of JSON rather than a new function in :mod:`bayesbuilding.models`:: + + FormulaModel( + mu="g[occ]*dt - fs[occ]*rad + C*(beta*DTint_dt + (1-beta)*DText_dt)", + sigma="s0[occ]", + lower=0.0, + ) + +Grammar +------- +- Operators: ``+ - * / **`` and unary ``+``/``-``. +- Numeric constants. +- Symbols: every name is either a prior (a key of the wrapper's + ``priors_dict``) or an input (a feature column). Inputs are whatever names + are not priors. Their order -- the column order of ``x`` -- is: continuous + drivers first, then categorical inputs (used as an index or as a + ``switch`` condition), each group in order of first appearance across + ``mu``, ``sigma`` and ``extras``. So ``x[:, 0]`` is the main driver, e.g. + ``dt`` in ``g[occ]*dt - fs[occ]*rad`` -> ``(dt, rad, occ)``. Pass + ``inputs`` explicitly to force another order. +- Indexing: ``g[occ]`` or ``g[occ, 1]`` indexes a vector/matrix prior by a + categorical input (cast to int) and/or integer constants. +- Functions: ``max(a, b)``, ``min(a, b)``, ``switch(cond, a, b)``, + ``sigmoid(a)``, ``sqrt(a)``, ``exp(a)``, ``log(a)``, ``abs(a)``. + +Anything else (attribute access, comparisons, arbitrary calls, lambdas...) is +rejected at construction time: expressions are parsed with :mod:`ast` and +walked against a whitelist, never passed to ``eval``. +""" + +import ast +import warnings + +import pymc as pm +import pytensor.tensor as pt + +FUNCTIONS = { + "max": (pm.math.maximum, 2), + "min": (pm.math.minimum, 2), + "switch": (pm.math.switch, 3), + "sigmoid": (pm.math.sigmoid, 1), + "sqrt": (pt.sqrt, 1), + "exp": (pt.exp, 1), + "log": (pt.log, 1), + "abs": (pt.abs, 1), +} + +_BINARY_OPERATORS = { + ast.Add: lambda a, b: a + b, + ast.Sub: lambda a, b: a - b, + ast.Mult: lambda a, b: a * b, + ast.Div: lambda a, b: a / b, + ast.Pow: lambda a, b: a**b, +} + +_UNARY_OPERATORS = { + ast.USub: lambda a: -a, + ast.UAdd: lambda a: a, +} + + +class FormulaError(ValueError): + """Raised when a formula uses a construct outside the allowed grammar.""" + + +def _parse(expr: str) -> ast.expr: + try: + tree = ast.parse(expr.strip(), mode="eval").body + except SyntaxError as e: + raise FormulaError(f"Invalid formula {expr!r}: {e.msg}") from e + _check(tree, expr) + return tree + + +def _check(node: ast.AST, expr: str): + """Walk ``node`` and raise FormulaError on anything outside the grammar.""" + if isinstance(node, ast.BinOp): + if type(node.op) not in _BINARY_OPERATORS: + raise FormulaError(f"Operator {type(node.op).__name__} not allowed in {expr!r}") + _check(node.left, expr) + _check(node.right, expr) + elif isinstance(node, ast.UnaryOp): + if type(node.op) not in _UNARY_OPERATORS: + raise FormulaError(f"Operator {type(node.op).__name__} not allowed in {expr!r}") + _check(node.operand, expr) + elif isinstance(node, ast.Constant): + if isinstance(node.value, bool) or not isinstance(node.value, (int, float)): + raise FormulaError(f"Only numeric constants allowed in {expr!r}") + elif isinstance(node, ast.Name): + if node.id in FUNCTIONS: + raise FormulaError(f"Function {node.id!r} used as a symbol in {expr!r}") + elif isinstance(node, ast.Subscript): + if not isinstance(node.value, ast.Name): + raise FormulaError(f"Only a symbol can be indexed in {expr!r}") + indices = node.slice.elts if isinstance(node.slice, ast.Tuple) else [node.slice] + for index in indices: + is_int = isinstance(index, ast.Constant) and type(index.value) is int + if not (isinstance(index, ast.Name) or is_int): + raise FormulaError( + f"Index must be an input name or an integer in {expr!r}" + ) + elif isinstance(node, ast.Call): + if not isinstance(node.func, ast.Name) or node.func.id not in FUNCTIONS: + name = getattr(node.func, "id", ast.unparse(node.func)) + raise FormulaError( + f"Function {name!r} not allowed in {expr!r}. " + f"Allowed: {sorted(FUNCTIONS)}" + ) + if node.keywords: + raise FormulaError(f"Keyword arguments not allowed in {expr!r}") + arity = FUNCTIONS[node.func.id][1] + if len(node.args) != arity: + raise FormulaError( + f"{node.func.id}() takes {arity} argument(s), " + f"got {len(node.args)} in {expr!r}" + ) + for arg in node.args: + _check(arg, expr) + else: + raise FormulaError(f"{type(node).__name__} not allowed in {expr!r}") + + +def _symbols(node: ast.AST) -> list[str]: + """Names in order of first appearance (functions excluded).""" + # ast.walk is breadth-first: sort by source position for textual order. + found = sorted( + ( + child + for child in ast.walk(node) + if isinstance(child, ast.Name) and child.id not in FUNCTIONS + ), + key=lambda child: (child.lineno, child.col_offset), + ) + return list(dict.fromkeys(child.id for child in found)) + + +def _index_symbols(node: ast.AST) -> set[str]: + """Names used inside a subscript: categorical inputs, cast to int.""" + names = set() + for child in ast.walk(node): + if isinstance(child, ast.Subscript): + indices = ( + child.slice.elts if isinstance(child.slice, ast.Tuple) else [child.slice] + ) + names |= {i.id for i in indices if isinstance(i, ast.Name)} + return names + + +def _switch_condition_symbols(node: ast.AST) -> set[str]: + """Names used directly as a ``switch`` condition (a flag input).""" + return { + child.args[0].id + for child in ast.walk(node) + if isinstance(child, ast.Call) + and child.func.id == "switch" + and isinstance(child.args[0], ast.Name) + } + + +def _evaluate(node: ast.AST, env: dict, categorical: dict): + if isinstance(node, ast.BinOp): + return _BINARY_OPERATORS[type(node.op)]( + _evaluate(node.left, env, categorical), + _evaluate(node.right, env, categorical), + ) + if isinstance(node, ast.UnaryOp): + return _UNARY_OPERATORS[type(node.op)](_evaluate(node.operand, env, categorical)) + if isinstance(node, ast.Constant): + return node.value + if isinstance(node, ast.Name): + return env[node.id] + if isinstance(node, ast.Subscript): + indices = node.slice.elts if isinstance(node.slice, ast.Tuple) else [node.slice] + resolved = tuple( + categorical[i.id] if isinstance(i, ast.Name) else i.value for i in indices + ) + target = env[node.value.id] + return target[resolved[0]] if len(resolved) == 1 else target[resolved] + if isinstance(node, ast.Call): + func = FUNCTIONS[node.func.id][0] + return func(*(_evaluate(arg, env, categorical) for arg in node.args)) + raise FormulaError(f"{type(node).__name__} not allowed") # unreachable + + +class FormulaModel: + """A PymcWrapper ``model_function`` defined by text formulas. + + :param mu: expression of the likelihood mean. + :param sigma: expression of the likelihood scale. + :param lower: lower truncation bound returned as the ``"lower"`` extra + (for a ``TruncatedNormal`` likelihood); ``None`` returns no ``"lower"``. + :param extras: further named expressions exposed as ``pm.Deterministic`` + for diagnostics. + :param inputs: ordered input names, i.e. the expected column order of + ``x``. Usually left out and derived with :meth:`bind` from the prior + names; when given, every other symbol is a prior. + """ + + def __init__( + self, + mu: str, + sigma: str, + lower: float | None = 0.0, + extras: dict[str, str] | None = None, + inputs: list[str] | tuple[str, ...] | None = None, + ): + self.mu = mu + self.sigma = sigma + self.lower = lower + self.extras = dict(extras or {}) + reserved = {"mu", "sigma", "lower"} & set(self.extras) + if reserved: + raise FormulaError(f"Reserved extras names: {sorted(reserved)}") + self._trees = {"mu": _parse(mu), "sigma": _parse(sigma)} + self._trees.update({name: _parse(e) for name, e in self.extras.items()}) + + self.symbols: list[str] = [] + for tree in self._trees.values(): + for name in _symbols(tree): + if name not in self.symbols: + self.symbols.append(name) + self.categorical_inputs = set().union( + *(_index_symbols(t) for t in self._trees.values()) + ) + self._flag_inputs = self.categorical_inputs | set().union( + *(_switch_condition_symbols(t) for t in self._trees.values()) + ) + self.inputs = None if inputs is None else tuple(inputs) + if self.inputs is not None: + self._check_inputs(self.inputs) + + def _check_inputs(self, inputs): + unknown = [name for name in inputs if name not in self.symbols] + if unknown: + raise FormulaError(f"Inputs {unknown} do not appear in the formula") + not_inputs = self.categorical_inputs - set(inputs) + if not_inputs: + raise FormulaError( + f"{sorted(not_inputs)} are used as indices, so they must be inputs" + ) + + def input_names(self, prior_names) -> tuple[str, ...]: + """Symbols that are not priors: continuous drivers first, then + categorical/flag inputs, each in order of first appearance.""" + names = [name for name in self.symbols if name not in set(prior_names)] + return tuple( + [n for n in names if n not in self._flag_inputs] + + [n for n in names if n in self._flag_inputs] + ) + + def bind(self, prior_names) -> "FormulaModel": + """Fix ``inputs`` from the prior names (keeping an explicit ``inputs`` + order if one was given); warns if a prior is not used by the formula. + Returns ``self`` for chaining.""" + prior_names = set(prior_names) + unused = sorted(prior_names - set(self.symbols)) + if unused: + # Not an error: the prior is still sampled, it just doesn't inform + # anything (legacy model functions silently ignored such priors). + warnings.warn(f"Priors {unused} are not used in the formula", stacklevel=2) + clash = sorted(prior_names & set(self.extras)) + if clash: + raise FormulaError(f"Extras {clash} have the same name as a prior") + inputs = self.input_names(prior_names) + if self.inputs is not None: + if set(self.inputs) != set(inputs): + raise FormulaError( + f"Explicit inputs {list(self.inputs)} do not match the " + f"non-prior symbols {list(inputs)}" + ) + inputs = self.inputs + self._check_inputs(inputs) + self.inputs = inputs + return self + + def __call__(self, x, variables: dict): + inputs = ( + self.inputs if self.inputs is not None else self.input_names(variables) + ) + missing = [ + name for name in self.symbols if name not in inputs and name not in variables + ] + if missing: + raise FormulaError(f"Symbols {missing} are neither priors nor inputs") + env = dict(variables) + env.update({name: x[:, i] for i, name in enumerate(inputs)}) + categorical = { + name: x[:, inputs.index(name)].astype("int64") + for name in self.categorical_inputs + } + values = { + name: _evaluate(tree, env, categorical) for name, tree in self._trees.items() + } + mu = values.pop("mu") + out = {"sigma": values.pop("sigma")} + if self.lower is not None: + out["lower"] = pt.constant(float(self.lower)) + out.update(values) + return mu, out + + def to_dict(self) -> dict: + spec = {"mu": self.mu, "sigma": self.sigma, "lower": self.lower} + if self.extras: + spec["extras"] = dict(self.extras) + if self.inputs is not None: + spec["inputs"] = list(self.inputs) + return spec + + @classmethod + def from_dict(cls, spec: dict) -> "FormulaModel": + return cls(**spec) + + def __repr__(self): + return f"FormulaModel(mu={self.mu!r}, sigma={self.sigma!r}, lower={self.lower!r})" diff --git a/bayesbuilding/legacy_formulas.py b/bayesbuilding/legacy_formulas.py new file mode 100644 index 0000000..497a389 --- /dev/null +++ b/bayesbuilding/legacy_formulas.py @@ -0,0 +1,152 @@ +"""Formula equivalents of the legacy functions of :mod:`bayesbuilding.models`. + +Frozen translation table, used to migrate configs that name a legacy model +(``"model": "heating_dt_occ_rad"``) to a formula (``"model": {"mu": ...}``), +see :func:`legacy_to_formula`. ``tests/test_formula_legacy_equivalence.py`` +checks each formula evaluates exactly like its legacy function. Do not add +entries: new models are written directly as formulas. + +``LEGACY_INPUTS`` gives the column order each legacy function indexes via +``x[:, i]``. +""" + +_SETBACK_RAD = "switch(is_heating, base[0] + {g}*max({tau} - text, 0) - fs*rad, base[1])" +_DT_CORE = "g[occ]*dt - fs[occ]*rad - alpha[occ]*elec_consumption" +_WALL = "C*(beta*DTint_dt + (1 - beta)*DText_dt)" + +LEGACY_FORMULAS: dict[str, dict] = { + "season_cp_heating_es_dt": {"mu": "g*max(dt - tau, 0) + base", "sigma": "sigma"}, + "season_cp_heating_es": {"mu": "g*max(tau - text, 0) + base", "sigma": "sigma"}, + "heating_es_dju": {"mu": "g*dju + base", "sigma": "sigma"}, + "season_cp_occ_cp_heating_es": { + "mu": "g[occ]*max(tau[occ] - text, 0) + baseline[occ]", + "sigma": "sigma[occ]", + }, + "season_cp_occ_cp_es_dt": { + "mu": "g[occ]*max(dt - tau[occ], 0) + baseline[occ]", + "sigma": "sigma[occ]", + }, + "heating_es_dju_rad": {"mu": "base + g*dju - fs*rad", "sigma": "sigma"}, + "season_cp_heating_es_rad": { + "mu": _SETBACK_RAD.format(g="g", tau="tau"), + "sigma": "sigma", + }, + "season_cp_heating_es_setback": { + "mu": "switch(is_heating, base[0] + g*max(tau - text, 0), base[1])", + "sigma": "sigma", + }, + "heating_es_dju_rad_occ": { + "mu": "base[occ] + g[occ]*dju - fs[occ]*rad", + "sigma": "sigma", + }, + "heating_es_dju_rad_occ_setback": { + "mu": "switch(occ, base[0] + g*dju - fs*rad, base[1])", + "sigma": "sigma", + }, + "heating_cp_occ_rad": { + "mu": "base[occ] + max(g[occ]*(tau[occ] - text) - fs[occ]*rad, 0)", + "sigma": "sqrt(s0[occ]**2 + (s1*sigmoid((tau[occ] - text)/1.5))**2)", + }, + "season_cp_heating_es_rad_g_by_period": { + "mu": _SETBACK_RAD.format(g="g[period]", tau="tau"), + "sigma": "sigma", + }, + "season_cp_heating_es_rad_tau_by_period": { + "mu": _SETBACK_RAD.format(g="g", tau="tau[period]"), + "sigma": "sigma", + }, + "season_cp_heating_es_rad_base_by_period": { + "mu": ( + "switch(is_heating, base[is_heating, period] + g*max(tau - text, 0)" + " - fs*rad, base[is_heating, period])" + ), + "sigma": "sigma", + }, + "season_cp_heating_es_rad_g_tau_by_period": { + "mu": _SETBACK_RAD.format(g="g[period]", tau="tau[period]"), + "sigma": "sigma", + }, + "heating_dt_occ_rad": {"mu": _DT_CORE, "sigma": "s0[occ]"}, + "heating_dt_occ_rad_multiroom": { + "mu": ( + "g[occ, 0]*dt_chambre1 + g[occ, 1]*dt_cuisine + g[occ, 2]*dt_sejour" + " - (fs[occ, 0] + fs[occ, 1] + fs[occ, 2])*rad" + " - alpha[occ]*elec_consumption" + ), + "sigma": "s0[occ]", + }, + "heating_dt_occ_rad_lag": { + "mu": f"{_DT_CORE} + h[occ]*dt_lag", + "sigma": "s0[occ]", + }, + "heating_dt_occ_rad_DTdt": {"mu": f"{_DT_CORE} + c*DT_dt", "sigma": "s0[occ]"}, + "heating_dt_occ_rad_DTdt_lags": { + "mu": ( + f"{_DT_CORE} + c0*DT_dt + c0*rho*DT_dt_lag1" + " + c0*rho**2*DT_dt_lag2 + c0*rho**3*DT_dt_lag3" + ), + "sigma": "s0[occ]", + }, + "heating_dt_occ_rad_DTdt_wall": { + "mu": f"{_DT_CORE} + {_WALL}", + "sigma": "s0[occ]", + }, + "heating_dt_occ_rad_DTdt_wall_Ci": { + "mu": f"{_DT_CORE} + {_WALL} + Ci*DTint_intraday", + "sigma": "s0[occ]", + }, + "heating_dt_occ_rad_DTdt_wall_Ci_radlag": { + "mu": f"{_DT_CORE} + {_WALL} + Ci*DTint_intraday - fm[occ]*rad_lag1", + "sigma": "s0[occ]", + }, +} + +_DT_INPUTS = ("dt", "rad", "elec_consumption", "occ") + +LEGACY_INPUTS: dict[str, tuple[str, ...]] = { + "season_cp_heating_es_dt": ("dt",), + "season_cp_heating_es": ("text",), + "heating_es_dju": ("dju",), + "season_cp_occ_cp_heating_es": ("occ", "text"), + "season_cp_occ_cp_es_dt": ("occ", "dt"), + "heating_es_dju_rad": ("dju", "rad"), + "season_cp_heating_es_rad": ("text", "rad", "is_heating"), + "season_cp_heating_es_setback": ("text", "is_heating"), + "heating_es_dju_rad_occ": ("dju", "rad", "occ"), + "heating_es_dju_rad_occ_setback": ("dju", "rad", "occ"), + "heating_cp_occ_rad": ("text", "rad", "occ"), + "season_cp_heating_es_rad_g_by_period": ("text", "rad", "is_heating", "period"), + "season_cp_heating_es_rad_tau_by_period": ("text", "rad", "is_heating", "period"), + "season_cp_heating_es_rad_base_by_period": ("text", "rad", "is_heating", "period"), + "season_cp_heating_es_rad_g_tau_by_period": ("text", "rad", "is_heating", "period"), + "heating_dt_occ_rad": _DT_INPUTS, + "heating_dt_occ_rad_multiroom": ( + "dt_chambre1", + "dt_cuisine", + "dt_sejour", + "rad", + "elec_consumption", + "occ", + ), + "heating_dt_occ_rad_lag": _DT_INPUTS + ("dt_lag",), + "heating_dt_occ_rad_DTdt": _DT_INPUTS + ("DT_dt",), + "heating_dt_occ_rad_DTdt_lags": _DT_INPUTS + + ("DT_dt", "DT_dt_lag1", "DT_dt_lag2", "DT_dt_lag3"), + "heating_dt_occ_rad_DTdt_wall": _DT_INPUTS + ("DTint_dt", "DText_dt"), + "heating_dt_occ_rad_DTdt_wall_Ci": _DT_INPUTS + + ("DTint_dt", "DText_dt", "DTint_intraday"), + "heating_dt_occ_rad_DTdt_wall_Ci_radlag": _DT_INPUTS + + ("DTint_dt", "DText_dt", "DTint_intraday", "rad_lag1"), +} + + +def legacy_to_formula(name: str) -> dict: + """Formula spec (``{"mu", "sigma", "lower"}``) equivalent to legacy model + ``name``.""" + try: + return {**LEGACY_FORMULAS[name], "lower": 0.0} + except KeyError: + raise KeyError( + f"No formula equivalent for legacy model {name!r}. " + f"Available: {list(LEGACY_FORMULAS)}" + ) diff --git a/bayesbuilding/models.py b/bayesbuilding/models.py index 229363c..c3d8d09 100644 --- a/bayesbuilding/models.py +++ b/bayesbuilding/models.py @@ -1,3 +1,11 @@ +"""Legacy forward-model functions -- frozen, do not add new models here. + +New models are declared as formulas (see :mod:`bayesbuilding.formula`), e.g. in +a candidate's JSON config. These functions are kept so that traces saved with +``PymcWrapper.save_model`` before formulas existed (which reference a function +of this module by name) can still be reloaded. +""" + import pymc as pm import pytensor.tensor as pt diff --git a/bayesbuilding/plotting.py b/bayesbuilding/plotting.py index 5f297c5..158c748 100644 --- a/bayesbuilding/plotting.py +++ b/bayesbuilding/plotting.py @@ -1480,3 +1480,101 @@ def compare_bars( plt.savefig(image_path, format="png", bbox_inches="tight") return plt.gcf() + + +def plot_candidate_metric_bars( + comparison: pd.DataFrame, + metric: str, + std_col: str = None, + n_std: float = 2.0, + ascending: bool = False, + thresholds: float | tuple[float, float] = None, + threshold_color: str = "#EF553B", + y_label: str = None, + title: str = None, + bar_color: str = "#2a78d6", +) -> go.Figure: + """Bar chart ranking candidates by one ``comparison_train`` metric + (descending by default), with symmetric error bars at ``+/- n_std`` times + the metric's own std/se column. + + ``comparison`` is the DataFrame returned by + :func:`bayesbuilding.training.compare_candidates` -- candidate name as index, one + column per metric (``elpd_loo``, ``elpd_loo_se``, ``r2_mean_d``, + ``r2_sd_d``, ...) -- **not** its ``.reset_index().T.reset_index()`` + transposed display form. + + ``std_col`` defaults to the naming convention + :func:`bayesbuilding.training.evaluate_candidate` produces: ``"{x}_mean_{suffix}"`` + -> ``"{x}_sd_{suffix}"`` (e.g. ``"r2_mean_d"`` -> ``"r2_sd_d"``), falling + back to ``"{metric}_se"``/``"{metric}_sd"`` for metrics with no ``_mean_`` + in their name (e.g. ``"elpd_loo"`` -> ``"elpd_loo_se"``). Pass it + explicitly to override. + + ``thresholds`` draws optional dashed reference lines (e.g. ASHRAE + Guideline 14-style M&V acceptance bounds): a single number (e.g. ``15`` + for a monthly CV(RMSE) <= 15% cap) draws one line; a pair (e.g. + ``(-5, 5)`` for a monthly NMBE band) draws both lines and shades the + region between them as the "passing" zone. No assumption is made about + which side of a single threshold is acceptable -- only a two-value pair + gets a shaded band. + + Single series (one bar per candidate, same metric) -- one hue, no legend, + per the "identity via color" rule not applying here (bars are distinguished + by their x position/label, not by color). + """ + if metric not in comparison.columns: + raise KeyError(f"{metric!r} not in comparison columns: {list(comparison.columns)}") + + if std_col is None: + candidate_cols = ( + [metric.replace("_mean_", "_sd_")] if "_mean_" in metric else [] + ) + [f"{metric}_se", f"{metric}_sd"] + std_col = next((c for c in candidate_cols if c in comparison.columns), None) + if std_col is None: + raise KeyError( + f"Could not infer a std/se column for {metric!r} (tried " + f"{candidate_cols}); pass std_col explicitly." + ) + elif std_col not in comparison.columns: + raise KeyError(f"{std_col!r} not in comparison columns: {list(comparison.columns)}") + + ranked = comparison[[metric, std_col]].sort_values(metric, ascending=ascending) + + fig = go.Figure( + go.Bar( + x=ranked.index.astype(str), + y=ranked[metric], + error_y=dict(type="data", array=n_std * ranked[std_col], visible=True), + marker_color=bar_color, + ) + ) + + if thresholds is not None: + threshold_values = ( + [thresholds] if isinstance(thresholds, (int, float)) else list(thresholds) + ) + if len(threshold_values) == 2: + low, high = sorted(threshold_values) + fig.add_hrect( + y0=low, y1=high, fillcolor="rgba(44,160,44,0.08)", line_width=0, layer="below" + ) + for value in threshold_values: + fig.add_hline( + y=value, + line_dash="dash", + line_color=threshold_color, + annotation_text=f"{value:g}", + annotation_position="right", + ) + + fig.update_layout( + template="plotly_white", + title=title or f"{metric} par modèle (± {n_std:g}σ)", + xaxis_title="Modèle", + yaxis_title=y_label or metric, + showlegend=False, + ) + return fig + + diff --git a/bayesbuilding/training.py b/bayesbuilding/training.py new file mode 100644 index 0000000..dcf02e8 --- /dev/null +++ b/bayesbuilding/training.py @@ -0,0 +1,254 @@ +"""Fit, score and compare the candidates of a :class:`~bayesbuilding.candidates. +BayesConfig`.""" + +import sys +from pathlib import Path + +import arviz as az +import numpy as np +import pandas as pd +import pymc as pm +import xarray as xr + +from bayesbuilding.candidates import BayesConfig, CandidateConfig, build_priors_dict +from bayesbuilding.wrapper import PymcWrapper, r2_score, resample_samples + + +def nmbe(y_true, y_pred) -> float: + """Normalized Mean Bias Error [%] (same definition as corrai's).""" + return float(np.sum(y_pred - y_true) / np.sum(y_true) * 100) + + +def cv_rmse(y_true, y_pred) -> float: + """Coefficient of Variation of the RMSE [%] (same definition as corrai's, + i.e. ASHRAE Guideline 14's, with n - 1 degrees of freedom).""" + n = len(y_true) + return float(np.sqrt(np.sum((y_true - y_pred) ** 2) / (n - 1)) / np.mean(y_true) * 100) + + +def mae(y_true, y_pred) -> float: + """Mean Absolute Error.""" + return float(np.mean(np.abs(y_true - y_pred))) + + +METRICS = {"r2": r2_score, "nmbe": nmbe, "cv_rmse": cv_rmse, "mae": mae} + + +def build_candidate_xy( + df: pd.DataFrame, target: str, candidate: CandidateConfig +) -> tuple[pd.DataFrame, pd.Series]: + """Select the columns of ``df`` a candidate needs, in its formula's input order. + + Raises KeyError early (naming the candidate and the missing keys) if + ``candidate.inputs`` doesn't map every input symbol of the formula, instead + of letting PyMC fail opaquely further down. + """ + required = candidate.model.inputs + missing = [k for k in required if k not in candidate.inputs] + if missing: + raise KeyError( + f"Candidate {candidate.name!r}: missing input mapping for {missing}" + ) + x = df[[candidate.inputs[k] for k in required]] + y = df[target] + return x, y + + +def fit_candidate( + candidate: CandidateConfig, + x: pd.DataFrame, + y: pd.Series, + artifact_path: Path | None, + priors: dict = None, +) -> PymcWrapper: + """Build, sample (prior + posterior) and persist a PymcWrapper for one candidate. + + The likelihood is a Normal truncated at ``candidate.model.lower`` (an energy + consumption cannot be negative), or a plain Normal when ``lower`` is None. + + ``priors`` defaults to ``candidate.priors`` (the config's declared priors). + Passing a different priors spec (e.g. from :func:`~bayesbuilding.candidates. + posterior_as_priors`) lets the same candidate be re-conditioned on a + different period, starting from an already-informed prior instead of the + config's original one -- a sequential Bayesian update. + + The wrapper is saved to ``artifact_path / "trace"`` unless ``artifact_path`` + is None. + """ + if candidate.model.lower is None: + likelihood, likelihood_params = pm.Normal, ["sigma"] + else: + likelihood, likelihood_params = pm.TruncatedNormal, ["sigma", "lower"] + wrapper = PymcWrapper( + model_function=candidate.model, + priors_dict=build_priors_dict( + priors if priors is not None else candidate.priors + ), + likelihood=likelihood, + likelihood_params=likelihood_params, + ) + wrapper.sample_prior( + samples=candidate.prior_draws, + x=x, + sample_kwargs={"random_seed": candidate.random_seed}, + ) + sample_kwargs = {"random_seed": candidate.random_seed} + if sys.platform == "win32": + # PyMC's multi-chain sampling multiprocesses via "spawn" on Windows, + # which re-imports the caller's top-level code in each child process -- + # this crashes scripts written to run cell-by-cell without an + # `if __name__ == "__main__":` guard. Run chains sequentially instead. + sample_kwargs["cores"] = 1 + wrapper.sample( + x=x, + y=y, + draws=candidate.draws, + tune=candidate.tune, + chains=candidate.chains, + sample_kwargs=sample_kwargs, + ) + if artifact_path is not None: + wrapper.save_model(Path(artifact_path) / "trace") + return wrapper + + +def sample_mu_and_observations( + wrapper: PymcWrapper, x: pd.DataFrame +) -> az.InferenceData: + """Posterior predictive of the deterministic ``mu``, the likelihood's noise + scale ``sigma`` and the noisy ``observations`` for arbitrary ``x``, in one + call. + + ``PymcWrapper.sample_posterior_predictive`` only samples ``observations``; + ``mu``/``sigma`` are needed to separate parametric from noise uncertainty + and to derive control-chart limits from the model's own (possibly + heteroscedastic) noise scale (see :mod:`bayesbuilding.control_charts`). + Does not mutate ``wrapper.traces["posterior"]``. + """ + with wrapper.model: + pm.set_data({"x": x}) + return pm.sample_posterior_predictive( + trace=wrapper.traces["sampling"], var_names=["mu", "sigma", "observations"] + ) + + +def flatten_regime_vars( + group: xr.Dataset, var_names: list[str] +) -> dict[str, np.ndarray]: + """Flatten each named variable of a posterior/prior xarray group to a 1D + array of draws, ready for e.g. a KDE. + + A variable indexed by a two-level regime -- shape 2 on a ``"{var}_dim_0"`` + dim, e.g. ``g[occ]`` -- is exploded into ``"{var}0"``/``"{var}1"``. A + ``var_names`` entry absent from ``group`` (e.g. a candidate whose priors + don't cover every var asked for) is silently skipped rather than raising. + """ + flattened = {} + for name in var_names: + if name not in group.data_vars: + continue + da = group[name] + dim = f"{name}_dim_0" + if dim in da.dims and da.sizes[dim] == 2: + for i in (0, 1): + flattened[f"{name}{i}"] = da.isel({dim: i}).to_numpy().reshape(-1) + else: + flattened[name] = da.to_numpy().reshape(-1) + return flattened + + +def evaluate_candidate( + wrapper: PymcWrapper, x_test: pd.DataFrame, y_test: pd.Series +) -> dict: + """LOO (elpd, se) and scores of a fitted candidate on ``(x_test, y_test)``. + + R2/NMBE/CV(RMSE)/MAE are computed both at native (daily) resolution and + aggregated to calendar weeks/months ("_d"/"_W"/"_ME" suffixes), as the mean + and std over posterior predictive sample paths. Weekly/monthly scores sum + each daily sample path (and y_test) into bins before scoring, rather than + resampling the driving variable and re-evaluating the model -- a nonlinear + forward model (e.g. max(driving_var - tau, 0)) isn't equivalent under the + two orderings. LOO is intrinsically pointwise (daily) and computed on the + training data. + + The posterior predictive is sampled once and reused for every resolution + and metric. + """ + loo = wrapper.get_loo_score() + results = { + "elpd_loo": float(loo.elpd), + "elpd_loo_se": float(loo.se), + } + + wrapper.sample_posterior_predictive(x=x_test) + post_trace = wrapper.traces["posterior"].posterior_predictive["observations"] + flattened_trace_full = np.array(post_trace).reshape(-1, post_trace.shape[-1]) + + for suffix, resample_rule in (("d", None), ("W", "W"), ("ME", "ME")): + if resample_rule is None: + flattened_trace, y_eval = flattened_trace_full, y_test + else: + flattened_trace, y_eval = resample_samples( + flattened_trace_full, y_test, resample_rule + ) + y_eval = np.asarray(y_eval) + for name, metric_func in METRICS.items(): + scores_array = np.array( + [metric_func(y_eval, sample) for sample in flattened_trace] + ) + results[f"{name}_mean_{suffix}"] = float(np.mean(scores_array)) + results[f"{name}_sd_{suffix}"] = float(np.std(scores_array)) + + return results + + +def compare_candidates(results: dict[str, dict]) -> pd.DataFrame: + """Assemble a comparison table across candidates, best elpd_loo first.""" + return pd.DataFrame(results).T.sort_values("elpd_loo", ascending=False) + + +def train_candidates( + config: BayesConfig, + df: pd.DataFrame, + out_dir: Path | str, + df_test: pd.DataFrame | None = None, + candidates: list[CandidateConfig] | None = None, + comparison_filename: str | None = "comparison_train.csv", +) -> tuple[dict[str, PymcWrapper], pd.DataFrame, dict[str, str]]: + """Fit and score every candidate, then build the comparison table. + + Each candidate is fitted on ``df`` and saved to ``out_dir//trace``, + then scored on ``df_test`` (in sample on ``df`` when None). A candidate + that fails (sampler crash, missing input column...) is recorded in the + returned errors and skipped instead of aborting the loop. + + :param candidates: subset to train, defaults to ``config.candidates``. + :param comparison_filename: file written in ``out_dir`` with the comparison + table, or None to write nothing. + :return: ``(wrappers, comparison, errors)``: fitted wrappers by name, the + :func:`compare_candidates` table, and error messages by name. + """ + out_dir = Path(out_dir) + candidates = config.candidates if candidates is None else candidates + wrappers, results, errors = {}, {}, {} + for candidate in candidates: + print(f"=== candidate: {candidate.name} ===") + try: + x, y = build_candidate_xy(df, config.target, candidate) + artifact_path = out_dir / candidate.name + artifact_path.mkdir(parents=True, exist_ok=True) + wrapper = fit_candidate(candidate, x, y, artifact_path) + if df_test is None: + x_eval, y_eval = x, y + else: + x_eval, y_eval = build_candidate_xy(df_test, config.target, candidate) + results[candidate.name] = evaluate_candidate(wrapper, x_eval, y_eval) + wrappers[candidate.name] = wrapper + except Exception as exc: + print(f" [failed: {exc}]") + errors[candidate.name] = str(exc) + + comparison = compare_candidates(results) if results else pd.DataFrame() + if comparison_filename is not None and not comparison.empty: + comparison.to_csv(out_dir / comparison_filename) + return wrappers, comparison, errors diff --git a/bayesbuilding/wrapper.py b/bayesbuilding/wrapper.py index e0a5596..1045e3d 100644 --- a/bayesbuilding/wrapper.py +++ b/bayesbuilding/wrapper.py @@ -3,6 +3,7 @@ from collections.abc import Callable from pathlib import Path import bayesbuilding.models as mods +from bayesbuilding.formula import FormulaModel import warnings import arviz as az @@ -144,9 +145,10 @@ class PymcWrapper: load_model(Path): Load the model traces and json file contained in the desired path and build the pymc model. - If the loaded model function is not in bayesbuilding.model, it must - be provided separately and the pymc model must be build using the build_model() - method. + A FormulaModel is rebuilt from its saved spec. A legacy model function + is looked up by name in bayesbuilding.models; if it is not there, it + must be provided separately and the pymc model must be build using the + build_model() method. plot_dist_comparison(var_names): Plot compare prior and posterior distributions of variables and observations var_names arguments filter the values to display. default is self.var_names @@ -203,7 +205,8 @@ def __repr__(self): return string_out - def save_model(self, dir_path: Path): + def save_model(self, dir_path: Path | str): + dir_path = Path(dir_path) if not os.path.exists(dir_path): os.makedirs(dir_path) @@ -211,8 +214,13 @@ def save_model(self, dir_path: Path): traces.to_netcdf((dir_path / f"{name}.nc").as_posix(), engine="h5netcdf") with open(dir_path / "config.json", "w", encoding="utf-8") as f: + # A FormulaModel is saved as its spec (rebuilt on load); a legacy + # model function is saved by name and looked up in + # bayesbuilding.models on load. + is_formula = hasattr(self.model_function, "to_dict") to_dump = { - "model_function": self.model_function, + "model_function": None if is_formula else self.model_function, + "model_spec": self.model_function.to_dict() if is_formula else None, "priors_dict": self.priors_dict, "likelihood": self.likelihood, "likelihood_params": self.likelihood_params, @@ -223,7 +231,8 @@ def save_model(self, dir_path: Path): to_dump, f, ensure_ascii=False, default=custom_serializer, indent=4 ) - def load_model(self, dir_path: Path): + def load_model(self, dir_path: Path | str): + dir_path = Path(dir_path) if not os.path.exists(dir_path): raise ValueError(f"Provided dir_path : {dir_path} not found") @@ -247,7 +256,10 @@ def load_model(self, dir_path: Path): if isinstance(self.likelihood, str): self.likelihood = getattr(pm, self.likelihood) - if self.model_function is not None: + model_spec = config_dict.get("model_spec") + if model_spec is not None: + self.model_function = FormulaModel.from_dict(model_spec) + elif self.model_function is not None: try: self.model_function = getattr(mods, self.model_function) except AttributeError: diff --git a/tests/test_candidates_training.py b/tests/test_candidates_training.py new file mode 100644 index 0000000..124ea66 --- /dev/null +++ b/tests/test_candidates_training.py @@ -0,0 +1,266 @@ +import json +import warnings + +import numpy as np +import pandas as pd +import pymc as pm +import pytest +import xarray as xr + +from bayesbuilding.candidates import ( + BayesConfig, + CandidateConfig, + build_priors_dict, + posterior_as_priors, +) +from bayesbuilding.formula import FormulaModel +from bayesbuilding.legacy_formulas import LEGACY_FORMULAS, legacy_to_formula +from bayesbuilding.training import ( + build_candidate_xy, + cv_rmse, + flatten_regime_vars, + mae, + nmbe, + train_candidates, +) +from bayesbuilding.wrapper import PymcWrapper +from tests.test_formula_legacy_equivalence import PRIOR_SHAPES + +PRIORS_DT = { + "g": {"dist": "Normal", "kwargs": {"mu": 100.0, "sigma": 10.0}}, + "tau": {"dist": "Normal", "kwargs": {"mu": 0.0, "sigma": 1.0}}, + "base": {"dist": "Normal", "kwargs": {"mu": 10.0, "sigma": 1.0}}, + "sigma": {"dist": "HalfNormal", "kwargs": {"sigma": 1.0}}, +} + +CONFIG_PAYLOAD = { + "config_id": 1, + "feature_pipe": {"resample": [["Resample", {"rule": "d"}]]}, + "target": "chauffage__Wh__Zone__FULL", + "candidates": [ + { + "name": "dt_change_point", + "model": {"mu": "g*max(dt - tau, 0) + base", "sigma": "sigma", "lower": 0.0}, + "inputs": {"dt": "dt__C__Zone__FULL"}, + "priors": PRIORS_DT, + "prior_draws": 500, + "draws": 1000, + "tune": 500, + "chains": 2, + "random_seed": 42, + }, + ], + "unused_candidates": [ + { + "name": "text_change_point", + "model": {"mu": "g*max(tau - text, 0) + base", "sigma": "sigma"}, + "inputs": {"text": "Text__C"}, + "priors": PRIORS_DT, + }, + ], +} + + +class TestConfig: + def test_from_json_round_trip(self, tmp_path): + path = tmp_path / "bayes_config.json" + path.write_text(json.dumps(CONFIG_PAYLOAD)) + + config = BayesConfig.from_json(path) + + assert config.config_id == 1 + dt_candidate = config.candidates[0] + assert isinstance(dt_candidate.model, FormulaModel) + assert dt_candidate.model.inputs == ("dt",) + assert dt_candidate.draws == 1000 + assert config.unused_candidates[0].draws == 4000 # default + assert config.get_candidate("text_change_point").model.inputs == ("text",) + + config.to_json(tmp_path / "copy.json") + assert BayesConfig.from_json(tmp_path / "copy.json").to_dict() == config.to_dict() + + def test_unused_candidates_is_optional(self): + payload = {k: v for k, v in CONFIG_PAYLOAD.items() if k != "unused_candidates"} + assert BayesConfig.from_dict(payload).unused_candidates == [] + + def test_legacy_model_name_is_translated_with_a_warning(self): + with pytest.warns(DeprecationWarning, match="legacy model name"): + candidate = CandidateConfig( + name="legacy", + model="season_cp_heating_es_dt", + inputs={"dt": "dt__C"}, + priors=PRIORS_DT, + ) + assert candidate.model.mu == LEGACY_FORMULAS["season_cp_heating_es_dt"]["mu"] + assert candidate.model.lower == 0.0 + + def test_candidate_errors_name_the_candidate(self): + with pytest.raises(ValueError, match="'bad'.*indices"): + CandidateConfig( + name="bad", + model={"mu": "g[k]*dt", "sigma": "sigma"}, + inputs={}, + priors={"g": {}, "sigma": {}, "k": {}}, + ) + + def test_build_priors_dict_injects_name_and_resolves_distribution(self): + priors = build_priors_dict(PRIORS_DT) + assert priors["g"][0] is pm.Normal + assert priors["g"][1] == {"name": "g", "mu": 100.0, "sigma": 10.0} + + +class TestBuildCandidateXY: + def test_selects_and_orders_columns_driver_first(self): + df = pd.DataFrame( + { + "target": [1.0, 2.0], + "occ_col": [0.0, 1.0], + "dt_col": [3.0, 4.0], + "rad_col": [5.0, 6.0], + "unrelated": [7.0, 8.0], + } + ) + candidate = CandidateConfig( + name="dt_occ", + model={"mu": "g[occ]*dt - fs*rad", "sigma": "s0[occ]"}, + inputs={"occ": "occ_col", "dt": "dt_col", "rad": "rad_col"}, + priors={"g": {}, "fs": {}, "s0": {}}, + ) + x, y = build_candidate_xy(df, "target", candidate) + assert list(x.columns) == ["dt_col", "rad_col", "occ_col"] + assert y.tolist() == [1.0, 2.0] + + def test_raises_on_missing_input_mapping(self): + df = pd.DataFrame({"target": [1.0], "dt_col": [2.0]}) + candidate = CandidateConfig( + name="incomplete", + model={"mu": "g*max(dt - tau, 0) + base", "sigma": "sigma"}, + inputs={}, + priors=PRIORS_DT, + ) + with pytest.raises(KeyError, match="incomplete"): + build_candidate_xy(df, "target", candidate) + + +@pytest.mark.parametrize("name", sorted(LEGACY_FORMULAS)) +def test_every_legacy_formula_builds_a_truncated_normal_wrapper(name): + """What fit_candidate builds, minus the (expensive) sampling.""" + priors = { + var: {"dist": "HalfNormal", "kwargs": {"sigma": 1.0, "shape": shape or None}} + for var, shape in PRIOR_SHAPES[name].items() + } + for spec in priors.values(): + if spec["kwargs"]["shape"] is None: + del spec["kwargs"]["shape"] + candidate = CandidateConfig( + name=name, model=legacy_to_formula(name), inputs={}, priors=priors + ) + wrapper = PymcWrapper( + model_function=candidate.model, + priors_dict=build_priors_dict(priors), + likelihood=pm.TruncatedNormal, + likelihood_params=["sigma", "lower"], + ) + assert "mu" in {rv.name for rv in wrapper.model.deterministics} + + +def test_posterior_as_priors_does_not_inject_mu_for_half_normal(): + """A HalfNormal (only a 'sigma' kwarg) must not come back with a spurious + 'mu' -- pm.HalfNormal doesn't accept one.""" + candidate = CandidateConfig( + name="c", + model={"mu": "g*dt", "sigma": "s0"}, + inputs={}, + priors={ + "g": {"dist": "Normal", "kwargs": {"mu": 100.0, "sigma": 10.0}}, + "s0": {"dist": "HalfNormal", "kwargs": {"sigma": 3000.0}}, + }, + ) + fake_posterior = xr.Dataset( + { + "g": (("chain", "draw"), [[95.0, 105.0, 100.0]]), + "s0": (("chain", "draw"), [[1.0, 2.0, 3.0]]), + } + ) + + class FakeTrace: + posterior = fake_posterior + + class FakeWrapper: + traces = {"sampling": FakeTrace()} + + updated = posterior_as_priors(candidate, FakeWrapper()) + + assert set(updated["g"]["kwargs"]) == {"mu", "sigma"} + assert set(updated["s0"]["kwargs"]) == {"sigma"} + + +def test_flatten_regime_vars_explodes_shape_2_regime_variable(): + group = xr.Dataset( + { + "sigma": (("chain", "draw"), [[1.0, 2.0, 3.0]]), + "g": (("chain", "draw", "g_dim_0"), [[[10.0, 20.0], [11.0, 21.0], [12.0, 22.0]]]), + } + ) + flattened = flatten_regime_vars(group, ["sigma", "g", "does_not_exist"]) + assert set(flattened) == {"sigma", "g0", "g1"} + np.testing.assert_allclose(flattened["g1"], [20.0, 21.0, 22.0]) + + +def test_metrics_match_corrai_definitions(): + # Values from corrai.base.metrics docstrings (y_pred, y_true order there). + y_true = np.array([100.0, 200.0, 300.0]) + y_pred = np.array([110.0, 190.0, 310.0]) + assert nmbe(y_true, y_pred) == pytest.approx(1.6666666666666667) + assert cv_rmse(y_true, y_pred) == pytest.approx(6.123724356957945) + assert mae(y_true, y_pred) == pytest.approx(10.0) + + +def test_train_candidates_fits_scores_and_records_failures(tmp_path): + rng = np.random.default_rng(1) + idx = pd.date_range("2024-01-01", periods=70, freq="D") + dt = rng.uniform(0, 15, 70) + df = pd.DataFrame( + {"dt__C": dt, "energy": 4.0 * np.maximum(dt - 3, 0) + 10 + rng.normal(0, 1, 70)}, + index=idx, + ) + priors = { + "g": {"dist": "HalfNormal", "kwargs": {"sigma": 10.0}}, + "tau": {"dist": "Normal", "kwargs": {"mu": 3.0, "sigma": 2.0}}, + "base": {"dist": "Normal", "kwargs": {"mu": 10.0, "sigma": 5.0}}, + "sigma": {"dist": "HalfNormal", "kwargs": {"sigma": 3.0}}, + } + small = {"prior_draws": 20, "draws": 60, "tune": 60, "chains": 1, "random_seed": 0} + config = BayesConfig.from_dict( + { + "target": "energy", + "candidates": [ + { + "name": "cp", + "model": {"mu": "g*max(dt - tau, 0) + base", "sigma": "sigma"}, + "inputs": {"dt": "dt__C"}, + "priors": priors, + **small, + }, + { + "name": "broken", + "model": {"mu": "g*max(dt - tau, 0) + base", "sigma": "sigma"}, + "inputs": {"dt": "missing_column"}, + "priors": priors, + **small, + }, + ], + } + ) + + with warnings.catch_warnings(): + warnings.simplefilter("ignore") + wrappers, comparison, errors = train_candidates(config, df, tmp_path) + + assert list(wrappers) == ["cp"] + assert list(errors) == ["broken"] + assert list(comparison.index) == ["cp"] + for col in ("elpd_loo", "elpd_loo_se", "r2_mean_d", "cv_rmse_sd_ME", "mae_mean_W"): + assert col in comparison.columns + assert (tmp_path / "comparison_train.csv").exists() + assert (tmp_path / "cp" / "trace" / "config.json").exists() diff --git a/tests/test_control_charts.py b/tests/test_control_charts.py new file mode 100644 index 0000000..6051b9f --- /dev/null +++ b/tests/test_control_charts.py @@ -0,0 +1,266 @@ +import numpy as np +import pandas as pd +import pytest + +from bayesbuilding.control_charts import ( + _contiguous_true_runs, + cusum_control_stats, + ewma_control_stats, + flatten_predictive_samples, + plot_ewma_chart, + standardized_residuals, + xbar_control_stats, +) + + +def test_flatten_predictive_samples_collapses_chain_and_draw(): + arr_3d = np.arange(2 * 3 * 4).reshape(2, 3, 4) # (chain, draw, time) + flat = flatten_predictive_samples(arr_3d) + assert flat.shape == (6, 4) + np.testing.assert_array_equal(flat, arr_3d.reshape(6, 4)) + + arr_2d = np.arange(5 * 4).reshape(5, 4) # already (samples, time) + flat_2d = flatten_predictive_samples(arr_2d) + assert flat_2d.shape == (5, 4) + np.testing.assert_array_equal(flat_2d, arr_2d) + + +def test_xbar_control_stats_residual_and_out_of_control_flag(): + mu_draws = np.tile([100.0, 200.0, 300.0], (4, 1)) # 4 draws, same mu each day + sigma_draws = np.full((4, 3), 2.0) + y_true = np.array([100.0, 210.0, 300.0]) # day 1: residual=10, others: 0 + + stats = xbar_control_stats(y_true, mu_draws, sigma_draws, L=1.96) + + np.testing.assert_allclose(stats["mu_hat"], [100.0, 200.0, 300.0]) + np.testing.assert_allclose(stats["sigma_hat"], [2.0, 2.0, 2.0]) + np.testing.assert_allclose(stats["residual"], [0.0, 10.0, 0.0]) + np.testing.assert_allclose(stats["ucl"], [1.96 * 2.0] * 3) + np.testing.assert_allclose(stats["lcl"], [-1.96 * 2.0] * 3) + np.testing.assert_array_equal(stats["out_of_control"], [False, True, False]) + + +def test_xbar_control_stats_limits_follow_heteroscedastic_sigma(): + # A model with a heteroscedastic noise term (e.g. heating_cp_occ_rad) + # should produce limits that widen/narrow with it, not a single flat + # sigma pooled over the whole period. + n_draws, n_days = 50, 3 + mu_draws = np.zeros((n_draws, n_days)) + sigma_draws = np.tile([1.0, 5.0, 1.0], (n_draws, 1)) + y_true = np.zeros(n_days) + + stats = xbar_control_stats(y_true, mu_draws, sigma_draws) + + assert stats["ucl"][1] > stats["ucl"][0] + assert stats["ucl"][1] > stats["ucl"][2] + np.testing.assert_allclose(stats["ucl"][0], stats["ucl"][2]) + + +def test_ewma_control_stats_matches_manual_recursion(): + mu_draws = np.zeros((1, 3)) + sigma_draws = np.ones((1, 3)) + y_true = np.array([2.0, 0.0, 0.0]) # residual = [2, 0, 0], sigma_hat = 1 + lam = 0.5 + + stats = ewma_control_stats(y_true, mu_draws, sigma_draws, lam=lam, L=1.0) + + z0 = lam * 2.0 + var0 = lam**2 * 1.0**2 + z1 = lam * 0.0 + (1 - lam) * z0 + var1 = lam**2 * 1.0**2 + (1 - lam) ** 2 * var0 + z2 = lam * 0.0 + (1 - lam) * z1 + var2 = lam**2 * 1.0**2 + (1 - lam) ** 2 * var1 + + np.testing.assert_allclose(stats["ewma"], [z0, z1, z2]) + np.testing.assert_allclose(stats["ucl"], np.sqrt([var0, var1, var2])) + np.testing.assert_allclose(stats["lcl"], -np.sqrt([var0, var1, var2])) + + +def test_ewma_control_stats_limits_widen_then_plateau(): + # Classic EWMA control-limit shape: tight near the first sample, widening + # towards the constant-sigma asymptote L*sigma*sqrt(lam/(2-lam)). + n_days = 200 + mu_draws = np.zeros((10, n_days)) + sigma_draws = np.full((10, n_days), 3.0) + y_true = np.zeros(n_days) + lam, L = 0.2, 1.96 + + stats = ewma_control_stats(y_true, mu_draws, sigma_draws, lam=lam, L=L) + + asymptotic_ucl = L * 3.0 * np.sqrt(lam / (2 - lam)) + assert stats["ucl"][0] < stats["ucl"][-1] + assert stats["ucl"][-1] == pytest.approx(asymptotic_ucl, rel=1e-3) + + +def test_ewma_control_stats_default_alarm_matches_out_of_control(): + # n_reset=None (default) must reproduce the old, memory-less behaviour + # exactly: alarm == the instantaneous out_of_control test, no resets ever. + mu_draws = np.zeros((1, 5)) + sigma_draws = np.ones((1, 5)) + y_true = np.array([0.0, 3.0, 0.0, -3.0, 0.0]) + + stats = ewma_control_stats(y_true, mu_draws, sigma_draws, lam=0.5, L=1.0) + + np.testing.assert_array_equal(stats["alarm"], stats["out_of_control"]) + assert not stats["reset_points"].any() + + +def test_ewma_control_stats_confirms_alarm_then_resets_after_transient_event(): + # residual jumps to 5 for 3 steps (enough to confirm an alarm after + # n_confirm=2 consecutive out-of-limit EWMA steps), then drops back to 0 + # for good -- enough consecutive in-control-error steps (n_reset=2) to + # confirm the event has genuinely cleared and reset the EWMA to 0. + mu_draws = np.zeros((1, 8)) + sigma_draws = np.ones((1, 8)) + y_true = np.array([0.0, 5.0, 5.0, 5.0, 0.0, 0.0, 0.0, 0.0]) + lam, L = 0.5, 1.0 + + stats = ewma_control_stats( + y_true, mu_draws, sigma_draws, lam=lam, L=L, + n_confirm=2, n_reset=2, reset_epsilon=0.5, + ) + + expected_alarm = [False, False, True, True, True, False, False, False] + np.testing.assert_array_equal(stats["alarm"], expected_alarm) + np.testing.assert_array_equal(stats["out_of_control"], expected_alarm) + + expected_reset = [False, False, False, False, False, True, False, False] + np.testing.assert_array_equal(stats["reset_points"], expected_reset) + + # at the reset step, the EWMA and its variance recursion restart exactly + # as they would at a fresh first point (z <- 0, var <- lam**2*sigma**2). + assert stats["ewma"][5] == 0.0 + assert stats["ucl"][5] == pytest.approx(lam * L) + + +def test_ewma_control_stats_persistent_drift_never_resets(): + # same alarm-triggering ramp, but the residual stays elevated forever + # afterwards -- a real sustained drift must never get reset away. + mu_draws = np.zeros((1, 8)) + sigma_draws = np.ones((1, 8)) + y_true = np.array([0.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0]) + + stats = ewma_control_stats( + y_true, mu_draws, sigma_draws, lam=0.5, L=1.0, + n_confirm=2, n_reset=2, reset_epsilon=0.5, + ) + + assert stats["alarm"][2:].all() + assert not stats["reset_points"].any() + + +def test_ewma_control_stats_short_excess_never_confirms_alarm(): + # n_confirm larger than the whole series: no run of consecutive + # out-of-limit EWMA steps can ever reach it, so the alarm never confirms + # regardless of the spike -- a short blip shouldn't be over-eagerly + # flagged as a confirmed alarm. + mu_draws = np.zeros((1, 5)) + sigma_draws = np.ones((1, 5)) + y_true = np.array([0.0, 5.0, 0.0, 0.0, 0.0]) + + stats = ewma_control_stats( + y_true, mu_draws, sigma_draws, lam=0.5, L=1.0, + n_confirm=10, n_reset=2, reset_epsilon=0.5, + ) + + assert not stats["alarm"].any() + assert not stats["out_of_control"].any() + assert not stats["reset_points"].any() + + +def test_cusum_control_stats_accumulates_a_sustained_upward_drift(): + # z = residual/sigma_hat = 1.0 every day, k=0.5 -> cusum_pos grows by 0.5 + # each day, cusum_neg stays clamped at 0 (never goes negative). + mu_draws = np.zeros((1, 4)) + sigma_draws = np.ones((1, 4)) + y_true = np.array([1.0, 1.0, 1.0, 1.0]) + + stats = cusum_control_stats(y_true, mu_draws, sigma_draws, k=0.5, h=5.0) + + np.testing.assert_allclose(stats["cusum_pos"], [0.5, 1.0, 1.5, 2.0]) + np.testing.assert_allclose(stats["cusum_neg"], [0.0, 0.0, 0.0, 0.0]) + assert not stats["out_of_control"].any() + + +def test_cusum_control_stats_accumulates_a_sustained_downward_drift(): + mu_draws = np.zeros((1, 4)) + sigma_draws = np.ones((1, 4)) + y_true = np.array([-1.0, -1.0, -1.0, -1.0]) + + stats = cusum_control_stats(y_true, mu_draws, sigma_draws, k=0.5, h=5.0) + + np.testing.assert_allclose(stats["cusum_pos"], [0.0, 0.0, 0.0, 0.0]) + np.testing.assert_allclose(stats["cusum_neg"], [0.5, 1.0, 1.5, 2.0]) + assert not stats["out_of_control"].any() + + +def test_cusum_control_stats_flags_out_of_control_once_h_is_exceeded(): + # z=2.0, k=0.5 -> cusum_pos grows by 1.5/day: 1.5, 3.0, 4.5, 6.0 -> crosses h=5 on day 4. + mu_draws = np.zeros((1, 4)) + sigma_draws = np.ones((1, 4)) + y_true = np.array([2.0, 2.0, 2.0, 2.0]) + + stats = cusum_control_stats(y_true, mu_draws, sigma_draws, k=0.5, h=5.0) + + np.testing.assert_array_equal(stats["out_of_control"], [False, False, False, True]) + + +def test_standardized_residuals_matches_residual_over_sigma_hat(): + mu_draws = np.tile([100.0, 200.0, 300.0], (4, 1)) # 4 draws, same mu each day + sigma_draws = np.full((4, 3), 2.0) + y_true = np.array([100.0, 210.0, 294.0]) # residual = [0, 10, -6] + + z = standardized_residuals(y_true, mu_draws, sigma_draws) + + np.testing.assert_allclose(z, [0.0, 5.0, -3.0]) + + +def test_standardized_residuals_follows_heteroscedastic_sigma(): + # Same absolute residual, but a day with a wider sigma should read as a + # smaller z -- sigma_hat isn't pooled into a single flat scale. + mu_draws = np.zeros((10, 2)) + sigma_draws = np.tile([1.0, 5.0], (10, 1)) + y_true = np.array([2.0, 2.0]) + + z = standardized_residuals(y_true, mu_draws, sigma_draws) + + assert z[0] == pytest.approx(2.0) + assert z[1] == pytest.approx(0.4) + + +def test_contiguous_true_runs_finds_each_block(): + mask = np.array([False, True, True, False, False, True, False]) + assert _contiguous_true_runs(mask) == [(1, 2), (5, 5)] + + +def test_plot_ewma_chart_without_reset_args_has_no_alarm_shading_or_markers(): + index = pd.date_range("2024-01-01", periods=5, freq="D") + measure = pd.Series([0.0, 5.0, 5.0, 5.0, 0.0], index=index) + mu_draws = np.zeros((1, 5)) + sigma_draws = np.ones((1, 5)) + + fig = plot_ewma_chart(measure, mu_draws, sigma_draws, lam=0.5, L=1.0) + + assert len(fig.layout.shapes) == 1 # just the y=0 center line, no alarm vrect + assert not any(trace.name == "reset EWMA" for trace in fig.data) + + +def test_plot_ewma_chart_with_reset_args_shades_alarm_and_marks_resets(): + index = pd.date_range("2024-01-01", periods=8, freq="D") + measure = pd.Series( + [0.0, 5.0, 5.0, 5.0, 0.0, 0.0, 0.0, 0.0], index=index + ) + mu_draws = np.zeros((1, 8)) + sigma_draws = np.ones((1, 8)) + + fig = plot_ewma_chart( + measure, mu_draws, sigma_draws, lam=0.5, L=1.0, + n_confirm=2, n_reset=2, reset_epsilon=0.5, + ) + + # the y=0 center line, plus one shaded vrect for the confirmed alarm run + assert len(fig.layout.shapes) == 2 + reset_trace = next(trace for trace in fig.data if trace.name == "reset EWMA") + assert list(reset_trace.x) == [index[5]] + + diff --git a/tests/test_formula.py b/tests/test_formula.py new file mode 100644 index 0000000..b058589 --- /dev/null +++ b/tests/test_formula.py @@ -0,0 +1,184 @@ +import json +import tempfile +from pathlib import Path + +import numpy as np +import pandas as pd +import pymc as pm +import pytensor.tensor as pt +import pytest + +from bayesbuilding.formula import FormulaError, FormulaModel +from bayesbuilding.wrapper import PymcWrapper + + +class TestParsing: + @pytest.mark.parametrize( + "expr", + [ + "__import__('os').system('echo')", + "g.__class__", + "eval('1')", + "(lambda: 1)()", + "g if dt else tau", + "g < dt", + "g[dt*2]", + "g[0.5]", + "max(g)", + "max(g, dt, tau)", + "max(a=g, b=dt)", + "'text'", + "[g, dt]", + "g @ dt", + "max", + "g +", + ], + ) + def test_rejects_constructs_outside_the_grammar(self, expr): + with pytest.raises(FormulaError): + FormulaModel(mu=expr, sigma="sigma") + + def test_rejects_reserved_extras(self): + with pytest.raises(FormulaError, match="Reserved"): + FormulaModel(mu="g*dt", sigma="sigma", extras={"sigma": "g"}) + + def test_symbols_in_order_of_first_appearance(self): + model = FormulaModel( + mu="g[occ]*dt - fs*rad + max(tau - dt, 0)", + sigma="s0[occ]", + extras={"heat": "g[occ]*dt + rad2"}, + ) + assert model.symbols == ["g", "occ", "dt", "fs", "rad", "tau", "s0", "rad2"] + assert model.categorical_inputs == {"occ"} + + def test_bind_derives_inputs_from_prior_names(self): + model = FormulaModel(mu="g[occ]*dt - fs*rad", sigma="s0[occ]") + model.bind(["g", "fs", "s0"]) + assert model.inputs == ("dt", "rad", "occ") + + def test_switch_condition_goes_after_drivers(self): + model = FormulaModel( + mu="switch(is_heating, base[0] + g*max(tau - text, 0) - fs*rad, base[1])", + sigma="sigma", + ).bind(["base", "g", "tau", "fs", "sigma"]) + assert model.inputs == ("text", "rad", "is_heating") + assert model.categorical_inputs == set() + + def test_bind_keeps_explicit_input_order(self): + model = FormulaModel(mu="g*dt + h*rad", sigma="s", inputs=["rad", "dt"]) + assert model.bind(["g", "h", "s"]).inputs == ("rad", "dt") + + def test_bind_rejects_explicit_inputs_missing_a_symbol(self): + model = FormulaModel(mu="g*dt + h*rad", sigma="s", inputs=["dt"]) + with pytest.raises(FormulaError, match="do not match"): + model.bind(["g", "h", "s"]) + + def test_multi_index_categorical(self): + model = FormulaModel(mu="g[occ, 1]*dt + base[is_heating, period]", sigma="s") + model.bind(["g", "base", "s"]) + assert model.categorical_inputs == {"occ", "is_heating", "period"} + + def test_bind_warns_on_unused_prior(self): + with pytest.warns(UserWarning, match="not used"): + FormulaModel(mu="g*dt", sigma="sigma").bind(["g", "sigma", "tau"]) + + def test_bind_rejects_index_that_is_a_prior(self): + with pytest.raises(FormulaError, match="indices"): + FormulaModel(mu="g[k]*dt", sigma="sigma").bind(["g", "k", "sigma"]) + + def test_explicit_inputs_must_appear(self): + with pytest.raises(FormulaError, match="do not appear"): + FormulaModel(mu="g*dt", sigma="sigma", inputs=["dt", "rad"]) + + def test_dict_round_trip(self): + model = FormulaModel( + mu="g*dt", sigma="sigma", lower=None, extras={"heat": "g*dt"} + ).bind(["g", "sigma"]) + spec = json.loads(json.dumps(model.to_dict())) + clone = FormulaModel.from_dict(spec) + assert clone.to_dict() == model.to_dict() + assert clone.inputs == ("dt",) + + +class TestEvaluation: + def test_evaluates_indexing_and_functions(self): + x = np.array([[10.0, 0.0], [2.0, 1.0]]) # columns: text, occ + variables = { + "g": pt.as_tensor_variable(np.array([1.0, 3.0])), + "tau": pt.as_tensor_variable(np.array([15.0, 5.0])), + "s0": pt.as_tensor_variable(np.array([0.5, 2.0])), + } + model = FormulaModel( + mu="g[occ]*max(tau[occ] - text, 0)", + sigma="s0[occ]", + extras={"heat": "sqrt(abs(-4.0))"}, + ).bind(variables) + mu, extras = model(pt.as_tensor_variable(x), variables) + np.testing.assert_allclose(mu.eval(), [5.0, 9.0]) + np.testing.assert_allclose(extras["sigma"].eval(), [0.5, 2.0]) + assert float(extras["lower"].eval()) == 0.0 + assert float(extras["heat"].eval()) == 2.0 + + def test_no_lower_when_none(self): + model = FormulaModel(mu="g*dt", sigma="sigma", lower=None) + variables = {"g": pt.constant(1.0), "sigma": pt.constant(1.0)} + _, extras = model(pt.as_tensor_variable(np.ones((3, 1))), variables) + assert "lower" not in extras + + def test_unbound_model_derives_inputs_at_call(self): + model = FormulaModel(mu="g*dt + base", sigma="sigma") + variables = {k: pt.constant(2.0) for k in ("g", "base", "sigma")} + mu, _ = model(pt.as_tensor_variable(np.array([[1.0], [3.0]])), variables) + np.testing.assert_allclose(mu.eval(), [4.0, 8.0]) + + +@pytest.fixture +def toy_data(): + rng = np.random.default_rng(42) + idx = pd.date_range("2024-01-01", periods=60, freq="D") + text = rng.uniform(-5, 20, 60) + occ = rng.integers(0, 2, 60).astype(float) + y = 3.0 * np.maximum(15 - text, 0) + np.where(occ == 1, 20, 10) + y = y + rng.normal(0, 1, 60) + x = pd.DataFrame({"text": text, "occ": occ}, index=idx) + return x, pd.Series(np.maximum(y, 0.1), index=idx, name="energy") + + +def test_wrapper_save_load_round_trip_with_formula(toy_data): + x, y = toy_data + priors = { + "g": (pm.HalfNormal, {"name": "g", "sigma": 5}), + "tau": (pm.Normal, {"name": "tau", "mu": 15, "sigma": 3}), + "base": (pm.Normal, {"name": "base", "mu": 15, "sigma": 10, "shape": 2}), + "sigma": (pm.HalfNormal, {"name": "sigma", "sigma": 3}), + } + model = FormulaModel( + mu="g*max(tau - text, 0) + base[occ]", sigma="sigma" + ).bind(priors) + wrapper = PymcWrapper( + model_function=model, + priors_dict=priors, + likelihood=pm.TruncatedNormal, + likelihood_params=["sigma", "lower"], + ) + wrapper.sample( + x=x[list(model.inputs)], + y=y, + draws=50, + tune=50, + chains=1, + sample_kwargs={"random_seed": 0, "cores": 1, "progressbar": False}, + ) + path = Path(tempfile.mkdtemp()) / "trace" + wrapper.save_model(path) + + config = json.loads((path / "config.json").read_text(encoding="utf-8")) + assert config["model_function"] is None + assert config["model_spec"]["mu"] == "g*max(tau - text, 0) + base[occ]" + + loaded = PymcWrapper() + loaded.load_model(path) + assert isinstance(loaded.model_function, FormulaModel) + assert loaded.model_function.inputs == ("text", "occ") # driver first + loaded.sample_posterior_predictive(x=x[list(model.inputs)]) + assert "observations" in loaded.traces["posterior"].posterior_predictive diff --git a/tests/test_formula_legacy_equivalence.py b/tests/test_formula_legacy_equivalence.py new file mode 100644 index 0000000..3458337 --- /dev/null +++ b/tests/test_formula_legacy_equivalence.py @@ -0,0 +1,140 @@ +import numpy as np +import pandas as pd +import pytensor.tensor as pt +import pytest + +import bayesbuilding.models as mods +from bayesbuilding.formula import FormulaModel +from bayesbuilding.legacy_formulas import ( + LEGACY_FORMULAS, + LEGACY_INPUTS, + legacy_to_formula, +) + +N = 50 +N_OCC, N_PERIOD, N_ROOM = 2, 3, 3 + +CATEGORICAL_LEVELS = {"occ": N_OCC, "is_heating": 2, "period": N_PERIOD} + +# Shape of each prior, per legacy model (scalar when absent). +PRIOR_SHAPES = { + "season_cp_heating_es_dt": {"g": (), "tau": (), "base": (), "sigma": ()}, + "season_cp_heating_es": {"g": (), "tau": (), "base": (), "sigma": ()}, + "heating_es_dju": {"g": (), "base": (), "sigma": ()}, + "season_cp_occ_cp_heating_es": { + k: (N_OCC,) for k in ("g", "tau", "baseline", "sigma") + }, + "season_cp_occ_cp_es_dt": {k: (N_OCC,) for k in ("g", "tau", "baseline", "sigma")}, + "heating_es_dju_rad": {"g": (), "fs": (), "base": (), "sigma": ()}, + "season_cp_heating_es_rad": { + "g": (), + "fs": (), + "tau": (), + "base": (2,), + "sigma": (), + }, + "season_cp_heating_es_setback": {"g": (), "tau": (), "base": (2,), "sigma": ()}, + "heating_es_dju_rad_occ": { + "g": (N_OCC,), + "fs": (N_OCC,), + "base": (N_OCC,), + "sigma": (), + }, + "heating_es_dju_rad_occ_setback": {"g": (), "fs": (), "base": (2,), "sigma": ()}, + "heating_cp_occ_rad": { + **{k: (N_OCC,) for k in ("base", "g", "tau", "fs", "s0")}, + "s1": (), + }, + "season_cp_heating_es_rad_g_by_period": { + "g": (N_PERIOD,), + "fs": (), + "tau": (), + "base": (2,), + "sigma": (), + }, + "season_cp_heating_es_rad_tau_by_period": { + "g": (), + "fs": (), + "tau": (N_PERIOD,), + "base": (2,), + "sigma": (), + }, + "season_cp_heating_es_rad_base_by_period": { + "g": (), + "fs": (), + "tau": (), + "base": (2, N_PERIOD), + "sigma": (), + }, + "season_cp_heating_es_rad_g_tau_by_period": { + "g": (N_PERIOD,), + "fs": (), + "tau": (N_PERIOD,), + "base": (2,), + "sigma": (), + }, +} +_DT = {k: (N_OCC,) for k in ("g", "fs", "alpha", "s0")} +PRIOR_SHAPES.update( + { + "heating_dt_occ_rad": _DT, + "heating_dt_occ_rad_multiroom": { + "g": (N_OCC, N_ROOM), + "fs": (N_OCC, N_ROOM), + "alpha": (N_OCC,), + "s0": (N_OCC,), + }, + "heating_dt_occ_rad_lag": {**_DT, "h": (N_OCC,)}, + "heating_dt_occ_rad_DTdt": {**_DT, "c": ()}, + "heating_dt_occ_rad_DTdt_lags": {**_DT, "c0": (), "rho": ()}, + "heating_dt_occ_rad_DTdt_wall": {**_DT, "C": (), "beta": ()}, + "heating_dt_occ_rad_DTdt_wall_Ci": {**_DT, "C": (), "beta": (), "Ci": ()}, + "heating_dt_occ_rad_DTdt_wall_Ci_radlag": { + **_DT, + "C": (), + "beta": (), + "Ci": (), + "fm": (N_OCC,), + }, + } +) + + +def test_tables_cover_the_same_models(): + assert set(LEGACY_FORMULAS) == set(LEGACY_INPUTS) == set(PRIOR_SHAPES) + + +@pytest.mark.parametrize("name", sorted(LEGACY_FORMULAS)) +def test_formula_matches_legacy_model(name): + rng = np.random.default_rng(0) + columns = {} + for col in LEGACY_INPUTS[name]: + if col in CATEGORICAL_LEVELS: + columns[col] = rng.integers(0, CATEGORICAL_LEVELS[col], N).astype(float) + else: + columns[col] = rng.normal(5.0, 3.0, N) + df = pd.DataFrame(columns) + variables = { + k: pt.as_tensor_variable(rng.uniform(0.1, 2.0, shape)) + for k, shape in PRIOR_SHAPES[name].items() + } + + legacy_mu, legacy_extras = getattr(mods, name)( + pt.as_tensor_variable(df[list(LEGACY_INPUTS[name])].to_numpy()), variables + ) + + formula = FormulaModel.from_dict(legacy_to_formula(name)).bind(variables) + assert set(formula.inputs) == set(LEGACY_INPUTS[name]) + # the main driver comes first (scripts plot against x.iloc[:, 0]) + assert formula.inputs[0] not in CATEGORICAL_LEVELS + mu, extras = formula( + pt.as_tensor_variable(df[list(formula.inputs)].to_numpy()), variables + ) + + np.testing.assert_allclose(mu.eval(), legacy_mu.eval(), rtol=1e-10) + np.testing.assert_allclose( + np.broadcast_to(extras["sigma"].eval(), N), + np.broadcast_to(legacy_extras["sigma"].eval(), N), + rtol=1e-10, + ) + assert float(extras["lower"].eval()) == float(legacy_extras["lower"].eval()) diff --git a/tests/test_models.py b/tests/test_models.py index ad74291..27fb8a9 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -42,3 +42,141 @@ def test_heating_cp_occ_rad_floors_heat_and_computes_sigma(): assert sigma[1] == pytest.approx(np.sqrt(s0**2 + (s1 * w1) ** 2)) # Colder day -> more heating -> higher w -> more noise added to the floor. assert sigma[0] > sigma[1] + + +def test_heating_dt_occ_rad_multiroom_sums_over_rooms(): + """Direct call (no PyMC involved) checking the room axis is summed, not + left dangling, and that occupation correctly selects a (n_obs, n_room) + slice of a (n_occ, n_room) parameter via fancy indexing on the first axis + only -- the failure mode this guards against is accidentally pairwise- + indexing occupation against room (like season_cp_heating_es_rad_base_by_ + period's base[is_heating, period]) instead of broadcasting one against + the other.""" + from bayesbuilding.models import heating_dt_occ_rad_multiroom + + n_obs, n_room = 5, 3 + x = np.zeros((n_obs, n_room + 3)) + x[:, :n_room] = np.arange(n_obs * n_room).reshape(n_obs, n_room) + x[:, n_room] = 1.0 # rad + x[:, n_room + 1] = 0.0 # elec_consumption + x[:, n_room + 2] = [0, 1, 0, 1, 0] # occupation + + g = np.array([[1.0, 2.0, 3.0], [10.0, 20.0, 30.0]]) + fs = np.zeros((2, n_room)) + variables_dict = { + "g": g, + "fs": fs, + "alpha": np.array([0.0, 0.0]), + "s0": np.array([1.0, 1.0]), + } + + mu, extras = heating_dt_occ_rad_multiroom(x, variables_dict) + + dt = x[:, :n_room] + occ = x[:, -1].astype(int) + expected = (g[occ] * dt).sum(axis=1) + + assert mu.shape == (n_obs,) + np.testing.assert_allclose(mu, expected) + + + +def test_heating_dt_occ_rad_lag_reheat_surcharge_sign(): + """Direct call (no PyMC) locking down the sign convention for h, worked + out from the underlying physics (see bayes.py docstring / plan): with a + thermal-mass recharge cost C/dt_step = 50 folded into g_fit=150, h=-50 on + top of a steady-state g=100, a stable day (dt == dt_lag == 10) must + reproduce the plain g*dt = 1000, while a reheat day coming out of a low + dt_lag=2 up to dt=10 must show the extra +400 surcharge (1400 total) -- + catches an accidental sign flip on h that unit shape/smoke tests alone + wouldn't.""" + from bayesbuilding.models import heating_dt_occ_rad_lag + + g_fit, h_fit = 150.0, -50.0 + x = np.array( + [ + [10.0, 0.0, 0.0, 1, 10.0], # stable day: dt == dt_lag + [10.0, 0.0, 0.0, 1, 2.0], # reheat day: dt_lag much lower than dt + ] + ) + variables_dict = { + "g": np.array([g_fit, g_fit]), + "fs": np.array([0.0, 0.0]), + "alpha": np.array([0.0, 0.0]), + "h": np.array([h_fit, h_fit]), + "s0": np.array([1.0, 1.0]), + } + + mu, extras = heating_dt_occ_rad_lag(x, variables_dict) + + np.testing.assert_allclose(mu, [1000.0, 1400.0]) + + + +def test_heating_dt_occ_rad_DTdt_lags_capacity_construction(): + """Direct call (no PyMC) locking down the c0 * rho**k geometric-decay + construction and the lag-column order: with c0=100 and rho=0.5, the + per-lag coefficients must be [100, 50, 25, 12.5] applied to DT_dt at + [lag0 (today), lag1, lag2, lag3] in that order -- catches an accidental + transpose of the lag columns or an off-by-one in the decay exponent that + a shape/smoke test alone wouldn't. Also locks down that every coefficient + stays positive (guaranteed by construction for c0 > 0, 0 < rho < 1), + unlike the earlier c0 + cumsum(d) version, whose unconstrained increments + let a higher lag's coefficient go negative -- an unphysical sign flip for + a fading capacitive effect.""" + from bayesbuilding.models import heating_dt_occ_rad_DTdt_lags + + c0_fit = 100.0 + rho_fit = 0.5 + dt_dt_lags = np.array([1.0, 2.0, 3.0, 4.0]) # lag0, lag1, lag2, lag3 + + x = np.array([[10.0, 0.0, 0.0, 1, *dt_dt_lags]]) + variables_dict = { + "g": np.array([0.0, 0.0]), + "fs": np.array([0.0, 0.0]), + "alpha": np.array([0.0, 0.0]), + "c0": c0_fit, + "rho": rho_fit, + "s0": np.array([1.0, 1.0]), + } + + mu, extras = heating_dt_occ_rad_DTdt_lags(x, variables_dict) + + c_expected = np.array([100.0, 50.0, 25.0, 12.5]) + assert (c_expected > 0).all() + expected = dt_dt_lags @ c_expected + + np.testing.assert_allclose(mu, [expected]) + + + +def test_heating_dt_occ_rad_DTdt_wall_Ci_adds_air_capacity_term(): + """Direct call (no PyMC) locking down the sign/construction of the new + Ci*DTint_intraday term: with every other coefficient zeroed out (g, fs, + alpha, C all 0, beta arbitrary since C=0 makes it irrelevant), mu must + equal exactly Ci * DTint_intraday -- catches an accidental sign flip or a + swapped column index for the new 7th input column.""" + from bayesbuilding.models import heating_dt_occ_rad_DTdt_wall_Ci + + Ci_fit = 8000.0 + DTint_intraday = np.array([1.5, -0.5]) + + x = np.array( + [ + [10.0, 0.0, 0.0, 1, 0.0, 0.0, DTint_intraday[0]], + [10.0, 0.0, 0.0, 0, 0.0, 0.0, DTint_intraday[1]], + ] + ) + variables_dict = { + "g": np.array([0.0, 0.0]), + "fs": np.array([0.0, 0.0]), + "alpha": np.array([0.0, 0.0]), + "C": 0.0, + "beta": 0.5, + "Ci": Ci_fit, + "s0": np.array([1.0, 1.0]), + } + + mu, extras = heating_dt_occ_rad_DTdt_wall_Ci(x, variables_dict) + + np.testing.assert_allclose(mu, Ci_fit * DTint_intraday) From 820b7affca75340205279172c3014a8b292e70ec Mon Sep 17 00:00:00 2001 From: BaptisteDE Date: Thu, 1 Oct 2026 16:05:45 +0200 Subject: [PATCH 2/2] =?UTF-8?q?=F0=9F=9A=A8=20lint?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- bayesbuilding/candidates.py | 5 +++- bayesbuilding/control_charts.py | 17 ++++++++--- bayesbuilding/formula.py | 31 ++++++++++++------- bayesbuilding/legacy_formulas.py | 4 ++- bayesbuilding/plotting.py | 16 ++++++---- bayesbuilding/training.py | 4 ++- tests/test_candidates_training.py | 20 ++++++++++--- tests/test_control_charts.py | 49 ++++++++++++++++++++++--------- tests/test_formula.py | 6 ++-- tests/test_models.py | 3 -- 10 files changed, 109 insertions(+), 46 deletions(-) diff --git a/bayesbuilding/candidates.py b/bayesbuilding/candidates.py index 67523df..065b199 100644 --- a/bayesbuilding/candidates.py +++ b/bayesbuilding/candidates.py @@ -11,7 +11,10 @@ "name": "dt_occ", "model": {"mu": "g[occ]*dt - fs[occ]*rad", "sigma": "s0[occ]", "lower": 0.0}, "inputs": {"dt": "", "rad": "", "occ": ""}, - "priors": {"g": {"dist": "HalfNormal", "kwargs": {"sigma": 10, "shape": 2}}, ...}, + "priors": { + "g": {"dist": "HalfNormal", "kwargs": {"sigma": 10, "shape": 2}}, + ... + }, "draws": 4000, "tune": 1000 } ], diff --git a/bayesbuilding/control_charts.py b/bayesbuilding/control_charts.py index 974e833..ba0beb6 100644 --- a/bayesbuilding/control_charts.py +++ b/bayesbuilding/control_charts.py @@ -236,7 +236,12 @@ def _out_of_limits(z, v): def cusum_control_stats( - y_true: pd.Series, mu_draws, sigma_draws, k: float = 0.5, h: float = 5.0, L: float = 1.96 + y_true: pd.Series, + mu_draws, + sigma_draws, + k: float = 0.5, + h: float = 5.0, + L: float = 1.96, ) -> dict: """Classical two-sided tabular CUSUM (Page) on the standardized residual ``z = residual / sigma_hat`` -- same ``mu_hat``/``sigma_hat`` point @@ -292,7 +297,6 @@ def cusum_control_stats( } - def _contiguous_true_runs(mask: np.ndarray) -> list[tuple]: """Positions ``(start_i, end_i)`` (inclusive) of each contiguous run of ``True`` in ``mask``. Simpler than @@ -429,7 +433,10 @@ def _plot_control_chart( align="left", xanchor="left", yanchor="top", - text=f"Points hors contrôle : {n_out}/{len(stat)} ({n_out / len(stat) * 100:.1f}%)", + text=( + f"Points hors contrôle : {n_out}/{len(stat)} " + f"({n_out / len(stat) * 100:.1f}%)" + ), ) fig.update_layout( @@ -552,7 +559,9 @@ def plot_cusum_chart( :func:`plot_xbar_chart`/:func:`plot_ewma_chart` qui suivent ``sigma_hat`` dans les unités brutes du résidu. """ - stats = cusum_control_stats(measure_ts, mu_prediction, sigma_prediction, k=k, h=h, L=L) + stats = cusum_control_stats( + measure_ts, mu_prediction, sigma_prediction, k=k, h=h, L=L + ) signed = stats["cusum_pos"] - stats["cusum_neg"] h_line = np.full_like(signed, stats["h"], dtype=float) return _plot_control_chart( diff --git a/bayesbuilding/formula.py b/bayesbuilding/formula.py index 845d53d..7700016 100644 --- a/bayesbuilding/formula.py +++ b/bayesbuilding/formula.py @@ -80,12 +80,16 @@ def _check(node: ast.AST, expr: str): """Walk ``node`` and raise FormulaError on anything outside the grammar.""" if isinstance(node, ast.BinOp): if type(node.op) not in _BINARY_OPERATORS: - raise FormulaError(f"Operator {type(node.op).__name__} not allowed in {expr!r}") + raise FormulaError( + f"Operator {type(node.op).__name__} not allowed in {expr!r}" + ) _check(node.left, expr) _check(node.right, expr) elif isinstance(node, ast.UnaryOp): if type(node.op) not in _UNARY_OPERATORS: - raise FormulaError(f"Operator {type(node.op).__name__} not allowed in {expr!r}") + raise FormulaError( + f"Operator {type(node.op).__name__} not allowed in {expr!r}" + ) _check(node.operand, expr) elif isinstance(node, ast.Constant): if isinstance(node.value, bool) or not isinstance(node.value, (int, float)): @@ -144,7 +148,9 @@ def _index_symbols(node: ast.AST) -> set[str]: for child in ast.walk(node): if isinstance(child, ast.Subscript): indices = ( - child.slice.elts if isinstance(child.slice, ast.Tuple) else [child.slice] + child.slice.elts + if isinstance(child.slice, ast.Tuple) + else [child.slice] ) names |= {i.id for i in indices if isinstance(i, ast.Name)} return names @@ -168,7 +174,9 @@ def _evaluate(node: ast.AST, env: dict, categorical: dict): _evaluate(node.right, env, categorical), ) if isinstance(node, ast.UnaryOp): - return _UNARY_OPERATORS[type(node.op)](_evaluate(node.operand, env, categorical)) + return _UNARY_OPERATORS[type(node.op)]( + _evaluate(node.operand, env, categorical) + ) if isinstance(node, ast.Constant): return node.value if isinstance(node, ast.Name): @@ -278,11 +286,11 @@ def bind(self, prior_names) -> "FormulaModel": return self def __call__(self, x, variables: dict): - inputs = ( - self.inputs if self.inputs is not None else self.input_names(variables) - ) + inputs = self.inputs if self.inputs is not None else self.input_names(variables) missing = [ - name for name in self.symbols if name not in inputs and name not in variables + name + for name in self.symbols + if name not in inputs and name not in variables ] if missing: raise FormulaError(f"Symbols {missing} are neither priors nor inputs") @@ -293,7 +301,8 @@ def __call__(self, x, variables: dict): for name in self.categorical_inputs } values = { - name: _evaluate(tree, env, categorical) for name, tree in self._trees.items() + name: _evaluate(tree, env, categorical) + for name, tree in self._trees.items() } mu = values.pop("mu") out = {"sigma": values.pop("sigma")} @@ -315,4 +324,6 @@ def from_dict(cls, spec: dict) -> "FormulaModel": return cls(**spec) def __repr__(self): - return f"FormulaModel(mu={self.mu!r}, sigma={self.sigma!r}, lower={self.lower!r})" + return ( + f"FormulaModel(mu={self.mu!r}, sigma={self.sigma!r}, lower={self.lower!r})" + ) diff --git a/bayesbuilding/legacy_formulas.py b/bayesbuilding/legacy_formulas.py index 497a389..85298b7 100644 --- a/bayesbuilding/legacy_formulas.py +++ b/bayesbuilding/legacy_formulas.py @@ -10,7 +10,9 @@ ``x[:, i]``. """ -_SETBACK_RAD = "switch(is_heating, base[0] + {g}*max({tau} - text, 0) - fs*rad, base[1])" +_SETBACK_RAD = ( + "switch(is_heating, base[0] + {g}*max({tau} - text, 0) - fs*rad, base[1])" +) _DT_CORE = "g[occ]*dt - fs[occ]*rad - alpha[occ]*elec_consumption" _WALL = "C*(beta*DTint_dt + (1 - beta)*DText_dt)" diff --git a/bayesbuilding/plotting.py b/bayesbuilding/plotting.py index 158c748..1e05621 100644 --- a/bayesbuilding/plotting.py +++ b/bayesbuilding/plotting.py @@ -1524,7 +1524,9 @@ def plot_candidate_metric_bars( by their x position/label, not by color). """ if metric not in comparison.columns: - raise KeyError(f"{metric!r} not in comparison columns: {list(comparison.columns)}") + raise KeyError( + f"{metric!r} not in comparison columns: {list(comparison.columns)}" + ) if std_col is None: candidate_cols = ( @@ -1537,7 +1539,9 @@ def plot_candidate_metric_bars( f"{candidate_cols}); pass std_col explicitly." ) elif std_col not in comparison.columns: - raise KeyError(f"{std_col!r} not in comparison columns: {list(comparison.columns)}") + raise KeyError( + f"{std_col!r} not in comparison columns: {list(comparison.columns)}" + ) ranked = comparison[[metric, std_col]].sort_values(metric, ascending=ascending) @@ -1557,7 +1561,11 @@ def plot_candidate_metric_bars( if len(threshold_values) == 2: low, high = sorted(threshold_values) fig.add_hrect( - y0=low, y1=high, fillcolor="rgba(44,160,44,0.08)", line_width=0, layer="below" + y0=low, + y1=high, + fillcolor="rgba(44,160,44,0.08)", + line_width=0, + layer="below", ) for value in threshold_values: fig.add_hline( @@ -1576,5 +1584,3 @@ def plot_candidate_metric_bars( showlegend=False, ) return fig - - diff --git a/bayesbuilding/training.py b/bayesbuilding/training.py index dcf02e8..80a31ee 100644 --- a/bayesbuilding/training.py +++ b/bayesbuilding/training.py @@ -23,7 +23,9 @@ def cv_rmse(y_true, y_pred) -> float: """Coefficient of Variation of the RMSE [%] (same definition as corrai's, i.e. ASHRAE Guideline 14's, with n - 1 degrees of freedom).""" n = len(y_true) - return float(np.sqrt(np.sum((y_true - y_pred) ** 2) / (n - 1)) / np.mean(y_true) * 100) + return float( + np.sqrt(np.sum((y_true - y_pred) ** 2) / (n - 1)) / np.mean(y_true) * 100 + ) def mae(y_true, y_pred) -> float: diff --git a/tests/test_candidates_training.py b/tests/test_candidates_training.py index 124ea66..ec09e10 100644 --- a/tests/test_candidates_training.py +++ b/tests/test_candidates_training.py @@ -40,7 +40,11 @@ "candidates": [ { "name": "dt_change_point", - "model": {"mu": "g*max(dt - tau, 0) + base", "sigma": "sigma", "lower": 0.0}, + "model": { + "mu": "g*max(dt - tau, 0) + base", + "sigma": "sigma", + "lower": 0.0, + }, "inputs": {"dt": "dt__C__Zone__FULL"}, "priors": PRIORS_DT, "prior_draws": 500, @@ -77,7 +81,9 @@ def test_from_json_round_trip(self, tmp_path): assert config.get_candidate("text_change_point").model.inputs == ("text",) config.to_json(tmp_path / "copy.json") - assert BayesConfig.from_json(tmp_path / "copy.json").to_dict() == config.to_dict() + assert ( + BayesConfig.from_json(tmp_path / "copy.json").to_dict() == config.to_dict() + ) def test_unused_candidates_is_optional(self): payload = {k: v for k, v in CONFIG_PAYLOAD.items() if k != "unused_candidates"} @@ -199,7 +205,10 @@ def test_flatten_regime_vars_explodes_shape_2_regime_variable(): group = xr.Dataset( { "sigma": (("chain", "draw"), [[1.0, 2.0, 3.0]]), - "g": (("chain", "draw", "g_dim_0"), [[[10.0, 20.0], [11.0, 21.0], [12.0, 22.0]]]), + "g": ( + ("chain", "draw", "g_dim_0"), + [[[10.0, 20.0], [11.0, 21.0], [12.0, 22.0]]], + ), } ) flattened = flatten_regime_vars(group, ["sigma", "g", "does_not_exist"]) @@ -221,7 +230,10 @@ def test_train_candidates_fits_scores_and_records_failures(tmp_path): idx = pd.date_range("2024-01-01", periods=70, freq="D") dt = rng.uniform(0, 15, 70) df = pd.DataFrame( - {"dt__C": dt, "energy": 4.0 * np.maximum(dt - 3, 0) + 10 + rng.normal(0, 1, 70)}, + { + "dt__C": dt, + "energy": 4.0 * np.maximum(dt - 3, 0) + 10 + rng.normal(0, 1, 70), + }, index=idx, ) priors = { diff --git a/tests/test_control_charts.py b/tests/test_control_charts.py index 6051b9f..8d02307 100644 --- a/tests/test_control_charts.py +++ b/tests/test_control_charts.py @@ -116,8 +116,14 @@ def test_ewma_control_stats_confirms_alarm_then_resets_after_transient_event(): lam, L = 0.5, 1.0 stats = ewma_control_stats( - y_true, mu_draws, sigma_draws, lam=lam, L=L, - n_confirm=2, n_reset=2, reset_epsilon=0.5, + y_true, + mu_draws, + sigma_draws, + lam=lam, + L=L, + n_confirm=2, + n_reset=2, + reset_epsilon=0.5, ) expected_alarm = [False, False, True, True, True, False, False, False] @@ -141,8 +147,14 @@ def test_ewma_control_stats_persistent_drift_never_resets(): y_true = np.array([0.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0, 5.0]) stats = ewma_control_stats( - y_true, mu_draws, sigma_draws, lam=0.5, L=1.0, - n_confirm=2, n_reset=2, reset_epsilon=0.5, + y_true, + mu_draws, + sigma_draws, + lam=0.5, + L=1.0, + n_confirm=2, + n_reset=2, + reset_epsilon=0.5, ) assert stats["alarm"][2:].all() @@ -159,8 +171,14 @@ def test_ewma_control_stats_short_excess_never_confirms_alarm(): y_true = np.array([0.0, 5.0, 0.0, 0.0, 0.0]) stats = ewma_control_stats( - y_true, mu_draws, sigma_draws, lam=0.5, L=1.0, - n_confirm=10, n_reset=2, reset_epsilon=0.5, + y_true, + mu_draws, + sigma_draws, + lam=0.5, + L=1.0, + n_confirm=10, + n_reset=2, + reset_epsilon=0.5, ) assert not stats["alarm"].any() @@ -195,7 +213,8 @@ def test_cusum_control_stats_accumulates_a_sustained_downward_drift(): def test_cusum_control_stats_flags_out_of_control_once_h_is_exceeded(): - # z=2.0, k=0.5 -> cusum_pos grows by 1.5/day: 1.5, 3.0, 4.5, 6.0 -> crosses h=5 on day 4. + # z=2.0, k=0.5 -> cusum_pos grows by 1.5/day: 1.5, 3.0, 4.5, 6.0 + # -> crosses h=5 on day 4. mu_draws = np.zeros((1, 4)) sigma_draws = np.ones((1, 4)) y_true = np.array([2.0, 2.0, 2.0, 2.0]) @@ -247,20 +266,22 @@ def test_plot_ewma_chart_without_reset_args_has_no_alarm_shading_or_markers(): def test_plot_ewma_chart_with_reset_args_shades_alarm_and_marks_resets(): index = pd.date_range("2024-01-01", periods=8, freq="D") - measure = pd.Series( - [0.0, 5.0, 5.0, 5.0, 0.0, 0.0, 0.0, 0.0], index=index - ) + measure = pd.Series([0.0, 5.0, 5.0, 5.0, 0.0, 0.0, 0.0, 0.0], index=index) mu_draws = np.zeros((1, 8)) sigma_draws = np.ones((1, 8)) fig = plot_ewma_chart( - measure, mu_draws, sigma_draws, lam=0.5, L=1.0, - n_confirm=2, n_reset=2, reset_epsilon=0.5, + measure, + mu_draws, + sigma_draws, + lam=0.5, + L=1.0, + n_confirm=2, + n_reset=2, + reset_epsilon=0.5, ) # the y=0 center line, plus one shaded vrect for the confirmed alarm run assert len(fig.layout.shapes) == 2 reset_trace = next(trace for trace in fig.data if trace.name == "reset EWMA") assert list(reset_trace.x) == [index[5]] - - diff --git a/tests/test_formula.py b/tests/test_formula.py index b058589..a85ad65 100644 --- a/tests/test_formula.py +++ b/tests/test_formula.py @@ -152,9 +152,9 @@ def test_wrapper_save_load_round_trip_with_formula(toy_data): "base": (pm.Normal, {"name": "base", "mu": 15, "sigma": 10, "shape": 2}), "sigma": (pm.HalfNormal, {"name": "sigma", "sigma": 3}), } - model = FormulaModel( - mu="g*max(tau - text, 0) + base[occ]", sigma="sigma" - ).bind(priors) + model = FormulaModel(mu="g*max(tau - text, 0) + base[occ]", sigma="sigma").bind( + priors + ) wrapper = PymcWrapper( model_function=model, priors_dict=priors, diff --git a/tests/test_models.py b/tests/test_models.py index 27fb8a9..a4eafdc 100644 --- a/tests/test_models.py +++ b/tests/test_models.py @@ -80,7 +80,6 @@ def test_heating_dt_occ_rad_multiroom_sums_over_rooms(): np.testing.assert_allclose(mu, expected) - def test_heating_dt_occ_rad_lag_reheat_surcharge_sign(): """Direct call (no PyMC) locking down the sign convention for h, worked out from the underlying physics (see bayes.py docstring / plan): with a @@ -112,7 +111,6 @@ def test_heating_dt_occ_rad_lag_reheat_surcharge_sign(): np.testing.assert_allclose(mu, [1000.0, 1400.0]) - def test_heating_dt_occ_rad_DTdt_lags_capacity_construction(): """Direct call (no PyMC) locking down the c0 * rho**k geometric-decay construction and the lag-column order: with c0=100 and rho=0.5, the @@ -149,7 +147,6 @@ def test_heating_dt_occ_rad_DTdt_lags_capacity_construction(): np.testing.assert_allclose(mu, [expected]) - def test_heating_dt_occ_rad_DTdt_wall_Ci_adds_air_capacity_term(): """Direct call (no PyMC) locking down the sign/construction of the new Ci*DTint_intraday term: with every other coefficient zeroed out (g, fs,