Compressor/backend/indicator_features/angle_bins.py

#!/usr/bin/env python3
"""Per-10-degree local statistics from every stored indicator cycle (read-only).

Each ``wave_sample_one`` row set is exactly one crankshaft revolution of chamber
pressure.  Map every raw sample to a crank angle (0-360, linear in sample index,
0 = keyphasor marker) and split it into 36 ten-degree bins.  On the **raw**
samples of each bin compute six simple numbers:

    mean      箱内平均压力
    range     箱内 max - min
    std       箱内标准差（波动）
    mad       箱内中位绝对偏差（稳健波动）
    slope     箱内压力对角度的线性斜率
    roughness 去掉箱内线性趋势后的残差 RMS（平滑 / 毛糙）

No resampling.  One long CSV row per (cycle, bin), so the file is large by
design and the caller should keep ``--per-day`` small.  Only running cycles
(``wave_file.rpm > 0``) are used.  Same-point before/after comparison is done by
``event_bins.py``.
"""
from __future__ import annotations

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

import numpy as np
import pymysql

BACKEND = Path(__file__).resolve().parents[1]
sys.path.insert(0, str(BACKEND))

from app.config import settings  # noqa: E402
from indicator_features.bins_common import BIN_METRICS, CYCLE_METRICS, FIELDNAMES  # noqa: E402
from indicator_features.indicator_features import CYLINDER_STAGE, otsu_threshold, smooth  # noqa: E402

CYLINDERS = ("一", "二", "三", "四", "五", "六")
SIDES = ("盖侧", "轴侧")
DEFAULT_BIN_DEG = 10.0


def connect() -> pymysql.connections.Connection:
    return pymysql.connect(
        host=settings.db_host, port=settings.db_port, user=settings.db_user,
        password=settings.db_password, database=settings.db_name, charset="utf8mb4",
        cursorclass=pymysql.cursors.DictCursor, connect_timeout=8, read_timeout=600,
    )


def point_name(unit: str, cylinder: str, side: str) -> str:
    return f"{unit}号机组{cylinder}缸压力{side}"


def fetch_file_rows(connection, point: str, start: str, end: str) -> list[dict]:
    """Running pressure files for one point in [start, end)."""
    with connection.cursor() as cursor:
        cursor.execute(
            """
            SELECT id, sample_time, rpm
            FROM wave_file
            WHERE point_name = %s AND measurement_type = '压力'
              AND rpm > 0
              AND sample_time >= %s AND sample_time < %s
            ORDER BY sample_time
            """,
            (point, start, end),
        )
        return list(cursor.fetchall())


def sample_per_day(rows: list[dict], per_day: int) -> list[dict]:
    if per_day <= 0:
        return rows
    by_day: dict = defaultdict(list)
    for row in rows:
        by_day[row["sample_time"].date()].append(row)
    picked: list[dict] = []
    for day in sorted(by_day):
        day_rows = by_day[day]
        if len(day_rows) <= per_day:
            picked.extend(day_rows)
            continue
        idx = np.linspace(0, len(day_rows) - 1, per_day).round().astype(int)
        picked.extend(day_rows[int(i)] for i in idx)
    return picked


def fetch_cycles(connection, file_ids: list[int], batch: int) -> dict[int, np.ndarray]:
    cycles: dict[int, list[tuple[int, float]]] = defaultdict(list)
    with connection.cursor() as cursor:
        for start in range(0, len(file_ids), batch):
            chunk = file_ids[start:start + batch]
            placeholders = ", ".join(["%s"] * len(chunk))
            cursor.execute(
                f"""
                SELECT wave_file_id, sample_index, CAST(signal_value AS FLOAT) AS sig
                FROM wave_sample_one
                WHERE wave_file_id IN ({placeholders})
                ORDER BY wave_file_id, sample_index
                """,
                tuple(chunk),
            )
            for row in cursor.fetchall():
                if row["sig"] is None:
                    continue
                cycles[int(row["wave_file_id"])].append((int(row["sample_index"]), float(row["sig"])))
    result: dict[int, np.ndarray] = {}
    for file_id, rows in cycles.items():
        rows.sort()
        result[file_id] = np.asarray([value for _, value in rows], dtype=float)
    return result


def _safe_ratio(numerator: float, denominator: float) -> float | None:
    return numerator / denominator if abs(denominator) > 1e-12 else None


