Change_Model.py

"""Factored two-head change model: one shared encoder, a land-cover head and a change head.

Why not just difference two maps. Classifying 2019 and 2025 independently and subtracting
compounds both maps' errors, and change is rare, so the false positives leaking out of the
large stable area can rival the real change. Measured on the Oregon run, only about 8% of the
2.38 M ha of candidate change reached 90% joint confidence. The change head here reads the
before and after embeddings together, so it never forms two independent maps and cannot
compound their errors: it learns the joint before-to-after signature of real conversion
against land that merely looked different across two annual snapshots.

Why the two heads train on different data. GLanCE has gold land-cover labels but is North
America wide and thin on Pacific Northwest fire and harvest. The Oregon reference points have
gold "what changed" labels from COLD, the NDVI trend test, LandTrendr and visual review, but
their from/to classes are the old classifier's own predictions. Using those as state targets
would feed that classifier's error straight back into the model meant to correct it. So the
loss is masked per example: each example contributes only to the head it has a trustworthy
label for, and the two tasks meet only through the shared encoder. The encoder is what
couples them, not the examples, which is why no example ever needs a label it does not have.

The four variants the notebook compares are all this one function with different head weights:
    state only   -> change_weight=0   (a plain land-cover classifier, the state baseline)
    change only  -> state_weight=0    (no GLanCE help, isolates what sharing buys)
    frozen       -> pretrained=<state-only model>, freeze_trunk=True, state_weight=0
    joint        -> both weights nonzero (the factored model itself)
"""

import time

import numpy as np
import pandas as pd
import torch
from sklearn.model_selection import GroupShuffleSplit, ShuffleSplit, StratifiedShuffleSplit
from sklearn.preprocessing import LabelEncoder, StandardScaler
from sklearn.utils.class_weight import compute_class_weight
from torch import nn
from torch.utils.data import DataLoader, TensorDataset

# Land-cover classes ordered from least to most vegetated. Used only by the map-differencing
# baseline, to turn a pair of class names into a greenness direction comparable to the change
# head's increase/decrease labels.
GREENNESS_ORDER = ['Water', 'Developed', 'Barren/Sparse', 'Herbaceous', 'Shrubs', 'Trees']


class _SharedTrunk(nn.Module):
    """Per-year encoder applied Siamese style, with the same weights, to both years.

    Sharing the weights is the point: a 2019 embedding and a 2025 embedding are the same kind
    of object, so encoding them with two different networks would let the model learn a
    year-specific shortcut instead of a representation in which real transitions are visible.
    """

    def __init__(self, input_dim, hidden_dims, latent_dim, dropout):
        super().__init__()
        layers = []
        prev = input_dim
        for h in list(hidden_dims) + [latent_dim]:
            layers += [nn.Linear(prev, h), nn.BatchNorm1d(h), nn.ReLU(), nn.Dropout(dropout)]
            prev = h
        self.net = nn.Sequential(*layers)

    def forward(self, x):
        return self.net(x)


class _FactoredChangeNet(nn.Module):
    """Shared trunk plus a one-year state head and a two-year change head.

    The change head sees [h_before, h_after, h_before - h_after]. The explicit difference is
    redundant in principle, since a linear layer over the concatenation could form it, but
    handing it over directly makes the change signal available in the first layer instead of
    something the model has to spend capacity discovering.
    """

    def __init__(self, input_dim, hidden_dims, latent_dim, change_hidden,
                 n_state_classes, n_change_classes, dropout):
        super().__init__()
        self.trunk = _SharedTrunk(input_dim, hidden_dims, latent_dim, dropout)
        self.state_head = nn.Linear(latent_dim, n_state_classes)

        layers = []
        prev = 3 * latent_dim
        for h in change_hidden:
            layers += [nn.Linear(prev, h), nn.BatchNorm1d(h), nn.ReLU(), nn.Dropout(dropout)]
            prev = h
        layers.append(nn.Linear(prev, n_change_classes))
        self.change_head = nn.Sequential(*layers)

    def forward_state(self, x):
        return self.state_head(self.trunk(x))

    def forward_change(self, x_before, x_after):
        h_before, h_after = self.trunk(x_before), self.trunk(x_after)
        return self.change_head(torch.cat([h_before, h_after, h_before - h_after], dim=1))


