Machine Learning Aide-mémoire

# Machine Learning: systems that learn patterns from data
# instead of being explicitly programmed.

# Three main paradigms:
# - Supervised:    learn f: X -> y from labeled (X, y) pairs
# - Unsupervised:  find structure in unlabeled X
# - Reinforcement: learn actions that maximize reward via trial-and-error

# Typical workflow:
# 1. Define problem & collect data
# 2. Preprocess & feature engineering
# 3. Split train / validation / test
# 4. Train model, tune hyperparameters
# 5. Evaluate on test set
# 6. Deploy & monitor

# The standard ML stack in Python
import numpy as np          # numerical arrays
import pandas as pd         # tabular data
import matplotlib.pyplot as plt  # plotting
import seaborn as sns       # statistical plots

# scikit-learn: classical ML (trees, SVM, linear models, preprocessing)
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score

# Deep learning (pick one):
import torch                # PyTorch
# import tensorflow as tf   # TensorFlow/Keras

# Install:
# pip install numpy pandas matplotlib seaborn scikit-learn torch

from sklearn.model_selection import train_test_split

X_train, X_test, y_train, y_test = train_test_split(
    X, y,
    test_size=0.2,      # 20% held out for test
    random_state=42,    # reproducibility
    stratify=y,         # keep class proportions (classification only)
)

# For a third validation set:
X_train, X_temp, y_train, y_temp = train_test_split(
    X, y, test_size=0.4, stratify=y, random_state=42)
X_val, X_test, y_val, y_test = train_test_split(
    X_temp, y_temp, test_size=0.5, stratify=y_temp, random_state=42)
# Now: 60% train, 20% val, 20% test

from sklearn.datasets import load_iris
from sklearn.model_selection import train_test_split
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import accuracy_score, classification_report

iris = load_iris()
X, y = iris.data, iris.target

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, stratify=y, random_state=42)

model = RandomForestClassifier(n_estimators=100, random_state=42)
model.fit(X_train, y_train)
y_pred = model.predict(X_test)

print("Accuracy:", accuracy_score(y_test, y_pred))
print(classification_report(y_test, y_pred, target_names=iris.target_names))

import numpy as np
import random
import torch

SEED = 42
np.random.seed(SEED)
random.seed(SEED)
torch.manual_seed(SEED)
if torch.cuda.is_available():
    torch.cuda.manual_seed_all(SEED)

# For sklearn, pass random_state to each estimator:
# RandomForestClassifier(random_state=SEED)

# Force deterministic operations (slower but fully reproducible):
torch.use_deterministic_algorithms(True)

# Optional: pin GPU non-determinism
import os
os.environ["CUBLAS_WORKSPACE_CONFIG"] = ":4096:8"

import pandas as pd
from sklearn.impute import SimpleImputer, KNNImputer

df = pd.DataFrame({"age": [25, np.nan, 35, 45],
                   "income": [50000, 60000, np.nan, 80000]})

# Drop rows with any missing values (use sparingly)
df_drop = df.dropna()

# Mean/median/mode imputation
imp = SimpleImputer(strategy="median")  # 'mean','most_frequent','constant'
df_imputed = pd.DataFrame(imp.fit_transform(df), columns=df.columns)

# KNN imputer (uses similar samples)
knn_imp = KNNImputer(n_neighbors=3)
df_knn = pd.DataFrame(knn_imp.fit_transform(df), columns=df.columns)

# Always fit imputers on TRAINING data only, then transform test
imp.fit(X_train)
X_train_imp = imp.transform(X_train)
X_test_imp  = imp.transform(X_test)   # uses train medians

from sklearn.preprocessing import StandardScaler, MinMaxScaler, RobustScaler

# Standardization: mean 0, std 1 (most common)
scaler = StandardScaler()
X_train_scaled = scaler.fit_transform(X_train)
X_test_scaled  = scaler.transform(X_test)   # use train mean/std!

# Min-max: scale to [0, 1]
minmax = MinMaxScaler()
X_mm = minmax.fit_transform(X_train)

# Robust to outliers: uses median & IQR
robust = RobustScaler()
X_robust = robust.fit_transform(X_train)

from sklearn.preprocessing import LabelEncoder, OneHotEncoder, OrdinalEncoder

# Label encode TARGET (classification): strings -> integers
le = LabelEncoder()
y = le.fit_transform(["cat", "dog", "cat", "bird"])  # [0, 1, 0, 2]
le.inverse_transform([0, 1, 2])                       # back to strings

# One-hot encode FEATURES (nominal, no order)
ohe = OneHotEncoder(sparse_output=False, handle_unknown="ignore")
X_ohe = ohe.fit_transform(df[["color"]])  # red/green/blue -> 3 cols

# Ordinal encode (ordered categories)
oe = OrdinalEncoder(categories=[["low", "medium", "high"]])
X_ord = oe.fit_transform(df[["size"]])

# For high-cardinality features, consider target encoding
# or just drop them if they're IDs.

from sklearn.utils import resample
from imblearn.over_sampling import SMOTE
from sklearn.utils.class_weight import compute_class_weight

# 1. Class weights (no resampling — preferred for trees/linear)
classes = np.unique(y_train)
weights = compute_class_weight("balanced", classes=classes, y=y_train)
class_weight = dict(zip(classes, weights))
model = RandomForestClassifier(class_weight=class_weight)

# 2. Oversample minority
df_majority = df[df.y == 0]
df_minority = df[df.y == 1]
df_minority_up = resample(df_minority, replace=True,
                          n_samples=len(df_majority), random_state=42)

# 3. SMOTE: synthesize new minority samples via k-nearest neighbors
smote = SMOTE(random_state=42)
X_res, y_res = smote.fit_resample(X_train, y_train)

# Always resample on TRAINING data only — never the test set!

import numpy as np
from sklearn.ensemble import IsolationForest

# Z-score: |x - mean| / std > 3 is an outlier
z_scores = np.abs((X - X.mean()) / X.std())
outliers = (z_scores > 3).any(axis=1)

# IQR method (more robust to outliers than z-score)
Q1, Q3 = X.quantile(0.25), X.quantile(0.75)
IQR = Q3 - Q1
outliers_iqr = ((X < (Q1 - 1.5 * IQR)) | (X > (Q3 + 1.5 * IQR))).any(axis=1)

# Isolation Forest (good for high-dimensional data)
iso = IsolationForest(contamination=0.05, random_state=42)
outliers_iso = iso.fit_predict(X)  # -1 = outlier, 1 = inlier

# Remove or cap (winsorize) outliers before training
X_clean = X[~outliers]

from sklearn.preprocessing import PolynomialFeatures

# Generate degree-2 features including interactions
poly = PolynomialFeatures(degree=2, interaction_only=False,
                          include_bias=False)
X_poly = poly.fit_transform(X)
# For X = [a, b], produces: [a, b, a^2, a*b, b^2]

# interaction_only=True keeps just [a, b, a*b] — no squares
poly_int = PolynomialFeatures(degree=2, interaction_only=True,
                              include_bias=False)
X_int = poly_int.fit_transform(X)

# Useful for linear models to capture non-linearity
# WARNING: feature count grows quickly — use with regularization

from sklearn.preprocessing import KBinsDiscretizer
import pandas as pd

# Equal-width bins (uniform spacing)
df["age_bin"] = pd.cut(df["age"], bins=5, labels=False)

# Equal-frequency bins (quantiles — same count per bin)
df["age_quantile"] = pd.qcut(df["age"], q=5, labels=False,
                              duplicates="drop")

# KBinsDiscretizer with strategy
kbd = KBinsDiscretizer(n_bins=5, encode="onehot-dense",
                        strategy="quantile")  # 'uniform','kmeans'
X_binned = kbd.fit_transform(X[["age"]])

# Custom bins with domain knowledge
df["age_group"] = pd.cut(df["age"],
                          bins=[0, 18, 35, 60, 100],
                          labels=["child","young","adult","senior"])

from sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer

corpus = [
    "machine learning is fun",
    "deep learning is a subset of machine learning",
    "python is great for data science",
]

# Bag of words (counts)
cv = CountVectorizer(max_features=1000, stop_words="english",
                     ngram_range=(1, 2))
X_counts = cv.fit_transform(corpus)

# TF-IDF: downweights common words
tfidf = TfidfVectorizer(max_features=1000, stop_words="english",
                        ngram_range=(1, 2), min_df=2)
X_tfidf = tfidf.fit_transform(corpus)

# For production, consider sentence-transformers for embeddings:
# from sentence_transformers import SentenceTransformer
# model = SentenceTransformer("all-MiniLM-L6-v2")
# embeddings = model.encode(corpus)

import pandas as pd

df["date"] = pd.to_datetime(df["date"])

# Extract calendar parts
df["year"]      = df["date"].dt.year
df["month"]     = df["date"].dt.month
df["day"]       = df["date"].dt.day
df["dayofweek"] = df["date"].dt.dayofweek    # 0=Mon
df["dayofyear"] = df["date"].dt.dayofyear
df["hour"]      = df["date"].dt.hour
df["is_weekend"] = df["date"].dt.dayofweek >= 5

