A small Transformer, the same seismic question

Francisco Plaza-Vega · STATSEI14 · 13 October 2026

Does a more flexible sequence encoder improve the probability of at least one catalogue M ≥ 5 event in the next seven days, using 30 past days from the same Chilean region?

Eight-minute discussion: fixed question (2 min), read the attention model (2 min), paired evidence (2 min), scientific interpretation (2 min). Follow the prepared outputs during the session; running all six fits is optional. The complete coding exercise remains follow-up work.

We compare a 16-unit LSTM and one compact Transformer encoder block using three catalogue channels. The ETAS-inspired channel is omitted here: architecture changes while representation and task remain fixed. This does not reproduce a published Transformer study.

Scope: the 2021–2025 test results have already been inspected. This is an exploratory retrospective comparison, not a new blind test or operational forecasting. The frozen catalogue contains revised event records; historical availability is not reconstructed.

1. Follow prepared results or make a new run

The participant notebook defaults to prepared pilot results. Change RUN_TRAINING to True for a separate fixed-protocol run. Keep EVALUATE_TEST=True only for the unchanged original comparison. For new seeds or a manually changed architecture, set EVALUATE_TEST=False; all cells can still run in order. The solved notebook was executed in a fresh local CPU kernel with training enabled; its outputs identify that run.

Locally, open this notebook inside the extracted course pack. In Colab, upload course-pack.zip when requested. The smaller data-pack.zip lacks the reference features and pilot results. No live catalogue download occurs.

If imports fail, run this in a separate cell and restart:

%pip install tensorflow==2.19.1 keras==3.15.1 numpy==2.1.3 pandas==2.3.3 scikit-learn==1.9.1 matplotlib

The environment files record the tested local versions. Local CPU timing does not certify Colab performance.

In [1]:
RUN_TRAINING = True  # True fits new models; False follows prepared pilot evidence.
SEEDS = [7, 17, 27]     # Predeclare three seeds before looking at their results.
EVALUATE_TEST = True   # False for every extension; True only for the original fixed rerun.
REFERENCE_ARCHITECTURE = True  # Set False before manually changing an architecture.
EXTENSION_NOTE = ""    # Describe any manual architecture change, such as omitting positions.
import os
os.environ["CUDA_VISIBLE_DEVICES"] = "-1"
os.environ["TF_CPP_MIN_LOG_LEVEL"] = "2"
os.environ["TF_NUM_INTRAOP_THREADS"] = "2"
os.environ["TF_NUM_INTEROP_THREADS"] = "2"
os.environ["OMP_NUM_THREADS"] = "2"
from pathlib import Path
import datetime, hashlib, json, platform, time, uuid, zipfile
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sklearn
from sklearn.preprocessing import StandardScaler
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import brier_score_loss, log_loss
import tensorflow as tf
from tensorflow import keras
from IPython.display import display
NOTEBOOK_STARTED = time.perf_counter()
try:
    tf.config.threading.set_intra_op_parallelism_threads(2)
    tf.config.threading.set_inter_op_parallelism_threads(2)
except RuntimeError:
    print("Restart the runtime to apply reference thread settings.")
print("Mode:", "NEW FIXED-PROTOCOL RUN" if RUN_TRAINING else "PREPARED PILOT RESULTS")
print("TensorFlow:", tf.__version__, "Keras:", keras.__version__, "Python:", platform.python_version())

Mode: NEW FIXED-PROTOCOL RUN
TensorFlow: 2.19.1 Keras: 3.15.1 Python: 3.12.14
In [2]:
ROOT = None
required_relative = ["config/course.json", "results/reference-run/features_raw.npz",
    "results/reference-run/window_metadata.csv", "results/reference-run/normalization.json",
    "results/transformer-pilot/protocol.json", "results/transformer-pilot/metrics.csv",
    "results/transformer-pilot/fits.csv", "results/transformer-pilot/predictions.csv",
    "results/transformer-pilot/runtime.json"]
for candidate in [Path.cwd(), *Path.cwd().parents]:
    if all((candidate / name).is_file() for name in required_relative):
        ROOT = candidate
        break
if ROOT is None:
    try:
        from google.colab import files
    except ImportError:
        raise FileNotFoundError("Open this notebook inside the complete extracted course-pack.zip.")
    print("Upload course-pack.zip, including reference data and prepared pilot results.")
    uploaded = files.upload()
    archives = [name for name in uploaded if name.lower().endswith(".zip")]
    if len(archives) != 1: raise ValueError("Upload exactly one course-pack.zip.")
    destination = (Path.cwd() / "transformer-course-data").resolve()
    destination.mkdir(exist_ok=True)
    with zipfile.ZipFile(archives[0]) as archive:
        for member in archive.infolist():
            if not (destination / member.filename).resolve().is_relative_to(destination):
                raise ValueError("Unsafe ZIP path.")
        archive.extractall(destination)
    candidates = [p.parents[1] for p in destination.rglob("config/course.json")]
    candidates = [p for p in candidates if all((p / name).is_file() for name in required_relative)]
    if len(candidates) != 1: raise ValueError("The complete pack must contain one course folder.")
    ROOT = candidates[0]
PILOT = ROOT / "results/transformer-pilot"
protocol = json.loads((PILOT / "protocol.json").read_text(encoding="utf-8-sig"))
config = json.loads((ROOT / "config/course.json").read_text(encoding="utf-8-sig"))
PILOT_SEEDS = [7, 17, 27]  # Immutable provenance of the prepared pilot.
FEATURES = ["log_count", "max_magnitude_excess", "has_event"]
assert protocol["seeds"] == PILOT_SEEDS and protocol["features"] == FEATURES
assert len(SEEDS) == 3 and len(set(SEEDS)) == 3 and all(isinstance(seed, int) for seed in SEEDS)
if SEEDS != PILOT_SEEDS or not REFERENCE_ARCHITECTURE:
    assert RUN_TRAINING and not EVALUATE_TEST, "Extensions require RUN_TRAINING=True and EVALUATE_TEST=False."