class _HeadView:
    """One head exposed as a standalone sklearn-style classifier.

    evaluate_model and compare_models call model.predict(x_test) and expect original labels
    back, so each head needs its own object with that interface. A view holds a reference to
    the parent classifier rather than a copy of the network, so both views and the parent
    pickle together as a single object graph.
    """

    def __init__(self, parent, kind):
        self.parent = parent
        self.kind = kind

    @property
    def classes_(self):
        return (self.parent.state_classes_ if self.kind == 'state'
                else self.parent.change_classes_)

    def predict(self, X):
        return (self.parent.predict_state(X) if self.kind == 'state'
                else self.parent.predict_change(X))

    def predict_proba(self, X):
        return (self.parent.predict_proba_state(X) if self.kind == 'state'
                else self.parent.predict_proba_change(X))


class _FactoredChangeClassifier:
    """Fitted factored model, with `.state` and `.change` views for the two heads.

    Inference runs on CPU so joblib.dump/load stays device agnostic, matching MLP.py.
    """

    def __init__(self, model, state_encoder, change_encoder, scaler, input_dim):
        self.model = model
        self.state_encoder = state_encoder
        self.change_encoder = change_encoder
        self.scaler = scaler
        self.input_dim = input_dim
        self.state_classes_ = state_encoder.classes_
        self.change_classes_ = change_encoder.classes_
        self.state = _HeadView(self, 'state')
        self.change = _HeadView(self, 'change')

    def _scale(self, X):
        return self.scaler.transform(np.asarray(X, dtype=np.float32)).astype(np.float32)

    def _split_pair(self, X):
        """Split an (n, 2*input_dim) paired matrix into its before and after halves."""
        arr = np.asarray(X, dtype=np.float32)
        assert arr.shape[1] == 2 * self.input_dim, (
            f'expected {2 * self.input_dim} paired columns (before then after), '
            f'got {arr.shape[1]}')
        return arr[:, :self.input_dim], arr[:, self.input_dim:]

    def state_logits(self, X):
        self.model.eval()
        with torch.no_grad():
            return self.model.forward_state(torch.from_numpy(self._scale(X))).numpy()

    def change_logits(self, X):
        before, after = self._split_pair(X)
        self.model.eval()
        with torch.no_grad():
            return self.model.forward_change(
                torch.from_numpy(self._scale(before)),
                torch.from_numpy(self._scale(after))).numpy()

    def predict_state(self, X):
        return self.state_encoder.inverse_transform(self.state_logits(X).argmax(axis=1))

    def predict_change(self, X):
        return self.change_encoder.inverse_transform(self.change_logits(X).argmax(axis=1))

    def predict_proba_state(self, X):
        return _softmax(self.state_logits(X))

    def predict_proba_change(self, X):
        return _softmax(self.change_logits(X))

    def embed(self, X):
        """Latent representation from the shared trunk, for inspecting what the encoder learned."""
        self.model.eval()
        with torch.no_grad():
            return self.model.trunk(torch.from_numpy(self._scale(X))).numpy()


def _softmax(logits):
    e = np.exp(logits - logits.max(axis=1, keepdims=True))
    return e / e.sum(axis=1, keepdims=True)


def _loader(dataset, batch_size, shuffle):
    """DataLoader that never yields a size-1 batch, which BatchNorm1d cannot handle."""
    batch_size = max(2, min(batch_size, len(dataset)))
    return DataLoader(dataset, batch_size=batch_size, shuffle=shuffle,
                      drop_last=shuffle and len(dataset) > batch_size)


def _endless(loader):
    """Yield batches forever, re-iterating the loader so each pass reshuffles."""
    while True:
        for batch in loader:
            yield batch


