Blame for src/phyphox_pipeline.py

"""Phyphox sensor pipeline for Human Activity Recognition.

Converts raw Phyphox accelerometer + gyroscope CSV exports into the
561-feature vector expected by the UCI HAR classifier.

Feature order matches features.txt exactly:
  1-200   : time-domain 3-axis signals  (5 signals × 40 features)
  201-265 : time-domain magnitudes      (5 signals × 13 features)
  266-502 : frequency-domain 3-axis     (3 signals × 79 features)
  503-554 : frequency-domain magnitudes (4 signals × 13 features)
  555-561 : angle features              (7 features)
"""

import io
import numpy as np
import pandas as pd
from scipy import signal as sp_signal
from scipy.stats import skew, kurtosis as sp_kurtosis

FS = 50        # target sampling rate Hz
WINDOW = 128   # samples per window (2.56 s)
STEP = 64      # hop size — 50% overlap
AR_ORDER = 4   # Burg AR model order

# 14 frequency band pairs (1-indexed, inclusive) applied per axis
BANDS = [
    (1, 8), (9, 16), (17, 24), (25, 32), (33, 40), (41, 48),
    (49, 56), (57, 64), (1, 16), (17, 32), (33, 48), (49, 64),
    (1, 24), (25, 48),
]

ACC_COLS = ["Time (s)", "X (m/s^2)", "Y (m/s^2)", "Z (m/s^2)"]
GYRO_COLS = ["Time (s)", "X (rad/s)", "Y (rad/s)", "Z (rad/s)"]


def _parse_csv(file_obj, expected_cols: list) -> pd.DataFrame:
    """Parse a Phyphox CSV export and validate required columns.

    Args:
        file_obj: file-like object (bytes or str) from Phyphox export
        expected_cols: list of required column names

    Returns:
        DataFrame with numeric data, NaN rows dropped

    Raises:
        ValueError: if columns are missing or file cannot be parsed
    """
    try:
        raw = file_obj.read()
        if isinstance(raw, bytes):
            raw = raw.decode("utf-8")
        df = pd.read_csv(io.StringIO(raw), float_precision="high")
    except Exception as exc:
        raise ValueError(f"Cannot parse CSV: {exc}") from exc

    df.columns = [c.strip('"').strip() for c in df.columns]
    missing = [c for c in expected_cols if c not in df.columns]
    if missing:
        raise ValueError(
            f"Missing columns {missing}. Found: {list(df.columns)}. "
            "Check you uploaded the correct Phyphox CSV (Accelerometer or Gyroscope)."
        )
    return df[expected_cols].apply(pd.to_numeric, errors="coerce").dropna()


def _butter_lp(data: np.ndarray, cutoff: float, fs: float = FS, order: int = 3) -> np.ndarray:
    """Zero-phase Butterworth low-pass filter applied along axis 0.

    Args:
        data: 1-D or 2-D array
        cutoff: cutoff frequency in Hz
        fs: sampling rate in Hz
        order: filter order

    Returns:
        Filtered array, same shape as input
    """
    b, a = sp_signal.butter(order, cutoff / (fs / 2.0), btype="low")
    # Use maximum safe padding — critical for low cutoffs (e.g. 0.3 Hz needs
    # ~167 samples to settle; scipy's default 9-sample pad is far too short).
    padlen = min(len(data) - 1, max(3 * int(fs / cutoff), 9))
    if data.ndim == 1:
        return sp_signal.filtfilt(b, a, data, padlen=padlen)
    return np.column_stack(
        [sp_signal.filtfilt(b, a, data[:, i], padlen=padlen) for i in range(data.shape[1])]
    )


def _median_filt(data: np.ndarray, k: int = 3) -> np.ndarray:
    """Median filter applied along axis 0.

    Args:
        data: 1-D or 2-D array
        k: kernel size (must be odd)

    Returns:
        Filtered array, same shape as input
    """
    if data.ndim == 1:
        return sp_signal.medfilt(data, kernel_size=k)
    return np.column_stack(
        [sp_signal.medfilt(data[:, i], kernel_size=k) for i in range(data.shape[1])]
    )


