Twenty Questions: Decision Trees

# Chapter 12: a different kind of model. A decision tree learns by asking questions, not by turning knobs.
# Plain NumPy: no gradients anywhere, so no tensor.py this time. Save gini, best_split, grow_tree,
# predict_proba and count_leaves as tree.py, with 'import numpy as np' at the top: the capstone uses them.
import numpy as np

ACTIVITIES = ["picnic", "kite", "read"]
FEATURES = ["sun", "rain", "wind"]


def make_days(n, seed):
    """Chapter 9's process: read if rainy, else kite if windy, else picnic; 1 day in 10 at random."""
    r = np.random.default_rng(seed)
    X = np.round(r.random((n, 3)), 2)
    rule = np.where(X[:, 1] > 0.5, 2, np.where(X[:, 2] > 0.6, 1, 0))
    labels = rule.copy()
    noisy = r.random(n) < 0.1
    labels[noisy] = r.integers(0, 3, noisy.sum())
    return X, labels, rule

X_train, y_train, _ = make_days(150, 9)            # Chapter 9's 150 days
X_val, y_val, rule_val = make_days(1000, 99)       # Chapter 11's validation days
X_test, y_test, rule_test = make_days(1000, 123)   # and its test days, looked at once at the end

# === 1. How mixed is a group? ======================================================
def gini(counts):
    """Chance that two days drawn at random (with replacement) from the group did different things."""
    p = counts / counts.sum()
    return 1.0 - np.sum(p ** 2)


def entropy(counts):
    """Chapter 3's cross-entropy of the group against its own proportions, in nats."""
    p = counts[counts > 0] / counts.sum()
    return float(-np.sum(p * np.log(p))) + 0.0         # + 0.0 turns -0.0 into 0.0

root = np.bincount(y_train, minlength=3)
print("150 days:", dict(zip(ACTIVITIES, root.tolist())))
print(f"gini {gini(root):.4f}   entropy {entropy(root):.4f} nats   (an even 3-way mix: {gini(np.ones(3)):.4f}, {np.log(3):.4f})")
for name, c in [("all one activity", np.array([50, 0, 0])), ("45 / 5 / 0", np.array([45, 5, 0]))]:
    print(f"  {name:17} gini {gini(c):.4f}   entropy {entropy(c):.4f}")

# === 2. The tree: find the best question, split, repeat (save this part as tree.py) ===
def best_split(X, y, n_classes, min_leaf):
    """Try every feature and every threshold between two neighbouring values. Score = the children's gini,
    weighted by size. Counting is done once per distinct value, so a split costs O(n) per feature."""
    best = None
    n = len(y)
    for f in range(X.shape[1]):
        values, which = np.unique(X[:, f], return_inverse=True)
        if len(values) < 2:
            continue
        table = np.bincount(which * n_classes + y, minlength=len(values) * n_classes).reshape(len(values), n_classes)
        left = np.cumsum(table, axis=0)[:-1]          # class counts with x <= each value (all but the last)
        right = left[-1] + table[-1] - left
        nl, nr = left.sum(1), right.sum(1)
        ok = (nl >= min_leaf) & (nr >= min_leaf)
        if not ok.any():
            continue
        score = (nl - (left ** 2).sum(1) / np.maximum(nl, 1) + nr - (right ** 2).sum(1) / np.maximum(nr, 1)) / n
        score[~ok] = np.inf
        i = int(np.argmin(score))
        if best is None or score[i] < best[0]:
            best = (float(score[i]), f, float((values[i] + values[i + 1]) / 2))
    return best


def grow_tree(X, y, n_classes, max_depth=None, min_leaf=1, depth=0):
    """A node is a dict. A leaf holds the class counts of its training examples; a question node
    holds the feature, the threshold, and two subtrees: x[f] <= t goes left."""
    counts = np.bincount(y, minlength=n_classes)
    node = {"counts": counts}
    if (max_depth is not None and depth >= max_depth) or counts.max() == len(y):
        return node
    split = best_split(X, y, n_classes, min_leaf)
    if split is None or split[0] >= gini(counts) - 1e-12:   # no question makes the groups purer
        return node
    _, f, t = split
    go_left = X[:, f] <= t
    node.update(feature=f, threshold=t,
                left=grow_tree(X[go_left], y[go_left], n_classes, max_depth, min_leaf, depth + 1),
                right=grow_tree(X[~go_left], y[~go_left], n_classes, max_depth, min_leaf, depth + 1))
    return node


def predict_proba(node, X, alpha=0.0):
    """Send every row down the tree; each gets its leaf's class proportions. alpha adds that many
    pretend examples of every class, so no class ever gets probability exactly 0."""
    out = np.zeros((len(X), len(node["counts"])))
    rows = [(node, np.arange(len(X)))]
    while rows:
        nd, idx = rows.pop()
        if "feature" not in nd:
            c = nd["counts"] + alpha
            out[idx] = c / c.sum()
            continue
        go_left = X[idx, nd["feature"]] <= nd["threshold"]
        rows += [(nd["left"], idx[go_left]), (nd["right"], idx[~go_left])]
    return out