if not REFERENCE_ARCHITECTURE:
    assert EXTENSION_NOTE.strip(), "Describe the manual architecture change in EXTENSION_NOTE."
if not RUN_TRAINING:
    assert SEEDS == PILOT_SEEDS and REFERENCE_ARCHITECTURE, "Prepared mode uses the original pilot configuration."
allowed_splits = ["validation", "test"] if EVALUATE_TEST else ["validation"]
pilot_metrics = pd.read_csv(PILOT / "metrics.csv").query("split in @allowed_splits").copy()
pilot_fits = pd.read_csv(PILOT / "fits.csv")
pilot_predictions = pd.read_csv(PILOT / "predictions.csv").query("split in @allowed_splits").copy()
pilot_runtime = json.loads((PILOT / "runtime.json").read_text(encoding="utf-8-sig"))
run_protocol = json.loads(json.dumps(protocol))
run_protocol.update({
    "purpose": "Demo rerun or declared validation-only extension",
    "seeds": SEEDS, "evaluate_test": EVALUATE_TEST,
    "reference_architecture_declared": REFERENCE_ARCHITECTURE,
    "extension_note": EXTENSION_NOTE or ("New predeclared seeds; unchanged architecture." if SEEDS != PILOT_SEEDS else "Original architecture and seeds."),
    "source_pilot_protocol_sha256": hashlib.sha256((PILOT / "protocol.json").read_bytes()).hexdigest(),
    "evaluation_sequence": "Fit all candidates on train/validation; freeze artifacts; " + ("evaluate previously inspected test once." if EVALUATE_TEST else "validation only; no test predictions or test metrics."),
})
if not REFERENCE_ARCHITECTURE:
    run_protocol["original_transformer_reference_specification"] = run_protocol.pop("transformer")
    run_protocol["actual_architecture"] = "See saved model configuration JSON files and EXTENSION_NOTE."
run_protocol["interpretation"] = "Exploratory retrospective; " + ("previously inspected test." if EVALUATE_TEST else "training and validation only; test not evaluated.")
assert config["lookback_days"] == 30 and config["horizon_days"] == 7
assert config["input_magnitude"] == 4.5 and config["target_magnitude"] == 5
for relative, expected in protocol["input_hashes"].items():
    assert hashlib.sha256((ROOT / relative).read_bytes()).hexdigest() == expected, f"Changed input: {relative}"
pilot_hashes_before = {p.name:hashlib.sha256(p.read_bytes()).hexdigest() for p in PILOT.iterdir() if p.is_file()}
print("Course folder:", ROOT)
print("Frozen inputs verified; original pilot completed:", pilot_runtime["completed_utc"])

Course folder: H:\Otros ordenadores\Mi PC\1_Proyectos\2026 - Statsei 14 -\statsei14-deep-learning
Frozen inputs verified; original pilot completed: 2026-10-07T14:01:48.551600+00:00

2. The comparison fixed before fitting

Fixed choice Setting
Input Log count, maximum magnitude excess and event presence for 30 days
Region 28–34°S, 69–74°W; input catalogue M ≥ 4.5
Target At least one catalogue M ≥ 5 event in the following seven days
Origins Mondays; target weeks crossing partition boundaries excluded
Train / validation / test 2001–2016 / 2017–2020 / 2021–2025
Training Adam 0.001, binary cross-entropy, batch 64, maximum 20 epochs, no shuffling
Stopping Validation loss, patience 3, restore best weights
Seeds 7, 17 and 27 for both architectures

There are 834 training weeks with 197 positives. More parameters do not automatically provide useful information from this sample.

The table below is existing pilot evidence, separate from a new run. Variation across seeds describes optimization variability, not uncertainty about future generalization.

In [3]:
display(pilot_metrics.groupby(["model", "split"])[["brier", "log_loss"]].agg(["mean", "std"]).round(6))
display(pilot_fits[["model", "seed", "parameter_count", "epochs", "best_epoch", "fit_seconds"]].round(3))
print("Original pilot fit seconds:", round(pilot_runtime["total_fit_seconds"], 3))
print("Original pilot Colab verified:", pilot_runtime["colab_execution_verified"])
brier log_loss
mean std mean std
model split
LSTM catalogue test 0.165413 0.000185 0.512942 0.000515
validation 0.183065 0.001354 0.550741 0.002846
Logistic catalogue test 0.164229 NaN 0.510317 NaN
validation 0.186611 NaN 0.560023 NaN
Training frequency test 0.165370 NaN 0.513204 NaN
validation 0.187690 NaN 0.562855 NaN
Transformer compact test 0.166367 0.000645 0.515602 0.001897
validation 0.188771 0.001369 0.565122 0.002935
model seed parameter_count epochs best_epoch fit_seconds
0 LSTM catalogue 7 1297 20 20 2.054
1 Transformer compact 7 2305 6 3 1.398
2 Transformer compact 17 2305 9 6 1.587
3 LSTM catalogue 17 1297 20 20 2.005
4 LSTM catalogue 27 1297 20 20 1.997
5 Transformer compact 27 2305 5 2 1.338 
Original pilot fit seconds: 10.379
Original pilot Colab verified: False

3. Identical windows and training-only scaling

The core practical constructs the supplied windows from the frozen catalogue. This notebook checks their fingerprints and boundaries, then takes only the three catalogue channels. The original pilot also reconstructed features and labels independently from the event records.