def _burg_ar(x: np.ndarray, order: int = AR_ORDER) -> np.ndarray:
    """Burg method autoregressive coefficients.

    Implements the standard Burg recursion with Levinson-Durbin update.
    Mean-centres the signal before fitting.

    Args:
        x: 1-D signal array
        order: AR model order

    Returns:
        Array of `order` AR coefficients [a1, a2, ..., ap]
    """
    x = np.asarray(x, dtype=np.float64)
    x = x - x.mean()
    N = len(x)
    ef = x.copy()
    eb = x.copy()
    a = np.zeros(order)

    for m in range(1, order + 1):
        f = ef[m:].copy()
        b = eb[m - 1: N - 1].copy()

        denom = np.dot(f, f) + np.dot(b, b) + 1e-12
        km = -2.0 * np.dot(f, b) / denom

        # Levinson-Durbin update of AR polynomial
        a_prev = a[:m - 1].copy()
        for j in range(m - 1):
            a[j] = a_prev[j] + km * a_prev[m - 2 - j]
        a[m - 1] = km

        # Update forward/backward prediction errors
        ef[m:] = f + km * b
        eb[m - 1: N - 1] = b + km * f

    return a


def _entropy(x: np.ndarray) -> float:
    """Normalised signal entropy via absolute-value probability distribution.

    Args:
        x: 1-D array

    Returns:
        Entropy value (>= 0)
    """
    total = np.abs(x).sum()
    if total < 1e-12:
        return 0.0
    p = np.abs(x) / total
    p = p[p > 0]
    return float(-np.sum(p * np.log(p)))


def _bands_energy(fft_mag: np.ndarray) -> np.ndarray:
    """Energy in each of the 14 UCI HAR frequency bands (1-indexed, inclusive).

    Args:
        fft_mag: FFT magnitude array, must contain at least 64 values

    Returns:
        Array of 14 energy values
    """
    m = fft_mag[:64]
    return np.array([float(np.sum(m[s - 1: e] ** 2)) for s, e in BANDS])


def _safe_corr(a: np.ndarray, b: np.ndarray) -> float:
    """Pearson correlation, returns 0.0 if either signal is constant.

    Args:
        a: first 1-D array
        b: second 1-D array

    Returns:
        Correlation coefficient in [-1, 1]
    """
    if a.std() < 1e-10 or b.std() < 1e-10:
        return 0.0
    r = np.corrcoef(a, b)[0, 1]
    return 0.0 if not np.isfinite(r) else float(r)


def _angle(u: np.ndarray, v: np.ndarray) -> float:
    """Angle in radians between two 3-D vectors.

    Args:
        u: first vector, shape (3,)
        v: second vector, shape (3,)

    Returns:
        Angle in radians, or 0.0 if either vector is zero
    """
    un, vn = np.linalg.norm(u), np.linalg.norm(v)
    if un < 1e-10 or vn < 1e-10:
        return 0.0
    return float(np.arccos(np.clip(np.dot(u, v) / (un * vn), -1.0, 1.0)))


def _t3ax(sig: np.ndarray) -> np.ndarray:
    """40 time-domain features from a 3-axis signal (N, 3).

    Order: mean×3, std×3, mad×3, max×3, min×3, sma,
           energy×3, iqr×3, entropy×3, arCoeff×12, correlation×3
    """
    N = len(sig)
    x, y, z = sig[:, 0], sig[:, 1], sig[:, 2]
    out = []

    out += [x.mean(), y.mean(), z.mean()]
    out += [x.std(), y.std(), z.std()]
    out += [
        float(np.median(np.abs(x - np.median(x)))),
        float(np.median(np.abs(y - np.median(y)))),
        float(np.median(np.abs(z - np.median(z)))),
    ]
    out += [x.max(), y.max(), z.max()]
    out += [x.min(), y.min(), z.min()]
    out += [float((np.abs(x) + np.abs(y) + np.abs(z)).sum() / N)]   # sma
    out += [float(np.sum(x ** 2) / N), float(np.sum(y ** 2) / N), float(np.sum(z ** 2) / N)]
    out += [
        float(np.percentile(x, 75) - np.percentile(x, 25)),
        float(np.percentile(y, 75) - np.percentile(y, 25)),
        float(np.percentile(z, 75) - np.percentile(z, 25)),
    ]
    out += [_entropy(x), _entropy(y), _entropy(z)]
    out += _burg_ar(x).tolist()
    out += _burg_ar(y).tolist()
    out += _burg_ar(z).tolist()
    out += [_safe_corr(x, y), _safe_corr(x, z), _safe_corr(y, z)]

    return np.array(out, dtype=np.float64)   # 40 values


