v0.6.0 — calibrated decoding (LowRankNoise)

The headline of 0.6 is a new per-voxel Gaussian likelihood option in
ParameterFitter, opt-in here and the default in 0.7. Same fixed points
as the existing sum-of-squares loss, much faster convergence under
heteroskedastic noise — which is to say, on every real fMRI dataset.

The second headline is calibrated decoding: braincoder.noise.LowRankNoise
(low-rank + diagonal noise covariance with cross-validated rank and
likelihood temperature) and braincoder.utils.calibration.

The rest of 0.6 is a substantial set of stability fixes, a new GP-prior
fitter, and a per-fold voxel-selection helper.

Default fitter changes (heads-up for downstream projects)

• ParameterFitter(noise_model={'ssq', 'gaussian'}) — new knob.
  • 'ssq' (default in 0.6): unweighted sum-of-squared-residuals.
    Equivalent to MLE under homoscedastic Gaussian noise. Same as 0.5.x.
  • 'gaussian' (new): per-voxel Gaussian negative log-likelihood with
    a free per-voxel σ²ᵥ. The θ-gradient is (1/σ²ᵥ)·∂SSQᵥ/∂θᵥ, so
    each voxel makes comparable progress under the shared Adam optimizer
    regardless of its noise level. Converges to the same per-voxel θ
    at high iterations — but does so much faster in practice on real
    data, where noise varies wildly across voxels.
• Synthetic verification: at max_n_iterations=20_000 both modes
  agree on per-voxel θ to <1e-3 normalized RMSE. At 2_000 iters
  Gaussian mode is at its asymptote; SSQ mode is still ~5× further
  out. See gp_prior_fmri/notes/analyses/classical_vs_ml_convergence/
  for the full convergence study.
• 0.7 will flip the default to 'gaussian'. If you need exact
  reproduction of pre-0.6 fits, pin noise_model='ssq' explicitly.
• ParameterFitter.r2_history_ / loss_history_ — the R² and loss
  trajectories of the last fit() call are now stored on the fitter.

Calibrated decoding: low-rank noise model with CV rank + temperature

• braincoder.noise.LowRankNoise (also braincoder.LowRankNoise) —
  noise covariance Ω = D (F Fᵀ + Ψ) D: k principal factors of the residual
  correlation plus diagonal uniquenesses (floored), Student-t (ν by ML) or
  Gaussian. Stored and evaluated in factored form (Woodbury + determinant
  lemma; no dense V × V inverse), so it scales to thousands of voxels.
• LowRankNoise.fit(model, data, paradigm, parameters, stimulus_range=grid, groups=runs) chooses the rank k and a likelihood temperature T
  (posterior ∝ likelihood^(1/T)) jointly by nested, group-wise CV within the
  training data on the decoding objective (held-out log posterior of the
  true stimulus). Fixed k= / temperature= also accepted. Encoding
  parameters are held fixed in the inner folds by default;
  inner_parameters= takes per-group parameters (dict) or a refit callable.
  noise.summary() / noise.cv_scores_ report the choice.
• get_stimulus_pdf(..., noise_model=, temperature=) — new optional
  arguments (also on EncodingRegressionModel). Defaults unchanged: without
  them the output is identical to 0.5.x. temperature also works with a
  dense omega=. New helper EncodingModel.predict_stimulus_grid.
• braincoder.noise.fit_lowrank_omega(residuals, k) → dense
  (omega, dof) for code that passes omega= / dof=.
• braincoder.utils.calibration — randomised PIT, HPD level of the
  truth, HDI coverage, log score vs uniform, marginalize,
  calibration_summary, plot_calibration; works on any
  get_stimulus_pdf DataFrame.
• Docs: new page docs/noise_models.rst, example
  examples/01_decoding_pipeline/04_lowrank_noise_calibration.py.
• Warnings documented there (from 7T single-trial numerosity data):
  ResidualFitter's one-factor Ω underestimates held-out noise (Mahalanobis
  ~4× too large; posteriors overconfident, all 38 participant × ROI worse
  than flat); lambd shrinkage is unsafe when voxels > trials; the
  pseudo-WWᵀ from init_pseudoWWT is uncentred (dominated by the
  common-mode response). ResidualFitter itself is unchanged.
• Tests: tests/test_noise_models.py (Woodbury vs dense, CV recovery of k,
  calibrated coverage, T ≈ 1 when correct and > 1 when tuning is
  misspecified, backward compatibility, all backends).

Hierarchical Bayesian / GP prior

A new hierarchical-Bayesian fitter for any encoding model, with a
Gaussian-Process prior over cortical geodesic distance.

• GeodesicGPPrior (braincoder.optimize.gp_prior): RBF kernel
  over a precomputed pairwise distance matrix. Adaptive jitter for
  Cholesky stability; classical MDS embedding so the kernel is PSD
  even on non-Euclidean graph distances (Schoenberg 1938); Cholesky
  and triangular solves in float64.
• BayesianParameterFitter (braincoder.optimize.bayesian_fitter):
  three-stage empirical-Bayes (classical → hyperparams via MLE → MAP),
  with options for type-II joint MAP (joint_hyperparams=True,
  co-fits hyperparameters in stage 3), shared lengthscale across
  priors (shared_lengthscale=True), and frozen Cholesky in MAP
  to keep TF's CholeskyGrad out of the backward graph.

Voxel selection / FDR utilities

• fit_r2_mixture / r2_fdr_threshold (braincoder.utils.stats) —
  2-Gaussian mixture on logit(R²) with a tail-FDR threshold for
  empirical-null voxel selection.
• r2_posterior_signal / r2_p_signal_threshold — per-voxel
  posterior P(signal | R²) and the R² at which it crosses a chosen
  probability. More forgiving than tail-FDR on small/degenerate ROIs.