# Cyclic encoding (preserves continuity: hour 23 -> 0)
df["hour_sin"] = np.sin(2 * np.pi * df["hour"] / 24)
df["hour_cos"] = np.cos(2 * np.pi * df["hour"] / 24)
df["month_sin"] = np.sin(2 * np.pi * df["month"] / 12)
df["month_cos"] = np.cos(2 * np.pi * df["month"] / 12)

# Time elapsed (useful as a trend feature)
df["days_since"] = (df["date"] - df["date"].min()).dt.days

from sklearn.feature_selection import (SelectKBest, f_classif,
    RFE, SelectFromModel)
from sklearn.ensemble import RandomForestClassifier

# 1. Univariate: keep top-k by statistical test
selector = SelectKBest(f_classif, k=10)
X_new = selector.fit_transform(X, y)
selected = selector.get_support(indices=True)

# 2. Recursive Feature Elimination (model-based, iterative)
rfe = RFE(RandomForestClassifier(random_state=42),
          n_features_to_select=10)
X_rfe = rfe.fit_transform(X, y)

# 3. Select from model importance
sfm = SelectFromModel(RandomForestClassifier(random_state=42),
                      threshold="median")
X_sfm = sfm.fit_transform(X, y)

# 4. Correlation filter — drop one of highly-correlated pairs
corr = pd.DataFrame(X).corr().abs()
to_drop = [c for c in corr.columns
           if any(corr[c] > 0.95)]
X_filtered = X.drop(columns=to_drop)

from sklearn.linear_model import LinearRegression, LogisticRegression

# Regression: predict a continuous value
lr = LinearRegression()
lr.fit(X_train, y_train)
y_pred = lr.predict(X_test)

# Inspect coefficients
print("Coef:", lr.coef_, "Intercept:", lr.intercept_)

# Classification: predict a class (uses softmax for multiclass)
clf = LogisticRegression(
    penalty="l2",            # 'l1' (Lasso), 'elasticnet', None
    C=1.0,                   # inverse of regularization strength
    max_iter=1000,
    multi_class="multinomial",
)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
y_proba = clf.predict_proba(X_test)  # class probabilities

from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier

# Single tree (high variance, easy to overfit)
tree = DecisionTreeClassifier(
    max_depth=5,              # limit depth to prevent overfitting
    min_samples_leaf=10,      # min samples per leaf
    random_state=42,
)
tree.fit(X_train, y_train)

# Random forest: ensemble of trees on bootstrapped samples
rf = RandomForestClassifier(
    n_estimators=200,         # number of trees
    max_depth=None,
    max_features="sqrt",      # features per split
    n_jobs=-1,                # use all CPUs
    random_state=42,
)
rf.fit(X_train, y_train)

# Feature importance (averaged across trees)
importances = pd.Series(rf.feature_importances_, index=X.columns)

from sklearn.ensemble import GradientBoostingClassifier, HistGradientBoostingClassifier
import lightgbm as lgb
import xgboost as xgb

# sklearn (slow on large data)
gbm = GradientBoostingClassifier(
    n_estimators=200, learning_rate=0.1, max_depth=3, random_state=42)
gbm.fit(X_train, y_train)

# sklearn's faster histogram-based version (like LightGBM)
hgbm = HistGradientBoostingClassifier(
    max_iter=200, learning_rate=0.1, max_leaf_nodes=31, random_state=42)
hgbm.fit(X_train, y_train)

# XGBoost (popular in competitions)
xgb_clf = xgb.XGBClassifier(
    n_estimators=200, learning_rate=0.1, max_depth=6,
    subsample=0.8, colsample_bytree=0.8, n_jobs=-1,
    eval_metric="logloss", random_state=42)
xgb_clf.fit(X_train, y_train, eval_set=[(X_val, y_val)], early_stopping_rounds=10)

# LightGBM (fastest on large data)
lgb_clf = lgb.LGBMClassifier(n_estimators=200, learning_rate=0.1, num_leaves=31)
lgb_clf.fit(X_train, y_train, eval_set=(X_val, y_val))

from sklearn.svm import SVC, SVR
from sklearn.preprocessing import StandardScaler

# SVC requires scaling!
scaler = StandardScaler()
X_train_s = scaler.fit_transform(X_train)
X_test_s  = scaler.transform(X_test)

# Kernels: 'linear', 'rbf' (Gaussian), 'poly', 'sigmoid'
svm = SVC(
    kernel="rbf",
    C=1.0,              # regularization (higher = fit harder)
    gamma="scale",      # kernel width
    probability=True,   # enable predict_proba (slower)
    random_state=42,
)
svm.fit(X_train_s, y_train)
y_pred = svm.predict(X_test_s)

# SVR for regression
svr = SVR(kernel="rbf", C=1.0, epsilon=0.1)
svr.fit(X_train_s, y_train)

from sklearn.neighbors import KNeighborsClassifier, KNeighborsRegressor

# KNN: classify by majority vote of k nearest training samples
knn = KNeighborsClassifier(
    n_neighbors=5,        # k — try odd values to break ties
    weights="uniform",    # or 'distance' (closer = more weight)
    metric="minkowski",   # 'euclidean','manhattan', etc.
    p=2,                  # 2=euclidean, 1=manhattan
    n_jobs=-1,
)
knn.fit(X_train_scaled, y_train)
y_pred = knn.predict(X_test_scaled)

# Choose k via cross-validation
from sklearn.model_selection import cross_val_score
for k in [1, 3, 5, 7, 9, 11]:
    score = cross_val_score(KNeighborsClassifier(n_neighbors=k),
                            X_train_scaled, y_train, cv=5).mean()
    print(f"k={k}: {score:.3f}")

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score

# Choose k via the elbow method or silhouette score
inertias, silhouettes = [], []
for k in range(2, 11):
    km = KMeans(n_clusters=k, n_init=10, random_state=42)
    labels = km.fit_predict(X)
    inertias.append(km.inertia_)
    silhouettes.append(silhouette_score(X, labels))

# Final model with chosen k
kmeans = KMeans(n_clusters=4, n_init=10, random_state=42)
clusters = kmeans.fit_predict(X)
centers = kmeans.cluster_centers_

# Predict cluster for new points
new_labels = kmeans.predict(X_new)

from sklearn.cluster import DBSCAN

# DBSCAN groups densely-connected points; marks sparse ones as noise (-1)
db = DBSCAN(
    eps=0.5,            # max distance between points in a cluster
    min_samples=5,      # min points to form a dense region
    metric="euclidean",
)
labels = db.fit_predict(X)
n_clusters = len(set(labels)) - (1 if -1 in labels else 0)
n_noise = list(labels).count(-1)

# Tune eps via the k-distance plot (k = min_samples)
from sklearn.neighbors import NearestNeighbors
nn = NearestNeighbors(n_neighbors=5).fit(X)
distances, _ = nn.kneighbors(X)
distances = np.sort(distances[:, -1])

# DBSCAN doesn't have predict() — refit to predict on new data
# Or use HDBSCAN for a more robust, parameter-light alternative

from sklearn.cluster import AgglomerativeClustering
from scipy.cluster.hierarchy import dendrogram, linkage
import matplotlib.pyplot as plt

# Compute linkage matrix (for dendrogram)
Z = linkage(X, method="ward")  # 'ward','complete','average','single'

# Plot dendrogram to choose number of clusters
plt.figure(figsize=(10, 5))
dendrogram(Z, truncate_mode="level", p=5)
plt.axhline(y=Z[-4, 2], color="r", linestyle="--")  # cut at 4 clusters
plt.show()

# Agglomerative clustering with chosen k
agg = AgglomerativeClustering(
    n_clusters=4,
    linkage="ward",   # minimizes variance
    metric="euclidean",
)
labels = agg.fit_predict(X)

from sklearn.mixture import GaussianMixture

# GMM: probabilistic clustering — each cluster is a Gaussian
gmm = GaussianMixture(
    n_components=4,
    covariance_type="full",  # 'spherical','tied','diag','full'
    n_init=10,
    random_state=42,
)
gmm.fit(X)

labels = gmm.predict(X)
probs = gmm.predict_proba(X)   # soft clustering (responsibilities)

# Generative: sample new data
X_sampled, y_sampled = gmm.sample(100)

# Choose n_components via BIC (lower is better)
bics = []
for k in range(1, 11):
    g = GaussianMixture(n_components=k, n_init=10, random_state=42).fit(X)
    bics.append(g.bic(X))

from sklearn.ensemble import IsolationForest
from sklearn.svm import OneClassSVM
from sklearn.neighbors import LocalOutlierFactor

# Isolation Forest: isolates anomalies with shallow trees
iso = IsolationForest(contamination=0.05, random_state=42)
labels = iso.fit_predict(X)  # -1 anomaly, 1 normal
scores = iso.decision_function(X)  # lower = more anomalous