X is weeks × 30 days × 3 features. Inputs are [origin − 30 days, origin); targets are [origin, origin + 7 days). Adjacent histories can overlap even though target weeks do not. Scaling uses training windows only. Event time does not certify the historical availability of revised records.

In [4]:
with np.load(ROOT / "results/reference-run/features_raw.npz", allow_pickle=False) as archive:
    assert archive["feature_names"].tolist() == FEATURES + ["log_etas_score"]
    X_raw = archive["X_raw"][:, :, :3].copy()
windows = pd.read_csv(ROOT / "results/reference-run/window_metadata.csv")
for column in ["origin", "input_start", "input_end", "target_end"]:
    windows[column] = pd.to_datetime(windows[column], utc=True)
y = windows.y.to_numpy(dtype="int32")
assert X_raw.shape == (len(windows), 30, 3) and np.isfinite(X_raw).all()
assert windows.origin.is_unique and windows.origin.is_monotonic_increasing
assert (windows.input_end == windows.origin).all()
assert ((windows.origin - windows.input_start).dt.days == 30).all()
assert ((windows.target_end - windows.origin).dt.days == 7).all()
assert (windows.origin.dt.dayofweek == 0).all()
assert (windows.target_end.iloc[:-1].to_numpy() <= windows.origin.iloc[1:].to_numpy()).all()
bounds = {name:pd.Timestamp(config[name], tz="UTC") for name in ["train_start", "validation_start", "test_start", "catalog_end"]}
masks = {
    "train": (windows.origin.ge(bounds["train_start"]) & windows.origin.lt(bounds["validation_start"]) & windows.target_end.le(bounds["validation_start"])).to_numpy(),
    "validation": (windows.origin.ge(bounds["validation_start"]) & windows.origin.lt(bounds["test_start"]) & windows.target_end.le(bounds["test_start"])).to_numpy(),
    "test": (windows.origin.ge(bounds["test_start"]) & windows.origin.lt(bounds["catalog_end"]) & windows.target_end.le(bounds["catalog_end"])).to_numpy(),
}
expected_sizes = {"train":(834,197), "validation":(208,52), "test":(260,54)}
assert (sum(mask.astype(int) for mask in masks.values()) <= 1).all()
for split, mask in masks.items():
    assert np.array_equal(mask, windows.split.eq(split).to_numpy())
    if split != "test" or EVALUATE_TEST:
        assert (int(mask.sum()), int(y[mask].sum())) == expected_sizes[split]
train, validation, test = masks["train"], masks["validation"], masks["test"]
scaler = StandardScaler().fit(X_raw[train].reshape(-1, 3))
reference_scaler = json.loads((ROOT / "results/reference-run/normalization.json").read_text())
np.testing.assert_allclose(scaler.mean_, reference_scaler["mean"][:3], rtol=1e-10, atol=1e-10)
np.testing.assert_allclose(scaler.scale_, reference_scaler["scale"][:3], rtol=1e-10, atol=1e-10)
X = scaler.transform(X_raw.reshape(-1, 3)).reshape(X_raw.shape).astype("float32")
print("Input shape:", X.shape, "training channel means:", scaler.mean_.round(4))
display(windows.loc[windows.split.isin(["train"] + allowed_splits)].groupby("split").agg(weeks=("y", "size"), positives=("y", "sum")))

Input shape: (1305, 30, 3) training channel means: [0.1256 0.055  0.1504]
weeks positives
split
test 260 54
train 834 197
validation 208 52

4. Two ways to combine the same past

The LSTM carries a learned state through the days. The Transformer lets each day combine information from other days in the same window.

The compact model uses projection to 16 features, fixed day positions, two attention heads, residual connections and normalization, a feed-forward block, pooling and a sigmoid output. A causal attention mask is unnecessary because all input days precede the forecast origin. Position encoding preserves day order.

This position layer has no learned parameters. It is registered for serialization; run its definition before loading a saved model in another process. Attention weights alone are not evidence of a causal seismic mechanism.

In [5]:
@keras.utils.register_keras_serializable(package="StatseiDemo")
class SinusoidalPosition(keras.layers.Layer):
    def __init__(self, length=30, width=16, **kwargs):
        super().__init__(**kwargs)
        self.length, self.width = length, width
        positions = np.arange(length, dtype=np.float32)[:, None]
        frequencies = np.exp(np.arange(0, width, 2, dtype=np.float32) * (-np.log(10000.0) / width))
        encoding = np.zeros((length, width), dtype=np.float32)
        encoding[:, 0::2] = np.sin(positions * frequencies)
        encoding[:, 1::2] = np.cos(positions * frequencies)
        self.encoding = encoding[None, :, :]
    def call(self, inputs):
        return inputs + tf.cast(tf.convert_to_tensor(self.encoding), inputs.dtype)
    def get_config(self):
        return {**super().get_config(), "length":self.length, "width":self.width}

5. Read the fitting loop, or run it after the session

Extensions: predeclare three new SEEDS with RUN_TRAINING=True and EVALUATE_TEST=False. For a manual architecture change, also set REFERENCE_ARCHITECTURE=False and describe it in EXTENSION_NOTE before editing the model code. The saved model configuration records the actual architecture. This is not an architecture search.

Input, layers, compile, fit and predict remain visible. The loops pair seeds and change architecture. Order alternates to reduce a systematic timing advantage from running second. All six candidates are retained; none is selected using test results.

The constant reference uses training prevalence. Logistic regression uses nine mean/max/last summaries of the same channels. Normalization and fitting remain restricted to training.

New outputs go to a timestamped folder under results/transformer-demo-run. The original pilot and core practical remain unchanged. If training is unavailable, leave RUN_TRAINING = False.