• plot_r2_mixture — diagnostic plot of the mixture fit + threshold.

Fisher information & expected uncertainty

• Analytical Fisher information for multivariate Student-t noise.
  EncodingModel.get_fisher_information(..., analytical=True) now
  accepts dof and applies the closed-form prefactor
  (ν + p) / (ν + p + 2) · Jᵀ Ω⁻¹ J. Previously this combination
  raised ValueError and forced users to the Monte-Carlo path; MC
  was both noisier (~12% RMS at n=200) and prone to int32 overflow in
  UnsortedSegmentSum at n × n_vox × n_stim ≈ 5e7. The analytical
  path is now exact for any ν and recovers the Gaussian formula as
  ν → ∞.
• EncodingModel.get_expected_uncertainty(stimuli, omega, dof, ...)
  — new convenience wrapper. Simulates n_simulations noisy responses
  per stimulus from (parameters, omega, dof), decodes each via
  get_stimulus_pdf, and returns a DataFrame indexed by stimulus
  value with mean_E, var_E, mean_error, mean_abs_error, n_sims.
  batch_stimuli lets you bound memory at large grid sizes.
• Tests: tests/test_fisher_information.py (9 tests): Student-t
  prefactor, dof → ∞ recovery of Gaussian, MC ↔ analytical agreement,
  and the FI ↔ 1/var_E relationship.
• Docs: new docs/fisher_information.rst page with the two-API
  story, the "is spikiness real?" diagnostic, and a worked example.

Cortical surface helpers

• braincoder.utils.cortex.geodesic_distance_matrix — Dijkstra-based
  pairwise geodesic distance on a triangulated mesh, restricted to a
  set of source vertices (typical use: project EPI-mask centroids to
  nearest vertex then run Dijkstra on the matched set).

Robustness fixes

• safe_cholesky(M, jitter=1e-4) (braincoder.utils.backend) —
  symmetrise + adaptive-jitter wrapper for ops.cholesky with a
  retry-on-NaN at 10× jitter. Wired into get_stimulus_pdf, Fisher
  information, both _simulate noise paths, and the residual fitter
  likelihood. Eliminates the Cholesky : Tensor had NaN values class
  of crashes that hit ~10–20% of subjects on production fits.
• ResidualFitter.fit(use_wwt=True) — when False, strips the
  σ²·WᵀW tuning-similarity term from Ω. Diagnostic; default
  unchanged.
• BayesianParameterFitter.fit_map uses a manual NaN-safe global-
  norm gradient clip instead of keras.optimizers.Adam(clipnorm=...),
  which under PyTorch would poison Adam's moment buffers with NaN if
  any single gradient went non-finite.
• ResidualFitter.fit no longer swallows every exception. It used
  to catch any exception, record cost=inf, back off the learning rate
  and carry on, so a fit that failed at every step (e.g. the old
  _lgamma under JAX) returned the initial Ω and dof as if it had
  converged. Now: programming errors (TypeError, which includes JAX
  tracer-conversion errors, plus AttributeError, NameError, …) are
  re-raised immediately; numerical failures (exceptions such as a
  Cholesky failure, non-finite cost or gradients) are still rolled back
  with a 10% learning-rate cut, but max_consecutive_failures=50 in a
  row raises RuntimeError, and so does a fit that never reaches a
  finite cost. A step with non-finite gradients is no longer applied
  (previously NaN gradients could poison the variables). Failed steps
  emit a RuntimeWarning and are counted in fitter.n_failed_steps.
  The stored best values now correspond to the parameters at which the
  best cost was evaluated (previously the post-update values were
  stored, one Adam step past the best cost).

Backend layer fixes (TF / JAX / torch parity)

• TF compute_jacobian uses tf.while_loop instead of pfor, which
  overflowed int32 indexing on large (trials × voxels × parameters) problems.

• Monte-Carlo Fisher information sums the log-probability before taking the
  gradient, so analytical=False works under JAX / torch (previously it
  raised on the non-scalar loss; TF summed it silently).

• JAX compute_gradients runs inside a StatelessScope so
  keras.Variable reads inside value_and_grad don't leak tracers.

• sample_mvt / sample_student_t derive independent sub-seeds for
  the χ² and normal components.

• _lgamma routes through differentiable backend-native
  lgamma (tf.math.lgamma / jax.scipy.special.gammaln /
  torch.lgamma) instead of the non-differentiable scipy fallback.
  Downstream impact: Student-t dof fitted on TF before this fix
  is biased low (the gradient was missing
  0.5[ψ((ν+k)/2) − ψ(ν/2)] per observation); on JAX the residual
  fit silently never ran. Refit Ω/dof for any decoding results that
  used method='t' with an older checkout.

• torch compute_gradients uses torch.autograd.grad instead of
  reading .grad off each variable.

• get_stimulus_pdf converts user data arrays to tensors before
  _likelihood (avoids numpy-on-MPS device-mismatch errors).

• safe_cholesky explicitly checks for NaN in the output, since some
  backends return NaN silently rather than raising.

Tests & CI

• New test files: test_gp_prior.py (+360 lines), test_safe_cholesky.py
  (+174), test_backend_multibackend.py (+376), additions to
  test_utils.py (+168).
• scikit-learn now declared as a runtime dep (used by
  fit_r2_mixture).
• TFP-dependent import test is skipped when tensorflow_probability is
  not installed.

Compatibility

• numpy.trapz (deprecated in NumPy 2.0) replaced with
  numpy.trapezoid via a shim.

Docs: https://braincoder-devs.github.io/ (What's new: https://braincoder-devs.github.io/whats_new.html, Noise models for decoding: https://braincoder-devs.github.io/noise_models.html)