Compressor/backend/diagnose_drift_running.py

from __future__ import annotations

import os
import sys
from time import monotonic

import numpy as np
import pandas as pd
import pymysql

sys.path.insert(0, os.path.dirname(os.path.abspath(__file__)))
from app.config import settings  # noqa: E402

BATCHES = (30, 31, 32)
CACHE = os.path.join(
    os.path.dirname(os.path.dirname(os.path.abspath(__file__))), "cache", "pks_minute_levels.csv"
)
SPEED_THRESHOLD = 300.0  # rpm: running gate


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.SSCursor,
        connect_timeout=settings.db_connect_timeout,
        read_timeout=3600,
        write_timeout=120,
    )


def wall_offset() -> int:
    connection = connect()
    try:
        with connection.cursor() as cursor:
            cursor.execute("SELECT sample_time, UNIX_TIMESTAMP(sample_time) AS u "
                           "FROM pks_long_sample LIMIT 1")
            wall, uni = cursor.fetchone()
    finally:
        connection.close()
    return int(pd.Timestamp(wall).value // 10 ** 9) - int(uni)


def load_minute_levels(offset: int) -> pd.DataFrame:
    if os.path.exists(CACHE):
        f = pd.read_csv(CACHE, parse_dates=["time"])
        print(f"   loaded minute levels from cache: {len(f):,} rows")
        return f
    t0 = monotonic()
    connection = connect()
    rows = []
    try:
        with connection.cursor() as cursor:
            cursor.execute(
                "SELECT import_batch_id AS b, UNIX_TIMESTAMP(sample_time) DIV 60 AS bk, "
                "COUNT(*) AS cnt, AVG(YSJ_1) AS lev, AVG(YSJ_41) AS rpm "
                "FROM pks_long_sample GROUP BY import_batch_id, bk "
                "ORDER BY import_batch_id, bk"
            )
            while True:
                chunk = cursor.fetchmany(200_000)
                if not chunk:
                    break
                rows.extend(chunk)
    finally:
        connection.close()
    f = pd.DataFrame(rows, columns=["b", "bk", "cnt", "lev", "rpm"])
    f["time"] = pd.to_datetime(f["bk"] * 60 + offset, unit="s")
    f[["b", "bk", "cnt", "lev", "rpm", "time"]].to_csv(CACHE, index=False)
    print(f"   aggregated minute levels: {len(f):,} rows in {monotonic()-t0:.0f}s")
    return f


def mad_std(x: np.ndarray) -> float:
    med = float(np.median(x))
    return 1.4826 * float(np.median(np.abs(x - med))), med


def flag_summary(values: np.ndarray, tau: float, minute_step: int = 1) -> dict:
    mask = np.abs(values) > tau
    if mask.size == 0:
        return {"count": 0, "frac": 0.0, "top": []}
    diff = np.diff(mask.astype(np.int8))
    starts = np.flatnonzero(diff == 1) + 1
    ends = np.flatnonzero(diff == -1) + 1
    if mask[0]:
        starts = np.concatenate(([0], starts))
    if mask[-1]:
        ends = np.concatenate((ends, [len(mask)]))
    if len(starts) == 0 or len(ends) == 0:
        return {"count": 0, "frac": round(float(mask.mean()), 4), "top": []}
    ends = ends[: len(starts)]
    durs = (ends - starts) * minute_step
    peaks = [float(np.max(np.abs(values[s:e]))) for s, e in zip(starts, ends)]
    order = np.argsort(peaks)[::-1][:5]
    return {
        "count": int(len(durs)),
        "durs_min": [int(x) for x in durs],
        "frac": round(float(mask.mean()), 4),
        "top": [{"peak": round(peaks[i], 2), "dur_min": int(durs[i])} for i in order],
    }


def build_run_blocks(times: pd.DatetimeIndex, running: np.ndarray) -> np.ndarray:
    """block id per minute; merge running minutes with gaps <= 6 min."""
    idx_run = np.flatnonzero(running)
    if len(idx_run) == 0:
        return np.full(len(times), -1, dtype=np.int64)
    block = np.full(len(times), -1, dtype=np.int64)
    bid = 0
    prev_t = None
    for i in idx_run:
        if prev_t is None or (times[i] - prev_t).total_seconds() > 360:
            bid += 1
        block[i] = bid
        prev_t = times[i]
    return block


def main() -> None:
    offset = wall_offset()
    f = load_minute_levels(offset)
    per: dict[int, pd.DataFrame] = {b: g.sort_values("time").set_index("time") for b, g in f.groupby("b")}

    result: dict = {}
    for b in BATCHES:
        df = per[b].copy()
        df["running"] = (df["rpm"] > SPEED_THRESHOLD).to_numpy()
        df["block"] = build_run_blocks(df.index, df["running"].to_numpy())
        dur = df.groupby("block")["lev"].transform("size")
        df["warmup"] = df.groupby("block").cumcount() < 120  # first 2h of each run block
        valid = df["running"] & ~df["warmup"] & df["lev"].notna()
        n_run = int(df["running"].sum())
        n_block = int((df["block"] >= 0).groupby(df["block"]).ngroups)

        run_levels = df.loc[df["running"], "lev"]
        rpm_run = df.loc[df["running"], "rpm"]

        # self-reference on valid (running, non-warmup) minutes
        ser = df["lev"].where(df["running"])
        local = ser.rolling("6h", min_periods=20).median()
        ref = ser.rolling("20d", min_periods=200).median()
        e_series = (local - ref).loc[valid].dropna()
        e = e_series.to_numpy()

        if len(e) == 0:
            self_info = None
        else:
            sig, med = mad_std(e)
            self_info = {
                "n_min": int(len(e)), "med": round(med, 3), "sigma": round(sig, 3),
                "tau4": round(4 * sig, 3), "max_abs": round(float(np.abs(e).max()), 3),
                "flag": {f"{k}s": flag_summary(e, k * sig) for k in (3, 4, 5)},
            }

        # plateau (steady level) per run block and drift of plateau across blocks
        plateau = (df.loc[valid, "lev"].groupby(df.loc[valid, "block"]).median())
        pstart = df.loc[df["running"]].groupby("block")["rpm"].apply(lambda s: s.index[0])
        prev_ref = plateau.rolling(30, min_periods=8).median().shift(1)
        dev = (plateau - prev_ref).dropna()
        drift = None
        if len(dev) > 3:
            sig_d, med_d = mad_std(dev.to_numpy())
            flagged = dev[np.abs(dev) > 4 * sig_d]
            drift = {
                "n_blocks": int(len(plateau)),
                "blocks_per_month": round(len(plateau) / 13.0, 1),
                "plateau_p50": round(float(plateau.median()), 2),
                "plateau_minmax": [round(float(plateau.min()), 2), round(float(plateau.max()), 2)],
                "dev_sigma": round(sig_d, 3),
                "dev_mad_med": round(med_d, 3),
                "blocks_over_4sigma": int(len(flagged)),
                "top_blocks": [
                    {"start": str(idx)[:16], "plateau": round(float(pl), 2),
                     "dev": round(float(dv), 2)}
                    for idx, (pl, dv) in flagged.head(5).items()
                ],
            }

        stats = {
            "n_run_minutes": n_run, "run_minutes_share": round(n_run / len(df), 4),
            "n_run_blocks": n_block, "blocks_over_24h": int((dur >= 1440).sum()),
            "rpm_run_p50": round(float(rpm_run.median()), 0),
            "run_level_p50": round(float(run_levels.median()), 2),
            "run_level_minmax": [round(float(run_levels.min()), 1), round(float(run_levels.max()), 1)],
        }
        result[str(b)] = {"gating": stats, "self_within_run": self_info, "plateau_drift": drift}
        print(f"[{b}] run_min={n_run:,} ({stats['run_minutes_share']:.1%}) blocks={n_block} "
              f"rpm_p50={stats['rpm_run_p50']:.0f} run_level_p50={stats['run_level_p50']}")
        if self_info:
            print(f"   within-run self: sigma={self_info['sigma']} tau4={self_info['tau4']} "
                  f"max={self_info['max_abs']} frac@4s={self_info['flag']['4s']['frac']} "
                  f"runs@4s={self_info['flag']['4s']['count']}")
        if drift:
            print(f"   plateau drift: blocks={drift['n_blocks']} dev_sigma={drift['dev_sigma']} "
                  f"over_4s={drift['blocks_over_4sigma']}")

    # cross-machine on common running minutes (all three machines running)
    frame = pd.DataFrame({str(b): per[b]["lev"] for b in BATCHES})
    frame["rpm"] = per[30]["rpm"]
    common = (frame["30"] > 0) & (frame["31"] > 0) & (frame["32"] > 0) & (frame["rpm"] > SPEED_THRESHOLD)
    sub = frame.loc[common, [str(b) for b in BATCHES]]
    cross = {}
    if len(sub) > 1000:
        for u in BATCHES:
            others = [str(o) for o in BATCHES if o != u]
            ref = sub[others].median(axis=1)
            d = (sub[str(u)] - ref).to_numpy()
            sig, med = mad_std(d)
            cross[str(u)] = {
                "n_min": int(len(d)), "med_bias": round(med, 3), "sigma": round(sig, 3),
                "tau4": round(4 * sig, 3), "max_abs": round(float(np.abs(d).max()), 3),
                "flag": {f"{k}s": flag_summary(d, k * sig) for k in (3, 4, 5)},
            }
            print(f"[cross u={u}] common_run_min={len(d):,} bias={med:.2f} sigma={sig:.3f} "
                  f"frac@4s={cross[str(u)]['flag']['4s']['frac']} runs@4s={cross[str(u)]['flag']['4s']['count']}")
    result["cross_common_running"] = cross

    print("\n===== JSON SUMMARY =====")
    print(__import__("json").dumps(result, indent=2, ensure_ascii=False))


if __name__ == "__main__":
    main()