use crate::estimate::EstimationError;
use crate::estimate::reml::FirthDenseOperator;
use crate::faer_ndarray::{
FaerCholesky, FaerEigh, FaerLinalgError, FaerSymmetricFactor, array1_to_col_matmut, fast_ab,
fast_av, fast_av_into,
};
use crate::linalg::sparse_exact::{factorize_sparse_spd, solve_sparse_spd};
use crate::linalg::utils::{StableSolver, boundary_hit_step_fraction};
use crate::matrix::{DesignMatrix, LinearOperator};
use crate::mixture_link::{InverseLinkJet as MixtureInverseLinkJet, logit_inverse_link_jet5};
use crate::solver::active_set;
use crate::types::{Coefficients, LinearPredictor, StandardLink};
use crate::types::{
GlmLikelihoodSpec, InverseLink, LikelihoodSpec, LinkFunction, MixtureLinkState, ResponseFamily,
SasLinkState, is_valid_tweedie_power,
};
use dyn_stack::{MemBuffer, MemStack};
use faer::sparse::SparseColMat;
use faer::sparse::linalg::matmul::{
SparseMatMulInfo, sparse_sparse_matmul_numeric, sparse_sparse_matmul_numeric_scratch,
sparse_sparse_matmul_symbolic,
};
use faer::sparse::{
SparseColMatMut, SparseColMatRef, SparseRowMat, SymbolicSparseColMat, SymbolicSparseColMatRef,
};
use faer::{Accum, Par, Side, Unbind, get_global_parallelism};
use log;
use ndarray::{Array1, Array2, ArrayView1, ArrayView2, ArrayView3, ShapeBuilder, Zip};
use rayon::iter::{IndexedParallelIterator, IntoParallelRefMutIterator, ParallelIterator};
use statrs::function::gamma::{digamma, ln_gamma};
use faer::linalg::cholesky::llt::factor::LltParams;
use faer::{Auto, Spec};
use std::collections::BTreeSet;
use std::collections::HashMap;
use std::sync::{Arc, Mutex, OnceLock};
pub use crate::solver::active_set::{ConstraintKktDiagnostics, LinearInequalityConstraints};
use crate::linalg::utils::{array_is_finite, inf_norm};
mod convergence;
mod damping;
mod edf;
mod gpu_dispatch;
mod loop_driver;
mod penalty;
mod pls_solver;
mod reweight;
mod state;
use convergence::effective_kkt_tolerance;
use damping::{
add_scaled_diagonal_to_upper_sparse, compute_lm_d2, update_scaled_diagonal_in_place,
};
pub use edf::StablePLSResult;
use edf::{
calculate_edf_from_sparse_factor, calculate_edf_with_penalty,
calculate_edfwithworkspace_from_factor, calculate_edfwithworkspace_with_penalty,
};
use penalty::{
KroneckerQsTransform, PirlsPenalty, WorkingCoordinateDesign, WorkingReparamTransform,
attach_penalty_shift,
};
use pls_solver::solve_penalized_least_squares_implicit;
pub use pls_solver::{GaussianFixedCache, SparseXtwxPrecomputed};
pub use reweight::runworking_model_pirls;
pub(crate) use state::array1_l2_norm;
pub use state::{
AdaptiveKktTolerance, ExportedLaplaceCurvature, FirthDiagnostics, HessianCurvatureKind,
PirlsCoordinateFrame, PirlsLinearSolvePath, PirlsResult, PirlsStatus,
WorkingModelIterationInfo, WorkingModelPirlsResult, WorkingState,
};
const GAMMA_SHAPE_MIN: f64 = 1e-8;
const GAMMA_SHAPE_MAX: f64 = 1e12;
const GAMMA_SHAPE_TARGET_TOL: f64 = 1e-12;
pub(super) const PIRLS_ETA_ABS_CAP: f64 = 40.0;
#[inline]
fn gamma_shape_score(shape: f64, target: f64) -> f64 {
shape.ln() - digamma(shape) - target
}
fn estimate_gamma_shape_from_eta(
y: ArrayView1<'_, f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<'_, f64>,
) -> f64 {
const EPS: f64 = 1e-12;
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let (weighted_target, total_weight) = (0..eta.len())
.into_par_iter()
.map(|i| {
let wi = priorweights[i].max(0.0);
if wi == 0.0 {
return (0.0_f64, 0.0_f64);
}
let yi = y[i].max(EPS);
let mui = eta[i].clamp(-700.0, 700.0).exp().max(EPS);
let ratio = yi / mui;
(wi * (ratio - ratio.ln() - 1.0), wi)
})
.reduce(
|| (0.0_f64, 0.0_f64),
|(t1, w1), (t2, w2)| (t1 + t2, w1 + w2),
);
if total_weight <= 0.0 {
return 1.0;
}
let target = (weighted_target / total_weight).max(0.0);
if target <= GAMMA_SHAPE_TARGET_TOL {
return GAMMA_SHAPE_MAX;
}
let discriminant = (target - 3.0) * (target - 3.0) + 24.0 * target;
let approx = ((3.0 - target) + discriminant.sqrt()) / (12.0 * target);
let mut lo = GAMMA_SHAPE_MIN;
let mut hi = approx.max(1.0);
while hi < GAMMA_SHAPE_MAX && gamma_shape_score(hi, target) > 0.0 {
hi = (hi * 2.0).min(GAMMA_SHAPE_MAX);
}
if gamma_shape_score(hi, target) > 0.0 {
return GAMMA_SHAPE_MAX;
}
for _ in 0..80 {
let mid = 0.5 * (lo + hi);
if gamma_shape_score(mid, target) > 0.0 {
lo = mid;
} else {
hi = mid;
}
if (hi - lo) <= GAMMA_SHAPE_TARGET_TOL * hi.max(1.0) {
break;
}
}
0.5 * (lo + hi)
}
#[derive(Clone, Debug)]
pub struct SparsePirlsDecision {
pub path: PirlsLinearSolvePath,
pub reason: &'static str,
pub p: usize,
pub nnz_x: usize,
pub nnz_xtwx_symbolic: Option<usize>,
pub nnz_s_lambda: usize,
pub nnz_h_est: Option<usize>,
pub density_h_est: Option<f64>,
}
fn fmt_opt_usize(v: Option<usize>) -> String {
v.map(|v| v.to_string()).unwrap_or_else(|| "na".to_string())
}
fn fmt_opt_f64(v: Option<f64>) -> String {
v.map(|v| format!("{v:.4}"))
.unwrap_or_else(|| "na".to_string())
}
impl SparsePirlsDecision {
fn path_str(&self) -> &'static str {
match self.path {
PirlsLinearSolvePath::DenseTransformed => "dense_transformed",
PirlsLinearSolvePath::SparseNative => "sparse_native",
}
}
fn format_fields(&self, path: &str) -> String {
format!(
"path={path} reason={} p={} nnz_x={} nnz_xtwx_symbolic={} nnz_s_lambda={} nnz_h_est={} density_h_est={}",
self.reason,
self.p,
self.nnz_x,
fmt_opt_usize(self.nnz_xtwx_symbolic),
self.nnz_s_lambda,
fmt_opt_usize(self.nnz_h_est),
fmt_opt_f64(self.density_h_est),
)
}
fn log_once(&self) {
let path = self.path_str();
let key = self.format_fields(path);
let repetition_count = pirls_decision_repetition_count(key.clone());
if repetition_count == 1 {
log::debug!("[pirls-path] {key}");
return;
}
if should_log_pirls_decision_summary(repetition_count) {
log::debug!(
"[pirls-path] repeated path={} reason={} count={} (suppressing identical decisions)",
path,
self.reason,
repetition_count,
);
}
}
}
fn pirls_decision_repetition_count(log_key: String) -> usize {
static PIRLS_DECISION_LOG_COUNTS: OnceLock<Mutex<HashMap<String, usize>>> = OnceLock::new();
let counts = PIRLS_DECISION_LOG_COUNTS.get_or_init(|| Mutex::new(HashMap::new()));
let mut counts = counts.lock().expect("pirls decision log counter poisoned");
let count = counts.entry(log_key).or_insert(0);
*count += 1;
*count
}
fn should_log_pirls_decision_summary(repetition_count: usize) -> bool {
repetition_count > 1 && repetition_count.is_power_of_two()
}
const SPARSE_NATIVE_MAX_H_DENSITY: f64 = 0.30;
#[derive(Clone, Debug)]
struct SparsePenaltyPattern {
upper_triplets: Vec<(usize, usize, f64)>,
nnz_upper: usize,
}
impl SparsePenaltyPattern {
fn from_dense_upper(matrix: &Array2<f64>, tol: f64) -> Self {
let p = matrix.nrows().min(matrix.ncols());
let mut upper_triplets = Vec::new();
for col in 0..p {
for row in 0..=col {
let value = matrix[[row, col]];
if value.abs() > tol {
upper_triplets.push((row, col, value));
}
}
}
let nnz_upper = upper_triplets.len();
Self {
upper_triplets,
nnz_upper,
}
}
}
#[derive(Clone, Debug)]
pub(crate) struct SparsePenalizedSystemStats {
pub(crate) nnz_xtwx_symbolic: usize,
pub(crate) nnz_s_lambda_upper: usize,
pub(crate) nnz_h_upper: usize,
pub(crate) density_upper: f64,
}
struct SparsePenalizedSystemCache {
xtwx_cache: SparseXtWxCache,
penalty_pattern: SparsePenaltyPattern,
h_upper_symbolic: SymbolicSparseColMat<usize>,
h_uppervalues: Vec<f64>,
h_upper_col_ptr: Vec<usize>,
h_upperrow_idx: Vec<usize>,
p: usize,
}
impl SparsePenalizedSystemCache {
fn new(
x: &SparseColMat<usize, f64>,
penalty_pattern: SparsePenaltyPattern,
) -> Result<Self, EstimationError> {
let xtwx_cache = SparseXtWxCache::new(x)?;
let p = x.ncols();
let h_upper_symbolic = build_penalized_symbolic(
p,
xtwx_cache.xtwx_symbolic.col_ptr(),
xtwx_cache.xtwx_symbolic.row_idx(),
&penalty_pattern.upper_triplets,
)?;
let h_uppervalues = vec![0.0; h_upper_symbolic.row_idx().len()];
Ok(Self {
xtwx_cache,
penalty_pattern,
h_upper_col_ptr: h_upper_symbolic.col_ptr().to_vec(),
h_upperrow_idx: h_upper_symbolic.row_idx().to_vec(),
h_upper_symbolic,
h_uppervalues,
p,
})
}
fn matches(
&self,
x: &SparseColMat<usize, f64>,
penalty_pattern: &SparsePenaltyPattern,
) -> bool {
self.xtwx_cache.matches(x)
&& self.penalty_pattern.nnz_upper == penalty_pattern.nnz_upper
&& self.penalty_pattern.upper_triplets == penalty_pattern.upper_triplets
}
fn stats(&self) -> SparsePenalizedSystemStats {
let upper_total = self.p.saturating_mul(self.p + 1) / 2;
SparsePenalizedSystemStats {
nnz_xtwx_symbolic: self.xtwx_cache.xtwx_symbolic.row_idx().len(),
nnz_s_lambda_upper: self.penalty_pattern.nnz_upper,
nnz_h_upper: self.h_upper_symbolic.row_idx().len(),
density_upper: if upper_total == 0 {
0.0
} else {
self.h_upper_symbolic.row_idx().len() as f64 / upper_total as f64
},
}
}
fn assemble_upper(
&mut self,
x: &SparseColMat<usize, f64>,
weights: &Array1<f64>,
ridge: f64,
precomputed_xtwx: Option<&SparseXtwxPrecomputed>,
) -> Result<SparseColMat<usize, f64>, EstimationError> {
if weights.len() != self.xtwx_cache.nrows {
crate::bail_invalid_estim!(
"weights length {} does not match design rows {}",
weights.len(),
self.xtwx_cache.nrows
);
}
let use_precomputed = match precomputed_xtwx {
Some(pre) => {
let col_ptr_ok =
pre.xtwx_symbolic_col_ptr.as_slice() == self.xtwx_cache.xtwx_symbolic.col_ptr();
let row_idx_ok =
pre.xtwx_symbolic_row_idx.as_slice() == self.xtwx_cache.xtwx_symbolic.row_idx();
let values_ok = pre.xtwxvalues.len() == self.xtwx_cache.xtwxvalues.len();
if col_ptr_ok && row_idx_ok && values_ok {
self.xtwx_cache.xtwxvalues.copy_from_slice(&pre.xtwxvalues);
true
} else {
log::warn!(
"[sparse-xtwx-cache] precomputed XᵀWX pattern mismatch; \
falling back to per-call recompute"
);
false
}
}
None => false,
};
if !use_precomputed {
self.xtwx_cache.compute_numeric(x, weights)?;
}
self.h_uppervalues.fill(0.0);
let mut cursor = self.h_upper_col_ptr[..self.p].to_vec();
let xtwx_col_ptr = self.xtwx_cache.xtwx_symbolic.col_ptr();
let xtwxrow_idx = self.xtwx_cache.xtwx_symbolic.row_idx();
for col in 0..self.p {
let start = xtwx_col_ptr[col];
let end = xtwx_col_ptr[col + 1];
for idx in start..end {
let row = xtwxrow_idx[idx];
if row <= col {
let cursor_idx = &mut cursor[col];
while *cursor_idx < self.h_upper_col_ptr[col + 1]
&& self.h_upperrow_idx[*cursor_idx] < row
{
*cursor_idx += 1;
}
if *cursor_idx >= self.h_upper_col_ptr[col + 1]
|| self.h_upperrow_idx[*cursor_idx] != row
{
crate::bail_invalid_estim!("penalized symbolic pattern missing XtWX entry");
}
self.h_uppervalues[*cursor_idx] += self.xtwx_cache.xtwxvalues[idx];
}
}
}
cursor.copy_from_slice(&self.h_upper_col_ptr[..self.p]);
for &(row, col, value) in &self.penalty_pattern.upper_triplets {
let cursor_idx = &mut cursor[col];
while *cursor_idx < self.h_upper_col_ptr[col + 1]
&& self.h_upperrow_idx[*cursor_idx] < row
{
*cursor_idx += 1;
}
if *cursor_idx >= self.h_upper_col_ptr[col + 1]
|| self.h_upperrow_idx[*cursor_idx] != row
{
crate::bail_invalid_estim!("penalized symbolic pattern missing penalty entry");
}
self.h_uppervalues[*cursor_idx] += value;
}
if ridge > 0.0 {
cursor.copy_from_slice(&self.h_upper_col_ptr[..self.p]);
for col in 0..self.p {
let cursor_idx = &mut cursor[col];
while *cursor_idx < self.h_upper_col_ptr[col + 1]
&& self.h_upperrow_idx[*cursor_idx] < col
{
*cursor_idx += 1;
}
if *cursor_idx >= self.h_upper_col_ptr[col + 1]
|| self.h_upperrow_idx[*cursor_idx] != col
{
crate::bail_invalid_estim!("penalized symbolic pattern missing diagonal entry");
}
self.h_uppervalues[*cursor_idx] += ridge;
}
}
Ok(SparseColMat::new(
self.h_upper_symbolic.clone(),
self.h_uppervalues.clone(),
))
}
}
fn build_penalized_symbolic(
p: usize,
xtwx_col_ptr: &[usize],
xtwxrow_idx: &[usize],
penalty_triplets: &[(usize, usize, f64)],
) -> Result<SymbolicSparseColMat<usize>, EstimationError> {
let mut cols: Vec<BTreeSet<usize>> = (0..p).map(|_| BTreeSet::new()).collect();
for col in 0..p {
cols[col].insert(col);
let start = xtwx_col_ptr[col];
let end = xtwx_col_ptr[col + 1];
for &row in &xtwxrow_idx[start..end] {
if row <= col {
cols[col].insert(row);
}
}
}
for &(row, col, _) in penalty_triplets {
if row > col || col >= p {
crate::bail_invalid_estim!(
"penalty sparse pattern must be upper-triangular within bounds"
);
}
cols[col].insert(row);
}
let mut col_ptr = Vec::with_capacity(p + 1);
let mut row_idx = Vec::new();
col_ptr.push(0);
for rows in cols {
row_idx.extend(rows.into_iter());
col_ptr.push(row_idx.len());
}
Ok(unsafe { SymbolicSparseColMat::new_unchecked(p, p, col_ptr, None, row_idx) })
}
pub trait WorkingModel {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError>;
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature_kind: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
assert!(core::mem::size_of_val(&curvature_kind) > 0);
self.update(beta)
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, curvature)
}
fn screen_candidate(
&mut self,
beta: &Coefficients,
arr: &Array1<f64>,
linear_predictor: &LinearPredictor,
curvature: HessianCurvatureKind,
) -> Result<CandidateEvaluation, EstimationError> {
assert!(arr.iter().all(|v| !v.is_nan()));
assert!(std::mem::size_of_val(linear_predictor) > 0);
self.update_candidate(beta, curvature)
.map(CandidateEvaluation::Full)
}
fn supports_observed_information_curvature(&self) -> bool {
false
}
}
#[derive(Debug, Clone)]
pub struct CandidateScreen {
pub penalized_objective: f64,
pub deviance: f64,
pub penalty_term: f64,
pub arithmetic_finite: bool,
}
pub enum CandidateEvaluation {
Screen(CandidateScreen),
Full(WorkingState),
}
impl CandidateEvaluation {
#[inline]
fn penalized_objective(&self, firth_bias_reduction: bool) -> f64 {
match self {
Self::Screen(s) => s.penalized_objective,
Self::Full(state) => {
let mut value = state.deviance + state.penalty_term;
if firth_bias_reduction && let Some(j) = state.jeffreys_logdet() {
value -= 2.0 * j;
}
value
}
}
}
#[inline]
fn arithmetic_finite(&self) -> bool {
match self {
Self::Screen(s) => s.arithmetic_finite,
Self::Full(state) => state.gradient.iter().all(|g| g.is_finite()),
}
}
#[inline]
fn into_full(self) -> Option<WorkingState> {
match self {
Self::Full(state) => Some(state),
Self::Screen(_) => None,
}
}
}
#[derive(Clone, Debug, PartialEq, Eq)]
pub(super) struct PirlsAcceptedStateCacheKey {
curvature: HessianCurvatureKind,
firth_active: bool,
beta_bits: Vec<u64>,
arrow_latent_bits: Option<Vec<u64>>,
}
impl PirlsAcceptedStateCacheKey {
fn requested(
beta: &Coefficients,
curvature: HessianCurvatureKind,
options: &WorkingModelPirlsOptions,
) -> Self {
Self::new(beta, curvature, options.firth_bias_reduction, options)
}
fn accepted(
beta: &Coefficients,
state: &WorkingState,
options: &WorkingModelPirlsOptions,
) -> Self {
Self::new(
beta,
state.hessian_curvature,
matches!(state.firth, FirthDiagnostics::Active { .. }),
options,
)
}
fn new(
beta: &Coefficients,
curvature: HessianCurvatureKind,
firth_active: bool,
options: &WorkingModelPirlsOptions,
) -> Self {
let arrow_latent_bits = options.arrow_schur.as_ref().map(|arrow_cfg| {
arrow_cfg.snapshot_t.as_ref()()
.iter()
.map(|value| value.to_bits())
.collect()
});
Self {
curvature,
firth_active,
beta_bits: beta.as_ref().iter().map(|value| value.to_bits()).collect(),
arrow_latent_bits,
}
}
}
#[derive(Clone, Copy)]
pub(crate) struct IntegratedWorkingInput<'a> {
pub quadctx: &'a crate::quadrature::QuadratureContext,
pub se: ArrayView1<'a, f64>,
pub mixture_link_state: Option<&'a MixtureLinkState>,
pub sas_link_state: Option<&'a SasLinkState>,
}
pub struct WorkingDerivativeBuffersMut<'a> {
c: &'a mut Array1<f64>,
d: &'a mut Array1<f64>,
dmu_deta: &'a mut Array1<f64>,
d2mu_deta2: &'a mut Array1<f64>,
d3mu_deta3: &'a mut Array1<f64>,
}
#[derive(Clone, Copy)]
struct WorkingBernoulliGeometry {
mu: f64,
weight: f64,
z: f64,
c: f64,
d: f64,
}
pub(crate) trait WorkingLikelihood {
fn irls_update(
&self,
y: ArrayView1<f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
integrated: Option<IntegratedWorkingInput<'_>>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) -> Result<(), EstimationError>;
fn loglik_deviance(
&self,
y: ArrayView1<f64>,
mu: &Array1<f64>,
priorweights: ArrayView1<f64>,
) -> Result<f64, EstimationError>;
}
impl WorkingLikelihood for GlmLikelihoodSpec {
fn irls_update(
&self,
y: ArrayView1<f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
integrated: Option<IntegratedWorkingInput<'_>>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) -> Result<(), EstimationError> {
match (&self.spec.response, &self.spec.link, integrated.is_some()) {
(ResponseFamily::Binomial, _, true) => {
let integ = integrated.unwrap();
update_glmvectors_integrated_by_family(
integ.quadctx,
y,
eta,
integ.se,
&self.spec,
priorweights,
mu,
weights,
z,
derivatives,
integ.mixture_link_state,
integ.sas_link_state,
)?;
Ok(())
}
(ResponseFamily::Binomial, link, false) => {
if matches!(link, InverseLink::Mixture(_)) {
crate::bail_invalid_estim!(
"BinomialMixture IRLS update requires explicit mixture link state"
.to_string(),
);
}
update_glmvectors(
y,
eta,
&self.spec.link,
priorweights,
mu,
weights,
z,
derivatives,
)?;
Ok(())
}
(ResponseFamily::Gaussian, _, _) => {
update_glmvectors(
y,
eta,
&InverseLink::Standard(StandardLink::Identity),
priorweights,
mu,
weights,
z,
None,
)?;
if let Some(phi) = self.scale.fixed_phi() {
if !(phi.is_finite() && phi > 0.0) {
crate::bail_invalid_estim!(
"Gaussian fixed dispersion phi must be finite and positive (got {})",
phi
);
}
if phi != 1.0 {
let inv_phi = 1.0 / phi;
weights.mapv_inplace(|w| w * inv_phi);
}
}
Ok(())
}
(ResponseFamily::Poisson, _, _) => {
write_poisson_log_working_state(y, eta, priorweights, mu, weights, z, derivatives);
Ok(())
}
(ResponseFamily::Tweedie { p }, _, _) => {
let p = *p;
write_tweedie_log_working_state(
y,
eta,
priorweights,
p,
fixed_glm_dispersion(self),
mu,
weights,
z,
derivatives,
)?;
Ok(())
}
(ResponseFamily::NegativeBinomial { theta }, _, _) => {
let theta = *theta;
write_negative_binomial_log_working_state(
y,
eta,
priorweights,
theta,
mu,
weights,
z,
derivatives,
)?;
Ok(())
}
(ResponseFamily::Beta { phi }, _, _) => {
let phi = *phi;
write_beta_logit_working_state(
y,
eta,
priorweights,
phi,
mu,
weights,
z,
derivatives,
)?;
Ok(())
}
(ResponseFamily::Gamma, _, _) => {
write_gamma_log_working_state(
y,
eta,
priorweights,
self.gamma_shape().unwrap_or(1.0),
mu,
weights,
z,
derivatives,
);
Ok(())
}
(ResponseFamily::RoystonParmar, _, _) => Err(EstimationError::InvalidInput(
"RoystonParmar is survival-specific and not a GLM IRLS family".to_string(),
)),
}
}
fn loglik_deviance(
&self,
y: ArrayView1<f64>,
mu: &Array1<f64>,
priorweights: ArrayView1<f64>,
) -> Result<f64, EstimationError> {
if matches!(self.spec.response, ResponseFamily::Tweedie { .. }) {
validate_tweedie_responses(&y, &priorweights)?;
}
Ok(calculate_deviance(y, mu, self, priorweights))
}
}
pub struct PirlsWorkspace {
pub wz: Array1<f64>,
pub eta_buf: Array1<f64>,
pub scaled_matrix: Array2<f64>, pub final_aug_matrix: Array2<f64>, pub rhs_full: Array1<f64>, pub working_residual: Array1<f64>,
pub weighted_residual: Array1<f64>,
pub delta_eta: Array1<f64>,
pub vec_buf_p: Array1<f64>,
sparse_penalized_system_cache: Option<SparsePenalizedSystemCache>,
pub factorization_scratch: MemBuffer,
pub perm: Vec<usize>,
pub perm_inv: Vec<usize>,
pub factorization_matrix: Array2<f64>,
pub weighted_xvalues: Vec<f64>,
pub weighted_x_chunk: Array2<f64>,
pub hessian_buf: Array2<f64>,
pub matvec_buf: Array1<f64>,
}
impl PirlsWorkspace {
pub fn new(n: usize, p: usize, idx: usize, idx2: usize) -> Self {
assert!(idx < usize::MAX);
assert!(idx2 < usize::MAX);
PirlsWorkspace {
wz: Array1::zeros(n),
eta_buf: Array1::zeros(n),
scaled_matrix: Array2::zeros((0, 0).f()),
final_aug_matrix: Array2::zeros((0, 0).f()),
rhs_full: Array1::zeros(0),
working_residual: Array1::zeros(n),
weighted_residual: Array1::zeros(n),
delta_eta: Array1::zeros(n),
vec_buf_p: Array1::zeros(p),
sparse_penalized_system_cache: None,
factorization_scratch: {
let par = faer::Par::Seq;
let req = faer::linalg::cholesky::llt::factor::cholesky_in_place_scratch::<f64>(
1,
par,
Spec::new(<LltParams as Auto<f64>>::auto()),
);
MemBuffer::new(req)
},
perm: vec![0; p],
perm_inv: vec![0; p],
factorization_matrix: Array2::zeros((0, 0)),
weighted_xvalues: Vec::new(),
weighted_x_chunk: Array2::zeros((0, 0).f()),
hessian_buf: Array2::zeros((0, 0).f()),
matvec_buf: Array1::zeros(n),
}
}
pub(super) fn add_dense_xtwx_signed(
weights: &Array1<f64>,
weighted_x_scratch: &mut Array2<f64>,
x: &Array2<f64>,
out: &mut Array2<f64>,
) {
*out = crate::solver::estimate::reml::assembly::xt_diag_x_dense_into(
x,
weights,
weighted_x_scratch,
);
}
fn ensure_sparse_penalty_cache(
&mut self,
x: &SparseColMat<usize, f64>,
s_lambda: &Array2<f64>,
) -> Result<(), EstimationError> {
let penalty_pattern = SparsePenaltyPattern::from_dense_upper(s_lambda, 1e-12);
let rebuild = match self.sparse_penalized_system_cache.as_ref() {
Some(cache) => !cache.matches(x, &penalty_pattern),
None => true,
};
if rebuild {
self.sparse_penalized_system_cache =
Some(SparsePenalizedSystemCache::new(x, penalty_pattern)?);
}
Ok(())
}
pub(crate) fn sparse_penalized_system_stats(
&mut self,
x: &SparseColMat<usize, f64>,
s_lambda: &Array2<f64>,
) -> Result<SparsePenalizedSystemStats, EstimationError> {
self.ensure_sparse_penalty_cache(x, s_lambda)?;
Ok(self.sparse_penalized_system_cache.as_ref().unwrap().stats())
}
pub(super) fn assemble_sparse_penalized_hessian(
&mut self,
x: &SparseColMat<usize, f64>,
weights: &Array1<f64>,
s_lambda: &Array2<f64>,
ridge: f64,
precomputed_xtwx: Option<&SparseXtwxPrecomputed>,
) -> Result<SparseColMat<usize, f64>, EstimationError> {
self.ensure_sparse_penalty_cache(x, s_lambda)?;
self.sparse_penalized_system_cache
.as_mut()
.unwrap()
.assemble_upper(x, weights, ridge, precomputed_xtwx)
}
}
#[derive(Clone, Debug)]
pub struct WorkingModelPirlsOptions {
pub max_iterations: usize,
pub convergence_tolerance: f64,
pub adaptive_kkt_tolerance: Option<AdaptiveKktTolerance>,
pub max_step_halving: usize,
pub min_step_size: f64,
pub firth_bias_reduction: bool,
pub coefficient_lower_bounds: Option<Array1<f64>>,
pub linear_constraints: Option<LinearInequalityConstraints>,
pub initial_lm_lambda: Option<f64>,
pub geodesic_acceleration: bool,
pub arrow_schur: Option<ArrowSchurInnerConfig>,
}
#[derive(Clone)]
pub struct ArrowSchurInnerConfig {
pub n_rows: usize,
pub latent_dim: usize,
pub n_beta: usize,
pub build: std::sync::Arc<
dyn Fn(&Array1<f64>) -> Option<crate::solver::arrow_schur::ArrowSchurSystem> + Send + Sync,
>,
pub solver_mode: Option<crate::solver::arrow_schur::ArrowSolverMode>,
pub streaming_chunk_size: Option<usize>,
pub trust_region_radius: f64,
pub block_offsets: Option<Arc<[std::ops::Range<usize>]>>,
pub apply_delta_t: std::sync::Arc<dyn Fn(&Array1<f64>) + Send + Sync>,
pub snapshot_t: std::sync::Arc<dyn Fn() -> Array1<f64> + Send + Sync>,
pub restore_t: std::sync::Arc<dyn Fn(&Array1<f64>) + Send + Sync>,
}
impl std::fmt::Debug for ArrowSchurInnerConfig {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("ArrowSchurInnerConfig")
.field("n_rows", &self.n_rows)
.field("latent_dim", &self.latent_dim)
.field("n_beta", &self.n_beta)
.field("solver_mode", &self.solver_mode)
.field("streaming_chunk_size", &self.streaming_chunk_size)
.field("trust_region_radius", &self.trust_region_radius)
.field(
"block_offsets",
&self.block_offsets.as_ref().map(|o| o.len()),
)
.finish_non_exhaustive()
}
}
fn restore_arrow_latent_if_needed(
options: &WorkingModelPirlsOptions,
snapshot: Option<Array1<f64>>,
) {
if let (Some(arrow_cfg), Some(snapshot)) = (options.arrow_schur.as_ref(), snapshot) {
arrow_cfg.restore_t.as_ref()(&snapshot);
}
}
pub(super) fn restore_pending_arrow_latent_if_needed(
options: &WorkingModelPirlsOptions,
pending_snapshot: &mut Option<Array1<f64>>,
) {
restore_arrow_latent_if_needed(options, pending_snapshot.take());
}
pub(super) fn commit_pending_arrow_latent(pending_snapshot: &mut Option<Array1<f64>>) {
drop(pending_snapshot.take());
}
pub(super) const FIXED_STABILIZATION_RIDGE: f64 = 1e-8;
pub(super) struct GamWorkingModel<'a> {
x_original: DesignMatrix,
coordinate_design: WorkingCoordinateDesign,
offset: Array1<f64>,
y: ArrayView1<'a, f64>,
priorweights: ArrayView1<'a, f64>,
penalty: PirlsPenalty,
workspace: PirlsWorkspace,
likelihood: GlmLikelihoodSpec,
link_kind: InverseLink,
firth_bias_reduction: bool,
lastmu: Array1<f64>,
lastweights: Array1<f64>,
lastz: Array1<f64>,
last_c: Array1<f64>,
last_d: Array1<f64>,
lasthessian_weights: Array1<f64>,
lasthessian_c: Array1<f64>,
lasthessian_d: Array1<f64>,
lasthessian_curvature: HessianCurvatureKind,
last_dmu_deta: Array1<f64>,
last_d2mu_deta2: Array1<f64>,
last_d3mu_deta3: Array1<f64>,
last_penalty_term: f64,
x_original_csr: Option<SparseRowMat<usize, f64>>,
covariate_se: Option<Array1<f64>>,
quadctx: crate::quadrature::QuadratureContext,
}
pub(super) struct GamModelFinalState {
likelihood: GlmLikelihoodSpec,
coordinate_frame: PirlsCoordinateFrame,
finalmu: Array1<f64>,
finalweights: Array1<f64>,
scoreweights: Array1<f64>,
finalz: Array1<f64>,
final_c: Array1<f64>,
final_d: Array1<f64>,
final_dmu_deta: Array1<f64>,
final_d2mu_deta2: Array1<f64>,
final_d3mu_deta3: Array1<f64>,
penalty_term: f64,
}
impl<'a> GamWorkingModel<'a> {
fn new(
x_transformed: Option<DesignMatrix>,
x_original: DesignMatrix,
coordinate_frame: PirlsCoordinateFrame,
offset: ArrayView1<f64>,
y: ArrayView1<'a, f64>,
priorweights: ArrayView1<'a, f64>,
penalty: PirlsPenalty,
workspace: PirlsWorkspace,
likelihood: GlmLikelihoodSpec,
link_kind: InverseLink,
firth_bias_reduction: bool,
transform: Option<WorkingReparamTransform>,
quadctx: crate::quadrature::QuadratureContext,
) -> Self {
let coordinate_design = match coordinate_frame {
PirlsCoordinateFrame::OriginalSparseNative => {
WorkingCoordinateDesign::OriginalSparseNative
}
PirlsCoordinateFrame::TransformedQs => {
if let Some(x_transformed) = x_transformed {
WorkingCoordinateDesign::TransformedExplicit {
x_csr: x_transformed.to_csr_cache(),
x_transformed,
}
} else {
WorkingCoordinateDesign::TransformedImplicit {
transform: transform.expect(
"TransformedQs PIRLS coordinate frame requires either x_transformed or qs",
),
}
}
}
};
let x_original_csr = x_original.to_csr_cache();
let n = match &coordinate_design {
WorkingCoordinateDesign::OriginalSparseNative => x_original.nrows(),
WorkingCoordinateDesign::TransformedExplicit { x_transformed, .. } => {
x_transformed.nrows()
}
WorkingCoordinateDesign::TransformedImplicit { .. } => x_original.nrows(),
};
GamWorkingModel {
x_original,
coordinate_design,
offset: offset.to_owned(),
y,
priorweights,
penalty,
workspace,
likelihood,
link_kind,
firth_bias_reduction,
lastmu: Array1::zeros(n),
lastweights: Array1::zeros(n),
lastz: Array1::zeros(n),
last_c: Array1::zeros(n),
last_d: Array1::zeros(n),
lasthessian_weights: Array1::zeros(n),
lasthessian_c: Array1::zeros(n),
lasthessian_d: Array1::zeros(n),
lasthessian_curvature: HessianCurvatureKind::Fisher,
last_dmu_deta: Array1::zeros(n),
last_d2mu_deta2: Array1::zeros(n),
last_d3mu_deta3: Array1::zeros(n),
last_penalty_term: 0.0,
x_original_csr,
covariate_se: None,
quadctx,
}
}
fn with_covariate_se(mut self, se: Array1<f64>) -> Self {
self.covariate_se = Some(se);
self
}
fn into_final_state(self) -> GamModelFinalState {
let GamWorkingModel {
coordinate_design,
lastmu,
lastweights,
lastz,
last_c: _,
last_d: _,
lasthessian_weights,
lasthessian_c,
lasthessian_d,
last_dmu_deta,
last_d2mu_deta2,
last_d3mu_deta3,
last_penalty_term,
..
} = self;
let coordinate_frame = match coordinate_design {
WorkingCoordinateDesign::OriginalSparseNative => {
PirlsCoordinateFrame::OriginalSparseNative
}
WorkingCoordinateDesign::TransformedExplicit { .. } => {
PirlsCoordinateFrame::TransformedQs
}
WorkingCoordinateDesign::TransformedImplicit { .. } => {
PirlsCoordinateFrame::TransformedQs
}
};
GamModelFinalState {
likelihood: self.likelihood.clone(),
coordinate_frame,
finalmu: lastmu,
finalweights: lasthessian_weights,
scoreweights: lastweights,
finalz: lastz,
final_c: lasthessian_c,
final_d: lasthessian_d,
final_dmu_deta: last_dmu_deta,
final_d2mu_deta2: last_d2mu_deta2,
final_d3mu_deta3: last_d3mu_deta3,
penalty_term: last_penalty_term,
}
}
fn transformed_matvec_into(&self, beta: &Coefficients, out: &mut Array1<f64>) {
self.transformed_matvec_array_into(beta.as_ref(), out);
}
fn transformed_matvec_array_into(&self, beta: &Array1<f64>, out: &mut Array1<f64>) {
match &self.coordinate_design {
WorkingCoordinateDesign::TransformedExplicit { x_transformed, .. } => {
if let Some(dense) = x_transformed.as_dense() {
fast_av_into(dense, beta, out);
return;
}
out.assign(&x_transformed.matrixvectormultiply(beta));
}
WorkingCoordinateDesign::TransformedImplicit { transform } => {
let beta_orig = transform.apply(beta);
if let Some(dense) = self.x_original.as_dense() {
fast_av_into(dense, &beta_orig, out);
} else {
out.assign(&self.x_original.apply(&beta_orig));
}
}
WorkingCoordinateDesign::OriginalSparseNative => {
out.assign(&self.x_original.matrixvectormultiply(beta));
}
}
}
fn transformed_transpose_matvec(&self, vec: &Array1<f64>) -> Array1<f64> {
match &self.coordinate_design {
WorkingCoordinateDesign::OriginalSparseNative => {
self.x_original.transpose_vector_multiply(vec)
}
WorkingCoordinateDesign::TransformedExplicit { x_transformed, .. } => {
x_transformed.transpose_vector_multiply(vec)
}
WorkingCoordinateDesign::TransformedImplicit { transform } => {
let xtv = self.x_original.transpose_vector_multiply(vec);
transform.apply_transpose(&xtv)
}
}
}
fn compute_xtwx_blas(
workspace: &mut PirlsWorkspace,
design: &DesignMatrix,
weights: &Array1<f64>,
) -> Result<Array2<f64>, EstimationError> {
match design {
DesignMatrix::Dense(x) if x.is_materialized_dense() => {
let p = x.ncols();
let x_dense = x.to_dense_arc();
if workspace.hessian_buf.nrows() != p || workspace.hessian_buf.ncols() != p {
workspace.hessian_buf = Array2::zeros((p, p).f());
} else {
workspace.hessian_buf.fill(0.0);
}
if crate::solver::gpu::cuda_selected() {
return crate::solver::gpu::pirls_gpu::weighted_crossprod_gpu(
x_dense.view(),
weights.view(),
)
.map_err(EstimationError::InvalidInput);
}
crate::gpu::log_backend_inventory_once();
let gpu_decision = crate::gpu::decide(
crate::gpu::GpuKernel::DenseXtWX,
crate::gpu::GpuEligibility::BackendNotCompiled,
);
gpu_decision
.require_supported()
.map_err(EstimationError::InvalidInput)?;
gpu_decision.log();
if weights.iter().any(|&w| w < 0.0) {
PirlsWorkspace::add_dense_xtwx_signed(
weights,
&mut workspace.weighted_x_chunk,
x_dense.as_ref(),
&mut workspace.hessian_buf,
);
} else {
PirlsWorkspace::add_dense_xtwx_signed(
weights,
&mut workspace.weighted_x_chunk,
x_dense.as_ref(),
&mut workspace.hessian_buf,
);
}
Ok(std::mem::take(&mut workspace.hessian_buf))
}
_ => crate::matrix::xt_diag_x_signed(
design,
crate::matrix::SignedWeightsView::from_array(weights),
)
.map(|h| h.to_dense())
.map_err(EstimationError::InvalidInput),
}
}
fn penalized_hessian(&mut self, weights: &Array1<f64>) -> Result<Array2<f64>, EstimationError> {
match &self.coordinate_design {
WorkingCoordinateDesign::TransformedExplicit { x_transformed, .. } => {
let mut h = Self::compute_xtwx_blas(&mut self.workspace, x_transformed, weights)?;
self.penalty.add_to_hessian(&mut h);
Ok(h)
}
WorkingCoordinateDesign::TransformedImplicit { transform } => {
let xtwx = Self::compute_xtwx_blas(&mut self.workspace, &self.x_original, weights)?;
let mut h = transform.conjugate_matrix(&xtwx);
self.penalty.add_to_hessian(&mut h);
Ok(h)
}
WorkingCoordinateDesign::OriginalSparseNative => {
let mut h =
Self::compute_xtwx_blas(&mut self.workspace, &self.x_original, weights)?;
self.penalty.add_to_hessian(&mut h);
Ok(h)
}
}
}
fn supports_observed_hessian_curvature(&self) -> bool {
supports_observed_hessian_curvature_for_likelihood(&self.likelihood, &self.link_kind)
}
fn update_hessian_curvature_arrays(
&mut self,
requested: HessianCurvatureKind,
) -> Result<HessianCurvatureKind, EstimationError> {
if requested == HessianCurvatureKind::Fisher || !self.supports_observed_hessian_curvature()
{
self.lasthessian_weights.assign(&self.lastweights);
self.lasthessian_c.assign(&self.last_c);
self.lasthessian_d.assign(&self.last_d);
return Ok(HessianCurvatureKind::Fisher);
}
compute_observed_hessian_curvature_arrays_into(
&self.likelihood,
&self.link_kind,
&self.workspace.eta_buf,
self.y,
&self.lastweights,
self.priorweights,
&mut self.lasthessian_weights,
&mut self.lasthessian_c,
&mut self.lasthessian_d,
)?;
Ok(HessianCurvatureKind::Observed)
}
fn sparse_penalized_hessian(
&mut self,
weights: &Array1<f64>,
ridge: f64,
) -> Result<SparseColMat<usize, f64>, EstimationError> {
let x_sparse = self.x_original.as_sparse().ok_or_else(|| {
EstimationError::InvalidInput(
"sparse-native PIRLS requires a sparse original design".to_string(),
)
})?;
let PirlsPenalty::Dense { s_transformed, .. } = &self.penalty else {
crate::bail_invalid_estim!(
"sparse-native PIRLS requires a dense transformed penalty matrix"
);
};
self.workspace.assemble_sparse_penalized_hessian(
x_sparse,
weights,
s_transformed,
ridge,
None,
)
}
fn screen_candidate_from_direction(
&mut self,
beta: &Coefficients,
direction: &Array1<f64>,
current_eta: &LinearPredictor,
) -> Result<CandidateScreen, EstimationError> {
let n = self.offset.len();
if self.workspace.eta_buf.len() != n {
self.workspace.eta_buf = Array1::zeros(n);
}
if self.workspace.delta_eta.len() != n {
self.workspace.delta_eta = Array1::zeros(n);
}
let mut delta_eta = std::mem::take(&mut self.workspace.delta_eta);
self.transformed_matvec_array_into(direction, &mut delta_eta);
Zip::from(&mut self.workspace.eta_buf)
.and(current_eta.as_ref())
.and(&delta_eta)
.par_for_each(|eta, &base, &d| *eta = base + d);
self.workspace.delta_eta = delta_eta;
if self.likelihood.scale.gamma_shape_is_estimated() {
let shape =
estimate_gamma_shape_from_eta(self.y, &self.workspace.eta_buf, self.priorweights);
self.likelihood = self.likelihood.clone().with_gamma_shape(shape);
}
let integrated = self.covariate_se.as_ref().map(|se| IntegratedWorkingInput {
quadctx: &self.quadctx,
se: se.view(),
mixture_link_state: self.link_kind.mixture_state(),
sas_link_state: self.link_kind.sas_state(),
});
match &self.link_kind {
InverseLink::Mixture(_)
| InverseLink::LatentCLogLog(_)
| InverseLink::Sas(_)
| InverseLink::BetaLogistic(_) => {
if let Some(integ) = integrated {
update_glmvectors_integrated_for_link(
integ.quadctx,
self.y,
&self.workspace.eta_buf,
integ.se,
&self.link_kind,
self.priorweights,
&mut self.lastmu,
&mut self.lastweights,
&mut self.lastz,
None,
)?;
} else {
update_glmvectors(
self.y,
&self.workspace.eta_buf,
&self.link_kind,
self.priorweights,
&mut self.lastmu,
&mut self.lastweights,
&mut self.lastz,
None,
)?;
}
}
InverseLink::Standard(_) => {
self.likelihood.irls_update(
self.y,
&self.workspace.eta_buf,
self.priorweights,
&mut self.lastmu,
&mut self.lastweights,
&mut self.lastz,
integrated,
None,
)?;
}
}
let deviance = self
.likelihood
.loglik_deviance(self.y, &self.lastmu, self.priorweights)?;
let penalty_term = self.penalty.shifted_quadratic(beta.as_ref());
let penalized_objective = deviance + penalty_term;
let arithmetic_finite = penalized_objective.is_finite()
&& self.workspace.eta_buf.iter().all(|v| v.is_finite())
&& self.lastmu.iter().all(|v| v.is_finite())
&& self.lastweights.iter().all(|v| v.is_finite());
Ok(CandidateScreen {
penalized_objective,
deviance,
penalty_term,
arithmetic_finite,
})
}
}
impl<'a> WorkingModel for GamWorkingModel<'a> {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
requested_curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
let n = self.offset.len();
if self.workspace.eta_buf.len() != n {
self.workspace.eta_buf = Array1::zeros(n);
}
if self.workspace.matvec_buf.len() != n {
self.workspace.matvec_buf = Array1::zeros(n);
}
let mut matvec_tmp = std::mem::take(&mut self.workspace.matvec_buf);
self.transformed_matvec_into(beta, &mut matvec_tmp);
self.workspace.eta_buf.assign(&self.offset);
self.workspace.eta_buf += &matvec_tmp;
self.workspace.matvec_buf = matvec_tmp;
if self.likelihood.scale.gamma_shape_is_estimated() {
let shape =
estimate_gamma_shape_from_eta(self.y, &self.workspace.eta_buf, self.priorweights);
self.likelihood = self.likelihood.clone().with_gamma_shape(shape);
}
let integrated = self.covariate_se.as_ref().map(|se| IntegratedWorkingInput {
quadctx: &self.quadctx,
se: se.view(),
mixture_link_state: self.link_kind.mixture_state(),
sas_link_state: self.link_kind.sas_state(),
});
match &self.link_kind {
InverseLink::Mixture(_) => {
if let Some(integ) = integrated {
update_glmvectors_integrated_for_link(
integ.quadctx,
self.y,
&self.workspace.eta_buf,
integ.se,
&self.link_kind,
self.priorweights,
&mut self.lastmu,
&mut self.lastweights,
&mut self.lastz,
Some(WorkingDerivativeBuffersMut {
c: &mut self.last_c,
d: &mut self.last_d,
dmu_deta: &mut self.last_dmu_deta,
d2mu_deta2: &mut self.last_d2mu_deta2,
d3mu_deta3: &mut self.last_d3mu_deta3,
}),
)?;
} else {
update_glmvectors(
self.y,
&self.workspace.eta_buf,
&self.link_kind,
self.priorweights,
&mut self.lastmu,
&mut self.lastweights,
&mut self.lastz,
Some(WorkingDerivativeBuffersMut {
c: &mut self.last_c,
d: &mut self.last_d,
dmu_deta: &mut self.last_dmu_deta,
d2mu_deta2: &mut self.last_d2mu_deta2,
d3mu_deta3: &mut self.last_d3mu_deta3,
}),
)?;
}
}
InverseLink::LatentCLogLog(_) | InverseLink::Sas(_) | InverseLink::BetaLogistic(_) => {
if let Some(integ) = integrated {
update_glmvectors_integrated_for_link(
integ.quadctx,
self.y,
&self.workspace.eta_buf,
integ.se,
&self.link_kind,
self.priorweights,
&mut self.lastmu,
&mut self.lastweights,
&mut self.lastz,
Some(WorkingDerivativeBuffersMut {
c: &mut self.last_c,
d: &mut self.last_d,
dmu_deta: &mut self.last_dmu_deta,
d2mu_deta2: &mut self.last_d2mu_deta2,
d3mu_deta3: &mut self.last_d3mu_deta3,
}),
)?;
} else {
update_glmvectors(
self.y,
&self.workspace.eta_buf,
&self.link_kind,
self.priorweights,
&mut self.lastmu,
&mut self.lastweights,
&mut self.lastz,
Some(WorkingDerivativeBuffersMut {
c: &mut self.last_c,
d: &mut self.last_d,
dmu_deta: &mut self.last_dmu_deta,
d2mu_deta2: &mut self.last_d2mu_deta2,
d3mu_deta3: &mut self.last_d3mu_deta3,
}),
)?;
}
}
InverseLink::Standard(_) => {
self.likelihood.irls_update(
self.y,
&self.workspace.eta_buf,
self.priorweights,
&mut self.lastmu,
&mut self.lastweights,
&mut self.lastz,
integrated,
Some(WorkingDerivativeBuffersMut {
c: &mut self.last_c,
d: &mut self.last_d,
dmu_deta: &mut self.last_dmu_deta,
d2mu_deta2: &mut self.last_d2mu_deta2,
d3mu_deta3: &mut self.last_d3mu_deta3,
}),
)?;
}
}
let mut firth = FirthDiagnostics::Inactive;
if self.firth_bias_reduction {
let (hat_diag, jeffreys_logdet) = match &self.coordinate_design {
WorkingCoordinateDesign::TransformedExplicit {
x_transformed,
x_csr,
} => {
if x_transformed.as_sparse().is_some() {
let csr = x_csr.as_ref().ok_or_else(|| {
EstimationError::InvalidInput(
"missing CSR cache for sparse transformed design".to_string(),
)
})?;
compute_jeffreys_pirls_diagnostics_sparse(
csr,
self.workspace.eta_buf.view(),
self.priorweights,
)?
} else {
let x_dense_cow = x_transformed.to_dense_cow();
compute_jeffreys_pirls_diagnostics(
x_dense_cow.view(),
self.workspace.eta_buf.view(),
self.priorweights,
)?
}
}
WorkingCoordinateDesign::TransformedImplicit { transform } => {
let x_t_dense =
fast_ab(&self.x_original.to_dense(), &transform.materialize_dense());
compute_jeffreys_pirls_diagnostics(
x_t_dense.view(),
self.workspace.eta_buf.view(),
self.priorweights,
)?
}
WorkingCoordinateDesign::OriginalSparseNative => {
if self.x_original.as_sparse().is_some() {
let csr = self.x_original_csr.as_ref().ok_or_else(|| {
EstimationError::InvalidInput(
"missing CSR cache for sparse original design".to_string(),
)
})?;
compute_jeffreys_pirls_diagnostics_sparse(
csr,
self.workspace.eta_buf.view(),
self.priorweights,
)?
} else {
let x_dense = self
.x_original
.try_to_dense_arc(
"Firth diagnostics require dense access to the original design",
)
.map_err(EstimationError::InvalidInput)?;
compute_jeffreys_pirls_diagnostics(
x_dense.view(),
self.workspace.eta_buf.view(),
self.priorweights,
)?
}
}
};
firth = FirthDiagnostics::Active {
jeffreys_logdet,
hat_diag: hat_diag.clone(),
};
ndarray::Zip::from(&mut self.lastz)
.and(&hat_diag)
.and(&self.lastweights)
.and(&self.lastmu)
.par_for_each(|zi, &hii, &wi, &mui| {
if wi > 0.0 {
*zi += hii * (0.5 - mui) / wi;
}
});
}
let z = &self.lastz;
ndarray::Zip::from(&mut self.workspace.weighted_residual)
.and(&mut self.workspace.working_residual)
.and(&self.workspace.eta_buf)
.and(z)
.and(&self.lastweights)
.par_for_each(|wr, r, &eta, &zi, &wi| {
let residual = eta - zi;
*r = residual;
*wr = residual * wi;
});
let mut gradient = self.transformed_transpose_matvec(&self.workspace.weighted_residual);
let score_norm = array1_l2_norm(&gradient);
let s_beta = self.penalty.shifted_gradient(beta.as_ref());
let s_beta_norm = array1_l2_norm(&s_beta);
gradient += &s_beta;
let hessian_curvature = self.update_hessian_curvature_arrays(requested_curvature)?;
self.lasthessian_curvature = hessian_curvature;
if self.workspace.matvec_buf.len() != n {
self.workspace.matvec_buf = Array1::zeros(n);
}
solver_hessian_weights_into(
&self.lasthessian_weights,
&self.lastweights,
&mut self.workspace.matvec_buf,
);
let solver_weights = std::mem::take(&mut self.workspace.matvec_buf);
let (penalized_hessian, sparsehessian, ridge_used) = if matches!(
self.coordinate_design,
WorkingCoordinateDesign::OriginalSparseNative
) {
let (h_sparse, _factor, ridge_used) =
ensure_sparse_positive_definitewithridge(|ridge| {
self.sparse_penalized_hessian(&solver_weights, ridge)
})?;
(Array2::zeros((0, 0)), Some(h_sparse), ridge_used)
} else {
let mut penalized_hessian = self.penalized_hessian(&solver_weights)?;
assert_symmetric_tol(&penalized_hessian, "PIRLS penalized Hessian", 1e-8);
let ridge_used = ensure_positive_definitewithridge(
&mut penalized_hessian,
"PIRLS penalized Hessian",
)?;
(penalized_hessian, None, ridge_used)
};
self.workspace.matvec_buf = solver_weights;
let deviance = self
.likelihood
.loglik_deviance(self.y, &self.lastmu, self.priorweights)?;
let log_likelihood = calculate_loglikelihood_omitting_constants(
self.y,
&self.lastmu,
&self.likelihood,
self.priorweights,
);
let mut penalty_term = self.penalty.shifted_quadratic(beta.as_ref());
let mut ridge_grad_norm = 0.0;
if ridge_used > 0.0 {
let ridge_penalty = ridge_used * beta.as_ref().dot(beta.as_ref());
penalty_term += ridge_penalty;
gradient.zip_mut_with(beta.as_ref(), |g, &b| *g += ridge_used * b);
ridge_grad_norm = ridge_used * array1_l2_norm(beta.as_ref());
}
self.last_penalty_term = penalty_term;
let gradient_natural_scale = score_norm + s_beta_norm + ridge_grad_norm;
Ok(WorkingState {
eta: LinearPredictor::new(std::mem::replace(
&mut self.workspace.eta_buf,
Array1::zeros(0),
)),
gradient,
hessian: match sparsehessian {
Some(h_sparse) => crate::linalg::matrix::SymmetricMatrix::Sparse(h_sparse),
None => crate::linalg::matrix::SymmetricMatrix::Dense(penalized_hessian),
},
log_likelihood,
deviance,
penalty_term,
firth,
ridge_used,
hessian_curvature,
gradient_natural_scale,
})
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
if !self.firth_bias_reduction {
return self.update_with_curvature(beta, curvature);
}
let firth_enabled = self.firth_bias_reduction;
self.firth_bias_reduction = false;
let result = self.update_with_curvature(beta, curvature);
self.firth_bias_reduction = firth_enabled;
result
}
fn screen_candidate(
&mut self,
beta: &Coefficients,
direction: &Array1<f64>,
current_eta: &LinearPredictor,
curvature: HessianCurvatureKind,
) -> Result<CandidateEvaluation, EstimationError> {
if self.firth_bias_reduction {
return self
.update_candidate(beta, curvature)
.map(CandidateEvaluation::Full);
}
self.screen_candidate_from_direction(beta, direction, current_eta)
.map(CandidateEvaluation::Screen)
}
fn supports_observed_information_curvature(&self) -> bool {
self.supports_observed_hessian_curvature()
}
}
const DENSE_OUTER_MAX_P: usize = 1024;
const DENSE_OUTER_PARALLEL_FLOP_THRESHOLD: u64 = 100_000;
enum XtWxBackend {
Dense(DenseOuterState),
Sparse(SparseSpGemmState),
}
struct DenseOuterState {
xtwx_dense: Array2<f64>,
thread_buffers: Vec<Array2<f64>>,
}
struct SparseSpGemmState {
wxvalues: Vec<f64>,
wx_tvalues: Vec<f64>,
sqrt_weights: Vec<f64>,
info: SparseMatMulInfo,
scratch: MemBuffer,
par: Par,
}
pub(crate) struct SparseXtWxCache {
xtwx_symbolic: SymbolicSparseColMat<usize>,
xtwxvalues: Vec<f64>,
nrows: usize,
ncols: usize,
nnz: usize,
x_col_ptr: Vec<usize>,
xrow_idx: Vec<usize>,
x_t_csc: SparseColMat<usize, f64>,
backend: XtWxBackend,
}
impl SparseXtWxCache {
fn new(x: &SparseColMat<usize, f64>) -> Result<Self, EstimationError> {
let x_t_csc =
x.as_ref().transpose().to_col_major().map_err(|_| {
EstimationError::InvalidInput("failed to transpose to CSC".to_string())
})?;
let (xtwx_symbolic, info) = sparse_sparse_matmul_symbolic(x_t_csc.symbolic(), x.symbolic())
.map_err(|_| {
EstimationError::InvalidInput("failed to build symbolic XtWX cache".to_string())
})?;
let xtwxvalues = vec![0.0; xtwx_symbolic.row_idx().len()];
let backend = if x.ncols() <= DENSE_OUTER_MAX_P {
XtWxBackend::Dense(DenseOuterState {
xtwx_dense: Array2::<f64>::zeros((x.ncols(), x.ncols())),
thread_buffers: Vec::new(),
})
} else {
let par = get_global_parallelism();
let scratch = MemBuffer::new(sparse_sparse_matmul_numeric_scratch::<usize, f64>(
xtwx_symbolic.as_ref(),
par,
));
XtWxBackend::Sparse(SparseSpGemmState {
wxvalues: vec![0.0; x.val().len()],
wx_tvalues: vec![0.0; x_t_csc.val().len()],
sqrt_weights: vec![0.0; x.nrows()],
info,
scratch,
par,
})
};
Ok(Self {
xtwx_symbolic,
xtwxvalues,
nrows: x.nrows(),
ncols: x.ncols(),
nnz: x.val().len(),
x_col_ptr: x.symbolic().col_ptr().to_vec(),
xrow_idx: x.symbolic().row_idx().to_vec(),
x_t_csc,
backend,
})
}
fn matches(&self, x: &SparseColMat<usize, f64>) -> bool {
if self.nrows != x.nrows() || self.ncols != x.ncols() || self.nnz != x.val().len() {
return false;
}
let sym = x.symbolic();
self.x_col_ptr.as_slice() == sym.col_ptr() && self.xrow_idx.as_slice() == sym.row_idx()
}
fn compute_numeric(
&mut self,
x: &SparseColMat<usize, f64>,
weights: &Array1<f64>,
) -> Result<(), EstimationError> {
if weights.len() != self.nrows {
crate::bail_invalid_estim!(
"weights length {} does not match design rows {}",
weights.len(),
self.nrows
);
}
match &mut self.backend {
XtWxBackend::Dense(state) => {
state.compute(self.x_t_csc.as_ref(), weights, self.nrows, self.ncols);
let col_ptr = self.xtwx_symbolic.col_ptr();
let row_idx = self.xtwx_symbolic.row_idx();
let dense = &state.xtwx_dense;
for col in 0..self.ncols {
let start = col_ptr[col];
let end = col_ptr[col + 1];
for idx in start..end {
let row = row_idx[idx];
if row <= col {
self.xtwxvalues[idx] = dense[[row, col]];
}
}
}
}
XtWxBackend::Sparse(state) => state.compute(
x,
self.x_t_csc.as_ref(),
weights,
self.ncols,
self.xtwx_symbolic.as_ref(),
&mut self.xtwxvalues,
),
}
Ok(())
}
}
impl DenseOuterState {
fn compute(
&mut self,
x_t: SparseColMatRef<'_, usize, f64>,
weights: &Array1<f64>,
n: usize,
p: usize,
) {
assert_eq!(self.xtwx_dense.dim(), (p, p));
self.xtwx_dense.fill(0.0);
if n == 0 || p == 0 {
return;
}
let xtwx_start = std::time::Instant::now();
let nnz_total = x_t.symbolic().row_idx().len() as u64;
let work = nnz_total
.saturating_mul(nnz_total)
.checked_div(n as u64)
.unwrap_or(u64::MAX);
let n_threads = rayon::current_num_threads();
let parallelize = n_threads > 1 && work >= DENSE_OUTER_PARALLEL_FLOP_THRESHOLD;
if !parallelize {
accumulate_outer_upper(&mut self.xtwx_dense, x_t, weights, 0..n);
log::info!(
"[STAGE] PIRLS dense XᵀWX assembly (serial) n={} p={} flops~{} elapsed={:.3}s",
n,
p,
(n as u64).saturating_mul((p as u64).saturating_mul(p as u64)),
xtwx_start.elapsed().as_secs_f64(),
);
return;
}
if self.thread_buffers.len() != n_threads {
self.thread_buffers
.resize_with(n_threads, || Array2::<f64>::zeros((p, p)));
}
let chunk = n.div_ceil(n_threads);
self.thread_buffers
.par_iter_mut()
.enumerate()
.for_each(|(t, buf)| {
buf.fill(0.0);
let start = t * chunk;
let end = (start + chunk).min(n);
if start < end {
accumulate_outer_upper(buf, x_t, weights, start..end);
}
});
for buf in &self.thread_buffers {
self.xtwx_dense += buf;
}
log::info!(
"[STAGE] PIRLS dense XᵀWX assembly (parallel, threads={}) n={} p={} flops~{} elapsed={:.3}s",
rayon::current_num_threads(),
n,
p,
(n as u64).saturating_mul((p as u64).saturating_mul(p as u64)),
xtwx_start.elapsed().as_secs_f64(),
);
}
}
impl SparseSpGemmState {
fn compute(
&mut self,
x: &SparseColMat<usize, f64>,
x_t: SparseColMatRef<'_, usize, f64>,
weights: &Array1<f64>,
p: usize,
xtwx_symbolic: SymbolicSparseColMatRef<'_, usize>,
xtwxvalues: &mut [f64],
) {
let n = x_t.ncols();
assert_eq!(weights.len(), n);
assert_eq!(self.sqrt_weights.len(), n);
assert!(
weights.iter().all(|&w| w.is_finite() && w >= 0.0),
"SparseSpGemmState::compute requires finite nonnegative PIRLS weights"
);
let sqrt_w = self.sqrt_weights.as_mut_slice();
for (dst, &w) in sqrt_w.iter_mut().zip(weights.iter()) {
*dst = w.sqrt();
}
let sqrt_w: &[f64] = sqrt_w;
let x_ref = x.as_ref();
for col in 0..p {
let rows = x_ref.row_idx_of_col_raw(col);
let xvals = x_ref.val_of_col(col);
let range = x_ref.col_range(col);
let dst = &mut self.wxvalues[range];
for ((d, &s), row) in dst.iter_mut().zip(xvals.iter()).zip(rows.iter()) {
*d = s * sqrt_w[row.unbound()];
}
}
for col in 0..n {
let w = sqrt_w[col];
let xvals = x_t.val_of_col(col);
let range = x_t.col_range(col);
let dst = &mut self.wx_tvalues[range];
for (d, &s) in dst.iter_mut().zip(xvals.iter()) {
*d = s * w;
}
}
let wx_ref = SparseColMatRef::new(x.symbolic(), &self.wxvalues[..]);
let wx_t_ref = SparseColMatRef::new(x_t.symbolic(), &self.wx_tvalues[..]);
let stack = MemStack::new(&mut self.scratch);
let xtwxmut = SparseColMatMut::new(xtwx_symbolic, xtwxvalues);
sparse_sparse_matmul_numeric(
xtwxmut,
Accum::Replace,
wx_t_ref,
wx_ref,
1.0,
&self.info,
self.par,
stack,
);
}
}
#[inline]
fn accumulate_outer_upper(
acc: &mut Array2<f64>,
x_t: SparseColMatRef<'_, usize, f64>,
weights: &Array1<f64>,
rows: std::ops::Range<usize>,
) {
assert_eq!(acc.nrows(), acc.ncols());
let p = acc.ncols();
let acc_data = acc
.as_slice_mut()
.expect("dense XᵀWX accumulator is row-major and contiguous");
for i in rows {
let w_i = weights[i].max(0.0);
if w_i == 0.0 {
continue;
}
let cols = x_t.row_idx_of_col_raw(i);
let vals = x_t.val_of_col(i);
let nnz_i = cols.len();
for jj in 0..nnz_i {
let j = cols[jj].unbound();
let wvj = w_i * vals[jj];
let row = &mut acc_data[j * p..j * p + p];
for kk in jj..nnz_i {
let k = cols[kk].unbound();
row[k] += wvj * vals[kk];
}
}
}
}
pub(super) fn compute_jeffreys_pirls_diagnostics_sparse(
x_design_csr: &SparseRowMat<usize, f64>,
eta: ArrayView1<f64>,
observation_weights: ArrayView1<f64>,
) -> Result<(Array1<f64>, f64), EstimationError> {
let n = x_design_csr.nrows();
let p = x_design_csr.ncols();
let mut x_dense = Array2::<f64>::zeros((n, p));
let xview = x_design_csr.as_ref();
for i in 0..n {
let vals = xview.val_of_row(i);
let cols = xview.col_idx_of_row_raw(i);
if cols.len() != vals.len() {
crate::bail_invalid_estim!(
"sparse row structure mismatch: column/value lengths differ"
);
}
for (idx, &col) in cols.iter().enumerate() {
x_dense[[i, col.unbound()]] = vals[idx];
}
}
compute_jeffreys_pirls_diagnostics(x_dense.view(), eta, observation_weights)
}
pub(super) fn compute_jeffreys_pirls_diagnostics(
x_design: ArrayView2<f64>,
eta: ArrayView1<f64>,
observation_weights: ArrayView1<f64>,
) -> Result<(Array1<f64>, f64), EstimationError> {
let op = FirthDenseOperator::build_with_observation_weights(
&x_design.to_owned(),
&eta.to_owned(),
observation_weights,
)?;
Ok((op.pirls_hat_diag(), op.jeffreys_logdet()))
}
fn ensure_positive_definitewithridge(
hess: &mut Array2<f64>,
label: &str,
) -> Result<f64, EstimationError> {
let ridge = if FIXED_STABILIZATION_RIDGE > 0.0 {
FIXED_STABILIZATION_RIDGE
} else {
0.0
};
if hess.cholesky(Side::Lower).is_ok() {
return Ok(0.0);
}
if ridge > 0.0 {
for i in 0..hess.nrows() {
hess[[i, i]] += ridge;
}
if hess.cholesky(Side::Lower).is_ok() {
log::debug!("{} stabilized with fixed ridge {:.1e}.", label, ridge);
return Ok(ridge);
}
}
if let Ok((evals, _)) = hess.eigh(Side::Lower) {
let min_eig = evals.iter().fold(f64::INFINITY, |a, &b| a.min(b));
return Err(EstimationError::HessianNotPositiveDefinite {
min_eigenvalue: min_eig,
});
}
Err(EstimationError::HessianNotPositiveDefinite {
min_eigenvalue: f64::NEG_INFINITY,
})
}
pub(super) fn solve_newton_direction_dense(
hessian: &Array2<f64>,
gradient: &Array1<f64>,
direction_out: &mut Array1<f64>,
) -> Result<(), EstimationError> {
solve_newton_direction_dense_with_factor(hessian, gradient, direction_out).map(|_| ())
}
pub(super) fn solve_direction_with_dense_factor(
factor: &FaerSymmetricFactor,
gradient: &Array1<f64>,
direction_out: &mut Array1<f64>,
) {
if direction_out.len() != gradient.len() {
*direction_out = Array1::zeros(gradient.len());
}
direction_out.assign(gradient);
let mut rhsview = array1_to_col_matmut(direction_out);
factor.solve_in_place(rhsview.as_mut());
direction_out.mapv_inplace(|v| -v);
}
pub(super) fn solve_newton_direction_dense_with_factor(
hessian: &Array2<f64>,
gradient: &Array1<f64>,
direction_out: &mut Array1<f64>,
) -> Result<Option<FaerSymmetricFactor>, EstimationError> {
let dense_solve_start = std::time::Instant::now();
let p = hessian.nrows();
if direction_out.len() != gradient.len() {
*direction_out = Array1::zeros(gradient.len());
}
if crate::solver::gpu::cuda_selected() {
let rhs = Array2::from_shape_vec((p, 1), gradient.to_vec()).map_err(|e| {
EstimationError::InvalidInput(format!("CUDA PIRLS RHS layout failed: {e}"))
})?;
let (solved, _) =
crate::solver::gpu::pirls_gpu::cholesky_solve_gpu(hessian.view(), rhs.view())
.map_err(EstimationError::InvalidInput)?;
direction_out.assign(&solved.column(0));
direction_out.mapv_inplace(|v| -v);
if array_is_finite(direction_out) {
log::info!(
"[STAGE] PIRLS dense newton solve backend=CUDA p={} flops~{} elapsed={:.3}s route=\"cuSOLVER potrf/potrs\"",
p,
(p as u64).saturating_mul((p as u64).saturating_mul(p as u64)) / 3,
dense_solve_start.elapsed().as_secs_f64(),
);
return Ok(None);
}
}
let cpu_route = String::from("CPU stable solver");
let factor = StableSolver::new("pirls newton direction")
.factorize(hessian)
.map_err(EstimationError::LinearSystemSolveFailed)?;
solve_direction_with_dense_factor(&factor, gradient, direction_out);
let validation_residual = {
let h_delta = hessian.dot(direction_out);
h_delta
.iter()
.zip(gradient.iter())
.map(|(h, g)| (h + g).abs())
.fold(0.0_f64, f64::max)
};
let g_inf = gradient.iter().map(|v| v.abs()).fold(0.0_f64, f64::max);
let rel = validation_residual / (1.0 + g_inf);
if !rel.is_finite() || rel > 1.0e-3 {
let rhs = gradient.mapv(|v| -v);
if let Some(pseudo) = StableSolver::new("pirls newton direction (pseudoinverse fallback)")
.solve_with_pseudoinverse_fallback(hessian, &rhs, 1.0e-10, 1.0e-3, 1.0e-10)
{
direction_out.assign(&pseudo);
log::info!(
"[STAGE] PIRLS dense newton solve backend=CPU p={} elapsed={:.3}s route=\"{} + pseudoinverse fallback (rel={:.3e} > 1e-3)\"",
p,
dense_solve_start.elapsed().as_secs_f64(),
cpu_route,
rel,
);
return Ok(Some(factor));
}
}
if array_is_finite(direction_out) {
log::info!(
"[STAGE] PIRLS dense newton solve backend=CPU p={} flops~{} elapsed={:.3}s route=\"{}\"",
p,
(p as u64).saturating_mul((p as u64).saturating_mul(p as u64)) / 3,
dense_solve_start.elapsed().as_secs_f64(),
cpu_route,
);
return Ok(Some(factor));
}
Err(EstimationError::LinearSystemSolveFailed(
FaerLinalgError::FactorizationFailed {
context: "PIRLS dense newton solve exhausted",
},
))
}
pub fn solve_newton_direction_implicit<F>(
apply_xtwx: F,
xtwx_diag: ArrayView1<'_, f64>,
dense_penalties: &[(f64, &Array2<f64>)],
op_penalties: &[(f64, &dyn crate::terms::penalty_op::PenaltyOp)],
gradient: &Array1<f64>,
direction_out: &mut Array1<f64>,
ridge: f64,
rel_tol: f64,
max_iter: usize,
) -> Result<(), EstimationError>
where
F: Fn(&Array1<f64>) -> Array1<f64>,
{
let p = gradient.len();
if xtwx_diag.len() != p {
crate::bail_invalid_estim!(
"solve_newton_direction_implicit: xtwx_diag length {} != gradient length {}",
xtwx_diag.len(),
p
);
}
for (_, s) in dense_penalties.iter() {
if s.nrows() != p || s.ncols() != p {
crate::bail_invalid_estim!(
"solve_newton_direction_implicit: dense penalty dim {}×{} != p={}",
s.nrows(),
s.ncols(),
p
);
}
}
for (_, op) in op_penalties.iter() {
if op.dim() != p {
crate::bail_invalid_estim!(
"solve_newton_direction_implicit: op penalty dim {} != p={}",
op.dim(),
p
);
}
}
if direction_out.len() != p {
*direction_out = Array1::zeros(p);
}
let pcg_start = std::time::Instant::now();
let mut precond_diag = xtwx_diag.to_owned();
if ridge > 0.0 {
precond_diag.mapv_inplace(|d| d + ridge);
}
for (lambda, s) in dense_penalties.iter() {
if *lambda == 0.0 {
continue;
}
for i in 0..p {
precond_diag[i] += *lambda * s[[i, i]];
}
}
for (lambda, op) in op_penalties.iter() {
if *lambda == 0.0 {
continue;
}
let d = op.diag();
for i in 0..p {
precond_diag[i] += *lambda * d[i];
}
}
let apply_h = |v: &Array1<f64>| -> Array1<f64> {
let mut hv = apply_xtwx(v);
if ridge > 0.0 {
hv.zip_mut_with(v, |h, &x| *h += ridge * x);
}
for (lambda, s) in dense_penalties.iter() {
if *lambda == 0.0 {
continue;
}
let sv = fast_av(s, v);
hv.scaled_add(*lambda, &sv);
}
for (lambda, op) in op_penalties.iter() {
if *lambda == 0.0 {
continue;
}
let mut sv = Array1::<f64>::zeros(p);
op.matvec(v.view(), sv.view_mut());
hv.scaled_add(*lambda, &sv);
}
hv
};
let solution =
crate::linalg::utils::solve_spd_pcg(apply_h, gradient, &precond_diag, rel_tol, max_iter)
.ok_or(EstimationError::LinearSystemSolveFailed(
FaerLinalgError::FactorizationFailed {
context: "PIRLS implicit PCG solve exhausted",
},
))?;
direction_out.assign(&solution);
direction_out.mapv_inplace(|v| -v);
if !array_is_finite(direction_out) {
return Err(EstimationError::LinearSystemSolveFailed(
FaerLinalgError::FactorizationFailed {
context: "PIRLS implicit PCG non-finite direction",
},
));
}
log::info!(
"[STAGE] PIRLS implicit (PCG) newton solve p={} dense_pens={} op_pens={} elapsed={:.3}s",
p,
dense_penalties.len(),
op_penalties.len(),
pcg_start.elapsed().as_secs_f64(),
);
Ok(())
}
pub(super) fn project_coefficients_to_lower_bounds(
beta: &mut Array1<f64>,
lower_bounds: &Array1<f64>,
) {
for i in 0..beta.len() {
let lb = lower_bounds[i];
if lb.is_finite() && beta[i] < lb {
beta[i] = lb;
}
}
}
pub(super) fn projected_gradient_norm(
gradient: &Array1<f64>,
beta: &Array1<f64>,
lower_bounds: Option<&Array1<f64>>,
) -> f64 {
let Some(lb) = lower_bounds else {
return gradient.dot(gradient).sqrt();
};
let mut sum_sq = 0.0;
for i in 0..gradient.len() {
let g = gradient[i];
if lb[i].is_finite() && g > 0.0 {
let slack = beta[i] - lb[i];
let scale = beta[i].abs().max(lb[i].abs()).max(1.0);
let tol = 1e-6 * scale + 1e-10;
if slack < tol {
continue;
}
}
sum_sq += g * g;
}
sum_sq.sqrt()
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub(super) enum PirlsSoftAccept {
NearStationaryPlateau,
BoundarySaturation,
RelativeBandPlateau,
}
#[derive(Clone, Copy, Debug)]
pub(super) enum SoftAcceptProgress {
Realized { dev_change: f64 },
Predicted {
predicted_reduction: f64,
current_penalized: f64,
},
}
#[inline]
pub(super) fn pirls_soft_acceptance(
state: &WorkingState,
projected_grad: f64,
progress: SoftAcceptProgress,
max_abs_eta: f64,
progress_tol: f64,
kkt_tol: f64,
) -> Option<PirlsSoftAccept> {
let objective_scale = state.deviance.abs().max(state.penalty_term.abs()).max(1.0);
let scaled_dev_tol = progress_tol * objective_scale;
let near_stationary_plateau = match progress {
SoftAcceptProgress::Realized { dev_change } => {
state.near_stationary_kkt(projected_grad, kkt_tol) && dev_change.abs() < scaled_dev_tol
}
SoftAcceptProgress::Predicted {
predicted_reduction,
current_penalized,
} => {
let reduction_noise_floor = current_penalized.abs().max(1.0) * 1e-12;
state.near_stationary_kkt(projected_grad, kkt_tol)
&& predicted_reduction.abs() <= reduction_noise_floor
}
};
if near_stationary_plateau {
return Some(PirlsSoftAccept::NearStationaryPlateau);
}
let dev_change = match progress {
SoftAcceptProgress::Realized { dev_change } => dev_change,
SoftAcceptProgress::Predicted { .. } => return None,
};
if max_abs_eta >= PIRLS_ETA_ABS_CAP * (1.0 - 1e-12) && dev_change.abs() < scaled_dev_tol {
return Some(PirlsSoftAccept::BoundarySaturation);
}
if projected_grad <= progress_tol.max(1e-6) * objective_scale
&& dev_change.abs() < scaled_dev_tol * 0.1
&& dev_change >= 0.0
{
return Some(PirlsSoftAccept::RelativeBandPlateau);
}
None
}
pub(super) fn constrained_stationarity_norm(
gradient: &Array1<f64>,
beta: &Array1<f64>,
lower_bounds: Option<&Array1<f64>>,
linear_constraints: Option<&LinearInequalityConstraints>,
) -> f64 {
if let Some(constraints) = linear_constraints {
let kkt = compute_constraint_kkt_diagnostics(beta, gradient, constraints);
return kkt
.primal_feasibility
.max(kkt.dual_feasibility)
.max(kkt.complementarity)
.max(kkt.stationarity);
}
projected_gradient_norm(gradient, beta, lower_bounds)
}
fn chunk_rows_for_nnz_count(n: usize, p: usize) -> usize {
const TARGET_BYTES: usize = 8 * 1024 * 1024;
const MIN_ROWS: usize = 256;
const MAX_ROWS: usize = 65_536;
if p == 0 {
return n.max(1);
}
(TARGET_BYTES / (p * 8))
.clamp(MIN_ROWS, MAX_ROWS)
.min(n.max(1))
}
fn count_dense_upper_nnz(matrix: &Array2<f64>, tol: f64) -> usize {
let p = matrix.nrows().min(matrix.ncols());
let mut nnz = 0usize;
for col in 0..p {
for row in 0..=col {
if matrix[[row, col]].abs() > tol {
nnz += 1;
}
}
}
nnz
}
fn estimate_sparse_native_decision(
workspace: &mut PirlsWorkspace,
x_original: &DesignMatrix,
s_lambda: &Array2<f64>,
coefficient_lower_bounds: Option<&Array1<f64>>,
linear_constraints_original: Option<&LinearInequalityConstraints>,
) -> SparsePirlsDecision {
let p = x_original.ncols();
let nnz_s_lambda = count_dense_upper_nnz(s_lambda, 1e-12);
let dense_reject = |reason: &'static str, nnz_x: usize| SparsePirlsDecision {
path: PirlsLinearSolvePath::DenseTransformed,
reason,
p,
nnz_x,
nnz_xtwx_symbolic: None,
nnz_s_lambda,
nnz_h_est: None,
density_h_est: None,
};
let has_finite_lower_bounds = coefficient_lower_bounds
.map(|lb| lb.iter().any(|bound| bound.is_finite()))
.unwrap_or(false);
if has_finite_lower_bounds || linear_constraints_original.is_some() {
return dense_reject("constraints_present", 0);
}
let x_sparse = if let Some(sparse) = x_original.as_sparse() {
sparse
} else {
let row_chunk_start = std::time::Instant::now();
let n = x_original.nrows();
let chunk = chunk_rows_for_nnz_count(n, x_original.ncols());
let mut nnz: usize = 0;
let mut chunks_processed = 0usize;
if chunk > 0 && n > 0 {
let mut start = 0;
while start < n {
let end = (start + chunk).min(n);
chunks_processed += 1;
match x_original.try_row_chunk(start..end) {
Ok(rows) => {
nnz = nnz.saturating_add(rows.iter().filter(|v| v.abs() > 1e-12).count());
}
Err(_) => {
nnz = nnz.saturating_add((end - start).saturating_mul(x_original.ncols()));
}
}
start = end;
}
}
log::info!(
"[STAGE] PIRLS row-chunk generation chunks={} n={} p={} nnz={} elapsed={:.3}s",
chunks_processed,
n,
x_original.ncols(),
nnz,
row_chunk_start.elapsed().as_secs_f64(),
);
return dense_reject("design_not_sparse", nnz);
};
let nnz_x = x_sparse.val().len();
match workspace.sparse_penalized_system_stats(x_sparse, s_lambda) {
Ok(stats) => SparsePirlsDecision {
path: if stats.density_upper <= SPARSE_NATIVE_MAX_H_DENSITY {
PirlsLinearSolvePath::SparseNative
} else {
PirlsLinearSolvePath::DenseTransformed
},
reason: if stats.density_upper <= SPARSE_NATIVE_MAX_H_DENSITY {
"sparse_native_eligible"
} else {
"penalized_hessian_too_dense"
},
p,
nnz_x,
nnz_xtwx_symbolic: Some(stats.nnz_xtwx_symbolic),
nnz_s_lambda: stats.nnz_s_lambda_upper,
nnz_h_est: Some(stats.nnz_h_upper),
density_h_est: Some(stats.density_upper),
},
Err(_) => dense_reject("sparse_stats_failed", nnz_x),
}
}
pub(super) fn should_use_sparse_native_pirls(
workspace: &mut PirlsWorkspace,
x_original: &DesignMatrix,
s_lambda: &Array2<f64>,
coefficient_lower_bounds: Option<&Array1<f64>>,
linear_constraints_original: Option<&LinearInequalityConstraints>,
) -> SparsePirlsDecision {
estimate_sparse_native_decision(
workspace,
x_original,
s_lambda,
coefficient_lower_bounds,
linear_constraints_original,
)
}
pub(crate) fn sparse_reml_penalized_hessian(
workspace: &mut PirlsWorkspace,
x: &SparseColMat<usize, f64>,
weights: &Array1<f64>,
s_lambda: &Array2<f64>,
ridge: f64,
precomputed_xtwx: Option<&SparseXtwxPrecomputed>,
) -> Result<SparseColMat<usize, f64>, EstimationError> {
workspace.assemble_sparse_penalized_hessian(x, weights, s_lambda, ridge, precomputed_xtwx)
}
pub(super) fn ensure_sparse_positive_definitewithridge<F>(
mut assemble: F,
) -> Result<
(
SparseColMat<usize, f64>,
crate::linalg::sparse_exact::SparseExactFactor,
f64,
),
EstimationError,
>
where
F: FnMut(f64) -> Result<SparseColMat<usize, f64>, EstimationError>,
{
let mut ridge = 0.0_f64;
for _ in 0..16 {
let h = assemble(ridge)?;
match factorize_sparse_spd(&h) {
Ok(factor) => return Ok((h, factor, ridge)),
Err(_) => {
ridge = if ridge == 0.0 {
FIXED_STABILIZATION_RIDGE
} else {
ridge * 10.0
};
}
}
}
Err(EstimationError::HessianNotPositiveDefinite {
min_eigenvalue: f64::NAN,
})
}
fn solve_subsystem_direction(
h_sub: ndarray::ArrayView2<f64>,
g_sub: ndarray::ArrayView1<f64>,
out: &mut Array1<f64>,
) -> Result<(), EstimationError> {
let n = g_sub.len();
if out.len() != n {
*out = Array1::zeros(n);
}
if let Ok(factor) = StableSolver::new("pirls bounded subsystem").factorize_any(&h_sub) {
out.assign(&g_sub);
let mut rhs = array1_to_col_matmut(out);
factor.solve_in_place(rhs.as_mut());
out.mapv_inplace(|v| -v);
if array_is_finite(out) {
return Ok(());
}
}
let diag_scale = (0..n)
.map(|i| h_sub[[i, i]].abs())
.fold(0.0_f64, f64::max)
.max(1.0);
let mut tau = 1e-8 * diag_scale;
let mut h_reg = h_sub.to_owned();
for _ in 0..12 {
for i in 0..n {
h_reg[[i, i]] = h_sub[[i, i]] + tau;
}
if let Ok(factor) = StableSolver::new("pirls bounded subsystem ridge").factorize(&h_reg) {
out.assign(&g_sub);
let mut rhs = array1_to_col_matmut(out);
factor.solve_in_place(rhs.as_mut());
out.mapv_inplace(|v| -v);
if array_is_finite(out) {
return Ok(());
}
}
tau *= 10.0;
}
let gnorm = g_sub.dot(&g_sub).sqrt();
if gnorm > 0.0 {
let scale = 1.0 / gnorm.max(diag_scale);
for i in 0..n {
out[i] = -g_sub[i] * scale;
}
return Ok(());
}
out.fill(0.0);
Ok(())
}
pub(super) fn linear_constraints_from_lower_bounds(
lower_bounds: &Array1<f64>,
) -> Option<LinearInequalityConstraints> {
LinearInequalityConstraints::from_per_coordinate_lower_bounds(lower_bounds)
}
pub(super) fn compute_constraint_kkt_diagnostics(
beta: &Array1<f64>,
gradient: &Array1<f64>,
constraints: &LinearInequalityConstraints,
) -> ConstraintKktDiagnostics {
active_set::compute_constraint_kkt_diagnostics(beta, gradient, constraints)
}
pub(super) fn select_active_set_release(
gradient: &Array1<f64>,
hd: &Array1<f64>,
active_idx: &[usize],
use_blands: bool,
) -> Option<usize> {
if use_blands {
for &i in active_idx {
let lambda_i = gradient[i] + hd[i];
let scale = gradient[i].abs().max(hd[i].abs()).max(1.0);
let tol = 64.0 * f64::EPSILON * scale;
if lambda_i < -tol {
return Some(i);
}
}
None
} else {
let mut worst = 0.0_f64;
let mut idx = None;
for &i in active_idx {
let lambda_i = gradient[i] + hd[i];
if lambda_i < worst {
worst = lambda_i;
idx = Some(i);
}
}
idx
}
}
pub(crate) fn solve_newton_directionwith_lower_bounds(
hessian: &Array2<f64>,
gradient: &Array1<f64>,
beta: &Array1<f64>,
lower_bounds: &Array1<f64>,
direction_out: &mut Array1<f64>,
active_hint: Option<&mut Vec<usize>>,
) -> Result<(), EstimationError> {
let p = gradient.len();
if lower_bounds.len() != p || beta.len() != p {
crate::bail_invalid_estim!(
"lower-bound size mismatch: beta={}, gradient={}, bounds={}",
beta.len(),
gradient.len(),
lower_bounds.len()
);
}
if direction_out.len() != p {
*direction_out = Array1::zeros(p);
}
direction_out.fill(0.0);
let has_active_hint = active_hint
.as_ref()
.map(|hint| !hint.is_empty())
.unwrap_or(false);
if !has_active_hint && solve_newton_direction_dense(hessian, gradient, direction_out).is_ok() {
let mut feasible = true;
for i in 0..p {
let lb = lower_bounds[i];
if lb.is_finite() && beta[i] + direction_out[i] < lb {
feasible = false;
break;
}
}
if feasible {
return Ok(());
}
}
let mut active = vec![false; p];
if let Some(hint) = active_hint.as_ref() {
for &idx in hint.iter() {
if idx < p {
active[idx] = true;
}
}
}
for i in 0..p {
let lb = lower_bounds[i];
if lb.is_finite() && gradient[i] > 0.0 {
let scale = beta[i].abs().max(lb.abs()).max(1.0);
let tol = 1e-6 * scale + 1e-10;
if beta[i] <= lb + tol {
active[i] = true;
}
}
}
const BLANDS_RULE_GRACE: usize = 2;
let blands_threshold = BLANDS_RULE_GRACE * (p + 1);
let max_iters = 8 * (p + 1);
let mut d_free = Array1::<f64>::zeros(p);
let mut h_ff_buf = Array2::<f64>::zeros((p, p));
let mut g_f_buf = Array1::<f64>::zeros(p);
for it in 0..max_iters {
let use_blands = it >= blands_threshold;
let free_idx: Vec<usize> = (0..p).filter(|&i| !active[i]).collect();
let active_idx: Vec<usize> = (0..p).filter(|&i| active[i]).collect();
direction_out.fill(0.0);
for &i in &active_idx {
let lb = lower_bounds[i];
if lb.is_finite() {
direction_out[i] = lb - beta[i];
}
}
if free_idx.is_empty() {
let hd = fast_av(hessian, direction_out);
if let Some(idx) = select_active_set_release(gradient, &hd, &active_idx, use_blands) {
active[idx] = false;
continue;
}
if let Some(hint) = active_hint {
hint.clear();
hint.extend((0..p).filter(|&i| active[i]));
}
return Ok(());
}
let n_free = free_idx.len();
{
let mut h_ff = h_ff_buf.slice_mut(ndarray::s![..n_free, ..n_free]);
let mut g_f = g_f_buf.slice_mut(ndarray::s![..n_free]);
for (ii, &i) in free_idx.iter().enumerate() {
let mut gi = gradient[i];
for &j in &active_idx {
gi += hessian[[i, j]] * direction_out[j];
}
g_f[ii] = gi;
for (jj, &j) in free_idx.iter().enumerate() {
h_ff[[ii, jj]] = hessian[[i, j]];
}
}
}
solve_subsystem_direction(
h_ff_buf.slice(ndarray::s![..n_free, ..n_free]),
g_f_buf.slice(ndarray::s![..n_free]),
&mut d_free,
)?;
for (ii, &i) in free_idx.iter().enumerate() {
direction_out[i] = d_free[ii];
}
let mut hit_idx: Option<usize> = None;
let mut best_alpha = 1.0_f64;
for &i in &free_idx {
let lb = lower_bounds[i];
if !lb.is_finite() {
continue;
}
let slack = beta[i] - lb;
let di = direction_out[i];
if let Some(alpha_i) = boundary_hit_step_fraction(slack, di, best_alpha) {
best_alpha = alpha_i;
hit_idx = Some(i);
}
}
if let Some(i_hit) = hit_idx {
for i in 0..p {
direction_out[i] *= best_alpha;
}
active[i_hit] = true;
continue;
}
let hd = fast_av(hessian, direction_out);
if let Some(idx) = select_active_set_release(gradient, &hd, &active_idx, use_blands) {
active[idx] = false;
continue;
}
if let Some(hint) = active_hint {
hint.clear();
hint.extend((0..p).filter(|&i| active[i]));
}
return Ok(());
}
let gnorm = gradient.dot(gradient).sqrt();
if gnorm > 0.0 {
let diag_scale = (0..p)
.map(|i| hessian[[i, i]].abs())
.fold(0.0_f64, f64::max)
.max(1.0);
let step_scale = 1.0 / diag_scale;
for i in 0..p {
let di = -gradient[i] * step_scale;
let lb = lower_bounds[i];
if lb.is_finite() && beta[i] + di < lb {
direction_out[i] = lb - beta[i];
} else {
direction_out[i] = di;
}
}
} else {
direction_out.fill(0.0);
}
if let Some(hint) = active_hint {
hint.clear();
}
Ok(())
}
pub(super) fn solve_newton_directionwith_linear_constraints(
hessian: &Array2<f64>,
gradient: &Array1<f64>,
beta: &Array1<f64>,
constraints: &LinearInequalityConstraints,
direction_out: &mut Array1<f64>,
active_hint: Option<&mut Vec<usize>>,
) -> Result<(), EstimationError> {
active_set::solve_newton_direction_with_linear_constraints(
hessian,
gradient,
beta,
constraints,
direction_out,
active_hint,
)
}
use loop_driver::assert_symmetric_tol;
pub(crate) use loop_driver::fit_model_for_fixed_rho_with_adaptive_kkt;
pub use loop_driver::{PenaltyConfig, PirlsConfig, PirlsProblem, fit_model_for_fixed_rho};
#[inline]
pub(super) fn standard_inverse_link_jet(
inverse_link: &InverseLink,
eta: f64,
) -> Result<MixtureInverseLinkJet, EstimationError> {
crate::mixture_link::inverse_link_jet_for_inverse_link(inverse_link, eta)
}
#[inline]
fn bernoulli_logit_geometry_from_jet(
eta_raw: f64,
eta_used: f64,
y: f64,
priorweight: f64,
jet: crate::mixture_link::LogitJet5,
zero_on_nonsmooth: bool,
) -> WorkingBernoulliGeometry {
let fisher = jet.d1;
let nonsmooth = eta_raw != eta_used || !fisher.is_finite() || fisher < 0.0;
let (c, d) = if nonsmooth && zero_on_nonsmooth {
(0.0, 0.0)
} else {
(priorweight * jet.d2, priorweight * jet.d3)
};
WorkingBernoulliGeometry {
mu: jet.mu,
weight: priorweight * fisher,
z: bernoulli_exact_working_response(eta_used, y, jet.mu, jet.d1),
c,
d,
}
}
#[inline]
fn bernoulli_geometry_from_jet(
eta_raw: f64,
eta_used: f64,
y: f64,
priorweight: f64,
jet: MixtureInverseLinkJet,
) -> WorkingBernoulliGeometry {
let mu = jet.mu;
let v = mu * (1.0 - mu);
let n0 = jet.d1 * jet.d1;
let fisher = if v.is_finite() && v > 0.0 {
n0 / v
} else {
0.0
};
let nonsmooth =
eta_raw != eta_used || !v.is_finite() || v <= 0.0 || !fisher.is_finite() || fisher < 0.0;
let (c, d) = if nonsmooth {
(0.0, 0.0)
} else {
let v1 = jet.d1 * (1.0 - 2.0 * mu);
let v2 = jet.d2 * (1.0 - 2.0 * mu) - 2.0 * jet.d1 * jet.d1;
let n1 = 2.0 * jet.d1 * jet.d2;
let n2 = 2.0 * (jet.d2 * jet.d2 + jet.d1 * jet.d3);
let numer1 = n1 * v - n0 * v1;
let c = priorweight * numer1 / (v * v);
let d = priorweight * ((n2 * v - n0 * v2) / (v * v) - 2.0 * numer1 * v1 / (v * v * v));
(c, d)
};
WorkingBernoulliGeometry {
mu,
weight: priorweight * fisher,
z: bernoulli_exact_working_response(eta_used, y, mu, jet.d1),
c,
d,
}
}
#[inline]
fn bernoulli_exact_working_response(eta: f64, y: f64, mu: f64, dmu_deta: f64) -> f64 {
if dmu_deta.is_finite() && dmu_deta > 0.0 {
let delta = (y - mu) / dmu_deta;
if delta.is_finite() {
return eta + delta;
}
}
eta
}
#[inline]
fn write_identityworking_state(
y: ArrayView1<f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) {
mu.assign(eta);
weights.assign(&priorweights);
z.assign(&y);
if let Some(derivs) = derivatives {
derivs.c.fill(0.0);
derivs.d.fill(0.0);
derivs.dmu_deta.fill(1.0);
derivs.d2mu_deta2.fill(0.0);
derivs.d3mu_deta3.fill(0.0);
}
}
#[inline]
fn write_poisson_log_working_state(
y: ArrayView1<f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) {
const MIN_MU: f64 = 1e-10;
const MIN_WEIGHT: f64 = 1e-12;
if let Some(derivs) = derivatives {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
let dmu_s = derivs
.dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = derivs
.d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = derivs
.d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
let c_s = derivs.c.as_slice_mut().expect("c must be contiguous");
let d_s = derivs.d.as_slice_mut().expect("d must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.zip(c_s.par_iter_mut())
.zip(d_s.par_iter_mut())
.enumerate()
.for_each(
|(i, (((((((mu_o, w_o), z_o), dmu_o), d2_o), d3_o), c_o), d_o))| {
let eta_raw = eta[i];
let eta_i = eta_raw.clamp(-700.0, 700.0);
let mu_i = eta_i.exp().max(MIN_MU);
*mu_o = mu_i;
let raw_weight = priorweights[i].max(0.0) * mu_i;
let floor_active = raw_weight > 0.0 && raw_weight <= MIN_WEIGHT;
*w_o = if raw_weight > 0.0 {
raw_weight.max(MIN_WEIGHT)
} else {
0.0
};
*z_o = eta_i + (y[i] - mu_i) / mu_i;
*dmu_o = mu_i;
*d2_o = mu_i;
*d3_o = mu_i;
if floor_active || eta_raw != eta_i {
*c_o = 0.0;
*d_o = 0.0;
} else {
*c_o = raw_weight;
*d_o = raw_weight;
}
},
);
} else {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.enumerate()
.for_each(|(i, ((mu_o, w_o), z_o))| {
let eta_i = eta[i].clamp(-700.0, 700.0);
let mu_i = eta_i.exp().max(MIN_MU);
*mu_o = mu_i;
let raw_weight = priorweights[i].max(0.0) * mu_i;
*w_o = if raw_weight > 0.0 {
raw_weight.max(MIN_WEIGHT)
} else {
0.0
};
*z_o = eta_i + (y[i] - mu_i) / mu_i;
});
}
}
#[inline]
fn write_gamma_log_working_state(
y: ArrayView1<f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
shape: f64,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) {
const MIN_MU: f64 = 1e-10;
if let Some(derivs) = derivatives {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
let dmu_s = derivs
.dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = derivs
.d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = derivs
.d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
let c_s = derivs.c.as_slice_mut().expect("c must be contiguous");
let d_s = derivs.d.as_slice_mut().expect("d must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.zip(c_s.par_iter_mut())
.zip(d_s.par_iter_mut())
.enumerate()
.for_each(
|(i, (((((((mu_o, w_o), z_o), dmu_o), d2_o), d3_o), c_o), d_o))| {
let eta_i = eta[i].clamp(-700.0, 700.0);
let mu_i = eta_i.exp().max(MIN_MU);
*mu_o = mu_i;
*w_o = priorweights[i].max(0.0) * shape;
*z_o = eta_i + (y[i] - mu_i) / mu_i;
*dmu_o = mu_i;
*d2_o = mu_i;
*d3_o = mu_i;
*c_o = 0.0;
*d_o = 0.0;
},
);
} else {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.enumerate()
.for_each(|(i, ((mu_o, w_o), z_o))| {
let eta_i = eta[i].clamp(-700.0, 700.0);
let mu_i = eta_i.exp().max(MIN_MU);
*mu_o = mu_i;
*w_o = priorweights[i].max(0.0) * shape;
*z_o = eta_i + (y[i] - mu_i) / mu_i;
});
}
}
pub const BETA_MU_EPS: f64 = 1.0e-12;
#[inline]
fn tweedie_log_weight_mu_power(mu: f64, p: f64) -> f64 {
mu.max(1.0e-300).powf(2.0 - p)
}
#[inline]
fn valid_negbin_theta(theta: f64) -> bool {
theta.is_finite() && theta > 0.0
}
#[inline]
fn valid_count_response(y: f64) -> bool {
y.is_finite() && y >= 0.0 && (y - y.round()).abs() <= 1e-9
}
fn validate_count_responses(
y: &ArrayView1<'_, f64>,
priorweights: &ArrayView1<'_, f64>,
family: &str,
) -> Result<(), EstimationError> {
for (i, (&yi, &wi)) in y.iter().zip(priorweights.iter()).enumerate() {
if wi > 0.0 && !valid_count_response(yi) {
crate::bail_invalid_estim!(
"{family} response must be a finite non-negative integer at positive-weight row {i}; got {yi}"
);
}
}
Ok(())
}
#[inline]
fn valid_beta_phi(phi: f64) -> bool {
phi.is_finite() && phi > 0.0
}
#[inline]
fn valid_beta_response(y: f64) -> bool {
y.is_finite() && y > 0.0 && y < 1.0
}
fn validate_beta_responses(
y: &ArrayView1<'_, f64>,
priorweights: &ArrayView1<'_, f64>,
) -> Result<(), EstimationError> {
for (i, (&yi, &wi)) in y.iter().zip(priorweights.iter()).enumerate() {
if wi > 0.0 && !valid_beta_response(yi) {
crate::bail_invalid_estim!(
"beta-regression response must be finite and strictly inside (0, 1) at positive-weight row {i}; got {yi}"
);
}
}
Ok(())
}
#[inline]
fn valid_tweedie_response(y: f64) -> bool {
y.is_finite() && y >= 0.0
}
fn validate_tweedie_responses(
y: &ArrayView1<'_, f64>,
priorweights: &ArrayView1<'_, f64>,
) -> Result<(), EstimationError> {
for (i, (&yi, &wi)) in y.iter().zip(priorweights.iter()).enumerate() {
if wi > 0.0 && !valid_tweedie_response(yi) {
crate::bail_invalid_estim!(
"Tweedie response must be finite and non-negative at positive-weight row {i}; got {yi}"
);
}
}
Ok(())
}
#[inline]
fn safe_beta_mu(mu: f64) -> f64 {
mu.clamp(BETA_MU_EPS, 1.0 - BETA_MU_EPS)
}
#[inline]
fn trigamma(mut x: f64) -> f64 {
if !(x.is_finite() && x > 0.0) {
return f64::NAN;
}
let mut acc = 0.0;
while x < 8.0 {
acc += 1.0 / (x * x);
x += 1.0;
}
let inv = 1.0 / x;
let inv2 = inv * inv;
acc + inv + 0.5 * inv2 + inv2 * inv / 6.0 - inv2 * inv2 * inv / 30.0
+ inv2 * inv2 * inv2 * inv / 42.0
- inv2 * inv2 * inv2 * inv2 * inv / 30.0
}
#[inline]
fn polygamma2(mut x: f64) -> f64 {
if !(x.is_finite() && x > 0.0) {
return f64::NAN;
}
let mut acc = 0.0;
while x < 8.0 {
acc -= 2.0 / (x * x * x);
x += 1.0;
}
let inv = 1.0 / x;
let inv2 = inv * inv;
let inv3 = inv2 * inv;
acc - inv2 - inv3 - 0.5 * inv2 * inv2 + inv3 * inv3 / 6.0 - inv2 * inv3 * inv3 / 6.0
+ 0.3 * inv2 * inv2 * inv3 * inv3
- 5.0 * inv2 * inv2 * inv2 * inv3 * inv3 / 6.0
}
#[inline]
fn polygamma3(mut x: f64) -> f64 {
if !(x.is_finite() && x > 0.0) {
return f64::NAN;
}
let mut acc = 0.0;
while x < 8.0 {
acc += 6.0 / (x * x * x * x);
x += 1.0;
}
let inv = 1.0 / x;
let inv2 = inv * inv;
let inv3 = inv2 * inv;
let inv4 = inv2 * inv2;
acc + 2.0 * inv3 + 3.0 * inv4 + 2.0 * inv4 * inv - inv4 * inv3 + 4.0 * inv4 * inv3 * inv2 / 3.0
- 3.0 * inv4 * inv3 * inv4
+ 10.0 * inv4 * inv4 * inv4 * inv
}
#[inline]
fn beta_logit_working_curvature_eta_derivatives(
prior_weight: f64,
phi: f64,
mu: f64,
q: f64,
a: f64,
b: f64,
trigamma_sum: f64,
) -> (f64, f64) {
let q_prime = q * (1.0 - 2.0 * mu);
let q_double_prime = q * (1.0 - 2.0 * mu) * (1.0 - 2.0 * mu) - 2.0 * q * q;
let psi2_diff = polygamma2(a) - polygamma2(b);
let psi3_sum = polygamma3(a) + polygamma3(b);
let phi_sq = phi * phi;
let q_sq = q * q;
let c = prior_weight * phi_sq * (2.0 * q * q_prime * trigamma_sum + q_sq * phi * q * psi2_diff);
let d = prior_weight
* phi_sq
* (2.0 * (q_prime * q_prime + q * q_double_prime) * trigamma_sum
+ 4.0 * q * q_prime * phi * q * psi2_diff
+ q_sq * (phi * q_prime * psi2_diff + phi_sq * q_sq * psi3_sum));
(c, d)
}
#[inline]
fn write_tweedie_log_working_state(
y: ArrayView1<f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
p: f64,
phi: f64,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) -> Result<(), EstimationError> {
const MIN_MU: f64 = 1e-10;
const MIN_WEIGHT: f64 = 1e-12;
if !is_valid_tweedie_power(p) {
crate::bail_invalid_estim!(
"Tweedie variance power must be finite and strictly between 1 and 2; got {p}"
);
}
if !(phi.is_finite() && phi > 0.0) {
crate::bail_invalid_estim!("Tweedie dispersion phi must be finite and > 0; got {phi}");
}
validate_tweedie_responses(&y, &priorweights)?;
let exponent = 2.0 - p;
if let Some(derivs) = derivatives {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
let dmu_s = derivs
.dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = derivs
.d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = derivs
.d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
let c_s = derivs.c.as_slice_mut().expect("c must be contiguous");
let d_s = derivs.d.as_slice_mut().expect("d must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.zip(c_s.par_iter_mut())
.zip(d_s.par_iter_mut())
.enumerate()
.for_each(
|(i, (((((((mu_o, w_o), z_o), dmu_o), d2_o), d3_o), c_o), d_o))| {
let eta_raw = eta[i];
let eta_i = eta_raw.clamp(-700.0, 700.0);
let clamp_active = eta_raw != eta_i;
let mu_i = eta_i.exp().max(MIN_MU);
*mu_o = mu_i;
let raw_weight =
priorweights[i].max(0.0) * tweedie_log_weight_mu_power(mu_i, p) / phi;
let floor_active = raw_weight > 0.0 && raw_weight <= MIN_WEIGHT;
*w_o = if raw_weight > 0.0 {
raw_weight.max(MIN_WEIGHT)
} else {
0.0
};
*z_o = eta_i + (y[i] - mu_i) / mu_i;
if clamp_active {
*dmu_o = 0.0;
*d2_o = 0.0;
*d3_o = 0.0;
} else {
*dmu_o = mu_i;
*d2_o = mu_i;
*d3_o = mu_i;
}
if floor_active || clamp_active {
*c_o = 0.0;
*d_o = 0.0;
} else {
*c_o = exponent * raw_weight;
*d_o = exponent * exponent * raw_weight;
}
},
);
} else {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.enumerate()
.for_each(|(i, ((mu_o, w_o), z_o))| {
let eta_i = eta[i].clamp(-700.0, 700.0);
let mu_i = eta_i.exp().max(MIN_MU);
*mu_o = mu_i;
let raw_weight =
priorweights[i].max(0.0) * tweedie_log_weight_mu_power(mu_i, p) / phi;
*w_o = if raw_weight > 0.0 {
raw_weight.max(MIN_WEIGHT)
} else {
0.0
};
*z_o = eta_i + (y[i] - mu_i) / mu_i;
});
}
Ok(())
}
#[inline]
fn write_negative_binomial_log_working_state(
y: ArrayView1<f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
theta: f64,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) -> Result<(), EstimationError> {
const MIN_MU: f64 = 1e-10;
const MIN_WEIGHT: f64 = 1e-12;
if !valid_negbin_theta(theta) {
crate::bail_invalid_estim!("negative-binomial theta must be finite and > 0; got {theta}");
}
validate_count_responses(&y, &priorweights, "negative-binomial")?;
if let Some(derivs) = derivatives {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
let dmu_s = derivs
.dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = derivs
.d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = derivs
.d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
let c_s = derivs.c.as_slice_mut().expect("c must be contiguous");
let d_s = derivs.d.as_slice_mut().expect("d must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.zip(c_s.par_iter_mut())
.zip(d_s.par_iter_mut())
.enumerate()
.for_each(
|(i, (((((((mu_o, w_o), z_o), dmu_o), d2_o), d3_o), c_o), d_o))| {
let eta_raw = eta[i];
let eta_i = eta_raw.clamp(-700.0, 700.0);
let mu_i = eta_i.exp().max(MIN_MU);
let denom = theta + mu_i;
let negbin_weight = if theta > mu_i {
mu_i / (1.0 + mu_i / theta)
} else {
theta / (1.0 + theta / mu_i)
};
let raw_weight = priorweights[i].max(0.0) * negbin_weight;
let floor_active = raw_weight > 0.0 && raw_weight <= MIN_WEIGHT;
*mu_o = mu_i;
*w_o = if raw_weight > 0.0 {
raw_weight.max(MIN_WEIGHT)
} else {
0.0
};
*z_o = eta_i + (y[i] - mu_i) / mu_i;
*dmu_o = mu_i;
*d2_o = mu_i;
*d3_o = mu_i;
if floor_active || eta_raw != eta_i {
*c_o = 0.0;
*d_o = 0.0;
} else {
*c_o = raw_weight * theta / denom;
*d_o = raw_weight * theta * (theta - mu_i) / (denom * denom);
}
},
);
} else {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.enumerate()
.for_each(|(i, ((mu_o, w_o), z_o))| {
let eta_i = eta[i].clamp(-700.0, 700.0);
let mu_i = eta_i.exp().max(MIN_MU);
let negbin_weight = if theta > mu_i {
mu_i / (1.0 + mu_i / theta)
} else {
theta / (1.0 + theta / mu_i)
};
let raw_weight = priorweights[i].max(0.0) * negbin_weight;
*mu_o = mu_i;
*w_o = if raw_weight > 0.0 {
raw_weight.max(MIN_WEIGHT)
} else {
0.0
};
*z_o = eta_i + (y[i] - mu_i) / mu_i;
});
}
Ok(())
}
#[inline]
fn write_beta_logit_working_state(
y: ArrayView1<f64>,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
phi: f64,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) -> Result<(), EstimationError> {
const MIN_WEIGHT: f64 = 1e-12;
if !valid_beta_phi(phi) {
crate::bail_invalid_estim!("beta-regression phi must be finite and > 0; got {phi}");
}
validate_beta_responses(&y, &priorweights)?;
if let Some(derivs) = derivatives {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
let dmu_s = derivs
.dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = derivs
.d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = derivs
.d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
let c_s = derivs.c.as_slice_mut().expect("c must be contiguous");
let d_s = derivs.d.as_slice_mut().expect("d must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.zip(c_s.par_iter_mut())
.zip(d_s.par_iter_mut())
.enumerate()
.for_each(
|(i, (((((((mu_o, w_o), z_o), dmu_o), d2_o), d3_o), c_o), d_o))| {
let eta_raw = eta[i];
let eta_i = eta_raw.clamp(-700.0, 700.0);
let jet = logit_inverse_link_jet5(eta_i);
let mu_i = safe_beta_mu(jet.mu);
let q = (mu_i * (1.0 - mu_i)).max(BETA_MU_EPS);
let yi = y[i];
let a = (mu_i * phi).max(BETA_MU_EPS);
let b = ((1.0 - mu_i) * phi).max(BETA_MU_EPS);
let score_mu = phi * (digamma(b) - digamma(a) + yi.ln() - (1.0 - yi).ln());
let trigamma_sum = trigamma(a) + trigamma(b);
let info_mu = phi * phi * trigamma_sum;
let prior_weight = priorweights[i].max(0.0);
let raw_weight = prior_weight * q * q * info_mu;
let floor_active = raw_weight > 0.0 && raw_weight <= MIN_WEIGHT;
*mu_o = mu_i;
*w_o = if raw_weight > 0.0 {
raw_weight.max(MIN_WEIGHT)
} else {
0.0
};
*z_o = eta_i + score_mu / (q * info_mu).max(MIN_WEIGHT);
*dmu_o = q;
*d2_o = q * (1.0 - 2.0 * mu_i);
*d3_o = q * (1.0 - 6.0 * q);
if floor_active || eta_raw != eta_i {
*c_o = 0.0;
*d_o = 0.0;
} else {
let (c_i, d_i) = beta_logit_working_curvature_eta_derivatives(
prior_weight,
phi,
mu_i,
q,
a,
b,
trigamma_sum,
);
*c_o = c_i;
*d_o = d_i;
}
},
);
} else {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.enumerate()
.for_each(|(i, ((mu_o, w_o), z_o))| {
let eta_i = eta[i].clamp(-700.0, 700.0);
let jet = logit_inverse_link_jet5(eta_i);
let mu_i = safe_beta_mu(jet.mu);
let q = (mu_i * (1.0 - mu_i)).max(BETA_MU_EPS);
let yi = y[i];
let a = (mu_i * phi).max(BETA_MU_EPS);
let b = ((1.0 - mu_i) * phi).max(BETA_MU_EPS);
let score_mu = phi * (digamma(b) - digamma(a) + yi.ln() - (1.0 - yi).ln());
let info_mu = phi * phi * (trigamma(a) + trigamma(b));
let raw_weight = priorweights[i].max(0.0) * q * q * info_mu;
*mu_o = mu_i;
*w_o = if raw_weight > 0.0 {
raw_weight.max(MIN_WEIGHT)
} else {
0.0
};
*z_o = eta_i + score_mu / (q * info_mu).max(MIN_WEIGHT);
});
}
Ok(())
}
#[inline]
pub fn update_glmvectors(
y: ArrayView1<f64>,
eta: &Array1<f64>,
inverse_link: &InverseLink,
priorweights: ArrayView1<f64>,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) -> Result<(), EstimationError> {
let link = inverse_link.link_function();
if matches!(link, LinkFunction::Logit)
&& inverse_link.mixture_state().is_none()
&& inverse_link.sas_state().is_none()
{
if let Some(derivs) = derivatives {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
let c_s = derivs.c.as_slice_mut().expect("c must be contiguous");
let d_s = derivs.d.as_slice_mut().expect("d must be contiguous");
let dmu_s = derivs
.dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = derivs
.d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = derivs
.d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.zip(c_s.par_iter_mut())
.zip(d_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.for_each(
|(i, (((((((mu_o, w_o), z_o), c_o), d_o), dmu_o), d2_o), d3_o))| {
let eta_raw = eta[i];
let eta_c = eta_raw.clamp(-700.0, 700.0);
let jet = logit_inverse_link_jet5(eta_c);
let geom = bernoulli_logit_geometry_from_jet(
eta_raw,
eta_c,
y[i],
priorweights[i],
jet,
true,
);
*mu_o = geom.mu;
*w_o = geom.weight;
*z_o = geom.z;
*c_o = geom.c;
*d_o = geom.d;
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
},
);
} else {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.enumerate()
.for_each(|(i, ((mu_o, w_o), z_o))| {
let eta_raw = eta[i];
let eta_c = eta_raw.clamp(-700.0, 700.0);
let jet = logit_inverse_link_jet5(eta_c);
let geom = bernoulli_logit_geometry_from_jet(
eta_raw,
eta_c,
y[i],
priorweights[i],
jet,
true,
);
*mu_o = geom.mu;
*w_o = geom.weight;
*z_o = geom.z;
});
}
return Ok(());
}
match link {
LinkFunction::Logit
| LinkFunction::Probit
| LinkFunction::CLogLog
| LinkFunction::Sas
| LinkFunction::BetaLogistic => {
let zero_on_nonsmooth = matches!(link, LinkFunction::Logit);
if let Some(derivs) = derivatives {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
let c_s = derivs.c.as_slice_mut().expect("c must be contiguous");
let d_s = derivs.d.as_slice_mut().expect("d must be contiguous");
let dmu_s = derivs
.dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = derivs
.d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = derivs
.d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.zip(c_s.par_iter_mut())
.zip(d_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.try_for_each(
|(
i,
(((((((mu_o, w_o), z_o), c_o), d_o), dmu_o), d2_o), d3_o),
)|
-> Result<(), EstimationError> {
let eta_used = eta_for_observed_hessian_jet(inverse_link, eta[i]);
if matches!(link, LinkFunction::Logit) {
let jet = logit_inverse_link_jet5(eta_used);
let geom = bernoulli_logit_geometry_from_jet(
eta[i],
eta_used,
y[i],
priorweights[i],
jet,
zero_on_nonsmooth,
);
*mu_o = geom.mu;
*w_o = geom.weight;
*z_o = geom.z;
*c_o = geom.c;
*d_o = geom.d;
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
} else {
let jet = standard_inverse_link_jet(inverse_link, eta_used)?;
let geom = bernoulli_geometry_from_jet(
eta[i],
eta_used,
y[i],
priorweights[i],
jet,
);
*mu_o = geom.mu;
*w_o = geom.weight;
*z_o = geom.z;
*c_o = geom.c;
*d_o = geom.d;
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
}
Ok(())
},
)?;
} else {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.enumerate()
.try_for_each(|(i, ((mu_o, w_o), z_o))| -> Result<(), EstimationError> {
let eta_used = eta_for_observed_hessian_jet(inverse_link, eta[i]);
if matches!(link, LinkFunction::Logit) {
let jet = logit_inverse_link_jet5(eta_used);
let geom = bernoulli_logit_geometry_from_jet(
eta[i],
eta_used,
y[i],
priorweights[i],
jet,
zero_on_nonsmooth,
);
*mu_o = geom.mu;
*w_o = geom.weight;
*z_o = geom.z;
} else {
let jet = standard_inverse_link_jet(inverse_link, eta_used)?;
let geom = bernoulli_geometry_from_jet(
eta[i],
eta_used,
y[i],
priorweights[i],
jet,
);
*mu_o = geom.mu;
*w_o = geom.weight;
*z_o = geom.z;
}
Ok(())
})?;
}
Ok(())
}
LinkFunction::Identity => {
write_identityworking_state(y, eta, priorweights, mu, weights, z, derivatives);
Ok(())
}
LinkFunction::Log => {
write_poisson_log_working_state(y, eta, priorweights, mu, weights, z, derivatives);
Ok(())
}
}
}
#[inline]
pub fn update_glmvectors_by_family(
y: ArrayView1<f64>,
eta: &Array1<f64>,
likelihood: &GlmLikelihoodSpec,
priorweights: ArrayView1<f64>,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
) -> Result<(), EstimationError> {
likelihood.irls_update(y, eta, priorweights, mu, weights, z, None, None)
}
fn integrated_inverse_link_from_family(
spec: &LikelihoodSpec,
mixture_link_state: Option<&MixtureLinkState>,
sas_link_state: Option<&SasLinkState>,
) -> Result<InverseLink, EstimationError> {
match (&spec.response, &spec.link) {
(ResponseFamily::Binomial, InverseLink::Standard(StandardLink::Logit))
| (ResponseFamily::Binomial, InverseLink::Standard(StandardLink::Probit))
| (ResponseFamily::Binomial, InverseLink::Standard(StandardLink::CLogLog)) => {
Ok(spec.link.clone())
}
(ResponseFamily::Binomial, InverseLink::Sas(_)) => {
let state = sas_link_state.ok_or_else(|| {
EstimationError::InvalidInput(
"Integrated BinomialSas update requires explicit SasLinkState".to_string(),
)
})?;
Ok(InverseLink::Sas(*state))
}
(ResponseFamily::Binomial, InverseLink::BetaLogistic(_)) => {
let state = sas_link_state.ok_or_else(|| {
EstimationError::InvalidInput(
"Integrated BinomialBetaLogistic update requires explicit SasLinkState"
.to_string(),
)
})?;
Ok(InverseLink::BetaLogistic(*state))
}
(ResponseFamily::Binomial, InverseLink::Mixture(_)) => {
let state = mixture_link_state.ok_or_else(|| {
EstimationError::InvalidInput(
"Integrated BinomialMixture update requires explicit MixtureLinkState"
.to_string(),
)
})?;
Ok(InverseLink::Mixture(state.clone()))
}
_ => Err(EstimationError::InvalidInput(format!(
"Integrated link-runtime update is not supported for likelihood (response={:?}, link={:?})",
spec.response, spec.link
))),
}
}
#[inline]
pub fn update_glmvectors_integrated_for_link(
quadctx: &crate::quadrature::QuadratureContext,
y: ArrayView1<f64>,
eta: &Array1<f64>,
se: ArrayView1<f64>,
inverse_link: &InverseLink,
priorweights: ArrayView1<f64>,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
) -> Result<(), EstimationError> {
let link = inverse_link.link_function();
if !matches!(
inverse_link,
InverseLink::Standard(StandardLink::Logit)
| InverseLink::Standard(StandardLink::Probit)
| InverseLink::Standard(StandardLink::CLogLog)
| InverseLink::LatentCLogLog(_)
| InverseLink::Sas(_)
| InverseLink::BetaLogistic(_)
| InverseLink::Mixture(_)
) {
crate::bail_invalid_estim!(
"Integrated link-runtime update is not supported for inverse link {:?}",
inverse_link
);
}
if let Some(derivs) = derivatives {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
let c_s = derivs.c.as_slice_mut().expect("c must be contiguous");
let d_s = derivs.d.as_slice_mut().expect("d must be contiguous");
let dmu_s = derivs
.dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = derivs
.d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = derivs
.d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.zip(c_s.par_iter_mut())
.zip(d_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.try_for_each(
|(i, (((((((mu_o, w_o), z_o), c_o), d_o), dmu_o), d2_o), d3_o))|
-> Result<(), EstimationError> {
let jet = if let InverseLink::LatentCLogLog(state) = inverse_link {
crate::families::lognormal_kernel::latent_cloglog_inverse_link_jet(
quadctx,
eta[i],
se[i].hypot(state.latent_sd),
)?
} else if matches!(inverse_link, InverseLink::Standard(StandardLink::Logit)) {
crate::quadrature::integrated_logit_inverse_link_jet_pirls(
quadctx, eta[i], se[i],
)?
} else {
crate::quadrature::integrated_inverse_link_jetwith_state(
quadctx,
link,
eta[i],
se[i],
inverse_link.mixture_state(),
inverse_link.sas_state(),
)?
};
let local_jet = MixtureInverseLinkJet {
mu: jet.mean,
d1: jet.d1,
d2: jet.d2,
d3: jet.d3,
};
let e = eta[i].clamp(-700.0, 700.0);
let geom = bernoulli_geometry_from_jet(
eta[i],
e,
y[i],
priorweights[i],
local_jet,
);
*mu_o = geom.mu;
*w_o = geom.weight;
*z_o = geom.z;
*c_o = geom.c;
*d_o = geom.d;
*dmu_o = local_jet.d1;
*d2_o = local_jet.d2;
*d3_o = local_jet.d3;
Ok(())
},
)?;
} else {
let mu_s = mu.as_slice_mut().expect("mu must be contiguous");
let weights_s = weights.as_slice_mut().expect("weights must be contiguous");
let z_s = z.as_slice_mut().expect("z must be contiguous");
mu_s.par_iter_mut()
.zip(weights_s.par_iter_mut())
.zip(z_s.par_iter_mut())
.enumerate()
.try_for_each(|(i, ((mu_o, w_o), z_o))| -> Result<(), EstimationError> {
let jet = if let InverseLink::LatentCLogLog(state) = inverse_link {
crate::families::lognormal_kernel::latent_cloglog_inverse_link_jet(
quadctx,
eta[i],
se[i].hypot(state.latent_sd),
)?
} else if matches!(inverse_link, InverseLink::Standard(StandardLink::Logit)) {
crate::quadrature::integrated_logit_inverse_link_jet_pirls(
quadctx, eta[i], se[i],
)?
} else {
crate::quadrature::integrated_inverse_link_jetwith_state(
quadctx,
link,
eta[i],
se[i],
inverse_link.mixture_state(),
inverse_link.sas_state(),
)?
};
let local_jet = MixtureInverseLinkJet {
mu: jet.mean,
d1: jet.d1,
d2: jet.d2,
d3: jet.d3,
};
let e = eta[i].clamp(-700.0, 700.0);
let geom = bernoulli_geometry_from_jet(eta[i], e, y[i], priorweights[i], local_jet);
*mu_o = geom.mu;
*w_o = geom.weight;
*z_o = geom.z;
Ok(())
})?;
}
Ok(())
}
#[inline]
pub fn update_glmvectors_integrated_by_family(
quadctx: &crate::quadrature::QuadratureContext,
y: ArrayView1<f64>,
eta: &Array1<f64>,
se: ArrayView1<f64>,
spec: &LikelihoodSpec,
priorweights: ArrayView1<f64>,
mu: &mut Array1<f64>,
weights: &mut Array1<f64>,
z: &mut Array1<f64>,
derivatives: Option<WorkingDerivativeBuffersMut<'_>>,
mixture_link_state: Option<&MixtureLinkState>,
sas_link_state: Option<&SasLinkState>,
) -> Result<(), EstimationError> {
let inverse_link =
integrated_inverse_link_from_family(spec, mixture_link_state, sas_link_state)?;
update_glmvectors_integrated_for_link(
quadctx,
y,
eta,
se,
&inverse_link,
priorweights,
mu,
weights,
z,
derivatives,
)
}
pub(crate) fn computeworkingweight_derivatives_from_eta(
likelihood: &GlmLikelihoodSpec,
inverse_link: &InverseLink,
eta: &Array1<f64>,
priorweights: ArrayView1<f64>,
) -> Result<
(
Array1<f64>,
Array1<f64>,
Array1<f64>,
Array1<f64>,
Array1<f64>,
),
EstimationError,
> {
let n = eta.len();
let mut c = Array1::<f64>::zeros(n);
let mut d = Array1::<f64>::zeros(n);
let mut dmu_deta = Array1::<f64>::zeros(n);
let mut d2mu_deta2 = Array1::<f64>::zeros(n);
let mut d3mu_deta3 = Array1::<f64>::zeros(n);
match &likelihood.spec.response {
ResponseFamily::Gaussian => {
dmu_deta.fill(1.0);
}
ResponseFamily::Poisson => {
const MIN_WEIGHT: f64 = 1e-12;
let c_s = c.as_slice_mut().expect("c must be contiguous");
let d_s = d.as_slice_mut().expect("d must be contiguous");
let dmu_s = dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
c_s.par_iter_mut()
.zip(d_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.try_for_each(
|(i, ((((c_o, d_o), dmu_o), d2_o), d3_o))| -> Result<(), EstimationError> {
let eta_used = eta[i].clamp(-700.0, 700.0);
let jet = standard_inverse_link_jet(inverse_link, eta_used)?;
let raw_weight = priorweights[i].max(0.0) * jet.mu;
let floor_active = raw_weight > 0.0 && raw_weight <= MIN_WEIGHT;
if eta[i] != eta_used || floor_active {
*c_o = 0.0;
*d_o = 0.0;
} else {
*c_o = raw_weight;
*d_o = raw_weight;
}
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
Ok(())
},
)?;
}
ResponseFamily::Tweedie { p } => {
let p = *p;
const MIN_WEIGHT: f64 = 1e-12;
if !is_valid_tweedie_power(p) {
crate::bail_invalid_estim!(
"Tweedie variance power must be finite and strictly between 1 and 2; got {p}"
);
}
let exponent = 2.0 - p;
let phi = fixed_glm_dispersion(likelihood);
if !(phi.is_finite() && phi > 0.0) {
crate::bail_invalid_estim!(
"Tweedie dispersion phi must be finite and > 0; got {phi}"
);
}
let c_s = c.as_slice_mut().expect("c must be contiguous");
let d_s = d.as_slice_mut().expect("d must be contiguous");
let dmu_s = dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
c_s.par_iter_mut()
.zip(d_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.try_for_each(
|(i, ((((c_o, d_o), dmu_o), d2_o), d3_o))| -> Result<(), EstimationError> {
let eta_raw = eta[i];
let eta_used = eta_raw.clamp(-700.0, 700.0);
let clamp_active = eta_raw != eta_used;
let jet = standard_inverse_link_jet(inverse_link, eta_used)?;
let raw_weight =
priorweights[i].max(0.0) * tweedie_log_weight_mu_power(jet.mu, p) / phi;
let floor_active = raw_weight > 0.0 && raw_weight <= MIN_WEIGHT;
if clamp_active || floor_active {
*c_o = 0.0;
*d_o = 0.0;
} else {
*c_o = exponent * raw_weight;
*d_o = exponent * exponent * raw_weight;
}
if clamp_active {
*dmu_o = 0.0;
*d2_o = 0.0;
*d3_o = 0.0;
} else {
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
}
Ok(())
},
)?;
}
ResponseFamily::NegativeBinomial { theta } => {
let theta = *theta;
const MIN_WEIGHT: f64 = 1e-12;
if !valid_negbin_theta(theta) {
crate::bail_invalid_estim!(
"negative-binomial theta must be finite and > 0; got {theta}"
);
}
let c_s = c.as_slice_mut().expect("c must be contiguous");
let d_s = d.as_slice_mut().expect("d must be contiguous");
let dmu_s = dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
c_s.par_iter_mut()
.zip(d_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.try_for_each(
|(i, ((((c_o, d_o), dmu_o), d2_o), d3_o))| -> Result<(), EstimationError> {
let eta_raw = eta[i];
let eta_used = eta_raw.clamp(-700.0, 700.0);
let jet = standard_inverse_link_jet(inverse_link, eta_used)?;
let denom = theta + jet.mu;
let negbin_weight = if theta > jet.mu {
jet.mu / (1.0 + jet.mu / theta)
} else {
theta / (1.0 + theta / jet.mu)
};
let raw_weight = priorweights[i].max(0.0) * negbin_weight;
let floor_active = raw_weight > 0.0 && raw_weight <= MIN_WEIGHT;
if eta_raw != eta_used || floor_active {
*c_o = 0.0;
*d_o = 0.0;
} else {
*c_o = raw_weight * theta / denom;
*d_o = raw_weight * theta * (theta - jet.mu) / (denom * denom);
}
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
Ok(())
},
)?;
}
ResponseFamily::Beta { phi } => {
let phi = *phi;
const MIN_WEIGHT: f64 = 1e-12;
if !valid_beta_phi(phi) {
crate::bail_invalid_estim!("beta-regression phi must be finite and > 0; got {phi}");
}
let c_s = c.as_slice_mut().expect("c must be contiguous");
let d_s = d.as_slice_mut().expect("d must be contiguous");
let dmu_s = dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
c_s.par_iter_mut()
.zip(d_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.for_each(|(i, ((((c_o, d_o), dmu_o), d2_o), d3_o))| {
let eta_raw = eta[i];
let eta_i = eta_raw.clamp(-700.0, 700.0);
let jet = logit_inverse_link_jet5(eta_i);
let mu_i = safe_beta_mu(jet.mu);
let q = (mu_i * (1.0 - mu_i)).max(BETA_MU_EPS);
let a = (mu_i * phi).max(BETA_MU_EPS);
let b = ((1.0 - mu_i) * phi).max(BETA_MU_EPS);
let trigamma_sum = trigamma(a) + trigamma(b);
let prior_weight = priorweights[i].max(0.0);
let raw_weight = prior_weight * q * q * phi * phi * trigamma_sum;
let floor_active = raw_weight > 0.0 && raw_weight <= MIN_WEIGHT;
if floor_active || eta_raw != eta_i {
*c_o = 0.0;
*d_o = 0.0;
} else {
let (c_i, d_i) = beta_logit_working_curvature_eta_derivatives(
prior_weight,
phi,
mu_i,
q,
a,
b,
trigamma_sum,
);
*c_o = c_i;
*d_o = d_i;
}
*dmu_o = q;
*d2_o = q * (1.0 - 2.0 * mu_i);
*d3_o = q * (1.0 - 6.0 * q);
});
}
ResponseFamily::Gamma => {
let dmu_s = dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
dmu_s
.par_iter_mut()
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.try_for_each(
|(i, ((dmu_o, d2_o), d3_o))| -> Result<(), EstimationError> {
let jet =
standard_inverse_link_jet(inverse_link, eta[i].clamp(-700.0, 700.0))?;
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
Ok(())
},
)?;
}
ResponseFamily::Binomial => {
let link = inverse_link.link_function();
let zero_on_nonsmooth = matches!(link, LinkFunction::Logit);
let c_s = c.as_slice_mut().expect("c must be contiguous");
let d_s = d.as_slice_mut().expect("d must be contiguous");
let dmu_s = dmu_deta
.as_slice_mut()
.expect("dmu_deta must be contiguous");
let d2_s = d2mu_deta2
.as_slice_mut()
.expect("d2mu_deta2 must be contiguous");
let d3_s = d3mu_deta3
.as_slice_mut()
.expect("d3mu_deta3 must be contiguous");
c_s.par_iter_mut()
.zip(d_s.par_iter_mut())
.zip(dmu_s.par_iter_mut())
.zip(d2_s.par_iter_mut())
.zip(d3_s.par_iter_mut())
.enumerate()
.try_for_each(
|(i, ((((c_o, d_o), dmu_o), d2_o), d3_o))| -> Result<(), EstimationError> {
let eta_used = match link {
LinkFunction::Logit => eta[i].clamp(-700.0, 700.0),
LinkFunction::Probit
| LinkFunction::CLogLog
| LinkFunction::Sas
| LinkFunction::BetaLogistic => eta[i].clamp(-30.0, 30.0),
LinkFunction::Log => eta[i].clamp(-700.0, 700.0),
LinkFunction::Identity => eta[i],
};
if matches!(link, LinkFunction::Logit) {
let jet = logit_inverse_link_jet5(eta_used);
let geom = bernoulli_logit_geometry_from_jet(
eta[i],
eta_used,
jet.mu,
priorweights[i],
jet,
zero_on_nonsmooth,
);
*c_o = geom.c;
*d_o = geom.d;
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
} else {
let jet = standard_inverse_link_jet(inverse_link, eta_used)?;
let geom = bernoulli_geometry_from_jet(
eta[i],
eta_used,
jet.mu,
priorweights[i],
jet,
);
*c_o = geom.c;
*d_o = geom.d;
*dmu_o = jet.d1;
*d2_o = jet.d2;
*d3_o = jet.d3;
}
Ok(())
},
)?;
}
ResponseFamily::RoystonParmar => {
crate::bail_invalid_estim!(
"RoystonParmar is survival-specific and not a GLM IRLS family"
);
}
}
Ok((c, d, dmu_deta, d2mu_deta2, d3mu_deta3))
}
#[derive(Clone, Copy, Debug)]
pub struct VarianceJet {
pub v: f64,
pub v1: f64,
pub v2: f64,
pub v3: f64,
pub v4: f64,
}
impl VarianceJet {
#[inline]
pub fn bernoulli(mu: f64) -> Self {
Self {
v: mu * (1.0 - mu),
v1: 1.0 - 2.0 * mu,
v2: -2.0,
v3: 0.0,
v4: 0.0,
}
}
#[inline]
pub fn poisson(mu: f64) -> Self {
Self {
v: mu,
v1: 1.0,
v2: 0.0,
v3: 0.0,
v4: 0.0,
}
}
#[inline]
pub fn gamma(mu: f64) -> Self {
Self {
v: mu * mu,
v1: 2.0 * mu,
v2: 2.0,
v3: 0.0,
v4: 0.0,
}
}
#[inline]
pub fn tweedie(mu: f64, p: f64) -> Self {
let mu = mu.max(1e-10);
Self {
v: mu.powf(p),
v1: p * mu.powf(p - 1.0),
v2: p * (p - 1.0) * mu.powf(p - 2.0),
v3: p * (p - 1.0) * (p - 2.0) * mu.powf(p - 3.0),
v4: p * (p - 1.0) * (p - 2.0) * (p - 3.0) * mu.powf(p - 4.0),
}
}
#[inline]
pub fn negative_binomial(mu: f64, theta: f64) -> Self {
let mu = mu.max(1e-10);
let inv_theta = if valid_negbin_theta(theta) {
1.0 / theta
} else {
f64::NAN
};
Self {
v: mu + mu * mu * inv_theta,
v1: 1.0 + 2.0 * mu * inv_theta,
v2: 2.0 * inv_theta,
v3: 0.0,
v4: 0.0,
}
}
#[inline]
pub fn gaussian() -> Self {
Self {
v: 1.0,
v1: 0.0,
v2: 0.0,
v3: 0.0,
v4: 0.0,
}
}
#[inline]
pub fn binomial_n(mu: f64) -> Self {
Self::bernoulli(mu)
}
#[inline]
pub fn beta(mu: f64, phi: f64) -> Self {
let scale = 1.0 / (1.0 + phi.max(1e-12));
let base = Self::bernoulli(mu);
Self {
v: base.v * scale,
v1: base.v1 * scale,
v2: base.v2 * scale,
v3: 0.0,
v4: 0.0,
}
}
}
const OBSERVED_HESSIAN_WEIGHT_FLOOR_FRAC: f64 = 1e-6;
const OBSERVED_HESSIAN_WEIGHT_ABS_FLOOR: f64 = 1e-12;
#[inline]
pub fn solver_hessian_weight_floor(fisher_weight: f64) -> f64 {
(fisher_weight.max(0.0) * OBSERVED_HESSIAN_WEIGHT_FLOOR_FRAC)
.max(OBSERVED_HESSIAN_WEIGHT_ABS_FLOOR)
}
pub fn outer_hessian_curvature_arrays(
hessian_weights: crate::matrix::SignedWeightsView<'_>,
fisher_weights: crate::matrix::PsdWeightsView<'_>,
c_array: &Array1<f64>,
d_array: &Array1<f64>,
eta: &Array1<f64>,
inverse_link: &InverseLink,
) -> (Array1<f64>, Array1<f64>, Array1<f64>) {
let hessian_view = hessian_weights.view();
let fisher_view = fisher_weights.view();
let n = hessian_view.len();
let mut w_out = Array1::<f64>::zeros(n);
let mut c_out = Array1::<f64>::zeros(n);
let mut d_out = Array1::<f64>::zeros(n);
for i in 0..n {
let floor = solver_hessian_weight_floor(fisher_view[i]);
let w = hessian_view[i];
let clamp_active = eta_clamp_active(inverse_link, eta[i]);
let w_below_floor = !(w.is_finite() && w > floor);
if w_below_floor {
w_out[i] = floor;
c_out[i] = 0.0;
d_out[i] = 0.0;
} else if clamp_active {
w_out[i] = w;
c_out[i] = 0.0;
d_out[i] = 0.0;
} else {
w_out[i] = w;
c_out[i] = c_array[i];
d_out[i] = d_array[i];
}
}
(w_out, c_out, d_out)
}
#[inline]
fn fixed_glm_dispersion(likelihood: &GlmLikelihoodSpec) -> f64 {
likelihood.fixed_phi().unwrap_or(1.0)
}
#[inline]
pub fn weight_family_for_glm_likelihood(likelihood: &GlmLikelihoodSpec) -> WeightFamily {
match &likelihood.spec.response {
ResponseFamily::Gaussian => WeightFamily::Gaussian,
ResponseFamily::Poisson => WeightFamily::Poisson,
ResponseFamily::Tweedie { p } => WeightFamily::Tweedie { p: *p },
ResponseFamily::NegativeBinomial { theta } => {
WeightFamily::NegativeBinomial { theta: *theta }
}
ResponseFamily::Beta { phi } => WeightFamily::Beta { phi: *phi },
ResponseFamily::Gamma => WeightFamily::Gamma,
ResponseFamily::Binomial => WeightFamily::Binomial,
ResponseFamily::RoystonParmar => WeightFamily::Gaussian,
}
}
#[inline]
fn weight_link_for_inverse_link(inverse_link: &InverseLink) -> WeightLink {
match inverse_link {
InverseLink::Standard(StandardLink::Identity) => WeightLink::Identity,
InverseLink::Standard(StandardLink::Log) => WeightLink::Log,
InverseLink::Standard(StandardLink::Logit) => WeightLink::Logit,
InverseLink::Standard(StandardLink::Probit)
| InverseLink::Standard(StandardLink::CLogLog)
| InverseLink::LatentCLogLog(_)
| InverseLink::Sas(_)
| InverseLink::BetaLogistic(_)
| InverseLink::Mixture(_) => WeightLink::Other,
}
}
#[inline]
fn supports_observed_hessian_curvature_for_likelihood(
likelihood: &GlmLikelihoodSpec,
inverse_link: &InverseLink,
) -> bool {
let spec = &likelihood.spec;
if matches!(spec.response, ResponseFamily::NegativeBinomial { .. }) {
return matches!(inverse_link, InverseLink::Standard(StandardLink::Log));
}
if matches!(spec.response, ResponseFamily::Gamma) {
return true;
}
if !matches!(spec.response, ResponseFamily::Binomial) {
return false;
}
matches!(
spec.link,
InverseLink::Standard(StandardLink::Probit)
| InverseLink::Standard(StandardLink::CLogLog)
| InverseLink::Sas(_)
| InverseLink::BetaLogistic(_)
| InverseLink::Mixture(_)
)
}
#[inline]
fn eta_for_observed_hessian_jet(inverse_link: &InverseLink, eta: f64) -> f64 {
match inverse_link {
InverseLink::Standard(StandardLink::Logit | StandardLink::Log) => eta.clamp(-700.0, 700.0),
InverseLink::Standard(StandardLink::Identity) => eta,
InverseLink::Standard(StandardLink::Probit) => eta.clamp(-6.0, 6.0),
InverseLink::Standard(StandardLink::CLogLog) | InverseLink::LatentCLogLog(_) => {
eta.clamp(-23.0, 3.0)
}
InverseLink::Sas(_) | InverseLink::BetaLogistic(_) | InverseLink::Mixture(_) => {
eta.clamp(-20.0, 20.0)
}
}
}
#[inline]
pub fn eta_clamp_active(inverse_link: &InverseLink, eta: f64) -> bool {
let clamped = eta_for_observed_hessian_jet(inverse_link, eta);
clamped != eta
}
fn solver_hessian_weights_into(
hessian_weights: &Array1<f64>,
fisher_weights: &Array1<f64>,
out: &mut Array1<f64>,
) {
if out.len() != hessian_weights.len() {
*out = Array1::<f64>::zeros(hessian_weights.len());
}
ndarray::Zip::from(out)
.and(hessian_weights)
.and(fisher_weights)
.par_for_each(|o, &w, &fw| {
let floor = solver_hessian_weight_floor(fw);
*o = if w.is_finite() && w > floor { w } else { floor };
});
}
fn compute_observed_hessian_curvature_arrays_into(
likelihood: &GlmLikelihoodSpec,
inverse_link: &InverseLink,
eta: &Array1<f64>,
y: ArrayView1<'_, f64>,
fisher_weights: &Array1<f64>,
priorweights: ArrayView1<'_, f64>,
hessian_weights: &mut Array1<f64>,
hessian_c: &mut Array1<f64>,
hessian_d: &mut Array1<f64>,
) -> Result<(), EstimationError> {
assert!(supports_observed_hessian_curvature_for_likelihood(
likelihood,
inverse_link
));
let n = eta.len();
if hessian_weights.len() != n {
*hessian_weights = Array1::<f64>::zeros(n);
}
if hessian_c.len() != n {
*hessian_c = Array1::<f64>::zeros(n);
}
if hessian_d.len() != n {
*hessian_d = Array1::<f64>::zeros(n);
}
let weight_family = weight_family_for_glm_likelihood(likelihood);
let weight_link = weight_link_for_inverse_link(inverse_link);
let phi = fixed_glm_dispersion(likelihood);
hessian_weights
.as_slice_mut()
.expect("hessian weights must be contiguous")
.par_iter_mut()
.zip(
hessian_c
.as_slice_mut()
.expect("hessian c must be contiguous")
.par_iter_mut(),
)
.zip(
hessian_d
.as_slice_mut()
.expect("hessian d must be contiguous")
.par_iter_mut(),
)
.enumerate()
.try_for_each(|(i, ((w_out, c_out), d_out))| -> Result<(), EstimationError> {
let eta_used = eta_for_observed_hessian_jet(inverse_link, eta[i]);
let jet =
crate::mixture_link::inverse_link_jet_for_inverse_link(inverse_link, eta_used)?;
let h4 = crate::mixture_link::inverse_link_pdfthird_derivative_for_inverse_link(
inverse_link, eta_used,
)?;
let (w_obs, c_obs, d_obs) = observed_weight_dispatch(
weight_family,
weight_link,
eta_used,
y[i],
jet.mu,
phi,
priorweights[i].max(0.0),
jet,
h4,
);
let fisher_weight = fisher_weights[i].max(0.0);
if !(w_obs.is_finite() && w_obs > 0.0) {
crate::bail_invalid_estim!(
"observed Hessian curvature is not positive finite at row {i}: observed={w_obs}, fisher={fisher_weight}"
);
}
if !c_obs.is_finite() || !d_obs.is_finite() {
crate::bail_invalid_estim!(
"observed Hessian curvature derivatives are non-finite at row {i}: c={c_obs}, d={d_obs}"
);
}
*w_out = w_obs;
*c_out = c_obs;
*d_out = d_obs;
Ok(())
})
}
pub(crate) fn compute_observed_hessian_curvature_arrays(
likelihood: &GlmLikelihoodSpec,
inverse_link: &InverseLink,
eta: &Array1<f64>,
y: ArrayView1<'_, f64>,
fisher_weights: &Array1<f64>,
priorweights: ArrayView1<'_, f64>,
) -> Result<(Array1<f64>, Array1<f64>, Array1<f64>), EstimationError> {
let n = eta.len();
let mut hessian_weights = Array1::<f64>::zeros(n);
let mut hessian_c = Array1::<f64>::zeros(n);
let mut hessian_d = Array1::<f64>::zeros(n);
compute_observed_hessian_curvature_arrays_into(
likelihood,
inverse_link,
eta,
y,
fisher_weights,
priorweights,
&mut hessian_weights,
&mut hessian_c,
&mut hessian_d,
)?;
Ok((hessian_weights, hessian_c, hessian_d))
}
#[inline]
pub fn observed_weight_noncanonical(
y: f64,
mu: f64,
h1: f64,
h2: f64,
h3: f64,
h4: f64,
vj: VarianceJet,
phi: f64,
pw: f64,
) -> (f64, f64, f64) {
let VarianceJet {
v,
v1,
v2,
v3,
v4: _,
} = vj;
let phi_v = phi * v;
let phi_v2 = phi * v * v;
let phi_v3 = phi * v * v * v;
let h1_sq = h1 * h1;
let w_f = h1_sq / phi_v;
let n0 = h1_sq; let n1 = 2.0 * h1 * h2; let n2 = 2.0 * (h2 * h2 + h1 * h3); let vd1 = h1 * v1; let vd2 = h2 * v1 + h1_sq * v2;
let c_f = (n1 * v - n0 * vd1) / phi_v2;
let numer_cf = n1 * v - n0 * vd1;
let dnumer_cf = n2 * v - n0 * vd2;
let d_f = (dnumer_cf * v - 2.0 * numer_cf * vd1) / (phi_v3);
let b_num = h2 * v - h1_sq * v1;
let b = b_num / phi_v2;
let b_eta_num =
h3 * v * v - 3.0 * h1 * h2 * v * v1 - h1_sq * h1 * v * v2 + 2.0 * h1_sq * h1 * v1 * v1;
let b_eta = b_eta_num / phi_v3;
let h1_cu = h1_sq * h1;
let h1_qu = h1_sq * h1_sq;
let db_eta_num = h4 * v * v + 2.0 * h3 * v * h1 * v1
- 3.0 * (h2 * h2 + h1 * h3) * v * v1
- 3.0 * h1 * h2 * (h1 * v1 * v1 + v * h1 * v2)
- 3.0 * h1_sq * h2 * v * v2
- h1_cu * (h1 * v1 * v2 + v * h1 * v3)
+ 6.0 * h1_sq * h2 * v1 * v1
+ 4.0 * h1_qu * v1 * v2;
let phi_v4 = phi_v3 * v;
let b_etaeta = (db_eta_num * v - 3.0 * b_eta_num * h1 * v1) / phi_v4;
let resid = y - mu;
let w_obs = w_f - resid * b;
let c_obs = c_f + h1 * b - resid * b_eta;
let d_obs = d_f + h2 * b + 2.0 * h1 * b_eta - resid * b_etaeta;
(pw * w_obs, pw * c_obs, pw * d_obs)
}
#[inline]
pub fn e_obs_from_jets(
y: f64,
mu: f64,
h1: f64,
h2: f64,
h3: f64,
h4: f64,
h5: f64,
vj: VarianceJet,
phi: f64,
pw: f64,
) -> f64 {
let VarianceJet { v, v1, v2, v3, v4 } = vj;
let q = phi * v;
let h1_sq = h1 * h1;
let h1_cu = h1_sq * h1;
let h1_qu = h1_sq * h1_sq;
let q1 = phi * v1 * h1;
let q2 = phi * (v1 * h2 + v2 * h1_sq);
let q3 = phi * (v1 * h3 + 3.0 * v2 * h1 * h2 + v3 * h1_cu);
let q4 = phi
* (v1 * h4 + 4.0 * v2 * h1 * h3 + 3.0 * v2 * h2 * h2 + 6.0 * v3 * h1_sq * h2 + v4 * h1_qu);
let t0 = h1 / q;
let t1 = (h2 - t0 * q1) / q;
let t2 = (h3 - 2.0 * t1 * q1 - t0 * q2) / q;
let t3 = (h4 - 3.0 * t2 * q1 - 3.0 * t1 * q2 - t0 * q3) / q;
let t4 = (h5 - 4.0 * t3 * q1 - 6.0 * t2 * q2 - 4.0 * t1 * q3 - t0 * q4) / q;
let w_f3 = h1 * t3 + 3.0 * h2 * t2 + 3.0 * h3 * t1 + h4 * t0;
let resid = y - mu;
let e_obs = w_f3 + h3 * t1 + 3.0 * h2 * t2 + 3.0 * h1 * t3 - resid * t4;
pw * e_obs
}
#[inline]
pub fn observed_weight_gaussian_log(y: f64, mu: f64, phi: f64, pw: f64) -> (f64, f64, f64) {
let inv_phi = pw / phi;
let w = inv_phi * mu * (2.0 * mu - y);
let c = inv_phi * mu * (4.0 * mu - y);
let d = inv_phi * mu * (8.0 * mu - y);
(w, c, d)
}
#[inline]
pub fn observed_weight_gaussian_inverse(y: f64, eta: f64, phi: f64, pw: f64) -> (f64, f64, f64) {
let eta2 = eta * eta;
let eta4 = eta2 * eta2;
let eta5 = eta4 * eta;
let eta6 = eta4 * eta2;
let ey = eta * y;
let inv_phi = pw / phi;
let w = inv_phi * (3.0 - 2.0 * ey) / eta4;
let c = inv_phi * 6.0 * (ey - 2.0) / eta5;
let d = inv_phi * 12.0 * (5.0 - 2.0 * ey) / eta6;
(w, c, d)
}
#[inline]
fn observed_weight_binomial_logit_from_jet(
n_trials: f64,
jet: MixtureInverseLinkJet,
pw: f64,
) -> (f64, f64, f64) {
let scale = pw * n_trials;
(scale * jet.d1, scale * jet.d2, scale * jet.d3)
}
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum WeightFamily {
Gaussian,
Binomial,
Poisson,
Tweedie { p: f64 },
NegativeBinomial { theta: f64 },
Beta { phi: f64 },
Gamma,
}
#[derive(Debug, Clone, Copy, PartialEq, Eq)]
pub enum WeightLink {
Identity,
Log,
Logit,
Inverse,
Other,
}
#[inline]
pub fn variance_jet_for_weight_family(family: WeightFamily, mu: f64) -> VarianceJet {
match family {
WeightFamily::Gaussian => VarianceJet::gaussian(),
WeightFamily::Binomial => VarianceJet::binomial_n(mu),
WeightFamily::Poisson => VarianceJet::poisson(mu),
WeightFamily::Tweedie { p } => VarianceJet::tweedie(mu, p),
WeightFamily::NegativeBinomial { theta } => VarianceJet::negative_binomial(mu, theta),
WeightFamily::Beta { phi } => VarianceJet::beta(mu, phi),
WeightFamily::Gamma => VarianceJet::gamma(mu),
}
}
pub fn observed_weight_dispatch(
family: WeightFamily,
link: WeightLink,
eta: f64,
y: f64,
mu: f64,
phi: f64,
prior_weight: f64,
jet: MixtureInverseLinkJet,
h4: f64,
) -> (f64, f64, f64) {
match (family, link) {
(WeightFamily::Gaussian, WeightLink::Log) => {
observed_weight_gaussian_log(y, mu, phi, prior_weight)
}
(WeightFamily::Gaussian, WeightLink::Inverse) => {
observed_weight_gaussian_inverse(y, eta, phi, prior_weight)
}
(WeightFamily::Binomial, WeightLink::Logit) => {
observed_weight_binomial_logit_from_jet(1.0, jet, prior_weight)
}
_ => {
let vj = variance_jet_for_weight_family(family, mu);
observed_weight_noncanonical(y, mu, jet.d1, jet.d2, jet.d3, h4, vj, phi, prior_weight)
}
}
}
#[derive(Clone)]
pub enum DirectionalWorkingCurvature {
Diagonal(Array1<f64>),
}
pub fn directionalworking_curvature_from_c_array(
c_array: &Array1<f64>,
hessian_weights: &Array1<f64>,
eta_direction: &Array1<f64>,
) -> DirectionalWorkingCurvature {
let mut w_direction = c_array * eta_direction;
for i in 0..w_direction.len() {
if hessian_weights[i] <= 0.0 || !w_direction[i].is_finite() {
w_direction[i] = 0.0;
}
}
DirectionalWorkingCurvature::Diagonal(w_direction)
}
const BINOMIAL_MU_EPS: f64 = 1e-12;
#[inline]
fn safe_mu_for_binomial(mu: f64) -> f64 {
mu.clamp(BINOMIAL_MU_EPS, 1.0 - BINOMIAL_MU_EPS)
}
#[inline]
fn xlogy(x: f64, y: f64) -> f64 {
if x == 0.0 { 0.0 } else { x * y.ln() }
}
#[inline]
fn log_gamma_stirling_correction(x: f64) -> f64 {
let inv = 1.0 / x;
let inv2 = inv * inv;
inv / 12.0 - inv * inv2 / 360.0 + inv * inv2 * inv2 / 1260.0
}
#[inline]
fn log_gamma_large_ratio(base: f64, delta: f64) -> f64 {
let ratio = delta / base;
delta * base.ln() + (base + delta - 0.5) * ratio.ln_1p() - delta
+ log_gamma_stirling_correction(base + delta)
- log_gamma_stirling_correction(base)
}
#[inline]
fn beta_log_normalizer(a: f64, b: f64, sum: f64) -> f64 {
let direct = ln_gamma(sum) - ln_gamma(a) - ln_gamma(b);
if direct.is_finite() {
return direct;
}
let small = a.min(b);
let large = a.max(b);
if small < 8.0 {
return log_gamma_large_ratio(large, small) - ln_gamma(small);
}
-xlogy(a, a / sum) - xlogy(b, b / sum)
+ 0.5 * (a.ln() + b.ln() - sum.ln() - (2.0 * std::f64::consts::PI).ln())
+ log_gamma_stirling_correction(sum)
- log_gamma_stirling_correction(a)
- log_gamma_stirling_correction(b)
}
#[inline]
fn poisson_unit_deviance(yi: f64, mui_c: f64) -> f64 {
xlogy(yi, yi / mui_c) - (yi - mui_c)
}
#[inline]
fn gamma_unit_deviance(yi_c: f64, mui_c: f64) -> f64 {
let ratio = yi_c / mui_c;
ratio - 1.0 - ratio.ln()
}
#[inline]
fn tweedie_unit_deviance(yi: f64, mui_c: f64, p: f64) -> f64 {
if !is_valid_tweedie_power(p) {
f64::NAN
} else if !valid_tweedie_response(yi) {
f64::NAN
} else if yi == 0.0 {
mui_c.powf(2.0 - p) / (2.0 - p)
} else {
yi.powf(2.0 - p) / ((1.0 - p) * (2.0 - p)) - yi * mui_c.powf(1.0 - p) / (1.0 - p)
+ mui_c.powf(2.0 - p) / (2.0 - p)
}
}
#[inline]
fn negative_binomial_unit_deviance(yi: f64, mui_c: f64, theta: f64) -> f64 {
if !valid_negbin_theta(theta) || !valid_count_response(yi) {
return f64::NAN;
}
let y_term = xlogy(yi, (yi * (theta + mui_c)) / (mui_c * (theta + yi)));
let theta_term = theta * ((theta + mui_c) / (theta + yi)).ln();
theta_term + y_term
}
#[inline]
fn beta_loglikelihood_full_unit(yi: f64, mui: f64, phi: f64) -> f64 {
if !valid_beta_phi(phi) || !valid_beta_response(yi) {
return f64::NAN;
}
let mui_c = safe_beta_mu(mui);
let a = (mui_c * phi).max(BETA_MU_EPS);
let b = ((1.0 - mui_c) * phi).max(BETA_MU_EPS);
beta_log_normalizer(a, b, phi) + phi * xlogy(mui_c, yi) + phi * xlogy(1.0 - mui_c, 1.0 - yi)
- yi.ln()
- (1.0 - yi).ln()
}
#[inline]
fn beta_unit_deviance(yi: f64, mui: f64, phi: f64) -> f64 {
if !valid_beta_response(yi) {
return f64::NAN;
}
beta_loglikelihood_full_unit(yi, yi, phi) - beta_loglikelihood_full_unit(yi, mui, phi)
}
#[inline]
pub fn calculate_deviance(
y: ArrayView1<f64>,
mu: &Array1<f64>,
likelihood: &GlmLikelihoodSpec,
priorweights: ArrayView1<f64>,
) -> f64 {
const EPS: f64 = 1e-8;
const MU_FLOOR: f64 = 1e-10;
match &likelihood.spec.response {
ResponseFamily::Binomial => {
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let total_residual: f64 = (0..y.len())
.into_par_iter()
.map(|i| {
let yi = y[i];
let mui_c = safe_mu_for_binomial(mu[i]);
let wi = priorweights[i];
let term1 = if yi > EPS {
yi * (yi.ln() - mui_c.ln())
} else {
0.0
};
let term2 = if yi < 1.0 - EPS {
(1.0 - yi) * ((1.0 - yi).ln() - (1.0 - mui_c).ln())
} else {
0.0
};
wi * (term1 + term2)
})
.sum();
2.0 * total_residual
}
ResponseFamily::Gaussian => {
let phi = likelihood.scale.fixed_phi().unwrap_or(1.0);
if !(phi.is_finite() && phi > 0.0) {
return f64::NAN;
}
let raw: f64 = ndarray::Zip::from(y)
.and(mu)
.and(priorweights)
.map_collect(|&yi, &mui, &wi| wi * (yi - mui) * (yi - mui))
.sum();
raw / phi
}
ResponseFamily::Poisson => {
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let total: f64 = (0..y.len())
.into_par_iter()
.map(|i| {
let yi = y[i];
let mui_c = mu[i].max(MU_FLOOR);
priorweights[i] * poisson_unit_deviance(yi, mui_c)
})
.sum();
2.0 * total
}
ResponseFamily::Tweedie { p } => {
let p = *p;
let phi = fixed_glm_dispersion(likelihood);
if !is_valid_tweedie_power(p) || !(phi.is_finite() && phi > 0.0) {
return f64::NAN;
}
if validate_tweedie_responses(&y, &priorweights).is_err() {
return f64::NAN;
}
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let total: f64 = (0..y.len())
.into_par_iter()
.map(|i| {
let yi = y[i];
let mui_c = mu[i].max(MU_FLOOR);
priorweights[i] * tweedie_unit_deviance(yi, mui_c, p) / phi
})
.sum();
2.0 * total
}
ResponseFamily::NegativeBinomial { theta } => {
let theta = *theta;
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let total: f64 = (0..y.len())
.into_par_iter()
.map(|i| {
let yi = y[i];
let mui_c = mu[i].max(MU_FLOOR);
priorweights[i] * negative_binomial_unit_deviance(yi, mui_c, theta)
})
.sum();
2.0 * total
}
ResponseFamily::Beta { phi } => {
let phi = *phi;
if !valid_beta_phi(phi) {
return f64::NAN;
}
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let total: f64 = (0..y.len())
.into_par_iter()
.map(|i| priorweights[i] * beta_unit_deviance(y[i], mu[i], phi))
.sum();
2.0 * total
}
ResponseFamily::Gamma => {
let shape = likelihood.gamma_shape().unwrap_or(1.0);
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let total: f64 = (0..y.len())
.into_par_iter()
.map(|i| {
let yi_c = y[i].max(EPS);
let mui_c = mu[i].max(MU_FLOOR);
priorweights[i] * shape * gamma_unit_deviance(yi_c, mui_c)
})
.sum();
2.0 * total
}
ResponseFamily::RoystonParmar => f64::NAN,
}
}
#[inline]
pub(crate) fn calculate_loglikelihood_omitting_constants(
y: ArrayView1<f64>,
mu: &Array1<f64>,
likelihood: &GlmLikelihoodSpec,
priorweights: ArrayView1<f64>,
) -> f64 {
const MU_FLOOR: f64 = 1e-10;
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let n = y.len();
match &likelihood.spec.response {
ResponseFamily::Gaussian => {
let phi = likelihood.scale.fixed_phi().unwrap_or(1.0);
if !(phi.is_finite() && phi > 0.0) {
return f64::NAN;
}
let inv_phi = 1.0 / phi;
(0..n)
.into_par_iter()
.map(|i| {
let resid = y[i] - mu[i];
-0.5 * priorweights[i] * resid * resid * inv_phi
})
.sum()
}
ResponseFamily::Binomial => (0..n)
.into_par_iter()
.map(|i| {
let mui_c = safe_mu_for_binomial(mu[i]);
priorweights[i] * (y[i] * mui_c.ln() + (1.0 - y[i]) * (1.0 - mui_c).ln())
})
.sum(),
ResponseFamily::Poisson => (0..n)
.into_par_iter()
.map(|i| {
let mui_c = mu[i].max(MU_FLOOR);
let log_term = if y[i] > 0.0 { y[i] * mui_c.ln() } else { 0.0 };
priorweights[i] * (log_term - mui_c)
})
.sum(),
ResponseFamily::Tweedie { p } => {
let p = *p;
let phi = fixed_glm_dispersion(likelihood);
if !is_valid_tweedie_power(p) || !(phi.is_finite() && phi > 0.0) {
return f64::NAN;
}
-0.5 * calculate_deviance(y, mu, likelihood, priorweights)
}
ResponseFamily::NegativeBinomial { theta } => {
let theta = *theta;
(0..n)
.into_par_iter()
.map(|i| {
if !valid_negbin_theta(theta) {
return f64::NAN;
}
let yi = y[i];
if !valid_count_response(yi) {
return f64::NAN;
}
let mui_c = mu[i].max(MU_FLOOR);
priorweights[i]
* (ln_gamma(yi + theta) - ln_gamma(theta) - ln_gamma(yi + 1.0)
+ theta * (theta.ln() - (theta + mui_c).ln())
+ xlogy(yi, mui_c)
- yi * (theta + mui_c).ln())
})
.sum()
}
ResponseFamily::Beta { phi } => {
let phi = *phi;
(0..n)
.into_par_iter()
.map(|i| {
if !valid_beta_phi(phi) {
return f64::NAN;
}
priorweights[i] * beta_loglikelihood_full_unit(y[i], mu[i], phi)
})
.sum()
}
ResponseFamily::Gamma => {
-0.5 * calculate_deviance(y, mu, likelihood, priorweights)
}
ResponseFamily::RoystonParmar => f64::NAN,
}
}
pub fn latent_predict_and_apply_ift_warm_start(
cache: &crate::solver::arrow_schur::ArrowFactorCache,
delta_beta: Option<ndarray::ArrayView1<'_, f64>>,
delta_gt: Option<ndarray::ArrayView1<'_, f64>>,
latent: &mut crate::terms::latent_coord::LatentCoordValues,
max_row_delta: f64,
) -> (usize, bool) {
use crate::solver::persistent_warm_start::{LatentIftOutcome, ift_warm_start_latent};
match ift_warm_start_latent(cache, delta_beta, delta_gt) {
LatentIftOutcome::Applied { delta_t } => {
let clamped = latent_apply_ift_warm_start(latent, &delta_t, max_row_delta);
(clamped, true)
}
LatentIftOutcome::Noop => (0, false),
}
}
pub fn latent_apply_ift_warm_start(
latent: &mut crate::terms::latent_coord::LatentCoordValues,
delta_t: &ndarray::Array1<f64>,
max_row_delta: f64,
) -> usize {
let n = latent.n_obs();
let d = latent.latent_dim();
assert_eq!(delta_t.len(), n * d);
let mut clamped_rows = 0_usize;
let mut applied = ndarray::Array1::<f64>::zeros(n * d);
for i in 0..n {
let mut row_norm_sq = 0.0_f64;
for a in 0..d {
let dv = delta_t[i * d + a];
row_norm_sq += dv * dv;
}
let row_norm = row_norm_sq.sqrt();
let scale = if row_norm > max_row_delta && row_norm > 0.0 {
clamped_rows += 1;
max_row_delta / row_norm
} else {
1.0
};
for a in 0..d {
applied[i * d + a] = scale * delta_t[i * d + a];
}
}
latent.retract_flat_delta(applied.view());
clamped_rows
}
use crate::linalg::low_rank_weight::LowRankWeight;
pub fn compute_xtwx_low_rank(
workspace: &mut PirlsWorkspace,
design: &DesignMatrix,
weight: &LowRankWeight<'_>,
) -> Result<Array2<f64>, EstimationError> {
let diag_owned = weight.diag.to_owned();
let mut xtwx = GamWorkingModel::compute_xtwx_blas(workspace, design, &diag_owned)?;
if weight.is_rank_zero() {
return Ok(xtwx);
}
weight
.add_low_rank_xtwx_correction(design, &mut xtwx)
.map_err(EstimationError::InvalidInput)?;
Ok(xtwx)
}
pub fn compute_xtwy_low_rank(
design: &DesignMatrix,
weight: &LowRankWeight<'_>,
y: &Array1<f64>,
) -> Result<Array1<f64>, EstimationError> {
weight
.xtw_y(design, y.view())
.map_err(EstimationError::InvalidInput)
}
pub fn dense_block_xtwx(
design: ArrayView2<'_, f64>,
fisher_blocks: ArrayView3<'_, f64>,
row_weights: Option<ArrayView1<'_, f64>>,
) -> Result<Array2<f64>, EstimationError> {
let n = design.nrows();
let k = design.ncols();
let shape = fisher_blocks.shape();
if shape.len() != 3 || shape[0] != n || shape[1] != shape[2] {
crate::bail_invalid_estim!(
"dense block Fisher shape mismatch: expected ({n}, p, p), got {shape:?}"
);
}
if let Some(w) = row_weights.as_ref() {
if w.len() != n {
crate::bail_invalid_estim!(
"dense block row weight length mismatch: expected {n}, got {}",
w.len()
);
}
if w.iter().any(|v| !v.is_finite() || *v < 0.0) {
crate::bail_invalid_estim!("dense block row weights must be finite and non-negative");
}
}
let p_out = shape[1];
let dim = k * p_out;
let mut out = Array2::<f64>::zeros((dim, dim));
for row in 0..n {
let rw = row_weights.as_ref().map(|w| w[row]).unwrap_or(1.0);
for a in 0..p_out {
for b in 0..p_out {
let wab = rw * fisher_blocks[[row, a, b]];
if !wab.is_finite() {
crate::bail_invalid_estim!(
"dense block Fisher entry ({row},{a},{b}) is not finite"
);
}
if wab == 0.0 {
continue;
}
let row_a = a * k;
let row_b = b * k;
for i in 0..k {
let xi = design[[row, i]];
if xi == 0.0 {
continue;
}
let scaled = wab * xi;
for j in 0..k {
out[[row_a + i, row_b + j]] += scaled * design[[row, j]];
}
}
}
}
}
for i in 0..dim {
for j in (i + 1)..dim {
let avg = 0.5 * (out[[i, j]] + out[[j, i]]);
out[[i, j]] = avg;
out[[j, i]] = avg;
}
}
Ok(out)
}
pub fn dense_block_xtwy(
design: ArrayView2<'_, f64>,
fisher_blocks: ArrayView3<'_, f64>,
response: ArrayView2<'_, f64>,
row_weights: Option<ArrayView1<'_, f64>>,
) -> Result<Array1<f64>, EstimationError> {
let n = design.nrows();
let k = design.ncols();
let shape = fisher_blocks.shape();
if shape.len() != 3 || shape[0] != n || shape[1] != shape[2] {
crate::bail_invalid_estim!(
"dense block Fisher shape mismatch: expected ({n}, p, p), got {shape:?}"
);
}
let p_out = shape[1];
if response.dim() != (n, p_out) {
crate::bail_invalid_estim!(
"dense block response shape mismatch: expected ({n}, {p_out}), got {}x{}",
response.nrows(),
response.ncols()
);
}
if let Some(w) = row_weights.as_ref()
&& w.len() != n
{
crate::bail_invalid_estim!(
"dense block row weight length mismatch: expected {n}, got {}",
w.len()
);
}
let mut out = Array1::<f64>::zeros(k * p_out);
for row in 0..n {
let rw = row_weights.as_ref().map(|w| w[row]).unwrap_or(1.0);
for a in 0..p_out {
let mut wy = 0.0_f64;
for b in 0..p_out {
let wab = rw * fisher_blocks[[row, a, b]];
if !wab.is_finite() {
crate::bail_invalid_estim!(
"dense block Fisher entry ({row},{a},{b}) is not finite"
);
}
wy += wab * response[[row, b]];
}
for i in 0..k {
out[a * k + i] += design[[row, i]] * wy;
}
}
}
Ok(out)
}
pub fn woodbury_gram_capacitance(
a_inv_uhat: &Array2<f64>,
vhat: &Array2<f64>,
) -> Result<Array2<f64>, EstimationError> {
LowRankWeight::gram_capacitance(a_inv_uhat, vhat).map_err(EstimationError::InvalidInput)
}
#[cfg(test)]
mod low_rank_weight_pirls_tests {
use super::{
DesignMatrix, LowRankWeight, PirlsWorkspace, compute_xtwx_low_rank, compute_xtwy_low_rank,
woodbury_gram_capacitance,
};
use crate::linalg::matrix::{LinearOperator, SignedWeightsView};
use ndarray::{Array2, array};
fn tiny_design() -> DesignMatrix {
let x = array![
[1.0, 0.5, -0.2],
[0.3, 1.2, 0.4],
[-0.1, 0.7, 1.0],
[0.6, -0.3, 0.8],
[0.2, 0.9, -0.5],
];
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(x))
}
#[test]
fn xtwx_low_rank_matches_diagonal_when_rank_zero() {
let design = tiny_design();
let d = array![1.0, 2.0, 0.5, 1.5, 0.8];
let u = Array2::<f64>::zeros((5, 0));
let v = Array2::<f64>::zeros((5, 0));
let weight = LowRankWeight::new(d.view(), u.view(), v.view()).unwrap();
let mut ws = PirlsWorkspace::new(5, 3, 0, 0);
let got = compute_xtwx_low_rank(&mut ws, &design, &weight).unwrap();
let want = design
.xt_diag_x_signed_op(SignedWeightsView::from_array(&d))
.unwrap();
let diff = (&got - &want).mapv(f64::abs).sum();
assert!(diff < 1e-12, "rank-0 path diverged from diagonal: {}", diff);
}
#[test]
fn xtwy_low_rank_matches_dense_reference() {
let design = tiny_design();
let d = array![1.0, 2.0, 0.5, 1.5, 0.8];
let u = array![
[0.1, -0.2],
[0.4, 0.3],
[-0.1, 0.5],
[0.2, 0.1],
[0.0, -0.3]
];
let v = array![[0.2, 0.1], [0.0, 0.4], [0.3, -0.2], [-0.1, 0.6], [0.5, 0.0]];
let weight = LowRankWeight::new(d.view(), u.view(), v.view()).unwrap();
let y = array![0.7, -1.2, 0.3, 0.9, -0.4];
let got = compute_xtwy_low_rank(&design, &weight, &y).unwrap();
let xdense = design.as_dense().unwrap().to_owned();
let mut w = Array2::<f64>::zeros((5, 5));
for i in 0..5 {
w[[i, i]] = d[i];
}
w += &u.dot(&v.t());
let want = xdense.t().dot(&w.dot(&y));
let diff: f64 = got
.iter()
.zip(want.iter())
.map(|(a, b)| (a - b).abs())
.sum();
assert!(diff < 1e-10, "xtwy_low_rank diverged: {}", diff);
}
#[test]
fn woodbury_capacitance_is_well_formed() {
let uhat = array![[0.5, 0.1], [-0.2, 0.7], [0.3, -0.4]];
let vhat = array![[0.1, 0.2], [0.6, -0.1], [-0.3, 0.4]];
let cap = woodbury_gram_capacitance(&uhat, &vhat).unwrap();
let want = {
let mut m = vhat.t().dot(&uhat);
for k in 0..2 {
m[[k, k]] += 1.0;
}
m
};
let diff: f64 = cap
.iter()
.zip(want.iter())
.map(|(a, b)| (a - b).abs())
.sum();
assert!(diff < 1e-12);
}
}
#[cfg(test)]
mod tests {
use super::loop_driver::default_beta_guess_external;
use super::reweight::madsen_lm_accept_factor;
use super::{
LinearInequalityConstraints, PenaltyConfig, PirlsConfig, PirlsLinearSolvePath,
PirlsProblem, PirlsWorkspace, bernoulli_geometry_from_jet, calculate_deviance,
compute_constraint_kkt_diagnostics, compute_observed_hessian_curvature_arrays,
fit_model_for_fixed_rho, select_active_set_release, should_log_pirls_decision_summary,
should_use_sparse_native_pirls, solve_newton_directionwith_linear_constraints,
solve_newton_directionwith_lower_bounds, update_glmvectors,
};
use crate::matrix::DesignMatrix;
use crate::mixture_link::InverseLinkJet as MixtureInverseLinkJet;
use crate::probability::standard_normal_quantile;
use crate::solver::active_set;
use crate::types::{
Coefficients, GlmLikelihoodSpec, InverseLink, LikelihoodSpec, LinkFunction,
LogSmoothingParamsView, ResponseFamily, StandardLink,
};
use approx::assert_relative_eq;
use faer::sparse::{SparseColMat, Triplet};
use ndarray::{Array1, Array2, ArrayView1, ArrayView2, ShapeBuilder, array};
#[test]
fn dense_workspace_xtwx_preserves_signed_observed_weights() {
let x = array![[1.0, 2.0], [3.0, -1.0], [-2.0, 4.0], [0.5, -3.0]];
let weights = array![2.0, -1.5, 0.25, -3.0];
let mut workspace = PirlsWorkspace::new(x.nrows(), x.ncols(), 0, 0);
let mut streamed = Array2::<f64>::zeros((x.ncols(), x.ncols()).f());
PirlsWorkspace::add_dense_xtwx_signed(
&weights,
&mut workspace.weighted_x_chunk,
&x,
&mut streamed,
);
let wx = Array2::from_shape_fn(x.raw_dim(), |(i, j)| weights[i] * x[[i, j]]);
let expected = x.t().dot(&wx);
for i in 0..x.ncols() {
for j in 0..x.ncols() {
assert_relative_eq!(streamed[[i, j]], expected[[i, j]], epsilon = 1e-12);
}
}
assert!(
streamed[[0, 0]] < 0.0,
"negative row weights must not be clipped through a sqrt(max(0,w)) Gram path"
);
}
fn calculate_scale(
beta: &Array1<f64>,
x: ArrayView2<f64>,
y: ArrayView1<f64>,
weights: ArrayView1<f64>,
offset: ArrayView1<f64>,
edf: f64,
link_function: LinkFunction,
) -> f64 {
match link_function {
LinkFunction::Logit
| LinkFunction::Probit
| LinkFunction::CLogLog
| LinkFunction::Sas
| LinkFunction::BetaLogistic
| LinkFunction::Log => 1.0,
LinkFunction::Identity => {
let mut fitted = x.dot(beta);
fitted += &offset;
let residuals = &y - &fitted;
let weighted_rss: f64 = weights
.iter()
.zip(residuals.iter())
.map(|(&w, &r)| w * r * r)
.sum();
let effective_n = y.len() as f64;
(weighted_rss / (effective_n - edf).max(1.0)).sqrt()
}
}
}
#[test]
fn madsen_lm_reject_trajectory_doubles_per_rejection() {
let mut loop_lambda = 1.0_f64;
let mut v = 2.0_f64;
let trajectory = (0..6)
.map(|_| {
loop_lambda *= v;
v *= 2.0;
loop_lambda
})
.collect::<Vec<_>>();
assert_eq!(
trajectory,
vec![2.0, 8.0, 64.0, 1024.0, 32_768.0, 2_097_152.0],
"Madsen rejection trajectory must double the multiplier each time"
);
}
#[test]
fn madsen_lm_accept_factor_matches_canonical_textbook_values() {
let cases: &[(f64, f64, &str)] = &[
(1.0, 1.0 / 3.0, "rho=1: floored at 1/3 (cube=1, 1-cube=0)"),
(0.75, 0.875, "rho=0.75: 1 - (0.5)^3 = 0.875 (slight shrink)"),
(0.5, 1.0, "rho=0.5: 1 - 0 = 1.0 (no change)"),
(
0.25,
1.125,
"rho=0.25: 1 - (-0.5)^3 = 1.125 (slight expand)",
),
];
for (rho, expected, why) in cases {
let got = madsen_lm_accept_factor(*rho);
assert!(
(got - expected).abs() < 1e-12,
"madsen_lm_accept_factor({rho}) = {got:.6}, expected {expected:.6} — {why}"
);
}
let small_positive = madsen_lm_accept_factor(1e-9);
assert!(
(small_positive - 2.0).abs() < 1e-6,
"rho ≈ 0⁺ must approach the 2.0 cap; got {small_positive:.6}"
);
assert_eq!(madsen_lm_accept_factor(-100.0), 2.0);
assert_eq!(madsen_lm_accept_factor(100.0), 1.0 / 3.0);
assert!(madsen_lm_accept_factor(0.99).is_finite());
assert!(madsen_lm_accept_factor(0.01) <= 2.0 + 1e-15);
assert!(madsen_lm_accept_factor(0.99) >= 1.0 / 3.0 - 1e-15);
}
#[test]
fn gaussian_scale_uses_offset_in_residuals() {
let x = array![[1.0], [2.0], [3.0]];
let beta = array![2.0];
let offset = array![10.0, 20.0, 30.0];
let y = array![12.0, 24.0, 36.0]; let w = Array1::ones(3);
let scale = calculate_scale(
&beta,
x.view(),
y.view(),
w.view(),
offset.view(),
0.0,
LinkFunction::Identity,
);
assert!(
scale.abs() < 1e-12,
"scale must be ~0 for exact fit with offset; got {}",
scale
);
}
#[test]
fn gaussian_scale_matchesweighted_sdwith_offset() {
let x = array![[1.0], [2.0], [4.0]];
let beta = array![1.5];
let offset = array![0.5, -1.0, 2.0];
let y = array![2.2, 2.0, 7.5];
let w = array![1.0, 2.0, 0.5];
let edf = 1.25;
let scale = calculate_scale(
&beta,
x.view(),
y.view(),
w.view(),
offset.view(),
edf,
LinkFunction::Identity,
);
let mut fitted = x.dot(&beta);
fitted += &offset;
let rss: f64 = w
.iter()
.zip(y.iter().zip(fitted.iter()))
.map(|(&wi, (&yi, &fi))| wi * (yi - fi).powi(2))
.sum();
let expected = (rss / ((y.len() as f64 - edf).max(1.0))).sqrt();
assert!(
(scale - expected).abs() < 1e-12,
"scale mismatch: got {}, expected {}",
scale,
expected
);
}
#[test]
fn kkt_diagnosticszero_for_strictly_feasible_stationary_point() {
let constraints = LinearInequalityConstraints {
a: array![[1.0, 0.0], [0.0, 1.0]],
b: array![0.0, 0.0],
};
let beta = array![1.0, 2.0];
let grad = array![0.0, 0.0];
let diag = compute_constraint_kkt_diagnostics(&beta, &grad, &constraints);
assert!(diag.primal_feasibility <= 1e-12);
assert!(diag.dual_feasibility <= 1e-12);
assert!(diag.complementarity <= 1e-12);
assert!(diag.stationarity <= 1e-12);
}
#[test]
fn kkt_diagnostics_capture_active_lower_bound_solution() {
let constraints = LinearInequalityConstraints {
a: array![[1.0, 0.0], [0.0, 1.0]],
b: array![0.0, 0.0],
};
let beta = array![0.0, 1.5];
let grad = array![2.0, 0.0];
let diag = compute_constraint_kkt_diagnostics(&beta, &grad, &constraints);
assert_eq!(diag.n_constraints, 2);
assert_eq!(diag.n_active, 1);
assert!(diag.primal_feasibility <= 1e-12);
assert!(diag.dual_feasibility <= 1e-12);
assert!(diag.complementarity <= 1e-12);
assert!(diag.stationarity <= 1e-10);
}
#[test]
fn linear_constraint_active_set_releases_positive_kkt_systemmultiplier() {
let hessian = array![[1.0]];
let gradient = array![-1.0];
let beta = array![0.0];
let constraints = LinearInequalityConstraints {
a: array![[1.0], [-1.0]],
b: array![0.0, -0.1],
};
let mut direction = Array1::zeros(1);
solve_newton_directionwith_linear_constraints(
&hessian,
&gradient,
&beta,
&constraints,
&mut direction,
None,
)
.expect("constrained Newton direction should solve");
assert!(
(direction[0] - 0.1).abs() <= 1e-10,
"expected step to upper bound (0.1), got {}",
direction[0]
);
}
#[test]
fn linear_constraint_active_set_ignores_near_tangential_inactiverows() {
let hessian = array![[1.0, 0.0], [0.0, 1.0]];
let gradient = array![-1.0, 0.0];
let beta = array![0.0, 0.0];
let constraints = LinearInequalityConstraints {
a: array![[-1e-16, 1.0]],
b: array![-1.0],
};
let mut direction = Array1::zeros(2);
solve_newton_directionwith_linear_constraints(
&hessian,
&gradient,
&beta,
&constraints,
&mut direction,
None,
)
.expect("near-tangential inactive row should not block the Newton step");
assert!(
(direction[0] - 1.0).abs() <= 1e-12,
"expected unconstrained x-step of 1.0, got {}",
direction[0]
);
assert!(
direction[1].abs() <= 1e-12,
"expected zero y-step, got {}",
direction[1]
);
}
#[test]
fn default_beta_guess_logit_uses_log_odds_prevalence() {
let y = array![0.0, 1.0, 1.0, 1.0];
let w = Array1::ones(4);
let beta =
default_beta_guess_external(3, LinkFunction::Logit, y.view(), w.view(), None, None);
let prevalence: f64 = (3.0 + 0.5) / (4.0 + 1.0);
let prevalence = prevalence.max(1e-6_f64).min(1.0_f64 - 1e-6_f64);
let expected = (prevalence / (1.0 - prevalence)).ln();
assert!((beta[0] - expected).abs() < 1e-12);
assert_eq!(beta[1], 0.0);
assert_eq!(beta[2], 0.0);
}
#[test]
fn default_beta_guess_probit_uses_standard_normal_quantile() {
let y = array![0.0, 1.0, 1.0, 1.0];
let w = Array1::ones(4);
let beta =
default_beta_guess_external(3, LinkFunction::Probit, y.view(), w.view(), None, None);
let prevalence: f64 = (3.0 + 0.5) / (4.0 + 1.0);
let prevalence = prevalence.max(1e-6_f64).min(1.0_f64 - 1e-6_f64);
let log_odds = (prevalence / (1.0 - prevalence)).ln();
let expected =
standard_normal_quantile(prevalence).expect("clamped prevalence must be valid");
assert!((expected - log_odds).abs() > 1e-3);
assert!((beta[0] - expected).abs() < 1e-12);
assert_eq!(beta[1], 0.0);
assert_eq!(beta[2], 0.0);
}
#[test]
fn sparse_native_decision_rejects_dense_design() {
let x = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(array![
[1.0, 0.0],
[0.0, 1.0]
]));
let s = array![[1.0, 0.0], [0.0, 1.0]];
let mut workspace = PirlsWorkspace::new(2, 2, 0, 0);
let decision = should_use_sparse_native_pirls(&mut workspace, &x, &s, None, None);
assert_eq!(decision.path, PirlsLinearSolvePath::DenseTransformed);
assert_eq!(decision.reason, "design_not_sparse");
}
fn fixed_gaussian_beta(
x: Array2<f64>,
y: Array1<f64>,
penalties: Vec<crate::smooth::BlockwisePenalty>,
rho: Array1<f64>,
) -> Array1<f64> {
let p = x.ncols();
let weights = Array1::<f64>::ones(y.len());
let offset = Array1::<f64>::zeros(y.len());
let specs: Vec<crate::estimate::PenaltySpec> = penalties
.iter()
.map(crate::estimate::PenaltySpec::from_blockwise_ref)
.collect();
let nulls = vec![0; specs.len()];
let (canonical, _) =
crate::construction::canonicalize_penalty_specs(&specs, &nulls, p, "prior mean test")
.expect("canonical penalties");
let config = PirlsConfig {
likelihood: GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Gaussian,
InverseLink::Standard(StandardLink::Identity),
)),
link_kind: InverseLink::Standard(StandardLink::Identity),
max_iterations: 20,
convergence_tolerance: 1e-12,
firth_bias_reduction: false,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let problem = PirlsProblem {
x,
offset: offset.view(),
y: y.view(),
priorweights: weights.view(),
covariate_se: None,
gaussian_fixed_cache: None,
};
let penalty = PenaltyConfig {
canonical_penalties: &canonical,
balanced_penalty_root: None,
reparam_invariant: None,
p,
coefficient_lower_bounds: None,
linear_constraints_original: None,
penalty_shrinkage_floor: None,
kronecker_factored: None,
};
let (fit, _) = fit_model_for_fixed_rho(
LogSmoothingParamsView::new(rho.view()),
problem,
penalty,
&config,
None,
)
.expect("fixed rho fit");
fit.beta_transformed.as_ref().clone()
}
#[test]
fn constant_prior_mean_centers_penalty() {
let x = Array2::<f64>::zeros((4, 1));
let y = Array1::<f64>::zeros(4);
let penalty = crate::smooth::BlockwisePenalty::ridge(0..1, 1.0)
.with_prior_mean(crate::estimate::CoefficientPriorMean::scalar(2.5));
let beta = fixed_gaussian_beta(x, y, vec![penalty], array![0.0]);
assert!((beta[0] - 2.5).abs() < 1e-10, "beta={beta:?}");
}
#[test]
fn functional_prior_mean_recovers_kernel_amplitude() {
let x = Array2::<f64>::zeros((5, 3));
let y = Array1::<f64>::zeros(5);
let metadata = array![2.0];
let alpha = 1.75;
let penalty = crate::smooth::BlockwisePenalty::ridge(0..3, 1.0).with_prior_mean(
crate::estimate::CoefficientPriorMean::functional(
metadata,
std::sync::Arc::new(move |a: &Array1<f64>| {
let t = a[0];
array![alpha, alpha * t, alpha * t * t]
}),
),
);
let beta = fixed_gaussian_beta(x, y, vec![penalty], array![0.0]);
let recovered_alpha = beta[0];
assert!((recovered_alpha - alpha).abs() < 1e-10, "beta={beta:?}");
assert!((beta[1] / 2.0 - alpha).abs() < 1e-10, "beta={beta:?}");
assert!((beta[2] / 4.0 - alpha).abs() < 1e-10, "beta={beta:?}");
}
#[test]
fn zero_prior_mean_matches_default_fixed_fit_bitwise() {
let x = array![[1.0, 0.0], [1.0, 1.0], [1.0, 2.0], [1.0, 3.0], [1.0, 4.0],];
let y = array![0.5, 1.0, 1.5, 2.0, 2.5];
let base_penalty = crate::smooth::BlockwisePenalty::ridge(0..2, 1.0);
let zero_penalty = crate::smooth::BlockwisePenalty::ridge(0..2, 1.0).with_prior_mean(
crate::estimate::CoefficientPriorMean::constant(Array1::zeros(2)),
);
let rho = array![0.25];
let beta_default =
fixed_gaussian_beta(x.clone(), y.clone(), vec![base_penalty], rho.clone());
let beta_zero = fixed_gaussian_beta(x, y, vec![zero_penalty], rho);
assert_eq!(beta_default.to_vec(), beta_zero.to_vec());
}
#[test]
fn pirls_decision_summary_logs_on_power_of_two_repetitions() {
assert!(!should_log_pirls_decision_summary(1));
assert!(should_log_pirls_decision_summary(2));
assert!(!should_log_pirls_decision_summary(3));
assert!(should_log_pirls_decision_summary(4));
assert!(!should_log_pirls_decision_summary(6));
assert!(should_log_pirls_decision_summary(8));
}
#[test]
fn sparse_native_decision_collects_sparse_stats_for_large_sparse_design() {
let triplets: Vec<_> = (0..300).map(|i| Triplet::new(i, i, 1.0)).collect();
let x = SparseColMat::try_new_from_triplets(300, 300, &triplets)
.expect("sparse identity should build");
let x = DesignMatrix::from(x);
let s = Array2::from_diag(&Array1::ones(300));
let mut workspace = PirlsWorkspace::new(300, 300, 0, 0);
let decision = should_use_sparse_native_pirls(&mut workspace, &x, &s, None, None);
assert_eq!(decision.path, PirlsLinearSolvePath::SparseNative);
assert_eq!(decision.reason, "sparse_native_eligible");
assert_eq!(decision.nnz_x, 300);
assert_eq!(decision.nnz_xtwx_symbolic, Some(300));
assert_eq!(decision.nnz_h_est, Some(300));
assert!(decision.density_h_est.expect("density") < 0.01);
}
#[test]
fn sparse_native_decision_allows_moderate_sparse_designs_below_old_width_gate() {
let triplets: Vec<_> = (0..64).map(|i| Triplet::new(i, i, 1.0)).collect();
let x = SparseColMat::try_new_from_triplets(64, 64, &triplets)
.expect("sparse identity should build");
let x = DesignMatrix::from(x);
let s = Array2::from_diag(&Array1::ones(64));
let mut workspace = PirlsWorkspace::new(64, 64, 0, 0);
let decision = should_use_sparse_native_pirls(&mut workspace, &x, &s, None, None);
assert_eq!(decision.path, PirlsLinearSolvePath::SparseNative);
assert_eq!(decision.reason, "sparse_native_eligible");
assert_eq!(decision.nnz_x, 64);
assert_eq!(decision.nnz_xtwx_symbolic, Some(64));
assert_eq!(decision.nnz_h_est, Some(64));
assert!(decision.density_h_est.expect("density") < 0.05);
}
#[test]
fn sparse_native_decision_rejects_finite_lower_bounds() {
let triplets: Vec<_> = (0..64).map(|i| Triplet::new(i, i, 1.0)).collect();
let x = SparseColMat::try_new_from_triplets(64, 64, &triplets)
.expect("sparse identity should build");
let x = DesignMatrix::from(x);
let s = Array2::from_diag(&Array1::ones(64));
let mut lower_bounds = Array1::from_elem(64, f64::NEG_INFINITY);
lower_bounds[0] = 0.0;
let mut workspace = PirlsWorkspace::new(64, 64, 0, 0);
let decision =
should_use_sparse_native_pirls(&mut workspace, &x, &s, Some(&lower_bounds), None);
assert_eq!(decision.path, PirlsLinearSolvePath::DenseTransformed);
assert_eq!(decision.reason, "constraints_present");
}
#[test]
fn sparse_penalized_assembly_matches_dense_diagonal_case() {
let triplets = vec![
Triplet::new(0, 0, 1.0),
Triplet::new(1, 1, 2.0),
Triplet::new(2, 2, 3.0),
];
let x = SparseColMat::try_new_from_triplets(3, 3, &triplets)
.expect("diagonal sparse matrix should build");
let weights = array![2.0, 3.0, 5.0];
let s_lambda = array![[4.0, 0.0, 0.0], [0.0, 6.0, 0.0], [0.0, 0.0, 8.0]];
let ridge = 1e-8;
let mut workspace = PirlsWorkspace::new(3, 3, 0, 0);
let assembled = super::sparse_reml_penalized_hessian(
&mut workspace,
&x,
&weights,
&s_lambda,
ridge,
None,
)
.expect("sparse penalized assembly should succeed");
let dense = DesignMatrix::from(x.clone()).to_dense();
let mut expected = dense.t().dot(&Array2::from_diag(&weights)).dot(&dense);
expected += &s_lambda;
for i in 0..3 {
expected[[i, i]] += ridge;
}
let actual = DesignMatrix::from(assembled).to_dense();
for i in 0..3 {
for j in 0..3 {
let target = if i <= j { expected[[i, j]] } else { 0.0 };
assert!(
(actual[[i, j]] - target).abs() < 1e-10,
"mismatch at ({}, {}): {} vs {}",
i,
j,
actual[[i, j]],
target
);
}
}
}
#[test]
fn pirls_result_stores_integrated_logit_derivative_jet() {
assert!(file!().ends_with(".rs"));
let x = array![[1.0], [1.0], [1.0], [1.0], [1.0]];
let y = array![0.0, 1.0, 0.0, 1.0, 1.0];
let w = Array1::ones(5);
let offset = Array1::zeros(5);
let rho = Array1::<f64>::zeros(1);
let covariate_se = array![0.9, 0.7, 0.8, 0.6, 0.75];
let rs = [array![[1.0]]];
let canonical: Vec<crate::construction::CanonicalPenalty> = rs
.iter()
.map(|r| {
let local = r.t().dot(r);
crate::construction::CanonicalPenalty {
root: r.clone(),
col_range: 0..r.ncols(),
total_dim: r.ncols(),
nullity: 0,
local,
prior_mean: Array1::zeros(r.ncols()),
positive_eigenvalues: Vec::new(),
op: None,
}
})
.collect();
let config = PirlsConfig {
likelihood: GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Binomial,
InverseLink::Standard(StandardLink::Logit),
)),
link_kind: InverseLink::Standard(StandardLink::Logit),
max_iterations: 100,
convergence_tolerance: 1e-8,
firth_bias_reduction: false,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let (fit, _) = fit_model_for_fixed_rho(
LogSmoothingParamsView::new(rho.view()),
PirlsProblem {
x: x.view(),
offset: offset.view(),
y: y.view(),
priorweights: w.view(),
covariate_se: Some(covariate_se.view()),
gaussian_fixed_cache: None,
},
PenaltyConfig {
canonical_penalties: &canonical,
balanced_penalty_root: None,
reparam_invariant: None,
p: 1,
coefficient_lower_bounds: None,
linear_constraints_original: None,
penalty_shrinkage_floor: None,
kronecker_factored: None,
},
&config,
Some(&Coefficients::new(array![0.0])),
)
.expect("integrated logit PIRLS fit");
let ctx = crate::quadrature::QuadratureContext::new();
for i in 0..y.len() {
let jet = crate::quadrature::integrated_inverse_link_jet(
&ctx,
LinkFunction::Logit,
fit.final_eta[i].clamp(-700.0, 700.0),
covariate_se[i],
)
.expect("logit integrated inverse-link jet should evaluate");
let expected = bernoulli_geometry_from_jet(
fit.final_eta[i],
fit.final_eta[i].clamp(-700.0, 700.0),
y[i],
w[i],
MixtureInverseLinkJet {
mu: jet.mean,
d1: jet.d1,
d2: jet.d2,
d3: jet.d3,
},
);
assert_relative_eq!(
fit.solve_dmu_deta[i],
jet.d1,
epsilon = 1e-9,
max_relative = 1e-9
);
assert_relative_eq!(
fit.solve_d2mu_deta2[i],
jet.d2,
epsilon = 1e-9,
max_relative = 1e-8
);
assert_relative_eq!(
fit.solve_d3mu_deta3[i],
jet.d3,
epsilon = 1e-8,
max_relative = 1e-7
);
assert_relative_eq!(
fit.solve_c_array[i],
expected.c,
epsilon = 1e-9,
max_relative = 1e-8
);
assert_relative_eq!(
fit.solve_d_array[i],
expected.d,
epsilon = 1e-8,
max_relative = 1e-7
);
}
}
#[test]
fn pure_logit_working_state_preserves_tail_fisher_mass() {
let y = array![1.0];
let eta = array![50.0];
let priorweights = array![1.0];
let inverse_link = InverseLink::Standard(StandardLink::Logit);
let mut mu = Array1::zeros(1);
let mut weights = Array1::zeros(1);
let mut z = Array1::zeros(1);
update_glmvectors(
y.view(),
&eta,
&inverse_link,
priorweights.view(),
&mut mu,
&mut weights,
&mut z,
None,
)
.expect("pure logit working state");
let jet = crate::mixture_link::logit_inverse_link_jet5(eta[0]);
assert!(jet.d1 > 0.0);
assert!(
(weights[0] - jet.d1).abs() < 1e-30,
"pure logit PIRLS weight should equal the stable tail formula at eta={}; got {} vs {}",
eta[0],
weights[0],
jet.d1
);
assert!(
(mu[0] - jet.mu).abs() < 1e-30,
"pure logit PIRLS mu mismatch at eta={}; got {} vs {}",
eta[0],
mu[0],
jet.mu
);
let expected_z = eta[0] + (y[0] - jet.mu) / jet.d1;
assert!(
(z[0] - expected_z).abs() < 1e-12,
"pure logit PIRLS z should preserve the exact working response at eta={}; got {} vs {}",
eta[0],
z[0],
expected_z
);
assert!(
(weights[0] * (z[0] - eta[0]) - (y[0] - jet.mu)).abs() < 1e-30,
"pure logit PIRLS score carrier should preserve y-mu at eta={}; got {} vs {}",
eta[0],
weights[0] * (z[0] - eta[0]),
y[0] - jet.mu
);
}
#[test]
fn noncanonical_binomial_working_state_clamps_saturating_standard_links() {
for link in [StandardLink::Probit, StandardLink::CLogLog] {
let y = array![1.0];
let eta = array![30.0];
let priorweights = array![1.0];
let inverse_link = InverseLink::Standard(link);
let mut mu = Array1::zeros(1);
let mut weights = Array1::zeros(1);
let mut z = Array1::zeros(1);
update_glmvectors(
y.view(),
&eta,
&inverse_link,
priorweights.view(),
&mut mu,
&mut weights,
&mut z,
None,
)
.expect("noncanonical binomial working state");
assert!(
mu[0] > 0.0 && mu[0] < 1.0,
"{link:?} working mu must stay inside (0,1) before variance evaluation; got {}",
mu[0]
);
assert!(
weights[0].is_finite() && weights[0] > 0.0,
"{link:?} working weight must remain positive finite at saturated eta; got {}",
weights[0]
);
assert!(
z[0].is_finite(),
"{link:?} working response must remain finite at saturated eta; got {}",
z[0]
);
}
}
#[test]
fn gamma_log_deviance_uses_gamma_formula() {
assert!(file!().ends_with(".rs"));
let y = array![2.0, 5.0];
let mu = array![1.0, 4.0];
let w = array![1.5, 0.75];
let dev = calculate_deviance(
y.view(),
&mu,
&GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Gamma,
InverseLink::Standard(StandardLink::Log),
)),
w.view(),
);
let expected = 2.0
* (1.5 * (2.0_f64 / 1.0 - 1.0 - (2.0_f64 / 1.0).ln())
+ 0.75 * (5.0_f64 / 4.0 - 1.0 - (5.0_f64 / 4.0).ln()));
assert_relative_eq!(dev, expected, epsilon = 1e-12, max_relative = 1e-12);
}
#[test]
fn gamma_log_observed_curvature_matches_shape_one_closed_form() {
assert!(file!().ends_with(".rs"));
let eta = array![0.2, -0.4];
let mu = eta.mapv(f64::exp);
let y = array![1.8, 0.7];
let w = array![2.0, 0.5];
let fisher = w.clone();
let (w_obs, c_obs, d_obs) = compute_observed_hessian_curvature_arrays(
&GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Gamma,
InverseLink::Standard(StandardLink::Log),
)),
&InverseLink::Standard(StandardLink::Log),
&eta,
y.view(),
&fisher,
w.view(),
)
.expect("gamma-log observed curvature should evaluate");
for i in 0..eta.len() {
let expected_w = w[i] * y[i] / mu[i];
assert_relative_eq!(w_obs[i], expected_w, epsilon = 1e-12, max_relative = 1e-12);
assert_relative_eq!(c_obs[i], -expected_w, epsilon = 1e-12, max_relative = 1e-12);
assert_relative_eq!(d_obs[i], expected_w, epsilon = 1e-12, max_relative = 1e-12);
}
}
#[test]
fn negative_binomial_log_observed_curvature_matches_size_theta_closed_form() {
assert!(file!().ends_with(".rs"));
let theta = 2.5;
let eta = array![0.2, -0.4, 1.1];
let mu = eta.mapv(f64::exp);
let y = array![0.0, 3.0, 8.0];
let w = array![2.0, 0.5, 1.25];
let fisher = Array1::from_iter(
mu.iter()
.zip(w.iter())
.map(|(&mu_i, &w_i)| w_i * theta * mu_i / (theta + mu_i)),
);
let (w_obs, c_obs, d_obs) = compute_observed_hessian_curvature_arrays(
&GlmLikelihoodSpec::canonical(LikelihoodSpec::negative_binomial_log(theta)),
&InverseLink::Standard(StandardLink::Log),
&eta,
y.view(),
&fisher,
w.view(),
)
.expect("negative-binomial-log observed curvature should evaluate");
for i in 0..eta.len() {
let denom = theta + mu[i];
let scale = w[i] * theta * (theta + y[i]);
let expected_w = scale * mu[i] / (denom * denom);
let expected_c = scale * mu[i] * (theta - mu[i]) / (denom * denom * denom);
let expected_d = scale * mu[i] * (theta * theta - 4.0 * theta * mu[i] + mu[i] * mu[i])
/ (denom * denom * denom * denom);
assert_relative_eq!(w_obs[i], expected_w, epsilon = 1e-12, max_relative = 1e-12);
assert_relative_eq!(c_obs[i], expected_c, epsilon = 1e-12, max_relative = 1e-12);
assert_relative_eq!(d_obs[i], expected_d, epsilon = 1e-12, max_relative = 1e-12);
}
}
#[test]
fn gamma_log_fit_profiles_shape_instead_of_fixing_one() {
let x = array![[1.0], [1.0], [1.0], [1.0], [1.0], [1.0]];
let y = array![0.8, 1.1, 1.7, 2.0, 2.6, 3.1];
let w = Array1::ones(y.len());
let offset = Array1::zeros(y.len());
let rho = array![0.0];
let rs = [array![[0.0]]];
let canonical: Vec<crate::construction::CanonicalPenalty> = rs
.iter()
.map(|r| {
let local = r.t().dot(r);
crate::construction::CanonicalPenalty {
root: r.clone(),
col_range: 0..r.ncols(),
total_dim: r.ncols(),
nullity: 0,
local,
prior_mean: Array1::zeros(r.ncols()),
positive_eigenvalues: Vec::new(),
op: None,
}
})
.collect();
let config = PirlsConfig {
likelihood: GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Gamma,
InverseLink::Standard(StandardLink::Log),
)),
link_kind: InverseLink::Standard(StandardLink::Log),
max_iterations: 100,
convergence_tolerance: 1e-8,
firth_bias_reduction: false,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let (result, _) = fit_model_for_fixed_rho(
LogSmoothingParamsView::new(rho.view()),
PirlsProblem {
x: x.view(),
offset: offset.view(),
y: y.view(),
priorweights: w.view(),
covariate_se: None,
gaussian_fixed_cache: None,
},
PenaltyConfig {
canonical_penalties: &canonical,
balanced_penalty_root: None,
reparam_invariant: None,
p: 1,
coefficient_lower_bounds: None,
linear_constraints_original: None,
penalty_shrinkage_floor: None,
kronecker_factored: None,
},
&config,
None,
)
.expect("gamma PIRLS fit");
let fitted_shape = result
.likelihood
.gamma_shape()
.expect("gamma fit should expose fitted shape");
let profiled_shape =
super::estimate_gamma_shape_from_eta(y.view(), &result.final_eta, w.view());
assert!(fitted_shape > 1.0, "shape should not stay fixed at one");
assert_relative_eq!(
fitted_shape,
profiled_shape,
epsilon = 1e-10,
max_relative = 1e-10
);
}
#[test]
fn poisson_cache_rehydration_preserves_log_derivatives() {
let x = array![[1.0], [1.0], [1.0], [1.0]];
let y = array![1.0, 2.0, 4.0, 8.0];
let w = Array1::ones(4);
let offset = Array1::zeros(4);
let rho = array![0.0];
let rs = [array![[1.0]]];
let canonical: Vec<crate::construction::CanonicalPenalty> = rs
.iter()
.map(|r| {
let local = r.t().dot(r);
crate::construction::CanonicalPenalty {
root: r.clone(),
col_range: 0..r.ncols(),
total_dim: r.ncols(),
nullity: 0,
local,
prior_mean: Array1::zeros(r.ncols()),
positive_eigenvalues: Vec::new(),
op: None,
}
})
.collect();
let config = PirlsConfig {
likelihood: GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Poisson,
InverseLink::Standard(StandardLink::Log),
)),
link_kind: InverseLink::Standard(StandardLink::Log),
max_iterations: 100,
convergence_tolerance: 1e-8,
firth_bias_reduction: false,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let (fit, _) = fit_model_for_fixed_rho(
LogSmoothingParamsView::new(rho.view()),
PirlsProblem {
x: x.view(),
offset: offset.view(),
y: y.view(),
priorweights: w.view(),
covariate_se: None,
gaussian_fixed_cache: None,
},
PenaltyConfig {
canonical_penalties: &canonical,
balanced_penalty_root: None,
reparam_invariant: None,
p: 1,
coefficient_lower_bounds: None,
linear_constraints_original: None,
penalty_shrinkage_floor: None,
kronecker_factored: None,
},
&config,
None,
)
.expect("poisson PIRLS fit");
let compacted = fit.compact_for_reml_cache();
let rehydrated = compacted
.rehydrate_after_reml_cache(
&DesignMatrix::from(x.clone()),
y.view(),
w.view(),
offset.view(),
&InverseLink::Standard(StandardLink::Log),
)
.expect("rehydration should succeed");
assert_eq!(fit.solve_c_array.len(), rehydrated.solve_c_array.len());
for i in 0..fit.solve_c_array.len() {
assert_relative_eq!(
fit.solve_c_array[i],
rehydrated.solve_c_array[i],
epsilon = 1e-12,
max_relative = 1e-12
);
assert_relative_eq!(
fit.solve_d_array[i],
rehydrated.solve_d_array[i],
epsilon = 1e-12,
max_relative = 1e-12
);
}
}
#[test]
fn linear_constraint_active_set_releases_stalewarm_boundary_hint() {
let hessian = array![[2.0]];
let gradient = array![0.0];
let beta = array![1e-9];
let constraints = LinearInequalityConstraints {
a: array![[1.0]],
b: array![0.0],
};
let mut direction = Array1::zeros(1);
let mut active_hint = vec![0];
solve_newton_directionwith_linear_constraints(
&hessian,
&gradient,
&beta,
&constraints,
&mut direction,
Some(&mut active_hint),
)
.expect("active-set solve should succeed");
assert_relative_eq!(direction[0], 0.0, epsilon = 1e-14);
let projected = &beta + &direction;
assert_relative_eq!(projected[0], beta[0], epsilon = 1e-14);
assert!(active_hint.is_empty());
}
#[test]
fn linear_constraint_active_set_releases_stalewarm_hint() {
let hessian = array![[1.0]];
let gradient = array![-1.0];
let beta = array![0.0];
let constraints = LinearInequalityConstraints {
a: array![[1.0], [-1.0]],
b: array![0.0, -0.1],
};
let mut direction = Array1::zeros(1);
let mut active_hint = vec![0];
solve_newton_directionwith_linear_constraints(
&hessian,
&gradient,
&beta,
&constraints,
&mut direction,
Some(&mut active_hint),
)
.expect("stale warm active-set hint should be releasable");
assert!(
(direction[0] - 0.1).abs() <= 1e-10,
"expected step to upper bound (0.1), got {}",
direction[0]
);
assert_eq!(active_hint, vec![1]);
}
#[test]
fn working_set_kkt_diagnostics_use_active_setmultipliers() {
let working_constraints = LinearInequalityConstraints {
a: array![[1.0, 0.0], [2.0, 0.0], [0.0, 1.0]],
b: array![0.0, 0.0, 0.0],
};
let x = array![0.0, 0.0];
let lambda_true = array![1.0, 0.5, 2.0];
let gradient = working_constraints.a.t().dot(&lambda_true);
let kkt = active_set::working_set_kkt_diagnostics_from_multipliers(
&x,
&gradient,
&working_constraints,
&lambda_true,
3,
)
.expect("working-set KKT diagnostics");
assert!(kkt.primal_feasibility <= 1e-12);
assert!(kkt.dual_feasibility <= 1e-12);
assert!(kkt.complementarity <= 1e-12);
assert!(kkt.stationarity <= 1e-12);
assert_eq!(kkt.n_active, 3);
}
#[test]
fn compress_activeworking_set_groups_near_collinearrows() {
let constraints = LinearInequalityConstraints {
a: array![
[0.0, 0.5, 0.0],
[0.0, 0.50000000000003, 0.0],
[1.0, 0.0, 0.0]
],
b: array![1e-8, 1.00000000000005e-8, 0.2],
};
let x = array![0.0, 0.0, 0.0];
let active = vec![0, 1, 2];
let compressed = active_set::compress_active_working_set(&x, &constraints, &active)
.expect("compress working set");
assert_eq!(compressed.constraints.a.nrows(), 2);
assert_eq!(compressed.groups.len(), 2);
assert!(
compressed.groups.iter().any(|g| g == &vec![0, 1]),
"near-collinear rows should be grouped together: {:?}",
compressed.groups
);
}
#[test]
fn lower_bound_active_set_releases_stalewarm_boundary_hint() {
let hessian = array![[2.0]];
let gradient = array![0.0];
let beta = array![1e-9];
let lower_bounds = array![0.0];
let mut direction = Array1::zeros(1);
let mut active_hint = vec![0];
solve_newton_directionwith_lower_bounds(
&hessian,
&gradient,
&beta,
&lower_bounds,
&mut direction,
Some(&mut active_hint),
)
.expect("lower-bound active-set solve should succeed");
assert_relative_eq!(direction[0], 0.0, epsilon = 1e-14);
let projected = &beta + &direction;
assert_relative_eq!(projected[0], beta[0], epsilon = 1e-14);
assert!(active_hint.is_empty());
}
#[test]
fn select_active_set_release_worst_violation_picks_most_negative() {
let gradient = array![-0.1, -0.5, -0.2];
let hd = array![0.0, 0.0, 0.0];
let active_idx = vec![0, 1, 2];
assert_eq!(
select_active_set_release(&gradient, &hd, &active_idx, false),
Some(1)
);
}
#[test]
fn select_active_set_release_blands_picks_lowest_index_with_negative_multiplier() {
let gradient = array![-0.1, -0.5, -0.2];
let hd = array![0.0, 0.0, 0.0];
let active_idx = vec![0, 1, 2];
assert_eq!(
select_active_set_release(&gradient, &hd, &active_idx, true),
Some(0)
);
}
#[test]
fn select_active_set_release_blands_deadband_ignores_round_off() {
let g = 1.0_f64;
let lambda_noise = -32.0 * f64::EPSILON * g; let gradient = array![g];
let hd = array![lambda_noise - g]; let active_idx = vec![0];
assert_eq!(
select_active_set_release(&gradient, &hd, &active_idx, true),
None,
"round-off-level multiplier must not trigger Bland's release"
);
let lambda_real = -128.0 * f64::EPSILON * g;
let hd = array![lambda_real - g];
assert_eq!(
select_active_set_release(&gradient, &hd, &active_idx, true),
Some(0)
);
}
#[test]
fn select_active_set_release_returns_none_when_kkt_satisfied() {
let gradient = array![0.5, 1.0, 0.0];
let hd = array![0.0, 0.0, 0.0];
let active_idx = vec![0, 1, 2];
assert_eq!(
select_active_set_release(&gradient, &hd, &active_idx, false),
None
);
assert_eq!(
select_active_set_release(&gradient, &hd, &active_idx, true),
None
);
}
#[test]
fn lower_bound_active_set_releases_stalewarm_hint() {
let hessian = array![[1.0]];
let gradient = array![-1.0];
let beta = array![0.0];
let lower_bounds = array![0.0];
let mut direction = Array1::zeros(1);
let mut active_hint = vec![0];
solve_newton_directionwith_lower_bounds(
&hessian,
&gradient,
&beta,
&lower_bounds,
&mut direction,
Some(&mut active_hint),
)
.expect("stale warm lower-bound hint should be releasable");
assert!(
(direction[0] - 1.0).abs() <= 1e-12,
"expected unconstrained step of 1.0 after releasing stale bound, got {}",
direction[0]
);
assert!(active_hint.is_empty());
}
}
#[cfg(test)]
mod root_cause_tests {
use super::*;
use crate::types::LogSmoothingParamsView;
use approx::assert_relative_eq;
use ndarray::{Array1, Array2, array};
fn capture_pirls_penalized_deviance<F, R>(run: F) -> (R, Vec<f64>)
where
F: FnOnce() -> R,
{
super::reweight::test_support::PIRLS_PENALIZED_DEVIANCE_TRACE.with(|trace| {
*trace.borrow_mut() = Some(Vec::new());
});
let result = run();
let captured = super::reweight::test_support::PIRLS_PENALIZED_DEVIANCE_TRACE
.with(|trace| trace.borrow_mut().take().unwrap());
(result, captured)
}
fn scalar_working_state(
beta: &Coefficients,
curvature: HessianCurvatureKind,
gradient: f64,
deviance: f64,
) -> WorkingState {
WorkingState {
eta: LinearPredictor::new(array![beta.as_ref()[0]]),
gradient: array![gradient],
hessian: crate::linalg::matrix::SymmetricMatrix::Dense(array![[1.0]]),
log_likelihood: 0.0,
deviance,
penalty_term: 0.0,
firth: FirthDiagnostics::Inactive,
ridge_used: 0.0,
hessian_curvature: curvature,
gradient_natural_scale: 0.0,
}
}
fn test_working_state(beta: &Coefficients, curvature: HessianCurvatureKind) -> WorkingState {
scalar_working_state(beta, curvature, 1.0, 1.0)
}
#[derive(Default)]
struct CandidateEvalFailureModel {
observed_updates: usize,
fisher_updates: usize,
observed_candidate_calls: usize,
fisher_candidate_calls: usize,
}
impl CandidateEvalFailureModel {
fn state(beta: &Coefficients, curvature: HessianCurvatureKind) -> WorkingState {
test_working_state(beta, curvature)
}
}
impl WorkingModel for CandidateEvalFailureModel {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
match curvature {
HessianCurvatureKind::Observed => self.observed_updates += 1,
HessianCurvatureKind::Fisher => self.fisher_updates += 1,
}
Ok(Self::state(beta, curvature))
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
match curvature {
HessianCurvatureKind::Observed => self.observed_candidate_calls += 1,
HessianCurvatureKind::Fisher => self.fisher_candidate_calls += 1,
}
Err(EstimationError::InvalidInput(format!(
"non-finite candidate evaluation under {curvature:?} curvature at beta={:.3e}",
beta.as_ref()[0],
)))
}
fn supports_observed_information_curvature(&self) -> bool {
true
}
}
#[derive(Default)]
struct PermanentCandidateErrorModel {
candidate_calls: usize,
}
impl WorkingModel for PermanentCandidateErrorModel {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(test_working_state(beta, curvature))
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
self.candidate_calls += 1;
Err(EstimationError::InvalidSpecification(format!(
"permanent candidate failure under {curvature:?} curvature at beta={:.3e}",
beta.as_ref()[0],
)))
}
}
#[derive(Default)]
struct FirthAcceptedStateFailureModel {
current_state_calls: usize,
candidate_state_calls: usize,
candidate_screen_calls: usize,
}
impl WorkingModel for FirthAcceptedStateFailureModel {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
if beta.as_ref()[0].abs() < 1e-12 {
self.current_state_calls += 1;
Ok(test_working_state(beta, curvature))
} else {
self.candidate_state_calls += 1;
Err(EstimationError::InvalidInput(format!(
"overflow while re-evaluating accepted candidate under {curvature:?} curvature at beta={:.3e}",
beta.as_ref()[0],
)))
}
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
self.candidate_screen_calls += 1;
let mut state = test_working_state(beta, curvature);
state.deviance = 0.5;
state.gradient = array![0.5];
Ok(state)
}
}
#[derive(Default)]
struct ActiveConstraintKktModel;
impl WorkingModel for ActiveConstraintKktModel {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(scalar_working_state(beta, curvature, 1.0, 0.0))
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(scalar_working_state(beta, curvature, 1.0, 0.0))
}
}
struct PlateauStatusModel {
gradient: f64,
current_deviance: f64,
candidate_deviance: f64,
}
impl PlateauStatusModel {
fn state(
beta: &Coefficients,
curvature: HessianCurvatureKind,
gradient: f64,
deviance: f64,
) -> WorkingState {
scalar_working_state(beta, curvature, gradient, deviance)
}
}
impl WorkingModel for PlateauStatusModel {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(Self::state(
beta,
curvature,
self.gradient,
self.current_deviance,
))
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(Self::state(
beta,
curvature,
self.gradient,
self.candidate_deviance,
))
}
}
struct LinearObjectivePlateauModel {
gradient: f64,
}
impl LinearObjectivePlateauModel {
fn state(&self, beta: &Coefficients, curvature: HessianCurvatureKind) -> WorkingState {
let deviance = 1.0 + self.gradient * beta[0];
scalar_working_state(beta, curvature, self.gradient, deviance)
}
}
impl WorkingModel for LinearObjectivePlateauModel {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(self.state(beta, curvature))
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(self.state(beta, curvature))
}
}
#[test]
fn projected_gradient_excludes_near_bound_kkt_forces() {
let gradient = array![0.5, 1e-4];
let beta = array![1e-6, 2.0];
let lower_bounds = array![0.0, f64::NEG_INFINITY];
let norm = projected_gradient_norm(&gradient, &beta, Some(&lower_bounds));
assert!(
norm < 0.01,
"projected gradient should exclude near-bound KKT force (beta=1e-6, lb=0), got {:.6e}",
norm
);
}
#[test]
fn bound_solver_treats_near_bound_positive_grad_as_active() {
let hessian = array![[2.0, 0.0], [0.0, 2.0]];
let gradient = array![1.0, 0.0];
let beta = array![1e-6, 5.0];
let lower_bounds = array![0.0, f64::NEG_INFINITY];
let mut direction = Array1::zeros(2);
let mut active_hint = vec![];
solve_newton_directionwith_lower_bounds(
&hessian,
&gradient,
&beta,
&lower_bounds,
&mut direction,
Some(&mut active_hint),
)
.expect("solve should succeed");
assert!(
active_hint.contains(&0),
"near-bound coeff with positive gradient should be in active set, got {:?}",
active_hint
);
assert!(
(direction[0] - (-1e-6)).abs() < 1e-14,
"direction should snap to bound (lb - beta = -1e-6), got {:.6e}",
direction[0]
);
}
#[test]
fn pirls_converges_at_active_linear_constraint_kkt_point() {
let mut model = ActiveConstraintKktModel;
let options = WorkingModelPirlsOptions {
max_iterations: 3,
convergence_tolerance: 1e-8,
adaptive_kkt_tolerance: None,
max_step_halving: 3,
min_step_size: 0.0,
firth_bias_reduction: false,
coefficient_lower_bounds: None,
linear_constraints: Some(LinearInequalityConstraints {
a: array![[1.0]],
b: array![0.0],
}),
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let summary =
runworking_model_pirls(&mut model, Coefficients::new(array![0.0]), &options, |_| {})
.expect("active-constraint KKT point should be accepted as converged");
assert_eq!(summary.status, PirlsStatus::Converged);
assert!(
summary.lastgradient_norm <= 1e-12,
"KKT-aware stationarity norm should vanish at the constrained optimum, got {:.6e}",
summary.lastgradient_norm
);
let kkt = summary
.constraint_kkt
.expect("linear constraint run should report KKT diagnostics");
assert!(kkt.primal_feasibility <= 1e-12);
assert!(kkt.dual_feasibility <= 1e-12);
assert!(kkt.complementarity <= 1e-12);
assert!(kkt.stationarity <= 1e-12);
}
#[test]
fn certifies_kkt_accepts_biobank_pathological_case() {
let n = 320_000usize;
let p = 20usize;
let g_norm = 1.465e-5;
let tol = 1e-6;
let state = WorkingState {
eta: LinearPredictor::new(Array1::zeros(n)),
gradient: Array1::zeros(p),
hessian: crate::linalg::matrix::SymmetricMatrix::Dense(Array2::zeros((p, p))),
log_likelihood: 0.0,
deviance: 1.0,
penalty_term: 0.0,
firth: FirthDiagnostics::Inactive,
ridge_used: 0.0,
hessian_curvature: HessianCurvatureKind::Fisher,
gradient_natural_scale: 1.0e3,
};
assert!(
state.certifies_kkt(g_norm, tol),
"scale-invariant certificate should accept biobank pathological case"
);
assert!(
!(g_norm < tol),
"this test must witness the failure of the old absolute test; \
otherwise it does not prove the fix"
);
}
#[test]
fn certifies_kkt_is_scale_invariant() {
let n = 1000usize;
let p = 10usize;
let tol = 1e-6;
let g_norm = 1.0;
let natural_scale = 5.0e6;
let mk_state = |g: Array1<f64>, ns: f64| WorkingState {
eta: LinearPredictor::new(Array1::zeros(n)),
gradient: g,
hessian: crate::linalg::matrix::SymmetricMatrix::Dense(Array2::zeros((p, p))),
log_likelihood: 0.0,
deviance: 0.0,
penalty_term: 0.0,
firth: FirthDiagnostics::Inactive,
ridge_used: 0.0,
hessian_curvature: HessianCurvatureKind::Fisher,
gradient_natural_scale: ns,
};
let base = mk_state(Array1::zeros(p), natural_scale);
let scaled = mk_state(Array1::zeros(p), natural_scale * 1000.0);
assert_eq!(
base.certifies_kkt(g_norm, tol),
scaled.certifies_kkt(g_norm * 1000.0, tol),
"KKT classification must be invariant under uniform F → c·F"
);
}
#[test]
fn certifies_kkt_accepts_under_either_bound() {
let n = 100usize;
let p = 5usize;
let tol = 1e-6;
let state_well_scaled = WorkingState {
eta: LinearPredictor::new(Array1::zeros(n)),
gradient: Array1::zeros(p),
hessian: crate::linalg::matrix::SymmetricMatrix::Dense(Array2::zeros((p, p))),
log_likelihood: 0.0,
deviance: 0.0,
penalty_term: 0.0,
firth: FirthDiagnostics::Inactive,
ridge_used: 0.0,
hessian_curvature: HessianCurvatureKind::Fisher,
gradient_natural_scale: 1.0e6,
};
assert!(state_well_scaled.certifies_kkt(0.99e-6 * (1.0 + 1.0e6), tol));
let state_unscaled = WorkingState {
eta: LinearPredictor::new(Array1::zeros(n)),
gradient: Array1::zeros(p),
hessian: crate::linalg::matrix::SymmetricMatrix::Dense(Array2::zeros((p, p))),
log_likelihood: 0.0,
deviance: 0.0,
penalty_term: 0.0,
firth: FirthDiagnostics::Inactive,
ridge_used: 0.0,
hessian_curvature: HessianCurvatureKind::Fisher,
gradient_natural_scale: 0.0,
};
assert!(state_unscaled.certifies_kkt(2.0e-6, tol));
}
#[test]
fn near_stationary_kkt_uses_ten_times_band() {
let n = 100usize;
let p = 4usize;
let tol = 1e-6;
let state = WorkingState {
eta: LinearPredictor::new(Array1::zeros(n)),
gradient: Array1::zeros(p),
hessian: crate::linalg::matrix::SymmetricMatrix::Dense(Array2::zeros((p, p))),
log_likelihood: 0.0,
deviance: 0.0,
penalty_term: 0.0,
firth: FirthDiagnostics::Inactive,
ridge_used: 0.0,
hessian_curvature: HessianCurvatureKind::Fisher,
gradient_natural_scale: 99.0,
};
assert!(state.near_stationary_kkt(9.9e-4, tol));
assert!(!state.near_stationary_kkt(2.0e-3, tol));
assert!(!state.certifies_kkt(9.9e-4, tol));
}
#[test]
fn newton_decrement_correction_upper_bounds_true_decrement() {
let lambda_min = 0.5_f64;
let lambda_lm = 0.25_f64;
let g = ndarray::array![1.0_f64, 1.0];
let true_decrement_sq = g[0].powi(2) / 2.0 + g[1].powi(2) / 0.5;
let damped_decrement_sq =
g[0].powi(2) / (2.0 + lambda_lm) + g[1].powi(2) / (0.5 + lambda_lm);
let correction = 1.0 + lambda_lm / lambda_min;
let upper_bound = damped_decrement_sq * correction;
assert!(
upper_bound >= true_decrement_sq,
"(1 + λ_lm/λ_min)·damped must upper-bound true decrement: \
upper={:.6} true={:.6}",
upper_bound,
true_decrement_sq,
);
assert!(
upper_bound <= 2.0 * true_decrement_sq,
"correction should not be wildly loose: upper={:.6} true={:.6}",
upper_bound,
true_decrement_sq,
);
}
#[test]
fn lm_gain_ratio_accepts_zero_step_at_stationarity() {
let current_penalized: f64 = 9e5;
let predicted_reduction: f64 = 5e-16;
let actual_reduction: f64 = -1e-14;
let noise_floor = current_penalized.abs().max(1.0) * 1e-14;
let rho = if predicted_reduction > noise_floor {
actual_reduction / predicted_reduction
} else if actual_reduction >= -noise_floor {
1.0 } else {
-1.0
};
assert!(
rho > 0.0,
"near-zero reductions should not hard-reject; rho={:.1}, pred={:.2e}, actual={:.2e}, noise={:.2e}",
rho,
predicted_reduction,
actual_reduction,
noise_floor
);
}
#[test]
fn candidate_evaluation_errors_respect_lm_exhaustion_budget() {
let mut model = CandidateEvalFailureModel::default();
let options = WorkingModelPirlsOptions {
max_iterations: 1,
convergence_tolerance: 1e-8,
adaptive_kkt_tolerance: None,
max_step_halving: 5,
min_step_size: 0.0,
firth_bias_reduction: false,
coefficient_lower_bounds: None,
linear_constraints: None,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let err = match runworking_model_pirls(
&mut model,
Coefficients::new(array![0.0]),
&options,
|_| {},
) {
Ok(_) => panic!("candidate evaluation failures should exhaust LM retries and surface"),
Err(err) => err,
};
match err {
EstimationError::PirlsDidNotConverge {
max_iterations,
last_change,
} => {
assert!(
max_iterations == options.max_iterations,
"expected LM exhaustion to surface as PIRLS non-convergence with screening cap"
);
assert!(last_change.is_finite() && last_change > 0.0);
}
other => {
panic!("expected PirlsDidNotConverge from candidate evaluation, got {other:?}")
}
}
assert_eq!(
model.observed_updates, 1,
"the PIRLS iteration should start on observed curvature once"
);
assert_eq!(
model.fisher_updates, 1,
"candidate failure should trigger exactly one observed->Fisher fallback"
);
assert_eq!(
model.observed_candidate_calls, 1,
"observed candidate evaluation should fail once before the Fisher fallback"
);
assert_eq!(
model.fisher_candidate_calls,
options.max_step_halving - 1,
"Fisher candidate evaluation must stop at the configured LM retry budget"
);
}
#[test]
fn permanent_candidate_errors_do_not_trigger_lm_retries() {
let mut model = PermanentCandidateErrorModel::default();
let options = WorkingModelPirlsOptions {
max_iterations: 1,
convergence_tolerance: 1e-8,
adaptive_kkt_tolerance: None,
max_step_halving: 5,
min_step_size: 0.0,
firth_bias_reduction: false,
coefficient_lower_bounds: None,
linear_constraints: None,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let err = match runworking_model_pirls(
&mut model,
Coefficients::new(array![0.0]),
&options,
|_| {},
) {
Ok(_) => panic!("permanent candidate failures should surface immediately"),
Err(err) => err,
};
match err {
EstimationError::InvalidSpecification(message) => {
assert!(
message.contains("permanent candidate failure"),
"expected permanent candidate failure, got {message}"
);
}
other => panic!("expected InvalidSpecification, got {other:?}"),
}
assert_eq!(
model.candidate_calls, 1,
"non-retriable candidate failures should not be re-evaluated under stronger damping"
);
}
#[test]
fn firth_candidate_reevaluation_respects_lm_retry_budget() {
let mut model = FirthAcceptedStateFailureModel::default();
let options = WorkingModelPirlsOptions {
max_iterations: 1,
convergence_tolerance: 1e-8,
adaptive_kkt_tolerance: None,
max_step_halving: 4,
min_step_size: 0.0,
firth_bias_reduction: true,
coefficient_lower_bounds: None,
linear_constraints: None,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let err = match runworking_model_pirls(
&mut model,
Coefficients::new(array![0.0]),
&options,
|_| {},
) {
Ok(_) => panic!("Firth candidate reevaluation failures should not loop indefinitely"),
Err(err) => err,
};
match err {
EstimationError::PirlsDidNotConverge {
max_iterations,
last_change,
} => {
assert_eq!(max_iterations, options.max_iterations);
assert!(last_change.is_finite() && last_change > 0.0);
}
other => panic!("expected PirlsDidNotConverge, got {other:?}"),
}
assert_eq!(model.current_state_calls, 1);
assert_eq!(
model.candidate_screen_calls, options.max_step_halving,
"screening pass should retry until the LM budget is exhausted"
);
assert_eq!(
model.candidate_state_calls, options.max_step_halving,
"Firth accepted-state reevaluation must stop at the configured LM retry budget"
);
}
#[test]
fn plateaued_accepted_step_does_not_report_converged_with_large_projected_gradient() {
let mut model = PlateauStatusModel {
gradient: 5e-5,
current_deviance: 1.0,
candidate_deviance: 1.0 - 1.25e-9,
};
let options = WorkingModelPirlsOptions {
max_iterations: 1,
convergence_tolerance: 1e-6,
adaptive_kkt_tolerance: None,
max_step_halving: 4,
min_step_size: 0.0,
firth_bias_reduction: false,
coefficient_lower_bounds: None,
linear_constraints: None,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let result =
runworking_model_pirls(&mut model, Coefficients::new(array![0.0]), &options, |_| {})
.expect("plateaued accepted step should still return a final state");
assert_eq!(
result.status,
PirlsStatus::MaxIterationsReached,
"projected gradient 5e-5 is well above the near-stationary band and must not be promoted to Converged/Stalled — the candidate step is accepted but the outer iteration counter must run out as MaxIterationsReached, not be silently re-classified"
);
}
#[test]
fn long_constrained_objective_plateau_reports_valid_stall() {
let mut model = LinearObjectivePlateauModel { gradient: -5e-5 };
let options = WorkingModelPirlsOptions {
max_iterations: 25,
convergence_tolerance: 1e-6,
adaptive_kkt_tolerance: None,
max_step_halving: 4,
min_step_size: 0.0,
firth_bias_reduction: false,
coefficient_lower_bounds: None,
linear_constraints: Some(LinearInequalityConstraints {
a: array![[1.0]],
b: array![-100.0],
}),
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let result =
runworking_model_pirls(&mut model, Coefficients::new(array![0.0]), &options, |_| {})
.expect("long constrained objective plateau should preserve the final state");
assert_eq!(
result.status,
PirlsStatus::StalledAtValidMinimum,
"a long monotone objective plateau under explicit constraints is a valid bounded stall, unlike the unconstrained one-step plateau guard above"
);
assert!(
result.iterations < options.max_iterations,
"the long-plateau certificate should exit before exhausting the whole iteration budget"
);
}
#[test]
fn rejected_noise_scale_step_requires_near_stationary_projected_gradient() {
let mut model = PlateauStatusModel {
gradient: 2e-5,
current_deviance: 1.0e6,
candidate_deviance: 1.0e6 + 1.0,
};
let options = WorkingModelPirlsOptions {
max_iterations: 1,
convergence_tolerance: 1e-6,
adaptive_kkt_tolerance: None,
max_step_halving: 1,
min_step_size: 0.0,
firth_bias_reduction: false,
coefficient_lower_bounds: None,
linear_constraints: None,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let result =
runworking_model_pirls(&mut model, Coefficients::new(array![0.0]), &options, |_| {})
.expect("noise-scale rejected step should still preserve the current state");
assert_eq!(
result.status,
PirlsStatus::LmStepSearchExhausted,
"projected gradient 2e-5 exceeds the near-stationary band and must hit the LM-exhaust exit, not be accepted after a noise-scale rejection or fall through to MaxIterationsReached"
);
}
fn assert_deviance_monotone(trace: &[f64], label: &str) {
assert!(
trace.len() >= 2,
"{}: expected at least 2 deviance recordings, got {}",
label,
trace.len()
);
for i in 1..trace.len() {
let prev = trace[i - 1];
let curr = trace[i];
let tol = 1e-8 * prev.abs() + 1e-12;
assert!(
curr <= prev + tol,
"{}: deviance increased at iteration {} -> {}: {:.12e} -> {:.12e} (delta = {:.3e})",
label,
i - 1,
i,
prev,
curr,
curr - prev,
);
}
}
#[test]
fn test_deviance_monotonicity_gaussian() {
assert!(file!().ends_with(".rs"));
let n = 20;
let mut x_data = Array2::<f64>::zeros((n, 2));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = i as f64 / (n - 1) as f64;
x_data[[i, 0]] = 1.0; x_data[[i, 1]] = t; y[i] = 3.0 + 2.0 * t + 0.3 * (((i * 17 + 5) % 11) as f64 / 11.0 - 0.5);
}
let w = Array1::ones(n);
let offset = Array1::zeros(n);
let rho = array![0.0]; let rs = [array![[0.0, 0.0], [0.0, 1.0]]];
let canonical: Vec<crate::construction::CanonicalPenalty> = rs
.iter()
.map(|r| {
let local = r.t().dot(r);
crate::construction::CanonicalPenalty {
root: r.clone(),
col_range: 0..r.ncols(),
total_dim: r.ncols(),
nullity: 0,
local,
prior_mean: Array1::zeros(r.ncols()),
positive_eigenvalues: Vec::new(),
op: None,
}
})
.collect();
let config = PirlsConfig {
likelihood: GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Gaussian,
InverseLink::Standard(StandardLink::Identity),
)),
link_kind: InverseLink::Standard(StandardLink::Identity),
max_iterations: 100,
convergence_tolerance: 1e-8,
firth_bias_reduction: false,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let (result, trace) = capture_pirls_penalized_deviance(|| {
fit_model_for_fixed_rho(
LogSmoothingParamsView::new(rho.view()),
PirlsProblem {
x: x_data.view(),
offset: offset.view(),
y: y.view(),
priorweights: w.view(),
covariate_se: None,
gaussian_fixed_cache: None,
},
PenaltyConfig {
canonical_penalties: &canonical,
balanced_penalty_root: None,
reparam_invariant: None,
p: 2,
coefficient_lower_bounds: None,
linear_constraints_original: None,
penalty_shrinkage_floor: None,
kronecker_factored: None,
},
&config,
None,
)
});
result.expect("Gaussian P-IRLS fit should succeed");
if trace.len() < 2 {
return;
}
assert_deviance_monotone(&trace, "Gaussian");
}
#[test]
fn test_deviance_monotonicity_logistic() {
assert!(file!().ends_with(".rs"));
let n = 30;
let mut x_data = Array2::<f64>::zeros((n, 2));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t = (i as f64 / (n - 1) as f64) * 4.0 - 2.0; x_data[[i, 0]] = 1.0;
x_data[[i, 1]] = t;
let eta = 0.5 + 1.5 * t;
let p = 1.0 / (1.0 + (-eta).exp());
let pseudo_random = ((i * 31 + 7) % 17) as f64 / 17.0;
y[i] = if pseudo_random < p { 1.0 } else { 0.0 };
}
let w = Array1::ones(n);
let offset = Array1::zeros(n);
let rho = array![0.0];
let rs = [array![[0.0, 0.0], [0.0, 1.0]]];
let canonical: Vec<crate::construction::CanonicalPenalty> = rs
.iter()
.map(|r| {
let local = r.t().dot(r);
crate::construction::CanonicalPenalty {
root: r.clone(),
col_range: 0..r.ncols(),
total_dim: r.ncols(),
nullity: 0,
local,
prior_mean: Array1::zeros(r.ncols()),
positive_eigenvalues: Vec::new(),
op: None,
}
})
.collect();
let config = PirlsConfig {
likelihood: GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Binomial,
InverseLink::Standard(StandardLink::Logit),
)),
link_kind: InverseLink::Standard(StandardLink::Logit),
max_iterations: 100,
convergence_tolerance: 1e-8,
firth_bias_reduction: false,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let (result, trace) = capture_pirls_penalized_deviance(|| {
fit_model_for_fixed_rho(
LogSmoothingParamsView::new(rho.view()),
PirlsProblem {
x: x_data.view(),
offset: offset.view(),
y: y.view(),
priorweights: w.view(),
covariate_se: None,
gaussian_fixed_cache: None,
},
PenaltyConfig {
canonical_penalties: &canonical,
balanced_penalty_root: None,
reparam_invariant: None,
p: 2,
coefficient_lower_bounds: None,
linear_constraints_original: None,
penalty_shrinkage_floor: None,
kronecker_factored: None,
},
&config,
None,
)
});
result.expect("Logistic P-IRLS fit should succeed");
assert_deviance_monotone(&trace, "Logistic");
}
#[test]
fn test_deviance_monotonicity_logistic_multiseed() {
let seeds: &[u64] = &[42, 137, 271, 314, 997];
let n = 25;
for &seed in seeds {
let mut x_data = Array2::<f64>::zeros((n, 3));
let mut y = Array1::<f64>::zeros(n);
for i in 0..n {
let t1 = (i as f64 / (n - 1) as f64) * 6.0 - 3.0;
let t2 =
((i as u64).wrapping_mul(seed).wrapping_add(13) % 100) as f64 / 100.0 - 0.5;
x_data[[i, 0]] = 1.0;
x_data[[i, 1]] = t1;
x_data[[i, 2]] = t2;
let eta = -0.3 + 1.0 * t1 + 0.8 * t2;
let p = 1.0 / (1.0 + (-eta).exp());
let hash = (i as u64)
.wrapping_mul(seed)
.wrapping_add(seed >> 2)
.wrapping_mul(2654435761);
let pseudo_uniform = (hash % 10000) as f64 / 10000.0;
y[i] = if pseudo_uniform < p { 1.0 } else { 0.0 };
}
let ones: f64 = y.iter().sum();
if ones < 1.0 {
y[0] = 1.0;
}
if ones > (n as f64 - 1.0) {
y[n - 1] = 0.0;
}
let w = Array1::ones(n);
let offset = Array1::zeros(n);
let rho = array![0.0, 0.0];
let rs = vec![
array![[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [0.0, 0.0, 0.0]],
array![[0.0, 0.0, 0.0], [0.0, 0.0, 0.0], [0.0, 0.0, 1.0]],
];
let canonical: Vec<crate::construction::CanonicalPenalty> = rs
.iter()
.map(|r| {
let local = r.t().dot(r);
crate::construction::CanonicalPenalty {
root: r.clone(),
col_range: 0..r.ncols(),
total_dim: r.ncols(),
nullity: 0,
local,
prior_mean: Array1::zeros(r.ncols()),
positive_eigenvalues: Vec::new(),
op: None,
}
})
.collect();
let config = PirlsConfig {
likelihood: GlmLikelihoodSpec::canonical(LikelihoodSpec::new(
ResponseFamily::Binomial,
InverseLink::Standard(StandardLink::Logit),
)),
link_kind: InverseLink::Standard(StandardLink::Logit),
max_iterations: 100,
convergence_tolerance: 1e-8,
firth_bias_reduction: false,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let (result, trace) = capture_pirls_penalized_deviance(|| {
fit_model_for_fixed_rho(
LogSmoothingParamsView::new(rho.view()),
PirlsProblem {
x: x_data.view(),
offset: offset.view(),
y: y.view(),
priorweights: w.view(),
covariate_se: None,
gaussian_fixed_cache: None,
},
PenaltyConfig {
canonical_penalties: &canonical,
balanced_penalty_root: None,
reparam_invariant: None,
p: 3,
coefficient_lower_bounds: None,
linear_constraints_original: None,
penalty_shrinkage_floor: None,
kronecker_factored: None,
},
&config,
None,
)
});
result.unwrap_or_else(|e| {
panic!("Logistic P-IRLS fit failed for seed {}: {:?}", seed, e)
});
assert_deviance_monotone(&trace, &format!("Logistic(seed={})", seed));
}
}
#[test]
fn solve_newton_direction_implicit_matches_dense_at_k500() {
use crate::terms::closed_form_operator::ClosedFormPenaltyOperator;
use crate::terms::penalty_op::PenaltyOp;
const K: usize = 500;
const D: usize = 4;
let mut state: u64 = 0xDEADBEEF_CAFEBABE;
let mut next = || {
state = state
.wrapping_mul(6364136223846793005)
.wrapping_add(1442695040888963407);
((state >> 11) as f64) / ((1u64 << 53) as f64)
};
let mut centers = Array2::<f64>::zeros((K, D));
for i in 0..K {
for j in 0..D {
centers[[i, j]] = next();
}
}
let op = std::sync::Arc::new(ClosedFormPenaltyOperator::new(
centers.view(),
2,
2,
1,
1.0,
None,
None,
0,
None,
));
let p = op.dim();
assert_eq!(p, K);
let s_dense = op.as_dense();
let mut xtwx = Array2::<f64>::zeros((p, p));
for i in 0..p {
for j in 0..=i {
let v = if i == j {
2.0 + ((i as f64) * 0.07).sin() * 0.3
} else {
(((i as f64 - j as f64) * 0.13).cos()) * 0.02 / (((i + 1) as f64).sqrt())
};
xtwx[[i, j]] = v;
xtwx[[j, i]] = v;
}
}
let xtwx_diag: Array1<f64> = (0..p).map(|i| xtwx[[i, i]]).collect();
let lambda = 0.1_f64;
let ridge = 0.0_f64;
let gradient = Array1::<f64>::from_shape_fn(p, |i| ((i as f64) * 0.31).sin());
let mut h_dense = xtwx.clone();
for i in 0..p {
for j in 0..p {
h_dense[[i, j]] += lambda * s_dense[[i, j]];
}
}
let mut dense_dir = Array1::<f64>::zeros(p);
super::solve_newton_direction_dense(&h_dense, &gradient, &mut dense_dir)
.expect("dense Newton solve should succeed on synthetic SPD");
let xtwx_for_closure = xtwx.clone();
let apply_xtwx = move |v: &Array1<f64>| -> Array1<f64> { xtwx_for_closure.dot(v) };
let op_pen: &dyn PenaltyOp = op.as_ref();
let mut implicit_dir = Array1::<f64>::zeros(p);
super::solve_newton_direction_implicit(
apply_xtwx,
xtwx_diag.view(),
&[],
&[(lambda, op_pen)],
&gradient,
&mut implicit_dir,
ridge,
1e-12,
4 * p,
)
.expect("implicit Newton solve should succeed on synthetic SPD");
let dense_norm: f64 = dense_dir.iter().map(|v| v * v).sum::<f64>().sqrt();
let mut diff_sq = 0.0_f64;
for i in 0..p {
let d = implicit_dir[i] - dense_dir[i];
diff_sq += d * d;
}
let rel = diff_sq.sqrt() / dense_norm.max(1e-300);
assert!(
rel < 1e-9,
"implicit-PCG vs dense-Cholesky Newton direction relative diff {} exceeds 1e-9",
rel
);
}
#[derive(Default)]
struct InnerFisherButObservedSpdAtMode {
observed_post_calls: usize,
}
impl WorkingModel for InnerFisherButObservedSpdAtMode {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
if curvature == HessianCurvatureKind::Observed {
self.observed_post_calls += 1;
}
Ok(scalar_working_state(beta, curvature, 0.0, 0.0))
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(scalar_working_state(beta, curvature, 0.0, 0.0))
}
fn supports_observed_information_curvature(&self) -> bool {
true
}
}
#[test]
fn exported_laplace_observed_exact_when_post_finalization_spd() {
let mut model = InnerFisherButObservedSpdAtMode::default();
let options = WorkingModelPirlsOptions {
max_iterations: 2,
convergence_tolerance: 1e-8,
adaptive_kkt_tolerance: None,
max_step_halving: 3,
min_step_size: 0.0,
firth_bias_reduction: false,
coefficient_lower_bounds: None,
linear_constraints: None,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let summary =
runworking_model_pirls(&mut model, Coefficients::new(array![0.0]), &options, |_| {})
.expect("converged scalar model should produce a result");
assert!(
matches!(
summary.exported_laplace_curvature,
ExportedLaplaceCurvature::ObservedExact
),
"post-convergence Observed-SPD must export ObservedExact, got {:?}",
summary.exported_laplace_curvature
);
assert!(
model.observed_post_calls >= 1,
"post-convergence finalization must call update_with_curvature(Observed) \
at least once to assert SPD inertia"
);
}
#[derive(Default)]
struct CanonicalSurrogateModel;
impl WorkingModel for CanonicalSurrogateModel {
fn update(&mut self, beta: &Coefficients) -> Result<WorkingState, EstimationError> {
self.update_with_curvature(beta, HessianCurvatureKind::Fisher)
}
fn update_with_curvature(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(scalar_working_state(beta, curvature, 0.0, 0.0))
}
fn update_candidate(
&mut self,
beta: &Coefficients,
curvature: HessianCurvatureKind,
) -> Result<WorkingState, EstimationError> {
Ok(scalar_working_state(beta, curvature, 0.0, 0.0))
}
}
#[test]
fn exported_laplace_surrogate_when_observed_unsupported() {
let mut model = CanonicalSurrogateModel;
let options = WorkingModelPirlsOptions {
max_iterations: 2,
convergence_tolerance: 1e-8,
adaptive_kkt_tolerance: None,
max_step_halving: 3,
min_step_size: 0.0,
firth_bias_reduction: false,
coefficient_lower_bounds: None,
linear_constraints: None,
initial_lm_lambda: None,
geodesic_acceleration: false,
arrow_schur: None,
};
let summary =
runworking_model_pirls(&mut model, Coefficients::new(array![0.0]), &options, |_| {})
.expect("canonical surrogate model should converge");
assert!(
matches!(
summary.exported_laplace_curvature,
ExportedLaplaceCurvature::ExpectedInformationSurrogate
),
"model that doesn't support observed information must export \
ExpectedInformationSurrogate (no silent ObservedExact relabel), \
got {:?}",
summary.exported_laplace_curvature
);
}
#[test]
fn dense_xtwx_signed_assembly_preserves_negative_weights() {
let x = array![[1.0, 2.0], [3.0, -1.0], [0.5, 4.0]];
let weights = array![2.0, -3.0, 0.25];
let mut chunk = Array2::<f64>::zeros((0, 0));
let mut got = Array2::<f64>::zeros((2, 2));
PirlsWorkspace::add_dense_xtwx_signed(&weights, &mut chunk, &x, &mut got);
let mut expected = Array2::<f64>::zeros((2, 2));
for i in 0..x.nrows() {
for a in 0..x.ncols() {
for b in 0..x.ncols() {
expected[[a, b]] += weights[i] * x[[i, a]] * x[[i, b]];
}
}
}
for (actual, expected) in got.iter().zip(expected.iter()) {
assert_relative_eq!(*actual, *expected, epsilon = 1e-12);
}
assert!(
got[[0, 0]] < 0.0,
"negative observed-Hessian weights must not be clipped away"
);
}
}
pub(crate) const PIRLS_CACHE_BYTE_BUDGET: usize = 128 * 1024 * 1024;