Implémentation de l’algorithme APS

import numpy as np
import polars as pl

# Polars display options
pl.Config.set_tbl_hide_column_data_types(True)
pl.Config.set_tbl_hide_dataframe_shape(True)
pl.Config.set_float_precision(2)

# Load preprocessed data
df = pl.read_parquet("../../data/forest_cover_clean.parquet")

print(f"Dataset shape: {df.shape}")
df.head()

from sklearn.model_selection import train_test_split

y = df.drop_in_place("Cover_Type")
X = df

# 60% train, 20% calibration, 20% test split
X_train, X_temp, y_train, y_temp = train_test_split(
    X, y, test_size=0.4, random_state=42
)
X_calib, X_test, y_calib, y_test = train_test_split(
    X_temp, y_temp, test_size=0.5, random_state=42
)

print(f"Train set size (60%): {X_train.shape[0]}")
print(f"Calibration set size (20%): {X_calib.shape[0]}")
print(f"Test set size (20%): {X_test.shape[0]}")

from sklearn.ensemble import HistGradientBoostingClassifier

model = HistGradientBoostingClassifier(
    loss="log_loss",
    learning_rate=0.1,
    max_iter=2000,
    max_depth=None,
    max_leaf_nodes=127,
    min_samples_leaf=50,
    l2_regularization=0.5,
    categorical_features="from_dtype",
    early_stopping=True,
    n_iter_no_change=20,
    validation_fraction=0.1,
    class_weight="balanced",  # Covertype classes are heavily imbalanced
    random_state=42,
)

model.fit(X_train, y_train);

from sklearn.metrics import classification_report

y_pred = model.predict(X_test)
print(classification_report(y_test, y_pred))

# 1. Predict probabilities on calibration set
y_calib_pred_proba = model.predict_proba(X_calib)
classes = model.classes_

# Helper to get sorting indices (descending probability)
# argsort returns indices that sort the array, we flip specifically for descending
sorted_indices = np.argsort(y_calib_pred_proba, axis=1)[:, ::-1]

# Sort the probabilities
sorted_probas = np.take_along_axis(y_calib_pred_proba, sorted_indices, axis=1)

# Compute cumulative probabilities
cum_probas = np.cumsum(sorted_probas, axis=1)

# 2. Compute nonconformity scores (Cumulative probability of true class)

# Map class names to indices
class_to_idx = {c: i for i, c in enumerate(classes)}
y_calib_indices = np.array([class_to_idx[yi] for yi in y_calib])

# Find the rank of the true class for each sample
# We compare the sorted indices to the true class index
# np.where returns (row_indices, rank_indices)
rows, ranks = np.where(sorted_indices == y_calib_indices[:, None])

# The score is the cumulative probability at the rank of the true class
nonconformity_scores = cum_probas[rows, ranks]

# 3. Compute quantile
alpha = 0.01
target_coverage = 1 - alpha
n = len(nonconformity_scores)
q_level = np.ceil((n + 1) * (1 - alpha)) / n
qhat = np.quantile(nonconformity_scores, q_level, method="linear")

f"APS calibration qhat (α={alpha}): {qhat:.6f}"

import matplotlib.pyplot as plt

plt.figure(figsize=(10, 5))
plt.hist(
    nonconformity_scores,
    bins=50,
    alpha=0.7,
    edgecolor="k",
    label="Cumulative Proba of True Class",
)
plt.axvline(
    qhat, color="red", linestyle="--", label=f"Quantile ({1 - alpha:.2f}) = {qhat:.6f}"
)
plt.xlabel("APS Score (Cumulative Probability)")
plt.ylabel("Frequency")
plt.title("Distribution of Nonconformity Scores on Calibration Set")
plt.legend()
plt.show()

# 4. Predict on test set
y_test_pred_proba = model.predict_proba(X_test)

# Sort probabilities (descending)
test_sorted_indices = np.argsort(y_test_pred_proba, axis=1)[:, ::-1]
test_sorted_probas = np.take_along_axis(y_test_pred_proba, test_sorted_indices, axis=1)
test_cum_probas = np.cumsum(test_sorted_probas, axis=1)

# 5. Build prediction sets
prediction_sets = []

for i in range(len(X_test)):
    # Find indices where cumulative proba is <= qhat
    # We must include the first class that crosses the threshold to ensure coverage
    # Logical trick: include class if the *previous* cumulative sum was < qhat
    # Or simply: take all classes where cum_probas <= qhat, plus the first one that exceeds it.

    # Identify which ranks are inside the quantile
    mask = test_cum_probas[i] <= qhat

    # If the first class is already > qhat, we must minimally include the top-1 class
    # Otherwise, we might get an empty set if qhat is very small (though unlikely with alpha=0.01)
    if not np.any(mask):
        included_ranks = np.array([0])  # Keep at least top 1
    else:
        # Get indices of True values
        true_indices = np.where(mask)[0]
        last_idx = true_indices[-1]

        # We generally include the one that crosses the boundary to be safe
        # (Though strictly it depends on exact definition, usually include boundary)
        included_ranks = (
            np.arange(last_idx + 1)  # +1 to include 0 to last_idx
            if last_idx < len(classes) - 1
            else np.arange(len(classes))
        )

    # Map back sorted ranks to original class names
    original_indices = test_sorted_indices[i, included_ranks]
    pset = classes[original_indices].tolist()
    prediction_sets.append(pset)

# Create results DataFrame
df_test_results = pl.DataFrame(
    {
        "True Class": y_test,
        "Prediction Set": prediction_sets,
    }
).with_columns(pl.col("Prediction Set").list.len().alias("Cardinality"))

# Add features for analysis
df_test_full = pl.concat([df_test_results, X_test], how="horizontal")
df_test_full.head()

df_test_full.describe()

coverage = (
    pl.col("Prediction Set")
    .list.contains(pl.col("True Class"))
    .mean()
    .mul(100)
    .round(2)
    .alias("Coverage")
)
avg_set_size = pl.col("Cardinality").mean().round(2).alias("Average Cardinality")
singleton_rate = (pl.col("Cardinality") == 1).mean().round(2).alias("Singleton Rate")

print(f"Coverage: {df_test_full.select(coverage).item()}% (target: 99%)")
print(f"Average Cardinality: {df_test_full.select(avg_set_size).item()} elements")
print(f"Singleton Rate: {df_test_full.select(singleton_rate).item()}%")

df_test_full.group_by("True Class").agg(
    coverage,
    avg_set_size,
    singleton_rate,
    pl.len().alias("Count"),
).sort("True Class")

df_test_full.group_by("Cardinality").agg(
    coverage,
    pl.len().alias("Count"),
).sort("Cardinality")