hackathon-ENGIN/severity_benchmark.py

662 lines
26 KiB
Python

"""Exploratory candidate-ranking benchmark for ENGIN fault severity.
EXPERIMENT STATUS: CANDIDATE SEARCH. Absolute end-to-end metrics from this
script do not describe ``final_pipeline.py`` because its frozen label model is
an older baseline. Its reusable transformers and winning severity candidate
are imported by the final pipeline. Use ``ml_polish_benchmark.py`` for final
end-to-end numbers.
The label model is frozen as combined spectral features + Logistic Regression.
Severity candidates compare raw spectra, engine-relative deviation features,
ordinal models, and fault-specific models. Every experiment uses complete held-
out engines and is evaluated on clean validation spectra and the same spectra
with 5% of measurements masked.
"""
from __future__ import annotations
import argparse
import json
from dataclasses import dataclass
from pathlib import Path
from typing import Iterable
import numpy as np
import pandas as pd
from sklearn.base import BaseEstimator
from sklearn.dummy import DummyClassifier
from sklearn.ensemble import ExtraTreesClassifier
from sklearn.linear_model import LogisticRegression, Ridge
from sklearn.metrics import accuracy_score
from sklearn.pipeline import make_pipeline as sklearn_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
from benchmark_grouped import (
FAULT_LABELS,
FREQ_COLS,
NOT_APPLICABLE,
SEVERITIES,
make_pipeline,
make_splits,
macro_f1,
ml_points,
raw_score,
validate_data,
)
from robustness_grouped import mask_spectral_cells, missing_scenario_name
DEFAULT_SEEDS = [7, 21, 42, 77, 123]
SEVERITY_TO_INT = {"male": 0, "srednie": 1, "duze": 2}
INT_TO_SEVERITY = np.asarray(SEVERITIES, dtype=object)
FEATURE_SETS = ("raw", "deviation", "all")
@dataclass(frozen=True)
class Candidate:
candidate_id: str
feature_set: str
estimator: str
include_fault_type: bool = False
per_fault: bool = False
CANDIDATES = [
Candidate("raw_extra_trees", "raw", "extra_trees"),
Candidate("deviation_extra_trees", "deviation", "extra_trees"),
Candidate(
"deviation_extra_trees_mf03",
"deviation",
"extra_trees_mf03",
),
Candidate("all_extra_trees", "all", "extra_trees", include_fault_type=True),
Candidate(
"all_extra_trees_leaf2",
"all",
"extra_trees_leaf2",
include_fault_type=True,
),
Candidate("all_logistic", "all", "logistic", include_fault_type=True),
Candidate("all_svc", "all", "svc", include_fault_type=True),
Candidate("all_ordinal_ridge", "all", "ordinal_ridge", include_fault_type=True),
Candidate(
"all_ordinal_logistic",
"all",
"ordinal_logistic",
include_fault_type=True,
),
Candidate("all_extra_trees_per_fault", "all", "extra_trees", per_fault=True),
]
CANDIDATE_BY_ID = {candidate.candidate_id: candidate for candidate in CANDIDATES}
class SeverityFeatures:
"""Create severity-oriented features without using validation labels.
Engine-relative values use a leave-one-cylinder-out median reference. Input
to ``transform`` must therefore contain every cylinder of each engine.
``diagnostic_label`` is optional and may contain only the model-predicted
label at validation or test time.
"""
def __init__(self, feature_set: str, include_fault_type: bool = False) -> None:
if feature_set not in FEATURE_SETS:
raise ValueError(f"Unknown feature_set={feature_set!r}")
self.feature_set = feature_set
self.include_fault_type = include_fault_type
def fit(
self, X: pd.DataFrame, y: Iterable[str] | None = None
) -> "SeverityFeatures":
self._validate(X)
spectra = X[FREQ_COLS].apply(pd.to_numeric, errors="coerce")
fallback = spectra.median(axis=0).to_numpy(dtype=float)
if np.isnan(fallback).any():
raise ValueError("A frequency column is entirely NaN in training.")
self.fallback_medians_ = fallback
return self
def transform(self, X: pd.DataFrame) -> np.ndarray:
self._validate(X)
if not hasattr(self, "fallback_medians_"):
raise RuntimeError("SeverityFeatures must be fitted before transform().")
raw = self._clean_spectra(X)
reference = self._leave_one_out_engine_median(raw, X["engine_id"].to_numpy())
relative = raw - reference
absolute_relative = np.abs(relative)
ratio_delta = raw / np.maximum(np.abs(reference), 1e-6) - 1.0
if self.feature_set == "raw":
features = raw
elif self.feature_set == "deviation":
features = np.hstack(
[relative, absolute_relative, ratio_delta, self._summary_features(raw, relative, ratio_delta)]
)
else:
features = np.hstack(
[raw, relative, absolute_relative, ratio_delta, self._summary_features(raw, relative, ratio_delta)]
)
if self.include_fault_type:
if "diagnostic_label" not in X.columns:
raise ValueError("diagnostic_label is required when include_fault_type=True")
labels = X["diagnostic_label"].to_numpy(dtype=object)
one_hot = np.column_stack([labels == label for label in FAULT_LABELS]).astype(float)
features = np.hstack([features, one_hot])
if np.isnan(features).any() or np.isinf(features).any():
raise ValueError("Severity features contain NaN or infinity.")
return features
def fit_transform(
self, X: pd.DataFrame, y: Iterable[str] | None = None
) -> np.ndarray:
return self.fit(X).transform(X)
def _clean_spectra(self, X: pd.DataFrame) -> np.ndarray:
spectra = X[FREQ_COLS].apply(pd.to_numeric, errors="coerce")
spectra = spectra.interpolate(axis=1, limit_direction="both")
spectra = spectra.fillna(pd.Series(self.fallback_medians_, index=FREQ_COLS))
return spectra.to_numpy(dtype=float)
@staticmethod
def _leave_one_out_engine_median(raw: np.ndarray, groups: np.ndarray) -> np.ndarray:
reference = np.empty_like(raw)
for engine_id in pd.unique(groups):
positions = np.flatnonzero(groups == engine_id)
if len(positions) < 2:
raise ValueError(
f"Engine {engine_id!r} has fewer than two cylinders in transform()."
)
engine_values = raw[positions]
for local_position, global_position in enumerate(positions):
keep = np.arange(len(positions)) != local_position
reference[global_position] = np.median(engine_values[keep], axis=0)
return reference
@staticmethod
def _summary_features(
raw: np.ndarray, relative: np.ndarray, ratio_delta: np.ndarray
) -> np.ndarray:
absolute_relative = np.abs(relative)
gradients = np.diff(raw, axis=1)
frequency = np.arange(raw.shape[1], dtype=float)
centered_frequency = frequency - frequency.mean()
slope = raw @ centered_frequency / np.sum(centered_frequency**2)
# ``np.trapezoid`` was added after NumPy 1.26, which is still allowed by
# requirements.txt. Keep the fresh-install path compatible with both
# NumPy 1.x and 2.x.
if hasattr(np, "trapezoid"):
spectral_auc = np.trapezoid(raw, axis=1)
else: # pragma: no cover - exercised only with NumPy 1.x
spectral_auc = np.trapz(raw, axis=1)
columns = [
raw.mean(axis=1),
raw.std(axis=1),
raw.min(axis=1),
raw.max(axis=1),
np.ptp(raw, axis=1),
spectral_auc,
raw.argmax(axis=1) / 20.0,
raw.argmin(axis=1) / 20.0,
slope,
np.abs(gradients).mean(axis=1),
np.abs(gradients).max(axis=1),
relative.mean(axis=1),
relative.std(axis=1),
relative.min(axis=1),
relative.max(axis=1),
absolute_relative.mean(axis=1),
absolute_relative.max(axis=1),
np.sqrt(np.mean(relative**2, axis=1)),
ratio_delta.mean(axis=1),
ratio_delta.std(axis=1),
]
for start, stop in ((0, 5), (5, 10), (10, 15), (15, 21)):
columns.extend(
[
raw[:, start:stop].mean(axis=1),
relative[:, start:stop].mean(axis=1),
absolute_relative[:, start:stop].mean(axis=1),
]
)
return np.column_stack(columns)
@staticmethod
def _validate(X: pd.DataFrame) -> None:
required = {"engine_id", *FREQ_COLS}
missing = sorted(required.difference(X.columns))
if missing:
raise ValueError(f"Missing columns: {missing}")
class OrdinalRidge(BaseEstimator):
"""Weighted ridge regression rounded to the three ordered severities."""
def __init__(self, alpha: float = 10.0) -> None:
self.alpha = alpha
def fit(self, X: np.ndarray, y: Iterable[str]) -> "OrdinalRidge":
y_array = np.asarray(list(y), dtype=object)
y_numeric = np.asarray([SEVERITY_TO_INT[value] for value in y_array], dtype=float)
self.scaler_ = StandardScaler().fit(X)
counts = pd.Series(y_array).value_counts()
weights = np.asarray([len(y_array) / counts[value] for value in y_array])
self.model_ = Ridge(alpha=self.alpha).fit(
self.scaler_.transform(X), y_numeric, sample_weight=weights
)
return self
def predict(self, X: np.ndarray) -> np.ndarray:
numeric = self.model_.predict(self.scaler_.transform(X))
levels = np.clip(np.rint(numeric), 0, 2).astype(int)
return INT_TO_SEVERITY[levels]
class OrdinalLogistic(BaseEstimator):
"""Two cumulative balanced logistic classifiers for ordered severity."""
def __init__(self, C: float = 1.0, random_state: int = 42) -> None:
self.C = C
self.random_state = random_state
def fit(self, X: np.ndarray, y: Iterable[str]) -> "OrdinalLogistic":
y_numeric = np.asarray([SEVERITY_TO_INT[value] for value in y], dtype=int)
self.scaler_ = StandardScaler().fit(X)
scaled = self.scaler_.transform(X)
self.threshold_models_: list[float | LogisticRegression] = []
for threshold in (0, 1):
target = (y_numeric > threshold).astype(int)
if np.unique(target).size == 1:
self.threshold_models_.append(float(target[0]))
else:
model = LogisticRegression(
C=self.C,
class_weight="balanced",
max_iter=5_000,
random_state=self.random_state + threshold,
).fit(scaled, target)
self.threshold_models_.append(model)
return self
def predict(self, X: np.ndarray) -> np.ndarray:
scaled = self.scaler_.transform(X)
probabilities = []
for model in self.threshold_models_:
if isinstance(model, float):
probabilities.append(np.full(len(X), model))
else:
probabilities.append(model.predict_proba(scaled)[:, 1])
above_male = probabilities[0]
above_srednie = np.minimum(probabilities[1], above_male)
levels = (above_male >= 0.5).astype(int) + (above_srednie >= 0.5).astype(int)
return INT_TO_SEVERITY[levels]
class PerFaultEstimator:
"""Fit a separate severity estimator for each predicted fault type."""
def __init__(self, estimator_name: str, random_state: int, n_jobs: int) -> None:
self.estimator_name = estimator_name
self.random_state = random_state
self.n_jobs = n_jobs
def fit(
self,
X: np.ndarray,
y: np.ndarray,
fault_labels: np.ndarray,
) -> "PerFaultEstimator":
self.global_model_ = make_estimator(
self.estimator_name, self.random_state, self.n_jobs
).fit(X, y)
self.models_: dict[str, object] = {}
for offset, fault_label in enumerate(FAULT_LABELS, start=1):
mask = fault_labels == fault_label
if not mask.any():
continue
if np.unique(y[mask]).size == 1:
model = DummyClassifier(strategy="most_frequent").fit(X[mask], y[mask])
else:
model = make_estimator(
self.estimator_name,
self.random_state + offset,
self.n_jobs,
).fit(X[mask], y[mask])
self.models_[fault_label] = model
return self
def predict(self, X: np.ndarray, fault_labels: np.ndarray) -> np.ndarray:
result = self.global_model_.predict(X).astype(object)
for fault_label, model in self.models_.items():
mask = fault_labels == fault_label
if mask.any():
result[mask] = model.predict(X[mask])
return result
def make_estimator(estimator_name: str, random_state: int, n_jobs: int):
if estimator_name == "extra_trees":
return ExtraTreesClassifier(
n_estimators=400,
class_weight="balanced",
max_features="sqrt",
min_samples_leaf=1,
n_jobs=n_jobs,
random_state=random_state,
)
if estimator_name == "extra_trees_leaf2":
return ExtraTreesClassifier(
n_estimators=400,
class_weight="balanced",
max_features=0.7,
min_samples_leaf=2,
n_jobs=n_jobs,
random_state=random_state,
)
if estimator_name == "extra_trees_mf03":
return ExtraTreesClassifier(
n_estimators=800,
class_weight="balanced",
max_features=0.3,
min_samples_leaf=1,
n_jobs=n_jobs,
random_state=random_state,
)
if estimator_name == "logistic":
return sklearn_pipeline(
StandardScaler(),
LogisticRegression(
C=1.0,
class_weight="balanced",
max_iter=5_000,
random_state=random_state,
),
)
if estimator_name == "svc":
return sklearn_pipeline(
StandardScaler(),
SVC(C=2.0, kernel="rbf", class_weight="balanced"),
)
if estimator_name == "ordinal_ridge":
return OrdinalRidge(alpha=10.0)
if estimator_name == "ordinal_logistic":
return OrdinalLogistic(C=1.0, random_state=random_state)
raise ValueError(f"Unknown estimator: {estimator_name}")
def parse_args() -> argparse.Namespace:
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--data", type=Path, default=Path("val.csv"))
parser.add_argument("--output-dir", type=Path, default=Path("severity_outputs"))
parser.add_argument("--seeds", nargs="+", type=int, default=DEFAULT_SEEDS)
parser.add_argument("--n-splits", type=int, default=5)
parser.add_argument("--missing-rate", type=float, default=0.05)
parser.add_argument("--mask-random-state", type=int, default=2026)
parser.add_argument("--model-random-state", type=int, default=42)
parser.add_argument("--n-jobs", type=int, default=-1)
parser.add_argument(
"--candidates",
nargs="+",
choices=sorted(CANDIDATE_BY_ID),
default=[candidate.candidate_id for candidate in CANDIDATES],
)
return parser.parse_args()
def prepare_labeled_frame(df: pd.DataFrame, labels: np.ndarray) -> pd.DataFrame:
result = df.copy()
result["diagnostic_label"] = labels
return result
def fit_candidate(
candidate: Candidate,
X_train_full: pd.DataFrame,
train_fault_mask: np.ndarray,
y_severity_fault: np.ndarray,
y_fault_label: np.ndarray,
random_state: int,
n_jobs: int,
):
transformer = SeverityFeatures(
candidate.feature_set, include_fault_type=candidate.include_fault_type
).fit(X_train_full)
train_features = transformer.transform(X_train_full)[train_fault_mask]
if candidate.per_fault:
estimator = PerFaultEstimator(candidate.estimator, random_state, n_jobs).fit(
train_features, y_severity_fault, y_fault_label
)
else:
estimator = make_estimator(candidate.estimator, random_state, n_jobs).fit(
train_features, y_severity_fault
)
return transformer, estimator
def predict_candidate(
candidate: Candidate,
transformer: SeverityFeatures,
estimator,
X_valid_full: pd.DataFrame,
predicted_labels: np.ndarray,
) -> np.ndarray:
features = transformer.transform(X_valid_full)
if candidate.per_fault:
return estimator.predict(features, predicted_labels)
return estimator.predict(features)
def main() -> None:
args = parse_args()
df_clean = pd.read_csv(args.data).reset_index(drop=True)
validate_data(df_clean, args.n_splits)
df_masked, mask_manifest = mask_spectral_cells(
df_clean, args.missing_rate, args.mask_random_state
)
masked_scenario = missing_scenario_name(args.missing_rate)
scenarios = {"clean": df_clean, masked_scenario: df_masked}
candidates = [CANDIDATE_BY_ID[candidate_id] for candidate_id in args.candidates]
args.output_dir.mkdir(parents=True, exist_ok=True)
mask_manifest.to_csv(args.output_dir / "mask_manifest.csv", index=False)
y_label = df_clean["label"].reset_index(drop=True)
y_severity = df_clean["severity"].reset_index(drop=True)
true_fault = y_label.isin(FAULT_LABELS).to_numpy()
all_run_rows: list[dict[str, object]] = []
all_fold_rows: list[dict[str, object]] = []
for cv_seed in args.seeds:
print(f"Running severity CV seed {cv_seed} ...", flush=True)
splits = make_splits(
df_clean,
splitter_name="stratified-group",
n_splits=args.n_splits,
random_state=cv_seed,
)
oof_label = {
scenario: np.full(len(df_clean), "", dtype=object) for scenario in scenarios
}
oof_severity = {
(candidate.candidate_id, scenario): np.full(
len(df_clean), NOT_APPLICABLE, dtype=object
)
for candidate in candidates
for scenario in scenarios
}
for fold, (train_idx, valid_idx) in enumerate(splits, start=1):
X_train_clean = df_clean.iloc[train_idx]
y_train_label = y_label.iloc[train_idx]
label_pipeline = make_pipeline(
"logistic", "combined", args.model_random_state + fold, args.n_jobs
).fit(X_train_clean, y_train_label)
scenario_predicted_labels: dict[str, np.ndarray] = {}
for scenario, scenario_df in scenarios.items():
predicted = label_pipeline.predict(scenario_df.iloc[valid_idx])
scenario_predicted_labels[scenario] = predicted
oof_label[scenario][valid_idx] = predicted
train_fault_mask = y_train_label.isin(FAULT_LABELS).to_numpy()
train_with_labels = prepare_labeled_frame(
X_train_clean, y_train_label.to_numpy(dtype=object)
)
y_severity_fault = (
y_severity.iloc[train_idx].to_numpy(dtype=object)[train_fault_mask]
)
y_fault_label = y_train_label.to_numpy(dtype=object)[train_fault_mask]
fold_true_fault = y_label.iloc[valid_idx].isin(FAULT_LABELS).to_numpy()
for candidate_offset, candidate in enumerate(candidates):
transformer, estimator = fit_candidate(
candidate=candidate,
X_train_full=train_with_labels,
train_fault_mask=train_fault_mask,
y_severity_fault=y_severity_fault,
y_fault_label=y_fault_label,
random_state=args.model_random_state + 1000 + 20 * fold + candidate_offset,
n_jobs=args.n_jobs,
)
for scenario, scenario_df in scenarios.items():
predicted_labels = scenario_predicted_labels[scenario]
valid_with_labels = prepare_labeled_frame(
scenario_df.iloc[valid_idx], predicted_labels
)
severity_pred = predict_candidate(
candidate,
transformer,
estimator,
valid_with_labels,
predicted_labels,
)
emitted = np.full(len(valid_idx), NOT_APPLICABLE, dtype=object)
predicted_fault = np.isin(predicted_labels, FAULT_LABELS)
emitted[predicted_fault] = severity_pred[predicted_fault]
oof_severity[(candidate.candidate_id, scenario)][valid_idx] = emitted
fold_label_f1 = macro_f1(
y_label.iloc[valid_idx], predicted_labels
)
fold_severity = float(
accuracy_score(
y_severity.iloc[valid_idx].to_numpy()[fold_true_fault],
emitted[fold_true_fault],
)
)
fold_score = raw_score(fold_label_f1, fold_severity)
all_fold_rows.append(
{
"cv_seed": cv_seed,
"fold": fold,
"scenario": scenario,
"candidate_id": candidate.candidate_id,
"valid_fault_rows": int(fold_true_fault.sum()),
"macro_f1": fold_label_f1,
"severity_accuracy_submission": fold_severity,
"raw_score": fold_score,
"ml_points": ml_points(fold_score),
}
)
for scenario in scenarios:
label_f1 = macro_f1(y_label, oof_label[scenario])
for candidate in candidates:
emitted = oof_severity[(candidate.candidate_id, scenario)]
severity_accuracy = float(
accuracy_score(
y_severity.to_numpy()[true_fault], emitted[true_fault]
)
)
score = raw_score(label_f1, severity_accuracy)
all_run_rows.append(
{
"cv_seed": cv_seed,
"scenario": scenario,
"candidate_id": candidate.candidate_id,
"feature_set": candidate.feature_set,
"estimator": candidate.estimator,
"include_fault_type": candidate.include_fault_type,
"per_fault": candidate.per_fault,
"macro_f1_oof": label_f1,
"severity_accuracy_submission_oof": severity_accuracy,
"raw_score_oof": score,
"ml_points_oof": ml_points(score),
}
)
oof_frame = df_clean[
["engine_id", "cylinder", "label", "severity"]
].copy()
oof_frame["predicted_label"] = oof_label[scenario]
oof_frame["predicted_severity_submission"] = emitted
oof_frame.to_csv(
args.output_dir
/ f"oof__{candidate.candidate_id}__seed_{cv_seed}__{scenario}.csv",
index=False,
)
run_frame = pd.DataFrame(all_run_rows)
fold_frame = pd.DataFrame(all_fold_rows)
summary_frame = (
run_frame.groupby(["scenario", "candidate_id"], sort=False)
.agg(
seeds=("cv_seed", "nunique"),
macro_f1_mean=("macro_f1_oof", "mean"),
severity_mean=("severity_accuracy_submission_oof", "mean"),
severity_std=("severity_accuracy_submission_oof", "std"),
severity_min=("severity_accuracy_submission_oof", "min"),
severity_max=("severity_accuracy_submission_oof", "max"),
raw_score_mean=("raw_score_oof", "mean"),
raw_score_std=("raw_score_oof", "std"),
raw_score_min=("raw_score_oof", "min"),
ml_points_mean=("ml_points_oof", "mean"),
ml_points_min=("ml_points_oof", "min"),
)
.reset_index()
.sort_values(
["scenario", "raw_score_mean", "severity_mean"],
ascending=[True, False, False],
)
)
run_frame.to_csv(args.output_dir / "severity_runs.csv", index=False)
fold_frame.to_csv(args.output_dir / "severity_folds.csv", index=False)
summary_frame.to_csv(args.output_dir / "severity_summary.csv", index=False)
metadata = {
"data": str(args.data),
"cv_seeds": args.seeds,
"n_splits": args.n_splits,
"splitter": "StratifiedGroupKFold",
"missing_rate": args.missing_rate,
"mask_random_state": args.mask_random_state,
"label_pipeline": "combined features + Logistic Regression (frozen)",
"candidates": args.candidates,
}
(args.output_dir / "metadata.json").write_text(
json.dumps(metadata, ensure_ascii=False, indent=2), encoding="utf-8"
)
display_columns = [
"scenario",
"candidate_id",
"severity_mean",
"severity_std",
"severity_min",
"raw_score_mean",
"raw_score_std",
"ml_points_mean",
]
print("\nSeverity ranking across CV seeds:")
print(
summary_frame[display_columns].to_string(
index=False, float_format="{:.4f}".format
)
)
print(f"\nDetailed outputs saved to: {args.output_dir.resolve()}")
if __name__ == "__main__":
main()