use crate::cache::{LoadSource, Session as CacheSession};
use crate::estimate::EstimationError;
use crate::solver::estimate::reml::unified::BarrierConfig;
use crate::solver::startup_stats::{
SeedRejection, StartupStats, format_no_seeds_passed, uniform_structural_key,
};
use ::opt::{
Arc as ArcOptimizer, ArcError, Bfgs, BfgsError, Bounds, FallbackPolicy as OptFallbackPolicy,
FirstOrderObjective, FirstOrderSample, FixedPoint, FixedPointError, FixedPointObjective,
FixedPointSample, FixedPointStatus, GradientTolerance, HessianFallbackPolicy,
HessianMaterialization, HessianOperator, HessianValue, MatrixFreeTrustRegion, MaxIterations,
ObjectiveEvalError, OperatorObjective, OperatorSample, OptimizationStatus, OptimizerObserver,
SecondOrderObjective, SecondOrderSample, Solution, StepInfo, Tolerance, ZerothOrderObjective,
};
use ndarray::{Array1, Array2, ArrayView2};
use std::sync::Arc;
use std::sync::atomic::{AtomicBool, AtomicU64, AtomicUsize, Ordering};
const OPERATOR_TRUST_RESTART_RADIUS_FLOOR: f64 = 1.0e-6;
fn outer_strategy_contract_panic(message: impl Into<String>) -> ! {
std::panic::panic_any(message.into())
}
#[derive(Clone, Debug)]
pub struct InnerProgressFeedback {
pub cap: Arc<AtomicUsize>,
pub accepted_iter: Arc<AtomicUsize>,
pub last_iters: Arc<AtomicUsize>,
pub last_converged: Arc<AtomicBool>,
pub ift_residual: Arc<AtomicU64>,
pub accept_rho: Arc<AtomicU64>,
}
impl InnerProgressFeedback {
fn snapshot(&self) -> Option<InnerProgressSnapshot> {
let iters = self.last_iters.load(Ordering::Relaxed);
if iters == 0 {
None
} else {
let residual_bits = self.ift_residual.load(Ordering::Relaxed);
let r = f64::from_bits(residual_bits);
let last_ift_residual = if r.is_finite() && r >= 0.0 {
Some(r)
} else {
None
};
let accept_rho_bits = self.accept_rho.load(Ordering::Relaxed);
let ar = f64::from_bits(accept_rho_bits);
let last_accept_rho = if ar.is_finite() && ar >= 0.0 {
Some(ar)
} else {
None
};
Some(InnerProgressSnapshot {
last_iters: iters,
last_converged: self.last_converged.load(Ordering::Relaxed),
last_ift_residual,
last_accept_rho,
})
}
}
}
#[derive(Clone, Copy, Debug)]
struct InnerProgressSnapshot {
last_iters: usize,
last_converged: bool,
last_ift_residual: Option<f64>,
last_accept_rho: Option<f64>,
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum OuterHessianMaterialization {
Unavailable,
RepeatedHvp,
BatchedHvp,
Explicit,
}
impl OuterHessianMaterialization {
fn is_available(self) -> bool {
!matches!(self, Self::Unavailable)
}
}
#[derive(Debug, Clone)]
pub enum OuterStrategyError {
OperatorShape { reason: String },
NonFiniteHessian { reason: String },
RhoBlockShape { reason: String },
}
impl std::fmt::Display for OuterStrategyError {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
OuterStrategyError::OperatorShape { reason }
| OuterStrategyError::NonFiniteHessian { reason }
| OuterStrategyError::RhoBlockShape { reason } => f.write_str(reason),
}
}
}
impl std::error::Error for OuterStrategyError {}
impl From<OuterStrategyError> for String {
fn from(err: OuterStrategyError) -> String {
err.to_string()
}
}
pub trait OuterHessianOperator: Send + Sync {
fn dim(&self) -> usize;
fn matvec(&self, v: &Array1<f64>) -> Result<Array1<f64>, String>;
fn apply_into(&self, v: &Array1<f64>, out: &mut Array1<f64>) -> Result<(), String> {
let result = self.matvec(v)?;
if result.len() != out.len() {
return Err(format!(
"outer Hessian operator matvec produced length {} but expected {}",
result.len(),
out.len()
));
}
out.assign(&result);
Ok(())
}
fn is_cheap_to_materialize(&self) -> bool {
false
}
fn materialization_capability(&self) -> OuterHessianMaterialization {
if self.is_cheap_to_materialize() {
OuterHessianMaterialization::RepeatedHvp
} else {
OuterHessianMaterialization::Unavailable
}
}
fn mul_mat(&self, factor: ArrayView2<'_, f64>) -> Result<Array2<f64>, String> {
use rayon::iter::{IntoParallelIterator, ParallelIterator};
let dim = self.dim();
if factor.nrows() != dim {
return Err(OuterStrategyError::OperatorShape {
reason: format!(
"outer Hessian operator factor row count mismatch: got {}, expected {}",
factor.nrows(),
dim
),
}
.into());
}
let m = factor.ncols();
let cols: Result<Vec<Array1<f64>>, String> = (0..m)
.into_par_iter()
.map(|j| {
let col = factor.column(j).to_owned();
let hv = self.matvec(&col)?;
if hv.len() != dim {
return Err(OuterStrategyError::OperatorShape {
reason: format!(
"outer Hessian operator matvec length mismatch: got {}, expected {}",
hv.len(),
dim
),
}
.into());
}
Ok(hv)
})
.collect();
let cols = cols?;
let mut out = Array2::<f64>::zeros((dim, m));
for (j, hv) in cols.into_iter().enumerate() {
out.column_mut(j).assign(&hv);
}
Ok(out)
}
fn materialize_dense(&self) -> Result<Array2<f64>, String> {
let dim = self.dim();
let identity = Array2::<f64>::eye(dim);
let mut dense = self.mul_mat(identity.view())?;
if dense.nrows() != dim || dense.ncols() != dim {
return Err(OuterStrategyError::OperatorShape {
reason: format!(
"outer Hessian operator mul_mat returned {}x{}, expected {}x{}",
dense.nrows(),
dense.ncols(),
dim,
dim
),
}
.into());
}
for row in 0..dim {
for col in (row + 1)..dim {
let sym = 0.5 * (dense[[row, col]] + dense[[col, row]]);
dense[[row, col]] = sym;
dense[[col, row]] = sym;
}
}
if !dense.iter().all(|value| value.is_finite()) {
return Err(OuterStrategyError::NonFiniteHessian {
reason: "outer Hessian dense materialization produced non-finite entries"
.to_string(),
}
.into());
}
Ok(dense)
}
}
struct RhoBlockAdditiveOuterHessian {
base: Arc<dyn OuterHessianOperator>,
rho_block: Array2<f64>,
dim: usize,
}
impl OuterHessianOperator for RhoBlockAdditiveOuterHessian {
fn dim(&self) -> usize {
self.dim
}
fn matvec(&self, v: &Array1<f64>) -> Result<Array1<f64>, String> {
if v.len() != self.dim {
return Err(OuterStrategyError::OperatorShape {
reason: format!(
"outer Hessian operator input length mismatch: got {}, expected {}",
v.len(),
self.dim
),
}
.into());
}
let mut out = self.base.matvec(v)?;
let k = self.rho_block.nrows();
if k > 0 {
let rho_v = v.slice(ndarray::s![..k]).to_owned();
let rho_out = self.rho_block.dot(&rho_v);
out.slice_mut(ndarray::s![..k]).scaled_add(1.0, &rho_out);
}
Ok(out)
}
fn mul_mat(&self, factor: ArrayView2<'_, f64>) -> Result<Array2<f64>, String> {
let mut out = self.base.mul_mat(factor)?;
let k = self.rho_block.nrows();
if k > 0 {
if k > out.nrows() {
return Err(OuterStrategyError::RhoBlockShape {
reason: format!(
"rho-block Hessian update shape mismatch: rho_block is {}x{}, mul_mat output has {} rows",
self.rho_block.nrows(),
self.rho_block.ncols(),
out.nrows()
),
}
.into());
}
let factor_top = factor.slice(ndarray::s![..k, ..]);
let rho_contrib = self.rho_block.dot(&factor_top);
out.slice_mut(ndarray::s![..k, ..])
.scaled_add(1.0, &rho_contrib);
}
Ok(out)
}
fn is_cheap_to_materialize(&self) -> bool {
self.base.is_cheap_to_materialize()
}
fn materialization_capability(&self) -> OuterHessianMaterialization {
self.base.materialization_capability()
}
}
pub(crate) const OUTER_HVP_MATERIALIZE_MAX_DIM: usize = 64;
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Derivative {
Analytic,
Unavailable,
}
#[derive(Clone, Copy, Debug, PartialEq)]
pub enum DeclaredHessianForm {
Dense,
Operator {
materialization: OuterHessianMaterialization,
estimated_materialization_cost: Option<f64>,
},
Either,
Unavailable,
}
impl DeclaredHessianForm {
pub const fn is_analytic(self) -> bool {
!matches!(self, DeclaredHessianForm::Unavailable)
}
pub const fn is_operator_only(self) -> bool {
matches!(self, DeclaredHessianForm::Operator { .. })
}
pub const fn is_dense_only(self) -> bool {
matches!(self, DeclaredHessianForm::Dense)
}
}
impl From<Derivative> for DeclaredHessianForm {
fn from(d: Derivative) -> Self {
match d {
Derivative::Analytic => DeclaredHessianForm::Either,
Derivative::Unavailable => DeclaredHessianForm::Unavailable,
}
}
}
const SMALL_OUTER_BFGS_MAX_PARAMS: usize = 8;
const SECOND_ORDER_GEOMETRY_PROBE_MAX_PARAMS: usize = 64;
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct OuterThetaLayout {
pub n_params: usize,
pub psi_dim: usize,
}
impl OuterThetaLayout {
pub const fn new(n_params: usize, psi_dim: usize) -> Self {
Self { n_params, psi_dim }
}
pub const fn rho_dim(&self) -> usize {
self.n_params.saturating_sub(self.psi_dim)
}
fn validate_capability(&self, context: &str) -> Result<(), EstimationError> {
if self.psi_dim > self.n_params {
return Err(EstimationError::RemlOptimizationFailed(format!(
"{context}: invalid outer theta layout (psi_dim={} exceeds n_params={})",
self.psi_dim, self.n_params
)));
}
Ok::<(), _>(())
}
fn validate_point_len(
&self,
theta: &Array1<f64>,
context: &str,
) -> Result<(), ObjectiveEvalError> {
if theta.len() != self.n_params {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: outer theta length mismatch: got {}, expected {} (rho_dim={}, psi_dim={})",
theta.len(),
self.n_params,
self.rho_dim(),
self.psi_dim
)));
}
Ok::<(), _>(())
}
fn validate_gradient_len(
&self,
gradient: &Array1<f64>,
context: &str,
) -> Result<(), ObjectiveEvalError> {
if gradient.len() != self.n_params {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: outer gradient length mismatch: got {}, expected {} (rho_dim={}, psi_dim={})",
gradient.len(),
self.n_params,
self.rho_dim(),
self.psi_dim
)));
}
Ok::<(), _>(())
}
fn validate_hessian_shape(
&self,
hessian: &Array2<f64>,
context: &str,
) -> Result<(), ObjectiveEvalError> {
if hessian.nrows() != self.n_params || hessian.ncols() != self.n_params {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: outer Hessian shape mismatch: got {}x{}, expected {}x{} (rho_dim={}, psi_dim={})",
hessian.nrows(),
hessian.ncols(),
self.n_params,
self.n_params,
self.rho_dim(),
self.psi_dim
)));
}
Ok::<(), _>(())
}
fn validate_efs_eval(&self, eval: &EfsEval, context: &str) -> Result<(), ObjectiveEvalError> {
if eval.steps.len() != self.n_params {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: outer EFS step length mismatch: got {}, expected {} (rho_dim={}, psi_dim={})",
eval.steps.len(),
self.n_params,
self.rho_dim(),
self.psi_dim
)));
}
if let Some(ref psi_gradient) = eval.psi_gradient
&& psi_gradient.len() != self.psi_dim
{
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: outer EFS psi-gradient length mismatch: got {}, expected {}",
psi_gradient.len(),
self.psi_dim
)));
}
if let Some(ref psi_indices) = eval.psi_indices {
if psi_indices.len() != self.psi_dim {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: outer EFS psi-index count mismatch: got {}, expected {}",
psi_indices.len(),
self.psi_dim
)));
}
if psi_indices.iter().any(|&idx| idx >= self.n_params) {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: outer EFS psi index out of range for n_params={}",
self.n_params
)));
}
}
Ok(())
}
}
#[derive(Clone, Debug)]
pub struct OuterCapability {
pub gradient: Derivative,
pub hessian: DeclaredHessianForm,
pub n_params: usize,
pub psi_dim: usize,
pub fixed_point_available: bool,
pub barrier_config: Option<BarrierConfig>,
pub prefer_gradient_only: bool,
pub disable_fixed_point: bool,
}
impl OuterCapability {
pub const fn theta_layout(&self) -> OuterThetaLayout {
OuterThetaLayout::new(self.n_params, self.psi_dim)
}
pub fn validate_layout(&self, context: &str) -> Result<(), EstimationError> {
self.theta_layout().validate_capability(context)
}
pub const fn all_penalty_like(&self) -> bool {
self.psi_dim == 0
}
pub const fn has_psi_coords(&self) -> bool {
self.psi_dim > 0
}
fn efs_plan_eligible(&self) -> bool {
self.fixed_point_available
&& !self.disable_fixed_point
&& self.all_penalty_like()
&& self.n_params > SMALL_OUTER_BFGS_MAX_PARAMS
}
fn hybrid_efs_plan_eligible(&self) -> bool {
self.fixed_point_available
&& !self.disable_fixed_point
&& self.has_psi_coords()
&& self.n_params > SMALL_OUTER_BFGS_MAX_PARAMS
}
fn declared_hessian_for_planning(&self) -> Derivative {
if self.hessian.is_analytic() {
Derivative::Analytic
} else {
Derivative::Unavailable
}
}
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum Solver {
Arc,
Bfgs,
Efs,
HybridEfs,
CompassSearch,
}
#[derive(Clone, Copy, Debug, PartialEq, Eq, Default)]
pub enum SolverClass {
#[default]
Primary,
AuxiliaryGradientFree,
}
#[inline]
fn effective_seed_budget(
requested_budget: usize,
solver: Solver,
risk_profile: crate::seeding::SeedRiskProfile,
screening_enabled: bool,
) -> usize {
let requested_budget = requested_budget.max(1);
let capped = match (solver, risk_profile) {
(Solver::Efs | Solver::HybridEfs, _) => 1,
(Solver::Arc, crate::seeding::SeedRiskProfile::Survival) => 1,
(Solver::Arc, crate::seeding::SeedRiskProfile::GeneralizedLinear) if screening_enabled => 1,
(Solver::Arc, crate::seeding::SeedRiskProfile::GeneralizedLinear) => 2,
(Solver::CompassSearch, _) => 1,
_ => requested_budget,
};
requested_budget.min(capped)
}
#[inline]
fn should_screen_seeds(
config: &OuterConfig,
solver: Solver,
generated_seed_count: usize,
seed_budget: usize,
) -> bool {
config.screening_cap.is_some()
&& generated_seed_count > seed_budget
&& matches!(
solver,
Solver::Arc | Solver::Bfgs | Solver::Efs | Solver::HybridEfs
)
}
#[inline]
fn expensive_unsuccessful_seed_limit(
solver: Solver,
risk_profile: crate::seeding::SeedRiskProfile,
) -> Option<usize> {
match (solver, risk_profile) {
(Solver::Efs | Solver::HybridEfs, _) => Some(1),
(Solver::Arc, crate::seeding::SeedRiskProfile::Survival) => Some(1),
(Solver::Arc, crate::seeding::SeedRiskProfile::GeneralizedLinear) => Some(2),
(Solver::CompassSearch, _) => Some(1),
_ => None,
}
}
const SEED_SCREENING_CASCADE_MULTIPLIERS: [usize; 3] = [1, 4, 16];
const SEED_SCREENING_UNCAPPED: usize = 0;
fn rank_seeds_with_screening(
obj: &mut dyn OuterObjective,
config: &OuterConfig,
context: &str,
seeds: &[Array1<f64>],
) -> Vec<Array1<f64>> {
let Some(screening_cap) = config.screening_cap.as_ref() else {
return seeds.to_vec();
};
let initial_cap = config.seed_config.screen_max_inner_iterations.max(1);
let previous_cap = screening_cap.swap(initial_cap, Ordering::Relaxed);
let cascade_caps = [
initial_cap.saturating_mul(SEED_SCREENING_CASCADE_MULTIPLIERS[0]),
initial_cap.saturating_mul(SEED_SCREENING_CASCADE_MULTIPLIERS[1]),
initial_cap.saturating_mul(SEED_SCREENING_CASCADE_MULTIPLIERS[2]),
SEED_SCREENING_UNCAPPED,
];
let mut ranked: Vec<(usize, f64)> = Vec::with_capacity(seeds.len());
let mut rejected = 0usize;
let mut final_cap_used = initial_cap;
let mut stages_consumed = 0usize;
let cascade_start = std::time::Instant::now();
for (stage, &cap) in cascade_caps.iter().enumerate() {
screening_cap.store(cap, Ordering::Relaxed);
ranked.clear();
rejected = 0;
for (idx, seed) in seeds.iter().enumerate() {
obj.reset();
screening_cap.store(cap, Ordering::Relaxed);
match obj.eval_screening_proxy(seed) {
Ok(cost) if cost.is_finite() => ranked.push((idx, cost)),
Ok(_) | Err(_) => rejected += 1,
}
}
final_cap_used = cap;
stages_consumed = stage + 1;
if !ranked.is_empty() {
if stage > 0 {
log::info!(
"[OUTER] {context}: seed screening cap escalated from {} to {} \
(initial cap was too shallow for this problem; {}/{} seeds ranked)",
initial_cap,
if cap == 0 {
"uncapped".to_string()
} else {
cap.to_string()
},
ranked.len(),
seeds.len(),
);
}
break;
}
}
screening_cap.store(previous_cap, Ordering::Relaxed);
obj.reset();
log::info!(
"[OUTER] {context}: seed screening cascade complete elapsed={:.3}s stages_used={} final_cap={} ranked={}/{}",
cascade_start.elapsed().as_secs_f64(),
stages_consumed,
if final_cap_used == 0 {
"uncapped".to_string()
} else {
final_cap_used.to_string()
},
ranked.len(),
seeds.len(),
);
if ranked.is_empty() {
log::info!(
"[OUTER] {context}: no finite seed cost even with full inner budget \
({} seeds, {} rejected, {} cascade stages tried); keeping heuristic order",
seeds.len(),
rejected,
stages_consumed,
);
return seeds.to_vec();
}
ranked.sort_by(|(idx_a, cost_a), (idx_b, cost_b)| {
cost_a.total_cmp(cost_b).then_with(|| idx_a.cmp(idx_b))
});
let mut ordered = Vec::with_capacity(seeds.len());
let mut seen = vec![false; seeds.len()];
for (idx, _) in ranked {
seen[idx] = true;
ordered.push(seeds[idx].clone());
}
for (idx, seed) in seeds.iter().enumerate() {
if !seen[idx] {
ordered.push(seed.clone());
}
}
log::debug!(
"[OUTER] {context}: seed screening ranked {}/{} candidates at cap={} \
(initial cap={}, stages used={}); rejected={}",
ordered.len() - rejected,
seeds.len(),
if final_cap_used == 0 {
"uncapped".to_string()
} else {
final_cap_used.to_string()
},
initial_cap,
stages_consumed,
rejected,
);
ordered
}
#[inline]
fn candidate_improves_best(candidate: &OuterResult, best: Option<&OuterResult>) -> bool {
match best {
None => true,
Some(best) if candidate.converged != best.converged => candidate.converged,
Some(best) => candidate.final_value < best.final_value,
}
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum HessianSource {
Analytic,
BfgsApprox,
EfsFixedPoint,
HybridEfsFixedPoint,
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum OuterEvalOrder {
ValueAndGradient,
ValueGradientHessian,
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub struct OuterPlan {
pub solver: Solver,
pub hessian_source: HessianSource,
}
pub(crate) const EFS_FIRST_ORDER_FALLBACK_MARKER: &str = "[outer-efs-first-order-fallback]";
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum FallbackPolicy {
Automatic,
Disabled,
}
impl std::fmt::Display for OuterPlan {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
write!(
f,
"solver={:?}, hessian_source={:?}",
self.solver, self.hessian_source
)
}
}
impl OuterPlan {
pub fn routing_log_line(&self) -> String {
let matrix_free = false;
format!(
"solver={:?};hessian={:?};matrix-free={}",
self.solver, self.hessian_source, matrix_free
)
}
}
pub fn plan(cap: &OuterCapability) -> OuterPlan {
use Derivative::*;
use HessianSource as H;
use Solver as S;
match (cap.gradient, cap.declared_hessian_for_planning()) {
(Analytic, Analytic) => OuterPlan {
solver: S::Arc,
hessian_source: H::Analytic,
},
(Analytic, Unavailable) if cap.efs_plan_eligible() => OuterPlan {
solver: S::Efs,
hessian_source: H::EfsFixedPoint,
},
(Unavailable, Unavailable) if cap.efs_plan_eligible() => OuterPlan {
solver: S::Efs,
hessian_source: H::EfsFixedPoint,
},
(Analytic, Unavailable) if cap.hybrid_efs_plan_eligible() => OuterPlan {
solver: S::HybridEfs,
hessian_source: H::HybridEfsFixedPoint,
},
(Unavailable, Unavailable) if cap.hybrid_efs_plan_eligible() => OuterPlan {
solver: S::HybridEfs,
hessian_source: H::HybridEfsFixedPoint,
},
(Analytic, Unavailable) => OuterPlan {
solver: S::Bfgs,
hessian_source: H::BfgsApprox,
},
(Unavailable, _) => OuterPlan {
solver: S::Bfgs,
hessian_source: H::BfgsApprox,
},
}
}
pub fn plan_with_class(cap: &OuterCapability, class: SolverClass) -> OuterPlan {
use Derivative::*;
if class == SolverClass::AuxiliaryGradientFree
&& cap.gradient == Unavailable
&& cap.declared_hessian_for_planning() == Unavailable
&& !cap.efs_plan_eligible()
&& !cap.hybrid_efs_plan_eligible()
{
return OuterPlan {
solver: Solver::CompassSearch,
hessian_source: HessianSource::BfgsApprox,
};
}
plan(cap)
}
pub fn log_plan(context: &str, cap: &OuterCapability, the_plan: &OuterPlan) {
let hess_warning = match the_plan.hessian_source {
HessianSource::BfgsApprox if cap.n_params > 0 => {
" [no Hessian: BFGS approximation]".to_string()
}
_ => String::new(),
};
let barrier_note = if cap.barrier_config.is_some() && cap.efs_plan_eligible() {
" [EFS with runtime barrier-curvature guard]"
} else {
""
};
let hybrid_note = if the_plan.solver == Solver::HybridEfs {
" [hybrid EFS(ρ) + preconditioned-gradient(ψ)]"
} else {
""
};
log::info!(
"[OUTER] {context}: n_params={}, gradient={:?}, hessian={:?} -> {} [{}]{hess_warning}{barrier_note}{hybrid_note}",
cap.n_params,
cap.gradient,
cap.hessian,
the_plan,
the_plan.routing_log_line(),
);
}
fn requests_immediate_first_order_fallback(message: &str) -> bool {
message.contains(EFS_FIRST_ORDER_FALLBACK_MARKER)
}
fn disable_fixed_point(cap: &OuterCapability) -> Option<OuterCapability> {
(!cap.disable_fixed_point && (cap.efs_plan_eligible() || cap.hybrid_efs_plan_eligible())).then(
|| {
let mut degraded = cap.clone();
degraded.disable_fixed_point = true;
degraded
},
)
}
fn automatic_fallback_attempts(cap: &OuterCapability) -> Vec<OuterCapability> {
let mut attempts = Vec::new();
if cap.gradient == Derivative::Analytic
&& matches!(plan(cap).solver, Solver::Efs | Solver::HybridEfs)
&& let Some(no_fp_cap) = disable_fixed_point(cap)
{
attempts.push(no_fp_cap.clone());
return attempts;
}
if matches!(plan(cap).solver, Solver::Arc) {
return attempts;
}
attempts
}
fn disabled_fallback_hybrid_efs_has_standalone_bfgs_primary(
cap: &OuterCapability,
config: &OuterConfig,
) -> bool {
config.solver_class == SolverClass::Primary
&& config.fallback_policy == FallbackPolicy::Disabled
&& cap.gradient == Derivative::Analytic
&& matches!(plan(cap).solver, Solver::HybridEfs)
}
fn primary_capability_for_config(
mut cap: OuterCapability,
config: &OuterConfig,
context: &str,
) -> OuterCapability {
if disabled_fallback_hybrid_efs_has_standalone_bfgs_primary(&cap, config) {
log::info!(
"[OUTER] {context}: HybridEFS requires the automatic first-order \
escape path for ψ coordinates; fallback is disabled, so routing the \
primary attempt to analytic-gradient BFGS"
);
cap.disable_fixed_point = true;
}
cap
}
pub struct OuterEval {
pub cost: f64,
pub gradient: Array1<f64>,
pub hessian: HessianResult,
pub inner_beta_hint: Option<Array1<f64>>,
}
impl OuterEval {
pub fn infeasible(n_params: usize) -> Self {
Self {
cost: f64::INFINITY,
gradient: Array1::zeros(n_params),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
}
}
}
pub enum HessianResult {
Analytic(Array2<f64>),
Operator(Arc<dyn OuterHessianOperator>),
Unavailable,
}
impl Clone for OuterEval {
fn clone(&self) -> Self {
Self {
cost: self.cost,
gradient: self.gradient.clone(),
hessian: self.hessian.clone(),
inner_beta_hint: self.inner_beta_hint.clone(),
}
}
}
impl std::fmt::Debug for OuterEval {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
f.debug_struct("OuterEval")
.field("cost", &self.cost)
.field("gradient", &self.gradient)
.field("hessian", &self.hessian)
.finish()
}
}
impl Clone for HessianResult {
fn clone(&self) -> Self {
match self {
Self::Analytic(h) => Self::Analytic(h.clone()),
Self::Operator(op) => Self::Operator(Arc::clone(op)),
Self::Unavailable => Self::Unavailable,
}
}
}
impl std::fmt::Debug for HessianResult {
fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
match self {
Self::Analytic(h) => f
.debug_tuple("Analytic")
.field(&format!("{}x{}", h.nrows(), h.ncols()))
.finish(),
Self::Operator(op) => f
.debug_tuple("Operator")
.field(&format!("dim={}", op.dim()))
.finish(),
Self::Unavailable => f.write_str("Unavailable"),
}
}
}
impl HessianResult {
pub fn unwrap_analytic(self) -> Array2<f64> {
match self {
HessianResult::Analytic(h) => h,
HessianResult::Operator(_) => {
outer_strategy_contract_panic(
"expected dense analytic Hessian but got HessianResult::Operator",
)
}
HessianResult::Unavailable => {
outer_strategy_contract_panic(
"expected analytic Hessian but got HessianResult::Unavailable",
)
}
}
}
pub fn is_analytic(&self) -> bool {
matches!(
self,
HessianResult::Analytic(_) | HessianResult::Operator(_)
)
}
pub fn into_option(self) -> Option<Array2<f64>> {
match self {
HessianResult::Analytic(h) => Some(h),
HessianResult::Operator(_) => None,
HessianResult::Unavailable => None,
}
}
pub fn dim(&self) -> Option<usize> {
match self {
HessianResult::Analytic(h) => Some(h.nrows()),
HessianResult::Operator(op) => Some(op.dim()),
HessianResult::Unavailable => None,
}
}
pub fn materialize_dense(&self) -> Result<Option<Array2<f64>>, String> {
match self {
HessianResult::Analytic(h) => Ok(Some(h.clone())),
HessianResult::Operator(op) => op.materialize_dense().map(Some),
HessianResult::Unavailable => Ok(None),
}
}
pub fn add_rho_block_dense(&mut self, rho_block: &Array2<f64>) -> Result<(), String> {
if rho_block.nrows() != rho_block.ncols() {
return Err(OuterStrategyError::RhoBlockShape {
reason: format!(
"rho-block Hessian update must be square, got {}x{}",
rho_block.nrows(),
rho_block.ncols()
),
}
.into());
}
match self {
HessianResult::Analytic(h) => {
if rho_block.nrows() > h.nrows() || rho_block.ncols() > h.ncols() {
return Err(OuterStrategyError::RhoBlockShape {
reason: format!(
"rho-block Hessian update shape mismatch: got {}x{}, outer Hessian is {}x{}",
rho_block.nrows(),
rho_block.ncols(),
h.nrows(),
h.ncols()
),
}
.into());
}
let k = rho_block.nrows();
let mut sl = h.slice_mut(ndarray::s![..k, ..k]);
sl += rho_block;
Ok(())
}
HessianResult::Operator(op) => {
let base = Arc::clone(op);
let dim = base.dim();
if rho_block.nrows() > dim {
return Err(OuterStrategyError::RhoBlockShape {
reason: format!(
"rho-block Hessian update dimension mismatch: got {}x{}, operator dim is {}",
rho_block.nrows(),
rho_block.ncols(),
dim
),
}
.into());
}
*self = HessianResult::Operator(Arc::new(RhoBlockAdditiveOuterHessian {
base,
rho_block: rho_block.clone(),
dim,
}));
Ok(())
}
HessianResult::Unavailable => Ok(()),
}
}
}
#[derive(Clone, Debug)]
pub struct EfsEval {
pub cost: f64,
pub steps: Vec<f64>,
pub beta: Option<Array1<f64>>,
pub psi_gradient: Option<Array1<f64>>,
pub psi_indices: Option<Vec<usize>>,
pub inner_hessian_scale: Option<f64>,
}
pub trait OuterObjective {
fn capability(&self) -> OuterCapability;
fn eval_cost(&mut self, rho: &Array1<f64>) -> Result<f64, EstimationError>;
fn eval_screening_proxy(&mut self, rho: &Array1<f64>) -> Result<f64, EstimationError> {
self.eval_cost(rho)
}
fn eval(&mut self, rho: &Array1<f64>) -> Result<OuterEval, EstimationError>;
fn eval_with_order(
&mut self,
rho: &Array1<f64>,
order: OuterEvalOrder,
) -> Result<OuterEval, EstimationError> {
match order {
OuterEvalOrder::ValueAndGradient | OuterEvalOrder::ValueGradientHessian => {
self.eval(rho)
}
}
}
fn eval_efs(&mut self, rho: &Array1<f64>) -> Result<EfsEval, EstimationError> {
Err(EstimationError::RemlOptimizationFailed(format!(
"EFS evaluation not implemented for this objective at rho_dim={}",
rho.len()
)))
}
fn reset(&mut self);
fn seed_inner_state(&mut self, beta: &Array1<f64>) -> Result<(), EstimationError>;
}
#[derive(serde::Serialize, serde::Deserialize)]
struct IteratePayload {
schema: u32,
rho: Vec<f64>,
#[serde(default)]
beta: Vec<f64>,
cost: f64,
eval_id: u64,
}
const ITERATE_PAYLOAD_SCHEMA: u32 = 2;
fn encode_iterate(
rho: &Array1<f64>,
beta: Option<&Array1<f64>>,
cost: f64,
eval_id: u64,
) -> Option<Vec<u8>> {
let p = IteratePayload {
schema: ITERATE_PAYLOAD_SCHEMA,
rho: rho.to_vec(),
beta: beta.map(|b| b.to_vec()).unwrap_or_default(),
cost,
eval_id,
};
serde_json::to_vec(&p).ok()
}
fn decode_iterate(bytes: &[u8], expected_rho_dim: usize) -> Option<IteratePayload> {
let p: IteratePayload = serde_json::from_slice(bytes).ok()?;
if p.schema != ITERATE_PAYLOAD_SCHEMA {
return None;
}
if p.rho.len() != expected_rho_dim {
return None;
}
if !p.rho.iter().all(|x| x.is_finite()) || !p.cost.is_finite() {
return None;
}
if !p.beta.iter().all(|x| x.is_finite()) {
return None;
}
Some(p)
}
#[derive(Debug)]
pub(crate) enum CacheSeedDecision {
ExactFinal {
rho: Array1<f64>,
beta: Vec<f64>,
final_value: f64,
iterations: usize,
prior_obj_display: f64,
},
Seed {
rho: Array1<f64>,
beta: Vec<f64>,
prior_obj_display: f64,
iteration: u64,
},
Discard {
reason: &'static str,
prior_obj_display: f64,
all_rho_finite: Option<bool>,
},
}
pub(crate) fn classify_cache_entry_for_outer(
loaded: &crate::cache::LoadedEntry,
expected_rho_dim: usize,
) -> CacheSeedDecision {
let entry = &loaded.entry;
let Some(payload) = decode_iterate(&entry.payload, expected_rho_dim) else {
return CacheSeedDecision::Discard {
reason: "payload-shape-mismatch",
prior_obj_display: entry.objective.unwrap_or(f64::NAN),
all_rho_finite: None,
};
};
let cached_rho = Array1::from_vec(payload.rho);
let prior_obj_display = entry.objective.unwrap_or(f64::NAN);
if matches!(entry.objective, Some(v) if !v.is_finite()) {
return CacheSeedDecision::Discard {
reason: "non-finite-payload",
prior_obj_display,
all_rho_finite: Some(cached_rho.iter().all(|v| v.is_finite())),
};
}
if !cached_rho.iter().all(|v| v.is_finite()) {
return CacheSeedDecision::Discard {
reason: "non-finite-payload",
prior_obj_display,
all_rho_finite: Some(false),
};
}
if loaded.source == LoadSource::Exact && entry.kind == crate::cache::EntryKind::Final {
return CacheSeedDecision::ExactFinal {
rho: cached_rho,
beta: payload.beta,
final_value: entry.objective.unwrap_or(payload.cost),
iterations: entry
.iteration
.unwrap_or(payload.eval_id)
.min(usize::MAX as u64) as usize,
prior_obj_display,
};
}
CacheSeedDecision::Seed {
rho: cached_rho,
beta: payload.beta,
prior_obj_display,
iteration: entry.iteration.unwrap_or(payload.eval_id),
}
}
pub(crate) fn cache_entry_would_help_outer(
loaded: &crate::cache::LoadedEntry,
expected_rho_dim: usize,
) -> bool {
matches!(
classify_cache_entry_for_outer(loaded, expected_rho_dim),
CacheSeedDecision::ExactFinal { .. } | CacheSeedDecision::Seed { .. }
)
}
struct CheckpointingObjective<'a> {
inner: &'a mut dyn OuterObjective,
session: Arc<CacheSession>,
mirror_sessions: Vec<Arc<CacheSession>>,
eval_counter: AtomicU64,
last_inner_beta: std::sync::Mutex<Option<Array1<f64>>>,
}
impl<'a> CheckpointingObjective<'a> {
fn new(
inner: &'a mut dyn OuterObjective,
session: Arc<CacheSession>,
mirror_sessions: Vec<Arc<CacheSession>>,
) -> Self {
Self {
inner,
session,
mirror_sessions,
eval_counter: AtomicU64::new(0),
last_inner_beta: std::sync::Mutex::new(None),
}
}
fn last_inner_beta(&self) -> Option<Array1<f64>> {
self.last_inner_beta.lock().ok().and_then(|g| g.clone())
}
fn note(&self, rho: &Array1<f64>, beta: Option<&Array1<f64>>, cost: f64) {
if !cost.is_finite() {
return;
}
if let Some(b) = beta {
if !b.iter().all(|v| v.is_finite()) {
return;
}
if let Ok(mut guard) = self.last_inner_beta.lock() {
*guard = Some(b.clone());
}
}
let i = self.eval_counter.fetch_add(1, Ordering::Relaxed);
if let Some(bytes) = encode_iterate(rho, beta, cost, i) {
self.session.checkpoint(&bytes, Some(cost), Some(i));
for mirror in &self.mirror_sessions {
mirror.checkpoint(&bytes, Some(cost), Some(i));
}
}
}
}
impl<'a> OuterObjective for CheckpointingObjective<'a> {
fn capability(&self) -> OuterCapability {
self.inner.capability()
}
fn eval_cost(&mut self, rho: &Array1<f64>) -> Result<f64, EstimationError> {
let v = self.inner.eval_cost(rho)?;
self.note(rho, None, v);
Ok(v)
}
fn eval_screening_proxy(&mut self, rho: &Array1<f64>) -> Result<f64, EstimationError> {
self.inner.eval_screening_proxy(rho)
}
fn eval(&mut self, rho: &Array1<f64>) -> Result<OuterEval, EstimationError> {
let r = self.inner.eval(rho)?;
self.note(rho, r.inner_beta_hint.as_ref(), r.cost);
Ok(r)
}
fn eval_with_order(
&mut self,
rho: &Array1<f64>,
order: OuterEvalOrder,
) -> Result<OuterEval, EstimationError> {
let r = self.inner.eval_with_order(rho, order)?;
self.note(rho, r.inner_beta_hint.as_ref(), r.cost);
Ok(r)
}
fn eval_efs(&mut self, rho: &Array1<f64>) -> Result<EfsEval, EstimationError> {
let r = self.inner.eval_efs(rho)?;
self.note(rho, None, r.cost);
Ok(r)
}
fn seed_inner_state(&mut self, beta: &Array1<f64>) -> Result<(), EstimationError> {
let result = self.inner.seed_inner_state(beta);
if result.is_ok()
&& beta.iter().all(|v| v.is_finite())
&& let Ok(mut guard) = self.last_inner_beta.lock()
{
*guard = Some(beta.clone());
}
result
}
fn reset(&mut self) {
self.inner.reset();
}
}
pub struct ClosureObjective<
S,
Fc,
Fe,
Fr = fn(&mut S),
Fefs = fn(&mut S, &Array1<f64>) -> Result<EfsEval, EstimationError>,
Feo = fn(&mut S, &Array1<f64>, OuterEvalOrder) -> Result<OuterEval, EstimationError>,
Fsp = fn(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
Fseed = fn(&mut S, &Array1<f64>) -> Result<(), EstimationError>,
> {
pub(crate) state: S,
pub(crate) cap: OuterCapability,
pub(crate) cost_fn: Fc,
pub(crate) eval_fn: Fe,
pub(crate) eval_order_fn: Option<Feo>,
pub(crate) reset_fn: Option<Fr>,
pub(crate) efs_fn: Option<Fefs>,
pub(crate) screening_proxy_fn: Option<Fsp>,
pub(crate) seed_fn: Option<Fseed>,
}
impl<S, Fc, Fe, Fr, Fefs, Feo, Fsp, Fseed> OuterObjective
for ClosureObjective<S, Fc, Fe, Fr, Fefs, Feo, Fsp, Fseed>
where
Fc: FnMut(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
Fe: FnMut(&mut S, &Array1<f64>) -> Result<OuterEval, EstimationError>,
Fr: FnMut(&mut S),
Fefs: FnMut(&mut S, &Array1<f64>) -> Result<EfsEval, EstimationError>,
Feo: FnMut(&mut S, &Array1<f64>, OuterEvalOrder) -> Result<OuterEval, EstimationError>,
Fsp: FnMut(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
Fseed: FnMut(&mut S, &Array1<f64>) -> Result<(), EstimationError>,
{
fn capability(&self) -> OuterCapability {
self.cap.clone()
}
fn eval_cost(&mut self, rho: &Array1<f64>) -> Result<f64, EstimationError> {
crate::solver::estimate::reml::runtime::record_current_outer_theta_for_ift(rho);
(self.cost_fn)(&mut self.state, rho)
}
fn eval_screening_proxy(&mut self, rho: &Array1<f64>) -> Result<f64, EstimationError> {
crate::solver::estimate::reml::runtime::record_current_outer_theta_for_ift(rho);
match self.screening_proxy_fn.as_mut() {
Some(f) => f(&mut self.state, rho),
None => (self.cost_fn)(&mut self.state, rho),
}
}
fn eval(&mut self, rho: &Array1<f64>) -> Result<OuterEval, EstimationError> {
crate::solver::estimate::reml::runtime::record_current_outer_theta_for_ift(rho);
(self.eval_fn)(&mut self.state, rho)
}
fn eval_with_order(
&mut self,
rho: &Array1<f64>,
order: OuterEvalOrder,
) -> Result<OuterEval, EstimationError> {
crate::solver::estimate::reml::runtime::record_current_outer_theta_for_ift(rho);
match self.eval_order_fn.as_mut() {
Some(f) => f(&mut self.state, rho, order),
None => (self.eval_fn)(&mut self.state, rho),
}
}
fn eval_efs(&mut self, rho: &Array1<f64>) -> Result<EfsEval, EstimationError> {
crate::solver::estimate::reml::runtime::record_current_outer_theta_for_ift(rho);
match self.efs_fn.as_mut() {
Some(f) => f(&mut self.state, rho),
None => Err(EstimationError::RemlOptimizationFailed(
"EFS evaluation not implemented for this objective".to_string(),
)),
}
}
fn seed_inner_state(&mut self, beta: &Array1<f64>) -> Result<(), EstimationError> {
if beta.is_empty() {
return Ok(());
}
if let Some(f) = self.seed_fn.as_mut() {
return f(&mut self.state, beta);
}
Err(EstimationError::InvalidInput(format!(
"cached inner beta has length {}, but this objective does not expose an inner-state seeding hook",
beta.len()
)))
}
fn reset(&mut self) {
if let Some(f) = self.reset_fn.as_mut() {
f(&mut self.state);
}
}
}
impl<S, Fc, Fe, Fr, Fefs, Feo, Fsp> ClosureObjective<S, Fc, Fe, Fr, Fefs, Feo, Fsp>
where
Fc: FnMut(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
Fe: FnMut(&mut S, &Array1<f64>) -> Result<OuterEval, EstimationError>,
Fr: FnMut(&mut S),
Fefs: FnMut(&mut S, &Array1<f64>) -> Result<EfsEval, EstimationError>,
Feo: FnMut(&mut S, &Array1<f64>, OuterEvalOrder) -> Result<OuterEval, EstimationError>,
Fsp: FnMut(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
{
pub fn with_seed_inner_state<Fseed>(
self,
seed_fn: Fseed,
) -> ClosureObjective<S, Fc, Fe, Fr, Fefs, Feo, Fsp, Fseed>
where
Fseed: FnMut(&mut S, &Array1<f64>) -> Result<(), EstimationError>,
{
ClosureObjective {
state: self.state,
cap: self.cap,
cost_fn: self.cost_fn,
eval_fn: self.eval_fn,
eval_order_fn: self.eval_order_fn,
reset_fn: self.reset_fn,
efs_fn: self.efs_fn,
screening_proxy_fn: self.screening_proxy_fn,
seed_fn: Some(seed_fn),
}
}
}
fn into_objective_error(context: &str, err: EstimationError) -> ObjectiveEvalError {
ObjectiveEvalError::recoverable(format!("{context}: {err}"))
}
fn finite_cost_or_error(context: &str, cost: f64) -> Result<f64, ObjectiveEvalError> {
if cost.is_finite() {
Ok(cost)
} else {
Err(ObjectiveEvalError::recoverable(format!(
"{context}: objective returned a non-finite cost"
)))
}
}
fn finite_outer_eval_or_error(
context: &str,
layout: OuterThetaLayout,
eval: OuterEval,
) -> Result<OuterEval, ObjectiveEvalError> {
layout.validate_gradient_len(&eval.gradient, context)?;
if !eval.cost.is_finite() {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: objective returned a non-finite cost"
)));
}
if !eval.gradient.iter().all(|v| v.is_finite()) {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: objective returned a non-finite gradient"
)));
}
match &eval.hessian {
HessianResult::Analytic(hessian) => {
layout.validate_hessian_shape(hessian, context)?;
if !hessian.iter().all(|v| v.is_finite()) {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: objective returned a non-finite Hessian"
)));
}
}
HessianResult::Operator(op) => {
if op.dim() != layout.n_params {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: outer Hessian operator dimension mismatch: got {}, expected {} (rho_dim={}, psi_dim={})",
op.dim(),
layout.n_params,
layout.rho_dim(),
layout.psi_dim
)));
}
}
HessianResult::Unavailable => {}
}
Ok(eval)
}
fn finite_outer_first_order_eval_or_error(
context: &str,
layout: OuterThetaLayout,
eval: OuterEval,
) -> Result<OuterEval, ObjectiveEvalError> {
layout.validate_gradient_len(&eval.gradient, context)?;
if !eval.cost.is_finite() {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: objective returned a non-finite cost"
)));
}
if !eval.gradient.iter().all(|v| v.is_finite()) {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: objective returned a non-finite gradient"
)));
}
Ok(eval)
}
fn validate_second_order_seed_hessian(
context: &str,
layout: OuterThetaLayout,
eval: &OuterEval,
) -> Result<(), ObjectiveEvalError> {
if layout.n_params > SECOND_ORDER_GEOMETRY_PROBE_MAX_PARAMS || !eval.hessian.is_analytic() {
return Ok(());
}
if matches!(
&eval.hessian,
HessianResult::Operator(op) if !op.materialization_capability().is_available()
) {
return Ok(());
}
let Some(hessian) = eval.hessian.materialize_dense().map_err(|message| {
ObjectiveEvalError::recoverable(format!(
"{context}: analytic outer Hessian materialization failed during second-order seed validation: {message}"
))
})?
else {
return Ok(());
};
layout.validate_hessian_shape(&hessian, context)?;
if !hessian.iter().all(|value| value.is_finite()) {
return Err(ObjectiveEvalError::recoverable(format!(
"{context}: analytic outer Hessian probe encountered non-finite entries"
)));
}
Ok(())
}
struct OuterFirstOrderBridge<'a> {
obj: &'a mut dyn OuterObjective,
layout: OuterThetaLayout,
outer_inner_cap: Option<InnerProgressFeedback>,
iter_count: usize,
g_norm_initial: Option<f64>,
last_g_norm: Option<f64>,
}
impl ZerothOrderObjective for OuterFirstOrderBridge<'_> {
fn eval_cost(&mut self, x: &Array1<f64>) -> Result<f64, ObjectiveEvalError> {
self.layout
.validate_point_len(x, "outer eval_cost failed")?;
let cost = self
.obj
.eval_cost(x)
.map_err(|err| into_objective_error("outer eval_cost failed", err))?;
finite_cost_or_error("outer eval_cost failed", cost)
}
}
impl FirstOrderObjective for OuterFirstOrderBridge<'_> {
fn eval_grad(&mut self, x: &Array1<f64>) -> Result<FirstOrderSample, ObjectiveEvalError> {
self.layout.validate_point_len(x, "outer eval failed")?;
if let Some(feedback) = self.outer_inner_cap.as_ref() {
let g_ratio = match (self.last_g_norm, self.g_norm_initial) {
(Some(g), Some(g0)) if g0 > 0.0 => Some(g / g0),
_ => None,
};
let snapshot = feedback.snapshot();
let accepted_iter = feedback.accepted_iter.load(Ordering::Relaxed);
let cap = first_order_inner_cap_schedule(accepted_iter, g_ratio, snapshot);
let prev = feedback.cap.swap(cap, Ordering::Relaxed);
if prev != cap {
let ratio_str = match g_ratio {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
};
let snap_str = match snapshot {
Some(s) => format!(
"last_iters={} converged={} ift_residual={} accept_rho={}",
s.last_iters,
s.last_converged,
match s.last_ift_residual {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
},
match s.last_accept_rho {
Some(r) => format!("{:.3}", r),
None => "n/a".to_string(),
},
),
None => "no-history".to_string(),
};
log::info!(
"[OUTER schedule] inner-PIRLS cap transition accepted_iter={} eval_count={} g_ratio={} {} prev={} new={} ({})",
accepted_iter,
self.iter_count,
ratio_str,
snap_str,
prev,
cap,
if cap == 0 { "uncapped" } else { "capped" }
);
}
}
let stage_start = std::time::Instant::now();
log::info!(
"[STAGE] outer eval start order=ValueAndGradient dim={} (first-order bridge, iter={})",
x.len(),
self.iter_count
);
let eval = self
.obj
.eval_with_order(x, OuterEvalOrder::ValueAndGradient)
.map_err(|err| into_objective_error("outer eval failed", err))?;
let eval = finite_outer_first_order_eval_or_error("outer eval failed", self.layout, eval)?;
let g_norm = eval.gradient.iter().map(|v| v * v).sum::<f64>().sqrt();
let gradient = eval.gradient;
if self.g_norm_initial.is_none() && g_norm.is_finite() && g_norm > 0.0 {
self.g_norm_initial = Some(g_norm);
}
if g_norm.is_finite() {
self.last_g_norm = Some(g_norm);
}
log::info!(
"[STAGE] outer eval end order=ValueAndGradient elapsed={:.3}s cost={:.6e} |g|={:.3e} (first-order bridge, iter={})",
stage_start.elapsed().as_secs_f64(),
eval.cost,
g_norm,
self.iter_count,
);
crate::solver::visualizer::record_outer_eval(eval.cost, g_norm);
self.iter_count = self.iter_count.saturating_add(1);
Ok(FirstOrderSample {
value: eval.cost,
gradient,
})
}
}
fn first_order_inner_cap_schedule(
iter_count: usize,
g_ratio: Option<f64>,
last: Option<InnerProgressSnapshot>,
) -> usize {
if matches!(g_ratio, Some(r) if r < 0.01) {
return 0;
}
if let Some(snap) = last {
let next = if snap.last_converged {
let mut margin = match snap.last_ift_residual {
Some(r) if r < 0.01 => 1usize,
Some(r) if r >= 0.10 => 4usize,
_ => 2usize,
};
if matches!(snap.last_accept_rho, Some(r) if r < 0.5) {
margin = margin.saturating_add(2);
}
snap.last_iters.saturating_add(margin)
} else {
let multiplier = if matches!(snap.last_accept_rho, Some(r) if r < 0.3) {
3
} else {
2
};
snap.last_iters
.saturating_mul(multiplier)
.max(snap.last_iters.saturating_add(4))
};
return next.clamp(3, 64);
}
match iter_count {
0 => 3,
1 => 5,
_ => 10,
}
}
#[cfg(test)]
mod outer_inner_cap_schedule_tests {
use super::{InnerProgressSnapshot, first_order_inner_cap_schedule};
fn snap(last_iters: usize, last_converged: bool) -> Option<InnerProgressSnapshot> {
Some(InnerProgressSnapshot {
last_iters,
last_converged,
last_ift_residual: None,
last_accept_rho: None,
})
}
fn snap_with_accept_rho(
last_iters: usize,
last_converged: bool,
accept_rho: f64,
) -> Option<InnerProgressSnapshot> {
Some(InnerProgressSnapshot {
last_iters,
last_converged,
last_ift_residual: None,
last_accept_rho: Some(accept_rho),
})
}
fn snap_with_residual(
last_iters: usize,
last_converged: bool,
residual: f64,
) -> Option<InnerProgressSnapshot> {
Some(InnerProgressSnapshot {
last_iters,
last_converged,
last_ift_residual: Some(residual),
last_accept_rho: None,
})
}
#[test]
fn snapshot_distinguishes_zero_residual_from_no_signal() {
use super::InnerProgressFeedback;
use crate::solver::estimate::reml::runtime::IFT_RESIDUAL_NO_SIGNAL_BITS;
use std::sync::Arc;
use std::sync::atomic::{AtomicBool, AtomicU64, AtomicUsize};
let make_feedback =
|iters: usize, converged: bool, residual_bits: u64| InnerProgressFeedback {
cap: Arc::new(AtomicUsize::new(0)),
accepted_iter: Arc::new(AtomicUsize::new(0)),
last_iters: Arc::new(AtomicUsize::new(iters)),
last_converged: Arc::new(AtomicBool::new(converged)),
ift_residual: Arc::new(AtomicU64::new(residual_bits)),
accept_rho: Arc::new(AtomicU64::new(
crate::solver::estimate::reml::runtime::IFT_RESIDUAL_NO_SIGNAL_BITS,
)),
};
let fb = make_feedback(5, true, IFT_RESIDUAL_NO_SIGNAL_BITS);
let snap = fb.snapshot().expect("iters > 0, snapshot present");
assert!(
snap.last_ift_residual.is_none(),
"sentinel must decode to None"
);
let fb = make_feedback(5, true, 0.0_f64.to_bits());
let snap = fb.snapshot().expect("iters > 0, snapshot present");
assert_eq!(
snap.last_ift_residual,
Some(0.0),
"residual of exactly 0.0 must round-trip as a real signal, \
not be confused with the no-signal sentinel",
);
let fb = make_feedback(5, true, 0.05_f64.to_bits());
let snap = fb.snapshot().expect("snapshot present");
assert_eq!(snap.last_ift_residual, Some(0.05));
let fb = make_feedback(0, false, IFT_RESIDUAL_NO_SIGNAL_BITS);
assert!(fb.snapshot().is_none());
}
#[test]
fn schedule_falls_back_to_iter_tier_without_feedback() {
assert_eq!(first_order_inner_cap_schedule(0, None, None), 3);
assert_eq!(first_order_inner_cap_schedule(1, None, None), 5);
assert_eq!(first_order_inner_cap_schedule(2, None, None), 10);
assert_eq!(first_order_inner_cap_schedule(20, None, None), 10);
}
#[test]
fn schedule_uses_last_iters_plus_margin_when_converged() {
assert_eq!(first_order_inner_cap_schedule(2, None, snap(4, true)), 6);
assert_eq!(first_order_inner_cap_schedule(5, None, snap(12, true)), 14);
}
#[test]
fn schedule_geometric_backoff_when_last_hit_cap() {
assert_eq!(first_order_inner_cap_schedule(2, None, snap(5, false)), 10);
assert_eq!(first_order_inner_cap_schedule(2, None, snap(1, false)), 5);
assert_eq!(first_order_inner_cap_schedule(2, None, snap(30, false)), 60);
}
#[test]
fn schedule_clamps_floor_and_ceiling() {
assert_eq!(first_order_inner_cap_schedule(2, None, snap(0, true)), 3);
assert_eq!(first_order_inner_cap_schedule(2, None, snap(100, true)), 64);
}
#[test]
fn schedule_uncaps_when_outer_converged() {
assert_eq!(first_order_inner_cap_schedule(0, Some(0.0001), None), 0);
assert_eq!(
first_order_inner_cap_schedule(0, Some(0.005), snap(4, true)),
0
);
assert_eq!(
first_order_inner_cap_schedule(20, Some(0.001), snap(50, false)),
0
);
}
#[test]
fn schedule_ignores_modest_g_ratio_decay() {
assert_eq!(
first_order_inner_cap_schedule(2, Some(0.30), snap(4, true)),
6
);
assert_eq!(
first_order_inner_cap_schedule(2, Some(0.05), snap(4, true)),
6
);
}
#[test]
fn schedule_uses_ift_residual_to_pick_margin() {
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_residual(4, true, 0.005)),
5
);
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_residual(4, true, 0.0001)),
5
);
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_residual(4, true, 0.05)),
6
);
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_residual(4, true, 0.20)),
8
);
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_residual(4, true, 0.80)),
8
);
let residuals = [0.001, 0.05, 0.30];
let caps: Vec<usize> = residuals
.iter()
.map(|&r| first_order_inner_cap_schedule(2, None, snap_with_residual(4, true, r)))
.collect();
for w in caps.windows(2) {
assert!(
w[0] <= w[1],
"ift-residual margin policy regressed monotonicity: {caps:?}"
);
}
}
#[test]
fn schedule_bumps_margin_on_poor_lm_accept_rho() {
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_accept_rho(4, true, 0.95)),
6
);
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_accept_rho(4, true, 0.5)),
6
);
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_accept_rho(4, true, 0.4)),
8
);
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_accept_rho(4, true, 0.1)),
8
);
assert_eq!(
first_order_inner_cap_schedule(2, None, snap_with_accept_rho(4, true, 0.49)),
8
);
}
#[test]
fn schedule_escalates_geometric_backoff_on_very_poor_accept_rho() {
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: false,
last_ift_residual: None,
last_accept_rho: Some(0.15),
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 12);
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: false,
last_ift_residual: None,
last_accept_rho: Some(0.4),
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 8);
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: false,
last_ift_residual: None,
last_accept_rho: Some(0.9),
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 8);
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: false,
last_ift_residual: None,
last_accept_rho: None,
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 8);
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: false,
last_ift_residual: None,
last_accept_rho: Some(0.3),
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 8);
}
#[test]
fn schedule_skips_lm_accept_rho_bump_when_signal_absent() {
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: true,
last_ift_residual: None,
last_accept_rho: None,
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 6);
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: true,
last_ift_residual: None,
last_accept_rho: Some(0.5),
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 6);
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: true,
last_ift_residual: None,
last_accept_rho: Some(1.0),
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 6);
}
#[test]
fn schedule_combines_ift_residual_and_lm_accept_rho() {
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: true,
last_ift_residual: Some(0.30),
last_accept_rho: Some(0.20),
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 10);
let snap = Some(InnerProgressSnapshot {
last_iters: 4,
last_converged: true,
last_ift_residual: Some(0.005),
last_accept_rho: Some(0.30),
});
assert_eq!(first_order_inner_cap_schedule(2, None, snap), 7);
}
}
struct OuterSecondOrderBridge<'a> {
obj: &'a mut dyn OuterObjective,
layout: OuterThetaLayout,
hessian_source: HessianSource,
materialize_operator_max_dim: usize,
eval_count: usize,
outer_inner_cap: Option<InnerProgressFeedback>,
g_norm_initial: Option<f64>,
last_g_norm: Option<f64>,
}
impl ZerothOrderObjective for OuterSecondOrderBridge<'_> {
fn eval_cost(&mut self, x: &Array1<f64>) -> Result<f64, ObjectiveEvalError> {
self.layout
.validate_point_len(x, "outer eval_cost failed")?;
let cost = self
.obj
.eval_cost(x)
.map_err(|err| into_objective_error("outer eval_cost failed", err))?;
finite_cost_or_error("outer eval_cost failed", cost)
}
}
impl FirstOrderObjective for OuterSecondOrderBridge<'_> {
fn eval_grad(&mut self, x: &Array1<f64>) -> Result<FirstOrderSample, ObjectiveEvalError> {
self.layout.validate_point_len(x, "outer eval failed")?;
if let Some(feedback) = self.outer_inner_cap.as_ref() {
let arc_iter = feedback.accepted_iter.load(Ordering::Relaxed);
let g_ratio = match (self.last_g_norm, self.g_norm_initial) {
(Some(g), Some(g0)) if g0 > 0.0 => Some(g / g0),
_ => None,
};
let snapshot = feedback.snapshot();
let cap = first_order_inner_cap_schedule(arc_iter, g_ratio, snapshot);
let prev = feedback.cap.swap(cap, Ordering::Relaxed);
if prev != cap {
let ratio_str = match g_ratio {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
};
let snap_str = match snapshot {
Some(s) => format!(
"last_iters={} converged={} ift_residual={} accept_rho={}",
s.last_iters,
s.last_converged,
match s.last_ift_residual {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
},
match s.last_accept_rho {
Some(r) => format!("{:.3}", r),
None => "n/a".to_string(),
},
),
None => "no-history".to_string(),
};
log::info!(
"[OUTER schedule] inner-PIRLS cap transition (ARC bridge) arc_iter={} g_ratio={} {} prev={} new={} ({})",
arc_iter,
ratio_str,
snap_str,
prev,
cap,
if cap == 0 { "uncapped" } else { "capped" }
);
}
}
let stage_start = std::time::Instant::now();
log::info!(
"[STAGE] outer eval start order=ValueAndGradient dim={}",
x.len()
);
let eval = self
.obj
.eval_with_order(x, OuterEvalOrder::ValueAndGradient)
.map_err(|err| into_objective_error("outer eval failed", err))?;
let eval = finite_outer_first_order_eval_or_error("outer eval failed", self.layout, eval)?;
self.eval_count += 1;
let g_norm = eval.gradient.iter().map(|v| v * v).sum::<f64>().sqrt();
if self.g_norm_initial.is_none() && g_norm.is_finite() && g_norm > 0.0 {
self.g_norm_initial = Some(g_norm);
}
if g_norm.is_finite() {
self.last_g_norm = Some(g_norm);
}
log::info!(
"[STAGE] outer eval end order=ValueAndGradient elapsed={:.3}s cost={:.6e} |g|={:.3e}",
stage_start.elapsed().as_secs_f64(),
eval.cost,
g_norm,
);
log::info!(
"[OUTER] eval#{n} (grad) cost={cost:.6e} |g|={gnorm:.3e} rho=[{rho}]",
n = self.eval_count,
cost = eval.cost,
gnorm = g_norm,
rho = x
.iter()
.map(|v| format!("{v:.3}"))
.collect::<Vec<_>>()
.join(","),
);
crate::solver::visualizer::record_outer_eval(eval.cost, g_norm);
Ok(FirstOrderSample {
value: eval.cost,
gradient: eval.gradient,
})
}
}
impl SecondOrderObjective for OuterSecondOrderBridge<'_> {
fn eval_hessian(&mut self, x: &Array1<f64>) -> Result<SecondOrderSample, ObjectiveEvalError> {
self.layout.validate_point_len(x, "outer eval failed")?;
if let Some(feedback) = self.outer_inner_cap.as_ref() {
let arc_iter = feedback.accepted_iter.load(Ordering::Relaxed);
let g_ratio = match (self.last_g_norm, self.g_norm_initial) {
(Some(g), Some(g0)) if g0 > 0.0 => Some(g / g0),
_ => None,
};
let snapshot = feedback.snapshot();
let cap = first_order_inner_cap_schedule(arc_iter, g_ratio, snapshot);
let prev = feedback.cap.swap(cap, Ordering::Relaxed);
if prev != cap {
let ratio_str = match g_ratio {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
};
let snap_str = match snapshot {
Some(s) => format!(
"last_iters={} converged={} ift_residual={} accept_rho={}",
s.last_iters,
s.last_converged,
match s.last_ift_residual {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
},
match s.last_accept_rho {
Some(r) => format!("{:.3}", r),
None => "n/a".to_string(),
},
),
None => "no-history".to_string(),
};
log::info!(
"[OUTER schedule] inner-PIRLS cap transition (ARC bridge) arc_iter={} g_ratio={} {} prev={} new={} ({})",
arc_iter,
ratio_str,
snap_str,
prev,
cap,
if cap == 0 { "uncapped" } else { "capped" }
);
}
}
let stage_start = std::time::Instant::now();
log::info!(
"[STAGE] outer eval start order=ValueGradientHessian dim={}",
x.len()
);
let eval = self
.obj
.eval_with_order(x, OuterEvalOrder::ValueGradientHessian)
.map_err(|err| into_objective_error("outer eval failed", err))?;
let eval = finite_outer_eval_or_error("outer eval failed", self.layout, eval)?;
self.eval_count += 1;
let g_norm = eval.gradient.iter().map(|v| v * v).sum::<f64>().sqrt();
if self.g_norm_initial.is_none() && g_norm.is_finite() && g_norm > 0.0 {
self.g_norm_initial = Some(g_norm);
}
if g_norm.is_finite() {
self.last_g_norm = Some(g_norm);
}
log::info!(
"[STAGE] outer eval end order=ValueGradientHessian elapsed={:.3}s cost={:.6e} |g|={:.3e}",
stage_start.elapsed().as_secs_f64(),
eval.cost,
g_norm,
);
log::info!(
"[OUTER] eval#{n} (hess) cost={cost:.6e} |g|={gnorm:.3e} rho=[{rho}]",
n = self.eval_count,
cost = eval.cost,
gnorm = g_norm,
rho = x
.iter()
.map(|v| format!("{v:.3}"))
.collect::<Vec<_>>()
.join(","),
);
let hessian = build_bridge_hessian_for_source(
self.hessian_source,
eval.hessian,
self.materialize_operator_max_dim,
)?;
Ok(SecondOrderSample {
value: eval.cost,
gradient: eval.gradient,
hessian,
})
}
}
struct OuterAcceptObserver {
feedback: InnerProgressFeedback,
}
impl OptimizerObserver for OuterAcceptObserver {
fn on_step_accepted(&mut self, info: &StepInfo) {
log::trace!(
"outer step accepted iter={} step_norm={:.3e} predicted_decrease={:.3e} actual_decrease={:.3e}",
info.iter,
info.step_norm,
info.predicted_decrease,
info.actual_decrease,
);
self.feedback.accepted_iter.fetch_add(1, Ordering::Relaxed);
}
}
struct OuterToOptHessianOperator(Arc<dyn OuterHessianOperator>);
impl HessianOperator for OuterToOptHessianOperator {
fn dim(&self) -> usize {
self.0.dim()
}
fn apply_into(&self, v: &Array1<f64>, out: &mut Array1<f64>) -> Result<(), ObjectiveEvalError> {
self.0
.apply_into(v, out)
.map_err(|message| ObjectiveEvalError::Fatal {
message: format!("outer Hessian operator apply_into failed: {message}"),
})
}
fn apply_mat(&self, x: ArrayView2<'_, f64>) -> Result<Array2<f64>, ObjectiveEvalError> {
self.0
.mul_mat(x)
.map_err(|message| ObjectiveEvalError::Fatal {
message: format!("outer Hessian operator mul_mat failed: {message}"),
})
}
fn materialization(&self) -> HessianMaterialization {
match self.0.materialization_capability() {
OuterHessianMaterialization::Unavailable => HessianMaterialization::Unavailable,
OuterHessianMaterialization::RepeatedHvp => HessianMaterialization::RepeatedHvp,
OuterHessianMaterialization::BatchedHvp => HessianMaterialization::BatchedHvp,
OuterHessianMaterialization::Explicit => HessianMaterialization::Explicit,
}
}
fn materialize_dense(&self) -> Result<Array2<f64>, ObjectiveEvalError> {
self.0
.materialize_dense()
.map_err(|message| ObjectiveEvalError::Fatal {
message: format!("outer Hessian operator materialization failed: {message}"),
})
}
}
fn hessian_result_to_value(hessian: HessianResult) -> HessianValue {
match hessian {
HessianResult::Analytic(h) => HessianValue::Dense(h),
HessianResult::Operator(op) => {
HessianValue::Operator(Arc::new(OuterToOptHessianOperator(op)))
}
HessianResult::Unavailable => HessianValue::Unavailable,
}
}
struct OuterOperatorBridge<'a> {
obj: &'a mut dyn OuterObjective,
layout: OuterThetaLayout,
outer_inner_cap: Option<InnerProgressFeedback>,
eval_count: usize,
g_norm_initial: Option<f64>,
last_g_norm: Option<f64>,
}
impl ZerothOrderObjective for OuterOperatorBridge<'_> {
fn eval_cost(&mut self, x: &Array1<f64>) -> Result<f64, ObjectiveEvalError> {
self.layout
.validate_point_len(x, "outer eval_cost failed")?;
let cost = self
.obj
.eval_cost(x)
.map_err(|err| into_objective_error("outer eval_cost failed", err))?;
finite_cost_or_error("outer eval_cost failed", cost)
}
}
impl FirstOrderObjective for OuterOperatorBridge<'_> {
fn eval_grad(&mut self, x: &Array1<f64>) -> Result<FirstOrderSample, ObjectiveEvalError> {
self.layout.validate_point_len(x, "outer eval failed")?;
let eval = self
.obj
.eval_with_order(x, OuterEvalOrder::ValueAndGradient)
.map_err(|err| into_objective_error("outer eval failed", err))?;
let eval = finite_outer_first_order_eval_or_error("outer eval failed", self.layout, eval)?;
let g_norm = eval.gradient.iter().map(|v| v * v).sum::<f64>().sqrt();
if self.g_norm_initial.is_none() && g_norm.is_finite() && g_norm > 0.0 {
self.g_norm_initial = Some(g_norm);
}
if g_norm.is_finite() {
self.last_g_norm = Some(g_norm);
}
Ok(FirstOrderSample {
value: eval.cost,
gradient: eval.gradient,
})
}
}
impl OperatorObjective for OuterOperatorBridge<'_> {
fn eval_value_grad_op(
&mut self,
x: &Array1<f64>,
) -> Result<OperatorSample, ObjectiveEvalError> {
self.layout.validate_point_len(x, "outer eval failed")?;
if let Some(feedback) = self.outer_inner_cap.as_ref() {
let g_ratio = match (self.last_g_norm, self.g_norm_initial) {
(Some(g), Some(g0)) if g0 > 0.0 => Some(g / g0),
_ => None,
};
let snapshot = feedback.snapshot();
let cap = first_order_inner_cap_schedule(self.eval_count, g_ratio, snapshot);
let previous_cap = feedback.cap.swap(cap, Ordering::Relaxed);
if previous_cap != cap {
log::trace!("outer operator bridge updated inner cap from {previous_cap} to {cap}");
}
}
let stage_start = std::time::Instant::now();
log::info!(
"[STAGE] outer eval start order=ValueGradientHessian dim={} (operator bridge)",
x.len(),
);
let eval = self
.obj
.eval_with_order(x, OuterEvalOrder::ValueGradientHessian)
.map_err(|err| into_objective_error("outer eval failed", err))?;
let eval = finite_outer_eval_or_error("outer eval failed", self.layout, eval)?;
self.eval_count += 1;
let g_norm = eval.gradient.iter().map(|v| v * v).sum::<f64>().sqrt();
if self.g_norm_initial.is_none() && g_norm.is_finite() && g_norm > 0.0 {
self.g_norm_initial = Some(g_norm);
}
if g_norm.is_finite() {
self.last_g_norm = Some(g_norm);
}
log::info!(
"[STAGE] outer eval end elapsed={:.3}s cost={:.6e} |g|={:.3e} (operator bridge)",
stage_start.elapsed().as_secs_f64(),
eval.cost,
g_norm,
);
Ok(OperatorSample {
value: eval.cost,
gradient: eval.gradient,
hessian: hessian_result_to_value(eval.hessian),
})
}
}
#[inline]
fn project_to_bounds(x: &Array1<f64>, bounds: Option<&(Array1<f64>, Array1<f64>)>) -> Array1<f64> {
match bounds {
Some((lower, upper)) => {
let mut out = x.clone();
for idx in 0..out.len() {
out[idx] = out[idx].clamp(lower[idx], upper[idx]);
}
out
}
None => x.clone(),
}
}
fn build_bridge_hessian_for_source(
source: HessianSource,
hessian: HessianResult,
materialize_operator_max_dim: usize,
) -> Result<Option<Array2<f64>>, ObjectiveEvalError> {
match source {
HessianSource::Analytic => match hessian {
HessianResult::Analytic(h) => Ok(Some(h)),
HessianResult::Operator(op)
if op.materialization_capability().is_available()
&& op.dim() <= materialize_operator_max_dim =>
{
op.materialize_dense()
.map(Some)
.map_err(|message| ObjectiveEvalError::Fatal {
message: format!(
"outer Hessian operator materialization failed: {message}"
),
})
}
HessianResult::Operator(op) => Err(ObjectiveEvalError::Fatal {
message: format!(
"outer plan declared HessianSource::Analytic but the runtime returned a \
non-materializable Hessian operator (dim={}, materialization={:?}); \
finite-difference Hessian estimation is not permitted on the analytic route",
op.dim(),
op.materialization_capability(),
),
}),
HessianResult::Unavailable => Err(ObjectiveEvalError::Fatal {
message: "outer plan declared HessianSource::Analytic but the runtime returned \
HessianResult::Unavailable; finite-difference Hessian estimation is \
not permitted on the analytic route"
.to_string(),
}),
},
HessianSource::BfgsApprox
| HessianSource::EfsFixedPoint
| HessianSource::HybridEfsFixedPoint => Ok(None),
}
}
struct OuterFixedPointBridge<'a> {
obj: &'a mut dyn OuterObjective,
layout: OuterThetaLayout,
barrier_config: Option<BarrierConfig>,
fixed_point_tolerance: f64,
consecutive_psi_zero_iters: usize,
}
impl OuterFixedPointBridge<'_> {
fn reject_nonstationary_tiny_psi_step(
&self,
step: &Array1<f64>,
psi_indices: Option<&[usize]>,
psi_gradient: Option<&Array1<f64>>,
cost: f64,
) -> Result<(), ObjectiveEvalError> {
let Some(psi_indices) = psi_indices else {
return Ok(());
};
let Some(psi_gradient) = psi_gradient else {
return Ok(());
};
let psi_step_inf = psi_indices
.iter()
.map(|&idx| step[idx].abs())
.fold(0.0_f64, f64::max);
let psi_grad_inf = psi_gradient.iter().map(|v| v.abs()).fold(0.0_f64, f64::max);
if psi_step_inf <= self.fixed_point_tolerance && psi_grad_inf > self.fixed_point_tolerance {
return Err(ObjectiveEvalError::recoverable(format!(
"{} HybridEFS ψ nonstationary: ||Δψ||∞={:.3e} <= tol={:.3e} \
but raw ||gψ||∞={:.3e} (rho_dim={}, psi_dim={}, n_params={}, cost={:.6e})",
EFS_FIRST_ORDER_FALLBACK_MARKER,
psi_step_inf,
self.fixed_point_tolerance,
psi_grad_inf,
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
cost,
)));
}
Ok(())
}
}
const MAX_EFS_BACKTRACK: usize = 8;
const EFS_NEGLIGIBLE_STEP: f64 = 1e-12;
const EFS_LINESEARCH_THRESHOLD: f64 = 0.5;
const EFS_COST_DESCENT_TOL: f64 = 1e-12;
const MAX_CONSECUTIVE_PSI_STAGNATION: usize = 1;
impl FixedPointObjective for OuterFixedPointBridge<'_> {
fn eval_step(&mut self, x: &Array1<f64>) -> Result<FixedPointSample, ObjectiveEvalError> {
self.layout.validate_point_len(x, "outer EFS eval failed")?;
let eval = match self.obj.eval_efs(x) {
Ok(eval) => eval,
Err(err @ EstimationError::GradientUnavailable { .. })
if requests_immediate_first_order_fallback(&err.to_string()) =>
{
log::warn!(
"[STAGE] EFS -> gradient fallback: gradient unavailable at \
fixed-point dispatch; retrying with fixed-point disabled \
(rho_dim={}, psi_dim={}, n_params={})",
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
);
return Err(ObjectiveEvalError::recoverable(format!(
"outer EFS eval failed: {err}"
)));
}
Err(err) => return Err(into_objective_error("outer EFS eval failed", err)),
};
self.layout
.validate_efs_eval(&eval, "outer EFS eval failed")?;
if !eval.cost.is_finite() {
return Err(ObjectiveEvalError::recoverable(
"outer EFS eval failed: objective returned a non-finite cost".to_string(),
));
}
if let Some((idx, value)) = eval.steps.iter().enumerate().find(|(_, v)| !v.is_finite()) {
let psi_indices = eval.psi_indices.as_deref();
let coord_kind = match psi_indices {
Some(indices) if indices.contains(&idx) => "ψ",
Some(_) => "ρ/τ",
None => "ρ",
};
return Err(ObjectiveEvalError::recoverable(format!(
"outer EFS eval failed: non-finite {coord_kind} step at coord {idx} \
(step[{idx}]={value}, rho_dim={}, psi_dim={}, n_params={}, cost={:.6e})",
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
eval.cost,
)));
}
if let Some(ref barrier_cfg) = self.barrier_config
&& let Some(ref beta) = eval.beta
{
const LOCAL_CONCENTRATION_RATIO: f64 = 0.1;
const BARRIER_CURVATURE_SATURATION: f64 = 1.0;
const BARRIER_CURVATURE_RELATIVE_THRESHOLD: f64 = 0.05;
if let Some(hessian_scale) = eval.inner_hessian_scale
&& hessian_scale.is_finite()
&& hessian_scale > 0.0
&& barrier_cfg.barrier_curvature_is_significant(
beta,
hessian_scale,
BARRIER_CURVATURE_RELATIVE_THRESHOLD,
)
{
return Err(ObjectiveEvalError::recoverable(format!(
"{} EFS barrier curvature significant relative to inner Hessian \
(rho_dim={}, psi_dim={}, n_params={}, cost={:.6e}, ref_diag={:.3e})",
EFS_FIRST_ORDER_FALLBACK_MARKER,
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
eval.cost,
hessian_scale,
)));
}
if barrier_cfg.barrier_curvature_locally_concentrated(
beta,
LOCAL_CONCENTRATION_RATIO,
BARRIER_CURVATURE_SATURATION,
) {
return Err(ObjectiveEvalError::recoverable(format!(
"{} EFS barrier curvature locally concentrated \
(rho_dim={}, psi_dim={}, n_params={}, cost={:.6e})",
EFS_FIRST_ORDER_FALLBACK_MARKER,
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
eval.cost,
)));
}
}
let status = FixedPointStatus::Continue;
let raw_step = Array1::from_vec(eval.steps);
let psi_indices = eval.psi_indices.clone();
self.reject_nonstationary_tiny_psi_step(
&raw_step,
psi_indices.as_deref(),
eval.psi_gradient.as_ref(),
eval.cost,
)?;
let max_step_abs = raw_step.iter().map(|s| s.abs()).fold(0.0_f64, f64::max);
let current_cost = eval.cost;
if self.fixed_point_step_converged(x, &raw_step, psi_indices.as_deref()) {
if psi_indices.is_some() {
self.consecutive_psi_zero_iters = 0;
}
return Ok(FixedPointSample {
value: current_cost,
step: raw_step,
status: FixedPointStatus::Stop,
});
}
if max_step_abs < EFS_NEGLIGIBLE_STEP {
if psi_indices.is_some() {
self.consecutive_psi_zero_iters = 0;
}
return Ok(FixedPointSample {
value: current_cost,
step: raw_step,
status,
});
}
if self.barrier_config.is_none() && max_step_abs < EFS_LINESEARCH_THRESHOLD {
if psi_indices.is_some() {
self.consecutive_psi_zero_iters = 0;
}
return Ok(FixedPointSample {
value: current_cost,
step: raw_step,
status,
});
}
if let Some(scaled) = self.efs_backtrack(x, &raw_step, current_cost, MAX_EFS_BACKTRACK)? {
if psi_indices.is_some() {
self.consecutive_psi_zero_iters = 0;
}
return Ok(FixedPointSample {
value: current_cost,
step: scaled,
status,
});
}
if let Some(psi_idx) = psi_indices.as_ref() {
let mut rho_only = raw_step.clone();
for &i in psi_idx {
rho_only[i] = 0.0;
}
let max_rho_abs = rho_only.iter().map(|s| s.abs()).fold(0.0_f64, f64::max);
if max_rho_abs >= EFS_NEGLIGIBLE_STEP
&& let Some(scaled) =
self.efs_backtrack(x, &rho_only, current_cost, MAX_EFS_BACKTRACK)?
{
self.consecutive_psi_zero_iters = self.consecutive_psi_zero_iters.saturating_add(1);
log::info!(
"[HYBRID-EFS] full-vector backtrack exhausted; ρ/τ-only step \
accepted. Consecutive ψ-zero iters = {}",
self.consecutive_psi_zero_iters,
);
if self.consecutive_psi_zero_iters >= MAX_CONSECUTIVE_PSI_STAGNATION {
log::info!(
"[STAGE] HybridEFS -> joint gradient (BFGS/L-BFGS) fallback: \
{} consecutive ψ-zero iterations after exhausted backtracking \
(rho_dim={}, psi_dim={}, n_params={}, cost={:.6e})",
self.consecutive_psi_zero_iters,
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
current_cost,
);
return Err(ObjectiveEvalError::recoverable(format!(
"{} HybridEFS ψ stagnation: {} consecutive iterations \
exhausted backtracking and zeroed ψ step \
(rho_dim={}, psi_dim={}, n_params={}, cost={:.6e})",
EFS_FIRST_ORDER_FALLBACK_MARKER,
self.consecutive_psi_zero_iters,
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
current_cost,
)));
}
return Ok(FixedPointSample {
value: current_cost,
step: scaled,
status,
});
}
log::info!(
"[STAGE] HybridEFS -> joint gradient fallback: ρ/τ-only step also \
failed all {} halvings (rho_dim={}, psi_dim={}, n_params={}, \
cost={:.6e})",
MAX_EFS_BACKTRACK,
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
current_cost,
);
return Err(ObjectiveEvalError::recoverable(format!(
"{} HybridEFS step rejected after {} halvings on full vector \
and {} halvings on ρ/τ-only fallback \
(rho_dim={}, psi_dim={}, n_params={}, cost={:.6e})",
EFS_FIRST_ORDER_FALLBACK_MARKER,
MAX_EFS_BACKTRACK,
MAX_EFS_BACKTRACK,
self.layout.rho_dim(),
self.layout.psi_dim,
self.layout.n_params,
current_cost,
)));
}
log::info!(
"[STAGE] EFS -> gradient fallback: no α ∈ {{1, …, 2^-{}}} decreased the \
cost (rho_dim={}, n_params={}, cost={:.6e})",
MAX_EFS_BACKTRACK,
self.layout.rho_dim(),
self.layout.n_params,
current_cost,
);
Err(ObjectiveEvalError::recoverable(format!(
"{} EFS step rejected after {} halvings on pure-ρ vector \
(rho_dim={}, n_params={}, cost={:.6e})",
EFS_FIRST_ORDER_FALLBACK_MARKER,
MAX_EFS_BACKTRACK,
self.layout.rho_dim(),
self.layout.n_params,
current_cost,
)))
}
}
impl OuterFixedPointBridge<'_> {
fn efs_backtrack(
&mut self,
x: &Array1<f64>,
raw_step: &Array1<f64>,
current_cost: f64,
max_halvings: usize,
) -> Result<Option<Array1<f64>>, ObjectiveEvalError> {
let cost_floor = current_cost + EFS_COST_DESCENT_TOL * current_cost.abs().max(1.0);
let mut alpha = 1.0_f64;
for bt in 0..=max_halvings {
let trial_step = raw_step * alpha;
let trial = x + &trial_step;
match self.obj.eval_cost(&trial) {
Ok(c) if c.is_finite() && c <= cost_floor => {
if bt > 0 {
log::debug!(
"[EFS] backtrack accepted at α=2^-{bt}={alpha:.4e} \
after {bt} halvings (cost: {current_cost:.6e} → {c:.6e})"
);
}
return Ok(Some(trial_step));
}
Ok(c) => {
log::trace!(
"[EFS] backtrack α=2^-{bt}={alpha:.4e}: trial cost {c:.6e} \
not below current {current_cost:.6e}, halving"
);
}
Err(err) => {
log::trace!(
"[EFS] backtrack α=2^-{bt}={alpha:.4e}: trial eval failed \
({err}), halving"
);
}
}
alpha *= 0.5;
}
Ok(None)
}
fn fixed_point_step_converged(
&self,
x: &Array1<f64>,
step: &Array1<f64>,
psi_indices: Option<&[usize]>,
) -> bool {
if x.len() != step.len() {
return false;
}
for idx in 0..step.len() {
let scale = match psi_indices {
Some(indices) if indices.contains(&idx) => x[idx].abs().max(1.0),
_ => 1.0,
};
let normalized = step[idx].abs() / scale;
if !normalized.is_finite() || normalized > self.fixed_point_tolerance {
return false;
}
}
true
}
}
enum CompassSearchOutcome {
Converged {
point: Array1<f64>,
cost: f64,
polls: usize,
},
BudgetExhausted {
point: Array1<f64>,
cost: f64,
polls: usize,
},
}
fn compass_search_outer(
obj: &mut dyn OuterObjective,
mut x: Array1<f64>,
mut best_cost: f64,
lower: ndarray::ArrayView1<'_, f64>,
upper: ndarray::ArrayView1<'_, f64>,
init_step: f64,
step_tol: f64,
max_polls: usize,
) -> CompassSearchOutcome {
for i in 0..x.len() {
x[i] = x[i].clamp(lower[i], upper[i]);
}
let mut step = init_step;
let mut polls: usize = 0;
while step > step_tol && polls < max_polls {
let mut improved = false;
'sweep: for i in 0..x.len() {
for &sign in &[1.0, -1.0] {
if polls >= max_polls {
break 'sweep;
}
polls += 1;
let candidate_i = (x[i] + sign * step).clamp(lower[i], upper[i]);
if (candidate_i - x[i]).abs() < step_tol {
continue;
}
let mut candidate = x.clone();
candidate[i] = candidate_i;
let probe = obj.eval_cost(&candidate).ok().filter(|v| v.is_finite());
if let Some(c) = probe
&& c < best_cost
{
x = candidate;
best_cost = c;
improved = true;
break 'sweep;
}
}
}
if !improved {
step *= 0.5;
}
}
if step <= step_tol {
CompassSearchOutcome::Converged {
point: x,
cost: best_cost,
polls,
}
} else {
CompassSearchOutcome::BudgetExhausted {
point: x,
cost: best_cost,
polls,
}
}
}
fn solution_into_outer_result(
solution: Solution,
converged: bool,
plan_used: OuterPlan,
) -> OuterResult {
OuterResult {
rho: solution.final_point,
final_value: solution.final_value,
iterations: solution.iterations,
final_grad_norm: solution.final_gradient_norm,
final_gradient: solution.final_gradient,
final_hessian: solution.final_hessian,
converged,
plan_used,
operator_trust_radius: None,
operator_stop_reason: None,
}
}
use crate::inference::diagnostics::format_top_abs as format_top_abs_components;
fn bfgs_line_search_failure_message(
context: &str,
solution: &Solution,
max_attempts: usize,
failure_reason: impl std::fmt::Debug,
) -> String {
let grad_norm = solution
.final_gradient_norm
.or_else(|| {
solution
.final_gradient
.as_ref()
.map(|gradient| gradient.iter().map(|v| v * v).sum::<f64>().sqrt())
})
.unwrap_or(f64::NAN);
let gradient_detail = solution
.final_gradient
.as_ref()
.map(|gradient| format_top_abs_components(gradient, "top_abs_gradient", 6))
.unwrap_or_else(|| "top_abs_gradient=<unavailable>".to_string());
format!(
"{context}: BFGS line search failed; reason={failure_reason:?} \
max_attempts={max_attempts} iterations={} final_value={:.6e} \
|g|={:.3e} func_evals={} grad_evals={} {} {}",
solution.iterations,
solution.final_value,
grad_norm,
solution.func_evals,
solution.grad_evals,
format_top_abs_components(&solution.final_point, "top_abs_rho", 6),
gradient_detail,
)
}
#[derive(Clone, Debug)]
struct OuterConfig {
tolerance: f64,
max_iter: usize,
bounds: Option<(Array1<f64>, Array1<f64>)>,
seed_config: crate::seeding::SeedConfig,
rho_bound: f64,
heuristic_lambdas: Option<Vec<f64>>,
initial_rho: Option<Array1<f64>>,
fallback_policy: FallbackPolicy,
screening_cap: Option<Arc<AtomicUsize>>,
screen_initial_rho: bool,
outer_inner_cap: Option<InnerProgressFeedback>,
solver_class: SolverClass,
operator_initial_trust_radius: Option<f64>,
arc_initial_regularization: Option<f64>,
objective_scale: Option<f64>,
bfgs_step_cap: Option<f64>,
bfgs_step_cap_psi: Option<f64>,
cache_session: Option<Arc<CacheSession>>,
cache_mirror_sessions: Vec<Arc<CacheSession>>,
}
impl Default for OuterConfig {
fn default() -> Self {
Self {
tolerance: 1e-5,
max_iter: 200,
bounds: None,
seed_config: crate::seeding::SeedConfig::default(),
rho_bound: 30.0,
heuristic_lambdas: None,
initial_rho: None,
fallback_policy: FallbackPolicy::Automatic,
screening_cap: None,
screen_initial_rho: false,
outer_inner_cap: None,
solver_class: SolverClass::Primary,
operator_initial_trust_radius: None,
arc_initial_regularization: None,
objective_scale: None,
bfgs_step_cap: None,
bfgs_step_cap_psi: None,
cache_session: None,
cache_mirror_sessions: Vec::new(),
}
}
}
pub struct OuterProblem {
n_params: usize,
gradient: Derivative,
hessian: DeclaredHessianForm,
prefer_gradient_only: bool,
disable_fixed_point: bool,
psi_dim: usize,
barrier_config: Option<BarrierConfig>,
tolerance: f64,
max_iter: usize,
bounds: Option<(Array1<f64>, Array1<f64>)>,
rho_bound: f64,
seed_config: crate::seeding::SeedConfig,
heuristic_lambdas: Option<Vec<f64>>,
initial_rho: Option<Array1<f64>>,
fallback_policy: FallbackPolicy,
screening_cap: Option<Arc<AtomicUsize>>,
screen_initial_rho: bool,
outer_inner_cap: Option<InnerProgressFeedback>,
solver_class: SolverClass,
operator_initial_trust_radius: Option<f64>,
arc_initial_regularization: Option<f64>,
objective_scale: Option<f64>,
bfgs_step_cap: Option<f64>,
bfgs_step_cap_psi: Option<f64>,
cache_session: Option<Arc<CacheSession>>,
cache_mirror_sessions: Vec<Arc<CacheSession>>,
}
impl OuterProblem {
pub fn new(n_params: usize) -> Self {
Self {
n_params,
gradient: Derivative::Unavailable,
hessian: DeclaredHessianForm::Unavailable,
prefer_gradient_only: false,
disable_fixed_point: false,
psi_dim: 0,
barrier_config: None,
tolerance: 1e-5,
max_iter: 200,
bounds: None,
rho_bound: 30.0,
seed_config: crate::seeding::SeedConfig::default(),
heuristic_lambdas: None,
initial_rho: None,
fallback_policy: FallbackPolicy::Automatic,
screening_cap: None,
screen_initial_rho: false,
outer_inner_cap: None,
solver_class: SolverClass::Primary,
operator_initial_trust_radius: None,
arc_initial_regularization: None,
objective_scale: None,
bfgs_step_cap: None,
bfgs_step_cap_psi: None,
cache_session: None,
cache_mirror_sessions: Vec::new(),
}
}
pub fn with_gradient(mut self, d: Derivative) -> Self {
self.gradient = d;
self
}
pub fn with_hessian(mut self, form: DeclaredHessianForm) -> Self {
self.hessian = form;
self
}
pub fn with_prefer_gradient_only(mut self, prefer_gradient_only: bool) -> Self {
self.prefer_gradient_only = prefer_gradient_only;
self
}
pub fn with_disable_fixed_point(mut self, disable: bool) -> Self {
self.disable_fixed_point = disable;
self
}
pub fn with_psi_dim(mut self, dim: usize) -> Self {
self.psi_dim = dim;
self
}
pub fn with_barrier(mut self, cfg: Option<BarrierConfig>) -> Self {
self.barrier_config = cfg;
self
}
pub fn with_tolerance(mut self, tol: f64) -> Self {
self.tolerance = tol;
self
}
pub fn with_max_iter(mut self, n: usize) -> Self {
self.max_iter = n;
self
}
pub fn with_bounds(mut self, lo: Array1<f64>, hi: Array1<f64>) -> Self {
self.bounds = Some((lo, hi));
self
}
pub fn with_rho_bound(mut self, b: f64) -> Self {
self.rho_bound = b;
self
}
pub fn with_seed_config(mut self, sc: crate::seeding::SeedConfig) -> Self {
self.seed_config = sc;
self
}
pub fn with_heuristic_lambdas(mut self, h: Vec<f64>) -> Self {
self.heuristic_lambdas = Some(h);
self
}
pub fn with_initial_rho(mut self, rho: Array1<f64>) -> Self {
self.initial_rho = Some(rho);
self
}
pub fn with_screening_cap(mut self, screening_cap: Arc<AtomicUsize>) -> Self {
self.screening_cap = Some(screening_cap);
self
}
pub fn with_screen_initial_rho(mut self, screen_initial_rho: bool) -> Self {
self.screen_initial_rho = screen_initial_rho;
self
}
pub fn with_outer_inner_cap(mut self, feedback: InnerProgressFeedback) -> Self {
self.outer_inner_cap = Some(feedback);
self
}
pub fn with_solver_class(mut self, class: SolverClass) -> Self {
self.solver_class = class;
self
}
pub fn with_operator_initial_trust_radius(mut self, radius: Option<f64>) -> Self {
self.operator_initial_trust_radius = sanitized_operator_trust_restart_radius(radius);
self
}
pub fn with_arc_initial_regularization(mut self, sigma: Option<f64>) -> Self {
self.arc_initial_regularization = sigma.filter(|v| v.is_finite() && *v > 0.0);
self
}
pub fn with_objective_scale(mut self, scale: Option<f64>) -> Self {
self.objective_scale = scale.filter(|v| v.is_finite() && *v > 0.0);
self
}
pub fn with_bfgs_step_cap(mut self, cap: Option<f64>) -> Self {
self.bfgs_step_cap = cap.filter(|v| v.is_finite() && *v > 0.0);
self
}
pub fn with_bfgs_step_cap_psi(mut self, cap: Option<f64>) -> Self {
self.bfgs_step_cap_psi = cap.filter(|v| v.is_finite() && *v > 0.0);
self
}
pub fn with_cache_session(mut self, session: Arc<CacheSession>) -> Self {
self.cache_session = Some(session);
self
}
pub fn with_cache_mirror_sessions(mut self, sessions: Vec<Arc<CacheSession>>) -> Self {
self.cache_mirror_sessions = sessions;
self
}
pub fn with_fallback_policy(mut self, policy: FallbackPolicy) -> Self {
self.fallback_policy = policy;
self
}
fn capability(&self) -> OuterCapability {
OuterCapability {
gradient: self.gradient,
hessian: self.hessian,
prefer_gradient_only: self.prefer_gradient_only,
disable_fixed_point: self.disable_fixed_point,
n_params: self.n_params,
psi_dim: self.psi_dim,
fixed_point_available: false,
barrier_config: self.barrier_config.clone(),
}
}
fn config(&self) -> OuterConfig {
OuterConfig {
tolerance: self.tolerance,
max_iter: self.max_iter,
bounds: self.bounds.clone(),
seed_config: self.seed_config,
rho_bound: self.rho_bound,
heuristic_lambdas: self.heuristic_lambdas.clone(),
initial_rho: self.initial_rho.clone(),
fallback_policy: self.fallback_policy,
screening_cap: self.screening_cap.clone(),
screen_initial_rho: self.screen_initial_rho,
outer_inner_cap: self.outer_inner_cap.clone(),
solver_class: self.solver_class,
operator_initial_trust_radius: self.operator_initial_trust_radius,
arc_initial_regularization: self.arc_initial_regularization,
objective_scale: self.objective_scale,
bfgs_step_cap: self.bfgs_step_cap,
bfgs_step_cap_psi: self.bfgs_step_cap_psi,
cache_session: self.cache_session.clone(),
cache_mirror_sessions: self.cache_mirror_sessions.clone(),
}
}
pub fn build_objective<S, Fc, Fe, Fr, Fefs>(
&self,
state: S,
cost_fn: Fc,
eval_fn: Fe,
reset_fn: Option<Fr>,
efs_fn: Option<Fefs>,
) -> ClosureObjective<S, Fc, Fe, Fr, Fefs>
where
Fc: FnMut(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
Fe: FnMut(&mut S, &Array1<f64>) -> Result<OuterEval, EstimationError>,
Fr: FnMut(&mut S),
Fefs: FnMut(&mut S, &Array1<f64>) -> Result<EfsEval, EstimationError>,
{
let mut cap = self.capability();
cap.fixed_point_available = efs_fn.is_some();
ClosureObjective {
state,
cap,
cost_fn,
eval_fn,
eval_order_fn: None,
reset_fn,
efs_fn,
screening_proxy_fn: None::<fn(&mut S, &Array1<f64>) -> Result<f64, EstimationError>>,
seed_fn: None::<fn(&mut S, &Array1<f64>) -> Result<(), EstimationError>>,
}
}
pub fn build_objective_with_eval_order<S, Fc, Fe, Feo, Fr, Fefs>(
&self,
state: S,
cost_fn: Fc,
eval_fn: Fe,
eval_order_fn: Feo,
reset_fn: Option<Fr>,
efs_fn: Option<Fefs>,
) -> ClosureObjective<S, Fc, Fe, Fr, Fefs, Feo>
where
Fc: FnMut(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
Fe: FnMut(&mut S, &Array1<f64>) -> Result<OuterEval, EstimationError>,
Feo: FnMut(&mut S, &Array1<f64>, OuterEvalOrder) -> Result<OuterEval, EstimationError>,
Fr: FnMut(&mut S),
Fefs: FnMut(&mut S, &Array1<f64>) -> Result<EfsEval, EstimationError>,
{
let mut cap = self.capability();
cap.fixed_point_available = efs_fn.is_some();
ClosureObjective {
state,
cap,
cost_fn,
eval_fn,
eval_order_fn: Some(eval_order_fn),
reset_fn,
efs_fn,
screening_proxy_fn: None::<fn(&mut S, &Array1<f64>) -> Result<f64, EstimationError>>,
seed_fn: None::<fn(&mut S, &Array1<f64>) -> Result<(), EstimationError>>,
}
}
pub fn build_objective_with_screening_proxy<S, Fc, Fe, Feo, Fr, Fefs, Fsp>(
&self,
state: S,
cost_fn: Fc,
eval_fn: Fe,
eval_order_fn: Feo,
reset_fn: Option<Fr>,
efs_fn: Option<Fefs>,
screening_proxy_fn: Fsp,
) -> ClosureObjective<S, Fc, Fe, Fr, Fefs, Feo, Fsp>
where
Fc: FnMut(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
Fe: FnMut(&mut S, &Array1<f64>) -> Result<OuterEval, EstimationError>,
Feo: FnMut(&mut S, &Array1<f64>, OuterEvalOrder) -> Result<OuterEval, EstimationError>,
Fr: FnMut(&mut S),
Fefs: FnMut(&mut S, &Array1<f64>) -> Result<EfsEval, EstimationError>,
Fsp: FnMut(&mut S, &Array1<f64>) -> Result<f64, EstimationError>,
{
let mut cap = self.capability();
cap.fixed_point_available = efs_fn.is_some();
ClosureObjective {
state,
cap,
cost_fn,
eval_fn,
eval_order_fn: Some(eval_order_fn),
reset_fn,
efs_fn,
screening_proxy_fn: Some(screening_proxy_fn),
seed_fn: None::<fn(&mut S, &Array1<f64>) -> Result<(), EstimationError>>,
}
}
pub fn run(
&self,
obj: &mut dyn OuterObjective,
context: &str,
) -> Result<OuterResult, EstimationError> {
let mut config = self.config();
let Some(session) = config.cache_session.clone() else {
return run_outer(obj, &config, context);
};
let key_hex = session.key().to_hex();
let short_key = &key_hex[..8.min(key_hex.len())];
let mut had_hit = false;
let mut cached_inner_beta: Option<Array1<f64>> = None;
if let Some(loaded) = session.try_load_with_source() {
match classify_cache_entry_for_outer(&loaded, self.n_params) {
CacheSeedDecision::ExactFinal {
rho,
beta: _beta_final,
final_value,
iterations,
prior_obj_display,
} => {
let cap = primary_capability_for_config(obj.capability(), &config, context);
let plan_used = plan_with_class(&cap, config.solver_class);
log::info!(
"[CACHE] final-hit key={}.. context={} rho_dim={} prior_obj={:.6e} iter={} action=skip-outer-validation",
short_key,
context,
rho.len(),
prior_obj_display,
iterations,
);
return Ok(OuterResult {
rho,
final_value,
iterations,
final_grad_norm: None,
final_gradient: None,
final_hessian: None,
converged: true,
plan_used,
operator_trust_radius: None,
operator_stop_reason: None,
});
}
CacheSeedDecision::Seed {
rho,
beta,
prior_obj_display,
iteration,
} => {
let beta_len = beta.len();
let beta_arr = if beta.is_empty() {
None
} else {
Some(Array1::from_vec(beta))
};
if config
.initial_rho
.as_ref()
.is_none_or(|initial| initial != rho)
{
log::info!(
"[CACHE] hit key={}.. context={} rho_dim={} beta_dim={} prior_obj={:.6e} iter={}",
short_key,
context,
rho.len(),
beta_len,
prior_obj_display,
iteration,
);
config.initial_rho = Some(rho);
config.screen_initial_rho = false;
had_hit = true;
} else {
log::info!(
"[CACHE] hit key={}.. context={} rho_dim={} beta_dim={} already-aligned prior_obj={:.6e}",
short_key,
context,
rho.len(),
beta_len,
prior_obj_display,
);
had_hit = true;
}
cached_inner_beta = beta_arr;
}
CacheSeedDecision::Discard {
reason: "payload-shape-mismatch",
..
} => {
log::info!(
"[CACHE] skip key={}.. context={} reason=payload-shape-mismatch n_params={}",
short_key,
context,
self.n_params,
);
}
CacheSeedDecision::Discard {
reason,
prior_obj_display,
all_rho_finite,
} => {
log::info!(
"[CACHE] skip key={}.. context={} reason={} prior_obj={:.6e} all_rho_finite={}",
short_key,
context,
reason,
prior_obj_display,
all_rho_finite.unwrap_or(false),
);
}
}
} else {
log::info!(
"[CACHE] miss key={}.. context={} reason=fresh-fingerprint n_params={}",
short_key,
context,
self.n_params,
);
}
let mut checkpointing = CheckpointingObjective::new(
obj,
Arc::clone(&session),
config.cache_mirror_sessions.clone(),
);
if let Some(beta) = cached_inner_beta.as_ref() {
match checkpointing.seed_inner_state(beta) {
Ok(()) => log::info!(
"[CACHE] beta-warm key={}.. context={} beta_dim={} action=installed",
short_key,
context,
beta.len(),
),
Err(err) => log::warn!(
"[CACHE] beta-warm key={}.. context={} beta_dim={} action=skip err={}",
short_key,
context,
beta.len(),
err,
),
}
}
let result = run_outer(&mut checkpointing, &config, context);
let final_beta = checkpointing.last_inner_beta();
if let Ok(result) = result.as_ref()
&& result.final_value.is_finite()
&& result.converged
&& let Some(bytes) = encode_iterate(
&result.rho,
final_beta.as_ref(),
result.final_value,
result.iterations as u64,
)
{
let saved = session.finalize(
&bytes,
Some(result.final_value),
Some(result.iterations as u64),
);
if saved {
log::info!(
"[CACHE] save key={}.. context={} final_obj={:.6e} iter={} resumed={}",
short_key,
context,
result.final_value,
result.iterations,
had_hit,
);
}
for mirror in &config.cache_mirror_sessions {
let mirror_saved = mirror.finalize(
&bytes,
Some(result.final_value),
Some(result.iterations as u64),
);
if mirror_saved {
let mirror_hex = mirror.key().to_hex();
log::info!(
"[CACHE] save key={}.. context={} mirror final_obj={:.6e} iter={}",
&mirror_hex[..8.min(mirror_hex.len())],
context,
result.final_value,
result.iterations,
);
}
}
}
result
}
}
#[derive(Clone, Debug)]
pub struct OuterResult {
pub rho: Array1<f64>,
pub final_value: f64,
pub iterations: usize,
pub final_grad_norm: Option<f64>,
pub final_gradient: Option<Array1<f64>>,
pub final_hessian: Option<Array2<f64>>,
pub converged: bool,
pub plan_used: OuterPlan,
pub operator_trust_radius: Option<f64>,
pub operator_stop_reason: Option<OperatorTrustRegionStopReason>,
}
impl OuterResult {
pub fn final_grad_norm_report(&self) -> String {
match self.final_grad_norm {
Some(g) => format!("{g:.3e}"),
None => "n/a".to_string(),
}
}
}
#[derive(Clone, Copy, Debug, PartialEq, Eq)]
pub enum OperatorTrustRegionStopReason {
Converged,
RejectFloor,
IterationBudget,
RoutingMismatch,
}
fn run_outer(
obj: &mut dyn OuterObjective,
config: &OuterConfig,
context: &str,
) -> Result<OuterResult, EstimationError> {
let cap = primary_capability_for_config(obj.capability(), config, context);
cap.validate_layout(context)?;
if let Some(initial_rho) = config.initial_rho.as_ref() {
cap.theta_layout()
.validate_point_len(initial_rho, "initial outer seed")
.map_err(|err| match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(format!("{context}: {message}"))
}
})?;
}
crate::solver::estimate::reml::runtime::clear_outer_ift_residual_energy_for_fit();
if cap.n_params == 0 {
let cost = obj.eval_cost(&Array1::zeros(0))?;
let the_plan = plan_with_class(&cap, config.solver_class);
return Ok(OuterResult {
rho: Array1::zeros(0),
final_value: cost,
iterations: 0,
final_grad_norm: Some(0.0),
final_gradient: None,
final_hessian: None,
converged: true,
plan_used: the_plan,
operator_trust_radius: None,
operator_stop_reason: None,
});
}
let fallback_attempts = match (config.fallback_policy, config.solver_class) {
(FallbackPolicy::Automatic, SolverClass::Primary) => automatic_fallback_attempts(&cap),
(FallbackPolicy::Automatic, SolverClass::AuxiliaryGradientFree)
| (FallbackPolicy::Disabled, _) => Vec::new(),
};
let mut attempts: Vec<OuterCapability> = Vec::with_capacity(1 + fallback_attempts.len());
attempts.push(cap.clone());
for degraded in fallback_attempts {
attempts.push(degraded);
}
let mut last_error: Option<EstimationError> = None;
for (attempt_idx, attempt_cap) in attempts.iter().enumerate() {
let the_plan = plan_with_class(attempt_cap, config.solver_class);
if attempt_idx > 0 {
log::debug!("[OUTER] {context}: primary plan failed; falling back to {the_plan}");
}
log_plan(context, attempt_cap, &the_plan);
obj.reset();
let mut arc_retries_left: u32 = if matches!(the_plan.solver, Solver::Arc) {
2
} else {
0
};
let mut retry_config: Option<OuterConfig> = None;
let mut prev_attempt_grad_norm: Option<f64> = None;
let outcome = loop {
let active_config_owned: OuterConfig =
retry_config.clone().unwrap_or_else(|| config.clone());
let active_config: &OuterConfig = &active_config_owned;
match run_outer_with_plan(obj, active_config, context, attempt_cap, &the_plan) {
Ok(result) => {
if result.converged
|| arc_retries_left == 0
|| matches!(
result.operator_stop_reason,
Some(OperatorTrustRegionStopReason::RejectFloor)
)
{
break Ok(result);
}
let Some(cur_grad_norm) = result.final_grad_norm else {
log::info!(
"[OUTER] {context}: ARC attempt exhausted budget at \
iter={} cost={:.6e} without a final gradient norm; \
falling through to degraded plan",
result.iterations,
result.final_value,
);
break Ok(result);
};
if let Some(prev_g) = prev_attempt_grad_norm {
let progressed = cur_grad_norm.is_finite()
&& prev_g.is_finite()
&& cur_grad_norm < 0.5 * prev_g;
if !progressed {
log::info!(
"[OUTER] {context}: ARC retry stalled at \
iter={} cost={:.6e} |g|={:.6e} (prev |g|={:.6e}); \
deterministic replay suspected, falling through \
to degraded plan",
result.iterations,
result.final_value,
cur_grad_norm,
prev_g,
);
break Ok(result);
}
}
let next_trust_radius =
sanitized_operator_trust_restart_radius(result.operator_trust_radius);
log::info!(
"[OUTER] {context}: ARC attempt exhausted budget at \
iter={} cost={:.6e} |g|={:.6e}; resuming from last \
rho + trust_radius={:?}, inner-PIRLS uncapped \
(objective caches wiped; operator-TR Cauchy/Newton \
state is not resumable)",
result.iterations,
result.final_value,
cur_grad_norm,
next_trust_radius,
);
let cap_feedback = active_config.outer_inner_cap.clone();
let mut next = active_config.clone();
prev_attempt_grad_norm = Some(cur_grad_norm);
next.initial_rho = Some(result.rho.clone());
next.operator_initial_trust_radius = next_trust_radius;
retry_config = Some(next);
arc_retries_left -= 1;
obj.reset();
if let Some(feedback) = cap_feedback.as_ref() {
feedback.cap.store(0, Ordering::Relaxed);
}
}
Err(e) => break Err(e),
}
};
match outcome {
Ok(result) => {
if result.converged || attempt_idx + 1 == attempts.len() {
if !result.converged {
log::warn!(
"[OUTER warning] {context}: final outer attempt returned without convergence \
(plan={the_plan}, iterations={}, final_value={:.6e}, |g|={})",
result.iterations,
result.final_value,
result.final_grad_norm_report(),
);
}
return Ok(result);
}
let message = format!(
"{context}: attempt {} (plan={the_plan}) exhausted without convergence",
attempt_idx + 1
);
log::debug!("[OUTER] {message}; trying degraded fallback plan");
last_error = Some(EstimationError::RemlOptimizationFailed(message));
}
Err(e) => {
log::debug!(
"[OUTER] {context}: attempt {} (plan={the_plan}) failed: {e}",
attempt_idx + 1
);
last_error = Some(e);
}
}
}
Err(last_error.unwrap_or_else(|| {
EstimationError::RemlOptimizationFailed(format!("all plan attempts exhausted ({context})"))
}))
}
fn outer_bounds(lo: &Array1<f64>, hi: &Array1<f64>) -> Result<Bounds, EstimationError> {
Bounds::new(lo.clone(), hi.clone(), 1e-6).map_err(|err| {
EstimationError::InvalidInput(format!("outer rho bounds are invalid: {err}"))
})
}
fn outer_tolerance(value: f64) -> Result<Tolerance, EstimationError> {
Tolerance::new(value)
.map_err(|err| EstimationError::InvalidInput(format!("outer tolerance is invalid: {err}")))
}
fn outer_gradient_tolerance(config: &OuterConfig) -> GradientTolerance {
let abs = config
.objective_scale
.map(|scale| config.tolerance.max(scale * 1.0e-9))
.unwrap_or(config.tolerance);
GradientTolerance {
abs,
rel_initial_grad: None,
rel_cost: Some(config.tolerance),
projected: true,
}
}
fn outer_max_iterations(value: usize) -> Result<MaxIterations, EstimationError> {
MaxIterations::new(value)
.map_err(|err| EstimationError::InvalidInput(format!("outer max_iter is invalid: {err}")))
}
fn sanitized_operator_trust_restart_radius(radius: Option<f64>) -> Option<f64> {
radius
.filter(|value| value.is_finite() && *value > 0.0)
.map(|value| value.max(OPERATOR_TRUST_RESTART_RADIUS_FLOOR))
}
fn bfgs_axis_step_caps(config: &OuterConfig, layout: OuterThetaLayout) -> Option<Array1<f64>> {
if config.bfgs_step_cap.is_none() && config.bfgs_step_cap_psi.is_none() {
return None;
}
let mut caps = Array1::from_elem(layout.n_params, f64::INFINITY);
if let Some(cap) = config.bfgs_step_cap {
for i in 0..layout.rho_dim() {
caps[i] = cap;
}
}
if let Some(cap) = config.bfgs_step_cap_psi {
for i in layout.rho_dim()..layout.n_params {
caps[i] = cap;
}
}
Some(caps)
}
fn run_outer_with_plan(
obj: &mut dyn OuterObjective,
config: &OuterConfig,
context: &str,
cap: &OuterCapability,
the_plan: &OuterPlan,
) -> Result<OuterResult, EstimationError> {
let mut seeds = {
let generated = crate::seeding::generate_rho_candidates(
cap.n_params,
config.heuristic_lambdas.as_deref(),
&config.seed_config,
);
if generated.is_empty() {
Vec::new()
} else {
generated
}
};
if let Some(initial_rho) = config.initial_rho.as_ref()
&& !seeds.iter().any(|seed| seed == initial_rho)
{
seeds.insert(0, initial_rho.clone());
}
if seeds.is_empty() {
return Err(EstimationError::RemlOptimizationFailed(format!(
"no seeds generated for outer optimization ({context})"
)));
}
let (lower, upper) = config.bounds.clone().unwrap_or_else(|| {
(
Array1::<f64>::from_elem(cap.n_params, -config.rho_bound),
Array1::<f64>::from_elem(cap.n_params, config.rho_bound),
)
});
crate::solver::estimate::reml::runtime::record_current_outer_rho_upper_bounds_for_ift(&upper);
let bounds_template = (lower, upper);
let mut projected_seeds = Vec::with_capacity(seeds.len());
for seed in seeds {
let projected = project_to_bounds(&seed, Some(&bounds_template));
if !projected_seeds.contains(&projected) {
projected_seeds.push(projected);
}
}
seeds = projected_seeds;
if seeds.is_empty() {
return Err(EstimationError::RemlOptimizationFailed(format!(
"no bounded seeds generated for outer optimization ({context})"
)));
}
let screening_enabled = config.screening_cap.is_some();
let seed_budget = effective_seed_budget(
config.seed_config.seed_budget,
the_plan.solver,
config.seed_config.risk_profile,
screening_enabled,
)
.min(seeds.len());
let explicit_initial_rho_owns_single_seed_budget = config.initial_rho.is_some()
&& seed_budget == 1
&& seeds.len() > 1
&& !config.screen_initial_rho;
if !explicit_initial_rho_owns_single_seed_budget
&& should_screen_seeds(config, the_plan.solver, seeds.len(), seed_budget)
{
seeds = rank_seeds_with_screening(obj, config, context, &seeds);
}
log::debug!(
"[OUTER] {context}: trying generated seeds directly (generated={}, budget={})",
seeds.len(),
seed_budget,
);
if seed_budget < config.seed_config.seed_budget.max(1) {
log::debug!(
"[OUTER] {context}: capped requested seed budget {} -> {} for {:?} ({:?})",
config.seed_config.seed_budget.max(1),
seed_budget,
the_plan.solver,
config.seed_config.risk_profile,
);
}
if seeds.len() > seed_budget {
log::debug!(
"[OUTER] {context}: trying up to {seed_budget}/{} generated seeds in heuristic order",
seeds.len(),
);
}
let mut best: Option<OuterResult> = None;
let mut rejection_reasons: Vec<(usize, &'static str, String)> = Vec::new();
let layout = cap.theta_layout();
let mut started_seeds = 0usize;
let expensive_seed_limit =
expensive_unsuccessful_seed_limit(the_plan.solver, config.seed_config.risk_profile);
let mut unsuccessful_expensive_seeds = 0usize;
let mut stopped_early_due_to_limit = false;
let mut seed_rejections: Vec<SeedRejection> = Vec::new();
let mut last_classified_reason_idx: usize = 0;
let mut structural_early_exit_key: Option<(
crate::families::custom_family::KktRefusalDiagnosis,
Option<String>,
)> = None;
const STRUCTURAL_EARLY_EXIT_MIN_COUNT: usize = 2;
'seed_attempts: for (seed_idx, seed) in seeds.iter().enumerate() {
if started_seeds == seed_budget {
break;
}
while last_classified_reason_idx < rejection_reasons.len() {
let (idx, phase, msg) = &rejection_reasons[last_classified_reason_idx];
seed_rejections.push(SeedRejection::from_message(*idx, phase, msg.clone()));
last_classified_reason_idx += 1;
}
if structural_early_exit_key.is_none() {
if let Some(key) =
uniform_structural_key(&seed_rejections, STRUCTURAL_EARLY_EXIT_MIN_COUNT)
{
log::warn!(
"[OUTER] {context}: structural early-exit after {} uniform CertRefused \
rejections (diagnosis={}, carrying-block={}); skipping remaining {} seed(s)",
seed_rejections.len(),
key.0.as_str(),
key.1.as_deref().unwrap_or("<unknown>"),
seeds.len().saturating_sub(seed_idx),
);
structural_early_exit_key = Some(key);
break;
}
}
crate::solver::estimate::reml::runtime::record_current_outer_iter_for_ift(0);
obj.reset();
let prewarm_start = std::time::Instant::now();
match crate::solver::estimate::reml::continuation::prime_outer_seed(
obj,
seed,
&bounds_template.1,
) {
Ok(summary) => {
if !summary.collapsed {
log::info!(
"[OUTER] {context}: continuation pre-warm seed {seed_idx} steps={} elapsed={:.3}s",
summary.steps_accepted,
prewarm_start.elapsed().as_secs_f64(),
);
}
}
Err(cf) => {
let msg = format!(
"continuation pre-warm refused before seed eval: {}",
cf.message()
);
log::warn!("[OUTER] {context}: rejecting seed {seed_idx} (continuation): {msg}");
rejection_reasons.push((seed_idx, "validation", msg));
continue 'seed_attempts;
}
}
let t_seed_start = std::time::Instant::now();
let seed_slot;
let result: Result<OuterResult, EstimationError> = match the_plan.solver {
Solver::Arc => {
let seed_eval = obj
.eval_with_order(seed, OuterEvalOrder::ValueGradientHessian)
.map_err(|err| into_objective_error("outer eval failed", err));
let seed_eval = match seed_eval {
Ok(seed_eval) => seed_eval,
Err(err) => {
let err = match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(message)
}
};
if requests_immediate_first_order_fallback(&err.to_string()) {
return Err(err);
}
log::warn!(
"[OUTER] {context}: rejecting seed {seed_idx} before solver start: {err}"
);
rejection_reasons.push((seed_idx, "validation", err.to_string()));
continue 'seed_attempts;
}
};
let seed_eval = finite_outer_eval_or_error("outer eval failed", layout, seed_eval)
.map_err(|err| match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(message)
}
});
let mut seed_eval = match seed_eval {
Ok(seed_eval) => seed_eval,
Err(err) => {
log::warn!(
"[OUTER] {context}: rejecting seed {seed_idx} before solver start: {err}"
);
rejection_reasons.push((seed_idx, "validation", err.to_string()));
continue 'seed_attempts;
}
};
validate_second_order_seed_hessian(context, layout, &seed_eval).map_err(|err| {
match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(message)
}
}
})?;
started_seeds += 1;
seed_slot = started_seeds;
let cheap_materializable_operator = matches!(
seed_eval.hessian,
HessianResult::Operator(ref op)
if op.materialization_capability().is_available()
&& op.dim() <= OUTER_HVP_MATERIALIZE_MAX_DIM
);
if cheap_materializable_operator {
if let HessianResult::Operator(op) = &seed_eval.hessian {
match op.materialize_dense() {
Ok(dense) => {
seed_eval.hessian = HessianResult::Analytic(dense);
}
Err(message) => {
let err = EstimationError::RemlOptimizationFailed(format!(
"outer Hessian operator materialization failed: {message}"
));
log::warn!(
"[OUTER] {context}: rejecting seed {seed_idx} before solver start: {err}"
);
rejection_reasons.push((seed_idx, "validation", err.to_string()));
continue 'seed_attempts;
}
}
}
}
if matches!(seed_eval.hessian, HessianResult::Operator(_)) {
log::debug!(
"[OUTER] {context}: analytic Hessian provided as Hv operator; \
routing to opt::MatrixFreeTrustRegion (Steihaug-Toint CG)"
);
let (lo, hi) = &bounds_template;
let bounds_obj = outer_bounds(lo, hi)?;
let grad_tol = outer_gradient_tolerance(config);
let max_iter = outer_max_iterations(config.max_iter)?;
let initial_op_sample = OperatorSample {
value: seed_eval.cost,
gradient: seed_eval.gradient.clone(),
hessian: hessian_result_to_value(seed_eval.hessian.clone()),
};
let bridge_obj = OuterOperatorBridge {
obj,
layout,
outer_inner_cap: config.outer_inner_cap.clone(),
eval_count: 0,
g_norm_initial: None,
last_g_norm: None,
};
let mut solver = MatrixFreeTrustRegion::new(seed.clone(), bridge_obj)
.with_bounds(bounds_obj)
.with_gradient_tolerance(grad_tol)
.with_max_iterations(max_iter)
.with_initial_sample(seed.clone(), initial_op_sample)
.with_cg_tolerance(0.5)
.with_hessian_fallback_policy(HessianFallbackPolicy::Error);
if let Some(feedback) = config.outer_inner_cap.as_ref() {
solver = solver.with_observer(OuterAcceptObserver {
feedback: feedback.clone(),
});
}
if let Some(r) = sanitized_operator_trust_restart_radius(
config.operator_initial_trust_radius,
) {
solver = solver.with_initial_trust_radius(r);
}
let mf_start = std::time::Instant::now();
let report = solver.run_report();
let mf_elapsed = mf_start.elapsed().as_secs_f64();
let final_radius = report.diagnostics.final_trust_radius;
log::info!(
"[OUTER summary] matrix-free TR finished status={:?} in {} iters \
elapsed={:.3}s final_value={:.6e} final_trust_radius={}",
report.status,
report.solution.iterations,
mf_elapsed,
report.solution.final_value,
match final_radius {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
},
);
match report.status {
OptimizationStatus::Converged
| OptimizationStatus::NumericallyConverged => {
let mut result =
solution_into_outer_result(report.solution, true, *the_plan);
result.operator_stop_reason =
Some(OperatorTrustRegionStopReason::Converged);
result.operator_trust_radius = final_radius;
Ok(result)
}
OptimizationStatus::MaxIterations => {
log::warn!(
"[OUTER warning] {context}: matrix-free TR hit max_iter={} at final_value={:.6e} |g|={:.3e} final_trust_radius={}",
config.max_iter,
report.solution.final_value,
report.solution.final_gradient_norm.unwrap_or(f64::NAN),
match final_radius {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
},
);
let mut result =
solution_into_outer_result(report.solution, false, *the_plan);
result.operator_stop_reason =
Some(OperatorTrustRegionStopReason::IterationBudget);
result.operator_trust_radius = final_radius;
Ok(result)
}
OptimizationStatus::TrustRegionRejectFloor => {
log::warn!(
"[OUTER warning] {context}: matrix-free TR reached trust-radius reject floor at final_value={:.6e} |g|={:.3e} final_trust_radius={}",
report.solution.final_value,
report.solution.final_gradient_norm.unwrap_or(f64::NAN),
match final_radius {
Some(r) => format!("{:.3e}", r),
None => "n/a".to_string(),
},
);
let mut result =
solution_into_outer_result(report.solution, false, *the_plan);
result.operator_stop_reason =
Some(OperatorTrustRegionStopReason::RejectFloor);
result.operator_trust_radius = final_radius;
Ok(result)
}
OptimizationStatus::ObjectiveFailed
| OptimizationStatus::NumericalFailure
| OptimizationStatus::LineSearchFailed => {
Err(EstimationError::RemlOptimizationFailed(format!(
"matrix-free TR solver failed with status={:?}",
report.status
)))
}
}
} else {
let hessian_source = the_plan.hessian_source;
let (lo, hi) = &bounds_template;
let bounds = outer_bounds(lo, hi)?;
let grad_tol = outer_gradient_tolerance(config);
let max_iter = outer_max_iterations(config.max_iter)?;
let objective = OuterSecondOrderBridge {
obj,
layout,
hessian_source,
materialize_operator_max_dim: OUTER_HVP_MATERIALIZE_MAX_DIM,
eval_count: 0,
outer_inner_cap: config.outer_inner_cap.clone(),
g_norm_initial: None,
last_g_norm: None,
};
let seed_hessian = build_bridge_hessian_for_source(
hessian_source,
seed_eval.hessian.clone(),
OUTER_HVP_MATERIALIZE_MAX_DIM,
)
.map_err(|err| match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(message)
}
})?;
let initial_sample = SecondOrderSample {
value: seed_eval.cost,
gradient: seed_eval.gradient.clone(),
hessian: seed_hessian,
};
let mut optimizer = ArcOptimizer::new(seed.clone(), objective)
.with_bounds(bounds)
.with_gradient_tolerance(grad_tol)
.with_max_iterations(max_iter)
.with_initial_sample(seed.clone(), initial_sample);
if let Some(sigma) = config.arc_initial_regularization {
optimizer = optimizer.with_initial_regularization(sigma);
}
if let Some(feedback) = config.outer_inner_cap.as_ref() {
optimizer = optimizer.with_observer(OuterAcceptObserver {
feedback: feedback.clone(),
});
}
if matches!(hessian_source, HessianSource::Analytic) {
optimizer = optimizer
.with_hessian_fallback_policy(HessianFallbackPolicy::Error)
.with_fallback_policy(OptFallbackPolicy::Never);
}
match optimizer.run() {
Ok(sol) => Ok(solution_into_outer_result(sol, true, *the_plan)),
Err(ArcError::MaxIterationsReached { last_solution, .. }) => {
log::warn!(
"[OUTER warning] {context}: ARC hit max_iter={} at final_value={:.6e} |g|={:.3e}",
config.max_iter,
last_solution.final_value,
last_solution.final_gradient_norm.unwrap_or(f64::NAN),
);
Ok(solution_into_outer_result(*last_solution, false, *the_plan))
}
Err(e) => Err(EstimationError::RemlOptimizationFailed(format!(
"Arc solver failed: {e:?}"
))),
}
}
}
Solver::Bfgs => {
if cap.gradient != Derivative::Analytic {
return Err(EstimationError::RemlOptimizationFailed(format!(
"{context}: outer BFGS requires an analytic gradient capability; \
no non-analytic fallback is available (plan={the_plan}, \
declared gradient={:?})",
cap.gradient,
)));
}
let seed_eval = obj
.eval_with_order(seed, OuterEvalOrder::ValueAndGradient)
.map_err(|err| into_objective_error("outer eval failed", err));
let seed_eval = match seed_eval {
Ok(seed_eval) => seed_eval,
Err(err) => {
let err = match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(message)
}
};
log::warn!(
"[OUTER] {context}: rejecting seed {seed_idx} before solver start: {err}"
);
rejection_reasons.push((seed_idx, "validation", err.to_string()));
continue 'seed_attempts;
}
};
let seed_eval = match finite_outer_first_order_eval_or_error(
"outer eval failed",
layout,
seed_eval,
)
.map_err(|err| match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(message)
}
}) {
Ok(eval) => eval,
Err(err) => {
log::warn!(
"[OUTER] {context}: rejecting seed {seed_idx} before solver start: {err}"
);
rejection_reasons.push((seed_idx, "validation", err.to_string()));
continue 'seed_attempts;
}
};
started_seeds += 1;
seed_slot = started_seeds;
let (lo, hi) = &bounds_template;
let bounds = outer_bounds(lo, hi)?;
let grad_tol = outer_gradient_tolerance(config);
let max_iter = outer_max_iterations(config.max_iter)?;
let objective = OuterFirstOrderBridge {
obj,
layout,
outer_inner_cap: config.outer_inner_cap.clone(),
iter_count: 0,
g_norm_initial: None,
last_g_norm: None,
};
let initial_sample = FirstOrderSample {
value: seed_eval.cost,
gradient: seed_eval.gradient.clone(),
};
let mut optimizer = Bfgs::new(seed.clone(), objective)
.with_initial_sample(seed.clone(), initial_sample)
.with_bounds(bounds)
.with_gradient_tolerance(grad_tol)
.with_max_iterations(max_iter);
if let Some(caps) = bfgs_axis_step_caps(config, layout) {
optimizer = optimizer.with_axis_step_caps(caps);
}
if let Some(feedback) = config.outer_inner_cap.as_ref() {
optimizer = optimizer.with_observer(OuterAcceptObserver {
feedback: feedback.clone(),
});
}
let bfgs_start = std::time::Instant::now();
let outcome = optimizer.run();
let bfgs_elapsed = bfgs_start.elapsed().as_secs_f64();
match &outcome {
Ok(sol) => log::info!(
"[OUTER summary] BFGS converged in {} iters elapsed={:.3}s final_value={:.6e}",
sol.iterations,
bfgs_elapsed,
sol.final_value
),
Err(BfgsError::MaxIterationsReached { last_solution }) => log::warn!(
"[OUTER summary] BFGS hit max_iter in {} iters elapsed={:.3}s final_value={:.6e}",
last_solution.iterations,
bfgs_elapsed,
last_solution.final_value
),
Err(BfgsError::LineSearchFailed {
last_solution,
max_attempts,
failure_reason,
}) => log::info!(
"[OUTER summary] BFGS line-search failed in {} iters elapsed={:.3}s final_value={:.6e} reason={:?} max_attempts={} |g|={:.3e}",
last_solution.iterations,
bfgs_elapsed,
last_solution.final_value,
failure_reason,
max_attempts,
last_solution.final_gradient_norm.unwrap_or(f64::NAN),
),
Err(e) => log::info!(
"[OUTER summary] BFGS failed elapsed={:.3}s err={:?}",
bfgs_elapsed,
e
),
}
match outcome {
Ok(sol) => Ok(solution_into_outer_result(sol, true, *the_plan)),
Err(BfgsError::MaxIterationsReached { last_solution }) => {
Ok(solution_into_outer_result(*last_solution, false, *the_plan))
}
Err(BfgsError::LineSearchFailed {
last_solution,
max_attempts,
failure_reason,
}) => {
if last_solution.final_value.is_finite()
&& last_solution.final_point.iter().all(|v| v.is_finite())
&& last_solution
.final_gradient
.as_ref()
.is_none_or(|g| g.iter().all(|v| v.is_finite()))
{
Ok(solution_into_outer_result(*last_solution, false, *the_plan))
} else {
Err(EstimationError::RemlOptimizationFailed(
bfgs_line_search_failure_message(
context,
&last_solution,
max_attempts,
failure_reason,
),
))
}
}
Err(BfgsError::ObjectiveFailed { message }) => {
Err(EstimationError::RemlOptimizationFailed(format!(
"BFGS solver failed: ObjectiveFailed {{ message: {message:?} }}"
)))
}
Err(e) => Err(EstimationError::RemlOptimizationFailed(format!(
"BFGS solver failed: {e:?}"
))),
}
}
Solver::Efs => {
let mut objective = OuterFixedPointBridge {
obj,
layout,
barrier_config: cap.barrier_config.clone(),
fixed_point_tolerance: config.tolerance,
consecutive_psi_zero_iters: 0,
};
match objective.eval_step(seed) {
Ok(_) => {}
Err(err) => {
let err = match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(message)
}
};
if requests_immediate_first_order_fallback(&err.to_string()) {
return Err(err);
}
log::warn!(
"[OUTER] {context}: rejecting seed {seed_idx} before solver start: {err}"
);
rejection_reasons.push((seed_idx, "validation", err.to_string()));
continue 'seed_attempts;
}
};
started_seeds += 1;
seed_slot = started_seeds;
let (lo, hi) = &bounds_template;
let bounds = outer_bounds(lo, hi)?;
let tol = outer_tolerance(config.tolerance)?;
let max_iter = outer_max_iterations(config.max_iter)?;
let mut optimizer = FixedPoint::new(seed.clone(), objective)
.with_bounds(bounds)
.with_tolerance(tol)
.with_max_iterations(max_iter);
match optimizer.run() {
Ok(sol) => Ok(solution_into_outer_result(sol, true, *the_plan)),
Err(FixedPointError::MaxIterationsReached { last_solution }) => {
log::warn!(
"[OUTER warning] {context}: EFS hit max_iter={} at final_value={:.6e} step_norm={:.3e}",
config.max_iter,
last_solution.final_value,
last_solution.final_gradient_norm.unwrap_or(f64::NAN),
);
Ok(solution_into_outer_result(*last_solution, false, *the_plan))
}
Err(e) => Err(EstimationError::RemlOptimizationFailed(format!(
"fixed-point solver failed: {e:?}"
))),
}
}
Solver::HybridEfs => {
let mut objective = OuterFixedPointBridge {
obj,
layout,
barrier_config: cap.barrier_config.clone(),
fixed_point_tolerance: config.tolerance,
consecutive_psi_zero_iters: 0,
};
match objective.eval_step(seed) {
Ok(_) => {}
Err(err) => {
let err = match err {
ObjectiveEvalError::Recoverable { message }
| ObjectiveEvalError::Fatal { message } => {
EstimationError::RemlOptimizationFailed(message)
}
};
if requests_immediate_first_order_fallback(&err.to_string()) {
return Err(err);
}
log::warn!(
"[OUTER] {context}: rejecting seed {seed_idx} before solver start: {err}"
);
rejection_reasons.push((seed_idx, "validation", err.to_string()));
continue 'seed_attempts;
}
};
started_seeds += 1;
seed_slot = started_seeds;
let (lo, hi) = &bounds_template;
let bounds = outer_bounds(lo, hi)?;
let tol = outer_tolerance(config.tolerance)?;
let max_iter = outer_max_iterations(config.max_iter)?;
let mut optimizer = FixedPoint::new(seed.clone(), objective)
.with_bounds(bounds)
.with_tolerance(tol)
.with_max_iterations(max_iter);
match optimizer.run() {
Ok(sol) => Ok(solution_into_outer_result(sol, true, *the_plan)),
Err(FixedPointError::MaxIterationsReached { last_solution }) => {
log::warn!(
"[OUTER warning] {context}: HybridEFS hit max_iter={} at final_value={:.6e} step_norm={:.3e}",
config.max_iter,
last_solution.final_value,
last_solution.final_gradient_norm.unwrap_or(f64::NAN),
);
Ok(solution_into_outer_result(*last_solution, false, *the_plan))
}
Err(e) => Err(EstimationError::RemlOptimizationFailed(format!(
"hybrid EFS solver failed: {e:?}"
))),
}
}
Solver::CompassSearch => {
let projected_seed = project_to_bounds(seed, Some(&bounds_template));
let seed_cost = obj.eval_cost(&projected_seed).map_err(|err| {
EstimationError::RemlOptimizationFailed(format!(
"aux direct-search seed cost failed ({context}): {err}"
))
})?;
if !seed_cost.is_finite() {
rejection_reasons.push((
seed_idx,
"validation",
format!("aux direct-search rejects non-finite seed cost ({seed_cost})"),
));
continue 'seed_attempts;
}
started_seeds += 1;
seed_slot = started_seeds;
let (lo, hi) = &bounds_template;
let outcome = compass_search_outer(
obj,
projected_seed,
seed_cost,
lo.view(),
hi.view(),
1.0,
config.tolerance,
config.max_iter,
);
match outcome {
CompassSearchOutcome::Converged { point, cost, polls } => Ok(OuterResult {
rho: point,
final_value: cost,
iterations: polls,
final_grad_norm: None,
final_gradient: None,
final_hessian: None,
converged: true,
plan_used: *the_plan,
operator_trust_radius: None,
operator_stop_reason: None,
}),
CompassSearchOutcome::BudgetExhausted { point, cost, polls } => {
log::warn!(
"[OUTER warning] {context}: compass search exhausted max_polls={} at best_cost={:.6e}",
config.max_iter,
cost,
);
Ok(OuterResult {
rho: point,
final_value: cost,
iterations: polls,
final_grad_norm: None,
final_gradient: None,
final_hessian: None,
converged: false,
plan_used: *the_plan,
operator_trust_radius: None,
operator_stop_reason: None,
})
}
}
}
};
let seed_elapsed = t_seed_start.elapsed().as_secs_f64();
match result {
Ok(candidate) => {
let candidate_converged = candidate.converged;
log::debug!(
"[outer-timing] seed {}/{} ({:?}): {:.3}s cost={:.6e} converged={}",
seed_slot,
seed_budget,
the_plan.solver,
seed_elapsed,
candidate.final_value,
candidate.converged,
);
if candidate_improves_best(&candidate, best.as_ref()) {
best = Some(candidate);
}
let quality_compare_remaining_gaussian_seeds = matches!(
config.seed_config.risk_profile,
crate::seeding::SeedRiskProfile::Gaussian
) && seed_budget > 1
&& started_seeds < seed_budget;
if best.as_ref().is_some_and(|b| b.converged)
&& !quality_compare_remaining_gaussian_seeds
{
break;
}
if !candidate_converged && matches!(expensive_seed_limit, Some(limit) if limit > 0)
{
unsuccessful_expensive_seeds += 1;
if let Some(limit) = expensive_seed_limit
&& unsuccessful_expensive_seeds >= limit
{
log::info!(
"[OUTER] {context}: stopping expensive multi-start after {} non-converged {:?} seed(s)",
unsuccessful_expensive_seeds,
the_plan.solver,
);
stopped_early_due_to_limit = true;
break;
}
}
}
Err(e) => {
if requests_immediate_first_order_fallback(&e.to_string()) {
return Err(e);
}
log::debug!(
"[outer-timing] seed {}/{} ({:?}): {:.3}s FAILED: {}",
seed_slot,
seed_budget,
the_plan.solver,
seed_elapsed,
e,
);
rejection_reasons.push((seed_idx, "solver", e.to_string()));
if let Some(limit) = expensive_seed_limit {
unsuccessful_expensive_seeds += 1;
if unsuccessful_expensive_seeds >= limit {
log::info!(
"[OUTER] {context}: stopping expensive multi-start after {} failed {:?} seed(s)",
unsuccessful_expensive_seeds,
the_plan.solver,
);
stopped_early_due_to_limit = true;
break;
}
}
}
}
}
best.ok_or_else(|| {
while last_classified_reason_idx < rejection_reasons.len() {
let (idx, phase, msg) = &rejection_reasons[last_classified_reason_idx];
seed_rejections.push(SeedRejection::from_message(*idx, phase, msg.clone()));
last_classified_reason_idx += 1;
}
let n_generated = seeds.len();
let n_screened = n_generated;
let n_exact_validated = seed_rejections.len() + started_seeds;
let stats = StartupStats::from_rejections(
n_generated,
n_screened,
n_exact_validated,
started_seeds,
&seed_rejections,
);
let structural = structural_early_exit_key
.clone()
.or_else(|| uniform_structural_key(&seed_rejections, 1));
let early_exit_note = if structural_early_exit_key.is_some() {
"early-exit triggered: every observed seed reported the same structural CertRefused"
.to_string()
} else if stopped_early_due_to_limit {
format!(
"stopped early after {unsuccessful_expensive_seeds} consecutive non-converged \
{:?} seed(s) (expensive_unsuccessful_seed_limit)",
the_plan.solver
)
} else {
String::new()
};
if started_seeds == 0 {
EstimationError::RemlOptimizationFailed(format_no_seeds_passed(
context,
&stats,
&seed_rejections,
structural.as_ref(),
&early_exit_note,
))
} else {
let header = format!(
"all {started_seeds} seed candidates failed ({context}); \
generated={}, screened={}, exact_validated={}, solver_started={}",
stats.generated, stats.screened, stats.exact_validated, stats.solver_started,
);
let body = format_no_seeds_passed(
context,
&stats,
&seed_rejections,
structural.as_ref(),
&early_exit_note,
);
EstimationError::RemlOptimizationFailed(format!("{header}\n{body}"))
}
})
}
#[cfg(test)]
mod tests {
use super::*;
use ::opt::FixedPointObjective;
use ndarray::array;
use std::sync::atomic::{AtomicUsize, Ordering};
use std::sync::{Arc, Mutex};
struct FailingSeedMaterializationOperator {
dim: usize,
}
impl OuterHessianOperator for FailingSeedMaterializationOperator {
fn dim(&self) -> usize {
self.dim
}
fn matvec(&self, v: &Array1<f64>) -> Result<Array1<f64>, String> {
Ok(v.clone())
}
fn is_cheap_to_materialize(&self) -> bool {
true
}
fn materialize_dense(&self) -> Result<Array2<f64>, String> {
Err("seed materialization failed".to_string())
}
}
#[test]
fn materialize_dense_uses_single_batched_mul_mat() {
struct BatchedOnlyHessian {
matrix: Array2<f64>,
matvec_calls: Arc<AtomicUsize>,
mul_mat_calls: Arc<AtomicUsize>,
rhs_columns: Arc<AtomicUsize>,
}
impl OuterHessianOperator for BatchedOnlyHessian {
fn dim(&self) -> usize {
self.matrix.nrows()
}
fn matvec(&self, v: &Array1<f64>) -> Result<Array1<f64>, String> {
self.matvec_calls.fetch_add(1, Ordering::Relaxed);
Ok(self.matrix.dot(v))
}
fn mul_mat(&self, factor: ArrayView2<'_, f64>) -> Result<Array2<f64>, String> {
self.mul_mat_calls.fetch_add(1, Ordering::Relaxed);
self.rhs_columns
.fetch_add(factor.ncols(), Ordering::Relaxed);
Ok(self.matrix.dot(&factor))
}
}
let matvec_calls = Arc::new(AtomicUsize::new(0));
let mul_mat_calls = Arc::new(AtomicUsize::new(0));
let rhs_columns = Arc::new(AtomicUsize::new(0));
let op = BatchedOnlyHessian {
matrix: array![[2.0, 0.25, -0.5], [0.5, 3.0, 1.0], [-0.25, 2.0, 4.0]],
matvec_calls: Arc::clone(&matvec_calls),
mul_mat_calls: Arc::clone(&mul_mat_calls),
rhs_columns: Arc::clone(&rhs_columns),
};
let dense = op
.materialize_dense()
.expect("batched dense materialization");
let expected = array![[2.0, 0.375, -0.375], [0.375, 3.0, 1.5], [-0.375, 1.5, 4.0]];
assert_eq!(dense, expected);
assert_eq!(
mul_mat_calls.load(Ordering::Relaxed),
1,
"dense materialization must batch all identity columns into one mul_mat call"
);
assert_eq!(
rhs_columns.load(Ordering::Relaxed),
3,
"the single batched materialization call must include every identity RHS"
);
assert_eq!(
matvec_calls.load(Ordering::Relaxed),
0,
"operators with batched mul_mat must not be probed column-by-column"
);
}
#[test]
fn plan_analytic_hessian_selects_arc() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 3,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn plan_prefer_gradient_only_does_not_hide_analytic_hessian() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 3,
psi_dim: 1,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: true,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn plan_survival_baseline_exact_hessian_selects_arc() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 3,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn plan_no_hessian_few_params_selects_bfgs() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 3,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn plan_no_hessian_many_params_selects_bfgs() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 12,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn plan_boundary_8_params_uses_bfgs() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: SMALL_OUTER_BFGS_MAX_PARAMS,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn plan_boundary_9_params_uses_bfgs() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: SMALL_OUTER_BFGS_MAX_PARAMS + 1,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn plan_efs_selected_for_penalty_like_many_params() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Efs);
assert_eq!(p.hessian_source, HessianSource::EfsFixedPoint);
}
#[test]
fn plan_penalty_like_without_fixed_point_stays_bfgs() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn plan_efs_not_selected_few_params_even_if_penalty_like() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 5,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn plan_efs_not_selected_with_analytic_hessian() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 20,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
}
#[test]
fn plan_efs_with_no_gradient_penalty_like_many_params() {
let cap = OuterCapability {
gradient: Derivative::Unavailable,
hessian: DeclaredHessianForm::Unavailable,
n_params: 20,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Efs);
assert_eq!(p.hessian_source, HessianSource::EfsFixedPoint);
}
#[test]
fn plan_efs_allowed_with_barrier_config() {
let barrier = BarrierConfig {
tau: 1e-6,
constrained_indices: vec![0, 1],
lower_bounds: vec![0.0, 0.0],
bound_signs: vec![1.0, 1.0],
};
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 0,
fixed_point_available: true,
barrier_config: Some(barrier),
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Efs);
assert_eq!(p.hessian_source, HessianSource::EfsFixedPoint);
}
#[test]
fn plan_efs_allowed_with_barrier_config_no_gradient() {
let barrier = BarrierConfig {
tau: 1e-6,
constrained_indices: vec![0],
lower_bounds: vec![0.0],
bound_signs: vec![1.0],
};
let cap = OuterCapability {
gradient: Derivative::Unavailable,
hessian: DeclaredHessianForm::Unavailable,
n_params: 20,
psi_dim: 0,
fixed_point_available: true,
barrier_config: Some(barrier),
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Efs);
assert_eq!(p.hessian_source, HessianSource::EfsFixedPoint);
}
#[test]
fn barrier_curvature_significant_blocks_efs_at_runtime() {
let barrier = BarrierConfig {
tau: 1e-6,
constrained_indices: vec![0],
lower_bounds: vec![0.0],
bound_signs: vec![1.0],
};
let beta_near = Array1::from_vec(vec![0.001]);
assert!(barrier.barrier_curvature_is_significant(&beta_near, 1.0, 0.01));
let beta_far = Array1::from_vec(vec![10.0]);
assert!(!barrier.barrier_curvature_is_significant(&beta_far, 1.0, 0.01));
}
#[test]
fn barrier_curvature_locally_concentrated_covers_both_failure_modes() {
let barrier = BarrierConfig {
tau: 1e-6,
constrained_indices: vec![0, 1],
lower_bounds: vec![0.0, 0.0],
bound_signs: vec![1.0, 1.0],
};
let mild_uniform = Array1::from_vec(vec![1.0e-2, 1.0e-2]);
assert!(!barrier.barrier_curvature_locally_concentrated(&mild_uniform, 0.1, 1.0));
let tight_uniform = Array1::from_vec(vec![1.0e-4, 1.0e-4]);
assert!(barrier.barrier_curvature_locally_concentrated(&tight_uniform, 0.1, 1.0));
assert!(!barrier.barrier_curvature_locally_concentrated(&tight_uniform, 0.1, 1.0e9));
let large_uniform = Array1::from_vec(vec![10.0, 10.0]);
assert!(!barrier.barrier_curvature_locally_concentrated(&large_uniform, 0.1, 1.0));
let imbalanced = Array1::from_vec(vec![1.0e-2, 1.0]);
assert!(barrier.barrier_curvature_locally_concentrated(&imbalanced, 0.1, 1.0));
assert!(!barrier.barrier_curvature_locally_concentrated(&imbalanced, 1.0e-3, 1.0e9));
let infeasible = Array1::from_vec(vec![-0.5, 1.0]);
assert!(barrier.barrier_curvature_locally_concentrated(&infeasible, 0.1, 1.0));
}
#[test]
fn hessian_result_unwrap_analytic() {
let h = Array2::<f64>::eye(3);
let result = HessianResult::Analytic(h.clone());
assert!(result.is_analytic());
let extracted = result.unwrap_analytic();
assert_eq!(extracted, h);
}
#[test]
#[should_panic(expected = "expected analytic Hessian")]
fn hessian_result_unwrap_unavailable_panics() {
let result = HessianResult::Unavailable;
result.unwrap_analytic();
}
#[test]
fn zero_params_selects_arc() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 0,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn hessian_result_into_option() {
let h = Array2::<f64>::eye(2);
let result = HessianResult::Analytic(h.clone());
assert_eq!(result.into_option(), Some(h));
let result = HessianResult::Unavailable;
assert_eq!(result.into_option(), None);
}
#[test]
fn closure_objective_delegates() {
let mut obj = ClosureObjective {
state: 42_i32,
cap: OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 1,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
},
cost_fn: |_: &mut i32, _: &Array1<f64>| Ok(1.0),
eval_fn: |_: &mut i32, _: &Array1<f64>| {
Ok(OuterEval {
cost: 1.0,
gradient: Array1::zeros(1),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
eval_order_fn: None::<
fn(&mut i32, &Array1<f64>, OuterEvalOrder) -> Result<OuterEval, EstimationError>,
>,
reset_fn: Some(|st: &mut i32| {
*st = 42;
}),
efs_fn: None::<fn(&mut i32, &Array1<f64>) -> Result<EfsEval, EstimationError>>,
screening_proxy_fn: None::<fn(&mut i32, &Array1<f64>) -> Result<f64, EstimationError>>,
seed_fn: None::<fn(&mut i32, &Array1<f64>) -> Result<(), EstimationError>>,
};
assert_eq!(obj.capability().n_params, 1);
assert_eq!(obj.eval_cost(&Array1::zeros(1)).unwrap(), 1.0);
}
#[test]
fn closure_objective_seed_inner_state_delegates_when_hook_present() {
let mut obj = ClosureObjective {
state: Vec::<f64>::new(),
cap: OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 1,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
},
cost_fn: |_: &mut Vec<f64>, _: &Array1<f64>| Ok(0.0),
eval_fn: |_: &mut Vec<f64>, _: &Array1<f64>| {
Ok(OuterEval {
cost: 0.0,
gradient: Array1::zeros(1),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
eval_order_fn: None::<
fn(
&mut Vec<f64>,
&Array1<f64>,
OuterEvalOrder,
) -> Result<OuterEval, EstimationError>,
>,
reset_fn: None::<fn(&mut Vec<f64>)>,
efs_fn: None::<fn(&mut Vec<f64>, &Array1<f64>) -> Result<EfsEval, EstimationError>>,
screening_proxy_fn: None::<
fn(&mut Vec<f64>, &Array1<f64>) -> Result<f64, EstimationError>,
>,
seed_fn: None::<fn(&mut Vec<f64>, &Array1<f64>) -> Result<(), EstimationError>>,
}
.with_seed_inner_state(|state: &mut Vec<f64>, beta: &Array1<f64>| {
state.extend(beta.iter().copied());
Ok(())
});
obj.seed_inner_state(&array![1.5, -2.0]).unwrap();
assert_eq!(obj.state, vec![1.5, -2.0]);
}
#[test]
fn hybrid_efs_backtracking_uses_half_step_after_first_rejection() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 12,
psi_dim: 1,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let mut obj = ClosureObjective {
state: (),
cap: cap.clone(),
cost_fn: |_: &mut (), theta: &Array1<f64>| {
let psi = theta[11];
let cost = if (psi - 0.0).abs() < 1e-12 {
1.0
} else if (psi - 0.5).abs() < 1e-12 {
0.5
} else {
2.0
};
Ok(cost)
},
eval_fn: |_: &mut (), theta: &Array1<f64>| {
Ok(OuterEval {
cost: theta[11].abs(),
gradient: Array1::zeros(theta.len()),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
eval_order_fn: None::<
fn(&mut (), &Array1<f64>, OuterEvalOrder) -> Result<OuterEval, EstimationError>,
>,
reset_fn: None::<fn(&mut ())>,
efs_fn: Some(|_: &mut (), theta: &Array1<f64>| {
let mut steps = vec![0.0; theta.len()];
steps[11] = 1.0;
Ok(EfsEval {
cost: 1.0,
steps,
beta: None,
psi_gradient: Some(array![1.0]),
psi_indices: Some(vec![11]),
inner_hessian_scale: None,
})
}),
screening_proxy_fn: None::<fn(&mut (), &Array1<f64>) -> Result<f64, EstimationError>>,
seed_fn: None::<fn(&mut (), &Array1<f64>) -> Result<(), EstimationError>>,
};
let mut bridge = OuterFixedPointBridge {
obj: &mut obj,
layout: cap.theta_layout(),
barrier_config: None,
fixed_point_tolerance: 1e-8,
consecutive_psi_zero_iters: 0,
};
let sample = bridge
.eval_step(&Array1::zeros(cap.n_params))
.expect("hybrid EFS step should backtrack cleanly");
assert_eq!(sample.status, FixedPointStatus::Continue);
assert_eq!(sample.step.len(), cap.n_params);
assert_eq!(sample.step[11], 0.5);
assert!(
sample
.step
.iter()
.enumerate()
.all(|(idx, &value)| idx == 11 || value == 0.0)
);
}
#[test]
fn run_bfgs_mode_aware_eval_skips_hessian_work() {
let seen_orders = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_initial_rho(array![1.0])
.with_max_iter(1);
let mut obj = problem.build_objective_with_eval_order(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0] * theta[0]),
|_: &mut (), _: &Array1<f64>| {
Err(EstimationError::InvalidInput(
"legacy eager eval should not run on BFGS".to_string(),
))
},
{
let seen_orders = Arc::clone(&seen_orders);
move |_: &mut (), theta: &Array1<f64>, order: OuterEvalOrder| {
seen_orders.lock().unwrap().push(order);
Ok(OuterEval {
cost: theta[0] * theta[0],
gradient: array![2.0 * theta[0]],
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "mode-aware bfgs first order")
.expect("BFGS should use the order-aware first-order bridge");
assert_eq!(result.plan_used.solver, Solver::Bfgs);
let seen_orders = seen_orders.lock().unwrap();
assert!(
!seen_orders.is_empty(),
"mode-aware eval hook should have been used"
);
assert!(
seen_orders
.iter()
.all(|order| *order == OuterEvalOrder::ValueAndGradient),
"BFGS should request only value+gradient, saw {seen_orders:?}"
);
}
#[test]
fn first_order_bridge_keeps_true_gradient_on_repeated_flat_cost() {
let eval_calls = Arc::new(AtomicUsize::new(0));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable);
let mut obj = problem.build_objective(
(),
|_: &mut (), _: &Array1<f64>| Ok(1000.0),
{
let eval_calls = Arc::clone(&eval_calls);
move |_: &mut (), _: &Array1<f64>| {
let call = eval_calls.fetch_add(1, Ordering::Relaxed);
let cost = match call {
0 => 999.9995,
1 => 999.9990,
_ => 999.9987,
};
Ok(OuterEval {
cost,
gradient: array![4.0],
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let mut bridge = OuterFirstOrderBridge {
obj: &mut obj,
layout: OuterThetaLayout::new(1, 0),
outer_inner_cap: None,
iter_count: 0,
g_norm_initial: None,
last_g_norm: None,
};
let first = FirstOrderObjective::eval_grad(&mut bridge, &array![0.0])
.expect("first flat-cost eval should expose the true gradient");
let second = FirstOrderObjective::eval_grad(&mut bridge, &array![0.0])
.expect("second flat-cost eval should expose the true gradient");
let third = FirstOrderObjective::eval_grad(&mut bridge, &array![0.0])
.expect("third flat-cost eval should expose the true gradient");
let fourth = FirstOrderObjective::eval_grad(&mut bridge, &array![0.0])
.expect("fourth flat-cost eval should expose the true gradient");
assert_eq!(first.gradient[0], 4.0);
assert_eq!(second.gradient[0], 4.0);
assert_eq!(third.gradient[0], 4.0);
assert_eq!(fourth.gradient[0], 4.0);
assert_eq!(bridge.last_g_norm, Some(4.0));
assert_eq!(eval_calls.load(Ordering::Relaxed), 4);
}
#[test]
fn outer_second_order_bridge_separates_first_and_second_order_requests() {
let seen_orders = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either);
let mut obj = problem.build_objective_with_eval_order(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0] * theta[0]),
|_: &mut (), _: &Array1<f64>| {
Err(EstimationError::InvalidInput(
"legacy eager eval should not run".to_string(),
))
},
{
let seen_orders = Arc::clone(&seen_orders);
move |_: &mut (), theta: &Array1<f64>, order: OuterEvalOrder| {
seen_orders.lock().unwrap().push(order);
Ok(OuterEval {
cost: theta[0] * theta[0],
gradient: array![2.0 * theta[0]],
hessian: match order {
OuterEvalOrder::ValueAndGradient => HessianResult::Unavailable,
OuterEvalOrder::ValueGradientHessian => {
HessianResult::Analytic(array![[2.0]])
}
},
inner_beta_hint: None,
})
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let mut bridge = OuterSecondOrderBridge {
obj: &mut obj,
layout: OuterThetaLayout::new(1, 0),
hessian_source: HessianSource::Analytic,
materialize_operator_max_dim: OUTER_HVP_MATERIALIZE_MAX_DIM,
eval_count: 0,
outer_inner_cap: None,
g_norm_initial: None,
last_g_norm: None,
};
let grad_sample =
FirstOrderObjective::eval_grad(&mut bridge, &array![1.0]).expect("grad eval");
assert_eq!(grad_sample.value, 1.0);
assert_eq!(grad_sample.gradient, array![2.0]);
let hess_sample =
SecondOrderObjective::eval_hessian(&mut bridge, &array![1.0]).expect("hessian eval");
assert_eq!(hess_sample.value, 1.0);
assert_eq!(hess_sample.gradient, array![2.0]);
assert_eq!(hess_sample.hessian, Some(array![[2.0]]));
let seen_orders = seen_orders.lock().unwrap();
assert!(
*seen_orders
== vec![
OuterEvalOrder::ValueAndGradient,
OuterEvalOrder::ValueGradientHessian
],
"second-order bridge should split first-order and second-order requests, saw {seen_orders:?}"
);
}
#[test]
fn analytic_route_unavailable_hessian_is_fatal() {
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either);
let mut obj = problem.build_objective_with_eval_order(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0] * theta[0]),
|_: &mut (), _: &Array1<f64>| {
Err(EstimationError::InvalidInput(
"legacy eager eval should not run".to_string(),
))
},
move |_: &mut (), theta: &Array1<f64>, _order: OuterEvalOrder| {
Ok(OuterEval {
cost: theta[0] * theta[0],
gradient: array![2.0 * theta[0]],
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let mut bridge = OuterSecondOrderBridge {
obj: &mut obj,
layout: OuterThetaLayout::new(1, 0),
hessian_source: HessianSource::Analytic,
materialize_operator_max_dim: OUTER_HVP_MATERIALIZE_MAX_DIM,
eval_count: 0,
outer_inner_cap: None,
g_norm_initial: None,
last_g_norm: None,
};
let err = SecondOrderObjective::eval_hessian(&mut bridge, &array![1.0])
.expect_err("Analytic route must reject Unavailable Hessian, not pass None to opt");
match err {
ObjectiveEvalError::Fatal { message } => {
assert!(
message.contains("HessianSource::Analytic") && message.contains("Unavailable"),
"fatal message should explain the analytic-route mismatch, saw: {message}"
);
}
ObjectiveEvalError::Recoverable { message } => panic!(
"Analytic-route Hessian violations must be Fatal (FD estimation is forbidden); \
got Recoverable: {message}"
),
}
}
#[test]
fn outer_config_default() {
let cfg = OuterConfig::default();
assert_eq!(cfg.tolerance, 1e-5);
assert_eq!(cfg.max_iter, 200);
assert_eq!(cfg.rho_bound, 30.0);
}
#[test]
fn plan_hybrid_efs_selected_for_psi_coords_many_params() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 1,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::HybridEfs);
assert_eq!(p.hessian_source, HessianSource::HybridEfsFixedPoint);
}
#[test]
fn plan_psi_without_fixed_point_stays_bfgs() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 1,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn plan_hybrid_efs_no_gradient_selected_for_psi_coords() {
let cap = OuterCapability {
gradient: Derivative::Unavailable,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 1,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::HybridEfs);
assert_eq!(p.hessian_source, HessianSource::HybridEfsFixedPoint);
}
fn cap_for_routing(
gradient: Derivative,
hessian: DeclaredHessianForm,
n_params: usize,
) -> OuterCapability {
OuterCapability {
gradient,
hessian,
n_params,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
}
}
#[test]
fn routing_analytic_analytic_stays_arc_at_biobank_scale() {
let cap = cap_for_routing(Derivative::Analytic, DeclaredHessianForm::Either, 6);
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn routing_analytic_analytic_stays_arc_at_dense_work_scale() {
let cap = cap_for_routing(Derivative::Analytic, DeclaredHessianForm::Either, 3);
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn routing_unavailable_hessian_routes_to_bfgs() {
let cap = cap_for_routing(Derivative::Analytic, DeclaredHessianForm::Unavailable, 8);
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn routing_explicit_prefer_gradient_only_does_not_override_exact_hessian() {
let mut cap = cap_for_routing(Derivative::Analytic, DeclaredHessianForm::Either, 6);
cap.prefer_gradient_only = true;
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn routing_log_line_arc_analytic_does_not_advertise_matrix_free() {
let p = OuterPlan {
solver: Solver::Arc,
hessian_source: HessianSource::Analytic,
};
let line = p.routing_log_line();
assert!(line.contains("solver=Arc"), "got {line}");
assert!(line.contains("hessian=Analytic"), "got {line}");
assert!(line.contains("matrix-free=false"), "got {line}");
}
#[test]
fn routing_log_line_bfgs_reports_no_matrix_free() {
let p = OuterPlan {
solver: Solver::Bfgs,
hessian_source: HessianSource::BfgsApprox,
};
let line = p.routing_log_line();
assert!(line.contains("solver=Bfgs"), "got {line}");
assert!(line.contains("hessian=BfgsApprox"), "got {line}");
assert!(line.contains("matrix-free=false"), "got {line}");
}
#[test]
fn routing_log_line_efs_reports_no_matrix_free() {
for source in [
HessianSource::EfsFixedPoint,
HessianSource::HybridEfsFixedPoint,
] {
let p = OuterPlan {
solver: Solver::Efs,
hessian_source: source,
};
assert!(
p.routing_log_line().contains("matrix-free=false"),
"{:?} should not advertise matrix-free",
source
);
}
}
#[test]
fn routing_custom_family_gamlss_stays_on_arc_when_both_derivs_analytic() {
let cap = cap_for_routing(Derivative::Analytic, DeclaredHessianForm::Either, 4);
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn routing_matern_iso_kappa_stays_on_arc_when_both_derivs_analytic() {
let cap = cap_for_routing(Derivative::Analytic, DeclaredHessianForm::Either, 5);
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn routing_matern_iso_large_kappa_dim_stays_on_arc_with_analytic_hessian() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 37,
psi_dim: 31,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn routing_marginal_slope_stays_on_arc_when_both_derivs_analytic() {
let cap = cap_for_routing(Derivative::Analytic, DeclaredHessianForm::Either, 3);
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn plan_hybrid_efs_not_selected_few_params() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 5,
psi_dim: 1,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn plan_exact_hvp_capability_selects_arc_even_when_fixed_point_is_available() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 64,
psi_dim: 16,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: true,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
assert_eq!(p.hessian_source, HessianSource::Analytic);
}
#[test]
fn plan_hybrid_efs_not_selected_with_analytic_hessian() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 20,
psi_dim: 1,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Arc);
}
#[test]
fn plan_pure_efs_not_hybrid_when_all_penalty_like() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Efs);
assert_eq!(p.hessian_source, HessianSource::EfsFixedPoint);
}
#[test]
fn automatic_fallbacks_preserve_analytic_hessian_for_arc_primary() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 12,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
assert_eq!(plan(&cap).solver, Solver::Arc);
let attempts = automatic_fallback_attempts(&cap);
assert!(
attempts.is_empty(),
"ARC primary must not lateral-demote to BFGS+BfgsApprox; \
ARC budget retries live in the runner",
);
}
#[test]
fn automatic_fallbacks_from_efs_prefer_analytic_bfgs_over_fd() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
assert_eq!(plan(&cap).solver, Solver::Efs);
let attempts = automatic_fallback_attempts(&cap);
assert!(!attempts.is_empty(), "EFS failure must have a fallback");
assert_eq!(attempts[0].gradient, Derivative::Analytic);
assert_eq!(attempts[0].hessian, DeclaredHessianForm::Unavailable);
assert!(attempts[0].disable_fixed_point);
assert_eq!(plan(&attempts[0]).solver, Solver::Bfgs);
assert!(
attempts.iter().all(|c| c.gradient == Derivative::Analytic),
"fallback cascade must stay on analytic-gradient attempts",
);
}
#[test]
fn automatic_fallbacks_from_hybrid_efs_prefer_analytic_bfgs_over_fd() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 2,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
assert_eq!(plan(&cap).solver, Solver::HybridEfs);
let attempts = automatic_fallback_attempts(&cap);
assert!(!attempts.is_empty());
assert_eq!(attempts[0].gradient, Derivative::Analytic);
assert!(attempts[0].disable_fixed_point);
assert_eq!(plan(&attempts[0]).solver, Solver::Bfgs);
}
#[test]
fn disabled_fallback_hybrid_efs_capability_routes_to_bfgs_primary() {
let trapped_cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 9,
psi_dim: 6,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
assert_eq!(plan(&trapped_cap).solver, Solver::HybridEfs);
let disabled_config = OuterConfig {
fallback_policy: FallbackPolicy::Disabled,
..OuterConfig::default()
};
let primary_cap = primary_capability_for_config(
trapped_cap.clone(),
&disabled_config,
"biobank exact adaptive",
);
assert!(primary_cap.disable_fixed_point);
assert_eq!(plan(&primary_cap).solver, Solver::Bfgs);
let pure_efs_cap = OuterCapability {
psi_dim: 0,
..trapped_cap.clone()
};
assert_eq!(plan(&pure_efs_cap).solver, Solver::Efs);
let pure_primary_cap =
primary_capability_for_config(pure_efs_cap.clone(), &disabled_config, "pure EFS");
assert!(!pure_primary_cap.disable_fixed_point);
assert_eq!(plan(&pure_primary_cap).solver, Solver::Efs);
let no_gradient_cap = OuterCapability {
gradient: Derivative::Unavailable,
..trapped_cap.clone()
};
assert_eq!(plan(&no_gradient_cap).solver, Solver::HybridEfs);
let no_gradient_primary_cap = primary_capability_for_config(
no_gradient_cap.clone(),
&disabled_config,
"gradient-unavailable hybrid EFS",
);
assert!(!no_gradient_primary_cap.disable_fixed_point);
assert_eq!(plan(&no_gradient_primary_cap).solver, Solver::HybridEfs);
let automatic_config = OuterConfig::default();
let automatic_cap = primary_capability_for_config(
trapped_cap.clone(),
&automatic_config,
"biobank exact adaptive",
);
assert!(!automatic_cap.disable_fixed_point);
assert_eq!(plan(&automatic_cap).solver, Solver::HybridEfs);
let automatic_attempts = automatic_fallback_attempts(&trapped_cap);
assert!(!automatic_attempts.is_empty());
assert!(automatic_attempts[0].disable_fixed_point);
assert_eq!(plan(&automatic_attempts[0]).solver, Solver::Bfgs);
}
#[test]
fn disabled_fallback_hybrid_efs_problem_uses_bfgs_without_calling_efs() {
let efs_calls = Arc::new(AtomicUsize::new(0));
let problem = OuterProblem::new(9)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_psi_dim(6)
.with_fallback_policy(FallbackPolicy::Disabled)
.with_initial_rho(Array1::zeros(9))
.with_max_iter(5);
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(0.5 * theta.dot(theta)),
|_: &mut (), theta: &Array1<f64>| {
Ok(OuterEval {
cost: 0.5 * theta.dot(theta),
gradient: theta.clone(),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
{
let efs_calls = Arc::clone(&efs_calls);
Some(move |_: &mut (), _: &Array1<f64>| {
efs_calls.fetch_add(1, Ordering::Relaxed);
Err(EstimationError::RemlOptimizationFailed(format!(
"{} synthetic biobank adaptive HybridEFS escape",
EFS_FIRST_ORDER_FALLBACK_MARKER,
)))
})
},
);
let result = problem
.run(&mut obj, "disabled fallback marker")
.expect("disabled-fallback HybridEFS-shaped problem should route directly to BFGS");
assert_eq!(result.plan_used.solver, Solver::Bfgs);
assert_eq!(
efs_calls.load(Ordering::Relaxed),
0,
"central primary-capability canonicalization should avoid the EFS hook entirely"
);
}
#[test]
fn automatic_fallbacks_without_gradient_stop_at_fixed_point_status() {
for (psi_dim, expected_solver) in [(0, Solver::Efs), (2, Solver::HybridEfs)] {
let cap = OuterCapability {
gradient: Derivative::Unavailable,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
assert_eq!(plan(&cap).solver, expected_solver);
assert!(
automatic_fallback_attempts(&cap).is_empty(),
"gradient-unavailable fixed-point capabilities must not fabricate a BFGS fallback",
);
}
}
#[test]
fn automatic_fallbacks_do_not_repeat_arc_when_fixed_point_is_irrelevant() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 15,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
assert_eq!(plan(&cap).solver, Solver::Arc);
let attempts = automatic_fallback_attempts(&cap);
assert!(
attempts.is_empty(),
"ARC primary with incidental fixed_point_available must not \
cascade through the EFS arm or lateral-demote to BFGS",
);
}
#[test]
fn plan_disable_fixed_point_forces_bfgs_even_when_efs_eligible() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 15,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: true,
};
let p = plan(&cap);
assert_eq!(p.solver, Solver::Bfgs);
assert_eq!(p.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn run_malformed_gradient_seed_surfaces_as_error() {
let problem = OuterProblem::new(2)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_initial_rho(Array1::zeros(2))
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
|_: &mut (), _: &Array1<f64>| Ok(0.0),
|_: &mut (), _: &Array1<f64>| {
Ok(OuterEval {
cost: 0.0,
gradient: Array1::zeros(1),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let err = problem
.run(&mut obj, "test gradient mismatch")
.expect_err("malformed analytic gradient must surface as error");
assert!(
matches!(err, EstimationError::RemlOptimizationFailed(_)),
"unexpected error variant: {err:?}",
);
}
#[test]
fn run_bfgs_ignores_malformed_hessian_payload() {
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_initial_rho(array![0.0])
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0] * theta[0]),
|_: &mut (), theta: &Array1<f64>| {
Ok(OuterEval {
cost: theta[0] * theta[0],
gradient: array![2.0 * theta[0]],
hessian: HessianResult::Analytic(array![[f64::NAN, 0.0], [0.0, 1.0]]),
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "bfgs should ignore malformed hessian payload")
.expect("valid first-order data should be enough for BFGS");
assert_eq!(result.plan_used.solver, Solver::Bfgs);
assert_eq!(result.plan_used.hessian_source, HessianSource::BfgsApprox);
}
#[test]
fn finite_outer_eval_reports_gradient_length_mismatch() {
let err = finite_outer_eval_or_error(
"test gradient mismatch",
OuterThetaLayout::new(2, 0),
OuterEval {
cost: 0.0,
gradient: Array1::zeros(1),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
},
)
.expect_err("gradient mismatch should be rejected");
let message = match err {
ObjectiveEvalError::Recoverable { message } | ObjectiveEvalError::Fatal { message } => {
message
}
};
assert!(
message.contains("outer gradient length mismatch"),
"unexpected error: {message}"
);
}
#[test]
fn run_with_initial_seed_still_considers_generated_candidates() {
let generated = crate::seeding::generate_rho_candidates(
1,
None,
&crate::seeding::SeedConfig::default(),
);
let valid_seed = generated
.first()
.expect("seed generator should yield at least one candidate")
.clone();
let expected_seed = valid_seed.clone();
let initial_seed = array![9.0];
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 1;
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_seed_config(seed_config)
.with_initial_rho(initial_seed)
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
{
let valid_seed = valid_seed.clone();
move |_: &mut (), theta: &Array1<f64>| {
if theta == valid_seed {
Ok(0.0)
} else {
Ok(f64::INFINITY)
}
}
},
move |_: &mut (), theta: &Array1<f64>| {
if theta == valid_seed {
Ok(OuterEval {
cost: 0.0,
gradient: Array1::zeros(1),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
} else {
Ok(OuterEval::infeasible(theta.len()))
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "generated seed should remain reachable")
.expect("generated seed should still be eligible when an initial seed is provided");
assert_eq!(result.rho, expected_seed);
}
#[test]
fn run_indefinite_analytic_seed_stays_on_arc() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 1;
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either)
.with_seed_config(seed_config)
.with_initial_rho(array![0.0])
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0] * theta[0]),
|_: &mut (), _: &Array1<f64>| {
Ok(OuterEval {
cost: 0.0,
gradient: array![0.0],
hessian: HessianResult::Analytic(array![[-1.0]]),
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "indefinite seed geometry")
.expect("indefinite analytic seed geometry should stay on the second-order plan");
assert_eq!(result.plan_used.solver, Solver::Arc);
assert_eq!(result.plan_used.hessian_source, HessianSource::Analytic);
}
#[test]
fn run_seed_materialization_failure_surfaces_arc_error_verbatim() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 1;
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either)
.with_seed_config(seed_config)
.with_initial_rho(array![0.0])
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0] * theta[0]),
|_: &mut (), _: &Array1<f64>| {
Ok(OuterEval {
cost: 0.0,
gradient: array![0.0],
hessian: HessianResult::Operator(Arc::new(
FailingSeedMaterializationOperator { dim: 1 },
)),
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let err = problem
.run(&mut obj, "seed materialization failure")
.expect_err(
"ARC primary must surface the materialization failure verbatim — \
no lateral demote to BFGS+BfgsApprox",
);
let msg = err.to_string();
assert!(
msg.contains("seed materialization failed"),
"error must propagate the underlying materialization message; got: {msg}"
);
}
#[test]
fn run_nonconverged_arc_stays_on_arc_after_budget_retry_ladder() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 1;
let (_d, session) = tmp_cache_session("nonconverged-arc-cache");
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either)
.with_seed_config(seed_config)
.with_initial_rho(array![5.0])
.with_max_iter(1)
.with_cache_session(Arc::clone(&session));
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0].powi(4)),
|_: &mut (), theta: &Array1<f64>| {
let x = theta[0];
Ok(OuterEval {
cost: x.powi(4),
gradient: array![4.0 * x.powi(3)],
hessian: HessianResult::Analytic(array![[12.0 * x.powi(2)]]),
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "nonconverged arc should stay on arc")
.expect(
"ARC ladder must surface the last non-converged ARC result rather than \
demoting to BFGS+BfgsApprox",
);
assert_eq!(
result.plan_used.solver,
Solver::Arc,
"ARC primary must not lateral-demote after budget exhaustion"
);
assert_eq!(
result.plan_used.hessian_source,
HessianSource::Analytic,
"analytic outer Hessian must be preserved across the budget-bump retry ladder"
);
assert!(
!result.converged,
"test fixture is engineered so the ladder cannot converge; \
converged=true would mean the fixture stopped exercising the ladder"
);
}
#[test]
fn candidate_selection_prefers_lower_cost_within_same_convergence_class() {
let nonconverged_hi = OuterResult {
rho: array![0.0],
final_value: 9.0,
iterations: 1,
final_grad_norm: Some(1.0),
final_gradient: None,
final_hessian: None,
converged: false,
plan_used: OuterPlan {
solver: Solver::Bfgs,
hessian_source: HessianSource::BfgsApprox,
},
operator_trust_radius: None,
operator_stop_reason: None,
};
let nonconverged_lo = OuterResult {
rho: array![1.0],
final_value: 1.0,
iterations: 1,
final_grad_norm: Some(1.0),
final_gradient: None,
final_hessian: None,
converged: false,
plan_used: OuterPlan {
solver: Solver::Bfgs,
hessian_source: HessianSource::BfgsApprox,
},
operator_trust_radius: None,
operator_stop_reason: None,
};
let converged = OuterResult {
rho: array![2.0],
final_value: 5.0,
iterations: 1,
final_grad_norm: Some(0.0),
final_gradient: None,
final_hessian: None,
converged: true,
plan_used: OuterPlan {
solver: Solver::Bfgs,
hessian_source: HessianSource::BfgsApprox,
},
operator_trust_radius: None,
operator_stop_reason: None,
};
assert!(candidate_improves_best(&nonconverged_hi, None));
assert!(candidate_improves_best(
&nonconverged_lo,
Some(&nonconverged_hi)
));
assert!(!candidate_improves_best(
&nonconverged_hi,
Some(&nonconverged_lo)
));
assert!(candidate_improves_best(&converged, Some(&nonconverged_lo)));
assert!(!candidate_improves_best(&nonconverged_lo, Some(&converged)));
}
#[test]
fn gaussian_multistart_compares_converged_seed_costs() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 2;
seed_config.risk_profile = crate::seeding::SeedRiskProfile::Gaussian;
let started = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_seed_config(seed_config)
.with_max_iter(4);
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(if theta[0] < -1.0 { 0.0 } else { 10.0 }),
{
let started = Arc::clone(&started);
move |_: &mut (), theta: &Array1<f64>| {
started.lock().unwrap().push(theta.clone());
Ok(OuterEval {
cost: if theta[0] < -1.0 { 0.0 } else { 10.0 },
gradient: array![0.0],
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "Gaussian quality multistart")
.expect("Gaussian multistart should compare both converged seeds");
let starts = started.lock().unwrap();
assert!(
starts.len() >= 2,
"Gaussian quality mode should not stop at the first converged seed"
);
assert!(
result.rho[0] < -1.0,
"lower-cost converged Gaussian seed should win"
);
assert_eq!(result.final_value, 0.0);
}
#[test]
fn run_starts_solver_with_direct_startup_eval() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 1;
let calls = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either)
.with_seed_config(seed_config)
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
{
let calls = Arc::clone(&calls);
move |_: &mut (), theta: &Array1<f64>| {
calls.lock().unwrap().push("cost");
Ok(theta[0] * theta[0])
}
},
{
let calls = Arc::clone(&calls);
move |_: &mut (), theta: &Array1<f64>| {
calls.lock().unwrap().push("eval");
Ok(OuterEval {
cost: theta[0] * theta[0],
gradient: array![2.0 * theta[0]],
hessian: HessianResult::Analytic(array![[2.0]]),
inner_beta_hint: None,
})
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
problem
.run(&mut obj, "solver should start from a direct startup eval")
.expect("analytic plans should start with a direct full evaluation");
let calls = calls.lock().unwrap();
let first_eval_idx = calls
.iter()
.position(|call| *call == "eval")
.expect("solver should eventually request a full eval");
assert!(
first_eval_idx == 0,
"startup should not perform a separate cost-screening pass first: {calls:?}"
);
}
#[test]
fn run_screening_reorders_expensive_generated_seeds_before_full_startup_eval() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.max_seeds = 4;
seed_config.seed_budget = 2;
seed_config.risk_profile = crate::seeding::SeedRiskProfile::GeneralizedLinear;
let screening_cap = Arc::new(AtomicUsize::new(0));
let valid_seed = crate::seeding::generate_rho_candidates(1, None, &seed_config)
.last()
.expect("seed generator should yield at least one candidate")
.clone();
let started = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either)
.with_seed_config(seed_config)
.with_screening_cap(Arc::clone(&screening_cap))
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
{
let valid_seed = valid_seed.clone();
move |_: &mut (), theta: &Array1<f64>| {
if theta == valid_seed {
Ok(0.0)
} else {
Ok(1000.0)
}
}
},
{
let valid_seed = valid_seed.clone();
let started = Arc::clone(&started);
move |_: &mut (), theta: &Array1<f64>| {
started.lock().unwrap().push(theta.clone());
if theta == valid_seed {
Ok(OuterEval {
cost: 0.0,
gradient: array![0.0],
hessian: HessianResult::Analytic(array![[1.0]]),
inner_beta_hint: None,
})
} else {
Ok(OuterEval::infeasible(theta.len()))
}
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "screening should reorder expensive seeds")
.expect("screened startup should reach the best generated seed");
assert_eq!(result.rho, valid_seed);
assert_eq!(
started.lock().unwrap().first().cloned(),
Some(valid_seed),
"screening should move the lowest-cost seed to the front before full startup eval",
);
assert_eq!(screening_cap.load(std::sync::atomic::Ordering::Relaxed), 0);
}
#[test]
fn initial_rho_with_single_seed_budget_skips_expensive_screening() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.max_seeds = 4;
seed_config.seed_budget = 1;
seed_config.risk_profile = crate::seeding::SeedRiskProfile::GeneralizedLinear;
let screening_cap = Arc::new(AtomicUsize::new(0));
let screening_calls = Arc::new(AtomicUsize::new(0));
let initial_seed = array![9.0];
let started = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either)
.with_seed_config(seed_config)
.with_screening_cap(Arc::clone(&screening_cap))
.with_initial_rho(initial_seed.clone())
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
{
let screening_calls = Arc::clone(&screening_calls);
move |_: &mut (), _theta: &Array1<f64>| {
screening_calls.fetch_add(1, Ordering::Relaxed);
Ok(0.0)
}
},
{
let started = Arc::clone(&started);
let initial_seed = initial_seed.clone();
move |_: &mut (), theta: &Array1<f64>| {
started.lock().unwrap().push(theta.clone());
if theta == initial_seed {
Ok(OuterEval {
cost: 0.0,
gradient: array![0.0],
hessian: HessianResult::Analytic(array![[1.0]]),
inner_beta_hint: None,
})
} else {
Ok(OuterEval::infeasible(theta.len()))
}
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "initial rho should be authoritative")
.expect("initial-rho startup should not spend seed-screening solves");
assert_eq!(result.rho, initial_seed);
assert_eq!(
screening_calls.load(Ordering::Relaxed),
0,
"explicit initial rho plus seed_budget=1 should skip screening"
);
assert_eq!(
started.lock().unwrap().first().cloned(),
Some(initial_seed),
"solver should start from the explicit initial rho"
);
assert_eq!(screening_cap.load(Ordering::Relaxed), 0);
}
#[test]
fn run_screening_reorders_bfgs_seeds_before_full_startup_eval() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 1;
seed_config.risk_profile = crate::seeding::SeedRiskProfile::Gaussian;
let screening_cap = Arc::new(AtomicUsize::new(0));
let initial_seed = array![9.0];
let valid_seed = crate::seeding::generate_rho_candidates(1, None, &seed_config)
.first()
.expect("seed generator should yield at least one candidate")
.clone();
let started = Arc::new(Mutex::new(Vec::new()));
let screening_calls = Arc::new(AtomicUsize::new(0));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_seed_config(seed_config)
.with_screening_cap(Arc::clone(&screening_cap))
.with_initial_rho(initial_seed)
.with_screen_initial_rho(true)
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
{
let valid_seed = valid_seed.clone();
let screening_calls = Arc::clone(&screening_calls);
move |_: &mut (), theta: &Array1<f64>| {
screening_calls.fetch_add(1, Ordering::Relaxed);
if theta == valid_seed {
Ok(0.0)
} else {
Ok(1000.0)
}
}
},
{
let valid_seed = valid_seed.clone();
let started = Arc::clone(&started);
move |_: &mut (), theta: &Array1<f64>| {
started.lock().unwrap().push(theta.clone());
if theta == valid_seed {
Ok(OuterEval {
cost: 0.0,
gradient: array![0.0],
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
} else {
Ok(OuterEval::infeasible(theta.len()))
}
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "BFGS screening should reorder expensive seeds")
.expect("screened BFGS startup should reach the best generated seed");
assert_eq!(result.plan_used.solver, Solver::Bfgs);
assert_eq!(result.rho, valid_seed);
assert_eq!(
started.lock().unwrap().first().cloned(),
Some(valid_seed),
"BFGS screening should move the lowest-cost seed to the front before full startup eval",
);
assert!(
screening_calls.load(Ordering::Relaxed) > 1,
"BFGS seed screening should rank candidates with cost-only probes first",
);
assert_eq!(screening_cap.load(Ordering::Relaxed), 0);
}
#[test]
fn screening_cap_survives_per_seed_reset_before_proxy_eval() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.max_seeds = 3;
seed_config.seed_budget = 1;
seed_config.risk_profile = crate::seeding::SeedRiskProfile::Gaussian;
let screening_cap = Arc::new(AtomicUsize::new(0));
let proxy_saw_cap = Arc::new(AtomicBool::new(false));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_seed_config(seed_config)
.with_screening_cap(Arc::clone(&screening_cap))
.with_max_iter(1);
let mut obj = problem.build_objective_with_screening_proxy(
(),
|_: &mut (), _: &Array1<f64>| Ok(0.0),
|_: &mut (), theta: &Array1<f64>| {
Ok(OuterEval {
cost: theta[0].abs(),
gradient: array![0.0],
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
|_: &mut (), theta: &Array1<f64>, _: OuterEvalOrder| {
Ok(OuterEval {
cost: theta[0].abs(),
gradient: array![0.0],
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
{
let screening_cap = Arc::clone(&screening_cap);
Some(move |_: &mut ()| {
screening_cap.store(0, Ordering::Relaxed);
})
},
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
{
let screening_cap = Arc::clone(&screening_cap);
let proxy_saw_cap = Arc::clone(&proxy_saw_cap);
move |_: &mut (), theta: &Array1<f64>| {
let cap = screening_cap.load(Ordering::Relaxed);
if cap > 0 {
proxy_saw_cap.store(true, Ordering::Relaxed);
Ok(theta[0].abs())
} else {
Err(EstimationError::RemlOptimizationFailed(
"screening proxy ran without an active cap".to_string(),
))
}
}
},
);
problem
.run(&mut obj, "screening cap reset regression")
.expect("screening cap should be restored after each per-seed reset");
assert!(
proxy_saw_cap.load(Ordering::Relaxed),
"screening proxy should observe a nonzero cap"
);
assert_eq!(screening_cap.load(Ordering::Relaxed), 0);
}
#[test]
fn rank_seeds_cascade_escalates_when_initial_cap_collapses_all() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 1;
seed_config.screen_max_inner_iterations = 3;
let screening_cap = Arc::new(AtomicUsize::new(0));
let initial_seed = array![5.0];
let valid_seed = crate::seeding::generate_rho_candidates(1, None, &seed_config)
.first()
.expect("seed generator should yield at least one candidate")
.clone();
let max_cap_seen = Arc::new(AtomicUsize::new(0));
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either)
.with_seed_config(seed_config)
.with_screening_cap(Arc::clone(&screening_cap))
.with_initial_rho(initial_seed.clone())
.with_screen_initial_rho(true)
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
{
let screening_cap = Arc::clone(&screening_cap);
let max_cap_seen = Arc::clone(&max_cap_seen);
let valid_seed = valid_seed.clone();
move |_: &mut (), theta: &Array1<f64>| {
let cap = screening_cap.load(Ordering::Relaxed);
max_cap_seen.fetch_max(cap, Ordering::Relaxed);
if cap > 0 && cap < 12 {
return Ok(f64::NAN);
}
if theta == valid_seed {
Ok(0.0)
} else {
Ok(1000.0)
}
}
},
{
let valid_seed = valid_seed.clone();
move |_: &mut (), theta: &Array1<f64>| {
if theta == valid_seed {
Ok(OuterEval {
cost: 0.0,
gradient: array![0.0],
hessian: HessianResult::Analytic(array![[1.0]]),
inner_beta_hint: None,
})
} else {
Ok(OuterEval::infeasible(theta.len()))
}
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
problem
.run(&mut obj, "cascade should escalate")
.expect("cascade should reach a finite cost at the 4× cap stage");
let max_cap = max_cap_seen.load(Ordering::Relaxed);
assert_eq!(
max_cap, 12,
"cascade should stop at the 4× cap stage; observed max cap = {max_cap}"
);
assert_eq!(
screening_cap.load(Ordering::Relaxed),
0,
"screening cap must be restored to its previous value after cascade"
);
}
#[test]
fn run_efs_skips_global_cost_screening() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.max_seeds = 6;
seed_config.seed_budget = 1;
let screening_calls = Arc::new(AtomicUsize::new(0));
let problem = OuterProblem::new(15)
.with_gradient(Derivative::Unavailable)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_seed_config(seed_config)
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
{
let screening_calls = Arc::clone(&screening_calls);
move |_: &mut (), _: &Array1<f64>| {
screening_calls.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
Ok(0.0)
}
},
|_: &mut (), theta: &Array1<f64>| Ok(OuterEval::infeasible(theta.len())),
None::<fn(&mut ())>,
Some(|_: &mut (), theta: &Array1<f64>| {
Ok(EfsEval {
cost: 0.0,
steps: vec![0.0; theta.len()],
beta: None,
psi_gradient: None,
psi_indices: None,
inner_hessian_scale: None,
})
}),
);
problem
.run(
&mut obj,
"EFS should not use a separate global cost-screening pass",
)
.expect("first generated EFS seed should be sufficient");
assert_eq!(
screening_calls.load(std::sync::atomic::Ordering::Relaxed),
0,
"EFS startup should not call eval_cost just to screen seeds"
);
}
#[test]
fn run_efs_skips_invalid_leading_seed_without_spending_budget() {
let generated = crate::seeding::generate_rho_candidates(
15,
None,
&crate::seeding::SeedConfig::default(),
);
let valid_seed = generated
.first()
.expect("seed generator should yield at least one candidate")
.clone();
let invalid_seed = Array1::from_elem(15, 9.0);
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 1;
let problem = OuterProblem::new(15)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_seed_config(seed_config)
.with_initial_rho(invalid_seed)
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
|_: &mut (), _: &Array1<f64>| Ok(0.0),
|_: &mut (), theta: &Array1<f64>| Ok(OuterEval::infeasible(theta.len())),
None::<fn(&mut ())>,
{
let valid_seed = valid_seed.clone();
Some(move |_: &mut (), theta: &Array1<f64>| {
if theta == valid_seed {
Ok(EfsEval {
cost: 0.0,
steps: vec![0.0; theta.len()],
beta: None,
psi_gradient: None,
psi_indices: None,
inner_hessian_scale: None,
})
} else {
Err(EstimationError::RemlOptimizationFailed(
"invalid EFS seed".to_string(),
))
}
})
},
);
let result = problem
.run(&mut obj, "efs generated seed should remain reachable")
.expect("invalid startup seeds should not consume the only EFS seed slot");
assert_eq!(result.rho, valid_seed);
assert_eq!(result.plan_used.solver, Solver::Efs);
}
#[test]
fn run_efs_runtime_fallback_marker_degrades_to_bfgs_immediately() {
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.seed_budget = 2;
let efs_calls = Arc::new(AtomicUsize::new(0));
let problem = OuterProblem::new(12)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_seed_config(seed_config)
.with_initial_rho(Array1::zeros(12))
.with_max_iter(5);
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(0.5 * theta.dot(theta)),
|_: &mut (), theta: &Array1<f64>| {
Ok(OuterEval {
cost: 0.5 * theta.dot(theta),
gradient: theta.clone(),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
{
let efs_calls = Arc::clone(&efs_calls);
Some(move |_: &mut (), _: &Array1<f64>| {
efs_calls.fetch_add(1, std::sync::atomic::Ordering::Relaxed);
Err(EstimationError::RemlOptimizationFailed(format!(
"{} synthetic runtime escape hatch",
EFS_FIRST_ORDER_FALLBACK_MARKER,
)))
})
},
);
let result = problem
.run(&mut obj, "efs runtime fallback marker")
.expect("runtime EFS escape hatch should degrade to BFGS");
assert_eq!(result.plan_used.solver, Solver::Bfgs);
assert_eq!(
efs_calls.load(std::sync::atomic::Ordering::Relaxed),
1,
"runtime fallback marker should abort the EFS attempt immediately"
);
}
#[test]
fn run_rejects_invalid_theta_layout() {
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_psi_dim(2)
.with_initial_rho(Array1::zeros(1))
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
|_: &mut (), _: &Array1<f64>| Ok(0.0),
|_: &mut (), _: &Array1<f64>| {
Ok(OuterEval {
cost: 0.0,
gradient: Array1::zeros(1),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let err = problem
.run(&mut obj, "test invalid layout")
.expect_err("invalid theta layout should fail cleanly");
assert!(
err.to_string().contains("invalid outer theta layout"),
"unexpected error: {err}"
);
}
#[test]
fn effective_seed_budget_caps_expensive_solver_retries() {
assert_eq!(
effective_seed_budget(
4,
Solver::Efs,
crate::seeding::SeedRiskProfile::GeneralizedLinear,
false,
),
1
);
assert_eq!(
effective_seed_budget(
4,
Solver::HybridEfs,
crate::seeding::SeedRiskProfile::Survival,
false,
),
1
);
assert_eq!(
effective_seed_budget(
3,
Solver::Arc,
crate::seeding::SeedRiskProfile::GeneralizedLinear,
true,
),
1
);
assert_eq!(
effective_seed_budget(
3,
Solver::Arc,
crate::seeding::SeedRiskProfile::Survival,
false,
),
1
);
assert_eq!(
effective_seed_budget(
3,
Solver::Bfgs,
crate::seeding::SeedRiskProfile::Survival,
false,
),
3
);
}
fn aux_cap_unavailable(n_params: usize) -> OuterCapability {
OuterCapability {
gradient: Derivative::Unavailable,
hessian: DeclaredHessianForm::Unavailable,
n_params,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
}
}
#[test]
fn plan_with_class_primary_is_identical_to_plan_for_unavailable_grad() {
let cap = aux_cap_unavailable(3);
assert_eq!(plan_with_class(&cap, SolverClass::Primary), plan(&cap));
}
#[test]
fn plan_with_class_aux_unavailable_routes_to_compass_search() {
let cap = aux_cap_unavailable(3);
let p = plan_with_class(&cap, SolverClass::AuxiliaryGradientFree);
assert_eq!(p.solver, Solver::CompassSearch);
}
#[test]
fn plan_with_class_aux_analytic_grad_defers_to_primary_plan() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 3,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan_with_class(&cap, SolverClass::AuxiliaryGradientFree);
assert_eq!(p.solver, Solver::Arc);
}
#[test]
fn plan_with_class_aux_efs_eligible_defers_to_primary() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Unavailable,
n_params: 12,
psi_dim: 0,
fixed_point_available: true,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let p = plan_with_class(&cap, SolverClass::AuxiliaryGradientFree);
assert_eq!(p.solver, Solver::Efs);
}
#[test]
fn automatic_fallback_never_includes_compass_search() {
let cap = OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Either,
n_params: 5,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
};
let attempts = automatic_fallback_attempts(&cap);
for attempt_cap in &attempts {
let p = plan_with_class(attempt_cap, SolverClass::Primary);
assert_ne!(p.solver, Solver::CompassSearch);
}
}
#[test]
fn compass_search_budget_accounts_for_single_seed() {
let b = effective_seed_budget(
8,
Solver::CompassSearch,
crate::seeding::SeedRiskProfile::Survival,
false,
);
assert_eq!(b, 1);
}
#[test]
fn run_aux_compass_projects_seed_before_seed_cost() {
let seen = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_solver_class(SolverClass::AuxiliaryGradientFree)
.with_bounds(array![0.0], array![1.0])
.with_initial_rho(array![2.0])
.with_max_iter(64);
let mut obj = problem.build_objective(
(),
{
let seen = Arc::clone(&seen);
move |_: &mut (), theta: &Array1<f64>| {
seen.lock().unwrap().push(theta.clone());
Ok((theta[0] - 2.0).powi(2))
}
},
|_: &mut (), _: &Array1<f64>| {
Err(EstimationError::InvalidInput(
"aux direct-search test should not call eval".to_string(),
))
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let result = problem
.run(&mut obj, "aux direct-search seed projection")
.expect("aux direct-search should evaluate the projected seed");
assert_eq!(result.plan_used.solver, Solver::CompassSearch);
assert_eq!(result.rho, array![1.0]);
assert_eq!(result.final_value, 1.0);
assert_eq!(
seen.lock().unwrap().first().cloned(),
Some(array![1.0]),
"aux direct-search must project the seed before evaluating its cost",
);
}
#[test]
fn run_arc_projects_seed_before_seed_validation_eval() {
let seen = Arc::new(Mutex::new(Vec::new()));
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.max_seeds = 1;
seed_config.seed_budget = 1;
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Either)
.with_bounds(array![0.0], array![1.0])
.with_initial_rho(array![2.0])
.with_seed_config(seed_config)
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok((theta[0] - 0.25).powi(2)),
{
let seen = Arc::clone(&seen);
move |_: &mut (), theta: &Array1<f64>| {
seen.lock().unwrap().push(theta.clone());
Ok(OuterEval {
cost: (theta[0] - 0.25).powi(2),
gradient: array![2.0 * (theta[0] - 0.25)],
hessian: HessianResult::Analytic(array![[2.0]]),
inner_beta_hint: None,
})
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
problem
.run(&mut obj, "arc seed projection")
.expect("arc should evaluate the projected seed");
assert_eq!(
seen.lock().unwrap().first().cloned(),
Some(array![1.0]),
"Arc must project the seed before validating the initial sample",
);
}
#[test]
fn run_bfgs_projects_seed_before_seed_validation_eval() {
let seen = Arc::new(Mutex::new(Vec::new()));
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.max_seeds = 1;
seed_config.seed_budget = 1;
let problem = OuterProblem::new(1)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_bounds(array![0.0], array![1.0])
.with_initial_rho(array![2.0])
.with_seed_config(seed_config)
.with_max_iter(1);
let mut obj = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok((theta[0] - 0.25).powi(2)),
{
let seen = Arc::clone(&seen);
move |_: &mut (), theta: &Array1<f64>| {
seen.lock().unwrap().push(theta.clone());
Ok(OuterEval {
cost: (theta[0] - 0.25).powi(2),
gradient: array![2.0 * (theta[0] - 0.25)],
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
}
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
problem
.run(&mut obj, "bfgs seed projection")
.expect("BFGS should evaluate the projected seed");
assert_eq!(
seen.lock().unwrap().first().cloned(),
Some(array![1.0]),
"BFGS must project the seed before validating the initial sample",
);
}
fn tmp_cache_session(label: &str) -> (tempfile::TempDir, Arc<CacheSession>) {
let dir = tempfile::tempdir().unwrap();
let store = crate::cache::WarmStartStore::open(
dir.path().to_path_buf(),
crate::cache::StoreOptions {
size_budget_bytes: 1024 * 1024,
ttl: std::time::Duration::from_secs(60),
},
)
.unwrap();
let mut fp = crate::cache::Fingerprinter::new();
fp.absorb_str(b"outer-test", label);
let key = fp.finalize();
(dir, Arc::new(CacheSession::open(store, key)))
}
#[test]
fn checkpointing_objective_persists_finite_evals() {
let (_d, session) = tmp_cache_session("ckpt-persist");
let problem = OuterProblem::new(1).with_gradient(Derivative::Unavailable);
let mut inner: ClosureObjective<_, _, _> = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0] * theta[0]),
|_: &mut (), _: &Array1<f64>| {
Err(EstimationError::InvalidInput("eval not used".into()))
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let mut wrapped = CheckpointingObjective::new(&mut inner, Arc::clone(&session), Vec::new());
assert!(session.try_load().is_none());
let v0 = wrapped.eval_cost(&array![3.0]).unwrap();
assert!((v0 - 9.0).abs() < 1e-12);
let on_disk = session.try_load().expect("first eval should checkpoint");
let payload = decode_iterate(&on_disk.payload, 1).expect("payload decodes");
assert!((payload.cost - 9.0).abs() < 1e-12);
assert_eq!(payload.rho, vec![3.0]);
let v1 = wrapped.eval_cost(&array![0.5]).unwrap();
assert!((v1 - 0.25).abs() < 1e-12);
let on_disk = session
.try_load()
.expect("improving eval should checkpoint");
let payload = decode_iterate(&on_disk.payload, 1).expect("payload decodes");
assert!((payload.cost - 0.25).abs() < 1e-12);
assert_eq!(payload.rho, vec![0.5]);
let v_inf = wrapped.eval_cost(&array![f64::NAN]);
match v_inf {
Ok(value) => assert!(!value.is_finite()),
Err(err) => assert!(!err.to_string().is_empty()),
}
let on_disk = session.try_load().expect("prior best preserved");
let payload = decode_iterate(&on_disk.payload, 1).expect("payload decodes");
assert!((payload.cost - 0.25).abs() < 1e-12);
}
#[test]
fn checkpointing_objective_rejects_wrong_dim_on_decode() {
let bytes = encode_iterate(&array![1.0, 2.0, 3.0], None, 0.5, 0).expect("encode");
assert!(decode_iterate(&bytes, 3).is_some());
assert!(decode_iterate(&bytes, 5).is_none());
}
#[test]
fn iterate_payload_round_trips_beta() {
let rho = array![10.0, -10.0, 5.0];
let beta = array![0.12, -0.34, 0.56, 7.89];
let bytes = encode_iterate(&rho, Some(&beta), 1.0, 7).expect("encode");
let decoded = decode_iterate(&bytes, rho.len()).expect("decode");
assert_eq!(decoded.rho, rho.to_vec());
assert_eq!(decoded.beta, beta.to_vec());
let ro_bytes = encode_iterate(&rho, None, 1.0, 7).expect("encode-rho-only");
let ro = decode_iterate(&ro_bytes, rho.len()).expect("decode-rho-only");
assert!(ro.beta.is_empty());
}
#[test]
fn note_persists_inner_beta_hint_from_eval() {
let (_d, session) = tmp_cache_session("note-persists-beta");
let problem = OuterProblem::new(1).with_gradient(Derivative::Unavailable);
let mut inner: ClosureObjective<_, _, _> = problem.build_objective(
(),
|_: &mut (), _: &Array1<f64>| Ok(1.0),
|_: &mut (), theta: &Array1<f64>| {
Ok(OuterEval {
cost: theta[0] * theta[0],
gradient: array![2.0 * theta[0]],
hessian: HessianResult::Unavailable,
inner_beta_hint: Some(array![1.5, 2.5, 3.5]),
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let mut wrapped = CheckpointingObjective::new(&mut inner, Arc::clone(&session), Vec::new());
let eval = wrapped.eval(&array![0.5]).expect("eval ok");
assert!((eval.cost - 0.25).abs() < 1e-12);
let on_disk = session
.try_load()
.expect("eval with finite β must persist a (ρ,β) checkpoint");
let payload = decode_iterate(&on_disk.payload, 1).expect("payload decodes");
assert_eq!(payload.beta, vec![1.5, 2.5, 3.5]);
let captured = wrapped.last_inner_beta().expect("β was captured");
assert_eq!(captured.to_vec(), vec![1.5, 2.5, 3.5]);
}
#[test]
fn note_rejects_nonfinite_inner_beta() {
let (_d, session) = tmp_cache_session("note-rejects-bad-beta");
let problem = OuterProblem::new(1).with_gradient(Derivative::Unavailable);
let mut inner: ClosureObjective<_, _, _> = problem.build_objective(
(),
|_: &mut (), _: &Array1<f64>| Ok(1.0),
|_: &mut (), theta: &Array1<f64>| {
Ok(OuterEval {
cost: theta[0] * theta[0],
gradient: array![2.0 * theta[0]],
hessian: HessianResult::Unavailable,
inner_beta_hint: Some(array![f64::NAN, 0.5]),
})
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let mut wrapped = CheckpointingObjective::new(&mut inner, Arc::clone(&session), Vec::new());
let eval = wrapped.eval(&array![0.5]).expect("eval ok");
assert!((eval.cost - 0.25).abs() < 1e-12);
assert!(
session.try_load().is_none(),
"non-finite β must abort the checkpoint write, not poison the cache",
);
assert!(
wrapped.last_inner_beta().is_none(),
"non-finite β must not be exposed via last_inner_beta()",
);
}
#[test]
fn classify_extracts_beta_from_v2_payload() {
let rho = array![1.0, 2.0];
let beta = array![10.0, 20.0, 30.0];
let payload = encode_iterate(&rho, Some(&beta), 1.0, 0).expect("encode");
let loaded = crate::cache::LoadedEntry {
entry: crate::cache::CachedEntry {
payload,
objective: Some(1.0),
iteration: Some(0),
kind: crate::cache::EntryKind::Checkpoint,
written_unix_secs: 0,
},
source: crate::cache::LoadSource::Preloaded,
};
let CacheSeedDecision::Seed {
beta: decoded_beta, ..
} = classify_cache_entry_for_outer(&loaded, 2)
else {
panic!("expected Seed decision");
};
assert_eq!(decoded_beta, beta.to_vec());
let payload = encode_iterate(&rho, None, 1.0, 0).expect("encode");
let loaded = crate::cache::LoadedEntry {
entry: crate::cache::CachedEntry {
payload,
objective: Some(1.0),
iteration: Some(0),
kind: crate::cache::EntryKind::Checkpoint,
written_unix_secs: 0,
},
source: crate::cache::LoadSource::Preloaded,
};
let CacheSeedDecision::Seed {
beta: decoded_beta, ..
} = classify_cache_entry_for_outer(&loaded, 2)
else {
panic!("expected Seed decision");
};
assert!(
decoded_beta.is_empty(),
"ρ-only payload must produce an empty beta so the dispatcher skips seed_inner_state"
);
}
#[test]
fn run_calls_seed_inner_state_with_cached_beta() {
struct RecordingObj {
seeded: Arc<Mutex<Option<Array1<f64>>>>,
eval_count: Arc<Mutex<usize>>,
}
impl OuterObjective for RecordingObj {
fn capability(&self) -> OuterCapability {
OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Dense,
n_params: 2,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
}
}
fn eval_cost(&mut self, theta: &Array1<f64>) -> Result<f64, EstimationError> {
Ok(theta.dot(theta))
}
fn eval(&mut self, theta: &Array1<f64>) -> Result<OuterEval, EstimationError> {
*self.eval_count.lock().unwrap() += 1;
Ok(OuterEval {
cost: theta.dot(theta),
gradient: 2.0 * theta,
hessian: HessianResult::Analytic(2.0 * Array2::<f64>::eye(theta.len())),
inner_beta_hint: None,
})
}
fn reset(&mut self) {}
fn seed_inner_state(&mut self, beta: &Array1<f64>) -> Result<(), EstimationError> {
*self.seeded.lock().unwrap() = Some(beta.clone());
Ok(())
}
}
let (_d, session) = tmp_cache_session("seed-inner-state-call");
let bytes = encode_iterate(&array![1.0, 2.0], Some(&array![7.5, 8.5, 9.5]), 5.0, 3)
.expect("encode");
session.checkpoint(&bytes, Some(5.0), Some(3));
let seeded: Arc<Mutex<Option<Array1<f64>>>> = Arc::new(Mutex::new(None));
let eval_count: Arc<Mutex<usize>> = Arc::new(Mutex::new(0));
let mut obj = RecordingObj {
seeded: Arc::clone(&seeded),
eval_count: Arc::clone(&eval_count),
};
let problem = OuterProblem::new(2)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_max_iter(1)
.with_cache_session(Arc::clone(&session));
match problem.run(&mut obj, "seed-inner-state-call") {
Ok(result) => assert!(result.final_value.is_finite()),
Err(err) => assert!(!err.to_string().is_empty()),
}
let observed = seeded.lock().unwrap().clone();
assert_eq!(
observed,
Some(array![7.5, 8.5, 9.5]),
"dispatcher must call seed_inner_state with the cached β before run_outer",
);
}
#[test]
fn run_skips_seed_inner_state_when_payload_has_no_beta() {
struct CountingObj {
seed_calls: Arc<Mutex<usize>>,
}
impl OuterObjective for CountingObj {
fn capability(&self) -> OuterCapability {
OuterCapability {
gradient: Derivative::Analytic,
hessian: DeclaredHessianForm::Dense,
n_params: 2,
psi_dim: 0,
fixed_point_available: false,
barrier_config: None,
prefer_gradient_only: false,
disable_fixed_point: false,
}
}
fn eval_cost(&mut self, theta: &Array1<f64>) -> Result<f64, EstimationError> {
Ok(theta.dot(theta))
}
fn eval(&mut self, theta: &Array1<f64>) -> Result<OuterEval, EstimationError> {
Ok(OuterEval {
cost: theta.dot(theta),
gradient: 2.0 * theta,
hessian: HessianResult::Analytic(2.0 * Array2::<f64>::eye(theta.len())),
inner_beta_hint: None,
})
}
fn reset(&mut self) {}
fn seed_inner_state(&mut self, beta: &Array1<f64>) -> Result<(), EstimationError> {
*self.seed_calls.lock().unwrap() += beta.len().max(1);
Ok(())
}
}
let (_d, session) = tmp_cache_session("seed-inner-state-skip");
let bytes = encode_iterate(&array![1.0, 2.0], None, 5.0, 3).expect("encode");
session.checkpoint(&bytes, Some(5.0), Some(3));
let seed_calls: Arc<Mutex<usize>> = Arc::new(Mutex::new(0));
let mut obj = CountingObj {
seed_calls: Arc::clone(&seed_calls),
};
let problem = OuterProblem::new(2)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_max_iter(1)
.with_cache_session(Arc::clone(&session));
match problem.run(&mut obj, "seed-inner-state-skip") {
Ok(result) => assert!(result.final_value.is_finite()),
Err(err) => assert!(!err.to_string().is_empty()),
}
assert_eq!(
*seed_calls.lock().unwrap(),
0,
"seed_inner_state must not fire when the cached payload carries no β",
);
}
#[test]
fn cache_entry_classifier_honors_finite_seeds_regardless_of_saturation() {
for rho_seed in [array![9.0, 0.0], array![10.0, -10.0], array![-10.0, 10.0]] {
let payload = encode_iterate(&rho_seed, None, 1.0, 0).expect("encode");
let loaded = crate::cache::LoadedEntry {
entry: crate::cache::CachedEntry {
payload,
objective: Some(1.0),
iteration: Some(0),
kind: crate::cache::EntryKind::Checkpoint,
written_unix_secs: 0,
},
source: crate::cache::LoadSource::Preloaded,
};
assert!(cache_entry_would_help_outer(&loaded, 2));
let CacheSeedDecision::Seed { rho, .. } = classify_cache_entry_for_outer(&loaded, 2)
else {
panic!(
"finite seed {:?} must be honored unchanged; the read-side clamp / \
all-saturated-discard branches were band-aids over the missing β cache",
rho_seed
);
};
assert_eq!(rho, rho_seed, "ρ must round-trip without reshaping");
}
}
#[test]
fn cache_entry_classifier_rejects_only_structural_failures() {
let payload = encode_iterate(&array![0.5, 0.5], None, 1.0, 0).expect("encode");
let loaded = crate::cache::LoadedEntry {
entry: crate::cache::CachedEntry {
payload,
objective: Some(f64::NAN),
iteration: Some(0),
kind: crate::cache::EntryKind::Checkpoint,
written_unix_secs: 0,
},
source: crate::cache::LoadSource::Preloaded,
};
assert!(matches!(
classify_cache_entry_for_outer(&loaded, 2),
CacheSeedDecision::Discard {
reason: "non-finite-payload",
..
}
));
let payload = encode_iterate(&array![0.5, 0.5], None, 1.0, 0).expect("encode");
let loaded = crate::cache::LoadedEntry {
entry: crate::cache::CachedEntry {
payload,
objective: Some(1.0),
iteration: Some(0),
kind: crate::cache::EntryKind::Checkpoint,
written_unix_secs: 0,
},
source: crate::cache::LoadSource::Preloaded,
};
assert!(matches!(
classify_cache_entry_for_outer(&loaded, 3),
CacheSeedDecision::Discard {
reason: "payload-shape-mismatch",
..
}
));
}
#[test]
fn exact_final_cache_hit_is_helpful_even_at_boundary() {
let payload = encode_iterate(&array![10.0, -10.0], None, 1.0, 3).expect("encode");
let loaded = crate::cache::LoadedEntry {
entry: crate::cache::CachedEntry {
payload,
objective: Some(1.0),
iteration: Some(3),
kind: crate::cache::EntryKind::Final,
written_unix_secs: 0,
},
source: crate::cache::LoadSource::Exact,
};
assert!(cache_entry_would_help_outer(&loaded, 2));
assert!(matches!(
classify_cache_entry_for_outer(&loaded, 2),
CacheSeedDecision::ExactFinal { iterations: 3, .. }
));
}
#[test]
fn checkpointing_objective_mirrors_checkpoints() {
let (_primary_dir, primary) = tmp_cache_session("ckpt-primary");
let (_mirror_dir, mirror) = tmp_cache_session("ckpt-mirror");
let problem = OuterProblem::new(1).with_gradient(Derivative::Unavailable);
let mut inner: ClosureObjective<_, _, _> = problem.build_objective(
(),
|_: &mut (), theta: &Array1<f64>| Ok(theta[0] * theta[0]),
|_: &mut (), _: &Array1<f64>| {
Err(EstimationError::InvalidInput("eval not used".into()))
},
None::<fn(&mut ())>,
None::<fn(&mut (), &Array1<f64>) -> Result<EfsEval, EstimationError>>,
);
let mut wrapped = CheckpointingObjective::new(
&mut inner,
Arc::clone(&primary),
vec![Arc::clone(&mirror)],
);
let value = wrapped.eval_cost(&array![4.0]).unwrap();
assert_eq!(value, 16.0);
let primary_payload =
decode_iterate(&primary.try_load().expect("primary checkpoint").payload, 1)
.expect("primary decode");
let mirror_payload =
decode_iterate(&mirror.try_load().expect("mirror checkpoint").payload, 1)
.expect("mirror decode");
assert_eq!(primary_payload.rho, vec![4.0]);
assert_eq!(mirror_payload.rho, vec![4.0]);
assert_eq!(primary_payload.cost, mirror_payload.cost);
}
#[test]
fn cached_rho_is_prepended_as_first_seed() {
let (_d, session) = tmp_cache_session("seed-prepend");
let payload = encode_iterate(&array![2.5], None, 0.25, 0).expect("encode");
session.checkpoint(&payload, Some(0.25), Some(0));
assert!(
session.try_load().is_some(),
"precondition: cache populated"
);
let seen: Arc<Mutex<Vec<Array1<f64>>>> = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_solver_class(SolverClass::AuxiliaryGradientFree)
.with_bounds(array![-5.0], array![5.0])
.with_initial_rho(array![-3.0]) .with_max_iter(8)
.with_cache_session(Arc::clone(&session));
let mut obj = problem.build_objective(
seen.clone(),
|seen: &mut Arc<Mutex<Vec<Array1<f64>>>>, theta: &Array1<f64>| {
seen.lock().unwrap().push(theta.clone());
Ok((theta[0] - 2.5).powi(2))
},
|_: &mut Arc<Mutex<Vec<Array1<f64>>>>, _: &Array1<f64>| {
Err(EstimationError::InvalidInput("eval not used".into()))
},
None::<fn(&mut Arc<Mutex<Vec<Array1<f64>>>>)>,
None::<
fn(
&mut Arc<Mutex<Vec<Array1<f64>>>>,
&Array1<f64>,
) -> Result<EfsEval, EstimationError>,
>,
);
match problem.run(&mut obj, "seed-prepend") {
Ok(result) => assert!(result.final_value.is_finite()),
Err(err) => assert!(!err.to_string().is_empty()),
}
let evals = seen.lock().unwrap();
let pos_cached = evals.iter().position(|r| (r[0] - 2.5).abs() < 1e-9);
let pos_initial = evals.iter().position(|r| (r[0] + 3.0).abs() < 1e-9);
assert!(
pos_cached.is_some(),
"cached rho must be evaluated; saw {:?}",
*evals
);
if let (Some(c), Some(i)) = (pos_cached, pos_initial) {
assert!(
c <= i,
"cached rho (idx {c}) must precede initial_rho (idx {i})",
);
}
}
#[test]
fn all_saturated_cached_rho_is_honored_as_seed() {
let (_d, session) = tmp_cache_session("all-saturated-honored");
let payload = encode_iterate(&array![10.0, -10.0], None, 1.0, 0).expect("encode");
session.checkpoint(&payload, Some(1.0), Some(0));
assert!(
session.try_load().is_some(),
"precondition: cache populated"
);
let seen: Arc<Mutex<Vec<Array1<f64>>>> = Arc::new(Mutex::new(Vec::new()));
let mut seed_config = crate::seeding::SeedConfig::default();
seed_config.max_seeds = 4;
seed_config.seed_budget = 1;
let problem = OuterProblem::new(2)
.with_gradient(Derivative::Analytic)
.with_hessian(DeclaredHessianForm::Unavailable)
.with_seed_config(seed_config)
.with_initial_rho(array![0.0, 0.0])
.with_rho_bound(10.0)
.with_max_iter(1)
.with_cache_session(Arc::clone(&session));
let mut obj = problem.build_objective(
seen.clone(),
|_: &mut Arc<Mutex<Vec<Array1<f64>>>>, theta: &Array1<f64>| Ok(theta.dot(theta)),
|seen: &mut Arc<Mutex<Vec<Array1<f64>>>>, theta: &Array1<f64>| {
seen.lock().unwrap().push(theta.clone());
Ok(OuterEval {
cost: theta.dot(theta),
gradient: theta.clone(),
hessian: HessianResult::Unavailable,
inner_beta_hint: None,
})
},
None::<fn(&mut Arc<Mutex<Vec<Array1<f64>>>>)>,
None::<
fn(
&mut Arc<Mutex<Vec<Array1<f64>>>>,
&Array1<f64>,
) -> Result<EfsEval, EstimationError>,
>,
);
match problem.run(&mut obj, "all-saturated-honored") {
Ok(result) => assert!(result.final_value.is_finite()),
Err(err) => assert!(!err.to_string().is_empty()),
}
let evals = seen.lock().unwrap();
assert!(
evals.iter().any(|rho| rho == array![10.0, -10.0]),
"cached saturated ρ must be evaluated unchanged under v2 (ρ, β) invariant; saw {:?}",
*evals
);
}
#[test]
fn exact_final_cache_hit_skips_outer_validation() {
let (_d, session) = tmp_cache_session("final-skip");
let payload = encode_iterate(&array![2.5], None, 0.25, 7).expect("encode");
session.finalize(&payload, Some(0.25), Some(7));
let seen: Arc<Mutex<Vec<Array1<f64>>>> = Arc::new(Mutex::new(Vec::new()));
let problem = OuterProblem::new(1)
.with_solver_class(SolverClass::AuxiliaryGradientFree)
.with_bounds(array![-5.0], array![5.0])
.with_initial_rho(array![-3.0])
.with_max_iter(8)
.with_cache_session(Arc::clone(&session));
let mut obj = problem.build_objective(
seen.clone(),
|seen: &mut Arc<Mutex<Vec<Array1<f64>>>>, theta: &Array1<f64>| {
seen.lock().unwrap().push(theta.clone());
Ok((theta[0] - 2.5).powi(2))
},
|_: &mut Arc<Mutex<Vec<Array1<f64>>>>, _: &Array1<f64>| {
Err(EstimationError::InvalidInput("eval not used".into()))
},
None::<fn(&mut Arc<Mutex<Vec<Array1<f64>>>>)>,
None::<
fn(
&mut Arc<Mutex<Vec<Array1<f64>>>>,
&Array1<f64>,
) -> Result<EfsEval, EstimationError>,
>,
);
let result = problem
.run(&mut obj, "final-skip")
.expect("final exact hit should return cached outer result");
assert_eq!(result.rho, array![2.5]);
assert_eq!(result.final_value, 0.25);
assert_eq!(result.iterations, 7);
assert!(result.converged);
assert!(
seen.lock().unwrap().is_empty(),
"exact final hit should not evaluate the outer objective"
);
}
}