センサーデータの前処理と特徴量設計｜Pythonで欠損・外れ値・時系列特徴量を扱う

"""Beginner tutorial: turn long-format sensor CSV into leakage-aware features."""
from __future__ import annotations

from pathlib import Path
import argparse
import json
import platform

import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import scipy
from scipy.signal import butter, sosfiltfilt

REQUIRED = {"timestamp", "sensor", "value", "unit"}
CONFIG = {
    "temperature": {"unit": "degC", "cutoff_hz": 0.010},
    "pressure": {"unit": "kPa", "cutoff_hz": 0.030},
    "motor_current": {"unit": "A", "cutoff_hz": 0.050},
}


def make_example(path: Path) -> pd.DataFrame:
    """Create editorial synthetic data; values are not a process recipe."""
    rng = np.random.default_rng(42)
    start = pd.Timestamp("2026-01-01 09:00:00", tz="UTC")
    seconds = np.arange(0, 2 * 60 * 60 + 1, 5)
    t = seconds / 60
    rows: list[tuple] = []
    signals = {
        "temperature": (45 + 0.025 * t + 3.0 * (t >= 55) + 0.35 * np.sin(t / 7), "degC", 0.20),
        "pressure": (210 + 4.5 * np.sin(t / 9) + 7.0 * (t >= 75), "kPa", 0.80),
        "motor_current": (8.0 + 0.35 * np.sin(t / 3) + 1.1 * ((t >= 85) & (t < 95)), "A", 0.12),
    }
    for sensor, (truth, unit, noise) in signals.items():
        values = truth + rng.normal(0, noise, len(t))
        for i in [330, 780, 1120]:
            values[i] += {"temperature": 8, "pressure": -30, "motor_current": 4}[sensor]
        for i, (sec, value) in enumerate(zip(seconds, values)):
            # Short gaps are intentionally inserted for interpolation practice.
            if (sensor == "temperature" and 500 <= i <= 502) or (sensor == "pressure" and 900 <= i <= 901):
                continue
            rows.append((start + pd.Timedelta(seconds=int(sec)), sensor, float(value), unit))
    frame = pd.DataFrame(rows, columns=["timestamp", "sensor", "value", "unit"])
    frame.to_csv(path, index=False)
    return frame


def validate(frame: pd.DataFrame) -> pd.DataFrame:
    if not REQUIRED <= set(frame.columns):
        raise ValueError(f"Required columns: {sorted(REQUIRED)}")
    out = frame[list(REQUIRED)].copy()
    if out[["timestamp", "sensor", "unit"]].isna().any().any():
        raise ValueError("Identity, timestamp, or unit is missing")
    out["timestamp"] = pd.to_datetime(out["timestamp"], utc=True, errors="raise")
    out["value"] = pd.to_numeric(out["value"], errors="raise")
    if not np.isfinite(out["value"]).all():
        raise ValueError("Sensor values must be finite")
    if out.duplicated(["sensor", "timestamp"]).any():
        raise ValueError("Duplicate timestamp within a sensor")
    unknown = set(out["sensor"]) - set(CONFIG)
    if unknown:
        raise ValueError(f"Unknown sensors; add reviewed settings: {sorted(unknown)}")
    for sensor, group in out.groupby("sensor"):
        units = set(group["unit"])
        if units != {CONFIG[sensor]["unit"]}:
            raise ValueError(f"Unit mismatch for {sensor}: {sorted(units)}")
        if len(group) < 100:
            raise ValueError(f"Too few observations for {sensor}")
    return out.sort_values(["sensor", "timestamp"]).reset_index(drop=True)