# One-Class SVM: learns a boundary around normal data
ocsvm = OneClassSVM(kernel="rbf", nu=0.05, gamma="scale")
labels = ocsvm.fit_predict(X_train_normal)

# Local Outlier Factor: density-based (good for local anomalies)
lof = LocalOutlierFactor(n_neighbors=20, contamination=0.05)
labels = lof.fit_predict(X)
# For new data: lof.negative_outlier_factor_ (lower = more anomalous)

from sklearn.linear_model import LinearRegression, Ridge, Lasso, ElasticNet

# Ordinary Least Squares
ols = LinearRegression()
ols.fit(X_train, y_train)

# Ridge (L2): shrinks coefficients, keeps all features
ridge = Ridge(alpha=1.0)
ridge.fit(X_train_scaled, y_train)

# Lasso (L1): shrinks some coefficients to 0 (feature selection)
lasso = Lasso(alpha=0.1, max_iter=10000)
lasso.fit(X_train_scaled, y_train)
print("Non-zero coefs:", (lasso.coef_ != 0).sum())

# ElasticNet: mix of L1 and L2
en = ElasticNet(alpha=0.1, l1_ratio=0.5)  # l1_ratio=1 -> pure Lasso
en.fit(X_train_scaled, y_train)

from sklearn.metrics import (mean_squared_error, mean_absolute_error,
    r2_score, mean_absolute_percentage_error)

y_pred = model.predict(X_test)

mse  = mean_squared_error(y_test, y_pred)         # sensitive to outliers
rmse = np.sqrt(mse)                                # same units as y
mae  = mean_absolute_error(y_test, y_pred)         # robust to outliers
r2   = r2_score(y_test, y_pred)                    # 1=perfect, 0=mean, <0=worse
mape = mean_absolute_percentage_error(y_test, y_pred)  # % error

print(f"RMSE: {rmse:.3f}")
print(f"MAE:  {mae:.3f}")
print(f"R2:   {r2:.3f}")
print(f"MAPE: {mape:.1%}")

# Adjusted R2 accounts for number of features
n, p = X_test.shape
adj_r2 = 1 - (1 - r2) * (n - 1) / (n - p - 1)

from sklearn.tree import DecisionTreeRegressor
from sklearn.ensemble import RandomForestRegressor, GradientBoostingRegressor

# Single tree (piecewise constant predictions)
dt = DecisionTreeRegressor(max_depth=5, min_samples_leaf=10,
                            random_state=42)
dt.fit(X_train, y_train)

# Random forest regressor (averages trees)
rf = RandomForestRegressor(n_estimators=200, max_features="sqrt",
                            n_jobs=-1, random_state=42)
rf.fit(X_train, y_train)

# Gradient boosting regressor
gb = GradientBoostingRegressor(n_estimators=200, learning_rate=0.1,
                                max_depth=3, random_state=42)
gb.fit(X_train, y_train)

# Quantile regression — predict ranges instead of points
from sklearn.ensemble import GradientBoostingRegressor
gbr_lower = GradientBoostingRegressor(loss="quantile", alpha=0.1)
gbr_upper = GradientBoostingRegressor(loss="quantile", alpha=0.9)
gbr_lower.fit(X_train, y_train)
gbr_upper.fit(X_train, y_train)

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import lasso_path

# Compute Lasso coefficients for many alpha values
alphas, coefs, _ = lasso_path(X_train_scaled, y_train,
                               alphas=np.logspace(-3, 1, 100))

# Plot: each line is a coefficient vs alpha
plt.figure(figsize=(8, 5))
for coef in coefs:
    plt.plot(alphas, coef)
plt.xscale("log")
plt.xlabel("alpha (regularization strength)")
plt.ylabel("coefficient")
plt.title("Lasso regularization path")
plt.axhline(0, color="k", lw=0.5)
plt.show()

import matplotlib.pyplot as plt
import seaborn as sns

residuals = y_test - y_pred

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# 1. Residuals vs predicted — should look random around 0
axes[0].scatter(y_pred, residuals, alpha=0.5)
axes[0].axhline(0, color="r", linestyle="--")
axes[0].set_xlabel("Predicted"); axes[0].set_ylabel("Residual")

# 2. Residual distribution — should be ~normal
sns.histplot(residuals, kde=True, ax=axes[1])
axes[1].set_xlabel("Residual")

plt.show()

# Patterns in residuals suggest model issues:
# - U-shape: missing non-linearity (add polynomial features)
# - Funnel: heteroscedasticity (try log-transforming y)
# - Skew: try transforming y or use a different loss

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (confusion_matrix, classification_report,
    roc_auc_score, precision_recall_curve)

clf = LogisticRegression(max_iter=1000)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
y_proba = clf.predict_proba(X_test)[:, 1]  # P(class=1)

# Confusion matrix: TN FP / FN TP
cm = confusion_matrix(y_test, y_pred)
print(cm)
print(classification_report(y_test, y_pred))

# Adjust threshold (default 0.5) to trade precision for recall
y_pred_custom = (y_proba > 0.3).astype(int)  # higher recall

# ROC-AUC: probability ranking quality (threshold-independent)
print("ROC-AUC:", roc_auc_score(y_test, y_proba))

from sklearn.metrics import (confusion_matrix, precision_score,
    recall_score, f1_score)

# Binary: confusion_matrix returns [[TN, FP], [FN, TP]]
tn, fp, fn, tp = confusion_matrix(y_test, y_pred).ravel()

accuracy  = (tp + tn) / (tp + tn + fp + fn)
precision = tp / (tp + fp)    # of predicted positives, how many correct
recall    = tp / (tp + fn)    # of actual positives, how many caught
f1        = 2 * precision * recall / (precision + recall)  # harmonic mean
specificity = tn / (tn + fp)  # true negative rate

# sklearn shortcuts
print("Precision:", precision_score(y_test, y_pred))
print("Recall:    ", recall_score(y_test, y_pred))
print("F1:        ", f1_score(y_test, y_pred))

# Multi-class: average='macro' (unweighted) or 'weighted' (by class size)
print("Macro F1:", f1_score(y_test, y_pred, average="macro"))

from sklearn.metrics import (roc_curve, precision_recall_curve,
    auc, RocCurveDisplay, PrecisionRecallDisplay)
import matplotlib.pyplot as plt

# ROC curve: TPR (recall) vs FPR at all thresholds
fpr, tpr, thresholds = roc_curve(y_test, y_proba)
roc_auc = auc(fpr, tpr)

# PR curve: precision vs recall (better for imbalanced data)
prec, rec, thresholds = precision_recall_curve(y_test, y_proba)
pr_auc = auc(rec, prec)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))
RocCurveDisplay.from_predictions(y_test, y_proba, ax=axes[0])
PrecisionRecallDisplay.from_predictions(y_test, y_proba, ax=axes[1])
axes[0].set_title(f"ROC (AUC={roc_auc:.3f})")
axes[1].set_title(f"PR (AUC={pr_auc:.3f})")
plt.show()

from sklearn.linear_model import LogisticRegression
from sklearn.metrics import (cohen_kappa_score, log_loss,
    classification_report, confusion_matrix, ConfusionMatrixDisplay)

# Multiclass: logistic regression uses softmax by default
clf = LogisticRegression(multi_class="multinomial", max_iter=1000)
clf.fit(X_train, y_train)
y_pred = clf.predict(X_test)
y_proba = clf.predict_proba(X_test)

# Per-class report
print(classification_report(y_test, y_pred))

# Confusion matrix (now NxN)
disp = ConfusionMatrixDisplay.from_predictions(
    y_test, y_pred, display_labels=clf.classes_,
    cmap="Blues", normalize="true")  # normalize per true class

# Multiclass metrics
print("Cohen's kappa:", cohen_kappa_score(y_test, y_pred))
print("Log loss:     ", log_loss(y_test, y_proba))

from sklearn.calibration import CalibratedClassifierCV, calibration_curve
import matplotlib.pyplot as plt

# Many classifiers output poorly-calibrated probabilities
# (e.g. Random Forest, SVM). Calibrate them.

# Method 1: isotonic regression (more flexible, needs more data)
# Method 2: sigmoid (Platt scaling) — parametric, works with less data
calibrated = CalibratedClassifierCV(rf, method="isotonic", cv=5)
calibrated.fit(X_train, y_train)
y_proba = calibrated.predict_proba(X_test)[:, 1]

# Calibration curve: predicted prob vs observed frequency
prob_true, prob_pred = calibration_curve(y_test, y_proba, n_bins=10)
plt.plot([0, 1], [0, 1], "k--")
plt.plot(prob_pred, prob_true, "s-")
plt.xlabel("Predicted probability"); plt.ylabel("Observed frequency")
plt.show()

from sklearn.cluster import KMeans
from sklearn.metrics import silhouette_score, silhouette_samples
import matplotlib.pyplot as plt

# Elbow method: inertia (sum of squared distances to centroid)
inertias = []
for k in range(1, 11):
    km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X)
    inertias.append(km.inertia_)

