{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "tp7-000",
   "metadata": {},
   "source": [
    "# Practical 7 — Decision trees and random forests\n",
    "\n",
    "**Objectives:** program impurity, split search, recursive fitting and tree traversal with\n",
    "NumPy; connect your implementation to scikit-learn; implement bootstrap aggregation;\n",
    "select and evaluate a tree and a random forest.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-001",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "from sklearn.datasets import load_breast_cancer, make_moons\n",
    "from sklearn.dummy import DummyClassifier, DummyRegressor\n",
    "from sklearn.ensemble import RandomForestClassifier, RandomForestRegressor\n",
    "from sklearn.linear_model import LogisticRegression\n",
    "from sklearn.metrics import (accuracy_score, ConfusionMatrixDisplay, precision_score,\n",
    "                             recall_score, f1_score, roc_auc_score, root_mean_squared_error,\n",
    "                             mean_absolute_error, r2_score)\n",
    "from sklearn.model_selection import GridSearchCV, StratifiedKFold, KFold, train_test_split\n",
    "from sklearn.pipeline import make_pipeline\n",
    "from sklearn.preprocessing import StandardScaler\n",
    "from sklearn.svm import SVC\n",
    "from sklearn.tree import DecisionTreeClassifier, DecisionTreeRegressor, plot_tree\n",
    "\n",
    "SEED = 42\n",
    "BLUE, RED = \"#2474b5\", \"#e30613\"\n",
    "plt.rcParams.update({\"figure.figsize\": (7, 4), \"font.size\": 11})\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp7-002",
   "metadata": {},
   "source": [
    "## 1. Learn the best split yourself — 15 min\n",
    "\n",
    "Implement `gini(y)` and `entropy(y)` using class proportions $p_k$:\n",
    "\n",
    "$$G(y)=1-\\sum_k p_k^2,\\qquad H(y)=-\\sum_{k:p_k>0}p_k\\log_2 p_k.$$\n",
    "\n",
    "Use `np.unique(y, return_counts=True)` to count classes. Check a pure group and\n",
    "an equally mixed binary group before using the eight observations below.\n",
    "\n",
    "Complete the four holes in `split_candidates`: the two child sizes, their weighted\n",
    "impurity $Q$, and the gain $I(y)-Q$, where\n",
    "\n",
    "$$Q=\\frac{n_L}{n}I(y_L)+\\frac{n_R}{n}I(y_R).$$\n",
    "\n",
    "The feature loop, candidate thresholds and routing mask are supplied. The thresholds\n",
    "are midpoints between distinct values, so both children contain observations.\n",
    "Each candidate is stored as a dictionary; a smaller score is better.\n",
    "\n",
    "1. Print each threshold, child sizes, score and gain, first with Gini, then entropy.\n",
    "2. Find the best threshold using `min(candidates, key=...)`.\n",
    "3. Which split would you keep? Do the criteria agree? Are the children pure?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-003",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_small = np.array([.5, 1.2, 1.8, 2.6, 3.4, 4.1, 4.8, 5.5])[:, None]\n",
    "y_small = np.array([0, 0, 0, 0, 1, 1, 0, 1])\n",
    "print(\"x:\", X_small[:, 0])\n",
    "print(\"y:\", y_small)\n",
    "\n",
    "# Provided helper: dispatch to your Gini or entropy function.\n",
    "def impurity(y, criterion=\"gini\"):\n",
    "    if criterion == \"gini\":\n",
    "        return gini(y)\n",
    "    if criterion == \"entropy\":\n",
    "        return entropy(y)\n",
    "    raise ValueError(\"Use gini or entropy\")\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-004",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: implement gini(y) and entropy(y) using NumPy.\n",
    "# TODO: test [0, 0, 0, 0] and [0, 0, 1, 1].\n",
    "\n",
    "def split_candidates(X, y, criterion=\"gini\"):\n",
    "    candidates = []\n",
    "    for j in range(X.shape[1]):\n",
    "        values = np.unique(X[:, j])\n",
    "        thresholds = (values[:-1] + values[1:]) / 2\n",
    "        for threshold in thresholds:\n",
    "            left = X[:, j] <= threshold  # True for observations sent to the left.\n",
    "            n_left = ...   # TODO: number of observations in the left child.\n",
    "            n_right = ...  # TODO: number of observations in the right child.\n",
    "            score = ...    # TODO: weighted child impurity, using impurity(..., criterion).\n",
    "            gain = ...     # TODO: parent impurity minus score.\n",
    "            candidates.append({\"feature\": j, \"threshold\": threshold,\n",
    "                               \"n_left\": n_left, \"n_right\": n_right,\n",
    "                               \"score\": score, \"gain\": gain})\n",
    "    return candidates\n",
    "\n",
    "# TODO: call split_candidates for each criterion and print the candidate values.\n",
    "# TODO: identify the candidate with the smallest score.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp7-007",
   "metadata": {},
   "source": [
    "## 2. Fit and traverse a small tree — 20 min\n",
    "\n",
    "Generate 80 random points in [0,6]² with a fixed seed. Class 1 occupies the upper-right\n",
    "corner: x1>3 and x2>3. Plot the data before fitting.\n",
    "\n",
    "Complete the skeleton below to build a binary tree with NumPy:\n",
    "\n",
    "1. At a node, call `split_candidates`, keep the smallest score, and split the data.\n",
    "2. Call `fit_manual` on each child. Stop at `max_depth`, at a pure node, or if no\n",
    "   positive-gain split is possible. These stopping rules are supplied.\n",
    "3. For prediction, follow the tests until a leaf, then return its class-1 proportion.\n",
    "\n",
    "The tree is a nested dictionary: leaves store `proba`; internal nodes also store\n",
    "`feature`, `threshold`, `left` and `right`. Classify with `proba > .5` (ties give 0).\n",
    "Compare depth 1 and depth 2, then follow the points (4,1) and (4,5).\n",
    "Why does the second split help? Why are learned thresholds not necessarily exactly 3?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-008",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng = np.random.default_rng(SEED)\n",
    "X_toy = rng.uniform(0, 6, size=(80, 2))\n",
    "y_toy = ((X_toy[:, 0] > 3) & (X_toy[:, 1] > 3)).astype(int)\n",
    "for label, color in [(0, BLUE), (1, RED)]:\n",
    "    plt.scatter(*X_toy[y_toy == label].T, color=color, label=f\"Class {label}\")\n",
    "plt.xlabel(\"x1\")\n",
    "plt.ylabel(\"x2\")\n",
    "plt.legend()\n",
    "plt.show()\n",
    "\n",
    "def show_regions(ax, score_function, X, y, title):\n",
    "    low, high = X.min(axis=0) - .6, X.max(axis=0) + .6\n",
    "    xx, yy = np.meshgrid(np.linspace(low[0], high[0], 120),\n",
    "                         np.linspace(low[1], high[1], 120))\n",
    "    p = score_function(np.c_[xx.ravel(), yy.ravel()]).reshape(xx.shape)\n",
    "    ax.contourf(xx, yy, p > .5, levels=[-.5, .5, 1.5], colors=[BLUE, RED], alpha=.15)\n",
    "    for label, color in [(0, BLUE), (1, RED)]:\n",
    "        ax.scatter(*X[y == label].T, color=color, edgecolors=\"white\", label=f\"Class {label}\")\n",
    "    ax.set(xlabel=\"$x_1$\", ylabel=\"$x_2$\", title=title)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-009",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "def fit_manual(X, y, max_depth=2, criterion=\"gini\", depth=0):\n",
    "    node = {\"proba\": float(np.mean(y))}\n",
    "    if depth >= max_depth or impurity(y, criterion) < 1e-12:\n",
    "        return node\n",
    "    candidates = split_candidates(X, y, criterion)\n",
    "    if not candidates:\n",
    "        return node\n",
    "    # TODO: best = candidate with the smallest score.\n",
    "    # If best[\"gain\"] <= 1e-12, return node.\n",
    "    # TODO: obtain feature j and threshold, then create the left boolean mask.\n",
    "    # TODO: add feature, threshold, left and right to node.\n",
    "    # The children are calls to fit_manual on the corresponding subsets, depth+1.\n",
    "    return node\n",
    "\n",
    "def predict_proba_manual(tree, X):\n",
    "    probabilities = []\n",
    "    for x in X:\n",
    "        node = tree\n",
    "        # TODO: while \"feature\" in node, follow the left or right child.\n",
    "        probabilities.append(node[\"proba\"])\n",
    "    return np.array(probabilities)\n",
    "\n",
    "# TODO: fit depth 1 and depth 2; print training accuracy and draw regions with show_regions.\n",
    "# TODO: predict for np.array([[4., 1.], [4., 5.]]).\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp7-012",
   "metadata": {},
   "source": [
    "## 3. Connect the implementation to scikit-learn — 10 min\n",
    "\n",
    "1. Fit `DecisionTreeClassifier(max_depth=2, random_state=SEED)` on `X_toy, y_toy`.\n",
    "   Draw it with `plot_tree` and compare its predictions with your manual tree.\n",
    "2. Fit a stump (`max_depth=1`) on `X_small, y_small`. Does its threshold agree with Section 1?\n",
    "3. On the noisy moons split below, compare depths 1, 3 and None. Print training\n",
    "   accuracy, validation accuracy and leaf count; draw the training-data regions.\n",
    "\n",
    "Does a better training score imply a better validation score? Use no scaling for\n",
    "these trees. The moons experiment is independent of the later benchmark.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-013",
   "metadata": {},
   "outputs": [],
   "source": [
    "X_moons, y_moons = make_moons(n_samples=240, noise=.25, random_state=SEED)\n",
    "X_fit, X_val, y_fit, y_val = train_test_split(\n",
    "    X_moons, y_moons, test_size=.3, stratify=y_moons, random_state=SEED)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-014",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: fit the simple sklearn tree on X_toy/y_toy, then plot_tree.\n",
    "# TODO: compare its predictions with predict_proba_manual(manual_tree, X_toy) > .5.\n",
    "# TODO: inspect the stump's tree_.threshold[0] on X_small/y_small.\n",
    "# TODO: compare depths 1, 3, None on X_fit/y_fit; print scores and leaf counts.\n",
    "# Use show_regions with lambda Z: model.predict_proba(Z)[:, 1].\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp7-017",
   "metadata": {},
   "source": [
    "## 4. Program bootstrap aggregation, then use a random forest — 15 min\n",
    "\n",
    "Keep the moons training/validation split. Complete the three holes below:\n",
    "\n",
    "1. Draw n=len(y) indices **with replacement** using `rng.integers(0, n, size=n)`.\n",
    "   For example, `np.random.default_rng(0).integers(0, 5, size=5)` gives `[4, 3, 2, 1, 1]`:\n",
    "   observation 1 appears twice and observation 0 is omitted.\n",
    "2. Fit one `DecisionTreeClassifier` on the indexed observations. The surrounding\n",
    "   loop and storage of trees and sample indices are supplied.\n",
    "3. Average class-1 probabilities over trees with `np.mean(..., axis=0)`. The helper\n",
    "   below extracts the right probability column, including single-class trees.\n",
    "   For one observation, probabilities [.2, .7, .9] give .6, hence class 1.\n",
    "\n",
    "Run the supplied comparison cell after completing the functions. Compare the\n",
    "validation scores and boundaries of bagging and a random forest. What extra\n",
    "randomness does the forest add, and when is the feature subset sampled?\n",
    "\n",
    "The ensemble predicts class 1 when the mean probability exceeds .5. This can\n",
    "produce a different result from a majority vote over hard class predictions.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-bagging-probability-helper",
   "metadata": {},
   "outputs": [],
   "source": [
    "def class_one_probability(model, X):\n",
    "    # predict_proba columns follow model.classes_, not a fixed column number.\n",
    "    column = np.flatnonzero(model.classes_ == 1)\n",
    "    if not len(column):\n",
    "        # A tree trained without class 1 assigns it probability zero.\n",
    "        return np.zeros(len(X))\n",
    "    # If class 1 is the only class, its single column contains ones.\n",
    "    return model.predict_proba(X)[:, column[0]]\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-018",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "def fit_bagging(X, y, B=40, seed=SEED):\n",
    "    rng = np.random.default_rng(seed)  # Reproducible bootstrap draws.\n",
    "    trees, samples = [], []\n",
    "    n = len(y)\n",
    "    for b in range(B):\n",
    "        indices = ...  # TODO: draw n indices in [0, n), with replacement.\n",
    "        tree = DecisionTreeClassifier(random_state=seed+b)\n",
    "        # TODO: fit tree on X[indices] and y[indices].\n",
    "        trees.append(tree)\n",
    "        samples.append(indices)  # Retain duplicates: this is the bootstrap sample.\n",
    "    return trees, samples\n",
    "\n",
    "def bagging_probability(X):\n",
    "    # Shape: (number of trees, number of observations).\n",
    "    probabilities = np.array([class_one_probability(tree, X) for tree in bagged_trees])\n",
    "    return ...  # TODO: average over trees, returning one probability per observation.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-bagging-comparison",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Run after completing fit_bagging and bagging_probability above.\n",
    "bagged_trees, samples = fit_bagging(X_fit, y_fit)\n",
    "\n",
    "# sklearn's forest additionally randomizes candidate features at every node.\n",
    "forest_moons = RandomForestClassifier(n_estimators=100, max_features=\"sqrt\", random_state=SEED).fit(X_fit, y_fit)\n",
    "# Duplicates reduce the number of distinct observations in each sample.\n",
    "print(\"Mean unique bootstrap observations:\", np.mean([len(np.unique(i)) for i in samples]), \"of\", len(y_fit))\n",
    "# Average probabilities, then apply the binary classification threshold.\n",
    "print(\"Bagging validation accuracy:\", accuracy_score(y_val, bagging_probability(X_val) > .5))\n",
    "print(\"Forest validation accuracy:\", forest_moons.score(X_val, y_val))\n",
    "# Plot predictions in input space using the already-provided helper.\n",
    "fig, axes = plt.subplots(1, 2, figsize=(10, 4))\n",
    "show_regions(axes[0], bagging_probability, X_fit, y_fit, \"Manual bootstrap aggregation\")\n",
    "show_regions(axes[1], lambda Z: class_one_probability(forest_moons, Z), X_fit, y_fit, \"Random forest\")\n",
    "plt.tight_layout()\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp7-021",
   "metadata": {},
   "source": [
    "## 5. Select models on development data — 20 min\n",
    "\n",
    "Use the breast-cancer dataset from TP4/TP5: malignant +1, benign -1, same stratified\n",
    "20% test split and seed 42. Keep the test set unused until Section 6.\n",
    "\n",
    "Compare five procedures using the same three CV folds and mean **ROC-AUC**:\n",
    "\n",
    "| Procedure | Parameters to search |\n",
    "|---|---|\n",
    "| DummyClassifier(strategy=\"prior\") | no parameters |\n",
    "| StandardScaler + LogisticRegression(max_iter=3000) | C in [.1, 1, 10] |\n",
    "| StandardScaler + SVC(kernel=\"rbf\") | C in [.1, 1, 10], gamma in [\"scale\", .01, .1] |\n",
    "| DecisionTreeClassifier | criterion in [\"gini\", \"entropy\"], max_depth in [3, None], min_samples_leaf in [1, 5] |\n",
    "| RandomForestClassifier(n_estimators=100) | max_features in [\"sqrt\", 1.0], min_samples_leaf in [1, 5] |\n",
    "\n",
    "1. Run `GridSearchCV` with `return_train_score=True`, `n_jobs=1` and fixed tree seeds.\n",
    "2. Print mean train/CV AUC, CV standard deviation and selected parameters for each model.\n",
    "   Select `winner` by CV only. For SVM, AUC uses `decision_function`; `probability=True` is unnecessary.\n",
    "3. Plot the selected forest's ten largest `feature_importances_` values using NumPy and matplotlib.\n",
    "\n",
    "Why keep scaling inside the logistic and SVM pipelines? What changes with\n",
    "`max_features=1.0`? Can correlated variables share importance? Is importance causal?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-022",
   "metadata": {},
   "outputs": [],
   "source": [
    "data = load_breast_cancer()\n",
    "X = data.data\n",
    "y = np.where(data.target == 0, 1, -1)\n",
    "X_dev, X_test, y_dev, y_test = train_test_split(\n",
    "    X, y, test_size=.2, stratify=y, random_state=SEED)\n",
    "cv = StratifiedKFold(3, shuffle=True, random_state=SEED)\n",
    "print(\"Development:\", X_dev.shape, \"test:\", X_test.shape)\n",
    "print(\"Development class counts:\", dict(zip(*np.unique(y_dev, return_counts=True))))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-023",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: define specs = {name: (estimator, parameter_grid)} for the five procedures.\n",
    "# TODO: run GridSearchCV on X_dev/y_dev and store the fitted searches in searches.\n",
    "# TODO: store results as a list of dictionaries; print the comparison and select winner by CV.\n",
    "# TODO: use np.argsort and plt.barh for the forest's ten largest MDI importances.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp7-026",
   "metadata": {},
   "source": [
    "## 6. Final test-performance table and error analysis — 10 min\n",
    "\n",
    "Freeze all choices from Section 5. Evaluate each selected procedure once on the\n",
    "same test set. Print a table containing parameters, ROC-AUC, accuracy, malignant\n",
    "precision, malignant recall and malignant F1 (`pos_label=1`). Use `decision_function` for SVM AUC; for the other models, use the +1\n",
    "probability column located with `classes_`. Use `predict` for class metrics.\n",
    "\n",
    "Display the CV winner's confusion matrix with labels [-1, +1], and list the test-row\n",
    "indices of its false negatives. What is a false negative here? Do not change\n",
    "parameters or decision thresholds after inspecting test errors. Keep the CV winner\n",
    "even if another model has a larger test score.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-027",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: build test_table (a list of dictionaries) for the five frozen procedures.\n",
    "# Use decision_function for SVM; otherwise use class_one_probability.\n",
    "# TODO: print parameters, AUC, accuracy, malignant precision, recall and F1.\n",
    "# TODO: plot the CV winner's confusion matrix and print its false-negative indices.\n",
    "# TODO: interpret differences between AUC, recall and accuracy.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "tp7-034",
   "metadata": {
    "tags": [
     "optional"
    ]
   },
   "source": [
    "## Optional — Trees and forests for regression\n",
    "\n",
    "Use the generated data below to predict a noisy continuous target.\n",
    "\n",
    "1. Fit a dummy regressor, a regression tree and a random forest. Select parameters\n",
    "   on development data with three-fold CV and `scoring=\"neg_root_mean_squared_error\"`.\n",
    "   Try tree `max_depth` in [2, 5, None] and `min_samples_leaf` in [1, 5]; for a\n",
    "   100-tree forest, try `min_samples_leaf` in [1, 5]. Select by the highest CV score\n",
    "   (the lowest RMSE). Keep the test data unused during selection.\n",
    "2. Print the selected parameters and positive CV RMSE, then freeze the choices.\n",
    "   Evaluate each selected regressor once on the test set: RMSE, MAE and R².\n",
    "3. Plot the fitted tree and forest predictions over a dense x-grid, together with\n",
    "   development observations and the noise-free function.\n",
    "\n",
    "What does a squared-error leaf predict? Why is a tree prediction piecewise constant?\n",
    "Does averaging trees remove all noise? What happens near the edges of the observed range?\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-review-031",
   "metadata": {},
   "outputs": [],
   "source": [
    "rng_reg = np.random.default_rng(SEED)\n",
    "X_reg = rng_reg.uniform(-3, 3, size=(240, 1))\n",
    "def regression_signal(x):\n",
    "    return np.sin(2 * x) + .4 * x\n",
    "\n",
    "y_reg = regression_signal(X_reg[:, 0]) + rng_reg.normal(0, .3, len(X_reg))\n",
    "X_reg_dev, X_reg_test, y_reg_dev, y_reg_test = train_test_split(\n",
    "    X_reg, y_reg, test_size=.25, random_state=SEED)\n",
    "cv_reg = KFold(3, shuffle=True, random_state=SEED)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "tp7-review-032",
   "metadata": {
    "tags": [
     "exercise"
    ]
   },
   "outputs": [],
   "source": [
    "# TODO: GridSearchCV for DummyRegressor, DecisionTreeRegressor and RandomForestRegressor.\n",
    "# TODO: select by development CV; print selected parameters and positive CV RMSE.\n",
    "# TODO: print final test RMSE, MAE and R² for each selected regressor.\n",
    "# TODO: plot the two fitted prediction curves, development points and true signal.\n"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (IML)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.10.20"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}