Compressor/backend/indicator_features/analyze.py

#!/usr/bin/env python3
"""Compare indicator features with three normalizations (read-only).

  * peer     : within the same stage and side at the same timestamp;
  * chamber  : cover side vs shaft side of the same cylinder;
  * self     : the cylinder own trailing baseline (median/MAD) -> z-score.

Absolute levels depend on the process, so only relative quantities are scored.
"""
from __future__ import annotations

import argparse
import csv
from collections import defaultdict

import numpy as np

DEAD_DISCHARGE = 0.20
PEER_FEATURES = ("p_discharge", "pressure_ratio", "n_comp", "discharge_frac", "efficiency", "indicated_work")
PROFILE = tuple(f"prof_{i:02d}" for i in range(12))


def to_float(value):
    try:
        result = float(value)
    except (TypeError, ValueError):
        return None
    return result if np.isfinite(result) else None


def load(path: str) -> list[dict]:
    rows = list(csv.DictReader(open(path, encoding="utf-8-sig")))
    for row in rows:
        for key in PEER_FEATURES + PROFILE + ("p_suction", "p_span"):
            row[key] = to_float(row.get(key))
    return rows


def healthy(row: dict, exclude: set | None = None) -> bool:
    if exclude and (row["cylinder"], row["side"]) in exclude:
        return False
    pd = row.get("p_discharge")
    ps = row.get("p_suction")
    span = row.get("p_span")
    return pd is not None and pd >= DEAD_DISCHARGE and ps and ps > 0.05 and span and span >= 0.3


def series_medians(rows: list[dict], feature: str) -> dict:
    by_series = defaultdict(list)
    for row in rows:
        value = row.get(feature)
        if value is not None:
            by_series[(row["cylinder"], row["side"])].append(value)
    return {key: float(np.median(values)) for key, values in by_series.items() if values}


def flag_channels(rows: list[dict]) -> dict:
    """Flag weak/dead channels relative to the stage-wide expected discharge level."""
    medians = series_medians(rows, "p_discharge")
    reference = np.median([v for v in medians.values() if v > DEAD_DISCHARGE]) if medians else 0.0
    flags = {}
    for key, value in medians.items():
        if value < DEAD_DISCHARGE:
            flags[key] = "dead"
        elif reference > 0 and value < 0.6 * reference:
            flags[key] = "weak"
        elif reference > 0 and value > 1.5 * reference:
            flags[key] = "high"
    return flags


def peer_deviations(rows: list[dict], exclude: set) -> dict:
    by_key = defaultdict(dict)
    for row in rows:
        if healthy(row, exclude):
            by_key[(row["stage"], row["side"], row["sample_time"])][row["cylinder"]] = row
    result = defaultdict(list)
    for (stage, side, _time), members in by_key.items():
        if len(members) < 2:
            continue
        for feature in PEER_FEATURES:
            values = [m[feature] for m in members.values() if m[feature] is not None]
            if len(values) < 2:
                continue
            med = float(np.median(values))
            if abs(med) < 1e-9:
                continue
            for cylinder, member in members.items():
                if member[feature] is not None:
                    result[(cylinder, side, feature)].append((member[feature] - med) / med)
    return {key: float(np.mean(values)) for key, values in result.items()}


def chamber_asymmetry(rows: list[dict], exclude: set) -> dict:
    by_key = defaultdict(dict)
    for row in rows:
        by_key[(row["cylinder"], row["sample_time"])][row["side"]] = row
    result = defaultdict(list)
    for (cylinder, _time), pair in by_key.items():
        if "盖侧" not in pair or "轴侧" not in pair:
            continue
        cover, shaft = pair["盖侧"], pair["轴侧"]
        if not (healthy(cover, exclude) and healthy(shaft, exclude)):
            continue
        for feature in ("p_discharge", "pressure_ratio", "n_comp"):
            if cover[feature] is not None and shaft[feature] is not None and abs(shaft[feature]) > 1e-9:
                result[(cylinder, feature)].append((cover[feature] - shaft[feature]) / abs(shaft[feature]))
        profile_cover = np.array([cover[name] for name in PROFILE], dtype=float)
        profile_shaft = np.array([shaft[name] for name in PROFILE], dtype=float)
        if np.all(np.isfinite(profile_cover)) and np.all(np.isfinite(profile_shaft)):
            result[(cylinder, "profile_l2")].append(float(np.linalg.norm(profile_cover - profile_shaft)))
    return {key: float(np.mean(values)) for key, values in result.items()}


def self_baseline(rows: list[dict], exclude: set) -> dict:
    by_series = defaultdict(list)
    for row in rows:
        if healthy(row, exclude):
            by_series[(row["cylinder"], row["side"])].append(row)
    result = {}
    for (cylinder, side), series in by_series.items():
        series.sort(key=lambda r: r["sample_time"])
        times = [r["sample_time"] for r in series]
        for feature in ("p_discharge", "pressure_ratio", "n_comp"):
            values = np.array([r[feature] if r[feature] is not None else np.nan for r in series], dtype=float)
            valid = values[np.isfinite(values)]
            if valid.size < 5:
                continue
            center = float(np.median(valid))
            mad = float(np.median(np.abs(valid - center)))
            scale = mad * 1.4826 or float(np.std(valid))
            if scale <= 1e-9:
                continue
            z = (values - center) / scale
            result[(cylinder, side, feature)] = {
                "center": center, "scale": scale,
                "first_half_z": float(np.nanmean(z[: len(z) // 2])),
                "second_half_z": float(np.nanmean(z[len(z) // 2:])),
                "max_abs_z": float(np.nanmax(np.abs(z))) if np.isfinite(z).any() else None,
                "last_time": times[-1],
            }
    return result


def main() -> None:
    parser = argparse.ArgumentParser(description="功图特征三种归一化对照")
    parser.add_argument("csv_path")
    args = parser.parse_args()
    rows = load(args.csv_path)
    flags = flag_channels(rows)
    exclude = set(flags)
    units = sorted({(r["cylinder"], r["side"]) for r in rows})
    print(f"总行 {len(rows)}；测点 {len(units)}")
    print(f"通道标记(排除对照)：{flags if flags else '无'}")

    print("\n=== 同级缸间偏差（同侧、同刻，相对中位数）===")
    for (cylinder, side, feature), value in sorted(peer_deviations(rows, exclude).items()):
        print(f"  {cylinder}缸{side:<3} {feature:<16} {value * 100:+.2f}%")

    print("\n=== 同缸双腔差（盖 - 轴）===")
    for (cylinder, feature), value in sorted(chamber_asymmetry(rows, exclude).items()):
        if feature == "profile_l2":
            print(f"  {cylinder}缸 profile L2 {value:.3f}")
        else:
            print(f"  {cylinder}缸 {feature:<16} {value * 100:+.2f}%")

    print("\n=== 本缸自身基线 z（前段 vs 后段，沿时间趋势）===")
    for (cylinder, side, feature), stat in sorted(self_baseline(rows, exclude).items()):
        print(f"  {cylinder}缸{side:<3} {feature:<16} 前段z={stat['first_half_z']:+.2f} 后段z={stat['second_half_z']:+.2f} max|z|={stat['max_abs_z']:.2f}")


if __name__ == "__main__":
    main()