plt.plot(range(1, 11), inertias, "o-")
plt.xlabel("k"); plt.ylabel("Inertia"); plt.title("Elbow method")
plt.show()

# Silhouette: -1 (bad) to 1 (well-clustered), 0 = on boundary
silhouettes = []
for k in range(2, 11):
    labels = KMeans(n_clusters=k, n_init=10, random_state=42).fit_predict(X)
    silhouettes.append(silhouette_score(X, labels))

from sklearn.metrics import silhouette_samples
import matplotlib.pyplot as plt
import numpy as np

k = 4
km = KMeans(n_clusters=k, n_init=10, random_state=42).fit(X)
labels = km.labels_
sil_values = silhouette_samples(X, labels)

fig, ax = plt.subplots(figsize=(8, 6))
y_lower = 10
for i in range(k):
    ith = np.sort(sil_values[labels == i])
    y_upper = y_lower + len(ith)
    ax.fill_betweenx(np.arange(y_lower, y_upper), 0, ith, alpha=0.7)
    ax.text(-0.05, (y_lower + y_upper) / 2, str(i))
    y_lower = y_upper + 10
ax.axvline(sil_values.mean(), color="red", linestyle="--")
ax.set_xlabel("Silhouette coefficient"); ax.set_ylabel("Cluster")
plt.show()

import hdbscan  # pip install hdbscan

# HDBSCAN: DBSCAN that automatically chooses eps via hierarchy
clusterer = hdbscan.HDBSCAN(
    min_cluster_size=5,   # min points to form a cluster
    min_samples=None,     # how conservative to be about noise
    metric="euclidean",
    cluster_selection_method="eom",  # 'leaf' for fine-grained
)
labels = clusterer.fit_predict(X)

# Soft clustering: probability of cluster membership
probs = clusterer.probabilities_

# Compare clusters across parameter settings
clusterer.condensed_tree_.plot()

from sklearn.metrics import (silhouette_score, davies_bouldin_score,
    calinski_harabasz_score, adjusted_rand_score, normalized_mutual_info_score)

# Internal metrics (no ground truth)
print("Silhouette:        ", silhouette_score(X, labels))      # higher = better
print("Davies-Bouldin:    ", davies_bouldin_score(X, labels))   # lower = better
print("Calinski-Harabasz: ", calinski_harabasz_score(X, labels)) # higher = better

# External metrics (when you have ground-truth labels)
true_labels = y  # e.g. the actual iris species
print("Adjusted Rand:    ", adjusted_rand_score(true_labels, labels))
print("NMI:              ", normalized_mutual_info_score(true_labels, labels))

# ARI and NMI: 0 = random, 1 = perfect match
# Both are invariant to label permutations (cluster 1 vs cluster 2 doesn't matter)

from sklearn.decomposition import PCA, UMAP  # umap-learn package
from sklearn.manifold import TSNE
from sklearn.cluster import KMeans

# Reduce dimensions BEFORE clustering high-dim data
# PCA: linear, fast, preserves global structure
pca = PCA(n_components=10, random_state=42).fit(X)
X_pca = pca.transform(X)

# UMAP: non-linear, fast, preserves local+global structure
reducer = UMAP(n_components=5, random_state=42)
X_umap = reducer.fit_transform(X)

# t-SNE: visualization only (no transform on new data)
X_tsne = TSNE(n_components=2, perplexity=30, random_state=42).fit_transform(X)

# Cluster on the embeddings
labels = KMeans(n_clusters=4, n_init=10, random_state=42).fit_predict(X_umap)

# Plot 2D embedding colored by cluster
plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=labels, cmap="tab10", s=10)
plt.show()

from sklearn.decomposition import PCA

# PCA finds orthogonal axes of maximum variance
pca = PCA(n_components=2, random_state=42)
X_pca = pca.fit_transform(X_scaled)   # scaling matters!

# Explained variance: how much info each PC retains
print("Explained variance ratio:", pca.explained_variance_ratio_)
print("Total retained:", pca.explained_variance_ratio_.sum())

# Choose n_components to retain 95% of variance
pca_full = PCA(n_components=0.95, random_state=42).fit(X_scaled)
print("Components for 95%:", pca_full.n_components_)
X_reduced = pca_full.transform(X_scaled)

# Reconstruct (lossy) from compressed representation
X_reconstructed = pca_full.inverse_transform(X_reduced)

from sklearn.manifold import TSNE
import matplotlib.pyplot as plt

# t-SNE is for 2D/3D visualization, NOT as preprocessing
tsne = TSNE(
    n_components=2,
    perplexity=30,         # try 5-50; ~ sqrt(n_samples) is a heuristic
    learning_rate="auto",
    n_iter=1000,
    init="pca",            # better than random init
    random_state=42,
)
X_tsne = tsne.fit_transform(X_scaled)

plt.scatter(X_tsne[:, 0], X_tsne[:, 1], c=y, cmap="tab10", s=10)
plt.colorbar(); plt.show()

from umap import UMAP  # pip install umap-learn

# UMAP: faster than t-SNE, preserves more global structure,
# and supports transform on new data (good for ML pipelines)
reducer = UMAP(
    n_components=2,
    n_neighbors=15,    # smaller = local, larger = global
    min_dist=0.1,      # how tightly points pack
    metric="euclidean",
    random_state=42,
)
X_umap = reducer.fit_transform(X_scaled)

# Can transform new data (unlike t-SNE)
X_new_umap = reducer.transform(X_new)

# Use in a pipeline as a preprocessor before clustering/classification
from sklearn.pipeline import Pipeline
pipe = Pipeline([
    ("umap", UMAP(n_components=5, random_state=42)),
    ("clf",  RandomForestClassifier(random_state=42)),
])

from sklearn.discriminant_analysis import LinearDiscriminantAnalysis

# LDA is SUPERVISED dimensionality reduction — uses class labels
# to find axes that best separate classes
lda = LinearDiscriminantAnalysis(n_components=2)
X_lda = lda.fit_transform(X_train, y_train)  # note: y_train passed in

# Also a classifier (Gaussian model per class)
y_pred = lda.predict(X_test)

# Max n_components = n_classes - 1
# For 3 classes, you get at most 2 LDA components

# Compare: PCA ignores class labels (unsupervised);
# LDA uses labels to maximize class separability

from sklearn.decomposition import TruncatedSVD

# PCA doesn't work on sparse matrices (it needs to center data)
# Use TruncatedSVD instead — works directly on sparse TF-IDF, etc.

from scipy.sparse import csr_matrix
X_sparse = csr_matrix(X_tfidf)  # e.g. from TfidfVectorizer

svd = TruncatedSVD(n_components=100, random_state=42)
X_reduced = svd.fit_transform(X_sparse)

# Often combined with TF-IDF for LSA (Latent Semantic Analysis)
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.decomposition import TruncatedSVD

vectorizer = TfidfVectorizer(max_features=10000)
lsa = Pipeline([
    ("tfidf", vectorizer),
    ("svd",   TruncatedSVD(n_components=100, random_state=42)),
])
X_topics = lsa.fit_transform(documents)

from sklearn.model_selection import (cross_val_score, cross_validate,
    StratifiedKFold, KFold, RepeatedStratifiedKFold)

# Basic 5-fold CV
scores = cross_val_score(model, X, y, cv=5, scoring="accuracy")
print(f"{scores.mean():.3f} +/- {scores.std():.3f}")

# Stratified (preserves class balance — use for classification)
cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)
scores = cross_val_score(model, X, y, cv=cv, scoring="f1_macro")

# Multiple metrics at once
results = cross_validate(model, X, y, cv=cv,
    scoring=["accuracy", "f1_macro", "roc_auc_ovr"],
    return_train_score=True)
print(results["test_f1_macro"].mean())

# Repeated CV for more stable estimates
rcv = RepeatedStratifiedKFold(n_splits=5, n_repeats=3, random_state=42)

from sklearn.model_selection import learning_curve
import matplotlib.pyplot as plt

train_sizes, train_scores, val_scores = learning_curve(
    model, X, y,
    train_sizes=np.linspace(0.1, 1.0, 10),
    cv=5, scoring="accuracy", n_jobs=-1, random_state=42,
)

train_mean = train_scores.mean(axis=1)
val_mean   = val_scores.mean(axis=1)

plt.plot(train_sizes, train_mean, "o-", label="Train")
plt.plot(train_sizes, val_mean,   "o-", label="Validation")
plt.xlabel("Training samples"); plt.ylabel("Accuracy")
plt.legend(); plt.show()

# Diagnose:
# - High train, low val  -> overfitting (high variance) -> more data or regularize
# - Both low and close   -> underfitting (high bias) -> more complex model/features
# - Gap closing with data -> more data helps

from sklearn.model_selection import validation_curve
import matplotlib.pyplot as plt

param_range = [1, 2, 4, 8, 16, 32, 64]
train_scores, val_scores = validation_curve(
    RandomForestClassifier(random_state=42),
    X, y,
    param_name="max_depth",
    param_range=param_range,
    cv=5, scoring="accuracy", n_jobs=-1,
)