def _even_steps(steps, n_batches):
    """Evenly spaced step indices at which one head contributes to the loss.

    Spreading a head's batches across the epoch rather than bunching them at the start keeps the
    two heads interleaved, so the shared encoder never spends a long run of steps seeing only one
    task and drifting toward it.
    """
    if not n_batches:
        return set()
    if n_batches >= steps:
        return set(range(steps))
    return set(np.linspace(0, steps - 1, n_batches).round().astype(int).tolist())


def _grouped_val_split(x, groups, val_fraction, random_state):
    if groups is not None:
        splitter = GroupShuffleSplit(n_splits=1, test_size=val_fraction, random_state=random_state)
        return next(splitter.split(x, groups=np.asarray(groups)))
    splitter = ShuffleSplit(n_splits=1, test_size=val_fraction, random_state=random_state)
    return next(splitter.split(x))


def _stratified_val_split(x, y, val_fraction, random_state):
    """Stratified split, falling back to a random one when a class is too small to stratify."""
    try:
        splitter = StratifiedShuffleSplit(n_splits=1, test_size=val_fraction,
                                          random_state=random_state)
        return next(splitter.split(x, y))
    except ValueError:
        print('  change validation split: a class is too small to stratify, '
              'falling back to a random split')
        splitter = ShuffleSplit(n_splits=1, test_size=val_fraction, random_state=random_state)
        return next(splitter.split(x))


def _class_weights(y_encoded, n_classes, device, balance):
    if not balance:
        return None
    present = np.unique(y_encoded)
    weights = np.ones(n_classes, dtype=np.float32)
    weights[present] = compute_class_weight('balanced', classes=present, y=y_encoded)
    return torch.tensor(weights, dtype=torch.float32, device=device)


