[BUG] LSS scoring/evaluation pipeline: score() double-transforms and returns a dict; evaluate() ignores the family; default metrics read wrong columns

[ChrisW09]

Describe the bug
Umbrella issue: the LSS evaluation/scoring pipeline is broken end-to-end — the individual bugs compound, so filing them together since the fixes overlap in lss_base.py, metrics/registry.py, and metrics/distributional.py.

1. score() returns a dict and double-applies parameter transforms (deeptab/models/lss_base.py:472-473). It calls self.predict(X) (transformed params) and feeds them to family.evaluate_nll, whose compute_loss applies the transforms again (e.g. softplus(softplus(raw)) for Normal's scale). Verified: reported NLL 1.4155 vs correct 1.5467 on the same model/data. It also returns evaluate_nll's dict, not a float — breaking sklearn scoring (cross_val_score, GridSearchCV). Fix: self.predict(X, raw=True) + return a scalar.
2. evaluate() always falls back to normal-family metrics (lss_base.py:412): getattr(self._task_model, "distribution_family", "normal") — no such attribute exists anywhere (the value lives on self.family_name). Verified: a fitted Gamma model returns the normal metric set (Gaussian CRPS etc.) unless the user passes distribution_family explicitly.
3. Default LSS metrics read the wrong parameter columns (metrics/registry.py:35-50, distributional.py:224-231, 482-487). Column order comes from each family's param_names, but the metrics assume [loc, scale, ...]:
   • studentt head outputs [df, loc, scale] → CRPS uses loc=df, scale=loc; StudentTLoss uses mu=df, scale=loc, df=scale — fully scrambled.
   • johnsonsu outputs [skew, shape, location, scale] → CRPS reads loc=skew, scale=shape.
   • zip outputs [pi, rate] → PoissonDeviance/RMSE read π (a probability) as the mean.
   • gamma outputs [shape, rate] → deviance/RMSE read shape as the mean (mean is shape/rate).
   • quantile (default [0.25, 0.5, 0.75]) → PinballLoss(quantile=0.5) with col=0 scores the 0.25-quantile at τ=0.5.
     Since evaluate() converts metric exceptions to NaN + warning, none of this surfaces.
4. load() drops distributional_kwargs (lss_base.py:589, save side core/serialization.py:322): only the family name is persisted, so a model trained with custom quantiles / tweedie power / MoG components reloads with a default-configured family — head size survives, so it fails silently at score/interpretation time.
5. Standalone build_model() raises (lss_base.py:141): num_classes=self.family.param_count but self.family is only assigned inside fit() — the documented "call .build_model() or set rebuild=True" workflow is impossible for LSS models (AttributeError).
6. LSS fit never seeds RNGs: unlike _FitMixin.fit (which calls set_seed(random_state)), LSS fit uses random_state only for the split — MLPLSS(random_state=42) is not reproducible while MLPRegressor(random_state=42) is.
7. Family-specific defects: Categorical/Dirichlet declare param_count = 1, so the head gets one output and the family is unusable for K > 1; their forward also softmaxes across the batch dimension (distributions/categorical.py:24, beta.py:63, base.py:122-130). MixtureOfGaussians defines no transforms, so predict(raw=False) returns raw weight/scale logits (negative "sigmas"). PoissonDistribution.evaluate_nll computes y·log(y/rate) which is NaN whenever a true count is 0 (poisson.py:50).

To Reproduce

import numpy as np, pandas as pd
from deeptab.models import MLPLSS
X = pd.DataFrame({"a": np.random.randn(80)}); y = np.random.rand(80) + 0.5

m = MLPLSS(); m.fit(X, y, family="normal", max_epochs=2, accelerator="cpu")
print(type(m.score(X, y)))                       # dict, not float
raw  = np.asarray(m.predict(X, raw=True))
print(m._task_model.family.evaluate_nll(y, raw)) # NLL differs from score()'s
print(hasattr(m._task_model, "distribution_family"))  # False -> evaluate() uses "normal"

g = MLPLSS(); g.fit(X, y, family="gamma", max_epochs=2, accelerator="cpu")
print(g.evaluate(X, y))                          # normal-family metrics for a gamma model

Expected behavior
score() returns a scalar computed on raw parameters; evaluate() uses self.family_name; registry metrics index columns via each family's param_names; save bundles persist distributional_kwargs; build_model() works standalone; LSS fit seeds like the other estimators.

Screenshots
n/a

Desktop (please complete the following information):

• OS: macOS (Darwin 25.5.0, arm64)
• Python version: 3.11.15
• deeptab Version: 2.0.0 (main @ 4e6a359)

Additional context
Items 1–3 verified by execution; 4–7 by code inspection. The GammaDeviance sign error and NegativeBinomial parameterization are filed separately as self-contained fixes.