plt.plot(param_range, train_scores.mean(axis=1), "o-", label="Train")
plt.plot(param_range, val_scores.mean(axis=1),   "o-", label="Validation")
plt.xlabel("max_depth"); plt.ylabel("Accuracy")
plt.legend(); plt.show()

# Sweet spot is where validation peaks before overfitting

from sklearn.metrics import classification_report

report = classification_report(y_test, y_pred, output_dict=True)

# Per-class metrics + macro/weighted averages
import pandas as pd
df_report = pd.DataFrame(report).T
print(df_report[["precision", "recall", "f1-score", "support"]])

# 'macro' avg     = unweighted mean across classes (treats all equally)
# 'weighted' avg  = support-weighted mean (reflects class imbalance)
# Focus on:
# - macro F1 if all classes matter equally (e.g. digit recognition)
# - weighted F1 if rare classes matter less
# - recall of the positive class if it's the costly one (fraud, disease)

# For imbalanced classification, also report:
# - Precision-Recall AUC
# - Brier score (probability calibration)

from sklearn.metrics import (mean_squared_error, mean_absolute_error,
    r2_score, explained_variance_score, max_error)

# RMSE — penalizes large errors; same units as y
rmse = mean_squared_error(y_test, y_pred, squared=False)

# MAE — robust to outliers; same units as y
mae = mean_absolute_error(y_test, y_pred)

# R^2 — fraction of variance explained (1.0 = perfect, 0 = predict mean)
r2 = r2_score(y_test, y_pred)

# Explained variance — like R^2 but ignores mean offset
evs = explained_variance_score(y_test, y_pred)

# Max error — single largest absolute error (worst-case)
max_err = max_error(y_test, y_pred)

# Custom asymmetric loss: penalize under-prediction more
def asymmetric_error(y_true, y_pred, under_penalty=2.0):
    error = y_true - y_pred
    return np.where(error > 0, under_penalty * error**2, error**2).mean()

from sklearn.model_selection import (KFold, StratifiedKFold,
    GroupKFold, TimeSeriesSplit, RepeatedKFold)

# Standard K-fold (regression or balanced classification)
kf = KFold(n_splits=5, shuffle=True, random_state=42)

# Stratified (classification — preserves class proportions)
skf = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

# Group K-fold (when samples are not independent — e.g. multiple rows per patient)
gkf = GroupKFold(n_splits=5)  # pass groups=patient_ids to CV

# Time series split (no shuffle — preserves temporal order)
tscv = TimeSeriesSplit(n_splits=5)

# Use with cross_val_score:
scores = cross_val_score(model, X, y, cv=gkf, groups=groups, scoring="f1")

from sklearn.model_selection import TimeSeriesSplit
import matplotlib.pyplot as plt
import numpy as np

tscv = TimeSeriesSplit(n_splits=5, test_size=100)
# gap parameter (sklearn 1.0+): prevent leakage from train to test
tscv = TimeSeriesSplit(n_splits=5, test_size=100, gap=10)

for fold, (train_idx, test_idx) in enumerate(tscv.split(X)):
    X_tr, X_te = X[train_idx], X[test_idx]
    y_tr, y_te = y[train_idx], y[test_idx]
    model.fit(X_tr, y_tr)
    score = model.score(X_te, y_te)
    print(f"Fold {fold}: train={len(train_idx)}, test={len(test_idx)}, score={score:.3f}")

# NEVER shuffle time-series data for CV — that leaks future into past
# Use gap to avoid leakage from autocorrelated near-future points

from sklearn.model_selection import GroupKFold, GroupShuffleSplit

# When the same entity appears in multiple rows (e.g. multiple X-rays per patient),
# random splits leak: the model memorizes patient-specific features
# instead of learning generalizable patterns.

# Solution: keep all rows of a group in the same split
gkf = GroupKFold(n_splits=5)
scores = cross_val_score(model, X, y, cv=gkf, groups=df["patient_id"])

# GroupShuffleSplit: random train/test split respecting groups
gss = GroupShuffleSplit(n_splits=5, test_size=0.2, random_state=42)
for train_idx, test_idx in gss.split(X, groups=df["patient_id"]):
    ...

# Always pass the groups argument to CV functions
# Common scenarios: patients, users, sessions, image batches

from sklearn.model_selection import (GridSearchCV, cross_val_score,
    StratifiedKFold)

# Outer loop: estimate generalization performance
# Inner loop: tune hyperparameters
# Without nesting, GridSearchCV's CV score is optimistic (it selected the best
# hyperparameters using the same data)

inner_cv = StratifiedKFold(n_splits=3, shuffle=True, random_state=42)
outer_cv = StratifiedKFold(n_splits=5, shuffle=True, random_state=42)

grid = GridSearchCV(
    RandomForestClassifier(random_state=42),
    param_grid={"max_depth": [3, 5, 10], "n_estimators": [50, 100]},
    cv=inner_cv, scoring="f1", n_jobs=-1,
)

nested_scores = cross_val_score(grid, X, y, cv=outer_cv, scoring="f1")
print(f"Nested CV F1: {nested_scores.mean():.3f} +/- {nested_scores.std():.3f}")

# Fit the final model on all data
grid.fit(X, y)
best_model = grid.best_estimator_

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import cross_val_score
from sklearn.impute import SimpleImputer

# RIGHT: wrap preprocessing in a Pipeline so each fold fits
# its own scaler/imputer on the training portion only
pipe = Pipeline([
    ("impute",  SimpleImputer(strategy="median")),
    ("scale",   StandardScaler()),
    ("model",   RandomForestClassifier(random_state=42)),
])

scores = cross_val_score(pipe, X, y, cv=5, scoring="f1")

# WRONG: fit scaler on the whole dataset, then CV
# (leaks test-fold statistics into training)
# scaler.fit(X)              # NO - includes future test folds
# X_scaled = scaler.transform(X)
# cross_val_score(model, X_scaled, y, cv=5)  # optimistic scores

# Pipelines also make deployment easier (one object, no leakage)

from sklearn.model_selection import GridSearchCV

# Exhaustive search over all combinations
param_grid = {
    "n_estimators": [50, 100, 200],
    "max_depth":    [3, 5, 10, None],
    "min_samples_leaf": [1, 5, 10],
}
# 3 * 4 * 3 = 36 combinations, * 5 folds = 180 fits

grid = GridSearchCV(
    RandomForestClassifier(random_state=42),
    param_grid=param_grid,
    cv=5,
    scoring="f1",
    n_jobs=-1,
    refit=True,         # refit best model on full data
    return_train_score=True,
)
grid.fit(X_train, y_train)

print("Best params: ", grid.best_params_)
print("Best CV F1:  ", grid.best_score_)
best_model = grid.best_estimator_

from sklearn.model_selection import RandomizedSearchCV
from scipy.stats import randint, loguniform

# Sample from distributions instead of enumerating all combos
# Much more efficient when the grid is large
param_dist = {
    "n_estimators":     randint(50, 500),     # discrete uniform
    "max_depth":        [3, 5, 10, None],
    "min_samples_leaf": randint(1, 20),
    "learning_rate":    loguniform(1e-3, 1e0),  # log-uniform
}

rand_search = RandomizedSearchCV(
    GradientBoostingClassifier(random_state=42),
    param_distributions=param_dist,
    n_iter=50,           # number of parameter settings sampled
    cv=5,
    scoring="f1",
    n_jobs=-1,
    random_state=42,
)
rand_search.fit(X_train, y_train)
print("Best params:", rand_search.best_params_)

import optuna

def objective(trial):
    # Suggest hyperparameters
    n_estimators = trial.suggest_int("n_estimators", 50, 500)
    max_depth    = trial.suggest_int("max_depth", 3, 20)
    lr           = trial.suggest_float("learning_rate", 1e-3, 1e0, log=True)
    min_samples  = trial.suggest_int("min_samples_leaf", 1, 20)

    model = GradientBoostingClassifier(
        n_estimators=n_estimators, max_depth=max_depth,
        learning_rate=lr, min_samples_leaf=min_samples,
        random_state=42,
    )
    # 5-fold CV F1
    score = cross_val_score(model, X_train, y_train, cv=5,
                            scoring="f1", n_jobs=-1).mean()
    return score  # Optuna maximizes by default

study = optuna.create_study(direction="maximize")
study.optimize(objective, n_trials=100, timeout=600)

print("Best trial:", study.best_trial.params)
print("Best F1:    ", study.best_trial.value)

from sklearn.experimental import enable_halving_search_cv  # noqa
from sklearn.model_selection import HalvingGridSearchCV, HalvingRandomSearchCV

# Allocate more resources to promising candidates iteratively
halving = HalvingGridSearchCV(
    RandomForestClassifier(random_state=42),
    param_grid={"max_depth": [3, 5, 10, None],
                "n_estimators": [50, 100, 200]},
    cv=5,
    scoring="f1",
    factor=3,           # eliminate bottom 2/3 each round
    aggressive_elimination=False,
    n_jobs=-1,
)
halving.fit(X_train, y_train)
print("Best params:", halving.best_params_)
print("Best score: ", halving.best_score_)