def _tmag(sig: np.ndarray) -> np.ndarray:
    """13 time-domain features from a 1-D magnitude signal.

    Order: mean, std, mad, max, min, sma, energy, iqr, entropy, arCoeff×4
    """
    N = len(sig)
    out = [
        float(sig.mean()),
        float(sig.std()),
        float(np.median(np.abs(sig - np.median(sig)))),
        float(sig.max()),
        float(sig.min()),
        float(np.abs(sig).sum() / N),            # sma (1-D)
        float(np.sum(sig ** 2) / N),             # energy
        float(np.percentile(sig, 75) - np.percentile(sig, 25)),
        _entropy(sig),
    ]
    out += _burg_ar(sig).tolist()
    return np.array(out, dtype=np.float64)   # 13 values


def _f3ax(sig: np.ndarray) -> np.ndarray:
    """79 frequency-domain features from a 3-axis signal (N, 3).

    Order: mean×3, std×3, mad×3, max×3, min×3, sma,
           energy×3, iqr×3, entropy×3,
           maxInds×3, meanFreq×3,
           (skewness, kurtosis)×3 interleaved,
           bandsEnergy×14 per axis (×3 axes = 42)
    """
    x, y, z = sig[:, 0], sig[:, 1], sig[:, 2]

    def _fft(s):
        return np.abs(np.fft.rfft(s))[:64]

    fx, fy, fz = _fft(x), _fft(y), _fft(z)
    bins = np.arange(1, 65, dtype=np.float64)   # 1-indexed bin numbers

    def _mfreq(fm):
        t = fm.sum()
        return float(np.dot(bins[:len(fm)], fm) / t) if t > 1e-12 else 0.0

    def _maxinds(fm):
        return float(np.argmax(fm) + 1)   # 1-indexed

    out = []
    out += [fx.mean(), fy.mean(), fz.mean()]
    out += [fx.std(), fy.std(), fz.std()]
    out += [
        float(np.median(np.abs(fx - np.median(fx)))),
        float(np.median(np.abs(fy - np.median(fy)))),
        float(np.median(np.abs(fz - np.median(fz)))),
    ]
    out += [fx.max(), fy.max(), fz.max()]
    out += [fx.min(), fy.min(), fz.min()]
    n = len(fx)
    out += [float((fx + fy + fz).sum() / n)]    # sma of FFT mags
    out += [float(np.sum(fx ** 2) / n), float(np.sum(fy ** 2) / n), float(np.sum(fz ** 2) / n)]
    out += [
        float(np.percentile(fx, 75) - np.percentile(fx, 25)),
        float(np.percentile(fy, 75) - np.percentile(fy, 25)),
        float(np.percentile(fz, 75) - np.percentile(fz, 25)),
    ]
    out += [_entropy(fx), _entropy(fy), _entropy(fz)]
    out += [_maxinds(fx), _maxinds(fy), _maxinds(fz)]
    out += [_mfreq(fx), _mfreq(fy), _mfreq(fz)]
    # skewness/kurtosis interleaved per axis (skX,kurX, skY,kurY, skZ,kurZ)
    out += [float(skew(fx)), float(sp_kurtosis(fx))]
    out += [float(skew(fy)), float(sp_kurtosis(fy))]
    out += [float(skew(fz)), float(sp_kurtosis(fz))]
    out += _bands_energy(fx).tolist()
    out += _bands_energy(fy).tolist()
    out += _bands_energy(fz).tolist()

    return np.array(out, dtype=np.float64)   # 79 values