In [6]:
models, histories, validation_probabilities = {}, {}, {}
fit_rows = []
new_run_output = None
freeze = None
if RUN_TRAINING:
    run_name = datetime.datetime.now(datetime.timezone.utc).strftime("%Y%m%dT%H%M%SZ") + "-" + uuid.uuid4().hex[:6]
    new_run_output = ROOT / "results/transformer-demo-run" / run_name
    (new_run_output / "models").mkdir(parents=True, exist_ok=False)
    (new_run_output / "histories").mkdir()
    base_rate = float(y[train].mean())
    summaries = np.concatenate([X_raw.mean(axis=1), X_raw.max(axis=1), X_raw[:, -1, :]], axis=1)
    summary_scaler = StandardScaler().fit(summaries[train])
    Z = summary_scaler.transform(summaries)
    logistic = LogisticRegression(C=1.0, max_iter=1000, random_state=20261013)
    logistic.fit(Z[train], y[train])
    validation_probabilities[("Training frequency", -1)] = np.full(int(validation.sum()), base_rate)
    validation_probabilities[("Logistic catalogue", -1)] = logistic.predict_proba(Z[validation])[:, 1]
    for seed_index, seed in enumerate(SEEDS):
        order = ["LSTM catalogue", "Transformer compact"] if seed_index % 2 == 0 else ["Transformer compact", "LSTM catalogue"]
        for model_name in order:
            keras.backend.clear_session()
            keras.utils.set_random_seed(seed)
            inputs = keras.Input(shape=(30, 3), name="past_30_days")
            if model_name == "LSTM catalogue":
                sequence_summary = keras.layers.LSTM(16, name="lstm")(inputs)
            else:
                x = keras.layers.Dense(16, name="daily_projection")(inputs)
                x = SinusoidalPosition(name="day_position")(x)
                attention = keras.layers.MultiHeadAttention(
                    num_heads=2, key_dim=8, dropout=0.0, name="self_attention"
                )(x, x)
                x = keras.layers.LayerNormalization(epsilon=1e-6, name="attention_norm")(
                    keras.layers.Add()([x, attention])
                )
                feed_forward = keras.layers.Dense(32, activation="relu", name="ff_expand")(x)
                feed_forward = keras.layers.Dense(16, name="ff_project")(feed_forward)
                x = keras.layers.LayerNormalization(epsilon=1e-6, name="ff_norm")(
                    keras.layers.Add()([x, feed_forward])
                )
                sequence_summary = keras.layers.GlobalAveragePooling1D(name="past_window_pooling")(x)
            output = keras.layers.Dense(1, activation="sigmoid", name="next_week_probability")(sequence_summary)
            model = keras.Model(inputs, output, name=model_name.lower().replace(" ", "_"))
            model.compile(optimizer=keras.optimizers.Adam(learning_rate=0.001), loss="binary_crossentropy")
            early_stopping = keras.callbacks.EarlyStopping(monitor="val_loss", patience=3, restore_best_weights=True)
            fit_started = time.perf_counter()
            history = model.fit(
                X[train], y[train], validation_data=(X[validation], y[validation]),
                epochs=20, batch_size=64, callbacks=[early_stopping], shuffle=False, verbose=0,
            )
            fit_seconds = time.perf_counter() - fit_started
            assert np.isfinite(np.asarray(list(history.history.values()))).all()
            probability = model.predict(X[validation], verbose=0).ravel()
            assert np.isfinite(probability).all() and ((0 <= probability) & (probability <= 1)).all()
            slug = model_name.lower().replace(" ", "_") + f"_seed_{seed}"
            model_path = new_run_output / "models" / f"{slug}.keras"
            model.save(model_path)
            model_path.with_suffix(".json").write_text(model.to_json(), encoding="utf-8")
            reloaded = keras.models.load_model(model_path)
            np.testing.assert_allclose(probability, reloaded.predict(X[validation], verbose=0).ravel(), rtol=1e-6, atol=1e-7)
            models[(model_name, seed)] = model_path
            histories[(model_name, seed)] = pd.DataFrame(history.history)
            histories[(model_name, seed)].rename_axis("epoch_index").to_csv(new_run_output / "histories" / f"{slug}.csv")
            validation_probabilities[(model_name, seed)] = probability
            fit_rows.append({"model":model_name, "seed":seed, "parameter_count":model.count_params(),
                "epochs":len(history.history["loss"]), "best_epoch":int(np.argmin(history.history["val_loss"])) + 1,
                "fit_seconds":fit_seconds, "model_file":str(model_path.relative_to(new_run_output))})
            print(model_name, "seed", seed, "epochs", fit_rows[-1]["epochs"],
                  "validation Brier", round(float(brier_score_loss(y[validation], probability)), 6))
    fits = pd.DataFrame(fit_rows)
    print("Six fits complete. New test predictions have not yet been computed.")
else:
    fits = pilot_fits.copy()
    print("Prepared mode: no models fitted and no new test predictions computed.")
display(fits[["model", "seed", "parameter_count", "epochs", "best_epoch", "fit_seconds"]].round(3))

WARNING:tensorflow:From C:\Users\Francisco Plaza\AppData\Local\Temp\statsei14-transformer-pilot-20261007\Lib\site-packages\keras\src\backend\common\global_state.py:82: The name tf.reset_default_graph is deprecated. Please use tf.compat.v1.reset_default_graph instead.


LSTM catalogue seed 7 epochs 20 validation Brier 0.181521

WARNING:tensorflow:5 out of the last 15 calls to <function TensorFlowTrainer.make_predict_function.<locals>.one_step_on_data_distributed at 0x0000024205A2E200> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.