def fit_factored_change_model(
        state_x, state_y, change_pair, change_y, state_groups=None,
        state_weight=1.0, change_weight=1.0, pretrained=None, freeze_trunk=False,
        latent_dim=64, hidden_dims=(128,), change_hidden=(64,), dropout=0.3,
        epochs=100, batch_size=256, change_batch_size=128, lr=1e-3, weight_decay=0.0,
        balance_classes=True, val_fraction=0.2, early_stopping_rounds=None,
        early_stopping_metric='change', epoch_basis='balanced',
        random_state=1234, device=None, verbose=False):
    """Train the factored model with a masked multi-task loss.

    Each optimizer step draws a batch from whichever sources are scheduled for that step and sums
    their losses, so a GLanCE row only ever touches the state head and a reference point only ever
    touches the change head. Gradients from both still flow back through the shared trunk, which is
    where the two tasks regularize each other.

    The two sources differ in size by more than an order of magnitude, which makes the definition
    of an epoch load-bearing rather than cosmetic. See `epoch_basis`.

    Args:
        state_x, state_y: GLanCE predictors (the 64 embedding bands) and land-cover labels.
        change_pair: Paired change predictors, (n, 128), before block then after block.
        change_y: Change class ids aligned to change_pair.
        state_groups: Glance_ID per state row, so the state validation split stays grouped and
            no segment leaks between fit and validation.
        state_weight, change_weight: Loss weight per head. Setting one to 0 switches that head
            off entirely; its parameters then stay at initialization and its predictions are
            meaningless, so only read the head you trained.
        pretrained: Optional already-fitted model whose trunk weights and scaler are reused.
            The scaler comes along deliberately: a trunk is only valid for inputs scaled the
            way it was trained.
        freeze_trunk: If True, the trunk takes no gradient, so only the heads train. Combined
            with `pretrained` this is the frozen-encoder variant.
        latent_dim, hidden_dims, change_hidden, dropout: Architecture.
        epochs, batch_size, change_batch_size, lr, weight_decay: Training hyperparameters
            (epochs is the cap when early stopping is on).
        balance_classes: Weight each head's loss by inverse class frequency. Both label sets
            are heavily imbalanced, stable and Trees especially, so this is on by default.
        val_fraction: Fraction held out from each source for the loss curves.
        early_stopping_rounds: Stop when the watched validation loss has not improved for this
            many epochs, and restore the best epoch's weights.
        early_stopping_metric: Which validation loss to watch, 'change', 'state' or 'total'.
            Defaults to 'change' because the change head is the point of the model and has
            far less data, so it overfits long before the state head.
        epoch_basis: What one epoch means when both heads are active.
            'balanced' (default) gives each head exactly one pass over its own data per epoch,
            by having the smaller source contribute on an evenly spaced subset of the steps.
            'max' runs both heads on every step and cycles the smaller loader to fill the epoch.
            This is not a cosmetic choice. Under 'max' the change source contributes on all
            ~149 steps an epoch of GLanCE takes, which is ~12 passes over ~1,900 change points,
            so the change head overfits roughly 12 times faster per epoch than it does when
            trained alone. Early stopping counts epochs, so its resolution is 12 times too
            coarse to catch that: on the Oregon run the joint model's change validation optimum
            fell at about epoch 1.45 and early stopping was forced to return epoch 1, well off
            the optimum, which made the joint variant look worse than the single-task ones for
            reasons that had nothing to do with sharing the encoder. 'balanced' makes an epoch
            mean the same thing for both heads and for every variant, which is what makes the
            comparison valid. Note it also reduces the change head's share of the gradient
            reaching the shared trunk, which is arguably right since it tracks the relative
            data volumes; `change_weight` is the lever if the change task should pull harder.
        random_state: Seed for reproducibility.
        device: Torch device, defaults to MPS when available else CPU.
        verbose: Print losses every 10 epochs.

    Returns:
        Tuple of (fitted classifier, training time in seconds). The classifier exposes
        `.state` and `.change` head views for evaluate_model / compare_models, plus
        `.history_`, `.state_history_` and `.change_history_` for plot_training_curve.
    """
    assert state_weight > 0 or change_weight > 0, 'at least one head must be trained'
    assert early_stopping_metric in ('change', 'state', 'total'), \
        "early_stopping_metric must be 'change', 'state' or 'total'"

    torch.manual_seed(random_state)
    np.random.seed(random_state)
    if device is None:
        device = 'mps' if torch.backends.mps.is_available() else 'cpu'

    state_arr = np.asarray(state_x, dtype=np.float32)
    pair_arr = np.asarray(change_pair, dtype=np.float32)
    input_dim = state_arr.shape[1]
    assert pair_arr.shape[1] == 2 * input_dim, (
        f'change_pair should have {2 * input_dim} columns (before then after), '
        f'got {pair_arr.shape[1]}')
    before_arr, after_arr = pair_arr[:, :input_dim], pair_arr[:, input_dim:]

    # Encode both label sets on the full data so every class is represented even if a
    # validation split happens to miss one.
    state_encoder, change_encoder = LabelEncoder(), LabelEncoder()
    state_enc = state_encoder.fit_transform(state_y)
    change_enc = change_encoder.fit_transform(change_y)
    n_state, n_change = len(state_encoder.classes_), len(change_encoder.classes_)

    state_fit, state_val = _grouped_val_split(state_arr, state_groups, val_fraction, random_state)
    change_fit, change_val = _stratified_val_split(pair_arr, change_enc, val_fraction, random_state)

    # One scaler for both heads. They feed the same trunk, so they must arrive in the same
    # units; fitting on the union of every training-side vector, both years included, covers
    # the whole input distribution the trunk will see.
    if pretrained is not None:
        scaler = pretrained.scaler
    else:
        scaler = StandardScaler().fit(np.vstack([
            state_arr[state_fit], before_arr[change_fit], after_arr[change_fit]]))

    def scaled(a):
        return scaler.transform(a).astype(np.float32)

    state_fit_ds = TensorDataset(torch.from_numpy(scaled(state_arr[state_fit])),
                                 torch.from_numpy(state_enc[state_fit].astype(np.int64)))
    state_val_ds = TensorDataset(torch.from_numpy(scaled(state_arr[state_val])),
                                 torch.from_numpy(state_enc[state_val].astype(np.int64)))
    change_fit_ds = TensorDataset(torch.from_numpy(scaled(before_arr[change_fit])),
                                  torch.from_numpy(scaled(after_arr[change_fit])),
                                  torch.from_numpy(change_enc[change_fit].astype(np.int64)))
    change_val_ds = TensorDataset(torch.from_numpy(scaled(before_arr[change_val])),
                                  torch.from_numpy(scaled(after_arr[change_val])),
                                  torch.from_numpy(change_enc[change_val].astype(np.int64)))

    state_loader = _loader(state_fit_ds, batch_size, True)
    change_loader = _loader(change_fit_ds, change_batch_size, True)
    state_val_loader = _loader(state_val_ds, batch_size, False)
    change_val_loader = _loader(change_val_ds, change_batch_size, False)

    model = _FactoredChangeNet(input_dim, hidden_dims, latent_dim, change_hidden,
                               n_state, n_change, dropout).to(device)
    if pretrained is not None:
        model.trunk.load_state_dict(pretrained.model.trunk.state_dict())
    if freeze_trunk:
        for param in model.trunk.parameters():
            param.requires_grad = False

    state_loss_fn = nn.CrossEntropyLoss(
        weight=_class_weights(state_enc, n_state, device, balance_classes))
    change_loss_fn = nn.CrossEntropyLoss(
        weight=_class_weights(change_enc, n_change, device, balance_classes))
    optimizer = torch.optim.Adam(
        [p for p in model.parameters() if p.requires_grad], lr=lr, weight_decay=weight_decay)

    assert epoch_basis in ('balanced', 'max'), "epoch_basis must be 'balanced' or 'max'"
    n_state_batches = len(state_loader) if state_weight else 0
    n_change_batches = len(change_loader) if change_weight else 0
    steps = max(n_state_batches, n_change_batches)

    # Which steps each head contributes on. Under 'balanced' each head makes exactly one pass
    # over its own data per epoch, so an epoch means the same thing for both heads and across
    # every variant; under 'max' both run every step and the smaller loader cycles.
    if epoch_basis == 'balanced':
        state_steps = _even_steps(steps, n_state_batches)
        change_steps = _even_steps(steps, n_change_batches)
    else:
        state_steps = set(range(steps)) if n_state_batches else set()
        change_steps = set(range(steps)) if n_change_batches else set()

    state_batches, change_batches = _endless(state_loader), _endless(change_loader)
    if verbose:
        print(f'{steps} steps/epoch; land cover on {len(state_steps)}, '
              f'change on {len(change_steps)} (epoch_basis={epoch_basis!r})')

    def _curve(ylabel):
        return {'train': [], 'val': [], 'xlabel': 'Epoch', 'ylabel': ylabel}

    history, state_history, change_history = (
        _curve('Total loss'), _curve('State loss'), _curve('Change loss'))
    best_val, best_state_dict, best_epoch, no_improve = float('inf'), None, -1, 0

    start_time = time.perf_counter()
    for epoch in range(epochs):
        model.train()
        if freeze_trunk:
            # requires_grad=False only freezes the trunk's weights. Its BatchNorm running
            # means and variances are buffers, not parameters, so in train mode they keep
            # drifting toward the change data and the "frozen" representation quietly moves.
            # Holding the trunk in eval mode freezes those statistics too (and disables its
            # dropout), which is what makes this variant a real fixed-encoder control.
            model.trunk.eval()
        run_state = run_change = 0.0
        for step in range(steps):
            if step not in state_steps and step not in change_steps:
                continue
            optimizer.zero_grad()
            loss = torch.zeros((), device=device)
            if step in state_steps:
                xb, yb = next(state_batches)
                s_loss = state_loss_fn(model.forward_state(xb.to(device)), yb.to(device))
                loss = loss + state_weight * s_loss
                run_state += s_loss.item()
            if step in change_steps:
                x1, x2, yb = next(change_batches)
                c_loss = change_loss_fn(
                    model.forward_change(x1.to(device), x2.to(device)), yb.to(device))
                loss = loss + change_weight * c_loss
                run_change += c_loss.item()
            loss.backward()
            optimizer.step()

        # Each head's mean is over the steps it actually ran on, not over every step, so the
        # reported losses stay comparable between variants and across epoch_basis settings.
        train_state = run_state / max(1, len(state_steps))
        train_change = run_change / max(1, len(change_steps))
        state_history['train'].append(train_state)
        change_history['train'].append(train_change)
        history['train'].append(state_weight * train_state + change_weight * train_change)

        model.eval()
        with torch.no_grad():
            s_total = sum(state_loss_fn(model.forward_state(xb.to(device)),
                                        yb.to(device)).item() * len(xb)
                          for xb, yb in state_val_loader)
            c_total = sum(change_loss_fn(model.forward_change(x1.to(device), x2.to(device)),
                                         yb.to(device)).item() * len(yb)
                          for x1, x2, yb in change_val_loader)
        val_state, val_change = s_total / len(state_val_ds), c_total / len(change_val_ds)
        state_history['val'].append(val_state)
        change_history['val'].append(val_change)
        val_total = state_weight * val_state + change_weight * val_change
        history['val'].append(val_total)

        watched = {'change': val_change, 'state': val_state, 'total': val_total}[
            early_stopping_metric]
        if early_stopping_rounds is not None:
            if watched < best_val - 1e-6:
                best_val, best_epoch, no_improve = watched, epoch, 0
                best_state_dict = {k: v.detach().cpu().clone()
                                   for k, v in model.state_dict().items()}
            else:
                no_improve += 1

        if verbose and (epoch + 1) % 10 == 0:
            print(f'epoch {epoch + 1}/{epochs}  state {train_state:.4f}/{val_state:.4f}  '
                  f'change {train_change:.4f}/{val_change:.4f}  (train/val)')

        if early_stopping_rounds is not None and no_improve >= early_stopping_rounds:
            if verbose:
                print(f'early stopping at epoch {epoch + 1}; best epoch {best_epoch + 1} '
                      f'({early_stopping_metric} val {best_val:.4f})')
            break
    train_time = time.perf_counter() - start_time

    model.to('cpu').eval()
    if early_stopping_rounds is not None and best_state_dict is not None:
        model.load_state_dict(best_state_dict)

    clf = _FactoredChangeClassifier(model, state_encoder, change_encoder, scaler, input_dim)
    clf.history_ = history
    clf.state_history_ = state_history
    clf.change_history_ = change_history
    clf.best_epoch_ = best_epoch if (early_stopping_rounds is not None and best_epoch >= 0) else None
    return clf, train_time


