Chapter 17: Stacking and Hill-Climbing

import numpy as np
from sklearn.linear_model import LogisticRegression

from kaggle_companion.activities._common import clean

SEED = 17
REPLICATES = 25                    # independent constructed datasets; every estimate is their mean
FEATURES, FAMILIES = 20, 4         # 20 features, 6 of them informative; 4 genuinely different model families
N_TRAIN, N_DEV, N_ASSESS = 400, 200, 5000


def generate(rng, n, coef):
    X = rng.normal(size=(n, FEATURES))
    y = (rng.random(n) < 1 / (1 + np.exp(-(X @ coef)))).astype(int)
    return X, y


def greedy_brier_weights(y_development, candidate_predictions, step=0.1, max_iterations=30):
    """The chapter's hill-climbing function: start from the best candidate, then repeatedly move `step` of the weight
    to whichever candidate lowers the development Brier loss most, and stop when nothing helps."""
    y = np.asarray(y_development, dtype=float)
    p = np.asarray(candidate_predictions, dtype=float)
    if p.ndim != 2 or p.shape[0] != len(y) or p.shape[1] == 0:
        raise ValueError('rows by candidate matrix required')
    if not np.isfinite(p).all() or ((p < 0) | (p > 1)).any():
        raise ValueError('finite binary probabilities required')
    if not np.isin(y, [0, 1]).all() or not 0 < step <= 1:
        raise ValueError('binary labels and step in (0,1] required')
    loss = lambda q: np.mean((q - y) ** 2)
    start = int(np.argmin(np.mean((p - y[:, None]) ** 2, axis=0)))
    weights = np.eye(p.shape[1])[start]
    current = p @ weights
    for _ in range(max_iterations):
        choices = (1 - step) * current[:, None] + step * p
        losses = np.mean((choices - y[:, None]) ** 2, axis=0)
        winner = int(np.argmin(losses))
        if losses[winner] >= loss(current) - 1e-12:
            break
        weights *= 1 - step
        weights[winner] += step
        current = p @ weights
    return weights


def brier(y, p):
    return float(np.mean((p - y) ** 2))


def run(library_size):
    rows = []
    for rep in range(REPLICATES):
        rng = np.random.default_rng(SEED * 100 + rep)
        coef = np.zeros(FEATURES)
        coef[:6] = rng.normal(0, 0.8, 6)
        (X, y), (Xd, yd), (Xa, ya) = [generate(rng, n, coef) for n in (N_TRAIN, N_DEV, N_ASSESS)]
        columns = [rng.choice(FEATURES, 5, replace=False) for _ in range(FAMILIES)]   # each family reads five features
        Pd, Pa = np.zeros((N_DEV, library_size)), np.zeros((N_ASSESS, library_size))
        for j in range(library_size):
            cols = columns[j % FAMILIES]
            rows_used = np.arange(N_TRAIN) if j < FAMILIES else rng.integers(0, N_TRAIN, N_TRAIN)  # extras: bootstrap refits
            model = LogisticRegression(max_iter=300).fit(X[rows_used][:, cols], y[rows_used])
            Pd[:, j], Pa[:, j] = model.predict_proba(Xd[:, cols])[:, 1], model.predict_proba(Xa[:, cols])[:, 1]
        best = int(np.argmin([brier(yd, Pd[:, j]) for j in range(library_size)]))   # best single candidate, chosen on development rows
        w_all = greedy_brier_weights(yd, Pd)
        w_start = greedy_brier_weights(yd, Pd[:, :FAMILIES])                        # the same search on the four original models
        rows.append({
            "best": (brier(yd, Pd[:, best]), brier(ya, Pa[:, best])),
            "equal": (brier(yd, Pd.mean(axis=1)), brier(ya, Pa.mean(axis=1))),
            "greedy": (brier(yd, Pd @ w_all), brier(ya, Pa @ w_all)),
            "greedy_start": (brier(yd, Pd[:, :FAMILIES] @ w_start), brier(ya, Pa[:, :FAMILIES] @ w_start)),
            "kept": int((w_all > 0).sum()),
        })
    pick = lambda name, i: np.array([r[name][i] for r in rows])
    gain_dev = pick("greedy_start", 0) - pick("greedy", 0)        # what the extra candidates add on the rows the search used
    gain_assess = pick("greedy_start", 1) - pick("greedy", 1)     # what they add on assessment rows
    return clean({
        "library_size": library_size, "candidates_kept": float(np.mean([r["kept"] for r in rows])),
        "development": {k: float(pick(k, 0).mean()) for k in ("best", "equal", "greedy", "greedy_start")},
        "assessment": {k: float(pick(k, 1).mean()) for k in ("best", "equal", "greedy", "greedy_start")},
        "added_gain_development": float(gain_dev.mean()), "added_gain_assessment": float(gain_assess.mean()),
        "per_replicate": {"development": gain_dev, "assessment": gain_assess},
        "replicates": REPLICATES, "development_rows": N_DEV, "assessment_rows": N_ASSESS,
    })

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

VALUES = [4, 12, 40, 100]
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.ch17 import draw, NCOLS, HEIGHT

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