SEP Occurrence Forecasting — SHAP Demonstration#
The archived 49 predictors are anonymous and have no timestamps or event IDs. Every result below is therefore sample-level teaching evidence, not an event-aware forecast claim or a physical attribution.
Runtime dependency check#
import importlib.util
import subprocess
import sys
from pathlib import Path
REQUIRED_RUNTIME = {}
COLAB_EXTRAS = {'shap': 'shap'}
missing_required = [
package for module, package in REQUIRED_RUNTIME.items()
if importlib.util.find_spec(module) is None
]
missing_extras = [
package for module, package in COLAB_EXTRAS.items()
if importlib.util.find_spec(module) is None
]
if missing_extras and "google.colab" in sys.modules:
subprocess.check_call(
[sys.executable, "-m", "pip", "install", "-q", *missing_extras]
)
missing_extras = []
if missing_required or missing_extras:
missing = ", ".join(missing_required + missing_extras)
raise RuntimeError(
f"Missing notebook dependencies: {missing}. Locally run "
"`uv sync --group notebooks`; in Colab restart the runtime if an "
"installation cell just changed the environment."
)
print("runtime dependency check passed")
runtime dependency check passed
%matplotlib inline
Imports and deterministic configuration#
These tools handle the archived samples, preprocessing, imbalance-aware metrics, and the framework-neutral tree analysis.
import json
import os
import random
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.metrics import (
accuracy_score,
average_precision_score,
balanced_accuracy_score,
classification_report,
confusion_matrix,
precision_recall_curve,
roc_auc_score,
)
from sklearn.model_selection import train_test_split
from sklearn.preprocessing import StandardScaler
SEED = 42
random.seed(SEED)
np.random.seed(SEED)
Resolve the immutable archive#
The four supplied files contain the archived training and test samples. Their checksums are verified before any analysis is performed.
import hashlib
import os
from pathlib import Path
from urllib.parse import quote
from urllib.request import urlopen
DATASET_ID = 'sep-curated'
DATASET_FILES = {'x_train.pkl': ('data/sep-curated/x_train.pkl', 'e809bf00498633f509a223d61f9b0006e6ed1803f6de22118bcf654f2ce8ba3b'), 'x_test.pkl': ('data/sep-curated/x_test.pkl', '1d0c5f84713d4fde34d567cdb62e9081c4d723f6fef9abd543137376350d5955'), 'y_train.pkl': ('data/sep-curated/y_train.pkl', 'd7aa048f6b081a9fb1fc00dde19872c0f67ae5b4c8620daa5984b679f9f9dbdc'), 'y_test.pkl': ('data/sep-curated/y_test.pkl', 'd44c5af108bab2b19f5f8082548282edd8aee89d57e15469516de1ca3f400ee5')}
def file_sha256(path):
digest = hashlib.sha256()
with Path(path).open("rb") as stream:
for chunk in iter(lambda: stream.read(1024 * 1024), b""):
digest.update(chunk)
return digest.hexdigest()
def resolve_dataset():
resolved = {}
override = os.getenv("HELIO_DATA_DIR")
cache_root = Path(
os.getenv("HELIO_DATA_CACHE", Path.home() / ".cache" / "helio-data-methods")
) / "datasets" / DATASET_ID
for filename, (relative_path, checksum) in DATASET_FILES.items():
candidates = []
if override:
root = Path(override).expanduser()
candidates.extend([root / DATASET_ID / filename, root / filename])
for root in [Path.cwd(), *Path.cwd().parents]:
candidates.append(root / relative_path)
target = cache_root / filename
candidates.append(target)
match = next(
(
candidate
for candidate in candidates
if candidate.is_file() and file_sha256(candidate) == checksum
),
None,
)
if match is None:
target.parent.mkdir(parents=True, exist_ok=True)
ref = os.getenv("HELIO_DATA_REF", "main")
url = (
"https://raw.githubusercontent.com/SavvasRaptis/helio-data-methods/"
f"{quote(ref, safe='')}/{quote(relative_path, safe='/')}"
)
try:
with urlopen(url, timeout=120) as response, target.open("wb") as output:
while chunk := response.read(1024 * 1024):
output.write(chunk)
except Exception as exc:
target.unlink(missing_ok=True)
raise RuntimeError(
f"Could not retrieve {DATASET_ID}/{filename}. Check network "
"access or set HELIO_DATA_DIR to the archived data directory."
) from exc
if file_sha256(target) != checksum:
target.unlink(missing_ok=True)
raise ValueError(
f"Checksum mismatch for {DATASET_ID}/{filename}; "
"the invalid download was removed."
)
match = target
resolved[filename] = match
return resolved
dataset_files = resolve_dataset()
print("verified dataset:", DATASET_ID)
for name in dataset_files:
print(f" {name} (checksum verified)")
verified dataset: sep-curated
x_train.pkl: data/sep-curated/x_train.pkl
x_test.pkl: data/sep-curated/x_test.pkl
y_train.pkl: data/sep-curated/y_train.pkl
y_test.pkl: data/sep-curated/y_test.pkl
Preserve the supplied test set and split training samples#
The supplied test set remains untouched. Validation samples are drawn only from the supplied training set, and preprocessing is fitted without using the test samples.
x_supplied_train = pd.read_pickle(dataset_files["x_train.pkl"]).to_numpy(dtype=np.float32)
x_test = pd.read_pickle(dataset_files["x_test.pkl"]).to_numpy(dtype=np.float32)
y_supplied_train = (
pd.read_pickle(dataset_files["y_train.pkl"]).to_numpy().reshape(-1).astype(np.int64)
)
y_test = pd.read_pickle(dataset_files["y_test.pkl"]).to_numpy().reshape(-1).astype(np.int64)
feature_names = np.asarray([f"anonymous feature {i}" for i in range(x_test.shape[1])])
train_indices, validation_indices = train_test_split(
np.arange(len(y_supplied_train)),
test_size=0.15,
random_state=SEED,
stratify=y_supplied_train,
)
x_train_raw = x_supplied_train[train_indices]
y_train = y_supplied_train[train_indices]
x_validation_raw = x_supplied_train[validation_indices]
y_validation = y_supplied_train[validation_indices]
scaler = StandardScaler().fit(x_train_raw)
x_train = scaler.transform(x_train_raw).astype(np.float32)
x_validation = scaler.transform(x_validation_raw).astype(np.float32)
x_test_scaled = scaler.transform(x_test).astype(np.float32)
counts = np.bincount(y_train, minlength=2)
majority_class = int(np.argmax(counts))
majority_prediction = np.full_like(y_test, majority_class)
class_weights = len(y_train) / (2.0 * np.maximum(counts, 1))
print(
f"train={len(y_train):,}, validation={len(y_validation):,}, "
f"supplied test={len(y_test):,}, positive prevalence={y_train.mean():.4f}"
)
print("class weights:", dict(enumerate(class_weights.round(3))))
train=13,846, validation=2,444, supplied test=1,811, positive prevalence=0.0125
class weights: {0: 0.506, 1: 40.017}
Establish the majority-class baseline#
Because non-SEP samples dominate the archive, the majority-class result is a useful reminder that accuracy by itself is not sufficient.
def classification_evidence(y_true, probability, label):
prediction = (probability >= 0.5).astype(np.int64)
evidence = {
"accuracy": float(accuracy_score(y_true, prediction)),
"balanced_accuracy": float(balanced_accuracy_score(y_true, prediction)),
"roc_auc": float(roc_auc_score(y_true, probability)),
"pr_auc": float(average_precision_score(y_true, probability)),
"confusion_matrix": confusion_matrix(y_true, prediction, labels=[0, 1]).tolist(),
}
print(label, json.dumps(evidence, indent=2))
print(classification_report(y_true, prediction, digits=3, zero_division=0))
return evidence
majority_probability = np.full(len(y_test), float(majority_class))
majority_evidence = classification_evidence(
y_test, majority_probability, "majority-class baseline"
)
majority-class baseline {
"accuracy": 0.9872998343456654,
"balanced_accuracy": 0.5,
"roc_auc": 0.5,
"pr_auc": 0.012700165654334622,
"confusion_matrix": [
[
1788,
0
],
[
23,
0
]
]
}
precision recall f1-score support
0 0.987 1.000 0.994 1788
1 0.000 0.000 0.000 23
accuracy 0.987 1811
macro avg 0.494 0.500 0.497 1811
weighted avg 0.975 0.987 0.981 1811
Train the weighted boosted-tree model#
XGBoost provides a nonlinear, non-neural comparison using the same split and class weighting as the other SEP demonstrations.
from xgboost import XGBClassifier
ROUNDS = 300 # Reduce to 50 or 100 for a quicker run.
# Define the boosted-tree model used for comparison.
model = XGBClassifier(
n_estimators=ROUNDS,
max_depth=6,
learning_rate=0.05,
subsample=0.8,
colsample_bytree=0.8,
objective="binary:logistic",
eval_metric="logloss",
scale_pos_weight=float(counts[0] / max(counts[1], 1)),
random_state=SEED,
n_jobs=2,
)
model.fit(x_train, y_train, eval_set=[(x_validation, y_validation)], verbose=False)
probabilities = model.predict_proba(x_test_scaled)[:, 1]
Summarize anonymous feature influence with SHAP#
SHAP summarizes how strongly each anonymous column affects this fitted model. Without physical feature names, the plot describes model sensitivity rather than a physical attribution.
import warnings
warnings.filterwarnings("ignore", message="IProgress not found.*")
import shap
background_count = min(500, len(x_train))
explain_count = min(300, len(x_test_scaled))
# Compute model-attribution values for a representative sample of test rows.
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(x_test_scaled[:explain_count])
if isinstance(shap_values, list):
shap_values = shap_values[-1]
importance = np.abs(np.asarray(shap_values)).mean(axis=0)
order = np.argsort(importance)[-12:]
fig, ax = plt.subplots(figsize=(8, 5))
ax.barh(feature_names[order], importance[order])
ax.set(title="Mean absolute SHAP value (anonymous columns)", xlabel="mean |SHAP|")
plt.tight_layout()
plt.show()
probabilities = model.predict_proba(x_test_scaled)[:, 1]
model_evidence = classification_evidence(y_test, probabilities, "model")
print(
"HELIO_RESULT "
+ json.dumps(
{
"explained_samples": explain_count,
"shap_shape": list(np.asarray(shap_values).shape),
"balanced_accuracy": model_evidence["balanced_accuracy"],
},
sort_keys=True,
)
)
model {
"accuracy": 0.9900607399226946,
"balanced_accuracy": 0.8232905359400837,
"roc_auc": 0.9775800019453361,
"pr_auc": 0.5670328619716795,
"confusion_matrix": [
[
1778,
10
],
[
8,
15
]
]
}
precision recall f1-score support
0 0.996 0.994 0.995 1788
1 0.600 0.652 0.625 23
accuracy 0.990 1811
macro avg 0.798 0.823 0.810 1811
weighted avg 0.990 0.990 0.990 1811
HELIO_RESULT {"balanced_accuracy": 0.8232905359400837, "explained_samples": 300, "shap_shape": [300, 49]}