def bin_features(signal: np.ndarray, bin_deg: float = DEFAULT_BIN_DEG) -> list[dict]:
    """Distribution / derivative / oscillation numbers for each angular bin (raw samples)."""
    p = np.asarray(signal, dtype=float)
    p = p[np.isfinite(p)]
    count = len(p)
    if count < 60:
        return []
    nbins = int(round(360.0 / bin_deg))
    angle = 360.0 * np.arange(count) / count
    features = []
    for b in range(nbins):
        low = b * bin_deg
        high = low + bin_deg
        mask = (angle >= low) & (angle < high)
        y = p[mask]
        x = angle[mask]
        if y.size < 3:
            continue
        center = float(np.mean(y))
        std = float(np.std(y))
        span_x = float(np.ptp(x))
        if span_x > 0:
            slope, intercept = np.polyfit(x, y, 1)
            residual = y - (slope * x + intercept)
        else:
            slope, residual = 0.0, y - center
        p10, p50, p90 = (float(v) for v in np.percentile(y, [10, 50, 90]))
        if std > 1e-12:
            z = (y - center) / std
            skew = float(np.mean(z ** 3))
            kurt = float(np.mean(z ** 4))
        else:
            skew = kurt = 0.0
        dy = np.diff(y)
        d2 = np.diff(dy) if dy.size > 1 else np.empty(0)
        if dy.size:
            abs_dy = np.abs(dy)
            d_abs_mean = float(np.mean(abs_dy))
            d_abs_max = float(np.max(abs_dy))
            d_rms = float(np.sqrt(np.mean(dy ** 2)))
            tv = float(np.sum(abs_dy))
            if d2.size:
                d2_rms = float(np.sqrt(np.mean(d2 ** 2)))
                d2_abs_max = float(np.max(np.abs(d2)))
                sign = np.sign(dy)
                sign[sign == 0] = 1
                zero_cross = int(np.count_nonzero(np.diff(sign) != 0))
            else:
                d2_rms = d2_abs_max = 0.0
                zero_cross = 0
        else:
            d_abs_mean = d_abs_max = d_rms = d2_rms = d2_abs_max = tv = 0.0
            zero_cross = 0
        n_extrema = int(np.count_nonzero(np.diff(np.sign(np.diff(y))) != 0)) if y.size > 2 else 0
        centered = y - center
        mc = np.sign(centered)
        mc[mc == 0] = 1
        mean_cross = int(np.count_nonzero(np.diff(mc) != 0)) if y.size > 1 else 0
        # Hjorth mobility / complexity (shape/irregularity)
        var_y = float(np.var(y))
        var_dy = float(np.var(dy)) if dy.size else 0.0
        var_d2 = float(np.var(d2)) if d2.size else 0.0
        mobility = float(np.sqrt(var_dy / var_y)) if var_y > 1e-18 else 0.0
        complexity = (float(np.sqrt(var_d2 / var_dy)) / mobility) if var_dy > 1e-18 and mobility > 1e-18 else 0.0
        features.append({
            "bin": b, "bin_start": low, "n_raw": int(y.size),
            "mean": center, "range": float(np.max(y) - np.min(y)), "std": std,
            "mad": float(np.median(np.abs(y - np.median(y)))), "slope": float(slope),
            "roughness": float(np.sqrt(np.mean(residual ** 2))),
            "p10": p10, "p50": p50, "p90": p90, "iqr": p90 - p10,
            "skew": skew, "kurt": kurt,
            "d_abs_mean": d_abs_mean, "d_abs_max": d_abs_max, "d_rms": d_rms,
            "d2_rms": d2_rms, "zero_cross": zero_cross, "n_extrema": n_extrema,
            "d2_abs_max": d2_abs_max, "tv": tv, "mean_cross": mean_cross,
            "mobility": mobility, "complexity": complexity,
        })
    return features


