#!/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()