Compressor/backend/indicator_features/scan_bins.py

#!/usr/bin/env python3
"""Full-timeline unsupervised change scan over the per-bin/-cycle features.

For every ``(point, bin, metric)`` the whole available time series is searched
for the single best two-segment split (a step change), scored as the between
segment mean difference over the pooled standard deviation (Cohen's d).  The
point-wide (all-bins) part is removed per metric so only the bins that move
*differently* rank high -- the same local/whole decomposition used by
``event_bins.py``, but now applied to the entire timeline without needing a
known event or a pre/post window.

Output: a ranked candidate list and a CSV.  Read-only.
"""
from __future__ import annotations

import argparse
import csv
from collections import defaultdict
from datetime import datetime
from pathlib import Path

import numpy as np

from bins_common import ANGLE_METRICS, BIN_METRICS, CYCLE_METRICS, VOLATILE, apply_normalization
from event_bins import EVENTS, load

MIN_DEG = 1.0  # minimum pooled scale for angle metrics


def best_split(times, values, min_days: int = 3):
    """Best two-segment step (mean difference / pooled std). Returns (score, index)."""
    n = len(values)
    if n < 8:
        return None, None
    x = np.asarray(values, dtype=float)
    csum = np.cumsum(x)
    csum2 = np.cumsum(x ** 2)
    total, total_sq = float(csum[-1]), float(csum2[-1])
    best_score, best_i = 0.0, None
    min_seg = max(5, n // 12)
    for i in range(min_seg, n - min_seg):
        if (times[i - 1] - times[0]).days < min_days or (times[-1] - times[i]).days < min_days:
            continue
        left, right = i, n - i
        sl, sl2 = float(csum[i - 1]), float(csum2[i - 1])
        mean_l, mean_r = sl / left, (total - sl) / right
        var_l = max(sl2 / left - mean_l ** 2, 0.0)
        var_r = max((total_sq - sl2) / right - mean_r ** 2, 0.0)
        pooled = float(np.sqrt(0.5 * (var_l + var_r)))
        if pooled <= 1e-12:
            continue
        score = abs(mean_r - mean_l) / pooled
        if score > best_score:
            best_score, best_i = score, i
    return (best_score, best_i) if best_i is not None else (None, None)


def daily_series(rows, getter):
    return sorted(((r["time"], getter(r)) for r in rows), key=lambda t: t[0])


def main() -> None:
    parser = argparse.ArgumentParser(description="全时间轴 分箱特征 无监督变点扫描")
    parser.add_argument("--unit", required=True, choices=("7", "8", "9"))
    parser.add_argument("--csv", nargs="+", required=True)
    parser.add_argument("--top", type=int, default=40)
    parser.add_argument("--z", type=float, default=5.0, help="最小 |Cohen's d| 阈值")
    parser.add_argument("--normalize", choices=("none", "span"), default="none",
                        help="span: 用本周期 p_span 归一再扫（压工况）")
    parser.add_argument("--out", default="")
    args = parser.parse_args()

    rows = []
    for path in args.csv:
        rows.extend(load(path))
    apply_normalization(rows, args.normalize)
    print(f"载入 {len(rows):,} 行，测点 {len({r['point_name'] for r in rows})}，归一 {args.normalize}")

    # records[(point,metric)][bin] = (time, score)
    bin_scores = defaultdict(dict)
    cycle_scores = defaultdict(list)
    points = sorted({(r["cylinder"], r["side"]) for r in rows})

    for cylinder, side in points:
        point_rows = [r for r in rows if r["cylinder"] == cylinder and r["side"] == side]
        # per-bin
        by_bin = defaultdict(list)
        for r in point_rows:
            by_bin[r["bin"]].append(r)
        for name in BIN_METRICS:
            for b, items in by_bin.items():
                series = daily_series(items, lambda r, nm=name: r[nm])
                series = [(t, v) for t, v in series if v is not None]
                times = [t for t, _ in series]
                values = [v for _, v in series]
                score, idx = best_split(times, values)
                if score is None:
                    continue
                bin_scores[(f"{cylinder}缸{side}", name)][b] = (score, times[idx], np.mean(values[:idx]), np.mean(values[idx:]))
        # whole-cycle
        for name in CYCLE_METRICS:
            if name in ANGLE_METRICS:
                continue
            series = daily_series(point_rows, lambda r, nm=name: r[nm])
            series = [(t, v) for t, v in series if v is not None]
            times = [t for t, _ in series]
            values = [v for _, v in series]
            score, idx = best_split(times, values)
            if score is not None:
                cycle_scores[f"{cylinder}缸{side}"].append((name, score, times[idx], np.mean(values[:idx]), np.mean(values[idx:])))

    records = []
    for (point, metric), per_bin in bin_scores.items():
        if len(per_bin) < 4:
            continue
        raw = np.asarray([v[0] for v in per_bin.values()], dtype=float)
        median_score = float(np.median(raw))
        for b, (score, split_time, before, after) in per_bin.items():
            local = score - median_score  # remove whole-point level shift
            if local >= args.z:
                records.append({
                    "family": "bin", "point": point, "bin": b,
                    "metric": metric, "score": local, "raw_score": score,
                    "split": split_time, "before": before, "after": after,
                    "volatile": metric in VOLATILE,
                })
    for point, items in cycle_scores.items():
        for name, score, split_time, before, after in items:
            if score >= args.z:
                records.append({
                    "family": "cycle", "point": point, "bin": "", "metric": name,
                    "score": score, "raw_score": score, "split": split_time,
                    "before": before, "after": after, "volatile": False,
                })

    records.sort(key=lambda r: -r["score"])
    print(f"候选 {len(records)} 条（分箱局部分数 >= {args.z}，整周期 >= {args.z}）")
    print(f"{'#':>3}  {'家族':<6} {'测点':<9} {'角度箱':<8} {'指标':<13} {'分数':>6} {'突变时间':<17} {'前':>9} {'后':>9}")
    for i, r in enumerate(records[: args.top], 1):
        binst = f"{r['bin'] * 10:g}-{r['bin'] * 10 + 10:g}°" if r["family"] == "bin" else "整周期"
        near = min(EVENTS, key=lambda e: abs((datetime.fromisoformat(e[1]) - r["split"]).days))
        gap = abs((datetime.fromisoformat(near[1]) - r["split"]).days)
        print(f"{i:>3}  {r['family']:<6} {r['point']:<9} {binst:<8} {r['metric']:<13} {r['score']:>6.1f} "
              f"{r['split']:%Y-%m-%d %H:%M} {r['before']:>9.4g} {r['after']:>9.4g}   ~{near[0]}号{near[1]}({gap}天)")

    if args.out and records:
        out = Path(args.out)
        out.parent.mkdir(parents=True, exist_ok=True)
        with out.open("w", encoding="utf-8-sig", newline="") as handle:
            writer = csv.DictWriter(handle, fieldnames=["family", "point", "bin", "metric", "score", "raw_score", "split", "before", "after", "volatile"])
            writer.writeheader()
            for r in records:
                writer.writerow({**r, "split": r["split"].isoformat()})
        print(f"明细 -> {out}")
    print(f"时间：{datetime.now():%Y-%m-%d %H:%M:%S}")


if __name__ == "__main__":
    main()