WARNING:tensorflow:5 out of the last 15 calls to <function TensorFlowTrainer.make_predict_function.<locals>.one_step_on_data_distributed at 0x0000024209583F60> triggered tf.function retracing. Tracing is expensive and the excessive number of tracings could be due to (1) creating @tf.function repeatedly in a loop, (2) passing tensors with different shapes, (3) passing Python objects instead of tensors. For (1), please define your @tf.function outside of the loop. For (2), @tf.function has reduce_retracing=True option that can avoid unnecessary retracing. For (3), please refer to https://www.tensorflow.org/guide/function#controlling_retracing and https://www.tensorflow.org/api_docs/python/tf/function for  more details.

Transformer compact seed 7 epochs 6 validation Brier 0.190349

Transformer compact seed 17 epochs 9 validation Brier 0.188069

LSTM catalogue seed 17 epochs 20 validation Brier 0.183624

LSTM catalogue seed 27 epochs 20 validation Brier 0.18405

Transformer compact seed 27 epochs 5 validation Brier 0.187895
Six fits complete. New test predictions have not yet been computed.
model seed parameter_count epochs best_epoch fit_seconds
0 LSTM catalogue 7 1297 20 20 2.132
1 Transformer compact 7 2305 6 3 1.532
2 Transformer compact 17 2305 9 6 1.704
3 LSTM catalogue 17 1297 20 20 2.068
4 LSTM catalogue 27 1297 20 20 2.073
5 Transformer compact 27 2305 5 2 1.443

6. Freeze, then evaluate the previously inspected test

A new run records all six model fingerprints, validation scores and the actual protocol. With EVALUATE_TEST=False, evaluation ends on validation: no test predictions or test metrics are computed or displayed. With True, the unchanged original comparison then computes test predictions. This keeps a changed candidate from being presented as the original tested candidate. It does not make the previously inspected test blind again.

Retain all seeds. A mean of per-seed scores is not the performance of an ensemble of their probabilities.

In [7]:
if RUN_TRAINING:
    validation_rows = []
    for (name, seed), probability in validation_probabilities.items():
        validation_rows.append({"model":name, "seed":seed, "split":"validation",
            "brier":float(brier_score_loss(y[validation], probability)),
            "log_loss":float(log_loss(y[validation], probability, labels=[0,1])),
            "n":int(validation.sum()), "positives":int(y[validation].sum()), "positive_rate":float(y[validation].mean())})
    pd.DataFrame(validation_rows).to_csv(new_run_output / "validation_metrics.csv", index=False)
    freeze = {
        "frozen_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(),
        "source_pilot_protocol_sha256":run_protocol["source_pilot_protocol_sha256"],
        "actual_run_protocol":run_protocol,
        "seeds":SEEDS, "features":FEATURES, "evaluate_test":EVALUATE_TEST,
        "models":{str(p.relative_to(new_run_output)):hashlib.sha256(p.read_bytes()).hexdigest() for p in models.values()},
        "validation_metrics_sha256":hashlib.sha256((new_run_output / "validation_metrics.csv").read_bytes()).hexdigest(),
        "selection":"All six predeclared fits retained; no test-based selection or tuning.",
        "test_scope":"Previously inspected test; exploratory retrospective." if EVALUATE_TEST else "Not evaluated: validation-only extension.",
    }
    (new_run_output / "fit_freeze.json").write_text(json.dumps(freeze, indent=2) + "\n", encoding="utf-8")
    evaluation_probabilities = {"validation": validation_probabilities}
    if EVALUATE_TEST:
        test_probabilities = {
            ("Training frequency", -1):np.full(int(test.sum()), base_rate),
            ("Logistic catalogue", -1):logistic.predict_proba(Z[test])[:, 1],
        }
        for key, model_path in models.items():
            frozen_model = keras.models.load_model(model_path)
            test_probabilities[key] = frozen_model.predict(X[test], verbose=0).ravel()
        evaluation_probabilities["test"] = test_probabilities
    metric_rows, prediction_frames = [], []
    for split, probabilities in evaluation_probabilities.items():
        mask = masks[split]
        for (name, seed), probability in probabilities.items():
            assert len(probability) == int(mask.sum()) and np.isfinite(probability).all()
            assert ((probability >= 0) & (probability <= 1)).all()
            metric_rows.append({"model":name, "seed":seed, "split":split,
                "brier":float(brier_score_loss(y[mask], probability)),
                "log_loss":float(log_loss(y[mask], probability, labels=[0,1])),
                "n":int(mask.sum()), "positives":int(y[mask].sum()), "positive_rate":float(y[mask].mean())})
            prediction_frames.append(pd.DataFrame({"model":name, "seed":seed, "split":split,
                "origin":windows.loc[mask,"origin"].to_numpy(), "y_true":y[mask], "probability":probability}))
    metrics = pd.DataFrame(metric_rows)
    predictions = pd.concat(prediction_frames, ignore_index=True)
    evidence_label = "New fixed-protocol demo run" if EVALUATE_TEST else "New validation-only run; see declared protocol"
else:
    metrics, predictions = pilot_metrics.copy(), pilot_predictions.copy()
    evidence_label = "Original prepared pilot" + ("" if EVALUATE_TEST else "; validation only")
assert not predictions.duplicated(["model", "seed", "split", "origin"]).any()
assert predictions.groupby(["split", "origin"]).y_true.nunique().eq(1).all()
print(evidence_label)
display(metrics.round(6))