def cycle_timing(signal: np.ndarray) -> dict:
    """Whole-cycle plateau / valve-event features (Otsu plateaus + derivative corners)."""
    p = np.asarray(signal, dtype=float)
    p = p[np.isfinite(p)]
    count = len(p)
    out = {name: None for name in CYCLE_METRICS}
    if count < 60:
        return out
    ps = smooth(p, 21)
    angle = 360.0 * np.arange(count) / count
    threshold = otsu_threshold(ps)
    low = ps[ps <= threshold]
    high = ps[ps > threshold]
    suction = float(np.median(low)) if low.size else float(np.min(ps))
    discharge = float(np.median(high)) if high.size else float(np.max(ps))
    high_mask = ps >= threshold
    out["p_suction"] = suction
    out["p_discharge"] = discharge
    out["pressure_ratio"] = _safe_ratio(discharge, suction)
    out["suction_frac"] = float(np.mean(~high_mask))
    out["discharge_frac"] = float(np.mean(high_mask))
    out["angle_min"] = float(angle[int(np.argmin(ps))])
    out["angle_max"] = float(angle[int(np.argmax(ps))])
    signed = np.sign(ps - threshold)
    signed[signed == 0] = 1
    rising = np.flatnonzero(np.diff(signed) > 0)
    falling = np.flatnonzero(np.diff(signed) < 0)
    if rising.size:
        out["cross_up"] = float(angle[rising[0]])
    if falling.size:
        out["cross_down"] = float(angle[falling[0]])
    dp = np.gradient(ps)
    out["max_rise"] = float(np.max(dp))
    out["max_fall"] = float(np.min(dp))
    out["max_rise_angle"] = float(angle[int(np.argmax(dp))])
    out["max_fall_angle"] = float(angle[int(np.argmin(dp))])
    return out


def main() -> None:
    parser = argparse.ArgumentParser(description="功图每10°分箱局部统计（只读，抽样）")
    parser.add_argument("--unit", required=True, choices=("7", "8", "9"))
    parser.add_argument("--cylinders", nargs="+", default=list(CYLINDERS), choices=CYLINDERS)
    parser.add_argument("--sides", nargs="+", default=list(SIDES), choices=SIDES)
    parser.add_argument("--start", required=True, help="YYYY-MM-DD")
    parser.add_argument("--end", required=True, help="YYYY-MM-DD（不含）")
    parser.add_argument("--per-day", type=int, default=4, help="每测点每天最多取几个运行文件")
    parser.add_argument("--bin-deg", type=float, default=DEFAULT_BIN_DEG, help="角度箱宽度（默认10°，可选5/15/20）")
    parser.add_argument("--batch", type=int, default=200, help="IN 查询批大小")
    parser.add_argument("--output", required=True)
    args = parser.parse_args()
    if not (0 < args.bin_deg <= 90):
        parser.error("--bin-deg 必须在 (0, 90] 之间")

    connection = connect()
    rows_out: list[dict] = []
    try:
        for cylinder in args.cylinders:
            for side in args.sides:
                point = point_name(args.unit, cylinder, side)
                files = fetch_file_rows(connection, point, args.start, args.end)
                if not files:
                    print(f"{point}: 无运行文件", flush=True)
                    continue
                picked = sample_per_day(files, args.per_day)
                cycles = fetch_cycles(connection, [int(r["id"]) for r in picked], args.batch)
                id_to_time = {int(r["id"]): r for r in picked}
                kept = 0
                for file_id, signal in cycles.items():
                    bins = bin_features(signal, args.bin_deg)
                    if not bins:
                        continue
                    meta = id_to_time[file_id]
                    finite = signal[np.isfinite(signal)]
                    base = {
                        "unit": args.unit,
                        "cylinder": cylinder,
                        "stage": CYLINDER_STAGE[cylinder],
                        "side": side,
                        "point_name": point,
                        "wave_file_id": file_id,
                        "sample_time": str(meta["sample_time"]),
                        "rpm": float(meta["rpm"] or 0.0),
                        "n_samples": int(finite.size),
                        "p_min": float(finite.min()),
                        "p_max": float(finite.max()),
                        "p_span": float(finite.max() - finite.min()),
                    }
                    timing = cycle_timing(signal)
                    for item in bins:
                        rows_out.append({**base, **item, **timing})
                    kept += 1
                days = len({r["sample_time"].date() for r in picked})
                print(f"{point}: 运行文件 {len(files)}，抽样 {len(picked)}（{days}天），出周期 {kept}", flush=True)
    finally:
        connection.close()

    if not rows_out:
        print("没有特征行，退出")
        return
    rows_out.sort(key=lambda r: (r["cylinder"], r["side"], r["sample_time"], r["bin"]))
    output = Path(args.output)
    output.parent.mkdir(parents=True, exist_ok=True)
    with output.open("w", encoding="utf-8-sig", newline="") as handle:
        writer = csv.DictWriter(handle, fieldnames=FIELDNAMES)
        writer.writeheader()
        writer.writerows(rows_out)
    periods = len({(r["point_name"], r["wave_file_id"]) for r in rows_out})
    print(f"周期 {periods:,}，特征行 {len(rows_out):,} -> {output}")
    print(f"时间：{datetime.now():%Y-%m-%d %H:%M:%S}")


if __name__ == "__main__":
    main()