def count_leaves(node):
    return 1 if "feature" not in node else count_leaves(node["left"]) + count_leaves(node["right"])


def tree_depth(node):
    return 0 if "feature" not in node else 1 + max(tree_depth(node["left"]), tree_depth(node["right"]))

# === 3. The first question, and a depth-2 tree ======================================
def split_score(X, y, f, t, measure):
    """The two groups' impurity after asking 'is x[f] <= t?', each weighted by its share of the days."""
    left = X[:, f] <= t
    return sum(part.mean() * measure(np.bincount(y[part], minlength=3)) for part in (left, ~left) if part.any())

score, f, t = best_split(X_train, y_train, 3, 1)
print(f"\nbest first question: is {FEATURES[f]} <= {t:.3f}?   children's gini {score:.4f} (was {gini(root):.4f})")
print("best question about each feature, by gini and by entropy:")
for g in range(3):
    values = np.unique(X_train[:, g])
    thresholds = (values[:-1] + values[1:]) / 2
    for measure in (gini, entropy):
        s = [split_score(X_train, y_train, g, thr, measure) for thr in thresholds]
        i = int(np.argmin(s))
        print(f"  {FEATURES[g]:5} {measure.__name__:8} <= {thresholds[i]:.3f}   children {s[i]:.4f}"
              f"   (parent {measure(root):.4f})")


def show(node, depth=0, answer=""):
    pad = "    " * depth + answer
    if "feature" not in node:
        c = node["counts"]
        print(f"{pad}{ACTIVITIES[int(c.argmax())]}  {c.tolist()}")
        return
    print(f"{pad}{FEATURES[node['feature']]} <= {node['threshold']:.3f}?")
    show(node["left"], depth + 1, "yes: ")
    show(node["right"], depth + 1, "no:  ")

small = grow_tree(X_train, y_train, 3, max_depth=2)
print("\nthe depth-2 tree:")
show(small)


def accuracy(node, X, y):
    return float(np.mean(predict_proba(node, X).argmax(1) == y))


def log_loss(node, X, y, alpha=1.0):
    return float(-np.mean(np.log(predict_proba(node, X, alpha)[np.arange(len(y)), y])))

print(f"train {accuracy(small, X_train, y_train):.1%}   validation {accuracy(small, X_val, y_val):.1%}"
      f"   (the rule itself: {np.mean(rule_val == y_val):.1%})")

# === 4. Deeper trees memorize ==========================================================
print("\nmax depth   leaves   train acc   val acc   val loss (alpha = 1)")
for d in [1, 2, 3, 4, 5, None]:
    tr = grow_tree(X_train, y_train, 3, max_depth=d)
    print(f"{str(d) if d else 'no limit':>9}   {count_leaves(tr):6}   {accuracy(tr, X_train, y_train):9.1%}"
          f"   {accuracy(tr, X_val, y_val):7.1%}   {log_loss(tr, X_val, y_val):.3f}")
full = grow_tree(X_train, y_train, 3)
print(f"grown until every leaf is pure: depth {tree_depth(full)}, {count_leaves(full)} leaves for 150 days")

print("\nor stop splitting groups smaller than min_leaf:")
for m in [1, 5, 10, 20, 40]:
    tr = grow_tree(X_train, y_train, 3, min_leaf=m)
    print(f"  min_leaf {m:3}: {count_leaves(tr):3} leaves   train {accuracy(tr, X_train, y_train):.1%}"
          f"   validation {accuracy(tr, X_val, y_val):.1%}")

# === 5. The final exam, against Chapter 11's network ===================================
X_more, y_more, _ = make_days(1500, 10)            # Chapter 11's 1,500 days
print(f"\ntest set (Chapter 11: rule {np.mean(rule_test == y_test):.1%}, best network on 150 days 88.4%, "
      "unregularized network on 1,500 days 90.8%):")
for name, tr in [("depth 2, 150 days", small), ("no limit, 150 days", full),
                 ("depth 2, 1,500 days", grow_tree(X_more, y_more, 3, max_depth=2)),
                 ("no limit, 1,500 days", grow_tree(X_more, y_more, 3))]:
    print(f"  {name:22} validation {accuracy(tr, X_val, y_val):.1%}   test {accuracy(tr, X_test, y_test):.1%}"
          f"   leaves {count_leaves(tr)}")

# === 6. Why there is no gradient =======================================================
day = X_val[:1].copy()
p0 = predict_proba(small, day)[0]
day[0, 1] += 1e-4                                   # nudge the rain a little, as in Chapter 4
print(f"\nnudging one day's rain by 0.0001 changes its probabilities by {np.abs(predict_proba(small, day)[0] - p0).max():.1f}")