def hampel(series: pd.Series, window: int = 25, n_sigma: float = 4.0):
    if window < 5 or window % 2 == 0 or n_sigma <= 0:
        raise ValueError("Hampel window must be odd >=5 and n_sigma >0")
    median = series.rolling(window, center=True, min_periods=window // 2).median()
    deviation = (series - median).abs()
    mad = deviation.rolling(window, center=True, min_periods=window // 2).median()
    threshold = n_sigma * 1.4826 * mad
    mask = (deviation > threshold) & threshold.gt(0) & series.notna()
    cleaned = series.mask(mask, median)
    return cleaned, mask


def filter_segments(series: pd.Series, fs_hz: float, cutoff_hz: float, order: int = 3) -> pd.Series:
    if not 0 < cutoff_hz < fs_hz / 2:
        raise ValueError("cutoff_hz must be between 0 and the Nyquist frequency")
    sos = butter(order, cutoff_hz, btype="low", fs=fs_hz, output="sos")
    result = pd.Series(np.nan, index=series.index, dtype=float)
    block = series.notna().ne(series.notna().shift()).cumsum()
    for _, part in series.dropna().groupby(block[series.notna()]):
        if len(part) >= 30:
            result.loc[part.index] = sosfiltfilt(sos, part.to_numpy())
    return result


def preprocess(frame: pd.DataFrame, freq="5s", max_interp="15s"):
    frame = validate(frame)
    starts = frame.groupby("sensor")["timestamp"].min()
    ends = frame.groupby("sensor")["timestamp"].max()
    start, end = starts.max(), ends.min()
    grid = pd.date_range(start.ceil(freq), end.floor(freq), freq=freq)
    if len(grid) < 100:
        raise ValueError("Common time range is too short")
    step_s = pd.Timedelta(freq).total_seconds()
    max_points = int(pd.Timedelta(max_interp).total_seconds() // step_s)
    if max_points < 1:
        raise ValueError("max_interp must be at least one grid interval")
    rows, quality = [], []
    for sensor, group in frame.groupby("sensor", sort=True):
        observed = group.set_index("timestamp")["value"].sort_index()
        union = observed.index.union(grid)
        regular = observed.reindex(union).sort_index().interpolate(
            method="time", limit=max_points, limit_area="inside"
        ).reindex(grid)
        cleaned, spike = hampel(regular)
        filtered = filter_segments(cleaned, 1.0 / step_s, CONFIG[sensor]["cutoff_hz"])
        sensor_rows = pd.DataFrame({
            "timestamp": grid, "sensor": sensor,
            "regularized": regular.to_numpy(), "cleaned": cleaned.to_numpy(),
            "filtered": filtered.to_numpy(), "spike_flag": spike.to_numpy(),
            "unit": CONFIG[sensor]["unit"],
        })
        rows.append(sensor_rows)
        quality.append({
            "sensor": sensor, "observations": len(group), "grid_points": len(grid),
            "interpolated_or_missing_points": int((~grid.isin(observed.index)).sum()),
            "remaining_missing_points": int(regular.isna().sum()),
            "spikes_replaced": int(spike.sum()), "cutoff_hz": CONFIG[sensor]["cutoff_hz"],
        })
    return pd.concat(rows, ignore_index=True), pd.DataFrame(quality)


def feature_table(long: pd.DataFrame, window="5min", step="1min") -> pd.DataFrame:
    wide = long.pivot(index="timestamp", columns="sensor", values="filtered").sort_index()
    window_n = int(pd.Timedelta(window) / (wide.index[1] - wide.index[0]))
    step_n = int(pd.Timedelta(step) / (wide.index[1] - wide.index[0]))
    if window_n < 3 or step_n < 1:
        raise ValueError("Window/step is too short")
    records = []
    # End-aligned, past-only windows: [end-window, end].
    for stop in range(window_n - 1, len(wide), step_n):
        part = wide.iloc[stop - window_n + 1:stop + 1]
        if part.isna().any().any():
            continue
        record = {"window_end": part.index[-1], "window_start": part.index[0]}
        x = np.arange(len(part), dtype=float)
        for sensor in wide.columns:
            y = part[sensor].to_numpy()
            record.update({
                f"{sensor}_mean": float(y.mean()), f"{sensor}_std": float(y.std(ddof=1)),
                f"{sensor}_min": float(y.min()), f"{sensor}_max": float(y.max()),
                f"{sensor}_slope_per_sample": float(np.polyfit(x, y, 1)[0]),
                f"{sensor}_delta": float(y[-1] - y[0]),
                f"{sensor}_rms": float(np.sqrt(np.mean(y ** 2))),
            })
        records.append(record)
    if not records:
        raise ValueError("No complete windows; inspect missing intervals")
    return pd.DataFrame(records)


def make_plots(long: pd.DataFrame, features: pd.DataFrame, out: Path):
    temp = long[long.sensor == "temperature"]
    fig, axes = plt.subplots(2, 1, figsize=(10, 6), sharex=True)
    axes[0].plot(temp.timestamp, temp.regularized, color="0.6", lw=1, label="regularized")
    axes[0].plot(temp.timestamp, temp.filtered, color="#087e8b", lw=1.6, label="filtered")
    marked = temp[temp.spike_flag]
    axes[0].scatter(marked.timestamp, marked.regularized, color="#d1495b", s=28, label="spike flag")
    axes[0].set(ylabel="Temperature / degC", title="Editorial synthetic sensor data")
    axes[0].legend(ncol=3, fontsize=8)
    axes[1].plot(features.window_end, features.temperature_mean, label="5 min trailing mean")
    axes[1].plot(features.window_end, features.temperature_std, label="5 min trailing SD")
    axes[1].set(xlabel="Time / UTC", ylabel="Feature value", title="Past-only window features")
    axes[1].legend(fontsize=8)
    fig.autofmt_xdate(); fig.tight_layout(); fig.savefig(out / "sensor_preprocessing.png", dpi=160)
    plt.close(fig)


def main():
    parser = argparse.ArgumentParser()
    parser.add_argument("--input", help="Long-format CSV; omit to create a synthetic example")
    parser.add_argument("--output", default="sensor_results")
    args = parser.parse_args()
    out = Path(args.output); out.mkdir(parents=True, exist_ok=True)
    if args.input:
        raw = pd.read_csv(args.input); mode = "user CSV"
    else:
        raw = make_example(out / "sensor_example.csv"); mode = "editorial synthetic data"
    long, quality = preprocess(raw)
    features = feature_table(long)
    long.to_csv(out / "regularized_signals.csv", index=False)
    quality.to_csv(out / "data_quality.csv", index=False)
    features.to_csv(out / "sensor_features.csv", index=False)
    make_plots(long, features, out)
    summary = {
        "mode": mode, "input_rows": len(raw), "regularized_rows": len(long),
        "feature_windows": len(features), "spikes_replaced": int(quality.spikes_replaced.sum()),
        "remaining_missing_points": int(quality.remaining_missing_points.sum()),
        "versions": {"python": platform.python_version(), "numpy": np.__version__,
                     "pandas": pd.__version__, "scipy": scipy.__version__},
    }
    (out / "summary.json").write_text(json.dumps(summary, ensure_ascii=False, indent=2))
    print(json.dumps(summary, ensure_ascii=False, indent=2))


if __name__ == "__main__":
    main()