New fixed-protocol demo run
model seed split brier log_loss n positives positive_rate
0 Training frequency -1 validation 0.187690 0.562855 208 52 0.250000
1 Logistic catalogue -1 validation 0.186611 0.560023 208 52 0.250000
2 LSTM catalogue 7 validation 0.181521 0.547512 208 52 0.250000
3 Transformer compact 7 validation 0.190349 0.568482 208 52 0.250000
4 Transformer compact 17 validation 0.188069 0.563823 208 52 0.250000
5 LSTM catalogue 17 validation 0.183624 0.551824 208 52 0.250000
6 LSTM catalogue 27 validation 0.184050 0.552887 208 52 0.250000
7 Transformer compact 27 validation 0.187895 0.563060 208 52 0.250000
8 Training frequency -1 test 0.165370 0.513204 260 54 0.207692
9 Logistic catalogue -1 test 0.164229 0.510317 260 54 0.207692
10 LSTM catalogue 7 test 0.165309 0.512480 260 54 0.207692
11 Transformer compact 7 test 0.165821 0.513814 260 54 0.207692
12 Transformer compact 17 test 0.167079 0.517592 260 54 0.207692
13 LSTM catalogue 17 test 0.165626 0.513497 260 54 0.207692
14 LSTM catalogue 27 test 0.165304 0.512850 260 54 0.207692
15 Transformer compact 27 test 0.166200 0.515401 260 54 0.207692

7. Read errors and probability forecasts

Brier score is squared probability error; log loss penalizes confident incorrect probabilities. Smaller is better. Pair models by seed and evaluation period instead of selecting each model's best seed.

Calibration uses the same ten fixed bins for every model and seed. Empty bins stay in the table with zero counts and missing means; only occupied bins are plotted. The table shows weeks per bin: a handful of high-probability weeks cannot establish precise calibration. The distribution panel reveals whether predictions span a narrow range.

In [8]:
display(metrics.groupby(["model", "split"])[["brier", "log_loss"]].agg(["mean", "std"]).round(6))
paired = metrics.loc[metrics.seed.isin(SEEDS)].pivot(index=["split", "seed"], columns="model", values=["brier", "log_loss"])
differences = pd.DataFrame({
    "Transformer minus LSTM: Brier":paired["brier"]["Transformer compact"] - paired["brier"]["LSTM catalogue"],
    "Transformer minus LSTM: log loss":paired["log_loss"]["Transformer compact"] - paired["log_loss"]["LSTM catalogue"],
})
display(differences.round(6))
print("Positive values mean higher Transformer error in that pair.")
print("Seed variation describes these fits, not generalization uncertainty.")
brier log_loss
mean std mean std
model split
LSTM catalogue test 0.165413 0.000185 0.512942 0.000515
validation 0.183065 0.001354 0.550741 0.002846
Logistic catalogue test 0.164229 NaN 0.510317 NaN
validation 0.186611 NaN 0.560023 NaN
Training frequency test 0.165370 NaN 0.513204 NaN
validation 0.187690 NaN 0.562855 NaN
Transformer compact test 0.166367 0.000645 0.515602 0.001897
validation 0.188771 0.001369 0.565122 0.002935
Transformer minus LSTM: Brier Transformer minus LSTM: log loss
split seed
test 7 0.000511 0.001334
17 0.001453 0.004095
27 0.000896 0.002551
validation 7 0.008828 0.020970
17 0.004444 0.011999
27 0.003844 0.010173 
Positive values mean higher Transformer error in that pair.
Seed variation describes these fits, not generalization uncertainty.
In [9]:
BIN_EDGES = np.asarray(protocol["calibration_bin_edges"], dtype=float)
calibration_rows = []
for (name, seed, split), group in predictions.loc[predictions.seed.isin(SEEDS)].groupby(["model", "seed", "split"]):
    grouped = group.assign(bin=pd.cut(group.probability, bins=BIN_EDGES, include_lowest=True))
    table = grouped.groupby("bin", observed=False).agg(
        n=("y_true", "size"), mean_probability=("probability", "mean"), observed_fraction=("y_true", "mean")
    ).reset_index()
    table["bin"] = table["bin"].astype(str)
    table["model"], table["seed"], table["split"] = name, seed, split
    calibration_rows.append(table)
calibration = pd.concat(calibration_rows, ignore_index=True)
colors = {"LSTM catalogue":"#087f78", "Transformer compact":"#d27b25"}
plot_split = "test" if EVALUATE_TEST else "validation"
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
axes[0].plot([0,1], [0,1], "--", color="gray", linewidth=1)
for (name, seed), group in calibration.loc[calibration.split.eq(plot_split)].groupby(["model", "seed"]):
    occupied_bins = group.loc[group.n.gt(0)]
    axes[0].plot(occupied_bins.mean_probability, occupied_bins.observed_fraction, "o-", alpha=0.65,
                 color=colors[name], label=f"{name}; seed {seed}")
    probability = predictions.loc[predictions.model.eq(name) & predictions.seed.eq(seed) & predictions.split.eq(plot_split), "probability"]
    axes[1].hist(probability, bins=BIN_EDGES, histtype="step", linewidth=1.5, alpha=0.7,
                 color=colors[name], label=f"{name}; seed {seed}")
axes[0].set(xlabel="Mean predicted probability", ylabel="Observed positive fraction",
            xlim=(0,1), ylim=(0,1), title=f"{plot_split.title()}: fixed bins")
