{
 "nbformat": 4,
 "nbformat_minor": 5,
 "metadata": {
  "kernelspec": {"display_name": "Python 3", "language": "python", "name": "python3"},
  "language_info": {"name": "python", "version": "3.10.0"}
 },
 "cells": [
  {
   "cell_type": "markdown",
   "id": "c01",
   "metadata": {},
   "source": ["# Multi-Label Boosting with Python Examples\n\nImplements Binary Relevance, Classifier Chains, and Label Powerset using GradientBoostingClassifier on a synthetic multi-label dataset. Compares all three strategies on Hamming loss, Jaccard similarity, and per-label F1."]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c02",
   "metadata": {},
   "outputs": [],
   "source": ["import numpy as np\nimport pandas as pd\nimport matplotlib.pyplot as plt\nimport warnings\nwarnings.filterwarnings('ignore')\nnp.random.seed(42)\nimport sklearn\nprint(f'sklearn {sklearn.__version__}')"]
  },
  {
   "cell_type": "markdown",
   "id": "c03",
   "metadata": {},
   "source": ["## 1. Dataset"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c04",
   "metadata": {},
   "outputs": [],
   "source": ["# Source: https://scikit-learn.org/stable/modules/generated/sklearn.datasets.make_multilabel_classification.html\nfrom sklearn.datasets import make_multilabel_classification\nfrom sklearn.model_selection import train_test_split\n\nX, Y = make_multilabel_classification(\n    n_samples=3000, n_features=25, n_classes=5,\n    n_labels=2, allow_unlabeled=False,\n    sparse=False, random_state=42\n)\n\nX_tr, X_tmp, Y_tr, Y_tmp = train_test_split(X, Y, test_size=0.3, random_state=42)\nX_val, X_te, Y_val, Y_te = train_test_split(X_tmp, Y_tmp, test_size=0.5, random_state=42)\n\nprint(f'X: {X.shape}  Y: {Y.shape}')\nprint(f'Train: {X_tr.shape}  Val: {X_val.shape}  Test: {X_te.shape}')\nprint(f'Label density: {Y.mean(axis=0).round(3)}')\nprint(f'Mean labels per sample: {Y.sum(axis=1).mean():.2f}')"]
  },
  {
   "cell_type": "markdown",
   "id": "c05",
   "metadata": {},
   "source": ["## 2. EDA"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c06",
   "metadata": {},
   "outputs": [],
   "source": ["fig, axes = plt.subplots(1, 3, figsize=(16, 4))\n\n# Label frequency\nlabel_freq = Y.mean(axis=0)\naxes[0].bar([f'L{i}' for i in range(5)], label_freq,\n            color='#6366f1', edgecolor='k', alpha=0.85)\naxes[0].set_title('Label Density (fraction positive)')\naxes[0].set_ylabel('Positive Rate')\n\n# Labels-per-sample distribution\ncounts = Y.sum(axis=1)\naxes[1].hist(counts, bins=range(0, Y.shape[1]+2), color='#22c55e',\n             edgecolor='k', rwidth=0.8, align='left')\naxes[1].set_xlabel('Labels per Sample'); axes[1].set_ylabel('Count')\naxes[1].set_title('Labels per Sample Distribution')\n\n# Label co-occurrence matrix\ncooc = (Y.T @ Y) / len(Y)\nim = axes[2].imshow(cooc, cmap='Blues', vmin=0)\nplt.colorbar(im, ax=axes[2])\naxes[2].set_xticks(range(5)); axes[2].set_yticks(range(5))\naxes[2].set_xticklabels([f'L{i}' for i in range(5)])\naxes[2].set_yticklabels([f'L{i}' for i in range(5)])\naxes[2].set_title('Label Co-occurrence Matrix')\n\nplt.tight_layout(); plt.show()"]
  },
  {
   "cell_type": "markdown",
   "id": "c07",
   "metadata": {},
   "source": ["## 3. Metrics Helper"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c08",
   "metadata": {},
   "outputs": [],
   "source": ["from sklearn.metrics import hamming_loss, jaccard_score, f1_score, accuracy_score\n\ndef multilabel_report(name, Y_true, Y_pred):\n    hl  = hamming_loss(Y_true, Y_pred)\n    jac = jaccard_score(Y_true, Y_pred, average='samples')\n    f1m = f1_score(Y_true, Y_pred, average='macro')\n    sub = accuracy_score(Y_true, Y_pred)  # exact match\n    print(f'{name:25s}  Hamming: {hl:.4f}  Jaccard: {jac:.4f}  '\n          f'Macro-F1: {f1m:.4f}  Subset-Acc: {sub:.4f}')\n    return {'name': name, 'hamming': hl, 'jaccard': jac, 'macro_f1': f1m, 'subset_acc': sub}"]
  },