# Much faster than full GridSearchCV — evaluates many candidates
# cheaply first, then invests more resources in the promising ones.

# 1. Tune on a VALIDATION set or via CV — never on the test set
# 2. Use the metric you actually care about (scoring=...)
# 3. Set random_state for reproducibility
# 4. Start coarse (randomized, few iterations), then refine around the best
# 5. Don't over-tune — small improvements may be noise

# Inspect full results to understand hyperparameter sensitivity
import pandas as pd
results = pd.DataFrame(grid.cv_results_)
cols = ["param_max_depth", "param_n_estimators",
        "mean_test_score", "std_test_score", "rank_test_score"]
print(results[cols].sort_values("rank_test_score").head(10))

# Refit the best model on the FULL training set (automatic with refit=True)
# Then evaluate ONCE on the held-out test set
test_score = grid.best_estimator_.score(X_test, y_test)
print(f"Test score: {test_score:.3f}")

from sklearn.ensemble import BaggingClassifier, RandomForestClassifier

# Bagging: train multiple models on bootstrap samples (with replacement)
bag = BaggingClassifier(
    estimator=DecisionTreeClassifier(),  # any base estimator
    n_estimators=100,
    max_samples=1.0,         # fraction of samples per bootstrap
    max_features=1.0,        # fraction of features per estimator
    bootstrap=True,
    n_jobs=-1,
    random_state=42,
)
bag.fit(X_train, y_train)

# Random Forest is bagging + random feature subsets at each split
# Out-of-bag (OOB) score: evaluate each tree on samples it didn't see
rf = RandomForestClassifier(n_estimators=200, oob_score=True,
                             random_state=42)
rf.fit(X_train, y_train)
print("OOB score:", rf.oob_score_)  # ~ validation accuracy for free

from sklearn.ensemble import (AdaBoostClassifier,
    GradientBoostingClassifier, HistGradientBoostingClassifier)

# AdaBoost: weight misclassified samples higher each round
ada = AdaBoostClassifier(
    estimator=DecisionTreeClassifier(max_depth=1),
    n_estimators=200, learning_rate=0.5, random_state=42,
)
ada.fit(X_train, y_train)

# Gradient Boosting: fit each new tree to the residuals of the previous
gb = GradientBoostingClassifier(
    n_estimators=200, learning_rate=0.1, max_depth=3,
    subsample=0.8,        # stochastic GB — adds randomness
    random_state=42,
)
gb.fit(X_train, y_train)

# HistGradientBoosting: faster, native NaN handling, large-data friendly
hgb = HistGradientBoostingClassifier(
    max_iter=200, learning_rate=0.1, max_leaf_nodes=31,
    l2_regularization=0.0, early_stopping=True,
    random_state=42,
)
hgb.fit(X_train, y_train)

from sklearn.ensemble import StackingClassifier, VotingClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.svm import SVC

# Voting: simple average of predictions (hard=majority, soft=probabilities)
voting = VotingClassifier(
    estimators=[
        ("rf",  RandomForestClassifier(n_estimators=100, random_state=42)),
        ("svc", SVC(probability=True, random_state=42)),
        ("lr",  LogisticRegression(max_iter=1000)),
    ],
    voting="soft",     # use predict_proba
    weights=[2, 1, 1], # weight each model
)

# Stacking: train a meta-model on base model predictions
stacking = StackingClassifier(
    estimators=[
        ("rf",  RandomForestClassifier(n_estimators=100, random_state=42)),
        ("svc", SVC(probability=True, random_state=42)),
        ("lr",  LogisticRegression(max_iter=1000)),
    ],
    final_estimator=LogisticRegression(),  # meta-model
    cv=5,            # use CV predictions to train meta-model (avoids leakage)
    n_jobs=-1,
)
stacking.fit(X_train, y_train)

# Voting:  simple average of predictions
# Stacking: meta-model trained on CV predictions of base models
# Blending: like stacking but meta-model trained on a single holdout set

# Blending (manual implementation):
# 1. Split train into train' and holdout
# 2. Train base models on train'
# 3. Get predictions on holdout
# 4. Train meta-model on (holdout_predictions, holdout_y)
# 5. Retrain base models on full train for final ensemble

from sklearn.model_selection import train_test_split
X_tr, X_hold, y_tr, y_hold = train_test_split(X_train, y_train, test_size=0.3,
                                               random_state=42)

# Train base models on X_tr
rf.fit(X_tr, y_tr); svc.fit(X_tr, y_tr); lr.fit(X_tr, y_tr)

# Get holdout predictions as new features
import numpy as np
holdout_preds = np.column_stack([
    rf.predict_proba(X_hold)[:, 1],
    svc.predict_proba(X_hold)[:, 1],
    lr.predict_proba(X_hold)[:, 1],
])
meta_model = LogisticRegression().fit(holdout_preds, y_hold)

import numpy as np
from scipy.stats import spearmanr

# An ensemble only helps if models make DIFFERENT errors
# Measure prediction correlation between base models
preds = np.column_stack([rf.predict_proba(X_test)[:, 1],
                          svc.predict_proba(X_test)[:, 1],
                          lr.predict_proba(X_test)[:, 1]])
corr = spearmanr(preds).correlation
print("Prediction correlation matrix:")
print(corr)

# Low correlation = diverse models = better ensemble
# High correlation = redundant models = no benefit

# Ways to encourage diversity:
# 1. Different algorithms (tree + linear + nearest-neighbor)
# 2. Different feature subsets
# 3. Different hyperparameters
# 4. Different training samples (bagging)
# 5. Different random seeds (weakest form of diversity)

import torch

# Tensors: like numpy arrays but on GPU and with autograd
x = torch.randn(3, 4, requires_grad=True)
y = torch.randn(4, 2, requires_grad=True)
z = x @ y               # matrix multiply
loss = z.sum()

# Automatic differentiation
loss.backward()         # populate .grad
print(x.grad.shape)     # torch.Size([3, 4])
print(y.grad.shape)     # torch.Size([4, 2])

# Move to GPU
device = torch.device("cuda" if torch.cuda.is_available() else "cpu")
x = x.to(device)

# Detach from graph (e.g. for logging)
numpy_array = z.detach().cpu().numpy()

import torch
import torch.nn as nn
import torch.optim as optim

model = nn.Sequential(
    nn.Linear(784, 128), nn.ReLU(),
    nn.Linear(128, 64),  nn.ReLU(),
    nn.Linear(64, 10),
).to(device)

criterion = nn.CrossEntropyLoss()
optimizer = optim.Adam(model.parameters(), lr=1e-3)

for epoch in range(epochs):
    model.train()                            # enable dropout/batchnorm
    for X_batch, y_batch in train_loader:
        X_batch, y_batch = X_batch.to(device), y_batch.to(device)

        optimizer.zero_grad()                # clear old gradients
        logits = model(X_batch)              # forward
        loss = criterion(logits, y_batch)    # compute loss
        loss.backward()                      # backprop
        optimizer.step()                     # update weights

    # Validation
    model.eval()                             # disable dropout/batchnorm
    with torch.no_grad():                    # no gradient tracking
        val_loss, correct = 0, 0
        for X_batch, y_batch in val_loader:
            logits = model(X_batch.to(device))
            val_loss += criterion(logits, y_batch.to(device))
            correct += (logits.argmax(1).cpu() == y_batch).sum().item()
    print(f"Epoch {epoch}: val_loss={val_loss/len(val_loader):.4f}, "
          f"val_acc={correct/len(val_loader.dataset):.4f}")

import torch.nn as nn

# Linear (fully connected)
nn.Linear(in_features=128, out_features=64)

# Activations
nn.ReLU()        # default for hidden layers
nn.GELU()        # smoother, used in transformers
nn.Sigmoid()     # outputs in (0, 1)
nn.Softmax(dim=1)# multi-class probabilities (built into CrossEntropyLoss)

# Convolution (images)
nn.Conv2d(in_channels=3, out_channels=32, kernel_size=3, padding=1)
nn.MaxPool2d(kernel_size=2)   # downsample by 2
nn.BatchNorm2d(32)            # stabilize training
nn.Dropout(p=0.5)             # regularization

# Recurrent (sequences)
nn.LSTM(input_size=64, hidden_size=128, num_layers=2, batch_first=True)

# Transformer (attention)
nn.TransformerEncoderLayer(d_model=512, nhead=8, batch_first=True)

import torch.optim as optim

# SGD with momentum — robust default
optim.SGD(model.parameters(), lr=0.01, momentum=0.9, weight_decay=1e-4)

# Adam — adaptive, often faster convergence
optim.Adam(model.parameters(), lr=1e-3, weight_decay=1e-4)

