1554 lines
67 KiB
Python
1554 lines
67 KiB
Python
# WashData - Home Assistant integration for appliance cycle monitoring via smart plugs.
|
||
# Copyright (C) 2026 Lukas Bandura
|
||
# SPDX-License-Identifier: AGPL-3.0-or-later
|
||
#
|
||
# This program is free software: you can redistribute it and/or modify
|
||
# it under the terms of the GNU Affero General Public License as published
|
||
# by the Free Software Foundation, either version 3 of the License, or
|
||
# (at your option) any later version.
|
||
#
|
||
# This program is distributed in the hope that it will be useful,
|
||
# but WITHOUT ANY WARRANTY; without even the implied warranty of
|
||
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
||
# GNU Affero General Public License for more details.
|
||
#
|
||
# You should have received a copy of the GNU Affero General Public License
|
||
# along with this program. If not, see <https://www.gnu.org/licenses/>.
|
||
"""Analysis module for heavy CPU tasks (offloaded to executor)."""
|
||
from __future__ import annotations
|
||
|
||
import logging
|
||
from typing import Any, Optional
|
||
|
||
import math
|
||
|
||
import numpy as np
|
||
|
||
from .const import (
|
||
DEFAULT_DTW_BANDWIDTH,
|
||
DEFAULT_DTW_MODE,
|
||
DEFAULT_PROFILE_MATCH_MAX_DURATION_RATIO,
|
||
DEFAULT_PROFILE_MATCH_MIN_DURATION_RATIO,
|
||
MATCH_CORR_WEIGHT,
|
||
MATCH_DDTW_DIST_SCALE,
|
||
MATCH_DTW_BLEND,
|
||
MATCH_DTW_DIST_SCALE,
|
||
MATCH_DTW_ENSEMBLE_W,
|
||
MATCH_DTW_REFINE_TOP_N,
|
||
MATCH_DTW_RESAMPLE_N,
|
||
MATCH_DURATION_SCALE,
|
||
MATCH_DURATION_SCALE_OVERRUN,
|
||
MATCH_PREFIX_MIN_POINTS,
|
||
MATCH_PREFIX_SHAPE_MAX_RATIO,
|
||
MATCH_DURATION_WEIGHT,
|
||
MATCH_DURATION_WEIGHT_IN_PROGRESS,
|
||
MATCH_MIN_RATIO_GRACE_S,
|
||
MATCH_ENERGY_SCALE,
|
||
MATCH_ENERGY_REF_MIN_CYCLES,
|
||
MATCH_ENERGY_WEIGHT,
|
||
MATCH_KEEP_MIN_SCORE,
|
||
MATCH_MAE_PEAK_FLOOR,
|
||
MATCH_MAE_REF_PEAK,
|
||
MATCH_MAE_SCALE,
|
||
MAX_ALIGN_GRID_POINTS,
|
||
STAGE4_INTEGRATED_ENERGY_DEVICE_TYPES,
|
||
)
|
||
from .signal_processing import integrate_wh
|
||
|
||
|
||
def stage4_energy_mode(device_type: str | None) -> str:
|
||
"""Return the Stage-4 ``energy_mode`` for a device type.
|
||
|
||
``"integrated"`` for device types in
|
||
``STAGE4_INTEGRATED_ENERGY_DEVICE_TYPES`` (washing machine / washer-dryer),
|
||
where same-duration temperature/spin variants make integrated energy the
|
||
right discriminator; ``"mean"`` (the historical default) otherwise. Single
|
||
source of truth for the gate, used by the manager, Playground and matching
|
||
tuner so all three stay consistent with the live matcher.
|
||
"""
|
||
return "integrated" if device_type in STAGE4_INTEGRATED_ENERGY_DEVICE_TYPES else "mean"
|
||
|
||
|
||
def _agreement(
|
||
observed: float, expected: float, scale: float, gaussian: bool = False
|
||
) -> float:
|
||
"""1.0 when observed==expected, decaying with the log-ratio / scale.
|
||
|
||
Two kernels over the same log-ratio. The default is Lorentzian/Cauchy-like
|
||
(``1/(1+|x|)``), which decays slowly and so keeps a badly-sized candidate in
|
||
contention. ``gaussian`` (``exp(-x^2/2)``) decays fast and separates much
|
||
harder on size.
|
||
|
||
**This is a discrimination choice, not a density fit** (register item 307).
|
||
Measured on the corpus, the within-profile duration log-residual is decidedly
|
||
NOT normal - excess kurtosis 16.65, |z|>3 at 2.64% against the 0.27% a normal
|
||
predicts, and stripping the cycles that cannot match their own profile at all
|
||
(item 304) makes it heavier still, not lighter. The Gaussian nonetheless wins
|
||
at cycle end because those extreme cycles are unwinnable either way, so the
|
||
sharper penalty costs nothing on them and buys separation on the bulk.
|
||
|
||
The same sharpness is why it must NOT be used mid-cycle: there the observed
|
||
duration is a prefix, necessarily far below the profile mean, and a fast kernel
|
||
crushes the correct long candidate. The slow tail is what keeps it alive.
|
||
"""
|
||
if observed <= 0 or expected <= 0 or scale <= 0:
|
||
return 0.0
|
||
ratio = np.log(observed / expected)
|
||
if gaussian:
|
||
return float(np.exp(-0.5 * (ratio / scale) ** 2))
|
||
return 1.0 / (1.0 + abs(ratio) / scale)
|
||
|
||
_LOGGER = logging.getLogger(__name__)
|
||
ALIGNMENT_CONTEXT_BUFFER = 50
|
||
|
||
def find_best_alignment(
|
||
current_power: list[float] | np.ndarray,
|
||
sample_power: list[float] | np.ndarray,
|
||
dt: float = 1.0, # pylint: disable=unused-argument
|
||
corr_weight: float = MATCH_CORR_WEIGHT,
|
||
) -> tuple[float, dict[str, float], int]:
|
||
"""Find Best Alignment using Coarse-to-Fine Search (CPU Bound)."""
|
||
|
||
curr = np.array(current_power)
|
||
ref = np.array(sample_power)
|
||
|
||
n_curr = len(curr)
|
||
n_ref = len(ref)
|
||
|
||
# Guard: cross-correlation crashes on empty or single-element arrays.
|
||
if n_curr < 2 or n_ref < 2:
|
||
return 0.0, {"corr": 0.0, "mae_score": 0.0}, 0
|
||
|
||
# 1. Coarse Alignment (Cross-Correlation)
|
||
# Downsample for speed if arrays are large
|
||
ds_factor = 1
|
||
if n_curr > 200:
|
||
ds_factor = int(n_curr / 100)
|
||
|
||
if ds_factor > 1:
|
||
c_coarse = curr[::ds_factor]
|
||
r_coarse = ref[::ds_factor]
|
||
else:
|
||
c_coarse = curr
|
||
r_coarse = ref
|
||
|
||
# Standardize
|
||
if np.std(c_coarse) > 1e-6:
|
||
c_norm = (c_coarse - np.mean(c_coarse)) / np.std(c_coarse)
|
||
else:
|
||
c_norm = c_coarse
|
||
|
||
if np.std(r_coarse) > 1e-6:
|
||
r_norm = (r_coarse - np.mean(r_coarse)) / np.std(r_coarse)
|
||
else:
|
||
r_norm = r_coarse
|
||
|
||
# Cross correlation
|
||
correlation = np.correlate(c_norm, r_norm, mode="full")
|
||
lags = np.arange(-len(r_norm) + 1, len(c_norm))
|
||
|
||
best_idx = int(np.argmax(correlation))
|
||
best_lag_coarse = lags[best_idx]
|
||
|
||
best_offset = best_lag_coarse * ds_factor
|
||
|
||
# 2. Fine Refinement
|
||
window = 10 * ds_factor
|
||
min_off = max(-len(ref) + 1, best_offset - window)
|
||
max_off = min(len(curr), best_offset + window)
|
||
|
||
best_mae = float("inf")
|
||
final_offset = best_offset
|
||
|
||
# ~0.2 n^2 element work on a complete cycle (audit MATCH-CORE-11). It is
|
||
# byte-identical to the old ``np.mean(np.abs(c_seg - r_seg))``: the same
|
||
# subtract, abs and pairwise ``add.reduce`` divided by the count, written into
|
||
# one reused buffer. Allocating three fresh n-element arrays per offset, and
|
||
# np.mean's Python wrapper, cost more than the arithmetic. There is no exact
|
||
# shortcut for an L1 distance over shifts, and a 2-D batch over offsets was
|
||
# measured slower (the overlap length differs per offset). Any other dtype
|
||
# keeps the old expression (np.mean sums an integer array through a cast).
|
||
buf = (
|
||
np.empty(min(n_curr, n_ref))
|
||
if curr.dtype == np.float64 and ref.dtype == np.float64
|
||
else None
|
||
)
|
||
for off in range(int(min_off), int(max_off) + 1):
|
||
# intersection
|
||
c_start = max(0, off)
|
||
c_end = min(n_curr, n_ref + off)
|
||
length = c_end - c_start
|
||
if length < 10:
|
||
continue
|
||
r_start = max(0, -off)
|
||
|
||
if buf is None:
|
||
mae = np.mean(np.abs(curr[c_start:c_end] - ref[r_start:r_start + length]))
|
||
else:
|
||
seg = buf[:length]
|
||
np.subtract(curr[c_start:c_end], ref[r_start:r_start + length], out=seg)
|
||
np.abs(seg, out=seg)
|
||
mae = np.add.reduce(seg) / length
|
||
if mae < best_mae:
|
||
best_mae = mae
|
||
final_offset = off
|
||
|
||
# Calculate Final Score metrics
|
||
off = final_offset
|
||
c_start = max(0, off)
|
||
c_end = min(n_curr, n_ref + off)
|
||
r_start = max(0, -off)
|
||
r_end = min(n_ref, n_curr - off)
|
||
|
||
if (c_end - c_start) < 5:
|
||
return 0.0, {"mae": float(best_mae)}, final_offset
|
||
|
||
c_final = curr[c_start:c_end]
|
||
r_final = ref[r_start:r_end]
|
||
|
||
mae = np.mean(np.abs(c_final - r_final))
|
||
|
||
# Correlation
|
||
if np.std(c_final) > 1e-6 and np.std(r_final) > 1e-6:
|
||
corr = np.corrcoef(c_final, r_final)[0, 1]
|
||
else:
|
||
corr = 0.0
|
||
|
||
# Scale-invariant MAE: express the error relative to the current cycle's
|
||
# peak and calibrate to the legacy behaviour at MATCH_MAE_REF_PEAK (see
|
||
# const.py). On a complete cycle the peak is common to every candidate, so
|
||
# ranking is unaffected; in prefix mode the compared slice, and so its peak,
|
||
# depends on each candidate's span (deep-dive 02, audit MR-09).
|
||
current_peak = float(np.max(np.abs(curr))) if curr.size else 0.0
|
||
scaled_mae = mae * MATCH_MAE_REF_PEAK / max(current_peak, MATCH_MAE_PEAK_FLOOR)
|
||
mae_score = MATCH_MAE_SCALE / (MATCH_MAE_SCALE + scaled_mae)
|
||
score = (corr_weight * max(0.0, corr)) + ((1.0 - corr_weight) * mae_score)
|
||
|
||
return float(score), {"mae": float(mae), "corr": float(corr)}, final_offset
|
||
|
||
def _dtw_lite_scalar(x: np.ndarray, y: np.ndarray, n: int, m: int, w: int) -> float:
|
||
"""Verbatim original scalar fill for :func:`compute_dtw_lite` — kept as the
|
||
correctness reference and automatic fallback on unexpected errors."""
|
||
prev_row = np.full(m + 1, float("inf"))
|
||
curr_row = np.full(m + 1, float("inf"))
|
||
prev_row[0] = 0
|
||
for i in range(1, n + 1):
|
||
center = int(i * (m / n))
|
||
start_j = max(1, center - w)
|
||
end_j = min(m, center + w + 1)
|
||
curr_row.fill(float("inf"))
|
||
val_x = x[i - 1]
|
||
for j in range(start_j, end_j + 1):
|
||
cost = abs(float(val_x - y[j - 1]))
|
||
m1 = prev_row[j]
|
||
m2 = curr_row[j - 1]
|
||
m3 = prev_row[j - 1]
|
||
if m1 < m2:
|
||
best_prev = m1 if m1 < m3 else m3
|
||
else:
|
||
best_prev = m2 if m2 < m3 else m3
|
||
curr_row[j] = cost + best_prev
|
||
prev_row[:] = curr_row[:]
|
||
return float(prev_row[m])
|
||
|
||
|
||
def compute_dtw_lite(
|
||
x: np.ndarray, y: np.ndarray, band_width_ratio: float = 0.1,
|
||
derivative: bool = False,
|
||
) -> float:
|
||
"""
|
||
Compute DTW distance with Sakoe-Chiba band constraint.
|
||
Optimized 1D DP implementation. O(N*W).
|
||
|
||
When ``derivative`` is True this warps on the first derivative (slope) of the
|
||
two curves (Derivative DTW): alignment is driven by shape/transitions rather
|
||
than absolute power level, which is robust to amplitude offset and scale.
|
||
|
||
The inner loop operates on Python-native float lists (converted via ``.tolist()``
|
||
once per row) to avoid per-element NumPy scalar boxing overhead. The results
|
||
for each row are written back as a single slice assignment. For the typical
|
||
matching case (n=m=200, band=0.1 → w=20, ~41 cells/row) this is ~1.9× faster
|
||
than the original element-by-element NumPy indexing loop. The anti-diagonal
|
||
vectorized fill from :func:`_dtw_cost_banded` is NOT used here
|
||
because its per-diagonal Python setup overhead dominates for small n (it is
|
||
2× *slower* than the scalar loop for n=200 — the opposite of its large-n
|
||
envelope-rebuild behaviour where it wins by 1.6–8×).
|
||
"""
|
||
if derivative:
|
||
x = np.gradient(np.asarray(x, dtype=float)) if len(x) > 1 else np.asarray(x, dtype=float)
|
||
y = np.gradient(np.asarray(y, dtype=float)) if len(y) > 1 else np.asarray(y, dtype=float)
|
||
n, m = len(x), len(y)
|
||
if n == 0 or m == 0:
|
||
return float("inf")
|
||
|
||
xf = np.asarray(x, dtype=float)
|
||
yf = np.asarray(y, dtype=float)
|
||
|
||
w = max(1, int(min(n, m) * band_width_ratio))
|
||
|
||
try:
|
||
# Precompute band bounds for all rows (eliminates per-row int/max/min calls).
|
||
i_idx = np.arange(1, n + 1, dtype=float)
|
||
centers = (i_idx * (m / n)).astype(np.intp)
|
||
start_js = np.maximum(1, centers - w)
|
||
end_js = np.minimum(m, centers + w + 1)
|
||
|
||
# Convert y to a plain Python list once so that inner-loop element access
|
||
# is native float retrieval rather than NumPy scalar unboxing.
|
||
ylist = yf.tolist()
|
||
|
||
prev_row = np.full(m + 1, np.inf)
|
||
curr_row = np.full(m + 1, np.inf)
|
||
prev_row[0] = 0.0
|
||
|
||
for i in range(n):
|
||
sj = int(start_js[i])
|
||
ej = int(end_js[i])
|
||
curr_row[:] = np.inf
|
||
val_x = float(xf[i])
|
||
|
||
# Convert the relevant prev_row slice to Python lists once per row.
|
||
# prev_prev[k] == prev_row[sj - 1 + k] (diagonal predecessor of cell j=sj+k)
|
||
# prev_curr[k] == prev_row[sj + k] (up-predecessor of cell j=sj+k)
|
||
prev_prev = prev_row[sj - 1 : ej].tolist() # length = ej - sj + 1
|
||
prev_curr = prev_row[sj : ej + 1].tolist() # length = ej - sj + 1
|
||
|
||
row_vals: list[float] = []
|
||
prev_j_val = np.inf # curr_row[sj - 1] — left predecessor, maintained locally
|
||
for y_val, pr_j1, pr_j in zip(ylist[sj - 1 : ej], prev_prev, prev_curr):
|
||
cost = abs(val_x - y_val)
|
||
# min(up=pr_j, left=prev_j_val, diag=pr_j1)
|
||
best = pr_j if pr_j < prev_j_val else prev_j_val
|
||
if pr_j1 < best:
|
||
best = pr_j1
|
||
prev_j_val = cost + best
|
||
row_vals.append(prev_j_val)
|
||
|
||
curr_row[sj : ej + 1] = row_vals # single slice write
|
||
|
||
prev_row, curr_row = curr_row, prev_row # swap without copy
|
||
|
||
return float(prev_row[m])
|
||
except Exception: # pylint: disable=broad-exception-caught
|
||
# The scalar reference is byte-identical (proven by tests), so degrade to it on
|
||
# any unexpected error rather than propagating out of the unguarded Stage-3
|
||
# refinement loop in compute_matches_worker. Mirrors compute_dtw_path.
|
||
_LOGGER.debug("compute_dtw_lite vectorized path failed; using scalar fallback", exc_info=True)
|
||
return _dtw_lite_scalar(xf, yf, n, m, w)
|
||
|
||
def dtw_lite_batch(x: np.ndarray, y: np.ndarray, band_width_ratio: float) -> np.ndarray:
|
||
"""``compute_dtw_lite`` for ``k`` equal-length pairs at once (rows of ``x``, ``y``).
|
||
|
||
The Stage-3 refines are all 200x200 with one band, so the ``top_n x 2`` DP
|
||
fills run as one row scan (audit MR-05): within a row, ``cell[j] = min(c[j],
|
||
local[j] + cell[j-1])`` with ``c = local + min(up, diag)`` unrolls to ``S[j] +
|
||
cummin(c[k] - S[k])`` over the row's prefix sums ``S``. ~8x faster for 10 DTWs;
|
||
the prefix sums reorder the additions, so results agree to ~1e-15 relative,
|
||
not bit for bit.
|
||
"""
|
||
xs = np.asarray(x, dtype=float)
|
||
ys = np.asarray(y, dtype=float)
|
||
k, n = xs.shape
|
||
m = ys.shape[1]
|
||
if n == 0 or m == 0:
|
||
return np.full(k, np.inf)
|
||
w = max(1, int(min(n, m) * band_width_ratio))
|
||
centers = (np.arange(1, n + 1, dtype=float) * (m / n)).astype(np.intp)
|
||
start_js = np.maximum(1, centers - w)
|
||
end_js = np.minimum(m, centers + w + 1)
|
||
prev = np.full((k, m + 1), np.inf)
|
||
prev[:, 0] = 0.0
|
||
zeros = np.zeros((k, 1))
|
||
for i in range(n):
|
||
a, b = int(start_js[i]), int(end_js[i])
|
||
local = np.abs(xs[:, i:i + 1] - ys[:, a - 1:b])
|
||
up_diag = np.minimum(prev[:, a:b + 1], prev[:, a - 1:b])
|
||
csum = np.cumsum(local, axis=1)
|
||
before = np.concatenate((zeros, csum[:, :-1]), axis=1)
|
||
row = csum + np.minimum.accumulate(up_diag - before, axis=1)
|
||
cur = np.full((k, m + 1), np.inf)
|
||
cur[:, a:b + 1] = row
|
||
prev = cur
|
||
return prev[:, m]
|
||
|
||
|
||
def _stage3_scores_batched(
|
||
pairs: list[tuple[np.ndarray, np.ndarray]],
|
||
current_peak: float,
|
||
*,
|
||
dtw_mode: str,
|
||
dtw_bandwidth: float,
|
||
l1_scale: float,
|
||
ddtw_scale: float,
|
||
ensemble_w: float,
|
||
) -> list[float]:
|
||
"""Stage-3 scores for ``scaled`` / ``ddtw`` / ``ensemble`` in one batched DTW.
|
||
|
||
``pairs`` are each candidate's (current, sample) series already on the
|
||
``MATCH_DTW_RESAMPLE_N`` grid; the arithmetic is `_dtw_component_score`'s.
|
||
"""
|
||
level = dtw_mode in ("scaled", "ensemble")
|
||
deriv = dtw_mode in ("ddtw", "ensemble")
|
||
xs: list[np.ndarray] = []
|
||
ys: list[np.ndarray] = []
|
||
if level:
|
||
xs.extend(a for a, _ in pairs)
|
||
ys.extend(b for _, b in pairs)
|
||
if deriv:
|
||
xs.extend(np.gradient(a) for a, _ in pairs)
|
||
ys.extend(np.gradient(b) for _, b in pairs)
|
||
dists = dtw_lite_batch(np.vstack(xs), np.vstack(ys), dtw_bandwidth)
|
||
peak = max(current_peak, MATCH_MAE_PEAK_FLOOR)
|
||
|
||
def _score(dist: float, scale: float) -> float:
|
||
scaled = (dist / MATCH_DTW_RESAMPLE_N) * MATCH_MAE_REF_PEAK / peak
|
||
return scale / (scale + scaled)
|
||
|
||
k = len(pairs)
|
||
if dtw_mode == "ensemble":
|
||
return [
|
||
ensemble_w * _score(float(dists[i]), l1_scale)
|
||
+ (1.0 - ensemble_w) * _score(float(dists[k + i]), ddtw_scale)
|
||
for i in range(k)
|
||
]
|
||
scale = ddtw_scale if deriv else l1_scale
|
||
return [_score(float(d), scale) for d in dists]
|
||
|
||
|
||
def _resample_to(arr: np.ndarray, n: int) -> np.ndarray:
|
||
"""Linearly resample a 1-D array to exactly ``n`` points over its index span.
|
||
|
||
Used to put the current cycle and a profile sample onto one common grid
|
||
before DTW so the Sakoe-Chiba band width and the distance normalisation mean
|
||
the same thing regardless of each series' native sampling cadence/length.
|
||
"""
|
||
a = np.asarray(arr, dtype=float)
|
||
length = len(a)
|
||
if length == 0:
|
||
return np.zeros(n)
|
||
if length == n:
|
||
return a
|
||
return np.interp(np.linspace(0.0, 1.0, n), np.linspace(0.0, 1.0, length), a)
|
||
|
||
|
||
def _dtw_component_score(
|
||
curr_arr: np.ndarray,
|
||
sample_arr: np.ndarray,
|
||
current_peak: float,
|
||
band: float,
|
||
derivative: bool,
|
||
scale: float,
|
||
curr_resampled: np.ndarray | None = None,
|
||
) -> float:
|
||
"""DTW similarity in [0,1] for one candidate: resample both series to a
|
||
common grid, warp (level or derivative), and express the distance relative
|
||
to the current peak (behaviour-neutral at MATCH_MAE_REF_PEAK)."""
|
||
a = curr_resampled if curr_resampled is not None else _resample_to(curr_arr, MATCH_DTW_RESAMPLE_N)
|
||
b = _resample_to(sample_arr, MATCH_DTW_RESAMPLE_N)
|
||
dtw_dist = compute_dtw_lite(a, b, band_width_ratio=band, derivative=derivative)
|
||
norm_dist = dtw_dist / MATCH_DTW_RESAMPLE_N
|
||
scaled = norm_dist * MATCH_MAE_REF_PEAK / max(current_peak, MATCH_MAE_PEAK_FLOOR)
|
||
return scale / (scale + scaled)
|
||
|
||
|
||
def _stage3_dtw_score(
|
||
curr_arr: np.ndarray,
|
||
sample_arr: np.ndarray,
|
||
current_peak: float,
|
||
*,
|
||
dtw_mode: str,
|
||
dtw_bandwidth: float,
|
||
l1_scale: float,
|
||
ddtw_scale: float,
|
||
ensemble_w: float,
|
||
curr_resampled: np.ndarray | None = None,
|
||
) -> float:
|
||
"""The DTW score for one candidate: the four-way ``dtw_mode`` branch of the
|
||
Stage-3 refinement.
|
||
|
||
Lifted verbatim out of ``compute_matches_worker`` (for the #364 Stage-6 prefix
|
||
pass, removed in 0.5.8). Behaviour-identical to the inlined version, including ``legacy`` mode's
|
||
``dtw_dist / len(curr_arr)`` normalisation. (The normalised distance it also
|
||
returned until audit MR-12 went to a ``dtw_dist`` candidate field nothing
|
||
read, always 0.0 under the default ``ensemble``.)
|
||
"""
|
||
if dtw_mode == "legacy":
|
||
# Original behaviour: raw sequences, distance / len(current),
|
||
# fixed absolute-watt scale (not peak-relative).
|
||
dtw_dist = compute_dtw_lite(curr_arr, sample_arr, band_width_ratio=dtw_bandwidth)
|
||
n_points = len(curr_arr)
|
||
norm_dist = (dtw_dist / n_points) if n_points > 0 else 999.0
|
||
return 1.0 / (1.0 + norm_dist / MATCH_DTW_DIST_SCALE)
|
||
if dtw_mode == "ensemble":
|
||
# Blend the level-based (L1) and shape-based (derivative) DTW
|
||
# scores; they are complementary signals.
|
||
s_l1 = _dtw_component_score(curr_arr, sample_arr, current_peak, dtw_bandwidth, False, l1_scale, curr_resampled=curr_resampled)
|
||
s_dd = _dtw_component_score(curr_arr, sample_arr, current_peak, dtw_bandwidth, True, ddtw_scale, curr_resampled=curr_resampled)
|
||
return ensemble_w * s_l1 + (1.0 - ensemble_w) * s_dd
|
||
# "scaled" (default) or "ddtw": resample both onto one grid so the
|
||
# band and normalisation are consistent, then express the distance
|
||
# relative to the current peak (behaviour-neutral at
|
||
# MATCH_MAE_REF_PEAK), mirroring the Stage-2 MAE treatment.
|
||
use_deriv = dtw_mode == "ddtw"
|
||
scale = ddtw_scale if use_deriv else l1_scale
|
||
return _dtw_component_score(
|
||
curr_arr, sample_arr, current_peak, dtw_bandwidth, use_deriv, scale, curr_resampled=curr_resampled
|
||
)
|
||
|
||
|
||
def _rank_key(candidate: dict[str, Any]) -> float:
|
||
"""Sort key for candidate ranking: a non-finite score ranks last, never first."""
|
||
try:
|
||
score = float(candidate.get("score", 0.0))
|
||
except (TypeError, ValueError, OverflowError):
|
||
return float("-inf")
|
||
return score if math.isfinite(score) else float("-inf")
|
||
|
||
|
||
def compute_matches_worker(
|
||
current_power: list[float],
|
||
current_duration: float,
|
||
snapshots: list[dict[str, Any]],
|
||
config: dict[str, Any]
|
||
) -> list[dict[str, Any]]:
|
||
"""Worker function to compute matches against snapshots."""
|
||
candidates: list[dict[str, Any]] = []
|
||
|
||
min_duration_ratio = config.get("min_duration_ratio", DEFAULT_PROFILE_MATCH_MIN_DURATION_RATIO)
|
||
max_duration_ratio = config.get("max_duration_ratio", DEFAULT_PROFILE_MATCH_MAX_DURATION_RATIO)
|
||
# Falls back to the OPTION default like every other key here. It used to be a
|
||
# literal 0.1, which silently ran a weaker Stage 3 than production for any
|
||
# caller that omitted the key - worth ~0.8pp of top-1, and it is what made
|
||
# the on-device matching tuner optimise against a pipeline nobody runs
|
||
# (register item 309).
|
||
dtw_bandwidth = config.get("dtw_bandwidth", DEFAULT_DTW_BANDWIDTH)
|
||
dtw_mode = config.get("dtw_mode", DEFAULT_DTW_MODE)
|
||
keep_min = float(config.get("keep_min_score", MATCH_KEEP_MIN_SCORE))
|
||
corr_weight = float(config.get("corr_weight", MATCH_CORR_WEIGHT))
|
||
dur_weight = float(config.get("duration_weight", MATCH_DURATION_WEIGHT))
|
||
if config.get("in_progress"):
|
||
dur_weight = float(
|
||
config.get("duration_weight_in_progress", MATCH_DURATION_WEIGHT_IN_PROGRESS)
|
||
)
|
||
en_weight = float(config.get("energy_weight", MATCH_ENERGY_WEIGHT))
|
||
dur_scale = float(config.get("duration_scale", MATCH_DURATION_SCALE))
|
||
dur_overrun_scale = float(
|
||
config.get("duration_overrun_scale", MATCH_DURATION_SCALE_OVERRUN)
|
||
)
|
||
en_scale = float(config.get("energy_scale", MATCH_ENERGY_SCALE))
|
||
# #400: while the cycle is still running, Stage 4 compares like with like. Both
|
||
# of its terms describe the cycle SO FAR; without this they are graded against
|
||
# each candidate's COMPLETE duration and energy, which makes a long program 40%
|
||
# through numerically indistinguishable from a finished short one. Opt-in (live
|
||
# match path only) so the final match at cycle end - where the whole-cycle
|
||
# figures are the right comparison - is byte-identical.
|
||
in_progress = bool(config.get("in_progress"))
|
||
min_ratio_grace_s = float(config.get("min_ratio_grace_s", MATCH_MIN_RATIO_GRACE_S))
|
||
# ...and Stages 2/3 score the SHAPE against the same truncated stretch, while the
|
||
# cycle is still clearly mid-run (MATCH_PREFIX_SHAPE_MAX_RATIO). On by default
|
||
# for a live match; `prefix_shape: False` turns it off for the A/B harnesses.
|
||
prefix_shape_on = in_progress and bool(config.get("prefix_shape", True))
|
||
prefix_shape_max_ratio = float(
|
||
config.get("prefix_shape_max_ratio", MATCH_PREFIX_SHAPE_MAX_RATIO)
|
||
)
|
||
|
||
curr_arr = np.array(current_power)
|
||
|
||
for item in snapshots:
|
||
name = item["name"]
|
||
profile_duration = item["avg_duration"]
|
||
sample_power = item["sample_power"]
|
||
|
||
# Duration Check. The lower bound waits out MATCH_MIN_RATIO_GRACE_S of a live
|
||
# match: early on it rejects every programme longer than the cycle is old.
|
||
if profile_duration > 0:
|
||
ratio = current_duration / profile_duration
|
||
lower_applies = not (
|
||
in_progress and current_duration < min_ratio_grace_s
|
||
)
|
||
if (lower_applies and ratio < min_duration_ratio) or ratio > max_duration_ratio:
|
||
continue
|
||
|
||
# Core Similarity. While the cycle is mid-run this compares it against the
|
||
# candidate truncated to the elapsed time, on a shared grid (#400) - the
|
||
# same pair Stage 3 then warps, so both shape stages ask one question.
|
||
# `sample` below stays the FULL template: Stage 4 takes its own prefix of it
|
||
# (analysis.prefix_mean), and truncating it here would truncate twice.
|
||
span_s = float(item.get("sample_span_s") or profile_duration or 0.0)
|
||
shape_pair = None
|
||
if (
|
||
prefix_shape_on
|
||
and span_s > 0
|
||
and current_duration <= span_s * prefix_shape_max_ratio
|
||
):
|
||
shape_pair = prefix_shape_arrays(
|
||
curr_arr, sample_power, current_duration, span_s
|
||
)
|
||
# The alignment offset is not kept: nothing read it, and it was in different
|
||
# units on the two paths (shared-grid index vs native index; audit MR-12).
|
||
if shape_pair is not None:
|
||
score, metrics, _offset = find_best_alignment(
|
||
shape_pair[0], shape_pair[1], 1.0, corr_weight=corr_weight
|
||
)
|
||
else:
|
||
score, metrics, _offset = find_best_alignment(
|
||
current_power, sample_power, 1.0, corr_weight=corr_weight
|
||
)
|
||
|
||
if score > keep_min:
|
||
candidates.append({
|
||
"name": name,
|
||
"score": score,
|
||
"metrics": metrics,
|
||
"profile_duration": profile_duration,
|
||
"current": current_power,
|
||
"sample": sample_power,
|
||
# Transient, popped after Stage 3: the truncated (current, template)
|
||
# pair Stage 2 scored, so Stage 3 warps the same thing.
|
||
"_shape_pair": shape_pair,
|
||
# True wall-clock span of `sample`, for prefix truncation (#364).
|
||
# Falls back to profile_duration so the other snapshot builders
|
||
# (devtools, playground) keep working unchanged.
|
||
"sample_span_s": float(item.get("sample_span_s") or profile_duration or 0.0),
|
||
})
|
||
|
||
candidates.sort(key=_rank_key, reverse=True)
|
||
|
||
# Stage 3: DTW Refinement on the top N candidates
|
||
if dtw_bandwidth > 0.0 and len(candidates) > 0:
|
||
# top-N, blend and the distance scales are config-overridable so the
|
||
# tuning harness can sweep them without editing constants; production
|
||
# uses the const defaults.
|
||
top_n = int(config.get("dtw_refine_top_n", MATCH_DTW_REFINE_TOP_N))
|
||
blend = float(config.get("dtw_blend", MATCH_DTW_BLEND))
|
||
to_refine = candidates[:top_n]
|
||
current_peak = float(np.max(curr_arr)) if curr_arr.size else 0.0
|
||
l1_scale = float(config.get("dtw_l1_scale", MATCH_DTW_DIST_SCALE))
|
||
ddtw_scale = float(config.get("dtw_ddtw_scale", MATCH_DDTW_DIST_SCALE))
|
||
ensemble_w = float(config.get("dtw_ensemble_w", MATCH_DTW_ENSEMBLE_W))
|
||
# Resample the current trace once — it's the same for every candidate.
|
||
curr_resampled = _resample_to(curr_arr, MATCH_DTW_RESAMPLE_N)
|
||
|
||
batched: list[float] | None = None
|
||
if dtw_mode in ("scaled", "ddtw", "ensemble") and to_refine:
|
||
# One batched DTW for every refine (audit MR-05: Stage 3 was ~89% of
|
||
# matcher CPU). Falls back to the per-candidate path on any error.
|
||
try:
|
||
grid_pairs = []
|
||
for cand in to_refine:
|
||
pair = cand.get("_shape_pair")
|
||
if pair is not None:
|
||
a = _resample_to(pair[0], MATCH_DTW_RESAMPLE_N)
|
||
b = _resample_to(pair[1], MATCH_DTW_RESAMPLE_N)
|
||
else:
|
||
a = curr_resampled
|
||
b = _resample_to(np.array(cand["sample"]), MATCH_DTW_RESAMPLE_N)
|
||
grid_pairs.append((a, b))
|
||
batched = _stage3_scores_batched(
|
||
grid_pairs, current_peak, dtw_mode=dtw_mode,
|
||
dtw_bandwidth=dtw_bandwidth, l1_scale=l1_scale,
|
||
ddtw_scale=ddtw_scale, ensemble_w=ensemble_w,
|
||
)
|
||
except Exception: # pylint: disable=broad-exception-caught
|
||
_LOGGER.debug("batched Stage-3 DTW failed; per-candidate path", exc_info=True)
|
||
batched = None
|
||
|
||
for idx, cand in enumerate(to_refine):
|
||
if batched is not None:
|
||
cand["original_score"] = float(cand["score"])
|
||
cand["score"] = float(blend * cand["score"] + (1.0 - blend) * batched[idx])
|
||
continue
|
||
pair = cand.get("_shape_pair")
|
||
if pair is not None:
|
||
warp_curr, sample_arr, cand_resampled = pair[0], pair[1], None
|
||
else:
|
||
warp_curr, sample_arr, cand_resampled = (
|
||
curr_arr, np.array(cand["sample"]), curr_resampled
|
||
)
|
||
|
||
dtw_score = _stage3_dtw_score(
|
||
warp_curr,
|
||
sample_arr,
|
||
current_peak,
|
||
dtw_mode=dtw_mode,
|
||
dtw_bandwidth=dtw_bandwidth,
|
||
l1_scale=l1_scale,
|
||
ddtw_scale=ddtw_scale,
|
||
ensemble_w=ensemble_w,
|
||
curr_resampled=cand_resampled,
|
||
)
|
||
|
||
cand["original_score"] = float(cand["score"])
|
||
cand["score"] = float(blend * cand["score"] + (1.0 - blend) * dtw_score)
|
||
|
||
candidates.sort(key=_rank_key, reverse=True)
|
||
|
||
# The truncated pair is scratch for the two shape stages; it must not reach the
|
||
# MatchResult ranking (numpy arrays, and the store/WS serialise that dict).
|
||
for cand in candidates:
|
||
cand.pop("_shape_pair", None)
|
||
|
||
# Final pass: blend in duration + energy agreement. Shape correlation alone
|
||
# cannot separate profiles that differ mainly in duration/energy (the main
|
||
# multi-program washing-machine failure mode), so nudge the score toward
|
||
# candidates whose expected duration/energy match the observed cycle.
|
||
# Sanitize the configured weights so the blended score stays a convex
|
||
# combination in [0, 1]: clamp negatives to 0 and, if duration+energy exceed
|
||
# 1.0, scale them down proportionally (shape then contributes 0) rather than
|
||
# letting shape_w go negative or the total exceed 1.
|
||
# Drop non-finite configured weights (NaN/inf) so de_sum, the normalized
|
||
# weights, and every candidate score stay finite.
|
||
dur_w = max(0.0, dur_weight) if np.isfinite(dur_weight) else 0.0
|
||
en_w = max(0.0, en_weight) if np.isfinite(en_weight) else 0.0
|
||
de_sum = dur_w + en_w
|
||
if de_sum > 1.0:
|
||
dur_w, en_w = dur_w / de_sum, en_w / de_sum
|
||
shape_w = max(0.0, 1.0 - dur_w - en_w)
|
||
if (dur_w > 0 or en_w > 0) and candidates and current_duration > 0:
|
||
# energy_mode: "mean" (default) compares whole-cycle mean power (W);
|
||
# "integrated" compares true integrated energy (mean x duration). Opt-in so
|
||
# the historical default is byte-for-byte preserved. See register item 99.
|
||
integrated = config.get("energy_mode", "mean") == "integrated"
|
||
own_energy = (
|
||
{str(s.get("name")): s.get("energy_ref") for s in snapshots}
|
||
if integrated and not in_progress else {}
|
||
)
|
||
cur_mean = float(np.mean(curr_arr))
|
||
cur_energy = cur_mean * current_duration if integrated else cur_mean
|
||
for cand in candidates:
|
||
prof_dur = float(cand.get("profile_duration") or 0.0)
|
||
if in_progress and prof_dur > 0 and current_duration > prof_dur:
|
||
# The cycle has outlasted this candidate: real evidence against it,
|
||
# penalised on the sharper scale. Below a candidate's duration the
|
||
# term is unchanged - suppressing the penalty there was measured and
|
||
# rejected (see MATCH_DURATION_SCALE_OVERRUN in const.py).
|
||
dur_ag = _agreement(current_duration, prof_dur, dur_overrun_scale)
|
||
else:
|
||
# Gaussian only on a COMPLETED cycle, where the observed duration
|
||
# is the real one. Mid-cycle it is a prefix and the sharp kernel
|
||
# costs -6.4pp at 60% elapsed (item 307).
|
||
dur_ag = _agreement(
|
||
current_duration, prof_dur, dur_scale, gaussian=not in_progress
|
||
)
|
||
sample = cand.get("sample") or []
|
||
if in_progress:
|
||
cand_mean, cand_span = prefix_mean(
|
||
sample,
|
||
current_duration,
|
||
float(cand.get("sample_span_s") or prof_dur or 0.0),
|
||
prof_dur,
|
||
)
|
||
else:
|
||
cand_mean = float(np.mean(sample)) if sample else 0.0
|
||
cand_span = prof_dur
|
||
# In integrated mode the candidate figure must cover the same stretch of
|
||
# time as cur_energy does, so a prefix mean is scaled by the elapsed
|
||
# duration and a whole-template mean by the candidate's own duration.
|
||
cand_energy = cand_mean * cand_span if integrated else cand_mean
|
||
if integrated and not in_progress:
|
||
# A COMPLETED cycle is graded against the median energy of the
|
||
# profile's own cycles, not mean(template) x duration: on a warped
|
||
# envelope that inherits the reference cycle's heating length (w7g).
|
||
# Not while running: rescaling the prefix by median/template gained
|
||
# 0.8-1.6pp top-1 at 25-75% elapsed, but the live matches it moved
|
||
# are the ones the end gates read (end_gate_eval --loo: 2 new ends
|
||
# > 5 min early, washer-dryer median lag +0.7 min).
|
||
own = own_energy_ws(own_energy.get(str(cand["name"])))
|
||
if own is not None:
|
||
cand_energy = own
|
||
en_ag = _agreement(cur_energy, cand_energy, en_scale)
|
||
cand["shape_score"] = float(cand["score"])
|
||
cand["score"] = float(
|
||
shape_w * cand["score"]
|
||
+ dur_w * dur_ag
|
||
+ en_w * en_ag
|
||
)
|
||
candidates.sort(key=_rank_key, reverse=True)
|
||
|
||
# A non-finite score (a NaN in an imported template, a NaN avg_duration) used to
|
||
# sort to rank 1 with a NaN margin that was never "ambiguous" (audit
|
||
# MATCH-CORE-05). It carries no evidence: drop it.
|
||
candidates[:] = [c for c in candidates if math.isfinite(_rank_key(c))]
|
||
|
||
# (Stage 6, the #364 prefix scores for the Smart-Termination guard, was removed
|
||
# in 0.5.8: Stages 2/3 already score a running cycle on each candidate's
|
||
# truncated curve, so the guard fired on 0 of 713 genuine ends and 0 of 7
|
||
# quiet split positives on the shipped matcher - devtools/prefix_guard_eval.py.)
|
||
return candidates
|
||
|
||
def prefix_mean(
|
||
sample: list[float] | np.ndarray,
|
||
current_duration: float,
|
||
sample_span_s: float,
|
||
profile_duration: float,
|
||
) -> tuple[float, float]:
|
||
"""Mean power of ``sample`` over its leading ``current_duration`` seconds, and
|
||
the span that mean covers - the Stage-4 like-for-like pair for a cycle that is
|
||
still running (#400).
|
||
|
||
The Stage-4 like-for-like pair, and the same pair Stage 5 uses to choose between
|
||
a group's members while a cycle is running - one definition, two callers.
|
||
|
||
Falls back to the whole template (and the candidate's own duration) when the
|
||
cycle has already outlasted it: that candidate has finished, so its total is
|
||
the honest comparison. Deliberately NOT ``_prefix_point_count``: that helper's
|
||
12-sample floor exists because Stages 2/3 *correlate* the prefix, while a mean
|
||
over a handful of leading samples is perfectly well defined - applying the
|
||
floor here would silently restore whole-template energy for the first few
|
||
percent of every cycle, which is exactly the window #400 is about.
|
||
"""
|
||
if len(sample) == 0:
|
||
return 0.0, profile_duration
|
||
if sample_span_s > 0 and current_duration < sample_span_s:
|
||
k = max(1, int(round(len(sample) * (current_duration / sample_span_s))))
|
||
return float(np.mean(sample[:k])), current_duration
|
||
return float(np.mean(sample)), profile_duration
|
||
|
||
|
||
def member_energy_reference(
|
||
members: list[tuple[list[float], list[float], float]],
|
||
) -> dict[str, Any] | None:
|
||
"""A profile's energy taken from its OWN cycles, for Stage 4.
|
||
|
||
``members`` are the ``(offsets, watts, duration)`` triples the envelope is built
|
||
from. Each member's energy is measured the way Stage 4 measures the cycle being
|
||
matched - mean of the linearly interpolated trace x duration - so both sides of
|
||
the agreement are the same quantity. ``{"n": members used, "median_wh": median
|
||
whole-cycle Wh}``, or None without a usable member.
|
||
|
||
Why (w7g): Stage 4 took a profile's energy as ``mean(template) x avg_duration``.
|
||
On a DTW-warped envelope that inherits the heating length of the cycle the
|
||
members were warped onto: > 20% off its own members' median on 17 of 85 washer
|
||
profiles in the corpus, 2.5x on one.
|
||
"""
|
||
totals: list[float] = []
|
||
for offsets, watts, duration in members:
|
||
t = np.asarray(offsets, dtype=float)
|
||
p = np.asarray(watts, dtype=float)
|
||
if t.size != p.size or t.size < 2:
|
||
continue
|
||
ok = np.isfinite(t) & np.isfinite(p)
|
||
t, p = t[ok], p[ok]
|
||
if t.size < 2:
|
||
continue
|
||
order = np.argsort(t, kind="mergesort")
|
||
t, p = t[order], p[order]
|
||
span = float(t[-1] - t[0])
|
||
if span <= 0:
|
||
continue
|
||
try:
|
||
dur = float(duration)
|
||
except (TypeError, ValueError, OverflowError):
|
||
dur = 0.0
|
||
if not math.isfinite(dur) or dur <= 0:
|
||
dur = span
|
||
energy_ws = integrate_wh(t, p) * 3600.0 / span * dur
|
||
if math.isfinite(energy_ws) and energy_ws > 0:
|
||
totals.append(energy_ws)
|
||
if not totals:
|
||
return None
|
||
return {"n": len(totals), "median_wh": round(float(np.median(totals)) / 3600.0, 3)}
|
||
|
||
|
||
def own_energy_ws(ref: dict[str, Any] | None) -> float | None:
|
||
"""The whole-cycle energy (W*s) Stage 4 expects from a profile's own cycles.
|
||
|
||
The median of a :func:`member_energy_reference`, or None (keep the template's
|
||
``mean x duration``) without one of at least ``MATCH_ENERGY_REF_MIN_CYCLES``
|
||
cycles or with a malformed one.
|
||
"""
|
||
if not isinstance(ref, dict):
|
||
return None
|
||
try:
|
||
if int(ref.get("n") or 0) < MATCH_ENERGY_REF_MIN_CYCLES:
|
||
return None
|
||
own_ws = float(ref["median_wh"]) * 3600.0
|
||
except (TypeError, ValueError, KeyError, OverflowError):
|
||
return None
|
||
return own_ws if math.isfinite(own_ws) and own_ws > 0 else None
|
||
|
||
|
||
def _prefix_point_count(
|
||
n_points: int, current_duration: float, sample_span_s: float
|
||
) -> int:
|
||
"""Leading template samples that cover ``current_duration`` seconds.
|
||
|
||
0 when the span is unknown/non-positive, when the elapsed time already covers
|
||
the whole template (then it is not a prefix), or when too few points remain to
|
||
judge. Fraction-of-array is the right operator because every snapshot flavour
|
||
is uniform in time over its own span (envelope: np.linspace; sample cycle:
|
||
resample_uniform at a fixed dt). (The group aggregate snapshot is gone, #400.)
|
||
"""
|
||
if n_points < MATCH_PREFIX_MIN_POINTS or sample_span_s <= 0 or current_duration <= 0:
|
||
return 0
|
||
k = int(round(n_points * (current_duration / sample_span_s)))
|
||
if k < MATCH_PREFIX_MIN_POINTS or k >= n_points:
|
||
return 0
|
||
return k
|
||
|
||
|
||
def prefix_shape_arrays(
|
||
curr_arr: np.ndarray,
|
||
sample: list[float] | np.ndarray,
|
||
current_duration: float,
|
||
sample_span_s: float,
|
||
) -> tuple[np.ndarray, np.ndarray] | None:
|
||
"""``(current, template)`` on one grid, with the template TRUNCATED to the
|
||
elapsed time - or None when it cannot be truncated meaningfully.
|
||
|
||
One definition of "the same stretch of both curves" for the live Stage-2/3
|
||
shape scoring (#400). (It also fed the #364 Stage-6 prefix guard, removed in
|
||
0.5.8; ``devtools/prefix_guard_eval.py`` keeps a copy of that rule.) Both
|
||
series go onto a shared grid so an index offset equals a time offset
|
||
regardless of the template's native cadence; the grid also honours the #388
|
||
OOM cap. The 12-sample floor is real here (unlike in ``prefix_mean``): these
|
||
arrays get correlated and warped, not averaged.
|
||
|
||
The grid is shared **between the two series**, not across candidates: ``k``
|
||
is the candidate's own truncated point count, so a longer template can be
|
||
scored on a finer grid than a shorter one. Capping the grid at ``k`` is what
|
||
stops ``arr[:k]`` being upsampled past the points it actually has, which
|
||
would invent template detail the recording never contained. (The one A/B of
|
||
dropping the ``k`` term was measured on the removed guard, not on matching.)
|
||
``devtools/dtw_ab_eval.py`` scores only complete cycles, which never take
|
||
the prefix path, so it cannot judge this function - use ``devtools/eval.py``.
|
||
"""
|
||
arr = np.asarray(sample, dtype=float)
|
||
k = _prefix_point_count(arr.size, current_duration, sample_span_s)
|
||
if k == 0:
|
||
return None
|
||
grid = int(min(curr_arr.size, k, MAX_ALIGN_GRID_POINTS))
|
||
if grid < MATCH_PREFIX_MIN_POINTS:
|
||
return None
|
||
return _resample_to(curr_arr, grid), _resample_to(arr[:k], grid)
|
||
|
||
|
||
def _dtw_cost_matrix_scalar(
|
||
x: np.ndarray, y: np.ndarray, n: int, m: int, w: int
|
||
) -> np.ndarray:
|
||
"""Reference (scalar) Sakoe-Chiba DTW cost-matrix fill. Kept verbatim as the
|
||
fallback for :func:`_dtw_cost_banded` so behavior can never regress."""
|
||
cost_matrix = np.full((n + 1, m + 1), float("inf"))
|
||
cost_matrix[0, 0] = 0
|
||
for i in range(1, n + 1):
|
||
center = i * (m / n)
|
||
start_j = max(1, int(center - w))
|
||
end_j = min(m, int(center + w) + 1)
|
||
for j in range(start_j, end_j + 1):
|
||
cost = abs(float(x[i - 1] - y[j - 1]))
|
||
cost_matrix[i, j] = cost + min(
|
||
cost_matrix[i - 1, j], cost_matrix[i, j - 1], cost_matrix[i - 1, j - 1]
|
||
)
|
||
return cost_matrix
|
||
|
||
|
||
class _BandedCostMatrix:
|
||
"""A Sakoe-Chiba DTW cost matrix that stores only its in-band cells.
|
||
|
||
Cells are kept per anti-diagonal ``d = i + j`` (the order the fill computes
|
||
them), each diagonal padded with one ``inf`` cell on either side. ``get``
|
||
returns exactly what the full ``(n+1) x (m+1)`` matrix held at ``(i, j)``:
|
||
the computed value in band, ``inf`` everywhere else (audit LIVE-20).
|
||
"""
|
||
|
||
__slots__ = ("_vals", "_first", "_count", "_start")
|
||
|
||
def __init__(self, vals: np.ndarray, first: list[int], count: list[int], start: list[int]) -> None:
|
||
self._vals = vals
|
||
self._first = first
|
||
self._count = count
|
||
self._start = start
|
||
|
||
def get(self, i: int, j: int) -> float:
|
||
d = i + j
|
||
k = i - self._first[d]
|
||
if 0 <= k < self._count[d]:
|
||
# .item(): a Python float, without boxing the whole array (a .tolist()
|
||
# copy of a 2000-point band costs ~4x the band itself).
|
||
return self._vals.item(self._start[d] + 1 + k)
|
||
return math.inf
|
||
|
||
|
||
def _dtw_cost_banded(
|
||
x: np.ndarray, y: np.ndarray, n: int, m: int, w: int
|
||
) -> _BandedCostMatrix:
|
||
"""Bit-identical vectorized fill of the scalar cost matrix, in-band cells only.
|
||
|
||
The DTW recurrence is sequential, but all cells on one anti-diagonal
|
||
(``i + j`` constant) depend only on the two before it, so each diagonal is
|
||
one vectorized NumPy update. The Sakoe-Chiba band, the per-row bounds
|
||
(``int`` truncation), the ``local + min(min(up, left), diag)`` recurrence and
|
||
the out-of-band ``inf`` cells match the scalar loop exactly, so every cell -
|
||
and the backtracked path - is identical. (#311 follow-up: this fill
|
||
dominates envelope rebuilds.)
|
||
|
||
Until audit LIVE-20 it filled a full ``(n+1) x (m+1)`` matrix and masked all
|
||
of every anti-diagonal down to the band. Both band edges are non-decreasing
|
||
in ``i``, so the in-band cells of a diagonal are one contiguous run of rows,
|
||
found by a binary search; and a run moves by at most one row from one
|
||
diagonal to the next, so one ``inf`` pad cell per side makes every
|
||
predecessor lookup a slice. A 2000 x 2000 pair at the default 20% band
|
||
stores ~2.5x fewer cells and fills ~4-8x faster.
|
||
"""
|
||
xf = np.asarray(x, dtype=float)
|
||
yf = np.asarray(y, dtype=float)
|
||
# Per-row band bounds, identical to the scalar start_j/end_j (int truncates
|
||
# toward zero, matching Python int()).
|
||
i_idx = np.arange(1, n + 1)
|
||
center = i_idx * (m / n)
|
||
lo = np.maximum(1, (center - w).astype(np.int64))
|
||
hi = np.minimum(m, (center + w).astype(np.int64) + 1)
|
||
# Rows in band on diagonal d: i + lo[i] <= d <= i + hi[i] (both sides strictly
|
||
# increasing in i), inside the matrix (1 <= i <= n, 1 <= d - i <= m).
|
||
d_all = np.arange(n + m + 1)
|
||
first = np.searchsorted(i_idx + hi, d_all, side="left") + 1
|
||
last = np.searchsorted(i_idx + lo, d_all, side="right")
|
||
first = np.maximum(first, np.maximum(1, d_all - m))
|
||
last = np.minimum(last, np.minimum(n, d_all - 1))
|
||
count = np.maximum(0, last - first + 1)
|
||
first[0], count[0] = 0, 1 # the origin (0, 0)
|
||
start = np.zeros(n + m + 2, dtype=np.int64)
|
||
np.cumsum(count + 2, out=start[1:])
|
||
vals = np.full(int(start[-1]), np.inf)
|
||
vals[1] = 0.0
|
||
first_l: list[int] = first.tolist()
|
||
count_l: list[int] = count.tolist()
|
||
start_l: list[int] = start.tolist()
|
||
|
||
def run(e: int, p: int, q: int) -> np.ndarray:
|
||
"""Cells of rows ``p..q`` on diagonal ``e`` (``inf`` outside its band)."""
|
||
fe, ce = first_l[e], count_l[e]
|
||
base = start_l[e] + 1 - fe
|
||
if p >= fe - 1 and q <= fe + ce:
|
||
return vals[base + p: base + q + 1]
|
||
out = np.full(q - p + 1, np.inf)
|
||
a, b = max(p, fe), min(q, fe + ce - 1)
|
||
if a <= b:
|
||
out[a - p: b - p + 1] = vals[base + a: base + b + 1]
|
||
return out
|
||
|
||
for d in range(2, n + m + 1):
|
||
c = count_l[d]
|
||
if c <= 0:
|
||
continue
|
||
p = first_l[d]
|
||
q = p + c - 1
|
||
# cells (i, d - i) for i = p..q: x[i - 1] against y[d - i - 1]
|
||
local = np.abs(xf[p - 1: q] - yf[d - q - 1: d - p][::-1])
|
||
best = np.minimum(
|
||
np.minimum(run(d - 1, p - 1, q - 1), run(d - 1, p, q)), # up, left
|
||
run(d - 2, p - 1, q - 1), # diag
|
||
)
|
||
s = start_l[d] + 1
|
||
np.add(local, best, out=vals[s: s + c])
|
||
return _BandedCostMatrix(vals, first_l, count_l, start_l)
|
||
|
||
|
||
def compute_dtw_path(
|
||
x: np.ndarray, y: np.ndarray, band_width_ratio: float = 0.1
|
||
) -> list[tuple[int, int]]:
|
||
"""
|
||
Compute DTW path with Sakoe-Chiba constraint.
|
||
Returns list of (x_index, y_index) tuples mapping X to Y.
|
||
"""
|
||
n, m = len(x), len(y)
|
||
if n == 0 or m == 0:
|
||
return []
|
||
|
||
# Pre-flight memory guard: the scalar fallback's cost matrix is (n+1)x(m+1)
|
||
# float64. An uncapped call from a 1 Hz long cycle can request >1 GB there.
|
||
# If the allocation would exceed ~80 MB, skip DTW and return an empty path so
|
||
# the caller falls back to linear interpolation (graceful degrade rather
|
||
# than OOM-killing Home Assistant, issue #388). The banded fill needs far
|
||
# less, but the cap stays where it was so which pairs get a path is unchanged
|
||
# (``_closed_end_path_exists`` mirrors it).
|
||
_DTW_CELL_BUDGET = 10_000_000 # 10 M cells x 8 B ≈ 80 MB
|
||
if (n + 1) * (m + 1) > _DTW_CELL_BUDGET:
|
||
_LOGGER.warning(
|
||
"DTW cost matrix %dx%d would need %.0f MB — skipping DTW refinement "
|
||
"(cap compute_envelope_worker inputs via MAX_ALIGN_GRID_POINTS to prevent this)",
|
||
n, m, (n + 1) * (m + 1) * 8 / 1e6,
|
||
)
|
||
return []
|
||
|
||
w = max(1, int(min(n, m) * band_width_ratio))
|
||
try:
|
||
cost = _dtw_cost_banded(x, y, n, m, w).get
|
||
except Exception: # pylint: disable=broad-exception-caught
|
||
full = _dtw_cost_matrix_scalar(x, y, n, m, w)
|
||
cost = lambda i, j: full[i, j] # noqa: E731
|
||
|
||
# Backtracking
|
||
if math.isinf(cost(n, m)):
|
||
# Endpoint is unreachable (e.g. Sakoe-Chiba band excluded it); no valid path.
|
||
return []
|
||
|
||
path: list[tuple[int, int]] = []
|
||
i, j = n, m
|
||
|
||
while i > 0 or j > 0:
|
||
# Record current zero-based coordinate before stepping back.
|
||
path.append((max(i - 1, 0), max(j - 1, 0)))
|
||
|
||
if i == 0:
|
||
j -= 1
|
||
elif j == 0:
|
||
i -= 1
|
||
else:
|
||
candidates_cost = [
|
||
(cost(i - 1, j), 0), # deletion (i-1)
|
||
(cost(i, j - 1), 1), # insertion (j-1)
|
||
(cost(i - 1, j - 1), 2) # match (both)
|
||
]
|
||
candidates_cost.sort(key=lambda item: item[0])
|
||
best_move = candidates_cost[0][1]
|
||
if best_move == 0:
|
||
i -= 1
|
||
elif best_move == 1:
|
||
j -= 1
|
||
else:
|
||
i -= 1
|
||
j -= 1
|
||
|
||
path.reverse()
|
||
|
||
return path
|
||
|
||
def compute_envelope_worker(
|
||
raw_cycles_data: list[tuple[list[float], list[float], Optional[float]]] | list[tuple[list[float], list[float]]],
|
||
dtw_bandwidth: float,
|
||
reference_mask: list[bool] | None = None,
|
||
) -> tuple[list[float], list[float], list[float], list[float], list[float], float] | None:
|
||
"""
|
||
Compute statistical envelope.
|
||
Args:
|
||
raw_cycles_data: list of (offsets, power_values, duration) tuples.
|
||
Duration may be None and is used to compute target_duration.
|
||
dtw_bandwidth: ratio.
|
||
reference_mask: optional per-cycle flags (parallel to raw_cycles_data).
|
||
When any entry is True, the robust reference curve is built from the
|
||
median of the flagged cycles only (e.g. user-verified "golden"
|
||
cycles), so trusted cycles define the shape every other cycle is
|
||
warped onto. Min/max/avg/std bands are still built from all cycles.
|
||
Returns:
|
||
(time_grid, min_curve, max_curve, avg_curve, std_curve, target_duration) or None.
|
||
|
||
The bands are the pointwise extremes of the DTW-**warped** members, so
|
||
only a consumer that re-derives the same warp
|
||
(:func:`align_trace_to_envelope`) can compare an observed trace against
|
||
them honestly.
|
||
"""
|
||
if not raw_cycles_data:
|
||
return None
|
||
normalized_curves: list[tuple[np.ndarray, np.ndarray, float]] = []
|
||
golden_flags: list[bool] = []
|
||
sampling_rates: list[float] = []
|
||
|
||
# 1. Pre-process input
|
||
for idx, curve in enumerate(raw_cycles_data):
|
||
# Unpack curve tuple: (offsets, values) or (offsets, values, duration)
|
||
# Backward compatible with 2-tuple (offsets, values) format
|
||
try:
|
||
offsets_list, values_list, *rest = curve
|
||
curve_duration = rest[0] if rest else None
|
||
except (ValueError, TypeError, OverflowError):
|
||
continue
|
||
|
||
if not offsets_list or not values_list:
|
||
continue
|
||
|
||
if len(offsets_list) != len(values_list):
|
||
min_len = min(len(offsets_list), len(values_list))
|
||
if min_len < 3:
|
||
continue
|
||
offsets_list = offsets_list[:min_len]
|
||
values_list = values_list[:min_len]
|
||
|
||
if len(offsets_list) < 3 or len(values_list) < 3:
|
||
continue
|
||
|
||
try:
|
||
offsets = np.asarray(offsets_list, dtype=float)
|
||
values = np.asarray(values_list, dtype=float)
|
||
except (TypeError, ValueError, OverflowError):
|
||
continue
|
||
|
||
# Drop paired entries where either coordinate is non-finite.
|
||
finite_mask = np.isfinite(offsets) & np.isfinite(values)
|
||
offsets = offsets[finite_mask]
|
||
values = values[finite_mask]
|
||
if len(offsets) < 3:
|
||
continue
|
||
|
||
# Stored offsets are rounded to 0.1s, so two readings less than 0.1s apart
|
||
# collapse onto the same offset. A single such duplicate must not discard the
|
||
# whole trace (#377): drop the duplicate sample(s) instead of the cycle. Only
|
||
# exact duplicates are collapsed here; a genuinely out-of-order (decreasing)
|
||
# offset - which sorted storage never produces - is left for the strict check
|
||
# below to reject, exactly as before.
|
||
if offsets.size > 1:
|
||
diffs = np.diff(offsets)
|
||
if np.any(diffs == 0):
|
||
keep = np.concatenate(([True], diffs != 0))
|
||
dropped = int((~keep).sum())
|
||
offsets = offsets[keep]
|
||
values = values[keep]
|
||
_LOGGER.debug(
|
||
"compute_envelope_worker: dropped %d duplicate sample offset(s) "
|
||
"from a cycle trace (0.1s offset rounding)",
|
||
dropped,
|
||
)
|
||
if len(offsets) < 3:
|
||
continue
|
||
|
||
if not np.all(np.diff(offsets) > 0):
|
||
continue
|
||
|
||
try:
|
||
dur = float(curve_duration) if curve_duration is not None else float(offsets[-1])
|
||
except (TypeError, ValueError, OverflowError):
|
||
continue
|
||
|
||
# Validate duration is positive and finite before appending.
|
||
if not (dur > 0 and np.isfinite(dur)):
|
||
continue
|
||
|
||
normalized_curves.append((offsets, values, dur))
|
||
golden_flags.append(bool(reference_mask[idx]) if reference_mask and idx < len(reference_mask) else False)
|
||
|
||
if len(offsets) > 1:
|
||
intervals = np.diff(offsets)
|
||
positive_intervals = intervals[intervals > 0]
|
||
if positive_intervals.size > 0:
|
||
sr = float(np.median(positive_intervals))
|
||
if np.isfinite(sr):
|
||
sampling_rates.append(sr)
|
||
if not normalized_curves:
|
||
return None
|
||
|
||
# 2. Reference Selection
|
||
# The grid is sized from the median duration. Input is (offsets, values, duration).
|
||
max_times = [float(dur) for _, _, dur in normalized_curves]
|
||
median_dur = float(np.median(max_times))
|
||
ref_idx = int(np.argmin([abs(t - median_dur) for t in max_times]))
|
||
|
||
target_duration = max_times[ref_idx]
|
||
avg_sample_rate = float(np.median(sampling_rates)) if sampling_rates else 2.0
|
||
|
||
# Ensure target_duration is valid for calculations
|
||
if not (target_duration > 0 and np.isfinite(target_duration)):
|
||
target_duration = 1.0 # Safe default
|
||
|
||
align_dt = avg_sample_rate
|
||
num_points = max(50, int(target_duration / align_dt))
|
||
if num_points > MAX_ALIGN_GRID_POINTS:
|
||
num_points = MAX_ALIGN_GRID_POINTS
|
||
align_dt = target_duration / num_points # re-derive so per-cycle grids inherit the cap
|
||
time_grid = np.linspace(0.0, target_duration, num_points)
|
||
|
||
# Robust reference curve: the pointwise MEDIAN across all cycles resampled
|
||
# onto the shared grid - a synthetic "medoid" that is not distorted by a
|
||
# single atypical cycle near the median duration and handles multi-mode
|
||
# profiles far better than picking one representative curve. Falls back to
|
||
# the single closest-to-median cycle when there are too few cycles for a
|
||
# stable median.
|
||
golden_indices = [i for i, g in enumerate(golden_flags) if g]
|
||
if golden_indices:
|
||
# Trusted "golden" cycles define the reference shape.
|
||
grid_curves = np.array(
|
||
[
|
||
np.interp(time_grid, normalized_curves[i][0], normalized_curves[i][1])
|
||
for i in golden_indices
|
||
]
|
||
)
|
||
ref_array = np.median(grid_curves, axis=0)
|
||
elif len(normalized_curves) >= 3:
|
||
grid_curves = np.array(
|
||
[np.interp(time_grid, offs, vals) for offs, vals, _ in normalized_curves]
|
||
)
|
||
ref_array = np.median(grid_curves, axis=0)
|
||
else:
|
||
ref_offsets, ref_values, _ = normalized_curves[ref_idx]
|
||
ref_array = np.interp(time_grid, ref_offsets, ref_values)
|
||
|
||
# 3. Resample & DTW: warp every cycle onto the robust reference.
|
||
resampled: list[np.ndarray] = []
|
||
|
||
for offsets, values, dur in normalized_curves:
|
||
this_dur = dur
|
||
# Cap this grid too, not just the reference one: a cycle far longer than the
|
||
# median would otherwise size its own grid past the cap and push the cost
|
||
# matrix over compute_dtw_path's budget, which silently drops the outlier
|
||
# back to plain interpolation. Capping keeps DTW alignment available for it.
|
||
this_num_points = min(MAX_ALIGN_GRID_POINTS, max(10, int(this_dur / align_dt)))
|
||
this_grid = np.linspace(0.0, this_dur, this_num_points)
|
||
this_array = np.interp(this_grid, offsets, values)
|
||
|
||
path = compute_dtw_path(this_array, ref_array, band_width_ratio=dtw_bandwidth)
|
||
|
||
if not path:
|
||
resampled.append(np.interp(time_grid, offsets, values))
|
||
continue
|
||
path_arr = np.array(path)
|
||
cand_indices = path_arr[:, 0]
|
||
ref_indices = path_arr[:, 1]
|
||
|
||
# Interpolate map
|
||
# Map ref indices (time_grid indices) to cand indices (this_grid indices)
|
||
# We assume monotonicity and filter duplicates by taking mean
|
||
|
||
# Simplified: Use numpy interp of indicies
|
||
# ref_indices are 0..N_ref
|
||
# cand_indices are 0..N_cand
|
||
# We need mapping: for ref_idx in 0..num_points, what is cand_idx?
|
||
|
||
# Since ref_indices in path are not strictly increasing (duplicates),
|
||
# we can't use them as 'x' for interp directly if strictness required.
|
||
# But we can sort/unique them.
|
||
|
||
# Sort by ref_index? Path is already sorted roughly.
|
||
# Handle duplicates: average candidate indices for same ref index.
|
||
unique_ref, inverse = np.unique(ref_indices, return_inverse=True)
|
||
# Computing mean candidate index for each unique ref index
|
||
# This is slow in python loop.
|
||
# Vectorized:
|
||
# np.bincount?
|
||
mean_cand_indices = np.zeros_like(unique_ref, dtype=float)
|
||
np.add.at(mean_cand_indices, inverse, cand_indices)
|
||
counts = np.bincount(inverse)
|
||
mean_cand_indices /= counts
|
||
|
||
# Now we have unique_ref -> mean_cand_indices
|
||
# Interpolate to full time_grid (0..num_points-1)
|
||
mapped_cand_indices = np.interp(
|
||
np.arange(num_points),
|
||
unique_ref,
|
||
mean_cand_indices,
|
||
left=0,
|
||
right=len(this_array)-1
|
||
)
|
||
|
||
# Now get values
|
||
mapped_times = mapped_cand_indices * (this_dur / (len(this_array)-1))
|
||
warped_values = np.interp(mapped_times, this_grid, this_array)
|
||
resampled.append(warped_values)
|
||
|
||
# 4. Compute Stats
|
||
stacked = np.vstack(resampled)
|
||
min_curve = np.min(stacked, axis=0)
|
||
max_curve = np.max(stacked, axis=0)
|
||
avg_curve = np.mean(stacked, axis=0)
|
||
std_curve = np.std(stacked, axis=0)
|
||
|
||
return (
|
||
time_grid.tolist(),
|
||
min_curve.tolist(),
|
||
max_curve.tolist(),
|
||
avg_curve.tolist(),
|
||
std_curve.tolist(),
|
||
float(target_duration)
|
||
)
|
||
|
||
|
||
def align_trace_to_envelope(
|
||
t_obs: list[float] | np.ndarray,
|
||
p_obs: list[float] | np.ndarray,
|
||
time_grid: list[float] | np.ndarray,
|
||
reference: list[float] | np.ndarray | None,
|
||
dtw_bandwidth: float,
|
||
) -> tuple[np.ndarray, bool]:
|
||
"""Map each observed sample's time onto an envelope's own time grid.
|
||
|
||
Uses the **same DTW warp** :func:`compute_envelope_worker` used to build that
|
||
envelope's bands, only reduced in the opposite direction (candidate index ->
|
||
mean reference index instead of reference index -> mean candidate index).
|
||
That matters: ``min``/``max`` are the pointwise extremes of the *warped*
|
||
member curves, so a consumer that re-derives the same warp sees every member
|
||
inside the band by construction, while a consumer that stretches time
|
||
proportionally does not.
|
||
|
||
A proportional stretch is not a cheap approximation of this, it is a worse
|
||
alignment: real programmes absorb their run-to-run duration variance in one
|
||
stretch of the cycle (a dishwasher's drying tail), not uniformly, so scaling
|
||
the whole axis *moves* the fixed-time features. Measured on the maintainer's
|
||
corpus (register item 324) the final heating block's placement error grew
|
||
from sd 2.5 min (no scaling at all) to sd 3.7 min under proportional scaling,
|
||
which is what fabricated the out-of-band ``spike``/``dip`` artifacts.
|
||
|
||
Args:
|
||
t_obs: observed sample offsets (seconds from cycle start, increasing).
|
||
p_obs: observed power values, parallel to ``t_obs``.
|
||
time_grid: the envelope's time grid.
|
||
reference: the curve to warp onto. Callers pass ``envelope["avg"]``; see
|
||
``ProfileStore._align_to_envelope`` for why that beats the build's
|
||
own pre-DTW reference. ``None`` forces the proportional stretch.
|
||
dtw_bandwidth: Sakoe-Chiba band ratio, same value the build used.
|
||
|
||
**Known limit, measured and accepted (register item 347).** The observed span
|
||
here is ``t_obs[-1]``, while ``_rebuild_envelope_sync`` builds each member
|
||
over ``manual_duration`` (when plausible for its trace, audit MATCH-EVAL-17)
|
||
or ``max(last_offset, stored_duration)``. Where the
|
||
stored duration runs past the last sample the build covered a slightly longer
|
||
span than this re-derivation does, so the warp is not bit-for-bit the build's
|
||
own. Measured over the 703 stored cycles in the maintainer's corpus it
|
||
affects **7 of them (1.0%)**, median gap 8.2 s and worst 201 s on a 14040 s
|
||
cycle (1.4%) - all well inside one grid step. Closing it means threading the
|
||
authoritative duration through ``compute_envelope_conformance`` and
|
||
``detect_cycle_artifacts``, whose signatures take points and no cycle, and
|
||
through their manager / ws_api / playground callers. Not done: that is a wide
|
||
change to alignment code whose current behaviour is measured (item 324), for
|
||
a sub-grid-step effect on 1% of cycles.
|
||
|
||
Returns:
|
||
``(envelope-space time per observed sample, used_dtw)``. ``used_dtw`` is
|
||
False when the warp was unavailable and the proportional stretch was
|
||
used instead, so callers can loosen any judgement they base on it.
|
||
Never raises: every failure degrades to the proportional stretch.
|
||
"""
|
||
t_arr = np.asarray(t_obs, dtype=float)
|
||
p_arr = np.asarray(p_obs, dtype=float)
|
||
tg = np.asarray(time_grid, dtype=float)
|
||
|
||
def _proportional() -> np.ndarray:
|
||
if tg.size < 2 or t_arr.size < 1:
|
||
return t_arr
|
||
obs_dur = float(t_arr[-1])
|
||
env_dur = float(tg[-1])
|
||
if obs_dur <= 0 or env_dur <= 0:
|
||
return np.clip(t_arr, tg[0], tg[-1])
|
||
return np.clip(t_arr * (env_dur / obs_dur), tg[0], tg[-1])
|
||
|
||
try:
|
||
ref = np.asarray(reference, dtype=float) if reference is not None else None
|
||
if (
|
||
ref is None
|
||
or ref.size != tg.size
|
||
or tg.size < 2
|
||
or t_arr.size < 2
|
||
or t_arr.size != p_arr.size
|
||
or dtw_bandwidth <= 0
|
||
):
|
||
return _proportional(), False
|
||
|
||
obs_dur = float(t_arr[-1])
|
||
env_dur = float(tg[-1])
|
||
if not (obs_dur > 0 and env_dur > 0):
|
||
return _proportional(), False
|
||
|
||
# Rebuild the per-cycle grid the same way the envelope build did: one
|
||
# point per grid step, so an index offset is a time offset on both axes
|
||
# and the Sakoe-Chiba band means the same span of minutes either side.
|
||
grid_dt = env_dur / (tg.size - 1)
|
||
if grid_dt <= 0:
|
||
return _proportional(), False
|
||
n_obs = int(min(MAX_ALIGN_GRID_POINTS, max(10, int(obs_dur / grid_dt))))
|
||
obs_grid = np.linspace(0.0, obs_dur, n_obs)
|
||
obs_array = np.interp(obs_grid, t_arr, p_arr)
|
||
|
||
path = compute_dtw_path(obs_array, ref, band_width_ratio=dtw_bandwidth)
|
||
if not path:
|
||
return _proportional(), False
|
||
|
||
path_arr = np.array(path)
|
||
obs_indices = path_arr[:, 0]
|
||
ref_indices = path_arr[:, 1]
|
||
# Average the reference indices a single observed index maps onto (DTW
|
||
# paths repeat indices wherever one axis is stretched).
|
||
unique_obs, inverse = np.unique(obs_indices, return_inverse=True)
|
||
mean_ref = np.zeros_like(unique_obs, dtype=float)
|
||
np.add.at(mean_ref, inverse, ref_indices)
|
||
mean_ref /= np.bincount(inverse)
|
||
|
||
env_idx = np.interp(
|
||
np.arange(n_obs), unique_obs, mean_ref, left=0, right=tg.size - 1
|
||
)
|
||
env_time_on_grid = env_idx * grid_dt
|
||
mapped = np.interp(t_arr, obs_grid, env_time_on_grid)
|
||
return np.clip(mapped, tg[0], tg[-1]), True
|
||
except Exception: # pylint: disable=broad-exception-caught
|
||
# Alignment is a comparison aid, never a correctness gate: degrade.
|
||
return _proportional(), False
|
||
|
||
|
||
_DTW_PATH_CELL_BUDGET = 10_000_000 # compute_dtw_path's matrix cap
|
||
|
||
|
||
def _closed_end_path_exists(n: int, m: int, band_width_ratio: float) -> bool:
|
||
"""Whether ``compute_dtw_path`` would return a path for an ``n x m`` pair.
|
||
|
||
Exactly its three failure cases without filling the matrix: an empty side,
|
||
the cell budget, and an end cell the Sakoe-Chiba band cannot reach. Each row's
|
||
band is ``[lo, hi]`` (the fill's own truncation); cells are reachable left to
|
||
right from the diagonal, so the end is reachable iff no row's band starts more
|
||
than one column past where the previous row's ends (row 0 is the origin).
|
||
"""
|
||
if n == 0 or m == 0 or (n + 1) * (m + 1) > _DTW_PATH_CELL_BUDGET:
|
||
return False
|
||
w = max(1, int(min(n, m) * band_width_ratio))
|
||
center = np.arange(1, n + 1) * (m / n)
|
||
lo = np.maximum(1, (center - w).astype(np.int64))
|
||
hi = np.minimum(m, (center + w).astype(np.int64) + 1)
|
||
prev_hi = np.concatenate(([0], hi[:-1]))
|
||
return bool(np.all(lo <= prev_hi + 1) and np.all(lo <= hi) and hi[-1] == m)
|
||
|
||
|
||
def verify_profile_alignment_worker(
|
||
current_power: list[float],
|
||
envelope_avg_curve: list[float],
|
||
envelope_time_grid: list[float],
|
||
dtw_bandwidth: float
|
||
) -> tuple[float, float, float]:
|
||
"""
|
||
Verify alignment of current trace against profile envelope.
|
||
Returns: (mapped_envelope_time, mapped_envelope_power, overlap_score)
|
||
"""
|
||
if not current_power or not envelope_avg_curve:
|
||
return 0.0, 9999.0, 0.0
|
||
|
||
curr = np.array(current_power)
|
||
ref = np.array(envelope_avg_curve)
|
||
|
||
# 1. Coarse Alignment
|
||
score, _, offset = find_best_alignment(curr, ref, 1.0)
|
||
|
||
# 2. Extract aligned segments
|
||
# Determine the mapping window.
|
||
|
||
# Symmetric context window: pad equally left and right of the coarse alignment.
|
||
half = ALIGNMENT_CONTEXT_BUFFER // 2
|
||
start_ref = max(0, offset - half)
|
||
end_ref = min(len(ref), offset + len(curr) + half)
|
||
|
||
if end_ref <= start_ref:
|
||
return 0.0, 9999.0, 0.0
|
||
|
||
ref_seg = ref[start_ref:end_ref]
|
||
curr_seg = curr
|
||
|
||
if offset < 0:
|
||
curr_seg = curr[-offset:]
|
||
|
||
# The DTW that used to run here was closed-end: its backtrack starts at the last
|
||
# cell, so the mapped index was ALWAYS the window's last index whenever a path
|
||
# existed - 34x the CPU and a 33-42 MiB matrix per tick for a constant (audit
|
||
# LIVE-03). Same answer, same fallback, no matrix: a path exists unless the
|
||
# window is empty, over the old cell budget, or the band cannot reach the end
|
||
# cell. The resulting position runs ALIGNMENT_CONTEXT_BUFFER // 2 steps ahead
|
||
# of the trace; measured harmless-to-helpful, so it is kept, not "fixed".
|
||
if _closed_end_path_exists(len(curr_seg), len(ref_seg), dtw_bandwidth):
|
||
mapped_idx = end_ref - 1
|
||
else:
|
||
# Fallback to linear mapping based on offset
|
||
mapped_idx = min(len(ref)-1, offset + len(curr) - 1)
|
||
mapped_idx = max(0, mapped_idx)
|
||
|
||
# Ensure sequences are non-empty before indexing
|
||
if not envelope_time_grid or len(ref) == 0:
|
||
return 0.0, 9999.0, 0.0
|
||
mapped_idx = min(mapped_idx, len(envelope_time_grid) - 1, len(ref) - 1)
|
||
|
||
mapped_time = float(envelope_time_grid[mapped_idx])
|
||
mapped_power = float(ref[mapped_idx])
|
||
|
||
return mapped_time, mapped_power, float(score)
|