axes[1].set(xlabel="Predicted probability", ylabel="Weeks per run", title="Probability distribution", xlim=(0,1))
axes[1].legend(fontsize=7)
fig.suptitle(evidence_label)
plt.tight_layout()
plt.show()
display(calibration.loc[calibration.split.eq(plot_split), ["model", "seed", "bin", "n", "mean_probability", "observed_fraction"]].round(4))
No description has been provided for this image
model seed bin n mean_probability observed_fraction
0 LSTM catalogue 7 (-0.001, 0.1] 0 NaN NaN
1 LSTM catalogue 7 (0.1, 0.2] 37 0.1718 0.2162
2 LSTM catalogue 7 (0.2, 0.3] 202 0.2295 0.2079
3 LSTM catalogue 7 (0.3, 0.4] 18 0.3344 0.1667
4 LSTM catalogue 7 (0.4, 0.5] 2 0.4180 0.0000
5 LSTM catalogue 7 (0.5, 0.6] 1 0.5245 1.0000
6 LSTM catalogue 7 (0.6, 0.7] 0 NaN NaN
7 LSTM catalogue 7 (0.7, 0.8] 0 NaN NaN
8 LSTM catalogue 7 (0.8, 0.9] 0 NaN NaN
9 LSTM catalogue 7 (0.9, 1.0] 0 NaN NaN
20 LSTM catalogue 17 (-0.001, 0.1] 0 NaN NaN
21 LSTM catalogue 17 (0.1, 0.2] 27 0.1740 0.2593
22 LSTM catalogue 17 (0.2, 0.3] 217 0.2311 0.2074
23 LSTM catalogue 17 (0.3, 0.4] 16 0.3291 0.1250
24 LSTM catalogue 17 (0.4, 0.5] 0 NaN NaN
25 LSTM catalogue 17 (0.5, 0.6] 0 NaN NaN
26 LSTM catalogue 17 (0.6, 0.7] 0 NaN NaN
27 LSTM catalogue 17 (0.7, 0.8] 0 NaN NaN
28 LSTM catalogue 17 (0.8, 0.9] 0 NaN NaN
29 LSTM catalogue 17 (0.9, 1.0] 0 NaN NaN
40 LSTM catalogue 27 (-0.001, 0.1] 0 NaN NaN
41 LSTM catalogue 27 (0.1, 0.2] 24 0.1828 0.1667
42 LSTM catalogue 27 (0.2, 0.3] 230 0.2317 0.2130
43 LSTM catalogue 27 (0.3, 0.4] 5 0.3440 0.0000
44 LSTM catalogue 27 (0.4, 0.5] 1 0.4061 1.0000
45 LSTM catalogue 27 (0.5, 0.6] 0 NaN NaN
46 LSTM catalogue 27 (0.6, 0.7] 0 NaN NaN
47 LSTM catalogue 27 (0.7, 0.8] 0 NaN NaN
48 LSTM catalogue 27 (0.8, 0.9] 0 NaN NaN
49 LSTM catalogue 27 (0.9, 1.0] 0 NaN NaN
60 Transformer compact 7 (-0.001, 0.1] 0 NaN NaN
61 Transformer compact 7 (0.1, 0.2] 54 0.1895 0.1667
62 Transformer compact 7 (0.2, 0.3] 194 0.2325 0.2216
63 Transformer compact 7 (0.3, 0.4] 12 0.3349 0.1667
64 Transformer compact 7 (0.4, 0.5] 0 NaN NaN
65 Transformer compact 7 (0.5, 0.6] 0 NaN NaN
66 Transformer compact 7 (0.6, 0.7] 0 NaN NaN
67 Transformer compact 7 (0.7, 0.8] 0 NaN NaN
68 Transformer compact 7 (0.8, 0.9] 0 NaN NaN
69 Transformer compact 7 (0.9, 1.0] 0 NaN NaN
80 Transformer compact 17 (-0.001, 0.1] 0 NaN NaN
81 Transformer compact 17 (0.1, 0.2] 2 0.1975 0.5000
82 Transformer compact 17 (0.2, 0.3] 249 0.2446 0.2048
83 Transformer compact 17 (0.3, 0.4] 9 0.3287 0.2222
84 Transformer compact 17 (0.4, 0.5] 0 NaN NaN
85 Transformer compact 17 (0.5, 0.6] 0 NaN NaN
86 Transformer compact 17 (0.6, 0.7] 0 NaN NaN
87 Transformer compact 17 (0.7, 0.8] 0 NaN NaN
88 Transformer compact 17 (0.8, 0.9] 0 NaN NaN
89 Transformer compact 17 (0.9, 1.0] 0 NaN NaN
100 Transformer compact 27 (-0.001, 0.1] 0 NaN NaN
101 Transformer compact 27 (0.1, 0.2] 0 NaN NaN
102 Transformer compact 27 (0.2, 0.3] 253 0.2570 0.2095
103 Transformer compact 27 (0.3, 0.4] 7 0.3127 0.1429
104 Transformer compact 27 (0.4, 0.5] 0 NaN NaN
105 Transformer compact 27 (0.5, 0.6] 0 NaN NaN
106 Transformer compact 27 (0.6, 0.7] 0 NaN NaN
107 Transformer compact 27 (0.7, 0.8] 0 NaN NaN
108 Transformer compact 27 (0.8, 0.9] 0 NaN NaN
109 Transformer compact 27 (0.9, 1.0] 0 NaN NaN

8. What does stopping tell us?

In the original pilot, the Transformer stopped after 5–9 epochs; restored best epochs were 2, 3 and 6. The LSTM's best epoch was 20 in all three runs, reaching the fixed budget. This describes the chosen configuration and stopping rule. It neither establishes faster Transformer convergence nor diagnoses a unique cause of error.

A shorter fit also used fewer epochs. Distinguish runtime, optimization behavior and predictive evidence. Do not extend training based on these test results while calling the same test untouched.