def _fmag(sig: np.ndarray) -> np.ndarray:
    """13 frequency-domain features from a 1-D magnitude signal.

    Order: mean, std, mad, max, min, sma, energy, iqr, entropy,
           maxInds, meanFreq, skewness, kurtosis
    """
    fm = np.abs(np.fft.rfft(sig))[:64]
    n = len(fm)
    bins = np.arange(1, n + 1, dtype=np.float64)
    total = fm.sum()
    out = [
        float(fm.mean()),
        float(fm.std()),
        float(np.median(np.abs(fm - np.median(fm)))),
        float(fm.max()),
        float(fm.min()),
        float(np.abs(fm).sum() / n),
        float(np.sum(fm ** 2) / n),
        float(np.percentile(fm, 75) - np.percentile(fm, 25)),
        _entropy(fm),
        float(np.argmax(fm) + 1),                                        # maxInds (1-indexed)
        float(np.dot(bins, fm) / total) if total > 1e-12 else 0.0,      # meanFreq
        float(skew(fm)),
        float(sp_kurtosis(fm)),
    ]
    return np.array(out, dtype=np.float64)   # 13 values


def _window_features(
    body_acc: np.ndarray,
    grav_acc: np.ndarray,
    body_jerk: np.ndarray,
    gyro: np.ndarray,
    gyro_jerk: np.ndarray,
) -> np.ndarray:
    """Extract all 561 features from one pre-processed window.

    Args:
        body_acc:  body linear acceleration (128, 3) m/s²
        grav_acc:  gravity component        (128, 3) m/s²
        body_jerk: body jerk                (127, 3) m/s³
        gyro:      angular velocity         (128, 3) rad/s
        gyro_jerk: gyro jerk                (127, 3) rad/s²

    Returns:
        1-D array of 561 features
    """
    # Magnitudes
    ba_mag = np.linalg.norm(body_acc, axis=1)
    ga_mag = np.linalg.norm(grav_acc, axis=1)
    bj_mag = np.linalg.norm(body_jerk, axis=1)
    gy_mag = np.linalg.norm(gyro, axis=1)
    gj_mag = np.linalg.norm(gyro_jerk, axis=1)

    parts = []

    # Time 3-axis (5 × 40 = 200)
    for sig in [body_acc, grav_acc, body_jerk, gyro, gyro_jerk]:
        parts.append(_t3ax(sig))

    # Time magnitudes (5 × 13 = 65)
    for mag in [ba_mag, ga_mag, bj_mag, gy_mag, gj_mag]:
        parts.append(_tmag(mag))

    # Freq 3-axis (3 × 79 = 237)
    for sig in [body_acc, body_jerk, gyro]:
        parts.append(_f3ax(sig))

    # Freq magnitudes (4 × 13 = 52)
    for mag in [ba_mag, bj_mag, gy_mag, gj_mag]:
        parts.append(_fmag(mag))

    # Angle features (7)
    ba_mean = body_acc.mean(axis=0)
    ga_mean = grav_acc.mean(axis=0)
    bj_mean = body_jerk.mean(axis=0)
    gy_mean = gyro.mean(axis=0)
    gj_mean = gyro_jerk.mean(axis=0)

    parts.append(np.array([
        _angle(ba_mean, ga_mean),
        _angle(bj_mean, ga_mean),
        _angle(gy_mean, ga_mean),
        _angle(gj_mean, ga_mean),
        _angle(np.array([1.0, 0.0, 0.0]), ga_mean),
        _angle(np.array([0.0, 1.0, 0.0]), ga_mean),
        _angle(np.array([0.0, 0.0, 1.0]), ga_mean),
    ]))

    result = np.concatenate(parts)
    assert result.shape == (561,), f"Feature count error: got {result.shape[0]}, expected 561"
    return result