# AdamW — proper weight decay (preferred over Adam for transformers)
optim.AdamW(model.parameters(), lr=1e-3, weight_decay=0.01)

# Learning rate scheduler — reduce LR over training
from torch.optim.lr_scheduler import CosineAnnealingLR, OneCycleLR
scheduler = CosineAnnealingLR(optimizer, T_max=epochs)
# or OneCycleLR for super-convergence (needs careful setup)

# In training loop:
optimizer.step()
scheduler.step()  # update LR each epoch (or batch for OneCycleLR)

# 1. Dropout — randomly zero out units during training
nn.Dropout(p=0.5)

# 2. Weight decay (L2 regularization) via optimizer
optim.AdamW(model.parameters(), lr=1e-3, weight_decay=0.01)

# 3. Early stopping — stop when validation loss stops improving
best_val, patience, counter = float("inf"), 5, 0
for epoch in range(epochs):
    val_loss = validate(model, val_loader)
    if val_loss < best_val:
        best_val, counter = val_loss, 0
        torch.save(model.state_dict(), "best.pt")  # save best
    else:
        counter += 1
        if counter >= patience:
            print("Early stopping at epoch", epoch)
            break

# 4. Data augmentation (especially for images/text)
# 5. Batch normalization (mild regularizing effect)
# 6. Reduce model capacity (fewer layers/units)

from sklearn.metrics import (precision_score, recall_score,
    f1_score, fbeta_score)

# Precision: TP / (TP + FP) — of predicted positives, how many are real
# Recall:    TP / (TP + FN) — of actual positives, how many we caught
# F1:        harmonic mean of precision and recall
# F-beta:    weighted (beta > 1 favors recall, < 1 favors precision)

precision = precision_score(y_test, y_pred)
recall    = recall_score(y_test, y_pred)
f1        = f1_score(y_test, y_pred)
f2        = fbeta_score(y_test, y_pred, beta=2.0)   # recall-weighted
f05       = fbeta_score(y_test, y_pred, beta=0.5)   # precision-weighted

# When to use which:
# - Precision matters when false positives are costly (spam filter)
# - Recall matters when false negatives are costly (cancer screening)
# - F2 when recall is more important; F0.5 when precision is more important

from sklearn.metrics import (roc_auc_score, average_precision_score,
    roc_curve, precision_recall_curve)

# ROC-AUC: probability that model ranks a random positive
#          above a random negative. Threshold-independent.
roc_auc = roc_auc_score(y_test, y_proba)

# PR-AUC (average precision): area under precision-recall curve.
# More informative than ROC-AUC on imbalanced data.
pr_auc = average_precision_score(y_test, y_proba)

# Multiclass ROC-AUC
roc_auc_ovr = roc_auc_score(y_test, y_proba_multi,
                            multi_class="ovr", average="macro")
roc_auc_ovo = roc_auc_score(y_test, y_proba_multi,
                            multi_class="ovo", average="macro")

# Find the optimal threshold (Youden's J = TPR - FPR)
fpr, tpr, thresholds = roc_curve(y_test, y_proba)
j_scores = tpr - fpr
best_threshold = thresholds[j_scores.argmax()]

from sklearn.metrics import (confusion_matrix,
    ConfusionMatrixDisplay)
import matplotlib.pyplot as plt

cm = confusion_matrix(y_test, y_pred)

# Display with sklearn
disp = ConfusionMatrixDisplay(confusion_matrix=cm,
                              display_labels=class_names)
fig, ax = plt.subplots(figsize=(8, 6))
disp.plot(cmap="Blues", ax=ax, values_format="d")  # or ".2%" for normalized
plt.title("Confusion matrix")
plt.show()

# Normalized per row (true class) — see per-class recall
cm_norm = cm.astype(float) / cm.sum(axis=1, keepdims=True)
disp = ConfusionMatrixDisplay(confusion_matrix=cm_norm,
                              display_labels=class_names)
disp.plot(cmap="Blues", values_format=".2%")

from sklearn.metrics import (mean_squared_error, mean_absolute_error,
    r2_score, median_absolute_error, max_error,
    explained_variance_score, d2_tweedie_score)

# RMSE (penalizes large errors, same units as y)
rmse = mean_squared_error(y_test, y_pred, squared=False)

# MAE (robust to outliers, same units as y)
mae = mean_absolute_error(y_test, y_pred)

# Median AE (extremely robust to outliers)
med_ae = median_absolute_error(y_test, y_pred)

# Max error (worst single prediction — useful for SLAs)
max_err = max_error(y_test, y_pred)

# R^2 (1=perfect, 0=mean prediction, <0=worse than mean)
r2 = r2_score(y_test, y_pred)

# Explained variance (like R^2 but ignores bias)
evs = explained_variance_score(y_test, y_pred)

# Sometimes the right metric isn't in sklearn — define your own.

# Example 1: Cost matrix (asymmetric costs)
def cost_weighted_error(y_true, y_pred, cost_fp=1, cost_fn=10):
    cm = confusion_matrix(y_true, y_pred)
    tn, fp, fn, tp = cm.ravel()
    return (cost_fp * fp + cost_fn * fn) / len(y_true)

# Example 2: Profit at K (top-K precision weighted by revenue)
def profit_at_k(y_true, y_proba, revenue_per_hit, k=100):
    top_k_idx = np.argsort(y_proba)[::-1][:k]
    hits = y_true[top_k_idx].sum()
    return hits * revenue_per_hit - k * cost_per_contact

# Example 3: Custom threshold optimization
def optimize_threshold_for_profit(y_true, y_proba, profit_matrix):
    best_t, best_profit = 0.5, -float("inf")
    for t in np.linspace(0.01, 0.99, 99):
        y_pred = (y_proba >= t).astype(int)
        profit = profit_matrix(y_true, y_pred)
        if profit > best_profit:
            best_t, best_profit = t, profit
    return best_t, best_profit

from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.impute import SimpleImputer
from sklearn.ensemble import RandomForestClassifier

pipe = Pipeline([
    ("impute", SimpleImputer(strategy="median")),
    ("scale",  StandardScaler()),
    ("clf",    RandomForestClassifier(random_state=42)),
])

# Use like any estimator
pipe.fit(X_train, y_train)
y_pred = pipe.predict(X_test)
score = pipe.score(X_test, y_test)

# Access steps
pipe.named_steps["clf"]
pipe[:-1].get_feature_names_out()  # see transformed feature names

# Pipelines prevent leakage: each CV fold fits its own preprocessing

from sklearn.compose import ColumnTransformer
from sklearn.preprocessing import (StandardScaler, OneHotEncoder)
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer

numeric_features     = ["age", "income", "score"]
categorical_features = ["city", "occupation"]

# Different preprocessing for different column types
preprocessor = ColumnTransformer(
    transformers=[
        ("num", Pipeline([
            ("impute", SimpleImputer(strategy="median")),
            ("scale",  StandardScaler()),
        ]), numeric_features),
        ("cat", Pipeline([
            ("impute", SimpleImputer(strategy="most_frequent")),
            ("onehot", OneHotEncoder(handle_unknown="ignore")),
        ]), categorical_features),
    ],
    remainder="drop",  # or 'passthrough' to keep other columns
)

full_pipe = Pipeline([
    ("preprocess", preprocessor),
    ("clf", RandomForestClassifier(random_state=42)),
])
full_pipe.fit(X_train, y_train)

from sklearn.model_selection import GridSearchCV

# Reference nested parameters with 'step__param' syntax
param_grid = {
    "preprocess__num__impute__strategy": ["mean", "median"],
    "clf__n_estimators": [50, 100, 200],
    "clf__max_depth":    [3, 5, 10, None],
}

grid = GridSearchCV(full_pipe, param_grid, cv=5, scoring="f1", n_jobs=-1)
grid.fit(X_train, y_train)
print("Best params:", grid.best_params_)

# Tune preprocessing AND model together — finds the best combo
# e.g. median imputation + 100 trees + max_depth=10

from sklearn.base import BaseEstimator, TransformerMixin
import numpy as np

class LogTransformer(BaseEstimator, TransformerMixin):
    def __init__(self, columns=None):
        self.columns = columns

    def fit(self, X, y=None):
        return self  # nothing to fit

    def transform(self, X):
        X = X.copy()
        if self.columns is None:
            return np.log1p(X)
        X[self.columns] = np.log1p(X[self.columns])
        return X

# Use in a pipeline
pipe = Pipeline([
    ("log",  LogTransformer(columns=["income", "price"])),
    ("scale", StandardScaler()),
    ("clf",   LogisticRegression()),
])

# For more complex logic, use FunctionTransformer or write a full class

import joblib

# Save the entire fitted pipeline
joblib.dump(full_pipe, "model.pkl")

# Later (e.g. in a web service)
loaded = joblib.load("model.pkl")
predictions = loaded.predict(new_data)

# Why joblib over pickle?
# - More efficient for numpy arrays (which sklearn models contain)
# - Handles objects with large numerical arrays via compression

# Version compatibility:
# - sklearn versions must match between save and load in production
# - Pin sklearn in requirements.txt: scikit-learn==1.4.2