def map_difference_labels(class_before, class_after, greenness=GREENNESS_ORDER):
    """Change direction implied by naively differencing two independently classified maps.

    This is the baseline the factored model exists to beat, so it is deliberately the honest
    version of the naive method rather than a straw man: same maps, same points, and the
    greenness ordering used to turn a class pair into a direction. What it structurally cannot
    do is distinguish abrupt from gradual, because a pair of annual snapshots carries no
    information about how the transition happened. That is a property of the method, not a
    shortcoming of this implementation, which is why the notebook compares it on the binary
    change-versus-stable question only.

    Args:
        class_before, class_after: Sequences of land-cover class names per point.
        greenness: Class names ordered from least to most vegetated.

    Returns:
        Array of 'stable', 'increase' or 'decrease' per point. Classes outside `greenness`
        are ranked equal, so an unrecognised pair reads as stable rather than inventing a
        direction.
    """
    rank = {name: i for i, name in enumerate(greenness)}
    before = np.array([rank.get(str(c), -1) for c in class_before])
    after = np.array([rank.get(str(c), -1) for c in class_after])
    out = np.full(len(before), 'stable', dtype=object)
    out[after > before] = 'increase'
    out[after < before] = 'decrease'
    return out


def binary_change(labels, stable_value='stable'):
    """Collapse any change labelling to the binary change-versus-stable gate.

    Args:
        labels: Sequence of labels, either change class ids or readable strings.
        stable_value: The value meaning "no change" in `labels`.

    Returns:
        Boolean array, True where the point is labelled as changed.
    """
    return np.asarray(labels) != stable_value


