[al4225]

Bug fixes (merged via fix-bug branch)

fix(posterior_mixture) — score collapsed pure-slab laplacemix as slab

When pi_null > 1 - tol_null_prior, the mixture collapsed to the null
component and the log posterior was computed using slab_log_density, but the
weight assignment logic treated it as null. The sign now matches the weight.

fix(em_helpers) — rescale mixsqp likelihood rows to avoid underflow

The raw dnorm evaluations could all underflow to zero for a row with very
small Shat, leaving mixsqp a zero row (which it rejects). Rows are now
log-space-shifted by rowMaxs before exponentiation so no row underflows.
rowMaxs now passes na.rm = TRUE (C2 fix; see below).

fix(iter_cache) — apply low-count mask to cached M-step sdmat

The cached sdmat in refresh_iter_cache was not masked by lowc_idx,
causing the M-step to see stale variance estimates for low-count wavelet
columns. The mask is now applied before caching. (Made moot by the P3
sdmat-cache removal described below, but was a correctness bug while present.)

fix(qnorm) — restore caller's RNG state after quantile normalization

mf_quantile_normalize called set.seed() without saving/restoring
.Random.seed, corrupting the caller's RNG stream in reproducibility-critical
paths (permutation tests, simulation scripts). State is now saved and restored
with on.exit.

fix(lowc) — ignore NA rows when computing low-count column medians

mf_lowcount_idx used colMedians without na.rm = TRUE. When missing-data
rows were present the median was NA, which propagated into the low-count mask
and silenced the entire wavelet level.

fix(adjust_covariates) — handle missing rows in OLS path

mf_adjust_for_covariates called the OLS path on the full n × T residual
matrix even when some rows were missing (NA), causing lm.fit to return
coefficients for the wrong design matrix. Missing rows are now excluded before
the OLS solve and re-inserted as NA in the output.

fix(cs_filter) — keep CS fields aligned under cs_min_lead_pip

When cs_min_lead_pip filtered out a credible set, the bookkeeping vectors
(cs_activity, cs_purity, cs_size) were trimmed but cs_index was not,
leaving them misaligned. All four vectors are now trimmed together.

fix(hmm) — guard slab loops for a single state and set G_t[1]

With a one-state HMM the inner loops over slab transitions were entered with
zero iterations, leaving G_t uninitialized. Added an explicit guard and
initializes G_t[1] to the null-component log probability.

fix(em) — remove P-scaling from mixture_null_weight in mixsqp M-step

The null-component penalty in mf_em_m_step_per_scale was multiplied by
p * idx_size (number of data points), making the prior increasingly
conservative as P grew. The penalty is now a fixed constant independent of P,
matching SuSiE's treatment of the null-component weight.

---

Algorithm refactor

refactor(wavelet_eb) — replace mixsqp M-step with exact EM responsibility update

mf_residualize_wavelet_eb previously ran a mixsqp sub-problem at each scale
for the covariate-adjustment prior. Ported to the closed-form exact EM update:

• New helper .wavelet_eb_normal_means fuses the E-step (posterior mean,
  variance) and responsibility accumulation in one sweep.
• M-step becomes colSums(responsibilities) / total, which is exact and
  requires no optimizer.
• Low-count (lowc_idx) wavelet columns are excluded from responsibility
  accumulation and sigma2 estimation, matching prior_groups semantics.
• Removes the ashr import from this path (postmean, postsd,
  set_data, get_fitted_g no longer needed).

Numerical result: fitted function correctly shrunk vs OLS (sin signal recovers
[-0.95, 0.97] vs OLS [-1.03, 1.04]). All 23 adjust_covariates tests pass.

---

Performance improvements

perf(bhat_shat,iter_cache) — closed-form marginal Bhat/Shat; remove sdmat cache (P2+P3)

P2: mf_compute_marginal_bhat_shat replaces T separate n × p allocations
with one crossprod and sufficient-statistic formulas (Sxx, Sxy, Syy).
At n=3000, p=2000, T=192 this is ~260× faster (22.9 s → 0.09 s per outcome).
All call sites updated.

P3: refresh_iter_cache.mf_individual no longer builds sdmat/log_sdmat
(previously 117 MB persistent, 252 MB per-iteration allocation).
mf_em_likelihood_per_scale rebuilds these on the fly in O(pK) when the cache
is NULL. Benchmark: M-step wall time 2.78 s → 2.18 s; mixture weights
identical to 1.4 × 10⁻¹⁷.

perf(ser_loglik) — replace X·Eb matrix multiply with cached XtR dot product

SER_posterior_e_loglik.mf_individual computed sum_n R[n,t] * (X Eb)[n,t]
via an explicit (X %*% Eb_m)[idx_m, ] — O(n p T) per (l, m) call.

The algebraic identity sum_n R[n,t]*(X Eb)[n,t] = sum_j Eb[j,t]*(X'R)[j,t]
reduces this to colSums(Eb_m * model$residuals[[m]]), which is O(p T) and
reuses the X'R cache already built by compute_residuals. At n = 229 this is
a ~229× per-call saving. Because the function is called up to
(1 + inner_em_steps) × L × M times per IBSS iteration the gain compounds.

Numerical diff: relative error ≤ 10⁻¹³ (verified via bench_identities).
247 core tests pass.

---

Correctness fixes bundled in edbf44b

• sigma2 lower bound (C1): update_variance_components now clamps sigma2
  to max(value, residual_variance_lowerbound %||% 1e-8) in both the
  per_outcome and per_scale paths. Prevents log(0) / NaN propagation when
  the bias-correction term briefly pushes the residual variance estimate below
  zero due to floating-point cancellation.
• zero-pw guard (C3): mf_compute_marginal_bhat_shat now handles predictors
  whose column norm is zero (e.g., constant post-centering). Sets Bhat = 0,
  Shat = 1 for those columns, consistent with mf_per_outcome_bhat_shat.
• Dead code removal: Variables p_full, cache_row_keep, cache_subset,
  cache_m_sdmat, cache_m_log_sdmat removed from .opv_mixsqp; they were
  always NULL after the P3 sdmat-cache removal.

---

Test coverage

• 247 core SER/IBSS tests pass after all changes.
• 23 adjust_covariates / wavelet-EB tests pass after the exact-EM refactor.
• End-to-end full-fit: PASS count 1577 → 1579; 22 pre-existing failures are
  reference-package version drift unrelated to these changes.