Chapter 64: Audio and Signal Competitions

import numpy as np
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_auc_score
from sklearn.model_selection import GroupKFold, KFold
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

from kaggle_companion.activities._common import clean

SEED = 64
REPLICATES = 5                  # independent constructed collections; every estimate is their mean
RECORDINGS, NEW_RECORDINGS, NEW_CHUNKS = 60, 100, 4
CHUNK, BANDS = 512, 24          # samples per chunk; log-spectrum bands used as features
CALL_AMP = 0.5                  # loudness of the call tone relative to unit-variance background noise
CALL_FREQ = {0: 0.16, 1: 0.24}  # call frequency of species 0 and 1, in cycles per sample
MODELS = ["flexible", "simple"]


def recording_background(rng):
    """Each recording has its own background: a smooth random frequency response (a stand-in for site and microphone)."""
    return rng.normal(0, 3.0, 4) / np.arange(1, 5) ** 0.5, rng.uniform(0, 2 * np.pi, 4)


def chunk_features(background, species, rng, n):
    """n chunks of one recording as log band powers. Background noise is shaped by the recording's response; each chunk
    carries the species' call tone with probability 0.5."""
    amp, phase = background
    freq = np.fft.rfftfreq(CHUNK)
    gain_db = sum(a * np.cos(2 * np.pi * (i + 1) * freq * 2 + p) for i, (a, p) in enumerate(zip(amp, phase)))
    noise = np.fft.irfft(np.fft.rfft(rng.normal(size=(n, CHUNK)), axis=1) * 10 ** (gain_db / 20), CHUNK, axis=1)
    t, window = np.arange(CHUNK), np.hanning(CHUNK)
    call_freq = CALL_FREQ[species] * (1 + rng.normal(0, 0.04, n))
    tone = CALL_AMP * np.sin(2 * np.pi * call_freq[:, None] * t + rng.uniform(0, 2 * np.pi, (n, 1))) * np.sqrt(window)
    signal = noise + (rng.random(n) < 0.5)[:, None] * tone
    power = np.abs(np.fft.rfft(signal * window, axis=1)) ** 2
    width = (CHUNK // 2) // BANDS
    return np.log(power[:, :BANDS * width].reshape(n, BANDS, width).mean(2) + 1e-9)


def collection(recordings, chunks, rng):
    """Recordings alternate between the two species; every chunk inherits its recording's species and id."""
    X, y, group = [], [], []
    for r in range(recordings):
        X.append(chunk_features(recording_background(rng), r % 2, rng, chunks))
        y += [r % 2] * chunks
        group += [r] * chunks
    return np.vstack(X), np.array(y), np.array(group)


def make_model(kind):
    if kind == "flexible":      # a high-capacity learner standing in for a CNN: it can memorize recordings
        return ExtraTreesClassifier(40, max_features=0.5, random_state=0, n_jobs=1)
    return make_pipeline(StandardScaler(), LogisticRegression(C=10, max_iter=300))


def pooled_auc(model_kind, X, y, splits):
    """Pool every held-out prediction and score one AUC over them."""
    oof = np.zeros(len(y))
    for fit, hold in splits:
        oof[hold] = make_model(model_kind).fit(X[fit], y[fit]).predict_proba(X[hold])[:, 1]
    return roc_auc_score(y, oof)


def run(chunks_per_recording):
    scores = {m: {"chunk_cv": [], "recording_cv": [], "new_recordings": []} for m in MODELS}
    for rep in range(REPLICATES):
        rng = np.random.default_rng(SEED * 100 + rep)
        X, y, group = collection(RECORDINGS, chunks_per_recording, rng)
        X_new, y_new, _ = collection(NEW_RECORDINGS, NEW_CHUNKS, rng)        # recordings the model never saw
        for m in MODELS:
            scores[m]["chunk_cv"].append(pooled_auc(m, X, y, KFold(5, shuffle=True, random_state=rep).split(X)))
            scores[m]["recording_cv"].append(pooled_auc(m, X, y, GroupKFold(5).split(X, y, group)))
            scores[m]["new_recordings"].append(roc_auc_score(y_new, make_model(m).fit(X, y).predict_proba(X_new)[:, 1]))
    return clean({
        "chunks_per_recording": chunks_per_recording,
        **{m: {k: float(np.mean(v)) for k, v in scores[m].items()} for m in MODELS},
        "per_replicate": scores,
        "recordings": RECORDINGS, "training_chunks": RECORDINGS * chunks_per_recording, "new_chunks": NEW_RECORDINGS * NEW_CHUNKS,
        "replicates": REPLICATES,
    })

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

VALUES = [1, 3, 8, 16]
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.ch64 import draw, NCOLS, HEIGHT

DEFAULT = 8
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()