def change_gate_table(models, x_test, y_test, stable_id=0, threshold=0.5,
                      extra_predictions=None):
    """Binary change-versus-stable metrics, the question the carbon accounting actually asks.

    Whether a pixel changed at all decides whether any area, and therefore any carbon, is
    counted for it, so the gate deserves its own scoring rather than being buried inside a
    five-class macro average dominated by the stable class. Recall is how much real change
    survives; precision is how much of the counted area is fabricated, which is the direction
    map counting is known to err in.

    Args:
        models: Dict of display name -> (model, train_time) or (model, train_time,
            x_test_override), the same spec Evaluation.compare_models takes. predict() gives
            the hard call; predict_proba() supplies the ranking score for AUC when available.
        x_test, y_test: Shared paired predictors and true change class ids.
        stable_id: The class id meaning "no change".
        threshold: Probability of change above which a point is called changed. The hard call
            comes from thresholding the summed change probability, not from collapsing the
            multiclass argmax, and the two are not the same thing. With three classes at
            P(stable)=0.4, P(increase)=0.35, P(decrease)=0.25, the argmax is stable while the
            point is 60% likely to have changed. Collapsing the argmax would therefore
            under-detect change relative to the very score the ROC curve is drawn from, which
            would both understate recall and make finer-grained models look worse than they
            are. Models with no predict_proba fall back to the collapsed argmax.
        extra_predictions: Dict of name -> boolean "changed" array, for baselines that emit a
            hard label and no score (the map difference). These get every metric except AUC,
            since a single hard call is one operating point, not a ranking.

    Returns:
        Tuple of (table, curves, points, predicted):
          table: DataFrame indexed by model, one row of gate metrics each.
          curves: Dict of name -> (false positive rate, true positive rate) for models with
            a score, ready to plot as ROC curves.
          points: Dict of name -> (false positive rate, true positive rate) single operating
            points, for the hard-label baselines.
          predicted: Dict of name -> boolean "changed" array, the exact hard call every metric
            above was computed from. Returned so a caller can decompose those same calls
            further (by land cover class, by where a point was originally sampled from, and
            so on) without re-deriving the threshold-versus-argmax logic a second time.
    """
    from sklearn.metrics import precision_score, recall_score, f1_score, roc_auc_score, roc_curve

    truth = binary_change(np.asarray(y_test), stable_value=stable_id)
    n_stable = int((~truth).sum())
    rows, curves, points, predicted_out = {}, {}, {}, {}

    def _add(name, predicted, score=None, train_time=np.nan):
        false_pos = int((predicted & ~truth).sum())
        row = {
            'Precision (change)': precision_score(truth, predicted, zero_division=0),
            'Recall (change)': recall_score(truth, predicted, zero_division=0),
            'F1 (change)': f1_score(truth, predicted, zero_division=0),
            'False positive rate': false_pos / n_stable if n_stable else np.nan,
            'ROC AUC': roc_auc_score(truth, score) if score is not None else np.nan,
            'Training time (s)': train_time,
        }
        rows[name] = row
        predicted_out[name] = predicted
        if score is not None:
            fpr, tpr, _ = roc_curve(truth, score)
            curves[name] = (fpr, tpr)
        else:
            points[name] = (row['False positive rate'], row['Recall (change)'])

    for name, spec in models.items():
        model, train_time = spec[0], spec[1]
        x_eval = spec[2] if len(spec) > 2 else x_test
        score = None
        if hasattr(model, 'predict_proba'):
            proba = np.asarray(model.predict_proba(x_eval))
            classes = list(getattr(model, 'classes_', []))
            if stable_id in classes:
                score = 1.0 - proba[:, classes.index(stable_id)]
        # Threshold the summed change probability when there is one, so the hard call and the
        # ROC curve answer the same question; otherwise fall back to collapsing the argmax.
        predicted = (score >= threshold if score is not None
                     else binary_change(model.predict(x_eval), stable_value=stable_id))
        _add(name, predicted, score, train_time)

    for name, predicted in (extra_predictions or {}).items():
        _add(name, np.asarray(predicted, dtype=bool))

    return pd.DataFrame(rows).T, curves, points, predicted_out


