SEP Occurrence Forecasting — Repeated Sample-Level Validation#

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 = {'xgboost': 'xgboost'}
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}

Repeated stratified validation within supplied training samples#

Repeated stratified folds show how sensitive the tree model is to the particular sample-level partition. Each fold preserves the strong class imbalance.

from sklearn.model_selection import RepeatedStratifiedKFold
from xgboost import XGBClassifier

N_SPLITS = 5  # Reduce to 2 for a quicker validation run.
N_REPEATS = 3  # Reduce to 1 for a quicker validation run.
ROUNDS = 200  # Reduce to 20 or 50 for a quicker validation run.
# Repeat the sample-level validation with the same class balance in each fold.
features = scaler.fit_transform(x_supplied_train).astype(np.float32)
folds = RepeatedStratifiedKFold(
    n_splits=N_SPLITS,
    n_repeats=N_REPEATS,
    random_state=SEED,
)
scores = []
for fold, (fold_train, fold_validation) in enumerate(
    folds.split(features, y_supplied_train), start=1
):
    fold_scaler = StandardScaler().fit(x_supplied_train[fold_train])
    fold_x_train = fold_scaler.transform(x_supplied_train[fold_train])
    fold_x_validation = fold_scaler.transform(x_supplied_train[fold_validation])
    fold_y_train = y_supplied_train[fold_train]
    fold_counts = np.bincount(fold_y_train, minlength=2)
    fold_model = XGBClassifier(
        n_estimators=ROUNDS,
        max_depth=6,
        learning_rate=0.05,
        subsample=0.8,
        colsample_bytree=0.8,
        eval_metric="logloss",
        scale_pos_weight=float(fold_counts[0] / max(fold_counts[1], 1)),
        random_state=SEED + fold,
        n_jobs=2,
    )
    fold_model.fit(fold_x_train, fold_y_train)
    fold_probability = fold_model.predict_proba(fold_x_validation)[:, 1]
    score = {
        "balanced_accuracy": balanced_accuracy_score(
            y_supplied_train[fold_validation], fold_probability >= 0.5
        ),
        "roc_auc": roc_auc_score(y_supplied_train[fold_validation], fold_probability),
        "pr_auc": average_precision_score(
            y_supplied_train[fold_validation], fold_probability
        ),
    }
    scores.append(score)
    print(f"fold {fold}:", score)
summary = {
    metric: {
        "mean": float(np.mean([score[metric] for score in scores])),
        "std": float(np.std([score[metric] for score in scores])),
    }
    for metric in scores[0]
}
print(json.dumps(summary, indent=2))
print("HELIO_RESULT " + json.dumps({"folds": len(scores), **summary}, sort_keys=True))
fold 1: {'balanced_accuracy': 0.7334602237414543, 'roc_auc': 0.9660348042262275, 'pr_auc': 0.5011370837209519}
fold 2: {'balanced_accuracy': 0.8086155997513984, 'roc_auc': 0.9876243008079552, 'pr_auc': 0.5723439974557973}
fold 3: {'balanced_accuracy': 0.800215319529633, 'roc_auc': 0.959422883007195, 'pr_auc': 0.5706807886420412}
fold 4: {'balanced_accuracy': 0.7761359242439176, 'roc_auc': 0.9451541733322213, 'pr_auc': 0.531038820978802}
fold 5: {'balanced_accuracy': 0.8375778069251008, 'roc_auc': 0.962683002646004, 'pr_auc': 0.59705881689583}
fold 6: {'balanced_accuracy': 0.8083048477315102, 'roc_auc': 0.9323259788688627, 'pr_auc': 0.520676086117499}
fold 7: {'balanced_accuracy': 0.7339263517712865, 'roc_auc': 0.9666330018645121, 'pr_auc': 0.4659807674226477}
fold 8: {'balanced_accuracy': 0.8016141383048894, 'roc_auc': 0.9664662577617383, 'pr_auc': 0.618327308568547}
fold 9: {'balanced_accuracy': 0.7877093489616898, 'roc_auc': 0.9594001379864591, 'pr_auc': 0.4779911259991}
fold 10: {'balanced_accuracy': 0.8738523241620355, 'roc_auc': 0.9661554091450146, 'pr_auc': 0.5714040495447025}
fold 11: {'balanced_accuracy': 0.757838719701678, 'roc_auc': 0.9444530764449969, 'pr_auc': 0.4872122944826352}
fold 12: {'balanced_accuracy': 0.8339263517712865, 'roc_auc': 0.9634400248601616, 'pr_auc': 0.6587634350719273}
fold 13: {'balanced_accuracy': 0.8124104414808525, 'roc_auc': 0.9727969551998908, 'pr_auc': 0.5609591054370808}
fold 14: {'balanced_accuracy': 0.7762913485522794, 'roc_auc': 0.9429251613001054, 'pr_auc': 0.469969841328524}
fold 15: {'balanced_accuracy': 0.8006815924547184, 'roc_auc': 0.9606511141269324, 'pr_auc': 0.515014038562412}
{
  "balanced_accuracy": {
    "mean": 0.7961706892722485,
    "std": 0.0364701993572829
  },
  "roc_auc": {
    "mean": 0.9597444187718852,
    "std": 0.013230413054756505
  },
  "pr_auc": {
    "mean": 0.5412371706818999,
    "std": 0.05567085833529893
  }
}
HELIO_RESULT {"balanced_accuracy": {"mean": 0.7961706892722485, "std": 0.0364701993572829}, "folds": 15, "pr_auc": {"mean": 0.5412371706818999, "std": 0.05567085833529893}, "roc_auc": {"mean": 0.9597444187718852, "std": 0.013230413054756505}}

How to read the validation result#

These folds estimate sensitivity to a sample-level partition. Without event identifiers, they cannot detect leakage between measurements from the same physical event.