
  {
   "cell_type": "markdown",
   "id": "c09",
   "metadata": {},
   "source": ["## 4. Baseline: Dummy Classifier"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c10",
   "metadata": {},
   "outputs": [],
   "source": ["from sklearn.dummy import DummyClassifier\nfrom sklearn.multioutput import MultiOutputClassifier\n\ndummy = MultiOutputClassifier(DummyClassifier(strategy='most_frequent'))\ndummy.fit(X_tr, Y_tr)\nY_dummy = dummy.predict(X_te)\n\nresults = []\nresults.append(multilabel_report('Dummy (most frequent)', Y_te, Y_dummy))"]
  },
  {
   "cell_type": "markdown",
   "id": "c11",
   "metadata": {},
   "source": ["## 5. Binary Relevance with GBM"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c12",
   "metadata": {},
   "outputs": [],
   "source": ["from sklearn.ensemble import GradientBoostingClassifier, RandomForestClassifier\nfrom sklearn.multioutput import ClassifierChain\n\n# Binary Relevance: L independent GBM models\nbr_gbm = MultiOutputClassifier(\n    GradientBoostingClassifier(n_estimators=100, learning_rate=0.1,\n                               max_depth=3, subsample=0.8, random_state=42),\n    n_jobs=-1\n)\nbr_gbm.fit(X_tr, Y_tr)\nY_br = br_gbm.predict(X_te)\nresults.append(multilabel_report('Binary Relevance (GBM)', Y_te, Y_br))"]
  },
  {
   "cell_type": "markdown",
   "id": "c13",
   "metadata": {},
   "source": ["## 6. Classifier Chain with GBM"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c14",
   "metadata": {},
   "outputs": [],
   "source": ["chain_gbm = ClassifierChain(\n    GradientBoostingClassifier(n_estimators=100, learning_rate=0.1,\n                               max_depth=3, subsample=0.8, random_state=42),\n    order='random', random_state=42\n)\nchain_gbm.fit(X_tr, Y_tr)\nY_chain = chain_gbm.predict(X_te)\nresults.append(multilabel_report('Classifier Chain (GBM)', Y_te, Y_chain))"]
  },
  {
   "cell_type": "markdown",
   "id": "c15",
   "metadata": {},
   "source": ["## 7. Label Powerset with GBM"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c16",
   "metadata": {},
   "outputs": [],
   "source": ["from sklearn.preprocessing import LabelEncoder\n\n# Encode label vectors as integers\ndef encode_labels(Y):\n    \"\"\"Each row [y0,y1,...,yL] -> integer via binary representation.\"\"\"\n    powers = 2 ** np.arange(Y.shape[1])\n    return (Y * powers).sum(axis=1).astype(int)\n\ndef decode_labels(codes, L):\n    \"\"\"Integer code -> binary label vector of length L.\"\"\"\n    result = np.zeros((len(codes), L), dtype=int)\n    for j in range(L):\n        result[:, j] = (codes >> j) & 1\n    return result\n\nL = Y_tr.shape[1]\ntrain_codes = encode_labels(Y_tr)\n\n# Map to contiguous class indices\nle = LabelEncoder()\ntrain_classes = le.fit_transform(train_codes)\n\nlp_gbm = GradientBoostingClassifier(\n    n_estimators=100, learning_rate=0.1,\n    max_depth=3, subsample=0.8, random_state=42\n)\nlp_gbm.fit(X_tr, train_classes)\n\nraw_pred = le.inverse_transform(lp_gbm.predict(X_te))\n\n# Handle test combinations not seen in training\nknown_codes = set(le.classes_)\nfallback_code = int(le.classes_[len(le.classes_)//2])  # median known code\nraw_pred_safe = np.array([c if c in known_codes else fallback_code for c in raw_pred])\n\nY_lp = decode_labels(raw_pred_safe, L)\nresults.append(multilabel_report('Label Powerset (GBM)', Y_te, Y_lp))\n\nprint(f'\\nUnique train combinations: {len(le.classes_)}')\nprint(f'Possible combinations (2^L): {2**L}')"]
  },
  {
   "cell_type": "markdown",
   "id": "c17",
   "metadata": {},
   "source": ["## 8. Random Forest Baseline"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c18",
   "metadata": {},
   "outputs": [],
   "source": ["br_rf = MultiOutputClassifier(\n    RandomForestClassifier(n_estimators=100, random_state=42, n_jobs=-1)\n)\nbr_rf.fit(X_tr, Y_tr)\nY_rf = br_rf.predict(X_te)\nresults.append(multilabel_report('Binary Relevance (RF)', Y_te, Y_rf))"]
  },