def false_positive_rate_by_origin(y_test, meta_test, predictions, stable_id=0,
                                  origin_col='sample_type'):
    """False positive rate among stable-labelled test points, split by where each point was
    originally drawn from, to separate a genuine result from a tautology.

    Every reference point starts out tagged sample_type: 'stable' if the deployed classifier's
    2019 and 2025 maps already agreed on that pixel's class at sampling time, 'candidate' if
    they disagreed. A point later labelled stable by COLD can come from either stratum, and the
    two are not interchangeable for scoring naive map differencing. On a sample_type=='stable'
    point, class_2019 == class_2025 by construction, so a method that reads those two columns
    and compares them cannot produce a false positive there: agreement was the sampling
    criterion, not something being predicted. Its only real test is the sample_type=='candidate'
    points, where the maps originally disagreed and were later found, independently, to be
    wrong about that. Measured on the Oregon run this table was built for, 96.6% of the
    stable-labelled test set was the tautological stratum, which is why naive map differencing's
    overall false positive rate looked competitive with the trained models despite being wrong
    on every single informative point.

    This dilution affects any metric built from the false positive count (precision, F1, the
    false positive rate itself); it does not affect recall, which is computed on the changed
    population and untouched by how the stable population was assembled.

    Args:
        y_test: True change class ids for the test set, aligned to meta_test.
        meta_test: Test metadata; must include `origin_col` (and whatever each entry in
            `predictions` needs to have already been computed from).
        predictions: Dict of name -> boolean "changed" array, aligned to y_test/meta_test.
        stable_id: The class id meaning "no change".
        origin_col: Column recording each point's original sampling stratum.

    Returns:
        DataFrame indexed by origin stratum, one false-positive-rate column per name in
        `predictions` plus an 'n' column, so a reader can see both the rate and how much of
        the stable population that stratum represents.
    """
    stable_mask = (np.asarray(y_test) == stable_id)
    origin = np.asarray(meta_test.loc[stable_mask, origin_col])
    strata = sorted(pd.unique(origin), key=lambda o: -int((origin == o).sum()))

    table = pd.DataFrame(index=strata)
    table['n'] = [int((origin == o).sum()) for o in strata]
    for name, changed in predictions.items():
        changed = np.asarray(changed, dtype=bool)[stable_mask]
        table[name] = [changed[origin == o].mean() for o in strata]
    return table


