diff --git a/benchmarking/baselines/block_bootstrap.py b/benchmarking/baselines/block_bootstrap.py index 3a2add4e..b3f15c5b 100644 --- a/benchmarking/baselines/block_bootstrap.py +++ b/benchmarking/baselines/block_bootstrap.py @@ -15,7 +15,7 @@ from __future__ import annotations import math -from dataclasses import dataclass +from dataclasses import dataclass, field from typing import TYPE_CHECKING import numpy as np @@ -32,6 +32,9 @@ # Per (cell, segment): the numerator and denominator of each segment's rho. _N_QUANTITIES = 4 _MIN_RESAMPLES_FOR_SPREAD = 2 +# A resample whose reference energy on either side falls below this share of the whole campaign's is a +# degenerate ratio (a calm spell or outage drawn over and over), not a draw from the estimator's spread. +_MIN_REFERENCE_FRACTION = 0.1 # With one block covering the whole campaign, every resample is that campaign: nothing varies. _MIN_BLOCKS_FOR_SPREAD = 2 # Normal -/+1 sigma percentiles, so (p84 - p16) / 2 is a sigma for a normal. @@ -74,10 +77,116 @@ class BootstrapResult: :param n_blocks: blocks drawn per resample (``ceil(campaign / block)``) :param cells: uncertainty per cell name, keyed as the caller keyed ``cell_membership`` + :param resamples: per cell, the resampled uplifts themselves (NaN where a resample was + degenerate), in draw order; empty when the bootstrap did not run. Two bootstraps run with the + same campaign span, timebase, block length, resample count and seed draw the same blocks, so + their resamples pair up draw by draw (see :func:`combine_estimates`). """ n_blocks: int cells: dict[str, CellUncertainty] + resamples: dict[str, npt.NDArray[np.float64]] = field(default_factory=dict) + + +@dataclass(frozen=True) +class CombinedEstimate: + """An inverse-variance blend of two estimates of the same uplift. + + :param estimate: ``weight_a * a + (1 - weight_a) * b`` + :param sigma: the blend's 1-sigma, including the two components' correlation + :param weight_a: the weight on the first component, in ``[0, 1]``; NaN when neither was usable + :param correlation: the correlation of the two components' paired resamples that ``sigma`` used + """ + + estimate: float + sigma: float + weight_a: float + correlation: float + + +# Correlation assumed between two components when their resamples cannot be paired: midway between +# independent and identical, so the blend neither claims the full inverse-variance gain nor none. +_DEFAULT_COMPONENT_CORRELATION = 0.5 + + +def combine_estimates( + a: tuple[float, float, npt.NDArray[np.float64]], + b: tuple[float, float, npt.NDArray[np.float64]], + *, + weighting: tuple[float, float, npt.NDArray[np.float64], npt.NDArray[np.float64]] | None = None, +) -> CombinedEstimate: + """Blend two ``(estimate, sigma, resamples)`` of the same quantity by minimum-variance weights. + + A component is usable when both its estimate and its sigma are finite. With two usable + components the weights are the minimum-variance ones for two correlated estimates, + ``wa ∝ max(sb² - r sa sb, 0)`` and ``wb ∝ max(sa² - r sa sb, 0)`` (a zero sigma takes all the + weight; both weights zero splits evenly), with ``r`` the Pearson correlation of the draw-paired + finite resamples, or :data:`_DEFAULT_COMPONENT_CORRELATION` when fewer than two pairs exist. + The blend's variance is ``wa² sa² + wb² sb² + 2 wa wb r sa sb``, which these weights keep at or + below the better component's: a noisier leg correlated at ``r >= sb / sa`` gets no weight rather + than widening the blend. One usable component is returned as is; none gives NaN throughout. + + ``weighting``, when given as ``(sigma_a, sigma_b, resamples_a, resamples_b)`` from another + bootstrap of the same two components (e.g. one with the on/off labels permuted), supplies the + weights instead, so they need not read the realised estimates' own spread; the blend's variance + and the reported correlation still come from ``a`` and ``b``. A ``weighting`` with a non-finite + sigma is ignored. + """ + est_a, sig_a, res_a = a + est_b, sig_b, res_b = b + usable_a = math.isfinite(est_a) and math.isfinite(sig_a) + usable_b = math.isfinite(est_b) and math.isfinite(sig_b) + nan = float("nan") + if usable_a and not usable_b: + return CombinedEstimate(estimate=est_a, sigma=sig_a, weight_a=1.0, correlation=nan) + if usable_b and not usable_a: + return CombinedEstimate(estimate=est_b, sigma=sig_b, weight_a=0.0, correlation=nan) + if not (usable_a and usable_b): + return CombinedEstimate(estimate=nan, sigma=nan, weight_a=nan, correlation=nan) + + correlation = _paired_correlation(res_a, res_b) + if weighting is not None and math.isfinite(weighting[0]) and math.isfinite(weighting[1]): + sig_a_w, sig_b_w, res_a_w, res_b_w = weighting + weight_a = _min_variance_weight(sig_a_w, sig_b_w, _paired_correlation(res_a_w, res_b_w)) + else: + weight_a = _min_variance_weight(sig_a, sig_b, correlation) + weight_b = 1.0 - weight_a + variance = weight_a**2 * sig_a**2 + weight_b**2 * sig_b**2 + 2.0 * weight_a * weight_b * correlation * sig_a * sig_b + return CombinedEstimate( + estimate=weight_a * est_a + weight_b * est_b, + sigma=math.sqrt(max(variance, 0.0)), + weight_a=weight_a, + correlation=correlation, + ) + + +def _paired_correlation(res_a: npt.NDArray[np.float64], res_b: npt.NDArray[np.float64]) -> float: + """Pearson correlation of two components' draw-paired finite resamples, or the default.""" + if len(res_a) and len(res_b) and len(res_a) != len(res_b): + msg = ( + f"the two components' resample arrays have different lengths ({len(res_a)} vs {len(res_b)}); " + f"they can only be paired draw by draw when both bootstraps ran with the same campaign span, " + f"timebase, block length, resample count and seed." + ) + raise ValueError(msg) + n_pairs = min(len(res_a), len(res_b)) + paired = np.isfinite(res_a[:n_pairs]) & np.isfinite(res_b[:n_pairs]) + if paired.sum() < _MIN_RESAMPLES_FOR_SPREAD: + return _DEFAULT_COMPONENT_CORRELATION + pa, pb = res_a[:n_pairs][paired], res_b[:n_pairs][paired] + if np.ptp(pa) == 0.0 or np.ptp(pb) == 0.0: + return _DEFAULT_COMPONENT_CORRELATION + r = float(np.corrcoef(pa, pb)[0, 1]) + return r if math.isfinite(r) else _DEFAULT_COMPONENT_CORRELATION + + +def _min_variance_weight(sig_a: float, sig_b: float, correlation: float) -> float: + """Weight on ``a`` that minimises the variance of a blend of two estimates correlated at ``correlation``.""" + if sig_a == 0.0 or sig_b == 0.0: + return 0.5 if sig_a == sig_b else float(sig_a == 0.0) + cross = correlation * sig_a * sig_b + wa, wb = max(sig_b**2 - cross, 0.0), max(sig_a**2 - cross, 0.0) + return 0.5 if wa + wb == 0.0 else wa / (wa + wb) def _nan_cells(names: list[str]) -> dict[str, CellUncertainty]: @@ -222,11 +331,12 @@ def bootstrap_ratio_uplift( # vary. Any sigma it returned would be float residue (~1e-15), not a real certainty. return BootstrapResult(n_blocks=n_blocks, cells=_fallback_only_cells(names, fallback=fallback)) - uplift = _uplift_from_totals(totals) + uplift = _uplift_from_totals(totals, reference_floor=_MIN_REFERENCE_FRACTION * prefix[n_records]) boot_weight = {name: _bootstrap_weight(min(on, off)) for name, (on, off) in counts.items()} return BootstrapResult( n_blocks=n_blocks, cells=_summarise(uplift, names=names, fallback=fallback, boot_weight=boot_weight), + resamples={name: uplift[:, i] for i, name in enumerate(names)}, ) @@ -280,17 +390,22 @@ def _prefix_sums( return np.concatenate([np.zeros((1, len(names), _N_QUANTITIES)), np.cumsum(doubled, axis=0)], axis=0) -def _uplift_from_totals(totals: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]: +def _uplift_from_totals( + totals: npt.NDArray[np.float64], *, reference_floor: npt.NDArray[np.float64] | None = None +) -> npt.NDArray[np.float64]: """Re-form ``rho_up / rho_base - 1`` per (resample, cell) from resampled sums. The degeneracy guards mirror the point estimate's, so a resample fails only where the point - estimate would have failed on the same rows. + estimate would have failed on the same rows; ``reference_floor`` (per cell, same last axis as + ``totals``) additionally fails a resample whose on or off reference sum falls below it. """ test_on, ref_on, test_off, ref_off = (totals[..., k] for k in range(_N_QUANTITIES)) nan = np.full(test_on.shape, np.nan) rho_up = np.divide(test_on, ref_on, out=nan.copy(), where=ref_on != 0) rho_base = np.divide(test_off, ref_off, out=nan.copy(), where=ref_off != 0) valid = np.isfinite(rho_base) & (rho_base != 0) & np.isfinite(rho_up) + if reference_floor is not None: + valid &= (ref_on >= reference_floor[..., 1]) & (ref_off >= reference_floor[..., 3]) return np.divide(rho_up, rho_base, out=nan.copy(), where=valid) - 1.0 @@ -328,7 +443,14 @@ def _summarise( def _blend(*, bootstrap: float, fallback: float, weight: float) -> float: - """Linear ramp ``weight*bootstrap + (1-weight)*fallback``, using whichever side is finite.""" + """Linear ramp ``weight*bootstrap + (1-weight)*fallback``, using whichever side is finite. + + A cell at the sparse end (``weight == 0``) reports the fallback alone, even when that is NaN: + with so few records every resample that sees the cell is the same handful, so the bootstrap's + std is float residue (~1e-17), which would read as near-certainty where there is none. + """ + if weight == 0.0: + return fallback if math.isfinite(bootstrap) and math.isfinite(fallback): return weight * bootstrap + (1.0 - weight) * fallback if math.isfinite(bootstrap): diff --git a/benchmarking/baselines/example_toggle_study.py b/benchmarking/baselines/example_toggle_study.py index 7990b118..955797ca 100644 --- a/benchmarking/baselines/example_toggle_study.py +++ b/benchmarking/baselines/example_toggle_study.py @@ -51,6 +51,8 @@ # 20 minutes on, 20 minutes off -> a 40-minute on/off cycle. DEFAULT_TOGGLE_PERIOD = pd.Timedelta(minutes=40) +# Cycles are counted from the Unix epoch, as a clock-aligned controller counts them. +DEFAULT_TOGGLE_DATUM = pd.Timestamp("1970-01-01", tz="UTC") # One upgrade for the toggle study: a flat +3% Cp change in region 2. TOGGLE_PROFILES: dict[str, list] = {"cp_plus_3pct": [ConstantCpChange(delta=0.03)]} @@ -199,6 +201,7 @@ def main( min_pre_months=MIN_PRE_MONTHS, campaign_months=campaign_months, toggle_period=toggle_period, + toggle_datum=DEFAULT_TOGGLE_DATUM, n_replicates=n_replicates, seed=0, ) diff --git a/benchmarking/baselines/old/inspect_naive.py b/benchmarking/baselines/old/inspect_naive.py index 9bad69c2..822321b8 100644 --- a/benchmarking/baselines/old/inspect_naive.py +++ b/benchmarking/baselines/old/inspect_naive.py @@ -171,6 +171,7 @@ def inspect_naive_run( campaign_months=campaign_months, n_replicates=n_replicates, toggle_period=toggle_period, + toggle_datum=pd.Timestamp("1970-01-01", tz="UTC"), seed=0, ) diff --git a/benchmarking/baselines/study_overnight_toggle.py b/benchmarking/baselines/study_overnight_toggle.py index cd5e9217..83918a43 100644 --- a/benchmarking/baselines/study_overnight_toggle.py +++ b/benchmarking/baselines/study_overnight_toggle.py @@ -26,6 +26,7 @@ MIN_PRE_MONTHS, ) from benchmarking.baselines.example_toggle_study import ( + DEFAULT_TOGGLE_DATUM, DEFAULT_TOGGLE_PERIOD, default_output_root, run_toggle_study, @@ -52,6 +53,7 @@ def main() -> None: min_pre_months=MIN_PRE_MONTHS, campaign_months=CAMPAIGN_MONTHS, toggle_period=DEFAULT_TOGGLE_PERIOD, + toggle_datum=DEFAULT_TOGGLE_DATUM, n_replicates=N_REPLICATES, seed=0, ) diff --git a/benchmarking/baselines/study_power_model_compare.py b/benchmarking/baselines/study_power_model_compare.py index c90fef78..270b670a 100644 --- a/benchmarking/baselines/study_power_model_compare.py +++ b/benchmarking/baselines/study_power_model_compare.py @@ -93,7 +93,7 @@ MIN_PRE_MONTHS, save_per_method_curve, ) -from benchmarking.baselines.example_toggle_study import DEFAULT_TOGGLE_PERIOD +from benchmarking.baselines.example_toggle_study import DEFAULT_TOGGLE_DATUM, DEFAULT_TOGGLE_PERIOD from benchmarking.baselines.hot_context import build_hot_v0_context from benchmarking.baselines.naive_ratio import NaiveRatioMethod from benchmarking.baselines.overnight_profiles import overnight_profiles @@ -171,6 +171,7 @@ def _toggle_study() -> StudyConfig: min_pre_months=MIN_PRE_MONTHS, campaign_months=TOGGLE_CAMPAIGN_MONTHS, toggle_period=DEFAULT_TOGGLE_PERIOD, + toggle_datum=DEFAULT_TOGGLE_DATUM, n_replicates=N_REPLICATES, seed=SEED, ) diff --git a/benchmarking/baselines/study_toggle_methods_compare.py b/benchmarking/baselines/study_toggle_methods_compare.py index c4e3320f..43e31523 100644 --- a/benchmarking/baselines/study_toggle_methods_compare.py +++ b/benchmarking/baselines/study_toggle_methods_compare.py @@ -61,7 +61,7 @@ DEFAULT_WTG_NUMBERS, MIN_PRE_MONTHS, ) -from benchmarking.baselines.example_toggle_study import DEFAULT_TOGGLE_PERIOD +from benchmarking.baselines.example_toggle_study import DEFAULT_TOGGLE_DATUM, DEFAULT_TOGGLE_PERIOD from benchmarking.baselines.hot_context import build_hot_v0_context from benchmarking.baselines.power_model import PowerModelMethod from benchmarking.baselines.toggle_specialist import ToggleSpecialistMethod @@ -171,6 +171,7 @@ def toggle_study() -> StudyConfig: min_pre_months=MIN_PRE_MONTHS, campaign_weeks=CAMPAIGN_WEEKS, toggle_period=DEFAULT_TOGGLE_PERIOD, + toggle_datum=DEFAULT_TOGGLE_DATUM, n_replicates=N_REPLICATES, seed=SEED, ) @@ -202,6 +203,8 @@ def _build_methods(out_dir: Path, *, era5_hourly_df: pd.DataFrame) -> list: return [ ToggleSpecialistMethod( columns=HOT_COLUMNS, + toggle_period=DEFAULT_TOGGLE_PERIOD, + toggle_datum=DEFAULT_TOGGLE_DATUM, conditions=("power",), rated_power_kw=HOT_RATED_POWER_KW, out_dir=out_dir / "toggle_specialist_runs", diff --git a/benchmarking/baselines/study_toggle_specialist_uncertainty.py b/benchmarking/baselines/study_toggle_specialist_uncertainty.py index 44bfc536..e2070fc2 100644 --- a/benchmarking/baselines/study_toggle_specialist_uncertainty.py +++ b/benchmarking/baselines/study_toggle_specialist_uncertainty.py @@ -49,7 +49,7 @@ DEFAULT_TURBINE_SUBSET, DEFAULT_WTG_NUMBERS, ) -from benchmarking.baselines.example_toggle_study import DEFAULT_TOGGLE_PERIOD +from benchmarking.baselines.example_toggle_study import DEFAULT_TOGGLE_DATUM, DEFAULT_TOGGLE_PERIOD from benchmarking.baselines.study_toggle_methods_compare import TOGGLE_PROFILES from benchmarking.baselines.toggle_specialist import DEFAULT_BLOCK_HOURS, ToggleSpecialistMethod from benchmarking.harness import ( @@ -99,6 +99,7 @@ def uncertainty_study(n_replicates: int) -> StudyConfig: min_pre_months=MIN_PRE_MONTHS_TOGGLE, campaign_weeks=CAMPAIGN_WEEKS, toggle_period=DEFAULT_TOGGLE_PERIOD, + toggle_datum=DEFAULT_TOGGLE_DATUM, n_replicates=n_replicates, seed=SEED, ) @@ -129,6 +130,8 @@ def build_methods(block_hours_grid: list[float], out_dir: Path | None = None) -> return [ ToggleSpecialistMethod( columns=HOT_COLUMNS, + toggle_period=DEFAULT_TOGGLE_PERIOD, + toggle_datum=DEFAULT_TOGGLE_DATUM, name=f"toggle_specialist_bl{block_hours:g}", conditions=("power",), rated_power_kw=HOT_RATED_POWER_KW, diff --git a/benchmarking/baselines/toggle_specialist.py b/benchmarking/baselines/toggle_specialist.py index 09544eb6..60a976fd 100644 --- a/benchmarking/baselines/toggle_specialist.py +++ b/benchmarking/baselines/toggle_specialist.py @@ -12,6 +12,36 @@ never conditions on the test turbine's post-treatment wind speed (design-note §3). It speaks the data source's own column names and has no wind_up dependency. +An optional **pairing filter** (``pairing_max_gap``) then drops any used row with no used row of the +other segment nearby, so a filter that removes one segment's rows in some conditions also removes the +other segment's rows from those conditions. Rows kept per segment after each selection stage are +reported as ``MethodOutput.selection_accounting`` and written to a selection CSV. A caller may also +flag rows through ``columns.exclude_block``: such a row removes its whole toggle cycle, so both +states lose the same conditions — a symmetric alternative to ``exclude_row`` for a selection that +would otherwise thin one state's rows in particular conditions. + +**Two legs, blended.** The test power is compared against two references, each estimated on +its own row selection, and the two estimates are blended by minimum-variance weights, headline and +per bin. The *reference leg* (``sum``) compares against the summed power of a **reference subset**: +the ``sum`` estimate runs on every non-empty subset of the available references and keeps the +subset whose noise is smallest (ties to the larger subset), so a reference that tracks badly, or is +screened out for long stretches, loses to the subsets without it; every subset's result is written +to a CSV. The *block leg* (``block_mean``) needs **no references at all**: the campaign is tiled +into fixed wall-clock toggle cycles of ``toggle_period`` and every used row of a cycle is +compared against the cycle's mean test power, so the on/off ratio cancels the cycle's wind level +exactly; a cycle whose used rows lack either state is dropped (the ``block`` selection +stage). The two legs' errors come from different places (reference mismatch versus within-cycle +drift), which is what the blend exploits; their bootstraps share their block draws, so the blend's +sigma includes the legs' correlation (:func:`benchmarking.baselines.block_bootstrap.combine_estimates`). +With no reference turbine at all the estimate is the block leg alone. + +**Decisions on permuted labels.** The two data-driven decisions, the subset choice and the blend +weights, are judged on a bootstrap of the same rows with the on/off labels shuffled within each +toggle cycle: the same noise, blind to how it happened to split between the states. Ranking on +the realised sigma would favour the alternative whose realised uplift came out low (the on/off +contrast noise is skewed), biasing the estimate down and under-reporting sigma. The reported uplift +and sigma are always the actual ones; only the decisions read the permuted bootstrap. + Every uplift — the headline and each power bin — comes with a non-optional 1-sigma uncertainty from a circular block bootstrap (:mod:`benchmarking.baselines.block_bootstrap`). It is computed after the uplift, from the uplift's own frozen row selection and bin assignment, and only when the uplift is @@ -20,25 +50,41 @@ Each run writes a per-run folder ``toggle_specialist___/`` (v0-style naming) under ``out_dir`` (a temp dir by default), holding a per-segment data-stats CSV, -a headline results CSV, and -- when ``save_plots`` -- three diagnostic plots (a test-vs-reference -scatter, a per-segment daily-ratio timeseries, and a per-segment used-data-coverage timeseries). -The rich stats let a human confirm the right data was received and interpreted: the headline uplift -is re-derivable from the stats CSV as ``rho = used_test_mwh / used_ref_total_mwh`` per segment. +a headline results CSV, the per-leg selection and reference-subset CSVs, and -- when +``save_plots`` -- four diagnostic plots (a test-vs-reference scatter, a per-segment daily-ratio +timeseries, a per-segment used-data-coverage timeseries, and a pairing-gap histogram). +The rich stats let a human confirm the right data was received and interpreted: each leg's uplift +is re-derivable from its rows of the stats CSV as ``rho = used_test_mwh / used_ref_total_mwh`` per +segment, and the headline is the blend of the legs at the weight the results CSV reports. + +**Bin power.** The power bins and ``power_band`` both read one per-row power, chosen by +``bin_power``. ``"baseline"`` (the default) is what the turbine makes untreated: ``rho_base * +ref_total`` on the reference leg, the mean of the cycle's baseline rows on the block leg; no +upgraded row enters it. ``"both_states"`` is the mean of the two states instead (``rho`` averaged +over the states; the two states' cycle means averaged), for a campaign in which neither state is +a true baseline. Either way a row's own state cannot move it between bins. """ from __future__ import annotations +import logging import tempfile from dataclasses import dataclass, replace +from itertools import combinations from pathlib import Path -from typing import TYPE_CHECKING +from typing import TYPE_CHECKING, Literal, NamedTuple import matplotlib.pyplot as plt import numpy as np import pandas as pd from matplotlib.ticker import PercentFormatter -from benchmarking.baselines.block_bootstrap import BootstrapResult, bootstrap_ratio_uplift +from benchmarking.baselines.block_bootstrap import ( + BootstrapResult, + CombinedEstimate, + bootstrap_ratio_uplift, + combine_estimates, +) from benchmarking.baselines.filtering import NormalOperationFilter from benchmarking.diagnostics import DiagnosticContext, stages, write_common_diagnostics, write_run_config from benchmarking.harness.conditions import condition_bins, energy_ratio_by_bin, validate_conditions @@ -47,6 +93,8 @@ from benchmarking.synthetic import ToggleSchedule if TYPE_CHECKING: + from collections.abc import Mapping + import numpy.typing as npt from benchmarking.synthetic import ColumnSchema @@ -68,6 +116,123 @@ # treatment cannot move a row between bins. Binning by the test turbine's ws/TI would condition on # post-treatment signals, which this method exists not to do (see the module docstring). _SUPPORTED_CONDITIONS: tuple[str, ...] = ("power",) +# The two legs the estimate blends (see the module docstring), in blend order: the first is the one +# the row-selection accounting and the reference-side diagnostics report. +_REFERENCE_LEG = "sum" +_BLOCK_LEG = "block_mean" +_COMBINED = "combined" +# The subset choice runs 2**n - 1 estimates; refuse past this many references. +_MAX_REFS_FOR_SUBSETS = 8 +_BIN_POWERS = ("baseline", "both_states") + +logger = logging.getLogger(__name__) + + +class _Selection(NamedTuple): + """The cumulative used-row mask after each selection stage.""" + + filtered: pd.Series + not_excluded: pd.Series + paired: pd.Series + blocked: pd.Series + + +@dataclass(frozen=True) +class _Estimate: + """One leg's (or the blend's) complete estimate: what the outputs, the labels and the blend read.""" + + mode: str + refs: list[str] + selection: _Selection + used: npt.NDArray[np.bool_] + ref_total: npt.NDArray[np.float64] + label: npt.NDArray[np.float64] + rho_base: float + rho_up: float + uplift: float + sigma_overall: float + per_bin: pd.DataFrame | None + boot: BootstrapResult | None + accounting: pd.DataFrame + diagnostics: pd.DataFrame + # The label-permuted bootstrap the decisions read; None for the blend (which decides nothing further). + decision_boot: BootstrapResult | None = None + + @property + def decision_sigma(self) -> float: + """The headline sigma the decisions are judged on: that of the label-permuted bootstrap.""" + return _cell_sigma(self.decision_boot, _OVERALL) + + +def permute_labels_within_blocks( + *, + upgraded: npt.NDArray[np.bool_], + baseline: npt.NDArray[np.bool_], + used: npt.NDArray[np.bool_], + ids: npt.NDArray[np.int64] | None, + rng: np.random.Generator, +) -> tuple[npt.NDArray[np.bool_], npt.NDArray[np.bool_]]: + """Shuffle the on/off labels among the used campaign rows of each block; unused rows keep theirs. + + Every block keeps its on and off counts, so the permuted campaign has the same structure as the + real one and differs only in which rows are called on. With ``ids`` None the labels are shuffled + over the whole campaign. + """ + up_p = upgraded.copy() + base_p = baseline.copy() + campaign_used = np.flatnonzero(used & (upgraded | baseline)) + if ids is None: + groups = [campaign_used] + else: + groups = [campaign_used[ids[campaign_used] == b] for b in np.unique(ids[campaign_used])] + for rows in groups: + if len(rows) < 2: # noqa: PLR2004 - nothing to shuffle + continue + order = rng.permutation(rows) + up_p[rows] = upgraded[order] + base_p[rows] = baseline[order] + return up_p, base_p + + +class PairedRows(NamedTuple): + """The rows of each segment that have a row of the other segment within the pairing gap.""" + + baseline: npt.NDArray[np.bool_] + upgraded: npt.NDArray[np.bool_] + + +def pair_within( + index: pd.DatetimeIndex, + *, + baseline: npt.NDArray[np.bool_], + upgraded: npt.NDArray[np.bool_], + max_gap: pd.Timedelta, +) -> PairedRows: + """Keep each flagged row only if the other segment has a flagged row within ``max_gap``, inclusive.""" + gap = gap_to_other_segment(index, baseline=baseline, upgraded=upgraded) + within = gap <= max_gap.total_seconds() + return PairedRows(baseline=baseline & within, upgraded=upgraded & within) + + +def gap_to_other_segment( + index: pd.DatetimeIndex, *, baseline: npt.NDArray[np.bool_], upgraded: npt.NDArray[np.bool_] +) -> npt.NDArray[np.float64]: + """Seconds from each flagged row to the nearest flagged row of the other segment; NaN if unflagged.""" + times = index.as_unit("ns").asi8 + gap = np.full(len(index), np.nan) + gap[baseline] = _nearest_gap_s(times[baseline], others=times[upgraded]) + gap[upgraded] = _nearest_gap_s(times[upgraded], others=times[baseline]) + return gap + + +def _nearest_gap_s(times: npt.NDArray[np.int64], *, others: npt.NDArray[np.int64]) -> npt.NDArray[np.float64]: + if not len(others): + return np.full(len(times), np.inf) + others = np.sort(others) + position = np.searchsorted(others, times) + after = others[np.minimum(position, len(others) - 1)] + before = others[np.maximum(position - 1, 0)] + return np.minimum(np.abs(after - times), np.abs(times - before)) / 1e9 def _infer_timebase(index: pd.DatetimeIndex) -> pd.Timedelta: @@ -78,6 +243,29 @@ def _infer_timebase(index: pd.DatetimeIndex) -> pd.Timedelta: return pd.Timedelta(np.median(np.diff(unique.to_numpy()))) +def block_ids(index: pd.DatetimeIndex, *, start: pd.Timestamp, block: pd.Timedelta) -> npt.NDArray[np.int64]: + """Index of the fixed wall-clock block of length ``block``, tiled forward from ``start``, holding each row.""" + return np.floor(np.asarray((index - start) / block, dtype=float)).astype(np.int64) + + +def blocks_with_both_states( + ids: npt.NDArray[np.int64], *, baseline: npt.NDArray[np.bool_], upgraded: npt.NDArray[np.bool_] +) -> npt.NDArray[np.bool_]: + """Per row: whether its block holds at least one ``baseline`` and one ``upgraded`` row.""" + both = set(np.unique(ids[baseline])) & set(np.unique(ids[upgraded])) + return np.isin(ids, list(both)) + + +def _block_sum(ids: npt.NDArray[np.int64], values: npt.NDArray[np.float64], mask: npt.NDArray[np.bool_]) -> pd.Series: + """Sum of ``values`` over the ``mask`` rows of each block, indexed by block id.""" + return pd.Series(values[mask]).groupby(ids[mask]).sum() + + +def _per_row(ids: npt.NDArray[np.int64], per_block: pd.Series) -> npt.NDArray[np.float64]: + """Broadcast a per-block value back onto rows; NaN for a row whose block has no value.""" + return per_block.reindex(ids).to_numpy(dtype=float) + + def _wide_column(scada_df: pd.DataFrame, *, turbine_col: str, value_col: str) -> pd.DataFrame: """Pivot long SCADA to a timestamp x turbine table of ``value_col`` (NaN where missing).""" tmp = scada_df[[turbine_col, value_col]].copy() @@ -127,9 +315,26 @@ class ToggleSpecialistMethod: period**, for which :data:`DEFAULT_BLOCK_HOURS` may span only a cycle or two. :param n_resamples: bootstrap resamples; block sums are precomputed, so this can be generous. :param bootstrap_seed: RNG seed for the bootstrap, so a reported sigma is reproducible. + :param pairing_max_gap: when set, a used row of one segment is kept only if the other segment has + a used row within this gap, inclusive (a gap of exactly ``pairing_max_gap`` is kept). Must be + a positive whole multiple of the timebase; ``None`` disables pairing. + :param toggle_period: **required** one toggle cycle (the campaign's own on + off duration): the + block length for the block leg, the label permutation and ``columns.exclude_block``. A + positive whole multiple of the timebase. It is the caller's to state, not inferred, and it + is independent of ``pairing_max_gap``. + :param toggle_datum: **required** any boundary of the toggle-cycle grid (the controller's own + cycle datum, not when the campaign happened to begin); the cycles of the block leg, the label + permutation and ``columns.exclude_block`` are tiled from it. A ``ToggleSchedule`` tiles its + cycles from its first timestamp, which must lie on this grid. + :param bin_power: the per-row power the bins and ``power_band`` read (see the module docstring): + ``"baseline"`` (default) or ``"both_states"``. + :param power_band: when set, ``(lo, hi]`` in kW: each leg keeps only the used rows whose bin + power lies in it, and estimates its uplift over those. ``None`` keeps every used row. """ columns: ColumnSchema + toggle_period: pd.Timedelta + toggle_datum: pd.Timestamp name: str = "toggle_specialist" out_dir: Path | None = None save_plots: bool = False @@ -139,11 +344,22 @@ class ToggleSpecialistMethod: block_hours: float = DEFAULT_BLOCK_HOURS n_resamples: int = 1000 bootstrap_seed: int = 0 + pairing_max_gap: pd.Timedelta | None = None + bin_power: Literal["baseline", "both_states"] = "baseline" + power_band: tuple[float, float] | None = None def __post_init__(self) -> None: - """Validate ``columns`` names every role this method reads, and the requested ``conditions``.""" + """Validate the column roles, ``conditions``, ``bin_power`` and ``power_band``.""" self.columns.require_roles(("active_power", "availability")) validate_conditions(self.conditions, supported=_SUPPORTED_CONDITIONS, method_name=self.name) + if self.bin_power not in _BIN_POWERS: + msg = f"{self.name}: bin_power must be one of {_BIN_POWERS}, got {self.bin_power!r}" + raise ValueError(msg) + if self.power_band is not None: + lo, hi = self.power_band + if not (np.isfinite(lo) and np.isfinite(hi) and lo < hi): + msg = f"{self.name}: power_band must be finite (lo, hi] with lo < hi, got {self.power_band}" + raise ValueError(msg) if "power" in self.conditions and self.rated_power_kw is None: msg = ( f"{self.name}: rated_power_kw is required when 'power' is in conditions — the power bin " @@ -160,22 +376,6 @@ def estimate(self, mi: MethodInput) -> MethodOutput: f"or toggle_df; this method has no prepost baseline to compare against." ) raise ValueError(msg) - - mi = restrict_to_campaign(mi) - wide = _wide_column(mi.scada_df, turbine_col=mi.turbine_col, value_col=self.columns.active_power) - test = mi.test_wtg - refs = mi.context.references_among(wide.columns) - if not refs: - msg = ( - f"no reference turbines available for test_wtg {test!r}: scada_df contains only " - f"{list(wide.columns)}. The toggle specialist method needs at least one reference turbine." - ) - raise ValueError(msg) - # Narrow to the campaign's turbines and blank the cells it says may not contribute, so the - # estimate and every diagnostic below see one consistent selection. - ref_set = set(refs) - wide = mi.context.mask_invalid(wide[[c for c in wide.columns if c == test or c in ref_set]]) - if self.columns.availability not in mi.scada_df.columns: msg = ( f"the availability column {self.columns.availability!r} (columns.availability) is not in " @@ -183,91 +383,290 @@ def estimate(self, mi: MethodInput) -> MethodOutput: ) raise ValueError(msg) + mi = restrict_to_campaign(mi) + wide_all = _wide_column(mi.scada_df, turbine_col=mi.turbine_col, value_col=self.columns.active_power) + test = mi.test_wtg timebase = self.timebase if self.timebase is not None else _infer_timebase(mi.scada_df.index) - rows = resolve_toggle(mi.upgrade_timing, wide.index) + self._check_pairing_gap(timebase) + rows = resolve_toggle(mi.upgrade_timing, wide_all.index) + block = self._toggle_period(timebase) + datum = self._toggle_datum(mi, wide_all.index) + + components: dict[str, _Estimate] = {} + subsets: pd.DataFrame | None = None + n_refs_available = len(mi.context.references_among(wide_all.columns)) + if n_refs_available: + components[_REFERENCE_LEG], subsets = self._choose_references( + mi, wide_all=wide_all, test=test, timebase=timebase, block=block, datum=datum, rows=rows + ) + else: + # Nothing for the reference leg to work with: the estimate is the block leg alone, which + # needs no references, rather than a failed estimate. + logger.warning( + "%s: no reference turbines available for test_wtg %r; the estimate is the %r leg alone", + self.name, + test, + _BLOCK_LEG, + ) + components[_BLOCK_LEG] = self._component( + mi, + mode=_BLOCK_LEG, + wide_all=wide_all, + test=test, + timebase=timebase, + block=block, + datum=datum, + rows=rows, + ) + est = ( + _combine(components[_REFERENCE_LEG], components[_BLOCK_LEG], upgraded=rows.upgraded, bins=self._bins()) + if _REFERENCE_LEG in components + else components[_BLOCK_LEG] + ) + + ref_set = set(est.refs) + wide = wide_all[[c for c in wide_all.columns if c == test or c in ref_set]] + stats = pd.concat( + [ + _segment_stats( + mi, + wide=wide_all[[c for c in wide_all.columns if c == test or c in set(component.refs)]], + used=component.used, + toggle_rows=rows, + ref_total=component.ref_total, + timebase=timebase, + active_power_col=self.columns.active_power, + ).assign(component=mode) + for mode, component in components.items() + ], + ignore_index=True, + ) + self._write_outputs( + mi, + wide=wide, + stats=stats, + est=est, + components=components, + timebase=timebase, + rows=rows, + subsets=subsets, + n_refs_available=n_refs_available, + ) + return MethodOutput( + p50_overall=float(est.uplift), + p50_by_condition=est.per_bin, + sigma_overall=est.sigma_overall, + uncertainty_diagnostics=_harness_diagnostics(est.diagnostics), + labeled_rows=self._labeled_rows(mi, wide=wide, test=test, used=est.used, rows=rows, label=est.label), + selection_accounting=est.accounting, + ) + + def _component( + self, + mi: MethodInput, + *, + mode: str, + wide_all: pd.DataFrame, + test: str, + timebase: pd.Timedelta, + block: pd.Timedelta | None, + datum: pd.Timestamp, + rows: ToggleRowSets, + refs: list[str] | None = None, + with_actual: bool = True, + ) -> _Estimate: + """Run one leg end to end: row selection, reference, uplift, bins and the two bootstraps. + + ``refs`` names the references the ``sum`` leg uses (the block leg compares the test turbine + against its own block, so references play no part, not even in the row filters: a reference + outage must not cost it rows). ``with_actual=False`` skips the actual bootstrap (a subset + that is only being ranked on its permuted sigma); the label-permuted bootstrap always runs. + """ + refs = [] if mode == _BLOCK_LEG else list(refs or []) + if mode == _REFERENCE_LEG and not refs: + msg = f"{self.name}: the {mode!r} leg needs at least one reference turbine" + raise ValueError(msg) + # Narrow to the mode's turbines and blank the cells the campaign says may not contribute, so + # the estimate and every diagnostic see one consistent selection. + ref_set = set(refs) + wide = mi.context.mask_invalid(wide_all[[c for c in wide_all.columns if c == test or c in ref_set]]) + baseline = rows.campaign_baseline test_pw = wide[test].to_numpy(dtype=float) - ref_total = wide[refs].sum(axis=1).to_numpy(dtype=float) - used = self._used_mask(mi, wide=wide, test=test, refs=refs, timebase=timebase).to_numpy() + campaign = baseline | rows.upgraded + ids = block_ids(wide.index, start=datum, block=block) if block is not None and campaign.any() else None + selection = self._selection( + mi, wide=wide, test=test, refs=refs, timebase=timebase, rows=rows, ids=ids, both_states=mode == _BLOCK_LEG + ) + used = selection.blocked.to_numpy() + ref_total = self._reference_total(mode=mode, wide=wide, refs=refs, test_pw=test_pw, used=used, ids=ids) rho_base = _rho(test_pw, ref_total, used & baseline) rho_up = _rho(test_pw, ref_total, used & rows.upgraded) + label = self._bin_power( + mode=mode, + test_pw=test_pw, + ref_total=ref_total, + baseline=used & baseline, + upgraded=used & rows.upgraded, + ids=ids, + rho_base=rho_base, + rho_up=rho_up, + ) + in_band = None + if self.power_band is not None: + lo, hi = self.power_band + in_band = used & (label > lo) & (label <= hi) + used = in_band + rho_base = _rho(test_pw, ref_total, used & baseline) + rho_up = _rho(test_pw, ref_total, used & rows.upgraded) + accounting = self._selection_accounting(selection, rows=rows, power_band=in_band) recoverable = np.isfinite(rho_base) and rho_base != 0 and np.isfinite(rho_up) uplift = rho_up / rho_base - 1.0 if recoverable else np.nan - rho_label = _rho_label(rho_base, rho_up) per_bin = ( self._conditional_frame( test_pw=test_pw, ref_total=ref_total, - rho_label=rho_label, + label=label, baseline=used & baseline, upgraded=used & rows.upgraded, - ) + ).rename(columns={"sum_actual": f"sum_actual_{mode}", "sum_counterfactual": f"sum_counterfactual_{mode}"}) if "power" in self.conditions else None ) # Uncertainty runs strictly after the uplift, off the same frozen row selection and bin # assignment, and only when there is a finite uplift to qualify. - membership = self._cell_membership(rho_label=rho_label, ref_total=ref_total, used=used) - boot = ( - self._bootstrap( + membership = self._cell_membership(label=label, used=used) + + def _run_bootstrap(up: npt.NDArray[np.bool_], base: npt.NDArray[np.bool_]) -> BootstrapResult: + return self._bootstrap( index=wide.index, test_pw=test_pw, ref_total=ref_total, used=used, - upgraded=rows.upgraded, - baseline=baseline, + upgraded=up, + baseline=base, membership=membership, timebase=timebase, ) - if np.isfinite(uplift) - else None - ) + + boot = _run_bootstrap(rows.upgraded, baseline) if with_actual and np.isfinite(uplift) else None + decision_boot = None + if np.isfinite(uplift): + # Same rows, span, blocks, seed and cells as the actual bootstrap; only the labels differ. + up_p, base_p = permute_labels_within_blocks( + upgraded=rows.upgraded, + baseline=baseline, + used=used, + ids=ids, + rng=np.random.default_rng(self.bootstrap_seed), + ) + decision_boot = _run_bootstrap(up_p, base_p) + uplift_by_cell = {_OVERALL: uplift} if per_bin is not None: per_bin["sigma_uplift"] = [_cell_sigma(boot, str(b)) for b in per_bin["condition_bin"]] + uplift_by_cell.update( + {str(b): float(u) for b, u in zip(per_bin["condition_bin"], per_bin["p50_uplift"], strict=True)} + ) diagnostics = _uncertainty_diagnostics( boot, membership=membership, upgraded=used & rows.upgraded, baseline=used & baseline, used=used, + uplift_by_cell=uplift_by_cell, ) - - stats = _segment_stats( - mi, - wide=wide, - used=used, - toggle_rows=rows, + diagnostics.insert(0, "component", mode) + return _Estimate( + mode=mode, refs=refs, - timebase=timebase, - active_power_col=self.columns.active_power, - ) - sigma_overall = _cell_sigma(boot, _OVERALL) - self._write_outputs( - mi, - wide=wide, - stats=stats, + selection=selection, used=used, + ref_total=ref_total, + label=label, rho_base=rho_base, rho_up=rho_up, uplift=uplift, - sigma_overall=sigma_overall, - n_refs=len(refs), - timebase=timebase, + sigma_overall=_cell_sigma(boot, _OVERALL), per_bin=per_bin, + boot=boot, + accounting=accounting, diagnostics=diagnostics, + decision_boot=decision_boot, ) - return MethodOutput( - p50_overall=float(uplift), - p50_by_condition=per_bin, - sigma_overall=sigma_overall, - uncertainty_diagnostics=diagnostics, - labeled_rows=self._labeled_rows( - mi, wide=wide, test=test, used=used, rows=rows, rho_label=rho_label, ref_total=ref_total - ), + + def _choose_references( + self, + mi: MethodInput, + *, + wide_all: pd.DataFrame, + test: str, + timebase: pd.Timedelta, + block: pd.Timedelta | None, + datum: pd.Timestamp, + rows: ToggleRowSets, + ) -> tuple[_Estimate, pd.DataFrame | None]: + """Run the ``sum`` leg on every non-empty reference subset; return the smallest-sigma one. + + The sigma ranked is the ``decision_sigma`` of the label-permuted bootstrap; the subsets skip + their actual bootstrap and the winner is re-run in full. Ties go to the larger subset, then + to the earlier in enumeration (the context's reference order). No finite sigma anywhere -> + the full set. The second value lists every subset's result with the winner marked. + """ + refs_all = mi.context.references_among(wide_all.columns) + if len(refs_all) > _MAX_REFS_FOR_SUBSETS: + msg = ( + f"{self.name}: {len(refs_all)} references would need {2 ** len(refs_all) - 1} subset estimates; " + f"pass at most {_MAX_REFS_FOR_SUBSETS} (the nearest or best-correlated) references" + ) + raise ValueError(msg) + subsets = [list(c) for k in range(1, len(refs_all) + 1) for c in combinations(refs_all, k)] + estimates = [ + self._component( + mi, + mode=_REFERENCE_LEG, + wide_all=wide_all, + test=test, + timebase=timebase, + block=block, + datum=datum, + rows=rows, + refs=refs, + with_actual=False, + ) + for refs in subsets + ] + + def _rank(i: int) -> tuple[float, int, int]: + sigma = estimates[i].decision_sigma + return (sigma if np.isfinite(sigma) else np.inf, -len(subsets[i]), i) + + best = min(range(len(subsets)), key=_rank) + estimates[best] = self._component( + mi, + mode=_REFERENCE_LEG, + wide_all=wide_all, + test=test, + timebase=timebase, + block=block, + datum=datum, + rows=rows, + refs=subsets[best], ) + table = pd.DataFrame( + { + "refs": [";".join(refs) for refs in subsets], + "n_refs": [len(refs) for refs in subsets], + "n_used_timestamps": [int(est.used.sum()) for est in estimates], + "uplift_frc": [est.uplift for est in estimates], + "uplift_sigma_frc": [est.sigma_overall for est in estimates], + "decision_sigma_frc": [est.decision_sigma for est in estimates], + "chosen": [i == best for i in range(len(subsets))], + } + ) + return estimates[best], table def _labeled_rows( self, @@ -277,8 +676,7 @@ def _labeled_rows( test: str, used: npt.NDArray[np.bool_], rows: ToggleRowSets, - rho_label: float, - ref_total: npt.NDArray[np.float64], + label: npt.NDArray[np.float64], ) -> pd.DataFrame: """Return the test turbine's own records, tagged with the labels this estimate was built from. @@ -295,20 +693,19 @@ def _on_test_rows(values: npt.NDArray[np.generic]) -> npt.NDArray[np.generic]: np.where(rows.upgraded, _UPGRADED, np.where(rows.campaign_baseline, _BASELINE, _EXCLUDED)) ) - # The bin label is the same reference-derived baseline power the uplift binned on, so a row + # The bin label is the same bin power the uplift binned on, so a row # cannot sit in one bin here and another there. Outside the outer edges pd.cut gives NaN, # which is carried through as "this row belongs to no bin" rather than clipped to an edge. - if "power" in self.conditions and np.isfinite(rho_label): + if "power" in self.conditions and np.isfinite(label).any(): assert self.rated_power_kw is not None # noqa: S101 - guaranteed by __post_init__ bins = condition_bins("power", rated_power_kw=self.rated_power_kw) - labeled["power_bin"] = _on_test_rows(np.asarray(pd.cut(rho_label * ref_total, bins=bins))) + labeled["power_bin"] = _on_test_rows(np.asarray(pd.cut(label, bins=bins))) return labeled def _cell_membership( self, *, - rho_label: float, - ref_total: npt.NDArray[np.float64], + label: npt.NDArray[np.float64], used: npt.NDArray[np.bool_], ) -> dict[str, npt.NDArray[np.bool_]]: """Which **used** records belong to each bootstrap cell: the headline, plus each power bin. @@ -318,11 +715,11 @@ def _cell_membership( """ used_idx = np.flatnonzero(used) membership: dict[str, npt.NDArray[np.bool_]] = {_OVERALL: np.ones(len(used_idx), dtype=bool)} - if "power" not in self.conditions or not np.isfinite(rho_label): + if "power" not in self.conditions or not np.isfinite(label[used_idx]).any(): return membership assert self.rated_power_kw is not None # noqa: S101 - guaranteed by __post_init__ bins = condition_bins("power", rated_power_kw=self.rated_power_kw) - assigned = pd.cut(rho_label * ref_total[used_idx], bins=bins) + assigned = pd.cut(label[used_idx], bins=bins) for category in assigned.categories: membership[str(category)] = np.asarray(assigned == category) return membership @@ -367,28 +764,27 @@ def _conditional_frame( *, test_pw: npt.NDArray[np.float64], ref_total: npt.NDArray[np.float64], - rho_label: float, + label: npt.NDArray[np.float64], baseline: npt.NDArray[np.bool_], upgraded: npt.NDArray[np.bool_], ) -> pd.DataFrame: - """Per-power-bin uplift: ``rho_up(b) / rho_base(b) - 1``, on bins of the mean operating point. + """Per-power-bin uplift: ``rho_up(b) / rho_base(b) - 1``, on bins of ``label``. Two decisions carry this, and both are needed: - **The bin label is** ``rho_label * ref_total`` (see :func:`_rho_label`): reference-derived and - state-neutral, so neither the upgrade nor which state is called baseline can move a row + **The bin label is the bin power** (:meth:`_bin_power`): a row's own state cannot move it between bins, and it is on the test turbine's own kW scale. **The denominator is the per-bin** ``rho_base(b)``, not the global one: the test-to-reference ratio varies with power, and a global denominator would read that structure as uplift. The price is that the per-bin numbers no longer aggregate exactly to ``p50_overall``, which is - deliberate and un-relevelled; ``sum_actual`` / ``sum_counterfactual`` expose the gap. + deliberate and un-relevelled; the leg's ``sum_actual_`` / ``sum_counterfactual_`` + expose the gap. Sparse bins report NaN with ``n_records = 0`` rather than being imputed. """ assert self.rated_power_kw is not None # noqa: S101 - guaranteed by __post_init__ bins = condition_bins("power", rated_power_kw=self.rated_power_kw) - label = rho_label * ref_total counterfactual = _per_bin_counterfactual( label=label, test_pw=test_pw, ref_total=ref_total, baseline=baseline, bins=bins ) @@ -396,59 +792,240 @@ def _conditional_frame( frame.insert(0, "condition", "power") return frame - def _used_mask( - self, mi: MethodInput, *, wide: pd.DataFrame, test: str, refs: list[str], timebase: pd.Timedelta - ) -> pd.Series: - """Complete-case timestamps that also pass downtime filtering on the test turbine and every reference. + def _bins(self) -> list[float] | None: + """Return the power bin edges when the method reports per bin; None otherwise.""" + if "power" not in self.conditions: + return None + assert self.rated_power_kw is not None # noqa: S101 - guaranteed by __post_init__ + return condition_bins("power", rated_power_kw=self.rated_power_kw) - Returns a bool Series on ``wide.index``. Every turbine (test and references) must be - available (counter >= a full period) and have finite power — a down turbine on either side - of the ratio is therefore excluded. The test turbine additionally goes through the shared - :class:`NormalOperationFilter` (the same downtime + finite-power logic the power model uses; - the stuck filter is left off here as the ratio sums raw power rather than fitting a model). - """ - turbines = [test, *refs] - complete = wide[turbines].notna().all(axis=1) + def _check_pairing_gap(self, timebase: pd.Timedelta) -> None: + gap = self.pairing_max_gap + if gap is None: + return + if gap <= pd.Timedelta(0) or gap % timebase != pd.Timedelta(0): + msg = ( + f"{self.name}: pairing_max_gap {gap} must be a positive whole multiple of the timebase " + f"{timebase}; a shorter gap can never reach a row of the other segment." + ) + raise ValueError(msg) - full = timebase.total_seconds() - avail = _wide_column(mi.scada_df, turbine_col=mi.turbine_col, value_col=self.columns.availability).reindex( - index=wide.index, columns=turbines - ) - all_available = (avail >= full).all(axis=1) + def _toggle_period(self, timebase: pd.Timedelta) -> pd.Timedelta: + """Return ``toggle_period`` after checking it tiles the timebase grid.""" + block = self.toggle_period + ratio = block / timebase + if block <= pd.Timedelta(0) or ratio != round(ratio): + msg = ( + f"{self.name}: toggle_period {block} must be a positive whole multiple of the timebase " + f"{timebase}; a period that does not tile the timebase grid cannot hold whole records." + ) + raise ValueError(msg) + return block - test_rows = mi.scada_df[mi.scada_df[mi.turbine_col] == test] - test_keep = ( - NormalOperationFilter( - active_power_col=self.columns.active_power, - availability_col=self.columns.availability, - apply_stuck_filter=False, + def _toggle_datum(self, mi: MethodInput, index: pd.DatetimeIndex) -> pd.Timestamp: + """Return ``toggle_datum`` after checking it is comparable with ``index`` and fits a ``ToggleSchedule``.""" + datum = pd.Timestamp(self.toggle_datum) + if (datum.tz is None) != (index.tz is None): + msg = f"{self.name}: toggle_datum {datum} and the data index must both be tz-aware or both naive" + raise ValueError(msg) + timing = mi.upgrade_timing + if isinstance(timing, ToggleSchedule) and len(index): + first = timing.start if timing.start is not None else index.min() + if (pd.Timestamp(first) - datum) % self.toggle_period != pd.Timedelta(0): + msg = ( + f"{self.name}: the ToggleSchedule tiles its cycles from {first}, which is not on the " + f"toggle_datum {datum} grid of {self.toggle_period}" + ) + raise ValueError(msg) + return datum + + def _bin_power( + self, + *, + mode: str, + test_pw: npt.NDArray[np.float64], + ref_total: npt.NDArray[np.float64], + baseline: npt.NDArray[np.bool_], + upgraded: npt.NDArray[np.bool_], + ids: npt.NDArray[np.int64] | None, + rho_base: float, + rho_up: float, + ) -> npt.NDArray[np.float64]: + """Per row, the power the bins and ``power_band`` read (``bin_power``); NaN where undefined.""" + if mode == _REFERENCE_LEG: + scale = rho_base if self.bin_power == "baseline" else 0.5 * (rho_base + rho_up) + return scale * ref_total + if ids is None: + return np.full(len(test_pw), np.nan) + + def _cycle_mean(mask: npt.NDArray[np.bool_]) -> npt.NDArray[np.float64]: + n = pd.Series(mask.astype(float)).groupby(ids).sum() + return _per_row(ids, _block_sum(ids, test_pw, mask) / n[n > 0]) + + if self.bin_power == "baseline": + return _cycle_mean(baseline) + return 0.5 * (_cycle_mean(baseline) + _cycle_mean(upgraded)) + + def _reference_total( + self, + *, + mode: str, + wide: pd.DataFrame, + refs: list[str], + test_pw: npt.NDArray[np.float64], + used: npt.NDArray[np.bool_], + ids: npt.NDArray[np.int64] | None, + ) -> npt.NDArray[np.float64]: + """Return the per-row reference power for ``mode``; NaN on rows no surviving block covers. + + The block leg gives every used row of a block the block's mean test power, so within a block + the on/off ratio cancels the block's wind level whatever the on/off count mix (the block mean + carries the uplift at ``n_on / n`` strength in both denominators, a second-order effect). + """ + if mode == _REFERENCE_LEG: + return wide[refs].sum(axis=1).to_numpy(dtype=float) + if ids is None or not used.any(): + return np.full(len(wide), np.nan) + block_test = _block_sum(ids, test_pw, used) + n_used = pd.Series(used.astype(float)).groupby(ids).sum() + return _per_row(ids, block_test / n_used.reindex(block_test.index)) + + def _selection( + self, + mi: MethodInput, + *, + wide: pd.DataFrame, + test: str, + refs: list[str], + timebase: pd.Timedelta, + rows: ToggleRowSets, + ids: npt.NDArray[np.int64] | None = None, + both_states: bool = True, + ) -> _Selection: + """Return the used-row mask after each selection stage, each a bool Series on ``wide.index``. + + ``ids`` are the toggle-cycle ids (``None`` when there are no campaign rows to tile). The + ``block`` stage drops every row of a block that carries a ``columns.exclude_block`` flag on + any of the turbines (both legs) or, when ``both_states`` (the block leg), lacks either state. + """ + turbines = [test, *refs] + filtered = self._filtered_mask(mi, wide=wide, test=test, refs=refs, timebase=timebase) + not_excluded = filtered & ~self._excluded(mi, turbines=turbines, index=wide.index) + in_segment = rows.campaign_baseline | rows.upgraded + kept = not_excluded.to_numpy() + if self.pairing_max_gap is not None: + paired = pair_within( + wide.index, + baseline=kept & rows.campaign_baseline, + upgraded=kept & rows.upgraded, + max_gap=self.pairing_max_gap, ) - .keep_mask(test_rows, timebase=timebase) - .reindex(wide.index, fill_value=False) + kept = kept & (~in_segment | paired.baseline | paired.upgraded) + paired_series = pd.Series(kept, index=wide.index) + block_flag = self._excluded(mi, turbines=turbines, index=wide.index, role="exclude_block").to_numpy() + if block_flag.any(): + if ids is None: + msg = ( + f"{self.name}: columns.exclude_block {self.columns.exclude_block!r} carries flags but there is " + f"no campaign row to tile into blocks." + ) + raise ValueError(msg) + kept = kept & ~np.isin(ids, np.unique(ids[block_flag])) + if ids is not None and both_states: + both = blocks_with_both_states(ids, baseline=kept & rows.campaign_baseline, upgraded=kept & rows.upgraded) + kept = kept & (~in_segment | both) + return _Selection( + filtered=filtered, + not_excluded=not_excluded, + paired=paired_series, + blocked=pd.Series(kept, index=wide.index), ) - return complete & all_available & test_keep & ~self._test_excluded(mi, test=test, index=wide.index) - def _test_excluded(self, mi: MethodInput, *, test: str, index: pd.DatetimeIndex) -> pd.Series: - """Boolean mask on *index*: the test turbine's caller-flagged rows to drop (empty when unused). + def _selection_accounting( + self, selection: _Selection, *, rows: ToggleRowSets, power_band: npt.NDArray[np.bool_] | None = None + ) -> pd.DataFrame: + """Rows kept per segment after each stage. - Reads the ``columns.exclude_row`` column of the test turbine's rows; references are never - excluded here (their special modes still carry information for the ratio). Absent column or - unset role -> nothing excluded. Reindex fills missing timestamps with ``False`` so an expanded - index never becomes an exclusion. NaN raises rather than coercing: ``astype(bool)`` reads a - missing flag as ``True`` and drops the row, the opposite of the safe default. + ``kept_fraction`` is relative to the segment's rows; ``stage_kept_fraction`` to the previous + stage, so it isolates which stage discriminates between segments. + """ + stages_masks = ( + ("segment", None), + ("filters", selection.filtered), + ("exclude_row", selection.not_excluded), + ("pairing", selection.paired), + ("block", selection.blocked), + *(() if power_band is None else (("power_band", pd.Series(power_band, index=selection.blocked.index)),)), + ) + records = [] + for segment, in_segment in ((_BASELINE, rows.campaign_baseline), (_UPGRADED, rows.upgraded)): + n_segment = int(in_segment.sum()) + previous = n_segment + for stage, mask in stages_masks: + n_kept = n_segment if mask is None else int((mask.to_numpy() & in_segment).sum()) + records.append( + { + "stage": stage, + "segment": segment, + "n_kept": n_kept, + "kept_fraction": n_kept / n_segment if n_segment else np.nan, + "stage_kept_fraction": n_kept / previous if previous else np.nan, + } + ) + previous = n_kept + return pd.DataFrame(records) + + def _filtered_mask( + self, mi: MethodInput, *, wide: pd.DataFrame, test: str, refs: list[str], timebase: pd.Timedelta + ) -> pd.Series: + """Timestamps at which the test turbine and every reference pass the same normal-operation filter. + + Returns a bool Series on ``wide.index``. One :class:`NormalOperationFilter` (the downtime + + finite-power logic the power model uses; the stuck filter is left off here as the ratio sums + raw power rather than fitting a model) is applied to each turbine's own rows, so a turbine + that is down or has no power on either side of the ratio excludes the timestamp. The same + filter on every turbine means any future change to it (e.g. the stuck filter) applies to + the references exactly as to the test turbine. """ - col = self.columns.exclude_row + normal = NormalOperationFilter( + active_power_col=self.columns.active_power, + availability_col=self.columns.availability, + apply_stuck_filter=False, + ) + keep = pd.Series(data=True, index=wide.index, dtype=bool) + for turbine in (test, *refs): + rows = mi.scada_df[mi.scada_df[mi.turbine_col] == turbine] + keep &= normal.keep_mask(rows, timebase=timebase).reindex(wide.index, fill_value=False) + # ``wide`` is the masked pivot: a cell the campaign context blanked must not count either. + keep &= wide[turbine].notna() + return keep + + def _excluded( + self, mi: MethodInput, *, turbines: list[str], index: pd.DatetimeIndex, role: str = "exclude_row" + ) -> pd.Series: + """Boolean mask on *index*: timestamps at which any of ``turbines`` carries a caller-set ``role`` flag. + + Reads the ``columns.`` column (``exclude_row`` or ``exclude_block``) of each listed + turbine's rows (the test turbine and the references the mode uses), so a flagged reference + row counts the way a reference outage does. Absent column or unset role -> nothing flagged. + Reindex fills missing timestamps with ``False`` so an expanded index never becomes an + exclusion. NaN raises rather than coercing: ``astype(bool)`` reads a missing flag as ``True`` + and drops the row, the opposite of the safe default. + """ + col: str | None = getattr(self.columns, role) + excluded = pd.Series(data=False, index=index, dtype=bool) if not col or col not in mi.scada_df.columns: - return pd.Series(data=False, index=index, dtype=bool) - test_rows = mi.scada_df[mi.scada_df[mi.turbine_col] == test] - flags = test_rows[col] - if flags.isna().any(): - msg = ( - f"exclude_row column {col!r} has {int(flags.isna().sum())} NaN value(s) for turbine {test!r}; " - f"it must be boolean with no missing values (fill unknown rows with False explicitly)" - ) - raise ValueError(msg) - return flags.astype(bool).reindex(index, fill_value=False) + return excluded + for turbine in turbines: + flags = mi.scada_df.loc[mi.scada_df[mi.turbine_col] == turbine, col] + if flags.isna().any(): + msg = ( + f"{role} column {col!r} has {int(flags.isna().sum())} NaN value(s) for turbine " + f"{turbine!r}; it must be boolean with no missing values (fill unknown rows with False explicitly)" + ) + raise ValueError(msg) + excluded |= flags.astype(bool).reindex(index, fill_value=False) + return excluded def _write_outputs( self, @@ -456,17 +1033,14 @@ def _write_outputs( *, wide: pd.DataFrame, stats: pd.DataFrame, - used: np.ndarray, - rho_base: float, - rho_up: float, - uplift: float, - sigma_overall: float, - n_refs: int, + est: _Estimate, + components: dict[str, _Estimate], timebase: pd.Timedelta, - per_bin: pd.DataFrame | None = None, - diagnostics: pd.DataFrame | None = None, + rows: ToggleRowSets, + subsets: pd.DataFrame | None = None, + n_refs_available: int | None = None, ) -> None: - """Write the data-stats CSV, the headline results CSV, the per-bin CSV and (optionally) the plots.""" + """Write the data-stats, results, per-bin, uncertainty, selection and reference-subsets CSVs and the plots.""" upgrade_start = toggle_upgrade_start(mi.upgrade_timing, wide.index) last_dt = wide.index.max() run_name = f"toggle_specialist_{mi.test_wtg}_{upgrade_start:%Y%m%d}_{last_dt:%Y%m%d}" @@ -487,11 +1061,17 @@ def _write_outputs( "test_wtg": mi.test_wtg, "mode": "toggle", "n_turbines": wide.shape[1], - "n_refs": n_refs, - "ratio_baseline": rho_base, - "ratio_upgraded": rho_up, - "uplift_frc": uplift, - "uplift_sigma_frc": sigma_overall, + "n_refs": len(est.refs), + "n_refs_available": len(est.refs) if n_refs_available is None else n_refs_available, + "refs_used": ";".join(est.refs), + "uplift_frc": est.uplift, + "uplift_sigma_frc": est.sigma_overall, + **{f"ratio_baseline_{mode}": component.rho_base for mode, component in components.items()}, + **{f"ratio_upgraded_{mode}": component.rho_up for mode, component in components.items()}, + f"weight_{_REFERENCE_LEG}": _overall_weight(est, components), + **{f"uplift_frc_{mode}": component.uplift for mode, component in components.items()}, + **{f"uplift_sigma_frc_{mode}": component.sigma_overall for mode, component in components.items()}, + **{f"decision_sigma_frc_{mode}": c.decision_sigma for mode, c in components.items()}, "block_hours": self.block_hours, "n_resamples": self.n_resamples, "n_used_timestamps_baseline": used_base, @@ -502,40 +1082,65 @@ def _write_outputs( ) results.to_csv(run_dir / f"{run_name}_results_{ts}.csv", index=False) - if per_bin is not None: - per_bin.to_csv(run_dir / f"{run_name}_by_power_bin_{ts}.csv", index=False) - if diagnostics is not None: - diagnostics.to_csv(run_dir / f"{run_name}_uncertainty_{ts}.csv", index=False) + if est.per_bin is not None: + est.per_bin.to_csv(run_dir / f"{run_name}_by_power_bin_{ts}.csv", index=False) + est.diagnostics.to_csv(run_dir / f"{run_name}_uncertainty_{ts}.csv", index=False) + if subsets is not None: + subsets.to_csv(run_dir / f"{run_name}_reference_subsets_{ts}.csv", index=False) + # Every leg's row selection, so the blend's two selections can be told apart. + pd.concat( + [component.accounting.assign(component=mode) for mode, component in components.items()], ignore_index=True + ).to_csv(run_dir / f"{run_name}_selection_{ts}.csv", index=False) if self.save_plots: + # Scatter and ratio plots of one leg's own rows against its own reference. + shown = components.get(_REFERENCE_LEG, components[_BLOCK_LEG]) _save_plots( run_dir / "plots", wide=wide, mi=mi, test=mi.test_wtg, - used=used, + used=shown.used, + ref_total=shown.ref_total, + reference_mode=shown.mode, timebase=timebase, active_power_col=self.columns.active_power, ) - if per_bin is not None: + if est.per_bin is not None: _save_per_bin_plot( run_dir / "plots" / stages.CONDITIONAL_UPLIFT / f"{mi.test_wtg}_per_bin_uplift.png", - per_bin=per_bin, + per_bin=est.per_bin, test=mi.test_wtg, active_power_col=self.columns.active_power, ) - self._write_shared_diagnostics(mi, run_dir=run_dir, wide=wide, timebase=timebase) + _save_pairing_gap_plot( + run_dir / "plots" / stages.FILTER / f"{mi.test_wtg}_pairing_gap.png", + index=wide.index, + selection=est.selection, + rows=rows, + timebase=timebase, + max_gap=self.pairing_max_gap, + test=mi.test_wtg, + ) + self._write_shared_diagnostics( + mi, run_dir=run_dir, wide=wide, used=est.used, timebase=timebase, turbines=[mi.test_wtg, *est.refs] + ) def _write_shared_diagnostics( - self, mi: MethodInput, *, run_dir: Path, wide: pd.DataFrame, timebase: pd.Timedelta + self, + mi: MethodInput, + *, + run_dir: Path, + wide: pd.DataFrame, + used: np.ndarray, + timebase: pd.Timedelta, + turbines: list[str], ) -> None: """Emit the shared cross-method diagnostics (coverage/curves/histograms) and the run config.""" # ``wide`` (a pivot) drops all-NaN timestamps, so align the masks to the full unique index # the DiagnosticContext uses (timestamps absent from ``wide`` are simply not used). index = pd.DatetimeIndex(pd.unique(mi.scada_df.index)).sort_values() - test, refs = mi.test_wtg, [c for c in wide.columns if c != mi.test_wtg] - used_series = self._used_mask(mi, wide=wide, test=test, refs=refs, timebase=timebase) - used = used_series.reindex(index, fill_value=False).to_numpy() + used = pd.Series(used, index=wide.index).reindex(index, fill_value=False).to_numpy() treated = resolve_toggle(mi.upgrade_timing, index).upgraded.astype(bool) ctx = DiagnosticContext( run_dir=run_dir, @@ -549,12 +1154,18 @@ def _write_shared_diagnostics( mode="toggle", era5_df=None, # the exclusion alone, so the plots show what the flag removed that downtime did not - excluded_ts=self._test_excluded(mi, test=test, index=index).to_numpy(), + excluded_ts=self._excluded(mi, turbines=turbines, index=index).to_numpy(), ) write_common_diagnostics(ctx) params = { "active_power_col": self.columns.active_power, "availability_col": self.columns.availability, + "pairing_max_gap": None if self.pairing_max_gap is None else str(self.pairing_max_gap), + "toggle_period": str(self.toggle_period), + "toggle_datum": str(self.toggle_datum), + "bin_power": self.bin_power, + "power_band": None if self.power_band is None else list(self.power_band), + "exclude_block_col": self.columns.exclude_block, } write_run_config(ctx, method_name=self.name, method_params=params) @@ -582,6 +1193,106 @@ def _per_bin_counterfactual( return rho_row * ref_total +def _combine(a: _Estimate, b: _Estimate, *, upgraded: npt.NDArray[np.bool_], bins: list[float] | None) -> _Estimate: + """Blend the two legs' estimates by minimum-variance weights, headline and per bin. + + The weights (and the correlation they need) come from the two legs' label-permuted bootstraps, + cell by cell; the blend's sigma comes from the actual ones. + + The blend's used rows are the union of the components' (a row contributed to at least one + estimate); its bin label is ``a``'s where ``a`` used the row and ``b``'s otherwise, so the two + components' bins, which share edges, are merged by name, and each bin's ``n_records`` counts + the used ``upgraded`` rows that merged label puts in it (what ``labeled_rows`` reproduces). Each + bin also carries both legs' energy sums and its weight on ``a``, so its blend is re-derivable. + Row-selection accounting, ``rho`` and the reference-side diagnostics are ``a``'s (the reference + leg); ``b``'s selection is written to the selection CSV. The diagnostics carry both legs' rows and + a ``combined`` row per cell with the blend's sigma, the weight on ``a`` and the correlation it used. + """ + if a.boot is not None and b.boot is not None and a.boot.n_blocks != b.boot.n_blocks: + # Both bootstraps come from the same method instance (same span, timebase, block length, + # resample count and seed), so their draws pair up; a block-count mismatch means they do not. + msg = f"component bootstraps drew different block counts ({a.boot.n_blocks} vs {b.boot.n_blocks})" + raise ValueError(msg) + + def _resamples(boot: BootstrapResult | None, cell: str) -> npt.NDArray[np.float64]: + return boot.resamples.get(cell, np.array([])) if boot is not None else np.array([]) + + def _cell(cell: str, est_a: float, sig_a: float, est_b: float, sig_b: float) -> CombinedEstimate: + weighting = ( + _cell_sigma(a.decision_boot, cell), + _cell_sigma(b.decision_boot, cell), + _resamples(a.decision_boot, cell), + _resamples(b.decision_boot, cell), + ) + return combine_estimates( + (est_a, sig_a, _resamples(a.boot, cell)), (est_b, sig_b, _resamples(b.boot, cell)), weighting=weighting + ) + + overall = _cell(_OVERALL, a.uplift, a.sigma_overall, b.uplift, b.sigma_overall) + blended_cells = {_OVERALL: overall} + + used = a.used | b.used + label = np.where(a.used, a.label, b.label) + per_bin = None + if a.per_bin is not None and b.per_bin is not None and bins is not None: + fa = a.per_bin.set_index(a.per_bin["condition_bin"].astype(str)) + fb = b.per_bin.set_index(b.per_bin["condition_bin"].astype(str)) + upgraded_bin = pd.Series(pd.cut(label[used & upgraded], bins=bins).astype(str)) + records = [] + for cell in fa.index: + ra, rb = fa.loc[cell], fb.loc[cell] + blend = _cell(cell, ra["p50_uplift"], ra["sigma_uplift"], rb["p50_uplift"], rb["sigma_uplift"]) + blended_cells[cell] = blend + records.append( + { + "condition": ra["condition"], + "condition_bin": ra["condition_bin"], + "p50_uplift": blend.estimate, + "n_records": int((upgraded_bin == cell).sum()), + **{f"{col}_{a.mode}": ra[f"{col}_{a.mode}"] for col in ("sum_actual", "sum_counterfactual")}, + **{f"{col}_{b.mode}": rb[f"{col}_{b.mode}"] for col in ("sum_actual", "sum_counterfactual")}, + f"weight_{a.mode}": blend.weight_a, + "sigma_uplift": blend.sigma, + } + ) + per_bin = pd.DataFrame(records) + + combined_rows = pd.DataFrame( + [ + { + "component": _COMBINED, + "condition": _OVERALL if cell == _OVERALL else "power", + "condition_bin": cell, + "uplift": blend.estimate, + "sigma": blend.sigma, + f"weight_{a.mode}": blend.weight_a, + "correlation": blend.correlation, + f"decision_sigma_{a.mode}": _cell_sigma(a.decision_boot, cell), + f"decision_sigma_{b.mode}": _cell_sigma(b.decision_boot, cell), + } + for cell, blend in blended_cells.items() + ] + ) + diagnostics = pd.concat([a.diagnostics, b.diagnostics, combined_rows], ignore_index=True) + + return _Estimate( + mode=_COMBINED, + refs=a.refs, + selection=a.selection, + used=used, + ref_total=a.ref_total, + label=label, + rho_base=a.rho_base, + rho_up=a.rho_up, + uplift=overall.estimate, + sigma_overall=overall.sigma, + per_bin=per_bin, + boot=None, + accounting=a.accounting, + diagnostics=diagnostics, + ) + + def _cell_sigma(boot: BootstrapResult | None, cell: str) -> float: """Return one cell's 1-sigma, or NaN when the bootstrap did not run or never saw that cell.""" if boot is None or cell not in boot.cells: @@ -589,6 +1300,15 @@ def _cell_sigma(boot: BootstrapResult | None, cell: str) -> float: return boot.cells[cell].sigma +def _harness_diagnostics(diagnostics: pd.DataFrame) -> pd.DataFrame: + """One row per ``(condition, condition_bin)``, each leg's columns suffixed ``_``.""" + keys = ["condition", "condition_bin"] + values = [c for c in diagnostics.columns if c not in (*keys, "component")] + wide = diagnostics.pivot(index=keys, columns="component", values=values).dropna(axis=1, how="all") # noqa: PD010 - unique keys, nothing to aggregate + wide.columns = [f"{value}_{component}" for value, component in wide.columns] + return wide.reset_index() + + def _uncertainty_diagnostics( boot: BootstrapResult | None, *, @@ -596,8 +1316,9 @@ def _uncertainty_diagnostics( upgraded: npt.NDArray[np.bool_], baseline: npt.NDArray[np.bool_], used: npt.NDArray[np.bool_], + uplift_by_cell: Mapping[str, float] | None = None, ) -> pd.DataFrame: - """Per-cell account of how the uncertainty was reached, keyed by ``(condition, condition_bin)``. + """Per-cell account of the estimate and how its uncertainty was reached, keyed by ``(condition, condition_bin)``. Carried through the harness seam uninterpreted, so an uncertainty model can be developed against a saved sweep rather than by re-running one. Emitted even when the bootstrap did not run: the @@ -615,9 +1336,11 @@ def _uncertainty_diagnostics( { "condition": _OVERALL if cell == _OVERALL else "power", "condition_bin": cell, + "uplift": uplift_by_cell.get(cell, nan) if uplift_by_cell is not None else nan, "n_upgraded_records": int((member & up_used).sum()), "n_baseline_records": int((member & base_used).sum()), "n_blocks": boot.n_blocks if boot is not None else 0, + "sigma": cell_boot.sigma if cell_boot is not None else nan, # Both components, not just the reported max: a blend rule can then be re-judged from # a saved sweep rather than by re-running one. "sigma_bootstrap": cell_boot.sigma_bootstrap if cell_boot is not None else nan, @@ -639,13 +1362,15 @@ def _rho(test_pw: npt.NDArray[np.float64], ref_total: npt.NDArray[np.float64], m return float(test_pw[mask].sum() / denom) -def _rho_label(rho_base: float, rho_up: float) -> float: - """Return the test-to-reference ratio used to *label* bins: the mean of the two states. - - State-neutral by construction, so relabelling which state is the baseline cannot move a row - between bins. Still a campaign-level scalar, so the upgrade cannot move a row either. - """ - return 0.5 * (rho_base + rho_up) +def _overall_weight(est: _Estimate, components: dict[str, _Estimate]) -> float: + """Return the headline blend's weight on the reference leg: 1 with that leg alone, 0 without it.""" + if _REFERENCE_LEG not in components: + return 0.0 + if _BLOCK_LEG not in components: + return 1.0 + d = est.diagnostics + row = d[(d["component"] == _COMBINED) & (d["condition_bin"] == _OVERALL)] + return float(row[f"weight_{_REFERENCE_LEG}"].iloc[0]) if len(row) else float("nan") def _segment_stats( @@ -654,24 +1379,25 @@ def _segment_stats( wide: pd.DataFrame, used: npt.NDArray[np.bool_], toggle_rows: ToggleRowSets, - refs: list[str], + ref_total: npt.NDArray[np.float64], timebase: pd.Timedelta, active_power_col: str, ) -> pd.DataFrame: """Build the per-segment (all/baseline/upgraded) diagnostics table.""" test = mi.test_wtg test_pw = wide[test].to_numpy(dtype=float) - ref_total = wide[refs].sum(axis=1).to_numpy(dtype=float) n_turbines = wide.shape[1] timebase_hours = timebase / pd.Timedelta(hours=1) row_rows = resolve_toggle(mi.upgrade_timing, mi.scada_df.index) row_power = mi.scada_df[active_power_col].to_numpy(dtype=float) + # Only the turbines the estimate used (the test and its chosen references) count as rows. + in_wide = mi.scada_df[mi.turbine_col].isin(wide.columns).to_numpy() ts_baseline = toggle_rows.campaign_baseline - row_baseline = row_rows.campaign_baseline + row_baseline = row_rows.campaign_baseline & in_wide ts_masks = {"all": np.ones(len(wide), dtype=bool), "baseline": ts_baseline, "upgraded": toggle_rows.upgraded} - row_masks = {"all": np.ones(len(mi.scada_df), dtype=bool), "baseline": row_baseline, "upgraded": row_rows.upgraded} + row_masks = {"all": in_wide, "baseline": row_baseline, "upgraded": row_rows.upgraded & in_wide} rows = [] for segment in _SEGMENTS: @@ -766,6 +1492,8 @@ def _save_plots( mi: MethodInput, test: str, used: np.ndarray, + ref_total: npt.NDArray[np.float64], + reference_mode: str, timebase: pd.Timedelta, active_power_col: str, ) -> None: @@ -776,11 +1504,10 @@ def _save_plots( The baseline is the strict campaign off-blocks the estimate used, so the plots never disagree with the headline. """ - refs = [c for c in wide.columns if c != test] toggle_rows = resolve_toggle(mi.upgrade_timing, wide.index) baseline_mask = toggle_rows.campaign_baseline test_pw = wide[test].to_numpy(dtype=float) - ref_total = wide[refs].sum(axis=1).to_numpy(dtype=float) + ref_label = f"{reference_mode} reference {active_power_col}" upgrade_start = toggle_upgrade_start(mi.upgrade_timing, wide.index) segments = ( ("baseline", used & baseline_mask, "C0"), @@ -795,7 +1522,7 @@ def _save_plots( if np.isfinite(rho) and seg.any(): x_max = float(np.nanmax(ref_total[seg])) ax.plot([0, x_max], [0, rho * x_max], color=color, linewidth=1.5) - ax.set_xlabel(f"sum of reference {active_power_col} [kW]") + ax.set_xlabel(f"{ref_label} [kW]") ax.set_ylabel(f"{active_power_col} @ {test} [kW]") ax.set_title(f"{test}: test vs reference-total power") ax.grid(visible=True, alpha=0.3) @@ -841,6 +1568,67 @@ def _save_plots( _save(fig, plots_dir / stages.FILTER / f"{test}_coverage_timeseries.png") +def _save_pairing_gap_plot( + path: Path, + *, + index: pd.DatetimeIndex, + selection: _Selection, + rows: ToggleRowSets, + timebase: pd.Timedelta, + max_gap: pd.Timedelta | None, + test: str, +) -> None: + """Histogram each used row's gap to the nearest used row of the other segment, before and after pairing. + + Gaps beyond the plotted range are piled into the last bin. + """ + step_min = timebase.total_seconds() / 60.0 + cap_steps = max(24, 3 * round(max_gap / timebase)) if max_gap is not None else 24 + edges = (np.arange(cap_steps + 2) - 0.5) * step_min + hours_per_row = timebase / pd.Timedelta(hours=1) + + fig, ax = plt.subplots(figsize=(9, 5)) + text = [] + drawn = [("before pairing" if max_gap is not None else "used rows", selection.not_excluded, "C0")] + if max_gap is not None: + drawn.append(("after pairing", selection.paired, "C1")) + for label, mask, color in drawn: + kept = mask.to_numpy() + gap_min = ( + gap_to_other_segment(index, baseline=kept & rows.campaign_baseline, upgraded=kept & rows.upgraded) / 60.0 + ) + gap_min = gap_min[~np.isnan(gap_min)] + ax.hist( + np.minimum(gap_min, cap_steps * step_min), + bins=edges, + histtype="stepfilled", + alpha=0.45, + color=color, + label=label, + ) + n_base = int((kept & rows.campaign_baseline).sum()) + n_up = int((kept & rows.upgraded).sum()) + text.append(f"{label}: {n_base * hours_per_row:.1f} h baseline, {n_up * hours_per_row:.1f} h upgraded") + if max_gap is not None: + ax.axvline( + max_gap / pd.Timedelta(minutes=1), + color="k", + linestyle=":", + label=f"pairing_max_gap {max_gap / pd.Timedelta(minutes=1):g} min", + ) + else: + text.append("no pairing filter configured") + ax.text(0.98, 0.95, "\n".join(text), transform=ax.transAxes, ha="right", va="top", fontsize=9) + ax.set_xlabel(f"gap to nearest used row of the other segment [min] (last bin: >= {cap_steps * step_min:g})") + ax.set_ylabel("used rows") + ax.set_yscale("log") + ax.set_title(f"{test}: gap from each used row to the opposite toggle state") + ax.grid(visible=True, alpha=0.3) + ax.legend(loc="upper right", bbox_to_anchor=(0.99, 0.83)) + fig.tight_layout() + _save(fig, path) + + def _save_per_bin_plot(path: Path, *, per_bin: pd.DataFrame, test: str, active_power_col: str) -> None: """Plot the per-power-bin uplift with each bin's used-record count underneath. diff --git a/benchmarking/campaigns/methods.py b/benchmarking/campaigns/methods.py index f708edca..bf88442f 100644 --- a/benchmarking/campaigns/methods.py +++ b/benchmarking/campaigns/methods.py @@ -8,7 +8,7 @@ from benchmarking.baselines.power_model import CURATED_ERA5_EXCLUDE, TUNED_MODEL_PARAMS, PowerModelMethod from benchmarking.baselines.toggle_specialist import ToggleSpecialistMethod from benchmarking.diagnostics.context import ERA5_UNLOCATED -from benchmarking.synthetic import HOT_COLUMNS +from benchmarking.synthetic import HOT_COLUMNS, ToggleSchedule if TYPE_CHECKING: from pathlib import Path @@ -52,9 +52,13 @@ def carried_forward_methods( """ methods: list[Method] = [NaiveRatioMethod(columns=HOT_COLUMNS, out_dir=out_dir / "naive_ratio", save_plots=True)] if spec.mode == "toggle": + assert isinstance(spec.upgrade_timing, ToggleSchedule) # noqa: S101 - what ``mode == "toggle"`` means methods.append( ToggleSpecialistMethod( columns=HOT_COLUMNS, + toggle_period=spec.upgrade_timing.period, + # a declared campaign's schedule tiles its cycles from its start + toggle_datum=spec.upgrade_timing.start, out_dir=out_dir / "toggle_specialist", save_plots=True, conditions=("power",), diff --git a/benchmarking/harness/method.py b/benchmarking/harness/method.py index 1c17ace4..a70897a2 100644 --- a/benchmarking/harness/method.py +++ b/benchmarking/harness/method.py @@ -120,6 +120,10 @@ class MethodOutput: reference the method ruled out, so it stays visible without dragging that headline. :param screen_passes: optional per-pass frame from a method's reference-validity screen, so an analyst can see, and disagree with, which references it ruled out and why. + :param selection_accounting: optional per-(``stage``, ``segment``) frame of rows kept after each + of a method's row-selection stages (``n_kept``, ``kept_fraction``, ``stage_kept_fraction``), + so a consumer can see how much data each stage removed and whether it treated the segments + evenly. The harness never interprets it. """ p50_overall: float @@ -129,6 +133,7 @@ class MethodOutput: labeled_rows: pd.DataFrame | None = None reference_uplifts: pd.DataFrame | None = None screen_passes: pd.DataFrame | None = None + selection_accounting: pd.DataFrame | None = None @runtime_checkable diff --git a/benchmarking/harness/replicates.py b/benchmarking/harness/replicates.py index 3b41570d..b3ece1e8 100644 --- a/benchmarking/harness/replicates.py +++ b/benchmarking/harness/replicates.py @@ -16,6 +16,7 @@ from typing import TYPE_CHECKING, Literal import numpy as np +import pandas as pd from benchmarking.harness.campaign import resolve_campaign_grid from benchmarking.harness.northing import north_scada @@ -25,7 +26,6 @@ from collections.abc import Iterator import numpy.typing as npt - import pandas as pd from benchmarking.harness.campaign import CampaignUnit from benchmarking.synthetic import ColumnSchema, SyntheticDataset, UpliftResult @@ -65,6 +65,8 @@ class StudyConfig: :param campaign_months: the campaign-length sweep grid, in months (exclusive with ``campaign_weeks``) :param campaign_weeks: the campaign-length sweep grid, in weeks (exclusive with ``campaign_months``) :param toggle_period: toggle on/off cycle length (toggle mode only) + :param toggle_datum: a boundary of the toggle-cycle grid (toggle mode only); treatment starts are + drawn only from cycle boundaries, as a clock-aligned controller toggles :param seed: RNG seed for the draws """ @@ -76,11 +78,15 @@ class StudyConfig: campaign_months: list[int] | None = None campaign_weeks: list[int] | None = None toggle_period: pd.Timedelta | None = None + toggle_datum: pd.Timestamp | None = None seed: int = 0 def __post_init__(self) -> None: - """Validate that exactly one campaign-length grid is set.""" + """Validate that exactly one campaign-length grid is set, and a toggle study's cycle grid.""" resolve_campaign_grid(campaign_months=self.campaign_months, campaign_weeks=self.campaign_weeks) + if self.mode == "toggle" and (self.toggle_period is None or self.toggle_datum is None): + msg = "toggle_period and toggle_datum are required for mode='toggle'" + raise ValueError(msg) @property def campaign_lengths(self) -> list[int]: @@ -157,7 +163,11 @@ def iter_replicates( :param rated_power_kw: turbine rating, passed to the generator and to the northing step """ subset = base_scada[base_scada[columns.turbine].isin(study.turbine_subset)] - candidates = _candidate_starts(subset.index, study.treatment_start_range) + candidates = _candidate_starts( + subset.index, + study.treatment_start_range, + cycle_grid=(study.toggle_datum, study.toggle_period) if study.mode == "toggle" else None, + ) rng = np.random.default_rng(study.seed) turbines = np.asarray(study.turbine_subset) @@ -230,22 +240,29 @@ def build_replicates( def _candidate_starts( - index: pd.DatetimeIndex, treatment_start_range: tuple[pd.Timestamp, pd.Timestamp] + index: pd.DatetimeIndex, + treatment_start_range: tuple[pd.Timestamp, pd.Timestamp], + *, + cycle_grid: tuple[pd.Timestamp | None, pd.Timedelta | None] | None = None, ) -> npt.NDArray[np.datetime64]: - """Return unique on-grid timestamps within the draw range (so draws land on real records).""" + """Return unique on-grid timestamps within the draw range (so draws land on real records). + + With ``cycle_grid`` ``(datum, period)``, only the records on a toggle-cycle boundary. + """ lo, hi = treatment_start_range unique = index.unique() candidates = unique[(unique >= lo) & (unique <= hi)] + if cycle_grid is not None: + datum, period = cycle_grid + candidates = candidates[(candidates - datum) % period == pd.Timedelta(0)] if len(candidates) == 0: - msg = f"no records in treatment_start_range {treatment_start_range}" + msg = f"no records in treatment_start_range {treatment_start_range} (on the toggle-cycle grid, if any)" raise ValueError(msg) return candidates.sort_values().to_numpy() def _upgrade_timing(study: StudyConfig, treatment_start: pd.Timestamp) -> pd.Timestamp | ToggleSchedule: if study.mode == "toggle": - if study.toggle_period is None: - msg = "toggle_period is required for mode='toggle'" - raise ValueError(msg) + assert study.toggle_period is not None # noqa: S101 - guaranteed by StudyConfig.__post_init__ return ToggleSchedule(period=study.toggle_period, start=treatment_start) return treatment_start diff --git a/benchmarking/synthetic/schema.py b/benchmarking/synthetic/schema.py index 075280cb..b4c0ad26 100644 --- a/benchmarking/synthetic/schema.py +++ b/benchmarking/synthetic/schema.py @@ -50,9 +50,15 @@ class ColumnSchema: :param ambient_temp: ambient temperature :param exclude_row: names a caller-supplied **boolean** column marking rows to drop from row selection (e.g. special operating modes the treatment cannot affect). Optional and off by - default; when set, a method that honours it excludes the *test* turbine's flagged rows - alongside the downtime filter. Must be all-``False`` where unknown (never NaN) so an expanded - time index never turns a gap into an exclusion; a method honouring the role raises on NaN. + default; when set, a method that honours it excludes the flagged timestamps alongside the + downtime filter. Which turbines' flags count is the method's: the toggle specialist counts + the test turbine's and its references' (a flagged reference row counts the way a reference + outage does). Must be all-``False`` where unknown (never NaN) so an expanded time index + never turns a gap into an exclusion; a method honouring the role raises on NaN. + :param exclude_block: names a caller-supplied **boolean** column like ``exclude_row``, but a + flagged row removes the whole fixed-length block (one toggle cycle) around it, for every + turbine, rather than its own timestamp alone. For a selection that must stay symmetric + between toggle states (both lose the same conditions). Same never-NaN contract. """ turbine: str @@ -67,6 +73,7 @@ class ColumnSchema: nacelle_position: str | None = None ambient_temp: str | None = None exclude_row: str | None = None + exclude_block: str | None = None def northed(self, role: str) -> str: """Return the derived north-calibrated column name for ``role``, e.g. ``northed_YawAngleMean``. diff --git a/tests/benchmarking/baselines/test_block_bootstrap.py b/tests/benchmarking/baselines/test_block_bootstrap.py index 2e30c95a..61dc8269 100644 --- a/tests/benchmarking/baselines/test_block_bootstrap.py +++ b/tests/benchmarking/baselines/test_block_bootstrap.py @@ -12,7 +12,7 @@ import pandas as pd import pytest -from benchmarking.baselines.block_bootstrap import bootstrap_ratio_uplift, relative_scatter +from benchmarking.baselines.block_bootstrap import bootstrap_ratio_uplift, combine_estimates, relative_scatter _TIMEBASE = pd.Timedelta(minutes=10) @@ -207,6 +207,19 @@ def test_no_records_gives_nan_cells(self) -> None: assert result.n_blocks == 0 assert np.isnan(result.cells["overall"].sigma) + def test_a_resample_of_almost_no_reference_energy_is_not_a_spread_sample(self) -> None: + """A calm spell drawn over and over divides by nearly nothing; those resamples are dropped, not believed.""" + case = _case(n=4032) + calm = np.arange(len(case["times"])) < int(0.9 * len(case["times"])) + case["ref_total"] = np.where(calm, 1e-3, case["ref_total"]) + case["test_power"] = np.where( + calm, 1e-3 * (1.0 + 0.5 * np.random.default_rng(1).standard_normal(len(calm))), case["test_power"] + ) + cell = _run(case).cells["overall"] + assert np.isfinite(cell.sigma) + assert cell.sigma < 1.0 + assert cell.frac_resamples_finite < 1.0 + def test_too_few_resamples_for_a_spread_gives_nan(self) -> None: assert np.isnan(_run(_case(), n_resamples=1).cells["overall"].sigma) @@ -380,3 +393,125 @@ def test_sigma_scales_down_smoothly_as_noise_falls(self) -> None: sigmas.append(_run(case).cells["overall"].sigma) assert sigmas == sorted(sigmas, reverse=True) assert sigmas[-1] < sigmas[0] / 10, "no floor is arresting the descent" + + +class TestResamplesAreExposed: + def test_each_cell_carries_its_resampled_uplifts(self) -> None: + result = _run(_case(), n_resamples=300) + samples = result.resamples["overall"] + assert samples.shape == (300,) + assert np.std(samples[np.isfinite(samples)], ddof=1) == pytest.approx(result.cells["overall"].sigma_bootstrap) + + def test_no_records_means_no_resamples(self) -> None: + empty = { + k: np.array([], dtype=bool if k in ("upgraded", "baseline") else float) + for k in ("test_power", "ref_total", "upgraded", "baseline") + } + result = bootstrap_ratio_uplift( + times=_timeline(0), + cell_membership={"overall": np.array([], dtype=bool)}, + campaign_start=pd.Timestamp("2020-01-01", tz="UTC"), + campaign_end=pd.Timestamp("2020-01-08", tz="UTC"), + timebase=_TIMEBASE, + block_hours=48.0, + n_resamples=100, + seed=0, + **empty, + ) + assert result.resamples == {} + + +class TestCombineEstimates: + """Inverse-variance blend of two estimates of the same quantity, honouring their resample correlation.""" + + @staticmethod + def _samples(rho: float, *, n: int = 2000, seed: int = 0) -> tuple[np.ndarray, np.ndarray]: + rng = np.random.default_rng(seed) + a = rng.standard_normal(n) + b = rho * a + np.sqrt(1.0 - rho**2) * rng.standard_normal(n) + return a, b + + def test_independent_components_give_the_naive_inverse_variance_sigma(self) -> None: + a, b = self._samples(0.0) + out = combine_estimates((0.02, 1.0, a), (0.04, 2.0, b)) + # The sampled correlation is ~0, not exactly 0, so the weights land near the naive 0.8 / 0.2. + assert out.weight_a == pytest.approx(0.8, abs=0.02) + assert out.estimate == pytest.approx(0.8 * 0.02 + 0.2 * 0.04, abs=1e-3) + assert out.sigma == pytest.approx(1.0 / np.sqrt(1.0 + 0.25), abs=0.03) + + def test_a_correlation_widens_the_combined_sigma(self) -> None: + a, b = self._samples(0.8) + out = combine_estimates((0.0, 1.0, a), (0.0, 1.0, b)) + assert out.correlation == pytest.approx(0.8, abs=0.03) + assert out.sigma == pytest.approx(np.sqrt((1 + 0.8) / 2), abs=0.03) + + def test_a_correlated_noisier_component_drops_out_rather_than_widening_the_blend(self) -> None: + # Minimum-variance weights for correlated estimates: with r >= sb/sa the noisier leg gets no weight. + a, b = self._samples(0.8) + out = combine_estimates((0.10, 2.0, 2.0 * a), (0.02, 1.0, b)) + assert out.weight_a == 0.0 + assert out.estimate == pytest.approx(0.02) + assert out.sigma == pytest.approx(1.0) + + def test_the_blend_is_never_worse_than_its_better_component(self) -> None: + for rho in (0.0, 0.3, 0.6, 0.9): + a, b = self._samples(rho, seed=1) + out = combine_estimates((0.0, 1.5, 1.5 * a), (0.0, 1.0, b)) + assert out.sigma <= 1.0 + 1e-9, rho + + def test_one_unusable_component_leaves_the_other_untouched(self) -> None: + a, b = self._samples(0.0) + out = combine_estimates((float("nan"), float("nan"), a), (0.03, 2.0, b)) + assert out.weight_a == 0.0 + assert out.estimate == 0.03 + assert out.sigma == 2.0 + + def test_neither_usable_is_nan(self) -> None: + out = combine_estimates((float("nan"), float("nan"), np.array([])), (float("nan"), float("nan"), np.array([]))) + assert np.isnan(out.estimate) + assert np.isnan(out.sigma) + + def test_missing_resamples_assume_a_half_correlation(self) -> None: + out = combine_estimates((0.0, 1.0, np.array([])), (0.0, 1.0, np.array([]))) + assert out.correlation == 0.5 + assert out.sigma == pytest.approx(np.sqrt((1 + 0.5) / 2)) + + def test_resamples_without_spread_assume_a_half_correlation(self) -> None: + # identical turbines: every resample of a leg reads exactly the same uplift + out = combine_estimates((0.0, 1.0, np.zeros(100)), (0.0, 1.0, np.linspace(-1, 1, 100))) + assert out.correlation == 0.5 + + +class TestASparseCellWithNoFallbackReportsNaN: + """Two records in a long campaign: every resample that sees them is identical, so the bootstrap's + std is float residue (~1e-17), and the per-record scatter cannot be measured either. That is no + uncertainty, not a tiny one.""" + + def test_two_records_report_nan_rather_than_float_residue(self) -> None: + times = _timeline(2016) + keep = np.zeros(2016, dtype=bool) + keep[[1000, 1001]] = True + on = np.zeros(2016, dtype=bool) + on[1000] = True + case = { + "times": times[keep], + "test_power": np.array([800.0, 750.0]), + "ref_total": np.array([1000.0, 1000.0]), + "upgraded": on[keep], + "baseline": ~on[keep], + "cell_membership": {"overall": np.ones(2, dtype=bool)}, + "campaign_start": times[0], + "campaign_end": times[-1], + "timebase": _TIMEBASE, + } + cell = _run(case).cells["overall"] + assert np.isnan(cell.sigma_fallback) + assert np.isnan(cell.sigma) + + +class TestCombineEstimatesRefusesUnpairedResamples: + """The correlation only means something when the two resample arrays came from the same block draws.""" + + def test_different_resample_counts_raise(self) -> None: + with pytest.raises(ValueError, match="resample"): + combine_estimates((0.0, 1.0, np.zeros(100)), (0.0, 1.0, np.zeros(99))) diff --git a/tests/benchmarking/baselines/test_toggle_end_to_end.py b/tests/benchmarking/baselines/test_toggle_end_to_end.py index 7fbe0043..b87f14b8 100644 --- a/tests/benchmarking/baselines/test_toggle_end_to_end.py +++ b/tests/benchmarking/baselines/test_toggle_end_to_end.py @@ -36,6 +36,7 @@ def test_naive_recovers_toggle_uplift(tmp_path) -> None: # noqa: ANN001 min_pre_months=24, campaign_months=[6], toggle_period=pd.Timedelta(minutes=40), + toggle_datum=pd.Timestamp("1970-01-01", tz="UTC"), n_replicates=1, seed=0, ) diff --git a/tests/benchmarking/baselines/test_toggle_specialist.py b/tests/benchmarking/baselines/test_toggle_specialist.py index a13417b6..9379b395 100644 --- a/tests/benchmarking/baselines/test_toggle_specialist.py +++ b/tests/benchmarking/baselines/test_toggle_specialist.py @@ -24,6 +24,8 @@ _daily_segment_ratio, _expected_per_day, _infer_timebase, + gap_to_other_segment, + pair_within, restrict_to_campaign, ) from benchmarking.harness.conditions import condition_bins @@ -41,6 +43,10 @@ # Larger than any test timebase's full period, so the (required) availability filter keeps # every row unless a test deliberately sets a lower value. _FULLY_AVAILABLE_SECS = 3600.0 +# One toggle cycle of the 20-minute ``ToggleSchedule`` most tests use (``toggle_period`` is required). +_CYCLE = pd.Timedelta(minutes=20) +# Every fixture's index starts here, on the cycle grid of every period the tests use. +_DATUM = pd.Timestamp("2020-01-01", tz="UTC") # The method reads active_power + availability from the schema; the other required roles are # unused by the ratio, so name them with placeholders. _COLUMNS = ColumnSchema( @@ -117,7 +123,7 @@ def test_toggle_specialist_shares_no_wind_up_code() -> None: class TestRecovery: def test_toggle_recovers_known_uplift(self) -> None: scada, schedule, _ = _toggle_case(uplift=0.03) - out = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert isinstance(out, MethodOutput) @@ -169,9 +175,9 @@ def _varying_rho_case( def _per_bin(scada: pd.DataFrame, schedule: ToggleSchedule) -> pd.DataFrame: """Run the method with power conditioning on and return its populated per-bin rows.""" - out = ToggleSpecialistMethod(columns=_COLUMNS, conditions=("power",), rated_power_kw=_RATED_KW).estimate( - MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) - ) + out = ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",), rated_power_kw=_RATED_KW + ).estimate(MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) assert out.p50_by_condition is not None frame = out.p50_by_condition return frame[frame["n_records"] > 0] @@ -236,7 +242,7 @@ def test_uplift_confined_to_high_power_shows_only_there(self) -> None: scada = _scada({"T1": test, "R1": ref1, "R2": ref2}, idx) per_bin = _per_bin(scada, schedule) - low_bins = per_bin[per_bin["sum_counterfactual"] > 0].iloc[:2] + low_bins = per_bin[per_bin["sum_counterfactual_sum"] > 0].iloc[:2] assert low_bins["p50_uplift"].abs().max() < 0.005 assert per_bin["p50_uplift"].max() == pytest.approx(0.10, abs=0.01) @@ -245,7 +251,7 @@ class TestConditionsConfiguration: def test_default_reports_no_conditions(self) -> None: # back-compat: existing callers get exactly today's behaviour and need no rating scada, schedule, _ = _toggle_case(uplift=0.03) - out = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert out.p50_by_condition is None @@ -254,33 +260,73 @@ def test_power_conditioning_does_not_move_the_headline(self) -> None: # the per-bin decomposition is additional reporting, never a change to the estimate scada, schedule, _ = _toggle_case(uplift=0.03) mi = MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) - plain = ToggleSpecialistMethod(columns=_COLUMNS).estimate(mi) + plain = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate(mi) conditioned = ToggleSpecialistMethod( - columns=_COLUMNS, conditions=("power",), rated_power_kw=_RATED_KW + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",), rated_power_kw=_RATED_KW ).estimate(mi) assert conditioned.p50_overall == plain.p50_overall def test_frame_is_labelled_with_the_power_condition(self) -> None: scada, schedule = _varying_rho_case(uplift=0.05) - out = ToggleSpecialistMethod(columns=_COLUMNS, conditions=("power",), rated_power_kw=_RATED_KW).estimate( - MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) - ) + out = ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",), rated_power_kw=_RATED_KW + ).estimate(MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) assert out.p50_by_condition is not None assert set(out.p50_by_condition["condition"]) == {"power"} - for col in ("condition_bin", "p50_uplift", "n_records", "sum_actual", "sum_counterfactual"): + for col in ("condition_bin", "p50_uplift", "n_records", "sigma_uplift", "weight_sum"): assert col in out.p50_by_condition.columns + for leg in ("sum", "block_mean"): + assert {f"sum_actual_{leg}", f"sum_counterfactual_{leg}"} <= set(out.p50_by_condition.columns) + assert "sum_actual" not in out.p50_by_condition.columns + + def test_each_bin_is_the_blend_of_its_legs_at_its_weight(self) -> None: + """Each leg's bin uplift is its own energy ratio; the reported bin is their blend at the bin's weight.""" + scada, schedule = _varying_rho_case(uplift=0.05, noise_frac=0.05) + out = ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",), rated_power_kw=_RATED_KW + ).estimate(MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) + assert out.p50_by_condition is not None + frame = out.p50_by_condition.dropna(subset=["p50_uplift", "weight_sum"]) + assert len(frame) + legs = { + leg: frame[f"sum_actual_{leg}"] / frame[f"sum_counterfactual_{leg}"] - 1.0 for leg in ("sum", "block_mean") + } + weight = frame["weight_sum"] + # a leg with no estimate in a bin carries no weight there + blended = np.where(weight > 0, weight * legs["sum"], 0.0) + np.where( + weight < 1, (1.0 - weight) * legs["block_mean"], 0.0 + ) + np.testing.assert_allclose(frame["p50_uplift"], blended) + + def test_a_single_leg_frame_carries_that_legs_sums(self) -> None: + scada, schedule = _block_case(uplift=0.04) + out = _estimate(scada, schedule, toggle_period=_BLOCK, conditions=("power",), rated_power_kw=_RATED_KW) + assert out.p50_by_condition is not None + assert {"sum_actual_block_mean", "sum_counterfactual_block_mean"} <= set(out.p50_by_condition.columns) def test_ws_condition_raises_citing_the_method_limit(self) -> None: with pytest.raises(ValueError, match="does not support"): - ToggleSpecialistMethod(columns=_COLUMNS, conditions=("ws",), rated_power_kw=_RATED_KW) + ToggleSpecialistMethod( + columns=_COLUMNS, + toggle_period=_CYCLE, + toggle_datum=_DATUM, + conditions=("ws",), + rated_power_kw=_RATED_KW, + ) def test_unknown_condition_raises(self) -> None: with pytest.raises(ValueError, match="unknown condition"): - ToggleSpecialistMethod(columns=_COLUMNS, conditions=("bogus",), rated_power_kw=_RATED_KW) + ToggleSpecialistMethod( + columns=_COLUMNS, + toggle_period=_CYCLE, + toggle_datum=_DATUM, + conditions=("bogus",), + rated_power_kw=_RATED_KW, + ) def test_power_without_a_rating_raises(self) -> None: with pytest.raises(ValueError, match="rated_power_kw"): - ToggleSpecialistMethod(columns=_COLUMNS, conditions=("power",)) + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",)) class TestPerBinSparseData: @@ -289,9 +335,9 @@ def test_bins_with_no_data_are_nan_not_imputed(self) -> None: # do about it, and an imputed prior would manufacture false confidence. Rating the turbine far # above the data's power range guarantees empty upper bins rather than hoping for them. scada, schedule = _varying_rho_case(uplift=0.05) - out = ToggleSpecialistMethod(columns=_COLUMNS, conditions=("power",), rated_power_kw=4000.0).estimate( - MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) - ) + out = ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",), rated_power_kw=4000.0 + ).estimate(MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) assert out.p50_by_condition is not None empty = out.p50_by_condition[out.p50_by_condition["n_records"] == 0] assert not empty.empty @@ -299,9 +345,9 @@ def test_bins_with_no_data_are_nan_not_imputed(self) -> None: def test_every_bin_is_represented(self) -> None: scada, schedule = _varying_rho_case(uplift=0.05) - out = ToggleSpecialistMethod(columns=_COLUMNS, conditions=("power",), rated_power_kw=_RATED_KW).estimate( - MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) - ) + out = ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",), rated_power_kw=_RATED_KW + ).estimate(MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) assert out.p50_by_condition is not None assert len(out.p50_by_condition) == len(condition_bins("power", rated_power_kw=_RATED_KW)) - 1 @@ -310,13 +356,18 @@ class TestPerBinDiagnostics: def test_per_bin_csv_is_written(self, tmp_path: Path) -> None: scada, schedule = _varying_rho_case(uplift=0.05) ToggleSpecialistMethod( - columns=_COLUMNS, conditions=("power",), rated_power_kw=_RATED_KW, out_dir=tmp_path + columns=_COLUMNS, + toggle_period=_CYCLE, + toggle_datum=_DATUM, + conditions=("power",), + rated_power_kw=_RATED_KW, + out_dir=tmp_path, ).estimate(MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) assert list(tmp_path.rglob("*_by_power_bin_*.csv")) def test_no_per_bin_csv_when_conditioning_is_off(self, tmp_path: Path) -> None: scada, schedule, _ = _toggle_case(uplift=0.03) - ToggleSpecialistMethod(columns=_COLUMNS, out_dir=tmp_path).estimate( + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, out_dir=tmp_path).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert not list(tmp_path.rglob("*_by_power_bin_*.csv")) @@ -325,6 +376,8 @@ def test_per_bin_plot_written_with_save_plots(self, tmp_path: Path) -> None: scada, schedule = _varying_rho_case(uplift=0.05) ToggleSpecialistMethod( columns=_COLUMNS, + toggle_period=_CYCLE, + toggle_datum=_DATUM, conditions=("power",), rated_power_kw=_RATED_KW, out_dir=tmp_path, @@ -342,7 +395,7 @@ def test_prepost_timestamp_raises(self) -> None: treated = np.asarray(idx >= upgrade) scada = _recovery_scada(idx, treated=treated, uplift=0.05) with pytest.raises(ValueError, match="toggle"): - ToggleSpecialistMethod(columns=_COLUMNS).estimate( + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=upgrade, turbine_col=_TURBINE_COL) ) @@ -359,13 +412,15 @@ def test_blank_availability_role_raises(self, blank: str) -> None: # A schema that leaves the availability role blank (empty or whitespace) would silently skip # downtime filtering, so construction must reject it. with pytest.raises(ValueError, match="availability"): - ToggleSpecialistMethod(columns=replace(_COLUMNS, availability=blank)) + ToggleSpecialistMethod( + columns=replace(_COLUMNS, availability=blank), toggle_period=_CYCLE, toggle_datum=_DATUM + ) def test_missing_availability_column_raises(self) -> None: scada, schedule, _ = _toggle_case() scada = scada.drop(columns=[_AVAIL_COL]) with pytest.raises(ValueError, match="availability"): - ToggleSpecialistMethod(columns=_COLUMNS).estimate( + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) @@ -378,7 +433,7 @@ def test_down_reference_timestamps_are_excluded(self) -> None: down = (corrupted[_TURBINE_COL] == "R1") & corrupted.index.isin([idx[3], idx[4]]) corrupted.loc[down, _POWER_COL] = 1e6 corrupted.loc[down, _AVAIL_COL] = 0.0 - out = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=corrupted, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert out.p50_overall == pytest.approx(0.05) @@ -388,7 +443,7 @@ class TestCompleteCase: def test_drops_timestamp_when_a_reference_is_nan(self) -> None: scada, schedule, _ = _toggle_case() idx = scada.index.unique().sort_values() - clean = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + clean = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) @@ -397,7 +452,7 @@ def test_drops_timestamp_when_a_reference_is_nan(self) -> None: corrupted = scada.copy() mask = (corrupted[_TURBINE_COL] == "R1") & (corrupted.index == idx[3]) corrupted.loc[mask, _POWER_COL] = np.nan - out = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=corrupted, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert out.p50_overall == pytest.approx(clean.p50_overall) @@ -408,12 +463,12 @@ def test_nan_at_unused_timestamp_does_not_change_estimate(self) -> None: # make idx[2] already unused (test NaN), then add a second NaN at the same timestamp base = scada.copy() base.loc[(base[_TURBINE_COL] == "T1") & (base.index == idx[2]), _POWER_COL] = np.nan - before = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + before = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=base, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) after_df = base.copy() after_df.loc[(after_df[_TURBINE_COL] == "R1") & (after_df.index == idx[2]), _POWER_COL] = np.nan - after = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + after = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=after_df, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert after.p50_overall == pytest.approx(before.p50_overall) @@ -428,7 +483,7 @@ def test_estimate_invariant_to_timebase(self) -> None: treated = np.asarray(treated_mask(idx, schedule)) scada = _recovery_scada(idx, treated=treated, uplift=0.04) results.append( - ToggleSpecialistMethod(columns=_COLUMNS) + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=schedule.period, toggle_datum=_DATUM) .estimate(MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) .p50_overall ) @@ -438,6 +493,8 @@ def test_override_timebase_used_for_mwh(self, tmp_path) -> None: # noqa: ANN001 scada, schedule, _ = _toggle_case() method = ToggleSpecialistMethod( columns=_COLUMNS, + toggle_period=pd.Timedelta(hours=1), + toggle_datum=_DATUM, out_dir=tmp_path, timebase=pd.Timedelta(minutes=30), ) @@ -452,7 +509,7 @@ def test_override_timebase_used_for_mwh(self, tmp_path) -> None: # noqa: ANN001 class TestActivePowerOnly: def test_extra_columns_ignored(self) -> None: scada, schedule, _ = _toggle_case() - plain = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + plain = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) @@ -460,21 +517,23 @@ def test_extra_columns_ignored(self) -> None: rng = np.random.default_rng(0) with_extra["ws"] = rng.normal(size=len(with_extra)) with_extra["rpm"] = rng.normal(size=len(with_extra)) - out = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=with_extra, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert out.p50_overall == pytest.approx(plain.p50_overall) class TestErrors: - def test_no_reference_raises(self) -> None: + def test_no_reference_is_the_block_leg_alone(self, caplog: pytest.LogCaptureFixture) -> None: idx = _index(40) schedule = ToggleSchedule(period=pd.Timedelta(minutes=20), start=idx[0]) scada = _scada({"T1": np.ones(len(idx))}, idx) - with pytest.raises(ValueError, match="reference"): - ToggleSpecialistMethod(columns=_COLUMNS).estimate( + with caplog.at_level("WARNING"): + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) + assert _legs(out) == {"block_mean"} + assert "no reference turbines" in caplog.text def test_no_used_baseline_returns_nan(self) -> None: scada, schedule, treated = _toggle_case() @@ -483,7 +542,7 @@ def test_no_used_baseline_returns_nan(self) -> None: off_ts = idx[~treated] is_baseline_test = (scada[_TURBINE_COL] == "T1") & scada.index.isin(off_ts) scada.loc[is_baseline_test, _POWER_COL] = np.nan - out = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert np.isnan(out.p50_overall) @@ -536,7 +595,9 @@ def _read_only_csv(folder: Path, kind: str) -> pd.DataFrame: class TestDiagnostics: def _run(self, tmp_path, *, save_plots: bool = False): # noqa: ANN001, ANN202 scada, schedule, _ = _toggle_case(uplift=0.06) - method = ToggleSpecialistMethod(columns=_COLUMNS, out_dir=tmp_path, save_plots=save_plots) + method = ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, out_dir=tmp_path, save_plots=save_plots + ) out = method.estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) @@ -546,26 +607,40 @@ def test_writes_both_csvs(self, tmp_path) -> None: # noqa: ANN001 self._run(tmp_path) stats = _read_only_csv(tmp_path, "data_stats") results = _read_only_csv(tmp_path, "results") - assert sorted(stats["segment"]) == ["all", "baseline", "upgraded"] + for _, leg in stats.groupby("component"): + assert sorted(leg["segment"]) == ["all", "baseline", "upgraded"] + assert set(stats["component"]) == {"sum", "block_mean"} assert len(results) == 1 def test_uplift_rederivable_from_stats(self, tmp_path) -> None: # noqa: ANN001 + """Each leg's uplift from its own stats rows; the headline is their blend at the reported weight.""" out, _ = self._run(tmp_path) - stats = _read_only_csv(tmp_path, "data_stats").set_index("segment") - rho_base = stats.loc["baseline", "used_test_mwh"] / stats.loc["baseline", "used_ref_total_mwh"] - rho_up = stats.loc["upgraded", "used_test_mwh"] / stats.loc["upgraded", "used_ref_total_mwh"] - assert (rho_up / rho_base - 1.0) == pytest.approx(out.p50_overall) + stats = _read_only_csv(tmp_path, "data_stats").set_index(["component", "segment"]) + results = _read_only_csv(tmp_path, "results").iloc[0] + leg_uplift = {} + for leg in ("sum", "block_mean"): + rho_base = ( + stats.loc[(leg, "baseline"), "used_test_mwh"] / stats.loc[(leg, "baseline"), "used_ref_total_mwh"] + ) + rho_up = stats.loc[(leg, "upgraded"), "used_test_mwh"] / stats.loc[(leg, "upgraded"), "used_ref_total_mwh"] + leg_uplift[leg] = rho_up / rho_base - 1.0 + assert leg_uplift[leg] == pytest.approx( + results[f"ratio_upgraded_{leg}"] / results[f"ratio_baseline_{leg}"] - 1.0 + ) + weight = results["weight_sum"] + assert weight * leg_uplift["sum"] + (1.0 - weight) * leg_uplift["block_mean"] == pytest.approx(out.p50_overall) def test_results_csv_matches_estimate(self, tmp_path) -> None: # noqa: ANN001 out, _ = self._run(tmp_path) results = _read_only_csv(tmp_path, "results").iloc[0] assert results["uplift_frc"] == pytest.approx(out.p50_overall) assert results["mode"] == "toggle" - assert results["n_refs"] == 2 + assert results["n_refs_available"] == 2 + assert results["n_refs"] == len(str(results["refs_used"]).split(";")) def test_stats_coverage_full_for_complete_data(self, tmp_path) -> None: # noqa: ANN001 self._run(tmp_path) - stats = _read_only_csv(tmp_path, "data_stats").set_index("segment") + stats = _read_only_csv(tmp_path, "data_stats").query("component == 'sum'").set_index("segment") assert stats.loc["all", "rows_data_coverage"] == pytest.approx(1.0) assert stats.loc["all", "used_data_coverage"] == pytest.approx(1.0) assert stats.loc["all", "n_used_timestamps"] == 40 @@ -594,12 +669,13 @@ def test_precampaign_dropped_from_baseline(self, tmp_path) -> None: # noqa: ANN schedule = ToggleSchedule(period=pd.Timedelta(minutes=20), start=start) treated = np.asarray(treated_mask(idx, schedule)) scada = _recovery_scada(idx, treated=treated, uplift=0.05) - ToggleSpecialistMethod(columns=_COLUMNS, out_dir=tmp_path).estimate( + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, out_dir=tmp_path).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) - baseline = _read_only_csv(tmp_path, "data_stats").set_index("segment").loc["baseline"] + stats = _read_only_csv(tmp_path, "data_stats") # the baseline class starts at/after the campaign start; the pre-campaign rows are excluded. - assert pd.Timestamp(baseline["first_timestamp"]) >= start + for first in stats.loc[stats["segment"] == "baseline", "first_timestamp"]: + assert pd.Timestamp(first) >= start # --- uncertainty (non-optional, computed after the uplift) -------------------------------------- @@ -614,7 +690,25 @@ def _noisy_toggle_case(n: int = 2016, *, uplift: float = 0.03, noise_frac: float return scada, schedule +_LEG_NAMES = ("sum", "block_mean", "combined") + + +def _long_diagnostics(out: MethodOutput) -> pd.DataFrame: + """The harness frame back in per-leg rows: one row per (component, condition, condition_bin).""" + wide = out.uncertainty_diagnostics + assert wide is not None + frames = [] + for leg in _LEG_NAMES: + renamed = {c: c.removesuffix(f"_{leg}") for c in wide.columns if c.endswith(f"_{leg}")} + if renamed: + leg_rows = wide[["condition", "condition_bin", *renamed]].rename(columns=renamed) + frames.append(leg_rows.dropna(subset=list(renamed.values()), how="all").assign(component=leg)) + return pd.concat(frames, ignore_index=True) + + def _estimate(scada: pd.DataFrame, schedule: ToggleSchedule, **kwargs: object) -> MethodOutput: + kwargs.setdefault("toggle_period", _CYCLE) + kwargs.setdefault("toggle_datum", _DATUM) return ToggleSpecialistMethod(columns=_COLUMNS, **kwargs).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) @@ -637,12 +731,20 @@ def test_every_reported_bin_carries_a_sigma(self) -> None: assert len(populated) > 0 assert populated["sigma_uplift"].notna().all() + def test_the_harness_frame_has_one_row_per_cell_and_no_reserved_column(self) -> None: + scada, schedule = _noisy_toggle_case() + diag = _estimate(scada, schedule, conditions=("power",), rated_power_kw=_RATED_KW).uncertainty_diagnostics + assert diag is not None + assert not diag.duplicated(["condition", "condition_bin"]).any() + assert not set(diag.columns) & {"sigma", "estimate", "truth"} + assert {"sigma_sum", "sigma_block_mean", "sigma_combined"} <= set(diag.columns) + def test_diagnostics_cover_the_headline_and_every_bin(self) -> None: scada, schedule = _noisy_toggle_case() out = _estimate(scada, schedule, conditions=("power",), rated_power_kw=_RATED_KW) - diag = out.uncertainty_diagnostics + diag = _long_diagnostics(out) assert diag is not None - assert (diag["condition"] == "overall").sum() == 1 + assert (diag["condition"] == "overall").sum() == diag["component"].nunique() assert set(diag.columns) >= { "condition", "condition_bin", @@ -652,38 +754,43 @@ def test_diagnostics_cover_the_headline_and_every_bin(self) -> None: "sigma_robust", "frac_resamples_finite", } - overall = diag[diag["condition"] == "overall"].iloc[0] + overall = diag[(diag["condition"] == "overall") & (diag["component"] == "sum")].iloc[0] assert overall["n_upgraded_records"] > 0 assert overall["n_baseline_records"] > 0 -class TestUncertaintyDoesNotChangeUplift: - def test_block_length_moves_the_bootstrap_and_nothing_else(self) -> None: - """The session's hard constraint, at the method's own seam. +class TestUncertaintyDoesNotChangeAnUplift: + """Every candidate's uplift is fixed before any bootstrap runs. - Asserted on the *bootstrap* component, not the reported sigma: the reported value is - ``max(bootstrap, fallback)``, and the fallback does not depend on block length — so on data - where it wins at both lengths it would mask the difference this test exists to check. - """ - scada, schedule = _noisy_toggle_case() - a = _estimate(scada, schedule, conditions=("power",), rated_power_kw=_RATED_KW, block_hours=6.0) - b = _estimate(scada, schedule, conditions=("power",), rated_power_kw=_RATED_KW, block_hours=96.0) - assert a.p50_overall == b.p50_overall - pd.testing.assert_series_equal( - a.p50_by_condition["p50_uplift"], - b.p50_by_condition["p50_uplift"], # type: ignore[index] - ) - boots = [ - out.uncertainty_diagnostics.set_index("condition_bin").loc["overall", "sigma_bootstrap"] # type: ignore[union-attr] - for out in (a, b) - ] - assert boots[0] != boots[1] + The bootstraps decide which reference subset is kept and how the legs are weighted, so the + headline may move with them; no subset's uplift and not the block leg's may. + """ - def test_uplift_matches_a_run_with_the_bootstrap_reduced_to_nothing(self) -> None: + @staticmethod + def _candidates(tmp_path: Path, **kwargs: object) -> tuple[pd.Series, float, pd.DataFrame]: scada, schedule = _noisy_toggle_case() - full = _estimate(scada, schedule, n_resamples=1000) - minimal = _estimate(scada, schedule, n_resamples=2) - assert full.p50_overall == minimal.p50_overall + out = _estimate(scada, schedule, out_dir=tmp_path, **kwargs) + subsets = _read_only_csv(tmp_path, "reference_subsets").set_index("refs")["uplift_frc"] + block = _read_only_csv(tmp_path, "results").loc[0, "uplift_frc_block_mean"] + diag = _long_diagnostics(out).set_index(["component", "condition_bin"]) + return subsets, block, diag + + def test_block_length_moves_the_bootstrap_and_no_uplift(self, tmp_path: Path) -> None: + """Asserted on the *bootstrap* sigma, not the reported one: the reported value is + ``max(bootstrap, fallback)`` and the fallback does not depend on block length. + """ + subsets_a, block_a, da = self._candidates(tmp_path / "a", block_hours=6.0) + subsets_b, block_b, db = self._candidates(tmp_path / "b", block_hours=96.0) + pd.testing.assert_series_equal(subsets_a, subsets_b) + assert block_a == block_b + for leg in ("sum", "block_mean"): + assert da.loc[(leg, "overall"), "sigma_bootstrap"] != db.loc[(leg, "overall"), "sigma_bootstrap"] + + def test_no_uplift_moves_when_the_bootstrap_is_reduced_to_nothing(self, tmp_path: Path) -> None: + subsets_full, block_full, _ = self._candidates(tmp_path / "full", n_resamples=1000) + subsets_min, block_min, _ = self._candidates(tmp_path / "minimal", n_resamples=2) + pd.testing.assert_series_equal(subsets_full, subsets_min) + assert block_full == block_min class TestUncertaintyRunsOnlyWhenThereIsAnUpliftToQualify: @@ -692,7 +799,7 @@ def test_a_nan_uplift_reports_a_nan_sigma(self) -> None: idx = _index(40) # every row treated -> no off rows -> rho_base is NaN scada = _recovery_scada(idx, treated=np.ones(len(idx), dtype=bool)) - out = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput( scada_df=scada, test_wtg="T1", @@ -708,7 +815,9 @@ def test_the_bootstrap_is_not_run_when_the_uplift_is_nan(self) -> None: """Diagnostics still report the counts: they are what explains why there was no answer.""" idx = _index(40) scada = _recovery_scada(idx, treated=np.ones(len(idx), dtype=bool)) - out = ToggleSpecialistMethod(columns=_COLUMNS, conditions=("power",), rated_power_kw=_RATED_KW).estimate( + out = ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",), rated_power_kw=_RATED_KW + ).estimate( MethodInput( scada_df=scada, test_wtg="T1", @@ -716,11 +825,13 @@ def test_the_bootstrap_is_not_run_when_the_uplift_is_nan(self) -> None: turbine_col=_TURBINE_COL, ) ) - diag = out.uncertainty_diagnostics + diag = _long_diagnostics(out) assert diag is not None - assert (diag["n_blocks"] == 0).all() - assert diag["frac_resamples_finite"].isna().all() - assert (diag["n_baseline_records"] == 0).all() + legs = diag[diag["component"] != "combined"] + assert (legs["n_blocks"] == 0).all() + # no resample ran, so no leg reports a finite fraction (the harness frame drops the empty column) + assert "frac_resamples_finite" not in legs or legs["frac_resamples_finite"].isna().all() + assert (legs["n_baseline_records"] == 0).all() class TestUncertaintyReproducibility: @@ -752,7 +863,8 @@ def test_an_uncertainty_csv_is_written(self, tmp_path: Path) -> None: assert len(written) == 1 frame = pd.read_csv(written[0]) assert "n_upgraded_records" in frame.columns - assert len(frame) == 7 # overall + six power bins + # overall + six power bins, per leg and for the blend + assert frame.groupby("component").size().to_dict() == {"sum": 7, "block_mean": 7, "combined": 7} # --- labeled_rows: the row selection the estimate actually used --------------------------------- @@ -792,7 +904,9 @@ def test_one_row_per_test_turbine_record(self) -> None: def test_used_reproduces_the_methods_own_used_mask(self) -> None: scada, schedule = _noisy_toggle_case() - method = ToggleSpecialistMethod(columns=_COLUMNS, conditions=("power",), rated_power_kw=_RATED_KW) + method = ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, conditions=("power",), rated_power_kw=_RATED_KW + ) mi = MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) out = method.estimate(mi) assert out.labeled_rows is not None @@ -801,13 +915,14 @@ def test_used_reproduces_the_methods_own_used_mask(self) -> None: wide = toggle_specialist._wide_column( # noqa: SLF001 mi_restricted.scada_df, turbine_col=_TURBINE_COL, value_col=_POWER_COL ) - expected = method._used_mask( # noqa: SLF001 + expected = method._selection( # noqa: SLF001 mi_restricted, wide=wide, test="T1", refs=[c for c in wide.columns if c != "T1"], timebase=pd.Timedelta(minutes=10), - ) + rows=resolve_toggle(mi_restricted.upgrade_timing, wide.index), + ).paired assert out.labeled_rows["used"].to_numpy().tolist() == expected.to_numpy().tolist() def test_a_downtime_row_is_labelled_unused(self) -> None: @@ -887,20 +1002,20 @@ def _swap_states(toggle_df: pd.DataFrame) -> pd.DataFrame: class TestReversalSymmetry: - """Which of two toggle states is *called* the baseline is a naming choice, not a measurement. + """With ``bin_power="both_states"``, which of two toggle states is *called* the baseline is a naming choice. - The bin label must therefore not depend on it. It used to: the label was the baseline state's - predicted power, so renaming the states shifted every label by the uplift and migrated rows - across bin edges -- moving per-bin estimates by an appreciable fraction of their own sigma for - no physical reason. + The bin label must then not depend on it. Under ``"baseline"`` it does by design (the label is + the baseline state's power): renaming the states shifts every label by the uplift and migrates + rows across bin edges, which is why a campaign with no true baseline state uses ``both_states``. """ def _pair(self) -> tuple[MethodOutput, MethodOutput]: scada, schedule = _varying_rho_case(n=600, uplift=0.05) index = pd.DatetimeIndex(pd.unique(scada.index)).sort_values() toggle_df = build_toggle_df(index, schedule) - forward = _estimate(scada, toggle_df, conditions=("power",), rated_power_kw=_RATED_KW) - reversed_ = _estimate(scada, _swap_states(toggle_df), conditions=("power",), rated_power_kw=_RATED_KW) + kwargs = {"conditions": ("power",), "rated_power_kw": _RATED_KW, "bin_power": "both_states"} + forward = _estimate(scada, toggle_df, **kwargs) + reversed_ = _estimate(scada, _swap_states(toggle_df), **kwargs) return forward, reversed_ def test_every_row_keeps_its_power_bin_when_the_states_are_swapped(self) -> None: @@ -953,7 +1068,7 @@ class TestExcludeRow: def test_flagged_test_rows_are_not_used(self) -> None: scada, schedule = _noisy_toggle_case() flagged = _flag(scada, turbine="T1", every=5) - out = ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=flagged, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert out.labeled_rows is not None @@ -969,42 +1084,52 @@ def test_excluding_corrupted_rows_restores_the_known_uplift(self) -> None: spoiled = corrupted[_EXCLUDE_COL].to_numpy(dtype=bool) & (corrupted[_TURBINE_COL] == "T1").to_numpy() corrupted.loc[spoiled, _POWER_COL] *= 0.5 - naive = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + naive = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=corrupted, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) - filtered = ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS).estimate( + filtered = ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=corrupted, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert filtered.p50_overall == pytest.approx(0.05, abs=1e-9) assert abs(naive.p50_overall - 0.05) > abs(filtered.p50_overall - 0.05) - def test_flagged_reference_rows_are_kept(self) -> None: - """References are never excluded here: their special modes still carry ratio information.""" + def test_flagged_reference_rows_leave_the_reference_leg_like_an_outage(self, tmp_path: Path) -> None: + """A flagged reference row (curtailed, iced, parked) removes the timestamp from the reference leg. + + The block leg uses no references, so it keeps the timestamp and the headline ``used`` stays. + """ scada, schedule = _noisy_toggle_case() + # flag both references, so the flags reach whichever subset the reference leg keeps flagged = _flag(scada, turbine="R1", every=3) - baseline = ToggleSpecialistMethod(columns=_COLUMNS).estimate( - MethodInput(scada_df=flagged, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) + flagged[_EXCLUDE_COL] = ( + flagged[_EXCLUDE_COL].to_numpy() | _flag(scada, turbine="R2", every=3)[_EXCLUDE_COL].to_numpy() ) - with_role = ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS).estimate( - MethodInput(scada_df=flagged, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) - ) - assert with_role.p50_overall == pytest.approx(baseline.p50_overall) - assert with_role.labeled_rows is not None - assert with_role.labeled_rows["used"].to_numpy().any() + out = ToggleSpecialistMethod( + columns=_EXCLUDE_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, out_dir=tmp_path + ).estimate(MethodInput(scada_df=flagged, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) + assert out.labeled_rows is not None + is_ref = (flagged[_TURBINE_COL] != "T1").to_numpy() + flagged_ts = flagged.index[is_ref & flagged[_EXCLUDE_COL].to_numpy(dtype=bool)].unique() + used = out.labeled_rows["used"].astype(bool) + assert used.loc[flagged_ts].all() + kept = _read_only_csv(tmp_path, "selection").set_index(["component", "stage", "segment"])["n_kept"] + for segment in ("baseline", "upgraded"): + assert kept[("sum", "exclude_row", segment)] < kept[("sum", "filters", segment)] + assert kept[("block_mean", "exclude_row", segment)] == kept[("block_mean", "filters", segment)] def test_nan_in_the_exclude_column_is_rejected(self) -> None: """NaN is not "excluded": the schema contract says the column is never NaN, so say so loudly.""" scada, schedule = _noisy_toggle_case() flagged = _flag(scada, turbine="T1", every=7, value=np.nan) with pytest.raises(ValueError, match=_EXCLUDE_COL): - ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS).estimate( + ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=flagged, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) def test_unset_role_excludes_nothing(self) -> None: scada, schedule = _noisy_toggle_case() flagged = _flag(scada, turbine="T1", every=5) - out = ToggleSpecialistMethod(columns=_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=flagged, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert out.labeled_rows is not None @@ -1014,7 +1139,7 @@ def test_unset_role_excludes_nothing(self) -> None: def test_absent_column_excludes_nothing(self) -> None: """The role names a column the frame does not carry: skip, do not raise.""" scada, schedule = _noisy_toggle_case() - out = ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS).estimate( + out = ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM).estimate( MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) ) assert out.labeled_rows is not None @@ -1029,9 +1154,9 @@ def test_exclusions_reach_the_shared_diagnostics(self, tmp_path: Path) -> None: """ scada, schedule = _noisy_toggle_case() flagged = _flag(scada, turbine="T1", every=5) - ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS, out_dir=tmp_path, save_plots=True).estimate( - MethodInput(scada_df=flagged, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) - ) + ToggleSpecialistMethod( + columns=_EXCLUDE_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, out_dir=tmp_path, save_plots=True + ).estimate(MethodInput(scada_df=flagged, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) run_dir = next(p for p in Path(tmp_path).iterdir() if p.is_dir()) names = {p.name for p in (run_dir / "plots").rglob("*.png")} assert "excluded_row_fraction.png" in names @@ -1039,9 +1164,9 @@ def test_exclusions_reach_the_shared_diagnostics(self, tmp_path: Path) -> None: def test_no_exclusion_plots_when_nothing_is_excluded(self, tmp_path: Path) -> None: scada, schedule = _noisy_toggle_case() - ToggleSpecialistMethod(columns=_COLUMNS, out_dir=tmp_path, save_plots=True).estimate( - MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) - ) + ToggleSpecialistMethod( + columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, out_dir=tmp_path, save_plots=True + ).estimate(MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL)) run_dir = next(p for p in Path(tmp_path).iterdir() if p.is_dir()) names = {p.name for p in (run_dir / "plots").rglob("*.png")} assert "ops_curves_excluded.png" not in names @@ -1075,7 +1200,7 @@ def _fixture() -> tuple[pd.DataFrame, ToggleSchedule]: @staticmethod def _estimate(scada: pd.DataFrame, **kwargs: object) -> float: - method = ToggleSpecialistMethod(columns=_COLUMNS) + method = ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM) return method.estimate(MethodInput(scada_df=scada, **kwargs)).p50_overall def test_only_offered_references_are_used(self) -> None: @@ -1107,3 +1232,551 @@ def test_rows_a_reference_may_not_contribute_are_dropped(self) -> None: ) untouched = self._estimate(scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) assert with_holes != pytest.approx(untouched) + + +# --- pairing filter (pairing_max_gap) ----------------------------------------------------------- + + +class TestPairWithin: + """``pair_within`` keeps a row only when the other segment has a row within the gap, inclusive.""" + + @staticmethod + def _one_each(gap_rows: int) -> tuple[pd.DatetimeIndex, np.ndarray, np.ndarray]: + index = _index(gap_rows + 1) + baseline = np.zeros(len(index), dtype=bool) + upgraded = np.zeros(len(index), dtype=bool) + baseline[0] = True + upgraded[-1] = True + return index, baseline, upgraded + + def test_a_gap_exactly_equal_to_the_limit_is_kept(self) -> None: + index, baseline, upgraded = self._one_each(gap_rows=2) + paired = pair_within(index, baseline=baseline, upgraded=upgraded, max_gap=pd.Timedelta(minutes=20)) + assert paired.baseline[0] + assert paired.upgraded[-1] + + def test_a_gap_one_row_beyond_the_limit_is_dropped(self) -> None: + index, baseline, upgraded = self._one_each(gap_rows=3) + paired = pair_within(index, baseline=baseline, upgraded=upgraded, max_gap=pd.Timedelta(minutes=20)) + assert not paired.baseline.any() + assert not paired.upgraded.any() + + @pytest.mark.parametrize("unit", ["s", "ms", "us", "ns"]) + def test_the_index_resolution_does_not_change_the_gap(self, unit: str) -> None: + """Real SCADA often arrives at microsecond resolution; a 10-min gap must still read as 10 min.""" + index, baseline, upgraded = self._one_each(gap_rows=1) + index = index.as_unit(unit) + assert gap_to_other_segment(index, baseline=baseline, upgraded=upgraded)[0] == pytest.approx(600.0) + paired = pair_within(index, baseline=baseline, upgraded=upgraded, max_gap=pd.Timedelta(minutes=5)) + assert not paired.baseline.any() + + def test_an_empty_other_segment_drops_everything(self) -> None: + index = _index(6) + baseline = np.ones(len(index), dtype=bool) + paired = pair_within( + index, baseline=baseline, upgraded=np.zeros(len(index), dtype=bool), max_gap=pd.Timedelta(hours=1) + ) + assert not paired.baseline.any() + + @pytest.mark.parametrize("gap_minutes", [10, 30, 60]) + def test_matches_v0_any_within_timedelta(self, gap_minutes: int) -> None: + """Parity with wind-up v0's ``_toggle_pairing_filter`` on a gappy toggle frame.""" + v0 = pytest.importorskip("wind_up_v0.main_analysis") + rng = np.random.default_rng(3) + full = _index(400) + index = full[rng.random(len(full)) > 0.3] # drop rows so gaps vary + state = rng.random(len(index)) > 0.5 + cols = {"ws": 1.0, "pw": 1.0, "wd": 1.0} + pre_df = pd.DataFrame(cols, index=index[~state]) + post_df = pd.DataFrame(cols, index=index[state]) + v0_pre, v0_post = v0._toggle_pairing_filter( # noqa: SLF001 - the v0 behaviour being matched + pre_df=pre_df, + post_df=post_df, + pairing_filter_method="any_within_timedelta", + pairing_filter_timedelta_seconds=gap_minutes * 60, + detrend_ws_col="ws", + test_pw_col="pw", + ref_wd_col="wd", + timebase_s=600, + ) + + paired = pair_within(index, baseline=~state, upgraded=state, max_gap=pd.Timedelta(minutes=gap_minutes)) + assert set(index[paired.baseline]) == set(v0_pre.index) + assert set(index[paired.upgraded]) == set(v0_post.index) + + +class TestPairingFilter: + """``pairing_max_gap`` on the method: validation, the rows it drops, and the selection accounting.""" + + @staticmethod + def _block_excluded_case() -> tuple[pd.DataFrame, ToggleSchedule, pd.DatetimeIndex]: + """A noisy toggle campaign whose upgraded test rows are all flagged over one 1-day block.""" + scada, schedule = _noisy_toggle_case() + idx = pd.DatetimeIndex(pd.unique(scada.index)) + treated = np.asarray(treated_mask(idx, schedule)) + block = (idx >= idx[500]) & (idx < idx[644]) + flagged = pd.Series(block & treated, index=idx) + out = scada.copy() + is_test = out[_TURBINE_COL] == "T1" + out[_EXCLUDE_COL] = False + out.loc[is_test, _EXCLUDE_COL] = flagged.reindex(out.index[is_test]).to_numpy() + return out, schedule, idx[block] + + @staticmethod + def _run(scada: pd.DataFrame, schedule: ToggleSchedule, **kwargs: object) -> MethodOutput: + kwargs.setdefault("toggle_period", _CYCLE) + kwargs.setdefault("toggle_datum", _DATUM) + return ToggleSpecialistMethod(columns=_EXCLUDE_COLUMNS, **kwargs).estimate( + MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) + ) + + @pytest.mark.parametrize("gap", [pd.Timedelta(0), pd.Timedelta(minutes=5), pd.Timedelta(minutes=25)]) + def test_rejects_a_gap_that_is_not_a_positive_timebase_multiple(self, gap: pd.Timedelta) -> None: + scada, schedule = _noisy_toggle_case(n=200) + with pytest.raises(ValueError, match="pairing_max_gap"): + _estimate(scada, schedule, pairing_max_gap=gap) + + def test_baseline_rows_far_from_any_used_upgraded_row_are_dropped(self) -> None: + scada, schedule, block = self._block_excluded_case() + out = self._run(scada, schedule, pairing_max_gap=pd.Timedelta(minutes=20)) + assert out.labeled_rows is not None + rows = out.labeled_rows + inner = block[6:-6] # more than 20 min from the block's edges + assert not rows.loc[inner, "used"].astype(bool).any() + outside = rows.index < block[0] + assert rows.loc[outside, "used"].astype(bool).all() + + def test_without_pairing_the_block_keeps_its_baseline_rows(self) -> None: + scada, schedule, block = self._block_excluded_case() + out = self._run(scada, schedule) + assert out.labeled_rows is not None + rows = out.labeled_rows.loc[block] + assert rows.loc[rows["segment"] == "baseline", "used"].astype(bool).all() + + def test_accounting_shows_each_stage_per_segment(self) -> None: + scada, schedule, block = self._block_excluded_case() + out = self._run(scada, schedule, pairing_max_gap=pd.Timedelta(minutes=20)) + assert out.selection_accounting is not None + kept = out.selection_accounting.set_index(["stage", "segment"])["n_kept"] + n_block = len(block) // 2 + assert kept["exclude_row", "upgraded"] == kept["filters", "upgraded"] - n_block + assert kept["exclude_row", "baseline"] == kept["filters", "baseline"] + assert kept["pairing", "baseline"] < kept["exclude_row", "baseline"] + assert kept["pairing", "upgraded"] <= kept["exclude_row", "upgraded"] + + def test_accounting_is_reported_without_pairing(self) -> None: + scada, schedule = _noisy_toggle_case(n=200) + out = _estimate(scada, schedule) + assert out.selection_accounting is not None + kept = out.selection_accounting.set_index(["stage", "segment"])["n_kept"] + assert kept["pairing", "baseline"] == kept["exclude_row", "baseline"] + + +@pytest.mark.parametrize("gap", [None, pd.Timedelta(minutes=20)], ids=["no-pairing", "pairing"]) +def test_pairing_gap_plot_is_written(tmp_path: Path, gap: pd.Timedelta | None) -> None: + """The gap histogram is drawn with or without a pairing filter, so a user can judge whether one is needed.""" + scada, schedule = _noisy_toggle_case(n=200) + _estimate(scada, schedule, out_dir=tmp_path, save_plots=True, pairing_max_gap=gap) + assert len(list(tmp_path.rglob("*_pairing_gap.png"))) == 1 + + +# --- the two legs, their blend and the decisions on permuted labels ----------------------------- + +_BLOCK = pd.Timedelta(minutes=40) # one full 20-min on / 20-min off cycle = 4 rows of 10 min + + +def _block_case( + n_blocks: int = 120, + *, + uplift: float = 0.04, + ramp: float = 0.0, + ref_scale: np.ndarray | None = None, + seed: int = 0, +) -> tuple[pd.DataFrame, ToggleSchedule]: + """A toggle campaign whose power level jumps wildly between 40-min blocks (two on/off cycles each). + + Test power is ``level_b * shape_i``, times ``1 + uplift`` when treated, where ``shape`` is a + within-block ramp of slope ``ramp`` per row (flat by default: a ramp in phase with the toggle + is exactly the drift the block leg cannot cancel). With ``ref_scale`` (one factor per block) a + reference ``R1 = ref_scale_b * level_b * shape_i`` is added: it tracks the test turbine inside + every block but its ratio to the test changes from block to block. + """ + rng = np.random.default_rng(seed) + idx = _index(4 * n_blocks) + schedule = ToggleSchedule(period=pd.Timedelta(minutes=20), start=idx[0]) + treated = np.asarray(treated_mask(idx, schedule)) + level = np.repeat(rng.uniform(200.0, 2000.0, n_blocks), 4) + shape = np.tile(1.0 + ramp * np.arange(4), n_blocks) + base = level * shape + test = np.where(treated, base * (1.0 + uplift), base) + turbines = {"T1": test} + if ref_scale is not None: + turbines["R1"] = np.repeat(ref_scale, 4) * base + return _scada(turbines, idx), schedule + + +def _two_reference_case(noise_frac: float = 0.3, seed: int = 5) -> tuple[pd.DataFrame, ToggleSchedule]: + """``_block_case`` with a perfect tracker ``R1`` plus a noisy tracker ``R2``: the union is worse than R1 alone.""" + scada, schedule = _block_case(uplift=0.04, ref_scale=np.full(120, 0.5), seed=seed) + rng = np.random.default_rng(seed) + r2 = scada[scada[_TURBINE_COL] == "R1"].copy() + r2[_TURBINE_COL] = "R2" + r2[_POWER_COL] = r2[_POWER_COL].to_numpy() * (1.0 + noise_frac * rng.standard_normal(len(r2))) + return pd.concat([scada, r2]), schedule + + +def _legs(out: MethodOutput) -> set[str]: + return set(_long_diagnostics(out)["component"]) + + +class TestBinPower: + """``bin_power`` picks the per-row power the bins and ``power_band`` read.""" + + @staticmethod + def _bins(uplift: float, bin_power: str) -> pd.Series: + out = _estimate( + *_varying_rho_case(n=600, uplift=uplift), + conditions=("power",), + rated_power_kw=_RATED_KW, + bin_power=bin_power, + ) + assert out.labeled_rows is not None + return out.labeled_rows["power_bin"].astype(str) + + def test_baseline_bins_do_not_move_with_the_upgraded_rows(self) -> None: + assert (self._bins(0.0, "baseline") == self._bins(0.10, "baseline")).all() + + def test_both_states_bins_move_with_the_upgraded_rows(self) -> None: + """The control for the test above: the same change does move the both-states label.""" + assert (self._bins(0.0, "both_states") != self._bins(0.10, "both_states")).any() + + def test_an_unknown_bin_power_is_refused(self) -> None: + with pytest.raises(ValueError, match="bin_power"): + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, bin_power="mean") # type: ignore[arg-type] + + +class TestPowerBand: + """``power_band`` keeps the used rows whose bin power lies in ``(lo, hi]``.""" + + @pytest.mark.parametrize("band", [(800.0, 400.0), (400.0, 400.0), (0.0, float("inf"))]) + def test_a_band_that_is_not_a_finite_interval_is_refused(self, band: tuple[float, float]) -> None: + with pytest.raises(ValueError, match="power_band"): + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE, toggle_datum=_DATUM, power_band=band) + + _BAND = (360.0, 840.0) # two whole 240 kW bins at 1200 kW rated + + @staticmethod + def _run(scada: pd.DataFrame, schedule: ToggleSchedule, **kwargs: object) -> MethodOutput: + return _estimate(scada, schedule, conditions=("power",), rated_power_kw=_RATED_KW, **kwargs) + + def test_every_used_row_lies_in_the_band(self) -> None: + out = self._run(*_varying_rho_case(n=2016, uplift=0.05), power_band=self._BAND) + assert out.labeled_rows is not None + used_bins = set(out.labeled_rows.loc[out.labeled_rows["used"].astype(bool), "power_bin"].astype(str)) + assert used_bins == {"(360.0, 600.0]", "(600.0, 840.0]"} + + def test_a_band_holding_every_row_changes_nothing(self) -> None: + scada, schedule = _noisy_toggle_case() + assert self._run(scada, schedule, power_band=(-1e9, 1e9)).p50_overall == pytest.approx( + self._run(scada, schedule).p50_overall, abs=1e-12 + ) + + def test_a_constant_uplift_reads_the_same_inside_the_band(self) -> None: + """The cut reads the reference-derived power, which the state cannot move, so it selects no state.""" + out = self._run(*_varying_rho_case(n=2016, uplift=0.05), power_band=self._BAND) + assert out.p50_overall == pytest.approx(0.05, abs=2e-3) + + def test_the_band_is_a_selection_stage_and_in_the_run_config(self, tmp_path: Path) -> None: + # the run config is written with the plots + self._run(*_varying_rho_case(n=2016, uplift=0.05), power_band=self._BAND, out_dir=tmp_path, save_plots=True) + stages = set(_read_only_csv(tmp_path, "selection")["stage"]) + assert "power_band" in stages + run_dir = next(p for p in tmp_path.iterdir() if p.is_dir()) + config_text = next(run_dir.glob("config_*.yaml")).read_text() + assert "power_band" in config_text + + +class TestTogglePeriod: + """``toggle_period`` is the caller's statement of one toggle cycle: required, never inferred.""" + + def test_it_is_required(self) -> None: + with pytest.raises(TypeError, match="toggle_period"): + ToggleSpecialistMethod(columns=_COLUMNS) # type: ignore[call-arg] + + def test_it_is_independent_of_the_pairing_gap(self) -> None: + scada, schedule = _block_case(uplift=0.04, ref_scale=np.full(120, 0.5)) + out = _estimate(scada, schedule, toggle_period=_BLOCK, pairing_max_gap=pd.Timedelta(minutes=10)) + assert out.p50_overall == pytest.approx(0.04, abs=1e-9) + + def test_a_block_off_the_timebase_grid_is_refused(self) -> None: + scada, schedule = _block_case(n_blocks=10, ref_scale=np.ones(10)) + with pytest.raises(ValueError, match="toggle_period"): + _estimate(scada, schedule, toggle_period=pd.Timedelta(minutes=25)) + + +class TestToggleDatum: + """``toggle_datum`` is the caller's statement of the cycle grid: required, never taken from the data.""" + + def test_it_is_required(self) -> None: + with pytest.raises(TypeError, match="toggle_datum"): + ToggleSpecialistMethod(columns=_COLUMNS, toggle_period=_CYCLE) # type: ignore[call-arg] + + def test_cycles_stay_on_the_datum_grid_when_the_first_record_is_missing(self) -> None: + """Tiling from the first record present would put every cycle across two power levels.""" + scada, schedule = _block_case(uplift=0.04) + scada = scada[scada.index != scada.index.min()] + out = _estimate(scada, schedule, toggle_period=_BLOCK) + assert out.p50_overall == pytest.approx(0.04, abs=1e-3) + + def test_a_schedule_off_the_datum_grid_is_refused(self) -> None: + scada, schedule = _block_case(n_blocks=10) + with pytest.raises(ValueError, match="toggle_datum"): + _estimate(scada, schedule, toggle_period=_BLOCK, toggle_datum=_DATUM + pd.Timedelta(minutes=10)) + + def test_a_naive_datum_on_tz_aware_data_is_refused(self) -> None: + scada, schedule = _block_case(n_blocks=10) + with pytest.raises(ValueError, match="tz-aware"): + _estimate(scada, schedule, toggle_period=_BLOCK, toggle_datum=_DATUM.tz_localize(None)) + + +class TestBlockLeg: + """The block leg: the test turbine against its own block's mean power, needing no references.""" + + def test_recovers_the_uplift_with_no_references(self) -> None: + scada, schedule = _block_case(uplift=0.04) + out = _estimate(scada, schedule, toggle_period=_BLOCK) + assert out.p50_overall == pytest.approx(0.04, abs=1e-9) + assert np.isfinite(out.sigma_overall) + assert _legs(out) == {"block_mean"} + + def test_ignores_a_reference_outage(self) -> None: + """A reference the leg does not use must not cost it rows; the estimate is then this leg alone.""" + scada, schedule = _block_case(uplift=0.04, ref_scale=np.full(120, 0.5)) + is_ref = (scada[_TURBINE_COL] == "R1").to_numpy() + scada.loc[is_ref, _POWER_COL] = np.nan + out = _estimate(scada, schedule, toggle_period=_BLOCK) + assert out.p50_overall == pytest.approx(0.04, abs=1e-9) + assert out.labeled_rows is not None + assert out.labeled_rows["used"].astype(bool).all() + assert _legs(out) == {"block_mean"} + + def test_a_block_with_one_state_is_dropped_and_accounted(self) -> None: + scada, schedule = _block_case(uplift=0.04) + idx = pd.DatetimeIndex(pd.unique(scada.index)) + block_rows = idx[40:44] # the 11th block: two on/off cycles + off_rows = block_rows[~np.asarray(treated_mask(block_rows, schedule))] + scada.loc[scada.index.isin(off_rows), _POWER_COL] = np.nan + out = _estimate(scada, schedule, toggle_period=_BLOCK) + assert out.labeled_rows is not None + # NaN (a timestamp absent from the pivot) must read as unused, so compare rather than cast + used = out.labeled_rows["used"] == True # noqa: E712 + assert not used.loc[block_rows].any() + assert used.sum() == len(idx) - 4 + kept = out.selection_accounting.set_index(["stage", "segment"])["n_kept"] + assert kept["block", "upgraded"] == kept["pairing", "upgraded"] - 2 + assert kept["block", "baseline"] == kept["pairing", "baseline"] + + def test_the_reference_leg_keeps_every_block_without_flags(self) -> None: + scada, schedule = _noisy_toggle_case(n=200) + out = _estimate(scada, schedule, pairing_max_gap=pd.Timedelta(minutes=20)) + kept = out.selection_accounting.set_index(["stage", "segment"])["n_kept"] + stages = ["segment", "filters", "exclude_row", "pairing", "block"] + assert list(out.selection_accounting["stage"].unique()) == stages + assert kept["block", "baseline"] == kept["pairing", "baseline"] + + +class TestBlend: + """The estimate blends the reference and block legs by minimum-variance weights.""" + + def test_recovers_the_uplift_when_both_legs_are_exact(self) -> None: + scada, schedule = _block_case(uplift=0.04, ref_scale=np.full(120, 0.5)) + out = _estimate(scada, schedule, toggle_period=_BLOCK) + assert out.p50_overall == pytest.approx(0.04, abs=1e-9) + assert _legs(out) == {"sum", "block_mean", "combined"} + + def test_used_rows_are_the_union_of_the_legs(self) -> None: + scada, schedule = _block_case(uplift=0.04, ref_scale=np.full(120, 0.5)) + is_ref = (scada[_TURBINE_COL] == "R1").to_numpy() + scada.loc[is_ref & (np.arange(len(scada)) % 7 == 0), _POWER_COL] = np.nan + out = _estimate(scada, schedule, toggle_period=_BLOCK) + used = out.labeled_rows["used"] == True # noqa: E712 - NaN (absent timestamp) must read as unused + assert used.all() + # The accounting is the reference leg's, which the reference outage did cost rows. + kept = out.selection_accounting.set_index(["stage", "segment"])["n_kept"] + assert kept["block", "baseline"] + kept["block", "upgraded"] < used.sum() + + def test_diagnostics_carry_both_legs_and_the_blend(self) -> None: + scada, schedule = _noisy_toggle_case(n=400) + out = _estimate(scada, schedule, pairing_max_gap=pd.Timedelta(minutes=20)) + assert _legs(out) == {"sum", "block_mean", "combined"} + diag = _long_diagnostics(out) + head = diag[(diag["component"] == "combined") & (diag["condition_bin"] == "overall")].iloc[0] + assert 0.0 <= head["weight_sum"] <= 1.0 + assert np.isfinite(head["decision_sigma_sum"]) + assert np.isfinite(head["decision_sigma_block_mean"]) + assert head["uplift"] == pytest.approx(out.p50_overall) + legs = diag[(diag["component"] != "combined") & (diag["condition_bin"] == "overall")] + assert np.isfinite(legs["uplift"]).all() + w = head["weight_sum"] + by_leg = legs.set_index("component")["uplift"] + assert out.p50_overall == pytest.approx(w * by_leg["sum"] + (1 - w) * by_leg["block_mean"]) + + def test_the_blend_is_the_weighted_legs_and_its_sigma_comes_from_the_actual_bootstraps( + self, tmp_path: Path + ) -> None: + scada, schedule = _noisy_toggle_case(n=600) + out = _estimate(scada, schedule, pairing_max_gap=pd.Timedelta(minutes=20), out_dir=tmp_path) + results = _read_only_csv(tmp_path, "results").iloc[0] + diag = _long_diagnostics(out) + head = diag[(diag["component"] == "combined") & (diag["condition_bin"] == "overall")].iloc[0] + legs = { + m: diag[(diag["component"] == m) & (diag["condition_bin"] == "overall")].iloc[0] + for m in ("sum", "block_mean") + } + w = head["weight_sum"] + assert out.p50_overall == pytest.approx( + w * results["uplift_frc_sum"] + (1 - w) * results["uplift_frc_block_mean"] + ) + sa, sb, r = legs["sum"]["sigma"], legs["block_mean"]["sigma"], head["correlation"] + assert sa == pytest.approx(results["uplift_sigma_frc_sum"]) + assert sb == pytest.approx(results["uplift_sigma_frc_block_mean"]) + variance = w**2 * sa**2 + (1 - w) ** 2 * sb**2 + 2 * w * (1 - w) * r * sa * sb + assert out.sigma_overall == pytest.approx(np.sqrt(variance)) + # The weights were judged on the permuted sigmas, which are not the reported ones. + assert results["decision_sigma_frc_sum"] != results["uplift_sigma_frc_sum"] + assert results["decision_sigma_frc_block_mean"] != results["uplift_sigma_frc_block_mean"] + + +_BLOCK_COL = "cycle_bad" +_BLOCK_COLUMNS = replace(_COLUMNS, exclude_block=_BLOCK_COL) + + +def _flag_block(scada: pd.DataFrame, *, turbine: str, at: pd.Timestamp) -> pd.DataFrame: + """Return a copy of ``scada`` with an all-False ``exclude_block`` column set on one ``turbine`` row.""" + out = scada.copy() + out[_BLOCK_COL] = False + out.loc[((out[_TURBINE_COL] == turbine) & (out.index == at)).to_numpy(), _BLOCK_COL] = True + return out + + +class TestExcludeBlock: + """``columns.exclude_block``: a flagged row removes its whole toggle cycle, in both legs.""" + + @staticmethod + def _run(scada: pd.DataFrame, schedule: ToggleSchedule, **kwargs: object) -> MethodOutput: + kwargs.setdefault("toggle_datum", _DATUM) + return ToggleSpecialistMethod(columns=_BLOCK_COLUMNS, **kwargs).estimate( + MethodInput(scada_df=scada, test_wtg="T1", upgrade_timing=schedule, turbine_col=_TURBINE_COL) + ) + + def test_a_flagged_test_row_removes_its_whole_block(self) -> None: + scada, schedule = _block_case(uplift=0.04, ref_scale=np.full(120, 0.5)) + idx = pd.DatetimeIndex(pd.unique(scada.index)) + block_5 = idx[20:24] + out = self._run(_flag_block(scada, turbine="T1", at=block_5[1]), schedule, toggle_period=_BLOCK) + assert out.labeled_rows is not None + used = out.labeled_rows["used"].astype(bool) + assert not used.loc[block_5].any() + assert used.drop(block_5).all() + assert out.p50_overall == pytest.approx(0.04, abs=1e-9) + acc = out.selection_accounting.set_index(["stage", "segment"])["n_kept"] + for segment in ("baseline", "upgraded"): + assert acc.loc[("block", segment)] == acc.loc[("pairing", segment)] - 2 + + def test_a_flagged_reference_row_removes_the_block_from_the_reference_leg(self) -> None: + """The block leg uses no references, so it keeps the block and the headline ``used`` stays.""" + scada, schedule = _block_case(uplift=0.04, ref_scale=np.full(120, 0.5)) + idx = pd.DatetimeIndex(pd.unique(scada.index)) + block_7 = idx[28:32] + out = self._run(_flag_block(scada, turbine="R1", at=block_7[3]), schedule, toggle_period=_BLOCK) + assert out.labeled_rows is not None + assert out.labeled_rows["used"].astype(bool).all() + acc = out.selection_accounting.set_index(["stage", "segment"])["n_kept"] + for segment in ("baseline", "upgraded"): + assert acc.loc[("block", segment)] == acc.loc[("pairing", segment)] - 2 + assert out.p50_overall == pytest.approx(0.04, abs=1e-9) + + def test_no_flags_change_nothing(self) -> None: + scada, schedule = _block_case(uplift=0.04, ref_scale=np.full(120, 0.5)) + scada[_BLOCK_COL] = False + with_role = self._run(scada, schedule, toggle_period=_BLOCK) + without = _estimate(scada.drop(columns=_BLOCK_COL), schedule, toggle_period=_BLOCK) + assert with_role.p50_overall == pytest.approx(without.p50_overall) + assert with_role.sigma_overall == pytest.approx(without.sigma_overall) + + +class TestReferenceSubsets: + """The reference leg takes the reference subset with the smallest label-permuted sigma.""" + + def test_picks_the_subset_with_the_smallest_permuted_sigma(self, tmp_path: Path) -> None: + scada, schedule = _two_reference_case() + out = _estimate(scada, schedule, out_dir=tmp_path) + results = _read_only_csv(tmp_path, "results") + assert results.loc[0, "refs_used"] == "R1" + assert results.loc[0, "n_refs"] == 1 + assert results.loc[0, "n_refs_available"] == 2 + assert out.p50_overall == pytest.approx(0.04, abs=1e-9) + subsets = _read_only_csv(tmp_path, "reference_subsets") + assert subsets["decision_sigma_frc"].notna().all() + chosen = subsets.loc[subsets["chosen"], "refs"].item() + assert chosen == subsets.loc[subsets["decision_sigma_frc"].idxmin(), "refs"] + assert chosen == "R1" # the perfect tracker is the least noisy on any basis + + def test_every_subset_is_written_and_only_the_winner_has_an_actual_sigma(self, tmp_path: Path) -> None: + scada, schedule = _two_reference_case() + _estimate(scada, schedule, out_dir=tmp_path) + subsets = _read_only_csv(tmp_path, "reference_subsets").set_index("refs") + assert sorted(subsets.index) == ["R1", "R1;R2", "R2"] + assert subsets["chosen"].sum() == 1 + assert (subsets["n_used_timestamps"] > 0).all() + assert subsets["uplift_sigma_frc"].notna().sum() == 1 + assert np.isfinite(subsets.loc["R1", "uplift_sigma_frc"]) + + def test_the_reported_sigma_is_the_actual_one_of_the_chosen_subset(self) -> None: + scada, schedule = _two_reference_case() + chosen = _estimate(scada, schedule) + only_r1 = _estimate(scada[scada[_TURBINE_COL] != "R2"], schedule) + assert chosen.p50_overall == pytest.approx(only_r1.p50_overall) + assert chosen.sigma_overall == pytest.approx(only_r1.sigma_overall) + + def test_the_permuted_sigma_measures_the_same_noise_as_the_actual_one(self, tmp_path: Path) -> None: + """For the noisy reference the two sigmas agree; they differ only in whether the on/off split is real.""" + scada, schedule = _two_reference_case() + _estimate(scada, schedule, out_dir=tmp_path / "both") + _estimate(scada[scada[_TURBINE_COL] != "R1"], schedule, out_dir=tmp_path / "r2") + permuted = _read_only_csv(tmp_path / "both", "reference_subsets").set_index("refs")["decision_sigma_frc"] + actual_r2 = _read_only_csv(tmp_path / "r2", "results").loc[0, "uplift_sigma_frc_sum"] + assert permuted["R2"] == pytest.approx(actual_r2, rel=0.3) + assert permuted["R1"] < permuted["R2"] + + +class TestPermutedLabels: + """``permute_labels_within_blocks``: the same campaign with the on/off labels shuffled per block.""" + + def test_every_block_keeps_its_on_off_counts(self) -> None: + rng = np.random.default_rng(0) + n = 40 + upgraded = np.arange(n) % 4 < 2 # two on, two off per block of four + baseline = ~upgraded + used = np.ones(n, dtype=bool) + used[[3, 10, 11]] = False + ids = np.arange(n) // 4 + up_p, base_p = toggle_specialist.permute_labels_within_blocks( + upgraded=upgraded, baseline=baseline, used=used, ids=ids, rng=rng + ) + assert not np.array_equal(up_p, upgraded) + assert np.array_equal(up_p & used, ~base_p & used) # every used row is exactly one state + assert np.array_equal(up_p[~used], upgraded[~used]) # unused rows untouched + for block in np.unique(ids): + rows = (ids == block) & used + assert (up_p & rows).sum() == (upgraded & rows).sum() + assert (base_p & rows).sum() == (baseline & rows).sum() + + def test_without_blocks_the_labels_are_shuffled_globally(self) -> None: + rng = np.random.default_rng(1) + upgraded = np.arange(12) % 2 == 0 + up_p, base_p = toggle_specialist.permute_labels_within_blocks( + upgraded=upgraded, baseline=~upgraded, used=np.ones(12, dtype=bool), ids=None, rng=rng + ) + assert up_p.sum() == upgraded.sum() + assert np.array_equal(up_p, ~base_p) diff --git a/tests/benchmarking/campaigns/test_placebo_end_to_end.py b/tests/benchmarking/campaigns/test_placebo_end_to_end.py index 2aebcb19..d10912a6 100644 --- a/tests/benchmarking/campaigns/test_placebo_end_to_end.py +++ b/tests/benchmarking/campaigns/test_placebo_end_to_end.py @@ -43,15 +43,15 @@ def placebo_scada(mode: str, *, wind_noise_sd: float, seed: int = 0) -> pd.DataFrame: - """Hourly SCADA over ``mode``'s placebo period: a shared wind signal plus per-turbine noise. + """10-minute SCADA over ``mode``'s placebo period: a shared wind signal plus per-turbine noise. :param mode: the placebo mode whose analysis period the frame spans :param wind_noise_sd: per-turbine wind-speed noise in m/s; 0 makes every turbine identical :param seed: RNG seed for the noise """ rng = np.random.default_rng(seed) - index = pd.date_range(*placebo_analysis_period(mode), freq="1h", tz="UTC", inclusive="left") - hours = np.arange(len(index), dtype=float) + index = pd.date_range(*placebo_analysis_period(mode), freq="10min", tz="UTC", inclusive="left") + hours = np.arange(len(index), dtype=float) / 6 # a slow seasonal swing plus a daily cycle, so the pre and post periods differ in wind resource shared_ws = 8.0 + 2.5 * np.sin(hours / (24 * 30)) + 1.5 * np.sin(hours / 12) frames = [] @@ -68,7 +68,7 @@ def placebo_scada(mode: str, *, wind_noise_sd: float, seed: int = 0) -> pd.DataF HOT_COLUMNS.wind_speed: ws, HOT_COLUMNS.wind_speed_sd: 0.8, HOT_COLUMNS.gen_rpm: 1400.0, - HOT_COLUMNS.availability: 3600.0, + HOT_COLUMNS.availability: 600.0, }, index=index, ) diff --git a/tests/benchmarking/harness/test_replicates.py b/tests/benchmarking/harness/test_replicates.py index 6aa1ad72..dad2d32b 100644 --- a/tests/benchmarking/harness/test_replicates.py +++ b/tests/benchmarking/harness/test_replicates.py @@ -42,6 +42,7 @@ def _study(mode: str = "prepost", n_replicates: int = 5, seed: int = 0) -> Study min_pre_months=12, campaign_months=[3, 6], toggle_period=pd.Timedelta(days=14), + toggle_datum=pd.Timestamp("2017-01-01", tz="UTC"), n_replicates=n_replicates, seed=seed, ) @@ -55,6 +56,7 @@ def _weeks_study(mode: str = "toggle", n_replicates: int = 5, seed: int = 0) -> min_pre_months=12, campaign_weeks=[1, 2, 4, 8], toggle_period=pd.Timedelta(days=14), + toggle_datum=pd.Timestamp("2017-01-01", tz="UTC"), n_replicates=n_replicates, seed=seed, ) @@ -106,6 +108,32 @@ def test_max_activity_months_raises_for_a_weeks_study(self) -> None: _ = _weeks_study().max_activity_months +class TestToggleCycleGrid: + """A toggle study draws treatment starts only on its cycle grid, so one datum serves every replicate.""" + + def test_a_toggle_study_needs_a_datum(self) -> None: + with pytest.raises(ValueError, match="toggle_datum"): + StudyConfig( + mode="toggle", + turbine_subset=["T1"], + treatment_start_range=(pd.Timestamp("2017-01-01", tz="UTC"), pd.Timestamp("2017-12-31", tz="UTC")), + min_pre_months=12, + campaign_weeks=[2], + toggle_period=pd.Timedelta(days=14), + n_replicates=1, + ) + + def test_every_drawn_start_is_a_cycle_boundary(self) -> None: + study = _weeks_study(n_replicates=8) + reps = build_replicates(_base_scada(), profile=PROFILE, study=study) + assert study.toggle_period is not None + assert study.toggle_datum is not None + for rep in reps: + start = pd.Timestamp(rep.treatment_start) + start = start if start.tz is not None else start.tz_localize("UTC") + assert (start - study.toggle_datum) % study.toggle_period == pd.Timedelta(0) + + def test_build_replicates_returns_n_replicate_records() -> None: reps = build_replicates(_base_scada(), profile=PROFILE, study=_study(n_replicates=5)) assert len(reps) == 5 diff --git a/tests/benchmarking/harness/test_scoring.py b/tests/benchmarking/harness/test_scoring.py index 6e23aa24..12b978bb 100644 --- a/tests/benchmarking/harness/test_scoring.py +++ b/tests/benchmarking/harness/test_scoring.py @@ -46,6 +46,7 @@ def _study(mode: str = "prepost", n_replicates: int = 4, seed: int = 0) -> Study min_pre_months=12, campaign_months=[3, 6], toggle_period=pd.Timedelta(days=14), + toggle_datum=pd.Timestamp("2017-01-01", tz="UTC"), n_replicates=n_replicates, seed=seed, )