Compressor/backend/indicator_features/indicator_features.py

#!/usr/bin/env python3
"""Phase-level features from one compressor indicator diagram cycle.

One stored cycle (wave_sample_one) is the pressure of a single cylinder chamber
over one crankshaft revolution (~60 ms at ~996 RPM).

Design goals:
  * robust two-plateau (suction P_s / discharge P_d) detection via Otsu split;
  * physical phase alignment (search the crank offset that maximises the PV loop
    area) so indicated work and a pressure-conversion efficiency are meaningful;
  * a normalised phase profile for peer / double-chamber / self-baseline review.

Compression staging (confirmed with the site):
  一级(1st) = 2/4/6 缸,  二级(2nd) = 3/5 缸,  三级(3rd) = 1 缸.
Absolute levels depend on the process, so features are meant to be compared.
"""
from __future__ import annotations

import numpy as np

CYLINDER_BORE_MM = {"一": 360.0, "二": 490.0, "三": 390.0, "四": 490.0, "五": 390.0, "六": 490.0}
CYLINDER_STAGE = {"一": "三级", "二": "一级", "三": "二级", "四": "一级", "五": "二级", "六": "一级"}
PISTON_STROKE_MM = 148.0
CONNECTING_ROD_LENGTH_MM = 460.0
CLEARANCE_VOLUME_L_BY_BORE = {490.0: 1.51, 390.0: 0.74, 360.0: 0.62}
PROFILE_BINS = 12
K_ISENTROPIC = 1.30


def smooth(values: np.ndarray, window: int) -> np.ndarray:
    window = max(1, int(window))
    if window <= 1 or len(values) < window:
        return values
    kernel = np.ones(window, dtype=float) / window
    return np.convolve(values, kernel, mode="same")


def otsu_threshold(values: np.ndarray, bins: int = 96) -> float:
    hist, edges = np.histogram(values, bins=bins)
    centers = (edges[:-1] + edges[1:]) / 2.0
    total = float(values.size)
    if total <= 0:
        return float(np.mean(values))
    sum_total = float(np.sum(centers * hist))
    sum_before = 0.0
    weight_before = 0.0
    best_variance = -1.0
    threshold = float(centers[0])
    for index in range(len(hist)):
        weight_before += float(hist[index])
        if weight_before == 0.0:
            continue
        weight_after = total - weight_before
        if weight_after == 0.0:
            break
        sum_before += centers[index] * hist[index]
        mean_before = sum_before / weight_before
        mean_after = (sum_total - sum_before) / weight_after
        variance = weight_before * weight_after * (mean_before - mean_after) ** 2
        if variance > best_variance:
            best_variance = variance
            threshold = float(centers[index])
    return threshold


def volume_from_angle(count: int, bore_mm: float) -> np.ndarray:
    crank_radius = PISTON_STROKE_MM / 2.0
    clearance = CLEARANCE_VOLUME_L_BY_BORE.get(bore_mm, 1.0)
    theta = np.linspace(0.0, 2.0 * np.pi, count, endpoint=False)
    travel = (
        crank_radius * (1.0 - np.cos(theta))
        + CONNECTING_ROD_LENGTH_MM
        - np.sqrt(CONNECTING_ROD_LENGTH_MM**2 - (crank_radius * np.sin(theta)) ** 2)
    )
    piston_area = np.pi * (bore_mm / 2.0) ** 2
    return clearance + piston_area * travel / 1_000_000.0


def _plateau_levels(p: np.ndarray) -> tuple[float, float, float]:
    threshold = otsu_threshold(p)
    low = p[p <= threshold]
    high = p[p > threshold]
    suction = float(np.median(low)) if low.size else float(np.min(p))
    discharge = float(np.median(high)) if high.size else float(np.max(p))
    return suction, discharge, threshold


def _align_offset(p: np.ndarray, dvolume: np.ndarray) -> tuple[int, float]:
    """Offset (in samples) that maximises the indicated work -sum(p dV)."""
    count = len(p)
    cross = np.real(np.fft.ifft(np.conj(np.fft.fft(p)) * np.fft.fft(dvolume)))
    work = -cross
    offset = int(np.argmax(work))
    return offset, float(work[offset])


def _polytropic_index(volume: np.ndarray, pressure: np.ndarray, low: float, high: float, rising: bool) -> float | None:
    dp = np.diff(pressure, prepend=pressure[0])
    band = (pressure > low + 0.15 * (high - low)) & (pressure < low + 0.85 * (high - low))
    direction = dp > 0 if rising else dp < 0
    mask = band & direction
    v = volume[mask]
    p = pressure[mask]
    valid = (v > 0) & (p > 0)
    if np.count_nonzero(valid) < 8:
        return None
    slope, _ = np.polyfit(np.log(v[valid]), np.log(p[valid]), 1)
    return float(-slope)


