use super::exact_eval_cache::*;
use super::family::*;
use super::gradient_paths::*;
use super::hessian_paths::*;
use super::install_flex::validate_spec;
use super::*;
#[cfg(test)]
mod tests {
use super::exact_kernel::{
DenestedCubicCell as ExactDenestedCubicCell, ExactCellBranch as ExactCellBranchShared,
LocalSpanCubic, branch_cell as branch_exact_cell, build_denested_partition_cells,
denested_cell_coefficient_partials as exact_denested_cell_coefficient_partials,
global_cubic_from_local as exact_global_cubic_from_local,
transformed_link_cubic as exact_transformed_link_cubic,
};
use super::*;
use crate::custom_family::{CustomFamily, ExactOuterDerivativeOrder};
use ndarray::array;
#[inline]
fn bernoulli_marginal_slope_probit_link() -> InverseLink {
InverseLink::Standard(StandardLink::Probit)
}
fn default_test_family() -> BernoulliMarginalSlopeFamily {
let empty_design = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((0, 0)),
));
BernoulliMarginalSlopeFamily {
y: Arc::new(Array1::zeros(0)),
weights: Arc::new(Array1::zeros(0)),
z: Arc::new(Array1::zeros(0)),
marginal_design: empty_design.clone(),
logslope_design: empty_design,
latent_measure: LatentMeasureKind::StandardNormal,
gaussian_frailty_sd: None,
base_link: InverseLink::Standard(crate::types::StandardLink::Logit),
score_warp: None,
link_dev: None,
policy: crate::resource::ResourcePolicy::default_library(),
cell_moment_lru: Arc::new(exact_kernel::CellMomentLruCache::new(1024)),
cell_moment_cache_stats: Arc::new(exact_kernel::CellMomentCacheStats::default()),
intercept_warm_starts: None,
auto_subsample_phase_counter: Arc::new(std::sync::atomic::AtomicUsize::new(0)),
auto_subsample_last_rho: Arc::new(std::sync::Mutex::new(None)),
}
}
fn empty_termspec() -> TermCollectionSpec {
TermCollectionSpec {
linear_terms: vec![],
random_effect_terms: vec![],
smooth_terms: vec![],
}
}
fn dummy_blockspec(p: usize, n_rows: usize) -> ParameterBlockSpec {
ParameterBlockSpec {
name: "dummy".to_string(),
design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::<f64>::zeros((n_rows, p)),
)),
offset: Array1::zeros(n_rows),
penalties: vec![],
nullspace_dims: vec![],
initial_log_lambdas: Array1::zeros(0),
initial_beta: Some(Array1::zeros(p)),
gauge_priority: 100,
jacobian_callback: None,
stacked_design: None,
stacked_offset: None,
}
}
fn dummy_block_state(beta: Array1<f64>, n_rows: usize) -> ParameterBlockState {
ParameterBlockState {
beta,
eta: Array1::zeros(n_rows),
}
}
fn flex_hessian_matvec_fixture(
n: usize,
) -> Result<
(
BernoulliMarginalSlopeFamily,
Vec<ParameterBlockState>,
BernoulliMarginalSlopeExactEvalCache,
Array1<f64>,
),
String,
> {
let z = Array1::from_iter((0..n).map(|i| {
let t = (i as f64 + 0.5) / n as f64;
(10.0 * t).sin() + 0.3 * (31.0 * t).cos()
}));
let y = Array1::from_iter((0..n).map(|i| if (i * 37 + 11) % 101 < 43 { 1.0 } else { 0.0 }));
let weights = Array1::from_iter((0..n).map(|i| 0.75 + 0.5 * ((i % 7) as f64) / 6.0));
let design = Array2::from_shape_fn((n, 2), |(row, col)| match col {
0 => 1.0,
1 => z[row],
_ => unreachable!(),
});
let cfg = DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
};
let score_prepared = build_score_warp_deviation_block_from_seed(&z, &cfg)?;
let q_seed = Array1::from_iter(z.iter().map(|zi| 0.05 + 0.2 * zi));
let link_seed = padded_deviation_seed(&q_seed, 1.0, 0.5);
let link_prepared = build_link_deviation_block_from_knots_design_seed_and_weights(
&link_seed, &q_seed, &cfg,
)?;
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(y),
weights: Arc::new(weights),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
design.clone(),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(design)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let marginal_beta = array![0.05, 0.08];
let logslope_beta = array![-0.15, 0.04];
let marginal_eta =
Array1::from_iter(z.iter().map(|zi| marginal_beta[0] + marginal_beta[1] * zi));
let logslope_eta =
Array1::from_iter(z.iter().map(|zi| logslope_beta[0] + logslope_beta[1] * zi));
let states = vec![
ParameterBlockState {
beta: marginal_beta,
eta: marginal_eta,
},
ParameterBlockState {
beta: logslope_beta,
eta: logslope_eta,
},
ParameterBlockState {
beta: Array1::zeros(score_prepared.block.design.ncols()),
eta: Array1::zeros(n),
},
ParameterBlockState {
beta: Array1::zeros(link_prepared.block.design.ncols()),
eta: Array1::zeros(n),
},
];
let mut cache = family.build_exact_eval_cache(&states)?;
cache.row_primary_hessians = family.build_row_primary_hessian_cache(&states, &cache)?;
let direction = Array1::from_iter((0..cache.slices.total).map(|j| {
let x = j as f64 + 1.0;
0.03 * x.sin() + 0.01 * (0.37 * x).cos()
}));
Ok((family, states, cache, direction))
}
fn assert_allclose_relative(actual: &Array1<f64>, expected: &Array1<f64>, tol: f64) {
assert_eq!(actual.len(), expected.len());
for (idx, (&a, &e)) in actual.iter().zip(expected.iter()).enumerate() {
let denom = a.abs().max(e.abs()).max(1.0);
let rel = (a - e).abs() / denom;
assert!(
rel <= tol,
"entry {idx}: actual={a:.17e}, expected={e:.17e}, rel={rel:.3e}, tol={tol:.3e}"
);
}
}
#[test]
fn cross_block_identifiability_anchor_wider_than_candidate_no_false_alias() {
use crate::linalg::matrix::{DenseDesignMatrix, DesignMatrix};
let n = 64usize;
let z = Array1::from_iter((0..n).map(|i| {
let t = (i as f64) / (n as f64 - 1.0);
-1.5 + 3.0 * t
}));
let weights = Array1::from_elem(n, 1.0);
let mut anchor_dense = ndarray::Array2::<f64>::zeros((n, 20));
for i in 0..n {
for j in 0..20 {
let x = (i as f64 * 0.13 + j as f64 * 0.91 + 1.0).sin();
let y = (i as f64 * 0.07 - j as f64 * 0.31 + 2.0).cos();
anchor_dense[[i, j]] = 0.5 * x + 0.5 * y;
}
}
let anchor_design = DesignMatrix::Dense(DenseDesignMatrix::from(anchor_dense.clone()));
let link_cfg = DeviationBlockConfig {
num_internal_knots: 2,
..DeviationBlockConfig::default()
};
let q0_seed = Array1::from_iter(z.iter().map(|zi| 0.05 + 0.4 * zi));
let link_seed = padded_deviation_seed(&q0_seed, 1.0, 0.5);
let mut link_prepared = build_link_deviation_block_from_knots_design_seed_and_weights(
&link_seed, &q0_seed, &link_cfg,
)
.expect("link-dev fixture");
let p_before = link_prepared.runtime.basis_dim();
assert!(p_before > 0, "fixture must have positive basis_dim");
use super::deviation_runtime::ParametricAnchorBlock;
install_compiled_flex_block_into_runtime(
&mut link_prepared,
&q0_seed,
&link_cfg,
&[(&anchor_design, ParametricAnchorBlock::Marginal)],
&[],
&weights,
)
.expect("wider anchor should not false-positive full alias");
let p_after = link_prepared.runtime.basis_dim();
assert!(
p_after > 0,
"new algorithm must keep strictly positive basis_dim when (I-P_A)C has rank > 0"
);
let new_design = link_prepared
.runtime
.design_at_training_with_residual(&q0_seed)
.expect("design after orthogonalisation");
assert_eq!(new_design.nrows(), n);
assert_eq!(new_design.ncols(), p_after);
let cross = anchor_dense.t().dot(&new_design);
let max_abs = cross.iter().fold(0.0_f64, |acc, &v| acc.max(v.abs()));
let anchor_norm = anchor_dense.iter().map(|v| v * v).sum::<f64>().sqrt();
let cand_norm = new_design.iter().map(|v| v * v).sum::<f64>().sqrt();
let scale = (anchor_norm * cand_norm).max(1.0);
assert!(
max_abs <= 1.0e-9 * scale,
"Aᵀ C̃ should be at noise floor; max|.|={max_abs:.3e}, scale={scale:.3e}",
);
}
#[test]
fn cross_block_identifiability_minimal_counterexample_keeps_orthogonal_complement() {
use crate::linalg::matrix::{DenseDesignMatrix, DesignMatrix};
let n = 64usize;
let z = Array1::from_iter((0..n).map(|i| {
let t = (i as f64) / (n as f64 - 1.0);
-1.5 + 3.0 * t
}));
let weights = Array1::from_elem(n, 1.0);
let link_cfg = DeviationBlockConfig {
num_internal_knots: 3,
..DeviationBlockConfig::default()
};
let q0_seed = Array1::from_iter(z.iter().map(|zi| 0.05 + 0.4 * zi));
let link_seed = padded_deviation_seed(&q0_seed, 1.0, 0.5);
let mut link_prepared = build_link_deviation_block_from_knots_design_seed_and_weights(
&link_seed, &q0_seed, &link_cfg,
)
.expect("link-dev fixture");
let p_before = link_prepared.runtime.basis_dim();
assert!(
p_before >= 2,
"fixture must have at least two basis columns"
);
let candidate_design = link_prepared
.runtime
.design(&q0_seed)
.expect("candidate training-row design");
let mut anchor_dense = ndarray::Array2::<f64>::zeros((n, 1));
anchor_dense
.column_mut(0)
.assign(&candidate_design.column(0));
let anchor_design = DesignMatrix::Dense(DenseDesignMatrix::from(anchor_dense.clone()));
use super::deviation_runtime::ParametricAnchorBlock;
install_compiled_flex_block_into_runtime(
&mut link_prepared,
&q0_seed,
&link_cfg,
&[(&anchor_design, ParametricAnchorBlock::Marginal)],
&[],
&weights,
)
.expect("minimal counterexample must keep p_before - 1 directions, not collapse to 0");
let p_after = link_prepared.runtime.basis_dim();
assert_eq!(
p_after,
p_before - 1,
"(I−P_A)C should keep exactly {} directions when one column of C is reproduced by A; got {}",
p_before - 1,
p_after,
);
let new_design = link_prepared
.runtime
.design_at_training_with_residual(&q0_seed)
.expect("design after orthogonalisation");
let cross = anchor_dense.t().dot(&new_design);
let max_abs = cross.iter().fold(0.0_f64, |acc, &v| acc.max(v.abs()));
let anchor_norm = anchor_dense.iter().map(|v| v * v).sum::<f64>().sqrt();
let cand_norm = new_design.iter().map(|v| v * v).sum::<f64>().sqrt();
let scale = (anchor_norm * cand_norm).max(1.0);
assert!(
max_abs <= 1.0e-9 * scale,
"AᵀC̃ should be at noise floor after residualisation; max|.|={max_abs:.3e}, scale={scale:.3e}",
);
}
#[test]
fn cross_block_identifiability_true_alias_returns_fully_aliased_outcome() {
use crate::linalg::matrix::{DenseDesignMatrix, DesignMatrix};
let n = 64usize;
let z = Array1::from_iter((0..n).map(|i| {
let t = (i as f64) / (n as f64 - 1.0);
-1.5 + 3.0 * t
}));
let weights = Array1::from_elem(n, 1.0);
let link_cfg = DeviationBlockConfig {
num_internal_knots: 2,
..DeviationBlockConfig::default()
};
let q0_seed = Array1::from_iter(z.iter().map(|zi| 0.05 + 0.4 * zi));
let link_seed = padded_deviation_seed(&q0_seed, 1.0, 0.5);
let mut link_prepared = build_link_deviation_block_from_knots_design_seed_and_weights(
&link_seed, &q0_seed, &link_cfg,
)
.expect("link-dev fixture");
let candidate_design = link_prepared
.runtime
.design(&q0_seed)
.expect("candidate training-row design");
let anchor_design = DesignMatrix::Dense(DenseDesignMatrix::from(candidate_design));
use super::deviation_runtime::ParametricAnchorBlock;
let outcome = install_compiled_flex_block_into_runtime(
&mut link_prepared,
&q0_seed,
&link_cfg,
&[(&anchor_design, ParametricAnchorBlock::Marginal)],
&[],
&weights,
)
.expect("true alias must produce a structured FullyAliased outcome");
match outcome {
FlexCompileOutcome::FullyAliased { reason } => {
assert!(
reason.contains("zero directions remaining"),
"expected FullyAliased reason mentioning 'zero directions remaining', got: {reason}",
);
}
FlexCompileOutcome::Reparameterised => {
panic!("expected FullyAliased outcome but got Reparameterised");
}
}
}
#[test]
fn cross_block_identifiability_flex_anchor_true_alias_returns_fully_aliased() {
let n = 64usize;
let z = Array1::from_iter((0..n).map(|i| {
let t = (i as f64) / (n as f64 - 1.0);
-1.5 + 3.0 * t
}));
let weights = Array1::from_elem(n, 1.0);
let link_cfg = DeviationBlockConfig {
num_internal_knots: 2,
..DeviationBlockConfig::default()
};
let q0_seed = Array1::from_iter(z.iter().map(|zi| 0.05 + 0.4 * zi));
let link_seed = padded_deviation_seed(&q0_seed, 1.0, 0.5);
let mut link_prepared = build_link_deviation_block_from_knots_design_seed_and_weights(
&link_seed, &q0_seed, &link_cfg,
)
.expect("link-dev fixture");
let flex_anchor_design = link_prepared
.runtime
.design(&q0_seed)
.expect("candidate training-row design");
let outcome = install_compiled_flex_block_into_runtime(
&mut link_prepared,
&q0_seed,
&link_cfg,
&[],
&[&flex_anchor_design],
&weights,
)
.expect("flex-anchor true alias must produce a structured FullyAliased outcome");
match outcome {
FlexCompileOutcome::FullyAliased { reason } => {
assert!(
reason.contains("zero directions remaining"),
"expected FullyAliased reason mentioning 'zero directions remaining', got: {reason}",
);
}
FlexCompileOutcome::Reparameterised => {
panic!(
"FlexEvaluation anchor whose span covers the candidate must yield FullyAliased; \
got Reparameterised, indicating the FlexEvaluation arm was silently skipped",
);
}
}
}
#[test]
fn cross_block_identifiability_outcome_fully_aliased_extracts_reason() {
let outcome = FlexCompileOutcome::FullyAliased {
reason: "candidate has zero directions remaining after residualisation".to_string(),
};
match outcome {
FlexCompileOutcome::FullyAliased { reason } => {
assert!(reason.contains("zero directions remaining"));
}
FlexCompileOutcome::Reparameterised => {
panic!("constructed FullyAliased; cannot pattern-match as Reparameterised")
}
}
}
#[test]
fn cross_block_identifiability_partial_alias_keeps_residual_rank() {
use crate::faer_ndarray::FaerQr;
use crate::linalg::matrix::{DenseDesignMatrix, DesignMatrix};
let n = 96usize;
let z = Array1::from_iter((0..n).map(|i| {
let t = (i as f64) / (n as f64 - 1.0);
-1.5 + 3.0 * t
}));
let weights = Array1::from_elem(n, 1.0);
let link_cfg = DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
};
let q0_seed = Array1::from_iter(z.iter().map(|zi| 0.03 + 0.42 * zi));
let link_seed = padded_deviation_seed(&q0_seed, 1.0, 0.5);
let mut link_prepared = build_link_deviation_block_from_knots_design_seed_and_weights(
&link_seed, &q0_seed, &link_cfg,
)
.expect("link-dev fixture");
let candidate_design = link_prepared
.runtime
.design(&q0_seed)
.expect("candidate training-row design");
let (q, _r) = candidate_design.qr().expect("thin QR of candidate");
let p_c = candidate_design.ncols();
let k_alias = (p_c / 2).max(1);
assert!(
k_alias < p_c,
"partial-alias test needs p_c > k_alias + 1, got p_c={p_c}, k_alias={k_alias}",
);
let p_a = k_alias + 1;
let mut anchor_dense = ndarray::Array2::<f64>::zeros((n, p_a));
for j in 0..k_alias {
anchor_dense.column_mut(j).assign(&q.column(j));
}
let mut extra = Array1::<f64>::zeros(n);
for i in 0..n {
extra[i] = ((i as f64 * 0.17).sin() + (i as f64 * 0.41).cos()) * 0.5;
}
let q_t_extra = q.t().dot(&extra);
let extra_orth = &extra - &q.dot(&q_t_extra);
anchor_dense.column_mut(k_alias).assign(&extra_orth);
let anchor_design = DesignMatrix::Dense(DenseDesignMatrix::from(anchor_dense.clone()));
use super::deviation_runtime::ParametricAnchorBlock;
let p_before = link_prepared.runtime.basis_dim();
install_compiled_flex_block_into_runtime(
&mut link_prepared,
&q0_seed,
&link_cfg,
&[(&anchor_design, ParametricAnchorBlock::Marginal)],
&[],
&weights,
)
.expect("partial alias must keep the surviving rank");
let p_after = link_prepared.runtime.basis_dim();
assert_eq!(
p_after,
p_before - k_alias,
"partial alias should drop exactly the {} aliased directions; got {} -> {}",
k_alias,
p_before,
p_after,
);
let new_design = link_prepared
.runtime
.design_at_training_with_residual(&q0_seed)
.expect("design after orthogonalisation");
let cross = anchor_dense.t().dot(&new_design);
let max_abs = cross.iter().fold(0.0_f64, |acc, &v| acc.max(v.abs()));
let anchor_norm = anchor_dense.iter().map(|v| v * v).sum::<f64>().sqrt();
let cand_norm = new_design.iter().map(|v| v * v).sum::<f64>().sqrt();
let scale = (anchor_norm * cand_norm).max(1.0);
assert!(
max_abs <= 1.0e-9 * scale,
"AᵀC̃ should be at noise floor; max|.|={max_abs:.3e}, scale={scale:.3e}",
);
}
#[test]
fn flex_hessian_matvec_matches_dense_hessian() {
assert!(file!().ends_with(".rs"));
let (family, states, cache, direction) =
flex_hessian_matvec_fixture(96).expect("flex Hv fixture");
let dense = family
.exact_newton_joint_hessian_dense_from_cache(&states, &cache)
.expect("dense Hessian");
let from_matvec = family
.exact_newton_joint_hessian_matvec_from_cache(&direction, &states, &cache)
.expect("parallel chunked Hv");
let from_dense = dense.dot(&direction);
assert_allclose_relative(&from_matvec, &from_dense, 1.0e-13);
}
#[test]
fn row_primary_third_trace_many_matches_single_direction_contracts() {
let (family, states, cache, direction_a) =
flex_hessian_matvec_fixture(24).expect("flex trace fixture");
let direction_b = Array1::from_iter((0..cache.slices.total).map(|j| {
let x = j as f64 + 0.5;
0.02 * (0.7 * x).sin() - 0.015 * (0.31 * x).cos()
}));
let r = cache.primary.total;
let gram = (0..r * r)
.map(|idx| {
let x = idx as f64 + 1.0;
0.03 * x.sin() + 0.01 * (0.17 * x).cos()
})
.collect::<Vec<_>>();
for row in 0..6 {
let row_ctx = BernoulliMarginalSlopeFamily::row_ctx(&cache, row);
let row_dirs = vec![
family
.row_primary_direction_from_flat(
row,
&cache.slices,
&cache.primary,
&direction_a,
)
.expect("row direction a"),
family
.row_primary_direction_from_flat(
row,
&cache.slices,
&cache.primary,
&direction_b,
)
.expect("row direction b"),
];
let many = family
.row_primary_third_trace_many_with_moments(
row, &states, &cache, row_ctx, &row_dirs, &gram,
)
.expect("many-direction row trace");
for (dir_idx, row_dir) in row_dirs.iter().enumerate() {
let third = family
.row_primary_third_contracted_recompute(row, &states, &cache, row_ctx, row_dir)
.expect("single-direction third contraction");
let single =
BernoulliMarginalSlopeFamily::row_primary_trace_contract(&third, &gram);
let denom = single.abs().max(many[dir_idx].abs()).max(1.0);
let rel = (single - many[dir_idx]).abs() / denom;
assert!(
rel <= 1.0e-12,
"row {row} dir {dir_idx}: many={:.17e} single={:.17e} rel={:.3e}",
many[dir_idx],
single,
rel
);
}
}
}
#[test]
fn bernoulli_margslope_warm_start_cache_persists_across_eval_cache_builds() {
let n = 256usize;
let z = Array1::from_iter((0..n).map(|i| {
let t = (i as f64 + 0.5) / n as f64;
(12.0 * t).sin() + 0.25 * (37.0 * t).cos()
}));
let y = Array1::from_iter((0..n).map(|i| if i % 3 == 0 { 1.0 } else { 0.0 }));
let weights = Array1::ones(n);
let cfg = DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
};
let score_prepared = build_score_warp_deviation_block_from_seed(&z, &cfg)
.expect("score-warp deviation block");
let q_seed = Array1::from_iter(z.iter().map(|zi| 0.1 + 0.25 * zi));
let link_seed = padded_deviation_seed(&q_seed, 1.0, 0.5);
let link_prepared = build_link_deviation_block_from_knots_design_seed_and_weights(
&link_seed, &q_seed, &cfg,
)
.expect("link-wiggle deviation block");
let cache = new_intercept_warm_start_cache(n);
let make_family =
|cache: Option<Arc<BernoulliInterceptWarmStartCache>>| BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::ones((n, 1)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::ones((n, 1)),
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
intercept_warm_starts: cache,
..default_test_family()
};
let marginal_eta = Array1::from_iter((0..n).map(|i| 0.15 * ((i as f64) * 0.001).sin()));
let slope_eta = Array1::from_iter((0..n).map(|i| 0.35 + 0.02 * ((i as f64) * 0.003).cos()));
let states = vec![
ParameterBlockState {
beta: array![0.0],
eta: marginal_eta,
},
ParameterBlockState {
beta: array![0.0],
eta: slope_eta,
},
ParameterBlockState {
beta: Array1::zeros(score_prepared.block.design.ncols()),
eta: Array1::zeros(n),
},
ParameterBlockState {
beta: Array1::zeros(link_prepared.block.design.ncols()),
eta: Array1::zeros(n),
},
];
let warm_family = make_family(Some(Arc::clone(&cache)));
let first = warm_family
.build_exact_eval_cache(&states)
.expect("first warm eval cache");
let nan_bits = f64::NAN.to_bits();
for slot in cache.intercept_value.iter() {
let bits = slot.load(Ordering::Relaxed);
let v = f64::from_bits(bits);
assert!(
v.is_finite(),
"cache slot should be populated with converged intercept after first build"
);
assert_ne!(bits, nan_bits);
}
let second = warm_family
.build_exact_eval_cache(&states)
.expect("second warm eval cache");
let cold_family = make_family(None);
let cold = cold_family
.build_exact_eval_cache(&states)
.expect("cold reference eval cache");
for ((warm_a, warm_b), cold_ctx) in first
.row_contexts
.iter()
.zip(second.row_contexts.iter())
.zip(cold.row_contexts.iter())
{
assert!((warm_a.intercept - cold_ctx.intercept).abs() < 1e-9);
assert!((warm_b.intercept - cold_ctx.intercept).abs() < 1e-9);
}
}
#[test]
fn bernoulli_margslope_flex_ll_early_exit_is_exact_or_provably_rejected() {
let n = 24_000usize;
let z = Array1::from_iter((0..n).map(|i| {
let t = (i as f64 + 0.5) / n as f64;
(10.0 * t).sin() + 0.2 * (29.0 * t).cos()
}));
let y = Array1::from_iter((0..n).map(|i| if (i * 17 + 3) % 7 >= 4 { 1.0 } else { 0.0 }));
let weights = Array1::ones(n);
let cfg = DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
};
let score_prepared = build_score_warp_deviation_block_from_seed(&z, &cfg)
.expect("score-warp deviation block");
let q_seed = Array1::from_iter(z.iter().map(|zi| 0.1 + 0.2 * zi));
let link_seed = padded_deviation_seed(&q_seed, 1.0, 0.5);
let link_prepared = build_link_deviation_block_from_knots_design_seed_and_weights(
&link_seed, &q_seed, &cfg,
)
.expect("link-wiggle deviation block");
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(y),
weights: Arc::new(weights),
z: Arc::new(z),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::ones((n, 1)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::ones((n, 1)),
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let states = vec![
ParameterBlockState {
beta: array![0.0],
eta: Array1::from_iter((0..n).map(|i| 0.08 * ((i as f64) * 0.002).sin())),
},
ParameterBlockState {
beta: array![0.0],
eta: Array1::from_iter((0..n).map(|i| 0.30 + 0.03 * ((i as f64) * 0.003).cos())),
},
ParameterBlockState {
beta: Array1::zeros(score_prepared.block.design.ncols()),
eta: Array1::zeros(n),
},
ParameterBlockState {
beta: Array1::zeros(link_prepared.block.design.ncols()),
eta: Array1::zeros(n),
},
];
let exact = family
.log_likelihood_only_with_options(&states, &BlockwiseFitOptions::default())
.expect("exact full FLEX ll");
let mut permissive = BlockwiseFitOptions::default();
permissive.early_exit_threshold = Some((-exact) + 1.0);
let accepted = family
.log_likelihood_only_with_options(&states, &permissive)
.expect("permissive threshold should compute full FLEX ll");
let rel = ((accepted - exact) / exact.abs().max(1.0)).abs();
assert!(
rel < 1e-10,
"accepted early-exit-enabled FLEX LL {accepted} differs from exact {exact} by rel {rel}"
);
let mut rejecting = BlockwiseFitOptions::default();
rejecting.early_exit_threshold = Some(1e-6);
let err = family
.log_likelihood_only_with_options(&states, &rejecting)
.expect_err("tight threshold should reject before the full FLEX row sweep");
assert!(
err.contains("line-search rejected early"),
"unexpected early-exit error: {err}"
);
}
#[test]
fn row_primary_fourth_contracted_rejects_bad_direction_lengths() {
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![1.0]),
weights: Arc::new(array![1.0]),
z: Arc::new(array![0.25]),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((1, 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((1, 0)),
)),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: Array1::zeros(0),
eta: array![0.15],
},
ParameterBlockState {
beta: Array1::zeros(0),
eta: array![0.2],
},
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(0, &block_states, None, true)
.expect("row context");
let bad_dir = array![1.0];
let good_dir = array![0.0, 1.0];
let err = family
.row_primary_fourth_contracted_recompute_ordered(
0,
&block_states,
&cache,
&row_ctx,
&bad_dir,
&good_dir,
)
.expect_err("bad direction length should be rejected before indexing");
assert!(
err.contains("direction lengths (1,2) != 2"),
"unexpected error: {err}"
);
}
fn base_spec(
y: Array1<f64>,
weights: Array1<f64>,
z: Array1<f64>,
) -> BernoulliMarginalSlopeTermSpec {
let n = y.len();
BernoulliMarginalSlopeTermSpec {
y,
weights,
z,
base_link: bernoulli_marginal_slope_probit_link(),
marginalspec: empty_termspec(),
logslopespec: empty_termspec(),
marginal_offset: Array1::zeros(n),
logslope_offset: Array1::zeros(n),
frailty: FrailtySpec::None,
score_warp: None,
link_dev: None,
latent_z_policy: LatentZPolicy::default(),
}
}
#[test]
fn bernoulli_marginal_link_map_zeroes_derivatives_on_clamped_tails() {
let link = bernoulli_marginal_slope_probit_link();
let lower = bernoulli_marginal_link_map(&link, -8.0).expect("lower tail map");
let upper = bernoulli_marginal_link_map(&link, 8.0).expect("upper tail map");
let lower_q = standard_normal_quantile(BERNOULLI_LINK_PROBABILITY_EPS).unwrap();
let upper_q = standard_normal_quantile(1.0 - BERNOULLI_LINK_PROBABILITY_EPS).unwrap();
assert_eq!(lower.mu, BERNOULLI_LINK_PROBABILITY_EPS);
assert_eq!(upper.mu, 1.0 - BERNOULLI_LINK_PROBABILITY_EPS);
assert!((lower.q - lower_q).abs() < 1e-12);
assert!((upper.q - upper_q).abs() < 1e-12);
assert_eq!([lower.mu1, lower.mu2, lower.mu3, lower.mu4], [0.0; 4]);
assert_eq!([upper.mu1, upper.mu2, upper.mu3, upper.mu4], [0.0; 4]);
assert_eq!([lower.q1, lower.q2, lower.q3, lower.q4], [0.0; 4]);
assert_eq!([upper.q1, upper.q2, upper.q3, upper.q4], [0.0; 4]);
}
#[test]
fn rigid_transformed_gradient_matches_negative_log_likelihood_derivative() {
let link = bernoulli_marginal_slope_probit_link();
let eta = 0.25;
let g = -0.15;
let z = 0.7;
let y = 1.0;
let weight = 1.3;
let probit_scale = 1.0;
let objective = |eta_value: f64, g_value: f64| {
let marginal = bernoulli_marginal_link_map(&link, eta_value).unwrap();
let kernel =
RigidProbitKernel::new(marginal.q, g_value, z, y, weight, probit_scale).unwrap();
-weight * kernel.logcdf
};
let marginal = bernoulli_marginal_link_map(&link, eta).expect("marginal map");
let kernel =
RigidProbitKernel::new(marginal.q, g, z, y, weight, probit_scale).expect("kernel");
let grad = rigid_transformed_gradient(marginal, &kernel);
let step = 1e-6;
let finite_eta = (objective(eta + step, g) - objective(eta - step, g)) / (2.0 * step);
let finite_g = (objective(eta, g + step) - objective(eta, g - step)) / (2.0 * step);
assert!(
(grad[0] - finite_eta).abs() < 1e-7,
"eta gradient {} != finite difference {}",
grad[0],
finite_eta
);
assert!(
(grad[1] - finite_g).abs() < 1e-7,
"g gradient {} != finite difference {}",
grad[1],
finite_g
);
assert!(
grad[0] < 0.0,
"y=1 probit nll should decrease as eta increases"
);
}
fn pair_distance(lhs: (f64, f64), rhs: (f64, f64)) -> f64 {
(lhs.0 - rhs.0).abs() + (lhs.1 - rhs.1).abs()
}
fn make_rigid_test_family(n: usize) -> BernoulliMarginalSlopeFamily {
let y: Array1<f64> =
Array1::from_iter((0..n).map(|i| if (i * 31 + 7) % 5 >= 3 { 1.0 } else { 0.0 }));
let weights: Array1<f64> =
Array1::from_iter((0..n).map(|i| 0.5 + ((i * 13 + 4) % 7) as f64 * 0.1));
let z: Array1<f64> = Array1::from_iter(
(0..n).map(|i| -1.5 + 3.0 * (((i * 17 + 5) % n) as f64 + 0.5) / (n as f64)),
);
let ones_col = Array2::from_shape_fn((n, 1), |_| 1.0);
BernoulliMarginalSlopeFamily {
y: Arc::new(y),
weights: Arc::new(weights),
z: Arc::new(z),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
ones_col.clone(),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(ones_col)),
..default_test_family()
}
}
fn rigid_block_states(
family: &BernoulliMarginalSlopeFamily,
q: f64,
b: f64,
) -> Vec<ParameterBlockState> {
let n = family.y.len();
vec![
ParameterBlockState {
beta: array![q],
eta: Array1::from_elem(n, q),
},
ParameterBlockState {
beta: array![b],
eta: Array1::from_elem(n, b),
},
]
}
#[test]
fn bernoulli_log_likelihood_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_rigid_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let baseline = family
.log_likelihood_only(&states)
.expect("baseline ll (no subsample)");
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full_mask = family
.log_likelihood_only_with_options(&states, &opts_full)
.expect("ll with mask=full");
let rel = ((with_full_mask - baseline) / baseline.abs().max(1.0)).abs();
assert!(
rel < 1e-12,
"subsample(mask=full) {} differs from baseline {} by rel {}",
with_full_mask,
baseline,
rel
);
}
#[test]
fn bernoulli_log_likelihood_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_rigid_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.log_likelihood_only_with_options(&states, &opts_half)
.expect("ll with mask=even");
let mut opts_even_unscaled = BlockwiseFitOptions::default();
opts_even_unscaled.outer_score_subsample = Some(Arc::new(
OuterScoreSubsample::with_uniform_weight(even_mask, m, 0, 1.0),
));
let raw_even_sum = family
.log_likelihood_only_with_options(&states, &opts_even_unscaled)
.expect("raw even-row ll sum");
let expected_scaled = (n as f64 / m as f64) * raw_even_sum;
let rel = ((scaled - expected_scaled) / expected_scaled.abs().max(1.0)).abs();
assert!(
rel < 1e-12,
"scaled {} != 2*even_sum {} (rel {})",
scaled,
expected_scaled,
rel
);
let baseline = family.log_likelihood_only(&states).expect("baseline ll");
let ht_rel = ((scaled - baseline) / baseline.abs().max(1.0)).abs();
assert!(
ht_rel < 0.05,
"Horvitz-Thompson scaled {} not near baseline {} (rel {})",
scaled,
baseline,
ht_rel
);
}
fn make_sigma_aware_test_family(n: usize) -> BernoulliMarginalSlopeFamily {
let mut family = make_rigid_test_family(n);
family.gaussian_frailty_sd = Some(0.7);
family
}
fn rel_diff_array1(a: &Array1<f64>, b: &Array1<f64>) -> f64 {
let mut max = 0.0f64;
for i in 0..a.len() {
let d = (a[i] - b[i]).abs() / b[i].abs().max(1.0);
if d > max {
max = d;
}
}
max
}
fn rel_diff_array2(a: &Array2<f64>, b: &Array2<f64>) -> f64 {
let mut max = 0.0f64;
for ((i, j), &av) in a.indexed_iter() {
let bv = b[[i, j]];
let d = (av - bv).abs() / bv.abs().max(1.0);
if d > max {
max = d;
}
}
max
}
#[test]
fn bernoulli_sigma_psi_terms_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_sigma_aware_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let specs = vec![dummy_blockspec(1, n), dummy_blockspec(1, n)];
let baseline = family
.sigma_exact_joint_psi_terms(&states, &specs)
.expect("baseline psi terms")
.expect("baseline some");
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full = family
.sigma_exact_joint_psi_terms_with_options(&states, &specs, &opts_full)
.expect("psi terms with full mask")
.expect("some");
let obj_rel = ((with_full.objective_psi - baseline.objective_psi)
/ baseline.objective_psi.abs().max(1.0))
.abs();
assert!(obj_rel < 1e-12, "objective_psi rel {}", obj_rel);
let score_rel = rel_diff_array1(&with_full.score_psi, &baseline.score_psi);
assert!(score_rel < 1e-12, "score_psi rel {}", score_rel);
let h_full = with_full
.hessian_psi_operator
.as_ref()
.expect("op")
.to_dense();
let h_baseline = baseline
.hessian_psi_operator
.as_ref()
.expect("op")
.to_dense();
let h_rel = rel_diff_array2(&h_full, &h_baseline);
assert!(h_rel < 1e-12, "hessian rel {}", h_rel);
}
#[test]
fn bernoulli_sigma_psi_terms_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_sigma_aware_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let specs = vec![dummy_blockspec(1, n), dummy_blockspec(1, n)];
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.sigma_exact_joint_psi_terms_with_options(&states, &specs, &opts_half)
.expect("scaled")
.expect("some");
let mut opts_raw = BlockwiseFitOptions::default();
opts_raw.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::with_uniform_weight(
even_mask, m, 0, 1.0,
)));
let raw = family
.sigma_exact_joint_psi_terms_with_options(&states, &specs, &opts_raw)
.expect("raw")
.expect("some");
let factor = n as f64 / m as f64;
let exp_obj = factor * raw.objective_psi;
let obj_rel = ((scaled.objective_psi - exp_obj) / exp_obj.abs().max(1.0)).abs();
assert!(obj_rel < 1e-12, "objective_psi rel {}", obj_rel);
let exp_score = &raw.score_psi * factor;
let score_rel = rel_diff_array1(&scaled.score_psi, &exp_score);
assert!(score_rel < 1e-12, "score_psi rel {}", score_rel);
let h_scaled = scaled.hessian_psi_operator.as_ref().expect("op").to_dense();
let h_raw = raw.hessian_psi_operator.as_ref().expect("op").to_dense();
let h_exp = &h_raw * factor;
let h_rel = rel_diff_array2(&h_scaled, &h_exp);
assert!(h_rel < 1e-12, "hessian rel {}", h_rel);
}
#[test]
fn bernoulli_sigma_psi_second_order_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_sigma_aware_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let baseline = family
.sigma_exact_joint_psisecond_order_terms(&states)
.expect("baseline")
.expect("some");
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full = family
.sigma_exact_joint_psisecond_order_terms_with_options(&states, &opts_full)
.expect("with full mask")
.expect("some");
let obj_rel = ((with_full.objective_psi_psi - baseline.objective_psi_psi)
/ baseline.objective_psi_psi.abs().max(1.0))
.abs();
assert!(obj_rel < 1e-12, "objective rel {}", obj_rel);
let score_rel = rel_diff_array1(&with_full.score_psi_psi, &baseline.score_psi_psi);
assert!(score_rel < 1e-12, "score rel {}", score_rel);
}
#[test]
fn bernoulli_sigma_psi_second_order_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_sigma_aware_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.sigma_exact_joint_psisecond_order_terms_with_options(&states, &opts_half)
.expect("scaled")
.expect("some");
let mut opts_raw = BlockwiseFitOptions::default();
opts_raw.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::with_uniform_weight(
even_mask, m, 0, 1.0,
)));
let raw = family
.sigma_exact_joint_psisecond_order_terms_with_options(&states, &opts_raw)
.expect("raw")
.expect("some");
let factor = n as f64 / m as f64;
let exp_obj = factor * raw.objective_psi_psi;
let obj_rel = ((scaled.objective_psi_psi - exp_obj) / exp_obj.abs().max(1.0)).abs();
assert!(obj_rel < 1e-12, "objective rel {}", obj_rel);
let exp_score = &raw.score_psi_psi * factor;
let score_rel = rel_diff_array1(&scaled.score_psi_psi, &exp_score);
assert!(score_rel < 1e-12, "score rel {}", score_rel);
}
#[test]
fn bernoulli_sigma_psihessian_directional_derivative_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_sigma_aware_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let dir = array![0.1, -0.2];
let baseline = family
.sigma_exact_joint_psihessian_directional_derivative(&states, &dir)
.expect("baseline")
.expect("some");
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full = family
.sigma_exact_joint_psihessian_directional_derivative_with_options(
&states, &dir, &opts_full,
)
.expect("with full")
.expect("some");
let rel = rel_diff_array2(&with_full, &baseline);
assert!(rel < 1e-12, "drift rel {}", rel);
}
#[test]
fn bernoulli_sigma_psihessian_directional_derivative_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_sigma_aware_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let dir = array![0.1, -0.2];
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.sigma_exact_joint_psihessian_directional_derivative_with_options(
&states, &dir, &opts_half,
)
.expect("scaled")
.expect("some");
let mut opts_raw = BlockwiseFitOptions::default();
opts_raw.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::with_uniform_weight(
even_mask, m, 0, 1.0,
)));
let raw = family
.sigma_exact_joint_psihessian_directional_derivative_with_options(
&states, &dir, &opts_raw,
)
.expect("raw")
.expect("some");
let factor = n as f64 / m as f64;
let exp = &raw * factor;
let rel = rel_diff_array2(&scaled, &exp);
assert!(rel < 1e-12, "drift rel {}", rel);
}
#[test]
fn bernoulli_psi_workspace_with_options_threads_subsample_to_first_order() {
use crate::custom_family::CustomFamilyBlockPsiDerivative;
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_sigma_aware_test_family(n);
let states = rigid_block_states(&family, 0.3, 0.4);
let specs = vec![dummy_blockspec(1, n), dummy_blockspec(1, n)];
let derivative_blocks: Vec<Vec<CustomFamilyBlockPsiDerivative>> = vec![
Vec::new(),
vec![CustomFamilyBlockPsiDerivative::new(
None,
Array2::zeros((0, 0)),
Array2::zeros((0, 0)),
None,
None,
None,
None,
)],
];
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xBEEF_CAFE,
)));
let direct = family
.sigma_exact_joint_psi_terms_with_options(&states, &specs, &opts_half)
.expect("direct sigma terms with options")
.expect("direct some");
let ws = family
.exact_newton_joint_psi_workspace_with_options(
&states,
&specs,
&derivative_blocks,
&opts_half,
)
.expect("workspace with options")
.expect("workspace some");
let psi_total: usize = derivative_blocks.iter().map(Vec::len).sum();
let sigma_psi = psi_total - 1;
let via_ws = ws
.first_order_terms(sigma_psi)
.expect("ws first_order_terms")
.expect("some");
assert_eq!(via_ws.objective_psi, direct.objective_psi);
let score_rel = rel_diff_array1(&via_ws.score_psi, &direct.score_psi);
assert!(score_rel == 0.0, "score_psi diverged: rel {}", score_rel);
let h_ws = via_ws
.hessian_psi_operator
.as_ref()
.expect("ws hessian op")
.to_dense();
let h_direct = direct
.hessian_psi_operator
.as_ref()
.expect("direct hessian op")
.to_dense();
let h_rel = rel_diff_array2(&h_ws, &h_direct);
assert!(h_rel == 0.0, "hessian diverged: rel {}", h_rel);
let ws_full = family
.exact_newton_joint_psi_workspace_with_options(
&states,
&specs,
&derivative_blocks,
&BlockwiseFitOptions::default(),
)
.expect("workspace full")
.expect("some");
let via_ws_full = ws_full
.first_order_terms(sigma_psi)
.expect("ws_full first_order_terms")
.expect("some");
assert!(
(via_ws.objective_psi - via_ws_full.objective_psi).abs() > 1e-9,
"subsample objective {} too close to full-data {} (subsample not threaded?)",
via_ws.objective_psi,
via_ws_full.objective_psi
);
let m_f = m as f64;
let n_f = n as f64;
let ratio_bound = (n_f / m_f) * 2.0 + 1.0;
let ratio = (via_ws.objective_psi.abs() + 1.0) / (via_ws_full.objective_psi.abs() + 1.0);
assert!(
ratio < ratio_bound && (1.0 / ratio) < ratio_bound,
"subsample/full ratio {} outside coarse bound {}",
ratio,
ratio_bound
);
}
fn build_test_link_deviation_block_from_seed(
seed: &Array1<f64>,
cfg: &DeviationBlockConfig,
) -> Result<DeviationPrepared, String> {
build_link_deviation_block_from_knots_design_seed_and_weights(seed, seed, cfg)
}
#[test]
fn score_warp_basis_smoothness_penalty_is_full_rank() {
let seed = array![-2.0, -1.0, 0.0, 1.0, 2.0];
let cfg = DeviationBlockConfig {
num_internal_knots: 5,
..DeviationBlockConfig::default()
};
let prepared = build_score_warp_deviation_block_from_seed(&seed, &cfg)
.expect("build smoothness-null-space-drop score-warp");
let max_order = cfg
.penalty_orders
.iter()
.copied()
.max()
.or(Some(cfg.penalty_order))
.filter(|o| *o > 0)
.expect("test config has a positive penalty order");
let (penalty, nullity) = prepared
.runtime
.integrated_derivative_penalty_with_nullity(max_order)
.expect("integrated penalty on transformed basis");
assert_eq!(
nullity, 0,
"smoothness-null-space-drop basis must have zero nullity at the max configured \
derivative order; got {nullity}"
);
use crate::faer_ndarray::FaerEigh;
let (evals, _) = penalty
.eigh(faer::Side::Lower)
.expect("eigendecomposition of transformed-basis penalty");
let evals_slice = evals
.as_slice()
.expect("contiguous transformed-basis penalty eigenvalues");
let threshold = crate::estimate::reml::unified::positive_eigenvalue_threshold(evals_slice);
let smallest = evals_slice.iter().copied().fold(f64::INFINITY, f64::min);
assert!(
smallest > threshold,
"smallest eigenvalue {smallest} of transformed-basis penalty must exceed positive \
threshold {threshold} — null space was not fully dropped"
);
}
#[test]
fn link_deviation_basis_smoothness_penalty_is_full_rank() {
let q = array![-2.0, -0.8, -0.1, 0.4, 1.3, 2.1];
let cfg = DeviationBlockConfig {
num_internal_knots: 5,
..DeviationBlockConfig::default()
};
let prepared = build_link_deviation_block_from_knots_design_seed_and_weights(&q, &q, &cfg)
.expect("build smoothness-null-space-drop link-deviation");
let max_order = cfg
.penalty_orders
.iter()
.copied()
.max()
.or(Some(cfg.penalty_order))
.filter(|o| *o > 0)
.expect("test config has a positive penalty order");
let (penalty, nullity) = prepared
.runtime
.integrated_derivative_penalty_with_nullity(max_order)
.expect("integrated penalty on transformed basis");
assert_eq!(
nullity, 0,
"smoothness-null-space-drop basis must have zero nullity at the max configured \
derivative order; got {nullity}"
);
use crate::faer_ndarray::FaerEigh;
let (evals, _) = penalty
.eigh(faer::Side::Lower)
.expect("eigendecomposition of transformed-basis penalty");
let evals_slice = evals
.as_slice()
.expect("contiguous transformed-basis penalty eigenvalues");
let threshold = crate::estimate::reml::unified::positive_eigenvalue_threshold(evals_slice);
let smallest = evals_slice.iter().copied().fold(f64::INFINITY, f64::min);
assert!(
smallest > threshold,
"smallest eigenvalue {smallest} of transformed-basis penalty must exceed positive \
threshold {threshold} — null space was not fully dropped"
);
}
#[test]
fn bernoulli_marginal_slope_rejects_nonprobit_base_link() {
let y = array![0.0, 1.0];
let weights = array![1.0, 1.0];
let z = array![-0.4, 0.9];
let design =
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(array![[1.0], [1.0]]));
let spec = BernoulliMarginalSlopeTermSpec {
y,
weights,
z,
base_link: InverseLink::Standard(StandardLink::Logit),
marginalspec: empty_termspec(),
logslopespec: empty_termspec(),
marginal_offset: Array1::zeros(2),
logslope_offset: Array1::zeros(2),
frailty: FrailtySpec::None,
score_warp: None,
link_dev: None,
latent_z_policy: LatentZPolicy::default(),
};
let err = validate_spec(design.to_dense().view(), &spec)
.expect_err("non-probit marginal-slope link should be rejected");
assert!(err.contains("requires link(type=probit)"));
let err = bernoulli_marginal_slope_eta_from_probability(
&InverseLink::Standard(StandardLink::Logit),
0.5,
"test logit inverse",
)
.expect_err("non-probit marginal-slope inverse should be rejected");
assert!(err.contains("requires link(type=probit)"));
}
fn expand_integer_weight_rows(
y: &Array1<f64>,
z: &Array1<f64>,
weights: &Array1<f64>,
) -> (Array1<f64>, Array1<f64>) {
let mut y_expanded = Vec::new();
let mut z_expanded = Vec::new();
for i in 0..y.len() {
let reps = weights[i] as usize;
assert!(
(weights[i] - reps as f64).abs() < 1e-12,
"test helper expects integer weights, got {}",
weights[i]
);
for _ in 0..reps {
y_expanded.push(y[i]);
z_expanded.push(z[i]);
}
}
(Array1::from_vec(y_expanded), Array1::from_vec(z_expanded))
}
#[test]
fn link_dev_without_score_warp_exposes_structural_derivative_lower_bounds() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link block initial beta")
.len();
let beta_link = Array1::from_iter((0..link_dim).map(|idx| 0.1 * (idx as f64 + 1.0)));
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(Array1::zeros(seed.len())),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(beta_link.clone(), seed.len()),
];
let slices = block_slices(&family);
assert!(slices.h.is_none(), "score-warp slice should be absent");
let link_slice = slices.w.as_ref().expect("link slice");
assert_eq!(
slices.marginal.len(),
0,
"zero-column marginal design should not contribute coefficient coordinates"
);
assert_eq!(
slices.logslope.len(),
0,
"zero-column logslope design should not contribute coefficient coordinates"
);
assert_eq!(
link_slice.start, 0,
"link-only coefficients should start at 0"
);
assert_eq!(link_slice.len(), link_dim);
let primary = primary_slices(&slices);
assert!(primary.h.is_none(), "primary h slice should be absent");
let primary_w = primary.w.as_ref().expect("primary link slice");
assert_eq!(primary_w.start, 2, "primary link slice should start at 2");
assert_eq!(primary.total, 2 + link_dim);
family
.build_exact_eval_cache(&block_states)
.expect("eval cache");
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(0, &block_states, None, true)
.expect("row context");
let (nll, grad, hess) = family
.compute_row_primary_gradient_hessian(0, &block_states, &primary, &row_ctx)
.expect("analytic flex eval");
assert!(nll.is_finite(), "neglog should be finite for link-dev-only");
assert!(
grad.iter().all(|v| v.is_finite()),
"gradient should be finite"
);
assert!(
hess.iter().all(|v| v.is_finite()),
"Hessian should be finite"
);
let dummy_spec = dummy_blockspec(link_dim, seed.len());
assert!(
family
.block_linear_constraints(&block_states, 1, &dummy_spec)
.expect("non-link constraint lookup")
.is_none(),
"non-link block should not expose auxiliary monotonicity constraints"
);
let constraints = family
.block_linear_constraints(&block_states, 2, &dummy_spec)
.expect("link constraint lookup")
.expect("link constraints");
assert_eq!(constraints.a.ncols(), link_dim);
assert_eq!(constraints.b.len(), constraints.a.nrows());
assert!(
constraints.a.nrows() >= link_dim,
"anchored link constraints should be expressed in raw derivative-control rows"
);
assert_eq!(
constraints.b,
Array1::<f64>::from_elem(
constraints.a.nrows(),
prepared.runtime.monotonicity_eps() - 1.0
)
);
}
#[test]
fn zero_deviation_intercept_fast_path_matches_denested_calibration() {
let seed = Array1::from_iter((0..25).map(|i| -2.4 + 4.8 * i as f64 / 24.0));
let score_prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 8,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let link_prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 8,
..DeviationBlockConfig::default()
},
)
.expect("build link-deviation block");
let n = seed.len();
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(Array1::zeros(n)),
weights: Arc::new(Array1::ones(n)),
z: Arc::new(seed.clone()),
gaussian_frailty_sd: Some(0.65),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((n, 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((n, 0)),
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
intercept_warm_starts: Some(new_intercept_warm_start_cache(n)),
..default_test_family()
};
let marginal_eta = 0.35;
let slope = -0.8;
let beta_h = Array1::zeros(score_prepared.runtime.basis_dim());
let beta_w = Array1::zeros(link_prepared.runtime.basis_dim());
let marginal = family
.marginal_link_map(marginal_eta)
.expect("marginal map");
let scale = family.probit_frailty_scale();
let rigid_a = rigid_prescale_intercept_from_marginal(marginal.q, slope, scale);
let (f_rigid, f_a_rigid, _) = family
.evaluate_denested_calibration_newton(
rigid_a,
marginal_eta,
slope,
Some(&beta_h),
Some(&beta_w),
)
.expect("denested zero-deviation calibration");
assert!(
f_rigid.abs() <= 5e-13,
"closed-form rigid intercept residual should be at machine epsilon, got {f_rigid}"
);
let analytic_deriv = rigid_prescale_intercept_derivative_abs(marginal.q, slope, scale);
assert!(
(f_a_rigid - analytic_deriv).abs() <= 5e-13,
"denested derivative {f_a_rigid} != analytic {analytic_deriv}"
);
let (a_fast, deriv_fast, fast_path) = family
.solve_row_intercept_base(0, marginal_eta, slope, Some(&beta_h), Some(&beta_w), None)
.expect("zero-deviation solve");
assert!(
fast_path,
"zero coefficients should take analytic fast path"
);
assert_eq!(a_fast, rigid_a);
assert_eq!(deriv_fast, analytic_deriv);
}
#[test]
fn exact_layout_ignores_dummy_beta_widths_for_empty_design_blocks() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let score_prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let link_prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(Array1::zeros(seed.len())),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::zeros(score_prepared.runtime.basis_dim()),
seed.len(),
),
dummy_block_state(Array1::zeros(link_prepared.runtime.basis_dim()), seed.len()),
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
assert_eq!(cache.slices.marginal.len(), 0);
assert_eq!(cache.slices.logslope.len(), 0);
assert_eq!(cache.slices.h.as_ref().expect("h slice").start, 0);
assert_eq!(
cache.slices.w.as_ref().expect("w slice").start,
score_prepared.runtime.basis_dim()
);
assert_eq!(
cache.slices.total,
score_prepared.runtime.basis_dim() + link_prepared.runtime.basis_dim()
);
assert_eq!(cache.primary.q, 0);
assert_eq!(cache.primary.logslope, 1);
assert_eq!(cache.primary.h.as_ref().expect("primary h").start, 2);
assert_eq!(
cache.primary.w.as_ref().expect("primary w").start,
2 + score_prepared.runtime.basis_dim()
);
}
#[test]
fn score_warp_block_exposes_structural_derivative_lower_bounds() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(Array1::zeros(seed.len())),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(Array1::zeros(score_dim), seed.len()),
];
let dummy_spec = dummy_blockspec(score_dim, seed.len());
let constraints = family
.block_linear_constraints(&block_states, 2, &dummy_spec)
.expect("constraint lookup")
.expect("score-warp constraints");
assert_eq!(constraints.a.ncols(), score_dim);
assert_eq!(constraints.b.len(), constraints.a.nrows());
assert!(
constraints.a.nrows() >= score_dim,
"anchored score-warp constraints should be expressed in raw derivative-control rows"
);
assert_eq!(
constraints.b,
Array1::<f64>::from_elem(
constraints.a.nrows(),
prepared.runtime.monotonicity_eps() - 1.0
)
);
}
#[test]
fn post_update_block_beta_rejects_infeasible_score_warp_step() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared.block.design.ncols();
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(Array1::zeros(seed.len())),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let current = Array1::<f64>::zeros(score_dim);
let mut proposed = current.clone();
proposed[0] = -128.0;
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(current.clone(), seed.len()),
];
let spec = dummy_blockspec(score_dim, seed.len());
let err = family
.post_update_block_beta(&block_states, 2, &spec, proposed.clone())
.expect_err("post-update must not repair infeasible constrained beta");
assert!(
err.contains("structural monotonicity") || err.contains("exact monotonicity"),
"unexpected error: {err}"
);
}
#[test]
fn structural_deviation_runtime_is_piecewise_cubic() {
let seed = array![-1.0, 0.0, 1.0];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
..DeviationBlockConfig::default()
},
)
.expect("structural deviation basis");
assert_eq!(prepared.runtime.degree(), 3);
assert_eq!(prepared.runtime.value_span_degree(), 3);
let has_cubic_curvature = prepared
.runtime
.span_c3()
.iter()
.any(|value| value.abs() > 1e-12);
assert!(
has_cubic_curvature,
"structural deviation basis must expose true cubic span coefficients"
);
}
#[test]
fn structural_deviation_runtime_is_c2_at_internal_breakpoints() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("structural deviation basis");
let dim = prepared.block.design.ncols();
let beta = Array1::from_iter((0..dim).map(|idx| 0.015 * (idx as f64 + 1.0)));
let n_spans = prepared.runtime.breakpoints().len().saturating_sub(1);
for span_idx in 1..n_spans {
let left_cubic = prepared
.runtime
.local_cubic_on_span(&beta, span_idx - 1)
.expect("left span cubic");
let right_cubic = prepared
.runtime
.local_cubic_on_span(&beta, span_idx)
.expect("right span cubic");
let knot = prepared.runtime.breakpoints()[span_idx];
assert!(
(left_cubic.evaluate(knot) - right_cubic.evaluate(knot)).abs() <= 1e-10,
"deviation value should be continuous at breakpoint {span_idx}"
);
assert!(
(left_cubic.first_derivative(knot) - right_cubic.first_derivative(knot)).abs()
<= 1e-10,
"deviation first derivative should be continuous at breakpoint {span_idx}"
);
assert!(
(left_cubic.second_derivative(knot) - right_cubic.second_derivative(knot)).abs()
<= 1e-10,
"deviation second derivative should be continuous at breakpoint {span_idx}"
);
}
}
#[test]
fn structural_deviation_rejects_noncubic_degree() {
let seed = array![-1.0, 0.0, 1.0];
let err = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
degree: 4,
..DeviationBlockConfig::default()
},
)
.expect_err("structural deviation block should reject non-cubic degree");
assert!(err.contains("degree must be 3"));
}
#[test]
fn deviation_runtime_replays_exact_training_design() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build deviation block");
let replayed = prepared.runtime.design(&seed).expect("replayed design");
let trained = prepared.block.design.to_dense();
assert_eq!(replayed.dim(), trained.dim());
for i in 0..replayed.nrows() {
for j in 0..replayed.ncols() {
assert!(
(replayed[[i, j]] - trained[[i, j]]).abs() <= 1e-10,
"training-basis replay mismatch at ({i},{j})"
);
}
}
}
#[test]
fn structural_constraints_match_exact_monotonicity_guard() {
let seed = array![-1.0, 0.0, 1.0, 2.0];
let prepared =
build_score_warp_deviation_block_from_seed(&seed, &DeviationBlockConfig::default())
.expect("build deviation block");
let constraints = prepared.runtime.structural_monotonicity_constraints();
let dim = constraints.a.ncols();
assert_eq!(dim, prepared.runtime.basis_dim());
assert_eq!(
constraints.a.nrows(),
3 * prepared.runtime.breakpoints().len().saturating_sub(1)
);
assert_eq!(
constraints.b,
Array1::<f64>::from_elem(
constraints.a.nrows(),
prepared.runtime.monotonicity_eps() - 1.0
)
);
let feasible = Array1::<f64>::zeros(dim);
prepared
.runtime
.monotonicity_feasible(&feasible, "feasible structural beta")
.expect("zero deviation should be feasible");
let d1_design = prepared
.runtime
.first_derivative_design(&seed)
.expect("derivative design");
let row_idx = (0..d1_design.nrows())
.find(|&idx| d1_design.row(idx).dot(&d1_design.row(idx)) > 0.0)
.expect("derivative design should include a nonzero row");
let derivative_row = d1_design.row(row_idx).to_owned();
let row_norm_sq = derivative_row.dot(&derivative_row);
let infeasible = derivative_row.mapv(|value| -2.0 * value / row_norm_sq);
assert!(
prepared
.runtime
.monotonicity_feasible(&infeasible, "infeasible structural beta")
.is_err()
);
}
#[test]
fn structural_constraints_are_quadratic_derivative_bernstein_controls() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build deviation block");
let constraints = prepared.runtime.structural_monotonicity_constraints();
let beta = Array1::from_iter((0..prepared.runtime.basis_dim()).map(|idx| {
let centered = idx as f64 - 0.5 * (prepared.runtime.basis_dim() as f64 - 1.0);
0.025 * centered
}));
let controls = constraints.a.dot(&beta);
let n_spans = prepared.runtime.breakpoints().len().saturating_sub(1);
for span_idx in 0..n_spans {
let cubic = prepared
.runtime
.local_cubic_on_span(&beta, span_idx)
.expect("local cubic");
let left = cubic.left;
let right = cubic.right;
let mid = 0.5 * (left + right);
let b0 = controls[3 * span_idx];
let b1 = controls[3 * span_idx + 1];
let b2 = controls[3 * span_idx + 2];
assert!(
(b0 - cubic.first_derivative(left)).abs() <= 1e-10,
"left Bernstein control should equal derivative at span start"
);
assert!(
(b2 - cubic.first_derivative(right)).abs() <= 1e-10,
"right Bernstein control should equal derivative at span end"
);
let midpoint_from_bernstein = 0.25 * b0 + 0.5 * b1 + 0.25 * b2;
assert!(
(midpoint_from_bernstein - cubic.first_derivative(mid)).abs() <= 1e-10,
"quadratic Bernstein controls should reconstruct derivative at span midpoint"
);
}
}
#[test]
fn deviation_penalties_are_integrated_function_penalties() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
penalty_order: 2,
penalty_orders: vec![1, 2, 3],
double_penalty: true,
..DeviationBlockConfig::default()
},
)
.expect("build deviation block");
let expected_orders = [1, 2, 3, 0];
assert_eq!(prepared.block.penalties.len(), expected_orders.len());
for ((penalty, &nullity), &order) in prepared
.block
.penalties
.iter()
.zip(prepared.block.nullspace_dims.iter())
.zip(expected_orders.iter())
{
let crate::solver::estimate::PenaltySpec::Dense(actual) = penalty else {
panic!("deviation penalties should be dense local Gram matrices");
};
let (expected, expected_nullity) = prepared
.runtime
.integrated_derivative_penalty_with_nullity(order)
.expect("integrated function penalty");
assert_eq!(nullity, expected_nullity);
assert_eq!(actual.dim(), expected.dim());
for i in 0..actual.nrows() {
for j in 0..actual.ncols() {
assert!(
(actual[[i, j]] - expected[[i, j]]).abs() <= 1e-10,
"penalty order {order} mismatch at ({i},{j}): got {}, expected {}",
actual[[i, j]],
expected[[i, j]]
);
}
}
}
let crate::solver::estimate::PenaltySpec::Dense(l2_penalty) = &prepared.block.penalties[1]
else {
panic!("deviation double penalty should be dense");
};
let mut max_identity_diff = 0.0_f64;
for i in 0..l2_penalty.nrows() {
for j in 0..l2_penalty.ncols() {
let identity = if i == j { 1.0 } else { 0.0 };
max_identity_diff = max_identity_diff.max((l2_penalty[[i, j]] - identity).abs());
}
}
assert!(
max_identity_diff > 1e-6,
"deviation double penalty must be integrated L2, not coefficient identity"
);
}
#[test]
fn local_cubic_span_reconstructs_deviation_exactly() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build deviation block");
let dim = prepared.block.design.ncols();
let beta = Array1::from_iter((0..dim).map(|idx| 0.025 * (idx as f64 + 1.0)));
let n_spans = prepared.runtime.breakpoints().len().saturating_sub(1);
let support_left = prepared.runtime.breakpoints()[0];
let support_right =
prepared.runtime.breakpoints()[prepared.runtime.breakpoints().len() - 1];
for span_idx in 0..n_spans {
let cubic = prepared
.runtime
.local_cubic_on_span(&beta, span_idx)
.expect("local cubic coefficients");
let left = cubic.left;
let right = cubic.right;
let x_eval = array![left, 0.5 * (left + right), right];
let value_design = prepared.runtime.design(&x_eval).expect("value design");
let d1_design = prepared
.runtime
.first_derivative_design(&x_eval)
.expect("first derivative design");
let d2_design = prepared
.runtime
.second_derivative_design(&x_eval)
.expect("second derivative design");
let expected = value_design.dot(&beta);
let expected_d1 = d1_design.dot(&beta);
let expected_d2 = d2_design.dot(&beta);
for i in 0..x_eval.len() {
let x = x_eval[i];
assert!(
(cubic.evaluate(x) - expected[i]).abs() < 1e-10,
"span {span_idx}, x={x:.6}: cubic value mismatch"
);
assert!(
(cubic.first_derivative(x) - expected_d1[i]).abs() < 1e-10,
"span {span_idx}, x={x:.6}: cubic first-derivative mismatch"
);
assert!(
(cubic.second_derivative(x) - expected_d2[i]).abs() < 1e-10,
"span {span_idx}, x={x:.6}: cubic second-derivative mismatch"
);
let selected = prepared
.runtime
.local_cubic_at(&beta, x)
.expect("lookup cubic at x");
if x < support_left || x > support_right {
assert!(selected.c1.abs() < 1e-12);
assert!(selected.c2.abs() < 1e-12);
assert!(selected.c3.abs() < 1e-12);
assert!((selected.evaluate(x) - expected[i]).abs() < 1e-10);
} else {
let expected_span_idx = if i == 0 && span_idx > 0 {
span_idx - 1
} else {
span_idx
};
let expected_cubic = prepared
.runtime
.local_cubic_on_span(&beta, expected_span_idx)
.expect("expected lookup cubic on span");
assert_eq!(selected, expected_cubic);
}
}
}
}
#[test]
fn basis_span_cubic_reconstructs_basis_column_exactly() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build deviation block");
let basis_idx = 0usize;
let support_left = prepared.runtime.breakpoints()[0];
let support_right =
prepared.runtime.breakpoints()[prepared.runtime.breakpoints().len() - 1];
let cubic = prepared
.runtime
.basis_span_cubic(0, basis_idx)
.expect("basis span cubic");
let x_eval = array![cubic.left, 0.5 * (cubic.left + cubic.right), cubic.right];
let design = prepared.runtime.design(&x_eval).expect("basis design");
let d1 = prepared
.runtime
.first_derivative_design(&x_eval)
.expect("basis d1 design");
for i in 0..x_eval.len() {
let x = x_eval[i];
assert!((cubic.evaluate(x) - design[[i, basis_idx]]).abs() < 1e-10);
assert!((cubic.first_derivative(x) - d1[[i, basis_idx]]).abs() < 1e-10);
let selected = prepared
.runtime
.basis_cubic_at(basis_idx, x)
.expect("basis cubic at x");
if x < support_left || x > support_right {
assert!(selected.c1.abs() < 1e-12);
assert!(selected.c2.abs() < 1e-12);
assert!(selected.c3.abs() < 1e-12);
assert!((selected.evaluate(x) - design[[i, basis_idx]]).abs() < 1e-10);
} else {
let expected_span_idx = 0;
let expected_cubic = prepared
.runtime
.basis_span_cubic(expected_span_idx, basis_idx)
.expect("expected basis span cubic");
assert_eq!(selected, expected_cubic);
}
}
}
#[test]
fn deviation_runtime_saturates_outside_support() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build deviation block");
let dim = prepared.block.design.ncols();
let beta = Array1::from_iter((0..dim).map(|idx| 0.02 * (idx as f64 + 1.0)));
let left = prepared.runtime.breakpoints()[0];
let right = prepared.runtime.breakpoints()[prepared.runtime.breakpoints().len() - 1];
let left_tail_near = prepared
.runtime
.local_cubic_at(&beta, left - 0.25)
.expect("left tail");
let left_tail_far = prepared
.runtime
.local_cubic_at(&beta, left - 3.0)
.expect("left far tail");
let right_tail_near = prepared
.runtime
.local_cubic_at(&beta, right + 0.25)
.expect("right tail");
let right_tail_far = prepared
.runtime
.local_cubic_at(&beta, right + 3.0)
.expect("right far tail");
for cubic in [
left_tail_near,
left_tail_far,
right_tail_near,
right_tail_far,
] {
assert!(cubic.c1.abs() < 1e-12);
assert!(cubic.c2.abs() < 1e-12);
assert!(cubic.c3.abs() < 1e-12);
}
assert!((left_tail_near.c0 - left_tail_far.c0).abs() < 1e-12);
assert!((right_tail_near.c0 - right_tail_far.c0).abs() < 1e-12);
}
#[test]
fn deviation_runtime_replays_the_exact_training_basis() {
let seed = array![-2.0, -1.0, -0.25, 0.25, 1.0, 2.0];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build deviation block");
let replayed = prepared
.runtime
.design(&seed)
.expect("replay anchored deviation design");
let trained = prepared.block.design.to_dense();
assert_eq!(replayed.dim(), trained.dim());
for i in 0..replayed.nrows() {
for j in 0..replayed.ncols() {
assert!(
(replayed[[i, j]] - trained[[i, j]]).abs() <= 1e-10,
"replayed anchored deviation design mismatch at ({i},{j})"
);
}
}
}
#[test]
fn denested_microcells_follow_score_and_link_breaks() {
let score_seed = array![-2.0, -1.0, 0.0, 1.0, 2.0];
let link_seed = array![-1.5, -0.5, 0.5, 1.5];
let score_prepared = build_score_warp_deviation_block_from_seed(
&score_seed,
&DeviationBlockConfig {
num_internal_knots: 3,
..DeviationBlockConfig::default()
},
)
.expect("build score warp block");
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 3,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let beta_h = Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.02 * (idx as f64 + 1.0)),
);
let beta_w = Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
);
let exact_cells_a0 = build_denested_partition_cells(
0.25,
0.9,
score_prepared
.runtime
.breakpoints()
.as_slice()
.expect("score breaks"),
link_prepared
.runtime
.breakpoints()
.as_slice()
.expect("link breaks"),
|z| score_prepared.runtime.local_cubic_at(&beta_h, z),
|u| link_prepared.runtime.local_cubic_at(&beta_w, u),
)
.expect("exact module microcells for a=0.25");
let exact_cells_a1 = build_denested_partition_cells(
0.55,
0.9,
score_prepared
.runtime
.breakpoints()
.as_slice()
.expect("score breaks"),
link_prepared
.runtime
.breakpoints()
.as_slice()
.expect("link breaks"),
|z| score_prepared.runtime.local_cubic_at(&beta_h, z),
|u| link_prepared.runtime.local_cubic_at(&beta_w, u),
)
.expect("exact module microcells for a=0.55");
assert!(
exact_cells_a0.len() >= score_prepared.runtime.breakpoints().len().saturating_sub(1),
"microcell partition should refine the score spans"
);
assert!(
exact_cells_a0
.windows(2)
.all(|w| (w[0].cell.right - w[1].cell.left).abs() <= 1e-12),
"microcells should tile the partition contiguously"
);
assert!(exact_cells_a0.first().unwrap().cell.left.is_infinite());
assert!(exact_cells_a0.last().unwrap().cell.right.is_infinite());
assert!(
exact_cells_a0
.iter()
.zip(exact_cells_a1.iter())
.any(|(lhs, rhs)| (lhs.cell.left - rhs.cell.left).abs() > 1e-10),
"changing the intercept should move at least one link-induced breakpoint"
);
}
#[test]
fn denested_microcell_eta_matches_direct_denested_formula() {
let score_seed = array![-2.0, -1.0, 0.0, 1.0, 2.0];
let link_seed = array![-1.5, -0.5, 0.5, 1.5];
let score_prepared = build_score_warp_deviation_block_from_seed(
&score_seed,
&DeviationBlockConfig {
num_internal_knots: 3,
..DeviationBlockConfig::default()
},
)
.expect("build score warp block");
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 3,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let beta_h = Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
);
let beta_w = Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.02 * (idx as f64 + 1.0)),
);
let a = 0.35;
let b = -0.7;
let cells = build_denested_partition_cells(
a,
b,
score_prepared
.runtime
.breakpoints()
.as_slice()
.expect("score breaks"),
link_prepared
.runtime
.breakpoints()
.as_slice()
.expect("link breaks"),
|z| score_prepared.runtime.local_cubic_at(&beta_h, z),
|u| link_prepared.runtime.local_cubic_at(&beta_w, u),
)
.expect("microcells");
for cell in &cells {
let z = exact_kernel::interval_probe_point(cell.cell.left, cell.cell.right)
.expect("finite microcell probe");
let h = score_prepared
.runtime
.design(&array![z])
.expect("score design")
.row(0)
.dot(&beta_h);
let link = link_prepared
.runtime
.design(&array![a + b * z])
.expect("link design")
.row(0)
.dot(&beta_w);
let expected = a + b * z + b * h + link;
assert!(
(cell.cell.eta(z) - expected).abs() < 1e-10,
"microcell eta should equal the direct de-nested predictor at z={z:.6}"
);
}
}
#[test]
fn local_cubic_global_transform_reconstructs_same_function() {
let cubic = exact_kernel::LocalSpanCubic {
left: -1.3,
right: 0.7,
c0: 0.4,
c1: -0.2,
c2: 0.15,
c3: -0.05,
};
let (g0, g1, g2, g3) = exact_global_cubic_from_local(LocalSpanCubic {
left: cubic.left,
right: cubic.right,
c0: cubic.c0,
c1: cubic.c1,
c2: cubic.c2,
c3: cubic.c3,
});
for &x in &[-1.3, -0.8, -0.1, 0.5, 0.7] {
let direct = cubic.evaluate(x);
let global = g0 + g1 * x + g2 * x * x + g3 * x * x * x;
assert!(
(direct - global).abs() < 1e-12,
"global cubic transform should preserve the span polynomial at x={x}"
);
}
}
#[test]
fn denested_branch_selection_uses_normalized_cell_coefficients() {
let affine = ExactDenestedCubicCell {
left: -1.0,
right: 1.0,
c0: 0.2,
c1: -0.4,
c2: 1e-13,
c3: -1e-13,
};
let quartic = ExactDenestedCubicCell {
c2: 2e-4,
c3: 1e-13,
..affine
};
let sextic = ExactDenestedCubicCell {
c2: 2e-4,
c3: 5e-3,
..affine
};
assert_eq!(
branch_exact_cell(affine).expect("affine branch"),
ExactCellBranchShared::Affine
);
assert_eq!(
branch_exact_cell(quartic).expect("quartic branch"),
ExactCellBranchShared::Quartic
);
assert_eq!(
branch_exact_cell(sextic).expect("sextic branch"),
ExactCellBranchShared::Sextic
);
}
#[test]
fn denested_cell_coefficient_partials_match_finite_differences() {
let score_span = exact_kernel::LocalSpanCubic {
left: -0.75,
right: 0.25,
c0: 0.08,
c1: -0.03,
c2: 0.02,
c3: -0.01,
};
let link_span = exact_kernel::LocalSpanCubic {
left: -0.6,
right: 0.9,
c0: -0.05,
c1: 0.04,
c2: -0.02,
c3: 0.015,
};
let a = 0.3;
let b = -0.7;
let eps = 1e-6;
let coeffs = |aa: f64, bb: f64| {
let (h0, h1, h2, h3) = exact_global_cubic_from_local(LocalSpanCubic {
left: score_span.left,
right: score_span.right,
c0: score_span.c0,
c1: score_span.c1,
c2: score_span.c2,
c3: score_span.c3,
});
let (d0, d1, d2, d3) = exact_transformed_link_cubic(
LocalSpanCubic {
left: link_span.left,
right: link_span.right,
c0: link_span.c0,
c1: link_span.c1,
c2: link_span.c2,
c3: link_span.c3,
},
aa,
bb,
);
[
aa + bb * h0 + d0,
bb + bb * h1 + d1,
bb * h2 + d2,
bb * h3 + d3,
]
};
let (dc_da, dc_db) = exact_denested_cell_coefficient_partials(
LocalSpanCubic {
left: score_span.left,
right: score_span.right,
c0: score_span.c0,
c1: score_span.c1,
c2: score_span.c2,
c3: score_span.c3,
},
LocalSpanCubic {
left: link_span.left,
right: link_span.right,
c0: link_span.c0,
c1: link_span.c1,
c2: link_span.c2,
c3: link_span.c3,
},
a,
b,
);
let plus_a = coeffs(a + eps, b);
let minus_a = coeffs(a - eps, b);
let plus_b = coeffs(a, b + eps);
let minus_b = coeffs(a, b - eps);
for j in 0..4 {
let fd_a = (plus_a[j] - minus_a[j]) / (2.0 * eps);
let fd_b = (plus_b[j] - minus_b[j]) / (2.0 * eps);
assert!(
(dc_da[j] - fd_a).abs() < 1e-6,
"dc/da mismatch at coefficient {j}: analytic={}, fd={fd_a}",
dc_da[j]
);
assert!(
(dc_db[j] - fd_b).abs() < 1e-6,
"dc/db mismatch at coefficient {j}: analytic={}, fd={fd_b}",
dc_db[j]
);
}
}
#[test]
fn observed_denested_partials_include_third_a_derivative_for_piecewise_cubic_link() {
let z = array![-0.8, 0.2, 1.1];
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let beta_w = Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
);
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 1.0]),
weights: Arc::new(array![1.0, 0.7, 1.3]),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let a = 0.35;
let b = 0.6;
let row = 1usize;
let obs = family
.observed_denested_cell_partials(row, a, b, None, Some(&beta_w))
.expect("observed denested partials");
let u_obs = a + b * z[row];
let link_span = link_prepared
.runtime
.local_cubic_at(&beta_w, u_obs)
.expect("local cubic at observed point");
let expected_daaa = exact_kernel::denested_cell_third_partials(link_span).0;
assert_eq!(obs.dc_daaa, expected_daaa);
assert!(
eval_coeff4_at(&obs.dc_daaa, z[row]).abs() > 1e-12,
"piecewise-cubic link spans should contribute a third a-derivative"
);
}
#[test]
fn pooled_probit_baseline_matches_expanded_integer_weight_fit() {
let y = array![0.0, 1.0, 0.0, 1.0];
let z = array![-1.5, -0.2, 0.4, 1.4];
let weights = array![25.0, 2.0, 1.0, 20.0];
let weighted = pooled_probit_baseline(&y, &z, &weights).expect("weighted baseline");
let unweighted =
pooled_probit_baseline(&y, &z, &Array1::ones(y.len())).expect("unweighted baseline");
let (y_expanded, z_expanded) = expand_integer_weight_rows(&y, &z, &weights);
let expanded =
pooled_probit_baseline(&y_expanded, &z_expanded, &Array1::ones(y_expanded.len()))
.expect("expanded baseline");
assert!(
pair_distance(expanded, unweighted) > 1e-2,
"test data should distinguish weighted from unweighted seeding"
);
assert!(
pair_distance(weighted, expanded) < 1e-8,
"weighted pilot baseline should match the expanded integer-weight fit"
);
}
#[test]
fn validate_spec_rejects_nonfinite_or_negative_weights() {
let data = Array2::<f64>::zeros((3, 0));
let y = array![0.0, 1.0, 0.0];
let z = array![-1.0, 0.0, 1.0];
let err = validate_spec(
data.view(),
&base_spec(y.clone(), array![1.0, f64::NAN, 1.0], z.clone()),
)
.expect_err("non-finite weights should be rejected");
assert!(err.contains("finite non-negative weights"));
let err = validate_spec(data.view(), &base_spec(y, array![1.0, -0.5, 1.0], z))
.expect_err("negative weights should be rejected");
assert!(err.contains("finite non-negative weights"));
}
#[test]
fn validate_spec_rejects_nonfinite_z_values() {
let data = Array2::<f64>::zeros((3, 0));
let err = validate_spec(
data.view(),
&base_spec(
array![0.0, 1.0, 0.0],
array![1.0, 1.0, 1.0],
array![-1.0, f64::INFINITY, 1.0],
),
)
.expect_err("non-finite z should be rejected");
assert!(err.contains("finite z values"));
}
#[test]
fn validate_spec_accepts_learnable_gaussian_shift_sigma() {
assert!(file!().ends_with(".rs"));
let data = Array2::<f64>::zeros((3, 0));
let mut spec = base_spec(
array![0.0, 1.0, 0.0],
array![1.0, 1.0, 1.0],
array![-1.0, 0.0, 1.0],
);
spec.frailty = FrailtySpec::GaussianShift { sigma_fixed: None };
validate_spec(data.view(), &spec).expect("learnable GaussianShift sigma should validate");
}
#[test]
fn signed_probit_helpers_handle_nonfinite_boundaries_explicitly() {
let (logcdf_pos, lambda_pos) = signed_probit_logcdf_and_mills_ratio(f64::INFINITY);
assert_eq!(logcdf_pos, 0.0);
assert_eq!(lambda_pos, 0.0);
let (logcdf_neg, lambda_neg) = signed_probit_logcdf_and_mills_ratio(f64::NEG_INFINITY);
assert_eq!(logcdf_neg, f64::NEG_INFINITY);
assert_eq!(lambda_neg, f64::INFINITY);
let (logcdf_nan, lambda_nan) = signed_probit_logcdf_and_mills_ratio(f64::NAN);
assert!(logcdf_nan.is_nan());
assert!(lambda_nan.is_nan());
}
#[test]
fn signed_probit_exact_derivative_helper_rejects_invalid_nonfinite_margins() {
assert_eq!(
signed_probit_neglog_derivatives_up_to_fourth(f64::INFINITY, 2.5)
.expect("+inf should use the zero tail"),
(0.0, 0.0, 0.0, 0.0)
);
let neg_inf_err = signed_probit_neglog_derivatives_up_to_fourth(f64::NEG_INFINITY, 2.5)
.expect_err("-inf should be rejected in the exact derivative path");
assert!(neg_inf_err.contains("non-finite signed margin"));
let nan_err = signed_probit_neglog_derivatives_up_to_fourth(f64::NAN, 2.5)
.expect_err("NaN should be rejected in the exact derivative path");
assert!(nan_err.contains("non-finite signed margin"));
}
#[test]
fn unary_neglog_phi_preserves_negative_infinity_and_nan_boundaries() {
assert_eq!(
unary_derivatives_neglog_phi(f64::INFINITY, 1.75),
[0.0, 0.0, 0.0, 0.0, 0.0]
);
assert_eq!(
unary_derivatives_neglog_phi(f64::NEG_INFINITY, 1.75),
[f64::INFINITY, f64::NEG_INFINITY, 1.75, 0.0, 0.0]
);
let nan_terms = unary_derivatives_neglog_phi(f64::NAN, 1.75);
assert!(nan_terms.iter().all(|value| value.is_nan()));
}
#[test]
fn flexible_family_routes_outer_derivatives_by_scale() {
let seed = array![-1.0, 0.0, 1.0];
let score_prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 3,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(3)),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
..default_test_family()
};
let specs = vec![
dummy_blockspec(1, 3),
dummy_blockspec(1, 3),
dummy_blockspec(2, 3),
];
assert_eq!(
family.exact_outer_derivative_order(&specs, &BlockwiseFitOptions::default()),
ExactOuterDerivativeOrder::Second
);
assert!(family.exact_newton_joint_psi_workspace_for_first_order_terms());
let n_large = 50_000;
let mut large_flex_family = family.clone();
large_flex_family.y = Arc::new(Array1::zeros(n_large));
large_flex_family.weights = Arc::new(Array1::ones(n_large));
large_flex_family.z = Arc::new(Array1::zeros(n_large));
let large_flex_specs = vec![
dummy_blockspec(1, n_large),
dummy_blockspec(1, n_large),
dummy_blockspec(2, n_large),
];
assert_eq!(
large_flex_family
.exact_outer_derivative_order(&large_flex_specs, &BlockwiseFitOptions::default()),
ExactOuterDerivativeOrder::Second
);
let (large_flex_gradient, large_flex_hessian) = custom_family_outer_derivatives(
&large_flex_family,
&large_flex_specs,
&BlockwiseFitOptions::default(),
);
assert_eq!(
large_flex_gradient,
crate::solver::outer_strategy::Derivative::Analytic
);
assert_eq!(
large_flex_hessian,
crate::solver::outer_strategy::DeclaredHessianForm::Either
);
let mut large_rigid_family = large_flex_family.clone();
large_rigid_family.score_warp = None;
let large_rigid_specs = vec![dummy_blockspec(1, n_large), dummy_blockspec(1, n_large)];
assert_eq!(
large_rigid_family
.exact_outer_derivative_order(&large_rigid_specs, &BlockwiseFitOptions::default()),
ExactOuterDerivativeOrder::Second
);
}
#[test]
fn exact_outer_order_stays_second_order_at_biobank_work_scale() {
use crate::custom_family::{
default_coefficient_hessian_cost, exact_outer_order_from_capability,
};
use crate::matrix::DesignMatrix;
use ndarray::Array2;
let small_design = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((10, 8)),
));
let small_specs: Vec<ParameterBlockSpec> = (0..2)
.map(|i| ParameterBlockSpec {
name: format!("block_{i}"),
design: small_design.clone(),
offset: ndarray::Array1::zeros(10),
penalties: (0..2)
.map(|_| crate::custom_family::PenaltyMatrix::Dense(Array2::zeros((8, 8))))
.collect(),
nullspace_dims: vec![0; 2],
initial_log_lambdas: ndarray::Array1::zeros(2),
initial_beta: None,
gauge_priority: 100,
jacobian_callback: None,
stacked_design: None,
stacked_offset: None,
})
.collect();
let small_cost = default_coefficient_hessian_cost(&small_specs);
assert_eq!(
exact_outer_order_from_capability(&small_specs, small_cost),
ExactOuterDerivativeOrder::Second
);
let big_design = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((5_000, 500)),
));
let big_specs: Vec<ParameterBlockSpec> = (0..2)
.map(|i| ParameterBlockSpec {
name: format!("block_{i}"),
design: big_design.clone(),
offset: ndarray::Array1::zeros(5_000),
penalties: (0..10)
.map(|_| crate::custom_family::PenaltyMatrix::Dense(Array2::zeros((500, 500))))
.collect(),
nullspace_dims: vec![0; 10],
initial_log_lambdas: ndarray::Array1::zeros(10),
initial_beta: None,
gauge_priority: 100,
jacobian_callback: None,
stacked_design: None,
stacked_offset: None,
})
.collect();
let big_cost = default_coefficient_hessian_cost(&big_specs);
assert_eq!(
exact_outer_order_from_capability(&big_specs, big_cost),
ExactOuterDerivativeOrder::Second
);
let ctn_specs = vec![ParameterBlockSpec {
name: "ctn".to_string(),
design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(Array2::zeros((
1, 1,
)))),
offset: ndarray::Array1::zeros(1),
penalties: (0..22)
.map(|_| crate::custom_family::PenaltyMatrix::Dense(Array2::zeros((1, 1))))
.collect(),
nullspace_dims: vec![0; 22],
initial_log_lambdas: ndarray::Array1::zeros(22),
initial_beta: None,
gauge_priority: 100,
jacobian_callback: None,
stacked_design: None,
stacked_offset: None,
}];
let ctn_cost: u64 = 400_000u64 * 300 * 300;
assert_eq!(
exact_outer_order_from_capability(&ctn_specs, ctn_cost),
ExactOuterDerivativeOrder::Second
);
use crate::custom_family::exact_outer_order_with_outer_hvp;
assert_eq!(
exact_outer_order_with_outer_hvp(&ctn_specs, ctn_cost, false),
ExactOuterDerivativeOrder::Second,
"exact outer order must not be cost-demoted"
);
let huge_k_specs = vec![ParameterBlockSpec {
name: "k".to_string(),
design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(Array2::zeros((
1, 1,
)))),
offset: ndarray::Array1::zeros(1),
penalties: (0..5_000)
.map(|_| crate::custom_family::PenaltyMatrix::Dense(Array2::zeros((1, 1))))
.collect(),
nullspace_dims: vec![0; 5_000],
initial_log_lambdas: ndarray::Array1::zeros(5_000),
initial_beta: None,
gauge_priority: 100,
jacobian_callback: None,
stacked_design: None,
stacked_offset: None,
}];
assert_eq!(
exact_outer_order_with_outer_hvp(&huge_k_specs, 0, true),
ExactOuterDerivativeOrder::Second,
"outer HVP support keeps the exact second-order declaration regardless of K"
);
}
#[test]
fn bernoulli_marginal_slope_coefficient_cost_uses_joint_coupled_formula() {
use crate::custom_family::default_coefficient_hessian_cost;
use crate::matrix::DesignMatrix;
let n = 1000usize;
let p_marg = 20usize;
let p_log = 8usize;
let marg_design = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((n, p_marg)),
));
let log_design = DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((n, p_log)),
));
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(Array1::zeros(n)),
weights: Arc::new(Array1::from_elem(n, 1.0)),
z: Arc::new(Array1::zeros(n)),
marginal_design: marg_design.clone(),
logslope_design: log_design.clone(),
..default_test_family()
};
let specs = vec![
ParameterBlockSpec {
name: "marginal".to_string(),
design: marg_design,
offset: Array1::zeros(n),
penalties: vec![crate::custom_family::PenaltyMatrix::Dense(Array2::zeros((
p_marg, p_marg,
)))],
nullspace_dims: vec![0],
initial_log_lambdas: Array1::zeros(1),
initial_beta: None,
gauge_priority: 100,
jacobian_callback: None,
stacked_design: None,
stacked_offset: None,
},
ParameterBlockSpec {
name: "logslope".to_string(),
design: log_design,
offset: Array1::zeros(n),
penalties: vec![crate::custom_family::PenaltyMatrix::Dense(Array2::zeros((
p_log, p_log,
)))],
nullspace_dims: vec![0],
initial_log_lambdas: Array1::zeros(1),
initial_beta: None,
gauge_priority: 100,
jacobian_callback: None,
stacked_design: None,
stacked_offset: None,
},
];
let p_total = (p_marg + p_log) as u64;
let expected_joint = (n as u64) * p_total * p_total;
let expected_block_diag = (n as u64) * ((p_marg * p_marg + p_log * p_log) as u64);
assert_eq!(family.coefficient_hessian_cost(&specs), expected_joint);
assert_eq!(
default_coefficient_hessian_cost(&specs),
expected_block_diag
);
assert!(family.coefficient_hessian_cost(&specs) > default_coefficient_hessian_cost(&specs));
}
#[test]
fn rigid_fast_path_matches_loglik_finite_differences() {
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![1.0]),
weights: Arc::new(array![1.2]),
z: Arc::new(array![0.3]),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(array![[
1.0
]])),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(array![[
1.0
]])),
..default_test_family()
};
let states_at = |q: f64, g: f64| {
vec![
ParameterBlockState {
beta: array![q],
eta: array![q],
},
ParameterBlockState {
beta: array![g],
eta: array![g],
},
]
};
let q = 0.4;
let g = 0.7;
let block_states = states_at(q, g);
let eval = family
.evaluate(&block_states)
.expect("rigid family evaluation");
let grad_q = match &eval.blockworking_sets[0] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
BlockWorkingSet::Diagonal { .. } => {
panic!("expected exact-newton marginal block")
}
};
let grad_g = match &eval.blockworking_sets[1] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
BlockWorkingSet::Diagonal { .. } => {
panic!("expected exact-newton log-slope block")
}
};
let hess_qq = match &eval.blockworking_sets[0] {
BlockWorkingSet::ExactNewton { hessian, .. } => match hessian {
SymmetricMatrix::Dense(h) => h[[0, 0]],
_ => panic!("expected dense marginal Hessian"),
},
BlockWorkingSet::Diagonal { .. } => {
panic!("expected exact-newton marginal block")
}
};
let hess_gg = match &eval.blockworking_sets[1] {
BlockWorkingSet::ExactNewton { hessian, .. } => match hessian {
SymmetricMatrix::Dense(h) => h[[0, 0]],
_ => panic!("expected dense log-slope Hessian"),
},
BlockWorkingSet::Diagonal { .. } => {
panic!("expected exact-newton log-slope block")
}
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("rigid exact eval cache");
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(0, &block_states, None, true)
.expect("rigid row context");
let (_, primary_grad, primary_hess) = family
.compute_row_primary_gradient_hessian(0, &block_states, &cache.primary, &row_ctx)
.expect("rigid exact row derivatives");
let expected_score_q =
BernoulliMarginalSlopeFamily::exact_newton_score_component_from_objective_gradient(
primary_grad[0],
);
let expected_score_g =
BernoulliMarginalSlopeFamily::exact_newton_score_component_from_objective_gradient(
primary_grad[1],
);
assert!(
(grad_q - expected_score_q).abs() < 1e-10,
"marginal gradient mismatch: fast={grad_q:.12e}, exact={expected_score_q:.12e}"
);
assert!(
(grad_g - expected_score_g).abs() < 1e-10,
"logslope gradient mismatch: fast={grad_g:.12e}, exact={expected_score_g:.12e}"
);
assert!(
(hess_qq - primary_hess[[0, 0]]).abs() < 1e-10,
"marginal Hessian mismatch: fast={hess_qq:.12e}, exact={:.12e}",
primary_hess[[0, 0]]
);
assert!(
(hess_gg - primary_hess[[1, 1]]).abs() < 1e-10,
"logslope Hessian mismatch: fast={hess_gg:.12e}, exact={:.12e}",
primary_hess[[1, 1]]
);
}
#[test]
fn w_only_gradient_hessian_finite_and_symmetric() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let beta_link = Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0)));
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(beta_link.clone(), seed.len()),
];
let slices = block_slices(&family);
assert!(slices.h.is_none(), "score-warp absent → no h slice");
let primary = primary_slices(&slices);
assert!(primary.h.is_none(), "primary h absent");
assert_eq!(primary.total, 2 + link_dim);
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let (_, grad, hess) = family
.compute_row_primary_gradient_hessian(row, &block_states, &primary, &row_ctx)
.unwrap_or_else(|e| {
panic!("row {row}: compute_row_primary_gradient_hessian failed: {e}")
});
assert_eq!(
grad.len(),
primary.total,
"row {row}: gradient length mismatch"
);
assert_eq!(
hess.dim(),
(primary.total, primary.total),
"row {row}: hessian shape mismatch"
);
assert!(
grad.iter().all(|v| v.is_finite()),
"row {row}: non-finite gradient entry: {grad:?}"
);
assert!(
hess.iter().all(|v| v.is_finite()),
"row {row}: non-finite hessian entry"
);
for a in 0..primary.total {
for b in 0..a {
let diff = (hess[[a, b]] - hess[[b, a]]).abs();
assert!(
diff < 1e-10,
"row {row}: hessian asymmetry at ({a},{b}): \
H[{a},{b}]={:.6e} vs H[{b},{a}]={:.6e}, diff={diff:.3e}",
hess[[a, b]],
hess[[b, a]]
);
}
}
}
}
#[test]
fn h_only_gradient_hessian_finite_and_symmetric() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let beta_score = Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0)));
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(beta_score, seed.len()),
];
let slices = block_slices(&family);
assert!(slices.w.is_none(), "link-dev absent → no w slice");
let primary = primary_slices(&slices);
assert!(primary.w.is_none(), "primary w absent");
assert_eq!(primary.total, 2 + score_dim);
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let (_, grad, hess) = family
.compute_row_primary_gradient_hessian(row, &block_states, &primary, &row_ctx)
.unwrap_or_else(|e| {
panic!("row {row}: compute_row_primary_gradient_hessian failed: {e}")
});
assert_eq!(
grad.len(),
primary.total,
"row {row}: gradient length mismatch"
);
assert_eq!(
hess.dim(),
(primary.total, primary.total),
"row {row}: hessian shape mismatch"
);
assert!(
grad.iter().all(|v| v.is_finite()),
"row {row}: non-finite gradient entry: {grad:?}"
);
assert!(
hess.iter().all(|v| v.is_finite()),
"row {row}: non-finite hessian entry"
);
for a in 0..primary.total {
for b in 0..a {
let diff = (hess[[a, b]] - hess[[b, a]]).abs();
assert!(
diff < 1e-10,
"row {row}: hessian asymmetry at ({a},{b}): \
H[{a},{b}]={:.6e} vs H[{b},{a}]={:.6e}, diff={diff:.3e}",
hess[[a, b]],
hess[[b, a]]
);
}
}
}
}
#[test]
fn w_only_exact_outer_directional_derivatives_are_present_and_finite() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let beta_link = Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0)));
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(beta_link, seed.len()),
];
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
let w_range = slices.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.15;
if w_range.len() > 1 {
dir_u[w_range.start + 1] = -0.07;
}
dir_v[w_range.start] = 0.09;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = 0.03;
}
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_u)
.expect("w-only third directional derivative")
.expect("w-only third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_v)
.expect("w-only fourth directional derivative")
.expect("w-only fourth directional derivative matrix");
assert_eq!(third.dim(), (total, total));
assert_eq!(fourth.dim(), (total, total));
assert!(third.iter().all(|value| value.is_finite()));
assert!(fourth.iter().all(|value| value.is_finite()));
let max_abs_third = third
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs()));
let max_abs_fourth = fourth
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs()));
assert!(
max_abs_third > 1e-10,
"expected nonzero w-only third directional derivative"
);
assert!(
max_abs_fourth > 1e-10,
"expected nonzero w-only fourth directional derivative"
);
for i in 0..total {
for j in 0..i {
assert!((third[[i, j]] - third[[j, i]]).abs() < 1e-8);
assert!((fourth[[i, j]] - fourth[[j, i]]).abs() < 1e-8);
}
}
}
#[test]
fn h_only_exact_outer_directional_derivatives_are_present_and_finite() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let beta_score = Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0)));
let scalar_design = || {
DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(Array2::from_elem(
(seed.len(), 1),
1.0,
)))
};
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: scalar_design(),
logslope_design: scalar_design(),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(seed.len(), 0.25),
},
ParameterBlockState {
beta: array![0.15],
eta: Array1::from_elem(seed.len(), 0.15),
},
ParameterBlockState {
beta: beta_score,
eta: Array1::zeros(seed.len()),
},
];
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[slices.marginal.start] = -0.35;
dir_u[slices.logslope.start] = 0.28;
let h_range = slices.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.12;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.06;
}
dir_v[slices.marginal.start] = 0.18;
dir_v[slices.logslope.start] = -0.22;
dir_v[h_range.start] = 0.07;
if h_range.len() > 1 {
dir_v[h_range.start + 1] = 0.05;
}
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_u)
.expect("h-only third directional derivative")
.expect("h-only third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_v)
.expect("h-only fourth directional derivative")
.expect("h-only fourth directional derivative matrix");
assert_eq!(third.dim(), (total, total));
assert_eq!(fourth.dim(), (total, total));
assert!(third.iter().all(|value| value.is_finite()));
assert!(fourth.iter().all(|value| value.is_finite()));
let max_abs_third = third
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs()));
let max_abs_fourth = fourth
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs()));
assert!(
max_abs_third > 1e-10,
"expected nonzero h-only third directional derivative"
);
assert!(
max_abs_fourth > 1e-10,
"expected nonzero h-only fourth directional derivative"
);
for i in 0..total {
for j in 0..i {
assert!((third[[i, j]] - third[[j, i]]).abs() < 1e-8);
assert!((fourth[[i, j]] - fourth[[j, i]]).abs() < 1e-8);
}
}
}
#[test]
fn h_only_row_primary_higher_order_contractions_are_finite_and_symmetric() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let beta_score = Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0)));
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(beta_score, seed.len()),
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = -0.35;
dir_u[cache.primary.logslope] = 0.28;
let h_range = cache.primary.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.12;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.06;
}
dir_v[cache.primary.q] = 0.18;
dir_v[cache.primary.logslope] = -0.22;
dir_v[h_range.start] = 0.07;
if h_range.len() > 1 {
dir_v[h_range.start + 1] = 0.05;
}
let mut max_abs_third = 0.0_f64;
let mut max_abs_fourth = 0.0_f64;
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_third_contracted_recompute failed: {e}")
});
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
&dir_v,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_fourth_contracted_recompute failed: {e}")
});
assert_eq!(third.dim(), (total, total));
assert_eq!(fourth.dim(), (total, total));
assert!(third.iter().all(|value| value.is_finite()));
assert!(fourth.iter().all(|value| value.is_finite()));
max_abs_third = max_abs_third.max(
third
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs())),
);
max_abs_fourth = max_abs_fourth.max(
fourth
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs())),
);
for i in 0..total {
for j in 0..i {
assert!((third[[i, j]] - third[[j, i]]).abs() < 1e-8);
assert!((fourth[[i, j]] - fourth[[j, i]]).abs() < 1e-8);
}
}
}
assert!(
max_abs_third > 1e-10,
"expected nonzero h-only third contraction"
);
assert!(
max_abs_fourth > 1e-10,
"expected nonzero h-only fourth contraction"
);
}
#[test]
fn w_only_row_primary_higher_order_contractions_are_finite_and_symmetric() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let beta_link = Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0)));
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(beta_link, seed.len()),
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = 0.4;
dir_u[cache.primary.logslope] = -0.3;
let w_range = cache.primary.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.15;
if w_range.len() > 1 {
dir_u[w_range.start + 1] = -0.07;
}
dir_v[cache.primary.q] = -0.2;
dir_v[cache.primary.logslope] = 0.25;
dir_v[w_range.start] = 0.09;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = 0.03;
}
let mut max_abs_third = 0.0_f64;
let mut max_abs_fourth = 0.0_f64;
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_third_contracted_recompute failed: {e}")
});
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
&dir_v,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_fourth_contracted_recompute failed: {e}")
});
assert_eq!(third.dim(), (total, total));
assert_eq!(fourth.dim(), (total, total));
assert!(third.iter().all(|value| value.is_finite()));
assert!(fourth.iter().all(|value| value.is_finite()));
max_abs_third = max_abs_third.max(
third
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs())),
);
max_abs_fourth = max_abs_fourth.max(
fourth
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs())),
);
for i in 0..total {
for j in 0..i {
assert!((third[[i, j]] - third[[j, i]]).abs() < 1e-8);
assert!((fourth[[i, j]] - fourth[[j, i]]).abs() < 1e-8);
}
}
}
assert!(
max_abs_third > 1e-10,
"expected nonzero w-only third contraction"
);
assert!(
max_abs_fourth > 1e-10,
"expected nonzero w-only fourth contraction"
);
}
#[test]
fn dual_flex_row_primary_higher_order_contractions_are_finite_and_symmetric() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = 0.7;
dir_u[cache.primary.logslope] = -0.2;
let h_range = cache.primary.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = cache.primary.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[cache.primary.q] = -0.4;
dir_v[cache.primary.logslope] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
let mut max_abs_third = 0.0_f64;
let mut max_abs_fourth = 0.0_f64;
for row in 0..z.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_third_contracted_recompute failed: {e}")
});
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
&dir_v,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_fourth_contracted_recompute failed: {e}")
});
assert_eq!(third.dim(), (total, total));
assert_eq!(fourth.dim(), (total, total));
assert!(third.iter().all(|value| value.is_finite()));
assert!(fourth.iter().all(|value| value.is_finite()));
max_abs_third = max_abs_third.max(
third
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs())),
);
max_abs_fourth = max_abs_fourth.max(
fourth
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs())),
);
for i in 0..total {
for j in 0..i {
assert!((third[[i, j]] - third[[j, i]]).abs() < 1e-8);
assert!((fourth[[i, j]] - fourth[[j, i]]).abs() < 1e-8);
}
}
}
assert!(
max_abs_third > 1e-10,
"expected nonzero dual-flex third contraction"
);
assert!(
max_abs_fourth > 1e-10,
"expected nonzero dual-flex fourth contraction"
);
}
#[test]
fn dual_flex_row_primary_higher_order_zero_direction_returns_zero() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let zero = Array1::<f64>::zeros(cache.primary.total);
for row in 0..z.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(row, &block_states, &cache, &row_ctx, &zero)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_third_contracted_recompute failed: {e}")
});
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&zero,
&zero,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_fourth_contracted_recompute failed: {e}")
});
assert!(
third.iter().all(|value| value.abs() <= 0.0),
"row {row}: expected zero third contraction for zero direction"
);
assert!(
fourth.iter().all(|value| value.abs() <= 0.0),
"row {row}: expected zero fourth contraction for zero directions"
);
}
}
#[test]
fn h_only_row_primary_higher_order_zero_direction_returns_zero() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let zero = Array1::<f64>::zeros(cache.primary.total);
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(row, &block_states, &cache, &row_ctx, &zero)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_third_contracted_recompute failed: {e}")
});
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&zero,
&zero,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_fourth_contracted_recompute failed: {e}")
});
assert!(
third.iter().all(|value| value.abs() <= 0.0),
"row {row}: expected zero h-only third contraction for zero direction"
);
assert!(
fourth.iter().all(|value| value.abs() <= 0.0),
"row {row}: expected zero h-only fourth contraction for zero directions"
);
}
}
#[test]
fn w_only_row_primary_higher_order_zero_direction_returns_zero() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let zero = Array1::<f64>::zeros(cache.primary.total);
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(row, &block_states, &cache, &row_ctx, &zero)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_third_contracted_recompute failed: {e}")
});
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&zero,
&zero,
)
.unwrap_or_else(|e| {
panic!("row {row}: row_primary_fourth_contracted_recompute failed: {e}")
});
assert!(
third.iter().all(|value| value.abs() <= 0.0),
"row {row}: expected zero w-only third contraction for zero direction"
);
assert!(
fourth.iter().all(|value| value.abs() <= 0.0),
"row {row}: expected zero w-only fourth contraction for zero directions"
);
}
}
#[test]
fn dual_flex_exact_outer_zero_direction_returns_zero() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let slices = block_slices(&family);
let zero = Array1::<f64>::zeros(slices.total);
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &zero)
.expect("dual-flex third directional derivative")
.expect("dual-flex third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &zero, &zero)
.expect("dual-flex fourth directional derivative")
.expect("dual-flex fourth directional derivative matrix");
assert!(
third.iter().all(|value| value.abs() <= 0.0),
"expected zero dual-flex third directional derivative for zero direction"
);
assert!(
fourth.iter().all(|value| value.abs() <= 0.0),
"expected zero dual-flex fourth directional derivative for zero directions"
);
}
#[test]
fn dual_flex_exact_outer_fourth_direction_swap_is_symmetric() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[slices.marginal.start] = 0.7;
dir_u[slices.logslope.start] = -0.2;
let h_range = slices.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = slices.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[slices.marginal.start] = -0.4;
dir_v[slices.logslope.start] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
let forward = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_v)
.expect("dual-flex fourth directional derivative")
.expect("dual-flex fourth directional derivative matrix");
let swapped = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_v, &dir_u)
.expect("dual-flex swapped fourth directional derivative")
.expect("dual-flex swapped fourth directional derivative matrix");
assert_eq!(forward.dim(), (total, total));
assert_eq!(swapped.dim(), (total, total));
for i in 0..total {
for j in 0..total {
assert!(
(forward[[i, j]] - swapped[[i, j]]).abs() < 1e-8,
"fourth directional derivative should be symmetric in direction arguments at ({i},{j})"
);
}
}
}
#[test]
fn dual_flex_row_primary_fourth_direction_swap_is_symmetric() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = 0.7;
dir_u[cache.primary.logslope] = -0.2;
let h_range = cache.primary.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = cache.primary.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[cache.primary.q] = -0.4;
dir_v[cache.primary.logslope] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
for row in 0..z.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let forward = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
&dir_v,
)
.unwrap_or_else(|e| panic!("row {row}: forward fourth contraction failed: {e}"));
let swapped = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_v,
&dir_u,
)
.unwrap_or_else(|e| panic!("row {row}: swapped fourth contraction failed: {e}"));
assert_eq!(forward.dim(), (total, total));
assert_eq!(swapped.dim(), (total, total));
for i in 0..total {
for j in 0..total {
assert!(
(forward[[i, j]] - swapped[[i, j]]).abs() < 1e-8,
"row {row}: fourth contraction should be symmetric in direction arguments at ({i},{j})"
);
}
}
}
}
#[test]
fn dual_flex_row_primary_higher_order_direction_sign_rules_hold() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = 0.7;
dir_u[cache.primary.logslope] = -0.2;
let h_range = cache.primary.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = cache.primary.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[cache.primary.q] = -0.4;
dir_v[cache.primary.logslope] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
let neg_dir_u = dir_u.mapv(|value| -value);
for row in 0..z.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction failed: {e}"));
let third_neg = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&neg_dir_u,
)
.unwrap_or_else(|e| panic!("row {row}: negated third contraction failed: {e}"));
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
&dir_v,
)
.unwrap_or_else(|e| panic!("row {row}: fourth contraction failed: {e}"));
let fourth_neg_u = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&neg_dir_u,
&dir_v,
)
.unwrap_or_else(|e| panic!("row {row}: negated-u fourth contraction failed: {e}"));
for i in 0..total {
for j in 0..total {
assert!(
(third_neg[[i, j]] + third[[i, j]]).abs() < 1e-8,
"row {row}: third contraction should be odd in its direction at ({i},{j})"
);
assert!(
(fourth_neg_u[[i, j]] + fourth[[i, j]]).abs() < 1e-8,
"row {row}: fourth contraction should be linear in dir_u sign at ({i},{j})"
);
}
}
}
}
#[test]
fn h_only_row_primary_fourth_direction_swap_is_symmetric() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = -0.35;
dir_u[cache.primary.logslope] = 0.28;
let h_range = cache.primary.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.12;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.06;
}
dir_v[cache.primary.q] = 0.18;
dir_v[cache.primary.logslope] = -0.22;
dir_v[h_range.start] = 0.07;
if h_range.len() > 1 {
dir_v[h_range.start + 1] = 0.05;
}
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let forward = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
&dir_v,
)
.unwrap_or_else(|e| panic!("row {row}: forward fourth contraction failed: {e}"));
let swapped = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_v,
&dir_u,
)
.unwrap_or_else(|e| panic!("row {row}: swapped fourth contraction failed: {e}"));
for i in 0..total {
for j in 0..total {
assert!(
(forward[[i, j]] - swapped[[i, j]]).abs() < 1e-8,
"row {row}: h-only fourth contraction should be symmetric in direction arguments at ({i},{j})"
);
}
}
}
}
#[test]
fn w_only_row_primary_fourth_direction_swap_is_symmetric() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = 0.4;
dir_u[cache.primary.logslope] = -0.3;
let w_range = cache.primary.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.15;
if w_range.len() > 1 {
dir_u[w_range.start + 1] = -0.07;
}
dir_v[cache.primary.q] = -0.2;
dir_v[cache.primary.logslope] = 0.25;
dir_v[w_range.start] = 0.09;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = 0.03;
}
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let forward = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
&dir_v,
)
.unwrap_or_else(|e| panic!("row {row}: forward fourth contraction failed: {e}"));
let swapped = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_v,
&dir_u,
)
.unwrap_or_else(|e| panic!("row {row}: swapped fourth contraction failed: {e}"));
for i in 0..total {
for j in 0..total {
assert!(
(forward[[i, j]] - swapped[[i, j]]).abs() < 1e-8,
"row {row}: w-only fourth contraction should be symmetric in direction arguments at ({i},{j})"
);
}
}
}
}
#[test]
fn h_only_row_primary_higher_order_direction_sign_rules_hold() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir = Array1::<f64>::zeros(total);
dir[cache.primary.q] = -0.35;
dir[cache.primary.logslope] = 0.28;
let h_range = cache.primary.h.as_ref().expect("h slice");
dir[h_range.start] = 0.12;
if h_range.len() > 1 {
dir[h_range.start + 1] = -0.06;
}
let neg_dir = dir.mapv(|value| -value);
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(row, &block_states, &cache, &row_ctx, &dir)
.unwrap_or_else(|e| panic!("row {row}: third contraction failed: {e}"));
let third_neg = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&neg_dir,
)
.unwrap_or_else(|e| panic!("row {row}: negated third contraction failed: {e}"));
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir,
&dir,
)
.unwrap_or_else(|e| panic!("row {row}: fourth contraction failed: {e}"));
let fourth_neg = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&neg_dir,
&neg_dir,
)
.unwrap_or_else(|e| {
panic!("row {row}: doubly-negated fourth contraction failed: {e}")
});
for i in 0..total {
for j in 0..total {
assert!(
(third_neg[[i, j]] + third[[i, j]]).abs() < 1e-8,
"row {row}: h-only third contraction should be odd at ({i},{j})"
);
assert!(
(fourth_neg[[i, j]] - fourth[[i, j]]).abs() < 1e-8,
"row {row}: h-only fourth contraction should be invariant under flipping both directions at ({i},{j})"
);
}
}
}
}
#[test]
fn w_only_row_primary_higher_order_direction_sign_rules_hold() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir = Array1::<f64>::zeros(total);
dir[cache.primary.q] = 0.4;
dir[cache.primary.logslope] = -0.3;
let w_range = cache.primary.w.as_ref().expect("w slice");
dir[w_range.start] = 0.15;
if w_range.len() > 1 {
dir[w_range.start + 1] = -0.07;
}
let neg_dir = dir.mapv(|value| -value);
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third = family
.row_primary_third_contracted_recompute(row, &block_states, &cache, &row_ctx, &dir)
.unwrap_or_else(|e| panic!("row {row}: third contraction failed: {e}"));
let third_neg = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&neg_dir,
)
.unwrap_or_else(|e| panic!("row {row}: negated third contraction failed: {e}"));
let fourth = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir,
&dir,
)
.unwrap_or_else(|e| panic!("row {row}: fourth contraction failed: {e}"));
let fourth_neg = family
.row_primary_fourth_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&neg_dir,
&neg_dir,
)
.unwrap_or_else(|e| {
panic!("row {row}: doubly-negated fourth contraction failed: {e}")
});
for i in 0..total {
for j in 0..total {
assert!(
(third_neg[[i, j]] + third[[i, j]]).abs() < 1e-8,
"row {row}: w-only third contraction should be odd at ({i},{j})"
);
assert!(
(fourth_neg[[i, j]] - fourth[[i, j]]).abs() < 1e-8,
"row {row}: w-only fourth contraction should be invariant under flipping both directions at ({i},{j})"
);
}
}
}
}
#[test]
fn dual_flex_exact_outer_direction_sign_rules_hold() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[slices.marginal.start] = 0.7;
dir_u[slices.logslope.start] = -0.2;
let h_range = slices.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = slices.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[slices.marginal.start] = -0.4;
dir_v[slices.logslope.start] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
let neg_dir_u = dir_u.mapv(|value| -value);
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_u)
.expect("dual-flex third directional derivative")
.expect("dual-flex third directional derivative matrix");
let third_neg = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &neg_dir_u)
.expect("dual-flex negated third directional derivative")
.expect("dual-flex negated third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_v)
.expect("dual-flex fourth directional derivative")
.expect("dual-flex fourth directional derivative matrix");
let fourth_neg_u = family
.exact_newton_joint_hessiansecond_directional_derivative(
&block_states,
&neg_dir_u,
&dir_v,
)
.expect("dual-flex negated-u fourth directional derivative")
.expect("dual-flex negated-u fourth directional derivative matrix");
assert_eq!(third.dim(), (total, total));
assert_eq!(third_neg.dim(), (total, total));
assert_eq!(fourth.dim(), (total, total));
assert_eq!(fourth_neg_u.dim(), (total, total));
for i in 0..total {
for j in 0..total {
assert!(
(third_neg[[i, j]] + third[[i, j]]).abs() < 1e-8,
"third directional derivative should be odd in its direction at ({i},{j})"
);
assert!(
(fourth_neg_u[[i, j]] + fourth[[i, j]]).abs() < 1e-8,
"fourth directional derivative should be linear in dir_u sign at ({i},{j})"
);
}
}
}
#[test]
fn dual_flex_exact_outer_fourth_double_sign_flip_is_invariant() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[slices.marginal.start] = 0.7;
dir_u[slices.logslope.start] = -0.2;
let h_range = slices.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = slices.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[slices.marginal.start] = -0.4;
dir_v[slices.logslope.start] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
let neg_dir_u = dir_u.mapv(|value| -value);
let neg_dir_v = dir_v.mapv(|value| -value);
let forward = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_v)
.expect("dual-flex fourth directional derivative")
.expect("dual-flex fourth directional derivative matrix");
let flipped = family
.exact_newton_joint_hessiansecond_directional_derivative(
&block_states,
&neg_dir_u,
&neg_dir_v,
)
.expect("dual-flex doubly-negated fourth directional derivative")
.expect("dual-flex doubly-negated fourth directional derivative matrix");
assert_eq!(forward.dim(), (total, total));
assert_eq!(flipped.dim(), (total, total));
for i in 0..total {
for j in 0..total {
assert!(
(forward[[i, j]] - flipped[[i, j]]).abs() < 1e-8,
"fourth directional derivative should be invariant under flipping both directions at ({i},{j})"
);
}
}
}
#[test]
fn dual_flex_exact_outer_third_direction_is_linear() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[slices.marginal.start] = 0.7;
dir_u[slices.logslope.start] = -0.2;
let h_range = slices.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = slices.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[slices.marginal.start] = -0.4;
dir_v[slices.logslope.start] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
let dir_sum = &dir_u + &dir_v;
let third_u = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_u)
.expect("dual-flex third directional derivative u")
.expect("dual-flex third directional derivative u matrix");
let third_v = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_v)
.expect("dual-flex third directional derivative v")
.expect("dual-flex third directional derivative v matrix");
let third_sum = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_sum)
.expect("dual-flex third directional derivative sum")
.expect("dual-flex third directional derivative sum matrix");
for i in 0..total {
for j in 0..total {
let expected = third_u[[i, j]] + third_v[[i, j]];
assert!(
(third_sum[[i, j]] - expected).abs() < 1e-8,
"third directional derivative should be linear in its direction at ({i},{j})"
);
}
}
}
#[test]
fn dual_flex_row_primary_third_direction_is_linear() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = 0.7;
dir_u[cache.primary.logslope] = -0.2;
let h_range = cache.primary.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = cache.primary.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[cache.primary.q] = -0.4;
dir_v[cache.primary.logslope] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
let dir_sum = &dir_u + &dir_v;
for row in 0..z.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third_u = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction u failed: {e}"));
let third_v = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_v,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction v failed: {e}"));
let third_sum = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_sum,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction sum failed: {e}"));
for i in 0..total {
for j in 0..total {
let expected = third_u[[i, j]] + third_v[[i, j]];
assert!(
(third_sum[[i, j]] - expected).abs() < 1e-8,
"row {row}: third contraction should be linear in its direction at ({i},{j})"
);
}
}
}
}
#[test]
fn h_only_row_primary_third_direction_is_linear() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = -0.35;
dir_u[cache.primary.logslope] = 0.28;
let h_range = cache.primary.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.12;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.06;
}
dir_v[cache.primary.q] = 0.18;
dir_v[cache.primary.logslope] = -0.22;
dir_v[h_range.start] = 0.07;
if h_range.len() > 1 {
dir_v[h_range.start + 1] = 0.05;
}
let dir_sum = &dir_u + &dir_v;
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third_u = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction u failed: {e}"));
let third_v = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_v,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction v failed: {e}"));
let third_sum = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_sum,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction sum failed: {e}"));
for i in 0..total {
for j in 0..total {
let expected = third_u[[i, j]] + third_v[[i, j]];
assert!(
(third_sum[[i, j]] - expected).abs() < 1e-8,
"row {row}: h-only third contraction should be linear at ({i},{j})"
);
}
}
}
}
#[test]
fn w_only_row_primary_third_direction_is_linear() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let cache = family
.build_exact_eval_cache(&block_states)
.expect("exact eval cache");
let total = cache.primary.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[cache.primary.q] = 0.4;
dir_u[cache.primary.logslope] = -0.3;
let w_range = cache.primary.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.15;
if w_range.len() > 1 {
dir_u[w_range.start + 1] = -0.07;
}
dir_v[cache.primary.q] = -0.2;
dir_v[cache.primary.logslope] = 0.25;
dir_v[w_range.start] = 0.09;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = 0.03;
}
let dir_sum = &dir_u + &dir_v;
for row in 0..seed.len() {
let row_ctx = family
.build_row_exact_context_with_stats_and_cell_cache(row, &block_states, None, true)
.unwrap_or_else(|e| panic!("row {row}: build_row_exact_context failed: {e}"));
let third_u = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_u,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction u failed: {e}"));
let third_v = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_v,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction v failed: {e}"));
let third_sum = family
.row_primary_third_contracted_recompute(
row,
&block_states,
&cache,
&row_ctx,
&dir_sum,
)
.unwrap_or_else(|e| panic!("row {row}: third contraction sum failed: {e}"));
for i in 0..total {
for j in 0..total {
let expected = third_u[[i, j]] + third_v[[i, j]];
assert!(
(third_sum[[i, j]] - expected).abs() < 1e-8,
"row {row}: w-only third contraction should be linear at ({i},{j})"
);
}
}
}
}
#[test]
fn dual_flex_exact_outer_fourth_first_direction_is_linear() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
let mut dir_w = Array1::<f64>::zeros(total);
dir_u[slices.marginal.start] = 0.7;
dir_u[slices.logslope.start] = -0.2;
let h_range = slices.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
let w_range = slices.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.08;
dir_v[slices.marginal.start] = -0.4;
dir_v[slices.logslope.start] = 0.3;
dir_v[h_range.start] = -0.03;
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
dir_w[slices.marginal.start] = 0.11;
dir_w[slices.logslope.start] = -0.09;
dir_w[h_range.start] = 0.04;
dir_w[w_range.start] = -0.05;
let dir_sum = &dir_u + &dir_v;
let fourth_u = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_w)
.expect("dual-flex fourth directional derivative u,w")
.expect("dual-flex fourth directional derivative u,w matrix");
let fourth_v = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_v, &dir_w)
.expect("dual-flex fourth directional derivative v,w")
.expect("dual-flex fourth directional derivative v,w matrix");
let fourth_sum = family
.exact_newton_joint_hessiansecond_directional_derivative(
&block_states,
&dir_sum,
&dir_w,
)
.expect("dual-flex fourth directional derivative (u+v),w")
.expect("dual-flex fourth directional derivative (u+v),w matrix");
for i in 0..total {
for j in 0..total {
let expected = fourth_u[[i, j]] + fourth_v[[i, j]];
assert!(
(fourth_sum[[i, j]] - expected).abs() < 1e-8,
"fourth directional derivative should be linear in its first direction at ({i},{j})"
);
}
}
}
#[test]
fn h_only_exact_outer_third_direction_is_linear() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
let h_range = slices.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.12;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.06;
}
dir_v[h_range.start] = 0.07;
if h_range.len() > 1 {
dir_v[h_range.start + 1] = 0.05;
}
let dir_sum = &dir_u + &dir_v;
let third_u = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_u)
.expect("h-only third directional derivative u")
.expect("h-only third directional derivative u matrix");
let third_v = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_v)
.expect("h-only third directional derivative v")
.expect("h-only third directional derivative v matrix");
let third_sum = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_sum)
.expect("h-only third directional derivative sum")
.expect("h-only third directional derivative sum matrix");
for i in 0..total {
for j in 0..total {
let expected = third_u[[i, j]] + third_v[[i, j]];
assert!(
(third_sum[[i, j]] - expected).abs() < 1e-8,
"h-only third directional derivative should be linear at ({i},{j})"
);
}
}
}
#[test]
fn w_only_exact_outer_third_direction_is_linear() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
let w_range = slices.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.15;
if w_range.len() > 1 {
dir_u[w_range.start + 1] = -0.07;
}
dir_v[w_range.start] = 0.09;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = 0.03;
}
let dir_sum = &dir_u + &dir_v;
let third_u = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_u)
.expect("w-only third directional derivative u")
.expect("w-only third directional derivative u matrix");
let third_v = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_v)
.expect("w-only third directional derivative v")
.expect("w-only third directional derivative v matrix");
let third_sum = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_sum)
.expect("w-only third directional derivative sum")
.expect("w-only third directional derivative sum matrix");
for i in 0..total {
for j in 0..total {
let expected = third_u[[i, j]] + third_v[[i, j]];
assert!(
(third_sum[[i, j]] - expected).abs() < 1e-8,
"w-only third directional derivative should be linear at ({i},{j})"
);
}
}
}
#[test]
fn h_only_exact_outer_direction_sign_rules_hold() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let slices = block_slices(&family);
let total = slices.total;
let mut dir = Array1::<f64>::zeros(total);
let h_range = slices.h.as_ref().expect("h slice");
dir[h_range.start] = 0.12;
if h_range.len() > 1 {
dir[h_range.start + 1] = -0.06;
}
let neg_dir = dir.mapv(|value| -value);
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir)
.expect("h-only third directional derivative")
.expect("h-only third directional derivative matrix");
let third_neg = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &neg_dir)
.expect("h-only negated third directional derivative")
.expect("h-only negated third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir, &dir)
.expect("h-only fourth directional derivative")
.expect("h-only fourth directional derivative matrix");
let fourth_neg = family
.exact_newton_joint_hessiansecond_directional_derivative(
&block_states,
&neg_dir,
&neg_dir,
)
.expect("h-only doubly-negated fourth directional derivative")
.expect("h-only doubly-negated fourth directional derivative matrix");
for i in 0..total {
for j in 0..total {
assert!(
(third_neg[[i, j]] + third[[i, j]]).abs() < 1e-8,
"h-only third directional derivative should be odd at ({i},{j})"
);
assert!(
(fourth_neg[[i, j]] - fourth[[i, j]]).abs() < 1e-8,
"h-only fourth directional derivative should be invariant under flipping both directions at ({i},{j})"
);
}
}
}
#[test]
fn w_only_exact_outer_direction_sign_rules_hold() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let slices = block_slices(&family);
let total = slices.total;
let mut dir = Array1::<f64>::zeros(total);
let w_range = slices.w.as_ref().expect("w slice");
dir[w_range.start] = 0.15;
if w_range.len() > 1 {
dir[w_range.start + 1] = -0.07;
}
let neg_dir = dir.mapv(|value| -value);
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir)
.expect("w-only third directional derivative")
.expect("w-only third directional derivative matrix");
let third_neg = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &neg_dir)
.expect("w-only negated third directional derivative")
.expect("w-only negated third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir, &dir)
.expect("w-only fourth directional derivative")
.expect("w-only fourth directional derivative matrix");
let fourth_neg = family
.exact_newton_joint_hessiansecond_directional_derivative(
&block_states,
&neg_dir,
&neg_dir,
)
.expect("w-only doubly-negated fourth directional derivative")
.expect("w-only doubly-negated fourth directional derivative matrix");
for i in 0..total {
for j in 0..total {
assert!(
(third_neg[[i, j]] + third[[i, j]]).abs() < 1e-8,
"w-only third directional derivative should be odd at ({i},{j})"
);
assert!(
(fourth_neg[[i, j]] - fourth[[i, j]]).abs() < 1e-8,
"w-only fourth directional derivative should be invariant under flipping both directions at ({i},{j})"
);
}
}
}
#[test]
fn h_only_exact_outer_fourth_direction_swap_is_symmetric() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
let h_range = slices.h.as_ref().expect("h slice");
dir_u[h_range.start] = 0.12;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.06;
}
dir_v[h_range.start] = 0.07;
if h_range.len() > 1 {
dir_v[h_range.start + 1] = 0.05;
}
let forward = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_v)
.expect("h-only fourth directional derivative")
.expect("h-only fourth directional derivative matrix");
let swapped = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_v, &dir_u)
.expect("h-only swapped fourth directional derivative")
.expect("h-only swapped fourth directional derivative matrix");
for i in 0..total {
for j in 0..total {
assert!(
(forward[[i, j]] - swapped[[i, j]]).abs() < 1e-8,
"h-only fourth directional derivative should be symmetric in direction arguments at ({i},{j})"
);
}
}
}
#[test]
fn w_only_exact_outer_fourth_direction_swap_is_symmetric() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
let w_range = slices.w.as_ref().expect("w slice");
dir_u[w_range.start] = 0.15;
if w_range.len() > 1 {
dir_u[w_range.start + 1] = -0.07;
}
dir_v[w_range.start] = 0.09;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = 0.03;
}
let forward = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_v)
.expect("w-only fourth directional derivative")
.expect("w-only fourth directional derivative matrix");
let swapped = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_v, &dir_u)
.expect("w-only swapped fourth directional derivative")
.expect("w-only swapped fourth directional derivative matrix");
for i in 0..total {
for j in 0..total {
assert!(
(forward[[i, j]] - swapped[[i, j]]).abs() < 1e-8,
"w-only fourth directional derivative should be symmetric in direction arguments at ({i},{j})"
);
}
}
}
#[test]
fn h_only_exact_outer_zero_direction_returns_zero() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_score_warp_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build score-warp block");
let score_dim = prepared
.block
.initial_beta
.as_ref()
.expect("score-warp initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..score_dim).map(|idx| 0.04 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
score_warp: Some(prepared.runtime.clone()),
..default_test_family()
};
let slices = block_slices(&family);
let zero = Array1::<f64>::zeros(slices.total);
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &zero)
.expect("h-only third directional derivative")
.expect("h-only third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &zero, &zero)
.expect("h-only fourth directional derivative")
.expect("h-only fourth directional derivative matrix");
assert!(
third.iter().all(|value| value.abs() <= 0.0),
"expected zero h-only third directional derivative for zero direction"
);
assert!(
fourth.iter().all(|value| value.abs() <= 0.0),
"expected zero h-only fourth directional derivative for zero directions"
);
}
#[test]
fn w_only_exact_outer_zero_direction_returns_zero() {
let seed = array![-1.5, -0.5, 0.0, 0.5, 1.5];
let prepared = build_test_link_deviation_block_from_seed(
&seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let link_dim = prepared
.block
.initial_beta
.as_ref()
.expect("link initial beta")
.len();
let block_states = vec![
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(array![0.0], seed.len()),
dummy_block_state(
Array1::from_iter((0..link_dim).map(|idx| 0.05 * (idx as f64 + 1.0))),
seed.len(),
),
];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(array![0.0, 1.0, 0.0, 1.0, 0.0]),
weights: Arc::new(Array1::ones(seed.len())),
z: Arc::new(seed.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
Array2::zeros((seed.len(), 0)),
)),
link_dev: Some(prepared.runtime.clone()),
..default_test_family()
};
let slices = block_slices(&family);
let zero = Array1::<f64>::zeros(slices.total);
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &zero)
.expect("w-only third directional derivative")
.expect("w-only third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &zero, &zero)
.expect("w-only fourth directional derivative")
.expect("w-only fourth directional derivative matrix");
assert!(
third.iter().all(|value| value.abs() <= 0.0),
"expected zero w-only third directional derivative for zero direction"
);
assert!(
fourth.iter().all(|value| value.abs() <= 0.0),
"expected zero w-only fourth directional derivative for zero directions"
);
}
#[test]
fn w_only_gradient_matches_loglik_finite_differences() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let beta_w = Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
);
let states_at = |q: f64, b: f64, bw: Array1<f64>| {
vec![
ParameterBlockState {
beta: array![q],
eta: Array1::from_elem(z.len(), q),
},
ParameterBlockState {
beta: array![b],
eta: Array1::from_elem(z.len(), b),
},
ParameterBlockState {
beta: bw,
eta: Array1::zeros(z.len()),
},
]
};
let q0 = 0.25;
let b0 = 0.6;
let block_states = states_at(q0, b0, beta_w.clone());
let eval = family.evaluate(&block_states).expect("family evaluation");
let grad_q = match &eval.blockworking_sets[0] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton q block"),
};
let grad_b = match &eval.blockworking_sets[1] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton b block"),
};
let grad_w0 = match &eval.blockworking_sets[2] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton w block"),
};
let fd = |which: &str, eps: f64| match which {
"q" => {
let plus = family
.log_likelihood_only(&states_at(q0 + eps, b0, beta_w.clone()))
.expect("ll plus q");
let minus = family
.log_likelihood_only(&states_at(q0 - eps, b0, beta_w.clone()))
.expect("ll minus q");
(plus - minus) / (2.0 * eps)
}
"b" => {
let plus = family
.log_likelihood_only(&states_at(q0, b0 + eps, beta_w.clone()))
.expect("ll plus b");
let minus = family
.log_likelihood_only(&states_at(q0, b0 - eps, beta_w.clone()))
.expect("ll minus b");
(plus - minus) / (2.0 * eps)
}
"w0" => {
let mut plus_w = beta_w.clone();
plus_w[0] += eps;
let mut minus_w = beta_w.clone();
minus_w[0] -= eps;
let plus = family
.log_likelihood_only(&states_at(q0, b0, plus_w))
.expect("ll plus w");
let minus = family
.log_likelihood_only(&states_at(q0, b0, minus_w))
.expect("ll minus w");
(plus - minus) / (2.0 * eps)
}
_ => panic!("unknown derivative target"),
};
let eps = 1e-5;
assert!((grad_q - fd("q", eps)).abs() < 2e-4);
assert!((grad_b - fd("b", eps)).abs() < 2e-4);
assert!((grad_w0 - fd("w0", eps)).abs() < 2e-4);
}
#[test]
fn h_only_gradient_matches_loglik_finite_differences() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
..default_test_family()
};
let beta_h = Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
);
let states_at = |q: f64, b: f64, bh: Array1<f64>| {
vec![
ParameterBlockState {
beta: array![q],
eta: Array1::from_elem(z.len(), q),
},
ParameterBlockState {
beta: array![b],
eta: Array1::from_elem(z.len(), b),
},
ParameterBlockState {
beta: bh,
eta: Array1::zeros(z.len()),
},
]
};
let q0 = 0.25;
let b0 = 0.6;
let block_states = states_at(q0, b0, beta_h.clone());
let eval = family.evaluate(&block_states).expect("family evaluation");
let grad_q = match &eval.blockworking_sets[0] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton q block"),
};
let grad_b = match &eval.blockworking_sets[1] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton b block"),
};
let grad_h0 = match &eval.blockworking_sets[2] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton h block"),
};
let fd = |which: &str, eps: f64| match which {
"q" => {
let plus = family
.log_likelihood_only(&states_at(q0 + eps, b0, beta_h.clone()))
.expect("ll plus q");
let minus = family
.log_likelihood_only(&states_at(q0 - eps, b0, beta_h.clone()))
.expect("ll minus q");
(plus - minus) / (2.0 * eps)
}
"b" => {
let plus = family
.log_likelihood_only(&states_at(q0, b0 + eps, beta_h.clone()))
.expect("ll plus b");
let minus = family
.log_likelihood_only(&states_at(q0, b0 - eps, beta_h.clone()))
.expect("ll minus b");
(plus - minus) / (2.0 * eps)
}
"h0" => {
let mut plus_h = beta_h.clone();
plus_h[0] += eps;
let mut minus_h = beta_h.clone();
minus_h[0] -= eps;
let plus = family
.log_likelihood_only(&states_at(q0, b0, plus_h))
.expect("ll plus h");
let minus = family
.log_likelihood_only(&states_at(q0, b0, minus_h))
.expect("ll minus h");
(plus - minus) / (2.0 * eps)
}
_ => panic!("unknown derivative target"),
};
let eps = 1e-5;
assert!((grad_q - fd("q", eps)).abs() < 2e-4);
assert!((grad_b - fd("b", eps)).abs() < 2e-4);
assert!((grad_h0 - fd("h0", eps)).abs() < 2e-4);
}
#[test]
fn flexible_denested_gradient_matches_loglik_finite_differences() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let beta_h = Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
);
let beta_w = Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
);
let states_at = |q: f64, b: f64, bh: Array1<f64>, bw: Array1<f64>| {
vec![
ParameterBlockState {
beta: array![q],
eta: Array1::from_elem(z.len(), q),
},
ParameterBlockState {
beta: array![b],
eta: Array1::from_elem(z.len(), b),
},
ParameterBlockState {
beta: bh,
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: bw,
eta: Array1::zeros(z.len()),
},
]
};
let q0 = 0.25;
let b0 = 0.6;
let block_states = states_at(q0, b0, beta_h.clone(), beta_w.clone());
let eval = family.evaluate(&block_states).expect("family evaluation");
let grad_q = match &eval.blockworking_sets[0] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton q block"),
};
let grad_b = match &eval.blockworking_sets[1] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton b block"),
};
let grad_h0 = match &eval.blockworking_sets[2] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton h block"),
};
let grad_w0 = match &eval.blockworking_sets[3] {
BlockWorkingSet::ExactNewton { gradient, .. } => gradient[0],
_ => panic!("expected exact-newton w block"),
};
let fd = |which: &str, eps: f64| match which {
"q" => {
let plus = family
.log_likelihood_only(&states_at(q0 + eps, b0, beta_h.clone(), beta_w.clone()))
.expect("ll plus q");
let minus = family
.log_likelihood_only(&states_at(q0 - eps, b0, beta_h.clone(), beta_w.clone()))
.expect("ll minus q");
(plus - minus) / (2.0 * eps)
}
"b" => {
let plus = family
.log_likelihood_only(&states_at(q0, b0 + eps, beta_h.clone(), beta_w.clone()))
.expect("ll plus b");
let minus = family
.log_likelihood_only(&states_at(q0, b0 - eps, beta_h.clone(), beta_w.clone()))
.expect("ll minus b");
(plus - minus) / (2.0 * eps)
}
"h0" => {
let mut plus_h = beta_h.clone();
plus_h[0] += eps;
let mut minus_h = beta_h.clone();
minus_h[0] -= eps;
let plus = family
.log_likelihood_only(&states_at(q0, b0, plus_h, beta_w.clone()))
.expect("ll plus h");
let minus = family
.log_likelihood_only(&states_at(q0, b0, minus_h, beta_w.clone()))
.expect("ll minus h");
(plus - minus) / (2.0 * eps)
}
"w0" => {
let mut plus_w = beta_w.clone();
plus_w[0] += eps;
let mut minus_w = beta_w.clone();
minus_w[0] -= eps;
let plus = family
.log_likelihood_only(&states_at(q0, b0, beta_h.clone(), plus_w))
.expect("ll plus w");
let minus = family
.log_likelihood_only(&states_at(q0, b0, beta_h.clone(), minus_w))
.expect("ll minus w");
(plus - minus) / (2.0 * eps)
}
_ => panic!("unknown derivative target"),
};
let eps = 1e-5;
assert!((grad_q - fd("q", eps)).abs() < 2e-4);
assert!((grad_b - fd("b", eps)).abs() < 2e-4);
assert!((grad_h0 - fd("h0", eps)).abs() < 2e-4);
assert!((grad_w0 - fd("w0", eps)).abs() < 2e-4);
}
#[test]
fn flexible_exact_outer_directional_derivatives_are_present_and_finite() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let beta_h = Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
);
let beta_w = Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
);
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: beta_h.clone(),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: beta_w.clone(),
eta: Array1::zeros(z.len()),
},
];
let slices = block_slices(&family);
let total = slices.total;
let mut dir_u = Array1::<f64>::zeros(total);
let mut dir_v = Array1::<f64>::zeros(total);
dir_u[slices.marginal.start] = 0.7;
dir_u[slices.logslope.start] = -0.2;
if let Some(h_range) = slices.h.as_ref() {
dir_u[h_range.start] = 0.1;
if h_range.len() > 1 {
dir_u[h_range.start + 1] = -0.05;
}
}
if let Some(w_range) = slices.w.as_ref() {
dir_u[w_range.start] = 0.08;
}
dir_v[slices.marginal.start] = -0.4;
dir_v[slices.logslope.start] = 0.3;
if let Some(h_range) = slices.h.as_ref() {
dir_v[h_range.start] = -0.03;
}
if let Some(w_range) = slices.w.as_ref() {
dir_v[w_range.start] = 0.06;
if w_range.len() > 1 {
dir_v[w_range.start + 1] = -0.02;
}
}
let third = family
.exact_newton_joint_hessian_directional_derivative(&block_states, &dir_u)
.expect("flex third directional derivative")
.expect("flex third directional derivative matrix");
let fourth = family
.exact_newton_joint_hessiansecond_directional_derivative(&block_states, &dir_u, &dir_v)
.expect("flex fourth directional derivative")
.expect("flex fourth directional derivative matrix");
assert_eq!(third.dim(), (total, total));
assert_eq!(fourth.dim(), (total, total));
assert!(third.iter().all(|value| value.is_finite()));
assert!(fourth.iter().all(|value| value.is_finite()));
let max_abs_third = third
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs()));
let max_abs_fourth = fourth
.iter()
.fold(0.0_f64, |acc, value| acc.max(value.abs()));
assert!(
max_abs_third > 1e-10,
"expected nonzero dual-flex third directional derivative"
);
assert!(
max_abs_fourth > 1e-10,
"expected nonzero dual-flex fourth directional derivative"
);
for i in 0..total {
for j in 0..i {
assert!((third[[i, j]] - third[[j, i]]).abs() < 1e-8);
assert!((fourth[[i, j]] - fourth[[j, i]]).abs() < 1e-8);
}
}
}
#[test]
fn flexible_evaluate_block_diagonals_match_joint_exact_oracle() {
let z = array![-1.1, -0.25, 0.35, 1.2];
let y = array![0.0, 1.0, 0.0, 1.0];
let weights = array![1.0, 0.8, 1.3, 0.7];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-1.8, -0.7, 0.1, 0.9, 1.7];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let marginal_x = array![[1.0, -0.4], [1.0, 0.2], [1.0, 0.7], [1.0, 1.1]];
let logslope_x = array![[1.0, 0.3], [1.0, -0.6], [1.0, 0.5], [1.0, -1.0]];
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
marginal_x.clone(),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
logslope_x.clone(),
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let beta_m = array![0.18, -0.07];
let beta_g = array![0.42, 0.05];
let beta_h = Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.006 * (idx as f64 + 1.0)),
);
let beta_w = Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| -0.004 * (idx as f64 + 1.0)),
);
let block_states = vec![
ParameterBlockState {
beta: beta_m.clone(),
eta: marginal_x.dot(&beta_m),
},
ParameterBlockState {
beta: beta_g.clone(),
eta: logslope_x.dot(&beta_g),
},
ParameterBlockState {
beta: beta_h,
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: beta_w,
eta: Array1::zeros(z.len()),
},
];
let eval = family.evaluate(&block_states).expect("block evaluation");
let joint_hessian = family
.exact_newton_joint_hessian(&block_states)
.expect("joint hessian result")
.expect("dense joint hessian");
let joint_gradient = family
.exact_newton_joint_gradient_evaluation(&block_states, &[])
.expect("joint gradient result")
.expect("joint gradient");
let slices = block_slices(&family);
let ranges = [
slices.marginal.clone(),
slices.logslope.clone(),
slices.h.clone().expect("score-warp block"),
slices.w.clone().expect("link-deviation block"),
];
assert!((eval.log_likelihood - joint_gradient.log_likelihood).abs() < 2e-12);
assert!(
(eval.log_likelihood
- family
.log_likelihood_only(&block_states)
.expect("log likelihood only"))
.abs()
< 2e-12
);
for (block_idx, range) in ranges.iter().enumerate() {
let (gradient, hessian) = match &eval.blockworking_sets[block_idx] {
BlockWorkingSet::ExactNewton { gradient, hessian } => (gradient, hessian),
_ => panic!("expected exact-newton block {block_idx}"),
};
let expected_gradient = joint_gradient.gradient.slice(s![range.clone()]);
for idx in 0..gradient.len() {
assert!(
(gradient[idx] - expected_gradient[idx]).abs() < 2e-10,
"gradient mismatch block {block_idx} idx {idx}: got {} expected {}",
gradient[idx],
expected_gradient[idx]
);
}
let dense_hessian = match hessian {
SymmetricMatrix::Dense(h) => h,
_ => panic!("expected dense hessian block {block_idx}"),
};
let expected_hessian = joint_hessian.slice(s![range.clone(), range.clone()]);
for i in 0..dense_hessian.nrows() {
for j in 0..dense_hessian.ncols() {
assert!(
(dense_hessian[[i, j]] - expected_hessian[[i, j]]).abs() < 2e-9,
"hessian mismatch block {block_idx} ({i},{j}): got {} expected {}",
dense_hessian[[i, j]],
expected_hessian[[i, j]]
);
}
}
}
}
#[test]
fn latent_z_normalization_accepts_finite_sample_gaussian_scores() {
let z = array![
-0.85, -0.12, 0.31, 1.04, -1.21, 0.56, 0.77, -0.44, 1.33, -0.09, 0.28, -0.67
];
let weights = Array1::from_elem(12, 1.0);
let (standardized, normalization) = standardize_latent_z_with_policy(
&z,
&weights,
"bernoulli-marginal-slope",
&LatentZPolicy::exploratory_fit_weighted(),
)
.expect("normalize z");
let replayed = normalization
.apply(&z, "bernoulli-marginal-slope replay")
.expect("replay normalized z");
let mean = standardized.sum() / standardized.len() as f64;
let var = standardized.iter().map(|v| v * v).sum::<f64>() / standardized.len() as f64;
assert_eq!(replayed, standardized);
assert!(mean.abs() < 1e-12);
assert!((var.sqrt() - 1.0).abs() < 1e-12);
}
#[test]
fn latent_z_normalization_rejects_extreme_non_gaussian_scores() {
let z = array![0.0, 0.0, 0.0, 0.0, 10.0, -10.0];
let weights = Array1::from_elem(6, 1.0);
let strict_policy = LatentZPolicy {
check_mode: LatentZCheckMode::Strict,
..LatentZPolicy::default()
};
let err = standardize_latent_z_with_policy(
&z,
&weights,
"bernoulli-marginal-slope",
&strict_policy,
)
.expect_err("expected non-gaussian rejection");
assert!(err.contains("approximately latent N(0,1)"));
}
#[test]
fn auto_latent_measure_uses_rank_int_calibration_for_bad_normal_diagnostics() {
let z = array![0.0, 0.0, 0.0, 0.0, 10.0, -10.0];
let weights = Array1::from_elem(6, 1.0);
let policy = LatentZPolicy {
check_mode: LatentZCheckMode::Off,
normalization: LatentZNormalizationMode::None,
latent_measure: LatentMeasureSpec::Auto { grid_size: 5 },
..LatentZPolicy::default()
};
let (measure, calibration) =
build_latent_measure_with_geometry(&z, &weights, &policy).expect("auto latent measure");
assert!(
matches!(measure, LatentMeasureKind::StandardNormal),
"bad-normal latent z must route through rank-INT to the standard-normal kernel"
);
match calibration {
LatentMeasureCalibration::RankInverseNormal(cal) => {
assert!(
!cal.sorted_z.is_empty(),
"rank-INT knot table must be non-empty"
);
assert_eq!(cal.sorted_z.len(), cal.weighted_cdf.len());
for w in cal.sorted_z.windows(2) {
assert!(w[0] < w[1], "sorted_z must be strictly increasing");
}
for w in cal.weighted_cdf.windows(2) {
assert!(
w[0] <= w[1],
"weighted_cdf must be non-decreasing (got {} -> {})",
w[0],
w[1]
);
}
assert!(
cal.post_mean.abs() < 0.5,
"rank-INT post-mean too far from 0: {}",
cal.post_mean
);
assert!(
cal.post_sd > 0.0 && cal.post_sd.is_finite(),
"rank-INT post-sd must be positive finite, got {}",
cal.post_sd
);
}
LatentMeasureCalibration::None => {
panic!("bad-normal latent z must produce a RankInverseNormal calibration")
}
}
}
#[test]
fn empirical_intercept_calibrates_marginal_probability() {
let nodes = vec![-2.0, -0.25, 0.5, 3.0];
let weights = vec![0.2, 0.3, 0.1, 0.4];
let target_q = -0.35;
let target_mu = normal_cdf(target_q);
let slope = 0.8;
let scale = 0.9;
let intercept = empirical_intercept_from_marginal(
target_mu, target_q, slope, scale, &nodes, &weights, None,
)
.expect("empirical intercept");
let calibrated = nodes
.iter()
.zip(weights.iter())
.map(|(&node, &weight)| weight * normal_cdf(intercept + scale * slope * node))
.sum::<f64>();
assert!((calibrated - target_mu).abs() <= 1e-10);
}
#[test]
fn skewed_rigid_empirical_grid_calibrates_marginal_probability() {
let z = array![-1.15, -1.05, -0.95, -0.8, -0.65, -0.45, 0.1, 0.9, 2.4, 4.7];
let weights = array![1.0, 1.0, 1.0, 0.9, 0.9, 0.8, 0.5, 0.35, 0.2, 0.1];
let grid = build_empirical_z_grid(&z, &weights, 7, "test skewed grid").expect("grid");
let target_q = 0.25;
let target_mu = normal_cdf(target_q);
let slope = 1.35;
let scale = 0.82;
let intercept = empirical_intercept_from_marginal(
target_mu,
target_q,
slope,
scale,
&grid.nodes,
&grid.weights,
None,
)
.expect("empirical intercept");
let calibrated = grid
.nodes
.iter()
.zip(grid.weights.iter())
.map(|(&node, &weight)| weight * normal_cdf(intercept + scale * slope * node))
.sum::<f64>();
assert!((calibrated - target_mu).abs() <= 1e-10);
}
#[test]
fn gaussian_rigid_intercept_miscalibrates_skewed_empirical_law() {
let nodes = vec![-0.95, -0.7, -0.45, -0.2, 0.4, 1.3, 3.1];
let weights = vec![0.28, 0.22, 0.17, 0.13, 0.1, 0.07, 0.03];
let target_q = -0.15;
let target_mu = normal_cdf(target_q);
let slope = 1.4;
let scale = 0.9;
let gaussian_intercept = rigid_intercept_from_marginal(target_q, slope, scale);
let gaussian_mu = nodes
.iter()
.zip(weights.iter())
.map(|(&node, &weight)| weight * normal_cdf(gaussian_intercept + scale * slope * node))
.sum::<f64>();
assert!(
(gaussian_mu - target_mu).abs() > 1e-3,
"skewed empirical law should not be calibrated by Gaussian identity"
);
let empirical_intercept = empirical_intercept_from_marginal(
target_mu, target_q, slope, scale, &nodes, &weights, None,
)
.expect("empirical intercept");
let empirical_mu = nodes
.iter()
.zip(weights.iter())
.map(|(&node, &weight)| weight * normal_cdf(empirical_intercept + scale * slope * node))
.sum::<f64>();
assert!((empirical_mu - target_mu).abs() <= 1e-10);
}
#[test]
fn empirical_intercept_recovers_from_deep_tail_warm_start() {
let nodes = vec![-2.0, -0.25, 0.5, 3.0];
let weights = vec![0.2, 0.3, 0.1, 0.4];
let target_q = -0.35;
let target_mu = normal_cdf(target_q);
let slope = 0.8;
let scale = 0.9;
let stale_warm_start = Some(-100.0_f64);
let intercept = empirical_intercept_from_marginal(
target_mu,
target_q,
slope,
scale,
&nodes,
&weights,
stale_warm_start,
)
.expect("empirical intercept must recover from deep-tail warm start");
let calibrated = nodes
.iter()
.zip(weights.iter())
.map(|(&node, &weight)| weight * normal_cdf(intercept + scale * slope * node))
.sum::<f64>();
assert!(
(calibrated - target_mu).abs() <= 1e-10,
"calibrated mu={calibrated} should match target_mu={target_mu} from any seed"
);
}
#[test]
fn empirical_intercept_recovers_from_far_right_warm_start() {
let nodes = vec![-2.0, -0.25, 0.5, 3.0];
let weights = vec![0.2, 0.3, 0.1, 0.4];
let target_q = -0.35;
let target_mu = normal_cdf(target_q);
let slope = 0.8;
let scale = 0.9;
let stale_warm_start = Some(50.0_f64);
let intercept = empirical_intercept_from_marginal(
target_mu,
target_q,
slope,
scale,
&nodes,
&weights,
stale_warm_start,
)
.expect("empirical intercept must recover from far-right warm start");
let calibrated = nodes
.iter()
.zip(weights.iter())
.map(|(&node, &weight)| weight * normal_cdf(intercept + scale * slope * node))
.sum::<f64>();
assert!(
(calibrated - target_mu).abs() <= 1e-10,
"calibrated mu={calibrated} should match target_mu={target_mu} from any seed"
);
}
#[test]
fn auto_latent_measure_preserves_standard_normal_fast_path() {
let n = 2001usize;
let z = Array1::from_iter((0..n).map(|idx| {
let p = (idx as f64 + 0.5) / n as f64;
standard_normal_quantile(p).expect("normal quantile")
}));
let weights = Array1::ones(n);
let policy = LatentZPolicy {
latent_measure: LatentMeasureSpec::Auto { grid_size: 17 },
..LatentZPolicy::default()
};
let (measure, calibration) =
build_latent_measure_with_geometry(&z, &weights, &policy).expect("measure");
assert!(matches!(measure, LatentMeasureKind::StandardNormal));
assert!(
matches!(calibration, LatentMeasureCalibration::None),
"well-conditioned standard-normal z must skip rank-INT calibration"
);
let slope = 0.7;
let scale = 0.85;
let target_q = 0.4;
assert_eq!(
rigid_intercept_from_marginal(target_q, slope, scale),
target_q * (1.0 + (scale * slope).powi(2)).sqrt()
);
}
#[test]
fn flexible_family_exposes_exact_newton_workspaces() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let score_prepared = build_score_warp_deviation_block_from_seed(
&z,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("score warp block");
let link_seed = array![-2.0, -0.5, 0.0, 0.5, 2.0];
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 4,
..DeviationBlockConfig::default()
},
)
.expect("link block");
let family =
BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..score_prepared.block.design.ncols()).map(|idx| 0.015 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: Array1::from_iter(
(0..link_prepared.block.design.ncols()).map(|idx| 0.01 * (idx as f64 + 1.0)),
),
eta: Array1::zeros(z.len()),
},
];
let specs = vec![
dummy_blockspec(1, z.len()),
dummy_blockspec(1, z.len()),
dummy_blockspec(score_prepared.block.design.ncols(), z.len()),
dummy_blockspec(link_prepared.block.design.ncols(), z.len()),
];
let derivative_blocks = vec![Vec::new(), Vec::new(), Vec::new(), Vec::new()];
assert!(
family
.exact_newton_joint_hessian_workspace(&block_states, &specs)
.expect("flex hessian workspace")
.is_some()
);
assert!(
family
.exact_newton_joint_psi_workspace(&block_states, &specs, &derivative_blocks)
.expect("flex psi workspace")
.is_some()
);
}
#[test]
fn sigma_exact_joint_psi_terms_returns_analytic_terms() {
let z = array![-0.8, 0.2, 1.1];
let y = array![0.0, 1.0, 1.0];
let weights = array![1.0, 0.7, 1.3];
let sigma = 0.7;
let make_family =
|sigma: f64| BernoulliMarginalSlopeFamily {
y: Arc::new(y.clone()),
weights: Arc::new(weights.clone()),
z: Arc::new(z.clone()),
gaussian_frailty_sd: Some(sigma),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
array![[1.0], [1.0], [1.0]],
)),
..default_test_family()
};
let family = make_family(sigma);
let block_states = vec![
ParameterBlockState {
beta: array![0.25],
eta: Array1::from_elem(z.len(), 0.25),
},
ParameterBlockState {
beta: array![0.6],
eta: Array1::from_elem(z.len(), 0.6),
},
];
let specs = vec![dummy_blockspec(1, z.len()), dummy_blockspec(1, z.len())];
let terms = family
.sigma_exact_joint_psi_terms(&block_states, &specs)
.expect("analytic sigma psi terms")
.expect("sigma terms present");
assert!(terms.objective_psi.is_finite());
assert_eq!(terms.score_psi.len(), 2);
assert!(terms.score_psi.iter().all(|value| value.is_finite()));
assert_eq!(
terms
.hessian_psi_operator
.as_ref()
.expect("sigma Hessian operator")
.to_dense()
.dim(),
(2, 2)
);
let second = family
.sigma_exact_joint_psisecond_order_terms(&block_states)
.expect("analytic second sigma terms")
.expect("second sigma terms present");
assert!(second.objective_psi_psi.is_finite());
assert_eq!(second.score_psi_psi.len(), 2);
let drift = family
.sigma_exact_joint_psihessian_directional_derivative(&block_states, &array![0.1, -0.2])
.expect("analytic sigma Hessian directional derivative")
.expect("sigma drift present");
assert_eq!(drift.dim(), (2, 2));
assert!(drift.iter().all(|value| value.is_finite()));
let tau = sigma.ln();
let eps = 1e-5;
let ll_plus = make_family((tau + eps).exp())
.log_likelihood_only(&block_states)
.expect("ll plus sigma");
let ll_minus = make_family((tau - eps).exp())
.log_likelihood_only(&block_states)
.expect("ll minus sigma");
let objective_fd = -(ll_plus - ll_minus) / (2.0 * eps);
assert!((terms.objective_psi - objective_fd).abs() < 1e-5);
}
fn make_block_psi_test_family(n: usize) -> BernoulliMarginalSlopeFamily {
let y: Array1<f64> =
Array1::from_iter((0..n).map(|i| if (i * 31 + 7) % 5 >= 3 { 1.0 } else { 0.0 }));
let weights: Array1<f64> =
Array1::from_iter((0..n).map(|i| 0.5 + ((i * 13 + 4) % 7) as f64 * 0.1));
let z: Array1<f64> = Array1::from_iter(
(0..n).map(|i| -1.5 + 3.0 * (((i * 17 + 5) % n) as f64 + 0.5) / (n as f64)),
);
let marginal_design = Array2::from_shape_fn((n, 1), |(i, _)| {
0.3 + 0.4 * (((i * 29 + 11) % n) as f64) / (n as f64)
});
let logslope_design = Array2::from_shape_fn((n, 1), |(i, _)| {
0.2 + 0.5 * (((i * 37 + 9) % n) as f64) / (n as f64)
});
BernoulliMarginalSlopeFamily {
y: Arc::new(y),
weights: Arc::new(weights),
z: Arc::new(z),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
marginal_design,
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
logslope_design,
)),
..default_test_family()
}
}
fn block_psi_test_block_states(
family: &BernoulliMarginalSlopeFamily,
m_beta: f64,
g_beta: f64,
) -> Vec<ParameterBlockState> {
let m_design = family.marginal_design.to_dense().to_owned();
let g_design = family.logslope_design.to_dense().to_owned();
let m_eta = m_design.dot(&array![m_beta]);
let g_eta = g_design.dot(&array![g_beta]);
vec![
ParameterBlockState {
beta: array![m_beta],
eta: m_eta,
},
ParameterBlockState {
beta: array![g_beta],
eta: g_eta,
},
]
}
fn block_psi_test_marginal_derivative_blocks(
n: usize,
) -> Vec<Vec<crate::custom_family::CustomFamilyBlockPsiDerivative>> {
let x_psi = Array2::from_shape_fn((n, 1), |(i, _)| {
0.4 + 0.3 * (((i * 41 + 13) % n) as f64) / (n as f64)
});
vec![
vec![crate::custom_family::CustomFamilyBlockPsiDerivative::new(
None,
x_psi,
Array2::zeros((1, 1)),
None,
None,
None,
None,
)],
Vec::new(),
]
}
fn block_psi_test_dual_derivative_blocks(
n: usize,
) -> Vec<Vec<crate::custom_family::CustomFamilyBlockPsiDerivative>> {
let x_psi_m = Array2::from_shape_fn((n, 1), |(i, _)| {
0.4 + 0.3 * (((i * 41 + 13) % n) as f64) / (n as f64)
});
let x_psi_g = Array2::from_shape_fn((n, 1), |(i, _)| {
0.2 + 0.5 * (((i * 43 + 17) % n) as f64) / (n as f64)
});
vec![
vec![crate::custom_family::CustomFamilyBlockPsiDerivative::new(
None,
x_psi_m,
Array2::zeros((1, 1)),
None,
None,
None,
None,
)],
vec![crate::custom_family::CustomFamilyBlockPsiDerivative::new(
None,
x_psi_g,
Array2::zeros((1, 1)),
None,
None,
None,
None,
)],
]
}
#[test]
fn bernoulli_psi_terms_from_cache_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let derivative_blocks = block_psi_test_marginal_derivative_blocks(n);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let baseline = family
.exact_newton_joint_psi_terms_from_cache(&states, &derivative_blocks, 0, &cache)
.expect("baseline psi terms")
.expect("some");
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full = family
.exact_newton_joint_psi_terms_from_cache_with_options(
&states,
&derivative_blocks,
0,
&cache,
&opts_full,
)
.expect("with full")
.expect("some");
let obj_rel = ((with_full.objective_psi - baseline.objective_psi)
/ baseline.objective_psi.abs().max(1.0))
.abs();
assert!(obj_rel < 1e-12, "objective_psi rel {}", obj_rel);
let score_rel = rel_diff_array1(&with_full.score_psi, &baseline.score_psi);
assert!(score_rel < 1e-12, "score_psi rel {}", score_rel);
}
#[test]
fn bernoulli_psi_terms_from_cache_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let derivative_blocks = block_psi_test_marginal_derivative_blocks(n);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.exact_newton_joint_psi_terms_from_cache_with_options(
&states,
&derivative_blocks,
0,
&cache,
&opts_half,
)
.expect("scaled")
.expect("some");
let mut opts_raw = BlockwiseFitOptions::default();
opts_raw.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::with_uniform_weight(
even_mask, m, 0, 1.0,
)));
let raw = family
.exact_newton_joint_psi_terms_from_cache_with_options(
&states,
&derivative_blocks,
0,
&cache,
&opts_raw,
)
.expect("raw")
.expect("some");
let factor = n as f64 / m as f64;
let exp_obj = factor * raw.objective_psi;
let obj_rel = ((scaled.objective_psi - exp_obj) / exp_obj.abs().max(1.0)).abs();
assert!(obj_rel < 1e-12, "objective_psi rel {}", obj_rel);
let exp_score = &raw.score_psi * factor;
let score_rel = rel_diff_array1(&scaled.score_psi, &exp_score);
assert!(score_rel < 1e-12, "score_psi rel {}", score_rel);
}
#[test]
fn bernoulli_psi_second_order_terms_from_cache_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let derivative_blocks = block_psi_test_dual_derivative_blocks(n);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let baseline = family
.exact_newton_joint_psisecond_order_terms_from_cache(
&states,
&derivative_blocks,
0,
1,
&cache,
)
.expect("baseline psi second-order")
.expect("some");
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full = family
.exact_newton_joint_psisecond_order_terms_from_cache_with_options(
&states,
&derivative_blocks,
0,
1,
&cache,
&opts_full,
)
.expect("with full")
.expect("some");
let obj_rel = ((with_full.objective_psi_psi - baseline.objective_psi_psi)
/ baseline.objective_psi_psi.abs().max(1.0))
.abs();
assert!(obj_rel < 1e-12, "objective rel {}", obj_rel);
let score_rel = rel_diff_array1(&with_full.score_psi_psi, &baseline.score_psi_psi);
assert!(score_rel < 1e-12, "score rel {}", score_rel);
}
#[test]
fn bernoulli_psi_second_order_terms_from_cache_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let derivative_blocks = block_psi_test_dual_derivative_blocks(n);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.exact_newton_joint_psisecond_order_terms_from_cache_with_options(
&states,
&derivative_blocks,
0,
1,
&cache,
&opts_half,
)
.expect("scaled")
.expect("some");
let mut opts_raw = BlockwiseFitOptions::default();
opts_raw.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::with_uniform_weight(
even_mask, m, 0, 1.0,
)));
let raw = family
.exact_newton_joint_psisecond_order_terms_from_cache_with_options(
&states,
&derivative_blocks,
0,
1,
&cache,
&opts_raw,
)
.expect("raw")
.expect("some");
let factor = n as f64 / m as f64;
let exp_obj = factor * raw.objective_psi_psi;
let obj_rel = ((scaled.objective_psi_psi - exp_obj) / exp_obj.abs().max(1.0)).abs();
assert!(obj_rel < 1e-12, "objective rel {}", obj_rel);
let exp_score = &raw.score_psi_psi * factor;
let score_rel = rel_diff_array1(&scaled.score_psi_psi, &exp_score);
assert!(score_rel < 1e-12, "score rel {}", score_rel);
}
#[test]
fn bernoulli_psihessian_directional_derivative_from_cache_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let derivative_blocks = block_psi_test_marginal_derivative_blocks(n);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let mut d_beta_flat = Array1::<f64>::zeros(slices.total);
d_beta_flat[slices.marginal.start] = 0.05;
d_beta_flat[slices.logslope.start] = -0.04;
let baseline = family
.exact_newton_joint_psihessian_directional_derivative_from_cache(
&states,
&derivative_blocks,
0,
&d_beta_flat,
&cache,
)
.expect("baseline")
.expect("some");
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full = family
.exact_newton_joint_psihessian_directional_derivative_from_cache_with_options(
&states,
&derivative_blocks,
0,
&d_beta_flat,
&cache,
&opts_full,
)
.expect("with full")
.expect("some");
let rel = rel_diff_array2(&with_full, &baseline);
assert!(rel < 1e-12, "drift rel {}", rel);
}
#[test]
fn bernoulli_psihessian_directional_derivative_from_cache_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let derivative_blocks = block_psi_test_marginal_derivative_blocks(n);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let mut d_beta_flat = Array1::<f64>::zeros(slices.total);
d_beta_flat[slices.marginal.start] = 0.05;
d_beta_flat[slices.logslope.start] = -0.04;
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.exact_newton_joint_psihessian_directional_derivative_from_cache_with_options(
&states,
&derivative_blocks,
0,
&d_beta_flat,
&cache,
&opts_half,
)
.expect("scaled")
.expect("some");
let mut opts_raw = BlockwiseFitOptions::default();
opts_raw.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::with_uniform_weight(
even_mask, m, 0, 1.0,
)));
let raw = family
.exact_newton_joint_psihessian_directional_derivative_from_cache_with_options(
&states,
&derivative_blocks,
0,
&d_beta_flat,
&cache,
&opts_raw,
)
.expect("raw")
.expect("some");
let factor = n as f64 / m as f64;
let exp = &raw * factor;
let rel = rel_diff_array2(&scaled, &exp);
assert!(rel < 1e-12, "drift rel {}", rel);
}
#[test]
fn bernoulli_psihessian_operator_from_cache_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let derivative_blocks = block_psi_test_marginal_derivative_blocks(n);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let mut d_beta_flat = Array1::<f64>::zeros(slices.total);
d_beta_flat[slices.marginal.start] = 0.05;
d_beta_flat[slices.logslope.start] = -0.04;
let baseline = family
.exact_newton_joint_psihessian_directional_derivative_operator_from_cache_with_options(
&states,
&derivative_blocks,
0,
&d_beta_flat,
&cache,
&BlockwiseFitOptions::default(),
)
.expect("baseline operator")
.expect("some");
let baseline_dense = baseline.to_dense();
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full = family
.exact_newton_joint_psihessian_directional_derivative_operator_from_cache_with_options(
&states,
&derivative_blocks,
0,
&d_beta_flat,
&cache,
&opts_full,
)
.expect("with full")
.expect("some");
let with_full_dense = with_full.to_dense();
let rel = rel_diff_array2(&with_full_dense, &baseline_dense);
assert!(rel < 1e-12, "operator drift rel {}", rel);
}
#[test]
fn bernoulli_psihessian_operator_from_cache_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let derivative_blocks = block_psi_test_marginal_derivative_blocks(n);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let mut d_beta_flat = Array1::<f64>::zeros(slices.total);
d_beta_flat[slices.marginal.start] = 0.05;
d_beta_flat[slices.logslope.start] = -0.04;
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.exact_newton_joint_psihessian_directional_derivative_operator_from_cache_with_options(
&states,
&derivative_blocks,
0,
&d_beta_flat,
&cache,
&opts_half,
)
.expect("scaled")
.expect("some");
let scaled_dense = scaled.to_dense();
let mut opts_raw = BlockwiseFitOptions::default();
opts_raw.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::with_uniform_weight(
even_mask, m, 0, 1.0,
)));
let raw = family
.exact_newton_joint_psihessian_directional_derivative_operator_from_cache_with_options(
&states,
&derivative_blocks,
0,
&d_beta_flat,
&cache,
&opts_raw,
)
.expect("raw")
.expect("some");
let raw_dense = raw.to_dense();
let factor = n as f64 / m as f64;
let exp = &raw_dense * factor;
let rel = rel_diff_array2(&scaled_dense, &exp);
assert!(rel < 1e-12, "operator drift rel {}", rel);
}
#[test]
fn bernoulli_jointhessian_directional_derivative_from_cache_subsample_full_equals_unsampled() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let mut d_beta_flat = Array1::<f64>::zeros(slices.total);
d_beta_flat[slices.marginal.start] = 0.05;
d_beta_flat[slices.logslope.start] = -0.04;
let baseline = family
.exact_newton_joint_hessian_directional_derivative_from_cache(
&states,
&d_beta_flat,
&cache,
)
.expect("baseline")
.expect("some");
let mut opts_full = BlockwiseFitOptions::default();
opts_full.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).collect(),
n,
0xDEADBEEF,
)));
let with_full = family
.exact_newton_joint_hessian_directional_derivative_from_cache_with_options(
&states,
&d_beta_flat,
&cache,
&opts_full,
)
.expect("with full")
.expect("some");
let rel = rel_diff_array2(&with_full, &baseline);
assert!(rel < 1e-12, "joint Hessian dH drift rel {}", rel);
}
#[test]
fn bernoulli_jointhessian_batched_directional_operators_match_single_direction_path() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let directions: Vec<Array1<f64>> = (0..3)
.map(|rep| {
let mut d = Array1::<f64>::zeros(slices.total);
d[slices.marginal.start] = 0.03 * (rep as f64 + 1.0);
d[slices.logslope.start] = -0.02 * (rep as f64 + 1.0);
d
})
.collect();
let mut opts = BlockwiseFitOptions::default();
opts.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..n).step_by(2).collect(),
n,
0xB47C,
)));
let batched = family
.exact_newton_joint_hessian_directional_derivative_operators_from_cache_with_options(
&states,
&directions,
&cache,
&opts,
)
.expect("batched operators");
assert_eq!(batched.len(), directions.len());
for (idx, direction) in directions.iter().enumerate() {
let single = family
.exact_newton_joint_hessian_directional_derivative_operator_from_cache_with_options(
&states, direction, &cache, &opts,
)
.expect("single operator")
.expect("single operator some")
.to_dense();
let batched_dense = batched[idx]
.as_ref()
.expect("batched operator some")
.to_dense();
let rel = rel_diff_array2(&batched_dense, &single);
assert!(rel < 1e-12, "batched operator {idx} drift rel {rel}");
}
}
fn make_flex_hvp_cache_test_family(
n: usize,
) -> (BernoulliMarginalSlopeFamily, Vec<ParameterBlockState>) {
let score_seed = Array1::linspace(-2.0, 2.0, n.max(6));
let link_seed = Array1::linspace(-1.8, 1.8, n.max(6));
let score_prepared = build_score_warp_deviation_block_from_seed(
&score_seed,
&DeviationBlockConfig {
num_internal_knots: 3,
..DeviationBlockConfig::default()
},
)
.expect("build score warp block");
let link_prepared = build_test_link_deviation_block_from_seed(
&link_seed,
&DeviationBlockConfig {
num_internal_knots: 3,
..DeviationBlockConfig::default()
},
)
.expect("build link deviation block");
let y: Array1<f64> =
Array1::from_iter((0..n).map(|i| if (i * 17 + 3) % 7 >= 4 { 1.0 } else { 0.0 }));
let weights: Array1<f64> =
Array1::from_iter((0..n).map(|i| 0.75 + ((i * 11 + 5) % 5) as f64 * 0.05));
let z: Array1<f64> =
Array1::from_iter((0..n).map(|i| -1.7 + 3.4 * (i as f64 + 0.5) / n as f64));
let marginal_x = Array2::from_shape_fn((n, 2), |(i, j)| {
if j == 0 {
1.0
} else {
-0.4 + 0.8 * ((i * 19 + 7) % n) as f64 / n as f64
}
});
let logslope_x = Array2::from_shape_fn((n, 2), |(i, j)| {
if j == 0 {
1.0
} else {
0.3 - 0.6 * ((i * 23 + 11) % n) as f64 / n as f64
}
});
let family = BernoulliMarginalSlopeFamily {
y: Arc::new(y),
weights: Arc::new(weights),
z: Arc::new(z.clone()),
gaussian_frailty_sd: Some(0.15),
marginal_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
marginal_x.clone(),
)),
logslope_design: DesignMatrix::Dense(crate::matrix::DenseDesignMatrix::from(
logslope_x.clone(),
)),
score_warp: Some(score_prepared.runtime.clone()),
link_dev: Some(link_prepared.runtime.clone()),
..default_test_family()
};
let beta_m = array![0.12, -0.04];
let beta_g = array![0.35, 0.03];
let beta_h = Array1::from_iter(
(0..score_prepared.runtime.basis_dim()).map(|idx| 0.0015 * (idx as f64 + 1.0)),
);
let beta_w = Array1::from_iter(
(0..link_prepared.runtime.basis_dim()).map(|idx| -0.001 * (idx as f64 + 1.0)),
);
let states = vec![
ParameterBlockState {
beta: beta_m.clone(),
eta: marginal_x.dot(&beta_m),
},
ParameterBlockState {
beta: beta_g.clone(),
eta: logslope_x.dot(&beta_g),
},
ParameterBlockState {
beta: beta_h,
eta: Array1::zeros(z.len()),
},
ParameterBlockState {
beta: beta_w,
eta: Array1::zeros(z.len()),
},
];
(family, states)
}
#[test]
fn bernoulli_flex_paired_subsample_ll_delta_sign_matches_full_ll() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let (family, old_states) = make_flex_hvp_cache_test_family(96);
let mut trial_states = old_states.clone();
trial_states[0].beta[0] += 0.015;
trial_states[0].eta += 0.015;
trial_states[1].beta[1] -= 0.01;
let logslope_col =
Array1::from_iter((0..96).map(|i| 0.3 - 0.6 * ((i * 23 + 11) % 96) as f64 / 96.0));
trial_states[1].eta.scaled_add(-0.01, &logslope_col);
let full_old = family
.log_likelihood_only(&old_states)
.expect("full old ll");
let full_trial = family
.log_likelihood_only(&trial_states)
.expect("full trial ll");
let mut opts = BlockwiseFitOptions::default();
opts.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
(0..96).step_by(2).collect(),
96,
0x5EED5EED,
)));
let sub_old = family
.log_likelihood_only_with_options(&old_states, &opts)
.expect("subsample old ll");
let sub_trial = family
.log_likelihood_only_with_options(&trial_states, &opts)
.expect("subsample trial ll");
let full_delta = full_trial - full_old;
let sub_delta = sub_trial - sub_old;
assert!(
full_delta.abs() > 1e-8,
"synthetic beta-pair should produce a non-degenerate full-LL delta: {full_delta}"
);
assert_eq!(
full_delta.is_sign_positive(),
sub_delta.is_sign_positive(),
"paired subsample LL delta sign ({sub_delta}) should match full LL delta sign ({full_delta})"
);
}
#[test]
fn bernoulli_flex_row_primary_hessian_cache_policy_materializes_aou_shape() {
let plan = decide_row_primary_hessian_cache(
195_780,
20,
BMS_ROW_PRIMARY_HESSIAN_EXPECTED_REUSE_PASSES,
16 * 1024 * 1024 * 1024,
0,
);
assert_eq!(plan.bytes, 626_496_000);
assert!(plan.materialize);
assert_eq!(
plan.reason,
RowPrimaryHessianCacheReason::ReuseAmortizesBuild
);
}
#[test]
fn bernoulli_flex_row_primary_hessian_cache_policy_streams_when_single_cache_exceeds_ram() {
let plan = decide_row_primary_hessian_cache(
195_780,
20,
BMS_ROW_PRIMARY_HESSIAN_EXPECTED_REUSE_PASSES,
2 * 1024 * 1024 * 1024,
0,
);
assert!(!plan.materialize);
assert_eq!(
plan.reason,
RowPrimaryHessianCacheReason::SingleCacheExceedsRamFraction
);
}
#[test]
fn bernoulli_flex_row_primary_hessian_cache_policy_streams_when_global_pin_exhausted() {
let plan = decide_row_primary_hessian_cache(
195_780,
20,
BMS_ROW_PRIMARY_HESSIAN_EXPECTED_REUSE_PASSES,
16 * 1024 * 1024 * 1024,
(15 * 1024 * 1024 * 1024) / 2,
);
assert!(!plan.materialize);
assert_eq!(
plan.reason,
RowPrimaryHessianCacheReason::GlobalPinExceedsRamFraction
);
}
#[test]
fn bernoulli_flex_row_primary_hessian_cache_policy_streams_low_reuse() {
let plan = decide_row_primary_hessian_cache(100, 4, 1, 16 * 1024 * 1024 * 1024, 0);
assert!(!plan.materialize);
assert_eq!(plan.reason, RowPrimaryHessianCacheReason::ReuseTooLow);
}
#[test]
fn bernoulli_flex_hvp_cache_matches_uncached_path_small_case() {
let (family, states) = make_flex_hvp_cache_test_family(12);
let mut cached = family
.build_exact_eval_cache(&states)
.expect("cached exact eval cache");
cached.row_primary_hessians = family
.build_row_primary_hessian_cache(&states, &cached)
.expect("row Hessian cache");
let uncached = BernoulliMarginalSlopeExactEvalCache {
slices: cached.slices.clone(),
primary: cached.primary.clone(),
row_contexts: cached.row_contexts.clone(),
row_cell_moments: None,
row_cell_moments_d15: crate::resource::RayonSafeOnce::new(),
row_cell_moments_d21: crate::resource::RayonSafeOnce::new(),
row_primary_hessians: RowPrimaryEvalCache::Empty,
rigid_third_full: crate::resource::RayonSafeOnce::new(),
rigid_fourth_full: crate::resource::RayonSafeOnce::new(),
};
let direction =
Array1::from_iter((0..cached.slices.total).map(|idx| 0.02 * ((idx % 5) as f64 - 2.0)));
let hv_cached = family
.exact_newton_joint_hessian_matvec_from_cache(&direction, &states, &cached)
.expect("cached Hv");
let hv_uncached = family
.exact_newton_joint_hessian_matvec_from_cache(&direction, &states, &uncached)
.expect("uncached Hv");
let rel = rel_diff_array1(&hv_cached, &hv_uncached);
assert!(rel < 5e-11, "cached Hv drift rel {rel}");
let dense_cached = family
.exact_newton_joint_hessian_dense_from_cache(&states, &cached)
.expect("cached dense Hessian");
let dense_uncached = family
.exact_newton_joint_hessian_dense_from_cache(&states, &uncached)
.expect("uncached dense Hessian");
let rel_dense = rel_diff_array2(&dense_cached, &dense_uncached);
assert!(
rel_dense < 5e-11,
"cached dense Hessian drift rel {rel_dense}"
);
let diag_cached = family
.exact_newton_joint_hessian_diagonal_from_cache(&states, &cached)
.expect("cached diag");
let diag_uncached = family
.exact_newton_joint_hessian_diagonal_from_cache(&states, &uncached)
.expect("uncached diag");
let rel_diag = rel_diff_array1(&diag_cached, &diag_uncached);
assert!(rel_diag < 5e-11, "cached diag drift rel {rel_diag}");
}
#[test]
fn bernoulli_flex_hvp_cache_timing_biobank_shape_pattern() {
let (family, states) = make_flex_hvp_cache_test_family(96);
let mut cached = family
.build_exact_eval_cache(&states)
.expect("cached exact eval cache");
cached.row_primary_hessians = family
.build_row_primary_hessian_cache(&states, &cached)
.expect("row Hessian cache");
let uncached = BernoulliMarginalSlopeExactEvalCache {
slices: cached.slices.clone(),
primary: cached.primary.clone(),
row_contexts: cached.row_contexts.clone(),
row_cell_moments: None,
row_cell_moments_d15: crate::resource::RayonSafeOnce::new(),
row_cell_moments_d21: crate::resource::RayonSafeOnce::new(),
row_primary_hessians: RowPrimaryEvalCache::Empty,
rigid_third_full: crate::resource::RayonSafeOnce::new(),
rigid_fourth_full: crate::resource::RayonSafeOnce::new(),
};
let directions: Vec<_> = (0..4)
.map(|rep| {
Array1::from_iter(
(0..cached.slices.total)
.map(|idx| 0.01 * (((idx * 13 + rep * 7) % 11) as f64 - 5.0)),
)
})
.collect();
let start_uncached = std::time::Instant::now();
for direction in &directions {
family
.exact_newton_joint_hessian_matvec_from_cache(direction, &states, &uncached)
.expect("uncached Hv");
}
let uncached_elapsed = start_uncached.elapsed();
let start_cached = std::time::Instant::now();
for direction in &directions {
family
.exact_newton_joint_hessian_matvec_from_cache(direction, &states, &cached)
.expect("cached Hv");
}
let cached_elapsed = start_cached.elapsed();
eprintln!("flex Hv cache timing: uncached={uncached_elapsed:?} cached={cached_elapsed:?}");
assert!(
cached_elapsed < uncached_elapsed,
"expected cached Hv loop to beat uncached: cached={cached_elapsed:?} uncached={uncached_elapsed:?}"
);
}
#[test]
fn bernoulli_jointhessian_directional_operator_matches_dense_small_case() {
let n = 17usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let mut d_beta_flat = Array1::<f64>::zeros(slices.total);
d_beta_flat[slices.marginal.start] = 0.05;
d_beta_flat[slices.logslope.start] = -0.04;
let dense = family
.exact_newton_joint_hessian_directional_derivative_from_cache_with_options(
&states,
&d_beta_flat,
&cache,
&BlockwiseFitOptions::default(),
)
.expect("dense dH")
.expect("dH present");
let operator = family
.exact_newton_joint_hessian_directional_derivative_operator_from_cache_with_options(
&states,
&d_beta_flat,
&cache,
&BlockwiseFitOptions::default(),
)
.expect("operator dH")
.expect("dH operator present");
let rel = rel_diff_array2(&operator.to_dense(), &dense);
assert!(rel < 1e-12, "operator dH rel {}", rel);
}
#[test]
fn bernoulli_jointhessian_second_directional_operator_matches_dense_small_case() {
let n = 17usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let mut d_beta_u = Array1::<f64>::zeros(slices.total);
d_beta_u[slices.marginal.start] = 0.05;
d_beta_u[slices.logslope.start] = -0.04;
let mut d_beta_v = Array1::<f64>::zeros(slices.total);
d_beta_v[slices.marginal.start] = -0.03;
d_beta_v[slices.logslope.start] = 0.02;
let dense = family
.exact_newton_joint_hessiansecond_directional_derivative_from_cache_with_options(
&states,
&d_beta_u,
&d_beta_v,
&cache,
&BlockwiseFitOptions::default(),
)
.expect("dense d2H")
.expect("d2H present");
let operator = family
.exact_newton_joint_hessiansecond_directional_derivative_operator_from_cache_with_options(
&states,
&d_beta_u,
&d_beta_v,
&cache,
&BlockwiseFitOptions::default(),
)
.expect("operator d2H")
.expect("d2H operator present");
let rel = rel_diff_array2(&operator.to_dense(), &dense);
assert!(rel < 1e-12, "operator d2H rel {}", rel);
}
#[test]
fn bernoulli_large_scale_outer_derivatives_keep_analytic_hessian_route() {
let n = 50_001usize;
let family = make_block_psi_test_family(n);
let specs = vec![dummy_blockspec(1, n), dummy_blockspec(1, n)];
let options = BlockwiseFitOptions::default();
let (gradient, hessian) =
crate::custom_family::custom_family_outer_derivatives(&family, &specs, &options);
assert_eq!(
gradient,
crate::solver::outer_strategy::Derivative::Analytic
);
assert_eq!(
hessian,
crate::solver::outer_strategy::DeclaredHessianForm::Either
);
}
#[test]
fn bernoulli_jointhessian_directional_derivative_from_cache_subsample_half_scales_correctly() {
use crate::families::marginal_slope_shared::OuterScoreSubsample;
let n = 200usize;
let family = make_block_psi_test_family(n);
let states = block_psi_test_block_states(&family, 0.15, 0.25);
let cache = family
.build_exact_eval_cache(&states)
.expect("exact eval cache");
let slices = &cache.slices;
let mut d_beta_flat = Array1::<f64>::zeros(slices.total);
d_beta_flat[slices.marginal.start] = 0.05;
d_beta_flat[slices.logslope.start] = -0.04;
let even_mask: Vec<usize> = (0..n).filter(|i| i % 2 == 0).collect();
let m = even_mask.len();
let mut opts_half = BlockwiseFitOptions::default();
opts_half.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::new(
even_mask.clone(),
n,
0xCAFE,
)));
let scaled = family
.exact_newton_joint_hessian_directional_derivative_from_cache_with_options(
&states,
&d_beta_flat,
&cache,
&opts_half,
)
.expect("scaled")
.expect("some");
let mut opts_raw = BlockwiseFitOptions::default();
opts_raw.outer_score_subsample = Some(Arc::new(OuterScoreSubsample::with_uniform_weight(
even_mask, m, 0, 1.0,
)));
let raw = family
.exact_newton_joint_hessian_directional_derivative_from_cache_with_options(
&states,
&d_beta_flat,
&cache,
&opts_raw,
)
.expect("raw")
.expect("some");
let factor = n as f64 / m as f64;
let exp = &raw * factor;
let rel = rel_diff_array2(&scaled, &exp);
assert!(rel < 1e-12, "joint Hessian dH HT rel {}", rel);
}
#[test]
fn auto_outer_subsample_two_phase_converges_to_full_data_optimum() {
let n = 35_000usize;
let family = make_block_psi_test_family(n);
let states: Vec<ParameterBlockState> = Vec::new();
let specs: Vec<ParameterBlockSpec> = Vec::new();
let deriv_blocks: Vec<Vec<crate::custom_family::CustomFamilyBlockPsiDerivative>> =
Vec::new();
let rho_dim = 3usize;
let mut opts = BlockwiseFitOptions::default();
opts.auto_outer_subsample = true;
let counter = || {
family
.auto_subsample_phase_counter
.load(std::sync::atomic::Ordering::SeqCst)
};
let mut distinct_calls = 0usize;
for step in 0..15usize {
let rho_step = Array1::<f64>::from_elem(rho_dim, step as f64 * 0.1);
family
.batched_outer_gradient_terms(
&states,
&specs,
&deriv_blocks,
&rho_step,
&opts,
None,
)
.expect("guard ok");
distinct_calls += 1;
assert_eq!(
counter(),
distinct_calls,
"distinct-ρ call {step}: counter should equal number of distinct ρ"
);
if step == 3 || step == 7 {
let rho_retry = Array1::<f64>::from_elem(rho_dim, step as f64 * 0.1);
family
.batched_outer_gradient_terms(
&states,
&specs,
&deriv_blocks,
&rho_retry,
&opts,
None,
)
.expect("guard ok on retry");
assert_eq!(
counter(),
distinct_calls,
"line-search retry at step {step} must NOT bump counter"
);
}
}
assert_eq!(counter(), 15, "final counter should be 15 distinct ρ");
let family_off = make_block_psi_test_family(n);
let opts_off = BlockwiseFitOptions {
auto_outer_subsample: false,
..BlockwiseFitOptions::default()
};
for step in 0..5 {
let rho_step = Array1::<f64>::from_elem(rho_dim, step as f64 * 0.1);
family_off
.batched_outer_gradient_terms(
&states,
&specs,
&deriv_blocks,
&rho_step,
&opts_off,
None,
)
.expect("guard ok off");
}
assert_eq!(
family_off
.auto_subsample_phase_counter
.load(std::sync::atomic::Ordering::SeqCst),
0,
"with auto_outer_subsample=false the counter must stay at 0"
);
}
}