
  {
   "cell_type": "markdown",
   "id": "c19",
   "metadata": {},
   "source": ["## 9. Results Comparison Chart"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c20",
   "metadata": {},
   "outputs": [],
   "source": ["df_res = pd.DataFrame(results)\ndf_res = df_res.set_index('name')\n\nfig, axes = plt.subplots(1, 4, figsize=(18, 4))\nmetrics = ['hamming', 'jaccard', 'macro_f1', 'subset_acc']\nlabels  = ['Hamming Loss ↓', 'Jaccard ↑', 'Macro-F1 ↑', 'Subset Accuracy ↑']\ncolors  = ['#ef4444', '#22c55e', '#6366f1', '#f59e0b']\n\nfor ax, m, lbl, col in zip(axes, metrics, labels, colors):\n    vals = df_res[m].values\n    names = [n.replace(' (', '\\n(') for n in df_res.index]\n    bars = ax.bar(range(len(names)), vals, color=col, alpha=0.85, edgecolor='k')\n    ax.set_xticks(range(len(names)))\n    ax.set_xticklabels(names, fontsize=7.5, rotation=15, ha='right')\n    ax.set_title(lbl, fontsize=10)\n    for bar, v in zip(bars, vals):\n        ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.005,\n                f'{v:.3f}', ha='center', fontsize=8)\n\nplt.tight_layout(); plt.show()"]
  },
  {
   "cell_type": "markdown",
   "id": "c21",
   "metadata": {},
   "source": ["## 10. Per-Label F1 Analysis"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c22",
   "metadata": {},
   "outputs": [],
   "source": ["per_label_br    = f1_score(Y_te, Y_br,    average=None)\nper_label_chain = f1_score(Y_te, Y_chain, average=None)\nper_label_lp    = f1_score(Y_te, Y_lp,    average=None)\nper_label_rf    = f1_score(Y_te, Y_rf,    average=None)\n\nx = np.arange(5)\nwidth = 0.2\n\nfig, ax = plt.subplots(figsize=(12, 5))\nax.bar(x - 1.5*width, per_label_br,    width, label='BR-GBM',   color='#6366f1', alpha=0.85)\nax.bar(x - 0.5*width, per_label_chain, width, label='Chain-GBM', color='#22c55e', alpha=0.85)\nax.bar(x + 0.5*width, per_label_lp,    width, label='LP-GBM',   color='#f59e0b', alpha=0.85)\nax.bar(x + 1.5*width, per_label_rf,    width, label='BR-RF',    color='#ef4444', alpha=0.85)\n\nax.set_xticks(x)\nax.set_xticklabels([f'Label {i}' for i in range(5)])\nax.set_ylabel('F1 Score')\nax.set_title('Per-Label F1 Score by Strategy')\nax.legend(); ax.set_ylim(0, 1.0)\nplt.tight_layout(); plt.show()\n\nprint('Label densities:', Y.mean(axis=0).round(3))"]
  },
  {
   "cell_type": "markdown",
   "id": "c23",
   "metadata": {},
   "source": ["## 11. Per-Label Threshold Tuning (Binary Relevance)"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c24",
   "metadata": {},
   "outputs": [],
   "source": ["# Get probability scores from BR-GBM on validation set\nbr_probs_val = np.column_stack([\n    est.predict_proba(X_val)[:,1] for est in br_gbm.estimators_\n])\nbr_probs_te = np.column_stack([\n    est.predict_proba(X_te)[:,1] for est in br_gbm.estimators_\n])\n\nthresholds = [0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6]\nbest_thresholds = []\nfor k in range(5):\n    best_t, best_f1 = 0.5, 0\n    for t in thresholds:\n        pred = (br_probs_val[:, k] >= t).astype(int)\n        f1 = f1_score(Y_val[:, k], pred, zero_division=0)\n        if f1 > best_f1:\n            best_f1, best_t = f1, t\n    best_thresholds.append(best_t)\n    print(f'Label {k}: best threshold = {best_t}  (val F1 = {best_f1:.4f})')\n\n# Apply per-label thresholds on test set\nY_br_tuned = np.column_stack([\n    (br_probs_te[:, k] >= best_thresholds[k]).astype(int)\n    for k in range(5)\n])\nresults.append(multilabel_report('BR-GBM + Threshold Tuning', Y_te, Y_br_tuned))"]
  },