def cycle_features(signal: np.ndarray, cylinder: str, smooth_window: int = 21) -> dict:
    signal = np.asarray(signal, dtype=float)
    signal = signal[np.isfinite(signal)]
    count = len(signal)
    if count < 60:
        return {}
    p = smooth(signal, smooth_window)
    p = np.where(np.isfinite(p), p, signal)
    suction, discharge, threshold = _plateau_levels(p)
    span = discharge - suction
    volume = volume_from_angle(count, CYLINDER_BORE_MM.get(cylinder, 490.0))
    dvolume = np.diff(volume, append=volume[0])

    offset, indicated_work = _align_offset(p, dvolume)
    aligned_volume = np.roll(volume, -offset)
    aligned_dvolume = np.roll(dvolume, -offset)

    # Isentropic ideal work (MPa*L = kJ) and pressure-conversion efficiency.
    v_suction = float(np.max(aligned_volume))
    ratio = (discharge / suction) if suction > 0 else None
    if ratio and ratio > 1:
        ideal = (K_ISENTROPIC / (K_ISENTROPIC - 1.0)) * suction * v_suction * (
            ratio ** ((K_ISENTROPIC - 1.0) / K_ISENTROPIC) - 1.0
        )
    else:
        ideal = None
    efficiency = (ideal / indicated_work) if ideal and indicated_work and indicated_work > 1e-9 else None

    # Phase boundaries by pressure crossing (aligned).
    high_mask = p >= threshold
    suction_frac = float(np.mean(~high_mask))
    discharge_frac = float(np.mean(high_mask))
    dp = np.diff(p)
    rising = dp[dp > 0]
    falling = dp[dp < 0]
    compression_slope = float(np.quantile(rising, 0.9)) if rising.size else 0.0
    expansion_slope = float(-np.quantile(falling, 0.1)) if falling.size else 0.0
    d2 = np.diff(p, 2)
    band = np.abs(p[1:-1] - threshold) <= 0.1 * span
    corner = float(np.mean(np.abs(d2[band]))) if band.any() else 0.0

    # Normalised phase profile (0..1) on the aligned cycle.
    idx = np.linspace(0, count - 1, PROFILE_BINS).round().astype(int)
    prof = (p[idx] - suction) / span if span else np.zeros(PROFILE_BINS)
    profile = {f"prof_{i:02d}": float(prof[i]) for i in range(PROFILE_BINS)}

    features = {
        "n_samples": count,
        "p_min": float(np.min(p)),
        "p_max": float(np.max(p)),
        "p_mean": float(np.mean(p)),
        "p_suction": suction,
        "p_discharge": discharge,
        "p_span": span,
        "pressure_ratio": ratio,
        "suction_frac": suction_frac,
        "discharge_frac": discharge_frac,
        "compression_slope": compression_slope,
        "expansion_slope": expansion_slope,
        "corner_sharpness": corner,
        "n_comp": _polytropic_index(aligned_volume, p, suction, discharge, rising=True),
        "n_exp": _polytropic_index(aligned_volume, p, suction, discharge, rising=False),
        "indicated_work": indicated_work,
        "ideal_work_isentropic": ideal,
        "efficiency": efficiency,
        "phase_offset": offset,
    }
    features.update(profile)
    return features


def normalized_profile(signal: np.ndarray, cylinder: str, bins: int = 180, smooth_window: int = 21) -> np.ndarray | None:
    """Process-invariant P-theta shape: PV-area-aligned, resampled, plateau-normalised.

    Returns (P - P_suction) / (P_discharge - P_suction) over one revolution, so
    level/load changes cancel and only the waveform shape remains.  Valve faults
    distort corners/tilts here even when scalar levels look normal.
    """
    signal = np.asarray(signal, dtype=float)
    signal = signal[np.isfinite(signal)]
    if len(signal) < 60:
        return None
    p = smooth(signal, smooth_window)
    p = np.where(np.isfinite(p), p, signal)
    suction, discharge, _ = _plateau_levels(p)
    span = discharge - suction
    if span <= 1e-9:
        return None
    volume = volume_from_angle(len(p), CYLINDER_BORE_MM.get(cylinder, 490.0))
    dvolume = np.diff(volume, append=volume[0])
    offset, _ = _align_offset(p, dvolume)
    rotated = np.roll(p, -offset)
    idx = np.linspace(0, len(rotated) - 1, bins)
    resampled = np.interp(idx, np.arange(len(rotated)), rotated)
    return (resampled - suction) / span