# For cross-language deployment, use ONNX:
# import skl2onnx
# onx = skl2onnx.convert_sklearn(full_pipe, ...)

import joblib
import pickle

# joblib — preferred for sklearn (efficient with numpy arrays)
joblib.dump(model, "model.joblib", compress=3)
model = joblib.load("model.joblib")

# pickle — standard library, works for any picklable object
with open("model.pkl", "wb") as f:
    pickle.dump(model, f)
with open("model.pkl", "rb") as f:
    model = pickle.load(f)

# Compress to save disk space (joblib supports it natively)
joblib.dump(model, "model.joblib", compress=("xz", 3))

# Load partially for inspection without full deserialize
# (advanced, not commonly needed)

# Save metadata alongside the model artifact
import joblib, json
from datetime import datetime

metadata = {
    "model_type": "RandomForestClassifier",
    "sklearn_version": sklearn.__version__,
    "python_version": sys.version,
    "features": list(X_train.columns),
    "target": "churn",
    "metrics": {"f1": 0.842, "roc_auc": 0.913},
    "trained_at": datetime.utcnow().isoformat(),
    "data_hash": hashlib.md5(pd.util.hash_pandas_object(X_train).values).hexdigest(),
    "git_commit": subprocess.check_output(["git", "rev-parse", "HEAD"]).decode().strip(),
}

joblib.dump({"model": model, "metadata": metadata}, "model_v1.2.0.joblib")

# Convention: semantic versioning
# v1.0.0 -> initial release
# v1.1.0 -> retrained on new data (backward compatible)
# v2.0.0 -> breaking change (new features, new preprocessing)

# pip install skl2onnx onnxruntime
from skl2onnx import convert_sklearn
from skl2onnx.common.data_types import FloatTensorType
import onnxruntime as rt

# Convert sklearn pipeline to ONNX
initial_type = [("float_input", FloatTensorType([None, X_train.shape[1]]))]
onnx_model = convert_sklearn(full_pipe, initial_types=initial_type,
                              target_opset=15)

with open("model.onnx", "wb") as f:
    f.write(onnx_model.SerializeToString())

# Run inference in any language (Python, C++, Java, JS, Go, ...)
sess = rt.InferenceSession("model.onnx")
input_name = sess.get_inputs()[0].name
predictions = sess.run(None, {input_name: X_test.astype(np.float32)})[0]

# Benefits: no Python at inference, fast, cross-platform

# After deployment, monitor for:
# 1. Data drift: input distribution changes
# 2. Concept drift: P(y|X) changes
# 3. Prediction drift: output distribution changes
# 4. Performance degradation: metric drops over time

# Simple data drift check: PSI (Population Stability Index)
def psi(expected, actual, bins=10):
    expected_pct = np.histogram(expected, bins=bins)[0] / len(expected)
    actual_pct   = np.histogram(actual,   bins=bins)[0] / len(actual)
    # Avoid log(0)
    expected_pct = np.clip(expected_pct, 1e-4, None)
    actual_pct   = np.clip(actual_pct,   1e-4, None)
    return np.sum((actual_pct - expected_pct) * np.log(actual_pct / expected_pct))

# PSI < 0.1: stable, 0.1-0.25: minor drift, > 0.25: significant drift
for col in X.columns:
    score = psi(X_train[col], X_new[col])
    if score > 0.25:
        print(f"Drift detected in {col}: PSI={score:.3f}")

# 1. Scheduled retraining (simplest)
#    Retrain weekly/monthly regardless of performance
schedule = "0 2 * * 0"  # every Sunday at 2am

# 2. Performance-triggered retraining
#    Retrain when a metric drops below threshold
if current_f1 < production_f1 * 0.95:  # 5% drop
    trigger_retrain()

# 3. Drift-triggered retraining
#    Retrain when PSI exceeds threshold
if max_psi > 0.25:
    trigger_retrain()

# 4. Online learning (incremental updates)
from sklearn.linear_model import SGDClassifier
model = SGDClassifier()
for X_batch, y_batch in stream:
    model.partial_fit(X_batch, y_batch, classes=[0, 1])

# Always:
# - Validate new model beats old on recent data before promoting
# - Shadow deploy: run new model in parallel, compare predictions
# - Roll back automatically if metrics regress
# - Keep training data fresh (drop old samples)

from sklearn.linear_model import Lasso, Ridge, ElasticNet
import numpy as np

# L1 (Lasso): penalty = alpha * sum(|coef|)
# Drives some coefficients to exactly 0 -> sparse -> feature selection
lasso = Lasso(alpha=0.1)
lasso.fit(X_train_scaled, y_train)
print("Non-zero coefs:", (lasso.coef_ != 0).sum(), "/", len(lasso.coef_))

# L2 (Ridge): penalty = alpha * sum(coef^2)
# Shrinks coefficients smoothly toward 0, keeps all features
ridge = Ridge(alpha=1.0)
ridge.fit(X_train_scaled, y_train)

# ElasticNet: mix of L1 and L2
# penalty = alpha * (l1_ratio * L1 + (1 - l1_ratio) * L2)
en = ElasticNet(alpha=0.1, l1_ratio=0.5)
en.fit(X_train_scaled, y_train)

# Rule of thumb:
# - L1: many features, only a few matter (sparse)
# - L2: many features, all matter a little
# - ElasticNet: correlated features (Lasso picks one randomly)

from sklearn.ensemble import RandomForestClassifier
from sklearn.tree import DecisionTreeClassifier

# Decision tree regularization hyperparameters
tree = DecisionTreeClassifier(
    max_depth=5,             # limit tree depth (strongest lever)
    min_samples_split=20,    # min samples to split a node
    min_samples_leaf=10,     # min samples per leaf
    max_leaf_nodes=20,       # cap total leaves
    max_features="sqrt",     # features considered per split
    ccp_alpha=0.01,          # cost-complexity pruning
)

# Random forest adds more regularization via:
rf = RandomForestClassifier(
    n_estimators=200,
    max_depth=None,
    min_samples_leaf=5,
    max_features="sqrt",     # random feature subset per split
    bootstrap=True,          # bagging
    max_samples=0.8,         # subsample ratio
)

# For boosting (XGBoost), also use:
# - learning_rate (smaller = more regularization)
# - subsample, colsample_bytree (randomness)
# - reg_lambda, reg_alpha (L2/L1 on leaf weights)
# - gamma (minimum loss reduction to split)
# - early_stopping_rounds

import torch.nn as nn

model = nn.Sequential(
    nn.Linear(784, 256),
    nn.BatchNorm1d(256),     # normalize activations -> faster, stable training
    nn.ReLU(),
    nn.Dropout(0.5),         # randomly zero 50% of activations during training

    nn.Linear(256, 128),
    nn.BatchNorm1d(128),
    nn.ReLU(),
    nn.Dropout(0.3),         # less dropout closer to output

    nn.Linear(128, 10),      # no dropout before output (or use very low)
)

# IMPORTANT: dropout is OFF during inference (model.eval())
# BatchNorm uses batch stats in training and running stats in eval

# Common patterns:
# - Higher dropout in big layers, lower in small layers
# - No dropout right after input or right before output
# - BatchNorm before activation (debated; try both)

import torch

best_val_loss = float("inf")
patience, counter = 5, 0
best_state = None

for epoch in range(max_epochs):
    train_loss = train_one_epoch(model, train_loader, optimizer, criterion)
    val_loss   = validate(model, val_loader, criterion)

    if val_loss < best_val_loss:
        best_val_loss = val_loss
        counter = 0
        best_state = {k: v.clone() for k, v in model.state_dict().items()}
    else:
        counter += 1
        if counter >= patience:
            print(f"Early stopping at epoch {epoch}")
            break

# Restore the best model
model.load_state_dict(best_state)

# For XGBoost/LightGBM, use early_stopping_rounds:
# model.fit(X_train, y_train,
#           eval_set=[(X_val, y_val)],
#           early_stopping_rounds=20)

# Images (torchvision)
import torchvision.transforms as T

train_transform = T.Compose([
    T.RandomResizedCrop(224),
    T.RandomHorizontalFlip(),
    T.ColorJitter(brightness=0.2, contrast=0.2, saturation=0.2),
    T.RandomRotation(10),
    T.ToTensor(),
    T.Normalize(mean=[0.485, 0.456, 0.406], std=[0.229, 0.224, 0.225]),
])

# Text (NLPAug or simple Python)
import random
def augment_text(text):
    # Synonym replacement, random deletion, back-translation
    words = text.split()
    if random.random() < 0.3:
        idx = random.randrange(len(words))
        words[idx] = synonym(words[idx])
    if random.random() < 0.2:
        words = [w for w in words if random.random() > 0.1]
    return " ".join(words)

# Tabular (SMOTE for imbalance, mixup for regression)
# mixup: X_mix = alpha*X1 + (1-alpha)*X2, y_mix = alpha*y1 + (1-alpha)*y2