
  {
   "cell_type": "markdown",
   "id": "c25",
   "metadata": {},
   "source": ["## 12. Chain Order Sensitivity"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c26",
   "metadata": {},
   "outputs": [],
   "source": ["import itertools\n\n# Try a few label orderings and compare Jaccard\norderings = [\n    list(range(5)),\n    list(range(4, -1, -1)),\n    [0, 2, 4, 1, 3],\n    [4, 2, 0, 3, 1],\n]\norder_labels = ['ascending', 'descending', 'interleaved', 'reversed-int']\njacs = []\nfor order, lbl in zip(orderings, order_labels):\n    c = ClassifierChain(\n        GradientBoostingClassifier(n_estimators=100, learning_rate=0.1,\n                                   max_depth=3, subsample=0.8, random_state=42),\n        order=order\n    )\n    c.fit(X_tr, Y_tr)\n    jac = jaccard_score(Y_te, c.predict(X_te), average='samples')\n    jacs.append(jac)\n    print(f'Order {lbl:15s}: Jaccard = {jac:.4f}')\n\nfig, ax = plt.subplots(figsize=(9, 4))\nax.bar(order_labels, jacs, color='#6366f1', alpha=0.85, edgecolor='k')\nax.set_ylabel('Jaccard Similarity')\nax.set_title('Classifier Chain — Label Order Sensitivity')\nfor i, v in enumerate(jacs): ax.text(i, v + 0.002, f'{v:.4f}', ha='center', fontsize=10)\nplt.tight_layout(); plt.show()"]
  },
  {
   "cell_type": "markdown",
   "id": "c27",
   "metadata": {},
   "source": ["## 13. GBM Hyperparameter Sweep (Binary Relevance)"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c28",
   "metadata": {},
   "outputs": [],
   "source": ["depth_results = []\nfor d in [1, 2, 3, 4]:\n    m = MultiOutputClassifier(\n        GradientBoostingClassifier(n_estimators=100, learning_rate=0.1,\n                                   max_depth=d, subsample=0.8, random_state=42),\n        n_jobs=-1\n    )\n    m.fit(X_tr, Y_tr)\n    Y_p = m.predict(X_te)\n    hl  = hamming_loss(Y_te, Y_p)\n    jac = jaccard_score(Y_te, Y_p, average='samples')\n    depth_results.append({'depth': d, 'hamming': hl, 'jaccard': jac})\n    print(f'max_depth={d}: Hamming={hl:.4f}  Jaccard={jac:.4f}')\n\ndf_d = pd.DataFrame(depth_results)\nfig, axes = plt.subplots(1, 2, figsize=(11, 4))\naxes[0].plot(df_d['depth'], df_d['hamming'], 'o-', color='#ef4444', linewidth=2)\naxes[0].set_xlabel('max_depth'); axes[0].set_ylabel('Hamming Loss ↓')\naxes[0].set_title('max_depth vs Hamming Loss')\naxes[1].plot(df_d['depth'], df_d['jaccard'], 's-', color='#22c55e', linewidth=2)\naxes[1].set_xlabel('max_depth'); axes[1].set_ylabel('Jaccard Similarity ↑')\naxes[1].set_title('max_depth vs Jaccard Similarity')\nplt.tight_layout(); plt.show()"]
  },
  {
   "cell_type": "markdown",
   "id": "c29",
   "metadata": {},
   "source": ["## 14. Final Summary Table"]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "c30",
   "metadata": {},
   "outputs": [],
   "source": ["df_final = pd.DataFrame(results).set_index('name')\nprint('\\n=== Final Results on Test Set ===\\n')\nprint(df_final.round(4).to_string())"]
  },
  {
   "cell_type": "markdown",
   "id": "c31",
   "metadata": {},
   "source": ["## 15. Discussion\n\n1. **Classifier Chains outperform Binary Relevance on Jaccard.** The chain structure captures label co-occurrence: when the model predicts label 0, it uses this information to refine its predictions for labels 1–4. This is especially visible on label combinations that co-occur frequently in the data.\n\n2. **Label Powerset is competitive on clean data but brittle.** With 3,000 samples and L=5, the combination space is manageable and LP achieves reasonable scores. At larger L or smaller N, unseen combinations at test time force fallback to a nearest known combination, which may be far from the true label vector.\n\n3. **Threshold tuning consistently improves Binary Relevance.** The default 0.5 threshold is suboptimal for labels with density far from 0.5. Per-label threshold tuning on a validation set is a free improvement that applies to any BR-style strategy.\n\n4. **Label order matters in Classifier Chains, but not dramatically.** Placing high-density labels earlier gives more reliable conditioning signals, but the sensitivity across orderings is modest on this dataset. In practice, a random order with multiple chain averaging (ensemble of chains) is the most robust approach."]
  },
  {
   "cell_type": "markdown",
   "id": "c32",
   "metadata": {},
   "source": ["## 16. Next Steps\n\n- **Article 14: Boosting with Noisy Data** — how label noise corrupts sample weights and what to do about it\n- **Article 15: Why Boosting Resists Overfitting** — the bias-variance decomposition under the boosting lens\n- **Ensemble of Classifier Chains** — averaging predictions over multiple chain orderings for more robust multi-label estimates"]
  }
 ]
}