Chapter 23: Sensors and HAR: The Cross-Subject Problem

import numpy as np
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.metrics import f1_score
from sklearn.model_selection import LeaveOneGroupOut, StratifiedKFold

from kaggle_companion.activities._common import clean

SEED = 23
DATASETS = 4                 # independent constructed cohorts; every estimate is their mean
DEV_SUBJECTS, NEW_SUBJECTS, WINDOWS = 8, 12, 24   # windows per subject and activity
HZ, T, STEP = 50, 100, 50    # 50 Hz sampling, 2-second windows, 50% overlap
# activity: (cadence in Hz, amplitude, vertical share, rhythm irregularity, gravity tilt in radians)
ACTIVITIES = {"walking": (1.8, 1.0, 0.8, 0.1, 0.0), "running": (2.8, 1.7, 0.8, 0.1, 0.0), "cycling": (1.4, 0.8, 0.3, 0.05, 0.0),
              "stairs": (1.6, 1.1, 0.8, 0.7, 0.0), "sitting": (0.0, 0.0, 0.0, 0.0, 0.0), "standing": (0.0, 0.0, 0.0, 0.0, 0.35)}
NAMES = list(ACTIVITIES)


def rotation(angle):
    c, s = np.cos(angle), np.sin(angle)
    return np.array([[1, 0, 0], [0, c, -s], [0, s, c]])


def make_subject(distinctiveness, rng):
    """A person: a cadence multiplier, an amplitude multiplier, a wrist orientation and a sensor bias."""
    d = distinctiveness
    return {"cadence": 1 + 0.15 * d * rng.normal(), "amplitude": np.exp(0.35 * d * rng.normal()),
            "rotation": rotation(0.45 * d * rng.normal()), "bias": 0.3 * d * rng.normal(size=3)}


def recording(activity, person, rng):
    """One continuous 6-channel recording (3 accelerometer + 3 gyroscope axes), cut into 50%-overlapping windows."""
    cadence, amp, vertical, irregular, tilt = ACTIVITIES[activity]
    length = (WINDOWS - 1) * STEP + T
    t = np.arange(length) / HZ
    f, a = cadence * person["cadence"], amp * person["amplitude"]
    phase = rng.uniform(0, 2 * np.pi) + np.cumsum(rng.normal(0, irregular * 0.15, length))   # irregular rhythm drifts the phase
    wave = np.sin(2 * np.pi * f * t + phase)
    acc = np.column_stack([a * (1 - vertical) * np.sin(2 * np.pi * f * t + phase + 1.0), 0.4 * a * wave ** 2, a * vertical * wave])
    acc = (acc + rotation(tilt) @ np.array([0, 0, 1.0])) @ person["rotation"].T + person["bias"]   # gravity, then the person's orientation
    gyro = 0.5 * a * np.column_stack([np.cos(2 * np.pi * f * t + phase), np.sin(2 * np.pi * f * t + phase), wave]) @ person["rotation"].T
    signal = np.column_stack([acc, gyro]) + rng.normal(0, 0.8, (length, 6))
    return np.stack([signal[i * STEP:i * STEP + T] for i in range(WINDOWS)])


def window_features(w):
    """Per channel: mean, standard deviation, energy, zero-crossing rate, dominant frequency and its magnitude; plus axis correlations."""
    mean, sd, energy = w.mean(1), w.std(1), (w ** 2).mean(1)
    zero_cross = ((w[:, 1:] - mean[:, None]) * (w[:, :-1] - mean[:, None]) < 0).mean(1)
    spectrum = np.abs(np.fft.rfft(w - mean[:, None], axis=1))
    peak, magnitude = spectrum[:, 1:].argmax(1) + 1, spectrum[:, 1:].max(1)
    corr = [np.einsum("nt,nt->n", w[:, :, i] - mean[:, None, i], w[:, :, j] - mean[:, None, j]) / (T * sd[:, i] * sd[:, j] + 1e-6)
            for i, j in ((0, 1), (0, 2), (1, 2))]
    return np.column_stack([mean, sd, energy, zero_cross, peak * HZ / T, magnitude / T] + [c[:, None] for c in corr])


def cohort(n_subjects, distinctiveness, rng):
    X, y, who = [], [], []
    for s in range(n_subjects):
        person = make_subject(distinctiveness, rng)
        for k, name in enumerate(NAMES):
            X.append(window_features(recording(name, person, rng)))
            y += [k] * WINDOWS
            who += [s] * WINDOWS
    return np.vstack(X), np.array(y), np.array(who)


def model():
    return ExtraTreesClassifier(n_estimators=60, random_state=0, n_jobs=1)


def macro_f1(y, pred, labels=None):
    return f1_score(y, pred, labels=labels, average="macro")


def one_cohort(distinctiveness, seed):
    rng = np.random.default_rng(seed)
    X, y, who = cohort(DEV_SUBJECTS, distinctiveness, rng)
    X_new, y_new, _ = cohort(NEW_SUBJECTS, distinctiveness, rng)      # different people, never seen in development
    random_pred, subject_pred = np.zeros(len(y), int), np.zeros(len(y), int)
    for fit, val in StratifiedKFold(5, shuffle=True, random_state=0).split(X, y):      # the chapter's mistake
        random_pred[val] = model().fit(X[fit], y[fit]).predict(X[val])
    for fit, val in LeaveOneGroupOut().split(X, y, who):                               # the chapter's fix
        subject_pred[val] = model().fit(X[fit], y[fit]).predict(X[val])
    new_pred = model().fit(X, y).predict(X_new)
    seated = [NAMES.index("sitting"), NAMES.index("standing")]
    return [macro_f1(y, random_pred), macro_f1(y, subject_pred), macro_f1(y_new, new_pred),
            macro_f1(y, random_pred, seated), macro_f1(y, subject_pred, seated), macro_f1(y_new, new_pred, seated)]


def run(distinctiveness):
    scores = np.array([one_cohort(distinctiveness, SEED * 1000 + i) for i in range(DATASETS)])
    mean = scores.mean(0)
    return clean({
        "distinctiveness": distinctiveness, "datasets": DATASETS, "dev_subjects": DEV_SUBJECTS, "new_subjects": NEW_SUBJECTS,
        "random": mean[0], "loso": mean[1], "new": mean[2],
        "sit_stand": {"random": mean[3], "loso": mean[4], "new": mean[5]},
        "per_dataset": {"random": scores[:, 0], "loso": scores[:, 1], "new": scores[:, 2]},
    })

import pandas as pd
from kaggle_companion.activities.ch23 import explain

VALUES = [0, 0.5, 1.0, 1.5]
results = {value: run(value) for value in VALUES}
table = pd.DataFrame({str(v): explain(results[v], v)['metrics'] for v in VALUES})
print(table.to_string())

import matplotlib.pyplot as plt
from kaggle_companion.activities.ch23 import draw, NCOLS, HEIGHT

DEFAULT = 1.0
fig, axes = plt.subplots(1, NCOLS, figsize=(6.6 if NCOLS == 1 else 10.4, HEIGHT), layout='constrained')
draw(axes, results[DEFAULT], DEFAULT)
plt.show()