In [10]:
if RUN_TRAINING:
    fig, axes = plt.subplots(1, 2, figsize=(11, 3.5))
    for (name, seed), history in histories.items():
        epochs = np.arange(1, len(history) + 1)
        for axis, metric in zip(axes, ["loss", "val_loss"]):
            axis.plot(epochs, history[metric], color=colors[name], alpha=0.65, label=f"{name}; {seed}")
    axes[0].set(title="Training loss", xlabel="Epoch", ylabel="Binary cross-entropy")
    axes[1].set(title="Validation loss", xlabel="Epoch", ylabel="Binary cross-entropy")
    axes[1].legend(fontsize=7)
    plt.tight_layout()
    plt.show()
else:
    print("Prepared mode: stopping evidence comes from the original pilot summaries.")
display(fits[["model", "seed", "epochs", "best_epoch", "fit_seconds"]].round(3))
No description has been provided for this image
model seed epochs best_epoch fit_seconds
0 LSTM catalogue 7 20 20 2.132
1 Transformer compact 7 6 3 1.532
2 Transformer compact 17 9 6 1.704
3 LSTM catalogue 17 20 20 2.068
4 LSTM catalogue 27 20 20 2.073
5 Transformer compact 27 5 2 1.443

9. Keep new evidence separate

A new run exports metrics, predictions, calibration, normalization, models, histories, protocol and runtime information into its own folder. The frozen pilot stays unchanged. Differences across environments are possible; we measure them rather than replacing the reference.

In [11]:
if RUN_TRAINING:
    metrics.to_csv(new_run_output / "metrics.csv", index=False)
    predictions.to_csv(new_run_output / "predictions.csv", index=False)
    fits.to_csv(new_run_output / "fits.csv", index=False)
    calibration.to_csv(new_run_output / "calibration.csv", index=False)
    differences.to_csv(new_run_output / "paired_differences.csv")
    windows.loc[windows.split.isin(["train"] + allowed_splits)].to_csv(new_run_output / "window_metadata.csv", index=False)
    normalization = {"fitted_split":"train", "features":FEATURES, "mean":scaler.mean_.tolist(), "scale":scaler.scale_.tolist()}
    (new_run_output / "normalization.json").write_text(json.dumps(normalization, indent=2) + "\n", encoding="utf-8")
    (new_run_output / "protocol.json").write_text(json.dumps(run_protocol, indent=2) + "\n", encoding="utf-8")
    merged = metrics.merge(pilot_metrics, on=["model", "seed", "split"], suffixes=("_new", "_pilot"), validate="one_to_one")
    maximum_difference = float(np.max(np.abs(merged[["brier_new", "log_loss_new"]].to_numpy() - merged[["brier_pilot", "log_loss_pilot"]].to_numpy()))) if len(merged) else None
    runtime = {
        "completed_utc":datetime.datetime.now(datetime.timezone.utc).isoformat(), "completed":True,
        "python":platform.python_version(), "tensorflow":tf.__version__, "keras":keras.__version__,
        "numpy":np.__version__, "pandas":pd.__version__, "scikit_learn":sklearn.__version__,
        "platform":platform.platform(), "device":"CPU", "configured_threads":2,
        "fit_seconds":float(fits.fit_seconds.sum()),
        "notebook_seconds_since_imports":time.perf_counter() - NOTEBOOK_STARTED,
        "max_abs_metric_difference_from_pilot":maximum_difference,
        "scope":run_protocol["interpretation"], "seeds":SEEDS, "evaluate_test":EVALUATE_TEST,
        "reference_architecture_declared":REFERENCE_ARCHITECTURE,
        "extension_note":run_protocol["extension_note"], "pilot_comparison_matched_rows":len(merged),
        "colab_execution_verified":False, "runtime_note":"This machine only; no separate Colab execution certified.",
    }
    (new_run_output / "runtime.json").write_text(json.dumps(runtime, indent=2) + "\n", encoding="utf-8")
    print("New output folder:", new_run_output)
    print("Total fit seconds:", round(runtime["fit_seconds"], 3))
    print("Matched pilot rows:", len(merged), "; maximum metric difference:", maximum_difference)
    if SEEDS != PILOT_SEEDS:
        print("New seeds: only shared baseline rows match the original pilot; no missing-seed scores inferred.")
else:
    runtime = pilot_runtime
    print("Prepared pilot mode wrote no run outputs.")
pilot_hashes_after = {p.name:hashlib.sha256(p.read_bytes()).hexdigest() for p in PILOT.iterdir() if p.is_file()}
assert pilot_hashes_before == pilot_hashes_after, "Prepared pilot files changed."
print("Original prepared pilot files are unchanged.")

New output folder: H:\Otros ordenadores\Mi PC\1_Proyectos\2026 - Statsei 14 -\statsei14-deep-learning\results\transformer-demo-run\20261007T144119Z-88594f
Total fit seconds: 10.951
Matched pilot rows: 16 ; maximum metric difference: 8.326672684688674e-17
Original prepared pilot files are unchanged.

Discuss the result, not the model name

  1. Did architectural flexibility improve probability errors against the LSTM and simple references?
  2. What information do daily regional summaries preserve or discard? Could another question require another representation?
  3. What do seed variation, stopping and calibration bin counts leave uncertain?
  4. What additional data and untouched evaluation would support a stronger claim?

The prepared pilot's compact Transformer did not improve probability errors. This is a useful result: flexibility must match the information and question available. It concerns one small architecture, dataset and protocol, not a ranking of model families. Inputs are catalogue summaries, not raw waveforms or fine spatial grids.

The core practical's ETAS-inspired channel explores a different decision: changing the representation. A structural ETAS–neural model would need a fitted intensity and an appropriate forecast/evaluation formulation; this notebook does not implement it.

Continue with the homework guide. Sources: Keras attention layer, Keras time-series Transformer example, and Farfán, Nicolis & Peralta (2025), a different spatial-grid task.