def plot_change_gate_roc(curves, points=None, title='Change versus stable gate',
                         ax=None, save_path=None):
    """ROC curves for the binary change gate, with hard-label baselines drawn as single points.

    A model that ranks pixels by how likely they are to have changed traces a curve, so the
    operating point can be moved after the fact. A rule that emits one hard call, such as
    differencing two maps, has exactly one operating point and no such freedom. Plotting them
    on the same axes shows that difference rather than hiding it behind a single number.

    Args:
        curves: Dict of name -> (fpr, tpr) from change_gate_table.
        points: Dict of name -> (fpr, tpr) single operating points.
        title: Plot title.
        ax: Optional existing Axes.
        save_path: If given, save the figure here.

    Returns:
        The matplotlib Axes.
    """
    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots(figsize=(7, 6))
    for name, (fpr, tpr) in curves.items():
        ax.plot(fpr, tpr, linewidth=2, label=name)
    for name, (fpr, tpr) in (points or {}).items():
        ax.scatter([fpr], [tpr], s=90, marker='D', zorder=5, edgecolors='black',
                   linewidths=0.6, label=f'{name} (single operating point)')
    ax.plot([0, 1], [0, 1], linestyle=':', color='grey', linewidth=1, label='chance')
    ax.set_xlabel('False positive rate (stable land counted as change)')
    ax.set_ylabel('True positive rate (real change detected)')
    ax.set_xlim(0, 1)
    ax.set_ylim(0, 1.02)
    ax.set_title(title)
    ax.legend(loc='lower right', fontsize=9)
    ax.grid(alpha=0.3)
    ax.figure.tight_layout()

    if save_path is not None:
        ax.figure.savefig(save_path, bbox_inches='tight', dpi=150)
    return ax