def process_phyphox_files(
    acc_file,
    gyro_file,
) -> tuple:
    """Convert Phyphox CSV exports to (n_windows, 561) feature array.

    Pipeline:
      1. Parse + validate both CSVs
      2. Convert accelerometer from m/s² to g (÷ 9.80665) to match UCI training units
      3. Interpolate onto common 50 Hz grid
      4. Segment: 128-sample windows, 64-sample hop (50% overlap)
      5. Per window: median filter → 20 Hz Butterworth → gravity separation
         at 0.3 Hz → jerk → magnitudes → 561 features

    Args:
        acc_file: file-like object — Phyphox Accelerometer CSV
                  (columns: Time (s), X (m/s^2), Y (m/s^2), Z (m/s^2))
        gyro_file: file-like object — Phyphox Gyroscope CSV
                   (columns: Time (s), X (rad/s), Y (rad/s), Z (rad/s))

    Returns:
        Tuple of:
          np.ndarray shape (n_windows, 561) — raw (un-normalised) features
          list[str]                          — warning messages

    Raises:
        ValueError: invalid format, wrong columns, or < 3 s of data
    """
    warnings: list = []

    acc_df = _parse_csv(acc_file, ACC_COLS)
    gyro_df = _parse_csv(gyro_file, GYRO_COLS)

    acc_t = acc_df["Time (s)"].values
    acc_xyz = acc_df[["X (m/s^2)", "Y (m/s^2)", "Z (m/s^2)"]].values / 9.80665  # m/s² → g
    gyro_t = gyro_df["Time (s)"].values
    gyro_xyz = gyro_df[["X (rad/s)", "Y (rad/s)", "Z (rad/s)"]].values

    t0 = max(acc_t[0], gyro_t[0])
    t1 = min(acc_t[-1], gyro_t[-1])
    duration = t1 - t0

    if duration < 3.0:
        raise ValueError(
            f"Recording is {duration:.2f} s — minimum 3 seconds required. "
            "Hold the phone still or walk for at least 3 seconds before exporting."
        )

    t_grid = np.arange(t0, t1, 1.0 / FS)
    am = (acc_t >= t0) & (acc_t <= t1)
    gm = (gyro_t >= t0) & (gyro_t <= t1)

    acc_50 = np.column_stack(
        [np.interp(t_grid, acc_t[am], acc_xyz[am, i]) for i in range(3)]
    )
    gyro_50 = np.column_stack(
        [np.interp(t_grid, gyro_t[gm], gyro_xyz[gm, i]) for i in range(3)]
    )

    n = min(len(acc_50), len(gyro_50))
    acc_50, gyro_50 = acc_50[:n], gyro_50[:n]

    n_windows = max(0, (n - WINDOW) // STEP + 1)
    if n_windows == 0:
        raise ValueError(
            f"Only {n} samples ({n / FS:.1f} s) after alignment — "
            f"need at least {WINDOW} samples ({WINDOW / FS:.1f} s)."
        )

    if duration > 60:
        warnings.append(f"Long recording ({duration:.0f} s) — {n_windows} windows extracted.")

    # Apply noise filters to the full signal before windowing.
    acc_50  = _butter_lp(_median_filt(acc_50), cutoff=20.0)
    gyro_50 = _butter_lp(_median_filt(gyro_50), cutoff=20.0)

    all_features = []
    dt = 1.0 / FS

    for start in range(0, n - WINDOW + 1, STEP):
        end = start + WINDOW
        aw = acc_50[start:end]   # (128, 3)
        gw = gyro_50[start:end]  # (128, 3)

        # Gravity separation: window mean as gravity estimate.
        # The UCI pipeline used a 0.3 Hz LP on a full continuous recording.
        # That filter needs ~167 samples to settle; a 10-second clip gives only
        # ~500 samples total, so the filter corrupts all but the central window.
        # The window mean is equivalent for symmetric activities (oscillations
        # cancel over a stride cycle) and exact for static activities.
        grav = np.tile(aw.mean(axis=0), (WINDOW, 1))
        body = aw - grav

        # Jerk: finite difference → (127, 3)
        body_jerk = np.diff(body, axis=0) / dt
        gyro_jerk = np.diff(gw, axis=0) / dt

        all_features.append(_window_features(body, grav, body_jerk, gw, gyro_jerk))

    return np.array(all_features), warnings