francop/Home-Assistant
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity
Files
0c420e445f2f93d4223574dba4b366014030cf96
Home-Assistant/custom_components/ha_washdata/ml/feature_extraction.py
T
francop 28a8cb98f6 Added Alexa Music
2026-07-17 10:12:15 -04:00

808 lines
32 KiB
Python
Raw Blame History

# 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/>.
"""NumPy-only runtime feature extraction for the embedded ML models.
This is the bridge between live cycle data and the embedded models: it computes
the exact ``FEATURE_COLUMNS`` each model expects, ported faithfully from the
``ml_washdata`` lab feature definitions (see each ``*_feature_contract.json``).
Three feature extractors are implemented:
- **Cycle-end detector** (``END_FEATURE_COLUMNS``): self-contained from a live
power series + profile expectation; call ``latest_end_event_features``.
- **Live-match commit confidence** (``LIVE_MATCH_FEATURE_COLUMNS``): requires the
match ranking from ``ProfileStore`` plus the observed prefix; call
``live_match_features``.
- **Hybrid cycle quality** (``QUALITY_FEATURE_COLUMNS``): requires the complete
cycle power trace plus profile/match context; call ``quality_features``.
All inputs are plain Python/NumPy (offset-seconds, watts), so this module has no
Home Assistant dependency and is unit-tested directly. It is only invoked when
the user opts into experimental ML models (see engine.py).
"""
from __future__ import annotations
import math
from typing import Sequence
import numpy as np
# ---------------------------------------------------------------------------
# Cycle-end detector
# ---------------------------------------------------------------------------
# Mirrors ml_washdata/wash_ml/end_detection.py. The test suite asserts this list
# equals the embedded model's FEATURE_COLUMNS so the two cannot drift.
END_FEATURE_COLUMNS = [
"elapsed_fraction",
"energy_fraction",
"energy_remaining_expected",
"power_before_ratio",
"drop_ratio",
"peak_seen_ratio",
"low_run_s_log",
"elapsed_log",
]
MIN_LOW_RUN_S = 45.0
Point = tuple[float, float]
def cumulative_energy_wh(points: Sequence[Point]) -> np.ndarray:
"""Trapezoidal cumulative energy (Wh) aligned to each reading, with gap handling.
Segments spanning sensor-outage gaps (larger than ``energy_gap_threshold_s``)
are zeroed out so energy does not inflate across outages, matching the behaviour
of ``signal_processing.integrate_wh`` used for stored ``energy_wh`` fields.
"""
from ..signal_processing import energy_gap_threshold_s # noqa: PLC0415
offsets = np.asarray([float(offset) for offset, _power in points], dtype=float)
powers = np.asarray([max(0.0, float(power)) for _offset, power in points], dtype=float)
if offsets.size < 2:
return np.zeros(offsets.size, dtype=float)
max_gap = energy_gap_threshold_s(offsets)
deltas = np.diff(offsets)
segment = (powers[:-1] + powers[1:]) / 2.0 * deltas / 3600.0
# Zero out segments that span a sensor outage gap to match integrate_wh behaviour.
segment[deltas > max_gap] = 0.0
return np.concatenate([[0.0], np.cumsum(segment)])
def profile_expectation(cycles_points: Sequence[Sequence[Point]]) -> dict[str, float] | None:
"""Median expected duration (s), energy (Wh) and peak (W) over a profile's cycles."""
durations: list[float] = []
energies: list[float] = []
peaks: list[float] = []
for points in cycles_points:
if len(points) < 2:
continue
offsets = [float(offset) for offset, _power in points]
duration = offsets[-1] - offsets[0]
if duration <= 0:
continue
durations.append(duration)
energies.append(float(cumulative_energy_wh(points)[-1]))
peaks.append(max((float(power) for _offset, power in points), default=0.0))
if not durations:
return None
return {
"duration": float(np.median(durations)),
"energy": float(np.median(energies)),
"peak": float(np.median(peaks)),
}
def profile_expectations(cycles: list[dict]) -> dict[str, dict[str, float]]:
"""Median duration/energy/peak per profile from stored cycle dicts.
The dict-based counterpart of :func:`profile_expectation` (which works from
decompressed traces): reads the ``duration``/``energy_wh``/``max_power``
scalar fields already stored on each cycle. Shared by on-device training
(``training_task``) and the ML suggestion engine so the "profile expectation"
definition lives in one place. Profiles with no usable duration are skipped;
missing energy/peak default to 500.
"""
stats: dict[str, dict[str, list[float]]] = {}
for c in cycles:
name = c.get("profile_name")
if not isinstance(name, str) or not name:
continue
s = stats.setdefault(name, {"d": [], "e": [], "p": []})
for key, field in (("d", "duration"), ("e", "energy_wh"), ("p", "max_power")):
v = c.get(field)
if isinstance(v, (int, float)) and not isinstance(v, bool):
s[key].append(float(v))
out: dict[str, dict[str, float]] = {}
for name, s in stats.items():
if not s["d"]:
continue
out[name] = {
"duration": float(np.median(s["d"])),
"energy": float(np.median(s["e"])) if s["e"] else 500.0,
"peak": float(np.median(s["p"])) if s["p"] else 500.0,
}
return out
def latest_end_event_features(
points: Sequence[Point],
expectation: dict[str, float],
*,
min_low_run_s: float = MIN_LOW_RUN_S,
) -> dict[str, float] | None:
"""Features for the most recent low-power run (the "is this the end?" moment).
Returns ``None`` when there is no qualifying low-power run yet (the cycle is
still clearly active), in which case the caller keeps the existing behavior.
Mirrors ``end_detection._cycle_events`` for a single (latest) event.
"""
if len(points) < 4 or not expectation:
return None
offsets = np.asarray([float(offset) for offset, _power in points], dtype=float)
powers = np.asarray([max(0.0, float(power)) for _offset, power in points], dtype=float)
start = float(offsets[0])
peak = max(float(expectation.get("peak") or 0.0), float(np.max(powers)), 1.0)
low_threshold = max(5.0, 0.02 * peak)
profile_duration = max(float(expectation.get("duration") or 0.0), 1.0)
profile_energy = max(float(expectation.get("energy") or 0.0), 1e-6)
cumulative = cumulative_energy_wh(points)
# Find the most recent contiguous low-power run.
count = len(points)
run_start: int | None = None
index = count - 1
while index >= 0 and powers[index] < low_threshold:
run_start = index
index -= 1
if run_start is None:
return None
run_end = count - 1
run_duration = float(offsets[run_end] - offsets[run_start])
if run_duration < min_low_run_s:
return None
event_time = float(offsets[run_start] - start)
elapsed = max(event_time, 1.0)
energy_so_far = float(cumulative[run_start])
window = powers[max(0, run_start - 4):run_start]
power_before = float(np.mean(window)) if window.size else float(powers[run_start])
running_peak = float(np.max(powers[: run_start + 1]))
return {
"elapsed_fraction": float(min(elapsed / profile_duration, 2.0)),
"energy_fraction": float(min(energy_so_far / profile_energy, 2.0)),
"energy_remaining_expected": float(max(0.0, 1.0 - energy_so_far / profile_energy)),
"power_before_ratio": float(min(power_before / peak, 2.0)),
"drop_ratio": float(np.clip((power_before - float(powers[run_start])) / peak, 0.0, 1.0)),
"peak_seen_ratio": float(min(running_peak / peak, 2.0)),
"low_run_s_log": float(math.log1p(max(0.0, run_duration))),
"elapsed_log": float(math.log1p(elapsed)),
}
# ---------------------------------------------------------------------------
# Live-match commit confidence
# ---------------------------------------------------------------------------
# Mirrors ml_washdata/wash_ml/live_matching.py COMMIT_FEATURE_COLUMNS.
# The test suite asserts this equals the embedded model's FEATURE_COLUMNS.
LIVE_MATCH_FEATURE_COLUMNS = [
"match_progress_top1",
"top1_distance",
"margin",
"distance_ratio",
"candidate_count_log",
"prefix_active_fraction",
"duration_ratio_top1",
"elapsed_log",
]
def live_match_features(
points: Sequence[Point],
elapsed_s: float,
top1_distance: float,
top2_distance: float | None,
top1_median_duration_s: float,
candidate_count: int,
) -> dict[str, float]:
"""Features for the live-match commit-confidence model.
Args:
points: Observed power readings (offset_s, watts) for the current prefix.
elapsed_s: Seconds elapsed since cycle start.
top1_distance: Blended RMSE+DTW shape distance to the top-1 candidate
prefix (as returned by the profile matcher).
top2_distance: Distance to the top-2 candidate; pass ``None`` or ``0.0``
when only one candidate is available (margin defaults to 1.0).
top1_median_duration_s: Expected (median) duration of the top-1 candidate
profile in seconds.
candidate_count: Number of candidate profiles on this device.
Returns a dict with exactly ``LIVE_MATCH_FEATURE_COLUMNS`` keys.
"""
elapsed = max(0.0, float(elapsed_s))
top1 = max(0.0, float(top1_distance))
top2_raw = float(top2_distance) if top2_distance is not None else 0.0
top2 = top2_raw if top2_raw > 1e-9 else top1 + 1.0
margin = max(0.0, top2 - top1)
dur = float(top1_median_duration_s)
progress = (elapsed / dur) if dur > 0 else 1.0
# prefix_active_fraction: fraction of prefix readings clearly above idle.
# Lab uses > 0.05 on a peak-normalised trace; equivalent here is > 5% of
# peak, with a 1 W floor so a cold trace never divides by near-zero.
if points:
powers = np.asarray([max(0.0, float(p)) for _, p in points], dtype=float)
peak = float(np.max(powers)) if powers.size else 0.0
active_thr = max(1.0, 0.05 * peak)
prefix_active_fraction = float(np.mean(powers > active_thr)) if powers.size else 0.0
else:
prefix_active_fraction = 0.0
return {
"match_progress_top1": float(min(progress, 2.0)),
"top1_distance": float(top1),
"margin": float(margin),
"distance_ratio": float(top1 / top2) if top2 > 1e-9 else 1.0,
"candidate_count_log": float(math.log1p(max(0, int(candidate_count)))),
"prefix_active_fraction": float(prefix_active_fraction),
"duration_ratio_top1": float(min(progress, 2.0)),
"elapsed_log": float(math.log1p(elapsed)),
}
# ---------------------------------------------------------------------------
# Remaining-time / progress regressor
# ---------------------------------------------------------------------------
# Feature columns for the on-device remaining-time regressor. Unlike the three
# classifier heads this model is a ``standardized_linear`` regressor whose target
# is the cycle completion fraction (elapsed / total_actual). There is no shipped
# baseline: the model exists only once on-device training promotes one over the
# naive elapsed/expected estimate (``elapsed_over_expected`` is deliberately the
# first column so the naive baseline is trivially recoverable). The same
# extractor runs at training time on synthesized prefixes and at inference on the
# live trace, so the columns cannot drift.
PROGRESS_FEATURE_COLUMNS = [
"elapsed_over_expected",
"energy_over_expected",
"mean_power_over_peak",
"recent_power_over_peak",
"tail_slope_norm",
"active_fraction",
"elapsed_log",
]
def progress_features(
points: Sequence[Point],
expectation: dict[str, float],
) -> dict[str, float] | None:
"""Features for the remaining-time regressor from a running-cycle prefix.
Args:
points: Observed prefix power readings (offset_s, watts).
expectation: Matched profile's median ``duration``/``energy``/``peak``
(as produced by :func:`profile_expectation`).
Returns a dict with exactly ``PROGRESS_FEATURE_COLUMNS`` keys, or ``None``
when there is too little data to characterise progress.
"""
pts = _clean_points(points)
if len(pts) < 4 or not expectation:
return None
offsets = np.asarray([o for o, _ in pts], dtype=float)
powers = np.asarray([p for _, p in pts], dtype=float)
elapsed = max(float(offsets[-1] - offsets[0]), 1.0)
exp_dur = max(float(expectation.get("duration") or 0.0), 1.0)
exp_energy = max(float(expectation.get("energy") or 0.0), 1e-6)
exp_peak = max(float(expectation.get("peak") or 0.0), 1.0)
energy_so_far = float(cumulative_energy_wh(pts)[-1])
active_thr = max(1.0, 0.05 * exp_peak)
active_mask = powers > active_thr
active = powers[active_mask]
mean_power = float(np.mean(active)) if active.size else 0.0
# Recent power: mean of the trailing ~5% of samples (min one sample).
tail_n = max(1, len(pts) // 20)
recent_power = float(np.mean(powers[-tail_n:]))
# Tail slope over the last quarter (W per sample), normalised by peak: a
# declining tail is a strong "near the end" signal.
quarter = max(2, len(pts) // 4)
tail = powers[-quarter:]
if tail.size >= 2:
x = np.arange(tail.size, dtype=float)
xm = x - float(np.mean(x))
denom = float(np.dot(xm, xm))
slope = float(np.dot(xm, tail - float(np.mean(tail))) / denom) if denom > 1e-9 else 0.0
else:
slope = 0.0
return {
"elapsed_over_expected": float(min(elapsed / exp_dur, 3.0)),
"energy_over_expected": float(min(energy_so_far / exp_energy, 3.0)),
"mean_power_over_peak": float(min(mean_power / exp_peak, 2.0)),
"recent_power_over_peak": float(min(recent_power / exp_peak, 2.0)),
"tail_slope_norm": float(np.clip(slope / exp_peak, -2.0, 2.0)),
"active_fraction": float(np.mean(active_mask)) if powers.size else 0.0,
"elapsed_log": float(math.log1p(elapsed)),
}
# ---------------------------------------------------------------------------
# Hybrid cycle quality
# ---------------------------------------------------------------------------
# Mirrors ml_washdata/wash_ml/hybrid_curve_quality.py HYBRID_FEATURE_COLUMNS.
# Order must match the embedded model exactly; the test suite asserts this.
QUALITY_FEATURE_COLUMNS = [
# profile / context
"duration_log_ratio",
"energy_log_ratio",
"peak_log_ratio",
"profile_distance",
"label_margin_positive",
"max_gap_ratio",
"low_power_gap_ratio",
"false_end_energy_ratio",
"sample_density_log",
"peak_density_log",
"local_spike_score",
"local_spike_rate",
"local_noise_score",
"leading_idle_ratio",
"trailing_idle_ratio",
"trimmed_duration_log_ratio",
"flag_pressure",
"shape_fit_penalty",
# trace shape
"shape_active_fraction",
"shape_early_energy_fraction",
"shape_late_energy_fraction",
"shape_mid_trough_depth",
"shape_peak_density",
"shape_max_step_drop",
"shape_max_step_rise",
"shape_active_cv",
"shape_autocorr_lag1",
"shape_derivative_sign_changes",
"shape_plateau_ratio",
"shape_tail_slope",
# availability
"has_trace",
]
_QUALITY_TRACE_LENGTH = 128
_IDLE_THRESHOLD_W = 2.0
_STOP_THRESHOLD_W = 2.0
_SHAPE_COLUMNS = [
"shape_active_fraction",
"shape_early_energy_fraction",
"shape_late_energy_fraction",
"shape_mid_trough_depth",
"shape_peak_density",
"shape_max_step_drop",
"shape_max_step_rise",
"shape_active_cv",
"shape_autocorr_lag1",
"shape_derivative_sign_changes",
"shape_plateau_ratio",
"shape_tail_slope",
]
def quality_features(
points: Sequence[Point],
profile_median_duration_s: float,
profile_median_energy_wh: float,
profile_median_peak_w: float,
profile_distance: float,
label_margin: float,
profile_fit_score: float,
flag_count: int,
*,
trace_length: int = _QUALITY_TRACE_LENGTH,
) -> dict[str, float]:
"""Features for the hybrid curve-quality model (problem/bad-cycle detector).
Args:
points: Complete cycle power trace (offset_s, watts).
profile_median_duration_s: Median duration of the matched profile (s).
profile_median_energy_wh: Median energy of the matched profile (Wh).
profile_median_peak_w: Median peak power of the matched profile (W).
profile_distance: Shape distance from the MatchResult to the assigned
profile envelope (higher = worse fit).
label_margin: Score margin between the top-1 and top-2 profile candidates
(positive = confident match; 0.0 when only one candidate exists).
profile_fit_score: Profile fit score in [0, 1] from the matcher.
flag_count: Number of detection/anomaly flags raised for this cycle by
the existing detector (early_power_dip, false_end_pause_seen, etc.).
Returns a dict with exactly ``QUALITY_FEATURE_COLUMNS`` keys.
"""
pts = _clean_points(points)
if len(pts) < 4:
return _no_trace_quality_features(
profile_distance=profile_distance,
label_margin=label_margin,
profile_fit_score=profile_fit_score,
flag_count=flag_count,
)
offsets = np.asarray([float(o) for o, _ in pts], dtype=float)
powers = np.asarray([float(p) for _, p in pts], dtype=float)
duration_s = max(float(offsets[-1] - offsets[0]), 1.0)
total_energy_wh = float(cumulative_energy_wh(pts)[-1])
max_power_w = float(np.max(powers))
# -- profile context ratios --
prof_dur = max(float(profile_median_duration_s), 1.0)
prof_energy = max(float(profile_median_energy_wh), 1e-6)
prof_peak = max(float(profile_median_peak_w), 1.0)
# -- sampling gap features --
intervals = np.diff(offsets)
usable = intervals[(intervals > 0) & (intervals < 3600)]
max_gap_s = float(np.max(usable)) if usable.size else 0.0
# -- low-power / false-end features --
low_gap_s = _longest_low_power_gap_s(pts, _STOP_THRESHOLD_W)
fe_energy_wh = _false_end_energy_wh(pts, _STOP_THRESHOLD_W)
# -- density features --
sample_count = len(pts)
peak_count = _power_peak_count_arr(powers)
# -- noise features (from raw points) --
noise = _trace_noise_features(pts)
# -- idle-padding features --
padding = _trace_padding_ratios(pts, offsets, powers, duration_s, _IDLE_THRESHOLD_W)
# -- trace shape descriptors (from resampled + trimmed trace) --
trace = _resample_to_length(pts, trace_length)
shape = _trace_shape_descriptors(trace) if trace is not None else {c: 0.0 for c in _SHAPE_COLUMNS}
trimmed_ratio = float(padding["trimmed_duration_ratio"])
trimmed_log = math.log(max(1e-6, trimmed_ratio)) # always <= 0
return {
"duration_log_ratio": _log_ratio(duration_s / prof_dur),
"energy_log_ratio": _log_ratio(total_energy_wh / prof_energy),
"peak_log_ratio": _log_ratio(max_power_w / prof_peak),
"profile_distance": float(profile_distance),
"label_margin_positive": float(max(0.0, float(label_margin))),
"max_gap_ratio": _safe_div(max_gap_s, duration_s),
"low_power_gap_ratio": _safe_div(low_gap_s, duration_s),
"false_end_energy_ratio": _safe_div(fe_energy_wh, max(total_energy_wh, 1e-6)),
"sample_density_log": math.log1p(_safe_div(sample_count * 60.0, duration_s)),
"peak_density_log": math.log1p(_safe_div(peak_count * 3600.0, duration_s)),
"local_spike_score": float(noise["local_spike_score"]),
"local_spike_rate": float(noise["local_spike_rate"]),
"local_noise_score": float(noise["local_noise_score"]),
"leading_idle_ratio": float(padding["leading_idle_ratio"]),
"trailing_idle_ratio": float(padding["trailing_idle_ratio"]),
"trimmed_duration_log_ratio": float(trimmed_log),
"flag_pressure": float(max(0, int(flag_count))),
"shape_fit_penalty": float(max(0.0, 1.0 - float(profile_fit_score))),
**shape,
"has_trace": 1.0,
}
def _no_trace_quality_features(
*,
profile_distance: float,
label_margin: float,
profile_fit_score: float,
flag_count: int,
) -> dict[str, float]:
"""Zero-valued quality features for cycles with no usable power trace."""
return {
"duration_log_ratio": 0.0,
"energy_log_ratio": 0.0,
"peak_log_ratio": 0.0,
"profile_distance": float(profile_distance),
"label_margin_positive": float(max(0.0, float(label_margin))),
"max_gap_ratio": 0.0,
"low_power_gap_ratio": 0.0,
"false_end_energy_ratio": 0.0,
"sample_density_log": 0.0,
"peak_density_log": 0.0,
"local_spike_score": 0.0,
"local_spike_rate": 0.0,
"local_noise_score": 0.0,
"leading_idle_ratio": 0.0,
"trailing_idle_ratio": 0.0,
"trimmed_duration_log_ratio": 0.0,
"flag_pressure": float(max(0, int(flag_count))),
"shape_fit_penalty": float(max(0.0, 1.0 - float(profile_fit_score))),
**{c: 0.0 for c in _SHAPE_COLUMNS},
"has_trace": 0.0,
}
# ---------------------------------------------------------------------------
# Shared helpers (ported from ml_washdata/wash_ml/features.py and
# ml_washdata/wash_ml/hybrid_curve_quality.py - NumPy only)
# ---------------------------------------------------------------------------
def _clean_points(points: Sequence[Point]) -> list[Point]:
"""Filter out non-finite readings, sort by offset, and deduplicate."""
clean: list[Point] = []
for offset_raw, power_raw in points:
offset = float(offset_raw)
power = float(power_raw)
if math.isfinite(offset) and math.isfinite(power):
clean.append((offset, max(0.0, power)))
clean.sort(key=lambda pt: pt[0])
deduped: list[Point] = []
for offset, power in clean:
if deduped and offset == deduped[-1][0]:
deduped[-1] = (offset, power)
else:
deduped.append((offset, power))
return deduped
def _longest_low_power_gap_s(points: list[Point], threshold_w: float) -> float:
"""Longest contiguous span below ``threshold_w`` (seconds)."""
longest = 0.0
current = 0.0
for i in range(1, len(points)):
prev_t, prev_p = points[i - 1]
curr_t, curr_p = points[i]
dt = curr_t - prev_t
if dt <= 0 or dt > 3600:
current = 0.0
continue
avg = (prev_p + curr_p) / 2.0
if avg < threshold_w:
current += dt
longest = max(longest, current)
else:
current = 0.0
return float(longest)
def _false_end_energy_wh(points: list[Point], threshold_w: float) -> float:
"""Energy accumulated during low-power pauses that were followed by more power.
These "false ends" indicate the cycle was interrupted but resumed. Returns
the maximum such pause energy (Wh); 0.0 if no false end occurred.
"""
false_energies: list[float] = []
in_pause = False
pause_energy = 0.0
for i in range(1, len(points)):
prev_t, prev_p = points[i - 1]
curr_t, curr_p = points[i]
dt = curr_t - prev_t
if dt <= 0 or dt > 3600:
in_pause = False
pause_energy = 0.0
continue
avg = (prev_p + curr_p) / 2.0
if avg < threshold_w:
in_pause = True
pause_energy += avg * (dt / 3600.0)
elif in_pause:
false_energies.append(pause_energy)
in_pause = False
pause_energy = 0.0
return float(max(false_energies) if false_energies else 0.0)
def _power_peak_count_arr(powers: np.ndarray) -> int:
"""Number of above-p75 power peaks (rising transitions through the p75 threshold)."""
if powers.size < 3:
return 0
threshold = max(float(np.percentile(powers, 75)), 10.0)
above = powers > threshold
transitions = np.diff(above.astype(int))
return int(np.sum(transitions == 1) + (1 if above[0] else 0))
def _trace_noise_features(points: list[Point]) -> dict[str, float]:
"""Local spike and noise floor metrics over the raw power trace.
Ported from ml_washdata/wash_ml/features.py ``trace_noise_features``.
Distinguishes narrow single-sample spikes from broad appliance phases.
"""
if len(points) < 5:
return {"local_spike_score": 0.0, "local_spike_rate": 0.0, "local_noise_score": 0.0}
powers = np.asarray([float(p) for _, p in points], dtype=float)
active = powers[powers > 0.5]
scale = float(np.percentile(active, 95)) if active.size else float(np.max(powers))
if not math.isfinite(scale) or scale <= 1e-6:
return {"local_spike_score": 0.0, "local_spike_rate": 0.0, "local_noise_score": 0.0}
normalized = np.clip(powers / scale, 0.0, 8.0)
spike_scores: list[float] = []
residuals: list[float] = []
n = normalized.size
for i, value in enumerate(normalized):
left = max(0, i - 2)
right = min(n, i + 3)
neighbors = np.concatenate([normalized[left:i], normalized[i + 1:right]])
if neighbors.size < 2:
continue
local_median = float(np.median(neighbors))
residual = abs(float(value - local_median))
residuals.append(residual)
left_nbr = float(normalized[i - 1]) if i > 0 else local_median
right_nbr = float(normalized[i + 1]) if i + 1 < n else local_median
shoulder = max(left_nbr, right_nbr)
narrow_jump = float(value - shoulder)
if value > 0.15 and (value - local_median) > 0.28 and narrow_jump > 0.18:
spike_scores.append(min(3.0, max(float(value - local_median), narrow_jump)))
spike_count = len(spike_scores)
return {
"local_spike_score": round(float(max(spike_scores, default=0.0)), 6),
"local_spike_rate": round(float(spike_count / max(1, n)), 6),
"local_noise_score": round(float(np.percentile(np.asarray(residuals, dtype=float), 95)) if residuals else 0.0, 6),
}
def _trace_padding_ratios(
pts: list[Point],
offsets: np.ndarray,
powers: np.ndarray,
duration_s: float,
idle_threshold_w: float,
) -> dict[str, float]:
"""Leading/trailing idle fractions and trimmed-duration ratio."""
active_mask = powers > idle_threshold_w
active_indexes = np.where(active_mask)[0]
if active_indexes.size == 0 or duration_s <= 0:
return {"leading_idle_ratio": 0.0, "trailing_idle_ratio": 0.0, "trimmed_duration_ratio": 1.0}
first_active_t = float(offsets[active_indexes[0]])
last_active_t = float(offsets[active_indexes[-1]])
start_t = float(offsets[0])
end_t = float(offsets[-1])
leading = max(0.0, first_active_t - start_t)
trailing = max(0.0, end_t - last_active_t)
trimmed = max(0.0, last_active_t - first_active_t)
return {
"leading_idle_ratio": float(leading / duration_s),
"trailing_idle_ratio": float(trailing / duration_s),
"trimmed_duration_ratio": float(trimmed / duration_s) if duration_s > 0 else 1.0,
}
def _resample_to_length(points: list[Point], length: int) -> np.ndarray | None:
"""Trim idle padding, resample to ``length`` points, and peak-normalise.
Ported from ml_washdata/wash_ml/hybrid_curve_quality.py ``_resample_trace``.
Returns ``None`` when the trace is too short to be useful.
"""
# Trim leading/trailing idle.
active_indexes = [i for i, (_, p) in enumerate(points) if p > _IDLE_THRESHOLD_W]
if not active_indexes:
return None
pad_s = 60.0
start_off = max(points[0][0], points[active_indexes[0]][0] - pad_s)
end_off = points[active_indexes[-1]][0] + pad_s
trimmed = [(o, p) for o, p in points if start_off <= o <= end_off]
if len(trimmed) < 2:
return None
offsets = np.asarray([float(o) for o, _ in trimmed], dtype=float)
powers = np.asarray([max(0.0, float(p)) for _, p in trimmed], dtype=float)
valid = np.isfinite(offsets) & np.isfinite(powers)
offsets = offsets[valid]
powers = powers[valid]
if offsets.size < 2 or offsets[-1] <= offsets[0]:
return None
grid = np.linspace(offsets[0], offsets[-1], length)
trace = np.interp(grid, offsets, powers)
active = trace[trace > 0.5]
scale = float(np.percentile(active, 95)) if active.size else float(np.max(trace))
if not math.isfinite(scale) or scale <= 1e-6:
scale = 1.0
return np.clip(trace / scale, 0.0, 5.0)
def _trace_shape_descriptors(trace: np.ndarray) -> dict[str, float]:
"""Deterministic, scale-robust shape descriptors over a normalised trace.
Ported from ml_washdata/wash_ml/hybrid_curve_quality.py
``_trace_shape_descriptors``. Every value is NumPy-computable at runtime.
"""
if trace is None or trace.size < 4:
return {c: 0.0 for c in _SHAPE_COLUMNS}
trace = np.asarray(trace, dtype=float)
length = trace.size
total = float(np.sum(trace))
active_mask = trace > 0.5
active = trace[active_mask]
quarter = max(1, length // 4)
early_energy = float(np.sum(trace[:quarter]))
late_energy = float(np.sum(trace[-quarter:]))
mid = trace[quarter: length - quarter]
active_level = float(np.median(active)) if active.size else 0.0
mid_trough_depth = 0.0
if mid.size and active_level > 1e-6:
mid_trough_depth = float(np.clip(1.0 - float(np.min(mid)) / active_level, 0.0, 1.0))
diffs = np.diff(trace)
return {
"shape_active_fraction": float(np.mean(active_mask)),
"shape_early_energy_fraction": _safe_div(early_energy, total),
"shape_late_energy_fraction": _safe_div(late_energy, total),
"shape_mid_trough_depth": float(mid_trough_depth),
"shape_peak_density": _shape_peak_density(trace),
"shape_max_step_drop": float(max(0.0, -float(np.min(diffs)))) if diffs.size else 0.0,
"shape_max_step_rise": float(max(0.0, float(np.max(diffs)))) if diffs.size else 0.0,
"shape_active_cv": _safe_div(float(np.std(active)), float(np.mean(active))) if active.size else 0.0,
"shape_autocorr_lag1": _autocorr_lag1(trace),
"shape_derivative_sign_changes": _safe_div(
float(np.sum(np.abs(np.diff(np.sign(diffs))) > 0)), float(diffs.size)
) if diffs.size else 0.0,
"shape_plateau_ratio": _plateau_ratio(trace, active_level),
"shape_tail_slope": _safe_div(float(trace[-1] - trace[-quarter]), float(quarter)),
}
def _shape_peak_density(trace: np.ndarray, prominence: float = 0.2) -> float:
"""Prominent local maxima per sample."""
if trace.size < 3:
return 0.0
peaks = sum(
1
for i in range(1, trace.size - 1)
if trace[i] > trace[i - 1] and trace[i] >= trace[i + 1] and trace[i] >= prominence
)
return float(peaks) / float(trace.size)
def _autocorr_lag1(trace: np.ndarray) -> float:
"""Lag-1 autocorrelation (smoothness indicator)."""
centered = trace - float(np.mean(trace))
denom = float(np.dot(centered, centered))
if denom <= 1e-9:
return 0.0
return float(np.dot(centered[:-1], centered[1:]) / denom)
def _plateau_ratio(trace: np.ndarray, active_level: float, band: float = 0.12) -> float:
"""Fraction of trace within ``band`` of the running active level."""
if active_level <= 1e-6:
return 0.0
within = np.abs(trace - active_level) <= band
return float(np.mean(within & (trace > 0.5)))
def _log_ratio(ratio: float) -> float:
"""log(ratio) clamped to finite; 0.0 for non-positive or non-finite inputs."""
if not math.isfinite(ratio) or ratio <= 0:
return 0.0
return float(math.log(max(1e-6, ratio)))
def _safe_div(numerator: float, denominator: float) -> float:
if not math.isfinite(numerator) or not math.isfinite(denominator) or abs(denominator) <= 1e-9:
return 0.0
return float(numerator / denominator)
Reference in New Issue View Git Blame Copy Permalink